在求解组合优化问题中，人们首先想到的是取整，即根据连续变量的优化问题求解关联组合优化问题的变量必须为整数的约束条件，按照连续变量的优化问题对其进行求解，对于得到的结果按照某种方法取整。方法简单，但得到的结果往往违反优化问题的约束条件或得到次优结果。在取整的基础上，有学者提出了分支法。分支法的本质是根据连续变量的优化问题求解，但每次得到的结果不满足整数约束，并不是简单的舍入运算，而是通过缩小变量的可行范围，逐渐接近问题的最优解。分支法可以有效地处理组合优化问题，但每次可行域缩小都需要求解一个连续的优化问题，计算量大。同时，由于采用连续变量的优化问题求解，对所求解的数学问题有严格的数学要求。

另一种求解组合优化问题的方法为智能化方法。常见的智能化方法有粒子群算法、遗传算法、模拟退火算法等。粒子群算法根据其迭代规则，如果变量为整数变量，需要在其迭代规则的基础上进行进一步讨论。遗传算法和模拟退火算法其初始变量随机生成，可以直接转化为相关的整数，不需要做额外的变化。也就是说，遗传算法和模拟退火算法都可以求解组合优化问题。遗传算法与模拟退火算法相对比，遗传算法需要选择、交叉、变异，而模拟退火算法仅需要生成新的解，采用Metropolis准则与原解进行比较并进行相应的跟新操作。模拟退火算法迭代过程简单，容易实现。

1 模拟退火算法理论

1.1 模拟退火算法的起源

美国物理学家 N.Metropolis 和同仁在1953年发表研究复杂系统、计算其中能量分布的文章，他们使用蒙特卡罗模拟法计算多分子系统中分子的能量分布。这相当于是本文所探讨之问题的开始，事实上，模拟退火中常常被提到的一个名词就是Metropolis准则。

美国IBM公司物理学家 S.Kirkpatrick、C. D. Gelatt 和 M. P. Vecchi 于1983年在《Science》上发表了一篇颇具影响力的文章：《以模拟退火法进行最优化（Optimization by Simulated Annealing）》。他们借用了Metropolis等人的方法探讨一种旋转玻璃态系统（spin glass system）时，发觉其物理系统的能量和一些组合最优（combinatorial optimization）问题（著名的旅行推销员问题TSP即是一个代表例子）的成本函数相当类似：寻求最低成本即似寻求最低能量。由此，他们发展出以 Metropolis方法为本的一套算法，并用其来解决组合问题等的寻求最优解。

Kirkpatrick等人受到Metropolis等人用蒙特卡罗模拟的启发而发明了“模拟退火”这个名词，因为它和物体退火过程相类似。寻找问题的最优解（最值）即类似寻找系统的最低能量。因此系统降温时，能量也逐渐下降，而同样意义地，问题的解也“下降”到最值。

从理论上说，它是一种全局最优算法。模拟退火算法以优化问题的求解与物理系统退火过程的相似性为基础，它利用Metropolis算法并适当地控制温度的下降过程来实现模拟退火，从而达到求解全局优化问题的目的.模拟退火算法是一种能应用到求最小值问题的优化过程。在此过程中，每一步更新过程的长度都与相应的参数成正比，这些参数扮演着温度的角色。与金属退火原理相类似，在开始阶段为了更快地最小化，温度被升得很高，然后才慢慢降温以求稳定。

1.2 物理退火过程

模拟退火算法的核心思想与热力学的原理极为类似。在高温下，液体的大量分子彼此之间进行着相对*移动。如果该液体慢慢地冷却，热能原子可动性就会消失。大量原子常常能够自行排列成行，形成一个纯净的晶体，该晶体在各个方向上都被完全有序地排列在几百万倍于单个原子的距离之内。对于这个系统来说，晶体状态是能量最低状态，而所有缓慢冷却的系统都可以自然达到这个最低能量状态。实际上，如果液态金属被迅速冷却，则它不会达到这一状态，而只能达到一种只有较高能量的多晶体状态或非结晶状态。因此，这一过程的本质在于缓慢地冷却，以争取足够的时间，让大量原子在丧失可动性之前进行重新分布，这是确保能量达到低能量状态所必需的条件。简单而言，如下图所示，物理退火过程由以下几部分组成：加温过程（左）、等温过程（中）和冷却过程（右）。

图 1 物理退火过程

加温过程 其目的是增强粒子的热运动，使其偏离平衡位置。当温度足够高时，固体将熔解为液体，从而消除系统原先可能存在的非均匀态，使随后进行的冷却过程以某一平衡态为起点。熔解过程与系统的能量增大过程相联系，系统能量也随温度的升高而增大。

等温过程 通过物理学的知识得知，对于与周围环境交换热量而温度不变的封闭系统，系统状态的自发变化总是朝着*能减小的方向进行：当*能达到最小时，系统达到平衡态。

冷却过程 其目的是使粒子的热运动减弱并逐渐趋于有序，系统能量逐渐下降，从而得到低能量的晶体结构。

1.3 模拟退火原理

模拟退火算法来源于固体退火原理，将固体加温至充分高，再让其徐徐冷却。加温时，固体内部粒子随温升变为无序状，内能增大：而徐徐冷却时粒子渐趋有序，在每个温度上都达到平衡态，最后在常温时达到基态，内能减为最小。模拟退火算法与物理退火过程相似性对比如表1 所示。

表 1 模拟退火算法与物理退火过程相似性对比

物理退火	模拟退火
粒子状态	解
能量最低态	最优解
熔解过程	设定初温
等温过程	Metropolis采样过程
冷却	控制参数下降
能量最低态	目标函数

1953年Metropolis提出了这样一个重要性采样的方法，即设从当前状态i生成新状态j，若新状态的内能小于状态i的内能即，则接受新状态j作为新的当前状态；否则，以概率接受状态j，其中k为Boltzmann常数，这就是通常所说的Metropolis准则，通常表示为：

1.4 模拟退火算法思想

在搜索区间随机游走(即随机选择点)，再利用Metropolis抽样准则，使随机游走逐渐收敛于局部最优解。而温度是Metropolis算法中的一个重要控制参数，可以认为这个参数的大小控制了随机过程向局部或全局最优解移动的快慢。一般采用上述蒙特卡洛准则判断解是否被接受，系统能量从状态E1到状态E2，状态E2被接受的概率p的计算公式为：

公式（2）表明：如果E2<E1,则系统接受此状态，否则以一个随机概率接受或丢弃此状态。

在实际问题中，一般假设k=1，于是kT可以成为一个参数。这样经过一定次数的迭代，系统会逐渐趋于一个稳定的分布状态。重点抽样时，新状态下如果向下，则接受(局部最优)；若向上(全局搜索)，则以一定概率接受。模拟退火算法从某个初始解出发，经过大量解的变换后，可以求得给定控制参数值时组合优化问题的相对最优解。然后减小控制参数T的值，重复执行Metropolis算法，就可以在控制参数T趋于零时，最终求得组合优化问题的整体最优解。控制参数的值必须缓慢衰减。温度是Metropolis算法的一个重要控制参数，模拟退火可视为递减控制参数T时Metroplis算法的迭代。开始时T值大，可以接受较差的恶化解；随着T的减小，只能接受较好的恶化解；最后在T趋于0时，就不再接受任何恶化解了。在无限高温时，系统立即均匀分布，接受所有提出的变换。T的衰减越小，到达终点的时间越长；但可使马尔可夫(Markov) 链减小，以使到达准平衡分布的时间变短。

模拟退火算法可以求解最小值问题也可以求解最大值问题，两种问题所相关的接受概率计算公式如下。

模拟退火算法可以分解为解空间、目标函数和初始解三部分。其基本思想如下：

(1) 初始化：初始温度T(充分大)，初始解状态S(是算法迭代的起点)，每个T值的迭代次数L

(2) 对k=1, …, L做第(3)至第6步：

(3) 产生新解S'

(4) 计算增量ΔT=C(S')-C(S)，其中C(S)为评价函数

(5) 若ΔT<0则接受S′作为新的当前解，否则以概率exp-ΔT/T接受S'作为新的当前解。

(6) 如果满足终止条件则输出当前解作为最优解，结束程序。

终止条件通常取为连续若干个新解都没有被接受时终止算法。

(7) T逐渐减少，且T≥0，然后转第2步。

模拟退火算法流程图如下所示，其解的产生和接受可以划分为如下四个步骤：

（1）第一步是由一个产生函数从当前解产生一个位于解空间的新解；为便于后续的计算和接受，减少算法耗时，通常选择由当前新解经过简单地变换即可产生新解的方法，如对构成新解的全部或部分元素进行置换、互换等，注意到产生新解的变换方法决定了当前新解的邻域结构，因而对冷却进度表的选取有一定的影响。

（2）第二步是计算与新解所对应的目标函数差。因为目标函数差仅由变换部分产生，所以目标函数差的计算最好按增量计算。事实表明，对大多数应用而言，这是计算目标函数差的最快方法。

图 2模拟退火算法流程图

（3）第三步是判断新解是否被接受，判断的依据是一个接受准则，最常用的接受准则是Metropolis准则: 若Δt'<0则接受S'作为新的当前解S，否则以概率exp(-Δt'/T)接受S′作为新的当前解S。

（4）第四步是当新解被确定接受时，用新解代替当前解，这只需将当前解中对应于产生新解时的变换部分予以实现，同时修正目标函数值即可。此时，当前解实现了一次迭代。可在此基础上开始下一轮试验。而当新解被判定为舍弃时，则在原当前解的基础上继续下一轮试验。

2 实例描述

2.1 TSP旅行商问题

2.1.1 问题描述

有个旅行商，他从中心城市A出发到另外n个城市（已给出城市的经纬度坐标）出售商品，最后在回到中心城市A。如果旅行商想在路上花费的时间最短，请帮他规划路线[4]。

2.1.2 解空间

解空间 S 可表为{1,2,3,……n，n+1 ，n+2}的所有固定起点和终点的循环排列集合，其中1和n+2都表示中心城市，{2~n+1}表示旅行商经过的城市。

[1 2 3 4 5 6 7 8 9 10]和[1 9 8 7 6 1 2 3 4 5 10] 表示的就是通过8个城市旅行商的两条不同路径。另外为了计算过程更方便，我们还应提前构建表示两两城市之间的距离矩阵d，如：dij,表示城市i到城市j的距离。

2.1.3 新解的产生

任选序号 u,v (u < v )交换u 与v 之间的顺序，此时的新路径为:

产生的路径差值为：

2.1.4 目标函数

求解路径的最小值minf(1……n+2)

为了模拟这个问题，创建了30个城市，其中有一个城市是中心城市，创建的城市数据如下：

表 2 TSP问题模拟数据表

TSP问题退火模拟源代码（使用Matlab实现）：

clc,clear
x=xlsread('题目数据.xlsx','题目数据','B3:B31');
y=xlsread('题目数据.xlsx','题目数据','C3:C31');
sj=[x y]; 
dl =[120.7015202,36.37423];
sj=[dl;sj;dl];
sj=sj*pi/180; 
d=ones(31);
for i =1:31 
for j=1:31
 mp=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)=6370*acos(mp);
end
end
long=inf;
path=[];
%巡航路径及长度初始解
for j=1:300 %求较好的初始解
path0=[1,1 + randperm(29),31]; temp=0;
for i =1:30
temp = temp + d(path0(i) ,path0(i+1));
end
 if temp < long
path = path0; long = temp;
end
end
T=1;%最高温度
t=0.1^ 30;%最低温度
L = 100*30;%Markov链长度
at =0.99;%控制参数T的衰减函数参数
while t<T
    for i=1:L
        c=2 + floor(29* rand(1,2)); %产生新解
        c=sort(c); c1=c(1);c2=c(2);%计算代价函数值的增量
        df =d(path(c1 - 1),path(c2)) +d(path(c1),path(c2 +1)) -d(path(c1 - 1),path(c1)) - d( path( c2) ,path(c2 +1));
        if df<0 %接受准则
        long= long +df;
        path= [path(1:c1 -1) ,path(c2: -1:c1),path(c2 +1:31)];
        elseif exp( -df/T) >= rand(1)
        path= [path(1:c1 -1) ,path(c2: -1:c1),path(c2 +1:31)];
        long = long+df;
        end
        T=T*at;
    end
end
disp((sprintf('最终旅程为%d',long)))
sj=sj*180/pi;
xx=sj(path,1);yy = sj( path,2);
plot(xx,yy,'- *') %画出巡航路径
text(120.7015202,36.37423,'中心城市','color','k')%标注数据中心点
for i=1:29
    text(x(i),y(i),num2str(i),'color','r')%标注各点
end

基于模拟退火算法求解旅行商（TSP）问题具体步骤如下：

步骤1设置算法参数：初始温度T0，降温系数q，结束温度Tend和链长L；设置代数计数器初始化，t=O。

步骤2初始解：使用randperm函数产生随机初始线路。

步骤3解变换生成新解:采用产生随机数的方法产生两个要交换的城市，用二领域变换法产生新的路径，确定新的可行解S'。

步骤4 Metropolis准则：计算增量△t'=C(S')-c(S)，其中C(S)为评价函数。根据Metropolis准则，如果增量△t'<0，则以S'接受新的路径;如果△t '>=0，则以概率exp(-At ' /T)接受新的路径。

步骤5降温：如果满足终止条件，则输出当前解作为最优解，结束程序。否则转步骤3继续迭代。

下图是模拟退火算法解决TSP问题的结果，再30个城市的地图上画出了一个最优路线：

图 3 模拟退火解决TSP问题结果

2.2 背包问题

2.2.1 问题描述

现有一小偷入室盗窃，可以偷窃的物品数量共有N件，每件物品的价值和相应的重量已知；假设物品的价值V=[v1,…,vN]，相应的重量为W=[w1,…,wN] 。小偷在盗窃时其背包能装下物品的最大重量为R，试问，此时小偷如何选取物品在不超过背包重量前提下，使得偷盗获取的总价值最大？

2.2.2 具体实现

1）输入：背包容量为R=269，物品数量N=10，对应价值分别为

V=v1,…,vN={55，10，48，5，4，50，8，61，85，87}，物品重量分别为W=w1,…,wN={95，4，60，32，23，72，80，64，65，46}；参数设定：迭代次数（马尔科夫链长度）t，温度T，af退火率，最低温度x，背包限制T。

2）产生一个随机的合法初始解

3）产生新解，随机选取物品, 若 i 不在背包中, 则将其直接放入背包中, 或同时从背包中随机取出另一物品 j ; 若 i已在背包中, 则将其取出, 并同时随机装入另一物品j。

4）计算背包价值与合法性：计算物品价值总和，并判断是否超重。

5）接受新解：超重，放弃此解；比当前解背包价值更优，直接接受新解；比当前解背包价值更劣，概率接受新解，概率p=exp(-∆t/T)其中∆t为背包价值差。

6）更新全局最优解

7）若达到平衡，则下降温度，重置平衡次数；

8）终止条件：背包价值达到最优解或者温度下降到x度以下。

2.2.3 结果展示

表 3 模拟退火解决背包问题结果表

温度	退火率	平衡次数	迭代次数
200	0.95	5	36
200	0.8	5	22
200	0.95	10	7
200	0.95	100	1

退火率的设置对于算法效率影响十分明显，温度下降过快虽然较早达到稳定，但是有时找不到最优解，温度下降过慢会导致耗时增多，经过多次实验得出，退火率选择0.8为最佳。

平衡次数设置越多，所需迭代次数越少，但是单次迭代时间变长，所以平衡次数设置10次为最佳。初始温度设置将会影响解的搜索范围，温度越高最终解的质量越优，但是耗时也相应变长，初始温度设置200为最佳。

3 总结

模拟退火算法适用范围广，求得全局最优解的可靠性高，算法简单，便于实现；该算法的搜索策略有利于避免搜索过程因陷入局部最优解的缺陷，有利于提高求得全局最优解的可靠性。模拟退火算法有十分强的鲁棒性，这是因为比起普通的优化搜索方法，它采用许多独特的方法和技术。主要有以下几个方面：

(1)以一定的概率接受恶化解。模拟退火算法在搜索策略上不仅引入了适当的随机因素，而且还引入了物理系统退火过程的自然机理。这种自然机理的引入，使模拟退火算法在迭代过程中不仅接受使目标函数值变“好”的点，而且还能够以一定的概率接受使目标函数值变“差”的点。迭代过程中出现的状态是随机产生的，并且不强求后一状态一定优于前一状态，接受概率随着温度的下降而逐渐减小。很多传统的优化算法往往是确定性的，从一个搜索点到另一个搜索点的转移有确定的转移方法和转移关系，这种确定性往往可能使得搜索点远达不到最优点，因而限制了算法的应用范围。而模拟退火算法以一种概率的方式来进行搜索，增加了搜索过程的灵活性。

(2)引进算法控制参数。引进类似于退火温度的算法控制参数，它将优化过程分成若干阶段，并决定各个阶段下随机状态的取舍标准，接受函数由Metropolis算法给出一个简单的数学模型。模拟退火算法有两个重要的步骤：一是在每个控制参数下，由前迭代点出发，产生邻近的随机状态，由控制参数确定的接受准则决定此新状态的取舍，并由此形成一定长度的随机Markov链；二是缓慢降低控制参数，提高接受准则，直至控制参数趋于零，状态链稳定于优化问题的最优状态，从而提高模拟退火算法全局最优解的可靠性。

(3)对目标函数要求少。传统搜索算法不仅需要利用目标函数值，而且往往需要目标函数的导数值等其他一些辅助信息才能确定搜索方向：当这些信息不存在时，算法就失效了。而模拟退火算法不需要其他的辅助信息，而只是定义邻域结构，在其邻域结构内选取相邻解，再用目标函数进行评估。

模拟退火算法已经应用在各个领域用来求解最优化问题，在确保一定要求的优化质量基础上，提高模拟退火算法的搜索效率，是未来对模拟退火算法改进的主要内容。

码农公寓

0 引言