启发式算法学习（二）：模拟退火算法SAA

算法简介

模拟退火算法（Simulated annealing algorithm，SAA）的思想借鉴于固体的退火原理，当固体的温度很高的时候，内能比较大，固体的内部粒子处于快速无序运动，当温度慢慢降低的过程中，固体的内能减小，粒子的慢慢趋于有序，最终，当固体处于常温时，内能达到最小，此时，粒子最为稳定。模拟退火算法便是基于这样的原理设计而成。

模拟退火算法从某一高温出发，在高温状态下计算初始解，然后以预设的邻域函数产生一个扰动量，从而得到新的状态，即模拟粒子的无序运动，比较新旧状态下的能量，即目标函数的解。如果新状态的能量小于旧状态，则状态发生转化；如果新状态的能量大于旧状态，则以一定的概率准则发生转化。当状态稳定后，便可以看作达到了当前状态的最优解，便可以开始降温，在下一个温度继续迭代，最终达到低温的稳定状态，便得到了模拟退火算法产生的结果。

算法过程

状态空间与邻域函数

状态空间也称为搜索空间，它由经过编码的可行解的集合所组成。而邻域函数应尽可能满足产生的候选解遍布全部状态空间。其通常由产生候选解的方式和候选解产生的概率分布组成。候选解一般按照某一概率密度函数对解空间进行随机采样获得，而概率分布可以为均匀分布、正态分布、指数分布等。

状态概率分布（Metropolis准则）

状态转移概率是指从一个状态转换成另一个状态的概率，模拟退火算法中一般采用Metropolis准则，具体如下：

$$f(x)=\left\lbrace\begin{array}{cll}
1 & , & E(x_{new}) < E(x_{old}) \\
exp(-\frac{E(x_{new})-E(x_{old})}{T}) & , & E(x_{new}) \ge E(x_{old})
\end{array}\right.$$

其与当前温度参数T有关，随温度的下降而减小。

冷却进度表

冷却进度表是指从某一高温状态T向低温状态冷却时的降温函数,设时刻的温度为T(t)，则经典模拟退火算法的降温方式为：

$$T(t)=\frac{T_{0}}{lg(1+t)}$$

而快速模拟退火算法的降温方式为：

$$T(t)=\frac{T_{0}}{1+t}$$

其他方法不再赘述。

初始温度

一般来说，初始温度越大，获得高质量解的几率越大，但是花费的时间也会随之增加，因此，初温的确定应该同时考虑计算效率与优化质量，常用的方法包括：

1.均匀抽样一组状态，以各状态目标值的方差为初温。

2.随机产生一组状态，确定两两状态间的最大目标值差，然后根据差值，利用一定的函数确定初温，如： $T_{0}=-\frac{\Delta_{max}}{Pr}$ ,其中Pr为初始接受概率。

3.根据经验公式给出。

循环终止准则

内循环（求解循环）终止准则：

1. 检验目标函数的均值是否稳定
2. 连续若干步的目标值变化较小
3. 按一定的步数进行抽样

外循环（降温循环）终止准则：

1. 设置终止温度
2. 设置外循环迭代次数
3. 算法搜索到的最优值连续若干步保持不变
4. 检验系统熵是否稳定

Python实现

实例函数： $f(x)=(x^{2}-5x)sin(x^2)$

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 # [-1,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')

Reference
https://www.imooc.com/article/30160
https://baike.baidu.com/item/模拟退火算法/355508?fr=aladdin
https://blog.csdn.net/google19890102/article/details/45395257
https://www.cnblogs.com/xxhbdk/p/9192750.htmlst=>start: Start op1=>operation: 随机生成初始解w，计算目标函数f(w) op2=>operation: 扰动产生新解w'，计算目标函数f(w') op3=>operation: 接受新解，w=w'，f(w)=f(w') op4=>operation: 按照Metropolis准则接受新解 op5=>operation: 缓慢降低温度，重置迭代次数 op6=>operation: 运算结束，返回最优解 cond1=>condition: f(w')更接近优化目标？ cond2=>condition: 是否达到迭代次数? cond3=>condition: 满足终止条件? ed=>end: End st->op1->op2->cond1(yes)->op3->cond2(no)->cond3(yes)->op6->ed cond1(no)->op4->cond2 cond2(yes)->op6->ed cond3(no)->op5(top)->op2{"scale":1,"line-width":2,"line-length":50,"text-margin":10,"font-size":12}