选择ssm框架做网站的好处/企业网站建设的流程
模拟退火算法(Simulated annealing, SA)是一种基于蒙特卡罗(Monte Carlo)思想设计的，常用于在较大的解空间中搜索近似全局最优解的优化算法。本文将从模拟退火算法的历史、形式、特点、改进和实现等方面出发，对模拟退火算法作出介绍。
一、模拟退火算法的诞生
在众多科学领域，组合优化问题广泛存在，其中不乏 NP 完全问题(Nondeterministic polynomial-time complete problem)，典型的如做决定版的旅行商问题(Decision version of the traveling salesman problem, TSP)。目前为止，所有已知的解决 NP 完全问题的算法的时间复杂度皆大于多项式时间。为处理这类问题，人们提出多种折衷的算法，包括近似算法(Approximation algorithm)、随机化算法(Randomized algorithm)、启发式算法(Heuristic algorithm)等，以求在可接受的时间内得到可接受的解。模拟退火算法即是近似算法、随机化算法中的一例。
美国物理学家 N. Metropolis 等人在 1953 年发表的《Equation of State Calculations by Fast Computing Machines》一文中，使用蒙特卡罗模拟法计算多分子系统中分子的能量分布。1983 年，物理学家 S. Kirkpatrick、C. D. Gelatt 和 M. P. Vecchi 在《Science》上发表了《Optimization by Simulated Annealing》。文章摘要中写道：「在统计力学(在有限温度下的热平衡中具有多个自由度的系统的行为)和多变量或组合优化问题(寻找给定多参数函数的最小值)之间，存在深刻且有用的联系。对固体退火的细节类比提供了一个优化复杂庞大系统的框架。这种与统计力学的联系带来了新信息，提供了一个研究传统优化问题和方法的全新视角。」文中指出 Metropolis 的程序可以被用在寻找更优解中，因为物理系统的能量和一些组合优化(Combinatorial optimization)问题中的成本函数相类似，而原子的随机微小移位可类比为优化问题中解的局部变动。由此，Kirkpatrick 等人以 Metropolis 的方法为基础发展出一套随机化算法，是为模拟退火算法。
二、由固体退火到模拟退火：实现形式
模拟退火算法来源于固体退火原理，将固体加温至充分高，再使其以足够慢的速度冷却，用原子或晶格空位的移动来释放内部残留应力，通过这些原子排列重组的过程来消除材料中的差排。加温时，固体内部粒子随温度升高变为无序状，内能增大，而缓慢冷却时粒子趋于有序，在每个温度都达到平衡态，按照物理规律，最终在常温时达到基态，内能减为最小。
常规的模拟退火算法类比这一过程进行实现，包括以下步骤：
( 1 ) 设定冷却进度表(Cooling schedule)，其由以下一组初始参数组成：
冷却开始时的温度
，为避免陷入局部最优，应使此温度较高；
控制温度衰减的函数，最简单可处理为
，其中
为一个设定好的常数，经典的降温函数还有
，其中
为迭代次数；
温度的终值
，为避免精确度过差，应使此温度较低；
Markov 链的长度
，即在每一温度下的最大迭代次数。
( 2 ) 在数学模型中找到解空间与目标函数，并生成初始解。
解空间(Feasible region)
：若所有可能解均为可行解，解空间即为可能解
的集合 ；若存在不可行解，则限定解空间
为所有可行解的集合，或通常地，允许包含所有可能解但在目标函数中使用罚函数(Penalty function)以最终排除不可行解。
目标函数
：通常是解空间到实数集的映射，类比固体退火中的能量(Energy)，常表示为
，在算法内部，解越优，目标函数值
应越小。
初始解
：模拟退火算法的健壮性(Robustness)较好，即最终解的求得不过分依赖于初始解的选取，故初始解通常随机选择，并得到相应的目标函数值
。
( 3 ) 新解的产生与接受，以及最优解的储存。
对当前解
进行随机扰动产生新解
，通常通过简单变换(如部分元素的互换)产生(所有可能产生的新解构成当前解的邻域)。然后计算相应的目标函数值
，并得到
。
在经典的模拟退火算法中，根据 Metropolis 准则接受新解。物理系统总趋向于能量最低，但分子热运动则趋向于破坏这种低能量的状态，故 Metropolis 根据 Gibbs 分布提出，采样时可认为认为粒子在温度 时接受新的内能更高的状态的几率为
，其中
为温度
时的内能，
为其改变量 ， 为 Boltzmann 常数 。在算法中体现为，若
，则新解
被接受；否则新解
按概率
被接受。
由于存在退化的可能，需将历史最优解记录为
，对应历史最佳目标函数值
。
% MATLAB 下用局部代码表示接受新解的过程
if E_new < E_current
E_current = E_new;
m_current = m_new;
if E_new < E_best
E_best = E_new;
m_best = m_new;
end
else
if rand < exp(-(E_new-E_current)./t)
E_current = E_new;
m_current = m_new;
else
m_new = m_current;
end
end
( 4 ) 在每一温度
下，以 Markov 链长度
重复步骤 ( 3 )，再按照冷却进度表的设定降低温度 ，直至达到终止温度
。
三、模拟退火算法的特性分析
由于每次迭代中按一定概率接受较差解，模拟退火算法具有较强的摆脱局部最优的能力，为说明这一特点，此处介绍爬山算法与之对比。爬山算法是一种启发式的局部搜索算法，为典型的贪心算法(Greedy algorithm)。其采取局部择优的方式，随机选择一个位置作为初始点，每次朝着领域中更优的位置移动，直到到达局部最优解，而根据初始点的位置，不一定得到全局最优解。
以函数
为例，其部分函数图像如图 1 所示，若爬山算法选取的初始解大于 1，则在搜索到 A 点这一局部最优点后结束搜索，无法到达全局最优点 B 点。
图 1 一个多极值的函数图像
经典的模拟退火算法按照 Metropolis 准则，以
的概率接受恶化解，在
较大时
亦较大，即在较高的初始温度下，算法容易接受恶化解，从而能够跳出局部最优解，这大大减少了算法对初始解的依赖；而通过缓慢地降低温度，温度趋于零值时不再接受恶化解，算法可逼近全局最优解。可以证明：在模拟退火算法执行期间，随着温度参数
减小，算法返回某个整体最优解的概率单调增加，返回某个非最优解的概率单调减少；且在足够多的扰动及迭代次数下，多项式时间内模拟退火算法渐近收敛于近似最优解集。
对于局部搜索算法的平均情况，时间复杂度以问题规模
的多项式为界，而最坏情况的时间复杂度是未知的。由于冷却进度表的控制作用，模拟退火算法的时间复杂度是
，其中
为迭代次数，
是
个 Markov 链中的最大长度，
是问题规模
的多项式函数。
综上所述，模拟退火算法最终解的质量、对初始解的依赖、时间复杂度均优于局部搜索算法。
四、模拟退火算法的改进
若要模拟退火算法以概率
找到全局最优，需要在冷却进度表中满足：足够高的初始温度
、足够慢的降温速度、足够低的终止温度
、每一温度下足够的扰动；这使得收敛速度通常较慢。针对这一问题，大量学者对模拟退火算法做出了改进，本文选取其中部分进行介绍。
L. Ingber 在发表于 1989 年的文章《Very Fast Simulated Re-annealing》中提出了非常快速的模拟退火(Very Fast Simulated Annealing, VFSA)。文中谈到，传统的 Boltzmann Annealing 以及 Harold Szu 等人在 1987 年提出的 Fast Annealing 都「太过缓慢」，「这提供了发展一种新算法的动力」。
Ingber 在自己的方法中，基于似 Cauchy 分布进行扰动以产生新解。方式如下：
考虑第
次迭代中解
的第
维分量
满足
,
并与随机变量
相运算，
,
.
其中
根据均匀分布(Uniform distribution)中的
产生，即
,
.
由此 Ingber 给出新的降温函数
,
式中：
为迭代次数，
为给定常数，
为状态空间的维数。上式可改写为：
.
在实际应用中，通常选择
，
或
。
若利用依据 Tsallis 在 1988 年给出的多维分形中的 Tsallis 熵推出的广义 Gibbs 分布，可给出新的接受概率计算公式：
,
其中
为实数。显然，当
时上式即化为
.
这就是经典形式的接受概率公式。
五、一个 MATLAB 下的模拟退火算法实例
模拟退火算法的通用性很强，可以解决多类优化问题，以下将展示用模拟退火算法在 MATLAB 下解决一个非线性规划问题的实例。
非线性规划问题是求解由一系列未知实函数组成的组方程和不等式(统称为约束)定义的最优化问题，伴随着一个要被最大化或最小化的目标函数，其中一些约束或目标函数是非线性的。由于没有单纯形法(Simplex algorithm)这样的通用方法，求解非线性规划问题比求解线性规划问题要困难得多。以下是一个典型的非线性规划问题：
MATLAB 下求解此问题的模拟退火算法代码为：
x_new=ones(3,1);
x_current=[];
x_best=[];
E_new=[];
E_current=inf;
E_best=inf;
x_new(1)=20+rand*10;
x_new(2)=x_new(1)-10;
x_new(3)=10+rand*26;
markov=10^4;
T=100;
Tf=3;
x_current=x_new;
while T>Tf
rate=1;
if T<37
markov=10^5;
rate=0.05;
elseif T<10
markov=10^6;
rate=0.001;
end
for i=1:markov
r=rand;
if r<=0.25
x_new(1)=x_current(1)+rand*10*rate;
x_new(2)=x_new(1)-10;
elseif r>0.25&r<=0.5
x_new(1)=x_current(1)-rand*10*rate;
x_new(2)=x_new(1)-10;
elseif r>0.5&r<=0.75
x_new(3)=x_current(3)+rand*26*rate;
elseif r>0.75
x_new(3)=x_current(3)-rand*26*rate;
end
if -x_new(1)+2*x_new(2)+2*x_new(3)>=0 & ...
x_new(1)+2*x_new(2)+2*x_new(3)<=72,
E_new=-x_new(1)*x_new(2)*x_new(3);
if E_new
E_current=E_new;
x_current=x_new;
if E_new
E_best=E_new;
x_best=x_new;
end
else
if rand
E_current=E_new;
x_current=x_new;
else
x_new=x_current;
end
end
end
end
T=T-6;
end
solution=vpa(xbest,10)
value=vpa(-Ebest,10)
执行上面的程序，得到如下输出：
solution =
22.60638267
12.60638267
12.09036798
value =
3445.570021
多次运行后，平均得到结果值为
，平均用时为
秒。
与此对比，此类问题中，以自变量个数为问题规模
，简单的蒙特卡罗方法的渐进时间复杂度为
，达到了指数时间，而模拟退火算法始终保持在多项式时间之内。
MATLAB 下求解此问题的蒙特卡罗方法代码为：
N=1000;
x10=[]; x20=[]; x30=[];
vmax=-inf;
x1=unifrnd(20,30,N,1);
x2=x1-10;
x3=unifrnd(-10,16,N,1);
for i=1:N
for j=1:N
for k=1:N
if -x1(i)+2*x2(j)+2*x3(k)>=0 & ...
x1(i)+2*x2(j)+2*x3(k)<=72 & ...
x1(i)-x2(j)==10,
v=x1(i)*x2(j)*x3(k);
if v>vmax,
vmax=v;
x10=x1(i);
x20=x2(j);
x30=x3(k);
end
end
end
end
end
x=[x10;x20;x30];
solution=vpa(x,10)
value=vpa(vmax,10)
执行上面的程序，得到如下输出：
solution =
22.57402185
12.57402185
12.1369649
value =
3445.031902
简单的蒙特卡罗方法多次运行后平均得到结果值为
，平均用时为
秒。可以看到测试中模拟退火算法较此效率提高了
倍。
六、另一个 MATLAB 下的模拟退火算法实例
背包问题(Knapsack problem)是一种组合优化的 NP 完全问题。问题可以描述为：给定一组物品，每种物品都有自己的重量和价格，在限定的总重量内，我们如何选择，使得物品的总价格最高。
定义如下：有
种物品，物品
的重量为
，价格为
；所有物品的重量和价格都是非负的；背包所能承受的最大重量为
。若限定每种物品只能选择 0 或 1 个，则称问题为 0 - 1 背包问题。各类复杂的背包问题总可以变换为 0 - 1 背包问题进行求解，0 - 1 背包问题可以用公式表示为：
在此例中
，
，
，
。MATLAB 下求解此问题的模拟退火算法代码为：
w=[2;5;18;3;2;5;10;4;11;7;14;6];
p=[5;10;13;4;3;11;13;10;8;16;7;4];
p=-p;
Weight=46;
num=size(p,1);
x_new=ones(1,num);
x_current=x_new;
x_best=x_new;
E_current=inf;
E_best=inf;
T0=97;
Tf=3;
T=T0;
markov=10000;
while T>Tf
soltemp=x_current;
for i=1:markov
soltemp=x_current;
random=rand;
if random<0.4
temprand=ceil(rand*12);
if x_new(1,temprand)==1
x_new(1,temprand)=0;
end
end
if random>=0.4||random<=0.90
temprand=ceil(rand*12);
if x_new(1,temprand)==1;
x_new(1,temprand)=0;
temprand=ceil(rand*12);
if x_new(1,temprand)==0
x_new(1,temprand)=1;
end
end
end
if random>0.90
temprand=ceil(rand*12);
if x_new(1,temprand)==0
x_new(1,temprand)=1;
end
end
if x_new*w<=Weight
E_new=x_new*p;
if E_new
E_current=E_new;
x_current=x_new;
if E_new
E_best=E_new;
x_best=x_new;
end
else
if rand
E_current=E_new;
x_current=x_new;
else
x_new=x_current;
end
end
end
end
T=T.*0.99;
end
solution=x_best
value=-E_best
执行上面的程序，得到如下输出：
solution =
1 1 0 0 1 1 1 1 1 1 0 0
value =
76
多次运行后，平均用时为
秒，稳定得到结果值为
；在本问题中，
是可通过穷举证明的最大值。
七、结语
模拟退火算法的提出，进一步丰富了算法的种类，提供了一种物理现象与算法相结合的视角，为众多科学领域带来了一种有力的工具。后人也从多个角度，不懈地改进此算法，或是将此算法与其他算法相结合。本文从模拟退火算法的历史、形式、特点、改进和其在 MATLAB 下的实现等方面出发，介绍了这一算法。本文使用的 MATLAB 版本为 R2017a, 64bit (maci64)，结果数据在搭载 2.3 GHz Intel Core i7 处理器的计算机上运行得出。
2018.9
题图：
Wolfgang Lettl. The Man who Puts Girls on Fire, 1967.
参考文献:
1. Metropolis N, Rosenbluth A W, Rosenbluth M N, et al. Equation of state calculations by fast computing machines [J]. The journal of chemical physics, 1953, 21 (6): 1087-1092.
2. Kirkpatrick S, Gelatt C D, Vecchi M P. Optimization by simulated annealing [J]. Science, 1983, 220 (4598): 671-680.
3. 谢云. 模拟退火算法综述 [J]. 计算机应用研究, 1998, 15 (5): 7-9.
4. 刘勇, 康立山, 陈毓屏. 非数值并行算法 (第一册) 模拟退火算法 [M]. 北京: 科学出版社, 1994.
5. Ingber L. Very fast simulated re-annealing [J]. Mathematical and Computer Modelling, 1989, 12 (8): 967-973.
6. 张霖斌, 姚振兴, 纪晨, 等. 快速模拟退火算法及应用 [J]. 石油地球物理勘探, 1997, 32 (5): 654-660.
7. 卓金武, 魏永生, 秦健, 等. MATLAB 在数学建模中的应用 [M]. 北京航空航天大学出版社, 2011.