模拟退火算法（SA）详解

article/2025/9/20 23:49:05

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

二.相关知识背景

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

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

$W=\frac{N!}{n_{0}!n_{1}!}(1)$

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

$S=KlnW\left ( 2 \right )$

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

$S=Kln\frac{N!}{n_{0}!n_{1}!}\left ( 3 \right )$

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

$S'=Kln\frac{N!}{\left (n _{0} -1\right )!\left (n _{1} +1\right )!}\left ( 4 \right )$

由(4)-(3)式得

$\Delta S=Kln\frac{n_{0}}{n_{1}+1}\rightarrow Kln\frac{n_{0}}{n_{1}}\left ( 5 \right )$

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

$\Delta S=\frac{\Delta Q}{T}=\frac{E_{1}-E_{0}}{T}\left ( 6 \right )$

(5)、(6)式联立得

$Boltzmann factor\Rightarrow \frac{n_{0}}{n_{1}}=e^{-\frac{E_{0}-E_{1}}{KT}}$

三.模拟退火过程

想象一下如果我们现在有下面这样一个函数，现在想求函数最小值的（全局）最优解。如果采用贪心策略，那么从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;ends_new = val(x_new);%求差值，这里是找最大值而非最小值，所以不是s_new - s_olddE = 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;elsecounter = counter + 1;end%当接受更差的解的概率太小时，若又长时间找不到更优的解，那么退出循环，结束算法if(counter > 10000)break;endendfunction [y] = val(x)y = x + 10 * sin(5 * x) + 7 * cos(4 * x);
endfunction draw_base()X = 0: 0.01:9;M = val(X);plot(X, M);
endfunction [y] = judge(dE, t)if(dE < 0)y = 1;elsed = exp(-(dE / t));if(d > rand)y = 1;elsey = 0;endend
end