模拟退火算法原理
退火这个词，其实是铁匠发明的。它的意思很简单，就是将铁匠炉烧热后，再把下边的火撤掉，让金属在炉子里边慢慢冷却。人们发现，这个缓慢的降温过程能消除金属内部的各种缺陷，使得其恢复能量最低的状态。后来，受到退火工艺的启发，研究人员把将退火的缓慢的降温思路推广开来，用于寻找复杂函数的全局最优值。举个简单的例子，氦原子在金属中空位附近运动时，其势能函数大概长这样：
举个简单的例子，氦原子在金属中空位附近运动时，其势能函数大概长这样：
把势能函数看成一个轨道，再把氦看成一个小球放到这个轨道上。那么，小球的高度就是它的重力势能。求这个函数的全局最小值，其实就是让小球滚到最深那个坑过程。
由于函数中存在非常多局域极小值（浅坑），常见的贪心算法（如梯度下降法，就是闭着眼睛沿坡往下滚）会陷入最近的局域极小值，很错过到全局最小值（深坑）。
那么，怎么才能让滚进深坑呢？
我们先晃一晃整个体系，让小球随机动起来。我们把晃动的剧烈程度称为温度T。在足够长的时间内。小球出现在i点的概率呈玻尔兹曼分布，与i点的高度(势能Ei)和温度T密切相关：
温度为0时，小球出现在全局最小值的概率为100%。
那么，把温度设为0不就行了？
当然不行。上面的概率成立的前提是提供无穷长的时间，使得小球能够遍历整个函数。
实际中，无穷长的时间当然是不现实的。小球从低概率区向高概率区转移，这个过程是需要消耗时间的。
高温下，转移速度快。但另一方面，高温下，小球落在各个坑里的概率差别比较小。
为了让小球能够尽快的转移到最低点，我们可以采取一个降温策略：
- 先设定一个较高的温度T，使得小球有足够的能量越过能垒，小球到处乱跑的概率比较大。
- 随机移动小球的位置，检查其移动前后的能量差 ΔE。
- 若 ΔE<0 ，小球能量降低，则接受此次移动。
- 若 ΔE>0，小球能量升高，有一定概率接受此次移动，概率为：exp(-ΔE/(KBT)) 。这个判定方法称为Metropolis判据。
- 逐渐降低温度T，使得小球滑向低处的概率逐渐增大，直到小球缓缓收敛到最低点，不再晃动。
通过这样一个逐渐降温的过程，小球最终有很大概率（并不是一定）能找到全局最低值。这个算法就是模拟退火算法。
模拟退火算法最大的优势在于其通用性：函数形式如何复杂我不管，我就一步步瞎跑，然后按Metropolis判据概率性的接受这一步就行。
流程图见下：
算例
计算函数f(x)=∑x^2(-20≤x≤20)的最小值，其中个体x的维数n= 10。 这是一个简单的平方和函数，只有一个极小点x=(0, 0, …0)， 理论最小值f(0，0, .，. 0)=0
MATLAB版求解
%%%%%%%%%%%%%%%%%%%%%%模拟退火算法解决函数极值%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%初始化%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear all; %清除所有变量
close all; %清图
clc; %清屏
D=10; %变量维数
Xs=20; %上限
Xx=-20; %下限
%%%%%%%%%%%%%%%%%%%%%%%%%%%冷却表参数%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
L = 200; %马可夫链长度
K = 0.998; %衰减参数
S = 0.01; %步长因子
T=100; %初始温度
YZ = 1e-8; %容差
P = 0; %Metropolis过程中总接受点
%%%%%%%%%%%%%%%%%%%%%%%%%%随机选点 初值设定%%%%%%%%%%%%%%%%%%%%%%%%%
PreX = rand(D,1)*(Xs-Xx)+Xx;
PreBestX = PreX;
PreX = rand(D,1)*(Xs-Xx)+Xx;
BestX = PreX;
%%%%%%%%%%%每迭代一次退火一次(降温), 直到满足迭代条件为止%%%%%%%%%%%%
deta=abs( func1( BestX)-func1(PreBestX));
while (deta > YZ) && (T>0.001)
T=K*T;
%%%%%%%%%%%%%%%%%%%%%在当前温度T下迭代次数%%%%%%%%%%%%%%%%%%%%%%
for i=1:L
%%%%%%%%%%%%%%%%%在此点附近随机选下一点%%%%%%%%%%%%%%%%%%%%%
NextX = PreX + S* (rand(D,1) *(Xs-Xx)+Xx);
%%%%%%%%%%%%%%%%%边界条件处理%%%%%%%%%%%%%%%%%%%%%%%%%%
for ii=1:D
if NextX(ii)>Xs | NextX(ii)<Xx
NextX(ii)=PreX(ii) + S* (rand *(Xs-Xx)+Xx);
end
end
%%%%%%%%%%%%%%%%%%%%%%%是否全局最优解%%%%%%%%%%%%%%%%%%%%%%
if (func1(BestX) > func1(NextX))
%%%%%%%%%%%%%%%%%%保留上一个最优解%%%%%%%%%%%%%%%%%%%%%
PreBestX = BestX;
%%%%%%%%%%%%%%%%%%%此为新的最优解%%%%%%%%%%%%%%%%%%%%%%
BestX=NextX;
end
%%%%%%%%%%%%%%%%%%%%%%%% Metropolis过程%%%%%%%%%%%%%%%%%%%
if( func1(PreX) - func1(NextX) > 0 )
%%%%%%%%%%%%%%%%%%%%%%%接受新解%%%%%%%%%%%%%%%%%%%%%%%%
PreX=NextX;
P=P+1;
else
changer = -1*(func1(NextX)-func1(PreX))/ T ;
p1=exp(changer);
%%%%%%%%%%%%%%%%%%%%%%%%接受较差的解%%%%%%%%%%%%%%%%%%%%
if p1 > rand
PreX=NextX;
P=P+1;
end
end
trace(P+1)=func1( BestX);
end
deta=abs( func1( BestX)-func1 (PreBestX));
end
disp('最小值在点:');
BestX
disp( '最小值为:');
func1(BestX)
figure
plot(trace(2:end))
xlabel('迭代次数')
ylabel('目标函数值')
title('适应度进化曲线')
%%%%%%%%%%%%%%%%%%%%%%%%%适应度函数%%%%%%%%%%%%%%%%%%%%%%%
function result=func1(x)
summ=sum(x.^2);
result=summ;
end
从图中可以看到迭代次数太多，是个体寻优，需要花费很多时间。相比其他群体寻优方法，模拟退火算法很差。
python版求解
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
# @Author: yudengwu(余登武）
# @Date : 2021/6/21
#@email:1344732766@qq.com
import numpy as np
import matplotlib.pyplot as plt
import matplotlib as mpl
import matplotlib; matplotlib.use('TkAgg')
mpl.rcParams['font.sans-serif'] = ['SimHei'] # 指定默认字体
mpl.rcParams['axes.unicode_minus'] = False # 解决保存图像是负号'-'显示为方块的问题
#=================初始化参数============
D=10 #变量维数
Xs=20 #上限
Xx=-20 #下限
#==冷却表参数==
L = 200 #马可夫链长度 #在温度为t情况下的迭代次数
K = 0.998 #衰减参数
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#降温
print(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))
plt.plot(trace,label='迭代曲线')
plt.xlabel('迭代次数')
plt.ylabel('目标函数值')
plt.show()
复杂约束求解：python版求解
在约束求解中，添加了惩罚项。
import numpy as np
import matplotlib.pyplot as plt
import matplotlib as mpl
import matplotlib; matplotlib.use('TkAgg')
mpl.rcParams['font.sans-serif'] = ['SimHei'] # 指定默认字体
mpl.rcParams['axes.unicode_minus'] = False # 解决保存图像是负号'-'显示为方块的问题
#=================初始化参数============
D=2 #变量维数 ,对应x,y
Xs=1 #搜索区间上限
Xx=-1 #搜索区间下限
#==冷却表参数==
L = 300 #马可夫链长度 #在温度为t情况下的迭代次数
K = 0.98 #衰减参数
S = 0.01 #步长因子
T=100 #初始温度
YZ = 1e-8 #容差(运行t-1,，t时刻的误差）
P = 0 #Metropolis过程中总接受点
eloss=0.1 #允许的惩罚项误差
#====随机选点 初值设定====
PreX = np.zeros(shape=(D,1))
PreX[0]=np.random.uniform(1,2,1)#X范围【1,2】
PreX[1]=np.random.uniform(-1,0,1)#y范围【-1,0】
PreBestX = PreX#t-1代的全局最优X
PreX = np.zeros(shape=(D,1))
PreX[0]=np.random.uniform(1,2,1)#X范围【1,2】
PreX[1]=np.random.uniform(-1,0,1)#y范围【-1,0】
BestX = PreX#t时刻的全局最优X
#==============目标函数=============
def func1(X):
A = 10
pi = np.pi
x = X[0]
y = X[1]
return 2 * A + x ** 2 - A * np.cos(2 * pi * x) + y ** 2 - A * np.cos(2 * pi * y)
#==========惩罚项===========
def calc_e(X):
"""计算个体的目惩罚项，X 的维度是 size * 2 """
ee = 0
"""计算第一个约束的惩罚项"""
e1 = X[0] + X[1] - 6
ee += max(0, e1)
"""计算第二个约束的惩罚项"""
e2 = 3 * X[0] - 2 * X[1] - 5
ee += max(0, e2)
return ee
#====每迭代一次退火一次(降温), 直到满足迭代条件为止===
deta=np.abs(func1(BestX)-func1(PreBestX))#前后能量差
trace=[]#记录
while (deta > YZ) and (T>0.1):#如果能量差大于允许能量差 或者温度大于阈值
T = K * T#降温
print(T)
# =在当前温度T下迭代次数=
for i in range(L): #
# =在此点附近随机选下一点==
NextX = PreX + S * (np.random.uniform(size=(D, 1)) * (Xs - Xx) + Xx)
# ===边界条件处理
while NextX[0] > 2 or NextX[0] < 1:#x属于【1,2】
NextX[0] =PreX[0] + S* (np.random.uniform() * (Xs - Xx) + Xx)
while NextX[1] > 0 or NextX[0] < -1: # y属于【-1,0】
NextX[1] = PreX[1] + S* (np.random.uniform() * (Xs - Xx) + Xx)
# ===是否全局最优解 ===
if (func1(BestX) > func1(NextX)) and (calc_e(BestX)<eloss) and (calc_e(NextX)<eloss):#两者都没有违反约束，下一代适应度更好时选下一代
# 保留上一个最优解
PreBestX = BestX
# 此为新的最优解
BestX = NextX
elif (calc_e(BestX)>eloss) and (calc_e(NextX)>eloss) and (func1(BestX) > func1(NextX)):#两者都违反约束，下一代适应度更好时（选下一代）
# 保留上一个最优解
PreBestX = BestX
# 此为新的最优解
BestX = NextX
elif (calc_e(BestX)>eloss) and (calc_e(NextX)<eloss) :#上一代惩罚项不满足约束，下一代惩罚项满足约束（选下一代）
# 保留上一个最优解
PreBestX = BestX
# 此为新的最优解
BestX = NextX
# ====Metropolis过程 当代温度T 前后个体判断===
if (func1(PreX) - func1(NextX) > 0) and (calc_e(PreX)<eloss) and(calc_e(NextX)<eloss) : # 两个惩罚项都满足约束，后一个的适应度更小（选后一个）
# 接受新解
PreX = NextX
P = P + 1
elif (calc_e(PreX)>eloss) and(calc_e(NextX)>eloss) and (func1(PreX) - func1(NextX) > 0) :#两个惩罚项都不满足时，后一个适应度比前一个小（选后一个）
# 接受新解
PreX = NextX
P = P + 1
elif (calc_e(PreX)>eloss) and(calc_e(NextX)<eloss) :#前一个惩罚项不满足，后一个满足（选后一个）
# 接受新解
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('最优值X\n',BestX)
print('最优目标函数值\n',func1(BestX))
print('最优值对应的惩罚项',calc_e(BestX))
plt.plot(trace,label='迭代曲线')
plt.xlabel('迭代次数')
plt.ylabel('目标函数值')
plt.show()
作者：电气-余登武
1186
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
