基于深度置信網絡與多目標粒子群算法的通用飛機機翼優化設計

吳涌釧，孫 剛，陶 俊

(復旦大學 航空航天系，上海 200433)

0 引 言

通用飛機是民用航空工業中重要的組成部分，在工業、農業、林業等領域有重要的應用，在醫療衛生、搶險救災、氣象探測等方面也有不可替代的作用。升力的提升可以提高通用飛機的載貨量、載油量，從而提升飛機的整體性能，因此，升力的優化對通用飛機的設計有重要意義。

國內外研究人員對機翼的設計進行了大量的研究。Benaouali[1]應用代理模型方法對ONERA M6 進行機翼外形優化。馬曉永[2]應用集成序列二次規劃（sequential quadratic programming, SQP）優 化 算 法 對NAX880 V2.0 輕中型高速渦扇公務機機翼進行多點優化設計研究。周志強[3]開展了基于遺傳算法的飛機機翼結構拓撲優化設計方法研究。郝璇[4]采用雷諾平均N-S（Reynolds average Navier-Stokes, RANS）方程實現阻力的精確求解并建立響應面模型，對不同升力系數、馬赫數的多個飛行狀態進行變彎度減阻優化。

對于傳統的機翼優化設計問題，如果在優化過程的每步迭代中都通過CFD 數值模擬來評估目標函數，會消耗大量的計算資源。為了提高優化效率，減少資源的消耗，可以使用代理模型取代CFD 在優化中的計算過程。反向傳播（back propagation, BP）神經網絡[5]是一種由誤差反向傳播算法訓練得到的多層前饋神經網絡，是當前被應用最廣泛的神經網絡模型之一。呂小龍等[6]利用本征正交分解 （proper orthogonal decomposition, POD）提取本征向量，并使用徑向基函數（radical basis function, RBF）方法進行插值，建立了具有較好反演精度和泛化性能的PODRBF 代理模型。此外，還有一種Kriging 模型[7]，是依據協方差函數對隨機過程/隨機場進行空間建模和預測（插值）的回歸算法。

在機翼的優化問題中，優化變量往往是高維數據，如果用原數據去進行訓練和優化，則會造成較大的計算資源消耗，增加計算時間成本。因此，許多研究會采用數據降維方法對原數據進行處理。獨立分量分析（independent component analysis, ICA）[8]基于信息理論，是使用最廣泛的降維技術之一，可用于將多元信號分離為加性子分量。POD[9]則可以把多維隨機過程進行低維近似描述并提取復雜隨機過程的本質特征，從而達到數據降維的目的，可用于外形氣動設計。線性識別分析（linear discriminant analysis,LDA）[10]使用統計學、模式識別和機器學習方法，在兩類物體或事件的特征中找到一個線性組合，從而能夠特征化或區分它們，所得的組合可用來為后續的分類做降維處理。

飛機通常在多個工況下飛行，這就使得機翼的氣動優化成為多目標、多約束問題。多目標多約束問題的求解依賴優化算法，為此，許多研究人員對優化算法進行了研究。多目標遺傳算法（multi-objective genetic algorithm, MOGA）[11]是一種模擬大自然中生物體進化規律的自然選擇和遺傳學機理的生物進化過程的計算模型，對于一個優化問題，使一定數量候選解（個體）的抽象表示（染色體）的種群向更好的解進化。模擬退火算法（simulated annealing, SA）[12-13]的思想來源于熱力學中對固體退火過程的模擬，其將溫度T模擬為控制參數，將內能模擬為目標函數值f，對于控制參數T的每一個取值，算法持續進行“產生—新解—判斷—接受或舍棄”的迭代過程，并逐步衰減T值，算法終止時的當前解記為所得近似最優解。多目標蟻群算法（multi-objective ant colony optimization, MOACO）[14-15]是模擬螞蟻的覓食行為與旅行商問題（traveling salesman problem，TSP）的相似性提出的，用螞蟻的路徑表示待優化問題的可行解，路徑較短的螞蟻釋放的信息素量較多，隨著時間的推進，較短的路徑上累積的信息素濃度逐漸增高，選擇該路徑的螞蟻個數也愈來愈多，最終蟻群集中到最優路徑，即優化問題的最優解。

綜上所述，機翼在多工況下的空氣動力學優化作為一個多目標、多約束的問題，具有重要意義。本文提出了一種基于深度置信網絡代理模型和多目標粒子群算法的多目標優化框架。首先，利用基于型函數/類函數變換參數化方法（class function/shape function transformation, CST）得到的幾何參數，利用CFD 模擬得到空氣動力參數，然后對機翼的幾何參數進行主成分分析（principal component analysis, PCA）降維處理，建立起機翼的數據庫，在此基礎上訓練深度置信網絡（deep belief network, DBN）代理模型，以準確預測機翼的氣動參數；其次，通過將優化算法與代理模型耦合，建立起多目標粒子群優化框架；最后，通過該框架對機翼進行多工況下的氣動優化，在不增加其阻力系數的情況下提升其升力系數。

1 數值模擬方法與驗證

1.1 RANS 方法

通過求解RANS方程，可得到機翼的氣動參數：

數值模擬中，空間離散采用Roe-FDS[16]方法，時間推進采用隱式LU-SGS[16]方法，湍流模型為方程SA （Spalart-Allmaras）[17]湍流模型。

1.2 算例驗證

為了驗證RANS 方法的精準度和可靠性，使用DLR翼身組合體進行驗證，圖1 展示了該翼身組合體的網格結構，網格單元數為6.34 × 106。

圖1 DLR-F6 翼身組合體網格Fig.1 Grid of DLR-F6 wing-body assembly

數值模擬計算結果如表1 所示?？梢钥闯觯瑪抵的M計算結果與實驗所得結果誤差較小，說明RANS 具有較高的精確度。

表1 DLR-F6 翼身組合體的數值模擬結果(Re = 3.0 × 106)Table 1 Numerical simulation results of the DLR-F6 wing-bodyassembly(Re = 3.0 × 106)

在對機翼進行優化時，需要對機翼的幾何數據進行處理，用盡量少的特征參數描述機翼的幾何特征，以降低計算量，同時必須具有足夠高的精度，以確保參數對幾何特性描述的準確性及計算的可信度。

本研究所優化的機翼半翼展長為8.305 m，翼根弦長為2.180 m，翼梢弦長為0.910 m。選取機翼在翼根、翼尖梢以及1/2 半展長處的剖面翼型（分別記為剖面翼型1、剖面翼型2、剖面翼型3），以其幾何數據作為機翼的幾何數據。圖2 為機翼的平面形狀圖，圖3 為選取的三個剖面的示意圖。

圖2 機翼的平面形狀示意圖Fig.2 Schematic figure of the plane shape of the wing

圖3 機翼三個剖面的示意圖Fig.3 Schematic figure of the three sections of the wing

圖4 原翼型曲線和CST 擬合翼型曲線Fig.4 Original airfoil curve and CST fitted airfoil curve

圖5 CST 參數化擬合翼型殘差Fig.5 Fitting residuals of the CST parameterized airfoil

使用CST 參數化方法[18]處理幾何數據。CST 方法是一種用類別函數和形狀函數表示翼型幾何形狀的參數化方法。類別函數表示的是要描述的幾何體的類型，形狀函數是在類別函數的基礎上對幾何特性進行修正。翼型用CST 參數化方法進行表示的公式如下：

2 優化問題

2.1 CST 參數化

型與原始翼型基本一致，殘差均小于1×10?2，說明采用的CST 方法在描述翼型幾何形狀方面具有很高的可靠性。經過CST 參數化后，機翼的幾何特性可以用42 個參數表示。

2.2 PCA（主成分分析）

參數化后的機翼幾何參數有42 個變量，這對于后續的代理模型訓練以及優化過程是偏高的，容易造成較大的計算資源損耗。為了讓計算更加高效便捷，對數據進行降維處理。

PCA[19]是一種基礎的數據降維方法，其作用是用較少的變量來代替較多的變量，并反映出原來多個變量的大部分信息，被廣泛應用在數據分析、機器學習等方面。

原機翼包含了3 個翼型的幾何數據，每個翼型數據為14 維，共42 維?，F對每一個翼型進行PCA 降維處理，當每個翼型取前10 個主成分時，累計貢獻率分別達到了93.9%、93.4%、94.04%，則可以認為處理后所得的30 維機翼幾何數據包含了絕大多數原機翼的數據。圖6 展示了三個剖面翼型的14 維參數貢獻率。

圖6 3 個剖面翼型的14 維參數的貢獻率Fig.6 Contribution rates of the 14 dimensional parameters of the three airfoils

2.3 優化問題描述

本研究要優化的是機翼在多種工作狀態下的升力系數。選取其在馬赫數為0.4 時的狀態作為優化工況1，在馬赫數為0.45 時的狀態作為優化工況2，迎角均為4°，工況1 和工況2 的雷諾數分別為Re1=9.83 × 106和Re2= 1.11 × 107。

基于5 套網格對原機翼開展了網格無關性驗證，5 套網格單元數分別為2.0 × 106、2.5 × 106、3.0 × 106、3.5 × 106、4.0 × 106，計算結果如表2 所示，其中，各參數的下標1 和2 分別對應工況1 和工況2。由計算結果可得，伴隨網格數增長，計算結果呈收斂趨勢，從計算精度和計算資源考慮，選用網格單元為3.5 ×106的網格作為計算網格，網格如圖7 所示。

表2 原機翼氣動參數在不同網格數量下計算結果Table 2 Simulation results of aerodynamic parameters oforiginal wing under different numbers of grids

圖7 機翼的網格Fig.7 Grid of the wing

以兩種工況下的機翼升力系數作為優化目標，約束條件為優化后機翼阻力系數小于或等于原機翼阻力系數。優化問題可以描述為：

3 優化方法

3.1 深度置信網絡

本研究使用的代理模型DBN[20]是一種概率生成模型，這種模型通過訓練調整神經元之間的權重，可以用來進行數據生成，常用于回歸分析、類型識別等。

深度置信網絡的基本構成元件是受限玻爾茲曼機（restricted Boltzmann machines, RBM）[21]，其典型網絡結構如圖8 所示。一個RBM 由兩層神經元構成，一層是可視層，由可視神經元組成，用于數據的輸入；一層是隱藏層，由隱藏神經元組成，用于特征檢測。RBM 兩層神經元之間對稱連接，同層之間的神經元相互獨立。

圖8 RBM 結構Fig.8 Structure of RBM

對于輸入樣本X=(x1,x2, ···,xn)，隱藏神經元開啟的概率為P(hi=1)=1/(1+e?wx)，從0 ～ 1 的均勻分布中抽取一個隨機值，與開啟概率進行比較，決定該神經元是否開啟。經過預訓練，可求得最佳的權重值，使得其產生訓練樣本的概率最大。

將多層RBM 進行串聯，以上一層RBM 的隱藏層作為下一層RBM 的可視層，在最后一層加上softmax 分類器，即形成了深度置信網絡。這種網絡將RBM 的訓練作為神經網絡的預訓練，在最后一層加入BP 算法。深度置信網絡典型結構如圖9 所示。

圖9 DBN 結構Fig.9 Structure of DBN

3.2 改進多目標粒子群優化算法

粒子群優化（particle swarm optimization, PSO）算法[22]是一種通過模擬鳥群捕食行為設計的優化算法。假設區域里只有一塊食物（即優化問題中的最優解），鳥群的任務是找到這個食物源。在整個搜索過程中，群體中的鳥會互相傳遞位置信息，判斷是否找到了食物源（最優解），同時將食物源（最優解）的信息傳遞給整個鳥群，最后，整個鳥群可以收斂在食物源周圍。

把鳥個體簡化為一個粒子，把粒子的位置作為目標問題的候選解，把候選解對應的目標函數值作為適應度，把粒子群中適應度最高/最低的粒子位置作為群體最優位置，把單個粒子歷史適應度最高的位置作為粒子歷史最優位置，粒子的飛行過程即對應個體的搜索過程。粒子的飛行速度以粒子歷史最優位置和群體最優位置為參考進行調整，此即為進化迭代過程。經過多輪迭代，可以求出問題的最優解。粒子的速度和位置的進化方程如下:

多目標粒子群優化（multi-objective particle swarm optimization, MOPSO）[23]算法旨在將只能用于單個對象的粒子群算法（PSO）應用于多目標。在其基礎上，Raquel[24]將粒子擁擠距離、儲存非支配解的倉庫等引入其中，來提高粒子群的多樣性。擁擠距離是估計一個解周圍其他解的密度的參數。對于每個目標函數，根據目標函數的值對外部精英解集中的解進行分類，然后計算與每個解相鄰的兩個粒子形成的立方體的平均邊長，如圖10 所示，計算的結果就是該解的擁擠距離，粒子的擁擠距離可以反映解的多樣性。該種算法即為擁擠距離多目標粒子群優化（crowding distance multi-objective particle swarm optimization,CDMOPSO）算法[25]。

圖10 擁擠距離示意圖Fig.10 Schematic figure of the crowding distance

為了提高收斂速度和全局搜索能力，引入了αstable 函數族[26]來改進標準的MOPSO 算法。在αstable 函數族中，穩定系數α是描述分布狀態的一個重要參數，它決定了所生成的隨機數的分布范圍。圖11 顯示了不同α值下的α-stable 分布。通過引入α-stable 分布，MOPSO 在迭代過程中生成隨機數，使粒子進行突變，避免過早陷入局部最優。該種算法即為α-stable 多目標粒子群優化（α-stable multi-objective particle swarm optimization, ASMOPSO）算法。多目標粒子群算法如圖12 所示。

圖11 不同α 值下的α-stable 分布示意圖Fig.11 Schematic figure of the α-stable distribution under different α

圖12 MOPSO 結構Fig.12 Structure of MOPSO

為了比較CDMOPSO 和ASMOPSO 的性能，采用目標函數ZDT3 進行測試。利用反世代距離（inverted generation distance, IGD）[27]對算法的性能進行了綜合評價。IGD 不僅可以評估多目標解集的收斂性，還可以評估解集的分布特性，見式（5）：

d(o,P) 是近似解集P中的個體與精確解集P*中的個體o之間的最小歐氏距離。該算法得到的最優解集用p表示。IGD 的值越低，優化結果越好，算法的性能也越好。ZDT3 函數可以描述為：

x是 一 個30 維 變 量，其 范 圍 為[0，1]，n= 30。ZDT3 由CDMOPSO 算法和ASMOPSO 算法獨立評估30 次。IGD 值的統計數據如表3 所示，可以看出ASMOPSO 的IGD 值 明 顯 低 于CDMOPSO， 因 此ASMOPSO 算法的性能優于CDMOPSO 算法。

表3 CDMOPSO 和ASMOPSO 算法所得IGD 的比較Table 3 Comparison of IGD values between the CDMOPSOalgorithm and ASMOPSO algorithm

3.3 EI 加點準則

在優化過程中，每次迭代均采用期望改進（expected improvement, EI）[28]加點準則對代理模型進行更新。將EI 作為代理優化的適應度函數，在每次迭代中將EI 最大的點添加到樣本中，在每次迭代中更新代理模型。多目標EI 的定義為

在每次迭代中，選擇EI 最大的點進行CFD 數值模擬，獲得空氣動力學參數，添加到樣本中，并更新替代模型。

4 計算結果與分析

4.1 代理模型訓練結果

基于原機翼的CST 參數，通過拉丁超立方采樣方法 （latin hypercube sampling, LHS）[29]，在原機翼參數[0.8, 1.2]的倍數范圍內，生成980 組CST 參數，并通過數值模擬獲得一個包含980 個機翼的幾何參數和氣動參數的數據庫。

DBN 由4 層RBM 組成，以幾何參數為輸入、氣動參數為輸出進行訓練。預訓練后，初始化權值，將最后一層RBM 的輸出作為最后兩層BP 神經網絡的輸入。

在數據庫中選取882 個機翼作為訓練樣本，98 個機翼作為測試樣本，將樣本輸入模型進行訓練，得到的代理模型記為代理模型1。在數據庫中隨機選取450 個機翼作為訓練樣本，50 個機翼作為測試樣本，為了降低后續的訓練和優化計算量，對機翼幾何數據進行PCA 降維處理，降維后的幾何參數為30 維，將PCA 降維處理后的樣本輸入模型進行訓練，得到的代理模型記為代理模型2。

圖13 為DBN 代理模型的訓練結果，其中Y為期望輸出，T為DBN 輸出。可以看出，代理模型1 的相關系數R為0.995 16，代理模型2 的相關系數R為0.974 02，均接近于1，兩個模型都取得了較為精確的結果。

圖13 DBN 代理模型的訓練結果Fig.13 Training results of the DBN surrogate model

圖14 和圖15 分別展示了測試樣本的計算值與訓練所得代理模型預測值的對比。從圖中可看出，絕大多數預測值與計算值基本吻合。

圖14 機翼氣動參數代理模型1 的預測值與計算值Fig.14 Surrogate model 1 prediction values and calculation values of the aerodynamic parameters of wings

圖15 機翼氣動參數代理模型2 的預測值與計算值Fig.15 Surrogate model 2 prediction values and calculation values of the aerodynamic parameters of wings

DBN 模型預測的升力系數和阻力系數的平均誤差如表4 所示。由表可以看出，代理模型1 的升力系數與阻力系數的平均誤差小于1%，代理模型2 平均誤差小于2%，說明DBN 模型的預測具有較高的精度。代理模型1 的誤差明顯小于代理模型2，由此可得，代理模型1 具有更高的精準度，原因是代理模型1 的訓練樣本更多且保留了更多原數據的信息。

表4 升力系數和阻力系數的平均預測誤差Table 4 Average prediction errors of the lift coefficients and thedrag coefficients

4.2 多目標粒子群優化結果

將訓練得到的無PCA 處理和有PCA 處理的代理模型嵌入多目標粒子群優化算法中，分別記為優化1 和優化2。

將粒子數設置為500，Pareto 解數設置為10，突變率設置為0.8。優化的終止條件是最近10 代的適應度值的最大差異均小于0.01%。在每次迭代中，如果滿足終止準則，則優化結束，否則迭代將繼續。

優化1 與優化2 分別經過91 步與102 步迭代后，收斂得到最優解。優化1 與優化2 的優化過程對比如圖16 所示，優化得到的Pareto 解集如圖17 所示。

圖16 優化1 與優化2 的優化過程Fig.16 Optimization process with/without PCA processing

圖17 優化1 與優化2 優化所得Pareto 解集Fig.17 Pareto solution set obtained by optimization with/without PCA processing

為了在每次迭代的Pareto 解集中選出一個最優解，現定義一個函數：

其中ω1、ω2是權重系數，在此處設為[0.5，0.5]。對Pareto 解集中的所有解計算出對應的fop值，fop值最大的解即為滿意解。

4.3 優化結果驗證

優化后機翼的幾何模型如圖18 所示。

圖18 優化機翼的幾何模型Fig.18 Geometric model of the optimized wing

對優化機翼進行數值模擬，得出優化后機翼在工況1 與工況2 下的升力和阻力系數計算值。原機翼和優化機翼的氣動參數在DBN 模型下的預測值及計算值如表5 所示。

表5 原機翼與優化機翼氣動參數的預測值與計算值Table 5 Prediction and calculation values of aerodynamicparameters of the original and optimized wings

由表可以看出，優化1 得到的優化機翼較原翼型的升力系數在工況1 下提升9.81%、工況2 下提升10.14%，阻力系數在工況1 下減少1.04%、工況2 下減少0.20%；優化2 得到的優化機翼較原翼型的升力系數在工況1 下提升9.44%、工況2 下提升9.71%，具有較明顯提升，阻力系數在工況1 下減少0.69%、工況2 下不變，優化效果較為理想。

可以看到，優化1 與優化2 的優化結果基本相近，為了進一步比較兩種優化的效率，表6 給出了以無代理模型時過程需要100 步迭代達到收斂計算，優化1、優化2 和無代理模型三種情況下的計算算例數量和計算時間。

表6 優化過程中 CFD 數值模擬算例數量和計算時間Table 6 Computation time and number of CFD simulation casesduring the optimization process

由表可知：無代理模型情況下優化過程中所需要數值模擬的算例數量遠遠大于有代理模型；有代理模型的情況下，優化2 的優化過程需要計算591 個算例，遠低于優化1 所需的1 082 個算例，而兩者所達到的優化效果是基本接近的。

5 結 論

本研究提出了一種基于DBN 的多目標優化框架，并應用于某通用飛機機翼的氣動優化設計，研究表明：

1）采用CST 方法對機翼的控制剖面進行參數化描述，在此基礎上，引入PCA 技術，有效地降低了機翼的設計變量。

2）建立了機翼的氣動數據庫，針對是否采用PCA 技術分別建立了DBN 代理模型，訓練結果表明所建立的DBN 代理模型對氣動參數具有很好的預測效果。

3）基于DBN 代理模型，采用MOPSO 算法對于是否采用PCA 技術分別開展了機翼的多目標氣動優化設計，結果表明兩種優化結果基本一致，在阻力系數不增加的前提下，升力系數明顯提升。然而，采用PCA 技術可大幅減少優化過程中所需的計算量。

在后續的研究工作中，有必要開展更為高效的降維技術研究，進一步減少設計變量，從而減少氣動優化設計過程中的計算量，提高優化效率。