模拟退火算法(1) - Python实现

• 抽象来源：模仿冶金过程中的退火原理。

• 核心思想：在冶金退火过程中，随着温度的下降，系统内部分子的平均动能逐渐降低，分子在自身位置附近的扰动能力也随之下降，即分子自身的搜索范围随着温度的下降而下降。利用该特性，我们可以对给定状态空间（待求解空间）内的某个状态产生函数（待求解函数）的最值进行求解。在高温状态下，由于分子的扰动能力较强，对较差状态（远离最值所对应的状态）的容忍性高，因此可以在给定状态空间内进行全局的随机搜索，从而有较高概率跳出局部极值。随着温度的逐渐下降，分子的扰动能力减弱，对较差状态的容忍性随之降低，导致此时的全局随机搜索能力下降，相应地对局部极值的搜索能力上升。综合整个退火过程，在理想情况下，最终解应该对应于给定状态空间内的最值。

• 迭代公式：
  系统温度为$T$时，出现能量差为$dE$的降温概率为：
  $$\begin{array}{}\text{(1)}& p(dE)=e^{dE/kT}\end{array}$$
  其中，$dE<0$，因此$0<p(dE)<1$。
  以该降温概率为基础，采用状态转移概率$p(\mathrm{\Delta }f)$来表示对较差状态的容忍性：
  $$\begin{array}{}\text{(2)}& p(\mathrm{\Delta }f)=e^{\mathrm{\Delta }f/T}\end{array}$$
  其中，$f$为状态产生函数，因此$\mathrm{\Delta }f=f(x_{new})-f(x_{old})$。常数$k$通过改变温度$T$的取值范围而被忽略，即式$\text{1}$中的$kT$等效于式$\text{2}$中的$T$。
  注意，在求取最值的过程中，存在两类可接受解：1）更优解 --- $BS$（较当前状态更加接近最值的状态）；2）容忍解 --- $TS$（较当前状态更加远离最值的解）。前者的接受概率1；后者的接受概率为状态转移概率$p(\mathrm{\Delta }f)$，且必须保证$p(\mathrm{\Delta }f)\in (0,1)$。
  • 求取最小值：
    $$\begin{array}{}\text{(3)}& \{\begin{array}{l}\mathrm{\Delta }f<0\to BS\to 1\\ \mathrm{\Delta }f\ge 0\to TS\to p(-\mathrm{\Delta }f)\end{array}\end{array}$$
  • 求取最大值：
    $$\begin{array}{}\text{(4)}& \{\begin{array}{l}\mathrm{\Delta }f>0\to BS\to 1\\ \mathrm{\Delta }f\le 0\to TS\to p(\mathrm{\Delta }f)\end{array}\end{array}$$

 

• Python代码实现：
  

  import numpy as np
  import matplotlib.pyplot as plt
  import random
  
  class 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_seed
          T = self.T_max
          while T >= self.T_min:
              for i in range(self.iterMax):
                  f1 = self.func(x1)
                  delta_x = random.random() * 2 - 1
                  if x1 + delta_x >= self.interval[0] and x1 + delta_x <= self.interval[1]:   # 将随机解束缚在给定状态空间内
                      x2 = x1 + delta_x
                  else:
                      x2 = x1 - delta_x
                  f2 = self.func(x2)
                  delta_f = f2 - f1
                  x1 = deal(x1, x2, delta_f, T)
              T *= self.rate
          self.x_solu = x1                                            # 提取最终退火解       
          
      def func(self, x):                                              # 状态产生函数 - 即待求解函数
          value = np.sin(x**2) * (x**2 - 5*x)
          return value
          
      def p_min(self, delta, T):                                      # 计算最小值时，容忍解的状态迁移概率
          probability = np.exp(-delta/T)
          return probability
          
      def p_max(self, delta, T):
          probability = np.exp(delta/T)                               # 计算最大值时，容忍解的状态迁移概率
          return probability
          
      def deal_min(self, x1, x2, delta, T):
          if delta < 0:                                               # 更优解
              return x2
          else:                                                       # 容忍解
              P = self.p_min(delta, T)
              if P > random.random(): return x2
              else: return x1
              
      def deal_max(self, x1, x2, delta, T):
          if delta > 0:                                               # 更优解
              return x2
          else:                                                       # 容忍解
              P = self.p_max(delta, T)
              if P > random.random(): return x2
              else: return x1
          
      def 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')

  笔者所用示例函数为 ：
  $$\begin{array}{}\text{(5)}& f(x)=(x^2-5x)sin(x^2)\end{array}$$
  
  

• 结果展示：
  
  
  

• 参考：https://blog.csdn.net/AI_BigData_wh/article/details/77943787?locationNum=2&fps=1