王 博,劉連生,韓紹程,祝世興
(1.中國民航大學基礎實驗中心,天津 300300;2.中國民航大學航空工程學院,天津 300300)
工程實踐和科學研究中很多決策及設計問題同時涉及兩個或以上的目標,這些問題稱為多目標優化(Multi-Objective Optimization,MOO)[1]。由于MOO 問題中各目標函數之間一般存在內在沖突,很難得到使所有目標同時達到最佳狀態的解決方案,求解過程中通常需要經過權衡得到一組相互折中的解集[2]。
近幾十年來,國內外學者針對MOO 問題開展了大量研究工作,涌現了眾多優化算法,主要可以分為兩類:基于偏好信息的方法和Pareto 最優解法[3]。第一類優化方法根據每個目標函數的重要性設定權重系數并計算它們的加權和,從而將MOO 問題轉化為單目標優化問題。此類方法較為簡單,但權重系數的選擇直接影響最終優化結果,實際應用中權重系數的設定帶有一定主觀性、隨機性,很難得到客觀準確的優化結果。第二類優化方法同時對多個子目標優化,通過搜索得到一組不存在優劣關系的非支配最優解(即Pareto 最優解)。該方法能夠獲得一組較高質量的最優解集,供決策者從中選擇最合適的解,相較第一類方法應用更廣泛[4]。
蝗蟲優化算法(Grasshopper Optimization Algorithm,GOA)是Saremi 等[5]受蝗蟲捕食過程中種群行為的啟發,在2016 年提出的一種新型啟發式智能優化算法。GOA 簡單,計算復雜度較低,具有較高的收斂速度。近年來,為提高GOA的性能眾多學者做了大量研究工作。Mafarja 等[6]利用動態種群進化(Evolutionary Population Dynamics,EPD)策略改進了GOA 的種群選擇機制,有效避免了算法過早收斂,提高算法的準確性。Luo 等[7]提出將高斯變異引入GOA 以增強局部搜索能力,利用Levy 飛行策略和反向學習機制(Levy flight and Opposition-Based Learning,LOBL)提高了算法的全局搜索能力和搜索效率。Taher 等[8]改進了GOA 的變異機制和種群更新策略,避免了局部最優。Mirjalili等[9]于2017年首次將GOA應用于多目標問題,提出了多目標蝗蟲優化算法(Multi-Objective Grasshopper Optimization Algorithm,MOGOA)。黃超等[10]利用曲線自適應策略調整控制參數,提高算法收斂速度,提出一種求解移動機器人路徑規劃問題的MOGOA。目前,多目標蝗蟲優化算法由于采用了標準GOA 的基本框架,繼承了蝗蟲優化算法在復雜度和收斂速度上的優勢,具有較好的應用前景;但由于GOA 的天然缺陷,限制了其優化性能。首先,由于種群采用優勢個體引導的更新機制,易造成群體進化聚集,算法容易陷入局部最優;其次,調整算法搜索能力的重要參數c只隨迭代次數被動變化,導致算法全局搜索和局部開發的不平衡,降低了收斂效率和解集質量。
針對上述問題,本文提出一種混合多目標蝗蟲優化算法(Hybrid Multi-Objective Grasshopper Optimization Algorithm,HMOGOA)。HMOGOA 利用Halton 序列進行種群初始化,增強了初始種群的多樣性以及分布均勻性;由差分變異算子引導種群更新,提高了算法全局搜索能力避免陷入局部最優,改善了解集分布特性;加入自適應權重因子,根據種群優化情況動態調整算法的全局搜索和局部開發能力,提高了算法優化效率和解集質量。采用ZDT、CEC09 系列函數對算法進行仿真實驗,通過與多目標蝗蟲優化、多目標粒子群(Multi-Objective Particle Swarm Optimization,MOPSO)、基于分解的多目標進化(Multi-Objective Evolutionary Algorithm based on Decomposition,MOEA/D)、非支配排序遺傳(Non-dominated Sorting Genetic Algorithm Ⅱ,NSGA Ⅱ)等算法對比,表明HMOGOA具有更好的優化性能和穩定性。
多目標蝗蟲優化算法是一種求解多目標優化問題的GOA。它模擬了蝗蟲的捕食行為,算法中每個搜索個體(即蝗蟲個體)的位置代表了問題的一種解決方案。算法通過判斷不同解之間的支配關系尋找Pareto解,并借助外部檔案集存儲優化過程中的臨時解。其優化過程分為探索和開發兩部分。在探索階段,搜索個體在整個解空間內進行快速隨機移動尋找所有可能的解;而在開發階段,搜索個體根據已有信息對部分區域進行精細搜索尋找最優解[11]。這恰好與蝗蟲在成蟲和幼蟲時期的運動特點相吻合。蝗蟲種群運動的數學模型[12]為:

其中:Xi是第i只蝗蟲在種群中的位置,Si、Gi、Ai分別代表第i只蝗蟲受到的種群間社交作用、重力以及風力的影響。種群間的社交作用是影響蝗蟲運動的最主要因素,可通過下述公式描述:

其中:dij為第i只與第j只蝗蟲之間的距離,為第i只與第j只蝗蟲之間的單位向量,函數s(r)為蝗蟲間的社會作用力。由于重力對蝗蟲種群的影響較小可以忽略不計,同時假設風向始終指向最佳個體所處的位置,蝗蟲位置的更新模型可表示為:

其中:ubd、lbd分別表示蝗蟲在d維度上的位置上界和下界表示當前種群中優勢引導個體的位置;參數c為線性縮減系數,它隨著迭代次數增加減小算法全局搜索能力的同時提高局部尋優能力,計算公式如下所示:

標準MOGOA流程如圖1所示。

圖1 標準MOGOA流程Fig.1 Flow chart of standard MOGOA
為提高MOGOA 的性能,需解決兩方面問題:1)種群更新過程中,缺少個體間信息的交換,進化容易集中在優勢個體的鄰域范圍內進行,導致算法全局搜索能力下降,算法過早收斂陷入局部最優。2)衰減系數c是GOA 求解MOO 問題過程中的一個重要參數,它與迭代次數呈線性反比關系,主要限制算法的搜索區域,平衡全局搜索能力和局部尋優能力。優化前期,c數值較大,促進個體在整個解空間內搜索所有可能的優勢解;優化后期,其數值減小,有利于種群在當前解的局部區域尋優,使最優解能夠收斂到理想Pareto 解集。通常MOO 的優化過程并不隨迭代次數線性變化,應當根據種群進化過程合理控制算法的搜索能力。
因此,本文從提高初始種群質量降低種群進化聚集的概率,加強進化過程中個體信息交換、提高種群多樣性,根據種群進化狀態平衡算法全局搜索和局部開發能力3個方面,提高MOGOA的優化性能,提出一種混合多目標蝗蟲優化算法。
算法初始階段蝗蟲種群的位置隨機產生,極易出現種群分布不均的情況,造成搜索個體的初始位置集中于局部極值點附近,導致算法陷入局部最優。Halton 序列是一種低差異化序列,能夠生成在搜索空間內均勻分布的點集[13]。由隨機法和Halton 序列生成的一組二維個體分布情況如圖2 所示。圖2(a)、(b)分別為使用兩種方法生成個體的散點分布圖,圖2(c)、(d)為個體數量的分布直方圖。通過對比可見,相較于隨機方法Halton序列產生的個體在二維空間上分布更均勻,沒有聚集情況,同時處于不同數值區間的個體數量更平均,表明Halton序列能夠生成分布均勻性更好的個體。為此,本文采用Halton 序列初始化種群,提高初始種群質量有效避免局部最優。

圖2 隨機法與Halton序列生成個體的分布情況Fig.2 Distribution of individuals generated by random method and Halton sequence

MOO 需要在優化過程中不斷尋找所有非支配解,組成非支配解集。由于多目標問題的解空間由多個分段連續區域組成,為使解集能夠逼近真實Pareto 最優前沿且分布更均勻,優化過程中搜索個體在決策空間中的搜索方向應指向真實Pareto 前沿,且可以在平行于真實Pareto 前沿的流行區域移動[15-16]。圖3 反映了MOGOA 優化過程中種群進化趨勢。圖中以ZDT1為目標函數,由10個初始父代解利用MOGOA 產生1 000 個新的子代解。可以明顯看出搜索個體由優勢個體引導沿不同軌跡向著目標位置移動,最終聚集到真實Pareto 前沿的部分區域,表明MOGOA 的種群更新策略易使算法陷入局部最優。針對該問題,在單目標蝗蟲優化算法中,文獻[6]中采取了動態種群進化策略,文獻[7]中利用了Levy 飛行和反向學習策略改進個體的更新機制[2]。分別將EPD 和LOBL策略應用于多目標蝗蟲優化算法,以ZDT1 為測試函數,由圖3 中的10 個父代解出發,產生1 000 個子代解,得到的種群進化情況如圖4、5 所示。圖4 中EPD 策略一定程度擴大了搜索個體在進化過程中的尋優范圍,但依然無法避免解聚集;而圖5 中LOBL 策略雖然增強了算法的全局搜索能力,但種群搜索范圍過大影響了算法的收斂性能。由此說明,EPD 和LOBL策略對MOGOA性能的改善有限。

圖3 MOGOA進化示意圖Fig.3 Schematic diagram of MOGOA evolution

圖4 基于EPD的MOGOA進化示意圖Fig.4 Schematic diagram of MOGOA evolution based on EPD

圖5 基于LOBL的MOGOA進化示意圖Fig.5 Schematic diagram of MOGOA evolution based on LOBL
由于差分算子有利于種群間信息交換,促進全局探索,改善種群分布和收斂精度,因此本文提出一種基于差分變異算子的種群更新策略。差分變異算子及個體更新公式如下:


式(9)中,α0、α1為縮放系數(α0,α1∈(0,2]),r為(0,1)內的隨機數,Xsi(t)表示以較大概率在解集稀疏區域選取的個體,Xc1i(t)、Xc2i(t)表示以較大概率在解集密集區域選取的個體,Xsi(t)表示由當前解集隨機選取的個體,上述4個種群個體互不相同。Xsi(t)、Xc1i(t)、Xc2i(t)組成的差分算子實現了種群個體間優勢信息的交換,能夠提高算法的全局尋優能力,有助于算法跳出局部最優,改善收斂精度。而隨機個體Xri(t)可以增加變異的隨機性,避免過早收斂。為驗證差分變異算子對MOGOA的改善作用,同樣選取圖3中的父代解作為初始解,以ZDT1為目標函數,利用基于差分變異算子的MOGOA生成1000個子代解,進化結果如圖6所示。由圖中可見,差分變異算子能夠引導種群在父代解構成的帶狀區域變異,可以改善解空間的整體分布性能。

圖6 基于差分變異算子的MOGOA進化示意圖Fig.6 Schematic diagram of MOGOA evolution based on differential evolution operators

多目標問題的優化是一個復雜的過程具有隨機性,MOGOA 求解過程中c是一個開環參數,單純通過它無法有效平衡全局搜索和局部尋優能力[18]。根據種群進化狀態調節參數c的衰減速度,實現閉環控制,有助于提高種群的搜索效率和收斂精度。通常種群進化過程中個體的分布特性是影響算法性能的一個重要因素[19]。本文以群體分布均勻性和分散廣度作為進化過程反饋信息,提出種群進化狀態動態評價指標(Evolution State Metric,ESM)。ESM定義如下:

其中:Ns為當前解的數量,M為目標函數的數量,表示當前Pareto 解集中相鄰解的歐氏距離表示的平均值表示理想Pareto 解集中使第m個目標函數取得最大值的解(即理想解空間第m維上的極值解),表示當前解集與理想解空間第m維上極值解的最短距離。式(13)中第一項表示當前解集相鄰解之間距離的標準差,體現了解集分布的均勻度,數值越小解集分布越均勻;第二項表示當前解集與理想解空間各維度上極值點的最短距離之和,體現了解集分布的廣度,其值越小說明解集分布范圍越廣。圖7所示為MOGOA針對目標函數ZDT4 運行三次得到的ESM 數值隨迭代次數變化曲線。可以看出ESM 是波動變化的,算法每次運行過程中其變化幅度與趨勢不同,說明種群進化狀態具有非線性和不確定性,根據進化過程調節算法搜索能力有助于提高算法性能。因此,本文建立基于ESM 的自適應權重因子wa,在種群更新公式中加入自適應權重因子,通過權重因子調整參數c的衰減速度平衡算法的全局搜索和局部尋優能力。
自適權重因子wa和種群個體更新公式如下:

其中,wa數值根據種群進化狀態變化。種群進化狀態不佳時ESM 數值較大,說明解分布較集中、均勻性不好,此時wa數值增加大可減緩c的衰減速度,促進搜索個體在更大區域尋找新解;反之ESM 數值較小,說明搜索個體在目標空間內分布均勻且分布范圍廣,wa數值減小,可以促進算法局部尋優,提高收斂精度和收斂速度。實現了根據算法尋優過程動態調整全局搜索和局部開發能力,提高了算法的優化效率,有效改善了Pareto最優解集的質量。

圖7 ESM變化曲線Fig.7 Variation curve of ESM
混合多目標蝗蟲優化算法的偽代碼流程如下:


本文為驗證HMOGOA的性能,以MOGOA、MOPSO、MOEA/D及NSGA Ⅱ作為對比算法,選取ZDT、CEC2009系列測試函數集中的兩目標函數(ZDT1、ZDT2、ZDT3、ZDT4)和三目標函數(UF8、UF9、UF10),使用軟件Matlab2016a進行仿真實驗測試。
實驗過程中,所有算法的相關參數均設置一致以保證公平性:種群規模N=100,外部檔案集規模NA=100,迭代次數Tmax=100。HMOGOA、MOGOA 參數c的取值范圍均設置為:cmax=1,cmin=0.000 5。其他算法參數設置如下:MOPSO算法慣性權重因子w=0.5,個體學習因子c1=2,社會學習因子c2=2,變異率pm=0.1;MOEA/D算法鄰域大小T=15,交叉概率pc=0.5;NSGA Ⅱ算法交叉概率pc=0.7,變異概率pm=0.02。為避免隨機誤差對測試結果的影響,針對每個測試函數所有算法均運行30次。
為分析比較不同優化算法的性能,選取兩種評價指標來量化分析算法的性能。具體性能指標[20]如下:
1)反世代距離(Inverted Generational Distance,IGD)。IGD 反映了算法求得的Pareto 解集P與真實Pareto 前沿之間的距離,它是評價算法收斂性的指標。其值越小說明所求解越逼近真實Pareto 前沿,表明算法的收斂性越好,解的精度越高。IGD定義如下:

其中,NPt為真實Pareto前沿中解的數量為真實Pareto前沿中的第i個解與解集P間的最短歐氏距離。
2)分布度量指標(Δ)。Δ體現了所求Pareto 最優前沿相對真實前沿的分布廣度及均勻度,它是評價Pareto 解分布特性的指標。Δ數值越小說明解的多樣性越高、分布范圍越廣、分布均勻度越好,非支配解能夠更好地覆蓋整個真實Pareto解集。Δ的定義如下:

其中:Np為所求解的個數,di為Pareto前沿中相鄰點間的距離,為di的平均值,df、dl為Pareto 最優前沿與真實前沿邊界點之間的最短距離。
五種算法針對多目標問題測試所得各項性能指標的均值(Mean)和方差(Var)如表1、2所示,同一測試問題中的最優數值加粗顯示。

表1 IGD實驗結果Tab.1 Results of IGD

表2 Δ實驗結果Tab.2 Results of Δ
由表1 可見,相較于MOPSO、MOEA/D、NSGA Ⅱ算法,HMOGOA 的IGD 數值在大部分測試函數上明顯優于此三種算法,只是在函數UF9 上稍差于MOPSO,但數值相近且明顯優于其他兩種算法,表明HMOGOA 的收斂性優于此三種算法。而對于MOGOA,HMOGOA 的收斂性能和收斂穩定性在全部測試問題上都體現出了明顯優勢。由此說明HMOGOA采取的基于差分變異算子的更新策略有助于引導種群向理想解空間中目標個體位置移動。
通過表2 對比分析五種算法的分布性能指標可以看出,由于HMOGOA 在進化初期提高了初始種群分布性,避免了算法前期陷入局部最優;差分變異算子能夠提高全局搜索能力,改善解集分布特性;同時采用自適應權重因子平衡全局和局部尋優能力,提高了最優解集的質量。因此,HMOGOA 在所有測試問題上的分布性能均明顯優于其他算法。
圖8 所 示 為HMOGOA、MOGOA、MOPSO、MOEA/D 的Pareto 最優前沿分布圖。通過圖8 可以看出,對于函數ZDT3,HMOGOA 求得的解相較其他三種算法能均勻分布在函數的4個波谷上,收斂精度更好。對于函數ZDT1、ZDT2,HMOGOA與其他算法相比能更均勻分布在整個真實Pareto 前沿面上。而對于三目標問題,HMOGOA 求得解的數量更多,能夠較好逼近真實Pareto前沿,且分布更均勻。

圖8 各算法對部分函數的Pareto前沿對比Fig.8 Comparison of Pareto optimal front for some benchmark functions obtained by various algorithms
為進一步比較HMOGOA 與其他算法性能之間的顯著優勢,本文利用Wilcoxon 秩和檢驗按照顯著性水平α=0.05,將HMOGOA 分別與其他四種算法進行對比分析。分析結果如表3 所示,“+”“-”“≈”分別表示HMOGOA 的性能指標明顯優于、劣于或近似于對比算法。由表3 可以看出,HMOGOA 的IGD、Δ指標在大部分測試函數上,相較于MOGOA、MOPSO、MOEA/D、NSGA Ⅱ具有明顯優勢,只是在個別測試函數上近似于MOPSO、MOEA/D。通過上述實驗分析可知,本文提出的HMOGOA 性能明顯優于MOGOA,算法具有較高的收斂精度和較好的分布性能,且穩定性較高。

表3 Wilcoxon秩和檢驗結果Tab.3 Results of Wilcoxon rank-sum test
本文為提高MOGOA 的優化性能,提出一種混合多目標蝗蟲優化算法。該算法利用Halton 序列初始化種群,使初始種群具有較好的多樣性和分布特性,減小了算法過早收斂的概率;設計基于差分變異算子的種群更新機制,促進種群向優勢個體移動的同時進行更大范圍尋優,實現種群間優勢個體的信息交換,提高了算法的全局尋優能力,有效避免算法陷入局部最優,提升了解集的分布特性,改善了收斂精度;引入自適應權重因子,根據種群進化狀態反饋信息動態調整算法的全局搜索和局部開發能力,實現算法尋優能力的閉環精細控制,提高了算法的優化效率和最優解集質量。通過與四種多目標優化算法的對比實驗表明,HMOGOA 的性能相較于MOGOA有明顯提升,算法具有較好的收斂精度、分布廣度、分布均勻性和穩定性。下一步研究的主要內容是算法對于求解高維多目標優化問題的性能。