现代优化算法之模拟退火

现代优化算法

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

启发式算法包含的算法很多，例如解决复杂优化问题的蚁群算法（Ant Colony Algorithms）。有些启发式算法是根据实际问题而产生的，如解空间分解、解空间的限制等；另一类算法是集成算法，这些算法是诸多启发式算法的合成。

现代优化算法解决组合优化问题，如 TSP（Traveling Salesman Problem）问题，QAP（Quadratic Assignment Problem）问题，JSP（Job-shop Scheduling Problem）问题等效果很好。

这一章讲解模拟退火的算法过程，之前也介绍过一些简单的模拟退火的思想，上次是基于ACM-ICPC的思想进行介绍的，这次是详细的计算推导过程。

模拟退火

模拟退火算法得益于材料的统计力学的研究成果。统计力学表明材料中粒子的不同结构对应于粒子的不同能量水平。在高温条件下，粒子的能量较高，可以自由运动和重新排列。在低温条件下，粒子能量较低。如果从高温开始，非常缓慢地降温（这个过程被称为退火），粒子就可以在每个温度下达到热平衡。当系统完全被冷却时，最终形成处于低能状态的晶体。

如果用粒子的能量定义材料的状态，Metropolis 算法用一个简单的数学模型描述了退火过程。假设材料在状态ii之下的能量为E(i)E(i)，那么材料在温度TT时从状态ii进入状态jj就遵循如下规律：

（1）如果E(j)≤E(i)E(j)\leq E(i)，接受该状态被转换。

（2）如果E(j)>E(i)E(j)>E(i)，则状态转换以如下概率被接受：

eE(i)−E(j)KT

e^{E(i)-E(j)\over KT}

其中KK是物理学中的波尔兹曼常数，TT是材料温度。

在某一个特定温度下，进行了充分的转换之后，材料将达到热平衡。这时材料处于状态ii的概率满足波尔兹曼分布：

PT(x=i)=e−E(i)KT∑j∈Se−E(j)KT

P_T(x=i)=\frac{e^{-E(i)\over KT}}{\sum_{j\in S}e^{-E(j)\over KT}}
其中 xx 表示材料当前状态的随机变量, SS 表示状态空间集合。
显然

limT→∞e−E(i)KT∑j∈Se−E(j)KT=1|S|

\lim_{T\rightarrow \infty} \frac{e^{-E(i)\over KT}}{\sum_{j\in S}e^{-E(j)\over KT}} = {1\over |S|}

其中|S||S|表示集合SS中状态的数量。这表明所有状态在高温下具有相同的概率。而当温度下降时，

limT→0e−E(i)−EminKT∑j∈Se−E(j)−EminKT=limT→0e−E(i)−EminKT∑j∈Smine−E(j)−EminKT+∑j∉Smine−E(j)−EminKT=limT→0e−E(i)−EminKT∑j∈Smine−E(j)−EminKT=⎧⎩⎨⎪⎪1|Smin|0 i∈Smin other

\lim_{T\rightarrow 0}{e^{-{E(i)-E_{min}\over KT}}\over \sum_{j\in S}e^{-{E(j)-E_{min}\over KT}}}=\lim_{T\rightarrow 0}{e^{-{E(i)-E_{min}\over KT}}\over \sum_{j\in S_{min}}e^{-{E(j)-E_{min}\over KT}}+\sum_{j\notin S_{min}}e^{-{E(j)-E_{min}\over KT}}}=\\ \lim_{T\rightarrow 0}{e^{-{E(i)-E_{min}\over KT}}\over \sum_{j\in S_{min}}e^{-{E(j)-E_{min}\over KT}}}=\left\{ \begin{aligned} {1\over |S_{min}|} &\ i\in S_{min} \\ 0& \ other \end{aligned} \right.

其中Emin=minj∈SE(j)E_{min} =\min_{j\in S}E(j)且Smin=i | E(i)=EminS_{min} ={i\ |\ E(i)=E_{min}}。
上式表明当温度降至很低时,材料会以很大概率进入最小能量状态。

假定我们要解决的问题是一个寻找最小值的优化问题。将物理学中模拟退火的思想应用于优化问题就可以得到模拟退火寻优方法。

考虑这样一个组合优化问题:优化函数为 f:x→R+f : x → R^+ ,其中 x∈Sx ∈ S ,它表示优化问题的一个可行解, R+=y|y∈R,y>0R^+ ={y|y∈R,y>0}，SS表示函数的定义域。N(x)⊆SN(x)⊆S表示xx 的一个邻域集合。

首先给定一个初始温度T0T_0和该优化问题的一个初始解x(0)x(0)，并由x(0)x(0)生成下一个解x′∈N(x(0))x^\prime ∈N(x(0))，是否接受x′x^\prime作为一个新解x(1)x(1)依赖于下面概率:

P(x(0)→x′)=⎧⎩⎨1e−f(x′)−f(x(0))T0 若f(x′)<f(x(0)) other

P(x(0)\rightarrow x^\prime)=\left\{ \begin{aligned} 1 &\ \ 若f (x^\prime)

换句话说，如果生成的解 x’ 的函数值比前一个解的函数值更小,则接受 x(1)=xx(1) = x’ 作为一个新解。
否则以概率 e−f(x′)−f(x(0))T0e^{− {f (x')− f (x(0))}\over T_0} T0 接受 xx’ 作为一个新解。

泛泛地说,对于某一个温度TiT_i和该优化问题的一个解x(k)x(k), 可以生成xx’。接受xx’ 作为下一个新解 x(k+1)x(k + 1) 的概率为:

P(x(k)→x′)=⎧⎩⎨1e−f(x′)−f(x(k))T0 若f(x′)<f(x(k)) other(1)

P(x(k)\rightarrow x^\prime)=\left\{ \begin{aligned} 1 &\ \ 若f (x^\prime)

在温度TiT_i下,经过很多次的转移之后,降低温度TiT_i ,得到Ti+1<TiT_{i+1} 。在Ti+1 下重复上述过程。因此整个优化过程就是不断寻找新解和缓慢降温的交替过程。最终的解是对该问题寻优的结果。

我们注意到,在每个TiT_i 下,所得到的一个新状态x(k+1)x(k +1)完全依赖于前一个状态 x(k)x(k) , 可以和前面的状态 x(0),…,x(k−1)x(0),\dots , x(k − 1) 无关,因此这是一个马尔可夫过程。使用马尔可夫过程对上述模拟退火的步骤进行分析,结果表明:从任何一个状态 x(k ) 生成 x’ 的概率,在N(x(k))N(x(k))中是均匀分布的,且新状态x’被接受的概率满足式(1),那么经过有限次的转换,在温度TiT_i下的平衡态 xix_i 的分布由下式给出:

Pi(xi)=e−f(x−i)T∑j∈Se−f(xi)Ti

P_i(x_i)={e^{-f(x-i)\over T}\over \sum_{j\in S}e^{-{f(x_i)\over T_i}}}

当温度TT降为0 0 时, xi x_i 的分布为:

P∗i=⎧⎩⎨⎪⎪1|Smin|0 i∈Smin other

P_i^*=\left\{ \begin{aligned} {1\over |S_{min}|} &\ i\in S_{min} \\ 0& \ other \end{aligned} \right.

并且

∑xi∈SminPi=1

\sum_{x_i\in S_{min}}P_i =1

这说明如果温度下降十分缓慢,而在每个温度都有足够多次的状态转移,使之在每一个温度下达到热平衡,则全局最优解将以概率 1 被找到。因此可以说模拟退火算法可以找到全局最优解。

在模拟退火算法中应注意以下问题:

(1)理论上,降温过程要足够缓慢,要使得在每一温度下达到热平衡。但在计算机实现中,如果降温速度过缓,所得到的解的性能会较为令人满意,但是算法会太慢, 相对于简单的搜索算法不具有明显优势。如果降温速度过快,很可能最终得不到全局最优解。因此使用时要综合考虑解的性能和算法速度,在两者之间采取一种折衷。

(2)要确定在每一温度下状态转换的结束准则。实际操作可以考虑当连续 m 次的转换过程没有使状态发生变化时结束该温度下的状态转换。最终温度的确定可以提前定为一个较小的值Te ,或连续几个温度下转换过程没有使状态发生变化算法就结束。

(3)选择初始温度和确定某个可行解的邻域的方法也要恰当。

应用举例

题目

已知敌方100 个目标的经度、纬度如表1 所示。

表1 经度和纬度数据表

经度	纬度	经度	纬度	经度	纬度	经度	纬度
53.7121	15.3046	51.1758	0.0322	46.3253	28.2753	30.3313	6.9348
56.5432	21.4188	10.8198	16.2529	22.7891	23.1045	10.1584	12.4819
20.1050	15.4562	1.9451	0.2057	26.4951	22.1221	31.4847	8.9640
26.2418	18.1760	44.0356	13.5401	28.9836	25.9879	38.4722	20.1731
28.2694	29.0011	32.1910	5.8699	36.4863	29.7284	0.9718	28.1477
8.9586	24.6635	16.5618	23.6143	10.5597	15.1178	50.2111	10.2944
8.1519	9.5325	22.1075	18.5569	0.1215	18.8726	48.2077	16.8889
31.9499	17.6309	0.7732	0.4656	47.4134	23.7783	41.8671	3.5667
43.5474	3.9061	53.3524	26.7256	30.8165	13.4595	27.7133	5.0706
23.9222	7.6306	51.9612	22.8511	12.7938	15.7307	4.9568	8.3669
21.5051	24.0909	15.2548	27.2111	6.2070	5.1442	49.2430	16.7044
17.1168	20.0354	34.1688	22.7571	9.4402	3.9200	11.5812	14.5677
52.1181	0.4088	9.5559	11.4219	24.4509	6.5634	26.7213	28.5667
37.5848	16.8474	35.6619	9.9333	24.4654	3.1644	0.7775	6.9576
14.4703	13.6368	19.8660	15.1224	3.1616	4.2428	18.5245	14.3598
58.6849	27.1485	39.5168	16.9371	56.5089	13.7090	52.5211	15.7957
38.4300	8.4648	51.8181	23.0159	8.9983	23.6440	50.1156	23.7816
13.7909	1.9510	34.0574	23.3960	23.0624	8.4319	19.9857	5.7902
40.8801	14.2978	58.8289	14.5229	18.6635	6.7436	52.8423	27.2880
39.9494	29.5114	47.5099	24.0664	10.1121	27.2662	28.7812	27.6659
8.0831	27.6705	9.1556	14.1304	53.7989	0.2199	33.6490	0.3980
1.3496	16.8359	49.9816	6.0828	19.3635	17.6622	36.9545	23.0265
15.7320	19.5697	11.5118	17.3884	44.0398	16.2635	39.7139	28.4203
6.9909	23.1804	38.3392	19.9950	24.6543	19.6057	36.9980	24.3992
4.1591	3.1853	40.1400	20.3030	23.9876	9.4030	41.1084	27.7149

我方有一个基地，经度和纬度为（70,40）（70,40）。假设我方飞机的速度为10001000 公里/小时。

我方派一架飞机从基地出发，侦察完敌方所有目标，再返回原来的基地。在敌方每一目标点的侦察时间不计，求该架飞机所花费的时间（假设我方飞机巡航时间可以充分长）。

这是一个旅行商问题，旅行社问题又是NP完全问题，目前没有已知的算法可以解决。我们依次给基地编号为11，敌方目标依次编号为2，3，…，1012，3，…， 101，最后我方基地再重复编号为 102102（这样便于程序中计算）。

距离矩阵D=(dij)102×102D = (d_{ij})_{102\times 102} ，其中dijd_{ij} 表示表示i,ji, j 两点的距离，i,j=1,2,…,102，i, j = 1,2,\dots ,102，这里DD为实对称矩阵。则问题抽象成：

求一个从点11 出发，走遍所有中间点，到达点102102 的一个最短路径。

上面问题中给定的是地理坐标（经度和纬度），我们必须求两点间的实际距离。设A, B两点的地理坐标分别为(x1,y1),(x2,y2)(x_1 ,y_1 ),(x_2,y_2 )，过 A, B两点的大圆的劣弧长即为两点的实际距离。以地心为坐标原点O，以赤道平面为XOYXOY 平面，以00 度经线圈所在的平面为XOZ XOZ 平面建立三维直角坐标系。则 A, B两点的直角坐标分别为：

A(Rcosx1cosy1,Rsinx1cosy1,Rsiny1)B(Rcosx2cosy2,Rsinx2cosy2,Rsiny2)

A(R\cos x_1\cos y_1,R\sin x_1\cos y_1,R\sin y_1)\\B(R\cos x_2\cos y_2,R\sin x_2\cos y_2,R\sin y_2)

其中R=6370R = 6370为地球半径。
A, B两点的实际距离：

d=R arccos(OA→⋅OB→|OA→|⋅|OB→|)

d=R\ \arccos({\vec{OA}\cdot \vec{OB}\over |\vec{OA}|\cdot |\vec{OB}|})

化简得

d=Rarccos[cos(x1−x2)cos y1cos y2+sin y1sin y2]

d=R\arccos[cos(x_1-x_2)cos\ y_1cos\ y_2+sin\ y_1sin\ y_2]

算法描述

求解的模拟退火算法描述如下：
（1）解空间
解空间SS可表为1,2,…,101,102{1,2,\dots ,101,102 }的所有固定起点和终点的循环排列集合，即

S={(π1,…,π102)|π1=1,(π2,…,π101为为2,3,…,101的循环排列),π102=102}

S = \{( π_1,\dots ,π_{102} ) | π_1=1, (π_2 ,\dots ,π_{101}为为{2,3,\dots ,101}的循环排列) , \pi_{102}=102\}

其中每一个循环排列表示侦察100 100 个目标的一个回路，πi=jπ_i = j表示在第ii次侦察 jj点，初始解可选为(1,2,…,102)(1,2,\dots ,102)，本文中我们使用 Monte Carlo 方法求得一个较好的初始解。

（2）目标函数
此时的目标函数为侦察所有目标的路径长度或称代价函数。我们要求

minf(π1,π2,…,π102)=∑i=1101dπiπi+1

\min f(\pi_1,\pi_2,\dots,\pi_{102})=\sum_{i=1}^{101}d_{\pi_i\pi_{i+1}}

而一次迭代由下列三步构成：
（3）新解的产生
① 2 变换法
任选序号u,v（u<v） u,v（u 交换uu与vv之间的顺序，此时的新路径为：

π1…πu−1πvπv+1…πu+1πuπv+1…π102

\pi_1\dots \pi_{u-1}\pi_v\pi_{v+1}\dots \pi_{u+1}\pi_u\pi_{v+1}\dots \pi_{102}

② 3 变换法
任选序号u,vu,v 和 ww ，将uu 和vv 之间的路径插到 ww 之后，（设u<v<wu ）对应的新路径为：

π1…πu−1πv+1…πwπu…πvπw+1…π102

\pi_1\dots \pi_{u-1}\pi_{v+1}\dots \pi_w\pi_u\dots \pi_{v}\pi_{w+1}\dots \pi_{102}

（4）代价函数差
对于2 变换法，路径差可表示为

Δf=(dπu−1πv+dπuπv+1)−(dπu−1πu+dπvπv+1)

\Delta f=(d_{\pi_{u-1}\pi_v} + d_{\pi_{u}\pi_{v+1}})-(d_{\pi_{u-1}\pi_u}+d_{\pi_{v}\pi_{v+1}})

（5）接受准则

P={1exp(−∇f/T) Δf<0 Δf≥0

P=\left\{ \begin{aligned} 1 &\ \Delta f

如果Δf<0Δf ，则接受新的路径。否则，以概率exp(−Δf/T)exp(−Δf /T) 接受新的路径，即若exp(−Δf/T)exp(−Δf /T)大于 0 到1之间的随机数则接受。

（6）降温
利用选定的降温系数αα 进行降温即：T←αTT ←αT ，得到新的温度，这里我们取
α=0.999α = 0.999。

（7）结束条件
用选定的终止温度e=10−30，判断退火过程是否结束。若e = 10^{−30}，判断退火过程是否结束。若T < e$，算法结束，输出当前状态。

MATLAB程序如下:

clc,clear
load sj.txt %加载敌方100 个目标的数据，数据按照表格中的位置保存在纯文本
文件sj.txt 中
x=sj(:,1:2:8);x=x(:);
y=sj(:,2:2:8);y=y(:);
sj=[x y];
d1=[70,40];
sj=[d1;sj;d1];
sj=sj*pi/180;
%距离矩阵d
d=zeros(102);
for i=1:101
for j=i+1:102
temp=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(temp);
end
end
d=d+d';
S0=[];Sum=inf;
rand('state',sum(clock));
for j=1:1000
S=[1 1+randperm(100),102];
temp=0;
for i=1:101
temp=temp+d(S(i),S(i+1));
end
if temp<Sum
S0=S;Sum=temp;
end
end
e=0.1^30;L=20000;at=0.999;T=1;
%退火过程
for k=1:L
%产生新解
c=2+floor(100*rand(1,2));
c=sort(c);
c1=c(1);c2=c(2);
%计算代价函数值
df=d(S0(c1-1),S0(c2))+d(S0(c1),S0(c2+1))-d(S0(c1-1),S0(c1))-d(S0(c2),S0(c2+1));
%接受准则
if df<0
S0=[S0(1:c1-1),S0(c2:-1:c1),S0(c2+1:102)];
Sum=Sum+df;
elseif exp(-df/T)>rand(1)
S0=[S0(1:c1-1),S0(c2:-1:c1),S0(c2+1:102)];
Sum=Sum+df;
end
T=T*at;
if T<e
break;
end
end
% 输出巡航路径及路径长度
S0,Sum

计算结果为 44 小时左右。其中的一个巡航路径如图所示。