clear
clc
a = 0.95
k = [5;10;13;4;3;11;13;10;8;16;7;4];
k = -k;	% 模擬退火演算法是求解最小值，故取負數
d = [2;5;18;3;2;5;10;4;11;7;14;6];
restriction = 46;
num = 12;
sol_new = ones(1,num);         % 生成初始解
E_current = inf;E_best = inf;  
% E_current是當前解對應的目標函式值（即揹包中物品總價值）；
% E_new是新解的目標函式值；
% E_best是最優解的
sol_current = sol_new; sol_best = sol_new;
t0=97; tf=3; t=t0;
p=1;

while t>=tf
	for r=1:100
		%產生隨機擾動
		tmp=ceil(rand.*num);
		sol_new(1,tmp)=~sol_new(1,tmp);
		
		%檢查是否滿足約束
		while 1
			q=(sol_new*d <= restriction);
			if ~q
                p=~p;	%實現交錯著逆轉頭尾的第一個1
                tmp=find(sol_new==1);
                if p
                    sol_new(1,tmp)=0;
                else
                    sol_new(1,tmp(end))=0;
                end
            else
                break
			end
		end
		
		% 計算揹包中的物品價值
		E_new=sol_new*k;
		if E_new<E_current
            E_current=E_new;
            sol_current=sol_new;
            if E_new<E_best
				% 把冷卻過程中最好的解儲存下來
                E_best=E_new;
                sol_best=sol_new;
            end
		else
            if rand<exp(-(E_new-E_current)./t)
                E_current=E_new;
                sol_current=sol_new;
            else
                sol_new=sol_current;
            end
		end
	end
	t=t.*a;
end

disp('最優解為：')
sol_best
disp('物品總價值等於：')
val=-E_best;
disp(val)
disp('揹包中物品重量是：')
disp(sol_best * d)