数学建模：模拟退火算法（SA）
Posted A-L-Kun
tags:
篇首语：本文由小常识网(cha138.com)小编为大家整理，主要介绍了数学建模：模拟退火算法（SA）相关的知识，希望对你有一定的参考价值。
模拟退火算法（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步。
其中有几个需要注意的点：
- 初始点的选取对算法结果有一定的影响，最好是多次运行对结果进行综合判断。
- 在算法运行初期，温度下降快，避免接受过多的差结果。当运行时间增加，温度下降减缓，以便于更快稳定结果。
- 当迭代次数增加到一定次数时，结果可能已经达到稳定，但是距离算法结束还有一段时间。在设计程序时应该加入适当的输出条件，满足输出条件即可结束程序。
可以大概分为这四个步骤：
- 第一步是由一个产生函数从当前解产生一个位于解空间的新解；为便于后续的计算和接受，减少算法耗时，通常选择由当前新解经过简单地变换即可产生新解的方法，如对构成新解的全部或部分元素进行置换、互换等，注意到产生新解的变换方法决定了当前新解的邻域结构，因而对冷却进度表的选取有一定的影响。
- 第二步是计算与新解所对应的目标函数差。因为目标函数差仅由变换部分产生，所以目标函数差的计算最好按增量计算。事实表明，对大多数应用而言，这是计算目标函数差的最快方法。
- 第三步是判断新解是否被接受,判断的依据是一个接受准则，最常用的接受准则是 Metropolis 准则: 若 ΔT < 0 则接受 S′ 作为新的当前解 S，否则以概率 P 接受 S′ 作为新的当前解 S。
- 第四步是当新解被确定接受时，用新解代替当前解，这只需将当前解中对应于产生新解时的变换部分予以实现，同时修正目标函数值即可。此时，当前解实现了一次迭代。可在此基础上开始下一轮试验。而当新解被判定为舍弃时，则在原当前解的基础上继续下一轮试验。
二、 实例分析
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)
最后的运行结果为：
本文来自博客园，作者：Steve_Anthony，转载请注明原文链接：https://www.cnblogs.com/liuzhongkun/p/17286030.html
模拟退火算法（SA）简介及Python实现
一、概述
模拟退火算法（Simulated Annealing，SA）是一种模拟物理退火过程而设计的优化算法。它的基本思想最早在1953年就被Metropolis提出，但直到1983年，Kirkpatrick等人才设计出真正意义上的模拟退火算法并进行应用。
模拟退火算法采用类似于物理退火的过程。先在一个高温状态下，然后逐渐退火，在每个温度下慢慢冷却，最终达到物理基态（相当于算法找到最优解）。模拟退火算法源于对固体退火过程的模拟，采用Metropolis准则，并用一组称为冷却进度表的参数控制算法的进程，使得算法在多项式时间里可以给出一个近似最优解。
模拟退火算法在某一初温下，伴随温数的不断下降，结合概率突跳特性在解空间中随机寻找目标函数的全局最优解。即在局部最优解的空间内能概率性地跳出并最终趋于全局最优。
二、模拟退火算法
一、算法的主要流程：
step1: 设定当前解（即为当前的最优解)：令
T
T
T =
T
0
T_0
T0 ，即开始退火的初始温度，随机生成一个初始解
x
0
x_0
x0 ，并计算相应的目标函数值E(
x
0
x_0
x0)。
step2: 产生新解与当前解差值：根据当前解
x
i
x_i
xi进行扰动，产生一个新解
x
j
x_j
xj，计算相应的目标函数值E(
x
j
x_j
xj),得到𝜟𝑬=𝑬(
x
j
x_j
xj)−𝑬(
x
i
x_i
xi)。
step3: 判断新解是否被接受 ：若𝜟𝑬<𝟎，则新解
x
j
x_j
xj被接受；若𝜟𝑬>𝟎 ，则新解
x
j
x_j
xj按概率e
−
(
𝑬
(
x
j
)
−
𝑬
(
x
i
)
)
T
i
-(𝑬(x_j)−𝑬(x_i)) \\over T_i
Ti−(E(xj)−E(xi))接受，
T
i
T_i
Ti为当前温度。
step4: 当新解被确定接受时：新解
x
j
x_j
xj被作为当前解
step5: 循环以上四个步骤：在温度
T
i
T_i
Ti下，重复k次的扰动和接受过程，接着执行下一步骤。
step6: 最后找到全局最优解：判断T是否已经达到终止温度
T
f
T_f
Tf，是，则终止算法；否，则转到循环步骤继续执行。
二、冷却进度表
退火过程由一组初始参数，即冷却进度表控制。它的目的是尽量使系统达到平衡，以使算法在有限的时间内逼近最优解。它包括：
1.控制温度参数的初值
T
0
T_0
T0
一般来说，只有足够大的 才能满足算法要求，但是由于不同问题处理的规模不同，所以这个“足够大”也是不同的，有的可能
T
0
T_0
T0 =100就好了，但是有的可能
T
0
T_0
T0 =1000还不够。实验表明：初温越大，获得全局最优解的机率越大，但花费的时间也会越长。
Metropolis准则：
在某个温度下固体分子从一个状态转移到另一个状态时，如果新状态的内能小，则无条件接受；如果新状态的内能大，则以一定的概率接受它。
2.马尔科夫链的长度
L
k
L_k
Lk：任意温度T的迭代次数
算法在马尔科夫链长度内持续进行”产生新解—判断—接受/舍弃”的迭代过程，对应着固体在某一恒定温度下趋于热平衡的过程。若在一定的温度下做无限次迭代，相应的马尔科夫链可以达到平稳分布概率。
马尔科夫链的选取还与温度控制参数
T
k
T_k
Tk的下降密切相关，缓慢下降可以避免过长的马尔科夫链。在控制参数的衰减函数已经选定的前提下，让每个取值都能够达到准平衡状态。
根据这一原则一般取
L
k
L_k
Lk=100N，其中N为问题的规模。
3. 控制参数T的终值
T
f
T_f
Tf（停止准则）：
合理的停止准则既能保证算法收敛于某一近似解,又能使最终解具有一定的全局性。通常可以根据迭代的次数或者终止的温度或者迭代过程在若干个相继的链中的解没有任何变化等条件来判断迭代的终止。
最终温度通常是0，但会耗掉许多模拟时间。温度趋近于0，其周围状态几乎是一样的。所以找寻一个低到可接受的温度即可。
4. 控制温度T的衰减函数(温度的更新)
不同退火方法的温度下降速度是不一样的,其中指数降温是最常用的一种退火策略,其温度变化很有规律,直接与参数相关,衰减函数：t
k
_k
k+
1
_1
1=α
t
k
t_k
tk，k=0，1，2… 其中α是一个接近1的常数。一般取0.5~0.99。该衰减函数对控制参数的衰减量是随算法进程递减的。
注意：小的衰减量可能导致算法进程迭代次数的增加,从而使算法进程接受更多的变换,从而访问更多的邻域,搜索更大范围的解空间,返回更高质量的最终解。
三、算法的优缺点
模拟退火算法的应用很广泛，可以高效地求解NP完全问题，如货郎担问题(Travelling Salesman Problem，简记为TSP)、最大截问题(Max Cut Problem)、0-1背包问题(Zero One Knapsack Problem)、图着色问题(Graph Colouring Problem)等等，但其参数难以控制，不能保证一次就收敛到最优值，一般需要多次尝试才能获得（大部分情况下还是会陷入局部最优值）。观察模拟退火算法的过程，具有以下主要优势：
- 迭代搜索效率高，并且可以并行化；
- 算法中有一定概率接受比当前解较差的解，因此一定程度上可以跳出局部最优；
- 算法求得的解与初始解状态S无关，因此有一定的鲁棒性；
- 具有渐近收敛性，已在理论上被证明是一种以概率l 收敛于全局最优解的全局优化算法。
三、python实现
代码我参考了别人的代码，如下
import numpy as np
# =================初始化参数===============
D = 10 # 变量维数
Xs = 20 # 上限
Xx = -20 # 下限
# ====冷却表参数====
L = 200 # 马可夫链长度 #在温度为t情况下的迭代次数
K = 0.95 # 衰减参数
S = 0.01 # 步长因子
T = 100 # 初始温度
YZ = 1e-7 # 容差
P = 0 # Metropolis过程中总接受点
# ====随机选点初值设定====
PreX = np.random.uniform(size=(D, 1)) * (Xs - Xx) + Xx
PreBestX = PreX # t-1代的全局最优X
PreX = np.random.uniform(size=(D, 1)) * (Xs - Xx) + Xx
BestX = PreX # t时刻的全局最优X
# ==============目标函数=============
def func1(x):
return np.sum([i ** 2 for i in x])
# ====每迭代一次退火一次(降温), 直到满足迭代条件为止===
deta = np.abs(func1(BestX) - func1(PreBestX)) # 前后能量差
trace = [] # 记录
while (deta > YZ) and (T > 0.1): # 如果能量差大于允许能量差 或者温度大于阈值
T = K * T # 降温
# ===在当前温度T下迭代次数====
for i in range(L): #
# ====在此点附近随机选下一点=====
NextX = PreX + S * (np.random.uniform(size=(D, 1)) * (Xs - Xx) + Xx)
# ===边界条件处理
for ii in range(D): # 遍历每一个维度
while NextX[ii] > Xs or NextX[ii] < Xx:
NextX[ii] = PreX[ii] + S * (np.random.random() * (Xs - Xx) + Xx)
# ===是否全局最优解 ===
if (func1(BestX) > func1(NextX)):
# 保留上一个最优解
PreBestX = BestX
# 此为新的最优解
BestX = NextX
# ====Metropolis过程====
if (func1(PreX) - func1(NextX) > 0): # 后一个比前一个好
# 接受新解
PreX = NextX
P = P + 1
else:
changer = -1 * (func1(NextX) - func1(PreX)) / T
p1 = np.exp(changer)
# 接受较差的解
if p1 > np.random.random():
PreX = NextX
P = P + 1
trace.append(func1(BestX))
deta = np.abs(func1(BestX) - func1(PreBestX)) # 修改前后能量差
print('最小值点\\n', BestX)
print('最小值\\n', func1(BestX))
以上是关于数学建模：模拟退火算法（SA）的主要内容，如果未能解决你的问题，请参考以下文章