曾明華，全 軻

（華東交通大學(xué)交通運(yùn)輸與物流學(xué)院，南昌330013）（*通信作者電子郵箱296575845@qq.com）

0 引言

雙層規(guī)劃（Bi-Level Programming，BLP）是一類具有主從遞階關(guān)系結(jié)構(gòu)的數(shù)學(xué)模型，它是將優(yōu)化問(wèn)題作為約束條件的極值問(wèn)題［1］。該數(shù)學(xué)模型已廣泛應(yīng)用于資源分配、交通網(wǎng)絡(luò)設(shè)計(jì)等實(shí)際問(wèn)題。文獻(xiàn)［2］證明了線性雙層規(guī)劃是NP-難問(wèn)題。事實(shí)上，大多數(shù)雙層規(guī)劃還包括更復(fù)雜的非線性雙層規(guī)劃問(wèn)題。求解雙層規(guī)劃的方法大致可分為傳統(tǒng)算法和智能優(yōu)化算法，但傳統(tǒng)算法依賴于目標(biāo)函數(shù)的可微性，不具有普遍的適用性。目前，對(duì)目標(biāo)函數(shù)要求不高的智能優(yōu)化算法已經(jīng)廣泛用于解決雙層規(guī)劃問(wèn)題。

目前，已經(jīng)有學(xué)者結(jié)合布谷鳥(niǎo)搜索（Cuckoo Search，CS）算法和粒子群優(yōu)化（Particle Swarm Optimization，PSO）算法并用于求解一些實(shí)際應(yīng)用問(wèn)題。文獻(xiàn)［3］將CS算法中的Lévy飛行和淘汰機(jī)制引入PSO算法中，形成混合優(yōu)化算法。文獻(xiàn)［4］將CS 算法引入動(dòng)態(tài)多種群（Dynamic Multi-Swarm，DMS）-PSO算法，提高了原算法的全局搜索能力，提出DMS-PSO-CS 算法。文獻(xiàn)［5］基于種群拓?fù)浣Y(jié)構(gòu)與粒子變異，同時(shí)加以CS 算法的偏好隨機(jī)游走變異策略，提出基于拓?fù)浣Y(jié)構(gòu)與粒子變異改進(jìn)的粒子群優(yōu)化算法。文獻(xiàn)［6］用CS算法中的Lévy飛行取代PSO 算法位置更新公式的隨機(jī)數(shù)，提出一種PSO-CS 算法。文獻(xiàn)［7］通過(guò)正交平方原理建立初始種群，動(dòng)態(tài)調(diào)整CS 參數(shù)以最小化CS 的固定步長(zhǎng)的影響，并將PSO 算法引入CS 算法中形成一種混合CSPSO 算法。文獻(xiàn)［8］引入QPSO 算法用于CS 算法的位置尋優(yōu)過(guò)程，以解決CS 算法在解決服務(wù)質(zhì)量（Quality of Service，QoS）組播路由問(wèn)題收斂速度慢、算法搜索效率低的問(wèn)題。

根據(jù)以往的研究，混合CSPSO 算法相較于傳統(tǒng)PSO 算法有更好的性能，但目前對(duì)QPSO 算法與CS 算法的混合算法研究較少，且尚未有學(xué)者將混合CSQPSO 算法用于求解雙層規(guī)劃。

本文結(jié)合量子行為粒子群優(yōu)化算法、布谷鳥(niǎo)搜索算法以及模擬退火（Simulated Annealing，SA）算法以尋找雙層規(guī)劃問(wèn)題最優(yōu)解。對(duì)CS算法設(shè)計(jì)了一種改進(jìn)動(dòng)態(tài)步長(zhǎng)Lévy飛行，并引入模擬退火算法中的Metropolis 準(zhǔn)則，進(jìn)而提出了一種改進(jìn)混合布谷鳥(niǎo)搜索量子行為粒子群優(yōu)化（Improved hybrid Cuckoo Search-based Quantum-behaved Particle Swarm Optimization，ICSQPSO）算法。該算法能有效增加粒子群多樣性，增強(qiáng)粒子群跳出局部最優(yōu)解的能力。對(duì)13 個(gè)復(fù)雜的雙層規(guī)劃的數(shù)值實(shí)驗(yàn)表明，所提出的算法非常有效，且大部分結(jié)果顯著優(yōu)于所參考的四種改進(jìn)粒子群算法和一種改進(jìn)遺傳算法。

1 雙層規(guī)劃問(wèn)題描述

雙層規(guī)劃數(shù)學(xué)模型一般可以表述如式（1）：

其中：x為上層決策變量，y為下層決策變量，，；X是設(shè)定了對(duì)上層決策變量的額外限制，例如上下界 。，R是實(shí)數(shù)集。

對(duì)于雙層規(guī)劃模型給出如下基本概念［9］。

1）約束區(qū)域：S={(x，y)|h(x，y) ≤0，g(x，y) ≤0}。

2）對(duì)于給定的x∈X，下層規(guī)劃可行域?yàn)椋?/p>

S(x)={y|g(x，y) ≤0}

3）S在上層決策空間的投影為：S(X)={x|?y，(x，y)∈S}。

4）對(duì)于每個(gè)x∈S(X)，下層合理反應(yīng)集為：

Q(x)={y|y∈arg min{f(x，y)，y∈S(x)}}

5）誘導(dǎo)域：IR={(x，y)|(x，y)∈S，y∈Q(x)}。

進(jìn)一步給出雙層規(guī)劃模型可行解和最優(yōu)解的概念：若點(diǎn)(x，y)∈IR，則稱點(diǎn)(x，y)為雙層規(guī)劃模型的可行解。若點(diǎn)(x*，y*)是雙層規(guī)劃模型的可行解，且對(duì)任意(x，y)∈IR有F(x*，y*)≤F(x，y)，則稱(x*，y*)為雙層規(guī)劃模型的最優(yōu)解。

2 求解雙層規(guī)劃問(wèn)題的ICSQPSO算法

2.1 標(biāo)準(zhǔn)算法

2.1.1 標(biāo)準(zhǔn)QPSO算法

PSO 算法是由Kennedy 等［10］在1995年提出的一種智能優(yōu)化算法，具有實(shí)現(xiàn)容易、收斂快等優(yōu)點(diǎn)；但是原始PSO 算法在優(yōu)化后期由于種群在搜索空間中多樣性的丟失容易早熟收斂，陷入局部最優(yōu)解。文獻(xiàn)［11］提出了一種量子行為粒子群優(yōu)化算法。在QPSO 算法中，由概率密度函數(shù)描述的束縛狀態(tài)的粒子可以以一定概率出現(xiàn)在整個(gè)可行搜索空間的任何區(qū)域，具有全局收斂性。

假設(shè)算法搜索區(qū)域的維度是D，粒子群粒子數(shù)目是N。代表第i個(gè)粒子在j維空間中的位置。pi=代表第i個(gè)粒子在j維空間中當(dāng)前搜索到的最優(yōu)位置，代表整個(gè)粒子群在j維空間中當(dāng)前搜索到的最優(yōu)位置，i=1，2，…，N，j=1，2，…，D。相較PSO 算法的粒子位置和速度更新模型，QPSO算法僅采用粒子位置更新模型，如式（2）～（6）所示：

其中：

其中：Pi(t)稱為吸引子，這個(gè)點(diǎn)在迭代過(guò)程中會(huì)不斷吸引粒子靠近，φ(t)和ui(t)是［0，1］上均勻分布的隨機(jī)數(shù)；β稱為收縮—擴(kuò)張系數(shù)，用于協(xié)調(diào)迭代過(guò)程中粒子全局和局部的搜索性能，常采用線性減小的方式，如式（4），M是最大迭代次數(shù)；t是當(dāng)前迭代次數(shù)；Li(t)表示當(dāng)前粒子位置與平均個(gè)體最好位置C(t)距離的絕對(duì)值。

2.1.2 標(biāo)準(zhǔn)CS算法

CS 算法是一種新興智能算法，由Yang 等［12］于2009 年提出，其主要思想是模擬布谷鳥(niǎo)借巢產(chǎn)卵的習(xí)性，利用鳥(niǎo)類的Lévy 飛行機(jī)制，通過(guò)偏好隨機(jī)游走搜索到最優(yōu)的鳥(niǎo)巢以孵化鳥(niǎo)蛋。該算法的優(yōu)點(diǎn)是不易陷入局部最優(yōu)且全局搜索能力強(qiáng)，但收斂速度較慢。文獻(xiàn)［12］假設(shè)了CS 算法中的3 種理想狀態(tài)：

1）每次布谷鳥(niǎo)只生產(chǎn)一枚卵，并隨機(jī)選擇一個(gè)鳥(niǎo)窩孵化；

2）在搜索鳥(niǎo)窩的過(guò)程中，保留擁有最高質(zhì)量卵的鳥(niǎo)窩到下一代；

3）每一代的鳥(niǎo)窩數(shù)量固定，鳥(niǎo)窩中的卵被宿主發(fā)現(xiàn)的概率是p。

基于上述第3個(gè)假設(shè)，CS算法通過(guò)式（7）將搜索的位置和路徑更新：

其中：zi(t)代表第i個(gè)鳥(niǎo)窩在第j維空間中第t代時(shí)的位置；α為步長(zhǎng)縮放因子，通常α為固定值，取值為0.01；⊕為點(diǎn)對(duì)點(diǎn)乘法；L(λ)為L(zhǎng)évy飛行的隨機(jī)搜索路徑。

針對(duì)第三個(gè)假設(shè)，利用偏好隨機(jī)游走機(jī)制生成新解，即用隨機(jī)數(shù)r∈[0，1]與p對(duì)比，若r＞p，則對(duì)zi按式（8）進(jìn)行隨機(jī)改變，反之不變：

其中：rand為[0，1]區(qū)間的隨機(jī)數(shù)；zq(t)和zk(t)是種群中隨機(jī)選取的兩個(gè)不同個(gè)體。文獻(xiàn)［13］已詳細(xì)介紹CS算法，本文不再贅述。

2.1.3 標(biāo)準(zhǔn)SA算法

模擬退火（SA）的思想起源于固體的退火過(guò)程，最早由Metropolis 于1953 年提出［14］。當(dāng)SA 算法用于最優(yōu)化問(wèn)題時(shí)，首先要設(shè)定一個(gè)初始溫度TS、終止溫度Tf以及溫度衰減系數(shù)?。溫度衰減函數(shù)采用常用線性函數(shù)，如式（9）：

其中T(t)為某時(shí)刻t的溫度。

在退火過(guò)程中溫度為T(t)的時(shí)刻，固體從狀態(tài)h轉(zhuǎn)變到另一個(gè)狀態(tài)l遵循如式（10）所示的Metropolis準(zhǔn)則：其中：K為波爾茲曼常數(shù)；E(l)和E(h)分別為在狀態(tài)l和h下固體的內(nèi)能，在最優(yōu)化問(wèn)題中，固體的內(nèi)能可以看作是目標(biāo)函數(shù)值。在Metropolis 準(zhǔn)則下，SA 算法一定能接受最優(yōu)解，同時(shí)也能以一定概率接受較差解。正因?yàn)橛羞@種機(jī)制，SA 算法不易陷入局部最優(yōu)解。

2.2 改進(jìn)的Lévy飛行

式（7）中的步長(zhǎng)縮放因子α是一個(gè)重要參數(shù)，在原始CS算法中，α是固定常數(shù)，可能會(huì)影響算法的收斂速度［15］。為了提高收斂速度，本文提出一種如式（11）所示的分段函數(shù)以動(dòng)態(tài)調(diào)整α，α隨著迭代次數(shù)的增加而減小。這與PSO 算法中隨著迭代次數(shù)增加而慣性權(quán)重減小的原理類似。在迭代初期，α足夠大以促使CS 算法搜索到更多可行解，保證了解的多樣性；在迭代末期，α較小，更好地微調(diào)可行解向最優(yōu)解靠攏。

其中：αmax和αmin分別是α的最大值和最小值。

2.3 約束條件處理

本文假設(shè)上層和下層函數(shù)分別存在m個(gè)不等式約束，采用罰函數(shù)法處理式（1）的約束條件。以下層規(guī)劃函數(shù)為例，可以根據(jù)式（12）計(jì)算粒子的適應(yīng)度f(wàn)itness(x，y)：

其中：S為足夠大的正常數(shù)，稱為懲罰因子。

2.4 ICSQPSO算法步驟

ICSQPSO 算法流程如圖1 所示，其詳細(xì)步驟包括上層規(guī)劃模型的求解算法（算法1）和下層規(guī)劃模型的求解算法（算法2）。算法1與算法2分別描述如下。

圖1 ICSQPSO算法流程Fig.1 Flowchart of ICSQPSO algorithm

算法1 上層規(guī)劃模型的求解算法。

1）初始化參數(shù)：最大迭代次數(shù)M1，迭代次數(shù)t1=0，初始溫度TS、終止溫度Tf以及溫度衰減系數(shù)?。

2）將上層規(guī)劃模型的解x*j(t1)（j=1，2，…，n1）代入下層規(guī)劃模型，調(diào)用算法2求得下層規(guī)劃模型的最優(yōu)解y*w(t1)（w=1，2，…，n2）。

3）將y*w(t1)代入上層規(guī)劃模型，調(diào)用算法2 求得上層規(guī)劃模型的最優(yōu)解x*j(t1)。

4）將x*j(t1)與y*w(t1)代入上層規(guī)劃的目標(biāo)函數(shù)，計(jì)算目標(biāo)函數(shù)值，按式（10）進(jìn)行Metropolis 準(zhǔn)則判斷，若接受，則記為Fbest。

5）按式（9）進(jìn)行退火降溫。置t1←t1+1。

6）若t1≥M1，T(t1)≤Tf，則轉(zhuǎn)至7）；否則，返回2）。

7）結(jié)束計(jì)算。輸出群體最優(yōu)解(x*j(t1)，y*w(t1))以及上層規(guī)劃和下層規(guī)劃的目標(biāo)函數(shù)值。

算法2 下層規(guī)劃模型的求解算法。

1）初始化參數(shù)：粒子群數(shù)目N，收縮-擴(kuò)張系數(shù)βmax和βmin，步長(zhǎng)縮放因子上下限αmax和αmin，最大迭代次數(shù)M2，迭代次數(shù)t2=0，懲罰因子S，發(fā)現(xiàn)概率p，每個(gè)決策變量的上下界；隨機(jī)初始化粒子群中粒子的位置zi(t2)；將第i個(gè)粒子的pi設(shè)為該粒子的當(dāng)前位置，并計(jì)算各個(gè)粒子的適應(yīng)度，設(shè)pg為初始粒子群中最佳粒子的位置。

2）粒子按式（7）進(jìn)行Lévy 飛行，并利用式（2）更新粒子的位置。

3）檢查更新完后的zi(t2)中各個(gè)決策變量大小是否超過(guò)上界或低于下界。若超過(guò)上界，則將其控制在上界邊界處；若低于下界，則將其控制在下界邊界處。

4）更新個(gè)體最優(yōu)解pi和群體最優(yōu)解pg。如果粒子i的適應(yīng)度優(yōu)于其pi的適應(yīng)度，則該粒子的pi更新為當(dāng)前位置zi(t2)；如果粒子i的適應(yīng)度優(yōu)于當(dāng)前pg的適應(yīng)度，則pg更新為當(dāng)前位置zi(t2)。

5）在一定概率p的條件下，利用式（8）對(duì)zi(t2)進(jìn)行隨機(jī)更新。置t2←t2+1。

6）若t2≥M2，則返回群體最優(yōu)解pg至算法1；否則，返回2）。

3 算例和算法性能分析

本文選取了13 個(gè)算例來(lái)驗(yàn)證ICSQPSO 算法的有效性。算例既包括線性雙層規(guī)劃，也包括高度非線性雙層規(guī)劃，還有分式規(guī)劃以及含有多個(gè)下層規(guī)劃的雙層規(guī)劃。其中算例1～4選自文獻(xiàn)［15］的例1、例2、例5、例8，算例5～7 選自文獻(xiàn)［16］的例7～9，算例8 選自文獻(xiàn)［17］的例11，算例9 選自文獻(xiàn)［18］的例2，算例10～13選自文獻(xiàn)［19］的例1～4。

ICSQPSO 算法參數(shù)設(shè)置：粒子群數(shù)目N=40，初始溫度TS=10 000，終止溫度Tf=0.36，溫度衰減系數(shù)?=0.95，收縮-擴(kuò)張系數(shù)最大值βmax=1，最小值βmin=0.5，步長(zhǎng)縮放因子最大值αmax=0.5，最小值αmin=0.01，最大迭代次數(shù)M1=200，M2=300，經(jīng)過(guò)數(shù)值實(shí)驗(yàn)測(cè)試選定ta=150，懲罰因子S=100 000，發(fā)現(xiàn)概率p=0.25。圖2 給出了算例1 的上層規(guī)劃目標(biāo)函數(shù)值的演進(jìn)過(guò)程，可以看出，ICSQPSO 算法在求解過(guò)程中快速收斂，算法在其他算例中也表現(xiàn)出相同的效果。

圖2 算例1目標(biāo)函數(shù)值的變化過(guò)程Fig.2 Changing process of objective function value of example 1

每個(gè)算例單獨(dú)運(yùn)行30次，記錄30次運(yùn)行中的上層規(guī)劃函數(shù)F(x，y)最優(yōu)值Fbest、下層規(guī)劃函數(shù)f(x，y)最優(yōu)值fbest、以及上下層規(guī)劃函數(shù)最優(yōu)值對(duì)應(yīng)的最優(yōu)解(x*，y*)。Fbest、fbest的 結(jié)果如表1 所 示，(x*，y*) 的 結(jié)果如表2 所 示，并 按(Fbest文獻(xiàn)-Fbest本文)/Fbest文獻(xiàn)計(jì)算了ICSQPSO 算法相比文獻(xiàn)［15-19］算法所得Fbest的減少（優(yōu)化）程度，計(jì)算結(jié)果如表1倒數(shù)第二列所示；同時(shí)，按相似方法計(jì)算比較了ICSQPSO 算法與文獻(xiàn)［15-19］的比較文獻(xiàn)算法的Fbest的減少（優(yōu)化）程度，如表1 最后一列所示。這里，所謂文獻(xiàn)［15-19］的比較文獻(xiàn)算法，是指文獻(xiàn)［15-19］用來(lái)與其所提出算法進(jìn)行計(jì)算結(jié)果比較的算法。

本文將ICSQPSO 算法的測(cè)試結(jié)果與文獻(xiàn)［15-19］算法進(jìn)行對(duì)比。用于對(duì)比的算法有：文獻(xiàn)［15］提出的PSO-CST（Particle Swarm Optimization with Chaos Searching Technique）算法；文獻(xiàn)［16］提出的IBPSO（Improved Bilevel Particle Swarm Optimization）算法；文 獻(xiàn)［17］提出的HPSOBLP（Bi-Level Programming based on Hierarchical Particle Swarm Optimization）算法；文獻(xiàn)［18］提出的IPSO-BLPP（Bi-Level Programming Problem based on Improved Particle Swarm Optimization）算法；文獻(xiàn)［19］提出的一種基于新編碼方式的改進(jìn)遺傳算法。

表1 目標(biāo)函數(shù)最優(yōu)值對(duì)比Tab.1 Comparison of objective function optimal value

表2 不同算法之間決策變量最優(yōu)解對(duì)比Tab.2 Comparison of optimal solutions of decision variables among different algorithms

從表1 和表2 可以看出，除了算例5 和9，本文提出的ICSQPSO 算法的計(jì)算結(jié)果均優(yōu)于文獻(xiàn)［15-19］的計(jì)算結(jié)果。其中對(duì)于復(fù)雜的非線性雙層規(guī)劃算例1，ICSQPSO 算法發(fā)現(xiàn)了一個(gè)優(yōu)秀的解（16.977 4，7.814 3，10.000 0，0.000 0），其對(duì)應(yīng)的上層規(guī)劃目標(biāo)函數(shù)最優(yōu)值Fbest=118.080 5，較文獻(xiàn)［15］的Fbest=232.521 9大幅減小了49.22%；對(duì)于高度非線性雙層規(guī)劃算例6，ICSQPSO 算法求得了一個(gè)出色的解（0.418 8，0.394 2，1.818 3），其上層規(guī)劃目標(biāo)函數(shù)最優(yōu)值Fbest=0.446 1，較文獻(xiàn)［16］的Fbest=1.066 5 減小了58.17%；對(duì)于包含多個(gè)下層規(guī)劃的高維非線性雙層規(guī)劃算例13，ICSQPSO 算法能求解出比文獻(xiàn)［19］的改進(jìn)遺傳算法更好的解，ICSQPSO算法求得Fbest=-0.712 8，而文獻(xiàn)［19］求得Fbest=0.355 3，這表明本文算法在求解多下層分式雙層規(guī)劃時(shí)表現(xiàn)良好。

綜上所述，對(duì)于大多數(shù)算例，本文提出的ICSQPSO 算法能找到比文獻(xiàn)［15-19］算法更好的解，能有效地求解線性雙層規(guī)劃、單下層分式線性雙層規(guī)劃、多下層線性雙層規(guī)劃以及復(fù)雜的非線性雙層規(guī)劃。

4 結(jié)語(yǔ)

針對(duì)PSO 算法求解雙層規(guī)劃容易陷入局部最優(yōu)解的問(wèn)題，從增加粒子群多樣性、增強(qiáng)粒子群跳出局部最優(yōu)解的能力的角度出發(fā)，基于SA 算法中的Metropolis 準(zhǔn)則，并將改進(jìn)CS算法的動(dòng)態(tài)步長(zhǎng)Lévy 飛行和偏好隨機(jī)游走機(jī)制嵌入QPSO 算法中，提出了一種ICSQPSO 算法。該算法具有保持種群多樣性與提高全局搜索能力的優(yōu)點(diǎn)：一方面，采用QPSO 算法代替PSO 算法，并引入Metropolis 準(zhǔn)則，在求解過(guò)程中既能接受好解也能以一定的概率接受壞解，增強(qiáng)全局尋優(yōu)能力；另一方面，改進(jìn)CS算法中的動(dòng)態(tài)步長(zhǎng)Lévy飛行能使粒子群在優(yōu)化中保持較高的多樣性，保證搜索廣度，偏好隨機(jī)游走機(jī)制幫助粒子跳出局部最優(yōu)解。對(duì)13 個(gè)具有典型代表性的算例進(jìn)行數(shù)值實(shí)驗(yàn)，結(jié)果表明，所提出的ICSQPSO 算法能得到顯著優(yōu)于文獻(xiàn)［15-19］算法的結(jié)果，對(duì)求解雙層規(guī)劃模型是非常有效的。