應(yīng)用篇-第12章-模擬退火算法的仿真與實(shí)現(xiàn)

模擬退火算法（Simulated

Annealing

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

模擬退火算法的仿真與實(shí)現(xiàn)與實(shí)現(xiàn)12.1.1物理退火過程

固體退火是先將固體加熱至熔化，再徐徐冷卻使之凝固成規(guī)整晶體的熱力學(xué)過程。退火過程是系統(tǒng)在每一溫度下達(dá)到平衡的過程可以用封閉系統(tǒng)的等溫過程來描述。由Boltzmann有序性原理，退火過程遵循熱平衡封閉系統(tǒng)的熱力學(xué)定律——自由能減少定律：對(duì)于與周圍環(huán)境交換能量而溫度保持不變的封閉系統(tǒng)，系統(tǒng)狀態(tài)的自發(fā)變化總是朝著自由能減少的方向進(jìn)行當(dāng)自由能達(dá)到最小值時(shí)，系統(tǒng)達(dá)到平衡態(tài)。

設(shè)E為系統(tǒng)的微觀狀態(tài)的能量，則系統(tǒng)處在狀態(tài)i的概率為：其中A是

無關(guān)的常數(shù)。(2)式Gibbs正則分布.。該分布給出溫度T時(shí)固體處于能量

的微觀態(tài)i的概率。顯然，固體處于能量較低的微觀態(tài)概率大，在溫度降低時(shí)，那些能量相對(duì)較低的微觀態(tài)最有可能出現(xiàn)。當(dāng)溫度趨于零時(shí)，固體基本上只能處于能量為最小值的基態(tài)。

12.1.2Metropolis準(zhǔn)則

固體在恒定溫度下達(dá)到平衡的過程可以用MonteCarl方法加以模擬，該方法雖然簡(jiǎn)單，但必須大量采樣才能得到比較精確的結(jié)果，因而計(jì)算量很大。1953年，Metropolis等提出了重要性采樣法，即以概率接受新狀態(tài)。具體而言：

在溫度為T，由當(dāng)前狀態(tài)i產(chǎn)生新狀態(tài)j，兩者的能量分別為

和，若<則接受新產(chǎn)狀態(tài)j為當(dāng)前狀態(tài)；否則，若概率

大于[0,1]區(qū)間內(nèi)的隨機(jī)數(shù)接受新狀態(tài)J為當(dāng)前狀態(tài)，若不成立則保留狀態(tài)i為當(dāng)前狀態(tài)。當(dāng)這種過程多次重復(fù)，即經(jīng)過大量遷移后系統(tǒng)趨于能量較低的平衡態(tài)，固體狀態(tài)的概率分布趨于(2)式的Gibbs分布。上述接受新狀態(tài)的準(zhǔn)則稱為Metropolis準(zhǔn)則，相應(yīng)的算法稱為Metropolis算法，這種算法的計(jì)算顯著減少。

12.1.3模擬退火算法介紹

模擬退火算法(SimulatedAnnealing，SA)最早由Kirkpatrick等應(yīng)用于組合優(yōu)化領(lǐng)域，它是基于Monte-Carlo迭代求解策略的一種隨機(jī)尋優(yōu)算法，其出發(fā)點(diǎn)是基于物理中固體物質(zhì)的退火過程與一般組合優(yōu)化問題之間的相似性。模擬退火算法從某一較高初溫出發(fā)，伴隨溫度參數(shù)的不斷下降,結(jié)合概率突跳特性在解空間中隨機(jī)尋找目標(biāo)函數(shù)的全局最優(yōu)解，即在局部最優(yōu)解能概率性地跳出并最終趨于全局最優(yōu)。模擬退火算法是一種通用的優(yōu)化算法，理論上算法具有概率的全局優(yōu)化性能,目前已在工程中得到了廣泛應(yīng)用，諸如VLSI、生產(chǎn)調(diào)度、控制工程、機(jī)器學(xué)習(xí)、神經(jīng)網(wǎng)絡(luò)、信號(hào)處理等領(lǐng)域。模擬退火算法是通過賦予搜索過程一種時(shí)變且最終趨于零的概率突跳性，從而可有效避免陷入局部極小并最終趨于全局最優(yōu)的串行結(jié)構(gòu)的優(yōu)化算法。12.1.3模擬退火算法介紹

模擬退火算法(SimulatedAnnealing，SA)最早由Kirkpatrick等應(yīng)用于組合優(yōu)化領(lǐng)域，它是基于Monte-Carlo迭代求解策略的一種隨機(jī)尋優(yōu)算法，其出發(fā)點(diǎn)是基于物理中固體物質(zhì)的退火過程與一般組合優(yōu)化問題之間的相似性。模擬退火算法從某一較高初溫出發(fā)，伴隨溫度參數(shù)的不斷下降,結(jié)合概率突跳特性在解空間中隨機(jī)尋找目標(biāo)函數(shù)的全局最優(yōu)解，即在局部最優(yōu)解能概率性地跳出并最終趨于全局最優(yōu)。模擬退火算法是一種通用的優(yōu)化算法，理論上算法具有概率的全局優(yōu)化性能,目前已在工程中得到了廣泛應(yīng)用，諸如VLSI、生產(chǎn)調(diào)度、控制工程、機(jī)器學(xué)習(xí)、神經(jīng)網(wǎng)絡(luò)、信號(hào)處理等領(lǐng)域。模擬退火算法是通過賦予搜索過程一種時(shí)變且最終趨于零的概率突跳性，從而可有效避免陷入局部極小并最終趨于全局最優(yōu)的串行結(jié)構(gòu)的優(yōu)化算法。12.1.4模擬退火算法要素

1.狀態(tài)空間與狀態(tài)產(chǎn)生函數(shù)（鄰域函數(shù)）

搜索空間也稱為狀態(tài)空間，它由經(jīng)過編碼的可行解的集合組成。

狀態(tài)產(chǎn)生函數(shù)（鄰域函數(shù)）應(yīng)盡可能保證產(chǎn)生候選解遍布全部解空間。通常由

倆部分組成，即產(chǎn)生候選解產(chǎn)生的概率分布。

候選解一般采用按照某一概率密度函數(shù)對(duì)解空間進(jìn)行隨機(jī)采樣來獲得。

概率分布可以是均勻分布，正太分布，指數(shù)分布等等。2.狀態(tài)轉(zhuǎn)移概率（接受概率）

狀態(tài)轉(zhuǎn)移概率是指從一個(gè)狀態(tài)向另一個(gè)狀態(tài)的轉(zhuǎn)移概率；

通俗的理解是接受一個(gè)新解為當(dāng)前解的概率；

它與當(dāng)前的溫度參數(shù)T有關(guān)，隨溫度參數(shù)T有關(guān)，隨溫度下降而減小。

一般采用Metropolis準(zhǔn)則。3.冷卻進(jìn)度表T(t)

冷卻進(jìn)度表是指從某一高溫狀態(tài)T(t)來表示，則經(jīng)典模擬退火算法的降溫方式為：而快速模擬退火算法的降溫方式為：這倆種方式都能夠使得模擬退火算法收斂于全局最小點(diǎn)。實(shí)驗(yàn)表明，初溫越大，獲得高質(zhì)量解的幾率越大，但花費(fèi)的計(jì)算時(shí)間將增加。因此，初溫的確定應(yīng)折衷考慮優(yōu)化質(zhì)量和優(yōu)化效率，常用方法包括：均勻抽樣一組狀態(tài)，以各狀態(tài)目標(biāo)值的方差初溫。隨機(jī)產(chǎn)生一組狀態(tài)，確定倆狀態(tài)間的最大目標(biāo)值差，然后依據(jù)差值，利用一定的函數(shù)確定初溫。比如，

其中為初始接受概率。利用經(jīng)驗(yàn)公式給出。5.內(nèi)循環(huán)終止準(zhǔn)則

或稱Metropolis抽樣穩(wěn)定準(zhǔn)則，用于決定在各溫度下產(chǎn)生候選解的數(shù)目。常用的抽樣穩(wěn)定準(zhǔn)則包括：

檢驗(yàn)?zāi)繕?biāo)函數(shù)的均值是否穩(wěn)定；

連續(xù)若干步的目標(biāo)值變化較??；

按一定的步數(shù)抽樣。

6.外循環(huán)終止準(zhǔn)則

即算法終止準(zhǔn)則，常用的包括：

設(shè)置終止溫度的閾值；

設(shè)置外循環(huán)迭代次數(shù)；

算法搜索到的最優(yōu)值連續(xù)若干步保持不變；

檢驗(yàn)系統(tǒng)熵是否穩(wěn)定。12.1.5模擬退火算法流程圖

見下圖：12.2.1模擬退火算法基本內(nèi)容初始化：初始溫度T(充分大)，初始解狀態(tài)S(是算法迭代的起點(diǎn))，每個(gè)T值的迭代次數(shù)L

對(duì)k=1，……，L做第(3)至第6步：

產(chǎn)生新解S1

計(jì)算增量

，其中為評(píng)價(jià)函數(shù)。

若則接受s1作為新的當(dāng)前解否則以概率，接受S1作為新的當(dāng)前

解。

如果滿足終止條件則輸出當(dāng)前解作為最優(yōu)解，結(jié)束程序。

終止條件通常取為連續(xù)若干個(gè)新解都沒有被接受時(shí)終止算法。

T逐漸減少，且T->0，然后轉(zhuǎn)第2步。

#include"stdafx.h"

/*MIN[f(x,y)=5sin(xy)+x^2+y^2]*/

#include<math.h>

#include<stdio.h>

#include<string.h>

#include<stdlib.h>

#include<time.h>

/*ObjectiveFunction*/

doubleobjectfunction(doublex,doubley)

{

doublez=0.0;

z=5.0*sin(x*y)+x*x+y*y;

returnz;

}

/*RandomNumberfrom0to1*/

doublernd()

{

doubler;

r=(double)rand()/RAND_MAX;

returnr;

}

intmain()

{

inti;

intmarkovlength=10000;//馬可夫鏈

doubledecayscale=0.95;//衰減參數(shù)

doublestepfactor=0.02;//Metropolis的步長(zhǎng)

doubletemperature=100;//初始溫度

doubletolerance=1e-8;//容差

doubleprex,nextx;

doubleprey,nexty;

doubleprebestx,prebesty;

doublebestx,besty;

doubleacceptpoints=0.0;//Metropolis過程中總接受點(diǎn)

doublexmax,ymax;time_tt;

srand((unsigned)time(&t));

printf("TheMaxNumberXis:\n");

scanf("%lf",&xmax);

printf("TheMaxNumberYis:\n");

scanf("%lf",&ymax);

printf("xmax=%f,ymax=%f\n",xmax,ymax);

prex=-xmax*rnd();

prey=-ymax*rnd();

prebestx=bestx=prex;

prebesty=besty=prey;

do

{

temperature*=decayscale;

acceptpoints=0.0;

for(i=0;i<markovlength;i++)

{

do

{//在此點(diǎn)附近隨機(jī)選下一點(diǎn)

nextx=prex+stepfactor*xmax*(rnd()-0.5);//隨機(jī)數(shù)在正負(fù)0.5之間

nexty=prey+stepfactor*ymax*(rnd()-0.5);

}

while(!(nextx>=-xmax&&nextx<=xmax&&nexty>=-ymax&&nexty<=ymax));

if(objectfunction(bestx,besty)>objectfunction(nextx,nexty))//是否全局最優(yōu)解

{

//保留上一個(gè)最優(yōu)解

prebestx=bestx;

prebesty=besty;

//此為新的最優(yōu)解

bestx=nextx;

besty=nexty;

}

//Metropolis過程

if(objectfunction(prex,prey)>objectfunction(nextx,nexty))

{

//接受,此處lastPoint即下一個(gè)迭代的點(diǎn)以新接受的點(diǎn)開始

prex=nextx;

prey=nexty;

acceptpoints++;

}

else

{

doublechange=(objectfunction(prex,prey)-objectfunction(nextx,nexty))/temperature;

if(exp(change)>rnd())

{

prex=nextx;

prey=nexty;

acceptpoints++;

}

//不接受,保存原解

}

}

}

while(fabs(objectfunction(bestx,besty)-objectfunction(prebestx,prebesty))>tolerance);

printf("abs=%g\n",fabs(objectfunction(bestx,besty)-objectfunction(prebestx,prebesty)));

printf("%f,%f,%f,%f\n",prex,prey,objectfunction(prex,prey),temperature);

printf("TheMinpointsis:%lf,%lf\n",bestx,besty);

printf("TheMindatais:%lf\n",objectfunction(bestx,besty));

getchar();

getchar();

}12.2.2模擬退火的算法描述

模擬退火其實(shí)也是一種貪心算法，但是它的搜索過程引入了隨機(jī)因素。模擬退火算法以一定的概率來接受一個(gè)比當(dāng)前解要差的解，因此有可能會(huì)跳出這個(gè)局部的最優(yōu)解，達(dá)到全局的最優(yōu)解。以圖1為例，模擬退火算法在搜索到局部最優(yōu)解A后，會(huì)以一定的概率接受到E的移動(dòng)。也許經(jīng)過幾次這樣的不是局部最優(yōu)的移動(dòng)后會(huì)到達(dá)D點(diǎn)，于是就跳出了局部最大值A(chǔ)。

圖1：關(guān)于爬山算法與模擬退火，有一個(gè)有趣的比喻：

爬山算法：兔子朝著比現(xiàn)在高的地方跳去。它找到了不遠(yuǎn)處的最高山峰。但是這座山不一定是珠穆朗瑪峰。這就是爬山算法，它不能保證局部最優(yōu)值就是全局最優(yōu)值。模擬退火：兔子喝醉了。它隨機(jī)地跳了很長(zhǎng)時(shí)間。這期間，它可能走向高處，也可能踏入平地。但是，它漸漸清醒了并朝最高方向跳去。這就是模擬退火。在模擬退火算法中應(yīng)注意以下問題：（1）理論上，降溫過程要足夠緩慢，要使得在每一溫度下達(dá)到熱平衡。但在計(jì)算機(jī)實(shí)現(xiàn)中，如果降溫速度過緩，所得到的解的性能會(huì)較為令人滿意，但是算法會(huì)太慢，相對(duì)于簡(jiǎn)單的搜索算法不具有明顯優(yōu)勢(shì)。如果降溫速度過快，很可能最終得不到全局最優(yōu)解。因此使用時(shí)要綜合考慮解的性能和算法速度，在兩者之間采取一種折衷。（2）要確定在每一溫度下狀態(tài)轉(zhuǎn)換的結(jié)束準(zhǔn)則。實(shí)際操作可以考慮當(dāng)連續(xù)m次的轉(zhuǎn)換過程沒有使?fàn)顟B(tài)發(fā)生變化時(shí)結(jié)束該溫度下的狀態(tài)轉(zhuǎn)換。最終溫度的確定可以提前定為一個(gè)較小的值eT，或連續(xù)幾個(gè)溫度下轉(zhuǎn)換過程沒有使?fàn)顟B(tài)發(fā)生變化算法就結(jié)束。

（3）選擇初始溫度和確定某個(gè)可行解的鄰域的方法也要恰當(dāng)。

模擬退火算法的參數(shù)控制問題：模擬退火算法的應(yīng)用很廣泛，可以求解NP完全問題，但其參數(shù)難以控制，其主要問題有以下三點(diǎn)：(1)溫度T的初始值設(shè)置問題。溫度T的初始值設(shè)置是影響模擬退火算法全局搜索性能的重要因素之一、初始溫度高，則搜索到全局最優(yōu)解的可能性大，但因此要花費(fèi)大量的計(jì)算時(shí)間；反之，則可節(jié)約計(jì)算時(shí)間，但全局搜索性能可能受到影響。實(shí)際應(yīng)用過程中，初始溫度一般需要依據(jù)實(shí)驗(yàn)結(jié)果進(jìn)行若干次調(diào)整。(2)退火速度問題:模擬退火算法的全局搜索性能也與退火速度密切相關(guān)。一般來說，同一溫度下的“充分”搜索(退火)是相當(dāng)必要的，但這需要計(jì)算時(shí)間。實(shí)際應(yīng)用中，要針對(duì)具體問題的性質(zhì)和特征設(shè)置合理的退火平衡條件。

(3)溫度管理問題:溫度管理問題也是模擬退火算法難以處理的問題之一。實(shí)際應(yīng)用中，由于必須考慮計(jì)算復(fù)雜度的切實(shí)可行性等問題，常采用如下所示的降溫方式：

式中k為正的略小于1.00的常數(shù)，t為降溫的次數(shù)。

main.mzuobiao=[0.370.750.450.760.710.070.420.590.320.60.30.670.620.670.20...0.350.270.940.820.370.610.420.60.390.530.40.630.50.980.68;0.910.870.850.750.720.740.710.690.640.640.590.590.550.550.5...0.450.430.420.380.270.260.250.230.190.190.130.080.040.020.85]plot(zuobiao(1,),zuobiao(2,),'g'),holdonplot(zuobiao(1,),zuobiao(2,))length=max(size(zuobiao));%求初始距離..zhixu=randperm(length)%隨機(jī)生成一個(gè)路線經(jīng)過點(diǎn)的順序temp=zuobiao(1,);newzuobiao(1,)=temp(zhixu);temp=zuobiao(2,);newzuobiao(2,)=temp(zhixu);newzuobiaof=juli(newzuobiao)%參數(shù)定義區(qū)--------------------------------------%初始溫度為10000tmax=100;t