《Nature-Inspired Metaheuristic Algorithms》——模拟退火算法 SIMULATED ANNEALING

模拟退火算法

1.1什么是退火？

退火从物理现象上来说就是物体逐渐降温的物理现象，物体的温度越低能量也就越低，当能量足够低的时候就会开始结晶，逐渐达到物体的最低能量状态。但是如果物体的温度下降的过快，物体不会结晶，而是处于一种次最低能量的状态，物体中仍然保持一定的能量没有进行释放。

如图进行退火：
①对物体进行加热（此时物体能量较低，但不是结晶态，因为能量还未达到最低）
②加热后的物体能量较高分子间运动活跃，然后对其进行足够缓慢地降温
③当物体温度缓慢下降到一定程度的时候，物体开始结晶而达到了最低能量状态。
退火完成。

1.2 模拟退火中得到的优化算法启发

现象:缓缓降温，使得物体分子在每一温度时，能够有足够时间找到安顿位置，则逐渐地，到最后可得到最低能态，系统最安稳。

由此在最基础的爬山算法中引入了温度T这个新参数，当温度高的时候，如果新的解更优则立即选择这个新的解，而即使新的解不是更优解仍然以较高的概率p（相对于T较低时的概率是较高的选择这个新的解（T->∞时，p->1）；当温度低的时候，如果这个解是更优解则依然立即选择，而这个解不是更优解则会以很小的概率（T->0时，p->0）去选择这个新解。

其中，p的计算公式为：

一般的，k取为1， ∆E=Snew - Sold，也即新的函数适应度值-之前的函数适应度值，T为温度。

个人感觉：这个退火现象和算法的关系最大的就是这个温度的引入以及温度的缓慢下降才能获取到最优解。

2.1 模拟退火算法的核心思想

模拟退火算法的核心思想是：从一个初始解开始，随机地在这个初始解附近进行跳跃，如果跳跃到的新点position对应的value更佳，那么就得到这个新的最优解；如果在这个新的解不比之前的最优解更优的话，就以p的概率接受这个解为。

e的负指数次幂小于1，因此可以十分方便地直接将p和范围为(0,> 1)的随机数做比较，若p大于随机数则接受该解。随着温度的降低，p的值也逐渐降低，于是接受更差的解的概率也就越小，退火过程就趋于稳定。

2.2 以小兔子爬山举例SA

普通的爬山算法：一只小兔子想要爬到山峰的顶端，但是他不知道哪里才是最高的山峰，每次他向左或者向右进行小距离的跳动，如果跳动之后的位置更高那么就来到这个位置，如果更低则回到原来的位置，并且再次开始向左向右的跳动。当小兔子爬到第一个山峰的时候就会陷入到一个局部最优，无论向左还是向右都只会回到原来这个位置。

加入随机干扰因素之后的爬山算法（也即SA）：同样是更优则接受，不优则保持原位，加入一个随机干扰也就是，即使这个新的位置不优，也有一定的概率接受这个新的解，那么小兔子就有可能突破这个第一个山峰一直走到最高那个山峰。

一直没有提到的温度T：温度T是的主要作用有两点，①是通过温度来判断是否结束迭代过程，每次迭代完成之后T=T*α（通常α=0.7~0.99）②通过温度计算即使非更优依然接受解的概率p=e^(-∆E/kT)随着温度T的降低，p的值也逐渐降低，于是接受更差的解的概率也就越小，退火过程就趋于稳定。

2.3 使用SA求解f(x)=x^2在[-2,2]之间的最小值

这个区间里，随机选择一个点， Xold=4*（rand-0.5) ，以及它对应的函数的数值S=f(Xold) ，然后以这个点作为起点。接下来产生一个随机扰动，扰动的大小可以自己设计一些参数来调节Xnew=Xold+k*(rand-0.5)，并计算出该点对应的函数值Snew=f(Xnew)
(rand 为0~1之间的随机数)

对两个数值求差值，就可以判断新的解是否更接近最优解。∆E=Snew - Sold，如果dE小于0，则说明新的接更加接近我们想要的解，即接受该解。

如果dE大于0，则说明新的解离我们想要的解更远，由于加入了随机干扰因素，仍然保持p的概率接受该解。

3.1 模拟退火算法代码

模拟退火的精髓，判断函数

function [y] = judge(dE, t)
  if(dE < 0)
    y = 1;
  else
    d = exp(-(dE / t));
    if(d > rand)
      y = 1;
    else
      y = 0;
    end
  end
end

函数传入两个参数，dE=Snew - Sold与温度t，然后返回是否接受该解。

例：使用模拟退火算法求解[0， 9]范围内函数f(x)=x+10sin5x+7cos4x的最大值。

clear;

%范围
section_l = 0;
section_h = 9;

%绘制函数图像
draw_base();

%初始温度，停止温度与降温系数
tmp = 1e5;
tmp_min = 1e-3;
alpha = 0.98;

%生成初始随机解
x_old = (section_h - section_l) * rand() + section_l;
x_new = x_old;
s_old = val(x_old);
s_new = s_old;

text_lt = text(0, 0, 'Init');

%计数器
counter = 0;

%退火的主体部分，一个循环
while(tmp > tmp_min)
  %随机扰动
  delta = (rand() - 0.5) * 3;
  x_new = x_old + delta;
  %扰动的值小于一半的区间范围时，可以用这种办法防止新值超出区间范围
  if(x_new < section_l || x_new > section_h)
    x_new = x_new - 2 * delta;
  end
 
  s_new = val(x_new);

  %求差值，这里是找最大值而非最小值，所以不是s_new - s_old
  dE = s_old - s_new;
        
  %判断
  j = judge(dE, tmp);
  if(j)
    s_old = s_new;
    x_old = x_new;
  end
          
  %只有当dE < 0的情况下才降温
  if(dE < 0)
    delete(text_lt);
    hold on, scatter(x_old, s_old);
    hold on, text_lt = text(0.3, 21, ['tmp: ', num2str(tmp)]);
    pause(0.01);
    %上面是绘图的代码，下面降温
    tmp = tmp * alpha;
  else
    counter = counter + 1;
  end
            
  %当接受更差的解的概率太小时，若又长时间找不到更优的解，那么退出循环，结束算法
  if(counter > 10000)
    break;
  end
    
end

function [y] = val(x)
    y = x + 10 * sin(5 * x) + 7 * cos(4 * x);
end

function draw_base()
    X = 0: 0.01:9;
    M = val(X);
    plot(X, M);
end

Pseudocode of Simulated Annealing Algorithm

SA的伪代码
[1]模拟退火算法学习笔记 知乎 李美哲https://zhuanlan.zhihu.com/p/33184423