基于改進NSGA-Ⅱ的載人潛水器多開孔耐壓結構優化

田震，劉峰，王萌，陳書培，李捷龍，趙彥凱

(1.哈爾濱工程大學 船舶工程學院， 黑龍江 哈爾濱 150001； 2.哈爾濱工程大學 煙臺研究院， 山東 煙臺 264000)

0 引言

耐壓結構是載人潛水器的重要系統之一，是其實現水下安全航行和完成相關任務的重要保證。耐壓結構所承受的海水壓力隨著其工作水深的增加而增大，在這種情況下，往往通過增加尺度以增強其結構性能，但這必將導致其質量的增加，從而對載人潛水器的總體性能、經濟性等方面造成不利影響。因此，耐壓結構的設計是一個典型的多目標優化問題，需要對質量和結構性能之間所存在的矛盾進行協調和平衡。在耐壓結構的設計中需要對眾多方案進行優選，在這種情況下，優化模型的精度和求解效率，以及優化算法的搜索能力和收斂性等都將直接影響到優化解的全局性和求解效率。

有限元方法已在耐壓結構的設計中得到了廣泛的應用，但對于需要反復迭代的優化過程來說，其仍存在計算量巨大的缺點[1]。近似模型對樣本點的輸入(設計變量)和輸出(響應)進行擬合，進一步將二者的表達方式作為替代模型用于分析，不僅保證了問題分析的精度和準確性，對研究效率的提升也有重要的影響。目前，響應面模型[2]、Kriging模型[3]、RBF[4]等近似模型已在耐壓結構的研究中得到了應用。而建立近似模型所需的樣本點依然采用有限元分析完成，參數化分析方法實現了耐壓結構自動建模與分析，可提高設計效率。由于所采用的軟件和平臺的不同，所采用的耐壓結構參數化分析方法也不盡相同，楊卓懿等[5]利用Apdl語言實現了環肋圓柱殼的參數化分析；宋保維等[6]在解決了基于VC語言的iSight軟件集成開發的問題的基礎上，實現了耐壓結構的自動設計，并將其應用于環肋耐壓結構的優化設計，提高了耐壓結構性能和設計效率；劉峰等[7]所研究的耐壓結構參數化設計流程，重點針對Abaqus軟件的二次開發開展了研究，進一步解決了Abaqus與多學科設計優化軟件iSight的集成與通信問題。網格尺寸的選擇對有限分析有重要的影響，網格劃分過大將影響計算精度，而網格劃分較小則會導致計算成本的提高，上述研究中對于網格劃分缺乏深入研究。此外，載人潛水器耐壓結構布置有多個較大的開孔，必將導致其設計變量的增加，從而加大了參數化分析和研究的難度。

在多目標優化算法中，改進非支配排序遺傳算法(NSGA-Ⅱ)已在耐壓結構優化中得到了應用，且取得了明顯的優化效果。朱學康等[8]圍繞具有離散特性的圓柱殼多目標優化進行了研究，并利用NSGA-Ⅱ進行求解；汪志強等[9]在考慮了綜合裕度的基礎上，建立了環肋圓柱殼的多目標優化模型，采用NSGA-Ⅱ算法進行了優化求解；何衍儒等[10]將浮重比和內部有效體積作為目標函數，利用NSGA-Ⅱ進行了耐壓結構多目標優化的求解。但NSGA-Ⅱ的全局搜索能力和快速收斂能力還存在一定的不足，導致這些不足的原因是模擬二進制交叉(SBX)算子的采用。正態分布交叉(NDX)算子的空間開發和探索能力與SBX算子相比有較大提升，將NDX算子應用于NSGA-Ⅱ算法中，可以有效提升算法解集的多樣性和收斂性[11-12]。

本文以載人潛水器多開孔耐壓結構為研究對象，在完成有限元模型的建立、計算網格的確定、參數化分析的實現以及設計變量靈敏度分析的基礎上，選擇擬合精度較高近似模型作為替代模型，進一步建立多開孔耐壓柱殼多目標優化模型。以NDX算子作為交叉算子對NSGA-Ⅱ進行改進，利用ZDT系列測試函數對改進算法INSGA-Ⅱ的改進效果進行測試。進行INSGA-Ⅱ與傳統NSGA-Ⅱ算法的多目標優化的求解對比，以對INSGA-Ⅱ算法的有效性進行驗證。

1 有限元分析方法及網格無關性分析

1.1 耐壓結構有限元模型

根據設計需求和設計規范[13]確定耐壓結構形式和尺寸，建立基于Abaqus軟件的有限元分析模型，如圖1所示。

圖1 耐壓結構有限元模型

對于計算載荷Pj，可換算成[14-15]

Pj=0.009 8hj

(1)

式中：hj為計算深度。

hj=Khjx

(2)

K=1.5表示耐壓結構安全系數，hjx表示極限深度。

hjx=hg/(0.85～0.90)

(3)

hg為潛水器工作下潛深度。

利用(1)式～(3)式計算得到Pj，根據規范[13]計算Mises應力σmax、軸向應力σ1、周向應力σ2和肋骨應力σ3的臨界值與許用應力σs之間的關系分別為[σmax]=σs、[σ1]=0.85σs、[σ2]=1.15σs和[σ3]=0.6σs，則上述應力的臨界值分別為785 MPa、667.25 MPa、902.75 MPa和471 MPa。

1.2 計算網格的確定

進行耐壓結構有限元模型的邊界設置時，將中縱剖面內節點的橫向線位移，以及與其環繞的2個坐標軸角位移等均設置為0，這是因為模型的左、右兩側對稱的原因。耐壓結構分析過程分為強度和穩定性兩部分內容，強度主要進行σmax、σ1、σ2和σ3等分析，穩定性完成臨界載荷Pcr的分析。在結構分析的過程中，模型的網格劃分越精細，計算精度越高，但計算時間也會增大。為此，需要進行網格無關性分析。將0.025 m、0.03 m、0.05 m、0.1 m、0.2 m 5種網格分別與0.02 m網格進行比較，強度分析過程中，主要計算σmax、σ1、σ2和σ3等，還輸出質量M。不同網格的強度和穩定性分析的計算結果和計算時間對比如表1所示。

表1 計算結果和計算時間對比

表1中，在強度校核階段，與0.02 m的網格相比，M的差別均不大，這是因為耐壓結構由若干個單元組成，M在計算過程中通過對所有單元求和得到，而單元則通過網格劃分得到，網格越小，單元越多，則單元所組成的結構越接近實際結構，求和得到的M值越接近實際值；對于σmax、σ1、σ2和σ3等而言，0.025 m網格計算結果與0.02 m網格的相對誤差不超過3%，0.03 m網格計算結果與0.02 m網格的相對誤差不超過5%，其他尺寸網格計算計算結果與0.02 m網格的相對誤差均超過7%，幾種網格的強度分析所耗時間花費差距不大，盡管0.03 m網格的誤差比0.025 m網格的誤差稍大，但誤差滿足計算要求，且在計算時間上更有優勢。在穩定性分析階段，不同尺寸網格計算得到的Pcr值均相差不大，但在計算時間方面相差較大。綜合考慮計算的精度和時間等兩個方面因素，強度校核采用0.03 m網格，穩定性分析采用0.1 m網格。

2 耐壓結構近似模型的建立

2.1 設計參數的確定

以M最小和Pcr最大作為耐壓結構優化求解的目標；約束為σmax、σ1、σ2和σ3；設計變量分別為上、下開孔圍壁的高度SG、XG，耐壓殼厚KH，上、下開孔圍壁厚SH、XH，肋骨的厚度LH、高度LG，觀察窗弧度JD，開孔加強厚度JQH和觀察窗位置WZ，耐壓結構圖和設計參數如圖2所示。

圖2 耐壓結構圖及設計參數

2.2 耐壓結構參數化分析

以耐壓結構的設計變量作為輸入，柱段長度、肋骨個數、柱段半徑為常量，Pcr、M、σmax、σ1、σ2、σ3等為輸出值。基于Abaqus軟件的耐壓結構參數化分析流程，需要實現對多個樣點自動建模和分析，需要對其進行二次開發，并解決其與多學科設計優化軟件iSight之間集成與通信等問題[16]。參數化分析的流程如圖3所示。

圖3 耐壓結構參數化分析流程

2.3 設計變量靈敏度分析

采用最優拉丁超立方方法(Opt LHD)作為實驗設計方法，選擇90組樣本點進行分析，設計變量的取值情況見表2。耐壓結構設計變量有10個，設計變量維數的降低可降低問題的分析難度，從而有利于計算時間的縮短和設計效率的提高。對90組樣本點進行靈敏度分析，可得到各設計變量對極限強度Pcr和質量M的影響比率如圖4所示。

表2 設計變量及其可行域

圖4 設計變量靈敏度分析

圖4中，對于M而言，KH的影響最大，呈正相關，其次為LH和LG，均呈正相關。對于Pcr而言，KH的影響比率最大，呈正相關，其次是LG和LH，均呈正相關。綜合考慮以上全部設計變量參數對M和Pcr的影響比率，將KH、XG、SH、XH、LG和LH等變量作為后續優化的設計變量，其余4個變量作為固定值，并用x1表示KH，用x2表示XG，用x3表示SH，用x4表示XH、用x5表示LG，用x6表示LH。

2.4 常用近似模型及精度分析

1)響應面模型(RSM)是一種多項式模型方法[17]，響應面模型的類型用階數進行區分，最高的階數為4階。4階RSM可用(4)式表示

(4)

(2)徑向基函數(RBF)神經網絡模型是一種高精度的多維空間非線性函數逼近技術[18]。當徑向基函數為高斯函數時，第i個隱藏層單元輸出的響應為

(5)

式中：i表示第i個隱藏層；ci表示中心寬度；σi表示單元實際寬度。則輸出層第j個輸出為

(6)

式中：f(x)表示RBF的函數值；w為權重值；下標i、j分別代表第i節點及該節點的第j個輸出。

(7)

式中：λi為未知的、需要待定的加權系數。

λi需要符合(8)式、(9)式，二者分別為無偏估計值和方差結果

(8)

(9)

式中：γ(xi,xj)和γ(xi,x0)分別表示xi與xj，以及x0之間距離為h時的參數半方差值。

基于耐壓結構參數化分析流程選擇樣本點進行分析，得到6個設計變量不同取值情況下響應值，利用近似模型對設計變量和響應進行擬合。近似模型能否對原有模型特性進行準確的表達，需要進行擬合精度分析，可決系數R2常被用于近似模型擬合精度的判斷，R2表達式如下：

(10)

2.5 近似模型的擬合對比及選擇

利用近似模型對樣本點的輸入和輸出進行擬合，近似模型的精度判斷結果如表3所示。

表3 近似模型擬合精度

表3中，4階RSM擬合得到的M、Pcr、σ3擬合最高，2階RSM擬合得到的σmax、σ2擬合精度最高，3階RSM擬合得到的σ1擬合精度最高。隨機選取50組樣本點將設計參數通過近似模型進行計算得到近似模型預測值，并與計算值即真實值進行比較，以展示近似模型對原有模型的特性能否做出準確描述，如圖5所示。

圖5 近似模型預測值與實際值誤差分析

圖5中，只有σ1的兩個樣本點相對誤差超過5%，σ1的其余樣本點、其他響應所有樣本點的相對誤差均小于5%，表明近似模型可對原有模型的特性做出準確描述，即近似模型具有較高的擬合精度。

3 改進非支配排序遺傳算法

(11)

式中：β表示隨機變量，對于每一維而言，都需要按照(12)式重新生成，

(12)

式中：u為隨機數，其取值范圍為[0，1]；ηc表示交叉系數，取值常數。

3.1 INSGA-Ⅱ算法

正態分布交叉算子(NDX)的空間開發能力與SBX算子相當，其空間搜索能力更強。因此，將NDX算子引入NSGA-Ⅱ算法中可提升算法的空間搜索能力[11-12]。NDX是在SBX的基礎上進行改進的，其用正態分布隨機變量A·|N(0,1)|代替SBX算子中的隨機變量β[23]。其中A為常數，其數值影響算法對空間開發和探索的概率。為使NDX算子具有與SBX算子完全一致的開發和探索概率，A取1.481。NDX算子[24]的表達式為

(13)

由于離散重組操作的引入，使新產生的子個體具有更廣泛的取值空間，算例的維度越高，算子搜索空間的提升越明顯。將NDX算子結合到NSGA-Ⅱ算法中，得到改進的非支配遺傳算法INSGA-Ⅱ，算法流程如圖6所示。

圖6 INSGA-Ⅱ算法流程圖

3.2 INSGA-Ⅱ算法性能測試

為測試INSGA-Ⅱ算法的性能，利用標準多目標測試函數ZDT1(連續性函數)見(14)式、ZDT3(非連續性函數)見(15)式進行測試[22]。

(14)

(15)

為更直觀展現算法的改進效果，以兩個測試函數的真實Pareto前沿為基準，將INSGA-Ⅱ算法與NSGA-Ⅱ算法的測試結果進行對比，對兩種算法的參數進行統一設置。其中，種群規模、最大進化代數、交叉概率、變異概率的取值分別為Pop=500、Gen=500、Pc=0.9、Pm=1/30。得到ZDT1和ZDT3的Pareto前沿如圖7、圖8所示。

圖7 算法在ZDT1上的Pareto前沿

圖8 算法在ZDT3上的Pareto前沿

通過圖7、圖8可以看出，與NSGA-Ⅱ算法相比，INSGA-Ⅱ得到的Pareto前沿更加貼近真實的前沿，其Pareto解集中的解的分布更加均勻，求解過程在收斂性方面也表現更優，表明INSGA-Ⅱ的求解能力較NSGA-Ⅱ有一定的提升。

4 耐壓結構多目標優化

4.1 多目標優化模型的建立

以M最小、Pcr最大作為優化目標，以σmax、σ1、σ2、σ3和Pcr為約束，x1、x2、x3、x4、x5和x6等變量為設計變量，結合樣本點擬合得到的近似模型，建立基于多目標的耐壓結構優化模型見(16)式：

(16)

4.2 多目標優化求解

分別采用INSGA-Ⅱ和NSGA-Ⅱ對(16)式進行求解，兩種算法的參數設置如表4所示，并將兩種算法求得的Pareto解集進行對比，如圖9所示。

表4 算法參數設置

圖9 兩種算法求得的Pareto解集對比

圖9中，NSGA-Ⅱ和INSGA-Ⅱ求解得到的Pareto解集的目標函數變化趨勢相同，即隨著M的增加Pcr也隨之增加，表明在相同算法求解的條件下，需要通過增加結構質量的方式提升其安全性。與NSGA-Ⅱ相比，在M相近的情況下，INSGA-Ⅱ的優化方案的Pcr值更大，即INSGA-Ⅱ求解得到方案性能更優，且INSGA-Ⅱ得到的Pareto解集分布范圍更大、均勻性和連續性更好。

4.3 部分優化方案對比

為驗證INSGA-Ⅱ求解耐壓結構多目標優化問題的改進效果，在兩個Pareto解集中分別選擇兩個結果相近的方案(方案1與方案3為一組，方案2與方案4為一組)進行兩兩對比，如表5所示。

由表5可知，與初始方案相比，4個方案的設計參數中，x5的變化較大，其他設計參數變化不明顯，且σmax、σ1、σ2、σ3等約束均是增加的；在滿足約束的條件下，方案1和方案3的M和Pcr較初始方案均有所降低；方案2和方案4的M均較初始方案降低的幅度要小于方案1和方案3，但Pcr的較初始方案有很大幅度提高，表明NSGA-Ⅱ和INSGA-Ⅱ均得到了優化的方案，由于不同方案優化側重或程度均有所不同，從而導致這些方案均不是最優解，這是因為在優化求解過程中，兩個優化目標M和Pcr之間是相互影響的，過度追求M的降低，可能導致Pcr的降低，在實際設計中，需要設計者根據相關要求做出選擇。與方案1相比，采用INSGA-Ⅱ算法求得的方案3，質量減少176.803 kg，Pcr增加1.81 MPa；與方案2相比，方案4的質量減少10.566 kg，Pcr增加0.635 MPa，因此，INSGA-Ⅱ求解得到方案更優，進一步驗證了INSGA-Ⅱ的有效性。

表5 多目標優化方案對比

5 結論

1)通過網格無關性分析確定了參數化分析的網格尺寸，實現了設計變量的篩選和降維，得到了滿足要求的近似模型，以上工作可協調計算精度和計算成本之間的矛盾，進一步提升了耐壓結構設計效率。

2)針對傳統NSGA-Ⅱ在求解多目標優化時所存在的解集分布性不足的問題，以NDX算子作為交叉算子對NSGA-Ⅱ進行了改進,選擇典型測試函數對改進算法INSGA-Ⅱ進行了測試。結果表明,INSGA-Ⅱ具有處理連續性問題和非連續問題的能力，且與傳統NSGA-Ⅱ算法相比，結果的收斂性和分布性均有提升，證明了INSGA-Ⅱ的可行性。

3)利用INSGA-Ⅱ進行了多開孔耐壓結構的多目標優化求解，選取部分方案與NSGA-Ⅱ得到的優化方案進行了對比，通過對比表明：INSGA-Ⅱ得到的優化方案優化效果更優，在耐壓結構性能提升方面更加顯著。