Kriging模型在機翼氣動外形優化中的應用

孫美建，詹 浩

(西北工業大學航空學院，陜西西安 710072)

0 引言

在三維氣動外形優化中，由于設計變量多以及解Euler/N-S方程計算量龐大，基于梯度信息的直接優化方法仍然是飛行器氣動外形優化設計的主要工具［1－3］，然而該方法比較依賴于初始模型，容易陷入局部最優。人們受自然界各種物理現象啟發發展了如遺傳算法、模擬退火算法、粒子群算法等隨機優化算法，與直接優化算法相比具有更好的全局性和魯棒性，但是其缺點是評價目標函數的次數遠遠高于直接優化算法。

基于代理模型的優化方法已經在翼型優化中得到廣泛應用［4－6］，用隨機優化算法與代理模型相結合，可大幅度減少流場計算的次數。本文采用改進的量子粒子群算法優化Kriging模型的相關模型參數，提高Kriging模型預測精度，并結合具有雙層結構的粒子群優化算法，發展了同樣適用于三維氣動外形優化的高效算法。針對具有44個設計變量的跨聲速機翼采用N-S方程流場求解器進行了優化驗證。

1 優化算法

1.1 改進的量子粒子群算法

2004年，Sun J等人提出具有量子行為的粒子群算法(Quantum-Behaved Particle Swarm Optimization，QPSO)［7］。QPSO中粒子搜尋的位置由概率密度函數確定，而不同于標準粒子群算法(Particle Swarm Optimization，PSO)［8］按照軌道的形式進行搜索，通常 QPSO表現出更優秀的全局搜索性能。QPSO算法搜索方程可以表述為:

本文采用進一步改進的量子粒子群算法(MQPSO):每個粒子計算完成后立即更新pbest、gbest、mbest，并根據pbest進行排序，每個粒子僅向優于自身pbest的粒子學習，而且排名越靠后的粒子向其他粒子學習的程度越高。由此通過提高算法對信息的利用率，加快收斂速度。

1.2 Kriging代理模型

Kriging代理模型［10］起源于地理空間統計學，是一種估計方差最小的無偏估計模型，通過相關函數的作用，具有局部估計的特點，它可以較好地預估未知點處函數值的分布情況。Kriging模型中響應值與設計變量之間的關系可以表示為:

其中F(β，x)為回歸模型由P個已知函數組合構成，

回歸模型是對設計空間的全局近似，為確定性部分，通常分成0階(常數)、1階(線性)和2階(二次多項式)三類模型，β為回歸模型參數。

z(x)是均值為0、方差為σ2的統計隨機過程，并且兩個插值點的協方差為:

其中R是點x(i)和x(j)的相關函數，本文選用高斯函數表示為:

未知點x0處的預測值(x0)和方差估計值通過如下形式給出:

通過極大似然估計確定相關模型參數θk:

相關模型參數的選取對代理模型的性能有重要的影響，Lophaven等開發的DACE-A Matlab Kriging Toolbox［11］采用模式搜索方法求解相關模型參數，但是該方法嚴重依賴于初始點的選取，容易陷入局部最優。游龍海等人［12－13］采用遺傳算法和粒子群算法對Kriging模型相關參數進行優化，本文采用MQPSO算法求解相關模型參數，以提高Kriging模型預測精度。

1.3 基于Kriging模型的雙層量子粒子群算法

代理模型方法可以根據已有的點預測未知點的適應值分布情況，當代理模型達到一定精度時，可代替真實求解用于工程優化。傳統的代理模型優化框架［14－16］如下:首先用一定數量的樣本點構造代理模型，優化找到代理模型的最優解并與真實值對比，若不滿足收斂要求則加入新計算的真實值更新代理模型，直到優化結束。該方法收斂速度快，但是嚴重依賴于代理模型預測的精度，對于復雜的高維度優化問題極容易陷入局部最優。

雖然代理模型預測精度有限，但是往往能較好地反映出真實解的大致情況。為了將代理模型獲得的信息充分用于優化，本文采用如下基于Kriging模型的雙層量子粒子群優化算法(KMQPSO)，該算法由兩層MQPSO算法構成，并以Kriging模型為聯系媒介。

首先在變量空間均勻產生具有Q個粒子的粒子池，通過CFD計算獲取適應值，然后從中選取適應值最好的M個粒子作為父級MQPSO的初始種群，對種群中每個粒子按照MQPSO更新位置的方法繁殖X個個體(式(2)中 φ，u為(0，1)的隨機數，所以該X 個個體位于不同的進化位置)，此時用已計算的粒子構造Kriging模型預測繁殖的粒子，選出適應值最優的粒子)，計算該粒子與父級粒子種群中最近的粒子的“歐氏距離”，據此確定一個合適的閥值半徑 R，以)為中心R為半徑確定一個變量區域，在該區域中嵌入種群規模為m迭代代數為n的子級MQPSO進行優化尋找代理模型的最優解，將該解作為正式進化的個體進行CFD適應值計算，完成后立即更新pbest、gbest、mbest和代理模型，直到N代父級MQPSO優化結束。

閥值半徑R由下式確定:

其中，Dim為變量維數，α為控制半徑的系數，以各粒子的子級區域能覆蓋父級的進化區域，而相互之間又不產生過多的重疊為宜。

KMQPSO算法中具有Q個粒子的粒子池用于獲得初始代理模型樣本點，可以通過實驗設計方法或約束由“歐氏距離”表述的相似度在解空間產生初始樣本點。根據適應值排序原則選擇M個粒子的方法可獲得優良的初始種群而且這些粒子較均勻地分布在解空間。粒子的繁殖策略可以為子級MQPSO篩選優良的局部區域，嵌套的子級優化則有助于找到代理模型的局部最優。通過粒子的繁殖篩選和局部優化，將要進行CFD計算的粒子精準地定位于代理模型的局部最優，通過迭代更新不斷提高代理模型在最優解附近的預測精度，從而收斂速度更快，而又能保持多樣性保證優化精度。

2 機翼氣動外形優化

2.1 參數化方法

本文采用Hicks-Henne解析函數法［17］參數化機翼剖面翼型，新翼型由基本翼型的彎度或厚度加擾動構成。本文為了減少設計變量的個數，取彎度與厚度設計變量各4個，表達如下:

以平面形狀確定的大展弦比跨聲速機翼為優化對象，沿展向取如圖1的5個設計剖面(即沿展向η=2z/b為0、0.14、0.34、0.65 和 1.0 的機翼剖面，z為展向位置，b為展長)，每個剖面由8個設計變量控制剖面翼型形狀，翼根翼剖面扭轉角固定為2.8°，其他4個剖面的扭轉角為設計變量，并繞1/4弦線扭轉，這樣整個機翼一共由44個設計變量決定。

圖1 機翼設計剖面選取Fig.1 Design sections on wing planform

2.2 多目標多約束處理

評價適應值方法的選取將對優化結果產生重要影響。通常采用加權、分效能曲線加權、幾何求積等方法［18］評價適應值。本文采用一種多目標非線性加權組合方法，當前評估機翼的某項性能為cnew，根據設計經驗設定cbad為不可接受解，cgood為滿意解，使用如下公式確定機翼的適應值:

對于幾何約束問題通常采用罰函數加權方法，但是該方法會增加優化的難度。由于QPSO中粒子出現的位置不是確定的，而是由概率密度函數確定，因此可以設定循環直到產生符合幾何約束條件的機翼。而對于與性能相關的約束如低頭力矩的約束，本文將其作為設計目標進行優化。

2.3 Kriging模型性能測試

以設計機翼的基準翼型為設計對象，采用上述的8個設計變量的參數化方法生成400個不同的翼型，采用N-S雷諾平均方程計算繞翼型流場，湍流模型為k-ω兩方程模型。選取圖2中的樣本點個數，剩下的作為測試樣本，分別采用1階、2階回歸模型的原始Kriging和1階回歸模型的改進的Kriging進行預測。結果表明對于原始Kriging模型，采用2階回歸模型具有較高的精度，而改進的Kriging模型則表現出更高的預測精度。

圖2 Kriging模型預測翼型阻力系數的平均誤差Fig.2 Average error of CDbetween forecast and CFD

Kriging模型中不同階次的回歸模型對樣本點的個數有不同的要求，設變量個數為n，構造0階的回歸模型至少要1個樣本，1階至少要n+1個，而2階至少要0.5(n+1)(n+2)個，44個設計變量將至少需要1035個樣本點才能構造2階的Kriging模型，這在機翼優化中將耗費的計算量是難以接受的。以上述的44個設計變量的機翼為設計對象，生成400個不同的機翼，采用N-S雷諾平均方程計算繞機翼流場，湍流模型為k-ω兩方程模型。選取圖3中的樣本點個數，剩下的作為測試樣本，分別采用1階回歸模型的原始Kriging和1階的改進Kriging進行預測，改進的Kriging模型同樣表現出良好的預測精度。

2.4 機翼氣動外形優化算例

以上述跨聲速機翼為優化對象，巡航Ma=0.785，Re=2.2 ×107，CL=0.54，多目標優化阻力系數 CD、低頭力矩系數Cm以及對翼根的彎矩系數Cmx(繞機體軸，Cmx對機翼結構重量有重要影響)，加權系數δ分別為1.0、0.15、0.12，限定各剖面翼型最大相對厚度以及機翼容積不減。

圖3 Kriging模型預測機翼阻力系數的平均誤差Fig.3 Average error of CDbetween forecast and CFD

設置KMQPSO算法中Q為200，M為18，N為20，X隨迭代從5線性增加到20，m為20，n隨迭代從20線性增加到80，子級區域半徑系數α隨迭代從0.4線性增加到1.5，這些參數隨迭代線性變化是為了在優化初期保持種群多樣性，而在后期加快收斂速度。當總計算次數達400次時，適應值變化的幅度已經很小，監視設計變量值的變化發現已經表現出收斂的形態，圖4是某個設計變量值收斂的情況，前200次計算是對粒子池進行的，之后在設計變量區域內進行了有效的搜索至收斂。

圖5為機翼優化前后各設計剖面外形變化情況，剖面最大相對厚度和機翼容積都沒有降低，而阻力、低頭力矩和翼根彎矩都得到了改善。由圖6機翼剖面的的壓力分布可以看出，設計機翼得到無激波或極弱激波形態的壓力分布，說明設計結果是理想的。

圖4 某變量收斂情況Fig.4 Optimization convergence of one design parameter

圖5 各剖面初始形狀(虛線)與設計形狀(實線)Fig.5 Wing section profiles(initial-dash，design-solid)

圖6 剖面壓力系數Fig.6 Section pressure distributions

表1 機翼優化結果Table 1 Wing optimization results

整個優化過程在單臺PC機上完成，耗時約6天，如果采用大型計算機進行并行和分布式計算，將可能在數小時之內完成一次優化，使設計周期顯著降低。

3 結論

采用改進的量子粒子群MQPSO優化Kriging模型的相關模型參數可顯著提高代理模型預測精度，與具有雙層結構的量子粒子群算法KMQPSO相結合，實現對每個參與CFD計算的粒子進行精準更新。高維度的跨聲速機翼優化算例表明該方法穩定高效，設計結果理想，可大幅度提高隨機算法在三維氣動外形優化中的工程實用性。

［1］JAMESON A.Aerodynamic-structural design studies of lowsweep transonic wings［R］.AIAA 2008-145，2008.

［2］MORRIS A M，ALLEN C B，RENDALL T C S.Aerodynamic optimisation of modern transport wing using efficient variable fidelity shape parameterisation［R］.AIAA 2009-1277.

［3］左英桃，高正紅，夏露.基于Euler方程和離散共軛方法的氣動外形優化設計［J］.應用力學學報，2009，26(1):22-27.

［4］DUVIGNEAU R，VISONNEAU M.Hybrid genetic algorithms and neural networks for fast CFD-based design［R］.AIAA 2002-5465，2002.

［5］KHURANA M S，WINARTO H，SINHA A K.Airfoil optimisation by swarm algorithm with mutation and artificial neural networks［R］.AIAA 2009-1278，2009.

［6］任慶祝，宋文萍.基于Kriging模型的翼型多目標氣動優化設計研究［J］.航空計算技術，2009，39(3):77-82.

［7］SUN J，FENG B，XU W B.Particle swarm optimization with particles having quantum behavior［A］.Proc of the IEEE Congress on Evolutionary Computation［C］，2004:325-331.

［8］KENNEDY J，EBERHART R C.Particle swarm optimization［A］.Proceedings of IEEE International Conference on Neural Networks Perth［C］.Australia，1995:1942-1948.

［9］孔慶琴，孫俊，須文波.基于QPSO的改進算法［J］.計算機工程與應用，2007，43(28):58-60.

［10］CRESSIE N A C.Statistics for spatial data［M］.New York:John Wiley ＆ Sons，1993.

［11］LOPHAVEN S N，NIELSEN H B，SNDERGAARD J.DACE:a matlab kriging toolbox［R］.Technical Report MM-TR-2002-12.Denmark:Technical University of Denmark，2002.

［12］游海龍，賈新章.基于遺傳算法的Kriging模型構造與優化［J］.計算機輔助設計與圖形學學報，2007，19(1):64-68.

［13］陳鵬，李劍，管濤.基于PSO的Kriging相關模型參數優化［J］.微電子學與計算機，2009，26(4):178-181.

［14］于向軍，張利輝，李春然，等.克里金模型及其在全局優化設計中的應用［J］.中國工程機械學報，2006，4(3):259-261.

［15］LONG T，LIU L，WANG J B，et al.Multi-objective multidisciplinary optimization of long-endurance UAV wing using surrogate model in Model Center［R］.AIAA 2008-5918.

［16］RAJAGOPAL S，GANGULI R.Multidisciplinary design optimization of an UAV wing using kriging based multi-objective genetic algorithm［R］.AIAA 2009-2219，2009.

［17］HICKS R M，HENNE P A.Wing design by numerical optimization［J］.Journal of Aircraft，1978，15(7):407-413.

［18］《飛機設計手冊》總編委會.民用飛機總體設計［M］.北京:航空工業出版社，2005:58-59.