基于梯度增強Kriging 方法的水下航行器結構優化設計

陳力銘，邱浩波，高亮

華中科技大學 機械科學與工程學院，湖北 武漢 430074

0 引 言

隨著船舶結構復雜程度的不斷提升，高精度的有限元仿真軟件已成為船舶結構分析的主流工具，而直接調用耗時的仿真分析程序進行優化已難以滿足船舶設計快速高效的需求，因此，發展基于代理模型的高效優化方法是船舶結構優化設計的重點研究領域之一[1]。所謂代理模型，是指對一系列通過仿真分析獲得的輸入/輸出樣本進行擬合得到的統計近似模型。在代理模型提供的信息引導下的優化可以極大地減少耗時的仿真次數，達到縮短設計周期、節省設計成本的目的。目前，經典的代理模型建模與優化方法已被廣泛用于船舶領域，例如對船舶板架強度和穩定性[2]、潛艇振動聲輻射[3]、水下結構物基座阻抗特性[4]等響應的快速預報，以及對深潛器的多球交接耐壓殼[5]、多用途船貨艙段[6]的結構優化和對海上發射船[7]、集裝箱船[8]等的船型優化。

考慮將梯度信息納入到用于構建代理模型的樣本中，通常可以顯著提高模型的精度[9]。而在眾多的梯度增強代理模型中，梯度增強Kriging (Gradientenhanced Kriging, GEK)模型[10]由于其對非線性問題具有良好的插值能力和獨特的概率模型屬性，所以在一些領域的應用效果明顯，例如飛機機翼/翼型的氣動優化設計[11]、懸臂梁與簡支梁的結構優化設計[12]、軸孔過盈配合的接觸應力分析[13]、核反應堆最大燃料溫度的不確定性量化[14]等。然而，在船舶優化領域，對實現高效梯度計算方法（例如伴隨法[15]）的工具開發尚不成熟，使用梯度增強Kriging模型分析的案例并不多見，因而，利用該模型的優勢進行設計優化的方法仍有待發展。

本文擬將基于梯度增強Kriging 代理模型的優化方法[16]擴展至梯度信息只能通過有限差分法獲取的場景，實現對船舶結構的高效優化設計。即提出基于縮減型梯度增強Kriging 模型的加點策略，僅對有實際改進的采樣位置進行梯度計算，以此減少仿真調用的次數。在加點過程中，首先采用多起點局部優化算法搜索改進期望函數的若干局部最優解作為候選加點位置，然后計算相應的近似駐點概率，并根據改進期望值和近似駐點概率值的一致程度來確定加點位置，達到提高優化效率的目的。最后，以某水下航行器的結構優化設計為對象，對所提方法進行驗證。

1 梯度增強Kriging 模型

1.1 標準型梯度增強Kriging 模型

Kriging 模型的基本假設為：將確定性響應y(x)視為高斯隨機過程Y(x)。不失一般性，描述為

式中，f(x) 為 包含若干回歸項的向量； β為待估計的回歸系數向量；平穩隨機過程Z(·)的均值為0，協方差Cov[Z(x),Z(x′)]=σ2Corr[Z(x),Z(x′)]，其中σ2為 過程方差，Corr[Z(x),Z(x′)]為相關函數。對于梯度增強Kriging，給定采樣矩陣X，其增廣后的樣本數據向量y～中既包含了響應也包含了梯度，記為

作為 (n+nk)×(n+nk) 對 稱矩陣的R～，可依照式中所示分塊建立。對于n×n的 子矩陣R，n×nk的子矩陣?R/?x以 及nk×nk的子矩陣?2R/?x2，其元素分別為響應之間、響應與偏導數之間和偏導數之間的相關函數。當采用高斯型相關函數時，具體表達式可詳見文獻[16]。

令式(3) 關于 β 和 σ2的偏導分別為0，可得到極大似然估計量分別如下：

由于超參數 θ沒有解析解，所以將式(6)和式（7）代回到式(3)，通過優化算法搜索似然函數最大值所對應的超參數作為估計值。為了簡化描述，直接給出梯度增強Kriging 在任意采樣點x處的預測式和均方預測誤差，并分別表示如下：

1.2 縮減型梯度增強Kriging 模型

考慮僅將一部分梯度數據包含到樣本數據向量中，可構建縮減型梯度增強Kriging 模型。針對提高構建梯度增強Kriging 模型效率的問題，Han等[17]建立了一系列考慮一組梯度信息的子模型，并加權求和；Chen 等[18]篩選出了與響應顯著相關的變量所對應的偏導數，將其包含到樣本中。本文僅考慮在部分樣本點處計算梯度，將其包含到樣本數據向量中，建立的縮減型梯度增強Kriging模型可視為文獻[17]所述子模型的普遍形式。

2 基于梯度增強Kriging 的優化方法

2.1 近似駐點概率

工程問題對應的真實函數通常是光滑的[19]，即函數在設計空間內處處可導。對于任意非單調的函數，若其全局最優不在設計空間的邊界上，則必然位于某駐點處。由于梯度增強Kriging 方法可以對梯度進行最優線性無偏預測，所以可用來估計任意點為真實函數的駐點概率。根據式（1）可知，對于梯度，有

對于求解過程，可參見文獻[16]，本文不再贅述。

由此，梯度增強Kriging 方法對偏導數的預測式和預測均方誤差可分別表示為：的情況類似（虛線對稱軸與y軸重合），故省略。

圖1 概率 Pl的幾何表示Fig. 1 Illustrations for the probability Pl

當選取接近0 的正實數 εl(l=1,···,k)作為區間邊界時，近似駐點概率(approximate probability of stationary point, APSP)可定義為

以保證近似駐點概率在整個設計空間內有連續定義。同樣地，對于縮減型梯度增強Kriging 模型，近似駐點概率的計算過程類似，僅需替換相應的預測式和均方預測誤差。

2.2 加點準則

對于梯度增強Kriging 代理模型的優化方法，經 典 的EGO （efficient global optimization）算 法[20]仍然適用。然而，此時利用梯度信息的情況仍限于建模階段，若樣本容量較少，梯度增強Kriging方法對響應的預測結果可能會出現一定程度的偏差，影響改進期望（expected improvement, EI）函數的判斷。這里，以使用基于梯度增強Kriging 模型的EGO 算法對如下一維函數進行全局優化為例進行說明。

通過拉丁超立方采樣[21]得到4 個初始樣本，計算各樣本點處的梯度，建立梯度增強Kriging 代理模型，其預測值和均方預測誤差如圖2 所示。

圖2 梯度增強Kriging 對一維函數值的預測Fig. 2 Prediction of the one-dimensional function by Gradient-enhanced Kriging

根據EGO 算法的原理，應選擇改進期望函數的最大值所在位置作為下一輪采樣點。然而，從圖3(a)中可見，當前改進期望函數的最大值位于菱形標識處，而真實函數的全局最優解位于遠離于此的星形標識處。因此，當前的加點效率因改進期望函數的偏離而受限。若同時考慮近似駐點概率值的大小，如圖3(b) 所示，則真實函數的全局最優解所對應的近似駐點概率值將高于改進期望最大值的所在位置，即有機會選擇到全局最優解附近加點，從而使優化效率得到提高。

圖3 改進期望和近似駐點概率在當前加點階段選擇更新點時的決策對比Fig. 3 Comparison of the expected improvement and APSP in selecting the infill point at the current infill stage

通常，改進期望函數具有多峰特性。Zhan 等[22]提出搜索改進期望函數的多個局部最優位置，以并行的方式加點。當不具備并行條件時，需要合理地選擇加點位置。本文提出了基于縮減型梯度增強Kriging 的優化方法，具體步驟如下：

1)使用拉丁超立方采樣確定采樣矩陣X，計算相應位置的響應與有限差分梯度，得到初始樣本S={X,y～}，建立縮減型梯度增強Kriging 模型。

2)使用拉丁超立方采樣產生Q個位置作為內點法[23]的初始點，對改進期望函數進行搜索，得到包含Q個局部最優解（允許存在重復）位置的集合H。

3)對集合H中元素按對應的改進期望值大小降序排列，得到排序后的集合Hsorted。

4)計算集合Hsorted中元素所對應的近似駐點概率。

5)計算集合Hsorted的 后(Q-1)個元素所對應的一致程度差值 Δq(q=2,···,Q)，計算式為

6)若所有差值 Δq均 不小于0，則在集合Hsorted的第1 個元素處計算響應；否則，在集合Hsorted中對應的最小差值 Δq的元素處計算響應。若新計算的響應小于當前樣本中的最優解，則通過有限差分計算當前位置的梯度；否則，不計算任何位置的梯度。

7)判斷是否滿足終止條件，若滿足，停止加點，輸出優化結果；否則，將新采樣的樣本響應和梯度補充到樣本集合中，更新梯度增強Kriging 模型，重復步驟2～步驟6。

3 算 例

水下航行器由于受到自身發動機工作的影響，外殼會發生振動并激勵外場的海水介質形成輻射聲場，所以分析水下航行器自身振動特性是研究其輻射聲場強度分布的基礎。航行器在水下的振動濕模態分析比其在空氣中的振動干模態分析更復雜：一方面，水下航行器航行時在周圍海水介質的作用下會產生流固耦合效應，無法對其濕模態進行解析求解； 另一方面，水下試驗難度大、成本高，對于全尺寸水下航行器進行振動試驗也很難完成。因此，利用有限元分析軟件對水下航行器進行濕模態分析，可以高精度地模擬其在水下振動時的流固耦合特性，從而為后續的結構優化設計提供分析手段。

本文以某水下航行器的簡化結構優化設計為例，驗證所提優化方法的可行性。為了模擬航行器在水中的懸浮狀態，不對其施加任何位移約束和力載荷，其在水下自由振動時的前6 階模態均為剛體模態（計算數值均接近于0），而將提高其第7 階固有頻率作為優化目標，各設計變量如圖4 所示。待優化的水下航行器基準形的平行舯體直徑φ =2.54 m，平行舯體長L=10.16 m，艏徑R1=4.52 m，艉徑R2=9.42 m，相應的取值范圍如表1所示。此外，水下航行器的材料參數包括：彈性模量2×1011Pa、泊松比0.3、密度7 800 kg/m3、海水密度1 025 kg/m3、水中聲速1 500 m/s。

圖4 設計變量示意圖Fig. 4 Diagram of the design variables

表1 設計變量的取值范圍Table 1 Ranges of the design variables

本文使用了ANSYS 的Modal 模塊來實現水下航行器的振動濕模態分析。首先，導入SOLIDWORKS 中創建的參數化幾何模型，然后，在水下航行器四周建立長方體形的流體域作為外界流場，其尺寸為水下航行器各端緣向外擴展3 m， 并雙選水下航行器和流體域建立一個整體件，以便后續劃分網格時二者公共面處的節點位置重合。水下航行器的有限元網格采取自動劃分方式，得到的網格尺寸為0.25 m，網格數37 388。流體域的有限元網格劃分方式為空間六面體占優，網格尺寸為0.5 m，網格數30 035。本文算例采用了單向流固耦合分析方法。分析時，流體域與水下航行器的公共面為流固耦合面，定義流體域外表面的邊界壓力為0，插入APDL 命令流，調用非對稱模態算法求解結構動力學有限元方程。由于每次仿真分析耗時約6 min，直接調用仿真程序進行優化的效率低下，所以本文算例符合基于縮減型梯度增強Kriging 模型的優化方法應用場景。

以下為優化設計流程：首先，設定初始采樣成本為20 次仿真調用。采用拉丁超立方采樣確定4×4的 初始采樣矩陣X，調用仿真得到4 組響應值，并通過有限差分法計算采樣位置的梯度。需注意，此時初始采樣樣本數據向量中共包含20 組響應和4 組梯度。 然后，運用基于縮減型梯度增強Kriging 模型的優化方法進行序列優化，終止條件為加點過程中調用的仿真次數超過40 次。對于本文算例，取Q=10， εl(l=1,···,4)等于樣本所對應維度偏導絕對值中最大值的千分之一。最后，作為參照，在相同樣本量的前提下，運用基于標準型梯度增強Kriging 模型的EGO 算法和經典的EGO算法進行序列優化，各方法優化過程的迭代曲線如圖5 所示。

圖5 水下航行器第7 階固有頻率優化過程的迭代曲線Fig. 5 Iterative curves of the optimization process of seventh-order natural frequency for the underwater vehicle

由圖5 可知，與基準形參考值相比，3 種優化方法在迭代終止時得到的設計方案一定程度上均提高了水下航行器的第7 階固有頻率，而基于縮減型梯度增強Kriging 模型的優化方法在效率方面最高。圖6 所示為水下航行器優化前后的仿真分析結果。優化后的水下航行器第7 階固有頻率提高了14.6%，從而驗證了所提出的優化方法在船舶優化領域的可行性。

圖6 水下航行器結構優化前后對比Fig. 6 Comparison of the underwater vehicle structural design before and after optimization

4 結 語

本文針對船舶優化領域中梯度信息只能通過有限差分法獲取的問題，提出了基于縮減型梯度增強Kriging 代理模型的優化方法，擴展了標準型梯度增強Kriging 方法的應用范圍，并以某水下航行器結構優化為例對所提方法進行了驗證。結果表明，該方法在船舶結構優化領域具有良好的應用前景。在基于梯度增強Kriging 方法的船舶結構優化領域，未來值得進一步研究的相關問題包括：

1） 探索合理的梯度/偏導數篩選方法，以緩解船舶結構優化中因設計變量的增多而導致計算量劇增的問題；同時，致力于開發船舶結構有限元分析所對應的高效梯度計算（伴隨方法）程序。

2） 探索基于梯度增強Kriging 方法的船舶結構優化過程中約束（例如強度、壽命等）昂貴、優化目標（例如重量、模態等）多的解決方法。

3） 探索基于梯度增強Kriging 方法的船舶結構優化過程中涉及離散（例如型材數量、規格等）/混合變量時的解決方法。