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

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

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

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

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

當優(yōu)化問題有m個目標函數(shù)、r個約束函數(shù)時, 其數(shù)學模型可表示為

(1)

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

定義1(Pareto支配)[6]對于式(1), 假設有設計變量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]對于式(1), 若不存在滿足x(3)x(4)的設計變量x(3), 則稱x(4)為非支配解或Pareto最優(yōu)解.

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

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

1.2 Kriging模型

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

(2)

(3)

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

(4)

預測值的方差和偏導值信息分別為

(5)

(6)

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

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

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

模型預測方差

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

2.1 局部領導粒子的更新策略

對于粒子群算法[11-13], 要完成粒子速度、 位置更新過程中的自學習部分, 就要合理地選擇一個局部領導者. 因此本文采用Pareto支配為主要準則, 基于Kriging模型設計一種快速支配的局部領導者更新策略.

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

(7)

(8)

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

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

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

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

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

(9)

通過

(10)

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

2.2 非劣質解儲備解集構建機制

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

非劣質解的選取情形如下：

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

2) 當儲備解集中非劣質解的數(shù)量未達到預設上限時, 將所有非支配粒子全部儲存到解集中;

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

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

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

(11)

再進行歸一化處理得到罰值Vcf和Vcg為

(12)

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

(13)

2.3 基于Kriging模型的加點策略

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

在維護后的儲備解中獲取每個非劣質解目標函數(shù)預測值的誤差函數(shù)sfk(x), 可得每個非劣質解的置信度評估函數(shù)CEj為

(14)

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

(15)

適應度函數(shù)的計算公式為

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

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

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

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

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

2) 建立Kriging代理模型;

3) 大誤差粒子剔除后進行支配關系比較, 非支配粒子進入儲備解集, 完成儲備解集的更新, 當儲備解集中非劣質解達到上限時, 通過動態(tài)擁擠距離的計算進行儲備解集的維護;

4) 找到粒子局部領導者與全局領導者, 然后進行粒子速度、 位置的更新, 同時利用構建的復雜程度函數(shù)進行變異處理完成粒子更新;

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

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

3 算法性能測試

3.1 性能指標IGD

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

(18)

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

3.2 仿真驗證及分析

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

表1 IGD性能指標

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

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

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

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

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

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

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