多目標(biāo)進(jìn)化算法和代理模型技術(shù)在氣動穩(wěn)健優(yōu)化設(shè)計(jì)中的應(yīng)用

李焦贊，高正紅

（西北工業(yè)大學(xué) 翼型葉柵空氣動力學(xué)國家級科技重點(diǎn)實(shí)驗(yàn)室，陜西 西安710072）

0 引 言

傳統(tǒng)的優(yōu)化設(shè)計(jì)模型中，設(shè)計(jì)變量、目標(biāo)函數(shù)和約束函數(shù)均為確定性的，不能夠反映工程設(shè)計(jì)和產(chǎn)品制造過程中以及服役環(huán)境等客觀存在的大量不確定性因素對設(shè)計(jì)方案的影響，使得傳統(tǒng)優(yōu)化設(shè)計(jì)的最優(yōu)解可能對于這些不確定性因素非常敏感。穩(wěn)健優(yōu)化設(shè)計(jì)是穩(wěn)健設(shè)計(jì)和優(yōu)化設(shè)計(jì)兩種方法的結(jié)合，就是通過調(diào)整設(shè)計(jì)變量的名義值和控制其偏差來保證設(shè)計(jì)最優(yōu)解的穩(wěn)健性。與傳統(tǒng)的確定性優(yōu)化問題相比，穩(wěn)健優(yōu)化的一個(gè)顯著特點(diǎn)是在優(yōu)化過程中需要不斷地進(jìn)行不確定性分析，應(yīng)用優(yōu)化的手段不斷減小設(shè)計(jì)目標(biāo)對不確定性因數(shù)的依賴程度，使設(shè)計(jì)結(jié)果對不確定性因數(shù)變得不敏感，即當(dāng)使用環(huán)境或產(chǎn)品本身的參數(shù)在一定的范圍內(nèi)變化時(shí)，產(chǎn)品的性能變化很?。?－2，4，9］。

根據(jù)穩(wěn)健設(shè)計(jì)各種不同應(yīng)用目的，可派生出兩種類型的優(yōu)化設(shè)計(jì)問題［3］，第一種為通過調(diào)整設(shè)計(jì)變量的選擇減小由于設(shè)計(jì)參數(shù)的變差引起系統(tǒng)輸出的方差；另一種為通過設(shè)計(jì)變量的選擇減小設(shè)計(jì)變量變差引起系統(tǒng)輸出的方差。應(yīng)用到飛行器外形氣動穩(wěn)健設(shè)計(jì)中，不確定性因素也就分為以下兩類考慮：一類為飛行狀態(tài)（設(shè)計(jì)點(diǎn)）的不確定性，通常考察馬赫數(shù)、迎角等的波動對氣動性能的影響；另一類為幾何參數(shù)的不確定性。作者于其他文獻(xiàn)中對幾何參數(shù)不確定性進(jìn)行了詳細(xì)研究，本文重點(diǎn)研究考慮飛行參數(shù)存在不確定性變化時(shí)的氣動優(yōu)化問題。

1 穩(wěn)健優(yōu)化設(shè)計(jì)模型

穩(wěn)健優(yōu)化一般分為三部分：（1）分析確認(rèn)設(shè)計(jì)輸入的不確定性因素并由此得到它們的概率分布狀況；（2）指認(rèn)不確定性分析方法用以考察不確定性在系統(tǒng)中的傳播，從而得到穩(wěn)健設(shè)計(jì)所需的目標(biāo)函數(shù)和約束等，通常用均值和方差表示；（3）優(yōu)化搜索穩(wěn)健設(shè)計(jì)目標(biāo)函數(shù)，并使其滿足約束條件。這里我們考察飛行環(huán)境參數(shù)如馬赫數(shù)、迎角存在不確定性波動時(shí)的情況，假定輸入變量的分布已知為均勻分布，于是（2）和（3）分析如下。

1.1 不確定性分析

在穩(wěn)健設(shè)計(jì)中通常用均值μ來表示產(chǎn)品所期望的性能指標(biāo)，用標(biāo)準(zhǔn)差σ來描述其穩(wěn)健性。已知設(shè)計(jì)問題的輸入?yún)?shù)并通過計(jì)算分析，求解產(chǎn)品性能隨機(jī)特征的過程稱為不確定性分析過程。穩(wěn)健優(yōu)化問題要進(jìn)行不確定分析，即在穩(wěn)健優(yōu)化問題中要計(jì)算輸出性能的均值μ和標(biāo)準(zhǔn)差σ。目前已提出多種簡化求解方法，如蒙特卡羅取樣法、基于泰勒級數(shù)展開、試驗(yàn)設(shè)計(jì)方法等。本文選取蒙特卡羅取樣法做不確定性均值和方差分析。取樣法的特點(diǎn)是要進(jìn)行大量的隨機(jī)取樣，以得到目標(biāo)的概率均值和方差。對于翼型穩(wěn)健設(shè)計(jì)，由于求解翼型氣動系數(shù)的計(jì)算量較大，它們的均值和方差更是難以求解，于是應(yīng)用代理模型技術(shù)代替流場求解，然后基于代理模型進(jìn)行不確定性分析。各不確定性指標(biāo)如升阻力力矩等的均值和標(biāo)準(zhǔn)差，可離散求解分析如下：

其中，CLi（x）、CDi（x）是關(guān)于x的函數(shù)，表示在x的波動范圍內(nèi)抽取n組隨機(jī)輸入分別對應(yīng)的升阻力系數(shù)；x是不確定性向量（幾何參數(shù)或其他不確定性參數(shù)），在設(shè)計(jì)點(diǎn)附近波動并服從某一分布。

1.2 多目標(biāo)穩(wěn)健優(yōu)化模型

針對二維翼型設(shè)計(jì)，本文根據(jù)穩(wěn)健設(shè)計(jì)的多目標(biāo)特性，提出包含四個(gè)優(yōu)化目標(biāo)的穩(wěn)健設(shè)計(jì)模型。選取RAE2822作為初始翼型，飛行環(huán)境參數(shù)的不確定性因素主要關(guān)注馬赫數(shù)（M＝0.73±0.05）和迎角（a＝2°±0.5°）在一定范圍內(nèi)變化，考察優(yōu)化目標(biāo)分別為最小化阻力的均值和最大化升力的均值，以求得到更好的氣動性能；同時(shí)最小化阻力和升力的標(biāo)準(zhǔn)差，以使翼型保持優(yōu)異的穩(wěn)健性。從數(shù)學(xué)的角度描述穩(wěn)健模型為：

并滿足以下約束條件：

設(shè)定優(yōu)化問題需滿足以下約束：翼型厚度不低于原始厚度即保持在12%弦長以上，對比升力阻力力矩的均值和標(biāo)準(zhǔn)差，優(yōu)化結(jié)果的氣動特性不低于原始RAE2822翼型的性能。優(yōu)化過程中，對約束條件進(jìn)行罰函數(shù)處理。

2 翼型穩(wěn)健優(yōu)化設(shè)計(jì)流程

把上文建立的穩(wěn)健優(yōu)化模型應(yīng)用到具體翼型優(yōu)化設(shè)計(jì)中，于是建立基于代理模型的穩(wěn)健優(yōu)化設(shè)計(jì)流程如圖1所示。將穩(wěn)健優(yōu)化應(yīng)用于氣動設(shè)計(jì)的一個(gè)障礙是所需的不確定性分析計(jì)算量太大，因此采用代理模型技術(shù)［6］來進(jìn)行不確定分析。下文詳細(xì)表述了翼型多目標(biāo)穩(wěn)健優(yōu)化設(shè)計(jì)框架中主要的相關(guān)技術(shù)：翼型參數(shù)化方法，流場求解技術(shù)，代理模型技術(shù)及多目標(biāo)優(yōu)化算法等。

圖1 基于代理模型的多目標(biāo)穩(wěn)健優(yōu)化設(shè)計(jì)流程圖Fig.1 Flow chart of multi－objective robust optimization design based on surrogate model

2.1 翼型幾何參數(shù)化方法

常用的翼型參數(shù)化方法有：B樣條法，Hicks－Henne型函數(shù)疊加法，PARSEC方法和CST方法等。針對本文的研究，我們選擇PARSEC方法這種以直觀描述翼型幾何參數(shù)為設(shè)計(jì)變量并且其個(gè)數(shù)較少的參數(shù)化方法作為翼型優(yōu)化的幾何外形生成方法。PARSEC方法最早由H.Sobieczky提出，研究表明PARSEC方法對氣動優(yōu)化問題能以較小的設(shè)計(jì)變量個(gè)數(shù)得到更好的設(shè)計(jì)結(jié)果，也有研究顯示PARSEC方法對已有外形的精確擬合有時(shí)會失敗。本文在優(yōu)化過程中，所有翼型的后緣厚度取為0，后緣點(diǎn)固定，最終設(shè)計(jì)變量共9個(gè)，如圖2所示，各參數(shù)的定義為前緣半徑rLE、上表面最大厚度位置（xup，zup）、上表面最大厚度處的曲率zxxup、下表面最大厚度位置（xlw，zlw）、下表面最大厚度處的曲率zxxlw、后緣方向角αTE、后緣夾角βTE等。

圖2 PARSEC翼型參數(shù)化Fig.2 Airfoil parameterization using PARSEC

2.2 CFD流場求解器

流場求解選擇課題組開發(fā)的求解守恒形式雷諾平均NS方程［6］，自動生成繞翼型的C型網(wǎng)格，如圖3所示。采用中心差分格式進(jìn)行空間離散，五步龍格庫塔法進(jìn)行時(shí)間推進(jìn)，選取B－L湍流模型，翼型上下表面作全湍流計(jì)算處理。

圖3 流場求解中的翼型網(wǎng)格圖Fig.3 Airfoil mesh used in flow computing

2.3 RBF神經(jīng)網(wǎng)絡(luò)代理模型

代理模型技術(shù)的應(yīng)用為大規(guī)模優(yōu)化設(shè)計(jì)的響應(yīng)求解提供了可能，特別是CFD流場求解方面。本文選擇RBF神經(jīng)網(wǎng)絡(luò)模型［7］作為優(yōu)化設(shè)計(jì)中使用的代理模型，訓(xùn)練代理模型時(shí)所用樣本點(diǎn)的試驗(yàn)設(shè)計(jì)抽樣方法采用拉丁超立方方法。

徑向函數(shù)是一類以待測點(diǎn)與中心點(diǎn)之間的歐氏距離為自變量的函數(shù)，以徑向函數(shù)為基函數(shù)，通過線性疊加構(gòu)造出來的模型即為徑向基函數(shù)模型。RBF神經(jīng)網(wǎng)絡(luò)模型在函數(shù)逼近方面有其特有的優(yōu)良性質(zhì)，其基本形式如下：

其中x為輸入向量，φ（·）為徑向基函數(shù)，ci為第i個(gè)隱層節(jié)點(diǎn)徑向基函數(shù)的中心，σi為該徑向基函數(shù)的寬度，λi為該節(jié)點(diǎn)與輸出節(jié)點(diǎn)的權(quán)值，Nh為隱層節(jié)點(diǎn)的數(shù)量，λ0為偏置值，‖x－ci‖為x到中心ci的徑向距離，這里范數(shù)通常取歐幾里德二次范數(shù)：

本文選典型的高斯函數(shù)作為徑向基函數(shù)：φ（r）＝，σ表示徑向基函數(shù)的寬度。

由于RBF神經(jīng)網(wǎng)絡(luò)中基函數(shù)的個(gè)數(shù)并不等于樣本的個(gè)數(shù)，中心也并不一定選擇在輸入樣本點(diǎn)上并且增加了偏置量，這使得神經(jīng)網(wǎng)絡(luò)不嚴(yán)格通過每個(gè)樣本，從而有效地抑制了數(shù)值噪聲帶來的影響，提高了代理模型對設(shè)計(jì)空間非樣本點(diǎn)做出準(zhǔn)確預(yù)測的泛化能力。訓(xùn)練代理模型的具體步驟如下：

（1）在設(shè)計(jì)空間內(nèi)應(yīng)用拉丁超立方方法選取樣本，樣本盡可能最大限度地反映整個(gè)設(shè)計(jì)空間的信息；

（2）采用原始NS模型對樣本進(jìn)行氣動特性計(jì)算，獲得樣本的目標(biāo)特性數(shù)據(jù)；

（3）訓(xùn)練RBF神經(jīng)網(wǎng)絡(luò)代理模型，通過樣本的帶入學(xué)習(xí)可以確定模型的權(quán)值系數(shù)等參數(shù)從而得到神經(jīng)網(wǎng)絡(luò)模型輸入輸出之間的映射關(guān)系。

2.4 多目標(biāo)優(yōu)化搜索方法

穩(wěn)健優(yōu)化設(shè)計(jì)是個(gè)典型的多目標(biāo)優(yōu)化問題，目前解決這一類問題最典型的方法是基于Pareto最優(yōu)解集的多目標(biāo)遺傳進(jìn)化算法，本文選擇NCGA多目標(biāo)遺傳算法［8］，也稱為鄰域培植遺傳算法。NCGA有效地繼承了前人的成果，交叉和變異方法與標(biāo)準(zhǔn)GA方法相同，但是引入了“近鄰交叉”的概念。近鄰交叉是指交叉不是在具有完全不同的遺傳因子信息的個(gè)體間進(jìn)行，而是在具有一定程度的類似性的個(gè)體之間進(jìn)行，這樣的操作使得子代個(gè)體與父代個(gè)體比較接近，有利于提高GA的全局探索能力。NCGA能夠生成具有多樣性的解，因?yàn)樵谀繕?biāo)空間里沒有重復(fù)的個(gè)體，比起多目標(biāo)進(jìn)化算法中較為典型的NSGA－II和SPEA2有更加容易生成更多樣性解的傾向，“近鄰交叉”也被確認(rèn)為在多峰值問題上可以得到好的結(jié)果。

找到Pareto前緣后，下一步工作就是從中選擇最優(yōu)的設(shè)計(jì)點(diǎn)。當(dāng)Pareto前緣中非劣解不很多時(shí)，設(shè)計(jì)者能容易地對比不同結(jié)果之間的差別并從中選出滿足設(shè)計(jì)要求的最優(yōu)解；但當(dāng)設(shè)計(jì)目標(biāo)較多且前緣解較多時(shí)，非劣解集就會比較復(fù)雜，設(shè)計(jì)者就需要一定的準(zhǔn)則來指導(dǎo)其選出滿足設(shè)計(jì)要求的最優(yōu)解，多準(zhǔn)則決策方法就是一個(gè)很好的指導(dǎo)準(zhǔn)則，本文用柱狀圖法表示了一部分非劣解與初始解的性能對比情況，從中可以比較直觀地挑選出滿足設(shè)計(jì)要求的結(jié)果，達(dá)到篩選非劣解的目的。

3 算例結(jié)果與分析

根據(jù)上文建立的翼型穩(wěn)健設(shè)計(jì)模型，選擇RAE2822為初始翼型進(jìn)行穩(wěn)健優(yōu)化設(shè)計(jì)實(shí)例，考察飛行狀態(tài)參數(shù)馬赫數(shù)和迎角存在不確定性波動時(shí)的翼型綜合性能最優(yōu)設(shè)計(jì)。設(shè)計(jì)狀態(tài)為：馬赫數(shù)0.73，雷諾數(shù)6500000，迎角2°，設(shè)計(jì)升力系數(shù)0.69。設(shè)定多目標(biāo)遺傳算法參數(shù)為種群數(shù)60個(gè)和進(jìn)化代數(shù)30代，應(yīng)用拉丁超立方方法產(chǎn)生500個(gè)樣本訓(xùn)練RBF神經(jīng)網(wǎng)絡(luò)代理模型，利用圖1的優(yōu)化框架進(jìn)行穩(wěn)健優(yōu)化設(shè)計(jì)。本文優(yōu)化問題定義四個(gè)目標(biāo)函數(shù)，于是最終結(jié)果不是單一的而是一組Pareto解，每一個(gè)設(shè)計(jì)結(jié)果對應(yīng)著四個(gè)目標(biāo)函數(shù)，其優(yōu)化結(jié)果的Pareto前緣如圖4所示，圖4（a）表示升力均值對應(yīng)阻力均值，圖4（b）表示升力標(biāo)準(zhǔn)差對應(yīng)阻力標(biāo)準(zhǔn)差，他們共同體現(xiàn)了Pareto前緣的組成情況。圖中用數(shù)字標(biāo)出選擇的設(shè)計(jì)結(jié)果在Pareto前緣中的位置，也同時(shí)標(biāo)出了初始翼型RAE2822均值和標(biāo)準(zhǔn)差對應(yīng)的位置，可以看出優(yōu)化進(jìn)行成功：我們得到了比初始翼型阻力更小（圖4（a）縱坐標(biāo)減小方向）并且升阻力標(biāo)準(zhǔn)差也顯著減?。▓D4（b）的縱橫坐標(biāo)減小方向）的設(shè)計(jì)結(jié)果，即優(yōu)化結(jié)果保持很好的穩(wěn)健性。從圖4（a）中也可以看出，一些結(jié)果的升力均值比初始翼型有些許的惡化，可能與約束罰函數(shù)處理的偏差有關(guān)。

圖4 Pareto前緣與初始翼型升阻力均值與方差對比Fig.4 Pareto frontier of mean and Std deviation and also the value of the airfoil RAE2822

為了直觀地表示Pareto解，采用柱狀圖的形式表示了設(shè)計(jì)結(jié)果與初始解的目標(biāo)性能，如圖5所示。其中每個(gè)設(shè)計(jì)結(jié)果由4個(gè)柱狀圖組成，每個(gè)柱狀圖代表結(jié)果的一個(gè)目標(biāo)值，柱狀圖的高低代表了目標(biāo)值的大小，每類目標(biāo)柱狀圖分別用初始翼型RAE2822的對應(yīng)目標(biāo)值做了歸一化處理以便進(jìn)行數(shù)值對比。對比柱狀圖的高低體現(xiàn)了Pareto解集的收斂情況，可以看出與2822相比，設(shè)計(jì)結(jié)果升阻力標(biāo)準(zhǔn)差有明顯的減小，體現(xiàn)了設(shè)計(jì)的穩(wěn)健性；阻力均值也有很大的減小，體現(xiàn)減阻設(shè)計(jì)的性能要求；而升力的均值保持不減小或些許增大，也體現(xiàn)了設(shè)計(jì)升力不減小的要求，于是可以認(rèn)為優(yōu)化設(shè)計(jì)是成功的。

圖5 翼型性能及穩(wěn)健性的柱狀圖對比Fig.5 Histogram of the airfoils＇performance and robustness comparison

表格1所示為圖5設(shè)計(jì)結(jié)果所對應(yīng)的具體數(shù)值，對基于代理模型的設(shè)計(jì)結(jié)果同時(shí)進(jìn)行了精確NS方程求解驗(yàn)算，對比了他們與初始翼型RAE2822的性能區(qū)別。不論是用代理模型優(yōu)化的結(jié)果還是用NS方程校驗(yàn)的結(jié)果都比初始翼型的性能優(yōu)異，特別是阻力系數(shù)，均值減小的同時(shí)方差也有明顯的下降，而升力系數(shù)在保持均值不減小的情況下方差有一些減小。圖6所示為其中一個(gè)設(shè)計(jì)結(jié)果與翼型RAE2822的外形對比，可以看出優(yōu)化結(jié)果的最大厚度位置后移，有效減弱了激波強(qiáng)度從而減小波阻?？傊€(wěn)健優(yōu)化算法找到了比初始翼型性能更優(yōu)并且更穩(wěn)健的設(shè)計(jì)解。

表1 設(shè)計(jì)翼型與初始翼型的性能及穩(wěn)健性對比Table 1 Performance and robustness comparison between original airfoil RAE2822and designs

圖6 設(shè)計(jì)翼型與初始翼型幾何外形對比Fig.6 Geometry profile comparison between design airfoils and the original airfoil RAE2822

4 結(jié) 論

本文建立了考慮飛行環(huán)境參數(shù)存在不確定性變化因素時(shí)的翼型穩(wěn)健設(shè)計(jì)優(yōu)化模型，為了權(quán)衡性能和穩(wěn)健性之間的沖突，選擇基于Pareto非劣解概念的NCGA多目標(biāo)遺傳算法做氣動外形優(yōu)化設(shè)計(jì)，并應(yīng)用RBF神經(jīng)網(wǎng)絡(luò)代理模型技術(shù)解決優(yōu)化過程中大量CFD流場求解的效率問題。最后通過優(yōu)化算例對比了優(yōu)化結(jié)果與初始翼型RAE2822之間不同的性能特征?？梢钥闯?，多目標(biāo)穩(wěn)健設(shè)計(jì)很好地均衡了性能優(yōu)化和穩(wěn)健性要求，在改善氣動性能均值的基礎(chǔ)上又保持了性能的穩(wěn)健性，克服了非設(shè)計(jì)狀態(tài)翼型性能惡化的問題，更符合飛行器實(shí)際飛行中環(huán)境參數(shù)存在波動特性對飛行性能的要求，達(dá)到了翼型穩(wěn)健設(shè)計(jì)的目標(biāo)。

［1］李焦贊，高正紅.氣動設(shè)計(jì)問題中確定性優(yōu)化與穩(wěn)健優(yōu)化的對比研究［J］.航空計(jì)算技術(shù)，2010，40（2）：28－31.

［2］CIPRIAN M，F(xiàn)ATEH N，PEDIRODA V，et al.Multi criteria decision making：a tool to take decisions in multi objective optimization［R］.May 7－10 2007，AIAA paper 2007－2875.

［3］王振國，陳小前，羅文彬，等.飛行器多學(xué)科設(shè)計(jì)優(yōu)化理論與應(yīng)用研究［M］.國防工業(yè)出版社，2006.

［4］丁繼峰，李為吉，張勇，等.基于響應(yīng)面的翼型穩(wěn)健設(shè)計(jì)研究［J］.空氣動力學(xué)學(xué)報(bào)，2007，25（1）：19－22.

［5］SOBIECZKY H.Parametric airfoils and wings［J］.Notes onNumericalFluidMechanics，1998，68：71－88.

［6］蘇偉.基于CFD技術(shù)和代理模型的氣動外形優(yōu)化設(shè)計(jì)方法研究［D］.西安：西北工業(yè)大學(xué)航空學(xué)院，2007.

［7］SU Wei，ZUO Yingtao，GAO Zhenghong.Preliminary aerodynamic shape optimization using genetic algorithm and neural network［R］.Sep.6－8 2006，AIAA paper 2006－7106.

［8］SHINYA W，TOMOLYULD H，MITSUNORI M.Neighborhood cultivation genetic algorithm for multi－objective optimization problems［A］.In the Proceedings of the 4th Asia－Pacific Conference on Simulated Evolution and learning［C］.Singapore：Nanyang Technological University，School of Electrical ＆ Electronic Engineering，2002：198－202.

［9］唐智禮.基于非確定性的翼型魯棒反設(shè)計(jì)方法［J］.南京航空航天大學(xué)學(xué)報(bào)，2007，39（1）：16－21.