王超，高正紅,*，張偉，夏露，黃江濤

1. 西北工業(yè)大學(xué) 航空學(xué)院，西安 710072 2. 中國(guó)空氣動(dòng)力研究與發(fā)展中心 計(jì)算空氣動(dòng)力研究所，綿陽(yáng) 621000

基于代理模型的優(yōu)化(Surrogate-Based Optimization，SBO)方法可以提高優(yōu)化設(shè)計(jì)效率，便于實(shí)現(xiàn)全局優(yōu)化搜索，因此在工程設(shè)計(jì)領(lǐng)域得到廣泛應(yīng)用[1-4]。對(duì)于確定的優(yōu)化問(wèn)題，代理模型利用設(shè)計(jì)空間中有限的樣本，通過(guò)一定的模型假設(shè)構(gòu)建設(shè)計(jì)變量與目標(biāo)函數(shù)的近似關(guān)系，并對(duì)設(shè)計(jì)空間中的未知點(diǎn)進(jìn)行預(yù)測(cè)。代理模型的預(yù)測(cè)精度是影響基于代理模型優(yōu)化效果的關(guān)鍵，工程設(shè)計(jì)中常采用具有描述強(qiáng)非線性問(wèn)題能力的代理模型，如徑向基函數(shù)(RBF)[5]、Kriging[6]、支持向量回歸(SVR)[7]等模型，這類代理模型的預(yù)測(cè)精度與設(shè)計(jì)空間內(nèi)樣本的密度密切相關(guān)[8]。對(duì)于工程優(yōu)化設(shè)計(jì)問(wèn)題，特別是在初步設(shè)計(jì)階段，通常需要選擇盡可能大的設(shè)計(jì)空間，以獲得具有全局最優(yōu)的設(shè)計(jì)結(jié)果。然而，設(shè)計(jì)空間的擴(kuò)大將直接導(dǎo)致樣本在空間中的分布密度急劇下降，嚴(yán)重影響代理模型的精度。增加樣本數(shù)目可以改善代理模型的精度，但是由此會(huì)帶來(lái)巨大的計(jì)算代價(jià)，從而大幅度降低代理模型的優(yōu)化效率。因此，如何構(gòu)建盡可能包含全局最優(yōu)設(shè)計(jì)點(diǎn)的有效設(shè)計(jì)空間，以及如何選擇用于重構(gòu)代理模型的有效樣本是基于代理模型優(yōu)化方法需要解決的關(guān)鍵問(wèn)題。

Jones等[9]于1998年提出了基于Kriging模型的高效全局優(yōu)化(Efficient Global Optimization，EGO)方法。EGO方法將優(yōu)化搜索與代理模型加點(diǎn)重構(gòu)相結(jié)合，根據(jù)Kriging模型對(duì)目標(biāo)函數(shù)預(yù)測(cè)的均值和方差，將目標(biāo)函數(shù)改善的期望(Expected Improvement，EI)作為Kriging模型優(yōu)化的適應(yīng)值函數(shù)。由于EGO方法同時(shí)考慮了代理模型的預(yù)測(cè)值以及預(yù)測(cè)的不確定性，因此具有全局優(yōu)化能力[10]。Jones[11]對(duì)Kriging模型不同的加點(diǎn)函數(shù)進(jìn)行對(duì)比研究時(shí)指出，EI加點(diǎn)準(zhǔn)則存在收斂緩慢，優(yōu)化效率低的問(wèn)題。為了進(jìn)一步提高EGO方法的全局性和高效性，Schonlau和Welch[12]將EI函數(shù)進(jìn)行了推廣，引入了一個(gè)調(diào)節(jié)算法全局和局部搜索能力的參數(shù)；Sasena[13]、韓忠華[14-15]等討論了不同的加點(diǎn)函數(shù)和約束處理方法；Soberster[16]、Ginsbourger[17-18]、Liu[19]等發(fā)展了多種并行加點(diǎn)策略。

盡管各類改進(jìn)的EGO方法提高了樣本的利用率，但是對(duì)于大設(shè)計(jì)空間內(nèi)的優(yōu)化問(wèn)題，EGO方法仍然難以實(shí)現(xiàn)全局性與高效性的統(tǒng)一。例如，對(duì)于氣動(dòng)外形優(yōu)化設(shè)計(jì)問(wèn)題，尤其是機(jī)翼/翼型等自由曲線曲面類的外形設(shè)計(jì)，如果增大設(shè)計(jì)空間，優(yōu)化過(guò)程中會(huì)生成許多流場(chǎng)分離甚至無(wú)法計(jì)算收斂的“非翼型”外形，給樣本帶來(lái)許多“噪聲”；另一方面，大的設(shè)計(jì)空間會(huì)導(dǎo)致代理模型的精度急劇下降[20]，代理優(yōu)化需要大量的加點(diǎn)和更新，從而導(dǎo)致無(wú)法在可接受的計(jì)算代價(jià)內(nèi)獲得理想的設(shè)計(jì)結(jié)果。因此，應(yīng)用EGO方法進(jìn)行氣動(dòng)優(yōu)化設(shè)計(jì)時(shí)，為了獲得全局最優(yōu)解，一般采用多輪優(yōu)化的策略，即以上一輪優(yōu)化的結(jié)果作為新的初始外形，重新劃定設(shè)計(jì)空間進(jìn)行優(yōu)化。例如，Zhang等[21]對(duì)NACA0012翼型進(jìn)行三輪優(yōu)化，達(dá)到了全局優(yōu)化的目的。但是，多輪優(yōu)化的方法會(huì)導(dǎo)致人工成本增加，優(yōu)化效率降低。

針對(duì)EGO方法在氣動(dòng)優(yōu)化設(shè)計(jì)中面臨的主要問(wèn)題，本文對(duì)基于Kriging代理模型優(yōu)化的加點(diǎn)方法和設(shè)計(jì)空間的構(gòu)建問(wèn)題進(jìn)行了研究。首先，針對(duì)EGO方法本身的收斂性問(wèn)題，本文采用先全局再局部的優(yōu)化思想，提出了一種混合加點(diǎn)方法，該方法通過(guò)引入EI閾值控制不同的加點(diǎn)準(zhǔn)則，以提高EGO方法在確定設(shè)計(jì)空間內(nèi)的收斂性。針對(duì)設(shè)計(jì)空間的構(gòu)建問(wèn)題，本文對(duì)比了擴(kuò)大設(shè)計(jì)變量范圍和多輪優(yōu)化兩種不同的設(shè)計(jì)空間構(gòu)建方法，分析了設(shè)計(jì)變量范圍對(duì)設(shè)計(jì)空間大小和樣本密度的影響，進(jìn)而提出了自適應(yīng)設(shè)計(jì)空間擴(kuò)展的高效代理模型優(yōu)化方法。自適應(yīng)設(shè)計(jì)空間擴(kuò)展的方法在動(dòng)態(tài)的設(shè)計(jì)空間中進(jìn)行優(yōu)化搜索，只在有潛力的維度擴(kuò)展設(shè)計(jì)空間邊界，通過(guò)構(gòu)建自適應(yīng)的設(shè)計(jì)空間高效地分配樣本。最后，通過(guò)AIAA氣動(dòng)優(yōu)化設(shè)計(jì)討論組(Aerodynamic Design Optimization Discussion Group，ADDG)[22]發(fā)布的標(biāo)準(zhǔn)翼型氣動(dòng)優(yōu)化設(shè)計(jì)算例對(duì)本文提出的自適應(yīng)設(shè)計(jì)空間擴(kuò)展的優(yōu)化方法進(jìn)行了驗(yàn)證。

1 基于Kriging模型的混合加點(diǎn)方法

1.1 Kriging模型

Kriging模型是一種基于統(tǒng)計(jì)理論的插值技術(shù)，對(duì)于一定數(shù)目的樣本集X=[x1x2…xn]T及其響應(yīng)值y=[y1y2…yn]T，Kriging模型假設(shè)目標(biāo)函數(shù)值與設(shè)計(jì)變量之間的真實(shí)關(guān)系可以表示為

y(x)=g(x)+z(x)

(1)

g(x)是回歸部分，是用來(lái)描述全局趨勢(shì)的函數(shù)，一般由多項(xiàng)式線性疊加表示，即

(2)

式中：fi(x)為低階多項(xiàng)式基函數(shù)，可為0階、一階或二階多項(xiàng)式；p為多項(xiàng)式基函數(shù)個(gè)數(shù)；βi為多項(xiàng)式回歸系數(shù)，一般通過(guò)廣義最小二乘法求得。

z(x)是一個(gè)靜態(tài)的隨機(jī)過(guò)程，用來(lái)描述真實(shí)樣本與趨勢(shì)函數(shù)之間的偏差，滿足：

(3)

(4)

(5)

Kriging模型在未知點(diǎn)xnew的預(yù)測(cè)值表示為已知訓(xùn)練樣本的線性組合，即

(6)

式中：c為線性加權(quán)系數(shù)。通過(guò)預(yù)測(cè)方差最小的無(wú)偏估計(jì)[23]方法，可以求得Kriging模型的預(yù)測(cè)值為

R-1(y-f(xnew)Tβ)

(7)

同時(shí)得到Kriging模型的預(yù)測(cè)方差為

σ2(xnew)=σz2[1+ζT(FTR-1F)-1ζ-rTR-1r]

(8)

式中：ζ=FTR-1r-f,F=[f(x1)f(x2)…f(xn)]T，R為樣本點(diǎn)之間的相關(guān)函數(shù)R(xi,xj)組成的相關(guān)矩陣，r為未知點(diǎn)xnew與樣本點(diǎn)之間的相關(guān)函數(shù)R(xnew,xi)組成的相關(guān)矩陣。

基于隨機(jī)過(guò)程的Kriging模型與其他代理模型最大的區(qū)別是它的輸出是一個(gè)隨機(jī)變量，可以同時(shí)給出預(yù)測(cè)的均值和方差。根據(jù)Kriging模型預(yù)測(cè)的均值和方差，通過(guò)一定的加點(diǎn)準(zhǔn)則[11]可以實(shí)現(xiàn)高效的優(yōu)化設(shè)計(jì)。

1.2 EGO方法

基于Kriging模型優(yōu)化最常用的是EGO方法[24-27]，又稱EI加點(diǎn)準(zhǔn)則?？紤]以下一般性優(yōu)化問(wèn)題：

(9)

(10)

假設(shè)當(dāng)前樣本中目標(biāo)函數(shù)的最小值為ymin，則Kriging模型在未知點(diǎn)處的預(yù)測(cè)值y相對(duì)于ymin的改善值為

I(x)=max(ymin-y,0)

(11)

根據(jù)Kriging模型對(duì)未知點(diǎn)預(yù)測(cè)的概率分布，可以計(jì)算任意未知點(diǎn)相對(duì)于當(dāng)前最優(yōu)值改善的期望。圖1給出了目標(biāo)函數(shù)改善的幾何意義，圖中紅色部分為目標(biāo)函數(shù)改善的概率，它反映了目標(biāo)函數(shù)改善可能性大小。對(duì)目標(biāo)改善值I(x)

求一階矩，可以得到目標(biāo)改善期望的函數(shù)為

(12)

式中：Φ為標(biāo)準(zhǔn)正態(tài)分布函數(shù)；φ為標(biāo)準(zhǔn)正態(tài)分布的密度函數(shù)。

綜上，對(duì)于式(9)中的原始優(yōu)化問(wèn)題，EGO方法的優(yōu)化模型為

(13)

EGO方法將EI函數(shù)作為Kriging模型優(yōu)化的適應(yīng)值函數(shù)，每次把搜索到EI函數(shù)的最大值點(diǎn)加入到樣本中，更新代理模型，直到滿足終止條件。由于EGO方法可以同時(shí)考慮Kriging模型的預(yù)測(cè)值及其不確定性，因此具有高效的全局優(yōu)化能力[9]。

用二維Rastrigin函數(shù)對(duì)EGO方法進(jìn)行測(cè)試， 其函數(shù)表達(dá)式為

f(x,y)=x2+y2-10cos(2πx)-

10cos(2πy)+20

(14)

圖2給出了二維Rastrigin函數(shù)圖，該函數(shù)在x,y∈[-1,1]的設(shè)計(jì)空間內(nèi)，存在多個(gè)局部最優(yōu)，其中全局最優(yōu)值為f(0,0)=0。

采用拉丁超立方取樣(Latin Hypercubic Sampling, LHS)方法[28]生成20個(gè)樣本，構(gòu)建初始Kriging模型，然后采用EGO方法對(duì)Kriging模型進(jìn)行優(yōu)化，優(yōu)化算法采用粒子群算法[29]，加點(diǎn)50次。圖3給出了目標(biāo)函數(shù)fmin的收斂歷程，圖4給出了每次搜索到EI函數(shù)最大值EImax的收斂歷程?？梢?jiàn)，在優(yōu)化后期，EI函數(shù)的最大值變得非常小，其量級(jí)已經(jīng)遠(yuǎn)小于目標(biāo)函數(shù)的最小值(如圖5所示)，此時(shí)EI函數(shù)已經(jīng)無(wú)法對(duì)代理優(yōu)化進(jìn)行有效地指引。圖6給出了初始樣本和加點(diǎn)樣本的位置，雖然EGO方法搜索到了最優(yōu)區(qū)域，但是由于EI函數(shù)逐漸失效的問(wèn)題，導(dǎo)致EGO方法在優(yōu)化的后期基本變?yōu)殡S機(jī)搜索，優(yōu)化收斂十分緩慢。

1.3 混合加點(diǎn)方法

針對(duì)EGO方法中EI函數(shù)在加點(diǎn)后期失效的問(wèn)題，本文提出混合(Hybrid)加點(diǎn)方法：當(dāng)EI函數(shù)的最大值小于設(shè)定的閾值CEI時(shí)，轉(zhuǎn)換為最小預(yù)測(cè)值(Minimum Prediction, MP)加點(diǎn)準(zhǔn)則。MP加點(diǎn)準(zhǔn)則直接將代理模型對(duì)目標(biāo)函數(shù)的預(yù)測(cè)值作為適應(yīng)值函數(shù)，具有很強(qiáng)的局部挖掘能力，收斂速度快[11]，因此可以彌補(bǔ)EI函數(shù)失效的問(wèn)題。本文構(gòu)建的Hybrid加點(diǎn)方法的適應(yīng)值函數(shù)表示為

Fitness=

(15)

式中：CEI=0.01ymin，即當(dāng)前樣本中目標(biāo)函數(shù)最小值改善的期望小于1%時(shí)，由EI加點(diǎn)準(zhǔn)則轉(zhuǎn)換為MP加點(diǎn)準(zhǔn)則。

使用Hybrid加點(diǎn)方法對(duì)二維Rastrigin函數(shù)重新進(jìn)行優(yōu)化，圖7給出了加點(diǎn)過(guò)程中最大EI值與最小目標(biāo)值的比值，可見(jiàn)當(dāng)加到第34個(gè)點(diǎn)時(shí)，EImax/ymin<0.01，此時(shí)EI準(zhǔn)則轉(zhuǎn)換為MP準(zhǔn)則。

圖8給出了Hybrid加點(diǎn)過(guò)程中不同適應(yīng)值函數(shù)搜索到的最優(yōu)值，其中左軸代表EI函數(shù)最大值，右軸代表預(yù)測(cè)目標(biāo)函數(shù)的最小值。

圖9給出了Hybrid加點(diǎn)方法與EGO方法中目標(biāo)函數(shù)的收斂歷程對(duì)比，可見(jiàn)Hybrid方法明顯優(yōu)于EGO方法，EI準(zhǔn)則向MP準(zhǔn)則的轉(zhuǎn)換大幅度加速了收斂進(jìn)程。圖10給出了Hybrid方法的加點(diǎn)位置，可見(jiàn)Hybrid方法首先采用EI準(zhǔn)則，對(duì)設(shè)計(jì)空間的最優(yōu)區(qū)域進(jìn)行全局探索，當(dāng)EI函數(shù)最大值小于CEI時(shí)，說(shuō)明Kriging模型在設(shè)計(jì)空間最優(yōu)區(qū)域內(nèi)已經(jīng)具有較高預(yù)測(cè)精度，此時(shí)轉(zhuǎn)換為MP準(zhǔn)則。由圖10可見(jiàn)，Hybrid方法中的MP樣本點(diǎn)均在全局最優(yōu)值附近，可以對(duì)代理模型進(jìn)行局部精細(xì)化挖掘，提高優(yōu)化的收斂速度。

2 基于混合加點(diǎn)方法的翼型優(yōu)化設(shè)計(jì)

分別采用EGO方法和Hybrid加點(diǎn)方法對(duì)RAE2822翼型進(jìn)行氣動(dòng)優(yōu)化設(shè)計(jì)。該算例為 ADODG[22]發(fā)布的第2個(gè)標(biāo)準(zhǔn)算例：跨聲速黏性條件下RAE2822翼型減阻優(yōu)化設(shè)計(jì)。翼型的設(shè)計(jì)狀態(tài)為：馬赫數(shù)Ma=0.734, 升力系數(shù)CL=0.824,雷諾數(shù)Re=6.5×105，優(yōu)化目標(biāo)為固定升力系數(shù)下阻力系數(shù)CD最小，約束條件為低頭力矩系數(shù)Cm≤0.92，同時(shí)翼型的面積A不減小。RAE2822翼型優(yōu)化的數(shù)學(xué)模型為

(16)

采用C型多塊結(jié)構(gòu)網(wǎng)格進(jìn)行CFD計(jì)算，RAE2822翼型網(wǎng)格的拓?fù)浣Y(jié)構(gòu)如圖11所示，c為翼型弦長(zhǎng)，網(wǎng)格的具體參數(shù)如表1所示。

采用CST方法[30]對(duì)翼型進(jìn)行參數(shù)化，翼型上下表面各12個(gè)設(shè)計(jì)變量，設(shè)計(jì)變量數(shù)目為24，給定每個(gè)設(shè)計(jì)變量的擾動(dòng)范圍，翼型的設(shè)計(jì)空間如圖12所示。

使用LHS方法在設(shè)計(jì)空間生成100個(gè)初始樣本，分別采用EGO和Hybrid方法添加200個(gè)樣本點(diǎn)，在Hybrid方法中，設(shè)定CEI=0.01CD_min，CD_min為樣本中的最小阻力系數(shù)。采用粒子群算法進(jìn)行代理模型優(yōu)化搜索，圖13給出了不同優(yōu)化

Table1ParametersofcomputationalgridofRAE2822airfoil

參數(shù)數(shù)值遠(yuǎn)場(chǎng)大小50網(wǎng)格規(guī)模601×213物面第一層網(wǎng)格距離5×10-6物面法向網(wǎng)格增長(zhǎng)率1.13翼型前緣網(wǎng)格距離0.001翼型后緣網(wǎng)格距離0.001

方法中翼型阻力系數(shù)的收斂歷程，表2給出了不同方法優(yōu)化的翼型與初始翼型的計(jì)算結(jié)果對(duì)比，可見(jiàn)Hybrid加點(diǎn)方法的設(shè)計(jì)結(jié)果優(yōu)于EGO方法，在滿足約束的條件下，Hybrid方法優(yōu)化翼型的阻力系數(shù)比RAE2822翼型降低了91.8 counts。

圖14給出了使用Hybrid方法優(yōu)化翼型與RAE2822翼型的外形對(duì)比，優(yōu)化翼型頭部變尖，最大厚度位置后移，后緣彎度增大。圖15給出了優(yōu)化翼型與RAE2822翼型的壓力系數(shù)Cp分布對(duì)比，圖16和圖17分別給出了優(yōu)化前后翼型流場(chǎng)的壓力云圖?？梢?jiàn)，優(yōu)化翼型的前緣吸力峰增大，后加載增加，壓力恢復(fù)平緩，激波得到消除。

Table2OptimizationresultsofEGOmethodandHybridinfillmethod

方法CLCDCmARAE28220.8240.020 31-0.092 70.077 9EGO0.8240.011 25-0.091 90.077 9混合加點(diǎn)方法0.8240.011 13-0.091 90.077 9

3 設(shè)計(jì)空間的擴(kuò)展

3.1 代理模型加點(diǎn)分析

為了進(jìn)一步探討設(shè)計(jì)空間對(duì)優(yōu)化結(jié)果的影響，首先對(duì)RAE2822翼型優(yōu)化算例中的設(shè)計(jì)變量進(jìn)行歸一化處理。圖18通過(guò)平行坐標(biāo)系顯示了Hybrid方法中設(shè)計(jì)空間各維度的邊界以及在優(yōu)化過(guò)程中所有的加點(diǎn)樣本，圖19給出了所有加點(diǎn)樣本對(duì)應(yīng)的翼型外形。統(tǒng)計(jì)發(fā)現(xiàn)，加點(diǎn)樣本中有98.5%的樣本的某一維(或多維)落在設(shè)計(jì)空間的邊界，其中最優(yōu)樣本也位于設(shè)計(jì)空間的邊界處。由此可見(jiàn)，第2節(jié)RAE2822翼型優(yōu)化算例中的設(shè)計(jì)空間不足以尋找全局最優(yōu)，因此需要對(duì)翼型的設(shè)計(jì)空間進(jìn)行擴(kuò)展。

3.2 擴(kuò)展設(shè)計(jì)空間

設(shè)計(jì)空間的構(gòu)建和樣本的選取是基于代理模型優(yōu)化的關(guān)鍵，為了分析不同的設(shè)計(jì)空間擴(kuò)展方法和樣本數(shù)目對(duì)設(shè)計(jì)結(jié)果的影響，本文對(duì)比以下兩種設(shè)計(jì)空間擴(kuò)展方法：

方法1直接擴(kuò)展空間。將原始設(shè)計(jì)空間每一維的變化范圍擴(kuò)大一倍，擴(kuò)展的翼型設(shè)計(jì)空間如圖20所示。該方法中，初始樣本數(shù)目為100，加點(diǎn)500次，樣本總數(shù)為原始空間優(yōu)化的2倍。

方法2兩輪優(yōu)化。以原始設(shè)計(jì)空間優(yōu)化的翼型為初始，進(jìn)行第2輪優(yōu)化，每一維設(shè)計(jì)變量的變化范圍與第1輪優(yōu)化保持一致，第2輪優(yōu)化翼型的設(shè)計(jì)空間如圖21所示。該方法中，初始樣本數(shù)目為100，加點(diǎn)200次，第2輪優(yōu)化使用的樣本總數(shù)與原始空間優(yōu)化相同。

使用Opt1代表方法1的優(yōu)化結(jié)果，Opt2代表方法2的優(yōu)化結(jié)果，兩種方法優(yōu)化阻力系數(shù)的收斂對(duì)比如圖22所示，優(yōu)化的翼型外形以及壓力分布對(duì)比如圖23和圖24所示?？梢?jiàn)，直接擴(kuò)展空間優(yōu)化的翼型與原始空間優(yōu)化的翼型差異很大，翼型上表面仍然存在明顯的激波。表3給出了兩種設(shè)計(jì)空間擴(kuò)展方法優(yōu)化翼型的計(jì)算結(jié)果對(duì)比，可以發(fā)現(xiàn)，直接擴(kuò)展空間的方法不如原始空間的設(shè)計(jì)結(jié)果，而第2輪優(yōu)化的結(jié)果優(yōu)于原始空間的設(shè)計(jì)結(jié)果。圖25給出了第2輪優(yōu)化中的加點(diǎn)樣本與設(shè)計(jì)空間邊界，可見(jiàn)最優(yōu)樣本落在設(shè)計(jì)空間內(nèi)部， 即第2輪優(yōu)化的設(shè)計(jì)空間包含了最優(yōu)點(diǎn)。

該算例表明，設(shè)計(jì)空間的構(gòu)建方法對(duì)代理模型優(yōu)化設(shè)計(jì)結(jié)果有重要的影響，在相同的計(jì)算成本下，直接擴(kuò)大設(shè)計(jì)空間的方法不如多輪優(yōu)化的方法。

方法CFD調(diào)用次數(shù)CDCmA初始翼型0.020 31-0.092 70.077 9原始空間優(yōu)化3000.011 13-0.091 90.077 9第2輪優(yōu)化3000.011 06-0.091 80.077 9直接擴(kuò)展空間6000.011 54-0.091 10.077 9

為了進(jìn)一步分析設(shè)計(jì)變量范圍對(duì)設(shè)計(jì)空間的影響，假設(shè)每一維設(shè)計(jì)變量的變化范圍均為R，維數(shù)為D，則設(shè)計(jì)空間的大小可以用空間的超體積表示為

V=RD

(17)

設(shè)計(jì)空間超體積隨著變量范圍的變化率為

(18)

可見(jiàn)，設(shè)計(jì)空間的大小隨變量范圍呈D-1次多項(xiàng)式增長(zhǎng)。對(duì)于RAE2822翼型優(yōu)化算例，如果每一維設(shè)計(jì)變量的范圍增加一倍，則設(shè)計(jì)空間擴(kuò)大為原來(lái)的224倍。

Kriging模型的精度與設(shè)計(jì)空間樣本的密度密切相關(guān)[20]，直接擴(kuò)大設(shè)計(jì)變量范圍會(huì)導(dǎo)致設(shè)計(jì)空間的暴增，從而帶來(lái)樣本嚴(yán)重的稀疏性，代理模型的精度大幅下降，如圖26所示。因此，如果直接擴(kuò)展設(shè)計(jì)空間，代理優(yōu)化難以收斂，無(wú)法在可接受的計(jì)算代價(jià)內(nèi)獲得理想的設(shè)計(jì)結(jié)果。

在兩輪優(yōu)化的方法中，以第1輪優(yōu)化的結(jié)果為初始外形，劃定與第1輪優(yōu)化同樣的變量范圍，則搜索空間只增加了一倍，所以第2輪優(yōu)化采用同樣的樣本得到了理想的設(shè)計(jì)結(jié)果。因此，在實(shí)際的翼型設(shè)計(jì)中，切實(shí)可行的方案是進(jìn)行多輪優(yōu)化而不是一次性給定足夠大的變量范圍。

4 自適應(yīng)設(shè)計(jì)空間擴(kuò)展方法

在基于代理模型的氣動(dòng)優(yōu)化設(shè)計(jì)方法中，設(shè)計(jì)者一般難以直接構(gòu)建合理的設(shè)計(jì)空間，如果設(shè)計(jì)空間太大，將導(dǎo)致有限的樣本難以滿足代理模型的精度要求，從而無(wú)法獲得理想的設(shè)計(jì)結(jié)果；如果設(shè)計(jì)空間太小，則需要多輪優(yōu)化，但是多輪優(yōu)化不僅增加了人工成本，而且降低了設(shè)計(jì)效率。針對(duì)代理模型氣動(dòng)優(yōu)化中設(shè)計(jì)空間的構(gòu)建問(wèn)題，本文提出自適應(yīng)設(shè)計(jì)空間擴(kuò)展的優(yōu)化方法，該方法不再固定設(shè)計(jì)空間的大小，設(shè)計(jì)空間在優(yōu)化過(guò)程中進(jìn)行動(dòng)態(tài)變化，只在有潛力的維度自適應(yīng)擴(kuò)展設(shè)計(jì)變量范圍，通過(guò)構(gòu)建自適應(yīng)的設(shè)計(jì)空間，提高樣本的利用效率。

在RAE2822翼型優(yōu)化算例中，加點(diǎn)樣本大部分都落到設(shè)計(jì)空間的邊界上，主要有兩方面的原因：一方面是因?yàn)樵O(shè)計(jì)空間太小，最優(yōu)解在設(shè)計(jì)空間之外，導(dǎo)致優(yōu)化搜索到了邊界；另一方面是因?yàn)榇砟Ｐ驮谶吔缣幍木忍?，?dǎo)致優(yōu)化搜索到了錯(cuò)誤的最優(yōu)。基于上述兩點(diǎn)原因，在進(jìn)行設(shè)計(jì)空間擴(kuò)展時(shí)，需要同時(shí)滿足邊界要求和精度要求。本文構(gòu)建的自適應(yīng)設(shè)計(jì)空間擴(kuò)展的代理模型氣動(dòng)優(yōu)化框架如圖27所示，具體步驟為

Step1設(shè)定初始的設(shè)計(jì)空間，即設(shè)計(jì)變量的個(gè)數(shù)以及每個(gè)設(shè)計(jì)變量的變化范圍。

Step2利用LHS方法在設(shè)計(jì)空間選取初始樣本，將初始樣本添加到樣本庫(kù)。

Step3根據(jù)樣本庫(kù)確定設(shè)計(jì)空間的邊界Xmin和Xmax，構(gòu)建當(dāng)前的設(shè)計(jì)空間。其中Xmin為樣本庫(kù)中每一維的最小值組成的向量，Xmax為樣本庫(kù)中每一維的最大值組成的向量。

Step4利用樣本庫(kù)中的樣本構(gòu)建代理模型。為了提高優(yōu)化效率，每添加10個(gè)樣本，重新訓(xùn)練Kriging模型的超參數(shù)。

Step5根據(jù)代理模型的加點(diǎn)準(zhǔn)則，優(yōu)化搜索得到Xbest。

Step6CFD評(píng)估Xbest，并將Xbest加入到樣本庫(kù)。

Step7判斷Xbest是否滿足設(shè)計(jì)空間擴(kuò)展條件。如果代理模型在Xbest處的預(yù)測(cè)值滿足精度要求，并且Xbest的某些維到達(dá)當(dāng)前設(shè)計(jì)空間的邊界，則生成自適應(yīng)樣本點(diǎn)Xadaptive，CFD評(píng)估Xadaptive，并將Xadaptive添加到樣本庫(kù)。

Xadaptive由Xbest通過(guò)在某些維度擴(kuò)展生成，其作用是擴(kuò)大設(shè)計(jì)空間的邊界。Xadaptive的計(jì)算方式為

(19)

式中：coefficient為設(shè)計(jì)空間擴(kuò)展系數(shù)，本文取0.1。

Step8判斷樣本總數(shù)N是否達(dá)到設(shè)定的最大樣本數(shù)目Nmax，如果N>Nmax，則終止程序，否則轉(zhuǎn)到Step 3。

采用本文構(gòu)建的自適應(yīng)設(shè)計(jì)空間擴(kuò)展的優(yōu)化方法對(duì)RAE2822翼型重新進(jìn)行優(yōu)化設(shè)計(jì)，初始樣本數(shù)為100，加點(diǎn)200次。以初始設(shè)計(jì)空間為參考對(duì)樣本每一維設(shè)計(jì)變量進(jìn)行歸一化處理，圖28給出了初始設(shè)計(jì)空間與最終設(shè)計(jì)空間中設(shè)計(jì)變量的邊界對(duì)比，對(duì)應(yīng)翼型外形邊界如圖29所示。由圖可見(jiàn)，初始設(shè)計(jì)空間在某些維度進(jìn)行了自適應(yīng)地?cái)U(kuò)展。表4給出了不同方法的優(yōu)化結(jié)果，自適應(yīng)設(shè)計(jì)空間的優(yōu)化結(jié)果與第2輪優(yōu)化結(jié)果基本一致，該方法通過(guò)一次優(yōu)化便達(dá)到了多輪優(yōu)化的效果，因此大幅提高了氣動(dòng)優(yōu)化的效率。

方法CFD調(diào)用次數(shù)CDCmA初始翼型0.020 31-0.092 70.077 9原始空間優(yōu)化3000.011 13-0.091 90.077 9第2輪優(yōu)化3000.011 06-0.091 80.077 9直接擴(kuò)展空間6000.011 54-0.091 10.077 9自適應(yīng)擴(kuò)展空間3000.011 05-0.091 90.077 9

圖30和圖31分別給出了使用不同方法優(yōu)化的翼型外形以及壓力分布對(duì)比，可見(jiàn)優(yōu)化的翼型外形以及壓力分布十分相似，說(shuō)明該算例的“最優(yōu)”翼型集中在初始翼型附近的一個(gè)區(qū)域。

為了對(duì)RAE2822翼型優(yōu)化結(jié)果的可靠性進(jìn)行驗(yàn)證，將本文的優(yōu)化結(jié)果與文獻(xiàn)[21]中的結(jié)果進(jìn)行對(duì)比，如表5所示?？梢园l(fā)現(xiàn)，由于網(wǎng)格和CFD求解器的不同，不同方法對(duì)初始翼型和優(yōu)化翼型計(jì)算的氣動(dòng)力系數(shù)存在一定的差異，因此難以直接衡量?jī)?yōu)化結(jié)果的優(yōu)劣。從翼型的外形與壓力分布來(lái)看，文獻(xiàn)[21]中優(yōu)化翼型的外形以及壓力分布的形態(tài)與本文的結(jié)果十分相似，二者均得到一個(gè)典型的跨聲速無(wú)激波壓力分布，如圖32和圖33所示。

與文獻(xiàn)[21]中基于代理模型優(yōu)化方法的計(jì)算代價(jià)進(jìn)行對(duì)比，本文的CFD計(jì)算次數(shù)為300，文獻(xiàn)[21]采用多輪優(yōu)化的方式，CFD計(jì)算次數(shù)為592。由此可以看出，本文提出的自適應(yīng)設(shè)計(jì)空間擴(kuò)展方法的優(yōu)化效率得到大幅提升。

表5 與文獻(xiàn)[21]中RAE2822翼型的優(yōu)化結(jié)果對(duì)比Table 5 Comparison with optimization results of RAE2822 airfoil in Ref. [21]

5 NACA0012翼型優(yōu)化設(shè)計(jì)

為了進(jìn)一步驗(yàn)證自適應(yīng)設(shè)計(jì)空間擴(kuò)展優(yōu)化方法的適用性，本節(jié)對(duì)NACA0012翼型進(jìn)行優(yōu)化設(shè)計(jì)。該算例為ADODG[22]發(fā)布的第1個(gè)標(biāo)準(zhǔn)算例：跨聲速無(wú)黏條件下NACA0012翼型的減阻設(shè)計(jì)。從文獻(xiàn)[31-36]的結(jié)果可以看出，該算例的優(yōu)化結(jié)果與初始翼型相比，變形尺度非常大，可以更好地考驗(yàn)自適應(yīng)設(shè)計(jì)空間擴(kuò)展方法的能力。

該算例的設(shè)計(jì)狀態(tài)為：Ma=0.85，迎角α=0°，優(yōu)化目標(biāo)為0°迎角下翼型阻力系數(shù)最小，約束條件為翼型每個(gè)位置處的厚度t不減小。NACA0012翼型優(yōu)化的數(shù)學(xué)模型為

(20)

由于該算例為0°迎角下對(duì)稱翼型的優(yōu)化，本文采用半模外形進(jìn)行計(jì)算和優(yōu)化，全模翼型的阻力為半模外形的2倍。采用O型多塊結(jié)構(gòu)網(wǎng)格進(jìn)行CFD計(jì)算，初始翼型網(wǎng)格的拓?fù)浣Y(jié)構(gòu)如圖34所示，網(wǎng)格參數(shù)的設(shè)置如表6所示。

采用CST方法對(duì)翼型進(jìn)行參數(shù)化，設(shè)計(jì)變量數(shù)目為24，給定每個(gè)設(shè)計(jì)變量的擾動(dòng)范圍，翼型的初始設(shè)計(jì)空間如圖35所示。

采用構(gòu)建的自適應(yīng)設(shè)計(jì)空間擴(kuò)展的優(yōu)化方法對(duì)NACA0012翼型進(jìn)行優(yōu)化，并與固定設(shè)計(jì)空間的優(yōu)化方法進(jìn)行對(duì)比。初始樣本數(shù)目為100，加點(diǎn)200次。圖36給出了歸一化后初始設(shè)計(jì)空間與最終設(shè)計(jì)空間的邊界對(duì)比，對(duì)應(yīng)的翼型外形邊界如圖37所示，可見(jiàn)自適應(yīng)設(shè)計(jì)空間的方法在翼型的后緣大幅度地?cái)U(kuò)展了設(shè)計(jì)變量的范圍。

Table6ParametersofNACA0012airfoilcomputationalgrid

參數(shù)數(shù)值遠(yuǎn)場(chǎng)大小50網(wǎng)格規(guī)模201×201物面第1層網(wǎng)格距離5×10-4物面法向網(wǎng)格增長(zhǎng)率1.05翼型前緣網(wǎng)格距離0.001翼型后緣網(wǎng)格距離0.001

圖38給出了自適應(yīng)設(shè)計(jì)空間與固定設(shè)計(jì)空間方法優(yōu)化阻力系數(shù)的收斂歷程，在相同的計(jì)算花費(fèi)下，自適應(yīng)設(shè)計(jì)空間的方法獲得了更優(yōu)的結(jié)果。表7給出了不同設(shè)計(jì)空間下翼型優(yōu)化結(jié)果對(duì)比，自適應(yīng)設(shè)計(jì)空間優(yōu)化的翼型阻力系數(shù)為55 counts，明顯小于固定設(shè)計(jì)空間的優(yōu)化結(jié)果。

圖39和圖40分別給出了固定設(shè)計(jì)空間與自適應(yīng)設(shè)計(jì)空間優(yōu)化的翼型外形以及壓力分布對(duì)比?？梢?jiàn)，自適應(yīng)設(shè)計(jì)空間優(yōu)化的翼型前緣和后緣更加飽滿， 翼型前緣半徑增加，后緣變形尺度增大；從壓力分布來(lái)看， 自適應(yīng)設(shè)計(jì)空間優(yōu)化的翼型前緣的吸力峰更高，翼型后緣通過(guò)連續(xù)地膨脹和壓縮使壓力得到緩和恢復(fù)，從而減弱了激波強(qiáng)度，降低了激波阻力。

方法CLCDNACA001200.047 1固定設(shè)計(jì)空間00.017 5自適應(yīng)設(shè)計(jì)空間00.005 5

從總壓損失的角度進(jìn)一步分析優(yōu)化翼型阻力系數(shù)減小的原因。在可壓縮條件下，總壓與靜壓的關(guān)系以及總壓系數(shù)表示為

(21)

(22)

式中：p0為總壓;p為靜壓;Ma為當(dāng)?shù)伛R赫數(shù);γ為氣體比熱比;p∞為遠(yuǎn)場(chǎng)靜壓;ρ∞為遠(yuǎn)場(chǎng)密度;v∞為自由來(lái)流速度。

NACA0012翼型以及兩種不同方法優(yōu)化翼型的流場(chǎng)總壓系數(shù)云圖如圖41～圖43所示，截取翼型表面的總壓系數(shù)如圖44所示。由圖可見(jiàn)，NACA0012翼型激波最強(qiáng)，在激波位置處總壓下降最明顯。圖45給出了翼型后緣一倍弦長(zhǎng)處(x/c=2)截面總壓系數(shù)的對(duì)比，可以發(fā)現(xiàn)，該截面處自適應(yīng)設(shè)計(jì)空間優(yōu)化翼型的總壓相對(duì)來(lái)流總壓損失最小，即能量損失最小，因此該翼型的阻力系數(shù)最小。

將本文的優(yōu)化結(jié)果與文獻(xiàn)[21]中結(jié)果進(jìn)行對(duì)比，圖46和圖47分別給出了文獻(xiàn)[21]通過(guò)多輪優(yōu)化方法得到的翼型外形和壓力分布，該結(jié)果與本文的優(yōu)化結(jié)果具有很強(qiáng)的一致性。表8給出了兩種不同方法優(yōu)化的結(jié)果對(duì)比，可以發(fā)現(xiàn)，本文的方法通過(guò)更少的CFD計(jì)算代價(jià)獲得了阻力更小的翼型。

該算例進(jìn)一步證實(shí)，自適應(yīng)設(shè)計(jì)空間擴(kuò)展優(yōu)化方法對(duì)于大尺度變形的氣動(dòng)優(yōu)化問(wèn)題具有很好的適應(yīng)性。

表8 與文獻(xiàn)[21]中NACA0012翼型的優(yōu)化結(jié)果對(duì)比

6 結(jié) 論

本文針對(duì)基于Kriging模型氣動(dòng)優(yōu)化中的加點(diǎn)方法以及設(shè)計(jì)空間的構(gòu)建問(wèn)題進(jìn)行了研究，提出了混合加點(diǎn)方法和自適應(yīng)設(shè)計(jì)空間擴(kuò)展的優(yōu)化方法，通過(guò)ADODG[22]標(biāo)準(zhǔn)翼型優(yōu)化算例對(duì)文中提出的方法進(jìn)行了驗(yàn)證。

1) 文中構(gòu)建的混合加點(diǎn)方法通過(guò)引入EI閾值控制EI和MP加點(diǎn)準(zhǔn)則，利用先全局再局部的優(yōu)化思想，可以提高EGO算法在確定設(shè)計(jì)空間內(nèi)的收斂性。

2) 設(shè)計(jì)空間的構(gòu)建方法對(duì)基于代理模型氣動(dòng)優(yōu)化的結(jié)果有重要的影響，如果設(shè)計(jì)空間的尺度偏于保守，則搜索空間受到限制，從而失去全局優(yōu)化能力；如果設(shè)計(jì)空間的尺度偏于激進(jìn)，將導(dǎo)致有限的樣本難以滿足代理模型的精度要求，從而無(wú)法獲得理想的設(shè)計(jì)結(jié)果。

3) 本文構(gòu)建的自適應(yīng)設(shè)計(jì)空間擴(kuò)展的方法打破了傳統(tǒng)代理模型優(yōu)化依賴人工經(jīng)驗(yàn)給定設(shè)計(jì)空間的局限，通過(guò)構(gòu)建自適應(yīng)的設(shè)計(jì)空間，實(shí)現(xiàn)了樣本的高效配置，很好地兼顧了氣動(dòng)優(yōu)化的高效性和全局性。通過(guò)ADODG標(biāo)準(zhǔn)翼型氣動(dòng)優(yōu)化設(shè)計(jì)算例證實(shí)，自適應(yīng)設(shè)計(jì)空間擴(kuò)展方法可以大幅提高氣動(dòng)優(yōu)化效率。