水下拖纜穩(wěn)態(tài)運動的多目標優(yōu)化研究

王沖霄，劉忠樂，文無敵，張志強，趙 苗

(海軍工程大學 兵器工程學院， 武漢 430033)

水下拖曳系統(tǒng)，是一種廣泛應(yīng)用于海洋監(jiān)測、海洋研究以及軍事等領(lǐng)域的水下探測裝置，在開發(fā)海洋的先進技術(shù)手段中，拖曳系統(tǒng)裝備具有極其重大的意義[1]。線列陣是水下拖曳系統(tǒng)的一個關(guān)鍵組成部分，其水動力特性直接影響整個水下航行器系統(tǒng)的快速性、操縱性及穩(wěn)定性，開展水下拖纜的水動力特性研究具有重要的理論意義和工程實用價值。近年來，基于響應(yīng)面的優(yōu)化方法已廣泛應(yīng)用在穩(wěn)健設(shè)計和多目標與多學科優(yōu)化設(shè)計的代理模型里[2]。這種近似模型技術(shù)是在初始數(shù)據(jù)集合基礎(chǔ)上構(gòu)造逼近目標函數(shù)和約束條件的方法，同時也為快速優(yōu)化和敏感性分析提供了一種高效的解決方法[3]。本文針對水下拖纜，在王飛[4-5]的穩(wěn)態(tài)運動求解和分析的基礎(chǔ)上，引入多目標優(yōu)化設(shè)計理論，將復(fù)雜求解方法所得結(jié)果進行回歸處理，建立起拖纜尾端拖曳深度和首端張力的二次響應(yīng)面模型，給出多目標優(yōu)化算法下的Pareto最優(yōu)解集，并分析了拖纜參數(shù)對尾端拖曳深度和首端張力的影響。

1 拖纜運動模型及數(shù)值計算

為簡化討論過程，本文僅討論穩(wěn)定海流下自由端拖纜的運動模型，其拖曳系統(tǒng)如圖1。

圖1 拖曳系統(tǒng)示意圖

其他的拖曳形式，例如尾端加載拖船、尾繩、水下拖體及其他的纜載設(shè)備，相比自由端僅在邊界條件上存在區(qū)別，此處不詳細討論。將拖纜視為理想的柔性圓形纜繩，由水下航行器(AUV)搭載，建立空間固定的慣性坐標系O-XYZ，單位矢量定義為(i,j,k)，附著在拖纜上的局部坐標系btn，單位矢量定義為(b,t,n)。軸t表示拖纜切向，方向為纜長s的增長方向；軸n表示拖纜的法向，處在軸t和軸t在OXY平面內(nèi)的投影所組成的平面內(nèi)，并垂直于軸t；軸b與軸n和軸t共同組成右手笛卡爾坐標系。定義歐拉角θ，φ為拖纜微元相對慣性坐標系的姿態(tài)角，θ為Otn平面偏離OX軸的角度，φ為軸t偏離OXY平面的角度，θ∈(-180°,180°],φ∈(-90°，90°]，兩個歐拉角均以逆時針方向為正方向。慣性坐標系和局部坐標系通過姿態(tài)角相關(guān)聯(lián)，轉(zhuǎn)換關(guān)系如下：

(1)

(2)

對于每個拖纜微元ds，在穩(wěn)定直航的狀態(tài)下，其重力、浮力、流體阻力的合力處于平衡狀態(tài)，有以下平衡方程：

(3)

式中:T代表拖纜張力，始終指向拖纜的切向，B和G分別代表拖纜單位長度的浮力和重力，D代表流體阻力。

單位矢量(b,t,n)對拖纜長度s微分用“′”表示，則有[6]：

(4)

則可得：

(5)

所以張力在局部坐標系可展開為下式：

(6)

將式(6)與重力、浮力和流體阻力代入平衡方程式(3)中，并在拖纜局部坐標系下沿各坐標軸方向展開，則平衡方程寫為如下標量形式[7]：

(7)

式中:w為拖纜單位長度在水中的質(zhì)量，表示為：w=(u-ρσ)g，其中u為拖纜單位長度質(zhì)量，ρ為流體密度，σ為橫截面積，d為拖纜直徑，ε為拖纜應(yīng)變，Ct和Cn分別為拖纜的切向和法向阻力系數(shù)，ut，ub，un為局部坐標系下的速度分量，由系統(tǒng)相對于水流的拖曳速度轉(zhuǎn)換到局部坐標系下得到，有：

(8)

式中:v為拖曳速度，J為海流，僅有水平面內(nèi)速度分量而無垂向分量。拖纜在慣性系下的坐標由下式給出：

(9)

在均勻海流狀態(tài)下，拖纜自由端即s=0處的張力和歐拉角變化率均為零，拖纜側(cè)向無作用力存在，有θ≡const，此時可將穩(wěn)態(tài)問題轉(zhuǎn)換到二維空間下求解，在自由端的邊界條件表示如下[4]：

(10)

式中ψ為AUV的航行艏向角，將以上初始值方程和控制方程(7)、(9)聯(lián)立，采用四階龍格庫塔方法進行積分計算，即可求出穩(wěn)態(tài)解。

2 近似模型的建立

正常情況下，拖曳系統(tǒng)均工作于穩(wěn)定狀態(tài)，而穩(wěn)態(tài)下的拖曳深度和拖纜首端張力是設(shè)計過程中兩個非常重要的指標，拖纜的密度、楊氏模量、阻力系數(shù)以及拖曳速度和流體密度，均會對此產(chǎn)生一定的影響。王飛[5]對穩(wěn)態(tài)拖曳各參數(shù)在一定的取值范圍下，進行了拖曳深度和拖纜首端張力的初步研究，但每個參數(shù)的變化，是在其余參數(shù)不變的基礎(chǔ)上進行的，而且最后僅對各參數(shù)對結(jié)果的影響進行了定性分析。以下將采用試驗設(shè)計的方法，在王飛給出的拖纜模型基礎(chǔ)上，構(gòu)建參數(shù)和響應(yīng)之間的近似模型，并得出尾端拖曳深度和首端張力的近似計算公式。

2.1 試驗設(shè)計

本文選擇拉丁超立方設(shè)計(Latin hypercube designs)作為試驗方法。拉丁超立方設(shè)計方法是一種優(yōu)秀的強調(diào)樣本點分布均勻的試驗設(shè)計方法[8]。該方法假設(shè)需要n個實驗設(shè)計點，則設(shè)計變量會被分成n等分，在每等分中選擇一個參數(shù)作為設(shè)計點。根據(jù)問題規(guī)模和復(fù)雜程度，樣本點的個數(shù)也應(yīng)該適當增多。通常情況下，對于5～10個變量的問題，樣本點數(shù)量建議取為1.5×(n+1) ×(n+2)/2個[9]，因此針對本文的6個參數(shù)變量，選擇42個樣本點進行試驗設(shè)計。拖纜模型和參數(shù)取值范圍采用文獻[5]給出的模型，參數(shù)如表1所示：

表1 拖纜參數(shù)

尾繩置于引導纜尾端，起到定深和穩(wěn)定的作用。針對引導纜的物理參數(shù)以及流體密度ρ和拖曳速度vy設(shè)定取值范圍如下：

108N/m2≤E≤1011N/m2

0.010≤Ct≤0.03

1.2≤Cn≤1.9

0.8 kg/m≤u≤1.2 kg/m

1 m/s≤vy≤3 m/s

1 020 kg/m3≤ρ≤1 030 kg/m3

2.2 二階多項式響應(yīng)面

二階多項式響應(yīng)面的數(shù)學表達式為：

(11)

采用復(fù)相關(guān)系數(shù)R2作為響應(yīng)面模型的誤差分析指標，定義如下：

(12)

采用多項式回歸技術(shù)對試驗設(shè)計的樣本點和響應(yīng)值進行最小二乘擬合，并求出待定系數(shù)，構(gòu)造出近似模型?；貧w模型表示如下：

yi=f(xi,θ)+εi，i=1,2,…，n

(13)

Z=33 539.7-258.38x1-2 728.3x2-1.26×10-9x3-

87.61x1x2-1.21×10-11x1x3+

0.19x1x4-27.43x1x5+12.28x1x6-

3.04×10-11x2x3-3.07x2x4+

121.57x2x5-0.47x2x6+1.22×10-12x3x4+

5.44×10-12x3x5-3.12×10-11x3x6-

0.19×10-3x4x5-0.21x4x6-52.96x5x6

(14)

F=-429 504.9-1 752.43x1-84 097.73x2-

9.14×10-9x3-842.3x4+3 016.54x5-

2.03x1x4-241.36x1x5+40.65x1x6+

3.1×10-9x2x3-52.68x2x4-5 144.58x2x5+

47 585.52x2x6+9.16×10-12x3x4+

5.12×10-11x3x5+1.54×10-11x3x6-

1.26x4x5+1.26x4x6-597.65x5x6

(15)

式中變量(x1,x2,x3,x4,x5,x6)分別代表參數(shù)中的Cn，Ct，E，ρ，u，vy。通過計算復(fù)相關(guān)系數(shù)來驗證響應(yīng)面模型的擬合精度。其中，Z的復(fù)相關(guān)系數(shù)為0.996 96，F(xiàn)的復(fù)相關(guān)系數(shù)為0.999 53。該響應(yīng)面模型擬合程度較好。

3 多目標遺傳算法

本文所探討的是如何實現(xiàn)尾端拖曳深度最大、首端張力最小的多目標優(yōu)化問題，則本優(yōu)化問題的數(shù)學模型可以描述為：

minF, maxZ

findCn,Ct,E,ρ,u,vy

s.t. 108N/m2≤E≤1011N/m2

0.010≤Ct≤0.03

1.2≤Cn≤1.9

0.8 kg/m≤u≤1.2 kg/m

1 m/s≤vy≤3 m/s

1 020 kg/m3≤ρ≤1 030 kg/m3

(16)

大多數(shù)情況下，上式的各個子目標往往是互相沖突的，其中一個子目標的優(yōu)化會帶來其他子目標的損失，所以多目標優(yōu)化問題的優(yōu)化解是一個解集，稱為Pareto最優(yōu)解集，解集中的元素稱為Pareto最優(yōu)解，Pareto最優(yōu)解集在目標函數(shù)空間中的像稱為Pareto前沿。在目前流行的優(yōu)化算法中，遺傳算法最適合求解多目標優(yōu)化問題的Pareto解集。本文采用遺傳算法中的NSGA-Ⅱ算法，基于ISIGHT平臺，完成整個優(yōu)化過程。NSGA-Ⅱ算法利用基于Pareto支配的排序方法將個體進行分層排序，并通過計算擁擠距離的方式對同一層級個體進行具體排序，具有較高的運算效率和較好的收斂速度[10]。本優(yōu)化問題的主要流程如下：

1) 初始化拖纜的二階多項式響應(yīng)面計算模型，調(diào)入所需參數(shù)，初始化進化過程的參數(shù)，并隨機生成初始種群N；

2) 將初始種群N中的個體依次賦值給拖纜計算模型，修改待優(yōu)化參數(shù)并運行該模型，輸出對應(yīng)的拖曳深度和張力值，判斷是否滿足約束條件，以懲罰系數(shù)來處理不滿足條件的個體，使其在進化過程中被淘汰；

3) 對種群中所有個體均以此方法進行操作，得到對應(yīng)的輸出值，從而完成種群的初始化工作；

4) 對初始種群N進行非支配排序并計算擁擠距離，得出每個個體的優(yōu)劣性指標；

5) 以t記錄進化次數(shù)，初始值為1，對父種群Nt執(zhí)行交叉、變異等進化操作，產(chǎn)生子代種群Na；

6) 對每個子代種群Na進行計算，得到對應(yīng)目標值；

7) 對父代種群Nt和子代種群Na進行非支配排序并計算擁擠距離；

8) 對種群進行更新，對上一步驟合并的種群進行修正，得到新的種群Nt+1，判斷是否達到最大代數(shù)，如果是則輸出非支配解集，否則繼續(xù)迭代。

4 計算結(jié)果

本優(yōu)化過程的種群規(guī)模為24，遺傳代數(shù)為50，交叉概率為0.9，所得深度和張力的Pareto前沿散點分布如圖2，選取3個優(yōu)化方案，方案1中Z最大，方案2中F最小，方案3位于中間，所得優(yōu)化結(jié)果如表2所示。Pareto前沿為非支配解集，最終方案的選擇依賴于實際情況和選擇者偏好[11]。

圖2 拖纜優(yōu)化的Pareto前沿散點分布

表2 Pareto解集中的三個優(yōu)化結(jié)果

本文同時在ISIGHT平臺上通過數(shù)據(jù)處理給出各參數(shù)對尾端拖曳深度和首端張力影響的Pareto柱狀圖，如圖3和圖4所示。圖3的第一條、第三條和第六條以及圖4的第四條和第六條表示該參數(shù)與響應(yīng)呈負相關(guān)關(guān)系，其余表示該參數(shù)與響應(yīng)呈正相關(guān)關(guān)系?？梢钥闯?，拖曳速度vy對尾端拖曳深度和張力的貢獻程度百分比最大，分別達到了-52.53%和62.56%；其余各參數(shù)中，拖纜密度u和法向阻力系數(shù)Cn對尾端拖曳深度貢獻程度百分比較大，分別為32.47%和-12.24%，切向阻力系數(shù)Ct和拖纜密度u對張力貢獻程度百分比較大，分別為22.18%和12.12%，而楊氏模量E和流體密度ρ對兩者的貢獻程度均可忽略不計。以上分析結(jié)果基本和文獻[5]所得結(jié)論一致。

圖3 參數(shù)對拖曳深度影響

圖4 參數(shù)對首端張力影響

5 結(jié)論

基于ISIGHT平臺實現(xiàn)了參數(shù)的多目標優(yōu)化并得到Pareto解集，為設(shè)計者提供了重要的設(shè)計信息；通過計算各參數(shù)對尾端拖曳深度和首端張力的影響，驗證了前人的結(jié)論。若拖纜模型和參數(shù)取值范圍有變，本文所得出的尾端拖曳深度和張力的近似計算公式相應(yīng)會發(fā)生變化。