用MATLAB實(shí)現(xiàn)模擬退火算法

模

擬

退

火

算

法

及

其

MATLAB實(shí)現(xiàn)模擬退火6.1算法基

本

理

論6.2算

法

的MATLAB

實(shí)

現(xiàn)6.3應(yīng)用實(shí)例模

擬

退

火

算

法

及其MATLAB第

6

章實(shí)

現(xiàn)簡

單

了

解

退

火

算

法

特

點(diǎn)介紹

模

擬退

火

前，

先

介

紹

爬山

算

法

。爬山算法是一種簡單的貪心搜索算法，

該算法每次從當(dāng)前解的臨近解空間中選擇一個(gè)最優(yōu)解作為當(dāng)前解，

直到達(dá)到

一

個(gè)

局

部

最

優(yōu)

解

。簡

單

了

解

退

火

算

法

特

點(diǎn)爬

山

算

法如圖所示：

假設(shè)C

點(diǎn)為當(dāng)前解，爬山算法搜索到A

點(diǎn)這個(gè)局部最優(yōu)解就會(huì)停止搜索，

因?yàn)樵贏

點(diǎn)無

論向那個(gè)方向小幅度移動(dòng)都不能得到更優(yōu)的解。模擬退火

算

法在搜索到局部最優(yōu)解A

后，會(huì)以一定的概率接受到E

的移動(dòng)。也許經(jīng)過幾次這樣的不是局部最優(yōu)的移動(dòng)后會(huì)

到

達(dá)D點(diǎn)，于是就跳

出

了局部最大值A(chǔ)。6.1

算法基本理論一、算法概述工程中許多實(shí)際優(yōu)化問題的目標(biāo)函數(shù)都是非凸的，

存在許多局部最優(yōu)解。求解全

局

優(yōu)

化問題的方法可

分

為

兩

類

：確定性方法和隨機(jī)性方法。確定性算法適用于求解具有一些特殊特征的問題，

而梯度法和一般的隨機(jī)搜索方法則沿著目標(biāo)函數(shù)下降方

向搜索，

因此常常陷入局部而非全局最優(yōu)解

。6.1

算法基本理論一

、

算

法

概

述模擬退火算法(

SA)

是一種通用概率算法。用來

在一個(gè)大的搜索空間內(nèi)尋找問題的最優(yōu)解。1953年

，Metropolis等提出了模擬退火的思想。1983

年，Kirkpatrick等將SA引入組合優(yōu)化領(lǐng)域。6.1

算法基本理論二

、

基

本

思

想退火

，俗稱

固體降

溫先把固體加熱至足夠高溫，使固體中所有粒子處于無序的狀態(tài)，

然后將溫度緩慢下降，

粒子漸漸有序，

這樣只要溫度上升得足夠高，冷卻過程足夠慢，則所

有粒子

最

終

會(huì)處于

最

低能

態(tài)

。算法試圖隨著控制參數(shù)T的降低，

使目標(biāo)函

數(shù)

值f

(

內(nèi)

能E)

也逐漸降低，

直至趨于全局最小

值

(

退

火

中

低

溫

時(shí)

的

最

低

能

量

狀

態(tài)

)

,

算

法

工

作

過

程

就

像

固

體

退

火

過

程

一

樣

。模擬退火退火解粒子狀態(tài)最優(yōu)解能

量

最

低的

狀

態(tài)目

標(biāo)

函

數(shù)

f內(nèi)能控制參數(shù)溫度T6.1

算法基本理論模擬退火算法的由來6.1

算法基本理論Metropolis準(zhǔn)

則—

以

概

率

接

受

新

狀

態(tài)若在溫度T,當(dāng)前狀態(tài)i(解1)

→新狀態(tài)

j(解2)隨機(jī)數(shù)，則仍接受狀態(tài)j

為當(dāng)前狀態(tài)；若不

成

立，則保

留

狀

態(tài)I

為當(dāng)前狀態(tài)

。若E;(目標(biāo)函數(shù)f?)<Ei(fi),若E;>E?,

概率則接受j為當(dāng)前狀態(tài)；大于(0,1)區(qū)間的6.1

算法基本理論Ej>Ei(

更差的解)時(shí)，O<P<1,P隨著T

的減小而減??；Ej-EiKTP

二

e當(dāng)前狀態(tài)的內(nèi)能新狀態(tài)的內(nèi)能溫度在高溫下，

可接受與當(dāng)前狀態(tài)能量差

較

大的新

狀

態(tài)；

在低溫下，

只接受與當(dāng)前狀態(tài)能量差較小的新狀態(tài)。當(dāng)初始溫度足

夠高時(shí)，

概

率P

接近于

1,所

以

當(dāng)前

解經(jīng)過擾動(dòng)產(chǎn)生的新解，

無論好壞，

基本都可以被接受為當(dāng)前解。即不受制于當(dāng)前解，

不會(huì)困在局部最優(yōu)解中，

可以遍及解空間的各個(gè)區(qū)域，

當(dāng)然也不會(huì)保持在最優(yōu)解處。隨著溫度降低，概率降低，

較差解被接受的次數(shù)減少，

當(dāng)前解逐漸停留到最優(yōu)解周圍。溫度達(dá)到終止溫度前，

概率足夠低，使得只有最優(yōu)解被接受，較差解都不接受。最優(yōu)解即為最后接受的當(dāng)前解。算

法

總

結(jié)6.1

算法基本理論6.1

算法基本理論三、算法其他參數(shù)的說明退火過程由一組初始參數(shù)，

即冷卻進(jìn)度表控

制

，它的核心是盡量使系統(tǒng)達(dá)到準(zhǔn)平衡，

以使算法在有限

的時(shí)間內(nèi)逼近最優(yōu)解。冷卻進(jìn)度表包括：1.控制參數(shù)的初值To:冷卻開始的溫度；2.控制參數(shù)T的衰減函數(shù)：

要將連續(xù)的降溫過程，

離散成降溫過程中的一系列溫度點(diǎn)，

衰減函數(shù)即計(jì)算這一系列溫度的表達(dá)式；3.控制參數(shù)T

的終值Tr(停止準(zhǔn)則);4.Markov

鏈的長度Lx:任意溫度T

的迭代次數(shù)。6.1

算法基本理論四、算法基本步驟1、令T

=To,隨機(jī)生成一個(gè)初始解xo,

并計(jì)算相應(yīng)的目標(biāo)函數(shù)值E(xo)

;2、

令T等于冷卻進(jìn)度表中的下一個(gè)值T:;3、

根據(jù)當(dāng)前解xi進(jìn)行擾動(dòng)，產(chǎn)生一個(gè)新解xj,計(jì)相應(yīng)

的目標(biāo)函數(shù)值E(xj),得到

△e=E(xj)-E(xi);4、△e

<0,

則新解x;被接受，作為新的當(dāng)前解；

△e

>0,則新解x;按概率接

受

；5、

在溫度T

下，重復(fù)Lg次的擾動(dòng)和接受過程，即執(zhí)行步

驟

(

3

)

、(4);6、

判斷T

是否已達(dá)到T,

是，則終止算法；否則轉(zhuǎn)到(2

)

繼

續(xù)

執(zhí)

行△y

>o平

Y計(jì)算概

率

與[o,1]隨

機(jī)

數(shù)之間的差

值差

值

大于ON

結(jié)束

，輸

出

當(dāng)前

解Tk+i

=aTk初

始

溫

度，

隨機(jī)

產(chǎn)

生

初

始

解

。T>Tr當(dāng)前解x擾動(dòng)

次數(shù)>Lp業(yè)

N擾動(dòng)，產(chǎn)生新解xi+1計(jì)算兩解的目標(biāo)函數(shù)差值△y接受

新

解

作

為

當(dāng)前

解NYYN6.1

算法基本理論四、算法基本步驟算法實(shí)質(zhì)分為兩層循環(huán)，

在任一溫度下隨機(jī)擾動(dòng)產(chǎn)生新解

，

計(jì)算

目標(biāo)函數(shù)值的變化，

決定

是

否接

受

。

由

于

算

法初始溫度比較高，這樣使E增大的新解在初始時(shí)也可能被

接受，因此能跳出局部極小值，

然后通過緩慢地降低溫度，

算法可能收斂到全局最優(yōu)解。雖然在低溫時(shí)接受函數(shù)已經(jīng)非常小了，但仍不排除有接受更差解得可能，

因此一般都會(huì)把退火過程中碰到的最好的可行解(歷史最優(yōu)解)也記錄下來，與終止算法前最后被

接受

解一并

輸出

。6.1

算法基本理論五

、

幾

點(diǎn)

說

明1、

新解

的

產(chǎn)

生要求盡可能地遍及解空間的各個(gè)區(qū)域，這樣，在某一

恒定溫度下，

不斷產(chǎn)生新解時(shí)，

就可能跳出局部最優(yōu)解。

2、

收

斂的一

般條

件

：·

初始溫

度

足

夠高

；·

熱平

衡時(shí)間足夠長

；●終止溫度

足

夠

低

；●

降溫

過程足

夠

緩

慢

；6.1

算法基本理論五

、

幾

點(diǎn)

說

明3

、參數(shù)的選

擇

：(1)初始值ToT?

越大越好，但為了減少計(jì)算量，要根據(jù)實(shí)際情況選擇；(2)控制參數(shù)T的衰減函數(shù)常用Te+1=aTx,a的取值范圍：

0.5~0.99;(3)Mark

ov鏈長度Lk=100n,n為問題規(guī)模。6.1

算法基本理論六

、

算

法

優(yōu)

缺

點(diǎn)優(yōu)點(diǎn)

：計(jì)算過程簡單，

通用，

魯棒性強(qiáng)，

適用于并行處理，可用于求解復(fù)雜的非線性優(yōu)化問題。缺點(diǎn)

：收斂速度慢，

執(zhí)行時(shí)間長，算法性能與初始值有關(guān)及參數(shù)敏感等缺點(diǎn)。6.

2算法的MATLAB

實(shí)現(xiàn)旅行商問題一

名商人要到n個(gè)不同的城市去推銷商品，每2個(gè)城市/和j

之間的距離為d,

如何選擇

一

條

路徑使得商人每

個(gè)

城市走

一

遍

后回到起點(diǎn)所走

的路徑最短。例

：有

5

2

座城市，

已知

每

座

城

市的

坐

標(biāo)，

求

每

個(gè)

城

市

走

一

遍

后回

到

起點(diǎn)，

所

走的

路

徑

最

短

。N結(jié)束，輸出當(dāng)前解Tk+1=0.99Tk猶動(dòng)次數(shù)>10000r

下擾動(dòng)，產(chǎn)生新解Xi+1△y

>oY計(jì)算概率與[o,1]

隨機(jī)

數(shù)之間的差值計(jì)算兩解的目標(biāo)函數(shù)差值A(chǔ)y(兩條路徑之差)初始溫度(93),隨

機(jī)產(chǎn)生初始解(1到

52的隨機(jī)排列)。隨機(jī)產(chǎn)生

0～1的數(shù)亞數(shù)

>

0

.

5二變換法

三變換法接受新解作為

當(dāng)前解T

>3當(dāng)前解x;差值大于ONY擾動(dòng)：YNYN1.TS

P

問題的解空間和初始解TSP問題的解空間s是遍訪每個(gè)城市恰好一次的所有回

路，

是所有城市的排列的集合。

即S={(c?,c?,…,cn)|(c?,c?,…,cn)為{1,2,…,n}的排列}Si:遍訪n

個(gè)城市

的一條

路

徑；Ci=j:

第i次訪問城市j。因?yàn)樽顑?yōu)解和初始狀態(tài)沒有強(qiáng)的依賴關(guān)系，所以

sO={1,2,…,n

}的隨機(jī)排列6.

2算法的MATLAB

實(shí)現(xiàn)一、算法設(shè)計(jì)步驟2.

新

解的產(chǎn)

生(擾動(dòng)

)(1)二變換法：任選序號u,v

(

設(shè)u<v<n),的訪問順序。(

2)三變換法：任選序號u,v,w

(

設(shè)u≤v<w),的路徑

插

到w

之

后

訪

問6.

2算法的MATLAB

實(shí)現(xiàn)一、算法設(shè)計(jì)步驟交

換u,v

之間將u,v

之

間else%否則，三變換i

ii

1

=

;

==ind3)

…i

d1(i

il(

i;nd1-ind2)==1)

p1=ind1;tmp2

=ind2;tmp3

=ind3;mndtea3dnran==inl102)d3dni0ndnd2ile=0h1while

t>=tffor

r

1

ngth%隨機(jī)產(chǎn)生0~1的數(shù)，若小于0.5,

則二變換i

d1

==i

)

_p

_d

(ind2);sol_new(ind2)=tmp1;w;n1olinn_d1l

nwse=n1olmndsted20;ind2hile1=0;iwnd5)ledrknaraMif=6.

2算法的MATLAB

實(shí)現(xiàn)一、算法設(shè)計(jì)步驟ind2

=ceil(rand.*amount);ind3

=ceil(rand.*amount);ind1=ceil(rand.*amount);ind2

=ceil(rand.*amount);6.

2算法的MATLAB

實(shí)現(xiàn)一、算法設(shè)計(jì)步驟if(ind1<ind2)&&(ind2<ind3)end%ind1<ind2<ind3tmplist1=sol_new((ind1+1):(ind2-1));%u

、v之間的城市sols_ol

((

):1(

)d);3-ind2+

)將:i

城市移到u后面tmplist1;

%u

、v之間的城市移到w后面end2%n3+iindnd(iinwwenen_elseif

(ind1<ind3)&&(ind3<ind2)

ind2=tmp3;ind3=tmp2;elseif

(ind2<ind1)&&(ind1<ind3)

ind1=tmp2;ind2=tmp1;sol_new((ind1+1):(ind1+ind3-ind2+1))=

…elseif

(ind2<ind3)&&(ind3<ind1)ind1=tmp2;ind2

=tmp3;ind3

=tmp1;elseif

(ind3<ind1)&&(ind1<ind2)ind1=tmp3;ind2

=tmp1;ind3

=tmp2;elseif

(ind3<ind2)&&(ind2<ind1)ind1=tmp3;ind2=tmp2;ind3

=tmp1;3.目標(biāo)

函

數(shù)訪問所有城市的路徑總長度：模擬退火算法求出目標(biāo)函數(shù)的最小值一、算法設(shè)計(jì)步驟6.

2算法的MATLAB

實(shí)現(xiàn)6.

2算法的MATLAB

實(shí)現(xiàn)一、算法設(shè)計(jì)步驟%計(jì)算目標(biāo)函數(shù)即內(nèi)能

1

-1)E_new=E_new+

…end%從第一個(gè)城市到最后一個(gè)城市的距離E_new

=E_new+

…mount;a0inefo_rEdist_matrix(sol_new(amount),sol_new(1));dist_matrix(sol_new(i),sol_new(i+1));4.

目

標(biāo)函數(shù)差計(jì)算變換前的解和變換后目標(biāo)函數(shù)的差值△c1=c(s')-c(s)5.Metropolis

接

受準(zhǔn)則以新解與當(dāng)前解的目標(biāo)函數(shù)差定義接受概率，

即一、算法設(shè)計(jì)步驟6.

2算法的MATLAB實(shí)現(xiàn)6.

2算法的MATLAB

實(shí)