基于新型多可信度代理模型的多目標(biāo)優(yōu)化方法

趙歡，高正紅，夏露

西北工業(yè)大學(xué) 航空學(xué)院，西安 710072

現(xiàn)代飛行器外形設(shè)計(jì)需要考慮復(fù)雜的飛行條件/作戰(zhàn)環(huán)境、多個(gè)飛行狀態(tài)和高的性能表現(xiàn)，設(shè)計(jì)要求越來(lái)越多且更加苛刻，已成為典型的多目標(biāo)復(fù)雜約束優(yōu)化問(wèn)題。比如高速旋翼翼型優(yōu)化不僅要考慮前飛、懸停和機(jī)動(dòng)等多個(gè)狀態(tài)下升阻特性和力矩特性等10 余項(xiàng)設(shè)計(jì)指標(biāo)要求；同時(shí)隨著飛行速度提高，要考慮的狀態(tài)范圍更寬（Ma=0.2～0.9），設(shè)計(jì)要求變高，各個(gè)指標(biāo)間的矛盾更加突出，設(shè)計(jì)空間也更加復(fù)雜，因此，需要采用復(fù)雜多目標(biāo)優(yōu)化模型。此外，隨著多學(xué)科協(xié)同設(shè)計(jì)開(kāi)始廣泛應(yīng)用于飛行器設(shè)計(jì)，如戰(zhàn)斗機(jī)氣動(dòng)/結(jié)構(gòu)/隱身/音爆等各個(gè)學(xué)科的設(shè)計(jì)要求以及單個(gè)學(xué)科內(nèi)亞、跨、超聲速多個(gè)狀態(tài)點(diǎn)的設(shè)計(jì)要求，會(huì)使得設(shè)計(jì)目標(biāo)和設(shè)計(jì)約束進(jìn)一步增多，面臨高維多目標(biāo)設(shè)計(jì)難題。目前以NSGA-Ⅲ為代表的經(jīng)典多目標(biāo)優(yōu)化算法僅對(duì)2～3 個(gè)目標(biāo)具有良好的搜索能力[1］，過(guò)多的目標(biāo)會(huì)導(dǎo)致算法的收斂性和多樣性急劇下降。這是因?yàn)槭紫入S著目標(biāo)個(gè)數(shù)m的增加，近似整個(gè)Pareto 前緣點(diǎn)的個(gè)數(shù)呈指數(shù)增加（m-1 維超曲面），巨大數(shù)量的Pareto 解導(dǎo)致高的計(jì)算花費(fèi)。其次非支配解比例指數(shù)增加（非支配解比例為1-2/2m），幾乎所有的粒子都變成非支配的，過(guò)高的非支配解比例會(huì)使得目前的多目標(biāo)進(jìn)化算法很難選擇具有優(yōu)勢(shì)的個(gè)體，種群的進(jìn)化壓力下降，甚至變?yōu)殡S機(jī)搜索算法，收斂性嚴(yán)重變差，而隨著目標(biāo)個(gè)數(shù)增加，種群數(shù)增大，使得進(jìn)化算法選擇非支配解和保持多樣性的計(jì)算復(fù)雜度劇烈增加[2］。最后，大于3 個(gè)目標(biāo)的Pareto 前緣可視化問(wèn)題也給設(shè)計(jì)決策者帶來(lái)了極大困難，選擇折衷解需要更多的時(shí)間。這類問(wèn)題也被統(tǒng)稱為高維多目標(biāo)優(yōu)化問(wèn)題。目前流行的多目標(biāo)算法對(duì)于3 個(gè)及3 個(gè)以上目標(biāo)問(wèn)題的擴(kuò)展性降低，即便是對(duì)于簡(jiǎn)單的Pareto 面。

自從2003 年進(jìn)化多目標(biāo)國(guó)際會(huì)議[3］召開(kāi)后，許多研究者對(duì)高維多目標(biāo)問(wèn)題進(jìn)行了深入研究，并提出了一系列可行的解決策略。針對(duì)第1 類問(wèn)題，即高維多目標(biāo)導(dǎo)致現(xiàn)有基于Pareto 支配的多目標(biāo)進(jìn)化算法搜索能力下降，收斂特性變差的問(wèn)題，可以通過(guò)以下方法解決：① 修改Pareto 支配方式以減少非支配解的數(shù)量[4］；② 使用不同的排序方法選擇非支配解[5-6］；③ 新的適應(yīng)值評(píng)價(jià)機(jī)理，而不是傳統(tǒng)的Pareto 支配[7］；④ 評(píng)價(jià)指標(biāo)方法的使用，比如超體積評(píng)價(jià)指標(biāo)[7］，但隨著目標(biāo)增加，計(jì)算復(fù)雜度指數(shù)增加；⑤ 將多目標(biāo)問(wèn)題轉(zhuǎn)化為大量的標(biāo)量化函數(shù)優(yōu)化問(wèn)題，例如使用權(quán)重總和方法、切比雪夫度量和法向邊界交叉（NBI），以及使用分解基于多目標(biāo)優(yōu)化算法（Multi-Objective Evolutionary Algorithm Based on Decomposition， MOEA/D）等將多目標(biāo)問(wèn)題轉(zhuǎn)化為多個(gè)單目標(biāo)優(yōu)化問(wèn)題[8-10］，然而生成真實(shí)Pareto前緣均勻分布或偏好分布點(diǎn)的權(quán)重向量卻很難選擇。針對(duì)第2 類問(wèn)題，即Pareto 前緣解集數(shù)量指數(shù)增加的難題，可以通過(guò)以下方法解決：① 引入偏好信息，使多目標(biāo)搜索集中在Pareto 前緣的一個(gè)小區(qū)域[11］；② 目標(biāo)降維方法，通過(guò)分析目標(biāo)之間的相關(guān)性，去掉冗余目標(biāo)，僅優(yōu)化目標(biāo)中最沖突的1 組以使其產(chǎn)生和原始問(wèn)題相同的Pareto前緣解。目前流行的一些目標(biāo)降維方法如下：① 支配關(guān)系守恒方法，即移除某個(gè)目標(biāo)后支配關(guān)系保持不變[12］；② 無(wú)監(jiān)督特征選擇方法[13］，通過(guò)測(cè)量目標(biāo)之間的相關(guān)性，將距離更遠(yuǎn)的目標(biāo)列為更沖突的一些；③ 主成分分析方法[14］，通過(guò)線性主 成 分 分 析（Principal Components Analysis， PCA）降維去掉冗余的目標(biāo)；④ Pareto 角搜索進(jìn)化[15］算法和幾何投影法等。目標(biāo)降維方法不僅能解決以上2 類問(wèn)題，還能解決可視化難題，因此是目前高維多目標(biāo)問(wèn)題的主要解決途徑。

盡管上述許多算法已經(jīng)在試圖解決高維多目標(biāo)優(yōu)化問(wèn)題，但有許多缺乏在實(shí)際工程問(wèn)題中的驗(yàn)證。其中，結(jié)合PCA 降維分析和多目標(biāo)進(jìn)化算法已經(jīng)得到了許多設(shè)計(jì)者的使用，并取得了良好的效果，比如超臨界翼型、旋翼翼型以及飛翼翼型多學(xué)科多目標(biāo)優(yōu)化設(shè)計(jì)等[16-17］。然而，目標(biāo)降維方法僅當(dāng)目標(biāo)中存在較強(qiáng)相關(guān)性的目標(biāo)，即多個(gè)冗余目標(biāo)時(shí)，才能起作用。并且一些線性降維方法如PCA 主要針對(duì)二階及二階以下相關(guān)性的目標(biāo)維數(shù)減少，缺少對(duì)更高階以及一般問(wèn)題的適應(yīng)性。Saxena 等[18］引入了PCA 和非線性主成分分析（Nonlinear Principal Component Analysis，NPCA）結(jié)合的降維方法，在一些測(cè)試函數(shù)上，相對(duì)于支配關(guān)系守恒降維方法，計(jì)算復(fù)雜度更低。可見(jiàn)，高維多目標(biāo)優(yōu)化問(wèn)題，即如何高效地找出高維空間中真實(shí)Pareto 最優(yōu)解，依然困擾著設(shè)計(jì)者，尤其對(duì)于包含很少甚至不含冗余目標(biāo)或相關(guān)性強(qiáng)的目標(biāo)集合，困難更大。

1 剛性共軸雙旋翼直升機(jī)槳葉翼型高維多目標(biāo)優(yōu)化難題

直升機(jī)需要兼顧前飛、懸停和機(jī)動(dòng)等各種性能要求。對(duì)單旋翼直升機(jī)而言，前行槳葉和后行槳葉需要同時(shí)提供拉力并使其合力作用于轉(zhuǎn)軸處，以維持直升機(jī)機(jī)身平衡，因此需要旋翼翼型兼顧低速高升力和高速低阻力特性。而對(duì)于剛性共軸雙旋翼直升機(jī)，在飛行時(shí)存在2 幅旋翼的對(duì)轉(zhuǎn)運(yùn)動(dòng)，槳盤兩側(cè)同時(shí)存在前行槳葉，合力自然作用于轉(zhuǎn)軸處，前行槳葉的高動(dòng)壓保證了足夠的拉力，后行槳葉可不用提供拉力；在直升機(jī)高速前飛時(shí)，可將后行槳葉卸載以避免大范圍反流區(qū)引起的阻力、噪聲激增等問(wèn)題，因此此類直升機(jī)可以大幅提高前飛速度。前飛速度的提高，使槳葉外側(cè)段翼型需要保持更高的阻力發(fā)散馬赫數(shù)和高速低阻特性，相對(duì)于單旋翼翼型，設(shè)計(jì)要求更加嚴(yán)苛。根據(jù)需求，此類高速旋翼翼型槳尖（0.9R～1.0R 處，R 為槳葉半徑）性能要求如下：

1）高速前行時(shí)，槳尖翼型工作狀態(tài)馬赫數(shù)超過(guò)0.87，波阻急劇增加，前飛需用功率增加，設(shè)計(jì)翼型應(yīng)具有高的零升阻力發(fā)散馬赫數(shù)，盡可能低的零升阻力系數(shù)，以及維持低的零升俯仰力矩系數(shù)。槳尖載荷要求可適當(dāng)降低。

2） 在 跨 聲 速 狀 態(tài) 范 圍 內(nèi)(0.8 ≤Ma ≤Madd0，Madd0為零升阻力發(fā)散馬赫數(shù))，激波位置和強(qiáng)度隨馬赫數(shù)等狀態(tài)變化非常敏感，設(shè)計(jì)翼型在此范圍內(nèi)應(yīng)該具有穩(wěn)健的氣動(dòng)表現(xiàn)。

3）懸停狀態(tài)下，要求翼型在中等馬赫數(shù)下提供高升阻比以減小需用功率。例如在設(shè)計(jì)馬赫數(shù)范圍內(nèi)(Ma=0.4～0.6)，要求具有高的升阻比。

4）機(jī)動(dòng)飛行時(shí)需要高的載荷系數(shù)，尤其是后行槳葉一側(cè)，在大迎角范圍內(nèi)極容易出現(xiàn)失速顫振現(xiàn)象以及升力特性的突然下跌，致使力矩變大，飛行包線減小。因此，在低馬赫數(shù)范圍內(nèi)(Ma=0.2～0.4），應(yīng)具有高的最大升力系數(shù)和最大升阻比，以提高其機(jī)動(dòng)性能。

5）在所有馬赫數(shù)范圍或整個(gè)飛行包線范圍內(nèi)，各剖面設(shè)計(jì)翼型應(yīng)具有極低的俯仰力矩系數(shù)以減小槳葉扭轉(zhuǎn)及控制系統(tǒng)的操縱載荷。

6）一些研究表明，在槳葉表面近似0.45R～1.00R 區(qū)域內(nèi)，層流區(qū)能穩(wěn)定存在，并且可以顯著提高升阻比。因此，這個(gè)區(qū)域內(nèi)翼型可以使用自然層流超臨界翼型設(shè)計(jì)。

以上這些指標(biāo)要求達(dá)10 余項(xiàng)之多，并通常是相互矛盾和制約的，面臨高維多目標(biāo)和復(fù)雜約束優(yōu)化難題。因此，發(fā)展高效的高維多目標(biāo)設(shè)計(jì)模型和設(shè)計(jì)方法，并設(shè)計(jì)滿足嚴(yán)苛設(shè)計(jì)要求的剛性共軸雙旋翼直升機(jī)翼型非常重要。

高維多目標(biāo)優(yōu)化問(wèn)題可以數(shù)學(xué)定義為

式中：定義了m 個(gè)目標(biāo){ f1(X)，f2(X)，…，fm(X) }，p 個(gè)不等式約束gi(X)以及q 個(gè)等式約束hj(X)的優(yōu)化問(wèn)題。代表n維設(shè)計(jì)變量。XL和XU分別為設(shè)計(jì)變量下邊界和上邊界。實(shí)際優(yōu)化中等式約束hj總可以轉(zhuǎn)化為不等式約束，即

式（2）為僅考慮設(shè)計(jì)點(diǎn)氣動(dòng)特性的確定性優(yōu)化設(shè)計(jì)模型，無(wú)法兼顧非設(shè)計(jì)點(diǎn)的氣動(dòng)特性，如來(lái)流馬赫數(shù)變化時(shí)的氣動(dòng)特性。

如果考慮飛行狀態(tài)等各種不確定源對(duì)翼型氣動(dòng)特性的影響，將式（2）模型直接修改為

式中：Ξ=[ ξ1，ξ2，…，ξi]表示各種可能的不確定源，如馬赫數(shù)和升力系數(shù)。本文假設(shè)所有的不確定源ξi為獨(dú)立的連續(xù)概率分布變量。很顯然，由于Ξ 為滿足一定概率分布的不確定變量，對(duì)式（3）模型直接優(yōu)化很難獲得在Ξ 可能的各種變化值后仍然保持穩(wěn)健氣動(dòng)特性的翼型。因此，式（3）模型很難考慮各種不確定源的影響，也無(wú)法滿足工程使用需求。為了使剛性共軸雙旋翼直升機(jī)翼型在高速狀態(tài)下保持穩(wěn)健的氣動(dòng)特性[19-21］，基于連續(xù)概率理論，將式（3）模型轉(zhuǎn)變?yōu)榉€(wěn)健多目標(biāo)優(yōu)化模型，即

式中：μfi(X，Ξ)和σfi(X，Ξ)分別代表目標(biāo)fi(i=1，2，…，m)基于不確定變量Ξ 影響的均值和方差；μgi(X，Ξ) 和σgi(X，Ξ) 分 別 代 約 束gi(i=1，2，…，p)基于不確定變量Ξ 影響的均值和標(biāo)準(zhǔn)差；μhj(X，Ξ)則代表約束hj( j=1，2，…，q)基于不確定變量Ξ 影響的均值；參數(shù)k 取值越大則約束可靠性越強(qiáng)，反之亦然。式（4）模型為確定性優(yōu)化模型式（2）在考慮多源不確定性影響后建立的穩(wěn)健優(yōu)化設(shè)計(jì)模型。在得到連續(xù)隨機(jī)變量Ξ的概率密度函數(shù)ρ(Ξ)后，可以通過(guò)積分方式得到目標(biāo)函數(shù)和約束函數(shù)的統(tǒng)計(jì)特性，目標(biāo)函數(shù)fi的均值和方差為

式中：ΞU和ΞL分別代表隨機(jī)變量Ξ 的變化上限和下限。

相對(duì)于傳統(tǒng)的確定性多目標(biāo)優(yōu)化，穩(wěn)健多目標(biāo)優(yōu)化引入了各個(gè)目標(biāo)表現(xiàn)和穩(wěn)健性，以及各個(gè)約束的可行穩(wěn)健性測(cè)量。因而目標(biāo)個(gè)數(shù)成倍增加，優(yōu)化可行域變窄且更加復(fù)雜，函數(shù)評(píng)估數(shù)量增加。目前一些經(jīng)典的多目標(biāo)算法如NSGA-Ⅲ等對(duì)于高維多目標(biāo)（目標(biāo)個(gè)數(shù)超過(guò)3 個(gè)）問(wèn)題，其搜索和收斂能力均下降，尤其一些冗余的目標(biāo)會(huì)影響算法的收斂效率和真實(shí)表現(xiàn)。因此，通過(guò)有效的目標(biāo)降維方法結(jié)合目前經(jīng)典的多目標(biāo)算法，將能有效解決高維多目標(biāo)難題。同時(shí)，基于前期發(fā)展的多可信度代理建模技術(shù)（Adaptive Multi-Fidelity Polynomial Chaos-Kriging， AMF-PCK），本文提出了高效的高維多目標(biāo)穩(wěn)健優(yōu)化平臺(tái)，以滿足工程使用。

2 自適應(yīng)多可信度多項(xiàng)式混沌-Kriging 代理模型

相 對(duì) 于Universal Kriging、Blind Kriging 和Dynamic Kriging， LAR-PC-Kriging （Least Angle Regression -Polynomial Chaos Expansion -Kriging）不僅具有更好的自適應(yīng)基函數(shù)選擇策略，同時(shí)由于多項(xiàng)式混沌展開(kāi)提供了完備的正交基，理論上多項(xiàng)式混沌對(duì)不同的問(wèn)題能提供更好的近似能力和更強(qiáng)的適應(yīng)性。然而，由于目前自適應(yīng)基函數(shù)選擇方法受樣本數(shù)的影響[22］，當(dāng)有限甚至非常少的樣本可利用時(shí)，LAR-PC-Kriging模型不能正確選擇一組最優(yōu)的多項(xiàng)式基函數(shù)，導(dǎo)致其近似能力可能比Universal Kriging 更差[23］。結(jié)合多可信度建模技術(shù)，本文提出了一種更加高效的代理建模途徑，即自適應(yīng)多可信度多項(xiàng)式混沌-Kriging （AMF-PCK）代理模型[24］。

AMF-PCK 在任意點(diǎn)X 的預(yù)測(cè)值為

式中：α0為縮放因子；yc(X)是對(duì)低可信度多項(xiàng)式混沌展開(kāi)（Polynomial Chaos Expansion， PCE）矯正修正項(xiàng)；C(X)為定義的第2 個(gè)高斯過(guò)程；zh(X )為零均值的靜態(tài)高斯隨機(jī)過(guò)程，其協(xié)方差為θh=[ θh1， θh2，…，θhn]T為 高 斯過(guò) 程 zh(X )的 模 型超參數(shù)。空間相關(guān)函數(shù)R(X(i)，X(j)；θh)為超參數(shù)θh和qh以及高可信度訓(xùn)練樣本Sh的函數(shù)。yc(X)的表達(dá)式為

其中：αci和Ac分別表示校正展開(kāi)多項(xiàng)式項(xiàng)系數(shù)和相應(yīng)的多項(xiàng)式項(xiàng)索引集合；y?l(X)為低可信度多項(xiàng)式混沌-Kriging（PC-Kriging）模型預(yù)測(cè)器；zl(X )為低可信度PC-Kriging 模型yl(X)的靜態(tài)高斯隨機(jī)過(guò)程；ψi(X)和ψci(X)分別為yl(X)的多項(xiàng)式混沌項(xiàng)和矯正展開(kāi)多項(xiàng)式混沌項(xiàng)；βi為低可信度多項(xiàng)式混沌項(xiàng)系數(shù)。多項(xiàng)式混沌矯正項(xiàng)集合Ac的勢(shì)應(yīng)該小于低可信度稀疏多項(xiàng)式集合A 的勢(shì)，即Ac?A，其中A={i1，i2，…，iM}為低可信度PCKriging 模型yl(X)中多項(xiàng)式項(xiàng)索引集合。通過(guò)使用低可信度樣本構(gòu)建基于集合Ac的勢(shì)應(yīng)該小于低可信度稀疏多項(xiàng)式集合A 的勢(shì)，即Ac?A。

通過(guò)使用基于最小角回歸的留一交叉驗(yàn)證（Leave One Out Cross-Validation-based Least Angle Regression， LOOCV-LAR）自 適 應(yīng) 基 函數(shù)選擇方法得到M 個(gè)最重要的基函數(shù)后，再使用一般高斯過(guò)程回歸分析可以得到低可信度PC-Kriging 預(yù)測(cè)器為[25］

式中：βij為基函數(shù)ψij(X)對(duì)應(yīng)的系數(shù)；rl、Rl、Ψ以及β分別為低可信度PC-Kriging 模型中自相關(guān)函數(shù)向量、樣本相關(guān)函數(shù)矩陣、多項(xiàng)式基函數(shù)矩陣以及多項(xiàng)式系數(shù)向量；Yl為低可信度響應(yīng)向量。通過(guò)代入式（8）～式（10）得

其中：α0以及{αci}是需要被精確訓(xùn)練的縮放因子。子集Ac應(yīng)該通過(guò)一個(gè)合適的標(biāo)準(zhǔn)進(jìn)行選擇。

本文利用LOOCV-LAR 算法選擇最相關(guān)的稀疏多項(xiàng)式基函數(shù)作為PC-Kriging 模型的趨勢(shì)函數(shù)，即PC-Kriging 模型中基函數(shù)已經(jīng)按照重要性進(jìn)行排序。越靠前的多項(xiàng)式項(xiàng)對(duì)輸出貢獻(xiàn)越大也最重要。因此，通過(guò)利用非常少的高可信度樣本僅僅矯正排序靠前的多項(xiàng)式系數(shù)，就能獲得對(duì)低可信度PC-Kriging 的有效矯正。為此，本文設(shè)計(jì)了一種新的基于LOOCV 的自適應(yīng)基函數(shù)選擇方法，即在選擇最相關(guān)的M項(xiàng)低可信度多項(xiàng)式項(xiàng)之后，再選擇最優(yōu)的矯正展開(kāi)項(xiàng)集 合{ψi1(X)，ψi2(X)，…，ψiMc(X) }(Mc≤M)，以建立最優(yōu)的AMF-PCK 模型，這個(gè)算法被定義為基 于 嵌 套 LOOCV 的 LAR 算 法（Nested LOOCV-LAR），可參考圖1 和文獻(xiàn)[24］。其中圖1 為自適應(yīng)選擇最相關(guān)的低可信度多項(xiàng)式集合和最優(yōu)的矯正多項(xiàng)式集合圖中：ML為候選低可信度多項(xiàng)式基函數(shù)項(xiàng)數(shù)；M*和M*c分別為選擇的最優(yōu)低可信度多項(xiàng)式項(xiàng)數(shù)和矯正多項(xiàng)式項(xiàng)數(shù)。

圖1 自適應(yīng)選擇最相關(guān)的低可信度多項(xiàng)式集合和最優(yōu)的矯正多項(xiàng)式集合Fig.1 Adaptive selection of the most relevant lowfidelity polynomial set and optimal correction polynomial set

在式（11）中的各項(xiàng)被確定后，通過(guò)標(biāo)準(zhǔn)Universal Kriging 擬合過(guò)程即可得到AMF-PCK 在未知點(diǎn)的預(yù)測(cè)模型。同Universal Kriging 模型一樣，AMF-PCK 滿足最小化均方誤差和無(wú)偏估計(jì)假設(shè)（Best Linear Unbiased Predictor，BLUP）。AMF-PCK 預(yù)測(cè)模型在未知點(diǎn)X的預(yù)測(cè)值為

式中：Mc(Mc≤M)表示校正展開(kāi)項(xiàng)中多項(xiàng)式項(xiàng)數(shù)量。

縮放因子以及校正展開(kāi)多項(xiàng)式系數(shù)向量α=[α0，αc1，αc2，αcMc]T可通過(guò)廣義最小 二乘方 法估計(jì)得到，即

式中：矩陣F中第1 列和其他列是通過(guò)順序代入高可信度觀測(cè)點(diǎn)集合Sh到(X)和ψci(X) (i=1，2，…，Mc)中計(jì)算得到。相關(guān)函數(shù)Rh中未知超參數(shù)θh和qh通過(guò)最大化似然函數(shù)得到。而AMF-PCK 預(yù)測(cè)模型在測(cè)試點(diǎn)X的估計(jì)均方誤差（Mean-Squared Error， MSE）為

式（14）提供了一個(gè)非常有用的MSE 估計(jì)，同Kriging（或高斯過(guò)程）模型本身的優(yōu)勢(shì)一樣，合理的MSE 估計(jì)能幫助AMF-PCK 代理模型自適應(yīng)地修訂預(yù)測(cè)精度以及有利于基于（多可信度）期望改進(jìn)（Expected Improvement， EI）的全局優(yōu)化。本文使用粒子群算法與基于梯度的稀疏非線性優(yōu)化器（Sparse Nonlinear Optimizer， SNOPT）相結(jié)合，以得到全局最優(yōu)的超參數(shù)θl、 ql、 θh和qh。其中θl、ql代表低可信度PC-Kriging 模型的超參數(shù)。優(yōu)化過(guò)程中優(yōu)化范圍定義為[10-4θini，104θini]，θini是超參數(shù)初始值。首先使用粒子群優(yōu)化算法搜索全部設(shè)計(jì)空間找出可能的全局最優(yōu)所在區(qū)域，然后使用基于梯度的SNOPT 利用粒子群優(yōu)化方法提供的初值進(jìn)行第2次搜索直至收斂。

3 基于新型AMF-PCK 模型的多目標(biāo)優(yōu)化方法

3.1 非線性目標(biāo)降維策略

以主成分分析（Principal Component Analysis， PCA）（ 本 征 正 交 分 解、奇 異 值 分 解、Karhunen-Loève 轉(zhuǎn)換等）為代表的線性目標(biāo)降維方法通過(guò)投影當(dāng)前高維數(shù)據(jù)到相關(guān)矩陣的特征向量方向，以找到數(shù)據(jù)的本質(zhì)結(jié)構(gòu)或主方向。PCA移除了當(dāng)前數(shù)據(jù)的二階相關(guān)性。然而對(duì)于更加復(fù)雜的相關(guān)結(jié)構(gòu)，如多峰高斯分布或非高斯分布等，PCA 難以找到最重要的主特征方向，使得降維失敗。一些文獻(xiàn)已經(jīng)提到了類似問(wèn)題[1，26］，并通過(guò)引入非線性降維途徑試圖解決，比如Correntropy-PCA（C-PCA）、最大方差展開(kāi)（Maximum Variance Unfolding， MVU）等。本文引入基于MVU的非線性目標(biāo)降維方法[27］，以處理一般的工程問(wèn)題。MVU 學(xué)習(xí)高維數(shù)據(jù)的相似性，在映射前后保持與臨近點(diǎn)局部歐式距離以及角度相等，進(jìn)而能有效學(xué)習(xí)高維數(shù)據(jù)中的低維流形，反映流形的維數(shù)。MVU 可視為自動(dòng)訓(xùn)練核矩陣K 的核主成分 分 析（Kernel Principal Component Analysis， KPCA）的一種，然而相對(duì)于KPCA，MVU 不包含待訓(xùn)練參數(shù)以及能保持映射前后的鄰接關(guān)系不變。基于MVU 的目標(biāo)降維算法為

1）首先通過(guò)經(jīng)典的多目標(biāo)進(jìn)化算法（Multi-Objective Evolutionary Algorithms， MOEAs），如基于分解的多目標(biāo)進(jìn)化算法對(duì)高維問(wèn)題進(jìn)行優(yōu)化搜索一些代數(shù)后得到近似Pareto 解集。假設(shè)每代種群數(shù)量為Np，目標(biāo)個(gè)數(shù)為M，得到的搜索多目標(biāo)結(jié)果集合為F=[ f1，f2，…，fNp]，其中列向量fi=[ fi1，fi2，…，fiM]T∈RM表示這一代第i 個(gè)粒子各個(gè)目標(biāo)的適應(yīng)值。

式中：ωij表示fi與fj是否是?-近鄰，如果是，則為1，否則為0。? 代表最近鄰居的個(gè)數(shù)，文獻(xiàn)[28］指 出 ? 通 常 取 為近 鄰 矩 陣Ω=首先通過(guò)輸入數(shù)據(jù)計(jì)算得出。上述優(yōu)化約束保證了高維空間數(shù)據(jù)映射到低維空間時(shí)近鄰的歐式距離不變。上述優(yōu)化問(wèn)題通常非常耗時(shí)，且不容易得到最優(yōu)解。通過(guò)學(xué)習(xí)KPCA 方法，引入核矩陣，則上述優(yōu)化問(wèn)題可轉(zhuǎn)化為

上述凸半正定優(yōu)化問(wèn)題可通過(guò)SeDuMi[29］求解。如果使用PCA 線性降維，則可直接計(jì)算協(xié)方差矩陣為

3）對(duì)得到的核矩陣K 或者R 進(jìn)行特征分解，即K=VΛVT。 其 中V 為 特 征 向 量 矩 陣V=[v1，v2，…，vM]，vi=[vi1，vi2，…，viM]T和相應(yīng)的特征 值 Λ=diag(λ1，λ2，…，λM) (λ1≥λ2≥…≥λM≥0)。

5）定義降維后的初始目標(biāo)集合F=?，對(duì)于 第1 個(gè) 重 大 方 向vj(j=1，2，…，m)：① 如 果fi(i=1，2，…，M)目標(biāo)對(duì)v1方向的貢獻(xiàn)中存在符號(hào)相反的一些，則選擇最正和最負(fù)貢獻(xiàn)對(duì)應(yīng)的2 個(gè)目標(biāo)進(jìn)入集合F；② 如果fi(i=1，2，…，M)目標(biāo)對(duì)vj(j≥2)方向的貢獻(xiàn)中存在符號(hào)相反的一些，記最正貢獻(xiàn)為vpj以及最負(fù)貢獻(xiàn)為vnj，如果|vnj|≥vpj≥θt|vnj|則 選 擇vpj和vnj對(duì) 應(yīng) 的2 個(gè) 方 向均進(jìn)入F，如果vpj＜θt|vnj|則選擇vnj對(duì)應(yīng)的方向進(jìn)入F，如果vpj≥|vnj|≥θtvpj則選擇vpj和vnj對(duì)應(yīng)的2 個(gè)方向均進(jìn)入F，如果|vnj|＜θtvpj則選擇vpj對(duì)應(yīng)的方向進(jìn)入F，本文取系數(shù)θt=0.8；③ 如果所有fi(i=1，2，…，M)對(duì)vj方向的貢獻(xiàn)符號(hào)相同，則選擇貢獻(xiàn)最大的2 個(gè)目標(biāo)進(jìn)入集合F，F(xiàn) 應(yīng)該滿足F ?F，記mv=Card(F )。

6）使用相關(guān)性分析對(duì)F 中的相關(guān)目標(biāo)（貢獻(xiàn)等同的）進(jìn)一步識(shí)別和縮減。將入選F 中的目標(biāo)按列按順序放在矩陣H 中，然后計(jì)算相關(guān)矩陣R=HHT/M。對(duì)于fi∈F 以及fj∈F，如果滿足sign(Rik)=sign(Rjk) (k=1，2，…，M)以及Rij≥Rθ，則說(shuō)明fi與fj的貢獻(xiàn)方向相同，則選擇其中對(duì)主導(dǎo)方向貢獻(xiàn)更大的一個(gè)，即如果ci≥cj則選擇fi，從F 中 剔 除fj，反 之 亦 然。其 中Rij表示R的 元素。Rθ的選擇與目標(biāo)間的冗余度有關(guān)，目標(biāo)越冗余Rθ越大，反之越小。文獻(xiàn)[28］使用Rθ=1.0-e1(1.0-M0.954/M)，其中M0.954指前M0.954個(gè)主成分占據(jù)了0.954 的總能量。

7）用新得到的集合F 替換舊集合F，重復(fù)步驟2～步驟6 直到連續(xù)2 次得到的集合F 一樣，則停止迭代，輸出降維目標(biāo)集合F*。

8）使用MOEAs 算法對(duì)新目標(biāo)集合F*選擇的目標(biāo)進(jìn)行搜索，得到Pareto 最優(yōu)解集合，折衷選擇最合適的設(shè)計(jì)，進(jìn)行評(píng)估分析。

3.2 變可信度偽期望改進(jìn)矩陣并行樣本填充方法

在外循環(huán)中，每次多目標(biāo)優(yōu)化迭代結(jié)束后，需要選擇新的樣本填充點(diǎn)。設(shè)計(jì)變量作為優(yōu)化粒子在不斷進(jìn)化，新樣本的設(shè)計(jì)變量值X應(yīng)該平衡全局探索與局部搜索的計(jì)算花費(fèi)，如多目標(biāo)期望改進(jìn)（EI）準(zhǔn)則，即多目標(biāo)高效全局優(yōu)化（Efficieng Global Optimization， EGO）。同時(shí)，每次迭代中，為了提高多可信度AMF-PCK 模型在特定空間內(nèi)的預(yù)測(cè)精度，新加點(diǎn)的可信度水平也需要同時(shí)被確定。多目標(biāo)EGO 準(zhǔn)則雖然從很早開(kāi)始已經(jīng)獲得了廣泛研究，比如基于歐式距離改進(jìn)的EI準(zhǔn)則[30］，基于最大最小距離改進(jìn)的EI[31］，基于超體積改進(jìn)的EI[32］，基于切比雪夫度量改進(jìn)的EI[33］等。然而，這些方法應(yīng)用到高維多目標(biāo)(m≥3)加點(diǎn)時(shí)，多目標(biāo)EI的計(jì)算是一個(gè)典型的高維積分問(wèn)題，計(jì)算花費(fèi)隨目標(biāo)個(gè)數(shù)增多劇烈增加。并且其每次只能添加1個(gè)樣本，面對(duì)高維多目標(biāo)問(wèn)題(m≥3)時(shí)，其有效性仍待驗(yàn)證[34］。Keane[30］提出使用蒙特卡洛模擬（MCS）近似高維積分，計(jì)算準(zhǔn)確性與所使用樣本數(shù)密切相關(guān)，其計(jì)算花費(fèi)仍然難以承擔(dān)。Singh 等[15］提出了一種快速計(jì)算基于超體積改進(jìn)EI和基于歐式距離改進(jìn)EI的方法，在較低維度時(shí)，其計(jì)算花費(fèi)有明顯降低。因此，針對(duì)以上問(wèn)題，提出了一種變可信度偽期望改進(jìn)矩陣（Variable-Fidelity Pseudo Expected Improvement Matrix， VF-PEIM）并行樣本填充方法。期望改進(jìn)矩陣（Expected Improvement Matrix， EIM）標(biāo) 準(zhǔn) 是 由Zhan 等[34］提 出 的一種高效多目標(biāo)EI 填充方法，其通過(guò)將上述傳統(tǒng)的多目標(biāo)高維EI函數(shù)積分問(wèn)題轉(zhuǎn)化為多個(gè)一維EI積分之和，每個(gè)EI函數(shù)即EIji代表針對(duì)第j個(gè)Pareto前緣點(diǎn)，第i個(gè)目標(biāo)的期望改進(jìn)。二維EIM 如圖2所示。本文提出使用VF-EI定義EIji為

圖2 二維EIMFig.2 Two-dimensional expected improvement matrix

式中：l 代表可信度水平。在定義了式（17）后，可將歐式距離改進(jìn)準(zhǔn)則[30］中的高維積分問(wèn)題轉(zhuǎn)化為多個(gè)一維積分求和，即

通過(guò)對(duì)式（19）單目標(biāo)優(yōu)化搜索可獲得1 個(gè)最有潛力改進(jìn)點(diǎn)位置，即

然而，式（20）每次只能添加1 個(gè)樣本點(diǎn)，這對(duì)于多目標(biāo)優(yōu)化問(wèn)題通常是非常低效的。為此，引入偽期望改進(jìn)矩陣（Pseudo Expected Improvement Matrix， PEIM）填充準(zhǔn)則[35］，在每個(gè)加點(diǎn)階段，通過(guò)迭代優(yōu)化選取若干個(gè)樣本點(diǎn)。其中，第1 個(gè)加點(diǎn)位置通過(guò)式（20）得到，記為{ X1，l1}。然后，計(jì)算PEI 準(zhǔn)則為

第Nc個(gè)樣本點(diǎn)位置為

式（21）～式（24）的優(yōu)化均是單目標(biāo)優(yōu)化問(wèn)題，因而相對(duì)多目標(biāo)優(yōu)化效率很高。VF-PEIM 通過(guò)定義每個(gè)修訂周期合適的加點(diǎn)數(shù)，即可有效地獲得全局最優(yōu)解，從而兼顧了多目標(biāo)優(yōu)化效率和搜索能力。

3.3 高維多目標(biāo)優(yōu)化框架

針對(duì)式（4）的多目標(biāo)穩(wěn)健優(yōu)化設(shè)計(jì)模型，每個(gè)目標(biāo)在進(jìn)化迭代過(guò)程中均需要進(jìn)行穩(wěn)健性σfi(X)和平均表現(xiàn)μfi(X)評(píng)估，這通常是非常耗時(shí)的。同時(shí)，因?yàn)橛械哪繕?biāo)通常在不同的設(shè)計(jì)狀態(tài)（針對(duì)同一設(shè)計(jì)狀態(tài)的阻力值和力矩值只需要1 次評(píng)估，而不同狀態(tài)下的阻力值則需要分別評(píng)估），也許會(huì)卷入不同的不確定因素，因此會(huì)包含多個(gè)樣本設(shè)計(jì)，即需要建立針對(duì)多個(gè)目標(biāo)和約束的多個(gè)代理模型。為此，本文基于AMF-PCK 模型建立同時(shí)考慮設(shè)計(jì)變量與隨機(jī)變量的聯(lián)合代理模型。利用3.1 節(jié)引入的非線性目標(biāo)降維策略[36］，建立了針對(duì)復(fù)雜工程優(yōu)化問(wèn)題的高維多目標(biāo)穩(wěn)健優(yōu)化框架，其流程如圖3 所示。其包含4 個(gè)主要步驟：① 針對(duì)式（4）的優(yōu)化模型定義設(shè)計(jì)變量和隨機(jī)變量，通過(guò)有效的抽樣方法獲取一定的樣本并評(píng)估其表現(xiàn)，建立針對(duì)不同目標(biāo)的代理模型；② 使用MOEA/D 多目標(biāo)優(yōu)化算法對(duì)初始包含多個(gè)目標(biāo)的設(shè)計(jì)問(wèn)題進(jìn)行優(yōu)化設(shè)計(jì)，并獲得收斂解；③ 使用降維方法對(duì)多目標(biāo)進(jìn)行降維分析，獲得降維后的目標(biāo)集合；③ 對(duì)降維后的目標(biāo)集合進(jìn)行重新優(yōu)化搜索，獲得收斂解；④ 分析Pareto 解，并利用提出的變可信度偽期望改進(jìn)矩陣（VFPEIM）并行樣本填充策略填充新的樣本，構(gòu)建新的代理模型，重復(fù)步驟②～步驟④，直到收斂。

圖3 基于AMF-PCK 的多目標(biāo)優(yōu)化流程Fig.3 Many-objective optimization flow chart based on AMF-PCK

3.4 數(shù)值算例測(cè)試

采用經(jīng)典的DTLZ5 函數(shù)測(cè)試所提出方法的優(yōu)化效果，具體函數(shù)表達(dá)式可參考文獻(xiàn)[1，26］。測(cè)試算例均使用2 層可信度函數(shù)構(gòu)建AMF-PCK 代理模型，即l=1。多目標(biāo)優(yōu)化算法均使用基于分解的多目標(biāo)進(jìn)化算法即MOEA/D，均使用基于懲罰的邊界交叉（Penalty-based boundary Intersection， PBI）聚合方法，鄰居個(gè)數(shù)均為20，粒子數(shù)均為300，進(jìn)化代數(shù)均為500。在并行加點(diǎn)子優(yōu)化中，均使用單目標(biāo)加強(qiáng)學(xué)習(xí)粒子群優(yōu)化（Comprehensive Learning Particle Swarm Optimizer， CLPSO）算法搜索新的加點(diǎn)位置。針對(duì)VF-PEIM 方法，使用拉丁超立方抽樣（LHS）方法在設(shè)計(jì)空間[0，1]D內(nèi)產(chǎn)生3D個(gè)初始高可信度樣本以及40D個(gè)初始低可信度樣本，D為設(shè)計(jì)變量個(gè)數(shù)。而針對(duì)PEIM 方法，使用Kriging 模型進(jìn)行代理以及使用LHS 方法選擇10D個(gè)高可信度點(diǎn)作為初始樣本。選擇轉(zhuǎn)世代距離（Generational Distance， GD）和反轉(zhuǎn)世代距離（Inverted Generational Distance， IGD）指標(biāo)評(píng)價(jià)2 種方法得到的非支配解集的收斂性和分布性，其中參考解是在真實(shí)的Pareto前緣上均勻地選取10 000個(gè)離散點(diǎn)得到。

圖4 為使用基于PEIM 和VF-PEIM 加點(diǎn)優(yōu)化方法分別執(zhí)行。在10 次迭代后，結(jié)果對(duì)比如圖4 （a）所示，基于AMF-PCK 代理模型的VFPEIM 加點(diǎn)方法在同樣的迭代次數(shù)下獲得更低的IGD 值，即Pareto 解集更好的均勻性和收斂性。并且基于VF-PEIM 方法到10 次迭代時(shí)，IGD 值已經(jīng)基本收斂，而PEIM 方法收斂較慢。比如在10 次迭代后，針對(duì)3 目標(biāo)DTLZ5 函數(shù)，PEIM 加點(diǎn)方法得到的IGD 值為0.047 104，而基于VFPEIM 方法得到的IGD 值為0.020 379。可見(jiàn)，相對(duì)于經(jīng)典PEIM 策略，提出的VF-PEIM 加點(diǎn)策略優(yōu)化效果和效率得到顯著改善。在10 次迭代后，2 種方法收斂的Pareto 解集與真實(shí)Pareto 解集對(duì)比如圖4（b）所示。可以看出，相對(duì)于經(jīng)典的PEIM方法，基于VF-PEIM 方法獲得了與真實(shí)Pareto 解更貼近和分布范圍更廣的解集。原因之一在于，在相同迭代次數(shù)下，AMF-PCK 能提供更準(zhǔn)確的預(yù)測(cè)精度，因而MOEA/D 方法能有效地搜索最優(yōu)Pareto 解集。結(jié)果驗(yàn)證了基于AMF-PCK 的VFPEIM 并行加點(diǎn)優(yōu)化方法的有效性和可靠性。

圖4 使用基于PEIM 和VF-PEIM 加點(diǎn)優(yōu)化方法分別執(zhí)行10 次迭代IGD 和得到的Pareto 解集對(duì)比（DTLZ5 函數(shù)）Fig.4 Comparison of IGD and Pareto solution set using infilling-based PEIM and VF-PEIM methods respectively based on 10 iterations （DTLZ5 function）

為了進(jìn)一步驗(yàn)證提出的方法在更多目標(biāo)優(yōu)化問(wèn)題中的適應(yīng)性，使用目標(biāo)個(gè)數(shù)M={3，4，5，6，7，8，9，10 }的DTLZ5 函 數(shù) 進(jìn) 行 分 析。基于PEIM 和VF-PEIM 的2 種方法所有設(shè)置同上，均執(zhí)行10 次加點(diǎn)迭代。基于PEIM 和VFPEIM 這2 種方法得到的結(jié)果對(duì)比如表1 所示。從表1 可以看出，隨著目標(biāo)個(gè)數(shù)增加，基于PEIM和VF-PEIM 加點(diǎn)優(yōu)化方法得到的GD 和IGD 值在M≥5 后均迅速增加，即Pareto 解集的收斂性和均勻性迅速變差。基于VF-PEIM 的加點(diǎn)優(yōu)化方法表現(xiàn)要優(yōu)于基于PEIM 加點(diǎn)優(yōu)化方法，即各個(gè)目標(biāo)下更低的GD 和IGD 值。

為了解決目標(biāo)個(gè)數(shù)M≥5 后收斂性和多樣性迅速變差的問(wèn)題，在迭代優(yōu)化過(guò)程中引入基于MVU 的非線性目標(biāo)降維策略，詳細(xì)過(guò)程如3.1 節(jié)所述，記為VF-PEIM- MVU 方法。該方法先使用基于VF-PEIM 加點(diǎn)優(yōu)化方法收斂后，對(duì)得到的高維多目標(biāo)（M≥5）Pareto 解集使用目標(biāo)降維方法篩掉冗余目標(biāo)，然后對(duì)降維后的目標(biāo)集合使用基于VE-PEIM 加點(diǎn)優(yōu)化方法進(jìn)行再次優(yōu)化。如圖3 所示執(zhí)行迭代優(yōu)化直至降維目標(biāo)集合和Pareto 解集收斂。從表1 可以看出，基于VFPEIM-MUV 方法在高維多目標(biāo)優(yōu)化問(wèn)題下得到的Pareto 解集GD 和IGD 指標(biāo)均有所改善，即Pareto 解集具有更好的收斂性和分布性，并明顯好于基于VF-PEIM 和PEIM 優(yōu)化加點(diǎn)方法得到的結(jié)果。從而驗(yàn)證了提出的方法在高維多目標(biāo)優(yōu)化問(wèn)題中的適應(yīng)性和有效性，以及其相對(duì)于經(jīng)典的PEIM 多目標(biāo)加點(diǎn)優(yōu)化方法產(chǎn)生的顯著改進(jìn)。

表1 不同代理優(yōu)化算法的Pareto 解集的評(píng)價(jià)指標(biāo)對(duì)比（DTLZ5 函數(shù)）Table 1 Comparison of evaluation indexes for Pareto solution set of different surrogate-based optimization algorithms （DTLZ5 function）

從圖3 也可以看出，相對(duì)于PEIM 和VFPEIM 方法，VF-PEIM-MVU 方法多了降維分析-多目標(biāo)再次優(yōu)化循環(huán)過(guò)程，可能帶來(lái)額外的計(jì)算花費(fèi)。但實(shí)際上由于降維分析過(guò)程非常快，而且基于降維后的目標(biāo)再次進(jìn)行多目標(biāo)優(yōu)化的計(jì)算花費(fèi)顯著低于高維多目標(biāo)直接優(yōu)化，因而圖3 所示的過(guò)程多出來(lái)的這一步（降維分析-多目標(biāo)再次優(yōu)化）并不會(huì)帶來(lái)計(jì)算花費(fèi)的顯著增加，反而加速了整個(gè)高維多目標(biāo)優(yōu)化問(wèn)題的收斂，以及縮短了搜索到相同Pareto 解集所需要的時(shí)間。

4 高速旋翼翼型高維多目標(biāo)優(yōu)化設(shè)計(jì)

7%厚度旋翼翼型位于槳葉槳尖剖面，前飛時(shí)前行馬赫數(shù)非常高，激波強(qiáng)度和激波位置隨馬赫數(shù)變化非常敏感。據(jù)此，建立了考慮高速馬赫數(shù)不確定性的多目標(biāo)穩(wěn)健優(yōu)化模型，即

式中：μ((Ma))和σ((Ma))分別代表馬赫數(shù)Ma服從均勻分布U( )0.85， 8.87 ，升 力 系 數(shù)Cl=0 時(shí)的均值和方差；ω1和ω2為強(qiáng)調(diào)性能表現(xiàn)和穩(wěn)健性的權(quán)值，ω1越大則越強(qiáng)調(diào)性能表現(xiàn)，ω2越大則越強(qiáng)調(diào)穩(wěn)健性，反之亦然；t和t0分別代表候選翼型最大厚度和最大厚度約束值，本文中取t0=7%；Re和Alpha 分別為流場(chǎng)雷諾數(shù)和翼型迎角。高速旋翼翼型要求盡可能低的俯仰力矩系數(shù)，以降低配平阻力。因此本文在式（25）中將一段范圍內(nèi)（Ma∈[0.85， 8.87]）的力矩系數(shù)絕對(duì)值最大值作為目標(biāo)進(jìn)行設(shè)計(jì)分析，以研究力矩系數(shù)與阻力系數(shù)等其他氣動(dòng)特性的相關(guān)性，為尋找最優(yōu)的多目標(biāo)穩(wěn)健設(shè)計(jì)翼型提供參考[37-38］。

基于圖3所示優(yōu)化流程，首先采用LHS抽樣方法在設(shè)計(jì)空間內(nèi)產(chǎn)生4n個(gè)初始高可信度樣本以及50n個(gè)初始低可信度樣本，n為設(shè)計(jì)變量數(shù)目。計(jì)算分析方法采用驗(yàn)證的基于剪切應(yīng)力輸運(yùn)（Shear Stress Transfer， SST）的數(shù)值轉(zhuǎn)捩模擬程序，計(jì)算雷諾數(shù)為Re=Ma×7.2×106，來(lái)流湍流度為0.5%。參數(shù)化方法為類函數(shù)/形函數(shù)（CST）方法，設(shè)計(jì)變量為d=24 個(gè)。依然采用第2 節(jié)提出的層流翼型網(wǎng)格生成方法分別生成不同網(wǎng)格量的數(shù)值網(wǎng)格，其中高、低可信度計(jì)算使用的網(wǎng)格量分別為126 360 和19 200。2 種網(wǎng)格近場(chǎng)分布如圖5 所示，x和z分別為翼型在流場(chǎng)中沿來(lái)流和垂直于來(lái)流方向的坐標(biāo)，c為弦長(zhǎng)。

圖5 高、低可信度CFD 數(shù)值模擬所使用的密網(wǎng)格和稀網(wǎng)格近場(chǎng)分布Fig.5 Distributions of fine and coarse grids used in highfidelity and low-fidelity CFD simulation

對(duì)這7 個(gè)目標(biāo)使用基于分解的多目標(biāo)進(jìn)化（MOEA/D）算法的VF-PEIM 并行加點(diǎn)準(zhǔn)則進(jìn)行優(yōu)化。MOEA/D 算法使用邊界交叉聚合方法，鄰居個(gè)數(shù)為20，粒子數(shù)為2 100，進(jìn)化代數(shù)為800。在加點(diǎn)子優(yōu)化中，均使用單目標(biāo)加強(qiáng)學(xué)習(xí)粒子群（CLPSO）算法搜索新的加點(diǎn)位置，粒子數(shù)為設(shè)計(jì)變量個(gè)數(shù)，算法停止標(biāo)準(zhǔn)為連續(xù)搜索到的加點(diǎn)距離變化小于10-8。其中每個(gè)子迭代中基于優(yōu)化的Pareto 解均選擇10 個(gè)新的高可信度樣本點(diǎn)添加到樣本庫(kù)中，即Nc=10 后停止子優(yōu)化過(guò)程。加點(diǎn)優(yōu)化迭代50次后，進(jìn)行評(píng)估獲得476個(gè)初步收斂的非支配解。詳細(xì)流程如圖3所示。對(duì)獲得的收斂解通過(guò)3.1 節(jié)陳述的降維方法進(jìn)行目標(biāo)降維分析，獲得核矩陣K∈R7×7的特征值和特征向量。對(duì)特征值按照降序的順序進(jìn)行排列，如表2 所示。可以發(fā)現(xiàn)特征值從e1到e7迅速衰減。本文按照θ=95%能量貢獻(xiàn)準(zhǔn)則進(jìn)行截?cái)噙x擇，保留前3 個(gè)重大特征值及主特征向量，其中對(duì)應(yīng)于前3 個(gè)特征值的主特征向量（列向量）按順序被放在表3，即V1～V3。在表格中，第1個(gè)列向量即最重要的主特征向量V1中第1個(gè)元素和第5 個(gè)元素分別為最負(fù)和最正的元素，其對(duì)應(yīng)的第1個(gè)和第5個(gè)最矛盾的目標(biāo)，且對(duì)主特征向量貢獻(xiàn)最大，因此第1個(gè)和第5個(gè)目標(biāo)進(jìn)入降維目標(biāo)集合F={f1，f5}。第1個(gè)目標(biāo)優(yōu)化高速阻力及阻力發(fā)散特性，要求翼型具有更低的彎度和靠后的最大厚度位置，而第5個(gè)目標(biāo)提高低速升力及升阻比特性，要求翼型具有更大的彎度以及更大的前緣半徑，因此其存在天然的矛盾。在第2個(gè)主特征向量V2中，第4個(gè)元素vp2和第6個(gè)元素vn2為分別為最正和最負(fù)的值，且滿足|vn2|＜0.8vp2，因此選擇vp2對(duì)應(yīng)的目標(biāo)f4進(jìn)入目標(biāo)集合F={f1，f5，f4}。在第3 個(gè)主特征向量V3中，第2個(gè)元素vn3和第6個(gè)元素vp3分別對(duì)應(yīng)最負(fù)和最正的值，且滿足vp3＜0.8|vn3|，因此選擇vn3對(duì)應(yīng)的目標(biāo)f2進(jìn)入目標(biāo)集合F={f1，f5，f4，f2}。至此，通過(guò)降維分析成功選取了4個(gè)目標(biāo)進(jìn)入降維目標(biāo)集合，略去了3個(gè)冗余目標(biāo)。

表2 核矩陣特征值Table 2 Eigenvalues of kernel matrix

在獲得以上降維目標(biāo)集合F 后，對(duì)降維目標(biāo)集合F 中的目標(biāo)對(duì)應(yīng)的目標(biāo)向量H={f1，f5，f4，f2}進(jìn)行相關(guān)性分析，計(jì)算目標(biāo)相關(guān)矩陣R。從表3 中可以看出這4 個(gè)目標(biāo)貢獻(xiàn)方向互不相同，其中Rθ=0.639 916，因此不再剔除任何目標(biāo)，如表4 所示。通過(guò)對(duì)這4 個(gè)目標(biāo)分析可以發(fā)現(xiàn)，第2 個(gè)目標(biāo)為力矩特性要求，第5 個(gè)目標(biāo)為降低中等升力系數(shù)下的阻力即提高其升阻比，在無(wú)激波狀態(tài)下其優(yōu)化潛力很小。因此通過(guò)將第2、第5 個(gè)目標(biāo)，以及冗余目標(biāo)變?yōu)榧s束，使原始7 目標(biāo)優(yōu)化問(wèn)題轉(zhuǎn)化為如下雙目標(biāo)優(yōu)化模型：

表3 對(duì)應(yīng)于特征值順序的特征向量（列向量）Table 3 Eigenvectors （column vector） for corresponding eigenvalues

表4 降維后目標(biāo)相關(guān)性分析Table 4 Correlation analysis of reduced objectives

式中：μ((Ma))和σ((Ma))分別代表馬赫數(shù)Ma服從均勻分布U( )0.85， 8.87 ，升力系數(shù)Cl=0 時(shí)的均值和方差；ω1和ω2為強(qiáng)調(diào)性能表現(xiàn)和穩(wěn)健性的權(quán)值參數(shù)，ω1越大則越強(qiáng)調(diào)性能表現(xiàn)，ω2越大則越強(qiáng)調(diào)穩(wěn)健性，反之亦然；ω3為第2 個(gè)目標(biāo)的權(quán)值系數(shù)；t和t0分別代表候選翼型最大厚度和最大厚度約束值，本文中取t0=7%。Cd0、Cd1、Cd2和Cd3分 別 代 表原 始 翼 型 在 設(shè) 計(jì) 狀態(tài)：①M(fèi)a=0.6，Cl=0.6，Alpha=4o；②Ma=0.5，Cl=0.6，Alpha=5o； ③Ma=0.4，Cl=0.9，Alpha=8o；④Ma=0.5，Cl=0.8，Alpha=6o下的阻力系數(shù)，而Cm0代表原始翼型在Ma∈[0.85， 8.87]范圍內(nèi)，Cl=0 時(shí)最大力矩系數(shù)的絕對(duì)值。通過(guò)將這5 個(gè)非主要沖突目標(biāo)轉(zhuǎn)變?yōu)橛行Ъs束，可有效縮減翼型設(shè)計(jì)空間和引導(dǎo)搜索方向，進(jìn)而提升多目標(biāo)優(yōu)化效率。

對(duì)以上兩目標(biāo)優(yōu)化問(wèn)題，目前許多流行的多目標(biāo)優(yōu)化算法均能取得良好的優(yōu)化效果。首先由LHS 方法產(chǎn)生4d個(gè)初始高可信度樣本以及50d個(gè)初始低可信度樣本，構(gòu)建AMF-PCK 代理模型。其中d為設(shè)計(jì)變量個(gè)數(shù)。本文使用經(jīng)典的MOEA/D 算法進(jìn)行優(yōu)化搜索，種群大小為200 并更新1 000 代，在每次搜索結(jié)束后采用提出的VF-PEIM 并行樣本填充策略選擇Nc=5 個(gè)新的高可信度樣本到樣本庫(kù)中。高、低可信度CFD 數(shù)值分析所使用的計(jì)算網(wǎng)格近場(chǎng)分布如圖5 所示。迭代60 次后收斂，進(jìn)行評(píng)估得到1 組收斂的Pareto 解 集，如 圖6 所 示，目 標(biāo)1（Obj1）和 目 標(biāo)2（Obj2）。圖7 給出了多目標(biāo)加點(diǎn)優(yōu)化過(guò)程中每次搜索完成加點(diǎn)時(shí)預(yù)測(cè)平均阻力系數(shù)絕對(duì)誤差（Absolute Error）隨迭代次數(shù)（Iteration Number）的變化。可以看出，阻力系數(shù)預(yù)測(cè)誤差從最初的接近10×10-4（10 counts）到第40 次迭代時(shí)已經(jīng)小于2 counts，代理模型預(yù)測(cè)精度基本收斂。

圖6 多目標(biāo)優(yōu)化得到的Pareto 前緣Fig.6 Pareto solutions obtained from multi-objective optimization

圖7 優(yōu)化過(guò)程中阻力均值絕對(duì)誤差Fig.7 Absolute errors of mean values of drag coefficients in optimization process

為了比較優(yōu)化效果，從收斂的Pareto 前緣中選擇3個(gè)有代表性的翼型進(jìn)行對(duì)比，如圖8所示，分別 記 為OPT0701、OPT0702 和OPT0703 翼 型。其中OPT0701 翼型是一種高速與低速特性比較均衡的外形。而OPT0702 和OPT0703 翼型更側(cè)重于高速特性，降低高速下的激波阻力并集中提高阻力發(fā)散馬赫數(shù)，以滿足高速旋翼翼型特殊設(shè)計(jì)需求。選擇的3 個(gè)翼型外形和經(jīng)典的OA407 外形對(duì)比如圖8 所示。可以看出，優(yōu)化的3 個(gè)翼型相對(duì)于OA407，彎度均有所降低，其中OPT0702 和OPT0703 彎度下降較多。OPT0701 相對(duì)于OA407 前緣半徑略有減少，最大相對(duì)厚度位置后移。OPT0702翼型前緣半徑以及前緣厚度相對(duì)于OPT0703 翼型變大，而OPT0703 翼型最大相對(duì)厚度位置相對(duì)于OPT0703和OA407均后移。

圖8 優(yōu)化翼型與對(duì)比翼型外形對(duì)比Fig.8 Comparison between optimized airfoils and comparative airfoil

圖9給出了零升狀態(tài)4個(gè)翼型壓力分布（Cp）曲線對(duì)比。OA407 上下表面均出現(xiàn)一道強(qiáng)激波，而OPT0701、OPT0702 以及OPT0703 翼型上表面均由雙弱激波代替，并且其下表面也出現(xiàn)了比OA407 更弱的一道激波。尤其OPT0703 翼型在上下表面的激波強(qiáng)度均要弱于OPT0702 和OPT0701翼型。零升阻力發(fā)散特性曲線如圖10所示，Cm表示俯仰力矩系數(shù)。OPT0701、OPT0702和OPT0703 翼型阻力發(fā)散馬赫數(shù)分別為0.861、0.869 和0.875，相對(duì)于OA407 翼型的0.854 阻力發(fā)散馬赫數(shù)，設(shè)計(jì)翼型阻力發(fā)散馬赫數(shù)分別提高了0.007、0.015 和 0.021。 并 且 OPT0702 和OPT0703 翼型零升阻力系數(shù)也顯著低于OA407翼型。在Ma=0.87 下，3 個(gè)設(shè)計(jì)翼型相對(duì)于OA407翼型阻力系數(shù)均降低超過(guò)20 counts。力矩特性曲線顯示設(shè)計(jì)翼型OPT0702 和OPT0703 力矩系數(shù)絕對(duì)值在馬赫數(shù)0.80～0.87 范圍內(nèi)均小于0.02，并明顯好于OA407 翼型的力矩特性，而OA407翼型力矩特性隨阻力發(fā)散后急劇變差。

圖9 零升狀態(tài)壓力分布對(duì)比Fig.9 Comparison of pressure distributions at Cl=0

圖10 氣動(dòng)力特性曲線對(duì)比Fig.10 Comparison of aerodynamic characteristics curves

7%厚度翼型位于槳葉尖部，（共軸雙旋翼直升機(jī)）高速前轉(zhuǎn)時(shí)馬赫數(shù)超過(guò)0.87，上下表面均非常容易出現(xiàn)強(qiáng)激波，因此提高其阻力發(fā)散馬赫數(shù)以及降低基礎(chǔ)阻力對(duì)減少飛行阻力具有重要價(jià)值。同時(shí)在機(jī)動(dòng)以及懸停時(shí)，巡航馬赫數(shù)降低，為了配合槳葉其他部分產(chǎn)生收益，希望其能提供較大的最大升力系數(shù)以及一定的升阻比。當(dāng)然，由于其位于槳尖部分較小的一段，因此其所能提供的載荷也非常有限。圖11～圖14 分別給出了設(shè)計(jì)翼型與OA407 翼型在Ma=0.3、Ma=0.4、Ma=0.5、Ma=0.6 下升力系數(shù)（Cl）隨迎角（Alpha）變化曲線，以及升阻比（K）隨升力系數(shù)（Cl）變化曲線對(duì)比。可以看出，設(shè)計(jì)的OPT0701翼型在Ma=0.3～0.6 范圍內(nèi)相比于OA407 翼型最大升力系數(shù)和最大升阻比均有不同程度提高。2 個(gè)高速翼型OPT0702 和OPT0703 在Ma=0.3～0.4 下升力特性和升阻比均相對(duì)于OA407有所降低，其中OPT0703 翼型相對(duì)于OPT0702翼型低速最大升力系數(shù)要低。到Ma=0.5～0.6時(shí)，OPT0702 翼型升力和升阻比特性基本與OA407 相 當(dāng)，OPT0703 翼 型 在Ma=0.6 下 也 提供了和OA407 翼型接近的升力和升阻比特性。因?yàn)橐獏f(xié)調(diào)旋翼翼型的多個(gè)氣動(dòng)特性，提高翼型高速阻力發(fā)散特性，翼型彎度降低，最大厚度位置后移，其低速特性（包括最大升力系數(shù)和最大升阻比）將會(huì)有所損失。并且為了維持穩(wěn)健的高速阻力發(fā)散特性，使得OPT0702 和OPT0703 翼型上表面較為平坦，大迎角時(shí)非常容易出現(xiàn)強(qiáng)激波或者分離流，導(dǎo)致阻力增加。綜合來(lái)看，OPT0701 翼型較OA407 翼型氣動(dòng)特性得到全面提升，包括高、低速特性，但高速特性提升有限；而OPT0703 翼型更強(qiáng)調(diào)高速特性，低速特性損失較大；OPT0702 翼型在保持一定高速特性的同時(shí)，低速特性損失較小。可見(jiàn)設(shè)計(jì)翼型要維持良好的低速特性（不差于OA407 翼型），高速特性將難以顯著提高，進(jìn)而難以滿足高速直升機(jī)特殊的設(shè)計(jì)需求。

圖11 Ma=0.3 下氣動(dòng)特性曲線對(duì)比Fig.11 Comparison of aerodynamic characteristics curves at Ma=0.3

圖12 Ma=0.4 下氣動(dòng)特性曲線對(duì)比Fig.12 Comparison of aerodynamic characteristics curves at Ma=0.4

圖13 Ma=0.5 下氣動(dòng)特性曲線對(duì)比Fig.13 Comparison of aerodynamic characteristics curves at Ma=0.5

圖14 Ma=0.6 下氣動(dòng)特性曲線對(duì)比Fig.14 Comparison of aerodynamic characteristics curves at Ma=0.6

高速直升機(jī)槳尖翼型幾乎不承擔(dān)載荷，設(shè)計(jì)翼型OPT0702 和OPT0703 的主要目的在于提高阻力發(fā)散馬赫數(shù)和減小激波阻力，以滿足槳尖翼型設(shè)計(jì)需求。結(jié)果也表明，此類翼型高速阻力發(fā)散特性與Ma=0.3～0.4 下升力特性矛盾性要強(qiáng)于Ma=0.5～0.6 下升力特性，這與上述目標(biāo)降維分析結(jié)果一致。一些典型的氣動(dòng)特性對(duì)比如表5 所示。其中：Madd0表示零升阻力發(fā)散馬赫數(shù)；Cd0表示零升阻力系數(shù)；Cm0表示零升力矩系數(shù)；Clmax表示最大升力系數(shù)；Kmax表示最大升阻比。可以看出，OPT0701 翼型所有指標(biāo)均較OA407 翼型有不同程度提高，而OPT0702 和OPT0703 翼 型 均 相 對(duì) 于OA407 翼型高速特性指標(biāo)提升明顯，而低速特性指標(biāo)有所損失。

表5 7%厚度翼型氣動(dòng)特性對(duì)比Table 5 Aerodynamic characteristics of 7% airfoils

5 結(jié) 論

隨著剛性共軸雙旋翼直升機(jī)飛行速度的提升，對(duì)高速旋翼翼型提出了更高的要求，包括更高的零升阻力發(fā)散馬赫數(shù)和更低的基礎(chǔ)零升阻力，中等速度高升阻比，低速高升力特性以及所有狀態(tài)下更好的力矩特性，由此引發(fā)了高維多目標(biāo)和復(fù)雜約束優(yōu)化的難題。

本文針對(duì)這一難題，建立了多目標(biāo)穩(wěn)健優(yōu)化方法，并發(fā)展了有效的多目標(biāo)降維策略，以獲得滿足指標(biāo)要求的最優(yōu)解。同時(shí)建立了基于多可信度代理模型AMF-PCK 的高效多目標(biāo)穩(wěn)健優(yōu)化方法，以及多目標(biāo)優(yōu)化加點(diǎn)方法使設(shè)計(jì)效率進(jìn)一步提高。利用所建立的高效多目標(biāo)穩(wěn)健優(yōu)化方法對(duì)7%厚度高速旋翼翼型進(jìn)行了優(yōu)化設(shè)計(jì)，設(shè)計(jì)翼型的氣動(dòng)表現(xiàn)與經(jīng)典的OA407 翼型在高、中、低馬赫數(shù)下進(jìn)行了全面比較。

結(jié)果表明，相對(duì)于經(jīng)典的OA407 翼型，設(shè)計(jì)的OPT0701 翼型較OA407 翼型高、中、低速氣動(dòng)特性指標(biāo)全面提升，但高速特性提升非常有限；而OPT0702 和OPT0703 翼型在低速特性損失不大的情況下，高速特性提升顯著且更加穩(wěn)健。3 個(gè)設(shè)計(jì)翼型所有狀態(tài)下的力矩特性也均好于OA407 翼型。設(shè)計(jì)結(jié)果再次說(shuō)明了旋翼翼型各個(gè)目標(biāo)特性強(qiáng)烈的矛盾性。考慮到剛性共軸雙旋翼直升機(jī)翼型特殊的設(shè)計(jì)要求，設(shè)計(jì)翼型相對(duì)于經(jīng)典的OA 系列翼型展示了顯著的優(yōu)勢(shì)。設(shè)計(jì)結(jié)果驗(yàn)證了提出的高維多目標(biāo)優(yōu)化方法的有效性和可靠性。