類比退火演算法——解決售貨員的難題

Simulation Annealing
    1982年,KirkPatrick將退火思想引入組合最佳化領域，提出一種解大規模組合最佳化問題的演算法，對NP完全組合最佳化問題尤其有效。這源於固體的退火過程，即先將溫度加到很高,再緩慢降溫(即退火)，使達到能量最低點。如果急速降溫(即為淬火)則不能達到最低點.。
   類比退火演算法是一種能應用到求最小值問題或基本先前的更新的學習過程（隨機或決定性的）。在此過程中，每一步更新過程的長度都與相應的參數成正比，這些參數扮演著溫度的角色。然後，與金屬退火原理相類似，在開始階段為了更快地最小化或學習，溫度被升得很高，然後才（慢慢）降溫以求穩定。
類比退火演算法是一種用於求解大規模最佳化問題的隨機搜尋演算法，它以最佳化問題求解過程與物理系統退火過程之間的相似性為基礎；最佳化的目標函數相當於金屬的內能；最佳化問題的自變數組合狀態空間相當於金屬的內能狀態空間；問題的求解過程就是找一個組合狀態，使目標函數值最小。利用Metropolis準則並適當地控制溫度的下降過程實現類比退火，從而達到在多項式時間內求解全域最佳化問題的目標。
   類比退火演算法來源於固體退火原理，將固體加溫至充分高，再讓其徐徐冷卻，加溫時，固體內部粒子隨溫升變為無序狀，內能增大，而徐徐冷卻時粒子漸趨有序，在每個溫度都達到平衡態，最後在常溫時達到基態，內能減為最小。根據Metropolis準則，粒子在溫度T時趨於平衡的機率為e-ΔE/(kT)，其中E為溫度T時的內能，ΔE為其改變數，k為Boltzmann常數。用固體退火類比組合最佳化問題，將內能E類比為目標函數值f，溫度T演化成控制參數t，即得到解組合最佳化問題的類比退火演算法：由初始解i和控制參數初值t開始，對當前解重複“產生新解→計算目標函數差→接受或捨棄”的迭代，並逐步衰減t值，演算法終止時的當前解即為所得近似最優解，這是基於蒙特卡羅迭代求解法的一種啟發學習法隨機搜尋過程。退火過程由冷卻進度表(Cooling
Schedule)控制，包括控制參數的初值t及其衰減因子Δt、每個t值時的迭代次數L和停止條件S。
類比退火演算法的模型
　　類比退火演算法可以分解為解空間、目標函數和初始解三部分。
　類比退火的基本思想:
　　(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步。
演算法對應動態示範圖：
類比退火演算法新解的產生和接受可分為如下四個步驟：
　　第一步是由一個產生函數從當前解產生一個位於解空間的新解；為便於後續的計算和接受，減少演算法耗時，通常選擇由當前新解經過簡單地變換即可產生新解的方法，如對構成新解的全部或部分元素進行置換、互換等，注意到產生新解的變換方法決定了當前新解的鄰域結構，因而對冷卻進度表的選取有一定的影響。
　　第二步是計算與新解所對應的目標函數差。因為目標函數差僅由變換部分產生，所以目標函數差的計算最好按增量計算。事實表明，對大多數應用而言，這是計算目標函數差的最快方法。
　　第三步是判斷新解是否被接受,判斷的依據是一個接受準則，最常用的接受準則是Metropo1is準則: 若Δt′<0則接受S′作為新的當前解S，否則以機率exp(-Δt′/T)接受S′作為新的當前解S。
　　第四步是當新解被確定接受時，用新解代替當前解，這隻需將當前解中對應於產生新解時的變換部分予以實現，同時修正目標函數值即可。此時，當前解實現了一次迭代。可在此基礎上開始下一輪實驗。而當新解被判定為捨棄時，則在原當前解的基礎上繼續下一輪實驗。
　　類比退火演算法與初始值無關，演算法求得的解與初始解狀態S(是演算法迭代的起點)無關；類比退火演算法具有漸近收斂性，已在理論上被證明是一種以機率l 收斂於全域最優解的全域最佳化演算法；類比退火演算法具有並行性。
類比退火演算法的簡單應用
　　作為類比退火演算法應用，討論貨郎擔問題(Travelling Salesman Problem，簡記為TSP)：設有n個城市，用數位1,…,n代表。城市i和城市j之間的距離為d(i，j) i, j=1,…,n．TSP問題是要找遍訪每個域市恰好一次的一條迴路，且其路徑總長度為最短.。
　　求解TSP的類比退火演算法模型可描述如下：
　　解空間 解空間S是遍訪每個城市恰好一次的所有迴路，是{1，……，n}的所有迴圈排列的集合，S中的成員記為(w1,w2 ,……，wn)，並記wn+1= w1。初始解可選為(1，……，n)
　　目標函數 此時的目標函數即為訪問所有城市的路徑總長度或稱為代價函數：
　　我們要求此代價函數的最小值。
　　新解的產生 隨機產生1和n之間的兩相異數k和m，若k，則將
　　(w1, w2 ,…，wk , wk+1 ,…，wm ,…，wn)
　　變為：
　　(w1, w2 ,…，wm , wm-1 ,…，wk+1 , wk ,…，wn).
　　如果是k>m，則將
　　(w1, w2 ,…，wk , wk+1 ,…，wm ,…，wn)
　　變為：
　　(wm, wm-1 ,…，w1 , wm+1 ,…，wk-1 ,wn , wn-1 ,…，wk).
　　上述變換方法可簡單說成是“逆轉中間或者逆轉兩端”。
　　也可以採用其他的變換方法，有些變換有獨特的優越性，有時也將它們交替使用，得到一種更好方法。
　　代價函數差 設將(w1, w2 ,……，wn)變換為(u1, u2 ,……，un), 則代價函數差為：
根據上述分析，可寫出用類比退火演算法求解TSP問題的偽程式：
Procedure TSPSA:
　begin
　　init-of-T; { T為初始溫度}
　　S={1，……，n}; {S為初始值}
　　termination=false;
　　while termination=false
　　　begin
　　　　for i=1 to L do
　　　　　　begin
　　　　　　　　generate(S′form S); { 從當前迴路S產生新迴路S′}
　　　　　　　　Δt:=f(S′))-f(S);{f(S)為路徑總長}
　　　　　　　　IF(Δt<0) OR (EXP(-Δt/T)>Random-of-[0,1])
　　　　　　　　S=S′;
　　　　　　　　IF the-halt-condition-is-TRUE THEN
　　　　　　　　termination=true;
　　　　　　End;
　　　　T_lower;
　　　End;
　End
類比退火演算法的應用很廣泛，可以較高的效率求解最大截問題(Max Cut Problem)、0-1背包問題(Zero One Knapsack Problem)、圖著色問題(Graph Colouring Problem)、調度問題(Scheduling Problem)等等。
類比退火演算法的參數控制問題
　　類比退火演算法的應用很廣泛，可以求解NP完全問題，但其參數難以控制，其主要問題有以下三點：
　　(1) 溫度T的初始值設定問題。
　　溫度T的初始值設定是影響類比退火演算法全域搜尋效能的重要因素之一、初始溫度高，則搜尋到全域最優解的可能性大，但因此要花費大量的計算時間；反之，則可節約計算時間，但全域搜尋效能可能受到影響。實際應用過程中，初始溫度一般需要依據實驗結果進行若干次調整。
　　(2) 退火速度問題。
　　類比退火演算法的全域搜尋效能也與退火速度密切相關。一般來說，同一溫度下的“充分”搜尋(退火)是相當必要的，但這需要計算時間。實際應用中，要針對具體問題的性質和特徵設定合理的退火平衡條件。
　　(3) 溫度管理問題。
　　溫度管理問題也是類比退火演算法難以處理的問題之一。實際應用中，由於必須考慮計算複雜度的切實可行性等問題，常採用如下所示的降溫方式：
T(t+1)＝k×T(t)
式中k為正的略小於1.00的常數，t為降溫的次數。

類比退火演算法的主要思想
就函數最小值問題來說，類比退火的主要思想是：在搜尋區間（二維平面中）隨機遊走（即隨機播放點），再以Metropolis抽樣準則，使隨機遊走逐漸收斂於局部最優解。而溫度即是Metropolis演算法中的一個重要控制參數，可以認為這個參數的大小控制了隨時過程向局部或全域最優解移動的快慢。
冷卻參數表、領域結構和新解產生器、接受準則和隨機數產生器（即Metropolis演算法）一起構成演算法的三大支柱。
l 重點抽樣與Metroplis演算法：
Metropolis是一種有效重點抽樣法，其演算法為：系統從能量一個狀態變化到另一個狀態時，相應的能量從E1變化到E2，機率為p = exp[ - (E2- E1)/kT ]。如果E2 < E1，系統接收此狀態，否則，以一個隨機的機率接收此或丟棄此狀態。這種經常一定次數的迭代，系統會逐漸趨於一引穩定的分布狀態。
重點抽樣時，新狀態下如果向下則接受（局部最優），若向上（全域搜尋），以一定機率接受。類比退火方法從某個初始解出發,經過大量解的變換後,可以求得給定控制參數值時組合最佳化問題的相對最優解。然後減小控制參數T的值，重複執行Metropolis演算法，就可以在控制參數T趨於零時，最終求得組合最佳化問題的整體最優解。控制參數的值必須緩慢衰減。
其中溫度是一個Metropolis的重要控制參數，類比退火可視為遞減控制參數什時Metroplis演算法的迭代。開始T值大，可能接受較差的惡化解，隨著T的減小，只能接受較好的惡化解，最後在T趨於0時，就不再接受任何惡化解了。
在無限高溫時，系統立即均勻分布，接受所有提出的變換。T的衰減越小，T到達終點的時間越長；但可使馬可夫鏈越小，到達准平衡分布的時間越短，
l 參數的選擇：
我們稱調整類比退火法的一系列重要參數為冷卻進度表。它控制參數T的初值及其衰減函數，對應的MARKOV鏈長度和停止條件，非常重要。
一個冷卻進度表應當規定下述參數：
1． 控制參數t的初值t0;
2． 控制參數t的衰減函數；
3． 馬爾可夫鏈的長度Lk。（即每一次隨機遊走過程，要迭代多少次，才能趨於一個準平衡分布，即一個局部收斂解位置）
4． 結束條件的選擇
有效冷卻進度表判據：
一．演算法的收斂：主要取決於衰減函數和馬可夫鏈的長度及停止準則的選擇
二．演算法的實驗效能：最終解的品質和CPU的時間
參數的選擇：
一）控制參數初值T0的選取
一般要求初始值t0的值要充分大，即一開始即處於高溫狀態，且Metropolis的接收率約為1。
二）衰減函數的選取　
衰減函數用於控制溫度的退火速度，一個常用的函數為：T(n + 1) = K*T(n)，其中K是一個非常接近於1的常數。
三）馬可夫鏈長度L的選取
原則是：在衰減參數T的衰減函數已選定的前提下，L應選得在控制參數的每一取值上都能恢複准平衡。
四）終止條件
有很多種終止條件的選擇，各種不同的條件對演算法的效能和解的品質有很大影響，本文只介紹一個常用的終止條件。即上一個最優解與最新的一個最優解的之差小於某個容差，即可停止此次馬爾可夫鏈的迭代。
使用類比退火法求函數f(x,y) = 5sin(xy) + x2 + y2的最小值
解：根據題意，我們設計冷卻表進度表為：
即初始溫度為100
衰減參數為0.95
馬可夫鏈長度為10000
Metropolis的步長為0.02
結束條件為根據上一個最優解與最新的一個最優解的之差小於某個容差。

 

 

 

售貨員的難題

Time Limit:1000MS  Memory Limit:65536K
Total Submit:153 Accepted:10

Description

某鄉有N個村莊(1<N<40),有一個售貨員,他要到各個村莊去售貨,各個村莊之間的路程S(0<S<10000)是已知的,且A村到B村與B村到A村的路大多不同.為了提高效率,他從商店出發到每個村莊一次,然後返回商店所在的村,假設商店所在的村莊為1,他不知道選擇什麼樣的路線才能使所走過的路程最短.請你幫他選擇一條最短的路.

Input

村莊N和各個村之間的路程(均是整數).

Output

最短的路程.

Sample Input

30 2 11 0 22 1 0

 

Sample Output

3 

 

 

CODE：

/*類比退火*//*AC代碼：0ms*/#include <iostream>#include <ctime>#include <cstdlib>using namespace std;const int MAX=41;const int RANN=1000;const int RUNN=50;const int INF=99999999;int map[MAX][MAX],rpath[RANN][MAX],min[RANN],N;void adjust(int x[],int rn)//rn為調整次數{int a,b;while(rn--){a=rand()%(N-1)+1;//調整的位置b=rand()%(N-1)+1;swap(x[a],x[b]);}}void get_map(){int i,j;for(i=1;i<=N;i++)for(j=1;j<=N;j++)scanf("%d",&map[i][j]);}void get_rpath()//產生RANN組初始數列{int i,j;for(i=0;i<RANN;i++){for(j=1;j<=N-1;j++)rpath[i][j]=j+1;adjust(rpath[i],N-1);//初始調整min[i]=INF;}}void swap(int &x,int &y)//交換x,y{int t=x;x=y;y=t;}int get_dis(int x[])//返迴路徑長度{int sum=0,i,p=1;for(i=1;i<=N-1;i++){sum+=map[p][x[i]];p=x[i];}sum+=map[p][1];return sum;}void numcpy(int x[],int y[]){for(int i=1;i<=N-1;i++)x[i]=y[i]; }void get_ans(){int i,j,t=N-1,temp,p[41];while(t--)//調整的範圍遞減{for(i=0;i<RANN;i++)//遍曆RanN組資料{for(j=0;j<RUNN;j++)//對於每組做RunN次{numcpy(p,rpath[i]);//因為每次都在原有的rpath[i]上改變adjust(p,t);temp=get_dis(p);if(temp<min[i]){numcpy(rpath[i],p);//更新rpath[i];min[i]=temp;}}}}int ans=min[0];for(i=1;i<RANN;i++)if(min[i]<ans)ans=min[i];printf("%d/n",ans);}int main (){srand(time(0));//以時間為隨機數種子產生隨機數while(scanf("%d",&N)!=EOF){get_map();get_rpath();get_ans();}return 0;}