徐圣冠，陳紅全，張加樂(lè)，高緩欽，賈雪松

（南京航空航天大學(xué)航空學(xué)院非定?？諝鈩?dòng)力學(xué)與流動(dòng)控制工信部重點(diǎn)實(shí)驗(yàn)室，江蘇南京 210016）

0 引 言

現(xiàn)代比較流行的隨機(jī)性啟發(fā)式算法如遺傳算法（genetic algorithms，GAs）、粒子群算法（particle swarm optimization，PSO）、模擬退火算法等（simulated annealing，SA），雖然具有全局搜索能力，但它們同時(shí)也存在收斂速度慢、需要函數(shù)評(píng)估次數(shù)多等缺點(diǎn)，對(duì)于函數(shù)評(píng)估比較耗時(shí)的工程問(wèn)題（如基于計(jì)算流體力學(xué)的優(yōu)化問(wèn)題），其計(jì)算量有時(shí)是不能接受的。另一類基于梯度的優(yōu)化方法又嚴(yán)重依賴于搜索起點(diǎn)，極易陷入局部最優(yōu)解。針對(duì)經(jīng)典優(yōu)化算法的不足，一種新的基于代理模型（surrogate model）優(yōu)化的方法為人們所青睞。它以一個(gè)代理模型的函數(shù)來(lái)代替耗時(shí)很高的精確模型的評(píng)估，能保證在一定精度的范圍內(nèi)快速得到最優(yōu)解。代理模型按照優(yōu)化過(guò)程中的構(gòu)造方式，又可以分為靜態(tài)代理模型和動(dòng)態(tài)代理模型。動(dòng)態(tài)代理模型由于在優(yōu)化過(guò)程中不斷優(yōu)化樣本空間，較靜態(tài)代理模型有更高的搜索效率和精度。Kriging代理模型由于其在獲得預(yù)測(cè)點(diǎn)預(yù)測(cè)值的同時(shí)能獲得此處預(yù)測(cè)值的標(biāo)準(zhǔn)方差，因此很適合構(gòu)造動(dòng)態(tài)代理模型。在此基礎(chǔ)上，Jones等提出了基于Kriging代理模型的高效全局優(yōu)化（EGO）算法，此算法以改善期望（EI）函數(shù)作為最優(yōu)加點(diǎn)策略，并用一組經(jīng)典的測(cè)試函數(shù)作了測(cè)試，發(fā)現(xiàn)其用較少的樣本點(diǎn)就能得到一定精度的全局最優(yōu)解。

經(jīng)典的EGO算法雖然能高效地獲得全局最優(yōu)解，但其精度仍存在進(jìn)一步改善的空間。這一方面是由于EGO算法是基于代理模型的，本身就存在一定的誤差；另一方面，目前沒(méi)有一個(gè)很好的收斂判斷準(zhǔn)則，導(dǎo)致EGO算法容易過(guò)早收斂，對(duì)于這一點(diǎn)，文獻(xiàn)［9］作了具體分析和研究，本文不再作具體討論。

本文圍繞最優(yōu)解精度進(jìn)一步改善的問(wèn)題，研究了三類面向精確最優(yōu)解的EGO算法。第一類算法考慮了與成熟的擬牛頓法和Powell法等局部?jī)?yōu)化方法的協(xié)同作用，第二類算法考慮在經(jīng)典EGO算法的改善期望函數(shù)中引入Kriging信任的作用。在對(duì)幾個(gè)常用檢驗(yàn)函數(shù)驗(yàn)證結(jié)果分析比較的基礎(chǔ)上提出了兩者優(yōu)缺點(diǎn)互補(bǔ)的第三類算法，并進(jìn)行了驗(yàn)證。結(jié)果表明，三類方法都能夠以相對(duì)較少的函數(shù)評(píng)估，得到比經(jīng)典EGO算法更為精確的全局最優(yōu)解。最后，本文選擇第三類算法并采用N-S方程流場(chǎng)求解器對(duì)一個(gè)具有10個(gè)設(shè)計(jì)變量的跨聲速翼型進(jìn)行了優(yōu)化驗(yàn)證，發(fā)現(xiàn)阻力系數(shù)在原有EGO優(yōu)化算法的基礎(chǔ)上有了進(jìn)一步的下降。

1 EGO算法的基本理論

1.1 Kriging代理模型

Kriging 代理模型是一種估計(jì)方差最小的無(wú)偏估計(jì)模型，最早由南非的Krige 在1951年提出。它可以表示為如下數(shù)學(xué)形式：

式中：f()為已知的退化模型；β為相應(yīng)的相關(guān)系數(shù)；()是一個(gè)均值為0、方差為的隨機(jī)函數(shù)。()的設(shè)計(jì)點(diǎn)x和x（、=1，2，…，，為樣本點(diǎn)個(gè)數(shù)）之間的協(xié)方差矩陣可以表示為

式中：是一個(gè)相關(guān)矩陣，它和數(shù)據(jù)點(diǎn)之間的空間分布相關(guān)。是一個(gè)對(duì)稱矩陣，其中一種常用的函數(shù)形式如下：

式中：θ≥0；0 ≤p≤2。θ值的大小反映了函數(shù)在第個(gè)坐標(biāo)軸上的非線性度，值越大表示非線性度越高。p值的大小反映了函數(shù)在第個(gè)坐標(biāo)軸上的光滑程度，其中常取p=2，此時(shí)函數(shù)被模擬為光滑函數(shù)。

式（6）中包含了新加入的預(yù)測(cè)點(diǎn)所帶來(lái)的誤差，通常這個(gè)值很小，可以忽略，于是得到一個(gè)簡(jiǎn)化的均方差MSE估計(jì)公式：

式（3）中的相關(guān)參數(shù)θ可由最大似然估計(jì)（MLE）在θ＞0的基礎(chǔ)上最大化下式得到：

式中：

1.2 改善期望函數(shù)EI

改善期望（expected improvement，EI）函數(shù)作為經(jīng)典EGO 算法的核心部分，是使EGO 算法具有全局性和高效性的重要保證。

它的概率為

進(jìn)一步得到期望為

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

從式（12）可以看出，無(wú)論是預(yù)測(cè)點(diǎn)值較小還是模型預(yù)測(cè)精度較低，都能使EI函數(shù)值變得較大，使最佳改善點(diǎn)偏向此處，這也是EGO算法具有高效性的原因所在。

1.3 經(jīng)典EGO算法的基本流程

經(jīng)典EGO的基本流程如下：

1）在設(shè)計(jì)空間中按照一定的試驗(yàn)設(shè)計(jì)方法（DOE）生成初始樣本點(diǎn)；

2）以這個(gè)初始樣本空間通過(guò)子優(yōu)化相關(guān)系數(shù)構(gòu)造Kriging代理模型；

3）以這個(gè)代理模型為基礎(chǔ)對(duì)其進(jìn)行最優(yōu)值搜索，得到此模型的最優(yōu)值；

4）以這個(gè)模型及模型的最優(yōu)值為基礎(chǔ)最優(yōu)化EI函數(shù)，獲得EI極大化時(shí)的最佳改善點(diǎn)；

5）判斷是否滿足終止條件，若滿足，對(duì)這個(gè)模型進(jìn)行最優(yōu)值搜索，輸出最優(yōu)值，優(yōu)化過(guò)程結(jié)束；若不滿足，則轉(zhuǎn)到第6步；

6）對(duì)這個(gè)最佳改善點(diǎn)用高精度函數(shù)評(píng)估得到改善點(diǎn)的響應(yīng)值，加入樣本空間回到第2步。

在本文中，試驗(yàn)設(shè)計(jì)方法采用拉丁超立方（LHS）采樣方法，初始樣本點(diǎn)個(gè)數(shù)全部取為5 個(gè)。對(duì)于EGO 算法中的兩個(gè)子優(yōu)化問(wèn)題，由于類似于函數(shù)優(yōu)化，本文則采用自適應(yīng)遺傳算法（AGA），并加入了自適應(yīng)變化的多點(diǎn)交叉技術(shù)和生存壓力函數(shù)來(lái)改善AGA 的全局尋優(yōu)能力，并采用了Open-MP 并行加速技術(shù)來(lái)加速整個(gè)優(yōu)化過(guò)程。

2 高效高精度EGO算法

基于Kriging 代理模型的經(jīng)典EGO 算法有其固有的近似性。要消除這種近似性，本文考慮在最終尋優(yōu)過(guò)程中不使用代理模型的替代方法。

2.1 EGO與局部?jī)?yōu)化協(xié)同算法

由于基于梯度等優(yōu)化算法具有強(qiáng)大的局部搜索能力?？紤]EGO 算法和此類局部?jī)?yōu)化算法協(xié)同，在EGO 算法收斂時(shí)用局部?jī)?yōu)化算法繼續(xù)搜尋，消除EGO算法中代理模型帶來(lái)的誤差，進(jìn)一步地搜尋得到更為精確的最優(yōu)解。本文選擇了比較成熟的需要計(jì)算梯度的擬牛頓法（BFGS）和不需要計(jì)算梯度的Powell法。

BFGS 法是目前為止不用Hessian 矩陣的算法中最有效的算法，被公認(rèn)為是最好的擬牛頓法。它是由Broyden、Fletcher、Goldfarb 和Shanno 于1970年各自從不同角度出發(fā)得到的。本文采用的有限內(nèi)存BFGS 方法是由Nocedal在1980年提出的，簡(jiǎn)記為L(zhǎng)BFGS 方法。此方法無(wú)需記憶矩陣，解決了存儲(chǔ)量大的問(wèn)題，非常適合解決大規(guī)模問(wèn)題，具體過(guò)程參見(jiàn)文獻(xiàn)［17］。

Powell法和其他局部?jī)?yōu)化方法相比，其最大的特點(diǎn)是不需要計(jì)算梯度。因?yàn)閷?duì)目標(biāo)函數(shù)的解析性質(zhì)要求不高，所以適用范圍比擬牛頓法等需要求解梯度的優(yōu)化算法更為廣泛。但由于不利用函數(shù)的解析性質(zhì)，因此其收斂也可能較慢，具體過(guò)程參見(jiàn)文獻(xiàn)［18］。

本文考慮EGO算法和局部?jī)?yōu)化算法的協(xié)同作用，其過(guò)程可描述為：當(dāng)EGO 算法收斂時(shí)，用LBFGS 或Powell 算法以EGO 算法得到的最優(yōu)解作為起點(diǎn)進(jìn)行最優(yōu)化搜索，直至收斂。為下文中描述方便，這兩種方法分別簡(jiǎn)記為EGO-LB和EGO-P。

2.2 改善最優(yōu)加點(diǎn)策略

與傳統(tǒng)的EGO算法不同，本文考慮在最終最優(yōu)解的搜尋過(guò)程中（見(jiàn)1.3節(jié)算法流程第5步），不是基于代理模型，而是直接基于樣本點(diǎn)尋優(yōu)，這樣就繞開(kāi)了最終代理模型帶來(lái)的近似性。但這樣做有一個(gè)明顯的缺點(diǎn)，就是要求全局最優(yōu)解對(duì)應(yīng)的點(diǎn)必須包含在最終的樣本空間里面。EI 最優(yōu)加點(diǎn)策略在計(jì)算后期由于全局最優(yōu)解附近改善點(diǎn)大量聚集，使得該區(qū)域的標(biāo)準(zhǔn)差變得非常小，在后期會(huì)使最佳改善點(diǎn)偏向外圍區(qū)域，造成全局最優(yōu)解對(duì)應(yīng)的最佳改善點(diǎn)很難被搜尋到，因此考慮對(duì)EI最優(yōu)加點(diǎn)策略進(jìn)行改善。

本文采用考慮KB作用的EI最優(yōu)加點(diǎn)策略。在優(yōu)化開(kāi)始階段采用EI最優(yōu)加點(diǎn)策略，在EGO算法收斂時(shí)，改用KB法作為最優(yōu)加點(diǎn)策略。當(dāng)KB得到的最佳改善點(diǎn)和前一步所得到的重合時(shí)返回EI，如此往復(fù)，直至收斂。為下文中描述方便，此算法簡(jiǎn)記為EGO-K。

2.3 數(shù)值驗(yàn)證

本文選用文獻(xiàn)中常用的測(cè)試不同復(fù)雜度的3個(gè)檢驗(yàn)函數(shù)。

1）Branin函數(shù)

式中：-5 ≤≤10，0 ≤≤15，該函數(shù)的最小值為-1.185 9。

2）Hartman3函數(shù)

0 ≤x≤1（=1，2，3）。該函數(shù)的最小值為-3.862 8。

3）Hartman6函數(shù)

0 ≤x≤1（=1，2，…，6）。該函數(shù)的最小值為-3.322 4。

采用上述發(fā)展的EGO、EGO-LB、EGO-P 和EGO-K 算法，對(duì)以上3 個(gè)函數(shù)進(jìn)行了數(shù)值實(shí)驗(yàn)，其結(jié)果如表1所示，表中數(shù)據(jù)從上往下依次為相對(duì)誤差和得到最優(yōu)解時(shí)的函數(shù)計(jì)算次數(shù)。

表1 數(shù)值驗(yàn)證結(jié)果及對(duì)比Tab.1 Test function results and comparisons

從表1中可以看出：

1）本文發(fā)展的EGO 算法比經(jīng)典的EGO 算法具有更高的搜索效率和精度，這是由于本文在子優(yōu)化過(guò)程中采用的是具有很強(qiáng)全局搜索能力的AGA 算法，而經(jīng)典EGO算法采用的是直接算法。

2）EGO 算法與局部?jī)?yōu)化算法協(xié)同來(lái)提高EGO算法精度的代價(jià)是高昂的，其函數(shù)計(jì)算次數(shù)要比原有的EGO算法多出一倍左右，并且會(huì)隨著問(wèn)題規(guī)模和問(wèn)題復(fù)雜程度的增加而增大。但相比于其他優(yōu)化算法，與局部?jī)?yōu)化算法協(xié)同的EGO算法還是有巨大優(yōu)勢(shì)的，如對(duì)于Branin 問(wèn)題，其他比如GAs 算法要計(jì)算300 多次才能得到函數(shù)值的最優(yōu)解。

3）從搜索精度上來(lái)看，EGO-LB、EGO-P 和EGO-K 算法中，除了EGO-K 中處理Hartman6 函數(shù)時(shí)最終還有點(diǎn)誤差外（但也已非常接近最優(yōu)解），其他幾種方法對(duì)于這3 個(gè)測(cè)試函數(shù)都能達(dá)到精確的全局最優(yōu)解。

4）從搜索效率上來(lái)看，EGO-LB 效率最差，需要額外計(jì)算的函數(shù)次數(shù)最多，這和它需要計(jì)算函數(shù)梯度有關(guān)。表1中EGO-K 的函數(shù)計(jì)算次數(shù)包括中間判斷初步收斂的次數(shù)，實(shí)際當(dāng)KB 啟動(dòng)時(shí)分別只經(jīng)過(guò)了7次、4 次和2 次就達(dá)到了全局最優(yōu)解。從圖1中Hartman6 函數(shù)收斂曲線中可以看到，當(dāng)KB 發(fā)生作用時(shí)只經(jīng)過(guò)了2 次函數(shù)計(jì)算就達(dá)到了最優(yōu)解，搜索效率非常高。EGO-P 計(jì)算效率介于EGO-LB 和EGO-K之間。

圖1 EGO/EGO-K優(yōu)化Hartmat6函數(shù)收斂過(guò)程Fig.1 Convergence process of Hrtman6 function with EGO/EGO-K algorithm

5）從比較中可以看出，采用考慮Kriging 信任的改善期望函數(shù)的EGO算法具有很高的搜索效率，但當(dāng)遇到高維復(fù)雜、初始模型構(gòu)造誤差較大的情況時(shí)，其精度不如與局部?jī)?yōu)化算法協(xié)同的算法，但后者的效率又大不如前者。既然兩者優(yōu)缺點(diǎn)互補(bǔ)，那么將兩種方案協(xié)同作用在理論上應(yīng)該可行。下一節(jié)討論兩者優(yōu)缺點(diǎn)互補(bǔ)的改進(jìn)算法。

2.3 改進(jìn)算法

在EGO算法和局部?jī)?yōu)化算法協(xié)同作用的過(guò)程中，改用考慮Kriging信任的改善期望函數(shù)，形成全局優(yōu)化的改進(jìn)算法，簡(jiǎn)記為EGO-KLB 和EGO-KP，分別代表本文的EGO-K 算法與LBFGS 法、Powell 法耦合的算法，程序流程如圖2所示。

圖2 EGO-KLB/EGO-KP程序流程Fig.2 Flow chart of EGO-KLB/EGO-KP

對(duì)Branin函數(shù)和Hartman3函數(shù)，兩種方案都能得到精確的最優(yōu)解，故這里取Hartman6 函數(shù)來(lái)驗(yàn)證算法。這里兩個(gè)判斷參數(shù)取為=10，=15，為2%精度范圍內(nèi)連續(xù)不變代數(shù)。

由表2可以看出：EGO-KLB 和EGO-KP 在得到Hartman6函數(shù)精確最優(yōu)解時(shí)，所用的函數(shù)計(jì)算次數(shù)分別比EGO-LB 和EGO-P 減少了10.5%和7.9%，比經(jīng)典EGO 算法獲得精確全局最優(yōu)解的函數(shù)計(jì)算次數(shù)分別減少了8.3%和31.0%。由此說(shuō)明，此方案在獲得精確最優(yōu)解的同時(shí)，能降低改進(jìn)EGO 算法和經(jīng)典EGO算法的函數(shù)計(jì)算次數(shù)，說(shuō)明此優(yōu)化算法在一定程度上是有效可行的。

表2 改進(jìn)算法數(shù)值驗(yàn)證結(jié)果Tab.2 Test function results of the improved EGO

3 跨音速翼型優(yōu)化算例

為了驗(yàn)證本文面向精確最優(yōu)解的EGO 算法的工程價(jià)值，選取了跨音速翼型阻力極小化優(yōu)化算例來(lái)進(jìn)行驗(yàn)證。優(yōu)化算法選取了精度和效率相對(duì)都比較高的EGO-KP法。

本文以RAE2822為參考翼型，設(shè)計(jì)要求是在來(lái)流馬赫數(shù)0.729、攻角為2.55°、雷諾數(shù)為6.5×10時(shí)，翼型的阻力最小，同時(shí)要求翼型的截面面積和升力系數(shù)在一定范圍內(nèi)不減小。此時(shí)該優(yōu)化問(wèn)題可以表示為

適應(yīng)度函數(shù)取為

升力系數(shù)約束罰函數(shù)為

面積約束罰函數(shù)為

翼型參數(shù)化方法采用CST 參數(shù)法。CST 參數(shù)法能以較少的參數(shù)表達(dá)完整的翼型信息，具有使用方便、直觀、控制參數(shù)可控等優(yōu)點(diǎn)。本文取了10 個(gè)設(shè)計(jì)參數(shù)。

流場(chǎng)解算器采用Fluent 耦合隱式N-S 方程求解器，湍流模型采用S-A全湍流模型，采用多重網(wǎng)格加速技術(shù)，計(jì)算網(wǎng)格量為36 024 個(gè)。本算例中取兩個(gè)判斷參數(shù)=80，=100（即當(dāng)代最優(yōu)值保持不變時(shí)，KB 開(kāi)始作用；代最優(yōu)值保持不變時(shí)，Powell 法開(kāi)始作用）。優(yōu)化結(jié)果和收斂曲線見(jiàn)表3和圖3，從中可以看出：EGO 算法在計(jì)算67 次流場(chǎng)以后，就收斂到了全局最優(yōu)解附近；在經(jīng)過(guò)80 次最優(yōu)值保持不變以后，改用KB 和Powell 法，可以看出KB 和Powell 法使最優(yōu)解在原有的基礎(chǔ)上再次提高了1.11%；同樣，EGOKP 算法的函數(shù)計(jì)算次數(shù)中也包括了判斷初步收斂的80次函數(shù)計(jì)算。

表3 CD優(yōu)化結(jié)果及對(duì)比Tab.3 Optimization results and comparison of CD

圖3 翼型優(yōu)化收斂曲線Fig.3 Convergence curve of the airfoil optimization

表4給出了本文最終優(yōu)化得到的翼型和RAE2822的對(duì)比。從表4中優(yōu)化結(jié)果可以看出：優(yōu)化翼型阻力系數(shù)比RAE2822翼型減小了11.38%，升力系數(shù)和面積略有下降，但都在約束允許范圍之內(nèi)。圖4給出了優(yōu)化后翼型與RAE2822翼型的外形對(duì)比，其中，為翼型弦長(zhǎng)。從圖5～7中可以看出：優(yōu)化后翼型表面的激波已基本消失，表明優(yōu)化得到的阻力減小主要源自激波減阻。

表4 優(yōu)化結(jié)果和參考翼型對(duì)比Tab.4 Optimization results and comparison

圖4 優(yōu)化翼型與初始翼型外形對(duì)比Fig.4 Sharp of the optimized airfoil and RAE2822

圖5 優(yōu)化翼型與初始翼型壓強(qiáng)系數(shù)對(duì)比Fig.5 Pressure coefficient of the optimized airfoil and RAE2822

圖6 RAE2822翼型壓強(qiáng)系數(shù)等值線Fig.6 Pressure coefficient contour line of RAE2822

圖7 優(yōu)化翼型壓強(qiáng)系數(shù)等值線Fig.7 Pressure coefficient contour line of the optimized airfoil

4 結(jié)束語(yǔ)

為了提高EGO 算法的優(yōu)化精度，本文發(fā)展了與Kriging 信任法、牛頓法、Powell 法等局部?jī)?yōu)化算法協(xié)同的優(yōu)化算法。整體上，本文所發(fā)展的改進(jìn)算法都能以相對(duì)較少的額外函數(shù)評(píng)估，得到比經(jīng)典EGO 算法更精確的全局最優(yōu)解。具體地，第1 類改進(jìn)的EGO與局部?jī)?yōu)化協(xié)同算法（EGO-LB 和EGO-P）得到的全局最優(yōu)解具有較高的精度，但效率相對(duì)較低；第2 類考慮Kriging 信任的改善期望函數(shù)的EGO 算法（EGO-K）具有很高的效率，但高維測(cè)試函數(shù)算例得到的最優(yōu)解，其精度不如第1 類算法；第3 類改進(jìn)算法（EGO-KLB 和EGO-KP）耦合了第1 類和第2 類算法的優(yōu)點(diǎn)，能以比經(jīng)典EGO 算法更少的函數(shù)評(píng)估，獲得更精確的全局最優(yōu)解。最后，把發(fā)展的算法應(yīng)用到具體的跨音速翼型優(yōu)化問(wèn)題中，優(yōu)化結(jié)果顯示，得到的阻力較原EGO 算法有進(jìn)一步的減小，優(yōu)化翼型上的激波基本消失（基本無(wú)波阻），展示了一定的工程價(jià)值。