本文借鉴1:模拟退火算法详细讲解(含实例python代码)
本文借鉴2:清风建模模拟退火讲解
一、退火的物理过程:
- 加温过程。其目的是增强粒子的热运动,使其偏离平衡位置。当温度足够高时,固体将熔为液体,从而消除系统原先存在的非均匀状态。
- 等温过程。对于与周围环境交换热量而温度不变的封闭系统,系统状态的自发变化总是朝自由能减少的方向进行的,当自由能达到最小时,系统达到平衡状态。
- 冷却过程。使粒子热运动减弱,系统能量下降,得到晶体结构。
总结:
先在一个高温状态下(设定初值),然后逐渐退火(控制温度T下降),在每个温度下慢慢冷却(迭代+Metropolis算法),最终达到物理基态(相当于算法找到最优解)。
其中Metropolis 准则是SA算法收敛于全局最优解的关键所在——Metropolis准则以一定的概率接受恶化解,这样就使算法跳离局部最优的陷阱。
二、退火算法的步骤:
求解函数最小值问题:
上面包含两个部分:Metropolis算法和退火过程,,分别对应内循环和外循环。
三.退火算法的过程:
在一切开始之前先在区间里随机取一个初值x(0)。
下面进入内循环。
1.内循环
内循环就是Metropolis算法:上面步骤的(3)至(5)步。即在每次温度下,迭代L次,寻找在该温度下能量的最小值(即最优解)。
- 扰动产生新解
迭代过程的新解x(n+1)是由x(n)随机扰动产生(形象理解就是在x(n)基础上向左或向右随机移动),并且在温度越高的时候,移动的幅度越大,这模拟了物理中粒子的无序运动。
-
Metropolis准则
(1)在固定温度下,整个迭代过程中温度不发生变化,能量发生变化,当前一个状态x(n)的能量大于后一个状态x(n+1)的能量时,说明状态x(n)的解没有状态x(n+1)的解好,所以接受状态x(n+1)。(这里说的能量即为实际问题中的目标函数值)(2)但是如果下一状态的能量比前一个状态的能量高时,该不该接受下一状态呢?在这里设置一个接受概率P,即如果下一状态的能量比前一个状态的能量高(即比前一个解差),则接受下一状态的概率为P,这就是Metropolis准则。
具体实现的时候,会随机生成一个在区间[0,1]上服从均匀分布的随机数r, 如果r<p则接受新解
- 关于概率P
(1)温度一定时, ∆E越小,概率越大,即目标函数相差越小接受的可能性越大。
(2)∆E一定时,温度越高,概率越大,即搜索前期温度较高时更有可能接受新解。
当迭代次数足够时,进入外循环
2.外循环
外循环就是退火过程,将固体达到较高的温度(初始温度T(0)),控制温度不断下降,当达到终止温度Tf时,冷却结束,即退火过程结束
- 冷却进度表
(1)初始温度最常取100。
(2)衰减函数最常用的是:
PYTHON代码实现
整个代码是我在浏览其他博友时看到的,并非我本人原创,只是我比较困惑内循环迭代为何是先随机生成L个点,每个迭代一次最终取最好的,而不是选取一个点迭代L次。
所以我按照后者改了一下(如果有明白为何按前者来做的朋友麻烦给我解答一下疑惑)
最后还画了能量随温度下降曲线和整个过程的一个动画图。
import math
from random import random
import matplotlib.pyplot as plt
import numpy as np
import matplotlib.animation as animation
# 目标函数
def func(x):
"""
函数优化问题: min f(x)=x^4-4*x^2-2,范围[-5,5]
"""
res = x ** 4 - 14 * x ** 2 + +10 * x + 4
return res
# SA类
class SA:
"""模拟退火算法的类"""
def __init__(self, func, iter=100, T0=100, Tf=0.01, alpha=0.95):
self.func = func
self.iter = iter # 内循环迭代次数
self.alpha = alpha # 降温系数
self.T0 = T0 # 初始温度
self.Tf = Tf # 温度终值
self.T = T0 # 当前温度
self.x = random() * 4 - 2 # 随机生成(-2,2)之间的一个值
self.most_best = []
self.history = {'x': [], 'f': [], 'T': [],'f1':[],'T1':[]} # 画图
def generate_new(self, x): # 扰动产生新解的过程
while True:
x_new = x + self.T * (random() - random())
if -5 <= x_new <= 5:
break # 重复得到新解,直到产生的新解满足约束条件
return x_new
# Metropolis准则
def Metrospolis(self, f, f_new):
if f_new <= f:
return 1 # 1表示接受,0表示拒绝
else:
p = math.exp((f - f_new) / self.T)
if random() < p: # random()生成一个[0,1)之间的随机数
return 1
else:
return 0
# 外循环迭代
def run(self):
while self.T > self.Tf: # 当温度大于终止温度的阈值
# 内循环迭代L次
for i in range(self.iter):
f = self.func(self.x) # 旧值
x_new = self.generate_new(self.x) # 产生新解
f_new = self.func(x_new) # 新值
if self.Metrospolis(f, f_new): # 判断是否接受新值
self.x = x_new # 如果接受新解,则把新解存入x
ft = self.func(self.x) # 迭代L次后最优值
xt = self.x
Tt = self.T
self.history['x'].append(xt)
self.history['f'].append(ft)
self.history['T'].append(Tt)
self.history['f1'].append(ft)
self.history['T1'].append(Tt)
# 温度按照一定的比例下降(冷却)
self.T = self.T * self.alpha
# 得到最优解
print(f"F={ft}, x={xt}")
# 运行程序
sa = SA(func)
sa.run()
# 能量变化
plt.plot(sa.history['T1'], sa.history['f1'])
plt.title('SA')
plt.xlabel('T')
plt.ylabel('f')
plt.gca().invert_xaxis()
plt.show()
# # 绘制随机点移动图像
def update_points(num):
"""
更新数据点
"""
point_ani.set_data(X[num], Y[num])
text_pt.set_text("T=%.3f" % T[num])
return point_ani, text_pt,
x = np.linspace(-5, 5, 100) # 生成[0,10]范围内,100个等距离的数组
y = func(x)
X = sa.history['x']
Y = sa.history['f']
T = sa.history['T']
fig = plt.figure(tight_layout=True)
plt.plot(x, y)
plt.xlabel('x')
plt.ylabel('f(x)')
point_ani, = plt.plot(X[0], Y[0], "ro")
plt.grid(ls="--")
text_pt = plt.text(4, 0.8, '', fontsize=10)
ani = animation.FuncAnimation(fig, update_points, np.arange(0, len(X)), interval=0.01, blit=True, repeat=False)
plt.show()