浅析模拟退火算法及其MATLAB程序实现

浅析模拟退火算法及其MATLAB程序实现

前言

现代优化算法是 80 年代初兴起的启发式算法。这些算法包括禁忌搜索（tabusearch），模拟退火（simulated annealing），遗传算法（genetic algorithms），人工神经网络（neural networks）。它们主要用于解决大量的实际应用问题。目前，这些算法在理论和实际应用方面得到了较大的发展。无论这些算法是怎样产生的，它们有一个共同的目标－求 NP-hard 组合优化问题的全局最优解。虽然有这些目标，但 NP-hard 理论限制它们只能以启发式的算法去求解实际问题。

本文主要讨论的是$模拟退火（simulatedannealing）$。
算法参见：《数学建模算法与应用》

一、模拟退火的基本原理

材料统计力学中的退火现象

模拟退火算法得益于材料的统计力学的研究成果。统计力学表明材料中粒子的不同结构对应于粒子的不同能量水平。在高温条件下，粒子的能量较高，可以自由运动和重新排列。在低温条件下，粒子能量较低。如果从高温开始，非常缓慢地降温（这个过程被称为退火），粒子就可以在每个温度下达到热平衡。当系统完全被冷却时，最终形成处于低能状态的晶体。如果用粒子的能量定义材料的状态，Metropolis 算法用一个简单的数学模型描述了退火过程。假设材料在状态i 之下的能量为 E(i) ，那么材料在温度T 时从状态i 进入状态 j 就遵循如下规律：
（1）如果 $E(j)\le E(i)，接受该状态被转换。（2）如果E(j)>E(i)$ ，则状态转换以如下概率被接受：$$e^{\frac{E(i)-E(j)}{KT}}$$

其中 K 是物理学中的波尔兹曼常数，T 是材料温度
概括来说，该物理现象的实质为：温度越高（即能量较高），系统越不稳定：温度越低（即能量较低），系统越稳定。
我们可以考虑如下极限问题：当温度降至很低时，材料会以很大概率进入最小能量状态（该结论有详细的理论推导，参见随机数学）。假定我们要解决的问题是一个寻找最小值的优化问题。将物理学中模拟退火的思想应用于优化问题就可以得到模拟退火寻优方法。

模拟退火

事实上，系统的能量常常形象为数学建模中一个$目标函数f(x)$，比如：成本、路程、收益$\dots \dots$
在对目标函数进行模拟退火之后，可以得到问题的优化解$x_{min}$。
通过对比材料统计力学中的模拟退火现象，我们可以得到模拟退火算法的如下规律：

（1）如果 $f(x_i)\le f(x_j)，接受该状态被转换。（2）如果f(x_j)>f(x_i)$ ，则状态转换以如下概率被接受：$$e^{\frac{f(x_i)-f(x_j)}{KT}}$$
这里将以最优路径为例详细介绍模拟退火法的具体步骤（有关图论问题，请参阅图论（一）基本概念：

$（1）$组装邻接矩阵D

$$D=\left(\begin{array}{ccccc}{d}_{11}& d_{12}& d_{13}& \cdots & d_{1n}\\ d_{21}& d_{22}& d_{23}& \cdots & d_{2n}\\ ⋮& ⋮& ⋮& \ddots & ⋮\\ d_{m1}& d_{m2}& d_{m3}& \cdots & d_{mn}\end{array}\right)$$
其中$d_{i,j}$表示从$2\Rightarrow 3$的距离。关于邻接矩阵，请参见邻接矩阵。

$（2）$解空间

解空间$S$可表示为 S = { ( π 1 , ⋯   , π n ) ∣ π 1 = 1 , ( π 2 , ⋯   , π n − 1 ) ， π n } S=\lbrace(\pi_1,\cdots,\pi_n)\lvert\pi_1=1,(\pi_2,\cdots,\pi_{n-1})，\pi_n\rbrace S={(π1​,⋯,πn​)∣π1​=1,(π2​,⋯,πn−1​)，πn​}其中$({\pi }_2,\cdots ,{\pi }_{n-1})$为节点$2,3,\cdots ,n-1$的随机组合，这里我们采用$({\pi }_1,\cdots ,{\pi }_n)$的初始路径顺序$${\pi }_1\cdots {\pi }_{u-1}{\pi }_u{\pi }_{u+1}\cdots {\pi }_{v-1}{\pi }_v{\pi }_{v+1}\cdots {\pi }_n\text{\hspace{1em}(1)}$$

$（3）$目标函数

$$f(x)=d_{min}(f({\pi }_1,{\pi }_2\cdots ,{\pi }_n))$$

$（4）$新解的产生

这里采用的是$2$变换法:任选序号$u,v$,交换$u$与$v$之间的顺序，变成$逆序$，此时新路径为：
$${\pi }_1\cdots {\pi }_{u-1}{\pi }_v{\pi }_{v-1}\cdots {\pi }_{u+1}{\pi }_u{\pi }_{v+1}\cdots {\pi }_n\text{\hspace{1em}(2)}$$

$（5）$代价函数差

$$\delta f=f(x_i)-f(x_j)=(d_{{\pi }_{u-1}{\pi }_v}+d_{{\pi }_u{\pi }_{v+1}})-(d_{{\pi }_{u-1}{\pi }_u}+d_{{\pi }_v{\pi }_{v+1}})\text{\hspace{1em}(3)}$$

$（6）$接受准则

$$P=\{\begin{array}{ll}1,& \text{δf<0}\\ e^{\frac{\delta f}{KT}},& \text{δf>0}\end{array}\text{\hspace{1em}(4)}$$

$（6）$降温

$$T=\alpha T\text{\hspace{1em}(5)}$$

$（7）$终止条件

理论上，$$\underset{T\to -\infty }{\lim }$$时，系统将有且仅有$一个最优解$，但实际上不可能有这种情况，于是将限定一个$e$,使$T<e$时跳出循环，得到一个相对最优解$x_{min}$。

二、MATLAB程序实现

1.组装邻接矩阵

代码如下（示例）：

clear
sj=load('Desktop\数据2.txt');
sj=[sj,sj(:,1)]';
n=size(sj,1);
% x=sj0(:,[1:2:8]);x=x(:);
% y=sj0(:,[2:2:8]);y=y(:);
% sj=[70,40;x,y;70,40];%[70,40]为起点和终点坐标
% sj=sj*pi/180;
path=[1:n,1];%路径初始化
long=0;%距离初始化
d=zeros(n);%邻接矩阵初始化
for i=1:n-1
    for j=i+1:n
%         d(i,j)=6370*acos(cos(sj(i,1)-sj(j,1))*cos(sj(i,2))*cos(sj(j,2))+sin(sj(i,2))*sin(sj(j,2)));
        d(i,j)=sqrt((sj(i,1)-sj(j,1))^2+(sj(i,2)-sj(j,2))^2);
    end
end
d=d+d';%获得邻接矩阵（对称阵）

2.退火与新解产生

代码如下（示例）：

T=1;e=0.1^300;dt=0.99;
while T>e
% c=sort(2+floor(100*rand(1,2)));
c=sort(2+floor((n-2)*rand(1,2)));
df=d(path(c(1)-1),path(c(2)))+d(path(c(1)),path(c(2)+1))-d(path(c(1)-1),path(c(1)))-d(path(c(2)),path(c(2)+1));
if df<0
    path=[path(1:c(1)-1),path(sort(c(1):c(2),'descend')),path(c(2)+1:n)];
elseif rand<=exp(-df/T)
     path=[path(1:c(1)-1),path(sort(c(1):c(2),'descend')),path(c(2)+1:n)];
end
T=T*dt;
end

for i=1:n-1
    long=long+d(path(i),path(i+1));
 end
 GreatX=sj(path,1);
 GreatY=sj(path,2);
 plot(GreatX,GreatY,'-*');
 path
 long