3 模拟退火算法

3 模拟退火算法

3.1 概述

模拟退火算法（Simulated Annealing Algorithm，SAA）是一种模拟物理退火过程而设计出来的优化算法，退火过程实际上可以理解成一个物体的温度缓慢从高温（能量高，活跃状态）降低到低温（能量低，稳定状态）的过程。

假设当前我们要求解目标函数\(f(x)\)的最小值，如果我们把\(f(x)\)的函数值理解为物体处于状态\(x\)时的能量，那么我们的目标就变成了找到使得粒子能量最低的状态。为了完整模拟物体退火（降温）过程，我们还需要为当前时刻定义一个温度\(T\)。于是，退火算法可以大体上阐述为：

• 首先，设置一个较高的起始温度\(T_0\)并选定一个随机初始状态\(x_0\)（通常\(f(x_0)\)都会比较大）。
• 然后，我们控制温度\(T_0\)缓慢下降至某个终止温度，在温度下降期间通过某种策略产生一系列的状态\(x_1,x_2,\cdots,x_n\)，并且这些状态原则上应该满足\(f(x_1)\ge f(x_2)\ge\cdots\ge f(x_n)\)。
• 最后，我们就得到了一个能量较低的状态\(x_n\)。

当然，这只是一种较为粗略的描述。

3.2 算法具体细节

初始温度的选取：

对于SSA，我们首先需要选定一个初始温度\(T_0\)，初始温度选取的不同将直接影响到最终解的质量。如果初始温度太高，我们会有更高的概率可以获得到更高质量的解，但代价是算法需要花费更长的时间来计算并搜索解。因此，在实际应用的时候，可以先做几组实验看看效果如何，然后在效率和解质量这两者进行权衡。

这里介绍两种经验公式确定初温的方法：首先，随机产生一组状态\(S\)。然后，计算这组状态中两个状态间的最大差值（正数），即

\[\Delta=\max_{s\in S}{f(s)}-\min_{s\in S}f(s) \]

最后，利用经验公式进行确定：

\[\begin{aligned} &(1)\;T_0=-\frac{\Delta}{\ln{p}},\quad p=0.8\\ &(2)\;T_0=K\Delta \end{aligned} \]

其中\(K\)应该是一个比较大的数。

退温策略：

这里介绍两种比较常用的降温策略：

（1）等差方式降温

\[T=T_0-k\Delta T_0\quad (k=0,1,2,\cdots) \]

其中\(k\)表示算法迭代的次数；\(\Delta T_0\)表示每次降温的大小。

（2）等比方式降温

\[T=\lambda^k T_0\quad(0<\lambda<1) \]

理论上，温度降低的速度越慢，获得高质量解的概率越大，但如果温度降得太慢，反而影响算法的效率。一种折中的做法是温度高的时候下降得慢些，温度低时下降得快些。上面介绍的两种降温策略都是按照固定规则进行降温，可以在这基础之上进行改进。

新状态的产生：

如何通过初始状态产生新的状态也是算法的一种重点，应遵循的原则是所产生的候选解（状态）尽可能覆盖了全部的解空间（状态空间）。新状态的产生通常是在当前状态的邻域内，以一定的概率方式抽样出一个新的状态。可供选择的概率分布很多，如均匀分布、高斯分布、指数分布等等，需要自己选择。

需要注意的是，新状态的产生方式需要根据问题的性质来决定，并且如何定义当前状态的邻域也是一个问题。这里列举两个不同的示例：

（1）状态空间是离散的。

如果对于\(1,2,\cdots,n\)这\(n\)个数字的每种排列都分别对应着一个代价，我们的目标是找到使得代价最小的排列。新状态产生函数可以定义为随机交换两个数字的位置，例如

\[123\cdots n\Rightarrow n23\cdots1 \]

类似的可以定义为随机反转某个子序列等。

（2）状态空间是连续的。

如果初始状态是\(x_0\)，那么我们可以从高斯分布中抽样出一个随机向量\(v\)，然后得到新的状态\(x=x_0+v\)。

状态转移策略：

新状态产生之后，如果无条件接受，那么算法就变成了在状态空间随机”游走“，显然是不行的。因此，我们需要指定一个合适的策略使得算法能够向好的方向发展（即函数值减少的方向）。一个常用的准则就是Metropolis准则，公式表示如下：

\[p=\left\{\begin{aligned} &1&\Delta f<0\\ &\exp({-\frac{\Delta f}{T}})&\Delta f\ge0 \end{aligned}\right. \]

其中\(T>0\)表示当前的温度；\(\Delta f=f(x_{new})-f(x_{old})\)；\(p\)表示新状态被接受的概率。

实际上在编程的时候可以将Metropolis准则写成更紧凑的形式：

\[\min\{1,\exp(-\frac{\Delta f}{T})\}\ge random[0,1] \]

如果上述不等式成立，则接受新的状态，实际上这是一种轮盘赌的方式。

也就是如果产生的新状态比当前状态的能量更低（函数值更小）就发生状态转移，进入下一次迭代。如果新状态的能量更高，就以一定的概率接受更差的状态，而不是直接舍弃，这有助于跳出局部最优，发现更好的解。通过公式可知当\(T\)比较大的时候，更有可能接受较差的解。这意味着当温度较高时，算法倾向于在全局范围内随机搜索，搜索能力强。但当温度较低时，算法表现得较为保守，倾向于接受更好的解。

显然这里还有一个问题：如果产生了一个更差的状态，并且不被接受应该怎么办。具体的做法有以下几种：

1. 不断循环抽样，直到产生的新状态被接受为止。
2. 状态不发生改变。
3. 设置一个抽样次数上限，如果达到次数上限后，新状态仍然未被接受，我们可以选择保留最后一次抽样的状态，或是不发生改变。

当然可能还有其他的做法。

停止条件：

1. 设置算法的中止温度。
2. 查看算法在连续多次迭代过程中，如果目标函数值已经不发生变化，或变化很小就停止迭代。

算法描述：

设置好初温\(t=t_0\)，一个随机初始状态\(s=s_0\)，迭代次数计数器\(k=0\)

while 终止条件未满足:
    while True:
        s_new = generate_new_state(s)
        if 满足Metropolis准则:
            s = s_new
            break
	t(k+1) = update(t(k))
    k = k + 1
输出历史最优解

3.3 算法实现与测试

模拟退火算法具有一定的随机性，为了寻找到较好的解，最好限制搜索空间的大小。并且在运行过程中记录历史最优解。

import numpy as np


class SA:
    """
    模拟退火算法：寻找函数的最小值
    """

    def __init__(self, current_temp, gamma, final_temp, x0, st_x, func):
        # 初始温度
        self.current_temp = current_temp
        # 温度采用线性递减策略
        self.gamma = gamma
        # 终止温度
        self.final_temp = final_temp
        # 空间维数
        self.n = x0.size
        # 约束x坐标的范围
        self.st_x = st_x
        # 目标函数
        self.func = func
        # 当前找到的解
        self.solution = x0
        # 当前解对应的函数值
        self.value = func(x0)
        # 历史最优解
        self.best_solution = x0
        self.best_value = func(x0)

    def solve(self):
        while self.current_temp > self.final_temp:
            accept = False
            offset = self.get_offset()
            next_solution = self.solution + offset
            next_solution[next_solution > self.st_x[1]] = self.st_x[1]
            next_solution[next_solution < self.st_x[0]] = self.st_x[0]
            if self.func(next_solution) < self.value:
                accept = True
            elif np.random.rand() < np.exp((self.value - self.func(next_solution)) / self.current_temp):
                accept = True
            if accept:
                self.solution = next_solution
                self.value = self.func(next_solution)
                if self.value < self.best_value:
                    self.best_solution = self.solution
                    self.best_value = self.value
            self.current_temp -= self.gamma
        return self.best_solution, self.best_value

    def get_offset(self):
        # 产生随机偏移量
        v = np.random.rand(self.n) - 0.5
        v *= self.current_temp / 100
        return v


def ackley(x):
    return - 20 * np.exp(-0.2 * np.sqrt((x * x).mean())) - np.exp(np.cos(2 * np.pi * x).mean()) + 20 + np.exp(1)


if __name__ == '__main__':
    # ackley函数的最小值
    for i in range(10):
        sa = SA(
            current_temp=500,
            gamma=0.1,
            final_temp=1,
            st_x=np.array([-5, 5]),
            x0=np.random.randn(2),
            func=ackley
        )
        solution, value = sa.solve()
        print(solution, value)

测试结果：