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

韓忠華

西北工業(yè)大學 航空學院 翼型葉柵空氣動力學國家級重點實驗室， 西安 710072

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

韓忠華*

西北工業(yè)大學 航空學院 翼型葉柵空氣動力學國家級重點實驗室， 西安 710072

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

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

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

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

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

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

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

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

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

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

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

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

1.1 代理模型問題的基本描述

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

(1)

對這n個設計方案進行數值分析(CFD或CSM)，從分析結果中得到相應的n個函數響應值：

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

(2)

1.2 Kriging模型及其預估值

Kriging模型是一種插值模型，其插值結果定義為已知樣本函數響應值的線性加權，即

(3)

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

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

(4)

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

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

(5)

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

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

(6)

最小，并滿足如下插值條件(或無偏差條件)：

(7)

采用拉格朗日乘數法，經過推導可證明最優(yōu)加權系數ω由式(8)描述的線性方程組(也稱為Kriging模型方程組)給出

(8)

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

(9)

式中：

(10)

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

(11)

通過分塊矩陣求逆，該模型可最終寫為

(12)

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

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

(13)

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

1.3 相關函數的選擇

在式(12) 所示的Kriging模型中，R和r的構造均涉及相關函數的選擇和計算。目前比較流行的一類相關函數為

(14)

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

(15)

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

實踐證明，在實際應用中一些其他的相關函數可近似滿足高斯假設，也具有很好的性能。例如，目前比較流行一類三次樣條函數：

(16)

該函數二階可導，在光滑性和魯棒性方面取得很好的折中。

1.4 模型參數訓練(或優(yōu)化)

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

(17)

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

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

(18)

將式(18)代入式(17)后，取對數，并忽略常數項，優(yōu)化問題轉化為使式(19)最大：

(19)

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

1.5 Kriging模型具體實現的幾個關鍵性問題

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

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

(20)

其次，必須合理地進行模型參數的優(yōu)化。合理的模型參數，可顯著提高Kriging模型近似精度，使之明顯優(yōu)于多項式響應面、徑向基函數等代理模型。但是，不恰當的模型參數值，可能使Kriging模型預測精度差，甚至建立失敗。

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

R′=R+γI

(21)

γ=(103+n)εM

(22)

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

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

R=LLT

(23)

(24)

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

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

FLOPS≈(1/3)n3+2n2

Memory size≈8×10-6n2MB

(25)

2 Kriging模型理論和算法研究新進展

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

這里結合本文的研究工作，重點介紹21世紀以來，在Kriging模型的研究和應用過程中，出現的具有代表性的3個研究新進展。

2.1 梯度增強型Kriging模型

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

實現引入梯度信息的最簡單做法是利用一階泰勒展開：

(26)

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

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

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

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

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

(27)

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

GEK模型對未知函數的預估值定義為所有抽樣函數值和偏導數值的線性加權，即

(28)

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

(29)

與傳統的Kriging模型類似，可推導出GEK模型對未知函數的預估值表達式為

(30)

(31)

(32)

其中：相關矩陣R、一階導數矩陣?R以及二階導數矩陣?2R定義為

(33)

相關矢量r及其一階導數?r定義為

(34)

GEK模型對預估值的均方差估計為

(35)

在構造GEK模型過程中，模型參數訓練方法見文獻[88，90]。

需要說明的是，GEK模型相關矩陣中還涉及相關函數的一階和二階導數。對于三次樣條函數，可推導出其一階導數為

(36)

(37)

且Rk和ξl的計算參見式(16);sign()為取符號運算。

同理，可推導出其二階導數表達式為

(38)

(39)

2.2 CoKriging模型

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

2010年，文獻[101]獨立提出了一種可用于建立變可信度代理模型的實用CoKriging模型(見文獻[14])。下面對該模型進行介紹。

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

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

(40)

式中：下標“1”和“2”分別代表高、低可信度，例如n1和n2分別代表高、低可信度樣本點數(假設n2?n1)。相應的目標函數或約束函數的響應值為

(41)

CoKriging模型預估值定義為

(42)

式中：λ1、λ2分別為對高、低可信度響應值的加權系數。假設存在分別與y1、y2對應的2個靜態(tài)隨機過程：

(43)

則設計空間不同位置處，隨機變量之間的協方差和交叉協方差定義為

(44)

(45)

(46)

且有

(47)

Cokriging模型預估值的均方差為

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

(48)

(49)

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

(50)

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

由于可合理假定高、低可信度數值分析結果具有相似的函數變化趨勢，為簡化起見，可假定

(51)

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

2.3 分層 Kriging 模型

首先，假設y1對應的隨機過程定義為

(52)

式中：β0為待定的自適應因子，它反映了高、低可信度函數之間的相關性；如果β0=0，說明高、低可信度函數之間沒有相關性。β0的引入，有效避免了引入低可信度分析后代理模型精度反而變差的風險。采用與Kriging模型類似的推導方法(詳細推導過程見文獻[13])，可得到HK模型的預估值為

(53)

式中：

(54)

預估值的均方差為

(55)

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

2.4 有關Kriging模型的其他研究新進展

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

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

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

Miny(x)

(56)

式中：y(x)、gi(x)為目標函數與約束函數；NC為約束個數；xu、xl為設計變量的上、下限。

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

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

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

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

圖2 傳統的梯度或無梯度優(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

為了更好地認識代理優(yōu)化算法，對圖3進行重新整理，將“建立代理模型和根據加點準則的子優(yōu)化”過程看作一個獨立模塊，得到圖4。比較圖4 和圖2，不難發(fā)現，除優(yōu)化器(Optimizer)模塊外，二者完全相同。這就是說，“建立代理模型和根據加點準則的子優(yōu)化”是代理優(yōu)化算法的優(yōu)化機制。

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

3.2 試驗設計與初始樣本選擇

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

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

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

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

3.3 優(yōu)化加點準則及其子優(yōu)化求解

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

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

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

(57)

3.3.1 最小化代理模型預測準則

最小化代理模型預測準則是最簡單、最直接，也是最早被采用的方法[15，106-110]。其原理是直接在代理模型上尋找目標函數的最小值。帶約束的子優(yōu)化問題數學模型為

(58)

圖6 一維問題MSP加點法則原理圖Fig.6 Schematic of MSP infill-sampling criterion for a 1D problem

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

(59)

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

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

(60)

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

(61)

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

3.3.2 改善期望準則

(62)

對于最小化問題，目標函數改善量I(x)定義為

(63)

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

E[I(x)]=

(64)

圖7 某一維問題EI函數分布示意圖 Fig.7 Schematic of EI function distribution for a 1D problem

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

MaxE[I(x)]

s.t.xl≤x≤xu

(65)

從而可得到新的樣本點x*。

(66)

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

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

(67)

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

s.t.xl≤x≤xu

(68)

從而得到新的樣本點x*。

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

3.3.3 改善概率準則

改善概率準則方法[31]實際上是EI方法派生的方法。與EI準則類似，該準則通過計算目標函數改善的概率，并將改善概率最大點作為新的樣本點。目標函數改善的概率為[106-107，109]

(69)

對于上述一維問題，PI函數分布如圖8所示。

圖8 某一維問題的PI函數分布示意圖 Fig.8 Schematic of PI function distribution for a 1D problem

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

s.t.xl≤x≤xu

(70)

得到新的樣本點x*。

3.3.4 均方差準則

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

s.t.xl≤x≤xu

(71)

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

s.t.xl≤x≤xu

(72)

得到新的樣本點x*。

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

3.3.5 置信下界準則

(73)

圖10 某一維問題的LCB函數分布示意圖Fig.10 Schematic of LCB function distribution for a 1D problem

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

(74)

從而得到新的樣本點x*。

文獻[106-107，113]均對上述優(yōu)化加點準則進行了比較，得到了一些具有指導性的建議。這里再補充幾點說明：

1) 除MSP準則適用于多項式響應面、徑向基函數等其他代理模型外，EI、PI、MSE、LCB準則一般只適用于Kriging模型以及其他具有誤差估計功能的代理模型。

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

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

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

5) 既可用于全局優(yōu)化，也可用于局部優(yōu)化問題。對于局部優(yōu)化問題，為了提高收斂性，可采用信賴域方法[107]和采用梯度增強型Kriging[113]。

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

優(yōu)化終止條件是優(yōu)化算法必須要具備的一個因素。目前國內外還沒有關于定義代理優(yōu)化算法優(yōu)化終止條件的文獻報道。于是，本文作者最近根據代理優(yōu)化的特點，并借鑒其他優(yōu)化算法的終止條件，定義了如下4種終止條件。

首先，可以根據樣本點之間的距離和目標函數值的差別進行定義：

(75)

(76)

其次,可根據代理模型在最優(yōu)點附近的近似精度來定義：

(77)

再次，可根據計算資源情況，設定調用精確數值模擬分析的最大次數：

N≤Nmax

(78)

例如針對氣動優(yōu)化設計問題，由于計算資源有限，一般希望在100～300次CFD計算以內，獲得最優(yōu)的外形，因此可以設定Nmax=100～300。

最后，可針對特定的加點準則進行定義。例如針對EI準則，可定義為

max(EI)≤EImin

(79)

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

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

代理模型最初作為一種近似方法，用于提高采用高可信度數值分析的工程設計問題的效率[76]。因此，人們往往針對特定的問題發(fā)展特定的代理優(yōu)化框架，例如文獻[20-21]中發(fā)展的翼型優(yōu)化設計方法，文獻[52，85]的機翼優(yōu)化設計方法，以及文獻[54]發(fā)展的結構優(yōu)化設計方法。但是，代理優(yōu)化算法如果作為一種通用的優(yōu)化算法框架，就需要經過優(yōu)化標準算例的驗證和確認，以評估算法的效率和精度，并對算法進行相應的改進。對于代理優(yōu)化算法的效率，一般通過目標函數計算次數來評估[107]；因為在實際工程設計問題中，調用高可信度數值模擬來計算目標函數往往是整個優(yōu)化流程中最費時的部分。目前，關于代理優(yōu)化算法的驗證和確認的研究工作開展得比較少。

4.1 標準測試函數優(yōu)化算例驗證

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

Rosenbrock函數算例[117]是測試優(yōu)化算法的經典算例之一。該函數最優(yōu)點位于一個狹長而平坦的峽谷內(如圖 11所示)。2維Rosenbrock函數算例優(yōu)化數學模型為

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

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

(80)

圖12給出了采用不同加點準則時，真實目標函數值的收斂情況?？梢?，對于單獨使用的加點準則，EI、局部EI和LCB效果較好，MSP、PI和MSE效果不佳。如果組合使用，“MSP+EI”和“MSP+LCB”的效果均優(yōu)于單獨使用的效果。圖13 給出了梯度增強型Kriging模型的優(yōu)化收斂情況?？梢娨胩荻刃畔⒑?，各種加點準則的收斂特性均得到改善；尤其是MSP準則的改善效果最明顯。表 1為采用代理優(yōu)化算法不同加點準則的優(yōu)化結果的比較，表明代理優(yōu)化算法不僅可以收斂到無約束問題的精確最優(yōu)解，而且具有很高的效率。

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

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

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

(81)

該函數在3維以上為多極值問題。

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

圖13 引入梯度信息后的代理優(yōu)化算法不同加點準則目標函數收斂歷程的比較(2維Rosenbrock函數算例)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)化結果的比較(2維Rosenbrock函數算例)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函數代理優(yōu)化Kriging模型和GEK模型的比較Fig.14 Comparison of Kriging and GEK models for 5-80 dimensional Rosenbrock functions surrogate optimization

G9函數算例是國際上用于測試全局優(yōu)化算法的一個經典算例。該算例含7個變量，4個約束(其中g1和g4為主動約束)，優(yōu)化數學模型為

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

4x6x7-10x6-8x7

(82)

該優(yōu)化問題是一個難度非常大的全局優(yōu)化問題[118]，原因在于：① 可行域很小，僅占整個設計空間的0.5%，屬于強約束問題；② 目標函數的變化幅度較大，從107量級到103量級。目前為止，國際上能得到的最優(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模型不同加點準則的優(yōu)化效率比較?？梢姡琈SP準則和LCB準則在未得到全局最優(yōu)解前滿足了優(yōu)化終止條件；EI和局部EI方法的結果均非常接近全局最優(yōu)；采用不同組合加點準則，只有“MSP+EI”收斂到了全局最優(yōu)。引入梯度信息后，代理優(yōu)化不同加點準則的優(yōu)化結果和優(yōu)化效率均得到顯著改善，但仍然只有EI、局部EI、“MSP+EI”準則獲得了全局最優(yōu)解。此外，代理優(yōu)化的結果不僅優(yōu)于直接遺傳算法優(yōu)化的結果，且調用真實函數分析次數小了1～2個量級。

表3給出了代理優(yōu)化結果與目前國際上能得到的最佳結果的比較，說明優(yōu)化結果已經非常接近國際上最佳結果，主動約束得到精確滿足。事實上，表中的國外最佳結果對g4約束是有所放寬的。這說明代理優(yōu)化算法針對復雜全局優(yōu)化問題具有很強的約束處理能力。

通過各種約束和無約束、局部和全局優(yōu)化的標準算例測試，充分說明代理優(yōu)化算法具有精確收斂到優(yōu)化問題的真實最優(yōu)解的能力。

4.2 氣動優(yōu)化設計標準算例確認

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

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

表3 7維G9函數優(yōu)化結果 (EI加點準則)Table 3 Optimum results of a 7D G9 function case(EI infill-sampling criterion)

MinCD

(83)

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

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

表4基準RAE2822翼型氣動性能網格收斂性

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翼型計算網格示意圖 Fig.15 Schematic of computational grid for RAE2822 airfoil

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

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

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)化翼型氣動性能網格收斂性

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

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

4.3 氣動/結構綜合優(yōu)化設計算例確認

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

MinW

(84)

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

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

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

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

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

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 機翼氣動/結構綜合優(yōu)化過程中目標函數收斂曲線Fig.19 Convergence history of objective function for wing aerostructural optimization

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

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

5 結 論

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

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

2) 代理優(yōu)化算法的優(yōu)化機制，是通過代理模型及其子優(yōu)化過程來產生新的樣本點。

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

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

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

6) 基于Kriging模型的代理優(yōu)化算法對于20維左右的工程設計問題已經基本成熟，在保證優(yōu)化質量的同時，能實現很高的優(yōu)化效率。

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

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

2) 大規(guī)模樣本數據問題。如果優(yōu)化過程中產生大量樣本點(2 000以上)，針對這些樣本點進行代理模型訓練的時間甚至超過精確數值分析的時間。

3) 對于梯度增強Kriging模型，隨著設計變量數增加，模型訓練的時間大大增加。

4) 復雜設計空間問題。主要針對設計空間呈超多極值、高度非線性、強烈振蕩等特性的復雜優(yōu)化問題。

5) 數值噪聲問題。樣本數據計算可能未充分收斂，這樣便帶有數值噪聲，從而對代理優(yōu)化的過程和結果產生重要影響。

經過文獻調研，作者認為今后在代理優(yōu)化方法研究方面還需要開展的工作如下：

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

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

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

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

5) 研究數據庫輔助的高效代理優(yōu)化算法。充分運用分析和優(yōu)化工作過程產生的大量歷史數據，甚至包括風洞試驗得到的試驗數據，對進一步的優(yōu)化設計形成強有力的數據支持。

6) 研究代理優(yōu)化算法在復雜系統多學科優(yōu)化設計問題的應用。包括具有大規(guī)模設計變量和大規(guī)模約束，且各學科之間具有復雜耦合關系的工程設計問題。

致 謝

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

[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] 余雄慶. 飛機總體多學科設計優(yōu)化的現狀與發(fā)展方向[J]. 南京航空航天大學學報, 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)星， 余雄慶, 等. 多學科設計優(yōu)化中常用代理模型的研究[J]. 計算力學學報, 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] 張健, 李為吉. 飛機多學科設計優(yōu)化中的近似方法分析[J]. 航空計算技術, 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模型的翼型氣動性能優(yōu)化設計[J]. 航空學報, 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] 許瑞飛， 宋文萍， 韓忠華. 改進Kriging模型在翼型氣動優(yōu)化設計中的應用研究[J]. 西北工業(yè)大學學報, 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)化設計研究[J]. 空氣動力學學報, 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] 白俊強, 王波, 孫智偉, 等. 基于松散式代理模型管理框架的亞音速機翼優(yōu)化設計方法研究[J]. 西北工業(yè)大學學報, 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模型在機翼氣動外形優(yōu)化中的應用[J]. 空氣動力學學報, 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模型的結構耐撞性優(yōu)化[J]. 南京航空航天大學學報, 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代理模型的結構形狀優(yōu)化方法研究[J]. 計算力學學報, 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] 楊軍, 劉勇瓊, 艾春安, 等. 改進Kriging模型在固沖發(fā)動機導彈氣動優(yōu)化設計中的應用[J]. 固體火箭技術, 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)化設計[J]. 航空計算技術, 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)星結構多目標優(yōu)化設計[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代理模型的高速列車頭型多目標優(yōu)化設計[J]. 中國科學, 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模型的渦輪葉片多學科設計優(yōu)化[J]. 航空動力學報, 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] 張科施, 韓忠華, 李為吉, 等. 基于近似技術的高亞聲速運輸機機翼氣動/結構優(yōu)化設計[J]. 航空學報, 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] 吳寬展， 劉學軍， 呂宏強. 超臨界翼型設計中的多響應代理模型[J]. 航空計算技術, 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] 劉俊. 基于代理模型的高效氣動優(yōu)化設計方法及應用[D]. 西安： 西北工業(yè)大學, 2015.

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

[114] 高月華， 王希誠. 基于Kriging代理模型的多點加點序列優(yōu)化方法[J]. 工程力學， 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] 劉俊, 宋文萍, 韓忠華, 等. 梯度增強的Kriging模型與Kriging模型在優(yōu)化設計中的比較研究[J]. 西北工業(yè)大學學報, 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] 張科施， 韓忠華， 李為吉， 等. 一種考慮氣動彈性的運輸機機翼多學科優(yōu)化方法[J]. 空氣動力學學報， 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.

韓忠華男， 博士， 教授， 博士生導師。主要研究方向： 代理模型理論與算法， 氣動與多學科優(yōu)化設計， 轉捩預測與自然層流翼型/機翼設計, 氣動數據庫技術。

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">網絡出版時間

時間：2016-03-291455

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

國家自然科學基金 (11272265)

*

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

韓忠華．Kriging模型及代理優(yōu)化算法研究進展J． 航空學報,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