在记录模拟退火之前,我们来看看它的启发灵感–(详情见WIKI)

“模拟退火”来自冶金学术语退火,是将材料加热后再经特定速率冷却的技术,目的是增大晶粒的体积,并且减少晶格中的缺陷,以改变材料的物理性质。材料中的原子原来会停留在使内能有局部最小值的位置,加热使能量变大,原子会离开原来位置,而随机在其他位置中移动。退火冷却时速度较慢,使得原子有较多可能可以找到内能比原先更低的位置。

模拟退火的原理也和金属退火的原理近似:我们将热力学的理论套用到统计学上,将搜寻空间内每一点想像成空气内的分子;分子的能量,就是它本身的动能;而搜寻空间内的每一点,也像空气分子一样带有“能量”,以表示该点对命题的合适程度。算法先以搜寻空间内一个任意点作起始:每一步先选择一个“邻居”,然后再计算从现有位置到达“邻居”的概率。

可以证明,模拟退火算法所得解依概率收敛到全局最优解。

模拟退火法可用于精确算法失效的高难度计算优化问题,虽然通常只能获得全局最优的近似,但对很多实际问题已经足够。

模拟退火直观收敛图:
模拟退火


模拟退火(SA)是一种优化算法,用来找全局最优。

如上面的灵感所说,模拟退火算法模仿的是金属退火过程,分三个阶段:加温、等温、冷却。

从物理学的角度来说,

  • 加温是将金属加热到高温,原子活动剧烈
  • 等温是保持温度,让原子有足够时间寻找低能位置
  • 冷却是缓慢降温,使原子逐渐稳定在低能状态

而在算法中,温度T的意义是接受差解的概率,能量对应目标函数值,状态代表当前的解。

算法详解

超参数

先设置参数:

  • 初始解:x0\mathbf{x}_0 (这个初始解是随机生成的,所以你运气够好的话可以一发入魂,直接碰到全局最优(bushi))
  • 初始温度:T0=ΔEavglnP0T_0 = \dfrac{-\Delta E_{\text{avg}}}{\ln P_0}
  • 降温系数:α\alpha(通常 0.8α0.990.8 \leq \alpha \leq 0.99
  • 马尔可夫链长度:LL
  • 终止温度:TminT_{\text{min}}

来讨论参数对算法的影响:

温度TT
TT 很大时,exp(ΔE/T)1\exp(-\Delta E/T) \approx 1,对全局的探索能力就很强,有利于跳出局部最优。
TT 很小时,exp(ΔE/T)0\exp(-\Delta E/T) \approx 0,差不多收敛了,找局部最优。

降温系数α\alpha

α\alpha较大,接近 1:说明降温慢,利于搜索,不过相应的计算时间就要长一些。
α\alpha较小:降温快,虽然计算快,但是容易错过全局最优,陷入局部最优。

马尔可夫链长度LL

LL太大:每个温度停留太久,效率低
LL太小:未达到热平衡就降温,效果差

超参的设置依据实际背景来,如果你不知道该怎么设置的话,记得老祖宗传下来的道理——中庸。


以下是对算法过程的口头叙述:

我们的目标函数

E(x)E(x)

xx是当前解向量。

初始温度(初始温度似乎有比较通用的设置,对数分母设置为ln(0.8)\ln (0.8))
这个公式比较常用:

T0=ΔEmaxlnP0T_0 = \frac{-\Delta E_{\text{max}}}{\ln P_0}

其中,ΔEmax\Delta E_{\text{max}} 是初始随机解中最大能量差估计值,P0P_0 是初始接受概率,通常把P0P_0设置为0.8。

先随机生成一个初始解,就像随机把你放在地形中的某个位置。然后基于当前解,通过某种策略生成一个"邻居"解。就像在你当前位置周围随机迈出一步

接着决定是否接受新解,这是算法的精髓所在。使用Metropolis准则来决定是否接受新解:

Metropolis 接受准则

从当前解 x\mathbf{x} 转移到新解 x\mathbf{x'} 的接受概率:

P={1ΔE<0eΔETΔE0P = \begin{cases} 1 & \text{, } \Delta E < 0 \\ e^{-\frac{\Delta E}{T}} & \text{, } \Delta E \geq 0 \end{cases}

其中:

ΔE=EnewEcurrent\Delta E = E_{\text{new}} - E_{\text{current}} 是能量变化;TT 是当前温度

如果新解更好(能量更低),总是接受;如果新解更差,以一定概率接受:P = exp(-ΔE/T)

然后是降温,常见的降温方式:

指数降温:

Tk+1=αTkT_{k+1} = \alpha \cdot T_k

对数降温:

Tk=T0ln(1+k)T_k = \dfrac{T_0}{\ln(1 + k)}

线性降温:

Tk=T0kδT_k = T_0 - k \cdot \delta

其中 δ\delta 是每次迭代的温度下降量.

降温速度挺重要的,太快容易陷入局部最优,太慢计算时间长。

终止条件有这些:

温度达到最小:

TkTminT_k \leq T_{\text{min}}

能量变化过小:

EbestEcurrentε|E_{\text{best}} - E_{\text{current}}| \leq \varepsilon

还有就是设定迭代次数,迭代完了自然就完了。

实际问题

SA算法解决问题时会根据问题来设置领域生成函数,对于连续优化、组合优化(比如TSP问题)等不同问题有不同的领域生成函数。具体的函数就根据问题背景来设置。

算法伪代码

初始化 x=x0,T=T0,k=0while T>Tmin dofor i=1 to L dox=generate_neighbor(x)ΔE=E(x)E(x)if ΔE0 or random()<exp(ΔE/T)x=xend ifend forT=αTk=k+1end while返回 xbest\begin{aligned} &\text{初始化 } \mathbf{x} = \mathbf{x}_0, T = T_0, k = 0 \\ &\text{while } T > T_{\text{min}} \text{ do} \\ &\quad \text{for } i = 1 \text{ to } L \text{ do} \\ &\quad \quad \mathbf{x}' = \text{generate\_neighbor}(\mathbf{x}) \\ &\quad \quad \Delta E = E(\mathbf{x}') - E(\mathbf{x}) \\ &\quad \quad \text{if } \Delta E \leq 0 \text{ or } \text{random}() < \exp(-\Delta E/T) \\ &\quad \quad \quad \mathbf{x} = \mathbf{x}' \\ &\quad \quad \text{end if} \\ &\quad \text{end for} \\ &\quad T = \alpha \cdot T \\ &\quad k = k + 1 \\ &\text{end while} \\ &\text{返回 } \mathbf{x}_{\text{best}} \end{aligned}