clear
clc
%生成初始解
sol_new1=0;%（1）解空间（初始解）
sol_new2=0;
sol_new3=0;
sol_new4=0;
sol_new5=0;
sol_new6=0;
sol_new7=0;
sol_new8=0;
sol_new9=0;
sol_new10=0;
sol_new11=0;
sol_new12=0;
sol_new13=0;
sol_new14=0;
sol_new15=0;
sol_new16=0;
sol_current1 = sol_new1;
sol_best1 = sol_new1;
sol_current2 = sol_new2;
sol_best2 = sol_new2;
sol_current3 = sol_new3;
sol_best3 = sol_new3;
sol_current4 = sol_new4;
sol_best4 = sol_new4;
sol_current5 = sol_new5;
sol_best5 = sol_new5;
sol_current6 = sol_new6;
sol_best6 = sol_new6;
sol_current7 = sol_new7;
sol_best7 = sol_new7;
sol_current8 = sol_new8;
sol_best8 = sol_new8;
sol_current9 = sol_new9;
sol_best9 = sol_new9;
sol_current10 = sol_new10;
sol_best10 = sol_new10;
sol_current11 = sol_new11;
sol_best11 = sol_new11;
sol_current12 = sol_new12;
sol_best12 = sol_new12;
sol_current13 = sol_new13;
sol_best13 = sol_new13;
sol_current14 = sol_new14;
sol_best14 = sol_new14;
sol_current15 = sol_new15;
sol_best15 = sol_new15;
sol_current16 = sol_new16;
sol_best16 = sol_new16;
E_current = inf;
E_best = inf;
rand('state',sum(clock)); %初始化随机数发生器
t=90; %初始温度
tf=89.9; %结束温度
a = 0.99; %温度下降比例
while t>=tf%（7）结束条件
for r=1:1000 %退火次数
%产生随机扰动（3）新解的产生
sol_new1=sol_new1+rand*0.2;
sol_new2=sol_new2+rand*0.2;
sol_new3=sol_new3+rand*0.2;
sol_new4=sol_new4+rand*0.2;
sol_new5=sol_new5+rand*0.2;
sol_new6=sol_new6+rand*0.2;
sol_new7=sol_new7+rand*0.2;
sol_new8=sol_new8+rand*0.2;
sol_new9=sol_new9+rand*0.2;
sol_new10=sol_new10+rand*0.2;
sol_new11=sol_new11+rand*0.2;
sol_new12=sol_new12+rand*0.2;
sol_new13=sol_new13+rand*0.2;
sol_new14=sol_new14+rand*0.2;
sol_new15=sol_new15+rand*0.2;
sol_new16=sol_new16+rand*0.2;
%检查是否满足约束
if 0.040208*sol_new4+0.040208*sol_new5+0.104208*sol_new7+1.430208*sol_new9-0.20979*sol_new10-0.20979*sol_new11+95.73021*sol_new12+1.430208*sol_new13+1.430208*sol_new14+29.73021*sol_new15+22.29944*sol_new16<=8043.321537 && -0.062691*sol_new4-0.062691*sol_new5-0.126691*sol_new7-1.452691*sol_new9+0.18731*sol_new10+0.18731*sol_new11-95.75269*sol_new12-1.452691*sol_new13-1.452691*sol_new14+29.75269*sol_new15+22.32192*sol_new16<=6422.337753 && 28.35672*sol_new8+64.75672*sol_new9+64.75672*sol_new13+64.75672*sol_new14>=118394.2 && 28.37923*sol_new8+64.77923*sol_new9+64.77923*sol_new13+64.77923*sol_new14<=116771.2 && sol_new1>=0 && sol_new2>=0 && sol_new3>=0 && sol_new4>=0 && sol_new5>=0 && sol_new6>=0 && sol_new7>=0 && sol_new8>=0 && sol_new9>=0 && sol_new10>=0 && sol_new11>=0 && sol_new12>=0 && sol_new13>=0 && sol_new14>=0 && sol_new15>=0 && sol_new16>=0
else
sol_new1=sol_new1+rand*0.2;
sol_new2=sol_new2+rand*0.2;
sol_new3=sol_new3+rand*0.2;
sol_new4=sol_new4+rand*0.2;
sol_new5=sol_new5+rand*0.2;
sol_new6=sol_new6+rand*0.2;
sol_new7=sol_new7+rand*0.2;
sol_new8=sol_new8+rand*0.2;
sol_new9=sol_new9+rand*0.2;
sol_new10=sol_new10+rand*0.2;
sol_new11=sol_new11+rand*0.2;
sol_new12=sol_new12+rand*0.2;
sol_new13=sol_new13+rand*0.2;
sol_new14=sol_new14+rand*0.2;
sol_new15=sol_new15+rand*0.2;
sol_new16=sol_new16+rand*0.2;
continue;
end
%退火过程
E_new=350000*sol_new1+6500*sol_new2+350000*sol_new3+205000*sol_new4+205000*sol_new5+11800+sol_new6+1000*sol_new7+8500*sol_new8+7600*sol_new9+6000*sol_new11+4600*sol_new12+8150*sol_new13+6150*sol_new14+6100*sol_new15+40000*sol_new16;%（2）目标函数
if E_new<E_current%（5）接受准则
E_current=E_new;
sol_current1=sol_new1;
sol_current2=sol_new2;
sol_current3=sol_new3;
sol_current4=sol_new4;
sol_current5=sol_new5;
sol_current6=sol_new6;
sol_current7=sol_new7;
sol_current8=sol_new8;
sol_current9=sol_new9;
sol_current10=sol_new10;
sol_current11=sol_new11;
sol_current12=sol_new12;
sol_current13=sol_new13;
sol_current14=sol_new14;
sol_current15=sol_new15;
sol_current16=sol_new16;
if E_new<E_best
%把冷却过程中最好的解保存下来
E_best=E_new;
sol_best1=sol_new1;
sol_best2=sol_new2;
sol_best3=sol_new3;
sol_best4=sol_new4;
sol_best5=sol_new5;
sol_best6=sol_new6;
sol_best7=sol_new7;
sol_best8=sol_new8;
sol_best9=sol_new9;
sol_best10=sol_new10;
sol_best11=sol_new11;
sol_best12=sol_new12;
sol_best13=sol_new13;
sol_best14=sol_new14;
sol_best15=sol_new15;
sol_best16=sol_new16;
end
else
if rand<exp(-(E_new-E_current)/t)%（4）代价函数差
E_current=E_new;
sol_current1=sol_new1;
sol_current2=sol_new2;
sol_current3=sol_new3;
sol_current4=sol_new4;
sol_current5=sol_new5;
sol_current6=sol_new6;
sol_current7=sol_new7;
sol_current8=sol_new8;
sol_current9=sol_new9;
sol_current10=sol_new10;
sol_current11=sol_new11;
sol_current12=sol_new12;
sol_current13=sol_new13;
sol_current14=sol_new14;
sol_current15=sol_new15;
sol_current16=sol_new16;
else
sol_new1=sol_current1;
sol_new2=sol_current2;
sol_new3=sol_current3;
sol_new4=sol_current4;
sol_new5=sol_current5;
sol_new6=sol_current6;
sol_new7=sol_current7;
sol_new8=sol_current8;
sol_new9=sol_current9;
sol_new10=sol_current10;
sol_new11=sol_current11;
sol_new12=sol_current12;
sol_new13=sol_current13;
sol_new14=sol_current14;
sol_new15=sol_current15;
sol_new16=sol_current16;
end
end
plot(r,E_best,'*')
hold on
end
t=t*a;%（6）降温
end
disp('最优解为：')
disp(sol_best1)
disp(sol_best2)
disp(sol_best3)
disp(sol_best4)
disp(sol_best5)
disp(sol_best6)
disp(sol_best7)
disp(sol_best8)
disp(sol_best9)
disp(sol_best10)
disp(sol_best11)
disp(sol_best12)
disp(sol_best13)
disp(sol_best14)
disp(sol_best15)
disp(sol_best16)
disp('目标表达式的最小值等于：')
disp(E_best)