经典算法-模拟退火算法的python实现

模拟退火算法基本思想

模拟退火算法来源于固体退火原理，将固体加温至充分高，再让其徐徐冷却。加温时，固体内部粒子随温度升高变为无序状，内能增大，而缓慢冷却时粒子又逐渐趋有序。从理论上讲，如果冷却过程足够缓慢，那么冷却中任一温度时固体都能达到热平衡，而冷却到低温时将达到这一低温下的内能最小状态。

LLM大模型相关文章：

大模型查询工具助手之股票免费查询接口

GPT实战系列-LangChain + ChatGLM3构建天气查询助手

GPT实战系列-大模型为我所用之借用ChatGLM3构建查询助手

GPT实战系列-简单聊聊LangChain

GPT实战系列-Baichuan2本地化部署实战方案

GPT实战系列-大话LLM大模型训练

GPT实战系列-ChatGLM3本地部署CUDA11+1080Ti+显卡24G实战方案

根据Metropolis准则，粒子在温度T时趋于平衡的概率为e(-ΔE/(kT))，其中E为温度T时的内能，ΔE为其改变量，k为Boltzmann常数。用固体退火模拟组合优化问题，将内能E模拟为目标函数值f，温度T演化成控制参数t，即得到解组合优化问题的模拟退火算法：由初始解i和控制参数初值t开始，对当前解重复“产生新解→计算目标函数差→接受或舍弃”的迭代，并逐步衰减t值，算法终止时的当前解即为所得近似最优解，这是基于蒙特卡罗迭代求解法的一种启发式随机搜索过程。退火过程由冷却进度表(Cooling Schedule)控制，包括控制参数的初值t及其衰减因子Δt、每个t值时的迭代次数L和停止条件S。

因为物理系统总是趋向于能量最低，而分子热运动则趋向于破坏这种低能量的状态，故而只需着重取贡献比较大的状态即可达到比较好的效果，因而1953年Metropolis提出了这样一个重要性采样的方法，即设从当前状态i生成新状态j。若新状态的内能小于状态i的内能(即Ej<Ei)，则接受新状态j作为新的当前状态；否则，以概率 $exp[\frac{-(E_j-E_i)}{KT}]$ 接受状态j，其中K为Boltzmann常数，这就是通常所说的Metropolis准则。

初始化算法参数

    self.interval = interval                                    # 给定状态空间 - 即待求解空间self.T_max = T_max                                          # 初始退火温度 - 温度上限self.T_min = T_min                                          # 截止退火温度 - 温度下限self.iterMax = iterMax                                      # 定温内部迭代次数self.rate = rate                                            # 退火降温速度

算法基本步骤

①令 $T=T_0$ ，即开始退火的初始温度，随机生成一个初始解工，并计算相应的目标函数值 $E(x_0)$ 。

②令Ｔ等于冷却进度表中的下一个值 $T_i$ 。

③根据当前 $x_i$ ，进行扰动(扰动方式可以参考后面的实例)，产生一个新解 $x_j$ 、计算应的目标函数值 $E(x_j)$ ，得到$ △E=E(x_j)一E(x_i)$。

④若△E<0，则新解 $x_j$ 被接受，作为新的当前解；若 $△ E > 0$ ，则新解 $x_j$ ，按概率 $exp(-△E/T_i)$ 接受， $T_i$ 为当前温度。

⑤在温度 $T_i$ 下，重复 $L_k$ 次的扰动和接受过程，即执行步骤③与④。

⑥判断T是否已到达 $T_j$ ，是，则终止算法；否，则转到步骤②继续执行。

		x1 = self.x_seedT = self.T_maxwhile T >= self.T_min:for i in range(self.iterMax):f1 = self.func(x1)delta_x = random.random() * 2 - 1if x1 + delta_x >= self.interval[0] and x1 + delta_x <= self.interval[1]: # 将随机解束缚在给定状态空间内x2 = x1 + delta_xelse:x2 = x1 - delta_xf2 = self.func(x2)delta_f = f2 - f1x1 = deal(x1, x2, delta_f, T)T *= self.rateself.x_solu = x1                                            # 提取最终退火解

算法实质分两层循环，在任一温度随机扰动产生新解，并计算目标函数值的变化，决定是否被接受。由于算法初始温度比较高，这样，使E增大的新解在初始时也可能被接受。因而能跳出局部极小值，然后通过缓慢地降低温度，算法就最终可能收敛到全局最优解。还有一点要说明的是，虽然在低温时接受函数已经非常小了，但仍不排除有接受更差的解的可能，因此一般都会把退火过程中碰到的最好的可行解(历史最优解)也记录下来，与终止算法前最后一个被接受解一并输出。

参数说明

退火过程由一组初始参数，即冷却进度表(cooling schedule) 控制，它的核心是尽量使系统达到准平衡，以使算法在有限的时间内逼近最优解。冷却进度表包括：

①控制参数的初值 $T_0$ ：冷却开始的温度。
②控制参数Ｔ的衰减函数：因计算机能够处理的都是离散数据，因此需要把连续的降温过程离散化成降温过程中的一系列温度点，衰减函数即计算这一系列温度的表达式。
③控制参数T的终值 $T_j$ ，(停止准则)。
④Markov链的长度 $L_k$ ：任一温度T的迭代次数。

参数的选择

(1)控制参数Ｔ的初值 $T_0$

求解全局优化问题的随机搜索算法一般都采用大范围的粗略搜索，与局部的精细搜索相结合的搜索策略。只有在初始的大范围搜索阶段找到全局最优解所在的区域，才能逐渐缩小搜索的范围.最终求出全局最优解。模拟退火算法是通过控制参数Ｔ的初值 $T_0$ 和其衰减变化过程，来实现大范围的粗略搜索和局部精细搜索。

一般来说，只有足够大的 $T_0$ 才能满足算法要求(但对不同的问题“足够大”的含义也不同，有的可能 $T_0$ =100就可以，有的则要1000)。在问题规模较大时，过小的 $T_0$ 往往导致算法难以跳出局部陷阱而达不到全局最优。但为了减少计算量， $T_0$ 不宜取得过大，而应与其他参数折中选取。

(2)控制参数Ｔ的衰减函数

衰减函数可以有多种形式，一个常用的衰减函数是
$T_{k+1} = \alpha T_k, k=0, 1, 2, ...$

其中.a是一个常数，可以取为0.5~0.99，它的取值决定了降温的过程。小的衰减量可能导致算法进程迭代次数的增加，从而使算法进程接受更多的变换，访问更多的邻域，搜索更大范围的解空间，返回更好的最终解。同时由于在 $T_k$ 值上已经达到准平衡，则在 $T_{k+1}$ 时只需少量的变换就可达到准平衡。这样就可选取较短长度的Markov链来减少算法时间。

(3) Markov链长度

选取原则：

在控制参数Ｔ的衰减函数已选定的前提下， $L_k$ 应能使在控制参数Ｔ的每一取值上达到准平衡。从经验上来说，对简单的情况可以令 $L_k$ =100n，n为问题规模。

算法停止准则：

对Metropolis准则中的接受函数 $exp[\frac{-(E_j-E_i)}{KT}]$ 分析可知，在Ｔ比较大的高温情况下，指数上的分母比较大，而这是一个负指数，所以整个接受函数可能会趋于1，即比当前解x，更差的新解工，也可能被接受，因此就有可能跳出局部极小而进行广域搜索，去搜索解空间的其他区域；而随着冷却的进行，T减小到一个比较小的值时，接受函数分母小了，整体也小了，即难以接受比当前解更差的解，也就是不太容易跳出当前的区域。如果在高温时，已经进行了充分的广域搜索，找到了可能存在最好解的区域，而在低温再进行足够的局部搜索，则可能最终找到全局最优了。因此，一般T，应设为一个足够小的正数，比如0.01~5，但这只是一个粗糙的经验，更精细的设置及其他的终止准则可以查阅文献。

设计要点

为了更好地实现模拟退火算法，还需要注意以下方面。

状态表达

上文已经提到过，SA算法中优化问题的一个解模拟了(或说可以想象为)退火过程中固体内部的一种粒子分布情况。这里状态表达即指实际问题的解(即状态)，如何以一种合适的数学形式被表达出来，它应当适用于SA的求解、又能充分表达实际问题，这需要仔细地设计。

可以参考遗传算法和禁忌搜索中编码的相关内容。常见的表达方式有：背包问题和指派问题的0-1编码， TSP问题和调度问题的自然数编码：还有用于连续函数优化的实数编码等。

新解的产生

新解产生机制的基本要求是能够尽量遍及解空间的各个区域，这样、在某一恒定温度不断产生新解时，就可能跳出当前区域的极小以搜索其他区域，这是模拟退火算法能够进行广域搜索的一个重要条件。

收敛的一般性条件

收敛到全局最优的一般性条件是：

①初始温度足够高：
②热平衡时间足够长；
③终止温度足够低；
④降温过程足够缓慢。但上述条件在应用中很难同时满足。

Python实现

函数： $f(x) = (x^2 - 5x)sin(x^2)$

import numpy as np
import matplotlib.pyplot as plt
import randomclass SA(object):def __init__(self, interval, tab='min', T_max=10000, T_min=1, iterMax=1000, rate=0.95):self.interval = interval                                    # 给定状态空间 - 即待求解空间self.T_max = T_max                                          # 初始退火温度 - 温度上限self.T_min = T_min                                          # 截止退火温度 - 温度下限self.iterMax = iterMax                                      # 定温内部迭代次数self.rate = rate                                            # 退火降温速度#############################################################self.x_seed = random.uniform(interval[0], interval[1])      # 解空间内的种子self.tab = tab.strip()                                      # 求解最大值还是最小值的标签: 'min' - 最小值；'max' - 最大值#############################################################self.solve()                                                # 完成主体的求解过程self.display()                                              # 数据可视化展示def solve(self):temp = 'deal_' + self.tab                                   # 采用反射方法提取对应的函数if hasattr(self, temp):deal = getattr(self, temp)else:exit('>>>tab标签传参有误："min"|"max"<<<')x1 = self.x_seedT = self.T_maxwhile T >= self.T_min:for i in range(self.iterMax):f1 = self.func(x1)delta_x = random.random() * 2 - 1if x1 + delta_x >= self.interval[0] and x1 + delta_x <= self.interval[1]:   # 将随机解束缚在给定状态空间内x2 = x1 + delta_xelse:x2 = x1 - delta_xf2 = self.func(x2)delta_f = f2 - f1x1 = deal(x1, x2, delta_f, T)T *= self.rateself.x_solu = x1                                            # 提取最终退火解def func(self, x):                                              # 状态产生函数 - 即待求解函数value = np.sin(x**2) * (x**2 - 5*x)return valuedef p_min(self, delta, T):                                      # 计算最小值时，容忍解的状态迁移概率probability = np.exp(-delta/T)return probabilitydef p_max(self, delta, T):probability = np.exp(delta/T)                               # 计算最大值时，容忍解的状态迁移概率return probabilitydef deal_min(self, x1, x2, delta, T):if delta < 0:                                               # 更优解return x2else:                                                       # 容忍解P = self.p_min(delta, T)if P > random.random(): return x2else: return x1def deal_max(self, x1, x2, delta, T):if delta > 0:                                               # 更优解return x2else:                                                       # 容忍解P = self.p_max(delta, T)if P > random.random(): return x2else: return x1def display(self):print('seed: {}\nsolution: {}'.format(self.x_seed, self.x_solu))plt.figure(figsize=(6, 4))x = np.linspace(self.interval[0], self.interval[1], 300)y = self.func(x)plt.plot(x, y, 'g-', label='function')plt.plot(self.x_seed, self.func(self.x_seed), 'bo', label='seed')plt.plot(self.x_solu, self.func(self.x_solu), 'r*', label='solution')plt.title('solution = {}'.format(self.x_solu))plt.xlabel('x')plt.ylabel('y')plt.legend()plt.savefig('SA.png', dpi=500)plt.show()plt.close()if __name__ == '__main__':SA([-5, 5], 'max')