现代优化算法 之 模拟退火

现代优化算法

现代优化算法是 80 年代初兴起的启发式算法。这些算法包括禁忌搜索（tabu search），模拟退火（simulated annealing），遗传算法（genetic algorithms），人工神经网络（neural networks）。它们主要用于解决大量的实际应用问题。目前，这些算法在理论和实际应用方面得到了较大的发展。无论这些算法是怎样产生的，它们有一个共同的目标-求 NP-hard 组合优化问题的全局最优解。虽然有这些目标，但 NP-hard 理论限制它们只能以启发式的算法去求解实际问题。

启发式算法包含的算法很多，例如解决复杂优化问题的蚁群算法（Ant Colony Algorithms）。有些启发式算法是根据实际问题而产生的，如解空间分解、解空间的限制等；另一类算法是集成算法，这些算法是诸多启发式算法的合成。

现代优化算法解决组合优化问题，如 TSP（Traveling Salesman Problem）问题，QAP（Quadratic Assignment Problem）问题，JSP（Job-shop Scheduling Problem）问题等效果很好。

这一章讲解模拟退火的算法过程，之前也介绍过一些简单的模拟退火的思想，上次是基于ACM-ICPC的思想进行介绍的，这次是详细的计算推导过程。

模拟退火

模拟退火算法得益于材料的统计力学的研究成果。统计力学表明材料中粒子的不同结构对应于粒子的不同能量水平。在高温条件下，粒子的能量较高，可以自由运动和重新排列。在低温条件下，粒子能量较低。如果从高温开始，非常缓慢地降温（这个过程被称为退火），粒子就可以在每个温度下达到热平衡。当系统完全被冷却时，最终形成处于低能状态的晶体。

如果用粒子的能量定义材料的状态，Metropolis 算法用一个简单的数学模型描述了退火过程。假设材料在状态$i$之下的能量为$E(i)$，那么材料在温度$T$时从状态$i$进入状态$j$就遵循如下规律：

• （1）如果$E(j)\leq E(i)$，接受该状态被转换。
• （2）如果$E(j)>E(i)$，则状态转换以如下概率被接受：
  $$e^{E(i)-E(j)\over KT}$$

其中$K$是物理学中的波尔兹曼常数，$T$是材料温度。

在某一个特定温度下，进行了充分的转换之后，材料将达到热平衡。这时材料处于状态$i$的概率满足波尔兹曼分布：

$$P_T(x=i)=\frac{e^{-E(i)\over KT}}{\sum_{j\in S}e^{-E(j)\over KT}}$$
其中 $x$ 表示材料当前状态的随机变量, $S$ 表示状态空间集合。
显然
$$\lim_{T\rightarrow \infty} \frac{e^{-E(i)\over KT}}{\sum_{j\in S}e^{-E(j)\over KT}} = {1\over |S|}$$

其中$|S|$表示集合$S$中状态的数量。这表明所有状态在高温下具有相同的概率。而当温度下降时，

$$\lim_{T\rightarrow 0}{e^{-{E(i)-E_{min}\over KT}}\over \sum_{j\in S}e^{-{E(j)-E_{min}\over KT}}}=\lim_{T\rightarrow 0}{e^{-{E(i)-E_{min}\over KT}}\over \sum_{j\in S_{min}}e^{-{E(j)-E_{min}\over KT}}+\sum_{j\notin S_{min}}e^{-{E(j)-E_{min}\over KT}}}=\lim_{T\rightarrow 0}{e^{-{E(i)-E_{min}\over KT}}\over \sum_{j\in S_{min}}e^{-{E(j)-E_{min}\over KT}}}=\left\{ \begin{aligned} {1\over |S_{min}|} &\ i\in S_{min} \0& \ other \end{aligned} \right.$$

其中$E_{min} =\min_{j\in S}E(j)$且$S_{min} ={i\ |\ E(i)=E_{min}}$。
上式表明当温度降至很低时,材料会以很大概率进入最小能量状态。

假定我们要解决的问题是一个寻找最小值的优化问题。将物理学中模拟退火的思想应用于优化问题就可以得到模拟退火寻优方法。

考虑这样一个组合优化问题:优化函数为 $f : x → R^+$ ,其中 $x ∈ S$ ,它表示优化问题的一个可行解, $R^+ ={y|y∈R,y>0}$，$S$表示函数的定义域。$N(x)?S$表示$x$ 的一个邻域集合。

首先给定一个初始温度$T_0$和该优化问题的一个初始解$x(0)$，并由$x(0)$生成下一个解$x^\prime ∈N(x(0))$，是否接受$x^\prime$作为一个新解$x(1)$依赖于下面概率:

$$P(x(0)\rightarrow x^\prime)=\left\{ \begin{aligned} 1 &\ \ 若f (x^\prime) < f (x(0))\e^{-{f(x^\prime)-f(x(0))}\over T_0} &\ \ other \end{aligned} \right.$$

换句话说，如果生成的解 x’ 的函数值比前一个解的函数值更小,则接受 $x(1) = x$’ 作为一个新解。
否则以概率 $e^{? {f (x‘)? f (x(0))}\over T_0}$ T0 接受 $x$’ 作为一个新解。

泛泛地说,对于某一个温度$T_i$和该优化问题的一个解$x(k)$, 可以生成$x$’。接受$x$’ 作为下一个新解 $x(k + 1)$ 的概率为:

$$P(x(k)\rightarrow x^\prime)=\left\{ \begin{aligned} 1 &\ \ 若f (x^\prime) < f (x(k))\e^{-{f(x^\prime)-f(x(k))}\over T_0} &\ \ other \end{aligned} \right.\tag{1}$$

在温度$T_i$下,经过很多次的转移之后,降低温度$T_i$ ,得到$T_{i+1} < T_i$ 。在Ti+1 下重复上述过程。因此整个优化过程就是不断寻找新解和缓慢降温的交替过程。最终的解是对该问 题寻优的结果。

我们注意到,在每个$T_i$ 下,所得到的一个新状态$x(k +1)$完全依赖于前一个状态 $x(k)$ , 可以和前面的状态 $x(0),\dots , x(k ? 1)$ 无关,因此这是一个马尔可夫过程。使用马 尔可夫过程对上述模拟退火的步骤进行分析,结果表明:从任何一个状态 x(k ) 生成 x’ 的 概率,在$N(x(k))$中是均匀分布的,且新状态x’被接受的概率满足式(1),那么经过有限次的转换,在温度$T_i$下的平衡态 $x_i$ 的分布由下式给出:

$$P_i(x_i)={e^{-f(x-i)\over T}\over \sum_{j\in S}e^{-{f(x_i)\over T_i}}}$$

当温度$T$降为$0$时, $x_i$ 的分布为:

$$P_i^*=\left\{ \begin{aligned} {1\over |S_{min}|} &\ i\in S_{min} \0& \ other \end{aligned} \right.$$

并且

$$\sum_{x_i\in S_{min}}P_i =1$$

这说明如果温度下降十分缓慢,而在每个温度都有足够多次的状态转移,使之在每一个 温度下达到热平衡,则全局最优解将以概率 1 被找到。因此可以说模拟退火算法可以找 到全局最优解。

在模拟退火算法中应注意以下问题:

(1)理论上,降温过程要足够缓慢,要使得在每一温度下达到热平衡。但在计算 机实现中,如果降温速度过缓,所得到的解的性能会较为令人满意,但是算法会太慢, 相对于简单的搜索算法不具有明显优势。如果降温速度过快,很可能最终得不到全局最 优解。因此使用时要综合考虑解的性能和算法速度,在两者之间采取一种折衷。

(2)要确定在每一温度下状态转换的结束准则。实际操作可以考虑当连续 m 次的 转换过程没有使状态发生变化时结束该温度下的状态转换。最终温度的确定可以提前定 为一个较小的值Te ,或连续几个温度下转换过程没有使状态发生变化算法就结束。

(3)选择初始温度和确定某个可行解的邻域的方法也要恰当。