现代优化算法 之 模拟退火

### 现代优化算法 现代优化算法是 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）问题等效果很好。 这一章讲解模拟退火的算法过程，之前也介绍过一些简单的[模拟退火的思想](http://blog.csdn.net/u013007900/article/details/43525733)，上次是基于ACM-ICPC的思想进行介绍的，这次是详细的计算推导过程。 ### 模拟退火 模拟退火算法得益于材料的统计力学的研究成果。统计力学表明材料中粒子的不同结构对应于粒子的不同能量水平。在高温条件下，粒子的能量较高，可以自由运动和重新排列。在低温条件下，粒子能量较低。如果从高温开始，非常缓慢地降温（这个过程被称为退火），粒子就可以在每个温度下达到热平衡。当系统完全被冷却时，最终形成处于低能状态的晶体。 如果用粒子的能量定义材料的状态，Metropolis 算法用一个简单的数学模型描述了退火过程。假设材料在状态i之下的能量为E(i)，那么材料在温度T时从状态i进入状态j就遵循如下规律： - （1）如果E(j)≤E(i)，接受该状态被转换。 - （2）如果E(j)>E(i)，则状态转换以如下概率被接受： eE(i)−E(j)KT 其中K是物理学中的波尔兹曼常数，T是材料温度。 在某一个特定温度下，进行了充分的转换之后，材料将达到热平衡。这时材料处于状态i的概率满足波尔兹曼分布： PT(x=i)=e−E(i)KT∑j∈Se−E(j)KT 其中 x 表示材料当前状态的随机变量, S 表示状态空间集合。 显然 limT→∞e−E(i)KT∑j∈Se−E(j)KT=1|S| 其中|S|表示集合S中状态的数量。这表明所有状态在高温下具有相同的概率。而当温度下降时， 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 其中Emin=minj∈SE(j)且Smin=i | E(i)=Emin。 上式表明当温度降至很低时,材料会以很大概率进入最小能量状态。 假定我们要解决的问题是一个寻找最小值的优化问题。将物理学中模拟退火的思想应用于优化问题就可以得到模拟退火寻优方法。 考虑这样一个组合优化问题:优化函数为 f:x→R+ ,其中 x∈S ,它表示优化问题的一个可行解, R+=y|y∈R,y>0，S表示函数的定义域。N(x)⊆S表示x 的一个邻域集合。 首先给定一个初始温度T0和该优化问题的一个初始解x(0)，并由x(0)生成下一个解x′∈N(x(0))，是否接受x′作为一个新解x(1)依赖于下面概率: P(x(0)→x′)=⎧⎩⎨1e−f(x′)−f(x(0))T0  若f(x′)<f(x(0))  other 换句话说，如果生成的解 x’ 的函数值比前一个解的函数值更小,则接受 x(1)=x’ 作为一个新解。 否则以概率 e−f(x′)−f(x(0))T0 T0 接受 x’ 作为一个新解。 泛泛地说,对于某一个温度Ti和该优化问题的一个解x(k), 可以生成x’。接受x’ 作为下一个新解 x(k+1) 的概率为: P(x(k)→x′)=⎧⎩⎨1e−f(x′)−f(x(k))T0  若f(x′)<f(x(k))  other(1) 在温度Ti下,经过很多次的转移之后,降低温度Ti ,得到Ti+1<Ti 。在Ti+1 下重复上述过程。因此整个优化过程就是不断寻找新解和缓慢降温的交替过程。最终的解是对该问 题寻优的结果。 我们注意到,在每个Ti 下,所得到的一个新状态x(k+1)完全依赖于前一个状态 x(k) , 可以和前面的状态 x(0),…,x(k−1) 无关,因此这是一个马尔可夫过程。使用马 尔可夫过程对上述模拟退火的步骤进行分析,结果表明:从任何一个状态 x(k ) 生成 x’ 的 概率,在N(x(k))中是均匀分布的,且新状态x’被接受的概率满足式(1),那么经过有限次的转换,在温度Ti下的平衡态 xi 的分布由下式给出: Pi(xi)=e−f(x−i)T∑j∈Se−f(xi)Ti 当温度T降为0时, xi 的分布为: P∗i=⎧⎩⎨⎪⎪1|Smin|0 i∈Smin other 并且 ∑xi∈SminPi=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）。假设我方飞机的速度为1000 公里/小时。 我方派一架飞机从基地出发，侦察完敌方所有目标，再返回原来的基地。在敌方每一目标点的侦察时间不计，求该架飞机所花费的时间（假设我方飞机巡航时间可以充分长）。 这是一个[旅行商问题](http://baike.baidu.com/link?url=uqnx_j59Qy_JKzQzfdQwnOiN0I3XzFGcxa_GfHci4WG53qVm02KwgZJddmEulnpl7a2TkxSlMzbhoukR_0UdTa)，旅行社问题又是NP完全问题，目前没有已知的算法可以解决。我们依次给基地编号为1，敌方目标依次编号为2，3，…，101，最后我方基地再重复编号为 102（这样便于程序中计算）。 距离矩阵D=(dij)102×102，其中dij 表示表示i,j两点的距离，i,j=1,2,…,102，这里D为实对称矩阵。则问题抽象成： 求一个从点1出发，走遍所有中间点，到达点102的一个最短路径。 上面问题中给定的是地理坐标（经度和纬度），我们必须求两点间的实际距离。设A, B两点的地理坐标分别为(x1,y1),(x2,y2)，过 A, B两点的大圆的劣弧长即为两点的实际距离。以地心为坐标原点O，以赤道平面为XOY平面，以0度经线圈所在的平面为XOZ平面建立三维直角坐标系。则 A, B两点的直角坐标分别为： A(Rcosx1cosy1,Rsinx1cosy1,Rsiny1)B(Rcosx2cosy2,Rsinx2cosy2,Rsiny2) 其中R=6370为地球半径。 A, B两点的实际距离： d=R arccos(OA→⋅OB→|OA→|⋅|OB→|) 化简得 d=Rarccos[cos(x1−x2)cos y1cos y2+sin y1sin y2] ### 算法描述 求解的模拟退火算法描述如下： （1）解空间 解空间S可表为1,2,…,101,102的所有固定起点和终点的循环排列集合，即 S={(π1,…,π102)|π1=1,(π2,…,π101为为2,3,…,101的循环排列),π102=102} 其中每一个循环排列表示侦察100个目标的一个回路，πi=j表示在第i次侦察 j点，初始解可选为(1,2,…,102)，本文中我们使用[ Monte Carlo 方法](http://baike.baidu.com/link?url=cW0fyu21-nMWHWfwe0oKhpk25p4ep2Nnk0CwlSz70ihI-OnnKxuHwlN4ej_MvO1saIOcqlKn4-j0L4bGCk0AZq)求得一个较好的初始解。 （2）目标函数 此时的目标函数为侦察所有目标的路径长度或称代价函数。我们要求 minf(π1,π2,…,π102)=∑i=1101dπiπi+1 而一次迭代由下列三步构成： （3）新解的产生 ① 2 变换法 任选序号u,v（u<v）交换u与v之间的顺序，此时的新路径为： π1…πu−1πvπv+1…πu+1πuπv+1…π102 ② 3 变换法 任选序号u,v 和 w，将u 和v之间的路径插到 w 之后，（设u<v<w）对应的新路径为： π1…πu−1πv+1…πwπu…πvπw+1…π102 （4）代价函数差 对于2 变换法，路径差可表示为 Δf=(dπu−1πv+dπuπv+1)−(dπu−1πu+dπvπv+1) （5）接受准则 P={1exp(−∇f/T) Δf<0 Δf≥0 如果Δf<0，则接受新的路径。否则，以概率exp(−Δf/T)接受新的路径，即若exp(−Δf/T)大于 0 到1之间的随机数则接受。 （6）降温 利用选定的降温系数α 进行降温即：T←αT，得到新的温度，这里我们取 α=0.999。 （7）结束条件 用选定的终止温度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 小时左右。其中的一个巡航路径如图所示。