基于多輸出高斯過程的超臨界翼型優化

吳寬展，劉學軍，＊，呂宏強

（1.南京航空航天大學計算機科學與技術學院，江蘇南京 210016；2.南京航空航天大學航空宇航學院，江蘇南京 210016）

基于多輸出高斯過程的超臨界翼型優化

吳寬展1，劉學軍1，＊，呂宏強2

（1.南京航空航天大學計算機科學與技術學院，江蘇南京 210016；2.南京航空航天大學航空宇航學院，江蘇南京 210016）

提出了一種基于多輸出高斯過程（MOGP）回歸模型和模擬退火算法相結合的翼型優化方法。采用拉丁超立方抽樣方法在設計空間內構造一系列樣本點，其中優化設計采用CST參數化方法對翼型的幾何外形進行參數化表示。通過CFD流體計算得到其響應值來建立初始的MOGP代理模型。以阻力最小化為設計目標，考慮面積、升力等約束條件。通過單點優化和多點優化試驗表明，發展的翼型優化設計方法達到了優化設計的目的，同時也說明基于MOGP模型的優化設計方法在氣動優化設計中的應用是可行的。

翼型優化設計；輸出高斯過程；CST參數化方法；模擬退火算法；超臨界翼型

0 引 言

隨著計算機計算能力的提升與計算技術的發展，CFD（計算流體動力學）軟件已經成為流場分析、流場計算、流場預測的重要工具。在翼型氣動優化過程中，由于需要在翼型設計空間大量調用高精度流場計算模擬，基于CFD的優化計算會受到時間、財力的限制。一個解決方案是采用代理模型［1］技術（又稱響應面技術）構建近似模型來代替復雜的CFD計算。近似模型描述了目標函數（輸出）與設計變量（輸入）之間對應關系。通過近似模型計算目標函數值節省了大量計算時間并且可以搜索更廣闊的設計空間。一些試驗表明基于代理模型的優化技術不僅高效而且可靠，已逐漸在工程上得到應用［2－5］。

翼型設計領域中常用的代理模型主要有多項式模型，人工神經網絡模型，Kriging模型（單輸出高斯回歸模型）［3－6］等。高斯過程回歸模型本質上是一種基于貝葉斯框架的回歸方法，有著嚴格的統計學理論基礎。高斯模型能很好地處理小樣本、非線性等復雜的問題。利用統計學方法來尋找輸入變量與輸出變量之間的條件分布，從而建立兩者的回歸方程。高斯回歸模型分為單輸出高斯過程（Single Output Gaussian Process，即Kriging）回歸模型［7］和多輸出高斯過程（Multiple Output Gaussian Process，MOGP）回歸模型［8－9］。單輸出模型對每個輸出端口分別建模，忽略了端口之間的潛在相關性。MOGP回歸模型則為各個輸出端口之間建立了協方差矩陣來模擬它們之間可能存在的依賴關系，每個端口都可以利用其它輸出端口來提升自己的預測性能。

在先前的工作中，我們驗證了MOGP模型相比Kriging模型在多響應之間存在顯著相關性時，能進一步提高代理模型在超臨界翼型設計中的預測精度［10］。本文將CST參數化方法、MOGP多響應代理模型與模擬退火算法相結合設計了一種翼型優化設計方法，并將該方法應用到超臨界翼型的優化設計中。通過單點和多點兩個翼型優化算例證明了提出方法的有效性。

1 多輸出高斯過程模型

MOGP設含有D個函數的集合，｛fd（x）｝Dd＝1，其中每一個函數都可以被表示成一個光滑核函數，｛Gd（x）｝Dd＝1，和一個隱含函數u（x）的卷積，即：

更一般的，可以認為每一個光滑核函數和隱函數都是由若干個子函數線性組合而成，

每一個子函數都擁有相同的協方差函數，且相互獨立。

多輸出高斯過程卷積的光滑核函數和隱函數的選擇有多種，但在實際應用中，一般采取高斯形式，使光滑核函數和隱函數擁有相同的高斯形式，其可行性已由Boyle和Frean證明［8］，光滑核函數表示為：

其中，Sd，q為依賴于輸出端口fd（x）和隱函數μq（x）的方差系數，P－1d為輸出端口的精度矩陣。

隱含過程的協方差函數可以表示為：

在上面的表達式中，Λq是關于隱函數的精度矩陣。于是我們可以得到任意兩個輸出端口的協方差表達式：

給定卷積形式之后，模型的似然函數表示為：

其中，y為所有輸出端口的列向量；Kf，f為卷積過程的協方差矩陣，它的每個元素都是式（5）計算得到的；Σ是一個獨立的高斯白噪聲過程矩陣，是為每個輸出端口所加的一個擾動；θ表示該模型的超參數集合，使用最大似然估計計算出θ的值，記為θ＾。

給定新的輸入集合X＊后，模型輸出的預測分布為：

由于式（7）的協方差矩陣Kf，f表示的是所有輸出數據的協方差，所以MOGP可以為多輸出建模，能夠模擬多輸出端口之間的相關性。

2 CST翼型參數化方法

CST方法［11－13］是近年提出的一種可以使用較少的設計變量來精確描述復雜幾何外形的翼型參數化方法。該方法通過類函數定義了要描述的一般的幾何外形，型函數用來描述給定類函數的特性，具有設計變量可調節、設計空間廣等優點。本文采用CST參數化方法對超臨界翼型進行建模。

用CST參數化方法表示翼型的數學公式為：

其中c為翼型弦長，Δzte為翼型的后緣厚度，C（x／c）為類函數，S（x／c）為型函數，分別用以下式子表示：

對于C（x／c），當N1、N2取不同的值時，可以定義不同分類的基本幾何形狀。本文中取N1、N2分別為0.5和1.0，方程（8）右邊的第一項和第二項分別表示翼型的圓前緣和尖后緣，以保證參數化在翼型的軸向極值處具有很好的數學特性。

在方程（10）中，多項式系數bi｛｝表示了對給定翼型的設計變量，稱之為Bernstein系數，用以生成更復雜更實用的翼型，可通過最小二乘法擬合翼型得到。

3 CFD分析方法

流場控制方程采用守恒形式的N-S方程。N-S方程空間離散采用有限體積方法，對流通量采用Roe通量差分格式，湍流模型采用了κ－ε模型。采用粘性網格，網格在整個流場范圍內的規模為387×104個單元，圖1表示了RAE2822翼型網格拓撲結構。圖2表示了該網格的Y＋分布。

為了驗證計算程序的準確性，圖3給出了RAE2822翼型馬赫數0.729，迎角2.31°狀態下計算與試驗壓力分布對比?？梢钥闯?，CFD計算與試驗結果吻合較好，能夠準確捕捉到上翼面的激波。

圖1 RAE2822的網格圖示Fig.1 Grid for RAE2822

圖2 翼型RAE2822網格的Y＋分布Fig.2 Y＋distribution of RAE2822

圖3 壓力分布比較Fig.3 Comparison of pressure coefficient distributions

4 模擬退火算法

模擬退火算法［14］（Simulated Annealing，SA）的思想來源于熱力學中對固體退火過程的模擬，將固體加溫至融化，再讓其徐徐冷卻凝固成晶體。加溫時，固體內部粒子隨者溫度升高處于不穩定狀態，內能增大，而溫度緩慢下降時，粒子內能減小，固體內部狀態趨于穩定。在實際應用中，將溫度t作模擬為控制參數，將內能模擬為目標函數值f，對于控制參數t的每一個取值，算法持續進行“產生——新解——判斷——接受或舍棄”的迭代過程，并逐步衰減t值，算法終止時的當前解記為所得近似最優解。其中，使用Metroplis準則作為判斷當前解是否為新解替換的標準。這是基于蒙特卡羅迭代算法的一種啟發式隨機搜索過程。模擬退火算法適合處理傳統搜索方法難以解決的高度復雜的非線性問題。本文采用模擬退火算法作為翼型設計的優化算法。

5 翼型優化系統

5.1 試驗設計

試驗設計的目的是選擇一些有代表性的樣本點從而利用少量試驗點得到較高精度的代理模型。試驗設計的方法有多種，但主要有正交試驗設計、均勻設計、中心復合設計和拉丁超立方設計等。拉丁超立方設計［15］（Latin Hypercube Sampling，LHS）是滿足空間試驗設計中一種較為優秀的試驗設計方法，現在已經獲得了廣泛的應用。因此本文選用拉丁超立方設計來選取所需要的樣本點。

5.2 翼型優化設計算法框架

基于MOGP多響應代理模型的翼型優化設計方法的步驟為：首先確立使用數學語言描述的優化目標，通過CST參數化方法確定合理的翼型設計空間（設計變量的個數及其取值范圍）；然后采用拉丁超立方抽樣在設計空間內選取合適的初始樣本點，并使用CFD軟件計算出初始翼型所對應的響應值（升力系數、阻力系數等）；根據初始樣本點信息構造MOGP模型；根據初始樣本點構造MOGP代理模型。對代理模型進行優化分析，翼型的幾何約束采用罰函數法［16］處理，即對違反約束的個體適應度施加一個懲罰項，使用模擬退火算法尋找最佳的設計變量取值，對搜索到的局部最優值，將它加入到初始樣本點中，重新構造模型，如此迭代計算，直到收斂到全局最優點。圖4是基于MOGP模型與模擬退火算法的優化流程圖。

圖4 優化流程圖Fig.4 Flow chart of optimization

6 算例分析

6.1 單點翼型優化設計

算例1。以超臨界翼型RAE2822為初始翼型，以面積、升力等作為約束條件，進行單目標減阻優化設計，設計狀態為迎角2.79°，馬赫數0.734，雷諾數6.5×106。優化問題可以描述為：

Minimize CD

Subject to CL＝0.7983

Area≥0.99Area0

在優化設計過程中，取12維Bernstein系數作為設計變量，在給定的設計變量約束范圍內，采用拉丁超立方抽樣選取121個初始試驗點，建立初始MOGP模型，通過模擬退火算法得到單次搜索最優解并加入原樣本集中，經過21次循環后得到最終最優解。表1給出了設計變量范圍與最終優化結果，表2給出了基準翼型和優化后翼型的氣動力參數比較。由表可知，在滿足約束條件的前提下，優化后的翼型阻力系數比基準翼型減小了31%。圖5和圖6分別給出了原始外形與最終優化得到的翼型的形狀以及壓力分布的比較。由圖可知，優化翼型有著較小的前緣半徑，并且翼型的激波有著明顯的減小。

表1 設計變量約束范圍及優化結果Table 1 Scope of design variable constraint and optimization result

表2 初始翼型和優化翼型氣動力參數比較Table 2 Comparison of airfoil character between initial and optimized airfoil

圖5 優化前后翼型形狀Fig.5 Comparison between initial and optimized airfoil configurations

圖6 優化前后翼型表面壓力分布比較Fig.6 Comparison of airfoil pressure distribution between initial and optimized airfoil

6.2 多點翼型優化設計

飛行器需要在一定的馬赫范圍內均有良好的性能。通常只有一個設計狀態的單點優化很難滿足要求，而多目標優化可以有效地將多個設計要求統一起來處理。這里我們對RAE2822翼型在馬赫數0.70和0.75進行兩點優化，并考察優化翼型在馬赫數范圍0.68至0.76內，阻力系數的發散情況。

對于多目標的優化問題，采用“統一目標函數法”［17－19］處理多目標函數。統一目標函數法，即將各個分目標函數統一到一個總的目標函數中去。在兩點翼型優化中，假設第一設計點和第二設計點的目標函數分別為F1（x）、F2（x），使用線性插值法將各個分目標函數規格化為：

其中，αi、βi為目標函數的最大值與最小值。為使各個分目標能均勻一致的趨向各自的最優目標，根據各個目標函數的重要性分別給予一個加權因子ω，于是，統一的目標函數可表示為：

本文中，目標函數F1（x）、F2（x）分別代表待優化翼型在低馬赫數與高馬赫數下的氣動性能參數，其權重是相同的，我們令加權因子ω1、ω2均為0.5。

算例2。設計狀態（1）：迎角2.79°，馬赫數0.70，雷諾數6.5×106；設計狀態（2）：迎角2.79°，馬赫數0.75，雷諾數6.5×106。

Minimize CD1，CD2

Subject to CL1＝0.7526

CL2＝0.7858

算例2采用了242個初始樣本點構造MOGP代理模型，表3、表4給出了優化結果。經過優化后，優化翼型在兩個設計狀態下，阻力系數比基準翼型分別減小了5.5%和23%。

表3 設計狀態1優化結果Table 3 Optimization result of design state 1

表4 設計狀態2優化結果Table 4 Optimization result of design state 2

圖7給出了優化前后的阻力發散特性。從圖7中可以看到，由于單點優化只是局部地優化了一個設計點（Ma＝0.73）的氣動性能，所以優化后的翼型雖然在設計點處和高馬赫區域取得了良好的減阻效果，卻較為嚴重地惡化了低馬赫數條件下翼型的阻力特性。而多目標優化在設計點Ma＝0.70（低馬赫數區域）和設計點Ma＝0.75（高馬赫數區域）進行了全局的優化，所以優化后的翼型在給定的馬赫數范圍0.68至0.76內都有良好的減阻效果，極大地提高了阻力發散特性。

圖7 翼型單點與多點優化阻力特性比較Fig.7 Drag character of airfoil at different optimization result

7 結 論

本文設計了一種基于機器學習算法和傳統優化算法相結合的優化方法，并用于超臨界翼型的優化設計。其核心部分為采用訓練后的MOPG模型作為氣動性能評估的代理模型，優化算法采用的是模擬退火算法。根據所采用的優化算例可以看出，基于MOGP模型的優化設計方法有效地改善了初始基準翼型的氣動性能。從優化設計壓力分布來看，所采用的優化方法能給出合理的翼型，并給出期望的壓力分布。優化結果表明，所發展翼型優化設計方法可以有效地應用于超臨界翼型的多目標氣動優化設計，有著較好的工程應用前景。

［1］ Pena F L，Casás V D，Gosset A，et al.A surrogate method based on the enhancement of low fidelity computational fluid dynamics approximations by artificial neural networks［J］.Computers and Fluids，2012，58：112-119.

［2］Hacioglu A.Fast evolutionary algorithm for airfoil design via neural network［J］.AIAA Journal，2007，45（9）：2196-2203.

［3］ Tao Y，Li Y H，Huang Y.Multi-point optimization of supercritical airfoil based on agent model［J］.Scientific Research，System Simulation Technology ＆Application，2012，（14）：77-81.

［4］Han Z，Grtz S，Zimmermann R.Improving variable-fidelity surrogate modeling via gradient enhanced kriging and a generalized hybrid bridge function［J］.Aerospace Science and Technology，2013，25（1）：177-189.

［5］ Jeong S，Murayama M，Yamamoto K.Efficient optimization design method using Kriging model［J］.Journal of Aircraft 2005，42（2）：413-420.

［6］ Ren Q Z，Song W P.Multi-objective aerodynamic design optimization for airfoil using Kriging model［J］.Aeronautical Computing Technique，2009，39（3）：78-82.（in Chinese）任慶祝，宋文萍.基于Kriging模型的翼型多目標氣動優化設計研究［J］.航空計算技術，2009，39（3）：78-82.

［7］Rasmussen C E，Williams C K I.Gaussian processes for machine learning［M］.Cambridge，MA：MIT Press，2006.

［8］ Boyle P，Frean M.Dependent gaussian processes［C］.Advances in Neural Information Processing Systems 17.Saul L K，Weiss Y，Bottou L，edit.Cambridge，MA：MIT Press，2005：217-224.

［9］ Alvarez M A，Lawrence N D.Computationally efficient convolved multiple output Gaussian processes［J］.Journal of Machine Learning Research，2011，（12）：1425-1466.

［10］Liu X，Zhu Q，Lu H.Modeling multi-response surfaces for airfoil design with multiple output gaussian process regression［J］.Journal of Aircraft，2014，51：740-747.doi：10.2514／1.C032465.

［11］Jacob C H，Willam A C，Gosset A.A parametric approach to supercritical airfoil design optimization［R］.AIAA 2007-62.

［12］Li J，Gao Z H，Huang J T，et al.Aerodynamic optimization system based on CST technique［J］.Acta Aerodynamica Sinica，2012，30（4）：443-449.（in Chinese）李靜，高正紅，黃江濤，等.基于CST參數化方法氣動優化設計研究［J］.空氣動力學學報，2012，30（4）：443-449.

［13］Kulfan B，Bussoletti J E.Fundamental parametric geometry representations for aircraft component shapes［R］.AIAA 2006-6498.

［14］Xia L，Chang Y X，Zhang L.The application of an improved simulated annealing algorithm to airfoil design［J］.Flight Dynamics，2008，26（1）：71-74.（in Chinese）夏露，常彥鑫，張龍.改進的模擬退火算法在翼型設計中的應用［J］.飛行力學，2008，26（1）：71-74.

［15］Bai J Q，Wang B，Sun Z W.The research of robust supercritical airfoil design optimization［J］.Acta Aerodynamica Sinica，2011，29（4）：459-463.（in Chinese）白俊強，王波，孫智偉.超臨界翼型穩健型優化設計研究［J］.空氣動力學學報，2011，29（4）：459-463.

［16］Zhang L，Chen H Q.Research on genetic algorithm for aerodynamic shape optimization based on CST［J］.Aeronautical Computing Technique，2011，41（6）：53-57.（in Chinese）張磊，陳紅全.基于CST參數化的翼型優化遺傳算法研究［J］.航空計算技術，2011，41（6）：53-57.

［17］Zhong X P，Ding J，Li W，et al.Robust airfoil optimization with multi-objective estimation of distribution algorithm［J］.Chinese Journal of Aeronautics，2008，21（4）：289-295.

［18］Ursache N M，Bressloff N W，Keane A J.Aircraft roll enhancement via multi-objective optimization using surrogate modeling［J］.AIAA Journal，2011，49（7）：1525-1541.

［19］Yu G，Li D.Optimization design of airfoil based on hybrid genetic algorithm and compound form method［J］.Sciece Technology and Engineering，2007，7（10）：2292-2295.

A supercritical airfoil design based on multi-output surrogate model

Wu Kuanzhan1，Liu Xuejun1，＊，Lyu Hongqiang2
（1.College of Computer Science and Technology，Nanjing University of Aeronautics and Astronautics，Nanjing 210016，China；2.College of Aerospace Engineering，Nanjing University of Aeronautics and Astronautics，Nanjing 210016，China）

An airfoil optimization method based on a MOGP（Muliple Output Gaussian Process）model and the simulated annealing algorithm is developed.The Latin hypercube method is used to construct a series of sample points in the design space，where the CST parameterization method is used to define the airfoil geometry.The initial MOGP model is established based on the response values obtained by CFD（Computational Fluid Dynamics）calculation.The proposed numerical test is to minimize the drag coefficient with consideration of some constrained conditions such as the area of airfoil，the lift coefficient and so on.Experimental results show that the drag coefficient is significantly reduced and the aera of airfoil and the lift coefficient remains，which is exactly meet the expectation of the optimization.Furthermore，the proposed method can be extended to 3Dcases straightforwardly and is applicable to practical engineering problems.

airfoil design；MOGP；CST method；simulated annealing；supercritical airfoil

V211.3

：Adoi：10.7638／kqdlxxb-2014.0059

2014-06-30；

2014-10-10

國家自然科學基金（61170152）；航空基金重點實驗室類（20151452021）

吳寬展（1988-），男，碩士，研究方向：機器學習.E-mail：514400413＠qq.com

劉學軍＊，教授，研究方向：計算機科學與技術.E-mail：xuejun.liu＠nuaa.edu.cn

吳寬展，劉學軍，呂宏強.基于多輸出高斯過程的超臨界翼型優化［J］.空氣動力學學報，2015，33（6）：728-732.

10.7638／kqdlxxb-2014.0059 Wu K Z，Liu X J，LV H Q.A supercritical airfoil design based on multi-output surrogate model［J］.Acta Aerodynamica Sinica，2015，33（6）：728-732.

0258-1825（2015）06-0728-06