526互联

数学建模（三）：模拟退火算法（SA）

发布时间 2023-04-04 12:50:33作者: Steve_Anthony

目录

模拟退火算法（SA）

一、 概述

1、 算法简介

模拟退火算法（simulated annealing，SA）来源于固体退火原理，是一种基于概率的算法。

模拟退火算法(SA)来源于固体退火原理，是一种基于概率的算法。将固体加温至充分高的温度，再让其徐徐冷却，加温时，固体内部粒子随温升变为无序状，内能增大，分子和原子越不稳定。而徐徐冷却时粒子渐趋有序，能量减少，原子越稳定。在冷却（降温）过程中，固体在每个温度都达到平衡态，最后在常温时达到基态，内能减为最小。

模拟退火算法从某一较高初温出发，伴随温度参数的不断下降,结合概率突跳特性在解空间中随机寻找目标函数的全局最优解，即在局部最优解能概率性地跳出并最终趋于全局最优。模拟退火算法是通过赋予搜索过程一种时变且最终趋于零的概率突跳性，从而可有效避免陷入局部极小并最终趋于全局最优的串行结构的优化算法。

2、 核心思想

模拟退火算法从某一较高初温出发，伴随温度参数的不断下降,结合一定的概率突跳特性在解空间中随机寻找目标函数的全局最优解，即在局部最优解能概率性地跳出并最终趋于全局最优。

这里的 “一定的概率” 的计算参考了金属冶炼的退火过程，这也是模拟退火算法名称的由来。将温度 T 当作控制参数，目标函数值 f 视为内能 E ，而固体在某温度 T 时的一个状态对应一个解 \(x_i\)，然后算法试图随着控制参数 T 的降低，使目标函数 f (内能 E )也逐渐降低，直至趋于全局最小值（退火中低温时的最低能量状态），就像金属退火过程一样。

关于爬山算法与模拟退火，有一个有趣的比喻：

  • 爬山算法：兔子朝着比现在高的地方跳去。它找到了不远处的最高山峰。但是这座山不一定是珠穆朗玛峰。这就是爬山算法，它不能保证局部最优值就是全局最优值。

  • 模拟退火：兔子喝醉了。它随机地跳了很长时间。这期间，它可能走向高处，也可能踏入平地。但是，它渐渐清醒了并朝最高方向跳去。这就是模拟退火。

3、 数学原理

从上面我们知道，会结合概率突跳特性在解空间中随机寻找目标函数的全局最优解，那么具体的更新解的机制是什么呢？如果新解比当前解更优，则接受新解，否则基于Metropolis准则判断是否接受新解。接受概率为：

假设开始状态在A，随着迭代次数更新到B局部最优解，这时发现更新到B时，能量比A要低，则说明接近最优解了，因此百分百转移，状态到达B后，发现下一步能量上升了，如果是梯度下降则是不允许继续向前的，而这里会以一定的概率跳出这个坑，这个概率和当前的状态、能量等都有关系，如果B最终跳出来了到达C，又会继续以一定的概率跳出来，直到到达D后，就会稳定下来。

4、 模拟退火的流程

(1) 初始化：初始温度 T (充分大)，初始解状态S(是算法迭代的起点)，每个 T 值的迭代次数 L

(2) 对 k=1, …, L 做第(3)至第（6）步：

(3) 产生新解 S′

(4) 计算增量 ΔT = C(S′)-C(S)，其中 C(S) 为目标函数， C(S) 相当于能量

(5) 若 ΔT<0 则接受 S′ 作为新的当前解，否则以概率 exp(-ΔT/T) 接受S′作为新的当前解.

(6) 如果满足终止条件则输出当前解作为最优解，结束程序。

(7) T 逐渐减少，且 T->0 ，然后转第2步。

其中有几个需要注意的点：

  • 初始点的选取对算法结果有一定的影响，最好是多次运行对结果进行综合判断。
  • 在算法运行初期，温度下降快，避免接受过多的差结果。当运行时间增加，温度下降减缓，以便于更快稳定结果。
  • 当迭代次数增加到一定次数时，结果可能已经达到稳定，但是距离算法结束还有一段时间。在设计程序时应该加入适当的输出条件，满足输出条件即可结束程序。

可以大概分为这四个步骤：

  1. 第一步是由一个产生函数从当前解产生一个位于解空间的新解；为便于后续的计算和接受，减少算法耗时，通常选择由当前新解经过简单地变换即可产生新解的方法，如对构成新解的全部或部分元素进行置换、互换等，注意到产生新解的变换方法决定了当前新解的邻域结构，因而对冷却进度表的选取有一定的影响。
  2. 第二步是计算与新解所对应的目标函数差。因为目标函数差仅由变换部分产生，所以目标函数差的计算最好按增量计算。事实表明，对大多数应用而言，这是计算目标函数差的最快方法。
  3. 第三步是判断新解是否被接受,判断的依据是一个接受准则，最常用的接受准则是 Metropolis 准则: 若 ΔT < 0 则接受 S′ 作为新的当前解 S，否则以概率 P 接受 S′ 作为新的当前解 S。
  4. 第四步是当新解被确定接受时，用新解代替当前解，这只需将当前解中对应于产生新解时的变换部分予以实现，同时修正目标函数值即可。此时，当前解实现了一次迭代。可在此基础上开始下一轮试验。而当新解被判定为舍弃时，则在原当前解的基础上继续下一轮试验。

二、 实例分析

1、 初始化参数

# -*- coding: utf-8 -*-
"""
Created on Mon Apr  3 19:17:28 2023

@author: steve
"""
from random import random
import math
import matplotlib.pyplot as plt


max_iter = 100  # 每一次温度降低的迭代次数
alpha = 0.99  # 降温系数
T_f = 0.01  # 温度的终值
T_n = 100  # 当前的温度，也是初始温度
x, y = [random() * 10 - 5 for i in range(max_iter)], [random() * 10 - 5 for i in range(max_iter)] # 进行数据的初始化
f = lambda x, y : (4 * x ** 2 - 2.1 * x ** 4 + x ** 6 / 3 + x * y - 4 * y ** 2 + 4 * y ** 4)  # 我们需要求的函数
result = {
    "f": [], 
    "t": []
    }  # 用来存放每一次下降的最优解

2、 Metrospolis 准则

def metrospolis(T_n, old, new):
    """
    进行 Metrospolis 准则的判断

    Parameters
    ----------
    T_n : int
        当前的温度.  
    old : double
        函数扰动之前的值.
    new : double
        函数扰动之后的值.

    Returns
    -------
    int
        是否需要重新寻找值.
    """
    if old >= new:
        return 1
    else:
        p = math.exp((old - new) / T_n)
        if random() < p:
            return 1
        else:
            return 0

3、 生成新的值

def generate_new(T_n, x, y):
    """
    其为扰动的过程

    Parameters
    ----------
    T_n : int
        当前的温度.
    x : double
        DESCRIPTION.
    y : double
        DESCRIPTION.

    Returns
    -------
    list
        返回生成的新值.
    """
    while True:
        x_new = x + T_n * (random() - random())
        y_new = y + T_n * (random() - random())
        if (-5 <= x_new <= 5) & (-5 <= y_new <= 5):  
            break                                  #重复得到新解，直到产生的新解满足约束条件
    return x_new, y_new

4、 获取最优值

def best(max_iter, f, x, y):
    """
    计算这个温度下的最优值

    Parameters
    ----------
    max_iter : int
        最大的迭代次数.
    f : function
        我们需要求的函数.
    x : double
        DESCRIPTION.
    y : double
        DESCRIPTION.

    Returns
    -------
    list
        返回最优值，以及它的索引.

    """
    min_val, min_inx = f(x[0], y[0]), 0
    for i in range(max_iter):
        val = f(x[i], y[i])
        if min_val > val:
            min_val, min_inx = val, i
            
    return [min_val, min_inx]

5、 主程序

def plot(result):
    plt.plot(result['t'], result['f'])
    plt.title('SA')
    plt.xlabel('t')
    plt.ylabel('f')
    plt.gca().invert_xaxis()
    plt.show()


def main(max_iter, alpha, T_f, T_n, x, y, f, result):
    count = 0  # 统计迭代了多少次
    while T_n > T_f:  # 外层循环，当前温度小于最低温度时，终止循环
        for i in range(max_iter):  # 内层循环
            x_new, y_new = generate_new(T_n, x[i], y[i])  # 产生新值
            if metrospolis(T_n, f(x[i], y[i]), f(x_new, y_new)):  # 将原值与新产生的值进行比较
                x[i] = x_new  # 如果接收新值，则存入数组中
                y[i] = y_new
        # 迭代 max_iter 后，记录该温度下的最优解
        [ft, _] = best(max_iter, f, x, y)
        result["f"].append(ft)
        result["t"].append(T_n)
        T_n = T_n * alpha  # 进行降温操作
        count += 1
        
    # 得到最优解
    [f_best, inx] = best(max_iter, f, x, y)
    print(f"F={f_best}, x_1={x[inx]}, x_2={y[inx]}")
    # 进行图像表示
    plot(result)

6、 总代码

# -*- coding: utf-8 -*-
"""
Created on Mon Apr  3 19:17:28 2023

@author: steve
"""
from random import random
import math
import matplotlib.pyplot as plt


def metrospolis(T_n, old, new):
    """
    进行 Metrospolis 准则的判断

    Parameters
    ----------
    T_n : int
        当前的温度.  
    old : double
        函数扰动之前的值.
    new : double
        函数扰动之后的值.

    Returns
    -------
    int
        是否需要重新寻找值.
    """
    if old >= new:
        return 1
    else:
        p = math.exp((old - new) / T_n)
        if random() < p:
            return 1
        else:
            return 0
        
        
def generate_new(T_n, x, y):
    """
    其为扰动的过程

    Parameters
    ----------
    T_n : int
        当前的温度.
    x : double
        DESCRIPTION.
    y : double
        DESCRIPTION.

    Returns
    -------
    list
        返回生成的新值.
    """
    while True:
        x_new = x + T_n * (random() - random())
        y_new = y + T_n * (random() - random())
        if (-5 <= x_new <= 5) & (-5 <= y_new <= 5):  
            break                                  #重复得到新解，直到产生的新解满足约束条件
    return x_new, y_new 
    
def best(max_iter, f, x, y):
    """
    计算这个温度下的最优值

    Parameters
    ----------
    max_iter : int
        最大的迭代次数.
    f : function
        我们需要求的函数.
    x : double
        DESCRIPTION.
    y : double
        DESCRIPTION.

    Returns
    -------
    list
        返回最优值，以及它的索引.

    """
    min_val, min_inx = f(x[0], y[0]), 0
    for i in range(max_iter):
        val = f(x[i], y[i])
        if min_val > val:
            min_val, min_inx = val, i
            
    return [min_val, min_inx]


def plot(result):
    plt.plot(result['t'], result['f'])
    plt.title('SA')
    plt.xlabel('t')
    plt.ylabel('f')
    plt.gca().invert_xaxis()
    plt.show()



def main(max_iter, alpha, T_f, T_n, x, y, f, result):
    count = 0  # 统计迭代了多少次
    while T_n > T_f:  # 外层循环，当前温度小于最低温度时，终止循环
        for i in range(max_iter):  # 内层循环
            x_new, y_new = generate_new(T_n, x[i], y[i])  # 产生新值
            if metrospolis(T_n, f(x[i], y[i]), f(x_new, y_new)):  # 将原值与新产生的值进行比较
                x[i] = x_new  # 如果接收新值，则存入数组中
                y[i] = y_new
        # 迭代 max_iter 后，记录该温度下的最优解
        [ft, _] = best(max_iter, f, x, y)
        result["f"].append(ft)
        result["t"].append(T_n)
        T_n = T_n * alpha  # 进行降温操作
        count += 1
        
    # 得到最优解
    [f_best, inx] = best(max_iter, f, x, y)
    print(f"F={f_best}, x_1={x[inx]}, x_2={y[inx]}")
    # 进行图像表示
    plot(result)
    


max_iter = 100  # 每一次温度降低的迭代次数
alpha = 0.99  # 降温系数
T_f = 0.01  # 温度的终值
T_n = 100  # 当前的温度，也是初始温度
x, y = [random() * 10 - 5 for i in range(max_iter)], [random() * 10 - 5 for i in range(max_iter)] # 进行数据的初始化
f = lambda x, y : (4 * x ** 2 - 2.1 * x ** 4 + x ** 6 / 3 + x * y - 4 * y ** 2 + 4 * y ** 4)  # 我们需要求的函数
result = {
    "f": [], 
    "t": []
    }  # 用来存放每一次下降的最优解

main(max_iter, alpha, T_f, T_n, x, y, f, result)

最后的运行结果为：