深度学习经典算法 | 模拟退火算法详解

模拟退火算法基本思想

 

 

其他一些参数的说明

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

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

算法基本步骤

 

 

几点说明

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

状态表达

上文已经提到过，SA算法中优化问题的一个解模拟了(或说可以想象为)退火过程中固体内部的一种粒子分布情况。这里状态表达即指实际问题的解(即状态)如何以一种合适的数学形式被表达出来，它应当适用于SA的求解、又能充分表达实际问题，这需要仔细地设计。可以参考遗传算法和禁忌搜索中编码的相关内容。常见的表达方式有：背包问题和指派问题的0-1编码， TSP问题和调度问题的自然数编码：还有用于连续函数优化的实数编码等。

新解的产生

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

收敛的一般性条件

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

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

参数的选择

(1)控制参数Ｔ的初值T。

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

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

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

衰减函数可以有多种形式，一个常用的衰减函数是

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

(3) Markov链长度

 

 

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')

 

结果展示

 

 

 

参考文献

《matlab在数学建模中的应用》