韓忠華

西北工業(yè)大學(xué) 航空學(xué)院 翼型葉柵空氣動(dòng)力學(xué)國(guó)家級(jí)重點(diǎn)實(shí)驗(yàn)室， 西安 710072

Kriging模型及代理優(yōu)化算法研究進(jìn)展

韓忠華*

西北工業(yè)大學(xué) 航空學(xué)院 翼型葉柵空氣動(dòng)力學(xué)國(guó)家級(jí)重點(diǎn)實(shí)驗(yàn)室， 西安 710072

代理模型方法由于能顯著提高工程優(yōu)化設(shè)計(jì)問(wèn)題的效率，在航空航天及其他領(lǐng)域得到了廣泛重視，并逐漸發(fā)展成為一類(lèi)優(yōu)化算法，本文稱(chēng)其為代理優(yōu)化(SBO)算法。在現(xiàn)有的代理模型方法中，如多項(xiàng)式響應(yīng)面、徑向基函數(shù)、神經(jīng)網(wǎng)絡(luò)、支持向量回歸、多變量插值/回歸、多項(xiàng)式混沌展開(kāi)等，源于地質(zhì)統(tǒng)計(jì)學(xué)的Kriging模型具有代表性，是一種非常具有應(yīng)用潛力的代理模型方法。以飛行器設(shè)計(jì)領(lǐng)域的優(yōu)化問(wèn)題為背景，介紹了Kriging代理模型及應(yīng)用于優(yōu)化設(shè)計(jì)的理論和算法的最新研究進(jìn)展。首先，概述了Kriging模型的基本理論和算法，并討論了影響Kriging模型魯棒性和效率的幾個(gè)關(guān)鍵性問(wèn)題。其次，回顧了Kriging模型理論和算法研究的3個(gè)最新研究進(jìn)展，包括梯度增強(qiáng)型Kriging、CoKriging和分層Kriging模型。而后，分析提煉了基于Kriging模型的代理優(yōu)化算法的優(yōu)化機(jī)制和優(yōu)化框架，給出了“優(yōu)化加點(diǎn)準(zhǔn)則”和“子優(yōu)化”的概念，并介紹了目前常用的幾種優(yōu)化加點(diǎn)準(zhǔn)則及其相應(yīng)子優(yōu)化問(wèn)題的求解與約束處理；同時(shí)，還介紹了最新提出的局部EI加點(diǎn)準(zhǔn)則以及代理優(yōu)化的終止條件。最后，介紹了代理優(yōu)化在標(biāo)準(zhǔn)測(cè)試函數(shù)算例驗(yàn)證、飛行器氣動(dòng)與多學(xué)科優(yōu)化設(shè)計(jì)典型算例確認(rèn)方面的研究進(jìn)展，并對(duì)當(dāng)前存在的一些關(guān)鍵科學(xué)問(wèn)題以及未來(lái)研究方向進(jìn)行了討論。

優(yōu)化方法； Kriging； 代理模型； 飛行器設(shè)計(jì)； 多學(xué)科設(shè)計(jì)優(yōu)化(MDO)

現(xiàn)代優(yōu)化方法在飛行器設(shè)計(jì)中的應(yīng)用可追溯到20世紀(jì)50年代末60年代初，最初應(yīng)用于結(jié)構(gòu)優(yōu)化設(shè)計(jì)[1-2]。大約20世紀(jì)70年代開(kāi)始被引入到氣動(dòng)優(yōu)化設(shè)計(jì)[3-4]，而后隨著多學(xué)科設(shè)計(jì)優(yōu)化(MDO)在20世紀(jì)80年代的興起[5-6]，逐漸成為航空航天各學(xué)科領(lǐng)域的研究熱點(diǎn)[7-8]。對(duì)于MDO的研究，由于涉及多個(gè)學(xué)科的大量設(shè)計(jì)變量和狀態(tài)變量，且各學(xué)科之間可能存在復(fù)雜耦合關(guān)系，使得直接調(diào)用各學(xué)科高可信度分析模型(如計(jì)算流體力學(xué)(CFD)或計(jì)算固體力學(xué)(CSM))時(shí)，面臨優(yōu)化難度大、優(yōu)化時(shí)間過(guò)長(zhǎng)的問(wèn)題。于是，代理模型方法應(yīng)運(yùn)而生，并逐漸得到重視，成為MDO研究的重要分支和關(guān)鍵技術(shù)[9-12]。

所謂代理模型，通常是指在分析和優(yōu)化設(shè)計(jì)過(guò)程中可替代那些比較復(fù)雜和費(fèi)時(shí)的數(shù)值分析的近似數(shù)學(xué)模型，也稱(chēng)為響應(yīng)面模型、近似模型或元模型[13-15]。代理模型方法不僅可大大提高優(yōu)化設(shè)計(jì)效率，而且可降低優(yōu)化難度，并有利于濾除數(shù)值噪聲和實(shí)現(xiàn)并行優(yōu)化設(shè)計(jì)。代理模型方法的雛形是20世紀(jì)70年代應(yīng)用于結(jié)構(gòu)優(yōu)化設(shè)計(jì)的多項(xiàng)式響應(yīng)面[16-17]；自20世紀(jì)80年代MDO興起后，代理模型方法成為了MDO研究的重要分支[18]；大約在20世紀(jì)末，代理模型方法被引入到氣動(dòng)優(yōu)化設(shè)計(jì)領(lǐng)域[19-21]。在代理模型方法研究方面，目前已發(fā)展了包括多項(xiàng)式響應(yīng)面 (RSM)[16-17]，Kriging模型[22-23]、徑向基函數(shù)(RBFs)[24-27]、神經(jīng)網(wǎng)絡(luò)(NN)[28-29]、支持向量回歸 (SVR)[30-31]、多變量插值和回歸 (MIR)[32-33]、多項(xiàng)式混沌展開(kāi)(PCE)[34-35]等多種代理模型方法[31，36-38]。文獻(xiàn)[31，36]對(duì)代理模型方法及其應(yīng)用的發(fā)展進(jìn)行了回顧。在多維空間樣本點(diǎn)合理選取研究方面，除了采用經(jīng)典試驗(yàn)設(shè)計(jì)方法外(例如全因子、中心組合、D-優(yōu)化等)，目前還發(fā)展了適用于計(jì)算機(jī)數(shù)值模擬實(shí)驗(yàn)的拉丁超立方 (LHS)[39]、正交設(shè)計(jì)、均勻設(shè)計(jì) (UD)[40]等現(xiàn)代試驗(yàn)設(shè)計(jì)方法。

在代理模型方法發(fā)展之初，其作用是作為近似模型，用于替代比較復(fù)雜和費(fèi)時(shí)的數(shù)值分析模型進(jìn)行優(yōu)化設(shè)計(jì)。對(duì)于這種“替代算法”，優(yōu)化結(jié)果很大程度上依賴(lài)于代理模型近似精度，且無(wú)法保證能收斂到真實(shí)最優(yōu)解。近年來(lái)，隨著研究的不斷深入，代理模型的作用發(fā)生了轉(zhuǎn)變，不再是簡(jiǎn)單替代，而是構(gòu)成了一種基于歷史數(shù)據(jù)來(lái)驅(qū)動(dòng)樣本點(diǎn)加入，并逼近局部或全局最優(yōu)解的優(yōu)化機(jī)制。在這種新的優(yōu)化機(jī)制下，通過(guò)合理構(gòu)造加點(diǎn)策略，可保證所產(chǎn)生的樣本點(diǎn)序列精確收斂于優(yōu)化問(wèn)題的真實(shí)最優(yōu)解。同時(shí)，復(fù)雜多維問(wèn)題的代理模型不必在整個(gè)設(shè)計(jì)空間內(nèi)具有高的近似精度，而只需在重要的區(qū)域特別是最優(yōu)解附近具有高的近似精度。總之，由于代理模型在優(yōu)化設(shè)計(jì)問(wèn)題中作用的轉(zhuǎn)變(以及一系列特定算法的產(chǎn)生)，使得基于代理模型的優(yōu)化方法正逐漸發(fā)展成為可與傳統(tǒng)梯度優(yōu)化和無(wú)梯度優(yōu)化具有相同功能的一類(lèi)通用優(yōu)化方法[38]。為了方便起見(jiàn)，本文稱(chēng)其為代理優(yōu)化(Surrogate-based Optimization， SBO)算法。

在現(xiàn)有的代理模型方法之中，源于地質(zhì)統(tǒng)計(jì)學(xué)的Kriging模型是一種具有代表性的方法。Kriging模型(也稱(chēng)為高斯隨機(jī)過(guò)程模型或克里金模型)的思想由南非采礦工程師Krige于1951年在其碩士論文中提出[22](Krige的生平可參見(jiàn)文獻(xiàn)[31])。在1963—1971年期間，法國(guó)數(shù)學(xué)家Matheron進(jìn)行了完善和發(fā)展，形成了完整的數(shù)學(xué)理論和模型[41](稱(chēng)為區(qū)域變量理論[42])。出于對(duì)Krige教授開(kāi)創(chuàng)性工作的尊重，他在1963年首次將該理論稱(chēng)為Krigeage[41](后來(lái)人們稱(chēng)其為Kriging模型)。

Kriging模型理論形成后，最初應(yīng)用于地質(zhì)儲(chǔ)量估計(jì)，后來(lái)發(fā)展成為一類(lèi)非常流行的地質(zhì)統(tǒng)計(jì)學(xué)插值方法。1989年，Sacks 教授等[23]將Kriging理論進(jìn)一步推廣到確定性計(jì)算機(jī)實(shí)驗(yàn)的設(shè)計(jì)和分析領(lǐng)域，并給出了一種較實(shí)用的Kriging算法(或Kriging模型)。由于Sacks 教授等里程碑式的研究工作，使得 Kriging模型不僅在地質(zhì)、水文、氣象、環(huán)境科學(xué)等自然科學(xué)領(lǐng)域得到應(yīng)用，也在航空航天、汽車(chē)等工程科學(xué)領(lǐng)域得到研究、發(fā)展和應(yīng)用。Kriging模型不僅能給出對(duì)未知函數(shù)的預(yù)估值，還能給出預(yù)估值的誤差估計(jì)，這是Kriging模型區(qū)別于其他代理模型的顯著特點(diǎn)。Kriging模型由于其對(duì)非線(xiàn)性函數(shù)的良好近似能力和獨(dú)特的誤差估計(jì)功能，正受到越來(lái)越多研究者的關(guān)注，是目前最具代表性和最具應(yīng)用潛力的代理模型方法之一。

在飛行器設(shè)計(jì)領(lǐng)域，隨著各學(xué)科高可信度數(shù)值模擬技術(shù)(如CFD或CSM)的發(fā)展和現(xiàn)代優(yōu)化算法的應(yīng)用，Kriging模型已經(jīng)逐步引起了該領(lǐng)域科研人員和工程設(shè)計(jì)人員的重視[43-45]?；贙riging模型的代理優(yōu)化算法，被嘗試應(yīng)用于氣動(dòng)優(yōu)化設(shè)計(jì)[46-52]、結(jié)構(gòu)優(yōu)化設(shè)計(jì)[53-54]、導(dǎo)彈設(shè)計(jì)[55-56]、航天飛行器設(shè)計(jì)[57]、高速列車(chē)外形優(yōu)化設(shè)計(jì)[58]以及多學(xué)科優(yōu)化設(shè)計(jì)[59-60]，并很可能發(fā)展成為該領(lǐng)域最實(shí)用和最流行的優(yōu)化方法。但是，與經(jīng)典的梯度優(yōu)化算法[61]和進(jìn)化算法(如遺傳算法[62]) 等優(yōu)化算法相比，基于Kriging模型的代理優(yōu)化算法雖然潛力巨大，但由于發(fā)展時(shí)間較短，有些算法還不夠成熟。目前所發(fā)表的關(guān)于Kriging模型的文獻(xiàn)分散在眾多不同學(xué)科領(lǐng)域，采用不同行業(yè)的科學(xué)語(yǔ)言進(jìn)行表述，缺乏系統(tǒng)性和通用性。特別在該航空航天領(lǐng)域，目前還沒(méi)有專(zhuān)門(mén)介紹Kriging模型理論、算法以及最新研究進(jìn)展的綜述性文獻(xiàn)報(bào)道，可供研究人員參考的經(jīng)驗(yàn)也很少。

在航空航天領(lǐng)域，由于各學(xué)科數(shù)值模擬技術(shù)的日益成熟，基于高可信度數(shù)值模擬技術(shù)的優(yōu)化設(shè)計(jì)方法及其在飛行器設(shè)計(jì)中的應(yīng)用，將繼續(xù)成為未來(lái)10～20年內(nèi)的研究熱點(diǎn)之一。例如最近NASA發(fā)布的2030年CFD愿景研究報(bào)告中[63]，多學(xué)科分析與優(yōu)化被列為需要重點(diǎn)發(fā)展的6大關(guān)鍵領(lǐng)域之一。因此，研究更高效、更實(shí)用的代理模型及其相應(yīng)的優(yōu)化理論和算法，無(wú)疑將是一項(xiàng)重要內(nèi)容。

本文以飛行器設(shè)計(jì)領(lǐng)域的優(yōu)化問(wèn)題為背景，系統(tǒng)介紹Kriging模型及其優(yōu)化設(shè)計(jì)方法的最新進(jìn)展。第1節(jié)概述了Kriging模型基本理論和算法；第2節(jié)對(duì)Kriging模型理論和算法研究的最新進(jìn)展進(jìn)行了綜述，并討論了影響建立Kriging模型的魯棒性和效率的幾個(gè)關(guān)鍵性問(wèn)題；第3節(jié)介紹了將Kriging模型應(yīng)用于優(yōu)化設(shè)計(jì)問(wèn)題的代理優(yōu)化算法最新進(jìn)展，包括優(yōu)化加點(diǎn)準(zhǔn)則和最新約束處理方法；第4節(jié)給出了代理優(yōu)化算法驗(yàn)證與確認(rèn)的最新進(jìn)展；最后，對(duì)該領(lǐng)域存在的若干關(guān)鍵科學(xué)問(wèn)題及未來(lái)發(fā)展方向進(jìn)行了討論。

1 Kriging模型基本理論和算法概述

Kriging模型存在不同變種，根據(jù)其采用的全局趨勢(shì)模型不同，可分為“簡(jiǎn)單Kriging”、“普通Kriging”和“泛Kriging”。這里以最常用的普通Kriging為例，對(duì)Kriging模型的基本理論和算法進(jìn)行簡(jiǎn)要介紹。

1.1 代理模型問(wèn)題的基本描述

應(yīng)用于飛行器優(yōu)化設(shè)計(jì)的數(shù)值分析方法(如CFD或CSM)，其共同特點(diǎn)是采用迭代方法求解所滿(mǎn)足的控制方程，最終獲得所關(guān)心的性能數(shù)據(jù)，如升力、阻力和力矩或應(yīng)力、變形等。對(duì)優(yōu)化設(shè)計(jì)問(wèn)題而言，可將數(shù)值分析程序或軟件看成一個(gè)“輸入-輸出”系統(tǒng)，即作為輸出量的目標(biāo)函數(shù)和約束函數(shù)可看成是設(shè)計(jì)變量的函數(shù)。如果能建立目標(biāo)函數(shù)和約束函數(shù)對(duì)設(shè)計(jì)變量的代理模型，就能對(duì)最優(yōu)點(diǎn)進(jìn)行快速預(yù)測(cè)，從而大大提高優(yōu)化效率。

(1)

對(duì)這n個(gè)設(shè)計(jì)方案進(jìn)行數(shù)值分析(CFD或CSM)，從分析結(jié)果中得到相應(yīng)的n個(gè)函數(shù)響應(yīng)值：

yS=[y(1)y(2)…y(n)]T=

(2)

1.2 Kriging模型及其預(yù)估值

Kriging模型是一種插值模型，其插值結(jié)果定義為已知樣本函數(shù)響應(yīng)值的線(xiàn)性加權(quán)，即

(3)

于是，只要能給出加權(quán)系數(shù)ω=[w(1)w(2)…w(n)]T的表達(dá)式，便可得到設(shè)計(jì)空間中任意設(shè)計(jì)方案的性能預(yù)估值。為了計(jì)算加權(quán)系數(shù)，Kriging模型引入統(tǒng)計(jì)學(xué)假設(shè)：將未知函數(shù)看成是某個(gè)高斯靜態(tài)隨機(jī)過(guò)程的具體實(shí)現(xiàn)。該靜態(tài)隨機(jī)過(guò)程定義為

Y(x)=β0+Z(x)

(4)

式中：β0為未知常數(shù)，也稱(chēng)為全局趨勢(shì)模型，代表Y(x)的數(shù)學(xué)期望值；Z(·)為均值為零、方差為σ2(σ2(x)≡σ2,?x)的靜態(tài)隨機(jī)過(guò)程。在設(shè)計(jì)空間不同位置處，這些隨機(jī)變量存在一定的相關(guān)性(或協(xié)方差)。該協(xié)方差可表述為

Cov[Z(x),Z(x′)]=σ2R(x,x′)

(5)

式中：R(x,x′)為 “相關(guān)函數(shù)”(只與空間距離有關(guān))，并滿(mǎn)足距離為零時(shí)等于1；距離無(wú)窮大時(shí)等于0；相關(guān)性隨著距離的增大而減小。

基于上述假設(shè)，Kriging模型尋找最優(yōu)加權(quán)系數(shù)ω，使得均方差

(6)

最小，并滿(mǎn)足如下插值條件(或無(wú)偏差條件)：

(7)

采用拉格朗日乘數(shù)法，經(jīng)過(guò)推導(dǎo)可證明最優(yōu)加權(quán)系數(shù)ω由式(8)描述的線(xiàn)性方程組(也稱(chēng)為Kriging模型方程組)給出

(8)

式中：i=1,2，…，n；μ為拉格朗日乘數(shù)。式(8)寫(xiě)成矩陣形式為

(9)

式中：

(10)

其中：R為“相關(guān)矩陣”，由所有已知樣本點(diǎn)之間的相關(guān)函數(shù)值組成；r為“相關(guān)矢量”，由未知點(diǎn)與所有已知樣本點(diǎn)之間的相關(guān)函數(shù)值組成。求解線(xiàn)性方程組式(9)，并代入式(3)可得Kriging模型預(yù)估值為

(11)

通過(guò)分塊矩陣求逆，該模型可最終寫(xiě)為

(12)

式中：β0=(FTR-1F)-1FTR-1yS；列向量Vkrig只與已知樣本點(diǎn)有關(guān)，可在模型訓(xùn)練(見(jiàn)1.4節(jié))結(jié)束后一次性計(jì)算并存儲(chǔ)。之后，預(yù)測(cè)任意x處的函數(shù)值只需要計(jì)算r(x)與Vkrig之間的點(diǎn)乘，其計(jì)算時(shí)間相比CFD或CSM分析而言完全可以忽略。此外，Kriging模型還能夠給出預(yù)估值的均方差估計(jì)：

σ2{1.0-rTR-1r+(1-FTR-1r)2/FTR-1F}

(13)

該均方差可用于指導(dǎo)加入新樣本點(diǎn)，以提高代理模型精度或逼近優(yōu)化問(wèn)題的最優(yōu)解(見(jiàn)第4節(jié))。

1.3 相關(guān)函數(shù)的選擇

在式(12) 所示的Kriging模型中，R和r的構(gòu)造均涉及相關(guān)函數(shù)的選擇和計(jì)算。目前比較流行的一類(lèi)相關(guān)函數(shù)為

(14)

1≤pk≤2；k=1,2，…,m

(15)

圖1 一維高斯指數(shù)函數(shù)示意圖Fig.1 Schematic of 1D Gaussian exponential function

實(shí)踐證明，在實(shí)際應(yīng)用中一些其他的相關(guān)函數(shù)可近似滿(mǎn)足高斯假設(shè)，也具有很好的性能。例如，目前比較流行一類(lèi)三次樣條函數(shù)：

(16)

該函數(shù)二階可導(dǎo)，在光滑性和魯棒性方面取得很好的折中。

1.4 模型參數(shù)訓(xùn)練(或優(yōu)化)

建立Kriging模型時(shí)，可以對(duì)其模型參數(shù)(或超參數(shù))進(jìn)行訓(xùn)練，從而大大增強(qiáng)代理模型的靈活性。一般采用“最大似然估計(jì)”或“交叉驗(yàn)證”[66]。最大似然估計(jì)是使式(17)函數(shù)值最大：

(17)

式中：超參數(shù)p僅針對(duì)高斯指數(shù)函數(shù)才存在；|R|為相關(guān)矩陣的行列式。對(duì)于式(17),β0和σ2的最優(yōu)值可解析地給出：

β0(θ,p)=(FTR-1F)-1FTR-1yS

(18)

將式(18)代入式(17)后，取對(duì)數(shù)，并忽略常數(shù)項(xiàng)，優(yōu)化問(wèn)題轉(zhuǎn)化為使式(19)最大：

(19)

由于無(wú)法解析地求出θ和p的最優(yōu)值，這里需要采用數(shù)值優(yōu)化算法，如擬牛頓方法等梯度優(yōu)化算法[61]或Hooke & Jeeves模式搜索[67]、Simplex單純形法[68]、遺傳算法[62]等無(wú)梯度算法。文獻(xiàn)[45]深入討論了Kriging模型超參數(shù)訓(xùn)練方法。

1.5 Kriging模型具體實(shí)現(xiàn)的幾個(gè)關(guān)鍵性問(wèn)題

采用上述理論編寫(xiě)一個(gè)可用的Kriging模型計(jì)算機(jī)程序并非難事,例如較早流行的DACE工具箱[69]就是用MATLAB編寫(xiě)，語(yǔ)句較少。但是，不同 Kriging模型程序的性能卻差異很大。根據(jù)作者多年成功開(kāi)發(fā)Kriging模型程序的經(jīng)驗(yàn)，這里介紹影響Kriging模型的效率和魯棒性的幾個(gè)關(guān)鍵問(wèn)題，供讀者參考。

首先，需要對(duì)樣本點(diǎn)進(jìn)行歸一化。所謂歸一化，也就是將樣本點(diǎn)各設(shè)計(jì)變量的取值均轉(zhuǎn)化到 某個(gè)統(tǒng)一的區(qū)間(例如[0,1]或[-1,1])上：

(20)

其次，必須合理地進(jìn)行模型參數(shù)的優(yōu)化。合理的模型參數(shù)，可顯著提高Kriging模型近似精度，使之明顯優(yōu)于多項(xiàng)式響應(yīng)面、徑向基函數(shù)等代理模型。但是，不恰當(dāng)?shù)哪Ｐ蛥?shù)值，可能使Kriging模型預(yù)測(cè)精度差，甚至建立失敗。

再次，需要在任何情況下保證Kriging模型相關(guān)矩陣R的正定性。在某些極限情況下，例如兩個(gè)樣本點(diǎn)無(wú)限靠近時(shí)，相關(guān)矩陣的兩行或兩列近似相等，出現(xiàn)奇異性。即便不出現(xiàn)上述極限情況，矩陣的條件數(shù)可能很大，使得求解Kriging線(xiàn)性方程組的誤差較大，導(dǎo)致模型預(yù)測(cè)精度降低甚至預(yù)測(cè)失敗。為了增加魯棒性，可采用所謂的“正則化”方法(或加入所謂的Nugget效應(yīng)[70])，在矩陣R對(duì)角線(xiàn)上加一個(gè)小的常數(shù)：

R′=R+γI

(21)

γ=(103+n)εM

(22)

其中：εM=2.22×10-16為雙精度時(shí)的機(jī)器零(不同的計(jì)算機(jī)可能不同，讀者可以自己測(cè)試)。另外，γ取值越大，Kriging模型越趨近于回歸模型, 從而可濾出數(shù)值噪聲[71]。

最后，需要采用合理的矩陣分解方法。就建立Kriging模型本身的計(jì)算量而言，最耗時(shí)部分是相關(guān)矩陣求逆。為了提高效率，一般不需要直接求逆，而是求出R-1與矢量的乘積(如R-1yS或R-1r)即可。一種比較有效的方法是采用LU(Lower-upper)分解；但由于相關(guān)矩陣對(duì)稱(chēng)正定，因此可采用效率更高的Cholesky分解：

R=LLT

(23)

(24)

式中：i,j=1,2,…,n。

理論分析和實(shí)踐表明，Cholesky分解的效率比LU分解提高1倍。對(duì)于具有n個(gè)樣本點(diǎn)的問(wèn)題，Cholesky分解的浮點(diǎn)計(jì)算次數(shù)(FLOPS)和內(nèi)存占用量約為

FLOPS≈(1/3)n3+2n2

Memory size≈8×10-6n2MB

(25)

2 Kriging模型理論和算法研究新進(jìn)展

經(jīng)過(guò)調(diào)研，認(rèn)為Kriging模型理論和算法的發(fā)展及其在航空航天工程設(shè)計(jì)領(lǐng)域的應(yīng)用共經(jīng)歷了4個(gè)階段：① 20世紀(jì)50年代初至70年代初，是Kriging模型思想的提出和理論形成階段[22，41-42]；② 20世紀(jì)70年代初至80年代后期，是Kriging模型理論和算法的發(fā)展和成熟階段[23]；③ 20世紀(jì)90年代初至20世紀(jì)末，是Kriging模型在航空航天領(lǐng)域初步應(yīng)用的階段[72-78]；特別是隨著MDO研究的興起，代理模型技術(shù)得到高度重視，Kriging模型成為最具代表性的代理模型之一；④ 21世紀(jì)是Kriging模型理論和算法的進(jìn)一步完善[43-45,71,79]，并在航空航天領(lǐng)域進(jìn)一步應(yīng)用的階段[80-81]。

這里結(jié)合本文的研究工作，重點(diǎn)介紹21世紀(jì)以來(lái)，在Kriging模型的研究和應(yīng)用過(guò)程中，出現(xiàn)的具有代表性的3個(gè)研究新進(jìn)展。

2.1 梯度增強(qiáng)型Kriging模型

梯度信息可用于提高Kriging模型精度[82-85]。而如果采用Adjoint方法[86]等快速求解梯度方法，還可提高建立Kriging模型的效率[87-88]。與傳統(tǒng)有限差分法相比，Adjoint方法的優(yōu)點(diǎn)在于只需求解一次流動(dòng)控制方程和伴隨方程[86]，便可獲得對(duì)所有設(shè)計(jì)變量的梯度，梯度計(jì)算效率與設(shè)計(jì)變量數(shù)無(wú)關(guān)。利用梯度信息來(lái)提高Kriging模型的精度，成為一種新的代理模型方法，稱(chēng)為梯度增強(qiáng)型Kriging (Gradient-Enhanced Kriging, GEK)模型

實(shí)現(xiàn)引入梯度信息的最簡(jiǎn)單做法是利用一階泰勒展開(kāi)：

(26)

式中：i=1,2，…,n；k=1,2，…,m。

將某樣本點(diǎn)處的m個(gè)偏導(dǎo)數(shù)值變成m個(gè)附加樣本點(diǎn)的函數(shù)值，然后再基于這n+n×m個(gè)樣本點(diǎn)建立Kriging模型。該方法缺點(diǎn)是：其近似精度與步長(zhǎng)Δxk密切相關(guān)；步長(zhǎng)太大，數(shù)值誤差較大；步長(zhǎng)太小，容易出現(xiàn)相關(guān)矩陣的奇異性。該方法是一種間接引入梯度信息方法[84]。

這里介紹一種直接引入梯度信息的方法。需要說(shuō)明的是，除了一階導(dǎo)數(shù)信息外，二階導(dǎo)數(shù)信息(Hessian)也可用于提高代理模型的精度[89]。由于獲得二階導(dǎo)數(shù)的計(jì)算代價(jià)較大，且對(duì)于提高代理模型精度的作用有限，因而下文中主要介紹如何在Kriging模型中引入一階偏導(dǎo)數(shù)信息。

對(duì)于一個(gè)具有m個(gè)設(shè)計(jì)變量的優(yōu)化問(wèn)題，首先假設(shè)對(duì)目標(biāo)函數(shù)(或狀態(tài)變量)y在n個(gè)抽樣位置，獲得n個(gè)函數(shù)值及n×m個(gè)偏導(dǎo)數(shù)值。則式(1)和式(2)中的抽樣位置及函數(shù)響應(yīng)值相應(yīng)地變?yōu)?/p>

S=[x(1)…x(n)x(1)…x(1)…

(27)

由式(27)可知，每一個(gè)偏導(dǎo)數(shù)信息都被看作是一個(gè)獨(dú)立的樣本信息。這種表述方法具有一般性：如果在某些樣本點(diǎn)處梯度信息不可用，便從上述樣本數(shù)據(jù)集的相應(yīng)位置上去除即可；如果沒(méi)有任何偏導(dǎo)數(shù)信息，GEK模型將退化為傳統(tǒng)的Kriging模型。因此，GEK模型是一種可考慮梯度信息的更一般化的Kriging模型。

GEK模型對(duì)未知函數(shù)的預(yù)估值定義為所有抽樣函數(shù)值和偏導(dǎo)數(shù)值的線(xiàn)性加權(quán)，即

(28)

Cov[Z(x(i)),Z(x(j))]=σ2R(x(i),x(j))

(29)

與傳統(tǒng)的Kriging模型類(lèi)似，可推導(dǎo)出GEK模型對(duì)未知函數(shù)的預(yù)估值表達(dá)式為

(30)

(31)

(32)

其中：相關(guān)矩陣R、一階導(dǎo)數(shù)矩陣?R以及二階導(dǎo)數(shù)矩陣?2R定義為

(33)

相關(guān)矢量r及其一階導(dǎo)數(shù)?r定義為

(34)

GEK模型對(duì)預(yù)估值的均方差估計(jì)為

(35)

在構(gòu)造GEK模型過(guò)程中，模型參數(shù)訓(xùn)練方法見(jiàn)文獻(xiàn)[88，90]。

需要說(shuō)明的是，GEK模型相關(guān)矩陣中還涉及相關(guān)函數(shù)的一階和二階導(dǎo)數(shù)。對(duì)于三次樣條函數(shù)，可推導(dǎo)出其一階導(dǎo)數(shù)為

(36)

(37)

且Rk和ξl的計(jì)算參見(jiàn)式(16);sign()為取符號(hào)運(yùn)算。

同理，可推導(dǎo)出其二階導(dǎo)數(shù)表達(dá)式為

(38)

(39)

2.2 CoKriging模型

CoKriging模型是20世紀(jì)70年代發(fā)展起來(lái)的一種更有效的地質(zhì)統(tǒng)計(jì)學(xué)插值模型[91-95]。在地質(zhì)統(tǒng)計(jì)學(xué)領(lǐng)域，為了提高對(duì)某個(gè)抽樣比較困難的量的預(yù)測(cè)精度，提出了采用更容易抽樣的量進(jìn)行輔助預(yù)測(cè)的Kriging模型，稱(chēng)為CoKriging模型[92]。2000年，Kennedy和O’Hagan[96]首次將CoKriging模型推廣應(yīng)用于工程科學(xué)領(lǐng)域，發(fā)展了一種采用低可信度計(jì)算機(jī)程序結(jié)果，來(lái)輔助預(yù)測(cè)高可信度計(jì)算機(jī)程序結(jié)果的CoKriging方法。目前國(guó)際上對(duì)CoKriging模型的研究主要集中在地質(zhì)統(tǒng)計(jì)學(xué)和數(shù)學(xué)統(tǒng)計(jì)學(xué)等領(lǐng)域，在航空航天等工程科學(xué)領(lǐng)域的研究也逐漸得到重視[14，97-101]。

2010年，文獻(xiàn)[101]獨(dú)立提出了一種可用于建立變可信度代理模型的實(shí)用CoKriging模型(見(jiàn)文獻(xiàn)[14])。下面對(duì)該模型進(jìn)行介紹。

對(duì)于一個(gè)具有m個(gè)設(shè)計(jì)變量的優(yōu)化問(wèn)題，在設(shè)計(jì)空間中同時(shí)采用高可信度分析(例如Navier-Stokes方程或密網(wǎng)格數(shù)值模擬)和低可信度分析(例如Euler方程或疏網(wǎng)格數(shù)值模擬)進(jìn)行抽樣，并建立所謂的變可信度代理模型。更多建立變可信度模型的方法請(qǐng)參見(jiàn)文獻(xiàn)[88]。變可信度代理模型在達(dá)到相同近似精度的條件下，可顯著提高建立代理模型的效率。

假設(shè)高、低可信度分析程序的抽樣位置為

(40)

式中：下標(biāo)“1”和“2”分別代表高、低可信度，例如n1和n2分別代表高、低可信度樣本點(diǎn)數(shù)(假設(shè)n2?n1)。相應(yīng)的目標(biāo)函數(shù)或約束函數(shù)的響應(yīng)值為

(41)

CoKriging模型預(yù)估值定義為

(42)

式中：λ1、λ2分別為對(duì)高、低可信度響應(yīng)值的加權(quán)系數(shù)。假設(shè)存在分別與y1、y2對(duì)應(yīng)的2個(gè)靜態(tài)隨機(jī)過(guò)程：

(43)

則設(shè)計(jì)空間不同位置處，隨機(jī)變量之間的協(xié)方差和交叉協(xié)方差定義為

(44)

(45)

(46)

且有

(47)

Cokriging模型預(yù)估值的均方差為

(FTR-1F)-1(φ-FTR-1r)

(48)

(49)

剩余參數(shù)θ(11)、θ(12)、θ(22)沒(méi)有解析最優(yōu)解，可以通過(guò)數(shù)值優(yōu)化方法求解式(50)的最大值得到

(50)

需要說(shuō)明的，上述CoKriging模型也可引入梯度信息，得到梯度增強(qiáng)CoKriging模型(GECK，見(jiàn)文獻(xiàn)[101]附錄部分)。

由于可合理假定高、低可信度數(shù)值分析結(jié)果具有相似的函數(shù)變化趨勢(shì)，為簡(jiǎn)化起見(jiàn)，可假定

(51)

式中：h為空間距離；ρ為避免出現(xiàn)矩陣奇異接近1的常數(shù)，可取ρ=0.999 9。

2.3 分層 Kriging 模型

首先，假設(shè)y1對(duì)應(yīng)的隨機(jī)過(guò)程定義為

(52)

式中：β0為待定的自適應(yīng)因子，它反映了高、低可信度函數(shù)之間的相關(guān)性；如果β0=0，說(shuō)明高、低可信度函數(shù)之間沒(méi)有相關(guān)性。β0的引入，有效避免了引入低可信度分析后代理模型精度反而變差的風(fēng)險(xiǎn)。采用與Kriging模型類(lèi)似的推導(dǎo)方法(詳細(xì)推導(dǎo)過(guò)程見(jiàn)文獻(xiàn)[13])，可得到HK模型的預(yù)估值為

(53)

式中：

(54)

預(yù)估值的均方差為

(55)

HK 模型的參數(shù)訓(xùn)練方法與第1節(jié)中Kriging模型基本相同，這里不再贅述。需要說(shuō)明的是，如果存在梯度信息，可采用2.1節(jié)類(lèi)似的方式引入梯度信息，形成梯度增強(qiáng)型分層Kriging模型。

2.4 有關(guān)Kriging模型的其他研究新進(jìn)展

限于篇幅，這里簡(jiǎn)要介紹關(guān)于Kriging模型的其他代表性研究進(jìn)展。文獻(xiàn)[71]采用最大似然估計(jì)對(duì)式(21)的正則化參數(shù)γ進(jìn)行優(yōu)化，發(fā)展了可濾除數(shù)值噪聲的Kriging模型。此外，傳統(tǒng) Kriging模型采用預(yù)設(shè)的全局趨勢(shì)模型,這樣的近似精度可能不是最佳的。于是，文獻(xiàn)[102]發(fā)展了一種采用貝葉斯理論框架來(lái)識(shí)別其趨勢(shì)模型的方法，稱(chēng)為“盲Kriging”模型。為了避免識(shí)別過(guò)程陷入局部最優(yōu)，文獻(xiàn)[103-104]分別采用遺傳算法和交叉驗(yàn)證來(lái)尋找Kriging模型最優(yōu)趨勢(shì)模型，發(fā)展了所謂的“動(dòng)態(tài)Kriging”模型。針對(duì)傳統(tǒng)Kriging模型無(wú)法處理具有隨機(jī)性的計(jì)算機(jī)實(shí)驗(yàn)的問(wèn)題，文獻(xiàn)[105]引入固有不確定性(Intrinsic Uncertainty) 和非固有不確定性(Extrinsic Uncertainty)的概念，并提出了1種所謂的“隨機(jī)Kriging模型”，可用于魯棒優(yōu)化設(shè)計(jì)問(wèn)題。

3 基于Kriging模型的代理優(yōu)化算法

第1節(jié)和第2節(jié)中介紹了如何建立Kriging代理模型。本節(jié)討論將其應(yīng)用于求解式(56)一般性?xún)?yōu)化問(wèn)題的研究新進(jìn)展：

Miny(x)

(56)

式中：y(x)、gi(x)為目標(biāo)函數(shù)與約束函數(shù)；NC為約束個(gè)數(shù)；xu、xl為設(shè)計(jì)變量的上、下限。

本節(jié)首先通過(guò)與傳統(tǒng)優(yōu)化算法比較，梳理提煉出了代理優(yōu)化算法的優(yōu)化機(jī)制和算法框架。其次，給出了“優(yōu)化加點(diǎn)準(zhǔn)則”和“子優(yōu)化”的概念，并針對(duì)常用優(yōu)化加點(diǎn)準(zhǔn)則，介紹了特定的約束處理和混合子優(yōu)化技術(shù)；同時(shí)，還介紹了1種改進(jìn)的加點(diǎn)準(zhǔn)則。最后，給出了優(yōu)化終止條件。

3.1 代理優(yōu)化算法框架

將代理模型用于求解優(yōu)化設(shè)計(jì)問(wèn)題，最初采用的方法是在建立具有合理精度的代理模型后，采用代理模型去替代比較費(fèi)時(shí)的精確數(shù)值模擬分析，進(jìn)行優(yōu)化設(shè)計(jì)。很顯然，這種替代算法很大程度上依賴(lài)于代理模型精度；如果精度過(guò)低，將可能導(dǎo)致優(yōu)化效果不佳甚至失敗。這種方法也被稱(chēng)為第1代方法。第1代方法盡管被人們廣為接受，但由于沒(méi)有形成獨(dú)立的優(yōu)化機(jī)制，因而無(wú)需討論其算法框架。后來(lái)，人們發(fā)展了根據(jù)一定的準(zhǔn)則加入新的樣本點(diǎn)，循環(huán)更新代理模型的第2代方法，并逐步發(fā)展成為今天的代理優(yōu)化算法。但直到目前，國(guó)際上還沒(méi)有討論代理優(yōu)化算法的通用算法框架的文獻(xiàn)報(bào)道。

于是，在文獻(xiàn)調(diào)研和研究實(shí)踐基礎(chǔ)上，本文對(duì)傳統(tǒng)優(yōu)化方法和代理優(yōu)化方法[15]的算法框架進(jìn)行了比較(圖2和圖3)。不難看出，代理優(yōu)化算法的基本原理和過(guò)程可描述為：① 對(duì)設(shè)計(jì)空間進(jìn)行試驗(yàn)設(shè)計(jì)抽樣并運(yùn)行精確數(shù)值模擬分析獲得響應(yīng)值，建立初始代理模型；② 基于代理模型，按照一定的優(yōu)化加點(diǎn)準(zhǔn)則(如最小化代理模型預(yù)測(cè)(MSP)、改善期望(EI)、改善概率(PI)、均方差(MSE)、置信下界(LCB)等)，采用傳統(tǒng)優(yōu)化算法求解相應(yīng)的子優(yōu)化問(wèn)題，以很小的計(jì)算代價(jià)，對(duì)最優(yōu)解進(jìn)行預(yù)測(cè)；③ 對(duì)這些預(yù)測(cè)的最優(yōu)解再次運(yùn)行精確數(shù)值分析，并將結(jié)果作為新的樣本數(shù)據(jù)添加到現(xiàn)有數(shù)據(jù)集中，不斷更新代理模型，直到所產(chǎn)生的樣本點(diǎn)序列收斂于局部或全局最優(yōu)解。

圖2 傳統(tǒng)的梯度或無(wú)梯度優(yōu)化框架Fig.2 Framework of conventional gradient-based and gradient-free optimizations

圖3 代理優(yōu)化算法框架Fig.3 Framework of SBO algorithm

圖4 按通用優(yōu)化算法重新整理的代理優(yōu)化算法框架Fig.4 Rearranged framework of SBO algorithm using general optimization algorithm

為了更好地認(rèn)識(shí)代理優(yōu)化算法，對(duì)圖3進(jìn)行重新整理，將“建立代理模型和根據(jù)加點(diǎn)準(zhǔn)則的子優(yōu)化”過(guò)程看作一個(gè)獨(dú)立模塊，得到圖4。比較圖4 和圖2，不難發(fā)現(xiàn)，除優(yōu)化器(Optimizer)模塊外，二者完全相同。這就是說(shuō)，“建立代理模型和根據(jù)加點(diǎn)準(zhǔn)則的子優(yōu)化”是代理優(yōu)化算法的優(yōu)化機(jī)制。

從代理優(yōu)化算法的原理不難看出，初始樣本的選擇、代理模型及其訓(xùn)練、優(yōu)化加點(diǎn)準(zhǔn)則及其子優(yōu)化求解，是代理優(yōu)化算法的3大要素。

3.2 試驗(yàn)設(shè)計(jì)與初始樣本選擇

代理優(yōu)化算法的第1步是選擇初始樣本點(diǎn)，并建立初始代理模型。雖然理論上可以像梯度優(yōu)化算法那樣只給出1個(gè)初始點(diǎn)，但針對(duì)全局優(yōu)化問(wèn)題，更好的方法是通過(guò)試驗(yàn)設(shè)計(jì)(DoE)選取1組樣本。

試驗(yàn)設(shè)計(jì)的思想是通過(guò)選取最少的樣本點(diǎn)，使獲取的關(guān)于未知設(shè)計(jì)空間的信息最大化。Giunta等在文獻(xiàn)[39]中將DoE方法分為兩類(lèi)：經(jīng)典試驗(yàn)設(shè)計(jì)和現(xiàn)代試驗(yàn)設(shè)計(jì)方法。經(jīng)典試驗(yàn)設(shè)計(jì)方法包括全因子設(shè)計(jì)、中心組合設(shè)計(jì)(CCD)、Box-behnken設(shè)計(jì)、D優(yōu)化(D-Optimal)方法等。經(jīng)典試驗(yàn)設(shè)計(jì)方法主要用于安排儀器實(shí)驗(yàn),并考慮到如何減小實(shí)驗(yàn)隨機(jī)誤差的影響?，F(xiàn)代試驗(yàn)設(shè)計(jì)方法包括擬蒙特卡羅方法、準(zhǔn)蒙特卡羅方法、拉丁超立方、正交試驗(yàn)設(shè)計(jì)(OAD)、哈默斯利序列采樣方法(HSS)。中國(guó)方開(kāi)泰教授發(fā)明的均勻設(shè)計(jì)[40]也屬于現(xiàn)代試驗(yàn)設(shè)計(jì)方法的范疇。現(xiàn)代試驗(yàn)設(shè)計(jì)主要采用“空間填充”的思想，用于安排確定性的計(jì)算機(jī)試驗(yàn)，其中尤以拉丁超立方和均勻設(shè)計(jì)方法比較流行。圖5給出了針對(duì)2個(gè)變量(V1和V2)選取25個(gè)樣本點(diǎn)的示意圖。

圖5 采用拉丁超立方和均勻設(shè)計(jì)針對(duì)2維設(shè)計(jì)空間選擇的25個(gè)樣本點(diǎn)示意圖Fig.5 Schematics of 25 sample points selected by Latin hypercube sampling and uniform design for a 2D design space

不同試驗(yàn)設(shè)計(jì)方法選取的樣本點(diǎn)不同，導(dǎo)致初始代理模型的近似精度不同，從而對(duì)代理優(yōu)化的效率有影響。同樣影響優(yōu)化效率的是初始樣本點(diǎn)數(shù)目。文獻(xiàn)[31，37]討論了樣本點(diǎn)數(shù)的選擇。對(duì)于傳統(tǒng)代理模型優(yōu)化方法，必須使代理模型具有足夠精度。因而一般初始樣本點(diǎn)數(shù)與后期增加的樣本點(diǎn)數(shù)的比值在2∶1以上。例如對(duì)于二次響應(yīng)面方法，對(duì)于m維問(wèn)題的初始樣本點(diǎn)數(shù)必須大于m(m+1)/2。而對(duì)于基于Kriging模型的代理優(yōu)化算法，初始樣本點(diǎn)數(shù)理論上不受設(shè)計(jì)空間維數(shù)的限制，且優(yōu)化效率對(duì)初始樣本點(diǎn)數(shù)的依賴(lài)也并不明顯。一般情況下初始樣本點(diǎn)數(shù)與后期增加的樣本點(diǎn)數(shù)之比在1∶2 以下。

3.3 優(yōu)化加點(diǎn)準(zhǔn)則及其子優(yōu)化求解

建立初始代理模型后，下一步是通過(guò)一定的法則循環(huán)選擇新的樣本點(diǎn)，直到優(yōu)化收斂。所謂“優(yōu)化加點(diǎn)準(zhǔn)則”[15,31,71,106]，是指如何由所建立代理模型去產(chǎn)生新的樣本點(diǎn)的法則或規(guī)則。所謂“子優(yōu)化”，是相對(duì)主優(yōu)化而言，指采用傳統(tǒng)優(yōu)化算法求解由加點(diǎn)準(zhǔn)則所確定的優(yōu)化問(wèn)題，得到新的樣本點(diǎn)的過(guò)程。主優(yōu)化加點(diǎn)循環(huán)中的每一步，都要進(jìn)行一次完整的子優(yōu)化迭代，直到子優(yōu)化收斂。但由于在子優(yōu)化中，無(wú)需訪問(wèn)精確數(shù)值分析，因此計(jì)算時(shí)間可以忽略。

針對(duì)基于Kriging模型的代理優(yōu)化算法，國(guó)際上已經(jīng)發(fā)展了多種加點(diǎn)準(zhǔn)則[31,106-107]，包括MSP準(zhǔn)則[15，108]、EI準(zhǔn)則[15，20，74]、PI準(zhǔn)則[106-107，109]、MSE準(zhǔn)則[106，109]、LCB準(zhǔn)則[75，106-107]。為了說(shuō)明這5種常見(jiàn)加點(diǎn)準(zhǔn)則的原理和子優(yōu)化問(wèn)題的建立，下面以某一維函數(shù)為例，采用4個(gè)樣本點(diǎn)S=[0 0.4 0.6 1.0]T建立Kriging模型。該一維函數(shù)來(lái)自文獻(xiàn)[31]，表達(dá)式為

y=(6x-2)2sin(12x-4)x∈[0,1]

(57)

3.3.1 最小化代理模型預(yù)測(cè)準(zhǔn)則

最小化代理模型預(yù)測(cè)準(zhǔn)則是最簡(jiǎn)單、最直接，也是最早被采用的方法[15，106-110]。其原理是直接在代理模型上尋找目標(biāo)函數(shù)的最小值。帶約束的子優(yōu)化問(wèn)題數(shù)學(xué)模型為

(58)

圖6 一維問(wèn)題MSP加點(diǎn)法則原理圖Fig.6 Schematic of MSP infill-sampling criterion for a 1D problem

對(duì)于全局優(yōu)化問(wèn)題，為了提高子優(yōu)化的精度，本文作者最近提出了一種遺傳算法與梯度優(yōu)化算法相結(jié)合的混合子優(yōu)化方法。將多種優(yōu)化算法獲得的目標(biāo)函數(shù)值最小的解作為式(58)子優(yōu)化的最終優(yōu)化解。對(duì)于遺傳算法，可采用文獻(xiàn)[111]提出的約束處理方法。在該方法中，遺傳算法的適應(yīng)度函數(shù)值為

(59)

對(duì)于局部?jī)?yōu)化問(wèn)題，可采用當(dāng)前最優(yōu)點(diǎn)作為初始點(diǎn)，直接采用BFGS、SQP等梯度優(yōu)化算法[61]或Hooke & Jeeves[67]、Simplex[68]等局部搜索算法進(jìn)行子優(yōu)化。另外，為了確保代理優(yōu)化算法收斂性，文獻(xiàn)[38，112]還給出了一種信賴(lài)域方法。即在初始點(diǎn)附近設(shè)定一個(gè)初始的信賴(lài)域半徑r0，建立代理模型后，子優(yōu)化只在信賴(lài)域內(nèi)進(jìn)行，即搜索范圍限制為

x∈[max(x0-r0,xl),min(x0+r0,xu)]

(60)

式中：x0為優(yōu)化初始點(diǎn)(或樣本點(diǎn)中滿(mǎn)足約束的最優(yōu)點(diǎn))。對(duì)子優(yōu)化結(jié)果，采用精確數(shù)值模擬分析并更新代理模型，同時(shí)調(diào)整信賴(lài)域半徑大小，直到整個(gè)優(yōu)化收斂。文獻(xiàn)[38，112]給出的信賴(lài)域半徑調(diào)整方法為

(61)

式中：質(zhì)量因子δ′為目標(biāo)函數(shù)的實(shí)際改善與預(yù)測(cè)改善的比值。δ′的合理計(jì)算是其中的關(guān)鍵；如果計(jì)算不合理，可能使信賴(lài)域半徑縮小過(guò)快，從而導(dǎo)致優(yōu)化出現(xiàn)早熟。因此，還有待于發(fā)展更好的質(zhì)量因子計(jì)算方法或其他提高收斂性的方法。

3.3.2 改善期望準(zhǔn)則

(62)

對(duì)于最小化問(wèn)題，目標(biāo)函數(shù)改善量I(x)定義為

(63)

I(x)的期望值為[15,74,107]

E[I(x)]=

(64)

圖7 某一維問(wèn)題EI函數(shù)分布示意圖 Fig.7 Schematic of EI function distribution for a 1D problem

通過(guò)求解最大化EI值的子優(yōu)化問(wèn)題：

MaxE[I(x)]

s.t.xl≤x≤xu

(65)

從而可得到新的樣本點(diǎn)x*。

(66)

這樣可以得到所謂的約束EI值為[109-110]

Ec[I(x)]=E[I(x)∩Gi≤0]=

(67)

于是，可采用3.3.1節(jié)中提出的混合子優(yōu)化方法，求解無(wú)約束子優(yōu)化問(wèn)題：

s.t.xl≤x≤xu

(68)

從而得到新的樣本點(diǎn)x*。

EI 方法理論上是一種全局優(yōu)化算法[74，107]，但優(yōu)化后期收斂較慢。于是，本文作者最近提出了一種所謂的“局部EI”方法。其基本思想是：優(yōu)化過(guò)程中，代理優(yōu)化的子優(yōu)化不必找到整個(gè)設(shè)計(jì)空間內(nèi)EI函數(shù)的最大值，而只需在當(dāng)前最優(yōu)點(diǎn)附近尋找一個(gè)較大的值即可。這種方法是受到梯度優(yōu)化算法的啟發(fā)。眾所周知，梯度優(yōu)化過(guò)程中每一步的線(xiàn)搜索過(guò)程不一定要精確進(jìn)行，只需搜索有限步數(shù)即可(滿(mǎn)足所謂的“Armijo準(zhǔn)則”)。因而不難想到, 對(duì)于EI最大值搜索也不需要精確進(jìn)行，而是在一定范圍內(nèi)搜索即可。實(shí)踐證明，這樣不僅能顯著提高收斂效率，而且能保持一定的全局優(yōu)化性能，對(duì)于實(shí)際的工程設(shè)計(jì)問(wèn)題可能是一種比較好的選擇。

3.3.3 改善概率準(zhǔn)則

改善概率準(zhǔn)則方法[31]實(shí)際上是EI方法派生的方法。與EI準(zhǔn)則類(lèi)似，該準(zhǔn)則通過(guò)計(jì)算目標(biāo)函數(shù)改善的概率，并將改善概率最大點(diǎn)作為新的樣本點(diǎn)。目標(biāo)函數(shù)改善的概率為[106-107，109]

(69)

對(duì)于上述一維問(wèn)題，PI函數(shù)分布如圖8所示。

圖8 某一維問(wèn)題的PI函數(shù)分布示意圖 Fig.8 Schematic of PI function distribution for a 1D problem

對(duì)于約束優(yōu)化問(wèn)題，可采用與約束EI類(lèi)似的方法進(jìn)行處理。于是，可采用3.3.1節(jié)提出的混合子優(yōu)化方法，求解下面的子優(yōu)化問(wèn)題

s.t.xl≤x≤xu

(70)

得到新的樣本點(diǎn)x*。

3.3.4 均方差準(zhǔn)則

均方差準(zhǔn)則(MSE)方法是直接運(yùn)用Kriging模型提供的均方差估計(jì)(或均方根誤差(RMSE),如圖 9所示)的一種改善代理模型全局精度的加點(diǎn)準(zhǔn)則。即采用遺傳算法等全局優(yōu)化算法，并結(jié)合局部?jī)?yōu)化算法，求解下面的子優(yōu)化問(wèn)題[106，109]：

s.t.xl≤x≤xu

(71)

從而得到新樣本點(diǎn)x*。對(duì)于帶約束的問(wèn)題，可引入滿(mǎn)足約束的概率，可求解下面的子優(yōu)化問(wèn)題

s.t.xl≤x≤xu

(72)

得到新的樣本點(diǎn)x*。

圖9 某一維問(wèn)題的均方根誤差函數(shù)分布示意圖Fig.9 Schematic of RMSE function distribution for a 1D problem

3.3.5 置信下界準(zhǔn)則

(73)

圖10 某一維問(wèn)題的LCB函數(shù)分布示意圖Fig.10 Schematic of LCB function distribution for a 1D problem

該準(zhǔn)則約束的處理與MSP準(zhǔn)則類(lèi)似。即可采用3.3.1節(jié)提出的混合子優(yōu)化方法求解下面的約束子優(yōu)化問(wèn)題

(74)

從而得到新的樣本點(diǎn)x*。

文獻(xiàn)[106-107，113]均對(duì)上述優(yōu)化加點(diǎn)準(zhǔn)則進(jìn)行了比較，得到了一些具有指導(dǎo)性的建議。這里再補(bǔ)充幾點(diǎn)說(shuō)明：

1) 除MSP準(zhǔn)則適用于多項(xiàng)式響應(yīng)面、徑向基函數(shù)等其他代理模型外，EI、PI、MSE、LCB準(zhǔn)則一般只適用于Kriging模型以及其他具有誤差估計(jì)功能的代理模型。

2) 除上述的5種比較流行的加點(diǎn)準(zhǔn)則，還可采用目標(biāo)搜索(Goal-seeking)[107]和條件下界等方法[31]。但這些準(zhǔn)則不是通用的加點(diǎn)準(zhǔn)則。

3) 各種加點(diǎn)準(zhǔn)則可組合使用，形成組合加點(diǎn)法則。組合加點(diǎn)不僅可克服單一準(zhǔn)則的缺陷，提高其魯棒性。另外，還可在每一步迭代加入任意點(diǎn)數(shù)的新樣本點(diǎn)，這就是所謂的多點(diǎn)加點(diǎn)法則[114]或并行加點(diǎn)法則[31]。

4) 不僅可用于單目標(biāo)優(yōu)化，也可用于多目標(biāo)優(yōu)化[20，31]。

5) 既可用于全局優(yōu)化，也可用于局部?jī)?yōu)化問(wèn)題。對(duì)于局部?jī)?yōu)化問(wèn)題，為了提高收斂性，可采用信賴(lài)域方法[107]和采用梯度增強(qiáng)型Kriging[113]。

3.4 代理優(yōu)化的終止條件

優(yōu)化終止條件是優(yōu)化算法必須要具備的一個(gè)因素。目前國(guó)內(nèi)外還沒(méi)有關(guān)于定義代理優(yōu)化算法優(yōu)化終止條件的文獻(xiàn)報(bào)道。于是，本文作者最近根據(jù)代理優(yōu)化的特點(diǎn)，并借鑒其他優(yōu)化算法的終止條件，定義了如下4種終止條件。

首先，可以根據(jù)樣本點(diǎn)之間的距離和目標(biāo)函數(shù)值的差別進(jìn)行定義：

(75)

(76)

其次,可根據(jù)代理模型在最優(yōu)點(diǎn)附近的近似精度來(lái)定義：

(77)

再次，可根據(jù)計(jì)算資源情況，設(shè)定調(diào)用精確數(shù)值模擬分析的最大次數(shù)：

N≤Nmax

(78)

例如針對(duì)氣動(dòng)優(yōu)化設(shè)計(jì)問(wèn)題，由于計(jì)算資源有限，一般希望在100～300次CFD計(jì)算以?xún)?nèi)，獲得最優(yōu)的外形，因此可以設(shè)定Nmax=100～300。

最后，可針對(duì)特定的加點(diǎn)準(zhǔn)則進(jìn)行定義。例如針對(duì)EI準(zhǔn)則，可定義為

max(EI)≤EImin

(79)

也就是說(shuō)當(dāng)子優(yōu)化求解得到的最大EI值小于某一個(gè)閥值時(shí)，優(yōu)化即可終止。例如，在氣動(dòng)減阻優(yōu)化中，可設(shè)定EImin=1.0×10-6，也就是說(shuō)如果最大EI值小于0.01 counts時(shí)，優(yōu)化即可終止。

4 基于Kriging模型的代理優(yōu)化算法驗(yàn)證和確認(rèn)

代理模型最初作為一種近似方法，用于提高采用高可信度數(shù)值分析的工程設(shè)計(jì)問(wèn)題的效率[76]。因此，人們往往針對(duì)特定的問(wèn)題發(fā)展特定的代理優(yōu)化框架，例如文獻(xiàn)[20-21]中發(fā)展的翼型優(yōu)化設(shè)計(jì)方法，文獻(xiàn)[52，85]的機(jī)翼優(yōu)化設(shè)計(jì)方法，以及文獻(xiàn)[54]發(fā)展的結(jié)構(gòu)優(yōu)化設(shè)計(jì)方法。但是，代理優(yōu)化算法如果作為一種通用的優(yōu)化算法框架，就需要經(jīng)過(guò)優(yōu)化標(biāo)準(zhǔn)算例的驗(yàn)證和確認(rèn)，以評(píng)估算法的效率和精度，并對(duì)算法進(jìn)行相應(yīng)的改進(jìn)。對(duì)于代理優(yōu)化算法的效率，一般通過(guò)目標(biāo)函數(shù)計(jì)算次數(shù)來(lái)評(píng)估[107]；因?yàn)樵趯?shí)際工程設(shè)計(jì)問(wèn)題中，調(diào)用高可信度數(shù)值模擬來(lái)計(jì)算目標(biāo)函數(shù)往往是整個(gè)優(yōu)化流程中最費(fèi)時(shí)的部分。目前，關(guān)于代理優(yōu)化算法的驗(yàn)證和確認(rèn)的研究工作開(kāi)展得比較少。

4.1 標(biāo)準(zhǔn)測(cè)試函數(shù)優(yōu)化算例驗(yàn)證

Jones等在1998年提出著名的EI方法時(shí)，只采用了4種比較簡(jiǎn)單的測(cè)試函數(shù)算例進(jìn)行驗(yàn)證[74，107]，包括2維Branin-Hoo函數(shù)、Goldstein-Price函數(shù)、3維和6維Hartman函數(shù)，且這些都是無(wú)約束問(wèn)題。Sasena[109]、 Forrester[31]、Parr[110]以及劉俊[106]等在驗(yàn)證約束處理方法時(shí)，也主要采用了2維Branin-Hoo函數(shù)等簡(jiǎn)單算例。但這些函數(shù)其實(shí)并不是優(yōu)化算法領(lǐng)域的最具挑戰(zhàn)性的例子，還不能全面驗(yàn)證代理優(yōu)化作為一種優(yōu)化算法的性能。最近，文獻(xiàn)[113，115]采用了更具代表性和挑戰(zhàn)性的Rosenbrock函數(shù)算例、Rastrigin函數(shù)算例、7維G9函數(shù)算例、5維Himmelblau算例、Pressure Vessel Design等單目標(biāo)優(yōu)化標(biāo)準(zhǔn)算例，以及FON和ZET多種多目標(biāo)優(yōu)化標(biāo)準(zhǔn)算例，對(duì)代理優(yōu)化算法的效率和精度進(jìn)行了系統(tǒng)地驗(yàn)證和確認(rèn)，而且掌握了代理優(yōu)化算法的一些特性，總結(jié)了一些最佳參數(shù)設(shè)置。為了不斷改進(jìn)代理優(yōu)化算法的性能，今后還需要采用更多的標(biāo)準(zhǔn)算例進(jìn)行測(cè)試(例如文獻(xiàn)[116]中給出了175種算例)。限于篇幅，下文僅給出了最具代表性的Rosenbrock函數(shù)和G9函數(shù)的最新測(cè)試結(jié)果。

Rosenbrock函數(shù)算例[117]是測(cè)試優(yōu)化算法的經(jīng)典算例之一。該函數(shù)最優(yōu)點(diǎn)位于一個(gè)狹長(zhǎng)而平坦的峽谷內(nèi)(如圖 11所示)。2維Rosenbrock函數(shù)算例優(yōu)化數(shù)學(xué)模型為

Minf(x1,x2)=(1.0-x1)2+

s.t.x1,x2∈[-2.0,2.0]

(80)

圖12給出了采用不同加點(diǎn)準(zhǔn)則時(shí)，真實(shí)目標(biāo)函數(shù)值的收斂情況?？梢?jiàn)，對(duì)于單獨(dú)使用的加點(diǎn)準(zhǔn)則，EI、局部EI和LCB效果較好，MSP、PI和MSE效果不佳。如果組合使用，“MSP+EI”和“MSP+LCB”的效果均優(yōu)于單獨(dú)使用的效果。圖13 給出了梯度增強(qiáng)型Kriging模型的優(yōu)化收斂情況。可見(jiàn)引入梯度信息后，各種加點(diǎn)準(zhǔn)則的收斂特性均得到改善；尤其是MSP準(zhǔn)則的改善效果最明顯。表 1為采用代理優(yōu)化算法不同加點(diǎn)準(zhǔn)則的優(yōu)化結(jié)果的比較，表明代理優(yōu)化算法不僅可以收斂到無(wú)約束問(wèn)題的精確最優(yōu)解，而且具有很高的效率。

圖11 二維Rosenbrock函數(shù)示意圖Fig.11 Schematic of 2D Rosenbrock function

圖12 代理優(yōu)化算法不同加點(diǎn)準(zhǔn)則收斂歷程的比較(2維 Rosenbrock函數(shù)算例)Fig.12 Comparison of SBO convergence history with different infill sampling criteria for a 2D Rosenbrock function case

為了驗(yàn)證代理優(yōu)化對(duì)更高維問(wèn)題的適用性，可以采用多維Rosenbrock函數(shù)：

(81)

該函數(shù)在3維以上為多極值問(wèn)題。

圖14給出了針對(duì)5～80維函數(shù)，基于Kriging和梯度增強(qiáng)型Kriging模型的代理優(yōu)化結(jié)果的比較?？梢钥闯?，對(duì)于Kriging模型，隨著設(shè)計(jì)變量數(shù)增加，需要的函數(shù)計(jì)算次數(shù)迅速增加；而引入梯度信息后，函數(shù)計(jì)算次數(shù)對(duì)設(shè)計(jì)變量數(shù)變得不敏感。如果用于氣動(dòng)優(yōu)化設(shè)計(jì)問(wèn)題，采用Adjoint方法計(jì)算梯度，則可大大提高代理優(yōu)化對(duì)高維問(wèn)題的適應(yīng)性。但需要說(shuō)明的是，隨著設(shè)計(jì)變量數(shù)增加，模型相關(guān)矩陣的階數(shù)增加(n個(gè)樣本，m個(gè)設(shè)計(jì)變量，矩陣階數(shù)為n+n×m)，從而導(dǎo)致模型訓(xùn)練的計(jì)算量大幅度增加。因此，發(fā)展更高效的模型訓(xùn)練方法或矩陣降階方法，成為今后一個(gè)值得關(guān)注的關(guān)鍵問(wèn)題。

圖13 引入梯度信息后的代理優(yōu)化算法不同加點(diǎn)準(zhǔn)則目標(biāo)函數(shù)收斂歷程的比較(2維Rosenbrock函數(shù)算例)Fig.13 Comparison of SBO convergence history incorporating gradient-information with different infill sampling criteria for a 2D Rosenbrock function case

表1 不同優(yōu)化算法優(yōu)化結(jié)果的比較(2維Rosenbrock函數(shù)算例)Table 1 Comparison of optimization results of different optimization algorithms for a 2D Rosenbrock function case

OptimizationalgorithmOptimalobjectivefunctionvalueNo.offunctionandgradientevaluationsTrue0NotapplicableKrigingMSP1.2358300fKrigingEI8.1681×10-6300fKrigingPI4.4276300fKrigingMSE3.4450×10-2300fKrigingLCB6.8738×10-5300fKrigingMSP+EI2.7596×10-7300fKrigingMSP+LCB1.4811×10-5300fKrigingMSP+MSE1.7529×10-3300fKriginglocalEI1.1496×10-8300fGEKMSP2.3888×10-10200f+200gGEKEI1.7859×10-5200f+200gGEKPI1.1485×10-3200f+200gGEKMSE3.5909200f+200gGEKLCB1.3395×10-5200f+200gGEKMSP+EI1.1157×10-11200f+200gGEKMSP+LCB1.7690×10-8200f+200gGEKMSP+MSE2.7855×10-9200f+200gGEKlocalEI3.7717×10-9200f+200g

Notes: f—function evaluation; g—gradient evaluation

圖14 5～80維Rosenbrock函數(shù)代理優(yōu)化Kriging模型和GEK模型的比較Fig.14 Comparison of Kriging and GEK models for 5-80 dimensional Rosenbrock functions surrogate optimization

G9函數(shù)算例是國(guó)際上用于測(cè)試全局優(yōu)化算法的一個(gè)經(jīng)典算例。該算例含7個(gè)變量，4個(gè)約束(其中g(shù)1和g4為主動(dòng)約束)，優(yōu)化數(shù)學(xué)模型為

Minf(x)=(x1-10)2+5(x2-12)2+

4x6x7-10x6-8x7

(82)

該優(yōu)化問(wèn)題是一個(gè)難度非常大的全局優(yōu)化問(wèn)題[118]，原因在于：① 可行域很小，僅占整個(gè)設(shè)計(jì)空間的0.5%，屬于強(qiáng)約束問(wèn)題；② 目標(biāo)函數(shù)的變化幅度較大，從107量級(jí)到103量級(jí)。目前為止，國(guó)際上能得到的最優(yōu)解為

x*={2.330 499,1.951 372,-0.477 541 4,

4.365 726,-0.624 487 0,1.038 161,

1.594 227}

f(x*)=680.630 057 3

表2中給出了采用Kriging和GEK模型不同加點(diǎn)準(zhǔn)則的優(yōu)化效率比較?？梢?jiàn)，MSP準(zhǔn)則和LCB準(zhǔn)則在未得到全局最優(yōu)解前滿(mǎn)足了優(yōu)化終止條件；EI和局部EI方法的結(jié)果均非常接近全局最優(yōu)；采用不同組合加點(diǎn)準(zhǔn)則，只有“MSP+EI”收斂到了全局最優(yōu)。引入梯度信息后，代理優(yōu)化不同加點(diǎn)準(zhǔn)則的優(yōu)化結(jié)果和優(yōu)化效率均得到顯著改善，但仍然只有EI、局部EI、“MSP+EI”準(zhǔn)則獲得了全局最優(yōu)解。此外，代理優(yōu)化的結(jié)果不僅優(yōu)于直接遺傳算法優(yōu)化的結(jié)果，且調(diào)用真實(shí)函數(shù)分析次數(shù)小了1～2個(gè)量級(jí)。

表3給出了代理優(yōu)化結(jié)果與目前國(guó)際上能得到的最佳結(jié)果的比較，說(shuō)明優(yōu)化結(jié)果已經(jīng)非常接近國(guó)際上最佳結(jié)果，主動(dòng)約束得到精確滿(mǎn)足。事實(shí)上，表中的國(guó)外最佳結(jié)果對(duì)g4約束是有所放寬的。這說(shuō)明代理優(yōu)化算法針對(duì)復(fù)雜全局優(yōu)化問(wèn)題具有很強(qiáng)的約束處理能力。

通過(guò)各種約束和無(wú)約束、局部和全局優(yōu)化的標(biāo)準(zhǔn)算例測(cè)試，充分說(shuō)明代理優(yōu)化算法具有精確收斂到優(yōu)化問(wèn)題的真實(shí)最優(yōu)解的能力。

4.2 氣動(dòng)優(yōu)化設(shè)計(jì)標(biāo)準(zhǔn)算例確認(rèn)

為了考查代理優(yōu)化算法在實(shí)際工程優(yōu)化問(wèn)題中的性能，文獻(xiàn)[119-120]開(kāi)展了針對(duì)標(biāo)準(zhǔn)氣動(dòng)優(yōu)化設(shè)計(jì)算例的確認(rèn)研究，并與同類(lèi)方法進(jìn)行了比較。2015年，AIAA航空宇航科學(xué)大會(huì)(ASM)專(zhuān)門(mén)設(shè)立了氣動(dòng)設(shè)計(jì)優(yōu)化討論組(ADODG)，并定義了5個(gè)標(biāo)準(zhǔn)算例，包括2個(gè)翼型設(shè)計(jì)、2個(gè)機(jī)翼設(shè)計(jì)和1個(gè)全機(jī)優(yōu)化設(shè)計(jì)算例。文獻(xiàn)[119-120]采用了前2個(gè)標(biāo)準(zhǔn)算例：NACA0012和RAE2822翼型減阻優(yōu)化設(shè)計(jì)。雖然NACA0012和RAE2822翼型已廣泛被國(guó)內(nèi)外研究人員作為展示氣動(dòng)優(yōu)化設(shè)計(jì)算法的算例，但這些算例的設(shè)計(jì)狀態(tài)和優(yōu)化模型缺乏統(tǒng)一的定義，從而無(wú)法在同一尺度下評(píng)判所采用的優(yōu)化算法的優(yōu)劣。AIAA氣動(dòng)設(shè)計(jì)優(yōu)化討論組定義的RAE2822翼型的設(shè)計(jì)狀態(tài)為：Ma=0.734，Re=6.5×106，優(yōu)化數(shù)學(xué)模型為

表2 不同優(yōu)化算法優(yōu)化結(jié)果的比較(7維G9函數(shù)算例)Table 2 Comparison of optimization results of different optimization algorithms for a 7D G9 function case

表3 7維G9函數(shù)優(yōu)化結(jié)果 (EI加點(diǎn)準(zhǔn)則)Table 3 Optimum results of a 7D G9 function case(EI infill-sampling criterion)

MinCD

(83)

式中：CL、CD和Cm分別為翼型的升力、阻力和俯仰力矩系數(shù)；A1和A0分別為優(yōu)化翼型和初始翼型的面積。ADODG還要求必須對(duì)基準(zhǔn)翼型和優(yōu)化翼型作網(wǎng)格收斂性研究，并且要證明優(yōu)化過(guò)程充分收斂。

文獻(xiàn)[120]采用CST(Class-function Shape-function Transformation)方法作為翼型外形參數(shù)化方法，并采用8階伯恩斯坦多項(xiàng)式，共18個(gè)設(shè)計(jì)變量。對(duì)初始翼型作了網(wǎng)格收斂性研究(見(jiàn)表4)，最后確定計(jì)算網(wǎng)格數(shù)量大約為13萬(wàn)(513×257)。如圖 15所示，網(wǎng)格遠(yuǎn)場(chǎng)邊界為90c(c為翼型弦長(zhǎng))；物面第1層網(wǎng)格高度為2.0×10-6。為了避免升力等式約束，文獻(xiàn)將升力約束改為升力不減，并固定迎角α=2.879 5°。

表4基準(zhǔn)RAE2822翼型氣動(dòng)性能網(wǎng)格收斂性

Table4GridconvergenceforbaselineRAE2822airfoilaerodynamicperformance

GridsizeCL/countsCD/counts192×9682.4221.66256×12882.4208.10320×16082.4201.85384×19282.4198.65448×22482.4196.53512×25682.4195.02576×28882.4194.09

圖15 RAE2822翼型計(jì)算網(wǎng)格示意圖 Fig.15 Schematic of computational grid for RAE2822 airfoil

表5給出了優(yōu)化翼型的氣動(dòng)性能及面積與基準(zhǔn)翼型的比較。在所有約束都得到精確滿(mǎn)足的前提下，阻力大大降低。與文獻(xiàn)[119]中梯度優(yōu)化結(jié)果的比較說(shuō)明，代理優(yōu)化結(jié)果明顯優(yōu)于梯度優(yōu)化結(jié)果。圖16給出了代理優(yōu)化過(guò)程中目標(biāo)函數(shù)阻力系數(shù)的收斂歷程。共進(jìn)行了3輪優(yōu)化，通過(guò)第2輪和第3輪優(yōu)化，阻力分別進(jìn)一步下降了0.5個(gè)counts和0.2個(gè)counts，說(shuō)明優(yōu)化已經(jīng)充分收斂。圖17給出了優(yōu)化前后翼型的幾何形狀及壓力系數(shù)Cp分布對(duì)比。優(yōu)化后上表面的強(qiáng)激波被削弱成兩道非常弱的激波，優(yōu)化效果顯著。按ADODG的要求還給出了優(yōu)化翼型的網(wǎng)格收斂性研究(表6)，發(fā)現(xiàn)優(yōu)化翼型最終阻力為102.98 counts，這個(gè)結(jié)果優(yōu)于國(guó)際同行的結(jié)果[121]。

表5優(yōu)化前后翼型的氣動(dòng)性能及面積對(duì)比

Table5Comparisonofaerodynamicperformanceandareaofbaselineandoptimumairfoils

AirfoilCL/countsCD/countsCmA1RAE282282.4195.00-0.092130.07794Optimizedbygradientmethod[119]82.4129.40-0.091900.07794OptimizedbySBO82.4104.29-0.088020.07794

圖16 RAE2822翼型減阻優(yōu)化收斂歷程 Fig.16 Convergence history of drag reduction optimization for RAE2822 airfoil

圖17 優(yōu)化前后翼型的壓力分布及幾何形狀比較Fig.17 Comparison of pressure distributions and geometries for baseline and optimum airfoils

表6優(yōu)化翼型氣動(dòng)性能網(wǎng)格收斂性

Table6Gridconvergenceforoptimumairfoilaerodynamicperformance

GridsizeCL/countsCD/counts192×9682.4127.75256×12882.4113.99320×16082.4108.92384×19282.4106.58448×22482.4105.21512×25682.4104.29576×28882.4103.62640×32082.4102.98

為了更進(jìn)一步確認(rèn)代理優(yōu)化算法在三維復(fù)雜外形氣動(dòng)設(shè)計(jì)中的優(yōu)化能力，下一步需要進(jìn)行ADODG定義的算例3[121]、算例4[122-123]和算例5[124]的NASA CRM標(biāo)準(zhǔn)機(jī)翼和全機(jī)模型的代理優(yōu)化設(shè)計(jì)研究，并與梯度優(yōu)化方法[122-123]的結(jié)果進(jìn)行比較。更多新的結(jié)果將可能在2017年的第3屆ADODG研討會(huì)上展示。

4.3 氣動(dòng)/結(jié)構(gòu)綜合優(yōu)化設(shè)計(jì)算例確認(rèn)

具有學(xué)科強(qiáng)耦合關(guān)系的氣動(dòng)/結(jié)構(gòu)綜合優(yōu)化設(shè)計(jì)，是飛機(jī)MDO的最典型問(wèn)題[125]。該算例可以用于確認(rèn)代理優(yōu)化算法在具有學(xué)科強(qiáng)耦合關(guān)系的飛行器精細(xì)化優(yōu)化設(shè)計(jì)中的適用性。文獻(xiàn)[81,126]采用代理優(yōu)化算法研究了某高亞聲速運(yùn)輸機(jī)翼身組合體(如圖 18所示)機(jī)翼氣動(dòng)/結(jié)構(gòu)綜合優(yōu)化。但受到當(dāng)時(shí)技術(shù)條件限制，只采用了8個(gè)氣動(dòng)和結(jié)構(gòu)學(xué)科的設(shè)計(jì)變量。最近，文獻(xiàn)作者又采用最新的代理優(yōu)化算法，得到了23個(gè)設(shè)計(jì)變量下的優(yōu)化設(shè)計(jì)結(jié)果。優(yōu)化數(shù)學(xué)模型為

MinW

(84)

式中：W為機(jī)翼質(zhì)量；L和L/D分別為組合體升力和升阻比；σmax和δmax分別為機(jī)翼結(jié)構(gòu)的最大應(yīng)力和變形。該優(yōu)化是對(duì)巡航狀態(tài)下的機(jī)翼進(jìn)行設(shè)計(jì)：通過(guò)調(diào)整迎角、機(jī)翼尖梢比、線(xiàn)性扭轉(zhuǎn)角以及上下蒙皮沿展向各段厚度(各10個(gè)變量)共23個(gè)設(shè)計(jì)變量，在滿(mǎn)足升力、升阻比以及最大應(yīng)力和翼尖最大變形約束的條件下，獲得最小的機(jī)翼結(jié)構(gòu)重量。該優(yōu)化的機(jī)翼氣動(dòng)/結(jié)構(gòu)分析模型參見(jiàn)文獻(xiàn)[81，126]：氣動(dòng)分析采用全速勢(shì)有黏/無(wú)黏迭代的方法；采用ANSYS軟件進(jìn)行結(jié)構(gòu)分析；用弱耦合法進(jìn)行靜氣動(dòng)彈性分析。

圖18 高亞聲速運(yùn)輸機(jī)翼身組合體外形圖 Fig.18 Schematic of a wing-body configuration for high-subsonic transport aircraft

表7給出了優(yōu)化機(jī)翼相對(duì)于基準(zhǔn)機(jī)翼的性能變化，圖19給出了代理模型優(yōu)化過(guò)程中目標(biāo)函數(shù)的收斂歷程。優(yōu)化以200個(gè)樣本點(diǎn)為基礎(chǔ)建立代理模型，加點(diǎn)準(zhǔn)則采用 “MSP+EI”組合方法?？梢钥闯觯趪?yán)格滿(mǎn)足升力、升阻比、最大應(yīng)力和變形約束的情況下，機(jī)翼重量降低了28.9%。這說(shuō)明所采用的代理優(yōu)化方法不僅優(yōu)化效率高(23個(gè)設(shè)計(jì)變量問(wèn)題只進(jìn)行不到300次的多學(xué)科耦合分析)、優(yōu)化效果良好，且能夠精確處理氣動(dòng)和結(jié)構(gòu)不同學(xué)科的約束。圖 20給出了優(yōu)化后的機(jī)翼變形圖。圖 21給出了優(yōu)化前后機(jī)翼上下蒙皮厚度尺寸的對(duì)比，優(yōu)化后的機(jī)翼厚度尺寸大大降低，從翼根到翼梢厚度變化更加平緩，優(yōu)化結(jié)果比較合理，優(yōu)化效果非常顯著。

通過(guò)該算例確認(rèn)，表明代理優(yōu)化方法在處理具有復(fù)雜耦合關(guān)系的飛行器復(fù)雜多學(xué)科精細(xì)化優(yōu)化設(shè)計(jì)問(wèn)題方面具有良好的發(fā)展前景。但是，隨著設(shè)計(jì)變量數(shù)的增加，優(yōu)化的計(jì)算量將會(huì)迅速增加。而引入氣動(dòng)/結(jié)構(gòu)耦合系統(tǒng)Adjoint梯度[127]，以及成本更低的低可信度分析，發(fā)展梯度增強(qiáng)型變可信度Kriging代理模型及其相應(yīng)的代理優(yōu)化算法，將是該領(lǐng)域未來(lái)重點(diǎn)發(fā)展的方向之一。代理優(yōu)化算法的研究也將繼續(xù)成為未來(lái) 10～20年內(nèi)飛行器精細(xì)優(yōu)化設(shè)計(jì)領(lǐng)域的研究熱點(diǎn)。

表7優(yōu)化前后機(jī)翼外形及性能對(duì)比

Table7Comparisonofshapeandperformanceofbaselineandoptimumwing

WingperformanceBaselineOptimizedTaperratio0.210.21Lineartwistangle/(°)-2.49-4.00Angleofattack/(°)0.821.03Wingweight(single)/kg2030.31311.1Lift/t61345.854003.4(≥54)Lift-dragratio27.5627.01(≥27)Maxequivalentstress/(108Pa)2.73272.7282Maxdeformationatwingtip/m0.9470.998(≤1)

圖19 機(jī)翼氣動(dòng)/結(jié)構(gòu)綜合優(yōu)化過(guò)程中目標(biāo)函數(shù)收斂曲線(xiàn)Fig.19 Convergence history of objective function for wing aerostructural optimization

圖20 優(yōu)化機(jī)翼變形前后對(duì)比Fig.20 Comparison of optimum wings (un-deformed and deformed)

圖21 優(yōu)化前后機(jī)翼蒙皮厚度對(duì)比 Fig.21 Comparison of skin thickness of baseline and optimum wings

5 結(jié) 論

本文介紹了Kriging模型的基本理論和算法，包括梯度增強(qiáng)型Kriging、CoKriging和分層Kriging等研究新進(jìn)展。從一種通用優(yōu)化算法的角度，闡述了基于Kriging模型的代理優(yōu)化算法的原理和算法框架，介紹了優(yōu)化加點(diǎn)準(zhǔn)則和子優(yōu)化的研究進(jìn)展,討論了常用的幾種優(yōu)化加點(diǎn)準(zhǔn)則、約束處理方法和優(yōu)化終止條件。介紹了代理優(yōu)化算法標(biāo)準(zhǔn)算例驗(yàn)證和確認(rèn)的研究進(jìn)展。

1) 基于Kriging模型的代理優(yōu)化算法，能夠精確收斂到優(yōu)化問(wèn)題的真實(shí)最優(yōu)解，且具有精確的約束處理能力。因而可以作為一種通用優(yōu)化算法，應(yīng)用于具有光滑、連續(xù)設(shè)計(jì)空間的任意優(yōu)化問(wèn)題。

2) 代理優(yōu)化算法的優(yōu)化機(jī)制，是通過(guò)代理模型及其子優(yōu)化過(guò)程來(lái)產(chǎn)生新的樣本點(diǎn)。

3) 加點(diǎn)準(zhǔn)則及其子優(yōu)化問(wèn)題求解是基于Kriging模型的代理優(yōu)化算法的關(guān)鍵。MSP準(zhǔn)則偏向局部探測(cè)；EI、LCB準(zhǔn)則在全局探索和局部探測(cè)之間取得很好折中；PI和MSE準(zhǔn)則偏重于全局探索，收斂速度很慢，不適合單獨(dú)使用。將MSP與EI、LCB或MSE同時(shí)混合使用，成為更加實(shí)用的方法；已有的算例測(cè)試表明，“MSP+EI”組合加點(diǎn)準(zhǔn)則的綜合性能最佳。

4) 基于Kriging模型的代理優(yōu)化算法既適用于局部?jī)?yōu)化問(wèn)題，又適用于全局最優(yōu)解。對(duì)于全局優(yōu)化問(wèn)題，優(yōu)化效率比遺傳算法可提高1～2個(gè)量級(jí)。

5) 對(duì)于梯度可以快速獲得的問(wèn)題，引入梯度信息來(lái)建立梯度增強(qiáng)Kriging模型，不僅可大幅度提高優(yōu)化效率，而且能改善優(yōu)化精度。

6) 基于Kriging模型的代理優(yōu)化算法對(duì)于20維左右的工程設(shè)計(jì)問(wèn)題已經(jīng)基本成熟，在保證優(yōu)化質(zhì)量的同時(shí)，能實(shí)現(xiàn)很高的優(yōu)化效率。

基于Kriging模型的代理優(yōu)化算法存在如下關(guān)鍵問(wèn)題值得進(jìn)一步研究：

1) 維度之咒。對(duì)于30個(gè)以上設(shè)計(jì)變量的高維度優(yōu)化問(wèn)題，優(yōu)化效率較低。

2) 大規(guī)模樣本數(shù)據(jù)問(wèn)題。如果優(yōu)化過(guò)程中產(chǎn)生大量樣本點(diǎn)(2 000以上)，針對(duì)這些樣本點(diǎn)進(jìn)行代理模型訓(xùn)練的時(shí)間甚至超過(guò)精確數(shù)值分析的時(shí)間。

3) 對(duì)于梯度增強(qiáng)Kriging模型，隨著設(shè)計(jì)變量數(shù)增加，模型訓(xùn)練的時(shí)間大大增加。

4) 復(fù)雜設(shè)計(jì)空間問(wèn)題。主要針對(duì)設(shè)計(jì)空間呈超多極值、高度非線(xiàn)性、強(qiáng)烈振蕩等特性的復(fù)雜優(yōu)化問(wèn)題。

5) 數(shù)值噪聲問(wèn)題。樣本數(shù)據(jù)計(jì)算可能未充分收斂，這樣便帶有數(shù)值噪聲，從而對(duì)代理優(yōu)化的過(guò)程和結(jié)果產(chǎn)生重要影響。

經(jīng)過(guò)文獻(xiàn)調(diào)研，作者認(rèn)為今后在代理優(yōu)化方法研究方面還需要開(kāi)展的工作如下：

1) 通過(guò)同時(shí)引入Adjoint梯度和低可信度分析，發(fā)展適用于大規(guī)模設(shè)計(jì)變量問(wèn)題的模型訓(xùn)練及高效代理優(yōu)化算法。

2) 研究適用于大規(guī)模樣本數(shù)據(jù)的數(shù)據(jù)分區(qū)技術(shù)及相應(yīng)的代理優(yōu)化算法。

3) 研究局部代理模型和全局代理模型相結(jié)合的方法，克服代理模型剛性和局部收斂精度差的問(wèn)題。

4) 研究新的加點(diǎn)準(zhǔn)則及其子優(yōu)化方法，例如研究可同時(shí)加任意N個(gè)樣本點(diǎn)的算法，提高并行優(yōu)化性能。

5) 研究數(shù)據(jù)庫(kù)輔助的高效代理優(yōu)化算法。充分運(yùn)用分析和優(yōu)化工作過(guò)程產(chǎn)生的大量歷史數(shù)據(jù)，甚至包括風(fēng)洞試驗(yàn)得到的試驗(yàn)數(shù)據(jù)，對(duì)進(jìn)一步的優(yōu)化設(shè)計(jì)形成強(qiáng)有力的數(shù)據(jù)支持。

6) 研究代理優(yōu)化算法在復(fù)雜系統(tǒng)多學(xué)科優(yōu)化設(shè)計(jì)問(wèn)題的應(yīng)用。包括具有大規(guī)模設(shè)計(jì)變量和大規(guī)模約束，且各學(xué)科之間具有復(fù)雜耦合關(guān)系的工程設(shè)計(jì)問(wèn)題。

致 謝

感謝德國(guó)宇航中心(DLR)Goertz博士、布倫瑞克工業(yè)大學(xué)Zimmermann博士、凱澤斯勞滕大學(xué)Gauger教授。同時(shí)也感謝荷蘭代爾夫特工業(yè)大學(xué)Dwight教授、美國(guó)愛(ài)荷華州立大學(xué)Leifsson教授。他們與作者就本文工作進(jìn)行了大量討論并給予了許多建議。本文撰寫(xiě)歷時(shí)2年多，感謝課題組的老師和同學(xué)們?cè)谶@期間對(duì)本文撰寫(xiě)的大力協(xié)助和支持。

[1] HERNANDEZ S. Structural optimization: 1960—2010 and beyond[J]. Computational Technology Reviews, 2010, 2: 177-222.

[2] SCHMIT L A. Structural design by systematic synthesis[C]//Second Conference on Electronic Computation. Pittsburg: ASCE, 1960: 105-132.

[3] HICKS R M, MUNNAN E M, VANDERPLAATS G N. An assessment of airfoil design by numerical optimization: NASA TM X-3092[R]. Washington, D.C.: NASA, 1974.

[4] HICKS R M, HENNE P A. Wing design by numerical optimization[J]. Journal of Aircraft, 1978, 15(7): 407-412.

[5] SOBIESZCZANSKI-SOBIESKI J, HAFTKA R T. Multidisciplinary aerospace design optimization: Survey of recent developments[J]. Structural Optimization, 1997, 14(1): 1-23.

[6] KROO I, ALTUS S, BRAUN R, et al. Multidisciplinary optimization methods for aircraft preliminary design: AIAA-1994-4325[R]. Reston: AIAA, 1994.

[7] AIAA MDO Technical Committee. Current state of the art: Multidisciplinary design optimization: AIAA White Paper[R]. Reston: AIAA, 1991.

[8] ZHANG K S, HAN Z H, LI W J, et al. Bilevel adaptive weighted sum method for multidisciplinary multi-objective optimization[J]. AIAA Journal, 2008, 46(10): 2611-2622.

[9] 余雄慶. 飛機(jī)總體多學(xué)科設(shè)計(jì)優(yōu)化的現(xiàn)狀與發(fā)展方向[J]. 南京航空航天大學(xué)學(xué)報(bào), 2008, 40(4): 417-426.

YU X Q. Multidisciplinary design optimization for aircraft conceptual and preliminary design: Status and directions[J]. Journal of Nanjing University of Aeronautics & Astronautics, 2008, 40(4): 417-426 (in Chinese).

[10] 穆雪峰， 姚衛(wèi)星， 余雄慶, 等. 多學(xué)科設(shè)計(jì)優(yōu)化中常用代理模型的研究[J]. 計(jì)算力學(xué)學(xué)報(bào), 2005, 22(5): 608-612.

MU X F, YAO W X, YU X Q, et al. A survey of surrogate models used in MDO[J]. Chinese Journal of Computational Mechanics, 2005, 22(5): 608-612 (in Chinese).

[11] VIANA F A C, SIMPSON T W, BALABANOV V, et al. Metamodeling in multidisciplinary design optimization: how far have we really come?[J]. AIAA Journal, 2014, 52(4): 670-690.

[12] 張健, 李為吉. 飛機(jī)多學(xué)科設(shè)計(jì)優(yōu)化中的近似方法分析[J]. 航空計(jì)算技術(shù), 2005, 35(3): 5-8.

ZHANG J, LI W J. Approximation methods analysis in multidisciplinary design optimization[J]. Aeronautical Computer Technique, 2005, 35(3): 5-8 (in Chinese) .

[13] HAN Z H, GOERTZ S. A hierarchical Kriging model for variable-fidelity surrogate modeling[J]. AIAA Journal, 2012, 50(3): 1885-1896.

[14] HAN Z H, ZIMMERMANN R, GOERTZ S. An alternative CoKriging model for variable-fidelity surrogate modeling[J]. AIAA Journal, 2012, 50(5): 1205-1210.

[15] HAN Z H, ZHANG K S. Surrogate-based optimization[M]. InTech Book, 2012: 343-362.

[16] SCHMIT L A, FARSHI B. Some approximation concepts for structural synthesis[J]. AIAA Journal, 1974, 12(5): 692-699.

[17] GIUNTA A A, WATSON L T. A Comparison of approximation modeling techniques: Polynomial versus interpolation models: AIAA-1998-4758[R]. Reston: AIAA, 1998.

[18] KOCH P N, SIMPSON T W, ALLEN J K, et al. Statistical approximations for multidisciplinary design optimization: the problem of the size[J]. Journal of Aircraft, 1999, 36(1): 275-286.

[19] SEVANT N E, BLOOR M I G, WILSON M J. Aerodynamic design of a flying wing using response surface methodology[J]. Journal of Aircraft, 2000, 37(4): 562-569.

[20] JEONG S, MURAYAMA M, YAMAMOTO K. Efficient optimization design method using Kriging model[J]. Journal of Aircraft, 2005, 42(2): 413-420.

[21] VAVALLE A, QIN N. Iterative response surface based optimization scheme for transonic airfoil design[J]. Journal of Aircraft, 2007, 44(2): 365-376.

[22] KRIGE D G. A statistical approach to some basic mine valuations problems on the witwatersrand[J]. Journal of the Chemical, Metallurgical and Mining Engineering Society of South Africa, 1951, 52(6): 119-139.

[23] SACKS J, WELCH W J, MITCHELL T J, et al. Design and analysis of computer experiments[J]. Statistical Science, 1989, 4(4): 409-423.

[24] POWELL M J D. Algorithms for approximation[M]. New York: Oxford University Press, 1987: 141-167.

[25] BUHMANN M D. Acta numerica[M]. New York: Cambridge University Press, 2000: 1-38.

[26] KRISHNAMURTHY T. Response surface approximation with augmented and compactly supported radial basis functions: AIAA-2003-1748[R]. Reston: AIAA, 2003.

[27] MULLUR A A, MESSAC A. Extended radial basis functions: more flexible and effective metamodeling[J]. AIAA Journal, 2005, 43(6): 1306-1315.

[28] PARK J, SANDBERG I W. Universal approximation using radial-basis-function networks[J]. Neural Computation, 1991, 3(2): 246-257.

[29] ELANAYAR S V T, SHIN Y C. Radial basis function neural network for approximation and estimation of nonlinear stochastic dynamic systems[J]. IEEE Transactions on Neural Networks, 1994, 5(4): 594-603.

[30] SMOLA A J, SCH?LKOPF B. A tutorial on support vector regression[J]. Statistics and Computing, 2004, 14(3): 199-222.

[31] FORRESTER A I J, KEANE A J. Recent advances in surrogate-based optimization[J]. Progress in Aerospace Sciences, 2009, 45(1): 50-79.

[32] WANG Q, MOIN P, IACCARINO G. A rational interpolation scheme with super-polynomial rate of convergence[J]. SIAM Journal of Numerical Analysis, 2010, 47(6): 4073-4097.

[33] WANG Q, MOIN P, IACCARINO G. A high-order multi-variate approximation scheme for arbitrary data sets[J]. Journal of Computational Physics, 2010, 229(18): 6343-6361.

[34] WIENER N. The homogeneous chaos[J]. American Journal of Mathematics, 1938, 60(4): 897-936.

[35] XIU D. Numerical methods for stochastic computations: A spectral method approach[M]. Princeton: Princeton University Press, 2010: 152.

[36] QUEIPO N V, HAFTKA R T, SHYY W, et al. Surrogate-based analysis and optimization[J]. Progress in Aerospace Sciences, 2005, 41(1): 1-28.

[37] FORRESTER A I J, SBESTER A, KEANE A. Engineering design via surrogate modeling: A practical guide[M]. Chichester: John Wiley & Sons, 2008.

[38] KEANE A J, NAIR P B. Computational approaches for aerospace design: The pursuit of excellence[M]. Chichester: John Wiley & Sons, 2005.

[39] GIUNTA A A, WOJTKIEWICZ J S F, ELDRED M S. Overview of modern design of experiments methods for computational simulations: AIAA-2003-649[R]. Reston: AIAA, 2003.

[40] FANG K T, LIN D, WINKER P, et al. Uniform design: Theory and application[J]. Technometrics, 2000, 42(3): 237-248.

[41] MATHERON G M. Principles of geostatistics[J]. Economic Geology, 1963, 58(8): 1246-1266.

[42] MATHERON G. Theory of regionalized variable and its applications[M]. Fontainebleau: Ecole des Mines, 1971.

[43] SIMPSON T W, MAUERY T M, KORTE J J, et al. Kriging models for global approximation in simulation-based multidisciplinary design optimization[J]. AIAA Journal, 2001, 39(12): 2233-2241.

[44] MARTIN J D, SIMPSON T W. Use of Kriging models to approximate deterministic computer models[J]. AIAA Journal, 2005, 43(4): 853-863.

[45] TOAL D J J, BRESSLOFF N W, KEAN A J. Kriging hyperparameter tuning strategies[J]. AIAA Journal, 2008, 46(5): 1240-1252.

[46] 王曉鋒, 席光. 基于Kriging模型的翼型氣動(dòng)性能優(yōu)化設(shè)計(jì)[J]. 航空學(xué)報(bào), 2005, 26(5): 545-549.

WANG X F, XI G. Aerodynamic optimization design for airfoil based on Kriging model[J]. Acta Aeronautica et Astronautica Sinica, 2005, 26(5): 545-549 (in Chinese).

[47] 許瑞飛， 宋文萍， 韓忠華. 改進(jìn)Kriging模型在翼型氣動(dòng)優(yōu)化設(shè)計(jì)中的應(yīng)用研究[J]. 西北工業(yè)大學(xué)學(xué)報(bào), 2010， 28(4)： 503-510.

XU R F, SONG W P, HAN Z H. Application of improved kriging-model-based optimization method in airfoil aerodynamic design[J]. Journal of Northwestern Polytechnical University, 2010， 28(4)： 503-510 (in Chinese).

[48] LIU J, HAN Z H, SONG W P. Efficient kriging-based optimization design of transonic airfoils: Some key issues: AIAA-2012-0967[R]. Reston: AIAA, 2012.

[49] 孫俊峰, 劉剛, 江雄, 等. 基于Kriging模型的旋翼翼型優(yōu)化設(shè)計(jì)研究[J]. 空氣動(dòng)力學(xué)學(xué)報(bào), 2013, 31(4): 437-441.

SUN J F, LIU G, JIANG X, et al. Research of rotor airfoil design optimization based on the Kriging model[J]. Acta Aerodynamica Sinica, 2013, 31(4): 437-441 (in Chinese).

[50] HAN Z H, LIU J, SONG W P, et al. Surrogate-based aerodynamic shape optimization with application to wind turbine airfoils: AIAA-2013-1108[R]. Reston: AIAA, 2013.

[51] 白俊強(qiáng), 王波, 孫智偉, 等. 基于松散式代理模型管理框架的亞音速機(jī)翼優(yōu)化設(shè)計(jì)方法研究[J]. 西北工業(yè)大學(xué)學(xué)報(bào), 2011, 29(4): 515-519.

BAI J Q, WANG B, SUN Z W, et al. Developing optimization design of subsonic wing with losse type of agent model[J]. Journal of Northwestern Polytechnical University, 2011, 29(4): 515-519 (in Chinese).

[52] 孫美建, 詹浩. Kriging模型在機(jī)翼氣動(dòng)外形優(yōu)化中的應(yīng)用[J]. 空氣動(dòng)力學(xué)學(xué)報(bào), 2011, 29(6): 759-764.

SUN M J, ZHAN H. Application of Kriging surrogate model for aerodynamic shape optimization of wing[J]. Acta Aerodynamica Sinica, 2011, 29(6): 759-764 (in Chinese).

[53] 何歡, 朱廣榮, 何成, 等. 基于Kriging模型的結(jié)構(gòu)耐撞性?xún)?yōu)化[J]. 南京航空航天大學(xué)學(xué)報(bào), 2014, 46(2): 297-303.

HE H, ZHU G R, HE C, et al. Crashworthiness optimization based on Kriging metamodeling[J]. Journal of Nanjing University of Aeronautics & Astronautics, 2014, 46(2): 297-303 (in Chinese).

[54] 劉克龍, 姚衛(wèi)星, 穆雪峰. 基于Kriging代理模型的結(jié)構(gòu)形狀優(yōu)化方法研究[J]. 計(jì)算力學(xué)學(xué)報(bào), 2006, 23(3): 344-347, 362.

LIU K L, YAO W X, MU X F. A method of structural shape optimization based on Kriging model[J]. Chinese Journal of Computational Mechanics, 2006, 23(3): 344-347, 362 (in Chinese).

[55] 楊軍, 劉勇瓊, 艾春安, 等. 改進(jìn)Kriging模型在固沖發(fā)動(dòng)機(jī)導(dǎo)彈氣動(dòng)優(yōu)化設(shè)計(jì)中的應(yīng)用[J]. 固體火箭技術(shù), 2013, 36(5): 590-593.

YANG J, LIU Y Q, AI C A, et al. Application of improved Kriging-model-based optimization method in solid rocket-ramjet missile’s aerodynamic design[J]. Journal of Solid Rocket Technology, 2013, 36(5): 590-593 (in Chinese).

[56] 孫凱軍, 宋文萍, 韓忠華. 基于Kriging模型的高超聲速舵面優(yōu)化設(shè)計(jì)[J]. 航空計(jì)算技術(shù), 2012, 42(2): 9-12.

SUN K J, SONG W P, HAN Z H. Optimization design of hypersonic control surface based on Kriging model[J]. Aeronautical Computing Technique, 2012, 42(2): 9-12 (in Chinese).

[57] 鄭侃, 廖文和, 張翔. 基于近似模型管理的微小衛(wèi)星結(jié)構(gòu)多目標(biāo)優(yōu)化設(shè)計(jì)[J]. 中國(guó)機(jī)械工程, 2012, 23(6): 655-659.

ZHENG K, LIAO W H, ZHANG X. Multi-objective optimization design for microsatellite structure based on approximation model management[J]. China Mechanical Engineering, 2012, 23(6): 655-659 (in Chinese).

[58] 姚拴寶， 郭迪龍， 孫振旭， 等. 基于Kriging代理模型的高速列車(chē)頭型多目標(biāo)優(yōu)化設(shè)計(jì)[J]. 中國(guó)科學(xué), 2013, 43(2): 186-200.

YAO S B, GUO D L, SUN Z X, et al. Multi-objective optimization of the streamlined head of high-speed trains based on the kriging model[J]. Science China, 2013, 43(2): 186-200 (in Chinese).

[59] 韓永志, 高行山, 李立州. 基于Kriging模型的渦輪葉片多學(xué)科設(shè)計(jì)優(yōu)化[J]. 航空動(dòng)力學(xué)報(bào), 2007, 22(7): 1055-1059.

HAN Y Z， GAO X S, LI L Z. Kriging model-based multidisciplinary design optimization for turbine blade[J]. Journal of Aerospace Power, 2007, 22(7): 1055-1059 (in Chinese).

[60] 張科施, 韓忠華, 李為吉, 等. 基于近似技術(shù)的高亞聲速運(yùn)輸機(jī)機(jī)翼氣動(dòng)/結(jié)構(gòu)優(yōu)化設(shè)計(jì)[J]. 航空學(xué)報(bào), 2006, 27(5): 810-815.

ZHANG K S, HAN Z H, LI W J, et al. Multidisciplinary aerodynamic/structural design optimization for high subsonic transport wing using approximation technique[J]. Acta Aeronautica et Astronautica Sinica, 2006, 27(5): 810-815 (in Chinese).

[61] SUN W Y, YUAN Y X. Optimization theory and methods: nonlinear programming[M]. New York: Springer Ebooks, 2006.

[62] HOLLAND J H. Adaptation in natural and artificial systems[M]. Ann Arbor: The University of Michigan Press, 1975.

[63] SLOTNICK J, KHODADOUST A, ALONSO J, et al. CFD vision 2030 study: A path to revolutionary computational aerosciences: NASA/CR-2014-218178[R]. Washington, D.C.: NASA Langley Research Center, 2014.

[64] 吳寬展， 劉學(xué)軍， 呂宏強(qiáng). 超臨界翼型設(shè)計(jì)中的多響應(yīng)代理模型[J]. 航空計(jì)算技術(shù), 2014, 44(4): 17-22.

WU K Z, LIU X J, LV H Q. Multi-output surrogate model in supercritical airfoil design[J]. Aeronautical Computing Technique, 2014, 44(4): 17-22 (in Chinese).

[65] ZIMMERMANN R. On the condition number anomaly of gaussian correlation matrices[J]. Linear Algebra and its Applications, 2015, 466: 521-526.

[66] VIANA F A, HAFTKA R T, STEFFEN V. Multiple surrogates: How cross-validation errors can help us to obtain the best predictor?[J]. Structural and Multidisciplinary Optimization, 2009, 39(4): 439-457.

[67] KOBLEWALIK J, OSBORNE M R. Methods for unconstrained optimization problems[M]. New York: Elsevier, 1968.

[68] NELDER J A, MEAD R. A simple method for function minimization[J]. The Computer Journal, 1965, 7(4): 308-313.

[69] LOPHAVEN S N, NIELSEN H B, S?NDERGAARD J. DACE—A MATLAB Kriging toolbox (version 2.0): IMM-REP-2002-12[R]. Lyngby: Informatics and Mathematical Modelling, Technical University of Denmark, 2002.

[70] YIN J, NG S H, NG K M. Kriging meta model with modified nugget-effect: The heteroscedastic variance case[J]. Computers Industrial Engineering, 2011, 61(3): 760-777.

[71] FORRESTER A I J, KEANE A J, BRESSLOFF N W. Design and analysis of noisy computer experiments[J]. AIAA Journal, 2006, 44(10): 2331-2339.

[72] ETMAN L F P. Design and analysis of computer experiments: The Method of Sacks et al： Engineering Mechanics report WFW 94-098[R]. Eindhoven, The Netherlands: Eindhoven University of Technology, 1994.

[73] OSIO I G, AMON C H. An engineering design methodology with multistage bayesian surrogate and optimal sampling[J]. Research in Engineering Design, 1996, 8(4): 189-206.

[74] JONES D R, SCHONLAU M, WELCH W J. Efficient global optimization of expensive black-box functions[J]. Journal of Global Optimization, 1998, 13(4): 455-492.

[75] COX D D, JOHN S. SDO: A statistical method for global optimization[C]//IEEE International Conference on Systems, Man & Cybernetics, 1992: 1241-1246.

[76] TORCZON V, TROSSET M W. Use approximation to accelerate engineering design optimization: AIAA-1998-4800[R]. Reston: AIAA, 1998.

[77] BOOKER A J. Design and analysis of computer experiments: AIAA-1998-4757[R]. Reston: AIAA, 1998.

[78] MECKESHEIMER M, BARTON R R, SIMPSON T, et al. Metamodeling of combined discrete/continuous responses[J]. AIAA Journal, 2001, 39(10): 1950-1959.

[79] CHEN S K, XIONG Y, CHE W. Multiresponse and multistage metamodeling approach for design optimization[J]. AIAA Journal, 2009, 47(1): 206-218.

[80] HAN Z H, ZHANG K S, SONG W P, et al. Optimization of active flow control over an airfoil using a surrogate-management framework[J]. Journal of Aircraft, 2010, 47(2): 603-612.

[81] ZHANG K S, HAN Z H, LI W J, et al. Coupled aerodynamic/structural optimization of a subsonic transport wing using a surrogate model[J]. Journal of Aircraft, 2008, 45(6): 2167-2171.

[82] CHUNG J J, ALONSO H S. Using gradients to construct cokriging approximation models for high-dimensional design optimization problems： AIAA-2002-0317[R]. Reston: AIAA, 2002.

[83] CHUNG J J, ALONSO H S. Design of a low-boom supersonic business jet using cokriging approximation models: AIAA-2002-5598[R]. Reston: AIAA, 2002.

[84] LIU W, BATILL S M. Gradient-enhanced response surface approximations using kriging models: AIAA-2002-5456[R]. Reston: AIAA, 2002.

[85] LAURENCEAU J, SAGAUT P. Building efficient response surfaces of aerodynamic functions with Kriging and Cokriging[J]. AIAA Journal, 2008, 46(2): 498-507.

[86] JAMESON A. Optimum aerodynamic design using cfd and control theory: AIAA-1995-1729[R]. Reston: AIAA, 1995.

[87] DWIGHT R, HAN Z H. Efficient uncertainty quantification using gradient enhanced Kriging: AIAA-2009-2276[R]. Reston: AIAA, 2009.

[88] HAN Z H, GOERTZ 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.

[89] YAMAZAKI W, RUMPFKEIL M P, MAVRIPLIS D J. Design optimization utilizing gradient/hessian enhanced surrogate model: AIAA-2010-4363[R]. Reston: AIAA, 2010.

[90] HAN Z H. Improving adjoint-based aerodynamic optimization via gradient-enhanced Kriging: AIAA-2012-0670[R]. Reston: AIAA, 2012.

[91] DAVID M. Geostatistical ore reserve estimation[M]. Amsterdam: Elsevier, 1977: 1-364.

[92] JOURNEL A G, HUIJBREGTS J C. Mining geostatistics[M]. New York: Academic Press, 1978: 1-600.

[93] MYERS D E. Matrix formulation of Cokriging[J]. Mathematical Geology, 1982, 14(3): 249-257.

[94] GOOVAERTS P. Ordinary CoKriging revisited[J]. Mathematical Geology, 1997, 30(1): 21-42.

[95] PAPRITZ A. Standardized vs customary ordinary CoKriging: some comments on the article “the geostatistical analysis of experiments at the landscape-scale” by T.F.A. Bishop and R.M., Lark[J]. Geoderma, 2008, 146(1): 391-396.

[96] KENNEDY M C, O’HAGAN A. Predicting the output from a complex computer code when fast approximations are available[J]. Biometrika, 2000, 87(1): 1-13.

[97] QIAN Z, WU C F J. Bayesian hierarchical modeling for integrating low-accuracy and high-accuracy experiments[J]. Technometrics, 2008, 50(2):192-204.

[98] FORRESTER A I J, SBESTER A, KEANE A J. Multi-fidelity optimization via surrogate modeling[J]. Proceedings of the Royal Society A, Mathematical, Physical and Engineering Sciences, 2007, 463(2088): 3251-3269.

[99] KUYA Y, TAKEDA K, ZHANG X, et al. Multifidelity surrogate modeling of experimental and computational aerodynamic data sets[J]. AIAA Journal, 49(2): 289-298.

[100] ZIMMERMANN R, HAN Z H. Simplified cross-correlation estimation for multi-fidelity surrogate cokriging models[J]. Advance and Applications in Mathematical Sciences, 2010, 7(2): 181-201.

[101] HAN Z H, ZIMMERMANN R, GOERTZ S. A new CoKriging method for variable-fidelity surrogate modeling of aerodynamic data: AIAA-2010-1225[R]. Reston: AIAA, 2010.

[102] JOSEPH V R, HUNG Y, SUDJIANTO A. Blind Kriging: A new method for developing metamodels[J]. Journal of Mechanical Design, 2008, 130(3): 350-353.

[103] ZHAO L, CHOI K K, LE I. Metamodeling method using dynamic Kriging for design optimization[J]. AIAA Journal, 2011, 49(9): 2034-2046.

[104] LIANG H Q, ZHU M, WU Z. Using cross-validation to design trend function in Kriging surrogate modeling[J]. AIAA Journal, 2014, 52(10): 2313-2327.

[105] ANKENMAN B E, NELSON B L, STAUM J. Stochastic Kriging for simulation metamodeling[J]. Operations Research, 2010, 58(2): 371-382.

[106] LIU J, HAN Z H, SONG W P. Comparison of infill sampling criteria in kriging-based aerodynamic optimization[C]//28th Congress of the International Council of the Aeronautical Sciences, 2012.

[107] JONES D R. A taxonomy of global optimization methods based on response surfaces[J]. Journal of Global Optimization，2001, 21(4): 345-383.

[108] BOOKER A J, DENNIS J J, FRANK P D, et al. A rigorous framework for optimization of expensive functions by surrogates[J]. Structural Optimization, 1998, 17(1): 1-13.

[109] SASENA M J, PAPALAMBROS P Y, GOOVAERTS P. Exploration of metamodeling sampling criteria for constrained global optimization[J]. Engineering Optimization, 2002, 34(34): 263-278.

[110] PARR J M, KEANE J A, FORRESTER I J, et al. Infill sampling criteria for surrogate-based optimization with constraint handling[J]. Engineering Optimization , 2012, 44(10): 1147-1166.

[111] DEB K. An efficient constraint handling method for genetic algorithms[J]. Computer Methods in Applied Mechanics and Engineering, 2010, 186(2-4): 311-338

[112] ONG Y S, NAIR P B, KEANE A J. Evolutionary optimization of computationally expensive problems via surrogate modeling[J]. AIAA Journal, 2003, 41(4): 687-696.

[113] 劉俊. 基于代理模型的高效氣動(dòng)優(yōu)化設(shè)計(jì)方法及應(yīng)用[D]. 西安： 西北工業(yè)大學(xué), 2015.

LIU J. Efficient surrogate-based optimization method and its application in aerodynamic design[D]. Xi’an: Northwestern Polytechnical University, 2015 (in Chinese).

[114] 高月華， 王希誠(chéng). 基于Kriging代理模型的多點(diǎn)加點(diǎn)序列優(yōu)化方法[J]. 工程力學(xué)， 2012， 29(4): 90-95.

GAO Y H, WANG X C. A sequential optimization method with multi-point sampling criterion based on kriging surrogate mode[J]. Engineering Mechanics, 2012, 29 (4): 90-95 (in Chinese).

[115] 劉俊, 宋文萍, 韓忠華, 等. 梯度增強(qiáng)的Kriging模型與Kriging模型在優(yōu)化設(shè)計(jì)中的比較研究[J]. 西北工業(yè)大學(xué)學(xué)報(bào), 2015, 3(3): 819-826.

LIU J, SONG W P, HAN Z H, et al. Comparative study of gek (gradient-enhanced Kriging) and Kriging when applied to design optimization[J]. Journal of Northwestern Polytechnical University, 2015, 3(3): 819-826 (in Chinese).

[116] JAMIL M, YANG X S. A literature survey of benchmark functions for global optimization problems[J]. Journal of Mathematical Modelling and Numerical Optimisation, 2013, 4(2): 150-194.

[117] ROSENBROCK H H. An automatic method for finding the greatest or least value of a function[J]. The Computer Journal, 1960, 3(3): 175-184.

[118] MICHALEWICZ Z. Genetic algorithm, numerical optimization, and constraints[C]//Proceedings of the Sixth International Conference on Genetic Algorithms, 1995: 151-158.

[119] TESFAHUNEGN Y A, KOZIEL S, GRAMANZINI J R, et al. Application of direct and surrogate-based optimization to two-dimensional benchmark aerodynamic problems: A comparative study: AIAA-2015-0265[R]. Reston: AIAA, 2015.

[120] ZHANG Y, HAN Z H, SHI L X, et al. Multi-round surrogate-based optimization for benchmark aerodynamic design problems: AIAA-2016-1545[R]. Reston: AIAA, 2016.

[121] REN J, THELEN A, AMRIT A, et al. Application of multi-fidelity optimization techniques to benchmark aerodynamic design problems: AIAA-2016-1542[R]. Reston: AIAA, 2016.

[122] LYU Z, KENWAY G KW, MARTINS J R R A. Aerodynamic shape optimization studies on the common research model wing benchmark[J]. AIAA Journal, 2015, 53(4): 968-985.

[123] KENWAY G K W, MARTINS J R R A. Multipoint aerodynamic shape optimization investigations of the common research model wing[J]. AIAA Journal, 2016, 54(1): 112-128.

[124] KENWAY G K W, MARTINS J R R A. Aerodynamic shape optimization of the CRM configuration including buffet-onset conditions: AIAA-2016-1294[R]. Reston: AIAA, 2016.

[125] LIEM R, KENWAY G K W, MARTINS J R R A. Multi-mission aircraft fuel burn minimization via multipoint aerostructural optimization[J]. AIAA Journal, 2015, 53(1): 104-122.

[126] 張科施， 韓忠華， 李為吉， 等. 一種考慮氣動(dòng)彈性的運(yùn)輸機(jī)機(jī)翼多學(xué)科優(yōu)化方法[J]. 空氣動(dòng)力學(xué)學(xué)報(bào)， 2008， 26(1): 1-7.

ZHANG K S, HAN Z H, LI W J, et al. A method of coupled aerodynamic/structural integration optimization for transport-wing design[J]. Acta Aerodynamica Sinica, 2008, 26(1):1-7 (in Chinese).

[127] KENWAY G K W, KENNEDY G J, MARTINS J R R A. Scalable parallel approach for high-fidelity steady-state aeroelastic analysis and derivative computations[J]. AIAA Journal, 2014, 52(5): 935-951.

韓忠華男， 博士， 教授， 博士生導(dǎo)師。主要研究方向： 代理模型理論與算法， 氣動(dòng)與多學(xué)科優(yōu)化設(shè)計(jì)， 轉(zhuǎn)捩預(yù)測(cè)與自然層流翼型/機(jī)翼設(shè)計(jì), 氣動(dòng)數(shù)據(jù)庫(kù)技術(shù)。

Tel.: 029-88492704

E-mail: hanzh@nwpu.edu.cn

*Correspondingauthor.Tel.：029-88492704E-mail:hanzh@nwpu.edu.cn

Krigingsurrogatemodelanditsapplicationtodesignoptimization:Areviewofrecentprogress

HANZhonghua*

NationalKeylaboratoryofScienceandTechnologyonAerodynamicDesignandResearch,SchoolofAeronautics,NorthwesternPolytechnicalUniversity,Xi’an710072,China

Overthepasttwodecades,surrogatemodelinghas

muchattentionfromtheresearchersintheareaofaerospacescienceandengineeringduetoitscapabilityofgreatlyimprovingtheefficiencyofdesignoptimizationwhenhigh-fidelitynumericalanalysisisemployed.Designoptimizationviasurrogatemodelsisintensivelyresearchedandeventuallyleadstoanewtypeofoptimizationalgorithmwhichiscalledsurrogate-basedoptimization(SBO).Amongtheavailablesurrogatemodels,suchaspolynomialresponsesurfacemodel,radial-basisfunctions,artificialneutralnetwork,support-vectorregression,multivariateinterpolationorregression,andpolynomialchaosexpansion,Krigingmodelisthemostrepresentativesurrogatemodelwhichhasgreatpotentialinengineeringdesignandoptimization.Inthecontextofaircraftdesign,thispaperreviewsthetheory,algorithmandrecentprogressforresearchesontheKrigingsurrogatemodel.First,thefundamentaltheoryandalgorithmofKrigingmodelarebrieflyreviewedandtheexperienceabouthowtoimprovetherobustnessandefficiencyispresented.Second,threemajorbreakthroughsofKrigingmodelinrecentyearsarereviewed,includinggradient-enhancedKriging,CoKrigingandhierarchicalKriging.Third,theoptimizationmechanismandframeworkofsurrogate-basedoptimizationusingKrigingmodelarediscussed.Inthemeanwhile,theconceptofinfill-samplingcriterionandsuboptimizationispresented.Fiveinfill-samplingcriteriaaswellasthededicatedconstrainthandlingmethodsaredescribed.Furthermore,thenewlydevelopedlocalEI(expectedimprovement)methodandterminationcriteriaforSBOareintroduced.Fourth,anumberoftestcasesincludingbenchmarkoptimizationproblemsaswellasaerodynamicandmultidisciplinarydesignoptimizationproblemsaregiventodemonstratetheexcellentperformanceandgreatpotentialofthesurrogate-basedoptimizationusingKrigingmodel.Atlast,thekeychallengesaswellasfuturedirectionsaboutthetheory,algorithmandapplicationsarediscussed.

optimizationmethod;Kriging;surrogatemodel;aircraftdesign;multidisciplinarydesignoptimization(MDO)

2016-01-05；Revised2016-01-27；Accepted2016-03-15；Publishedonline2016-03-291455

URL：www.cnki.net/kcms/detail/11.1929.V.20160329.1455.004.html

NationalNaturalScienceFoundationofChina(11272265)

2016-01-05；退修日期2016-01-27；錄用日期2016-03-15； < class="emphasis_bold">網(wǎng)絡(luò)出版時(shí)間

時(shí)間：2016-03-291455

www.cnki.net/kcms/detail/11.1929.V.20160329.1455.004.html

國(guó)家自然科學(xué)基金 (11272265)

*

.Tel.：029-88492704E-mailhanzh@nwpu.edu.cn

韓忠華．Kriging模型及代理優(yōu)化算法研究進(jìn)展J． 航空學(xué)報(bào),2016,37(11):3197-3225.HANZH.KrigingsurrogatemodelanditsapplicationtodesignoptimizationareviewofrecentprogressJ.ActaAeronauticaetAstronauticaSinica,2016,37(11):3197-3225.

http://hkxb.buaa.edu.cnhkxb@buaa.edu.cn

10.7527/S1000-6893.2016.0083

V211.3

A

1000-6893(2016)11-3197-29