模擬退火算法及其Matlab實(shí)現(xiàn)

1、模擬退火算法及其 Matlab實(shí)現(xiàn)模擬退火算法 (Simulated Annealing algorithm ,簡(jiǎn)稱(chēng) SA)是柯克帕垂克 (S. Kirkpatrick ) 于1982年受熱力學(xué)中的固體退火過(guò)程與組合優(yōu)化問(wèn)題求解之間的某種“相似性”所啟發(fā)而 提出的，用于求解大規(guī)模組合優(yōu)化問(wèn)題的一種具有全局搜索功能的隨機(jī)性近似算法。與求解線(xiàn)性規(guī)劃的單純形法、Karmarkar投影尺度法，求解非線(xiàn)性規(guī)劃的最速下降法、Newton法、共軛梯度法，求解整數(shù)規(guī)劃的分支定界法、割平面法等經(jīng)典的優(yōu)化算法相比，模擬退火算法在很大程度上不受制于優(yōu)化問(wèn)題的具體形式和結(jié)構(gòu)，具有很強(qiáng)的適應(yīng)性和魯棒性，因而也具有廣泛的

2、應(yīng)用價(jià)值。模擬退火算法源于對(duì)固體退火過(guò)程的模擬；采用Metropolis接受準(zhǔn)則；并用一組稱(chēng)為冷卻進(jìn)度表的參數(shù)來(lái)控制算法進(jìn)程，使得算法在多項(xiàng)式時(shí)間里給出一個(gè)近似最優(yōu)解。固體退火過(guò)程的物理現(xiàn)象和統(tǒng)計(jì)性質(zhì)是模擬退火算法的物理背景；Metropolis接受準(zhǔn)則使算法能夠跳離局部最優(yōu)的“陷阱”，是模擬退火算法能夠獲得整體最優(yōu)解的關(guān)鍵；而冷卻進(jìn)度表的合理 選擇是算法應(yīng)用的關(guān)鍵。1物理退火過(guò)程物理中的固體退火是先將固體加熱至熔化，再徐徐冷卻，使之凝固成規(guī)整晶體的熱力 學(xué)過(guò)程。在加熱固體時(shí)，固體粒子的熱運(yùn)動(dòng)不斷增加，隨著溫度的升高，粒子與其平衡位置 的偏離越來(lái)越大，當(dāng)溫度升至溶解溫度后，固體的規(guī)則性被徹底破

3、壞，固體溶解為液體，粒子排列從較有序的結(jié)晶態(tài)轉(zhuǎn)變?yōu)闊o(wú)序的液態(tài)，這個(gè)過(guò)程稱(chēng)為溶解。溶解過(guò)程的目的是消除系統(tǒng)中原先可能存在的非均勻狀態(tài)，使隨后進(jìn)行的冷卻過(guò)程以某一平衡態(tài)為始點(diǎn)。溶解過(guò)程與系統(tǒng)的熵增過(guò)程相聯(lián)系，系統(tǒng)能量也隨溫度的升高而增大。冷卻時(shí)，液體粒子的熱運(yùn)動(dòng)漸漸減弱，隨著溫度的徐徐降低，粒子運(yùn)動(dòng)漸趨有序。當(dāng) 溫度降至結(jié)晶溫度后， 粒子運(yùn)動(dòng)變?yōu)閲@晶體格點(diǎn)的微小振動(dòng)，液體凝固成固體的晶態(tài)，這個(gè)過(guò)程稱(chēng)為退火。退火過(guò)程之所以必須“徐徐”進(jìn)行，是為了使系統(tǒng)在每一溫度下都達(dá)到平 衡態(tài)，最終達(dá)到固體的基態(tài)(圖 1-1)。退火過(guò)程中系統(tǒng)的熵值(衡量不能利用的熱能數(shù)量) 不斷減少，系統(tǒng)能量也隨溫度降低趨于最小

4、值。冷卻時(shí)，若急劇降低溫度，則將引起淬火效應(yīng)，即固體只能冷凝為非均勻的亞穩(wěn)態(tài)，系統(tǒng)能量也不會(huì)達(dá)到最小值。退火過(guò)程中系統(tǒng)在每一溫度下達(dá)到平衡態(tài)的過(guò)程，可以用封閉系統(tǒng)的等溫過(guò)程來(lái)描述。根據(jù)玻爾茲曼(Boltzmann )有序性原理，退火過(guò)程遵循應(yīng)用于熱平衡封閉系統(tǒng)的熱力學(xué)定 律自由能減少定律：“對(duì)于與周?chē)h(huán)境交換熱量而溫度保持不變的封閉系統(tǒng)，系統(tǒng)狀態(tài)的自發(fā)變化總是朝 著自由能減少的方向進(jìn)行，當(dāng)自由能達(dá)到最小值時(shí)，系統(tǒng)達(dá)到平衡態(tài)”。系統(tǒng)的自由能F = E -TS，其中E是系統(tǒng)的內(nèi)能，T是系統(tǒng)溫度，S是系統(tǒng)的熵。設(shè)i和j是恒溫系統(tǒng)的兩個(gè)狀態(tài)，即Fi = E-TS和Fj = Ej -TSj，而=Fj -

5、Fi =(Ej -EJ-TQ -S)八E-T ：S若系統(tǒng)狀態(tài)由i自發(fā)變化到j(luò)，則應(yīng)有 F ：0。顯然，能量減少(AE ： 0 )與熵增加（.S . 0）有利于自發(fā)變化。因此任一恒定溫度下， 系統(tǒng)狀態(tài)從非平衡自發(fā)變化到平衡態(tài), 都是能量和熵競(jìng)爭(zhēng)的結(jié)果，溫度決定著這兩個(gè)因素的相對(duì)權(quán)重。在高溫下，熵占統(tǒng)治地位，有利于變化的方向就是熵增加的方向，因而顯出粒子的無(wú)序狀態(tài)，而低溫對(duì)應(yīng)于低熵， 低溫下能量占優(yōu)勢(shì)，能量減少的方向有利于自發(fā)變化，因而得到有序（低熵）和低能的晶體結(jié)構(gòu)。2 Metropolis 算法固體在恒定溫度下達(dá)到熱平衡的過(guò)程可以蒙特卡羅（Monte Carlo）方法進(jìn)行模擬。MonteCar

6、lo方法的特點(diǎn)是算法簡(jiǎn)單，但必須大量采樣才能得到比較精確的結(jié)果，因而計(jì)算量很大。從物理系統(tǒng)傾向于能量較低的狀態(tài)，而熱運(yùn)動(dòng)又妨礙它準(zhǔn)確落入最低態(tài)的物理圖像出發(fā)，采樣時(shí)著重取那些有重要貢獻(xiàn)的狀態(tài)，則可以較快地得到較好的結(jié)果。1953年，Metropolis等提出了固體在恒定溫度下達(dá)到熱平衡的重要性采樣法。他們用下述方法產(chǎn)生固體的狀態(tài)序列：先給定以粒子相對(duì)位置表征的初始狀態(tài) i，作為固體的當(dāng)前狀態(tài)，該狀態(tài)的能量是 Ei。 然后用攝動(dòng)裝置使隨機(jī)選取的某個(gè)粒子的位移隨機(jī)地產(chǎn)生一微小變化，得到一個(gè)新?tīng)顟B(tài)j，新?tīng)顟B(tài)的能量是 Ej。如果Ej ： Ei，則該新?tīng)顟B(tài)就作為“重要”狀態(tài)。如果Ej Ei，則考慮到熱運(yùn)

7、動(dòng)的影響，該新?tīng)顟B(tài)是否是“重要”狀態(tài)，要根據(jù)固體處于該狀態(tài)的概率來(lái)判斷。固體處于狀態(tài)j和狀態(tài)i的概率的比值等于相應(yīng)Boltzmann因子由統(tǒng)計(jì)物理的知識(shí)，溫度為T(mén)的固體處于能量為Ej的微觀(guān)態(tài)i的概率為C exp_EikT ，其中exp-EikT稱(chēng)為玻爾茲曼（Boltzmann ）因子,k為Boltzmann常數(shù)，T為絕對(duì)溫度，C是與Ei無(wú)關(guān)的常數(shù)，上述概率分布也稱(chēng)為吉布斯（Gibbs）正則分布。容易知道，當(dāng)固體處于能量較低的微觀(guān)態(tài)的概率較大；在溫度降低時(shí)，那些能量相比最低的微觀(guān)態(tài)最有可能出現(xiàn)；當(dāng)溫度趨于零時(shí)，固體只能處于能量為最小值的基態(tài)上。的比值，即r = exp(E - Ei j)kT(1

8、.1)則r ：1。用隨機(jī)數(shù)發(fā)生器產(chǎn)生區(qū)間0,1）上服從均勻分布的隨機(jī)數(shù)，若r，則新?tīng)顟B(tài)j作為重要狀態(tài)，并以j取代i成為當(dāng)前狀態(tài)；否則舍去狀態(tài)j，仍以i為當(dāng)前狀態(tài)。重復(fù)以上新?tīng)顟B(tài)的產(chǎn)生過(guò)程。在大量遷移（固體狀態(tài)的變換稱(chēng)為遷移）后，系統(tǒng)趨于能量較低的平衡狀態(tài)，固體狀態(tài)的概率分布趨于（1.1）式的Gibbs正則分布。由（1.1）式可知，高溫下可接受與當(dāng)前狀態(tài)能差較大的新?tīng)顟B(tài)為重要狀態(tài)，而在低溫下只 能接受與當(dāng)前狀態(tài)能差較小的新?tīng)顟B(tài)為重要狀態(tài)。這與不同溫度下熱運(yùn)動(dòng)的影響完全一致， 在溫度趨于零時(shí)，就不能接受任何成立Ej a E時(shí)的新?tīng)顟B(tài)j 了。上述接受新?tīng)顟B(tài)的準(zhǔn)則稱(chēng)為Metropolis準(zhǔn)則，相應(yīng)的算

9、法被稱(chēng)為Metropolis算法。.3模擬退火算法對(duì)固體退火過(guò)程的研究給人們以新的啟示。1982年，Kirkpatrick等首先意識(shí)到固體退火過(guò)程與組合優(yōu)化問(wèn)題之間存在的類(lèi)似性，Metropolis等對(duì)固體在恒定溫度下達(dá)到熱平衡過(guò)程的模擬也給他們以啟迪：應(yīng)該把 Metropolis準(zhǔn)則引入到優(yōu)化過(guò)程中來(lái)。最終他們得到一種對(duì)Metropolis 算法進(jìn)行迭代的組合優(yōu)化算法，因這種算法模擬固體退火過(guò)程，故稱(chēng) 之為“模擬退火算法”。在模擬退火算法中，設(shè)優(yōu)化問(wèn)題的控制參數(shù)為t時(shí)的一個(gè)解乂及其非負(fù)目標(biāo)函數(shù) f （乂）分別與固體的在某一溫度下的一個(gè)微觀(guān)狀態(tài)i及其能量E等價(jià)，設(shè)隨著算法進(jìn)程遞減其值的控制參數(shù)

10、t相當(dāng)固體退火過(guò)程中的溫度的角色，則對(duì)于控制參數(shù)t的每一取值，算法持續(xù)進(jìn)行“產(chǎn)生新解一判斷一接受 /舍棄”的迭代過(guò)程就對(duì)應(yīng)著固體在某一恒定溫度下趨于熱平衡 的過(guò)程，也就是執(zhí)行了一次 Metropolis 算法。與Metropolis 算法從某一初始狀態(tài)出發(fā)，通 過(guò)計(jì)算系統(tǒng)的時(shí)間演化過(guò)程，求出系統(tǒng)最終達(dá)到的狀態(tài)相似，模擬退火算法從某個(gè)初始解出發(fā)，經(jīng)過(guò)大量解的變換后，可以求得給定控制參數(shù)值t時(shí)組合優(yōu)化問(wèn)題的相對(duì)最優(yōu)解xipt。然后減少控制參數(shù)t的值，重復(fù)執(zhí)行 Metropolis 算法，就可以在控制參數(shù) t趨于零時(shí)，最終 求得組合優(yōu)化問(wèn)題的整體最優(yōu)解。由于固體退火必須“徐徐”降溫，才能使固體在每一

11、溫度下都達(dá)到熱平衡，最終趨于能量最小的基態(tài)， 控制參數(shù)的值也必須緩慢衰減， 才能確保模擬 退火算法最終趨于組合優(yōu)化問(wèn)題的整體最優(yōu)解集。模擬退火算法用 Metropolis算法產(chǎn)生組合優(yōu)化問(wèn)題解的序列，并由與Metropolis準(zhǔn)則對(duì)應(yīng)的轉(zhuǎn)移概率P :t t、R（K 二 Xj）11 f（x；） f（x：）exp（）,f （x：）蘭 f （乂） f（x：） f （x：）（1.2）確定是否接受從當(dāng)前解x：到新解X；的轉(zhuǎn)移。式（1.2）中的r r 表示控制參數(shù)。開(kāi)始讓 t取較大的值（與固體的熔解溫度相對(duì)應(yīng)），在進(jìn)行足夠多的轉(zhuǎn)移后，緩慢減少t的值（與“徐徐”降溫相對(duì)應(yīng)），如此重復(fù)，直至滿(mǎn)足某個(gè)停止準(zhǔn)則是

12、算法終止。因此，模擬退火算法可視為 遞減控制參數(shù)值時(shí) Metropolis算法的迭代。圖1.1和圖1.2描述了固體退火過(guò)程與模擬退火 算法之間的相似性。溫度Tf基態(tài)能量圖1-1固體退火過(guò)程示意圖參數(shù)t,相對(duì)最優(yōu)解X爲(wèi)參數(shù)參數(shù) 切相對(duì)最優(yōu)解x,1opt參數(shù)“新解一判斷一接受/舍棄“新解一判斷一接受 /舍棄下降下降Metropolis算法（內(nèi)迭代）Metropolis算法（內(nèi)迭代）初始解xo目標(biāo)值f（x0）參數(shù)t0參數(shù)*下降參數(shù)f相對(duì)最優(yōu)解X爲(wèi)新解一判斷一接受 /舍棄Metropolis算法（內(nèi)迭代）圖1-2模擬退火算法流程示意圖設(shè)tk =T（tkj）表示Metropolis算法第k次迭代時(shí)控制參

13、數(shù)t的值，T t表示控制參數(shù)更新函數(shù)，tf表示終止溫度。表示Metropolis算法第k次迭代時(shí)產(chǎn)生的變換個(gè)數(shù)。下面 最小化目標(biāo)函數(shù)f（x）為例，給出模擬退火算法的具體操作步驟:(1)設(shè)置初始溫度to、終止溫度tf及控制參數(shù)更新函數(shù) T t ; 隨機(jī)產(chǎn)生初始解Xo，以此作為當(dāng)前最優(yōu)點(diǎn) x)pt =Xo，計(jì)算目標(biāo)函數(shù)值 f(Xopt)；(3)對(duì)當(dāng)前最優(yōu)點(diǎn)作一隨機(jī)變動(dòng)，產(chǎn)生一新解xk，計(jì)算新解的目標(biāo)函數(shù)值f (x)，并計(jì)算目標(biāo)函數(shù)值增量.：f = f(X)- f(Xopt)； 若.: 0 ,則接受該新解為當(dāng)前最優(yōu)點(diǎn)，人pt = x ;若厶-0 ,則以概率p的方式接受該新解為當(dāng)前最優(yōu)點(diǎn)；若k Lk，則kk 1，轉(zhuǎn)；tk二T tk二。設(shè)置令循環(huán)計(jì)數(shù)器初值k = 0；(6)若t _Tf，則轉(zhuǎn)(3)；若t : Tf，則