模拟退火算法(1) - Python实现
- 抽象来源：模仿冶金过程中的退火原理。
- 核心思想：在冶金退火过程中，随着温度的下降，系统内部分子的平均动能逐渐降低，分子在自身位置附近的扰动能力也随之下降，即分子自身的搜索范围随着温度的下降而下降。利用该特性，我们可以对给定状态空间（待求解空间）内的某个状态产生函数（待求解函数）的最值进行求解。在高温状态下，由于分子的扰动能力较强，对较差状态（远离最值所对应的状态）的容忍性高，因此可以在给定状态空间内进行全局的随机搜索，从而有较高概率跳出局部极值。随着温度的逐渐下降，分子的扰动能力减弱，对较差状态的容忍性随之降低，导致此时的全局随机搜索能力下降，相应地对局部极值的搜索能力上升。综合整个退火过程，在理想情况下，最终解应该对应于给定状态空间内的最值。
- 迭代公式：
系统温度为$T$时，出现能量差为$dE$的降温概率为：
\begin{equation}\label{eq_1}
p(dE) = e^{dE/kT}
\end{equation}
其中，$dE<0$，因此$0<p(dE)<1$。
以该降温概率为基础，采用状态转移概率$p(\Delta f)$来表示对较差状态的容忍性：
\begin{equation}\label{eq_2}
p(\Delta f) = e^{\Delta f/T}
\end{equation}
其中，$f$为状态产生函数，因此$\Delta f = f(x_{new}) - f(x_{old})$。常数$k$通过改变温度$T$的取值范围而被忽略，即式\ref{eq_1}中的$kT$等效于式\ref{eq_2}中的$T$。
注意，在求取最值的过程中，存在两类可接受解：1）更优解 --- $BS$（较当前状态更加接近最值的状态）；2）容忍解 --- $TS$（较当前状态更加远离最值的解）。前者的接受概率1；后者的接受概率为状态转移概率$p(\Delta f)$，且必须保证$p(\Delta f)\in (0, 1)$。- 求取最小值：
\begin{equation}
\begin{cases}
\Delta f < 0 \rightarrow BS \rightarrow 1\\
\Delta f \geq 0 \rightarrow TS \rightarrow p(-\Delta f)
\end{cases}
\end{equation} - 求取最大值：
\begin{equation}
\begin{cases}
\Delta f > 0 \rightarrow BS \rightarrow 1\\
\Delta f \leq 0 \rightarrow TS \rightarrow p(\Delta f)
\end{cases}
\end{equation}
- 求取最小值：
- Python代码实现：
View Code1 import numpy as np 2 import matplotlib.pyplot as plt 3 import random 4 5 class SA(object): 6 7 def __init__(self, interval, tab='min', T_max=10000, T_min=1, iterMax=1000, rate=0.95): 8 self.interval = interval # 给定状态空间 - 即待求解空间 9 self.T_max = T_max # 初始退火温度 - 温度上限 10 self.T_min = T_min # 截止退火温度 - 温度下限 11 self.iterMax = iterMax # 定温内部迭代次数 12 self.rate = rate # 退火降温速度 13 ############################################################# 14 self.x_seed = random.uniform(interval[0], interval[1]) # 解空间内的种子 15 self.tab = tab.strip() # 求解最大值还是最小值的标签: 'min' - 最小值；'max' - 最大值 16 ############################################################# 17 self.solve() # 完成主体的求解过程 18 self.display() # 数据可视化展示 19 20 def solve(self): 21 temp = 'deal_' + self.tab # 采用反射方法提取对应的函数 22 if hasattr(self, temp): 23 deal = getattr(self, temp) 24 else: 25 exit('>>>tab标签传参有误："min"|"max"<<<') 26 x1 = self.x_seed 27 T = self.T_max 28 while T >= self.T_min: 29 for i in range(self.iterMax): 30 f1 = self.func(x1) 31 delta_x = random.random() * 2 - 1 32 if x1 + delta_x >= self.interval[0] and x1 + delta_x <= self.interval[1]: # 将随机解束缚在给定状态空间内 33 x2 = x1 + delta_x 34 else: 35 x2 = x1 - delta_x 36 f2 = self.func(x2) 37 delta_f = f2 - f1 38 x1 = deal(x1, x2, delta_f, T) 39 T *= self.rate 40 self.x_solu = x1 # 提取最终退火解 41 42 def func(self, x): # 状态产生函数 - 即待求解函数 43 value = np.sin(x**2) * (x**2 - 5*x) 44 return value 45 46 def p_min(self, delta, T): # 计算最小值时，容忍解的状态迁移概率 47 probability = np.exp(-delta/T) 48 return probability 49 50 def p_max(self, delta, T): 51 probability = np.exp(delta/T) # 计算最大值时，容忍解的状态迁移概率 52 return probability 53 54 def deal_min(self, x1, x2, delta, T): 55 if delta < 0: # 更优解 56 return x2 57 else: # 容忍解 58 P = self.p_min(delta, T) 59 if P > random.random(): return x2 60 else: return x1 61 62 def deal_max(self, x1, x2, delta, T): 63 if delta > 0: # 更优解 64 return x2 65 else: # 容忍解 66 P = self.p_max(delta, T) 67 if P > random.random(): return x2 68 else: return x1 69 70 def display(self): 71 print('seed: {}\nsolution: {}'.format(self.x_seed, self.x_solu)) 72 plt.figure(figsize=(6, 4)) 73 x = np.linspace(self.interval[0], self.interval[1], 300) 74 y = self.func(x) 75 plt.plot(x, y, 'g-', label='function') 76 plt.plot(self.x_seed, self.func(self.x_seed), 'bo', label='seed') 77 plt.plot(self.x_solu, self.func(self.x_solu), 'r*', label='solution') 78 plt.title('solution = {}'.format(self.x_solu)) 79 plt.xlabel('x') 80 plt.ylabel('y') 81 plt.legend() 82 plt.savefig('SA.png', dpi=500) 83 plt.show() 84 plt.close() 85 86 87 if __name__ == '__main__': 88 SA([-5, 5], 'max')
笔者所用示例函数为 ：
\begin{equation}
f(x) = (x^2 - 5x)sin(x^2)
\end{equation}
- 结果展示：
- 参考：https://blog.csdn.net/AI_BigData_wh/article/details/77943787?locationNum=2&fps=1
