基于新型多可信度代理模型的多目標優化方法

趙歡，高正紅，夏露

西北工業大學 航空學院，西安 710072

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

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

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

1 剛性共軸雙旋翼直升機槳葉翼型高維多目標優化難題

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

1）高速前行時，槳尖翼型工作狀態馬赫數超過0.87，波阻急劇增加，前飛需用功率增加，設計翼型應具有高的零升阻力發散馬赫數，盡可能低的零升阻力系數，以及維持低的零升俯仰力矩系數。槳尖載荷要求可適當降低。

2） 在 跨 聲 速 狀 態 范 圍 內(0.8 ≤Ma ≤Madd0，Madd0為零升阻力發散馬赫數)，激波位置和強度隨馬赫數等狀態變化非常敏感，設計翼型在此范圍內應該具有穩健的氣動表現。

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

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

5）在所有馬赫數范圍或整個飛行包線范圍內，各剖面設計翼型應具有極低的俯仰力矩系數以減小槳葉扭轉及控制系統的操縱載荷。

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

以上這些指標要求達10 余項之多，并通常是相互矛盾和制約的，面臨高維多目標和復雜約束優化難題。因此，發展高效的高維多目標設計模型和設計方法，并設計滿足嚴苛設計要求的剛性共軸雙旋翼直升機翼型非常重要。

高維多目標優化問題可以數學定義為

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

式（2）為僅考慮設計點氣動特性的確定性優化設計模型，無法兼顧非設計點的氣動特性，如來流馬赫數變化時的氣動特性。

如果考慮飛行狀態等各種不確定源對翼型氣動特性的影響，將式（2）模型直接修改為

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

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

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

相對于傳統的確定性多目標優化，穩健多目標優化引入了各個目標表現和穩健性，以及各個約束的可行穩健性測量。因而目標個數成倍增加，優化可行域變窄且更加復雜，函數評估數量增加。目前一些經典的多目標算法如NSGA-Ⅲ等對于高維多目標（目標個數超過3 個）問題，其搜索和收斂能力均下降，尤其一些冗余的目標會影響算法的收斂效率和真實表現。因此，通過有效的目標降維方法結合目前經典的多目標算法，將能有效解決高維多目標難題。同時，基于前期發展的多可信度代理建模技術（Adaptive Multi-Fidelity Polynomial Chaos-Kriging， AMF-PCK），本文提出了高效的高維多目標穩健優化平臺，以滿足工程使用。

2 自適應多可信度多項式混沌-Kriging 代理模型

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

AMF-PCK 在任意點X 的預測值為

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

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

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

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

其中：α0以及{αci}是需要被精確訓練的縮放因子。子集Ac應該通過一個合適的標準進行選擇。

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

圖1 自適應選擇最相關的低可信度多項式集合和最優的矯正多項式集合Fig.1 Adaptive selection of the most relevant lowfidelity polynomial set and optimal correction polynomial set

在式（11）中的各項被確定后，通過標準Universal Kriging 擬合過程即可得到AMF-PCK 在未知點的預測模型。同Universal Kriging 模型一樣，AMF-PCK 滿足最小化均方誤差和無偏估計假設（Best Linear Unbiased Predictor，BLUP）。AMF-PCK 預測模型在未知點X的預測值為

式中：Mc(Mc≤M)表示校正展開項中多項式項數量。

縮放因子以及校正展開多項式系數向量α=[α0，αc1，αc2，αcMc]T可通過廣義最小 二乘方 法估計得到，即

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

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

3 基于新型AMF-PCK 模型的多目標優化方法

3.1 非線性目標降維策略

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

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

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

上述凸半正定優化問題可通過SeDuMi[29］求解。如果使用PCA 線性降維，則可直接計算協方差矩陣為

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

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

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

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

8）使用MOEAs 算法對新目標集合F*選擇的目標進行搜索，得到Pareto 最優解集合，折衷選擇最合適的設計，進行評估分析。

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

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

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

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

通過對式（19）單目標優化搜索可獲得1 個最有潛力改進點位置，即

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

第Nc個樣本點位置為

式（21）～式（24）的優化均是單目標優化問題，因而相對多目標優化效率很高。VF-PEIM 通過定義每個修訂周期合適的加點數，即可有效地獲得全局最優解，從而兼顧了多目標優化效率和搜索能力。

3.3 高維多目標優化框架

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

圖3 基于AMF-PCK 的多目標優化流程Fig.3 Many-objective optimization flow chart based on AMF-PCK

3.4 數值算例測試

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

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

圖4 使用基于PEIM 和VF-PEIM 加點優化方法分別執行10 次迭代IGD 和得到的Pareto 解集對比（DTLZ5 函數）Fig.4 Comparison of IGD and Pareto solution set using infilling-based PEIM and VF-PEIM methods respectively based on 10 iterations （DTLZ5 function）

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

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

表1 不同代理優化算法的Pareto 解集的評價指標對比（DTLZ5 函數）Table 1 Comparison of evaluation indexes for Pareto solution set of different surrogate-based optimization algorithms （DTLZ5 function）

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

4 高速旋翼翼型高維多目標優化設計

7%厚度旋翼翼型位于槳葉槳尖剖面，前飛時前行馬赫數非常高，激波強度和激波位置隨馬赫數變化非常敏感。據此，建立了考慮高速馬赫數不確定性的多目標穩健優化模型，即

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

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

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

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

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

在獲得以上降維目標集合F 后，對降維目標集合F 中的目標對應的目標向量H={f1，f5，f4，f2}進行相關性分析，計算目標相關矩陣R。從表3 中可以看出這4 個目標貢獻方向互不相同，其中Rθ=0.639 916，因此不再剔除任何目標，如表4 所示。通過對這4 個目標分析可以發現，第2 個目標為力矩特性要求，第5 個目標為降低中等升力系數下的阻力即提高其升阻比，在無激波狀態下其優化潛力很小。因此通過將第2、第5 個目標，以及冗余目標變為約束，使原始7 目標優化問題轉化為如下雙目標優化模型：

表3 對應于特征值順序的特征向量（列向量）Table 3 Eigenvectors （column vector） for corresponding eigenvalues

表4 降維后目標相關性分析Table 4 Correlation analysis of reduced objectives

式中：μ((Ma))和σ((Ma))分別代表馬赫數Ma服從均勻分布U( )0.85， 8.87 ，升力系數Cl=0 時的均值和方差；ω1和ω2為強調性能表現和穩健性的權值參數，ω1越大則越強調性能表現，ω2越大則越強調穩健性，反之亦然；ω3為第2 個目標的權值系數；t和t0分別代表候選翼型最大厚度和最大厚度約束值，本文中取t0=7%。Cd0、Cd1、Cd2和Cd3分 別 代 表原 始 翼 型 在 設 計 狀態：①Ma=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下的阻力系數，而Cm0代表原始翼型在Ma∈[0.85， 8.87]范圍內，Cl=0 時最大力矩系數的絕對值。通過將這5 個非主要沖突目標轉變為有效約束，可有效縮減翼型設計空間和引導搜索方向，進而提升多目標優化效率。

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

圖6 多目標優化得到的Pareto 前緣Fig.6 Pareto solutions obtained from multi-objective optimization

圖7 優化過程中阻力均值絕對誤差Fig.7 Absolute errors of mean values of drag coefficients in optimization process

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

圖8 優化翼型與對比翼型外形對比Fig.8 Comparison between optimized airfoils and comparative airfoil

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

圖9 零升狀態壓力分布對比Fig.9 Comparison of pressure distributions at Cl=0

圖10 氣動力特性曲線對比Fig.10 Comparison of aerodynamic characteristics curves

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

圖11 Ma=0.3 下氣動特性曲線對比Fig.11 Comparison of aerodynamic characteristics curves at Ma=0.3

圖12 Ma=0.4 下氣動特性曲線對比Fig.12 Comparison of aerodynamic characteristics curves at Ma=0.4

圖13 Ma=0.5 下氣動特性曲線對比Fig.13 Comparison of aerodynamic characteristics curves at Ma=0.5

圖14 Ma=0.6 下氣動特性曲線對比Fig.14 Comparison of aerodynamic characteristics curves at Ma=0.6

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

表5 7%厚度翼型氣動特性對比Table 5 Aerodynamic characteristics of 7% airfoils

5 結 論

隨著剛性共軸雙旋翼直升機飛行速度的提升，對高速旋翼翼型提出了更高的要求，包括更高的零升阻力發散馬赫數和更低的基礎零升阻力，中等速度高升阻比，低速高升力特性以及所有狀態下更好的力矩特性，由此引發了高維多目標和復雜約束優化的難題。

本文針對這一難題，建立了多目標穩健優化方法，并發展了有效的多目標降維策略，以獲得滿足指標要求的最優解。同時建立了基于多可信度代理模型AMF-PCK 的高效多目標穩健優化方法，以及多目標優化加點方法使設計效率進一步提高。利用所建立的高效多目標穩健優化方法對7%厚度高速旋翼翼型進行了優化設計，設計翼型的氣動表現與經典的OA407 翼型在高、中、低馬赫數下進行了全面比較。

結果表明，相對于經典的OA407 翼型，設計的OPT0701 翼型較OA407 翼型高、中、低速氣動特性指標全面提升，但高速特性提升非常有限；而OPT0702 和OPT0703 翼型在低速特性損失不大的情況下，高速特性提升顯著且更加穩健。3 個設計翼型所有狀態下的力矩特性也均好于OA407 翼型。設計結果再次說明了旋翼翼型各個目標特性強烈的矛盾性。考慮到剛性共軸雙旋翼直升機翼型特殊的設計要求，設計翼型相對于經典的OA 系列翼型展示了顯著的優勢。設計結果驗證了提出的高維多目標優化方法的有效性和可靠性。