基于梯度增強Kriging 方法的水下航行器結(jié)構(gòu)優(yōu)化設(shè)計

陳力銘，邱浩波，高亮

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

0 引 言

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

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

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

1 梯度增強Kriging 模型

1.1 標準型梯度增強Kriging 模型

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

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

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

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

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

1.2 縮減型梯度增強Kriging 模型

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

2 基于梯度增強Kriging 的優(yōu)化方法

2.1 近似駐點概率

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

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

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

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

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

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

2.2 加點準則

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

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

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

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

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

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

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

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

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

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

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

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

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

3 算 例

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

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

圖4 設(shè)計變量示意圖Fig. 4 Diagram of the design variables

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

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

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

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

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

圖6 水下航行器結(jié)構(gòu)優(yōu)化前后對比Fig. 6 Comparison of the underwater vehicle structural design before and after optimization

4 結(jié) 語

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

1） 探索合理的梯度/偏導(dǎo)數(shù)篩選方法，以緩解船舶結(jié)構(gòu)優(yōu)化中因設(shè)計變量的增多而導(dǎo)致計算量劇增的問題；同時，致力于開發(fā)船舶結(jié)構(gòu)有限元分析所對應(yīng)的高效梯度計算（伴隨方法）程序。

2） 探索基于梯度增強Kriging 方法的船舶結(jié)構(gòu)優(yōu)化過程中約束（例如強度、壽命等）昂貴、優(yōu)化目標（例如重量、模態(tài)等）多的解決方法。

3） 探索基于梯度增強Kriging 方法的船舶結(jié)構(gòu)優(yōu)化過程中涉及離散（例如型材數(shù)量、規(guī)格等）/混合變量時的解決方法。