### 現代優化算法 現代優化算法是 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 小時左右。其中的一個巡航路徑如圖所示。