模拟退火学习

　　模拟退火算法（Simulated Annealing Algorithm)是一种随机类全局优化方法，其最初思想由Metropolis在1953年提出，Kirkpatrick在1983年成功将其应用在组合优化问题中，它来源于热力学中固体物质的退火冷却过程，将固体加热至充分高，在让其徐徐冷却，加热时，固体内部温度升高，粒子变为无序状，内能增大，而徐徐冷却时粒子渐渐有序，在每个温度都达到平衡态，最后在常温时达到基态，温度最低。根据Metropolis准则，粒子在温度T时趋于平衡的概率为，其中，E为温度T时的内能，deltaE为其该变量，k为波尔兹曼(Boltzmann)常数。

　　对于一个目标函数f(X),求其最小值的问题。设f_k、f_k+1分别为目标函数在第k次和第k+1次的迭代值，即f_k=f(X_k),f_k+1=f(X_k+1).若f_k>f_k+1,则接受X_k+1为当前点，作为下一次迭代的初始值继续迭代，直到满足收敛结束条件。若f_k<f_k+1,则X_k+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;

　　模拟退火算法的关键是新的迭代点的随机产生以及对其的接受模式。新的迭代点X_k+1可在X_k附近某个范围内随机产生，范围不能过小，否则会陷入局部极值点。

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('多元函数优化过程')