模拟退火算法(Simulated Annealing Algorithm)是一种随机类全局优化方法,其最初思想由Metropolis在1953年提出,Kirkpatrick在1983年成功将其应用在组合优化问题中,它来源于热力学中固体物质的退火冷却过程,将固体加热至充分高,在让其徐徐冷却,加热时,固体内部温度升高,粒子变为无序状,内能增大,而徐徐冷却时粒子渐渐有序,在每个温度都达到平衡态,最后在常温时达到基态,温度最低。根据Metropolis准则,粒子在温度T时趋于平衡的概率为,其中,E为温度T时的内能,deltaE为其该变量,k为波尔兹曼(Boltzmann)常数。
对于一个目标函数f(X),求其最小值的问题。设fk、fk+1分别为目标函数在第k次和第k+1次的迭代值,即fk=f(Xk),fk+1=f(Xk+1).若fk>fk+1,则接受Xk+1为当前点,作为下一次迭代的初始值继续迭代,直到满足收敛结束条件。若fk<fk+1,则Xk+1可能被接受也可能不被接受,接受的概率为Boltzmann概率p:
其中,T为控制参数,在迭代过程中,T必须缓慢减少,T变化太快,会使优化陷入局部极值点。
一元函数优化举例:同样计算函数f(x)=x*cos(5*pi*x)+3.5在区间[-1,2.5]上的最大值
先在工作区建立上面函数的一个m文件并保存,命名为fun_sigv.m:
function y=fun_sigv(x) y=x.*cos(5*pi*x)+3.5;
模拟退火算法的关键是新的迭代点的随机产生以及对其的接受模式。新的迭代点Xk+1可在Xk附近某个范围内随机产生,范围不能过小,否则会陷入局部极值点。
opt_minmax=1; %目标优化类型:1最大化、-1最小化 sub=-1; %变量取值下限 up=2.5; %变量取值上限 delt=(up-sub)/5; yita=0.99; trace=[]; %模拟退火迭代性能跟踪器 k_total=3000; %迭代总次数 tx=sub:.01:up; y=fun_sigv(tx); T=max(y)-min(y); %模拟温度初始化 plot(tx,y) xlabel('x') ylabel('y') title('一元函数优化结果') hold on x0=sub+(up-sub)*rand; f0=fun_sigv(x0); %随机产生初始点 k=0; plot(x0,f0,'ro','linewidth',2) %在函数图像上标出初始点位置 while k<k_total x1=x0+(2*rand-1)*delt; x1=min(x1,up); x1=max(x1,sub); %在当前点附近随机产生下一个迭代点位置,并保证在所考查区域内 f1=fun_sigv(x1); if opt_minmax*f1>opt_minmax*f0 %迭代点优于当前点,接受迭代结果并设置为当前点 x0=x1; f0=f1; elseif rand<exp(opt_minmax*(f1-f0)/T) %迭代点劣于当前点,以boltzmann概率接受迭代结果并设置为当前点 x0=x1; f0=f1; end T=yita*T; %以缓慢速度衰减温度T k=k+1; trace=[trace;f0]; %跟踪模拟退火的迭代优化过程 [x0,f0] end plot(x0,f0,'r*','linewidth',2) figure plot(trace(:),'r-*') hold on xlabel('迭代次数') ylabel('目标函数优化情况') title('一元函数优化过程')
多元函数优化举例:
多峰的Shubert为:
求f(x,y)在[-10,10]x[-10,10]上的最大值。
同样新建 fun_mutv函数为:
function my=fun_mutv(x,y) t1=zeros(size(x)); t2=t1; for i=1:5 t1=t1+i*cos((i+1)*x+i); t2=t2+i*cos((i+1)*y+i); end my=t1.*t2;
clc clear all opt_minmax=1; %目标优化类型:1最大化、-1最小化 sub=-10; %变量取值下限 up=10; %变量取值上限 delt=(up-sub)/5; yita=0.99; trace=[]; %模拟退火迭代性能跟踪器 k_total=3000; %迭代总次数 [tx,ty]=meshgrid(sub:.1:up); tz=fun_mutv(tx,ty); T=max(max(tz))-min(min(tz)); %模拟温度初始化 mesh(tx,ty,tz) xlabel('x') ylabel('y') zlabel('z') title('多元函数优化结果') hold on x0=[sub+(up-sub)*rand,sub+(up-sub)*rand]; f0=fun_mutv(x0(1),x0(2)); %随机产生初始点 k=0; plot3(x0(1),x0(2),f0,'ko','linewidth',2) %在函数图像上标出初始点位置 while k<k_total x1=[x0(1)+(2*rand-1)*delt,x0(2)+(2*rand-1)*delt]; x1=[min(x1(1),up),min(x1(2),up)]; x1=[max(x1(1),sub),max(x1(2),sub)]; %在当前点附近随机产生下一个迭代点位置,并保证在所考查区域内 f1=fun_mutv(x1(1),x1(2)); if opt_minmax*f1>opt_minmax*f0 %迭代点优于当前点,接受迭代结果并设置为当前点 x0=x1; f0=f1; elseif rand<exp(opt_minmax*(f1-f0)/T) %迭代点劣于当前点,以boltzmann概率接受迭代结果并设置为当前点 x0=x1; f0=f1; end T=yita*T; %以缓慢速度衰减温度T k=k+1; trace=[trace;f0]; %跟踪模拟退火的迭代优化过程 [x0,f0] end plot3(x0(1),x0(2),f0,'k*','linewidth',2) figure plot(trace(:),'r-*') hold on xlabel('迭代次数') ylabel('目标函数优化情况') title('多元函数优化过程')