本節的目的是記錄以下學習和掌握模擬退火(Simulated Annealing，簡稱SA算法)過程。模擬退火算法是一種通用概率算法，用來在一個大的搜尋空間內尋找命題的最優解。這里分別使用隨機模擬退火算法和確定性模擬退火算法，在MATLAB平臺上進行編程，以尋找一個6-單元全連接網絡的能量最小化模型。 參考書籍：Richard O. Duda, Peter E. Hart, David G. Stork 著《模式分類》 ##一、技術論述 **1.隨機方法** 學習在構造模式分類器中起著中心的作用。一個通常的做法是先假設一個單參數或多參數的模型，然后根據訓練樣本來估計各參數的取值。當模型相當簡單并且低維時，可以采用解析的方法，如求函數導數，來顯式求解方程以獲得最優參數。如果模型相對復雜一些，則可以通過計算局部的導數而采用梯度下降算法來解，如人工神經網絡或其他一些最大似然方法。而對于更復雜的模型，經常會出現許多局部極值，上述方法的效果往往不盡人意。 如果一個問題越復雜，或者先驗知識和訓練樣本越少，我們對能夠自動搜索可行解的復雜搜索算法的依賴性就越強，如基于參數搜索的隨即方法。一個通常的做法是使搜索朝著預期最優解的區域前進，同時允許一定程度的隨機擾動，以發現更好的解。 **2.隨機搜索** 這里主要研究在多個候選解中搜索最優解的方法。假設給定多個變量si,i=1,2,…,N，其中每個變量的數值都取兩個離散值之一（如-1和1）。優化問題是這樣描述的：確定N個si的合適取值，時下述能量函數（又稱為代價函數）最小：  其中w_ij是一個實對稱矩陣，取值可正可負，其中令到自身的反饋權為零（即w_ii=0），這是因為非零w_ii只是在能量函數E上增加一個與si無關的常數，不影響問題的本質。這個優化問題可以用網絡和節點的方式表示，如下圖所示，其中節點之間的鏈接對應每個權值w_ij。  **3.貪心算法的局限性** 如上所述，對于求解有N個變量si的能量E最小化問題，除非N的取值很小，否則往往無法直接求解，因為其構型數目高達N^2。如果使用貪心算法搜索最優的構型，需要先隨機選取每個節點的起始狀態，然后順序考查每個節點從而計算與之相聯系的si=1狀態和si=-1狀態的能量，最后選取能夠降低能量的狀態遷移。這種判斷只用到了直接與之相連的具有非零權值的相鄰節點。 這種貪心搜索算法成功找到最優解的可能性很小，因為系統常常會陷入局部能量最小值，或根本就不收斂，因此需要選擇其他搜索方法。 **4.模擬退火** 在熱力學中，固體的退火過程主要由以下三部分組成： 1. 加溫過程。其目的是增強粒子的熱運動，使其偏離平衡位置。當溫度足夠高時，固體將熔解為液體，從而消除系統原先可能存在的非均勻態，使隨后進行的冷卻過程以某一平衡態為起點。熔解過程實際是系統的熵增過程，系統能量也隨溫度的升高而增大。 2. 等溫過程。物理學的知識告訴我們，對于與周圍環境交換熱量而溫度不變的封閉系統，系統狀態的自發變化總是朝自由能減少的方向進行，當自由能達到最小時，系統達到平衡態。 3. 冷卻過程。其目的是使粒子的熱運動減弱并漸趨有序，系統能量逐漸下降，從而得到低能的晶體結構。 Metropolis等在1953年提出了重要性采樣法，即以概率大小接受新狀態。具體而言，在溫度T時, 由當前狀態i產生新狀態j，兩者的能量分別為Ei和Ej，若Ej小于Ei，則接受新狀態j為當前狀態；否則，計算概率p(?E)：  若p(?E)大于[0,1]區間內的隨機數，則仍舊接受新狀態j為當前狀態；若不成立則保留i為當前狀態，其中k為玻爾茲曼常數，T為系統溫度。上述重要性采樣過程通常稱為Metropolis準則： 1. 在高溫下可接受與當前狀態能量差較大的新狀態； 2. 而在低溫下基本只接受與當前能量差較小的新狀態； 3. 且當溫度趨于零時，就不能接受比當前狀態能量高的新狀態。 1983年Kirkpatrick 等意識到組合優化與物理退火的相似性，并受到Metropolis 準則的啟迪，提出了模擬退火算法。模擬退火算法是基于Monte Carlo 迭代求解策略的一種隨機尋優算法，其出發點是基于物理退火過程與組合優化之間的相似性，模擬退火方法由某一較高初溫開始，利用具有概率突跳特性的Metropolis抽樣策略在解空間中進行隨機搜索，伴隨溫度的不斷下降，重復抽樣過程，最終得到問題的全局最優。對比貪心算法，模擬退火算法主要的優勢在于它使系統有可能從局部最小處跳出。 對于一個優化問題： 把優化問題的求解過程與統計熱力學的熱平衡問題進行類比，通過模擬高溫物體退火過程的方法，試圖找到優化問題的全局最優或近似全局最優解； 允許隨著參數的調整，目標函數可以偶爾向能量增加的方向發展（對應于能量有時上升），以利于跳出局部極小的區域，隨著假想溫度的下降（對應于物體的退火），系統活動性降低，最終穩定在全局最小所在的區域。 **5.兩種模擬退火算法** 兩種模擬退火算法，即隨機模擬退火和和確定性模擬退火算法的實現步驟如下所示：  隨機模擬退火算法收斂很慢，部分原因在于其中搜索的全部的構型空間的離散本質，即構型空間是一個N維超立方體。每一次搜索軌跡都只能沿著超立方體的一條邊，狀態只能落在超立方體的頂點上，因此失去了完整的梯度信息。而梯度信息是可以用超立方體內部的連續狀態值提供的。一種更快的方法就是以下的確定性模擬退火算法：  ##二、實驗結果討論 構造一個6-單元全連接網絡，能量函數使用公式：  其中網絡的連接權值矩陣如下：  設計步驟主要包括以下幾個部分：? 編寫程序[E, s_out] = RandomSimulatedAnnealing(T_max, time_max, c, s, w)，實現以上算法1所述的隨機模擬退火算法。這里需要設定以下參數： T_max=10，T(m+1) =c*T(m)，c=0.9，進行實驗，能量隨溫度下降次數的變化曲線如圖2所示（由于模擬退火算法所得到的結果有一定的隨機性，因此以下步驟均執行四次算法進行觀察），四次所得到的最終構型s如圖3所示。 改變參數：初始溫度：T_max=5，T(m+1) =c*T(m)，c=0.5，進行實驗，能量隨溫度下降次數的變化曲線如圖4所示，四次所得到的最終構型s如圖5所示。 編寫程序[E, s_out] = DeterministicAnnealing(T_max, time_max, c, s, w)，實現以上算法2所述的確定性模擬退火算法。這里需要設定以下參數： T_max=10，T(m+1) =c*T(m)，c=0.9，進行實驗，能量隨溫度下降次數的變化曲線如圖6所示，四次所得到的最終構型s如圖7所示。 改變參數：初始溫度：T_max=5，T(m+1) =c*T(m)，c=0.5，進行實驗，能量隨溫度下降次數的變化曲線如圖8所示，四次所得到的最終構型s如圖9所示。 結論：圖2、3給出了多次隨機模擬退火算法的運行結果，可以看到構型s不一定完全一樣；能量函數E的波形在經過若干次逐漸遞減的震蕩后基本收斂到全局最小值-19。當改變T(1)=5,c=0.5時，從圖4中可觀察到能量函數E的波形極速下降，并達到較小的值，中間少了一個溫度漸變和震蕩調整的過程。 圖6、7給出多次確定性模擬退火算法的運行結果，且每次得到的最終構型s均一致；能量函數E的波形在經過平緩遞減后收斂到全局最小值-19，不出現隨機模擬退火中劇烈震蕩的情況。當改變T(1)=5,c=0.5時，從圖8中可觀察到能量函數E的波形同樣呈現出極速下降的態勢，并達到較小的值，中間少了一個溫度漸變和調整的過程。 綜合以上的實驗結果，我們發現隨機退火和確定性退火均能給出相似的最終解，但對于一些大規模的現實問題，隨機模擬退火的運行速度很慢，而相比之下確定性退火算法要快很多，有時可以快2~3個數量級。 另外，初溫T(1)和溫度下降系數c的選擇對算法性能也有很大影響。初溫的確定應折衷考慮優化質量和優化效率，常用方法包括以下幾種： * 均勻抽樣一組狀態，以各狀態目標值的方差為初溫； * 隨機產生一組狀態，確定兩兩狀態間的最大目標值差|?max|，然后依據差值，利用一定的函數確定初溫。譬如T(1)=-?/ln? p_r，其中p_r為初始接受概率； * 利用經驗公式給出。 模擬退火算法設計中包括三個重要的函數：狀態產生函數、狀態接受函數、溫度更新函數；同時在程序設計時，需遵循內循環終止準則、外循環終止準則。這些環節的設計將決定模擬退火算法的優化性能。 ##三、實驗結果     ##四、簡單代碼實現 ~~~ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % 隨機模擬退火函數 % 輸入參數： % T_max：初始溫度 % time_max：最大迭代次數 % c：溫度下降比率 % s：初始構型 % w：權值矩陣 % 輸出參數： % E：能量變化矩陣 % s_out：經過算法計算后的構型 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function [E, s_out] = RandomSimulatedAnnealing(T_max, time_max, c, s, w) [x, y] = size(s); time = 1; % 迭代次數 T(time) = T_max; % 初始溫度設置 while (time < (time_max + 1)) % (T(time) > T_min) for i = 1:1000 num = ceil(rand(1) * y); % 選擇產生一個1到y之間的隨機數 for j = 1:y e(j) = w(num, j) * s(num) * s(j); end Ea(time) = -1 / 2 * sum(e); Eb(time) = -Ea(time); if Eb(time) < Ea(time) s(num) = -s(num); elseif (exp(-(Eb(time) - Ea(time)) / T(time)) > rand()) s(num) = -s(num); else s(num) = s(num); end end % 計算能量E E(time) = 0; for it = 1:6 for jt = 1:6 E(time) = E(time) + w(it, jt) * s(it) * s(jt); end end E(time) = E(time) * (-0.5); s_out(time,:) = s; % 每次形成的構型 time = time + 1; T(time) = T(time - 1) * c; end ~~~ ~~~ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % 確定性模擬退火函數 % 輸入參數： % T_max：初始溫度 % time_max：最大迭代次數 % c：溫度下降比率 % s：初始構型 % w：權值矩陣 % 中間函數： % tanh(l / T)：響應函數，該函數有一個隱含的重新規格化的作用 % 輸出參數： % E：能量變化矩陣 % s_out：經過算法計算后的構型 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function [E, s_out] = DeterministicAnnealing(T_max, time_max, c, s, w) [x, y] = size(s); time = 1; % 迭代次數 T(time) = T_max; % 初始溫度設置 while (time < (time_max + 1)) num = ceil(rand(1) * y); % 選擇產生一個1到y之間的隨機數 for j = 1:y e(j) = w(num, j) * s(j); end l(time) = sum(e); s(num) = tanh(l(time) / T(time)); % 計算能量E E(time) = 0; for it = 1:6 for jt = 1:6 E(time) = E(time) + w(it, jt) * s(it) * s(jt); end end E(time) = E(time) * (-0.5); s_out(time,:) = s; % 每次形成的構型 time = time + 1; T(time) = T(time - 1) * c; end ~~~ ~~~ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %模擬退火算法實驗 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% clear; close all; % 網絡的連接權值矩陣 w = [ 0 5 -3 4 4 1;... 5 0 -1 2 -3 1;... -3 -1 0 2 2 0;... 4 2 2 0 3 -3;... 4 -3 2 3 0 5;... 1 1 0 -3 5 0]; num = 6; % 總共產生6個數 s_in = rand(1,num); % 生成1和-1的隨機序列 s_in(s_in > 0.5) = 1; s_in(s_in < 0.5) = -1; disp(['初始構型S為：',num2str(s_in)]); % 以下是隨機模擬退火算法 T_max = 10; % 初始溫度設置 time_max = 100; % 最大迭代次數 c = 0.9; % 溫度變化比率 [E1, s_out1] = RandomSimulatedAnnealing(T_max, time_max, c, s_in, w); subplot(221),plot(E1);grid on; title(['T(1) = ',num2str(T_max),',c = ',num2str(c),',隨機模擬退火算法能量變化曲線']); disp(['T(1) = 10,c = 0.9,隨機模擬退火算法最終構型S為：',num2str(s_out1(time_max,:))]); T_max = 5; % 初始溫度設置 time_max = 100; % 最大迭代次數 c = 0.5; % 溫度變化比率 [E2, s_out2] = RandomSimulatedAnnealing(T_max, time_max, c, s_in, w); subplot(222),plot(E2);grid on; title(['T(1) = ',num2str(T_max),',c = ',num2str(c),',隨機模擬退火算法能量變化曲線']); disp(['T(1) = 5,c = 0.5,隨機模擬退火算法最終構型S為：',num2str(s_out2(time_max,:))]); % 以下是確定性模擬退火算法 T_max = 10; % 初始溫度設置 time_max = 100; % 最大迭代次數 c = 0.9; % 溫度變化比率 [E3, s_out3] = DeterministicAnnealing(T_max, time_max, c, s_in, w); subplot(223),plot(E3);grid on; title(['T(1) = ',num2str(T_max),',c = ',num2str(c),',確定性模擬退火算法能量變化曲線']); disp(['T(1) = 10,c = 0.9,確定性模擬退火算法最終構型S為：',num2str(s_out3(time_max,:))]); T_max = 5; % 初始溫度設置 time_max = 100; % 最大迭代次數 c = 0.5; % 溫度變化比率 [E4, s_out4] = DeterministicAnnealing(T_max, time_max, c, s_in, w); subplot(224),plot(E4);grid on; title(['T(1) = ',num2str(T_max),',c = ',num2str(c),',確定性模擬退火算法能量變化曲線']); disp(['T(1) = 5,c = 0.5,確定性模擬退火算法最終構型S為：',num2str(s_out4(time_max,:))]); ~~~ 參考鏈接：[http://www.cnblogs.com/growing/archive/2010/12/16/1908255.html](http://www.cnblogs.com/growing/archive/2010/12/16/1908255.html)