%多處理器調度仿真程序%%%%%%%%%%%%%%%%%%%%%%
%           圖10.4 17任務                %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function msp()
clc;
N=17; %任務數
P=3;  %處理器個數
M=6; %時間段數
%連接矩陣賦值%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
C0=zeros(N,N);
C0(1,4)=1;C0(1,5)=1;C0(1,6)=1;
C0(2,7)=1;C0(3,11)=1;
C0(4,8)=1;
C0(5,8)=1;
C0(6,9)=1;
C0(7,10)=1;
C0(8,12)=1;C0(8,13)=1;
C0(9,14)=1;
C0(10,15)=1;
C0(11,16)=1;C0(11,17)=1;
%參數賦值%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%A=7.5,F=3,C=3,D=1,E=1%%
A=7.5;
B=0;
C=3;
D=1;
E=1; 
F=3; %拉格朗日參數
T0=9;  %初始溫度
beta=0.9; %溫度迭代參數
delta_1=10^(-5);
num=3000;
diedai=50000;
t=1;%能量函數自變量
V=0.9*rand(N,M);%初始化
%主降溫程序%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for x=1:num
    %迭代%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    for y=1:diedai
        tempV=V;
        V_1=NewV(A,B,C,D,E,F,P,N,M,T0,C0,V);
                 energy1(t)=energy(A,B,C,D,E,F,P,N,M,C0,V_1);
                 ene=energy1(t) %能量函數
        V=V_1;
        if sum(sum(abs(tempV-V_1),1))<delta_1
            break
        end
                 t=t+1;
    end
    T0=beta*T0; %降溫
    [ch,c3,VF]=check(P,C0,V);
    if ch==0
        if c3==0
            disp('ok')
            zlV=VF
            for a=1:M
                tt=1;
                for i=1:N
                    if zlV(i,a)==1
                        mspf(a,tt)=i;
                        tt=tt+1;
                    end
                end  %調度結果輸出
            end
        end
        diaodu=mspf
        break
    end
end
save mspf
if x>num
    disp('miss')
%     [ch,c3,VF]=check(P,C0,V)
end
 plot(energy1) %能量變化繪制


%更新函數%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function New_V=NewV(A,B,C,D,E,F,P,N,M,T0,C0,V)
% i=unidrnd(N,1,1);
a=unidrnd(M,1,1); %隨機選擇更新神經元
for i=1:N
%     for a=1:M
        D1=2*sum(V(i,:))-2*V(i,a);
        temp=sum(V(:,a));
%        D2=2*(temp-P)*temp+(temp-P)^2;
        if a==1
            D3=sum(C0(:,i));
        else
            tc=C0(:,i)';
            tv=V(:,1:a-1);
            ta=(1-sum(tv,2)).^2;
            D3=tc*ta;
        end
        D4=2*(sum(V(i,:))-1);
        D5=1-2*V(i,a);
        D6=2*sum(sum(V(:,a))-P);
        d_V=A*D1+C*D3+D*D4+E*D5+F*D6;
        V(i,a)=0.5*(1+tanh(-d_V/(2*T0)));%神經元更新
%     end
end
for i=1:N
    for a=1:M
        if V(i,a)<0.01
            V(i,a)=0;
        else
            if V(i,a)>0.99
                V(i,a)=1;
            end
        end
    end
end     
New_V=V;

%能量函數%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function ET=energy(A,B,C,D,E,F,P,N,M,C0,V)
E1=sumsqr(sum(V,2))-sumsqr(V);
%E2=(sum(V,1)-P).^2*sum(V,1)';
E3=0;
for i=2:N
    E3=E3+sum(C0(1:i-1,i))*V(i,1);
end
for a=2:M
    for i=2:N
        E3=E3+V(i,a)*( C0(1:i-1,i)'*(1-sum(V(1:i-1,1:a-1),2)).^2 );
    end
end
E4=sum((1-sum(V,2)).^2);
E5=sum(sum(V,1))-sumsqr(V);
E6=sumsqr(sum(V,1)-P);
ET=A*E1+C*E3+D*E4+E*E5+F*E6;

%約束檢驗函數%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [check,c3,VF]=check(P,C0,V)
[N,M]=size(V);
for i=1:N
    for a=1:M
        if V(i,a)<0.7
            V(i,a)=0;
        else 
            if V(i,a)>0.3
                V(i,a)=1;
            end
        end
    end
end
VF=V;
check1=sumsqr(sum(V,2)-1);%每個任務在一個時間段執行
c2=max(sum(V,1));
if c2>P
    check2=1;
else
    check2=0;
end  %每個時間段任務數不超過P
E3=0;
for i=2:N
    E3=E3+sum(C0(:,i))*V(i,1);
end
for a=2:M
    for i=2:N
        E3=E3+V(i,a)*( C0(:,i)'*(1-sum(V(:,1:a-1),2)).^2 );
    end
end
check3=N-sum(sum(V,1)); %總完成任務數
check=check1+check2+check3;
c3=E3;%調度順序檢驗