水下拖纜穩態運動的多目標優化研究

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

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

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

1 拖纜運動模型及數值計算

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

圖1 拖曳系統示意圖

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

(1)

(2)

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

(3)

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

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

(4)

則可得：

(5)

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

(6)

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

(7)

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

(8)

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

(9)

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

(10)

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

2 近似模型的建立

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

2.1 試驗設計

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

表1 拖纜參數

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

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 二階多項式響應面

二階多項式響應面的數學表達式為：

(11)

采用復相關系數R2作為響應面模型的誤差分析指標，定義如下：

(12)

采用多項式回歸技術對試驗設計的樣本點和響應值進行最小二乘擬合，并求出待定系數，構造出近似模型。回歸模型表示如下：

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)分別代表參數中的Cn，Ct，E，ρ，u，vy。通過計算復相關系數來驗證響應面模型的擬合精度。其中，Z的復相關系數為0.996 96，F的復相關系數為0.999 53。該響應面模型擬合程度較好。

3 多目標遺傳算法

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

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)

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

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

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

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

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

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

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

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

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

4 計算結果

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

圖2 拖纜優化的Pareto前沿散點分布

表2 Pareto解集中的三個優化結果

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

圖3 參數對拖曳深度影響

圖4 參數對首端張力影響

5 結論

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