【数模/启发式算法】模拟退火算法

简介

模拟退火算法（Simulate Anneal Arithmetic，SAA）是一种通用概率演算法，用来在一个大的搜寻空间内找寻命题的最优解。模拟退火是S.Kirkpatrick, C.D.Gelatt和M.P.Vecchi在1983年所发明。而V.Černý在1985年也独立发明此演算法。

模拟退火算法是解决TSP问题的有效方法之一。

符号说明

符号	含义
$n$	粒子个数
$T0$	初始温度
$T_f$	最终温度
$\alpha$	降温速率
$T_i$	第i次迭代时的温度
$po{s}_i$	第i次迭代时自变量取值

核心思想

模拟退火来自冶金学的专有名词退火。

假设有一温度很高的物体，其在自然降温过程中会发生什么？当温度较高时，物体分子的热运动速率较大，分子较为活跃；当温度逐步下降时，物体分子的热运动速率较小，分子较为不活跃。

我们将上述过程带入到一个函数求极值问题中，并假设自变量即为分子。那么，温度较高时自变量取值发生变化的幅度较大；温度较低时自变量取值发生变化的幅度较小。另外，当自变量取到一个较好的取值时我们无条件接受取值；而当自变量取到一个不那么好甚至较差的取值时我们以一定的概率接受取值。因为，取值较差的自变量也有可能经过变化变化到最优解附近。

概率接受非较好取值，保证了算法跳出局部最优解避免陷入局部最优。

流程图

文章使用到的测试函数

一元函数：$$y=-11sin(x)-7cos(5x),x\in [-3,3]$$参考最小值为：$min(y)=-17.4905$

二元函数：$$y=x_1^2+x_2^2-x_1{x}_2-10x_1-4x_2+60$$参考最小值为：$min(y)=8.0$

基本原理

自变量变化控制

降温公式：$$T_{i+1}=\alpha T_i$$其中，$\alpha \in (0,1)$通常取0.9或0.99，即温度衰减很慢。$\alpha$取值越接近1（不等与1），温度衰减越慢、求解精度越高，但求解速度越久。

使用温度控制自变量变化步长：$$po{s}_{i+1}=po{s}_i+(Rand-0.5)\ast T_i$$其中，$Rand\in [0,1]$为随机数，$Rand-0.5$保证自变量变化正负均匀。
注意：此处的$po{s}_{i+1}$并不一定会被采纳，因为是概率取值。

概率接受非最优取值

注意：如果一个自变量取值较优，则会无条件接受。

概率接受非较好取值，保证了算法跳出局部最优解避免陷入局部最优。是算法的核心部分，其取值要求：变量距离较优解越远概率越小，温度越低概率越小。

概率计算公式：$$P=e^{\frac{△val}{T_i}}$$其中，$△val=oldval-newval$即自变量$po{s}_i$取值减去$po{s}_{i+1}$取值，$△val$由于$po{s}_{i+1}$的取值较$po{s}_i$的取值较差，因此$△val<0$。

$△val<0$保证了：变量距离较优解越远概率越小，温度越低概率越小。

模拟退火算法代码

Python版本：

import math
import random
from copy import deepcopy


def func(x):
    # 一元函数测试
    # return -11 * math.sin(x[0]) - 7 * math.cos(5 * x[0])
    # 二元函数测试
    return x[0] ** 2 + x[1] ** 2 - x[0] * x[1] - 10 * x[0] - 4 * x[1] + 60


class SA:
    def __init__(self, func, n, var, iter=100, T0=100, Tf=0.01, alpha=0.99, lb = None, ub = None):
        """ 默认寻找最小值，以及对应自变量取值。 若需要寻找最大值，对目标函数乘-1，并对最终结果乘-1即可
        :param func: type:函数。 des: 目标函数
        :param n: type:int。 des:原子个数
        :param var: type:int。 des: 自变量种类数，即目标函数是几元函数
        :param iter: type:int。 des:迭代次数
        :param T0: type:float。 des:初始温度
        :param Tf: type:float。 des:终止温度
        :param alpha: type:float。 des:温度衰减系数
        :param lb:type:list。 des: 每种自变量的下界
        :param ub:type:list。 des: 每种自变量的上界
        """
        if lb is None:
            lb = []
        if ub is None:
            ub = []

        # 初始化目标函数
        self.func = func
        # 初始化并行执行粒子个数
        self.n = n
        # 初始化自变量个数。即几元函数
        self.var = var
        # 初始化迭代次数
        self.iter = iter
        # 初始化温度衰减系数
        self.alpha = alpha
        # 初始化初始温度
        self.T0 = T0
        # 初始化结束温度
        self.Tf = Tf
        # 初始化边界
        self.lb = lb
        self.ub = ub

        # 为并行的粒子申请空间
        self.pos = [[0] * var for _ in range(n)]
        ## 初始化粒子位置
        for i in range(n):
            for j in range(var):
                self.pos[i][j] = lb[j] + random.random() * (ub[j] - lb[j])

        # 初始化最优解位置
        self.x = self.pos[0]
        for i in range(1, n):
            # 函数默认寻找最小值位置
            if self.func(self.pos[i]) < self.func(self.x):
                self.x = deepcopy(self.pos[i])


    def Metropolis(self, new_val, val):
        if new_val < val:
            return True
        else:
            # Metropolis准则
            p = math.exp((val - new_val)/self.T0)
            if random.random() < p:
                return True
            else:
                return False


    def run(self):
        # 迭代iter次
        for k in range(self.iter):
            # 温度从T0 -> Tf
            while self.T0 > self.Tf:
                # 每次迭代n个粒子
                for i in range(self.n):
                    # 为新解申请空间
                    new_pos = [0] * self.var
                    ## 新解是在原解的基础上变动来的
                    for j in range(self.var):
                        # 随机数-0.5保证正负分布均匀； *self.T0保证温度越高越活跃
                        new_pos[j] = self.pos[i][j] + (random.random() - 0.5) * self.T0
                        # 越界检查
                        if new_pos[j] > self.ub[j]:
                            new_pos[j] = self.ub[j]
                        if new_pos[j] < self.lb[j]:
                            new_pos[j] = self.lb[j]

                    new_val = self.func(new_pos)
                    val = self.func(self.pos[i])

                    # 更新全局最优解
                    if new_val < self.func(self.x):
                        self.x = deepcopy(new_pos)

                    # 判断是否更新解
                    if self.Metropolis(new_val, val):
                        self.pos[i] = deepcopy(new_pos)

                self.T0 = self.T0 * self.alpha

        print(f'函数最小值为{self.func(self.x)}')
        print(f'最小值对应的一组x取值为{self.x}')


"""一元函数测试"""
# sa = PSO(func, 50, 1, 0.9, 100, [-3], [3])


"""二元函数测试"""
sa = SA(func, 40, 2, 100, 100, 0.01, 0.9, [-15, -15], [15, 15])

sa.run()

二元函数测试结果：

函数最小值为8.000000028277022
最小值对应的一组x取值为[7.999866574862623, 6.000055456680725]