基于流形結構重建的多目標氣動優(yōu)化算法

宋超，李偉斌，周鑄，劉紅陽，藍慶生

中國空氣動力研究與發(fā)展中心 計算空氣動力研究所，綿陽 621000

飛行器氣動外形的優(yōu)化設計必須考慮多個目標，例如，民用飛機需要綜合考慮巡航升阻比、阻力發(fā)散、力矩特性以及抖振邊界等問題[1]，直升機旋翼需要在懸停、前飛和機動等多種運動狀態(tài)下有優(yōu)良的高速超臨界性能、低速大迎角氣動性能，以及接近于零的力矩系數(shù)[2]。然而，多個目標往往是互相矛盾的，一個目標的提升可能是以另一個目標的妥協(xié)為代價[3]，因此需要一種全面衡量所有設計目標的系統(tǒng)設計方法。

Fonseca和Fleming[4]利用Pareto最優(yōu)解集概念，提出了多目標遺傳算法(Multi Objective Genetic Algorithm, MOGA)。隨后，Srinivas和Deb[5]提出了非劣分層遺傳算法(Non-dominated Sorting Genetic Algorithm, NSGA)。該算法基于非支配排序思想設計，能夠有效解決多目標設計問題，得到均勻分布的解。Deb等對NSGA進行改進，并提出了基于精英保留策略的非劣排序進化算法(NSGA2)[6]，與NSGA相比，NSGA2的計算復雜度有所降低，并且種群的分布更加均勻。此類傳統(tǒng)的多目標進化算法已經(jīng)廣泛應用于飛行器的多目標優(yōu)化設計中[7-11]，并取得了一定的成果。然而，演化多目標優(yōu)化算法具有一些固有缺點，比如當種群接近收斂時，盲目交叉、變異操作會使本已近似收斂的種群偏離實際的Pareto解集(Pareto Set, PS)，給算法性能造成不良影響[12]。

分布估計算法(Estimation of Distribution Algorithm, EDA)[13]通過建立解集分布的概率模型，然后在該模型上隨機產(chǎn)生新的子代個體，從而實現(xiàn)種群進化，避免了傳統(tǒng)的交叉變異操作。EDA憑借其高效優(yōu)化思想，已經(jīng)成了多目標優(yōu)化領域的研究熱點之一[14-15]。Zhang等[16]利用設計空間的Pareto解集是一個分段連續(xù)的k維流形(k=m-1，m為設計目標個數(shù))這一規(guī)律[17]，發(fā)展了基于規(guī)則模型的多目標EDA。如何建立解集分布的概率模型，是EDA的核心問題。一些建模方法，如主成分分析(Principal Component Analysis, PCA)[18]等需要對種群進行聚類，過程復雜且準確度不高。詹煒[19]利用局部線性嵌入(Locally Linear Embedding, LLE)方法發(fā)掘Pareto前沿的流形結構，進行了衛(wèi)星星座的多目標優(yōu)化設計，取得了較好的效果。Yang[20]等采用LLE和免疫克隆算法來發(fā)掘Pareto前沿結構。目前，基于流形結構的EDA在氣動多目標優(yōu)化問題中的應用還在初步研究階段。本文提出一種利用流形結構重建的分布估計算法，將其用于多目標氣動優(yōu)化問題，可顯著提高優(yōu)化效率，實現(xiàn)多目標氣動優(yōu)化問題的高效求解。

1 算法介紹

分布估計算法首先建立一個描述解集分布的概率模型，在此基礎上對概率模型進行采樣，得到新的解集分布，如此反復進行，實現(xiàn)種群的進化，直到終止條件[21]。本文采用流形結構重建的方法來建立概率模型，并對目標空間的Pareto前沿進行擴展以實現(xiàn)概率模型的采樣，實現(xiàn)了基于流形結構重建的多目標優(yōu)化算法。

1.1 流形結構重建算法

優(yōu)化算法迭代過程中的種群逐漸形成清晰的Pareto前沿，高維設計空間中的Pareto解集存在一定的相關性和規(guī)律性。根據(jù)卡羅需-庫恩-塔克條件(Karush-Kuhn-Tucker condition， KKT)[16]，設計空間的Pareto解集是一個分段連續(xù)的低維流形。若獲得低維流形表示，可挖掘出其中蘊含的信息，進而指導優(yōu)化算法的演化。另一方面，由于低維流形的簡單性，一定程度上解決了高維設計變量帶來的諸多問題，如維數(shù)災難。

已有的利用流形結構的多目標優(yōu)化算法，借助了流形學習算法來發(fā)掘其流形結構。大多數(shù)流形學習算法，如前面提到的PCA、LLE等方法，都只是將高維空間中的樣本數(shù)據(jù)點投影到低維流形上，投影關系并沒有具體的解析函數(shù)表達式。具體地，針對氣動優(yōu)化問題，雖然可以采用流形學習算法將設計空間中的樣本點投影到低維的流形空間，但在低維流形空間中的點不具有物理意義，直接利用這些信息是困難的。即低維流形空間中的任意坐標不能顯式地對應于設計空間中的點，因此也就無從求得對應的氣動力。

x=g(y)=xi+J(yi)(y-yi)+R(y,yi)

(1)

式中：J(yi)為函數(shù)g(y)在點yi的雅克比矩陣；R(y,yi)為余項。略去余項并假設xi有z個鄰域點，則有

(2)

(3)

至此得到了點yi的雅克比矩陣J(yi)。

對于流形空間Sk中的任意一點y，解出其在高維空間Dn中坐標x的過程稱為流形結構重建。對應于氣動優(yōu)化問題，即給定目標空間中的一點，求其設計空間中的坐標。流形結構重建的具體步驟如下：

步驟2 找出ys在高維空間Dn中的坐標xs,xs的鄰域點xs,1,xs,2,…，xs,z以及這些鄰域點在Sk中的坐標ys,1,ys,2, …，ys,z。

步驟4 點y在Dn中的坐標x≈xs+Qs(y-ys)。

值得注意的是，以上流形結構重建算法建立的高維流形空間與對應的低維表示空間的顯式映射關系，是基于流形結構的連續(xù)性、局部等距性等一系列假設的，更詳細的約束條件及過程推導見文獻[22]。

1.2 目標空間Pareto前沿擴展

利用1.1節(jié)的流形結構重建算法，可以得到低維目標空間到高維設計空間的映射關系。通過在目標空間Pareto前沿的擴展就可以得到對應的設計空間的解集分布，用于指導算法進化。生成新個體的算法如圖1所示。設ζ為流形結構，將種群個體看作是一個隨機變量ξ的獨立觀測量，且ε是一個符合正態(tài)分布的隨機噪聲向量。種群個體在ζ附近隨機分布為

ξ=ζ+ε

(4)

這一方法也被Zhang等[16]、詹煒[19]采用。本文直接在目標空間對Pareto前沿進行“擴展”，一方面極大簡化了算法，另一方面目標函數(shù)空間中的流形結構有清晰的物理意義，能夠指導算法演化。

一般地，假設希望達到比當前最優(yōu)Pareto前沿上目標函數(shù)更小的值，隨機向量ε取值為

(5)

式中：上標i表示向量的第i個分量；N(yi)為均值為1、方差為1的正態(tài)分布函數(shù)；κ為推進因子。

圖1 流形結構擴展示意圖Fig.1 Illustration of manifold structure extension

1.3 基于流形結構重建的多目標優(yōu)化流程

本文采用流形結構重建算法估計解集的概率分布，利用Pareto前沿的流形分布規(guī)律，用于指導多目標優(yōu)化算法的演化，以下稱該算法為MR-EDA。記算法種群為

Pop(t)=[x1,x2,…，xN]

(6)

式中：t為當前迭代步數(shù)；N為種群大小。個體適應度值F(x)=[f1(x),f2(x),…，fm(x)]Τ，m為目標個數(shù)。采用MR-EDA法的優(yōu)化流程如圖2所示，具體如下：

步驟1 產(chǎn)生初始種群Pop(0)，計算種群中的每個個體的適應度值F(x)。

步驟2 對當前種群進行非支配排序，選取前N1個個體，組成新的種群Q(t)。針對Q(t)，利用1.2節(jié)所述的流形結構重建算法，生成N2個個體，形成種群P(t)。其中N=N1+N2。

步驟3 若滿足停止條件，迭代結束。否則，Pop(t)=Q(t)+P(t)(t=t+1)，返回步驟1繼續(xù)進行迭代。

圖2 MR-EDA優(yōu)化流程圖Fig.2 Flowchart of MR-EDA optimazation

2 解析函數(shù)測試算例

ZDT測試集是目前廣泛采用的兩目標測試函數(shù)。本文選取了ZDT1、ZDT2、ZDT3作為測試函數(shù)，其中ZDT1的最優(yōu)前沿為凸函數(shù)，ZDT2的最優(yōu)前沿為凹函數(shù)，ZDT3的最優(yōu)前沿為分段凹函數(shù)。測試函數(shù)的定義如表1所示，其中自變量維數(shù)m均取為30。

表1 測試函數(shù)定義Table 1 Definitions of test functions

注：上標1表示矢量的第1個維度的值。

為評價算法性能，引入IGD(Inverted Generational Distance)指標[23]。該指標首先需要均勻分布在真實Pareto最優(yōu)前沿上的點集S*，設S為優(yōu)化算法得到的Pareto非支配解集，則IGD定義為

(7)

式中：d(v,S)為點v到解集P所有解中的最小距離；|S*|為樣本集S*中的點數(shù)。IGD指標能同時度量收斂性和多樣性。越小的IGD值表示優(yōu)化算法得到的非支配解越逼近真實最優(yōu)解，同時非支配解的均勻性也越好。

本文分別采用MR-EDA與NSGA2進行解析函數(shù)的優(yōu)化測試。設計變量數(shù)目為30，初始種群數(shù)目為100，最大迭代步數(shù)設為200。NSGA2的交叉概率為0.8，變異概率為0.3。MR-EDA中流形重構種群數(shù)目N1=N/2。推進因子κ的取值過小可能影響收斂效率，取值過大則不滿足流形重構的條件，經(jīng)測試κ=0.05取得較好的效果。因此，本文算例中推進因子κ均取為0.05。圖3(a)給出了測試函數(shù)ZDT1的優(yōu)化結果，分別給出了ZDT1函數(shù)的真實最優(yōu)Pareto前沿，NSGA2及本文算法MR-EDA得到的最優(yōu)Pareto前沿。如圖3(a)所示，本文所提算法優(yōu)化得到的Pareto最優(yōu)前沿與真實前沿幾乎重合，并且在Pareto解的分布更加均勻。經(jīng)過相同迭代步數(shù)的NSGA2得到的最優(yōu)Pareto前沿與真實Pareto前沿還有一定的距離。圖3(b)給出了ZDT2測試函數(shù)的優(yōu)化結果。類似于ZDT1優(yōu)化結果，本文算法得到的ZDT2函數(shù)最優(yōu)Pareto前沿相比NSGA2收斂更好，前者基本與真實Pareto前沿重合。圖3(c)為ZDT3函數(shù)的測試結果，經(jīng)過200步的迭代，本文算法及NSGA2均未收斂到真實Pareto前沿，但本文算法明顯地得到了更好的Pareto解集。

圖3 3種測試優(yōu)化算例的Pareto前沿Fig.3 Pareto fronts obtained by three test optimization examples

圖4給出了3種測試優(yōu)化算例的IGD指標收斂過程。從圖中可以看出，本文算法的IGD指標收斂更快，對于ZDT1與ZDT2，其IGD分別在約150步和110步時幾乎收斂到0。對于ZDT3，本文算法約50步迭代后的IGD值已經(jīng)跟NS-GA2迭代200步得到的值相當。對于具有不同特征的Pareto前沿的解析函數(shù)，優(yōu)化結果表明本文算法的優(yōu)化效率明顯高于NSGA2。

圖4 3種測試優(yōu)化算例的IGD指標收斂過程Fig.4 Convergence histories of IGD obtained by three test optimization examples

3 氣動優(yōu)化算例

選擇低雷諾數(shù)翼型SD7032為基準進行多目標氣動優(yōu)化設計。翼型參數(shù)化采用CST (Class function/Shape function Transformation)方法[24]，翼型上下表面型函數(shù)階數(shù)均取為5，類函數(shù)指數(shù)N1=0.5，N2=1.0，設計變量數(shù)目為12。以SD7032的CST參數(shù)為基準，設計變量變化范圍為±30%。分別考慮2目標和3目標的翼型優(yōu)化設計問題，其中2目標優(yōu)化設計問題定義為

(8)

施加力矩系數(shù)約束為-0.12≤Cm1≤-0.09，-0.12≤Cm2≤-0.09，并約束翼型最大厚度tmax≥9.0%，其中下標1、2分別表示設計狀態(tài)1與設計狀態(tài)2的氣動參數(shù)。約束采用罰函數(shù)方式處理。3目標優(yōu)化問題在2目標優(yōu)化問題基礎上增加1個設計點，定義為

(9)

采用與2目標問題相同方式處理約束。翼型氣動性能求解采用XFOIL，XFOIL在低速低雷諾數(shù)范圍內(nèi)能較精確高效預測翼型氣動特性[25]。

同樣分別利用MR-EDA和NSGA2算法進行優(yōu)化。對于2目標優(yōu)化，NSGA2初始種群數(shù)目N=100，最大迭代步數(shù)設為100，交叉概率0.8，變異概率0.3。MR-EDA中流形重構種群數(shù)目N1=N/2，其余參數(shù)與NSGA2一致。3目標優(yōu)化問題種群數(shù)目增加為300，其余參數(shù)不變。

圖5給出了低雷諾數(shù)翼型的2目標優(yōu)化Pareto前沿的收斂過程，其中圖5(a)、圖5(b)分別給出了MR-EDA、NSGA2的Pareto前沿收斂過程，圖中的迭代步數(shù)均為5、10、20、100。從圖中可以看出，僅迭代5步，本文所提MR-EDA得到的Pareto前沿已經(jīng)形成了較為規(guī)則的分布，在之后的迭代中，Pareto前沿分布更加趨向于規(guī)則分布。迭代20步與迭代100步的Pareto前沿基本重合，說明了該算法具有很高的收斂效率。作為對比，NSGA2收斂更為緩慢，在相同的迭代步數(shù)下，Pareto前沿收斂性稍差。盡管本文算法最大迭代步數(shù)設置為100，但實際上由于算法的收斂效率高，在經(jīng)過約20步的迭代后，Pareto前沿幾乎與NSGA2迭代100步得到的最優(yōu)Pareto前沿幾乎重合。

圖5 2目標翼型優(yōu)化Pareto前沿收斂過程Fig.5 Convergence histories of Pareto front of 2-objective airfoil optimization

圖6(a)給出了3目標優(yōu)化迭代20步時的Pareto前沿。從圖中可以明顯看出，MR-EDA的收斂速度明顯快于NSGA2。迭代100步得到的Pareto前沿在圖6(b)中給出，2種算法都已經(jīng)收斂到相當?shù)乃健７謩e對比圖6(a)和圖6(b) 中的Pareto前沿，對于MR-EDA迭代20步的解與迭代100步的解已經(jīng)十分接近，即本文算法減少了約80%的計算量，說明了本文算法具有高效收斂的特點。

選取3目標優(yōu)化最優(yōu)Pareto解集中的一個翼型進行氣動性能分析，所選翼型在雷諾數(shù)為3.6×105時的功率因子約為117。圖7給出了該優(yōu)化翼型與SD7032翼型的形狀。優(yōu)化后的翼型前緣半徑減小，前半部分厚度減小而后半部分厚度增大，上表面曲率過渡更平緩。優(yōu)化翼型厚度約為9.13%，在約束范圍內(nèi)。圖8給出了該翼型在設計點2處的壓力分布，并與SD7032翼型的壓力分布進行了比較。從圖中看出，優(yōu)化后的翼型轉捩點后移，在后半段負壓幅值更高，這就使得升力增加，同時低頭力矩增大。

表2給出了優(yōu)化翼型與SD7032在3個設計點的氣動性能。從表中看出，優(yōu)化翼型的升力系數(shù)在3個設計點下都有所提升，增加約8%，同時阻力均有所減小。功率因子在3個設計點分別增加了約20.2%、31.7%、27.8%，優(yōu)化效果明顯。力矩系數(shù)的在設計點均明顯減小，但未超出約束范圍。圖9給出了該優(yōu)化翼型與SD7032翼型功率因子隨雷諾數(shù)變化的關系。從圖中可以看出在整個優(yōu)化區(qū)間內(nèi)，功率因子均有明顯提升，這就保證了翼型性能在一定速度范圍內(nèi)平緩過渡。

圖6 3目標翼型優(yōu)化得到的最優(yōu)Pareto前沿Fig.6 Pareto front of 3-objective optimal airfoil design problems

圖7 3目標優(yōu)化得到的Pareto解集中的翼型形狀與SD7032翼型Fig.7 Optimal airfoil shape in Pereto set obtained by 3-objective design and SD7032 airfoil

圖8 3目標優(yōu)化得到的Pareto最優(yōu)解壓力分布與SD7032翼型壓力分布(Re=3.6×105, α=4°)Fig.8 Pressure coefficient optimal airfoil in Pereto set obtained by 3-objective optimal and SD7032 airfoil (Re=3.6×105, α=4°)

表2 優(yōu)化翼型與SD7032氣動性能比較Table 2 Aerodynamic performance of optimal airfoil and SD7032 airfoil

圖9 3目標優(yōu)化得到的Pareto最優(yōu)解與SD7032的翼型功率因子隨Re變化曲線(α=4°)Fig.9 Variation curves of PI with Re for optimal airfoil in Pareto set obtained by 3-objective design and SD7032 airfoil (α=4°)

4 結 論

本文提出了一種采用流形結構重建的分布估計方法，應用于多目標優(yōu)化算法。通過解析函數(shù)及翼型氣動多目標優(yōu)化算例，驗證的本文方法的高效收斂特性，有以下結論：

1) 本文算法充分利用了目標函數(shù)空間Pareto前沿的流形結構，以此建立了解集的分布概率模型，用于指導算法演化，避免了傳統(tǒng)進化算法盲目交叉變異的缺點，顯著提升了優(yōu)化效率。

2) 通過解析函數(shù)算例及2目標、3目標氣動優(yōu)化算例，比較了NSGA2與本文所提算法。結果表明本文算法收斂更快，其中翼型3目標優(yōu)化設計能夠減小約80%的計算量。

3) 流形結構信息能夠很好地指導算法演化，充分利用流形結構信息能夠極大提高多目標設計問題的優(yōu)化效率。