模擬退火算法原理及matlab源代碼

1、模擬退火算法模擬退火算法是一種通用的隨機(jī)搜索算法，是局部搜索算法的擴(kuò)展。它的思想是再1953年由metropolis提出來(lái)的，到1983年由kirkpatrick等人成功地應(yīng)用在組合優(yōu)化問(wèn)題中。模擬退火算法來(lái)源于固體退火原理，將固體加溫至充分高，再讓其徐徐冷卻，加溫時(shí)，固體內(nèi)部粒子隨溫升變?yōu)闊o(wú)序狀，內(nèi)能增大，而徐徐冷卻時(shí)粒子漸趨有序，在每個(gè)溫度都達(dá)到平衡態(tài)，最后在常溫時(shí)達(dá)到基態(tài)，內(nèi)能減為最小。根據(jù)Metropolis準(zhǔn)則，粒子在溫度T時(shí)趨于平衡的概率為e-AE/(kT),其中E為溫度T時(shí)的內(nèi)能，AE為其改變量，k為Boltzmann常數(shù)。用固體退火模擬組合優(yōu)化問(wèn)題，將內(nèi)能E模擬為目標(biāo)函數(shù)值f,

2、溫度T演化成控制參數(shù)t,即得到解組合優(yōu)化問(wèn)題的模擬退火算法：由初始解i和控制參數(shù)初值t開(kāi)始，對(duì)當(dāng)前解重復(fù)“產(chǎn)生新解-計(jì)算目標(biāo)函數(shù)差一接受或舍棄”的迭代，并逐步衰減t值，算法終止時(shí)的當(dāng)前解即為所得近似最優(yōu)解，這是基于蒙特卡羅迭代求解法的一種啟發(fā)式隨機(jī)搜索過(guò)程。退火過(guò)程由冷卻進(jìn)度表(CoolingSchedule)控制，包括控制參數(shù)的初值t及其衰減因子At、每個(gè)t值時(shí)的迭代次數(shù)L和停止條件So模擬退火算法新解的產(chǎn)生和接受可分為如下四個(gè)步驟：第一步是由一個(gè)產(chǎn)生函數(shù)從當(dāng)前解產(chǎn)生一個(gè)位于解空間的新解；為便于后續(xù)的計(jì)算和接受，減少算法耗時(shí)，通常選擇由當(dāng)前新解經(jīng)過(guò)簡(jiǎn)單地變換即可產(chǎn)生新解的方法，如對(duì)構(gòu)成新解的

3、全部或部分元素進(jìn)行置換、互換等，注意到產(chǎn)生新解的變換方法決定了當(dāng)前新解的鄰域結(jié)構(gòu)，因而對(duì)冷卻進(jìn)度表的選取有一定的影響。第二步是計(jì)算與新解所對(duì)應(yīng)的目標(biāo)函數(shù)差。因?yàn)槟繕?biāo)函數(shù)差僅由變換部分產(chǎn)生，所以目標(biāo)函數(shù)差的計(jì)算最好按增量計(jì)算。事實(shí)表明，對(duì)大多數(shù)應(yīng)用而言，這是計(jì)算目標(biāo)函數(shù)差的最快方法。第三步是判斷新解是否被接受，判斷的依據(jù)是一個(gè)接受準(zhǔn)則，最常用的接受準(zhǔn)則是Metropo1is準(zhǔn)則：若At0則接受S作為新白當(dāng)前解S,否則以概率exp(-“/T)接受S作為新白當(dāng)前解So第四步是當(dāng)新解被確定接受時(shí)，用新解代替當(dāng)前解，這只需將當(dāng)前解中對(duì)應(yīng)于產(chǎn)生新解時(shí)的變換部分予以實(shí)現(xiàn)，同時(shí)修正目標(biāo)函數(shù)值即可。此時(shí)，當(dāng)前解

4、實(shí)現(xiàn)了一次迭代。可在此基礎(chǔ)上開(kāi)始下一輪試驗(yàn)。而當(dāng)新解被判定為舍棄時(shí)，則在原當(dāng)前解的基礎(chǔ)上繼續(xù)下一輪試驗(yàn)。模擬退火算法與初始值無(wú)關(guān)，算法求得的解與初始解狀態(tài)S(是算法迭代的起點(diǎn))無(wú)關(guān)；模擬退火算法具有漸近收斂性，已在理論上被證明是一種以概率l收斂于全局最優(yōu)解的全局優(yōu)化算法；模擬退火算法具有并行性。模擬退火算法的步驟：(1)初始化：初始溫度T(充分大)，初始解狀態(tài)S(是算法迭代的起點(diǎn))，每個(gè)T值的迭代次數(shù)L(2)對(duì)k=1,，L做第(3)至第6步：(3)產(chǎn)生新解S(4)計(jì)算增量久=C(S)-C(S),其中C(S)為評(píng)價(jià)函數(shù)(5)若At0,然后轉(zhuǎn)第2步。退火算法解非線性方程組Matlab程序clear

5、,clc應(yīng)是退火算法的主程序，它需要調(diào)用的函數(shù)是%數(shù)(1)nonLinearSumError1:計(jì)算非線性方程組總誤差的函數(shù)%數(shù)(2)newSolver1:在一組解的鄰域產(chǎn)生另一組解%數(shù)(3)isSolution:驗(yàn)證方程是否得解煩置初始值i=0;T=10001;j=0;%i:同一溫度下?tīng)顟B(tài)轉(zhuǎn)移次數(shù);T:溫度；j:下降溫度precision=0.1;x1Group=1;%x1Group:可能解的組數(shù)x1N=4;%|營(yíng)線性方程組的元數(shù)x1=round(-0.5+rand(x1Group,x1N)*20);%隨機(jī)生成-1010之間的初解errorHold=Inf;xHold=0;%x1=-751-

6、3;i=0;whilei200i=i+1;j=0；T=T-50;%退火whilej200j=j+1;functionError1=nonLinearSumError1(x1);%計(jì)算x1的誤差x2=newSolver1(x1,functionError1,-10,1,10);%在x1的鄰域生成新一組解x2functionError2=nonLinearSumError1(x2);%計(jì)算x2的誤差%僉查方程是否得解solution1,minError1,isTrue1=isSolution(x1,functionError1,precision);solution2,minError2,isTr

7、ue2=isSolution(x2,functionError2,precision);ifisTrue1=1方程得解functionError1solutiourni,jreturnelseifisTrue2=1方程得解solution2functionError2i,jreturnend%x1%x2iffunctionError2-functionError10x1=x2;%x2比x1好，用x2取代x1elseiferrorHold-functionError20%x1=xHold;elsep_x2x1=exp(-log(functionError2-functionError1)/T);

8、屈態(tài)轉(zhuǎn)移概率，注意：誤差取對(duì)數(shù)，因?yàn)橐獾姆蔷€性方程組比較復(fù)雜，阿能解的一點(diǎn)偏差會(huì)引起方程很大的變化。所以通過(guò)取對(duì)數(shù)縮小差距。ifrand(1)p_x2x1%狀態(tài)轉(zhuǎn)移xHold=x1;%hHold:把比較好的解保留下來(lái)errorHold=functionError1;%比較好的解對(duì)應(yīng)的誤差x1=x2;endendendendsolution1functionError1solution2functionError2函數(shù)(1):計(jì)算待解方程的絕對(duì)總誤差functionfuntionError=nonLinearSumError1(X)%方程的解是-7,5,1,-3funtionError=.ab

9、s(X(:,1).A2-sin(X(:,2),A3)+X(:,3).A2-exp(X(:,4)-50.566253390821)+.abs(X(:,1)A3+X(:,2)A2-X(:,4)A2+327)+.abs(cos(X(:,1)A4)+X(:,2)A4-X(:,3)A3-624.679868769613)+.abs(X(:,1)A4-X(:,2)A3+2.AX(:,3)-X(:,4)A4-2197)；函數(shù)(2):在x1的領(lǐng)域產(chǎn)生一組新的解%newSolver1根據(jù)x1的誤差給出一個(gè)新的可能解xfunctionx2=newSolver1(x1,x1Error,leftBound,dista

10、nce,rightBound)%parameter=leftBound,distance,rightBound%leftBound:解空間的左邊界，distance:可能解的間隔，rightBound:解空間的右邊界%軍空間是指在一個(gè)坐標(biāo)軸上解的左右邊界和解之間的間隔x1Group,x1N=size(x1);%x1Group:x1的行數(shù)，x1N:方程的元數(shù)%round(-0.5+rand(x1Group,x1N)*2)ifx1Error=30%在解空間上移動(dòng)1格x2=x1+round(-0.5+rand(x1Group,x1N)*2)*distance;k=x2rightBound;%防止新解

11、越過(guò)右邊界x2(:,k)=rightBound;elseifx1Error30&x1Error=100%在解空間上移動(dòng)3格以下x2=x1+round(-0.5+rand(x1Group,x1N)*6)*distance;k=x2rightBound;x2(:,k)=rightBound;elseifx1Error100&x1Error=1000%在解空間上移動(dòng)9格以下x2=x1+round(-0.5+rand(x1Group,x1N)*20)*distance;k=x2rightBound;x2(:,k)=rightBound;elseifx1Error1000&x1Error=10000%在解空間上移動(dòng)20格以下x2=x1+round(-0.5+rand(x1Group,x1N)*40)*distance;k=x2rightBound;x2(:,k)=rightBound;elseifx1Error10000%在解空間上移動(dòng)30格以下x2=x1+round(-0.5+rand(x1Group,x1N)*60)*distance;k=x2rightBound;x2(:,k)=rightBound;endifx1=x2x2=round(-0.5+rand(x1Group,x1N)*20);end函數(shù)(3):喈ij斷方程