高嘉 任亞明

【關(guān)鍵詞】組合優(yōu)化問題;模擬退火;分支界定

【中圖分類號(hào)】O221.4 【文獻(xiàn)標(biāo)識(shí)碼】A 【文章編號(hào)】1674-0688（2021）05-0066-03

0 引言

在優(yōu)化領(lǐng)域中，根據(jù)變量性質(zhì)的不同大體可以分為兩類：一類是包含連續(xù)變量的優(yōu)化問題;另一類是包含整數(shù)變量的優(yōu)化問題（也可稱之為組合優(yōu)化問題）。組合優(yōu)化問題的目標(biāo)是從組合問題的可行解集中求出最優(yōu)解，組合優(yōu)化往往涉及排序、分類等問題，它是優(yōu)化領(lǐng)域的一個(gè)重要分支。

在求解組合優(yōu)化問題中，人們首先想到的是取整的方法，即互聯(lián)組合優(yōu)化問題變量必須為整數(shù)的約束條件，按照連續(xù)變量的優(yōu)化問題對(duì)其進(jìn)行求解，對(duì)于得到的結(jié)果按照某種方法取整。該方法簡(jiǎn)單，但是所得到的結(jié)果往往會(huì)違背優(yōu)化問題的約束條件或者得到次優(yōu)的結(jié)果。在取整的基礎(chǔ)上，分支界定方法被提出，分支界定方法的本質(zhì)也是按照連續(xù)變量的優(yōu)化問題對(duì)其進(jìn)行求解，不過在每次得到的結(jié)果不滿足整數(shù)約束時(shí)，并不是進(jìn)行簡(jiǎn)單的取整操作，而是通過縮小變量可行域的方法，逐步逼近問題的最優(yōu)解。分支界定方法可以有效處理組合優(yōu)化問題，但是其每次縮小可行域就要進(jìn)行一次連續(xù)優(yōu)化問題的求解[1]，因此計(jì)算量大。同時(shí)，由于采取連續(xù)變量的優(yōu)化問題對(duì)其進(jìn)行求解，對(duì)于所求解的數(shù)學(xué)問題有這嚴(yán)格的數(shù)學(xué)要求，例如函數(shù)必須連續(xù)且必須為凸，這樣才可以保證所求的解為問題的全局最優(yōu)解，如果是多極值問題，無法保證結(jié)果為問題的全局最優(yōu)解，解的狀態(tài)與問題的初始值密切相關(guān)，而初始值一般是隨機(jī)給定，因此分支界定方法的使用受到較多的約束與限制。

另一種求解組合優(yōu)化問題的方法為智能化方法。常見的智能化方法有粒子群算法[2]、遺傳算法[3]、模擬退火算法[4]等。粒子群算法根據(jù)其迭代規(guī)則，如果變量為整數(shù)變量，需要在其迭代規(guī)則的基礎(chǔ)上進(jìn)行進(jìn)一步討論。遺傳算法和模擬退火算法其初始變量隨機(jī)生成，可以直接轉(zhuǎn)化為相關(guān)的整數(shù)，不需要做額外的變化。也就是說，遺傳算法和模擬退火算法都可以求解組合優(yōu)化問題。遺傳算法與模擬退火算法相對(duì)比，遺傳算法需要選擇、交叉、變異，而模擬退火算法僅需要生成新的解，采取Metropolis準(zhǔn)則與原解進(jìn)行比較并進(jìn)行相應(yīng)的跟新操作。模擬退火算法迭代過程簡(jiǎn)單，容易實(shí)現(xiàn)，因此本文采取模擬退火算法求解組合優(yōu)化問題。

1 模擬退火算法的基本原理及編程實(shí)現(xiàn)

模擬退火算法通過觀察自然科學(xué)中固體退火過程類推而來，最早的思想是由N. Metropolis等人于1953年提出[5]。我們首先給出模擬退火算法的具體迭代步驟，然后逐條對(duì)其進(jìn)行相應(yīng)的解釋。

模擬退火算法的具體迭代步驟如下。

（1）設(shè)置模擬退火算法初始溫度T，降溫速率q，結(jié)束溫度Tend，Metropolis準(zhǔn)則鏈長(zhǎng)L。

（2）隨機(jī)生成組合優(yōu)化問題的隨機(jī)初始解S1。

（3）將S1帶入組合優(yōu)化問題的目標(biāo)函數(shù)，求解目標(biāo)函數(shù)值f? ?1。

（4）設(shè)置K=0。

（5）隨機(jī)生成新的組合優(yōu)化問題的解S2;將S2帶入組合優(yōu)化問題的目標(biāo)函數(shù)，求解目標(biāo)函數(shù)值f? ? 2。

（6）Metropolis準(zhǔn)則，比較f? ?1和f? ?2更新S1。

（7）更新k=k+1;假如k

（8）更新溫度T=T*q，假如T<Tend，則停止迭代輸出結(jié)果;否則轉(zhuǎn)步驟（4）。

1.1 隨機(jī)生成組合優(yōu)化問題的解的隨機(jī)生成

組合優(yōu)化問題的隨機(jī)解生成代碼如下所示。

S1=round（LB+（UB-LB）*rand（1，N））? （1）

公式（1）中，S1為隨機(jī)生成組合優(yōu)化問題的隨機(jī)解，LB和UB分別代表與預(yù)先給定的變量的下限和上限向量，rand為Matlab函數(shù)，表示生成0～1的均勻分布隨機(jī)數(shù)。round（）為Matlab函數(shù)，表示按照指定的小數(shù)位數(shù)進(jìn)行四舍五入運(yùn)算的結(jié)果，通過round操作可以保證組合優(yōu)化問題的初始解為整數(shù)，而且是在給出的上下限范圍之內(nèi)。

1.2 模擬退火算法的Metropolis準(zhǔn)則

Metropolis準(zhǔn)則的制定在于以一定的概率跳出當(dāng)前的最優(yōu)解，即以一定的概率接收新的解，即使該新解目標(biāo)函數(shù)評(píng)價(jià)低于已經(jīng)存在的最優(yōu)解，其目的是跳出局部最優(yōu)解，使得算法具有全局搜索的能力，這也是模擬退火算法最為核心的概念[6]。

S1表示當(dāng)前解，其對(duì)應(yīng)的目標(biāo)函數(shù)值為f? ?（S1）;S2表示新生成的解，其對(duì)應(yīng)的目標(biāo)函數(shù)值為f? ?（S2）。組合優(yōu)化問題的解由S1變?yōu)镾2的接受概率P：

P=1，f （S2）f （S1） （2）

根據(jù)公式（2）的描述可知，當(dāng)f? ?（S2）f? （S1）時(shí)模擬退火算法不會(huì)立刻將新生成的解S2拋棄，其實(shí)進(jìn)行相應(yīng)的概率操作，首先利用均勻分布函數(shù)生成一個(gè)0～1的服從均勻分布概率的隨機(jī)數(shù)，然后將這個(gè)隨機(jī)數(shù)和e^（-（f （S2）-f （S1））/T）相比較，假如生成的隨機(jī)數(shù)小于e^（-（f （S2）-f （S1））/T）則接受組合優(yōu)化問題的解由S1變?yōu)镾2，否則拋棄生成的新解S2。具體的偽代碼如下所示。

隨機(jī)生成新解S2;

? S1評(píng)價(jià)指標(biāo)=以S1為變量帶入目標(biāo)函數(shù)，返回的目標(biāo)函數(shù)值;

? S2評(píng)價(jià)指標(biāo)=以S2為變量帶入目標(biāo)函數(shù)，返回的? 目標(biāo)函數(shù)值;

評(píng)價(jià)指標(biāo)差=S2評(píng)價(jià)指標(biāo)-S1評(píng)價(jià)指標(biāo)

if評(píng)價(jià)指標(biāo)差<0

S1=S2;

S1評(píng)價(jià)指標(biāo)=S2評(píng)價(jià)指標(biāo);

else

if exp（-dC/溫度）>0-1之間的均勻隨機(jī)數(shù)

S1=S2;

S1評(píng)價(jià)指標(biāo)=S2評(píng)價(jià)指標(biāo);

end

end

通過對(duì)Metropolis準(zhǔn)則的分析我們發(fā)現(xiàn)，只有當(dāng)f （S2）>f （S1）時(shí)，Metropolis準(zhǔn)則才會(huì)以e^（-（f （S2）-f （S1））/T）的概率接受新解，此時(shí)-（f （S2）-f （S1））/T的值一定為負(fù)，因此e^（-（f （S2）-f （S1））/T）的值一定是0～1的開區(qū)間。同時(shí)我們發(fā)現(xiàn)e^（-（f （S2）-f （S1））/T）的大小還與當(dāng)前時(shí)刻溫度T有關(guān)，根據(jù)模擬退火算法的具體迭代步驟可知，隨著迭代次數(shù)的增加，溫度T是逐漸變小，我們假設(shè)f （S2）-f （S1）大小不變的情況下，隨著溫度T變小e^（-（f （S2）-f （S1））/T）的值也在不斷變小趨近于0。因此，在迭代的初期，模擬退火算法以較高的概率跳出當(dāng)前最優(yōu)解，其目的是保證算法的全局收斂能力，在迭代的末期，模擬退火算法以較低的概率跳出當(dāng)前最優(yōu)解，其目的是在最優(yōu)解附件搜索提高算法結(jié)果的精度。

2 仿真

2.1 仿真算例

在這一節(jié)中，我們采取的算例是來自2019年“電工杯數(shù)學(xué)建模大賽”A題，電源規(guī)劃第二問。具體的算例描述如下：IEEE-RTS單階段電源擴(kuò)展規(guī)劃（目標(biāo)函數(shù)只考慮機(jī)組投資費(fèi)用），IEEE-RTS系統(tǒng)由32臺(tái)發(fā)電機(jī)組構(gòu)成，總裝機(jī)容量為3 405 MW，峰值負(fù)荷為2 850 MW。以2019年為基準(zhǔn)年，假設(shè)2030年系統(tǒng)峰值負(fù)荷增長(zhǎng)30%，規(guī)劃2030年當(dāng)年增裝機(jī)組的類型和數(shù)量（見表1）。

根據(jù)問題的描述，我們建立如下的數(shù)學(xué)模型：

minF=x1*220+x2*90+x3*60+x4*50

s.t. 5≥x1≥0，x1≥為整數(shù)

11≥x2≥0，x2≥為整數(shù) （3）

17≥x3≥0，x3≥為整數(shù)

21≥x4≥0，x4≥為整數(shù)

x1*250+x2*100+x3*65+x4*50+3 405≥4 446

其中，變量x1、x2、x3和x4分別代表表2中機(jī)組類型1、2、3和4的新建機(jī)組數(shù)量。當(dāng)然根據(jù)工程建設(shè)的要求，我們需要以最小的工程造價(jià)為目標(biāo)函數(shù)。

2.2 仿真分析

設(shè)置模擬退火算法初始溫度T=150，降溫速率q=0.9，結(jié)束溫度Tend=1e-06，Metropolis準(zhǔn)則鏈長(zhǎng)L=1 000，運(yùn)行模擬退火算法得到結(jié)果：[x1，x2，x3，x4]=[3 3 0 0]，目標(biāo)函數(shù)為930，其目標(biāo)函數(shù)變化曲線如圖1所示。

為了驗(yàn)證模擬退火算法得到結(jié)果的準(zhǔn)確性，我們采取分支界定算法對(duì)問題（3）進(jìn)行求解，得到的變量的解和最優(yōu)目標(biāo)函數(shù)值為[x1，x2，x3，x4]=[3 1 3 0]，目標(biāo)函數(shù)為930。很明顯，模擬退火算法和分支界定算法得到的最優(yōu)目標(biāo)函數(shù)均為930，但是兩者得到的解不相同。這說明了組合優(yōu)化問題（3）有多個(gè)解。為了求解組合優(yōu)化問題（3）的解，我們定義矩陣RESULT保存迭代過程中最優(yōu)解及其函數(shù)值的信息，前4列保存解，第5列保存該解所對(duì)應(yīng)的目標(biāo)函數(shù)，矩陣RESULT的初始值=[迭代步驟（2）隨機(jī)生成的初始解S1，S1對(duì)應(yīng)的目標(biāo)函數(shù)]，對(duì)模擬退火算法每次迭代的Metropolis準(zhǔn)則進(jìn)行改進(jìn)，具體的偽代碼如下所示。

隨機(jī)生成新解S2;

S1評(píng)價(jià)指標(biāo)=以S1為變量帶入目標(biāo)函數(shù)，返回的? 目標(biāo)函數(shù)值;

S2評(píng)價(jià)指標(biāo)=以S2為變量帶入目標(biāo)函數(shù)，返回的? 目標(biāo)函數(shù)值;

評(píng)價(jià)指標(biāo)差= S2評(píng)價(jià)指標(biāo)- S1評(píng)價(jià)指標(biāo)

if評(píng)價(jià)指標(biāo)差<0

S1=S2;

S1評(píng)價(jià)指標(biāo)=S2評(píng)價(jià)指標(biāo);

end

if評(píng)價(jià)指標(biāo)差>0

if exp（-dC/溫度）>0-1之間的均勻隨機(jī)數(shù)

S1=S2;

S1評(píng)價(jià)指標(biāo)=S2評(píng)價(jià)指標(biāo);

end

end

if評(píng)價(jià)指標(biāo)差==0

flag=0;

if 矩陣RESULT第一行最后一列>S2評(píng)價(jià)指標(biāo)

矩陣RESULT置為空矩陣;

RESULT=[S2，S2評(píng)價(jià)指標(biāo)];

end

if矩陣RESULT第一行最后一列==S2評(píng)價(jià)指標(biāo)

for i=1：矩陣RESULT的行數(shù)

if矩陣RESULT第i行前4列與S2相同

置變量flag為1，跳出當(dāng)前循環(huán);

end

end

if flag==0

RESULT=[RESULT;[S2，S2評(píng)價(jià)指標(biāo)]];

end

end

end

現(xiàn)在根據(jù)我們改進(jìn)過的模擬退火算法，再次執(zhí)行程序，得到結(jié)果見表2。在表2中，我們得到了3組最優(yōu)解，其最優(yōu)目標(biāo)函數(shù)值均為930。說明這3組解均為原問題（3）的最優(yōu)解。

3 總結(jié)

本文利用模擬退火算法求解組合優(yōu)化問題，仿真結(jié)果證明了模擬退火算法的有效性，通過進(jìn)一步分析與比較仿真結(jié)果與分支界定算法仿真結(jié)果，發(fā)現(xiàn)仿真算例存在多個(gè)極值點(diǎn)，我們?cè)谠谢A(chǔ)上進(jìn)一步改進(jìn)模擬退火算法的Metropolis準(zhǔn)則，同時(shí)輸出組合優(yōu)化調(diào)度問題的多個(gè)極值點(diǎn)。

參 考 文 獻(xiàn)

[1]龔純，王正林.精通MATLAB最優(yōu)化計(jì)算[M].北京：電子工業(yè)出版社，2009.

[2]吳辰斌，李海明，劉棟，等.一種改進(jìn)型粒子群優(yōu)化算法在電力系統(tǒng)經(jīng)濟(jì)負(fù)荷分配中的應(yīng)用[J].電力系統(tǒng)保護(hù)與控制，2016（10）：44-48.

[3]何大闊，王福利，毛志忠.基于改進(jìn)遺傳算法的電力系統(tǒng)經(jīng)濟(jì)負(fù)荷分配[J].控制與決策，2007（2）：230-232.

[4]王述紅，朱寶強(qiáng)，王鵬宇.模擬退火聚類算法在結(jié)構(gòu)面產(chǎn)狀分組中的應(yīng)用[J].東北大學(xué)學(xué)報(bào)（自然科學(xué)版），2020，41（9）：1328-1333.

[5]鄭小虎，鮑勁松，馬清文，等.基于模擬退火遺傳算法的紡紗車間調(diào)度系統(tǒng)[J].紡織學(xué)報(bào)，2020，41（6）：36-41.

[6]陳科勝，鮮思東，郭鵬.求解旅行商問題的自適應(yīng)升溫模擬退火算法[J].控制理論與應(yīng)用，2021，38（2）：245-

254.