模拟退火算法（SA）详解

一.算法来源：模拟退火算法(SA)源于固体退火原理，是一种基于概率的算法。将固体加温至充分高的温度，再让其徐徐冷却，加温时，固体内部粒子随温升变为无序状，内能增大，分子和原子越不稳定。而徐徐冷却时粒子渐趋有序，能量减少，原子越稳定。在冷却（降温）过程中，固体在每个温度都达到平衡态，最后在常温时达到基态，内能减为最小。

---

二.相关知识背景

玻尔兹曼分布推导（可跳过）

假设我们有一个包含宏观粒子总数为N、环境温度为T、定容条且达到热力平衡的封闭系统。假定此处只有、粒子是有效的，它们分别包含了、个微粒。 其中N=+。总微观状态数为

微观状态数(W)和熵(S)的关系在统计物理内的定义：

其中K为普朗克常数，依据（2）式则当前系统熵为：

若给予一定能量使1个处于状态的粒子激发到状态，则此时系统熵为：

 由(4)-(3)式得

 

熵(S)与系统从外界吸收的热量(Q)在热力学中的关系为:

(5)、(6)式联立得

​​​​​​​

---

 三.模拟退火过程

想象一下如果我们现在有下面这样一个函数，现在想求函数最小值的（全局）最优解。如果采用贪心策略，那么从A点开始试探，如果函数值继续减少，那么试探过程就会继续。而当到达点B时，显然我们的探求过程就结束了。最终我们只能找打一个局部最后解B，显然不是我们想要的全局最优解C。相交于普通的贪心算法，模拟退火算法以一定的概率来接受一个比当前解要差的解，因此有可能会跳出这个局部的最优解，达到全局的最优解。

将函数值视为能量，模拟退火算法使第n次迭代状态与n+1次迭代状态进行比较。定义系统由x(n)变为x(n+1)的接受概率P为：

当E(n+1)>E(n)时， Metropolis接受准则使模拟退火算法具有一定爬坡能力，且这种爬坡能力会随着温度降低逐渐降低，最终趋于稳定。

用固体退火模拟组合优化问题，将内能E模拟为目标函数值f，温度T演化成控制参数t，即得到解组合优化问题的模拟退火算法：由初始解i和控制参数初值t开始，对当前解重复“产生新解→计算目标函数差→接受或舍弃”的迭代，并逐步衰减t值，算法终止时的当前解即为所得近似最优解，退火过程由冷却进度表(Cooling Schedule)控制，包括控制参数的初值t及其衰减因子Δt、每个t值时的迭代次数L和停止条件Tf。而温度的作用就是来计算转移概率P的。当温度每次下降后，转移概率也发生变化，因此在所有温度下迭代L次的结果也都是不相同的。在每个温度下迭代L次来寻找当前温度下的最优解，然后降低温度继续寻找，直到到达终止温度，即转移概率P接近于0.

四.应用示例（matlab）

求解：x+10sin(5x)+7cos(4x)在[0,9]区间最大值

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

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