clear
clc
close all
global M w c n
%w=[92,4,43,83,84,68,92,82,6,44,32,18,56,83,25,96,70,48,14,58];c=[44,46,90,72,91,40,75,35,8,54,78,40,77,15,61,17,75,29,75,63];M=878
%w=[95,4,60,32,23,72,80,62,65,46];c=[55,10,47,5,4,50,8,61,85,87];M=269
c=[220,208,198,192,180,180,165,162,160,158,155,130,125,122,120,118,115,110,105,101,100,100,98,96,95,90,88,82,80,77,75,73,72,70,69,66,...
    65,63,60,58,56,50,30,20,15,10,8,5,3,1];
w=[80,82,85,70,72,70,66,50,55,25,50,55,40,48,50,32,22,60,30,32,40,38,35,32,25,28,30,22,50,30,45,30,60,50,20,65,20,25,30,10,...
    20,25,15,10,10,10,4,4,2,1];
M=1000;
n=length(w);
x=round(rand(1,n));
x=correct(x);
f=x*c';
m=w*x';
T=500;%T0=15比較好 
T0=T;
a=0.95;
L=40*n;
e=0.001;
N=1;
X1(N,1:n+1)=[x f];
F1(N,1)=f;
str1='程序運行正在進行,請稍等......';
disp(str1)
tic
cnt=0;
while T>e
       for k=1:L
        i=round(rand*(n-1))+1;
        if x(i)==0
            if m+w(i)<=M
                x(i)=1;
                f=f+c(i);
                m=m+w(i);
            else
                j=round(rand*(n-1))+1;
                while x(j)~=1
                    j=round(rand*(n-1))+1;
                end
                df=c(i)-c(j);
                dm=w(i)-w(j);
                if m+dm<=M
                    if df>0 | exp(df/T)>rand
                        x(i)=1;
                        x(j)=0;
                        f=f+df;
                        m=m+dm;
                    end
                end
            end
        else
            j=round(rand*(n-1))+1;
            while x(j)~=0
                j=round(rand*(n-1))+1;
            end
            df=c(j)-c(i);
            dm=w(j)-w(i);
            if m+dm<=M
                if df>0 | exp(df/T)>rand
                    x(i)=0;
                    x(j)=1;
                    f=f+df;
                    m=m+dm;
                end
            end
        end
       end%與for循環對應
    N=N+1;
    X1(N,1:n)=x;    F1(N,1)=x*c';    wt(N,1)=x*w';
    T=T*a;
end
[fmax,indx]=max(F1);
x=X1(indx,1:n)
plot(F1)
xlabel('退火次數 N')
ylabel('目標函數值 f')
toc
[wt F1];
m=x*w'
fmax