一種基于Kriging模型加點(diǎn)策略的多目標(biāo)粒子群優(yōu)化算法

陳 靜, 唐傲天, 劉 震, 徐 森, 曹曉聰

(1. 吉林大學(xué) 汽車工程學(xué)院, 長(zhǎng)春 130022; 2. 一汽轎車銷售有限公司, 長(zhǎng)春 130013)

在工程實(shí)踐和研究中存在許多使多個(gè)目標(biāo)在設(shè)計(jì)區(qū)域上盡可能最優(yōu)的多目標(biāo)優(yōu)化問題(multi objective optimization problems, MOP), 以及求解或?qū)嶒?yàn)成本很高的多目標(biāo)優(yōu)化問題(expensive multi objective optimization problems, EMOPs). 近年來, 由于粒子群優(yōu)化算法(particle swarm optimization, PSO)在多目標(biāo)優(yōu)化問題中具有協(xié)調(diào)多個(gè)目標(biāo)的作用, 因此多目標(biāo)粒子群優(yōu)化(multi objective particle swarm optimization, MOPSO)被廣泛應(yīng)用于工程實(shí)踐中[1], 目前已取得了許多研究成果. 孫小強(qiáng)等[2]采用聚類算法裁剪非支配解, 提高了解的分布性； 郭玉潔等[3]提出了一種雙種群協(xié)同多目標(biāo)粒子群優(yōu)化算法, 通過雙種群協(xié)同進(jìn)化策略擴(kuò)大搜索空間, 提高算法的全局搜索能力, 并結(jié)合Lévy飛行保證了種群多樣性； 韓紅桂等[4]提出了一種可提高多目標(biāo)粒子群算法優(yōu)化解的多樣性和收斂性的MOPSO, 提高了計(jì)算速率. 但上述算法在收斂性方面, 特別是解決高維多目標(biāo)問題或模型求解復(fù)雜問題上還有待改進(jìn). 劉華鎣等[5]采用網(wǎng)格擁擠距離與網(wǎng)格密度相結(jié)合的策略選取全局極值, 提出了一種基于空間劃分樹的多目標(biāo)粒子群優(yōu)化算法, 在保持種群多樣性的同時(shí)提高了計(jì)算結(jié)果的準(zhǔn)確性, 但該方法的估值仍有一定誤差. 針對(duì)以上問題, 本文將多目標(biāo)粒子群優(yōu)化算法進(jìn)行改進(jìn), 針對(duì)MOPSO提出一種基于Kriging模型的變異模式和基于Kriging模型的加點(diǎn)策略, 并將改進(jìn)算法運(yùn)用到多目標(biāo)的Pareto最優(yōu)解測(cè)試中, 取得了較好的效果.

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

1.1 多目標(biāo)優(yōu)化問題的數(shù)學(xué)描述

當(dāng)優(yōu)化問題有m個(gè)目標(biāo)函數(shù)、r個(gè)約束函數(shù)時(shí), 其數(shù)學(xué)模型可表示為

(1)

其中：m>1為目標(biāo)空間的維數(shù)；x=(x1,x2,…,xo)為o維設(shè)計(jì)變量；Ωo為x的可行解空間；f(x)為目標(biāo)函數(shù)；gj(x)≤0為不等式約束條件. 由于MOP中各目標(biāo)存在沖突關(guān)系, 無法使所有目標(biāo)都實(shí)現(xiàn)最優(yōu), 因此只能協(xié)調(diào)考慮各目標(biāo), 尋求折中方案找到綜合結(jié)果最優(yōu)的解.

定義1(Pareto支配)[6]對(duì)于式(1), 假設(shè)有設(shè)計(jì)變量x(1),x(2), 滿足fi(x(1))≤fi(x(2))(i=1,2,…,m), 且?k∈{1,2…,m}, 使得fk(x(1))

定義2(Pareto最優(yōu)解)[6]對(duì)于式(1), 若不存在滿足x(3)x(4)的設(shè)計(jì)變量x(3), 則稱x(4)為非支配解或Pareto最優(yōu)解.

定義3(Pareto最優(yōu)解集)[6]可行解空間Ωo中所有Pareto最優(yōu)解構(gòu)成的集合稱為Pareto最優(yōu)解集, 記為Ps={x(5)|?x(6)x(5)}.

定義4(Pareto前沿)[6]Pareto最優(yōu)解集中所有Pareto最優(yōu)解目標(biāo)函數(shù)值形成的曲線(面)稱為Pareto前沿, 記為PF={f(x)|x∈Ps}.

1.2 Kriging模型

Kriging[7]是一種基于高斯過程理論的元建模技術(shù), 其不僅能精準(zhǔn)預(yù)測(cè)未試驗(yàn)點(diǎn)的仿真響應(yīng)值, 還能給出預(yù)測(cè)不確定性的度量, 在高維度、 非線性數(shù)學(xué)擬合方面表現(xiàn)優(yōu)異, 對(duì)黑箱問題的解決有突出優(yōu)勢(shì)[8]. 設(shè)N次取樣得到樣本點(diǎn){x(1),x(2),…,x(N)}的仿真結(jié)果為{y(1),y(2),…,y(N)}, 則Kriging模型可由下式表示：

(2)

(3)

其中:o為設(shè)計(jì)變量x的個(gè)數(shù);θk,pk為待定系數(shù),k=1,2,…,o, 其取值分別決定Ωo第o維的權(quán)重和平順程度. 應(yīng)用隨機(jī)過程理論可得Kriging模型在未知點(diǎn)x處的預(yù)測(cè)值[10]：

(4)

預(yù)測(cè)值的方差和偏導(dǎo)值信息分別為

(5)

(6)

其中：R為樣本點(diǎn)的N×N相關(guān)函數(shù)矩陣, 組成元素為

Rij=Corr[z(x(i)),z(x(j))];

r為未知點(diǎn)x與樣本點(diǎn)間N×1的相關(guān)系數(shù)向量, 其中ri=Corr[z(x),z(x(i))];y是樣本點(diǎn)響應(yīng)值的N×1階向量;1為各元素均為1的N×1列向量; 估算值

模型預(yù)測(cè)方差

2 基于Kriging模型的多目標(biāo)粒子群優(yōu)化算法

2.1 局部領(lǐng)導(dǎo)粒子的更新策略

對(duì)于粒子群算法[11-13], 要完成粒子速度、 位置更新過程中的自學(xué)習(xí)部分, 就要合理地選擇一個(gè)局部領(lǐng)導(dǎo)者. 因此本文采用Pareto支配為主要準(zhǔn)則, 基于Kriging模型設(shè)計(jì)一種快速支配的局部領(lǐng)導(dǎo)者更新策略.

用樣本點(diǎn)及其目標(biāo)函數(shù)值建立Kriging模型, 可獲取目標(biāo)函數(shù)f(x)=(f1(x),f2(x),…,fm(x))及約束函數(shù)g(x)=(g1(x),g2(x),…,gr(x)). 參考式(5)得到目標(biāo)函數(shù)與約束函數(shù)預(yù)測(cè)值的誤差函數(shù)分別為

(7)

(8)

其中：k={1,2,…,m}表示第k個(gè)目標(biāo)函數(shù)響應(yīng)值;t={1,2,…,r}表示第t個(gè)約束函數(shù)響應(yīng)值.

對(duì)粒子各歷史位置的sfk(x)和sgt(x)值進(jìn)行支配關(guān)系比較, 非支配粒子歷史位置即為局部領(lǐng)導(dǎo)者. 在進(jìn)行支配關(guān)系比較前, 先根據(jù)所有粒子歷史軌跡的目標(biāo)函數(shù)、 約束函數(shù)預(yù)測(cè)值的誤差函數(shù)進(jìn)行預(yù)處理, 剔除誤差較大的粒子歷史位置, 以提高搜索的準(zhǔn)確性及效率. 若此時(shí)粒子的各歷史位置互不支配, 則將各位置的誤差函數(shù)進(jìn)行比較, 選擇誤差最小的歷史位置作為局部領(lǐng)導(dǎo)者.

大誤差粒子歷史位置剔除原則如下：

1) 根據(jù)式(7),(8)計(jì)算得到種群中所有粒子的sfk(x)和sgt(x)值, 先分別進(jìn)行歸一化處理, 再依次選取綜合罰值Vi,m>0.95與Vi,r>0.975的粒子歷史位置進(jìn)行剔除;

2) 對(duì)多維度sfk(x),sgt(x)進(jìn)行歸一化處理, 綜合考慮得到各維度的罰值分別為

(9)

通過

(10)

得到綜合罰值Vi,m與Vi,r[14]. 其中：i表示粒子的歷史位置;sfi,k表示目標(biāo)函數(shù)為k時(shí), 粒子當(dāng)前位置的誤差函數(shù);sfmin,k與sfmax,k分別表示目標(biāo)函數(shù)為k時(shí), 粒子所有離歷史位置誤差函數(shù)的最小值與最大值;sgi,t,sgmax,t,sgmin,t定義參考sfi,k,sfmin,k,sfmax,k;ak與at分別為粒子第i個(gè)歷史位置的目標(biāo)函數(shù)和約束函數(shù)綜合罰值中各維度罰值的權(quán)重.

2.2 非劣質(zhì)解儲(chǔ)備解集構(gòu)建機(jī)制

2.2.1 儲(chǔ)備解集的維護(hù) 隨著迭代過程的進(jìn)行, 越來越多的非劣質(zhì)解通過支配關(guān)系比較被找到, 并儲(chǔ)存在儲(chǔ)備解集中. 當(dāng)非劣質(zhì)解數(shù)量達(dá)到儲(chǔ)備解集的預(yù)設(shè)值上限時(shí), 即需對(duì)儲(chǔ)備解集進(jìn)行實(shí)時(shí)維護(hù), 即動(dòng)態(tài)擁擠距離維護(hù), 使種群能更準(zhǔn)確地向PF搜索, 保證PF的合理性.

非劣質(zhì)解的選取情形如下：

1) 在迭代過程中, 首先根據(jù)各粒子的誤差函數(shù)進(jìn)行大誤差粒子的剔除, 再進(jìn)行支配關(guān)系比較選擇非支配解(非劣質(zhì)解)進(jìn)入儲(chǔ)備解集;

2) 當(dāng)儲(chǔ)備解集中非劣質(zhì)解的數(shù)量未達(dá)到預(yù)設(shè)上限時(shí), 將所有非支配粒子全部?jī)?chǔ)存到解集中;

3) 當(dāng)儲(chǔ)備解集中非劣質(zhì)解的數(shù)量達(dá)到預(yù)設(shè)上限時(shí), 需根據(jù)動(dòng)態(tài)擁擠距離方法求解各非劣質(zhì)解的擁擠距離值(crowding distance, CD)并排序, 去除掉CD值最小的一個(gè)粒子, 重復(fù)上述操作, 直到非劣質(zhì)解的數(shù)量達(dá)到儲(chǔ)備解集上限要求.

2.2.2 基于Kriging模型的變異機(jī)制 由上述非劣質(zhì)解儲(chǔ)備解的維護(hù)機(jī)制可知, 隨著迭代次數(shù)的增加, 獲取的非支配粒子數(shù)量也大幅度增加. 雖然動(dòng)態(tài)擁擠距離可對(duì)儲(chǔ)備解集進(jìn)行維護(hù), 但也會(huì)丟失部分最優(yōu)解信息. 此外, 種群粒子的社會(huì)學(xué)習(xí)方式易導(dǎo)致算法出現(xiàn)早熟現(xiàn)象, 可能避開全局最優(yōu)而使所有粒子陷入局部最優(yōu). 為盡量避免種群算法的早熟現(xiàn)象, 保證求解精度, 本文利用Kriging模型響應(yīng)中的偏導(dǎo)值信息, 針對(duì)粒子群算法提出一種新的變異模式.

建立目標(biāo)函數(shù)和約束函數(shù)的復(fù)雜程度函數(shù)cf(x)和cg(x)為

(11)

再進(jìn)行歸一化處理得到罰值Vcf和Vcg為

(12)

在對(duì)所有粒子的位置和速度進(jìn)行更新時(shí), 獲取所有粒子約束函數(shù)的復(fù)雜程度函數(shù)值cg(x),對(duì)cg(x)=0的粒子進(jìn)行變異處理; 同時(shí)對(duì)cg(x)≠0的粒子, 利用本文定義的罰值Vcg, 取Vcg<0.1的粒子進(jìn)行變異. 粒子變異公式為

(13)

2.3 基于Kriging模型的加點(diǎn)策略

2.3.1 多目標(biāo)EIMe準(zhǔn)則及加點(diǎn) 在獲取初始樣本點(diǎn)后, 若進(jìn)一步提高模型準(zhǔn)確性, 則需獲取新的樣本點(diǎn)并進(jìn)行更多的試驗(yàn)設(shè)計(jì)(design of experiment, DOE)實(shí)驗(yàn), 以保證Pareto解集的前沿性. 但加點(diǎn)過程通常是盲目、 隨機(jī)的, 需進(jìn)行大量DOE實(shí)驗(yàn)才能達(dá)到期望的結(jié)果準(zhǔn)確性. 為解決該問題, 在保證模型結(jié)果準(zhǔn)確性的同時(shí)提高計(jì)算效率, 本文結(jié)合期望提高準(zhǔn)則(expected improvement, EI), 提出一種基于Kriging模型的多目標(biāo)粒子群優(yōu)化算法的加點(diǎn)策略, 使DOE、 Kriging代理模型、 MOPSO這三者間形成信息的協(xié)調(diào).

在維護(hù)后的儲(chǔ)備解中獲取每個(gè)非劣質(zhì)解目標(biāo)函數(shù)預(yù)測(cè)值的誤差函數(shù)sfk(x), 可得每個(gè)非劣質(zhì)解的置信度評(píng)估函數(shù)CEj為

(14)

得到上述信息后, 獲取每個(gè)非劣質(zhì)解基于Euler距離的期望提高準(zhǔn)則值(expected improvement matrix, EIMe)[9]為

(15)

適應(yīng)度函數(shù)的計(jì)算公式為

其中hkt為系數(shù),k={1,2,…,m},t={1,2,…,r}.

2.3.2 收斂準(zhǔn)則 當(dāng)實(shí)現(xiàn)加點(diǎn)策略后, 新的樣本點(diǎn)將被加入到初始點(diǎn)集中重新計(jì)算出新的Kriging模型, 并開始新一輪迭代過程. 只有當(dāng)Kriging模型達(dá)到一定的模型精度后才能完成迭代, 輸出Pareto前沿解集. 因此, 本文采用平均相對(duì)誤差檢驗(yàn)法驗(yàn)證模型的精度, 平均相對(duì)誤差的計(jì)算公式為

2.4 基于Kriging模型的粒子群優(yōu)化算法求解多目標(biāo)優(yōu)化問題流程

本文改進(jìn)的基于Kriging模型加點(diǎn)策略的多目標(biāo)粒子群優(yōu)化算法(KMOPSO)步驟如下：

1) 得到目標(biāo)問題, 定義其設(shè)計(jì)變量、 目標(biāo)和約束函數(shù), 在初始設(shè)計(jì)空間中對(duì)復(fù)雜系統(tǒng)進(jìn)行拉丁超立方試驗(yàn)設(shè)計(jì)生成初始樣本點(diǎn), 同時(shí)初始化種群;

2) 建立Kriging代理模型;

3) 大誤差粒子剔除后進(jìn)行支配關(guān)系比較, 非支配粒子進(jìn)入儲(chǔ)備解集, 完成儲(chǔ)備解集的更新, 當(dāng)儲(chǔ)備解集中非劣質(zhì)解達(dá)到上限時(shí), 通過動(dòng)態(tài)擁擠距離的計(jì)算進(jìn)行儲(chǔ)備解集的維護(hù);

4) 找到粒子局部領(lǐng)導(dǎo)者與全局領(lǐng)導(dǎo)者, 然后進(jìn)行粒子速度、 位置的更新, 同時(shí)利用構(gòu)建的復(fù)雜程度函數(shù)進(jìn)行變異處理完成粒子更新;

5) 對(duì)維護(hù)后儲(chǔ)備解集中的非劣質(zhì)解進(jìn)行置信度評(píng)估函數(shù)的獲取, 并得到近似Pareto前沿解集, 根據(jù)式(15)的期望準(zhǔn)則從儲(chǔ)備解集中尋找新的樣本點(diǎn)加入原樣本點(diǎn)集中;

6) 重復(fù)步驟1)～5), 直到滿足收斂準(zhǔn)則時(shí)完成迭代, 找到最優(yōu)解.

3 算法性能測(cè)試

3.1 性能指標(biāo)IGD

針對(duì)提出的KMOPSO算法, 本文用綜合指標(biāo)(inverted generational distance, IGD)評(píng)價(jià)算法, 并采用多目標(biāo)測(cè)試函數(shù)ZDT1,ZDT2,ZDT3,ZDT6[15]進(jìn)行驗(yàn)證對(duì)比. IGD值表示算法計(jì)算出的近似與真實(shí)Pareto前沿解集之間的距離, IGD值大則表明近似Pareto前沿解集有較差的收斂性和多樣性. IGD值的計(jì)算公式為

(18)

其中：P為均勻分布在PF上的點(diǎn)集, |P|為P的種群數(shù)量;Q為算法獲取的最優(yōu)Ps;d(v,Q)為P中個(gè)體v到種群Q的最小歐氏距離.

3.2 仿真驗(yàn)證及分析

下面對(duì)本文算法與多目標(biāo)人工蜂群算法(multi-objective artificial bee colony algorithm, MOABC)和多目標(biāo)遺傳算法(non-dominated sorting genetic algorithm Ⅱ, NSGAⅡ)進(jìn)行對(duì)比實(shí)驗(yàn). 算法參數(shù)設(shè)置如下： 種群數(shù)量為200, 儲(chǔ)備解集大小為300, 學(xué)習(xí)因子c1=c2=1.5, 慣性權(quán)重ω=0.73, 變量維度設(shè)為10, 并分別獨(dú)立進(jìn)行10次計(jì)算, 所得IGD性能指標(biāo)列于表1.

表1 IGD性能指標(biāo)

由表1可見, 本文提出的KMOPSO算法在上述4個(gè)測(cè)試函數(shù)驗(yàn)證下, 其IGD指標(biāo)的均值與標(biāo)準(zhǔn)差均小于另外兩種對(duì)比算法, 即KMOPSO算法在高維問題下所得Pareto前沿解集最接近真實(shí)解集. 這是由于KMOPSO提出的加點(diǎn)策略更接近目標(biāo), 且粒子的變異模式也能更好地進(jìn)行全局搜索.

圖1 3種算法在測(cè)試函數(shù)ZDT1下返回的非支配前沿解集Fig.1 Non-dominated frontier solution sets returned by three algorithms on test function ZDT1

仿真結(jié)果如圖1～圖4所示, 其中包括ZDT1,ZDT2,ZDT3和ZDT6的非支配前沿解. 圖1為測(cè)試函數(shù)ZDT1下3種算法返回的非支配解集. 由圖1(A)可見, 使用MOABC算法得到的結(jié)果收斂性較差; 由圖1(C)可見, KMOPSO算法的多樣性比NSGAⅡ算法差, 但Pareto鋒面的形狀更好. 由圖2～圖4可見, 使用KMOPSO算法獲得的結(jié)果具有理想的收斂性與多樣性.

圖2 3種算法在測(cè)試函數(shù)ZDT2下返回的非支配前沿解集Fig.2 Non-dominated frontier solution sets returned by three algorithms on test function ZDT2

圖3 3種算法在測(cè)試函數(shù)ZDT3下返回的非支配前沿解集Fig.3 Non-dominated frontier solution sets returned by three algorithms on test function ZDT3

圖4 3種算法在測(cè)試函數(shù)ZDT6下返回的非支配前沿解集Fig.4 Non-dominated frontier solution sets returned by three algorithms on test function ZDT6

綜上可見, 本文提出的基于Kriging模型加點(diǎn)策略的粒子群優(yōu)化算法, 首先利用Kriging模型的響應(yīng)信息, 改進(jìn)了粒子的變異模式以更好地解決粒子“早熟”問題; 然后在迭代過程的儲(chǔ)備解集中選取非劣質(zhì)解進(jìn)入樣本點(diǎn)集, 極大提高了計(jì)算結(jié)果的準(zhǔn)確性; 最后提出大誤差粒子剔除原則, 使支配關(guān)系比較的結(jié)果更精確. 仿真實(shí)驗(yàn)結(jié)果表明, KMOPSO算法可用于解決高維度、 非線性多目標(biāo)問題, 且可利用較少的樣本點(diǎn)得到更精確的代理模型, 在提高收斂速度的同時(shí)能較好地逼近Pareto前端.