基于聚類的新型多目標(biāo)分布估計(jì)算法及其應(yīng)用

張秀杰, 高肖霞, 張 虎, 趙 杰

(1. 哈爾濱工業(yè)大學(xué)機(jī)電工程學(xué)院, 黑龍江 哈爾濱 150001; 2. 哈爾濱工業(yè)大學(xué)航天學(xué)院, 黑龍江 哈爾濱 150001)

0 引 言

實(shí)際工程中存在著大量的具有多約束、多變量以及非線性等性質(zhì)的復(fù)雜多目標(biāo)優(yōu)化問(wèn)題(multi-objective optimization problem, MOP)。由于大多數(shù)情況下MOP中各子目標(biāo)之間相互沖突,不存在一個(gè)最優(yōu)解使所有子目標(biāo)同時(shí)達(dá)到最優(yōu)。因此,不同于單目標(biāo)優(yōu)化問(wèn)題只有一個(gè)或者若干個(gè)孤立的最優(yōu)解,MOP具有大量的對(duì)于所有目標(biāo)都可以接受的折衷解,即Pareto最優(yōu)解。所有Pareto最優(yōu)解組成的集合稱為Pareto解集(pareto set, PS),Pareto解集投影到目標(biāo)空間獲得的目標(biāo)向量的集合稱為Pareto前沿(Pareto front, PF)。并且連續(xù)MOP的PS和PF的結(jié)構(gòu)具有規(guī)則特性,即根據(jù)Karush-Kuhn-Tucker條件,在寬松的條件下,具有m個(gè)目標(biāo)的連續(xù)MOP的PS(或PF)的結(jié)構(gòu)是一個(gè)m-1維的分段連續(xù)的流形。對(duì)于一個(gè)MOP,由于不可能求解出其所有的Pareto最優(yōu)解,因此在求解過(guò)程中,決策者往往希望獲得一個(gè)有限數(shù)目的逼近解的集合(逼近解集),其對(duì)應(yīng)的目標(biāo)向量(構(gòu)成逼近前沿)越靠近PF越好(收斂性),并且沿著PF分布越廣泛以及越均勻越好(多樣性)。

由于傳統(tǒng)的確定性優(yōu)化技術(shù)不能較好地對(duì)復(fù)雜的MOP進(jìn)行求解,因此基于自然啟發(fā)搜索的全局優(yōu)化算法——演化算法(evolutionary algorithm, EA)成為了解決MOP流行的方法。多目標(biāo)演化算法(multi-objective evolutionary algorithm, MOEA)具有良好的并行性、魯棒性,而且其求解不依賴于問(wèn)題特性、通用性強(qiáng),并且單次運(yùn)行就可獲得MOP的Pareto解集的逼近,近年來(lái)得到了蓬勃發(fā)展[1]。

在EA當(dāng)中,包含有個(gè)體重組和環(huán)境選擇兩個(gè)重要組成部分。個(gè)體重組用于產(chǎn)生新解,而環(huán)境選擇則負(fù)責(zé)為下一代挑選有效的新解。在MOEA中,根據(jù)新解產(chǎn)生的方式,重組算子可以粗略地劃分為兩大類,即基于遺傳的重組算子和基于模型的重組算子。基于遺傳的MOEA應(yīng)用傳統(tǒng)的重組算子(例如:模擬二進(jìn)制交叉[2]、多項(xiàng)式變異[3]等)產(chǎn)生新解。基于模型的MOEA采用概率模型描述群體中個(gè)體的分布,并通過(guò)建立的模型采樣產(chǎn)生新個(gè)體,多元高斯模型、Bayesian網(wǎng)絡(luò)、流形學(xué)習(xí)、密度估計(jì)等常用的機(jī)器學(xué)習(xí)方法被廣泛應(yīng)用于群體建模[4]。目前已有的MOEA大多數(shù)采用的是基于遺傳的新解產(chǎn)生方法,但是基于模型的MOEA也得到了學(xué)者們?cè)絹?lái)越多的關(guān)注,近些年變得流行的多目標(biāo)分布估計(jì)算法(multi-objective estimation of distribution algorithm, MEDA)就是一個(gè)重要的代表[5]。

分布估計(jì)算法(estimation of distribution algorithm,EDA)[6]是EA中一類特別的方法。EDA并不采用傳統(tǒng)的交叉變異等遺傳操作,取而代之,它從所挑選的有效解中明確地提取全局統(tǒng)計(jì)信息,基于提取的統(tǒng)計(jì)信息建立一個(gè)有效解后驗(yàn)概率分布模型,進(jìn)而從建立的模型中抽樣產(chǎn)生新解。在基于遺傳的MOEA中,遺傳操作可能會(huì)破壞種群強(qiáng)模式的建立,因此種群個(gè)體朝向最優(yōu)解的移動(dòng)方向非常難以預(yù)測(cè)。然而,MEDA能夠預(yù)測(cè)PF的位置或模式,或者預(yù)測(cè)在搜索空間中有效的搜索方向。通過(guò)調(diào)整搜索使其沿著發(fā)掘或預(yù)測(cè)的有效的搜索方向,就能夠較好地產(chǎn)生有效解。學(xué)者們已經(jīng)提出了各種MEDAs,并且這些算法顯示出了良好的性能。

雖然MEDA已經(jīng)得到了越來(lái)越多學(xué)者的關(guān)注和研究,但是依然存在著不足,典型的有:算法中沒(méi)有充分考慮MOP的規(guī)則特性,種群中的異常解處理不正確,種群多樣性容易丟失,以及過(guò)多的計(jì)算開銷用于構(gòu)建最優(yōu)的種群模型[4,7]。針對(duì)上述不足,本文提出一種基于聚類的新型多目標(biāo)分布估計(jì)算法(clustering based MEDA, CEDA)。在CEDA中的每一代,首先利用聚類算法發(fā)掘種群中個(gè)體的分布結(jié)構(gòu),然后基于結(jié)構(gòu)信息,為每一個(gè)個(gè)體構(gòu)建一個(gè)多元高斯模型(multivariate Gaussian model, MGM),基于此模型,抽樣產(chǎn)生新解。

1 多目標(biāo)分布估計(jì)算法

EDA已經(jīng)被大量地應(yīng)用于MOP的求解。文獻(xiàn)[8]提出了一種基于混合的多目標(biāo)迭代密度估計(jì)演化算法(multi-objective mixture-based iterated density estimation evolutionary algorithm, MIEDA),用于求解連續(xù)和離散優(yōu)化問(wèn)題,MIEDA被認(rèn)為是其他MEDAs的基準(zhǔn)算法。文獻(xiàn)[9]采用基于支配的框架并且使用K-means聚類算法進(jìn)行建模,設(shè)計(jì)了一種多目標(biāo)分層Bayesian優(yōu)化算法(multi-objective hierarchical Bayesian optimization algorithm, mohBOA)。文獻(xiàn)[10]提出了一種延伸的緊湊遺傳算法(extended compact genetic algorithm, ECGA)以解決可擴(kuò)展的欺騙問(wèn)題。文獻(xiàn)[11]將EDA集成到基于分解的MOEA框架中,提出了一種混合局部搜索元啟發(fā)式方法的基于分解的EDA。文獻(xiàn)[12]構(gòu)建基于高斯過(guò)程的逆模型(inverse models)將所有已發(fā)掘的非支配解從目標(biāo)空間映射到?jīng)Q策空間,通過(guò)對(duì)目標(biāo)空間抽樣產(chǎn)生新解,提出了基于逆模型的MOEA(inverse models multi-objective evolutionary algorithm,IM-MOEA)。

為了利用連續(xù)MOP的規(guī)則特性提高M(jìn)EDA的求解性能,學(xué)者們又提出了多種基于規(guī)則特性的MEDA。文獻(xiàn)[7]提出了一種基于規(guī)則模型的MEDA,即(regularity model-based multi-objective estimation of distribution algorithm, RM-MEDA),其使用局部主成分分析方法對(duì)有效解建立概率分布模型,并通過(guò)概率模型抽樣產(chǎn)生新個(gè)體。在此之后,又設(shè)計(jì)了一種基于概率模型的MOEA[13],在決策空間和目標(biāo)空間同時(shí)建立概率模型逼近PS和PF。受到RM-MEDA思想的啟發(fā),出現(xiàn)了一系列變種的RM-MEDA,例如,基于冗余類消減的MEDA[14],帶有局部學(xué)習(xí)策略的MEDA[15],基于流形學(xué)習(xí)的演化多目標(biāo)優(yōu)化算法等[16]。

目前為止,在已有的MEDAs中,大多數(shù)在設(shè)計(jì)過(guò)程中并沒(méi)有充分地考慮MOP的規(guī)則特性,并且建模過(guò)程中,只是運(yùn)用較少個(gè)數(shù)的高斯模型描述有效解的分布,往往丟棄了異常解。實(shí)際上,運(yùn)用較少個(gè)數(shù)的高斯模型抽樣產(chǎn)生新解,新解大量的集中在模型的中心位置(均值)附近,多樣性不足,并且異常解可能代表著新的有效區(qū)域,需要開展搜索。另外,在已有的基于規(guī)則特性的MEDAs中,大部分借鑒RM-MEDA的思想,而RM-MEDA建模復(fù)雜,在搜索的早期階段種群的多樣性保持不好,并且難以設(shè)定主成分的個(gè)數(shù)。為了改善前述問(wèn)題,本文充分考慮MOP的規(guī)則特性,計(jì)劃基于聚類運(yùn)用較多數(shù)目的簡(jiǎn)單的高斯模型逼近種群結(jié)構(gòu),進(jìn)而抽樣產(chǎn)生新解,從而降低算法結(jié)構(gòu)的復(fù)雜性,增強(qiáng)算法的易用性,并提高算法產(chǎn)生多樣解的能力。

2 基于聚類的新型多目標(biāo)分布估計(jì)算法

2.1 算法框架

CEDA采用凝聚層次聚類(agglomerative hierarchical clustering, AHC)算法[17]發(fā)掘種群結(jié)構(gòu),CEDA的基本框架如算法1所示。

算法1CEDA基本框架

步驟1初始化種群P={x1,x2,…,xN}和控制概率β,設(shè)置演化代數(shù)t=0。

步驟2主循環(huán)。

步驟2.1設(shè)置一個(gè)空的外部文檔A=?。

步驟2.2對(duì)種群進(jìn)行聚類,

{LC1,…,LCk}=AHC(P,K)

步驟2.3構(gòu)建一個(gè)全局類GC。

步驟2.4分別計(jì)算局部類LCk和全局類GC的協(xié)方差矩陣Σk(k=1,…,K)和ΣGC。

步驟2.5新解產(chǎn)生:對(duì)于每一個(gè)個(gè)體xi∈P,i=1,2,…,N。

步驟2.5.1為個(gè)體xi選擇一個(gè)協(xié)方差矩陣Σi。

步驟2.5.2產(chǎn)生新個(gè)體yi=SolGen(Σi,xi)。

步驟2.5.3保留新解A=A∪{yi}。

步驟2.6環(huán)境選擇:更新種群P=EnvSel(A∪P)。

步驟2.7令t=t+1。

步驟2.8如果t>T,算法結(jié)束,輸出P;否則轉(zhuǎn)向步驟2。

步驟3停機(jī):輸出P。

在算法1中,N代表種群大小,K為AHC中定義的最大聚類個(gè)數(shù),T為最大演化代數(shù),GC和LCk分別代表全局類和第K個(gè)局部類,Σki為xi所在的局部類的協(xié)方差矩陣,β表示利用LCk建立抽樣池的概率(稱之為重組控制概率),rand()生成一個(gè)在[0,1]區(qū)間內(nèi)均勻分布的隨機(jī)數(shù)。

在CEDA算法的每一代中,首先利用AHC將種群劃分為K個(gè)局部類(步驟2.2),并從每一個(gè)局部類中隨機(jī)抽取一個(gè)個(gè)體共同構(gòu)建一個(gè)全局類(步驟2.3)。然后計(jì)算出全局類和所有局部類的協(xié)方差矩陣ΣGC和Σk(k=1,2,…,K)(步驟2.4)。緊接著為每個(gè)個(gè)體xi確定一個(gè)協(xié)方差矩陣Σi,該協(xié)方差矩陣分別以β和1-β的概率設(shè)置為Σk或ΣGC(步驟2.5.1),并且由xi和Σi構(gòu)成高斯模型抽樣產(chǎn)生新個(gè)體yi(步驟2.5.2),將yi存于外部檔案A中(步驟2.5.3)。最后基于A和P利用環(huán)境選擇機(jī)制更新種群P(步驟2.6)。以下小節(jié)對(duì)CEDA的細(xì)節(jié)進(jìn)行詳細(xì)介紹。

2.2 凝聚層次聚類算法

利用AHC算法將種群P中包含的N個(gè)個(gè)體,即P={x1,x2,…,xN},劃分到K個(gè)類中的原理如算法2所示。

算法2AHC(P,K)

步驟1將種群P中每個(gè)個(gè)體看作一個(gè)類。

步驟2循環(huán)。

步驟2.1計(jì)算每?jī)蓚€(gè)不同的類之間的歐氏距離;

步驟2.2找出距離最小的兩個(gè)類合并成新的類;

步驟2.3判斷是否滿足終止條件(聚類個(gè)數(shù)是否大于K),滿足則終止,輸出最終的聚類結(jié)果,否則轉(zhuǎn)至步驟2.1。

AHC首先將每一個(gè)個(gè)體視為一個(gè)類,然后利用一系列機(jī)制合并不同類,直至種群聚類個(gè)數(shù)不大于K。在CEDA的AHC算法中利用組間平均聯(lián)接算法(group average linkage algorithm)定義兩個(gè)類之間的距離。關(guān)于AHC算法的細(xì)節(jié)內(nèi)容參考文獻(xiàn)[17]。

2.3 新解產(chǎn)生

算法1中步驟2.5.2產(chǎn)生新個(gè)體,此過(guò)程包含基于多元高斯模型的抽樣和多項(xiàng)式變異,其細(xì)節(jié)如算法3所示。

算法3SolGen(Σi,xi)

步驟1利用平方根法分解協(xié)方差矩陣Σi得到一個(gè)下三角矩陣Λ,并且Σi=ΛΛT。

步驟2產(chǎn)生向量v=(v1,v2,…,vn)T,其中vj～N(0,I),j=1,2,…,n服從高斯分布。

步驟4修復(fù)該試驗(yàn)解

步驟5對(duì)試驗(yàn)解進(jìn)行變異

y?

aj和bj代表第j個(gè)變量的上下邊界;pm為變異概率,ηm為變異指數(shù),r=rand()。

對(duì)于種群中的每個(gè)個(gè)體,首先基于協(xié)方差矩陣為其建立多元高斯模型并抽樣產(chǎn)生一個(gè)初始試驗(yàn)解(第1～3行)。然后對(duì)試驗(yàn)解進(jìn)行修補(bǔ),保證其可行性(第4行),緊接著對(duì)試驗(yàn)解進(jìn)行變異操作以增強(qiáng)解的多樣性(第5行),最后再次對(duì)試驗(yàn)解進(jìn)行邊界修補(bǔ)確保其可行(第6行)。

2.4 環(huán)境選擇

每一代當(dāng)新解產(chǎn)生完畢之后,運(yùn)用EnvSel(A∪P)從A∪P中選出優(yōu)秀個(gè)體進(jìn)入下一代演化過(guò)程。CEDA采用(S-metric selection based evolutionary multi-objective optimization algorithm, SMS-EMOA)[18]中提出的基于超體積指標(biāo)的環(huán)境選擇方法。超體積指標(biāo)是已知的唯一一個(gè)為“Pareto兼容(Pareto compliant)”的一元指標(biāo)[19],基于超體積指標(biāo)的環(huán)境選擇方法在求解具有復(fù)雜PF的MOP時(shí)展現(xiàn)出了良好的性能[20]。EnvSel(A∪P)的具體過(guò)程如算法4所示。

算法4EnvSel(A∪P)

步驟1利用快速非支配排序方法對(duì)A∪P中的個(gè)體進(jìn)行排序

{B1,B2,…,BL}=Fast_Nondominated_Sort(A∪P)

步驟2復(fù)制較優(yōu)的個(gè)體到輔助種群P′中

步驟3如果l>1

循環(huán):當(dāng)|P′|>N

步驟3.2將x*從P′中移除,P′=P′{x*}。

步驟4如果l=1

循環(huán):當(dāng)|P′|>N

步驟4.2將x*從P′中移除,P′=P′{x*}。

步驟5將P′賦給P,P=P′。

步驟6輸出P。

首先將當(dāng)前種群P和外部文檔A合并成一個(gè)新的種群,并且利用(nondominated sorting genetic algorithm Ⅱ, NSGA-Ⅱ)[21]中提出的快速非支配排序方法將新種群中的個(gè)體劃分到L個(gè)不同的非支配前沿{B1,B2,…,BL}當(dāng)中。然后根據(jù)排序的結(jié)果,將新種群中的較優(yōu)個(gè)體復(fù)制到一個(gè)輔助種群P′當(dāng)中,直到|P′|≥N。如果P′當(dāng)中包含有多個(gè)非支配前沿(即l>1),則逐個(gè)移除第l前沿中d(x,P′)最大的個(gè)體,直到|P′|=N,d(x,P′)指的是在P′中支配x的點(diǎn)的個(gè)數(shù);否則如果l=1,逐個(gè)移除P′中超體積貢獻(xiàn)度Δφ(x,P′)最小的個(gè)體,直到|P′|=N,超體積貢獻(xiàn)度Δφ的計(jì)算方法參考文獻(xiàn)[18]。最后,將P′賦給P,作為下一代的種群。

2.5 算法說(shuō)明

關(guān)于CEDA,有如下說(shuō)明:

(1) 文獻(xiàn)[22]指出在MOEA中,相似個(gè)體重組,能提高產(chǎn)生新解的質(zhì)量。產(chǎn)生此效果的原因是增強(qiáng)了算法的局部搜索,暗含地利用了MOP的規(guī)則特性。與此類似,本文的CEDA采用鄰近個(gè)體為每個(gè)當(dāng)前個(gè)體構(gòu)建高斯模型逼近種群結(jié)構(gòu)進(jìn)而抽樣產(chǎn)生新解,這一機(jī)制同樣增強(qiáng)了算法的局部搜索,充分地考慮了MOP的規(guī)則特性,理應(yīng)也能更好地產(chǎn)生高質(zhì)量的新解。

(2) 與RM-MEDA中利用局部主成分分析方法提取流形結(jié)構(gòu),然后抽樣產(chǎn)生新解的方式相比,CEDA中的基于聚類建立高斯模型抽樣產(chǎn)生新解的方式更簡(jiǎn)單易用。并且,在進(jìn)化的早期階段,PS的流形結(jié)構(gòu)尚未顯現(xiàn),種群需要良好多樣性,但是RM-MEDA的新解產(chǎn)生方式限制了新解的產(chǎn)生方向,不利于產(chǎn)生多樣化的解,而CEDA利用完全協(xié)方差矩陣抽樣產(chǎn)生新解,能從各個(gè)方向生成新解,更好地維護(hù)新解的多樣性。

(3) MIEDA[8]等傳統(tǒng)的為每一個(gè)類構(gòu)建一個(gè)高斯模型進(jìn)行抽樣的方式,其產(chǎn)生的新解大量地分布在均值向量附近,新解的多樣性不夠,而CEDA為每個(gè)種群個(gè)體以自身為均值建立一個(gè)高斯模型抽樣產(chǎn)生新解,實(shí)際上是為每個(gè)個(gè)體添加一個(gè)高斯擾動(dòng),此種方式能產(chǎn)生更為多樣化的解。

(4) 在為個(gè)體構(gòu)建高斯模型的時(shí)候,如果對(duì)于每個(gè)個(gè)體都計(jì)算協(xié)方差矩陣,則需要大量的建模計(jì)算開銷。為了解決此問(wèn)題,CEDA中使同一類中的個(gè)體共享相同的協(xié)方差矩陣進(jìn)行建模從而大大降低建模計(jì)算開銷。之所以能夠進(jìn)行此策略是因?yàn)橄嗨频膫€(gè)體理應(yīng)具有相近的高斯模型,并且近似的高斯模型就已經(jīng)滿足算法的要求,沒(méi)有必要花費(fèi)大量的計(jì)算開銷建立精確的模型。

(5) 不同于以往丟棄異常解的建模方式,CEDA中為每個(gè)個(gè)體建立一個(gè)高斯模型進(jìn)行抽樣,實(shí)際上就是充分地考慮了異常解代表的解空間區(qū)域,因此能更好地加強(qiáng)對(duì)解空間的搜索。

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

3.1 標(biāo)準(zhǔn)測(cè)試實(shí)例和性能度量指標(biāo)

為了測(cè)試CEDA的性能,首先利用標(biāo)準(zhǔn)測(cè)試題對(duì)其進(jìn)行測(cè)試。大多數(shù)工程中的MOP具有復(fù)雜的PF結(jié)構(gòu),因此CEDA算法理應(yīng)適用于求解此類具有復(fù)雜PF結(jié)構(gòu)的MOPs。本文利用具有復(fù)雜PF和PS結(jié)構(gòu)的6道標(biāo)準(zhǔn)測(cè)試題GLT1～GLT6對(duì)CEDA進(jìn)行測(cè)試。其中,GLT1～GLT4為雙目標(biāo)測(cè)試問(wèn)題,GLT5～GLT6為三目標(biāo)測(cè)試問(wèn)題。(Gu F, Liu H, Tan K, GLT)測(cè)試題的具體細(xì)節(jié)參考文獻(xiàn)[20]。

為了評(píng)估算法的性能,運(yùn)用兩個(gè)常用的性能指標(biāo),即反世代距離(inverted generational distance, IGD)[7,23]和超體積(hypervolume, HV)[24]度量算法獲得的逼近前沿的效果。IGD和HV是兩個(gè)均能夠綜合評(píng)價(jià)獲得的逼近前沿的收斂性與多樣性的性能指標(biāo)。并且IGD值越小，HV值越大，代表算法所求得的逼近前沿的收斂性和多樣性越好。

在接下來(lái)的實(shí)驗(yàn)中,計(jì)算HV指標(biāo)值時(shí),各測(cè)試題的參考點(diǎn)取值為:GLT1取r=(2,2)T,GLT2取r=(2,11)T,GLT3取r=(2,2)T,GLT4取r=(2,3)T,GLT5-GLT6取r=(2,2,2)T。

3.2 對(duì)比算法及參數(shù)設(shè)置

選取4種典型的MOEAs,即NSGA-Ⅱ[21]、SMS-EMOA[18]、RM-MEDA[7]和(transform each individual objective function into the one multi-objective evolutionary algorithm based on decomposition, TMOEA/D)[25],與CEDA開展對(duì)比實(shí)驗(yàn)。NSGA-Ⅱ是一種基于支配的MOEA,SMS-EMOA是一種基于指標(biāo)的MOEA,TMOEA/D是一種用于解決具有復(fù)雜PF形狀的MOPs的基于分解的MOEA,RM-MEDA是一種基于規(guī)則特性的MEDA,這幾種算法涵蓋了當(dāng)前主流的MOEA類型。為了保證對(duì)比的公平性,所有對(duì)比算法的參數(shù)均通過(guò)前期的實(shí)驗(yàn)進(jìn)行了系統(tǒng)地優(yōu)化,對(duì)比實(shí)驗(yàn)中采用最佳的參數(shù)組合。所有算法均由Matlab進(jìn)行實(shí)現(xiàn),并且在同一臺(tái)計(jì)算機(jī)上運(yùn)行,具體的算法參數(shù)設(shè)置如下:

(1) 公共參數(shù)

(2) NSGA-Ⅱ參數(shù)

模擬二進(jìn)制交叉參數(shù)Pc=0.9,ηc=20;多項(xiàng)式變異算子控制參數(shù)pm=1/n,ηm=20。

(3) SMS-EMOA參數(shù)

模擬二進(jìn)制交叉參數(shù)Pc=0.9,ηc=20;多項(xiàng)式變異算子控制參數(shù)pm=1/n,ηm=20。

(4) RM-MEDA參數(shù)

聚類個(gè)數(shù)(principal component analysis, PCA)=5;局部主成分分析最大迭代次數(shù)為50;擴(kuò)展采樣率為0.25。

(5) TMOEA/D參數(shù)

鄰居大小NS=30;第一搜索階段演化代數(shù)T1=T/10;第二搜索階段演化代數(shù)T2=αT,α={0.01,0.02,…,0.1,0.1,0.1,0.15};差分演化交叉算子控制參數(shù)F=0.5,CR=1。

(6) CEDA參數(shù)

重組控制概率β=0.9;最大聚類數(shù)目K=5;多項(xiàng)式變異算子控制參數(shù)pm=1/n,ηm=20。

為了在實(shí)驗(yàn)中獲得統(tǒng)計(jì)置信的結(jié)論,每種算法對(duì)每道測(cè)試題獨(dú)立運(yùn)行33次,基于統(tǒng)計(jì)指標(biāo)值(均值和標(biāo)準(zhǔn)差)進(jìn)行算法性能的比較。在比較表格中,關(guān)于某一道測(cè)試題,各算法對(duì)其統(tǒng)計(jì)運(yùn)算獲得的指標(biāo)值的均值進(jìn)行升序(IGD指標(biāo))或降序(HV指標(biāo))排序,排序結(jié)果在表格的方括號(hào)中顯示,并且每種算法對(duì)GLT測(cè)試集的計(jì)算性能排序的平均值(平均秩)也列在表格中。對(duì)于每道測(cè)試題,各算法獲得的平均指標(biāo)值中最優(yōu)值用深灰色背景表示,次優(yōu)值用淺灰色背景表示。另外,當(dāng)CEDA與任一算法進(jìn)行比較的時(shí)候,執(zhí)行在5%顯著性水平的Wilcoxon秩和檢驗(yàn)觀測(cè)差異的顯著性。“?”,“§”和“≈”表示CEDA求解某問(wèn)題的性能在5%顯著水平上是優(yōu)于、劣于以及相似于比較算法對(duì)于該問(wèn)題的求解能力。

3.3 實(shí)驗(yàn)結(jié)果

3.3.1 質(zhì)量指標(biāo)

表1給出了NSGA-Ⅱ、SMS-EMOA、RM-MEDA、TMOEA/D以及CEDA算法分別獨(dú)立計(jì)算GLT測(cè)試集33次各自獲得的逼近前沿的HV和IGD值的統(tǒng)計(jì)結(jié)果。

表1 NSGA-Ⅱ、SMS-EMOA、RM-MEDA、TMOEA/D以及CEDA分別獨(dú)立計(jì)算GLT測(cè)試集33次所得的逼近前沿的IGD和HV指標(biāo)值的統(tǒng)計(jì)結(jié)果(均值(標(biāo)準(zhǔn)差)[排序])

從表1中可以看出,通過(guò)演化300代,與對(duì)比算法進(jìn)行比較,在12個(gè)指標(biāo)值中,CEDA獲得了8個(gè)最優(yōu)和2個(gè)次優(yōu)平均指標(biāo)值。根據(jù)統(tǒng)計(jì)顯著性檢驗(yàn),相對(duì)于NSGA-Ⅱ、SMS-EMOA、RM-MEDA和TMOEA/D,在與每種算法的12次比較中,CEDA分別獲得了12、11、10和7個(gè)明顯較優(yōu)的平均指標(biāo)值。另外,平均秩的值表明在求解GLT測(cè)試集時(shí),性能從最優(yōu)到最劣的算法分別是CEDA、TMOEA/D、RM-MEDA、SMS-EMOA、NSGA-Ⅱ。

3.3.2 搜索效率

圖1繪制出了NSGA-Ⅱ、SMS-EMOA、RM-MEDA、TMOEA/D以及CEDA算法分別獨(dú)立計(jì)算GLT測(cè)試集33次各自獲得的平均IGD指標(biāo)值的演化曲線。從圖中可以看出對(duì)于GLT2-GLT3和GLT5-GLT6的求解中,CEDA均在最少的演化代數(shù)內(nèi)獲得了最低的平均IGD指標(biāo)值。對(duì)于GLT1,CEDA的求解性能劣于RM-MEDA和TMOEA/D。對(duì)于GLT4,CEDA劣于TMOEA/D獲得了次優(yōu)的求解效果?？傮w而言,與其他4種相比,CEDA在求解GLT測(cè)試集的演化過(guò)程中收斂速度最快并且能夠維持最好的種群多樣性。

圖1 平均IGD指標(biāo)值演化曲線Fig.1 Evolution of the mean IGD metric values

3.3.3 結(jié)果可視化

圖2繪制出了統(tǒng)計(jì)比較中性能最好的兩種算法CEDA和TMOEA/D分別獨(dú)立計(jì)算GLT測(cè)試集33次各自獲得的全部最終逼近前沿(如圖2(a)和2(b)),以及分別獲得中位IGD指標(biāo)值時(shí)對(duì)應(yīng)的代表性逼近前沿(如圖2(c)和2(d))。從圖2(a)和2(b)可以看出,在33次獨(dú)立運(yùn)算中,在求解GLT1和GLT4時(shí),CEDA獲得的逼近前沿還有部分并未收斂到PFs上,但是在求解GLT2、GLT3、GLT5、GLT6時(shí),CEDA獲得的全部逼近前沿都能穩(wěn)定的收斂到PFs上,而且覆蓋整個(gè)PFs。然而,TMOEA/D求解GLT4獲得的逼近前沿并沒(méi)有全部收斂到PF上,并且求解GLT5和GLT6時(shí),其獲得的逼近前沿并未完全覆蓋整個(gè)的PFs。從圖2(c)和2(d)可以觀察到,TMOEA/D求解GLT3和GLT4時(shí),獲得的代表性的前沿雖然最終都能收斂到PFs,但并不能完全覆蓋PFs,求解GLT5和GLT6時(shí),得到的代表性前沿中仍存在一些個(gè)體沒(méi)有完全收斂到PFs,且前沿分布的均勻性也并不理想。與TMOEA/D相比,CEDA對(duì)于GLT2～GLT6獲得的代表性的前沿均具有更好的收斂性和多樣性。

圖2 TMOEA/D和CEDA獲得的逼近前沿Fig.2 Approximated fronts obtained by TMOEA/D and CEDA

根據(jù)上述的質(zhì)量指標(biāo)、搜索效率以及結(jié)果可視化,可以推斷出相對(duì)于NSGA-Ⅱ、SMS-EMOA、RM-MEDA和TMOEA/D,CEDA算法對(duì)于GLT測(cè)試集具有最佳的求解性能。

3.4 參數(shù)靈敏度分析

3.4.1 交配限制概率

在CEDA中,采用重組控制概率β維護(hù)算法演化過(guò)程中勘探與開發(fā)之間的平衡。為了分析β對(duì)算法性能的影響,采用不同β值(β=0.5,0.6,0.7,0.8,0.9)構(gòu)造CEDA算法對(duì)GLT測(cè)試集進(jìn)行求解,算法的其他參數(shù)與第3.2節(jié)設(shè)置相同。每種帶有不同β值的算法對(duì)每道測(cè)試題進(jìn)行22次獨(dú)立運(yùn)算,獲得的逼近前沿的IGD指標(biāo)值的均值和標(biāo)準(zhǔn)差如圖3所示。

圖3 重組控制概率(β)分析Fig.3 Mating limit probability (β) analysis

由圖3可以看出,當(dāng)求解GLT1,GLT3和GLT4時(shí),不同的β值獲得的平均IGD值有明顯的差異,然而對(duì)其他測(cè)試題求解時(shí),不同的β值卻得到了相似的平均IGD值。但是總體來(lái)說(shuō),當(dāng)β=0.9時(shí),CEDA對(duì)于GLT1-GLT3以及GLT5-GLT6均具有較好的求解效果,因此說(shuō)明算法的性能對(duì)于β的取值并不十分敏感。

3.4.2 聚類數(shù)目

CEDA中采用AHC方法發(fā)掘種群結(jié)構(gòu)。為了分析AHC中的最大聚類數(shù)目K對(duì)CEDA性能的影響,采用不同K值(K=4,5,7,10,20)構(gòu)造CEDA算法對(duì)GLT測(cè)試集進(jìn)行求解,算法中的其他參數(shù)與第3.2節(jié)設(shè)置相同。每種帶有不同K值的算法對(duì)每道測(cè)試題進(jìn)行22次獨(dú)立運(yùn)算,獲得的逼近前沿的IGD指標(biāo)值的均值和標(biāo)準(zhǔn)差如圖4所示。

圖4 聚類數(shù)目(K)分析Fig.4 Number of cluster (K) analysis

由圖4可以看出,當(dāng)求解GLT1～GLT4時(shí),不同K值的CEDA獲得的平均IGD值有顯著的差異,然而對(duì)于GLT5～GLT6求解時(shí),不同的K值卻得到了相近的平均IGD值?？傮w來(lái)說(shuō),當(dāng)K=5,CEDA對(duì)于不同的測(cè)試題均能獲得較小的平均IGD值,因此說(shuō)明CEDA的性能對(duì)于最大聚類數(shù)目的取值也并不十分敏感。

4 工程應(yīng)用

4.1 優(yōu)化模型

齒輪減速器是原動(dòng)機(jī)和工作機(jī)之間獨(dú)立的閉式傳動(dòng)裝置,用來(lái)降低轉(zhuǎn)速和增大轉(zhuǎn)矩,以滿足工作需要。其無(wú)須聯(lián)軸器和適配器,結(jié)構(gòu)緊湊。負(fù)載分布在行星齒輪上,因而承載能力比一般斜齒輪減速機(jī)高,滿足小空間高扭矩輸出的需要,廣泛應(yīng)用于大型礦山、鋼鐵、化工、港口、環(huán)保等領(lǐng)域。雖然齒輪減速器應(yīng)用廣泛,但長(zhǎng)期以來(lái)減速器的設(shè)計(jì)僅由設(shè)計(jì)人員依靠相關(guān)的資料、文獻(xiàn),以及多年的經(jīng)驗(yàn)完成,因而不僅效率低,而且可能造成人力、物力和財(cái)力的浪費(fèi),因此目前需要找到一種快速有效的方法來(lái)優(yōu)化設(shè)計(jì)齒輪減速器。齒輪減速器的優(yōu)化設(shè)計(jì)實(shí)際上是一個(gè)多目標(biāo)、多約束、多耦合的優(yōu)化問(wèn)題,普通的算法難以較好地求解此問(wèn)題。本文以此問(wèn)題為實(shí)例檢驗(yàn)CEDA求解復(fù)雜工程優(yōu)化問(wèn)題的效果。文獻(xiàn)[26-27]中建立了某輕型飛機(jī)的齒輪減速器簡(jiǎn)易模型如圖5所示。

該MOP的設(shè)計(jì)目標(biāo)是使減速器的體積和軸2所承受的應(yīng)力最小,并且滿足輪齒的彎曲應(yīng)力、接觸應(yīng)力、軸的扭轉(zhuǎn)變形以及應(yīng)力等約束。該問(wèn)題的數(shù)學(xué)模型描述為

圖5 齒輪減速器結(jié)構(gòu)簡(jiǎn)圖Fig.5 Gear Reducer structure diagram

(1)

(2)

式中,x1為齒寬;x2為齒輪模數(shù);x3為小齒輪齒數(shù);x4為軸承1之間的距離;x5為軸承2之間的距離;x6為軸1的直徑;x7為軸2的直徑;g1為齒的彎曲應(yīng)力約束;g2為齒的接觸應(yīng)力約束;g3、g4為軸的變形約束;g5、g6、g7為基于空間尺寸限制和經(jīng)驗(yàn)約束;g8、g9為由經(jīng)驗(yàn)得到的設(shè)計(jì)軸的要求;g10、g11為軸的應(yīng)力約束;g12至g25為7個(gè)變量的上下邊界。

4.2 實(shí)驗(yàn)設(shè)計(jì)與結(jié)果分析

利用NSGA-Ⅱ、SMS-EMOA、RM-MEDA、TMOEA/D以及CEDA對(duì)齒輪減速器優(yōu)化設(shè)計(jì)模型進(jìn)行求解。經(jīng)過(guò)參數(shù)優(yōu)化,運(yùn)算過(guò)程中的參數(shù)設(shè)置如表2所示,其余設(shè)計(jì)與第3.2節(jié)相同。每種算法對(duì)模型獨(dú)立運(yùn)算33遍,利用超體積HV指標(biāo)值度量每一次獲得的逼近前沿的效果。其中,計(jì)算HV值時(shí)取參考點(diǎn)r=[6 600,1 600]T。

5種算法對(duì)齒輪減速器優(yōu)化設(shè)計(jì)模型獨(dú)立運(yùn)行33次獲得的HV指標(biāo)值的箱型圖對(duì)比結(jié)果如圖6所示(圖6(a)為原圖,圖6(b)為局部放大圖)。從圖6中可以看出CEDA獲得了最大的中位HV指標(biāo)值和最小的四分位距,從而說(shuō)明CEDA對(duì)于齒輪減速器優(yōu)化設(shè)計(jì)模型能穩(wěn)定地求解出具有良好多樣性和收斂性的解。

表2 NSGA-Ⅱ、SMS-EMOA、RM-MEDA、TMOEA/D及CEDA算法 求解齒輪減速器優(yōu)化設(shè)計(jì)模型時(shí)參數(shù)設(shè)置

圖7(圖7(a)為原圖,圖7(b)為局部放大圖)繪制出了NSGA-Ⅱ、SMS-EMOA、RM-MEDA、TMOEA/D以及CEDA算法分別獨(dú)立計(jì)算齒輪減速器優(yōu)化設(shè)計(jì)模型33次各自獲得的平均HV指標(biāo)值的演化曲線。從圖7中可以看出CEDA在最小的演化代數(shù)內(nèi)獲得了最高的平均HV指標(biāo)值。也就是說(shuō)與其他4種相比,CEDA在演化過(guò)程中收斂速度最快并且能夠維持最好的種群多樣性。

圖8為分別利用NSGA-Ⅱ和CEDA對(duì)齒輪減速器優(yōu)化設(shè)計(jì)模型求解時(shí),獨(dú)立運(yùn)算33次各自獲得的全部逼近前沿以及中位IGD指標(biāo)值對(duì)應(yīng)的代表性逼近前沿。從圖8(a)和圖8(c)可以看出,CEDA獲得的全部逼近前沿均能夠穩(wěn)定地收斂,并且與NSGA-Ⅱ獲得的逼近前沿相比,其覆蓋面更廣。從圖8(b)和圖8(d)可以看出,相對(duì)于NSGA-Ⅱ,CEDA獲得了更為寬廣和均勻的代表性的逼近前沿。

從圖6～圖8的分析中可以推斷出CEDA算法對(duì)于齒輪減速器優(yōu)化設(shè)計(jì)模型具有優(yōu)異的求解性能。

圖6 NSGA-Ⅱ、SMS-EMOA、RM-MEDA、TMOEA/D以及CEDA算法對(duì)齒輪減速器優(yōu)化設(shè)計(jì)模型分別獨(dú)立運(yùn)算33次獲得的HV指標(biāo)值的箱線圖 Fig.6 Boxplot of HV metric values obtained by NSGA-П、SMS-EMOA、RM-MEDA、TMOEA/D and CEDA, respectively,over 33 independent runs on the application model

圖7 平均HV指標(biāo)值演化曲線Fig.7 Evolution of the mean HV metric values

圖8 NSGA-Ⅱ和CEDA獲得的逼近前沿Fig.8 Approximated fronts obtained by NSGA-Ⅱ and CEDA

5 結(jié) 論

本文設(shè)計(jì)了一種CEDA算法。在CEDA中,首先利用凝聚層次聚類算法將種群劃分為若干個(gè)局部類,從每一個(gè)局部類中隨機(jī)選擇一個(gè)個(gè)體構(gòu)成一個(gè)全局類,然后為每個(gè)個(gè)體構(gòu)建一個(gè)高斯模型(此高斯模型的均值為個(gè)體本身,協(xié)方差矩陣為個(gè)體所在局部類的協(xié)方差矩陣或者是全局類的協(xié)方差矩陣)去逼近種群結(jié)構(gòu),并抽樣產(chǎn)生新個(gè)體。此新解產(chǎn)生方法充分地考慮了多目標(biāo)優(yōu)化問(wèn)題的規(guī)則特性,其本質(zhì)是為每個(gè)個(gè)體添加了一個(gè)外部擾動(dòng),能改善已有的大部分分布估計(jì)算法中存在的異常個(gè)體處理不合理,種群多樣性容易丟失的問(wèn)題。并且處于同類中的個(gè)體共享協(xié)方差矩陣用于建模,極大地降低了建模的計(jì)算開銷。

以具有復(fù)雜Pareto前沿和復(fù)雜Pareto解集形狀的多目標(biāo)優(yōu)化測(cè)試題為求解對(duì)象,將CEDA與典型的MOEAs進(jìn)行了對(duì)比實(shí)驗(yàn)。實(shí)驗(yàn)結(jié)果表明,CEDA對(duì)于此類問(wèn)題具有優(yōu)異的求解性能。將CEDA算法應(yīng)用于齒輪減速器優(yōu)化設(shè)計(jì)中,結(jié)果表明,CEDA算法同樣能夠快速有效地求解此類復(fù)雜的實(shí)際工程問(wèn)題。

[1] ZHOU A, QU B Y, LI H, et al. Multiobjective evolutionary algorithms: a survey of the state of the art[J]. Swarm & Evolutionary Computation, 2011, 1(1):32-49.

[2] DEB K, BEYER H G. Self-adaptive genetic algorithms with simulated binary crossover[J].Evolutionary Computation,2001,9(2):197-221.

[3] SCHAER J D. Multiple objective optimization with vector evaluated genetic algorithms[C]∥Proc.of the International Conference on Genetic Algorithms and Their Applications, 1985: 93-100.

[4] MARTí L, GRIMME C, KERSCHKE P, et al. Averaged Hausdorff approximations of pareto fronts based on multiobjective estimation of distribution algorithms[C]∥Proc.of the Companion Publication of the Annual Conference on Genetic and Evolutionary Computation, 2015:1427-1428.

[5] PELIKAN M, SASTRY K, GOLDBERG D E. Multiobjective estimation of distribution algorithms[C]∥Proc.of the Scalable Optimization Via Probabilistic Modeling, 2006:223-248.

[7] ZHANG Q, ZHOU A, JIN Y. RM-MEDA: a regularity model-based multiobjective estimation of distribution algorithm[J]. IEEE Trans.on Evolutionary Computation, 2008, 12(1): 41-63.

[8] BOSMAN P A, THIERENS D. Multi-objective optimization with diversity preserving mixture-based iterated density estimation evolutionary algorithms[J]. International Journal of Approximate Reasoning, 2002, 31(3): 259-289.

[9] PELIKAN M, SASTRY K, GOLDBERG D E. Multiobjective hBOA, clustering, and scalability[C]∥Proc.of the 7th Annual Conference on Genetic and Evolutionary Computation, 2005: 663-670.

[10] SASTRY K, GOLDBERG D E, PELIKAN M. Limits of scalability of multiobjective estimation of distribution algorithms[C]∥Proc.of the IEEE Congress on Evolutionary Computation, 2005: 2217-2224.

[11] SHIM V A, TAN K C, CHEONG C Y. A hybrid estimation of distribution algorithm with decomposition for solving the multiobjective multiple traveling salesman problem[J]. IEEE Trans.on Systems, Man, and Cybernetics, Part C(Applications and Reviews), 2012, 42(5): 682-691.

[12] CHENG R, JIN Y, NARUKAWA K, et al. A multiobjective evolutionary algorithm using Gaussian process-based inverse modeling[J]. IEEE Trans.on Evolutionary Computation, 2015, 19(6): 838-856.

[13] ZHOU A, ZHANG Q, JIN Y. Approximating the set of Pareto-optimal solutions in both the decision and objective spaces by an estimation of distribution algorithm[J]. IEEE Trans.on Evolutionary Computation, 2009, 13(5): 1167-1189.

[14] WANG Y, XIANG J, CAI Z. A regularity model-based multiobjective estimation of distribution algorithm with reducing redundant cluster operator[J]. Applied Soft Computing, 2012, 12(11): 3526-3538.

[15] LI Y, XU X, LI P, et al. Improved RM-MEDA with local learning[J]. Soft Computing, 2014, 18(7): 1383-1397.

[16] LI K, KWONG S. A general framework for evolutionary multiobjective optimization via manifold learning[J]. Neurocomputing, 2014, 146(1): 65-74.

[17] XU R,WUNSCH D. Clustering[M]. New Jersey: Wiley,2008.

[18] BEUME N, NAUJOKS B, EMMERICH M. SMS-EMOA: multiobjective selection based on dominated hypervolume[J]. European Journal of Operational Research, 2007, 181(3): 1653-1669.

[19] ZITZLER E, THIELE L, LAUMANNS M, et al. Performance assessment of multiobjective optimizers: an analysis and review[J]. IEEE Trans.on Evolutionary Computation, 2003, 7(2): 117-132.

[20] ZHANG H, ZHOU A, SONG S, et al. A self-organizing multiobjective evolutionary algorithm[J], IEEE Trans.on Evolutionary Computation, 2016, 20(5): 792-806.

[21] DEB K, PRATAP A, AGARWAL S, et al. A fast and elitist multiobjective genetic algorithm: NSGA-Ⅱ[J]. IEEE Trans.on Evolutionary Computation, 2002, 6(2): 182-197.

[22] JIN Y, SENDHOFF B. Connectedness, regularity and the success of local search in evolutionary multi-objective optimization[C]∥Proc.of the IEEE Congress on Evolutionary Computation, 2003: 1910-1917.

[23] ZHOU A, ZHANG Q, JIN Y, et al. A model-based evolutionary algorithm for bi-objective optimization[C]∥Proc.of the IEEE Congress on Evolutionary Computation, 2005: 2568-2575.

[24] ZITZLER E, THIELE L. Multiobjective evolutionary algorithms: a comparative case study and the strength Pareto approach[J].IEEE Trans.on Evolutionary Computation,1999,3(4):257-271.

[25] LIU H L, GU F, CHEUNG Y. T-MOEA/D: MOEA/D with objective transform in multi-objective problems[C]∥Proc.of the International Conference on Information Science and Management Engineering, 2010: 282-285.

[26] FARHANG-MEHR A, AZARM S. Entropy-based multi-objective genetic algorithm for design optimization[J]. Structural & Multidisciplinary Optimization, 2002, 24(5): 351-361.

[27] AZARM S, TITS A, FAN M. Tradeoff driven optimization-based design of mechanical systems[C]∥Proc.of the Symposium on Multidisciplinary Analysis and Optimization, 1989: 551-558.