经典模拟退火算法介绍（含实例）

遇见繁星之夜~

已于 2024-05-31 10:48:55 修改

阅读量2.3k

点赞数 26

文章标签：模拟退火算法算法机器学习启发式算法 python

于 2024-05-28 15:29:35 首次发布

本文链接：https://blog.csdn.net/ygxdxm/article/details/139267899

版权

经典模拟退火算法

简述

模拟退火算法（simulated annealing，SA），最早的思想是由N．Metropolis等人于1953年提出的。

模拟退火算法的思想来源于对固体退火降温过程的模拟。即将固体加温至充分高,再让其徐徐冷却。在加热固体时，固体中原子的热运动不断增强,内能增大,随着温度的不断升高,固体的长程有序被彻底破坏,固体内部粒子随温度的升高而变为无序状态；冷却时，粒子逐渐趋于有序，在每个温度下都达到平衡状态；最后在常温下达到基态，同时内能也减为最小。根据 Metropolis准则，粒子在温度T时趋于平衡的几率为 $e^{-\frac{ΔE}{kT} }$ ，其中 $E$ 为温度 $T$ 时的内能， $Δ E$ 为其改变量， $k$ 为Boltzmann常数。用固体退火模拟组合优化问题，将内能 $E$ 模拟为目标函数值 $f$ ，温度T演化成控制参数 $t$ 。

模拟退火算法用于组合优化问题的出发点是基于物理中固体物质的退火过程与一般优化问题的相似性。

物理退火	优化问题
物质状态	解
能量最低的物质状态	最优解
退火过程	求解过程
温度	控制参数
能量	目标函数
等温过程	Metropolis抽样过程

算法流程

在这里插入图片描述

SA算法的基本思想是从选定的初始解开始，借助于控制参数t递减时产生的一系列Markov链中，利用一个新解产生装置和接受准则，重复进行“产生新解 —计算目标函数差—判断是否接受新解—接受或舍弃新解”，不断对当前解迭代，从而使目标函数得到最优结果的执行过程。由于固体退火必须缓慢降温，才能使固体在每一温度下都达到热平衡，最终趋于平衡状态。因此，控制参数 $t$ 的值必须缓慢衰减,才能确保模拟退火算法最终趋于优化问题的整体最优解。其求解步骤如下:

1°从可行解空间中任选一初始状态 $m_0$ 计算其目标函数值 $E(m_0)$ ,并选择初始控制温度 $T_0$ 和马尔可夫链的长度 $L$ ；
2°在可行解空间中产生一个随机扰动,用状态产生函数产生一个新状态 $m_1$ ,计算其目标函数值 $E(m_1)$ ；
3°根据接受准则判断是否接受：如果 $E(m_1)<E(m_0)$ ，则接受新状态 $m_1$ 为当前状态,否则按 Metropolis准则判决是否接受 $m_1$ ,若接受,则令当前状态等于 $m_1$ ;若不接受,则令当前状态等于 $m_0$ ；
4°根据某个收敛准则，判断抽样过程是否终止，是则转5°，否则转2°；（内循环终止准则,马尔可夫链长度判断）
5°按照某个温度冷却方案降低控制温度 $T$
6°根据某个收敛准则，判断退火过程是否终止,是则转7°,否则转2°；（外循环终止准则，是否收敛）
7°当前解作为最优解输出。

《模拟退火算法改进综述及参数探究》论文提到的传统模拟退火算法的步骤，描述得更清楚：

设目标函数为 $y = f (x)$ ，传统模拟退火算法流程如下：
Step1：选定初始控制温度 $T_0$ ，马氏链长度 $L_0$ ，在可行解空间中随机选取一个初始解 $i_0$ ，此时，最优解 $i=i_0$ ，迭代次数 $k = 0$ ，降温函数（即控制参数衰减函数） $T_k=h(k)$ ；
Step2：产生一次随机扰动，在可行解空间中得到一个新解 $ｊ$ ；

Step3：判断是否接受新解，判断准则为Metropolis准则：

(i)若 $f (i) >= f (j)$ ，则接受新解 $j$ ，此时最优解 $i = j$ ，
(ii)若 $f (i) < f (j)$ ，则依概率接受新解 $j$ ，即 $e^{\frac{f(i)-f(j)}{T_k} } >random(0,1]$ 时，接受新解 $j$ ，此时最优解 $i = j$ ，否则，拒绝 $j$ ，此时最优解仍为 $i$ ；
Step4：重复执行 $L_0$ 次Step2和Step3，得到链长为 $L_0$ 的马氏过程下的一个最优解；

Step5：判断是否满足停止准则，若满足则输出最优解，算法停止，否则执行Step6；

Step6：迭代次数ｋ＝ｋ＋１，最优解更新为Step4得到的解，温度函数变成 $T_{k+1}$ ，马氏链长度变为 $L_{k+1}$ ，回到Step2

冷却进度表

冷却进度表规定了控制参数 $T$ 的初值 $T_s$ 、控制参数 $T$ 的衰减函数、控制参数 $T$ 的终值 $T_e$ （停止准则）和Markov链的长度 $L_k$ .

降温函数

常见的降温函数有：
$T_k=\frac{T_s}{k} \tag{1}$

$T_k=T_s\alpha^{k} \tag{2}$

其中 $\alpha$ 为温度衰减率，一般取 $0.7\le \alpha \le 1.0$ .
$T_k=T_se^{-ck} \tag{3}$
其中c为常数。

初温、末温

冷却进度表中的另一关键参数就是初末温的设置．根据模拟退火算法解的收敛性理论可知，当初温足够高，降温过程足够慢时，算法才可以收敛到全局最优解．温差越大，解越精确，但相应的运行时间也会增加。

马氏链长度

在模拟退火算法中，马尔可夫链的长度指的是在某个特定温度下进行的迭代次数或状态转移的次数。在每一温度下，可以用一个马氏链来刻画状态空间中的状态转移，其中的转移概率就是每一个当前解会转移到新解的概率．在每一温度下，重复多次产生新解并去判断是否转移，而产生新解的次数就是对应的马氏链的长度。

迭代次数：内层循环（即在同一温度下的迭代）的迭代次数构成了马尔可夫链的长度。在每一次迭代中，算法会尝试对当前解进行微小变动以探索可能更优的解。
扰动次数：马尔可夫链的长度决定了在每个温度下系统状态转换（扰动）的次数。这个长度可以影响算法的搜索能力和最终结果的质量。太短的链可能导致搜索不充分，而太长则可能导致计算量过大。

接受准则(Metropolis准则)

模拟退火的概念是根据熔融金属中粒子的统计力学与复杂的组合最优化问题的求解过程的相似性而提出来的．统计力学的研究表明，在温度 $T$ 、粒子停留在状态 $r$ 满足波兹曼（Boltzmann）概率分布.
$P(E(r))=\frac{1}{Z(T)} e^{\frac{-E(r)}{k_{b}T } }$
式中： $E (r)$ 为 $r$ 状态的能量； $k_b＞0$ 为Boltzmann 常数； $T$ 为绝对温度； $Z (T)$ 为概率分布的标准化因子.

在温度 $T$ 时,根据金属内部粒子的状态选择规律,Metropolis等提出了以下准则：设初始状态 $i$ 为当前状态,该状态的能量为 $E_i$ ,然后对该状态作一随机微扰得到1个新状态 $j$ ,新状态的能量为 $E_j$ ．如果 $E_j＜ E_i$ , $j$ 就为重要状态；如果 $E_j> E_i$ 考虑到热运动的影响, $j$ 是否为重要状态要根据固体处该状态的概率来判断．固体处于状态 $i$ 和 $j$ 的概率的比值等于相应Boltzmann因子的比值即
$e^{\frac{E_i-E_j}{k_{b}T } } ,0<r<1$
用随机数据发生器产生一个[0,1]区间的随机数 $ξ$ ,若 $r ＞ ξ$ ,则新状态 $j$ 为重要状态,就以 $j$ 取代 $i$ 成为当前状态,否则仍以 $i$ 为当前状态,再重复以上新状态的产生过程．在大量固体状态变换后,系统趋于能量较低平衡状态,固体状态的概率分布趋于 Boltzmann概率分布。从r的表达式来看，T越小(T>0),r的值越小，接受状态 $j$ 的概率就越低。

经典模拟退火算法求解TSP问题

问题描述

在前面的文章遗传算法求解TSP问题介绍过，利用遗传算法来求解TSP问题。

TSP（traveling salesman problem）问题是典型的NP完全问题，可描述为：已知n个城市相互之间的距离，某一旅行商从某个城市出发访问每个城市一次且仅一次，最后回到出发城市，如何安排才能使所走路线最短。该问题可描述为搜索自然子集 $X=\left \{ 1,2,3,...,n \right \}$ (其中每个元素表示城市的编号)的一个排列 $\rho (X)=\left \{ V_{1},V_{2},...,V_{n} \right \}$ ，使得
$T_{d} =\sum_{i=1}^{n-1}d(V_{i},V_{i+1})+ d(V_{n},V_{1})$
取最小值，其中 $d(V_{i},V_{i+1})$ 表示城市 $V_{i}$ 到城市 $V_{i+1}$ 的距离。现提供14个城市的位置坐标如下。

城市编号	X坐标	Y坐标	城市编号	X坐标	Y坐标
1	16.47	96.10	8	17.20	96.29
2	16.47	94.44	9	16.30	97.38
3	20.09	92.54	10	14.05	98.12
4	22.39	93.37	11	16.53	97.38
5	25.23	97.24	12	21.52	95.59
6	22.00	96.05	13	19.41	97.13
7	20.47	97.02	14	20.09	92.55

模拟退火算法流程

冷却进度表参数设置

降温函数

降温函数选择幂函数，设置温度衰减率，模拟退火降温过程
$T_k=T_s\alpha^{k}$
这里 $\alpha$ 取0.9.
初温、末温

初温设置为1000，末温设置为0.001.
马氏链长度

这里选择固定的马氏链长度，均取为500，即每个温度下马氏链长度均为500.

解的定义

类似于之前遗传算法中解的定义。假设有5个城市，编号分别为1，2，3，4，5，那么 $S_{1}$ =|2|5|1|3|4|就是一个可行解，可以利用随机数发生器生成初始解，当前问题中有14个城市，以此类推即可。程序中将城市的完整的访问顺序定义为一条路径，每条路径拥有总长度。

产生新解的方法

依然以5个城市进行举例，假设城市编号为1，2，3，4，5，若当前可接受解为 $S_{1}$ =|2|5|1|3|4|，那么在它基础之上产生新解，利用随机数发生器生成两个不同的[1,5]之间的随机整数，假设r1=1，r2=3，直接交换两个位置的元素即可，新解即为 $S_{2}$ =|1|5|2|3|4|。

Metropolis准则

假设评价路径总距离的函数为 $f (S)$ ,记原来的可接受解的路径长度为 $f(S_1)$ ,新解的路径长度为 $f(S_2)$ ，设 $\Delta f(s)=f(S_2)-f(s_1)$ ,那么定义接受准则如下：
$P=\left\{\begin{matrix} 1,\Delta f(s)\le 0\\ e^{\frac{-\Delta f(s)}{k_{b}T } },\Delta f(s)>0 \end{matrix}\right.$

当 $\Delta f(s)\le 0$ 时，完全接受新解，即令 $S_1 = S_2$ ;当 $\Delta f(s)\gt 0$ 时，以概率 $e^{\frac{-\Delta f(s)}{k_{b}T } }$ 的情况接受新解。

python实现算法流程

好久没写过python代码了，将就着用吧🤣

import random
import math
import copy


class Individual:
    """ 定义解的形式，包含路径和路径距离 """
    def __init__(self, m_route, m_obj_val):
        self.route = m_route
        self.obj_val = m_obj_val


class CityPos:
    """ 城市坐标定义 """
    def __init__(self, pos_x, pos_y):
        self.x = pos_x
        self.y = pos_y


class Sa_Tsp:
    """ 使用模拟退火算法求解tsp问题的类 """

    def __init__(self, t_start, t_end, t_factor, markov_l):
        self.best_solution = None
        self.accept_solution = None
        self.city_num = None
        self.distance_info = None
        self.t_start = t_start
        self.t_end = t_end
        self.t_factor = t_factor
        self.markov_l = markov_l

    def open_route_date(self, file_path):
        """ 从文件中加载各个城市的位置信息并计算城市之间的距离 """
        pos_list = []
        # 读取文件数据
        with open(file_path, "r", encoding="utf-8") as file:
            # 逐行读取数据
            for line in file:
                x, y = line.strip().split()
                pos_list.append(CityPos(float(x), float(y)))

        self.city_num = len(pos_list)

        self.distance_info = [[0.0 for j in range(self.city_num)] for i in range(self.city_num)]

        for i in range(self.city_num):
            for j in range(self.city_num):
                if i < j:
                    self.distance_info[i][j] = self.cal_distance(pos_list[i], pos_list[j])
                elif i > j:
                    self.distance_info[i][j] = self.distance_info[j][i]

    @staticmethod
    def cal_distance(pos1, pos2):
        """计算两点之间的距离"""
        res = 0.0
        res += pow(pos1.x - pos2.x, 2)
        res += pow(pos1.y - pos2.y, 2)
        return math.sqrt(res)

    def cal_obj_val(self, route):
        """ 计算一条路径的距离 """
        obj_val = 0.0
        for i in range(self.city_num - 1):
            city_idx_pre = route[i] - 1
            city_idx_next = route[i + 1] - 1
            obj_val += self.distance_info[city_idx_pre][city_idx_next]

        city_idx_start = route[0] - 1
        city_idx_end = route[self.city_num - 1] - 1
        obj_val += self.distance_info[city_idx_start][city_idx_end]
        return obj_val

    def product_new_solution(self):
        """ 产生新解 """
        cur_route = copy.deepcopy(self.accept_solution.route)

        # 随机产生两个位置
        pos1 = random.randint(0, self.city_num - 1)
        pos2 = random.randint(0, self.city_num - 1)

        while pos1 == pos2:
            pos2 = random.randint(0, self.city_num - 1)

        pos1_val = cur_route[pos1]
        cur_route[pos1] = cur_route[pos2]
        cur_route[pos2] = pos1_val

        return Individual(cur_route, self.cal_obj_val(cur_route))

    def metropolis(self, cur_solution):
        """ metropolis接受准则 """
        if self.accept_solution.obj_val > cur_solution.obj_val:
            return True
        else:
            sub_obj_val = cur_solution.obj_val - self.accept_solution.obj_val
            accept_pro = math.exp(-sub_obj_val / self.t_start)
            if accept_pro > random.random():
                return True
            return False

    def sa(self):
        """ 模拟退火算法主流程 """
        base_route = []

        for i in range(self.city_num):
            base_route.append(i + 1)
        random.shuffle(base_route)

        # 初始化可接受解
        self.accept_solution = Individual(base_route, self.cal_obj_val(base_route))
        self.best_solution = copy.deepcopy(self.accept_solution)

        while self.t_start > self.t_end:

            for i in range(self.markov_l):
                cur_solution = self.product_new_solution()

                # 根据接受准则决定当前解的去留
                if self.metropolis(cur_solution):
                    self.accept_solution = cur_solution

                # 直接比较当前可接受解和最优解的目标函数值
                if self.best_solution.obj_val > self.accept_solution.obj_val:
                    self.best_solution = copy.deepcopy(self.accept_solution)

            # 降温
            self.t_start *= self.t_factor

        self.out_put_best_solution()

    def out_put_best_solution(self):
        print("best solution and obj_val:")
        best_route = self.best_solution.route
        for i in best_route:
            print(i, end='->')
        print(best_route[0])
        print(self.best_solution.obj_val)


if __name__ == '__main__':
    my_sa_tsp = Sa_Tsp(1000, 1e-3, 0.9, 500)
    # 自行编写城市数据
    my_file_path = "C://xxx//xxx//SA_tsp//city_pos_data1.txt"
    my_sa_tsp.open_route_date(my_file_path)
    my_sa_tsp.sa()

运行结果

python程序运行结果：

best solution and obj_val:
10->9->11->8->13->7->12->6->5->4->3->14->2->1->10
29.340520066994227

matlab仿真结果：

从运行结果上来看，模拟退火算法处理该问题得到的最短路径依然为29.3405，随着退火过程的不断进行，逐渐向最优解靠近，并最终找到最优解。

遇见繁星之夜~

关注

26
点赞
踩
28

收藏

觉得还不错? 一键收藏
0
评论
经典模拟退火算法介绍（含实例）

模拟退火算法的思想来源于对固体退火降温过程的模拟。即将固体加温至充分高,再让其徐徐冷却。在加热固体时，固体中原子的热运动不断增强,内能增大,随着温度的不断升高,固体的长程有序被彻底破坏,固体内部粒子随温度的升高而变为无序状态；冷却时，粒子逐渐趋于有序，在每个温度下都达到平衡状态；最后在常温下达到基态，同时内能也减为最小。
复制链接

扫一扫