转载自:https://blog.csdn.net/pangel18/article/details/52349840
模拟退火用于处理最优化问题,可以求出当目标函数取得最小值时的决策变量的值。
在编写程序时需要根据具体问题设计算法,算法描述为:
(1)解空间(初始解)
(2)目标函数
(3)新解的产生
① 2 变换法
② 3 变换法
(4)代价函数差
(5)接受准则
(6)降温
(7)结束条件
下面MATLAB程序用于求解非线性规划:
min f(x)=x1^2+x2^2+8
st.
x1^2-x2>=0
-x1-x2^2+2=0
x1,x2>=0
- clear
- clc
- %生成初始解
- sol_new2=1;%(1)解空间(初始解)
- sol_new1=2-sol_new2^2;
- sol_current1 = sol_new1;
- sol_best1 = sol_new1;
- sol_current2 = sol_new2;
- sol_best2 = sol_new2;
- 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_new2=sol_new2+rand*0.2;
- sol_new1=2-sol_new2^2;
- %检查是否满足约束
- if sol_new1^2-sol_new2>=0 && -sol_new1-sol_new2^2+2==0 && sol_new1>=0 &&sol_new2>=0
- else
- sol_new2=rand*2;
- sol_new1=2-sol_new2^2;
- continue;
- end
- %退火过程
- E_new=sol_new1^2+sol_new2^2+8;%(2)目标函数
- if E_new<E_current%(5)接受准则
- E_current=E_new;
- sol_current1=sol_new1;
- sol_current2=sol_new2;
- if E_new<E_best
- %把冷却过程中最好的解保存下来
- E_best=E_new;
- sol_best1=sol_new1;
- sol_best2=sol_new2;
- end
- else
- if rand<exp(-(E_new-E_current)/t)%(4)代价函数差
- E_current=E_new;
- sol_current1=sol_new1;
- sol_current2=sol_new2;
- else
- sol_new1=sol_current1;
- sol_new2=sol_current2;
- end
- end
- plot(r,E_best,'*')
- hold on
- end
- t=t*a;%(6)降温
- end
- disp('最优解为:')
- disp(sol_best1)
- disp(sol_best2)
- disp('目标表达式的最小值等于:')
- disp(E_best)
- 最优解为:
- 1.0038
- 0.9981
- 目标表达式的最小值等于:
- 10.0038
司守奎的算法大全给出了利用模拟退火求解TSP问题的MATLAB算法:
- clc,clear
- load sj.txt %加载敌方100 个目标的数据,数据按照表格中的位置保存在纯文本
- 文件sj.txt 中
- x=sj(:,1:2:8);x=x(:);
- y=sj(:,2:2:8);y=y(:);
- sj=[x y];
- d1=[70,40];
- sj=[d1;sj;d1];
- sj=sj*pi/180;
- %距离矩阵d
- d=zeros(102);
- for i=1:101
- for j=i+1:102
- temp=cos(sj(i,1)-sj(j,1))*cos(sj(i,2))*cos(sj(j,2))+sin(sj(i,2))*sin(sj(j,2));
- d(i,j)=6370*acos(temp);
- end
- end
- d=d+d';
- S0=[];Sum=inf;
- rand('state',sum(clock));
- for j=1:1000
- S=[1 1+randperm(100),102];
- temp=0;
- -276-
- for i=1:101
- temp=temp+d(S(i),S(i+1));
- end
- if temp<Sum
- S0=S;Sum=temp;
- end
- end
- e=0.1^30;L=20000;at=0.999;T=1;
- %退火过程
- for k=1:L
- %产生新解
- c=2+floor(100*rand(1,2));
- c=sort(c);
- c1=c(1);c2=c(2);
- %计算代价函数值
- df=d(S0(c1-1),S0(c2))+d(S0(c1),S0(c2+1))-d(S0(c1-1),S0(c1))-d(S0(c2),S0(c2+1));
- %接受准则
- if df<0
- S0=[S0(1:c1-1),S0(c2:-1:c1),S0(c2+1:102)];
- Sum=Sum+df;
- elseif exp(-df/T)>rand(1)
- S0=[S0(1:c1-1),S0(c2:-1:c1),S0(c2+1:102)];
- Sum=Sum+df;
- end
- T=T*at;
- if T<e
- break;
- end
- end
- % 输出巡航路径及路径长度
- S0,Sum