基于目標(biāo)空間分解和連續(xù)變異的多目標(biāo)粒子群算法

錢小宇，葛洪偉，蔡明

（1. 江南大學(xué) 輕工過程先進(jìn)控制教育部重點(diǎn)實(shí)驗(yàn)室，江蘇 無錫 214122; 2. 江南大學(xué) 物聯(lián)網(wǎng)工程學(xué)院，江蘇 無錫 214122; 3. 江南大學(xué) 信息化建設(shè)與管理中心，江蘇 無錫 214122）

粒子群算法自1995年被Kennedy和Eberhart提出后[1]，由于其簡(jiǎn)單高效，收斂速度快，逐漸在優(yōu)化算法中脫穎而出。隨之粒子群算法被Coelloran應(yīng)用到多目標(biāo)優(yōu)化問題上[2-3]，得到了各界的廣泛認(rèn)可。后來更多學(xué)者對(duì)多目標(biāo)粒子群算法(MOPSO)進(jìn)行了更深入的研究：Raquel等[4]使用了擁擠距離作為適應(yīng)值對(duì)粒子進(jìn)行排序；Li等[5]通過整合粒子的GMR(global margin ranking)和粒子密度信息的方法來對(duì)粒子進(jìn)行排序，能高效快速地選擇出pbest和gbest粒子；除了對(duì)排序選擇的方法進(jìn)行優(yōu)化外，有的利用聚合函數(shù)來對(duì)多個(gè)目標(biāo)函數(shù)進(jìn)行處理[6]；有的使用了目標(biāo)空間分解的方法來保障粒子最優(yōu)解的多樣性[7]；以及將MOPSO和其他優(yōu)化算法進(jìn)行交叉使用，如教與學(xué)方法[8]、差分方法[9]等, 通過一定的比例調(diào)用MOPSO。

上述的這些方法不斷地對(duì)MOPSO的多樣性和全局收斂進(jìn)行改進(jìn)和優(yōu)化，雖然取得了較好的研究成果，但是這些方法的改進(jìn)策略只對(duì)MOPSO某一性能效果進(jìn)行改進(jìn)，其他的性能指標(biāo)并沒有同時(shí)得到很好的改進(jìn)。例如：文獻(xiàn)[7]中，多樣性方面得到很大的提升，但是其收斂性方面仍有很大不足。所以，本文借鑒Pareto支配強(qiáng)度[10]和連續(xù)變異的方法[11]對(duì)基于目標(biāo)空間分解方法的多目標(biāo)粒子群算法[7]進(jìn)行優(yōu)化，在給每個(gè)子區(qū)域分配粒子時(shí)，刪除沒有歸屬的粒子，給沒分到粒子的區(qū)域初始化新的粒子，增加獲得較優(yōu)粒子的可能性。這樣既從整體和局部?jī)煞矫嫣嵘肆Ｗ拥氖諗啃裕瑫r(shí)多樣性也得到了一定的優(yōu)化提升，這種新的算法被稱為基于目標(biāo)空間分解和連續(xù)變異的多目標(biāo)粒子群優(yōu)化算法(multi-objective particle swarm optimization based on decomposition and continuous mutation ，MOPSO/DC)。

1 多目標(biāo)優(yōu)化問題

1.1 基本概念

式中： x =(x1,x2,···,xn) 為 n 維 決策變量； m 為目標(biāo)函數(shù)的個(gè)數(shù)； g ( x) 函 數(shù)為目標(biāo)函數(shù)的 q 個(gè)不等式約束； h ( x) 為 目標(biāo)函數(shù)的 p 個(gè)等式約束。所有這些滿足條件的決策變量用集合 ? 表 示，Y={F(x)|x∈?}為目標(biāo)空間。接下來介紹4個(gè)關(guān)于多目標(biāo)問題的重要定義。

定義1 Pareto支配。解 d , e ∈? ， d 支 配 e，記為 d ? e，滿足下面的兩個(gè)關(guān)系式：

定義2 Pareto最優(yōu)。如果 x 是Pareto最優(yōu)解，則在 ? 中 ， ? ?z∈? 使 z ? x 成立。

定義3 Pareto最優(yōu)解集(P S )：

1.2 粒子群優(yōu)化

粒子群算法中的粒子由速度信息和位置信息組成，更新公式為

式中：k指粒子群中第k個(gè)粒子；t為當(dāng)前迭代次數(shù)；W為權(quán)衡局部搜索和全局搜索的參數(shù)，W∈[0.1,0.9]； C1和 C2為學(xué)習(xí)因子，其大小都為2；R1和 R2都是[0，1]之間的隨機(jī)數(shù)；P指當(dāng)前粒子最好位置pbest；G指引導(dǎo)粒子gbest的位置；指在第t次迭代中第k個(gè)粒子的速度；指在第t次迭代中第k個(gè)粒子的位置。

2 MOPSO/DC算法

2.1 目標(biāo)空間分解

把 目標(biāo) 空間 Y 分 成 M 個(gè) 子 區(qū) 域 Y1，Y2，···，YM。令 j ∈ { 1,2,···,M}，對(duì)任一給定的第j個(gè)子區(qū)域，每個(gè)目標(biāo)函數(shù)在所有目標(biāo)函數(shù)中所占有的權(quán)重所組成的向量(a1,a2,···,am) 定 義為該子區(qū)域的中心向量 Aj。當(dāng)目標(biāo)函數(shù)個(gè)數(shù) m =2 時(shí) ，則第 j 子區(qū)域的中心向量 Aj表示為當(dāng) m =3 時(shí)，進(jìn)行兩層循環(huán)，令 k1為第1層(外層)循環(huán)變量， k2為第2層(內(nèi)層)循環(huán)變量， k1從 0取到h， k2從0取到 h - k1， 每 次 循 環(huán) 令其中h和M滿足：當(dāng) 時(shí)h取最小值，則第u個(gè)子區(qū)域的中心向量為當(dāng)m>3 時(shí) ，進(jìn)行 m - 1 層 循環(huán)， k1為第一層循環(huán)變量(最外層)，依次由外向內(nèi)，令 ki為第i層循環(huán)變量，km-1為第 m - 1 層 (最內(nèi)層)循環(huán)變量， k1從0取到h，k2從0取到 h - k1， ki從 0取到 ( h -k1-k2- · ··ki)，最里層循環(huán)變量km-1從0取到 ( h -k1-k2-···km-2)，每次循環(huán)令生個(gè)中心向量，這些向量的下標(biāo)為若M ， 則 從 下 標(biāo) u = 2 開始，以為步長(zhǎng)依次刪除M 個(gè)下標(biāo)對(duì)應(yīng)的中心向量，然后將剩下的M個(gè)中心向量依次作為 Aj，j=1, 2, ···, M。通過每個(gè)子區(qū)域的中心向量找出該子區(qū)域的T個(gè)相鄰的子區(qū)域, 參考指標(biāo)為兩個(gè)中心向量的余弦值。

2.2 粒子的分類與更新

在對(duì)粒子進(jìn)行分配時(shí)，通過參考點(diǎn)確定每個(gè)粒子的方向向量[7]。令參考點(diǎn)為 R ( r1,r2,···,rm)，其中為目標(biāo)函數(shù)個(gè)數(shù)。例如：如圖1所示， f1和 f2表示兩個(gè)目標(biāo)函數(shù)，m=2，則圖中粒子C的方向向量為從點(diǎn)R(r1,r2)指向點(diǎn)C的向量X，X=C-R。通過比較粒子方向向量和所有子區(qū)域中心向量的余弦值，確定該粒子屬于最大余弦值所對(duì)應(yīng)的中心，向量子區(qū)域。

圖 1 粒子的方向向量Fig. 1 Direction vector of the particle

在對(duì)粒子進(jìn)行分配時(shí)會(huì)有以下兩種特殊情況：

1)有的子區(qū)域粒子數(shù)大于子區(qū)域的容量Vol時(shí)，通過適應(yīng)值進(jìn)行取舍，適應(yīng)值公式為式中：T是Pareto支配強(qiáng)度，其數(shù)值為當(dāng)前粒子支配的粒子數(shù)，在適應(yīng)值計(jì)算中加入了Pareto支配強(qiáng)度T，增強(qiáng)每個(gè)子區(qū)域中的粒子趨向真實(shí)PF的能力；參數(shù)a表示支配強(qiáng)度對(duì)適應(yīng)值的影響程度，，M指目標(biāo)空間子區(qū)域數(shù)；CD為擁擠距離，通過每個(gè)目標(biāo)函數(shù)值對(duì)粒子進(jìn)行排序，序列兩端粒子在當(dāng)前目標(biāo)函數(shù)中的擁擠距離設(shè)為5，其他粒子在當(dāng)前目標(biāo)函數(shù)中的擁擠距離為在序列中該粒子前后兩粒子的目標(biāo)函數(shù)值之差的絕對(duì)值，最后將求出的該粒子在每個(gè)目標(biāo)函數(shù)中計(jì)算的擁擠距離之和作為當(dāng)前該粒子的擁擠距離CD。

計(jì)算該子區(qū)域包含的粒子的適應(yīng)值，從大到小排序，選擇序列中前30%的粒子，再?gòu)倪@些粒子中選擇離所在子區(qū)域中心向量最近的Vol個(gè)粒子，多余的劣質(zhì)粒子刪除。

2)當(dāng)子區(qū)域粒子數(shù)少于該子區(qū)域容量Vol時(shí)，由于情況1)中已經(jīng)對(duì)劣質(zhì)的粒子進(jìn)行刪除，為了盡可能保證整體的優(yōu)越性，所以重新初始化該子區(qū)域所缺數(shù)目的粒子作為該區(qū)域的粒子，增加獲得較優(yōu)粒子的可能性，增加粒子多樣性。然后計(jì)算新粒子的目標(biāo)函數(shù)值。

2.3 連續(xù)變異操作

當(dāng)只采用一種方法進(jìn)行變異時(shí)，不能兼顧全局粒子和局部粒子的特性。在MOPSO/DC中提出差分+柯西[12]+高斯[13]連續(xù)的變異策略，公式為式中：X為當(dāng)前引導(dǎo)粒子gbest的位置；X1和X2是從EPOP中隨機(jī)選出的兩個(gè)不同粒子的位置。這樣在變異時(shí)可保留全局粒子的一部分性質(zhì)，起到一定信息交流的作用，增加了多樣性。t為當(dāng)前迭代的次數(shù)，gmax為最大迭代次數(shù)， g ( 1)=2，g(t+1)=和分別是對(duì)gbest的位置進(jìn)行變異操作之后產(chǎn)生的新粒子。每次變異后讓產(chǎn)生的新粒子和當(dāng)前的gbest進(jìn)行比較，選擇支配權(quán)優(yōu)先的作為gbest，然后進(jìn)行后續(xù)變異操作，直到3次變異結(jié)束，最終確定引導(dǎo)粒子gbest。

高斯變異步長(zhǎng)較短，能很好地吸取局部粒子性質(zhì)，柯西變異具有相對(duì)大的步長(zhǎng)，能進(jìn)行較大范圍的變異，具有全局范圍變異的特性，這樣能很好地產(chǎn)生全新的粒子。觀地看出，柯西分布具有較寬的變異范圍，粒子具有較大范圍的變異。用這3種方法進(jìn)行變異操作，既增強(qiáng)了

圖 2 標(biāo)準(zhǔn)高斯和分布柯西分布對(duì)比Fig. 2 Comparison of standard Gaussian and Cauchy distributions

gbest的引導(dǎo)能力，同時(shí)gbest吸取局部或全局粒子的性質(zhì)，促進(jìn)了粒子之間信息交流共享，增加了多樣性。

2.4 MOSPO/DC步驟

1)初始化2N個(gè)粒子的位置，這些粒子的位置集合記為POP，初始化N個(gè)速度，這些速度集合記為V，設(shè)目標(biāo)函數(shù)的個(gè)數(shù)為m，當(dāng)前迭代次數(shù)為t，初始化最大迭代次數(shù)gmax、每個(gè)子區(qū)域容量Vol∈[1,3]、鄰域個(gè)數(shù)T以及子區(qū)域數(shù)目M。

2)進(jìn)行目標(biāo)空間分解操作。

3)計(jì)算粒子群POP目標(biāo)函數(shù)值，選出每個(gè)目標(biāo)函數(shù)的最小值作為參考點(diǎn)R用于計(jì)算每個(gè)粒子的方向向量，方向向量的計(jì)算參考圖1。

4)進(jìn)行粒子的分類與更新操作，接著將當(dāng)前所有子區(qū)域粒子的位置信息存入集合EPOP中，最后更新參考點(diǎn)R，并重新計(jì)算每個(gè)粒子的方向向量，清空POP備用。

5)gbest的選擇：產(chǎn)生[0，1]之間的隨機(jī)數(shù) r1,當(dāng)r1>0.8時(shí)，隨機(jī)從EPOP集合中選擇一個(gè)粒子作為當(dāng)前粒子的gbest；否則，在該粒子所在區(qū)域的鄰域內(nèi)操作。首先計(jì)算該粒子鄰域中每個(gè)子區(qū)域的中心向量和該子區(qū)域中粒子方向向量的余弦值，然后選出余弦值最大的子區(qū)域中的粒子作為當(dāng)前粒子的gbest。過多地使用變異會(huì)讓算法變得更加隨機(jī)，降低了算法的效率，所以產(chǎn)生[0,1]之間的隨機(jī)數(shù)r2,當(dāng)r2<0.6時(shí)，對(duì)選出的引導(dǎo)粒子位置進(jìn)行連續(xù)變異操作，確定最終的引導(dǎo)粒子gbest。否則，不對(duì)gbest進(jìn)行變異操作。

6)pbest的選擇：產(chǎn)生[0,1]之間的隨機(jī)數(shù)r3，當(dāng)r3>0.8時(shí)，隨機(jī)從EPOP中選擇一個(gè)粒子；否則，從該粒子的所在區(qū)域的鄰域中隨機(jī)選擇一個(gè)粒子。讓當(dāng)前粒子和這個(gè)隨機(jī)選出的粒子進(jìn)行比較，選出支配權(quán)優(yōu)先的作為當(dāng)前粒子的pbest。

7)通過上述選出的gbest和pbest，以及粒子速度更新式(5)和位置更新式(6)產(chǎn)生下一代新的粒子群體，這些新粒子位置集合記為NPOP，下一代新的速度集合重新覆蓋集合V。

8)如果當(dāng)前迭代次數(shù)t大于最大迭代次數(shù)gmax，則循環(huán)結(jié)束，輸出EPOP作為最優(yōu)解集，否則把NPOP和EPOP合并放入POP中，然后跳轉(zhuǎn)至3)，繼續(xù)循環(huán)。

3 實(shí)驗(yàn)分析

為了測(cè)試所提出的MOPSO/DC算法的性能，將其與目前較流行的MOPSO/D[7]、MOPSOTL[8]、MOPSODE[9]和NNIA[14]算法進(jìn)行對(duì)比，實(shí)驗(yàn)中同樣選擇了文獻(xiàn)[7]中的6個(gè)測(cè)試函數(shù)，依次為F1，F(xiàn)2···F6。為了定量比較5種算法的性能，采用以下3種廣泛使用的性能指標(biāo)：HV[15]，測(cè)量了粒子收斂于真實(shí)PF的效果以及最優(yōu)解集的多樣性，所獲得值越大越好；GD[16]，指算法獲得的PF到真實(shí)PF的距離，數(shù)值越小越接近真實(shí)PF，效果越好；IGD[15]，顯示真實(shí)PF到算法所得到PF的距離，值越小越能說明算法得到的PF上的點(diǎn)可以分散地均勻地收斂于真實(shí)PF上。

MOPSO/DC參數(shù)設(shè)置：當(dāng)目標(biāo)函數(shù)個(gè)數(shù)m=2時(shí)， M =500 。 當(dāng) m =3 時(shí)， M =1 035。V o l=1，N=M Vol， T =[0.1N]，粒子維度n=10，最大循環(huán)次數(shù)gmax=1 000，一共進(jìn)行了30次實(shí)驗(yàn)求各性能指標(biāo)的均值和標(biāo)準(zhǔn)差。其他算法參數(shù)請(qǐng)參考文獻(xiàn)[7]。

MOPSO/DC與其他算法關(guān)于這6個(gè)測(cè)試函數(shù)的3種性能指標(biāo)對(duì)比結(jié)果見表1，表中給出了通過MOPSO/DC和其他4個(gè)算法獲得關(guān)于這6個(gè)問題的IGD、GD和HV的平均值和標(biāo)準(zhǔn)差。從表1中可以看出：1)MOPSO/DC獲得的IGD平均值遠(yuǎn)小于其他4種算法；2)MOPSO/DC獲得的GD平均值都是最佳的，除F3函數(shù)外，不過關(guān)于F3的GD均值，本算法也優(yōu)于MOPSO/D算法[7]；3)通過MOPSO/DC獲得的HV值總體來說比其他算法獲得的HV值好(除F1函數(shù))，說明MOPSO/DC的收斂性以及解的多樣性也相當(dāng)好。整體來說，MOPSO/DC的這3個(gè)性能指標(biāo)優(yōu)于其他4個(gè)算法，對(duì)MOPSO/D算法進(jìn)行了很好的改進(jìn)。

表 1 各算法效果指標(biāo)對(duì)比Table 1 Comparison of the effect of each algorithm

MOPSO/DC算法的仿真圖如圖3，結(jié)果顯示該算法得到的PF和真實(shí)PF完全重合。

通過這些對(duì)比實(shí)驗(yàn)得出結(jié)論：MOPSO/DC對(duì)原算法[7]的收斂性和多樣性進(jìn)行了很好地優(yōu)化，增強(qiáng)了粒子趨向真實(shí)PF的效果，并且具有很好的穩(wěn)定性。

圖 3 測(cè)試函數(shù)的仿真圖Fig. 3 Test function's simulations

4 結(jié)束語

本文提出的算法采用目標(biāo)空間分解與子區(qū)域粒子更新和對(duì)引導(dǎo)粒子gbest的位置進(jìn)行連續(xù)變異的方法來提高粒子多樣性和收斂性，對(duì)原算法的不足之處進(jìn)行了很好地彌補(bǔ)和優(yōu)化。通過和4個(gè)不同種類的優(yōu)化算法對(duì)不同種測(cè)試函數(shù)進(jìn)行多種性能的測(cè)試對(duì)比實(shí)驗(yàn)證明，MOPSO/DC整體效果最好，具有很高收斂性，多樣性方面也有提升。總之，MOPSO/DC對(duì)多目標(biāo)優(yōu)化問題的多個(gè)性能指標(biāo)同時(shí)進(jìn)行了很好的優(yōu)化。未來的研究重點(diǎn)是，將該算法更好地適應(yīng)更多目標(biāo)的優(yōu)化問題，增強(qiáng)其通用性，用于解決實(shí)際問題。