Kriging模型在機(jī)翼氣動(dòng)外形優(yōu)化中的應(yīng)用

孫美建，詹 浩

(西北工業(yè)大學(xué)航空學(xué)院，陜西西安 710072)

0 引言

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

基于代理模型的優(yōu)化方法已經(jīng)在翼型優(yōu)化中得到廣泛應(yīng)用［4－6］，用隨機(jī)優(yōu)化算法與代理模型相結(jié)合，可大幅度減少流場(chǎng)計(jì)算的次數(shù)。本文采用改進(jìn)的量子粒子群算法優(yōu)化Kriging模型的相關(guān)模型參數(shù)，提高Kriging模型預(yù)測(cè)精度，并結(jié)合具有雙層結(jié)構(gòu)的粒子群優(yōu)化算法，發(fā)展了同樣適用于三維氣動(dòng)外形優(yōu)化的高效算法。針對(duì)具有44個(gè)設(shè)計(jì)變量的跨聲速機(jī)翼采用N-S方程流場(chǎng)求解器進(jìn)行了優(yōu)化驗(yàn)證。

1 優(yōu)化算法

1.1 改進(jìn)的量子粒子群算法

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

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

1.2 Kriging代理模型

Kriging代理模型［10］起源于地理空間統(tǒng)計(jì)學(xué)，是一種估計(jì)方差最小的無(wú)偏估計(jì)模型，通過(guò)相關(guān)函數(shù)的作用，具有局部估計(jì)的特點(diǎn)，它可以較好地預(yù)估未知點(diǎn)處函數(shù)值的分布情況。Kriging模型中響應(yīng)值與設(shè)計(jì)變量之間的關(guān)系可以表示為:

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

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

z(x)是均值為0、方差為σ2的統(tǒng)計(jì)隨機(jī)過(guò)程，并且兩個(gè)插值點(diǎn)的協(xié)方差為:

其中R是點(diǎn)x(i)和x(j)的相關(guān)函數(shù)，本文選用高斯函數(shù)表示為:

未知點(diǎn)x0處的預(yù)測(cè)值(x0)和方差估計(jì)值通過(guò)如下形式給出:

通過(guò)極大似然估計(jì)確定相關(guān)模型參數(shù)θk:

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

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

代理模型方法可以根據(jù)已有的點(diǎn)預(yù)測(cè)未知點(diǎn)的適應(yīng)值分布情況，當(dāng)代理模型達(dá)到一定精度時(shí)，可代替真實(shí)求解用于工程優(yōu)化。傳統(tǒng)的代理模型優(yōu)化框架［14－16］如下:首先用一定數(shù)量的樣本點(diǎn)構(gòu)造代理模型，優(yōu)化找到代理模型的最優(yōu)解并與真實(shí)值對(duì)比，若不滿足收斂要求則加入新計(jì)算的真實(shí)值更新代理模型，直到優(yōu)化結(jié)束。該方法收斂速度快，但是嚴(yán)重依賴于代理模型預(yù)測(cè)的精度，對(duì)于復(fù)雜的高維度優(yōu)化問(wèn)題極容易陷入局部最優(yōu)。

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

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

閥值半徑R由下式確定:

其中，Dim為變量維數(shù)，α為控制半徑的系數(shù)，以各粒子的子級(jí)區(qū)域能覆蓋父級(jí)的進(jìn)化區(qū)域，而相互之間又不產(chǎn)生過(guò)多的重疊為宜。

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

2 機(jī)翼氣動(dòng)外形優(yōu)化

2.1 參數(shù)化方法

本文采用Hicks-Henne解析函數(shù)法［17］參數(shù)化機(jī)翼剖面翼型，新翼型由基本翼型的彎度或厚度加擾動(dòng)構(gòu)成。本文為了減少設(shè)計(jì)變量的個(gè)數(shù)，取彎度與厚度設(shè)計(jì)變量各4個(gè)，表達(dá)如下:

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

圖1 機(jī)翼設(shè)計(jì)剖面選取Fig.1 Design sections on wing planform

2.2 多目標(biāo)多約束處理

評(píng)價(jià)適應(yīng)值方法的選取將對(duì)優(yōu)化結(jié)果產(chǎn)生重要影響。通常采用加權(quán)、分效能曲線加權(quán)、幾何求積等方法［18］評(píng)價(jià)適應(yīng)值。本文采用一種多目標(biāo)非線性加權(quán)組合方法，當(dāng)前評(píng)估機(jī)翼的某項(xiàng)性能為cnew，根據(jù)設(shè)計(jì)經(jīng)驗(yàn)設(shè)定cbad為不可接受解，cgood為滿意解，使用如下公式確定機(jī)翼的適應(yīng)值:

對(duì)于幾何約束問(wèn)題通常采用罰函數(shù)加權(quán)方法，但是該方法會(huì)增加優(yōu)化的難度。由于QPSO中粒子出現(xiàn)的位置不是確定的，而是由概率密度函數(shù)確定，因此可以設(shè)定循環(huán)直到產(chǎn)生符合幾何約束條件的機(jī)翼。而對(duì)于與性能相關(guān)的約束如低頭力矩的約束，本文將其作為設(shè)計(jì)目標(biāo)進(jìn)行優(yōu)化。

2.3 Kriging模型性能測(cè)試

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

圖2 Kriging模型預(yù)測(cè)翼型阻力系數(shù)的平均誤差Fig.2 Average error of CDbetween forecast and CFD

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

2.4 機(jī)翼氣動(dòng)外形優(yōu)化算例

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

圖3 Kriging模型預(yù)測(cè)機(jī)翼阻力系數(shù)的平均誤差Fig.3 Average error of CDbetween forecast and CFD

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

圖5為機(jī)翼優(yōu)化前后各設(shè)計(jì)剖面外形變化情況，剖面最大相對(duì)厚度和機(jī)翼容積都沒有降低，而阻力、低頭力矩和翼根彎矩都得到了改善。由圖6機(jī)翼剖面的的壓力分布可以看出，設(shè)計(jì)機(jī)翼得到無(wú)激波或極弱激波形態(tài)的壓力分布，說(shuō)明設(shè)計(jì)結(jié)果是理想的。

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

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

圖6 剖面壓力系數(shù)Fig.6 Section pressure distributions

表1 機(jī)翼優(yōu)化結(jié)果Table 1 Wing optimization results

整個(gè)優(yōu)化過(guò)程在單臺(tái)PC機(jī)上完成，耗時(shí)約6天，如果采用大型計(jì)算機(jī)進(jìn)行并行和分布式計(jì)算，將可能在數(shù)小時(shí)之內(nèi)完成一次優(yōu)化，使設(shè)計(jì)周期顯著降低。

3 結(jié)論

采用改進(jìn)的量子粒子群MQPSO優(yōu)化Kriging模型的相關(guān)模型參數(shù)可顯著提高代理模型預(yù)測(cè)精度，與具有雙層結(jié)構(gòu)的量子粒子群算法KMQPSO相結(jié)合，實(shí)現(xiàn)對(duì)每個(gè)參與CFD計(jì)算的粒子進(jìn)行精準(zhǔn)更新。高維度的跨聲速機(jī)翼優(yōu)化算例表明該方法穩(wěn)定高效，設(shè)計(jì)結(jié)果理想，可大幅度提高隨機(jī)算法在三維氣動(dòng)外形優(yōu)化中的工程實(shí)用性。

［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方程和離散共軛方法的氣動(dòng)外形優(yōu)化設(shè)計(jì)［J］.應(yīng)用力學(xué)學(xué)報(bào)，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模型的翼型多目標(biāo)氣動(dòng)優(yōu)化設(shè)計(jì)研究［J］.航空計(jì)算技術(shù)，2009，39(3):77-82.

［7］SUN J，F(xiàn)ENG 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ìn)算法［J］.計(jì)算機(jī)工程與應(yīng)用，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模型構(gòu)造與優(yōu)化［J］.計(jì)算機(jī)輔助設(shè)計(jì)與圖形學(xué)學(xué)報(bào)，2007，19(1):64-68.

［13］陳鵬，李劍，管濤.基于PSO的Kriging相關(guān)模型參數(shù)優(yōu)化［J］.微電子學(xué)與計(jì)算機(jī)，2009，26(4):178-181.

［14］于向軍，張利輝，李春然，等.克里金模型及其在全局優(yōu)化設(shè)計(jì)中的應(yīng)用［J］.中國(guó)工程機(jī)械學(xué)報(bào)，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］《飛機(jī)設(shè)計(jì)手冊(cè)》總編委會(huì).民用飛機(jī)總體設(shè)計(jì)［M］.北京:航空工業(yè)出版社，2005:58-59.