模拟退火算法

~hello world~

已于 2023-03-07 22:56:40 修改

阅读量2.1k

点赞数 1

分类专栏：模型算法文章标签：模拟退火算法算法机器学习

于 2023-02-23 21:27:28 首次发布

本文链接：https://blog.csdn.net/qq_50086023/article/details/129109165

版权

模型算法专栏收录该内容

9 篇文章 2 订阅

订阅专栏

1、模拟退火算法介绍

模拟退火算法 (Simulate Anneal Arithmetic，SAA) 是一种通用概率演算法，用来在一个大的搜寻空间内找寻命题的最优解。模拟退火的原理和金属退火的原理近似：将热力学的理论套用到统计学上，将搜寻空间内每一点想像成空气内的分子；分子的能量，就是它本身的动能；而搜寻空间内的每一点，也像空气分子一样带有“能量”，以表示该点对命题的合适程度。演算法先以搜寻空间内一个任意点作起始：每一步先选择一个“邻居”，然后再计算从现有位置到达“邻居”的概率。模拟退火算法是通过赋予搜索过程一种时变且最终趋于零的概率突跳性，从而可有效避免陷入局部极小并最终趋于全局最优的串行结构的优化算法。

2、模拟退火算法的理解(图解)

假如有一个函数 y=f(x)，其图行如下图所示。
在这里插入图片描述

现在，想找到这个函数的最大值。模拟退火算法是这么认为的：先在 x 的定义域内，取一个起始点 x=i，如图红色虚线，得到 y=f(i) 。接着让红色线往右移动一个单位，然后判定：如果 f（i + 1） > f（i），则接收该移动。可以简单理解为发现移动到的新点求出来的值更大，则接收这个移动。现在，红色线在不断重复上述操作的过程中，会到达函数的第一个极大值。如下图：
在这里插入图片描述
现在，再进行上述步骤，会发现 f（i + 1） < f（i），如果按照上述规则那么这里将不再接受移动了。可是，我们知道只要这个点不是真正的极大值，但电脑是不会知道的，那么该怎么做呢？模拟退火算法的核心就在这里，它不会放弃这个新的 f(i + 1)，它允许 f(i + 1) 可以比 f(i) 小，但并不是始终允许，而是看经过一段时间观察 f(i + 1) 能不能蜕变成“最大值”。因此对于 f(i + 1) < f(i) ，会执行如下操作：如果 f(i + 1) < f(i)，以一定概率接收该移动，并随时间推移概率降低。而这个概率，就是参考了金属冶炼的退火过程，被称为Metropolis准则。

3、Metropolis准则

在金属冶炼的退火过程中，在温度为 T 时，出现能量差为 dE 的粒子的降温的概率为P(dE)，表示为：在这里插入图片描述
其中: k 是波尔兹曼常数，值为 k=1.3806488(13)×10⁻²³，exp表示自然指数，dE<0， dE/kT<0，p(dE)函数的取值范围是(0,1)，满足概率密度函数的定义。其实这条公式更直观意思就是：温度越高，出现一次能量差为p(dE)的降温的概率就越大；温度越低，则出现降温的概率就越小。
　　
在实际问题中，假定当前可行解为 x，迭代更新后的解为 x_new，那么对应的“能量差”定义为： Δf = f(x_new) − f(x)，则其对应的“一定概率”为：
在这里插入图片描述
注：在实际问题中，可以设定k=1。因为kT可以等价于一个参数 T。

4、模拟退火算法模型

4.1 模拟退火的基本思想

(1) 初始化：初始温度T(充分大)，初始解状态S(是算法迭代的起点)， 每个T值的迭代次数L               
(2) 对k=1，……，L做第(3)至第6步：
(3) 产生新解S′
(4) 计算增量Δt′=C(S′)-C(S)，其中C(S)为评价函数||优化目标
(5) 若Δt′<0则接受S′作为新的当前解，否则以概率exp(-Δt′/T)接受S′作为新的当前解.
(6) 如果满足终止条件则输出当前解作为最优解，结束程序。终止条件通常取为连续若干个新解都没有被接受时终止算法。
(7) T逐渐减少，且T->0，然后转第2步。

4.2 算法对应动态演示图

在这里插入图片描述
模拟退火算法新解的产生和接受可分为如下四个步骤：
第一步是由一个产生函数从当前解产生一个位于解空间的新解
第二步是计算与新解所对应的目标函数差。
第三步是判断新解是否被接受，判断的依据是一个接受准则，最常用的接受准则是Metropolis准则: 若Δt′<0则接受S′作为新的当前解S，否则以概率exp(-Δt′/T)接受S′作为新的当前解S。
第四步是当新解被确定接受时，用新解代替当前解，这只需将当前解中对应于产生新解时的变换部分予以实现，同时修正目标函数值即可。此时，当前解实现了一次迭代。可在此基础上开始下一轮试验。而当新解被判定为舍弃时，则在原当前解的基础上继续下一轮试验。

模拟退火算法与初始值无关，算法求得的解与初始解状态S(是算法迭代的起点)无关；模拟退火算法具有渐近收敛性，已在理论上被证明是一种以概率l收敛于全局最优解的全局优化算法。

5、模拟退火算法的参数

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

（1）控制参数的初值T0：冷却开始的温度。
（2）控制参数T的衰减函数：因计算机能够处理的都是离散数据，因此把连续的降温过程离散化成降温过程中的一系列温度点，衰减函数即计算这一系列温度的表达式。
（3）控制参数T的终值Tf (停止准则)。
（4）Markov链的长度Lk：任一温度T 的迭代次数。

5.1 参数设置

模拟退火算法的应用很广泛，可以求解NP完全问题，但其参数难以控制，其主要问题有以下三点：
（1）控制参数T的初值T0：求解全局优化问题的随机搜索算法一般都采用大范围的粗略搜索与局部的精细搜索相结合的搜索策略。只有在初始的大范围搜索阶段找到全局最优解所在的区域，才能逐步缩小搜索的范围，最终求出全局最优解。在问题规模较大时，过小的T0往往导致算法难以跳出局部陷阱而达不到全局最优。一般为100℃。【初始温度高，则搜索到全局最优解的可能性大，但因此要花费大量的计算时间；反之，则可节约计算时间，但全局搜索性能可能受到影响。实际应用过程中，初始温度一般需要依据实验结果进行若干次调整。】
（2）控制参数 T 的衰减函数：衰减函数可以有多种形式，一个常用的衰减函数是：
在这里插入图片描述
其中，小k为降温的次数，α 是一个常数，可以取为0.5~0.99，它的取值决定了降温过程的快慢。【为了保证较大的搜索空间，α一般取接近于1的值，如0.95、0.9】

（3）退火速度：Markov链长度。Markov链长度的选取原则是：在控制参数 T 的衰减函数已经选定的前提下，L_k 应能使在控制参数T的每一取值上达到准平衡。从经验上说，对简单的情况可以令L_k =100n，n为问题规模。【循环次数增加必定带来计算开销的增大，实际应用中，要针对具体问题的性质和特征设置合理的退火平衡条件。】

（4）控制参数T的终值 T_f (停止准则)。或者若干个相继的Mapkob链解没有产生变化，就是连续好多个解都没变化。算法停止准则：对Metropolis准则中的接受函数在这里插入图片描述
分析可知，在T比较大的高温情况，指数上的分母比较大，而这是一个负指数，所以整个接受函数可能会趋于1，即比当前解 x_i 更差的新解 x_j 也可能被接受，因此就有可能跳出局部极小而进行广域搜索，去搜索解空间的其他区域；而随着冷却的进行，T减小到一个比较小的值时，接收函数分母小了整体也小了，即较难接受比当前解更差的解，也就不太容易跳出当前的区域。如果在高温时，已经进行了充分的广域搜索，找到了可能存在的最好解的区域，而在低温再进行足够的局部搜索，则可能最终找到全局最优了。因此，一般终止温度 T_j 应设为足够小的正数，比如0.01~0.05。

完整代码

# 模拟退火算法

import matplotlib.pyplot as plt
import math
import random
 
"""
求解的目标表达式为：y = 10 * math.sin(5 * x) + 7 * math.cos(4 * x)  x belongs to (0,10)
"""
 
 
def main():
    plot_obj_func()
    T_init = 100  # 初始最大温度
    alpha = 0.90  # 降温系数
    T_min = 1e-3  # 最小温度，即退出循环条件
    T = T_init
    x = random.random() * 10  # 初始化x，生成一个 0 到 10 之间的随机小数，并将其赋值给变量 x。
    y = 10 * math.sin(5 * x) + 7 * math.cos(4 * x)
    results = []  # 存x，y
    while T > T_min:
        x_best = x
        # y_best = float('-inf')  # 设置成这个有可能会陷入局部最优，不一定全局最优
        y_best = y  # 设置成这个收敛太快了，令人智熄
        flag = 0  # 用来标识该温度下是否有新值被接受
        # 每个温度迭代50次，找最优解
        for i in range(50):
            delta_x = random.random() - 0.5  # 自变量进行波动
            # 自变量变化后仍要求在[0,10]之间
            if 0 < (x + delta_x) < 10:
                x_new = x + delta_x
            else:
                x_new = x - delta_x
            y_new = 10 * math.sin(5 * x_new) + 7 * math.cos(4 * x_new)
            # 要接受这个y_new为当前温度下的理想值，要满足
            # 1y_new>y_old
            # 2math.exp(-(y_old-y_new)/T)>random.random()
            # 以上为找最大值，要找最小值就把>号变成<
            if (y_new > y or math.exp(-(y - y_new) / T) > random.random()):
                flag = 1  # 有新值被接受
                x = x_new
                y = y_new
                if y > y_best:
                    x_best = x
                    y_best = y
        if flag:
            x = x_best
            y = y_best
        results.append((x, y))
        T *= alpha # 降温
 
    print('最优解 x:%f,y:%f' % results[-1])
 
    plot_final_result(results)
    plot_iter_curve(results)
 
 
 
# 看看我们要处理的目标函数
def plot_obj_func():
    """y = 10 * math.sin(5 * x) + 7 * math.cos(4 * x)"""
    X1 = [i / float(10) for i in range(0, 100, 1)] # 创建列表 X1，其中包含了从0到99的所有整数，得到的每个整数都除以10得到的浮点数。
    Y1 = [10 * math.sin(5 * x) + 7 * math.cos(4 * x) for x in X1]
    plt.plot(X1, Y1)
    plt.show()
 
 
# 看看最终的迭代变化曲线
def plot_iter_curve(results):
    X = [i for i in range(len(results))]
    Y = [results[i][1] for i in range(len(results))]
    plt.plot(X, Y)
    plt.show()
 
def plot_final_result(results):
    X1 = [i / float(10) for i in range(0, 100, 1)]
    Y1 = [10 * math.sin(5 * x) + 7 * math.cos(4 * x) for x in X1]
    plt.plot(X1, Y1)
    plt.scatter(results[-1][0], results[-1][1], c='r', s=10)
    plt.show()
 
if __name__ == '__main__':
    # for i in range(100):
    main()