置信域優(yōu)化法求解船舶浮態(tài)方程

張 弢

（巴斯大學(xué) 理學(xué)院 英國 巴斯）

引 言

船舶的浮態(tài)平衡問題是靜力學(xué)中最基本的問題，也是任何浮體在設(shè)計過程中首先考慮的問題之一，中外學(xué)者多年來均對該問題進行了深入研究。傳統(tǒng)的船舶靜力學(xué)原理基于等體積假設(shè)，即體積微幅傾斜時，浮體旋轉(zhuǎn)軸必定通過當前水線面的漂心位置。趙曉非［1］根據(jù)靜力學(xué)原理建立平衡方程并使用牛頓法迭代求解，但每次迭代需要計算包含水線面慣性矩等參數(shù)的雅可比矩陣。桑松［2］在此基礎(chǔ)上利用最小功原理將該方法擴展到適用于任何浮體的變軸橫傾穩(wěn)性上，并推導(dǎo)出浮式結(jié)構(gòu)空間穩(wěn)性臂計算公式。此外，孫承猛［3］使用有限差分法近似計算雅可比矩陣，避免了復(fù)雜的公式推導(dǎo)。求方程的根其實是求函數(shù)極值的特殊情況，可以轉(zhuǎn)化為數(shù)值優(yōu)化中求函數(shù)最小值的問題。因此，馬坤［4］基于非線性規(guī)劃法設(shè)定懲罰函數(shù)來求目標函數(shù)的最小值，并將排水體積和浮心表示為3個自變量——吃水、橫傾角和縱傾角的函數(shù)，進而避免計算水線面慣性矩等參數(shù)。之后，陸叢紅［5］使用約束優(yōu)化的遺傳算法來求解浮態(tài)方程；金寧［6］在此基礎(chǔ)上進一步優(yōu)化遺傳算法，使目標方程收斂更快；張維英［7］則使用差分進化法來求解優(yōu)化問題。

其實，浮態(tài)平衡方程是二階可導(dǎo)的非線性方程組，使用非約束優(yōu)化（比如線搜索法，置信域法）就可以求解；遺傳算法、差分進化法等卻更適合求解約束優(yōu)化的全局最小值問題，并且由于不需要計算梯度，目標方程可以是非連續(xù)不可微的函數(shù)，因此被廣泛應(yīng)用在生物和經(jīng)濟學(xué)領(lǐng)域（若用于解決船舶浮態(tài)平衡問題，則比牛頓法和置信域法的超線性收斂速度要慢很多）。

本文前4節(jié)將建立浮態(tài)平衡的數(shù)學(xué)模型并論述兩種計算排水量和浮心的通用方法，第5和第6節(jié)介紹置信域Dogleg算法（又稱“狗腿算法”），并通過MATLAB算例驗證該方法用于求解浮態(tài)方程的實際效果。

1 坐標系

圖1 正浮狀態(tài)下的坐標系

2 歐拉角和線性變換

剛體和隨剛體運動的局部坐標系OX′Y′Z′相對于靜止的全局坐標系OXYZ的方位可以用3個歐拉角來表示，但同樣的方位可以用不同組合的3個歐拉角表示。因此，船舶工程中通常將繞Z′軸的轉(zhuǎn)動稱為首搖（γ），繞Y′軸的轉(zhuǎn)動稱為縱搖（β），繞X′軸的轉(zhuǎn)動稱為橫搖（α），并規(guī)定順序為首搖-縱搖-橫搖。在船舶靜力學(xué)中，首搖并不會改變船體的穩(wěn)性，因此歐拉角可以省略為2個。在本文中，歐拉角的正負號按照右手螺旋法則定義，即當船首下沉時縱搖為正，右舷下沉時橫搖為正。

船體從正浮狀態(tài)（OX′Y′Z′與OXYZ重合）縱搖，可以看作是3空間中的線性變換，用矩陣P表示：

船體從正浮狀態(tài)橫搖的線性變換用矩陣R表示為：

根據(jù)定義的歐拉角順序（即先繞Y′軸縱搖再繞X′軸橫搖），合并后的線性變換表示為矩陣T：

式中：下標G表示基底為的線性變換，下標L表示基底為的線性變換。因Y′軸與Y軸在正浮時重合，所以［P］G=［P］L。從式（1）中可以看出先以Y′軸縱搖再以X′軸橫搖等效于先以X軸橫搖再以Y軸縱搖。由于此線性變換只發(fā)生在3空間內(nèi)，因此也可以理解為局部坐標不變進行基底變換，求得的矩陣相同。至此，船體傾斜后的型線坐標（x，y，z）可由式（2）求出：

式中：（x0，y0，z0）為船體正浮的型線坐標。

3 浮態(tài)平衡方程

船舶靜力學(xué)的模型通常只涉及到垂蕩、橫搖、縱搖3個自由度，如圖2根據(jù)牛頓定律可以表示為由3個未知量以及3個等式組成的方程組參見式（3）：

其中：

水密度ρ、船體排水量m為已知量，水線面的Z軸坐標WLz、縱傾角β和橫傾角α是未知量，重心坐標CG是β和α的函數(shù)。由此，式（2）便可表示為該式中：（CGx0CGy0CGz0）為船體正浮時的重心坐標；下標x、y和z表示全局參考系OXYZ下的坐標。除特別說明以外，之后所有的坐標均參考全局坐標系OXYZ。

圖2 浮態(tài)平衡

排水體積V和浮心坐標CB是WLz、β以及α的函數(shù)。本文所用方法不基于等體積假設(shè)，因此坐標系原點無需定于水線面漂心處，而是定于如圖1和圖2所示的船尾中心基線處。這樣做的好處是：

（1）和船體型值表的參考系一致，不需進行額外的坐標轉(zhuǎn)換。

（2）對船長遠大于船寬的普通船舶來說，縱傾對浮心變化的影響遠遠大于橫傾，因此將縱傾的旋轉(zhuǎn)軸設(shè)在船尾可以減小浮心的變化速度。用數(shù)值優(yōu)化方法求解時，合理的縮放目標函數(shù)可以提高收斂的穩(wěn)定性。

（3）旋轉(zhuǎn)軸的平移并不會影響平衡后的橫傾和縱傾角度，只會影響水線面高度WLz。而此模型中，因水線面不受約束，故可在全局坐標系Z軸上任意平移，而對最終計算結(jié)果沒有影響。

4 浮心和排水體積

船體模型通過式（2）線性變換后與水線面相交，水線面切割模型（通常由離散的多邊形組成，見下頁圖3、圖4）后得到水下的部分，其算法的原理是找到水線面和多邊形的交點并舍棄水線以上部分，再按順序重新連接交點和多邊形在水下的頂點。接下來，計算式（3）中的排水體積V和浮心坐標CB主要有面元法［8］和切片法。

4.1 面元法

如下頁圖3所示，船體表面可劃分成很多平面單元，對水下單元的表面積分計算其浮力、浮力矩和體積矩等，然后用斯托克斯定理（或格林定理）將表面積分轉(zhuǎn)換為單元輪廓的封閉曲線積分，見式（4）。由于單元的輪廓為多邊形，轉(zhuǎn)換為曲線積分以后，只需對每條邊的直線方程積分即可。通過這個原理對所有水下單元的積分求和，即可得到總的排水體積和浮心。例如，當式 （4） 中curl F·n取常數(shù)1時，表面積分為面積；當curl F·n取z·nz時，表面積分為體積；當curl F·n取0.5·z2·nz時，表面積分為對Z軸的體積矩。

圖3 面元模型（Heimdal小型工程船）

圖4 切片模型（Heimdal小型工程船）

式中：curl F為向量場F的旋度；n（nx，ny，nz）為平面單元的法向量；dr為曲線的切向量。

4.2 切片法

將船體切割成很多與船長方向垂直的切片（見圖4），計算每個切片水下區(qū)域的面積Aw和形心dCB，然后使用式（5）和式（6）進行數(shù)值積分（如梯形法）求出排水體積和浮心，其中Lwl為水線長。

文獻［4］按照上述面元法的思路，以格林定理來求Aw和dCB，但本文給出一種更簡潔的表述方法，并可直接應(yīng)用于三維平面。在3中有任意兩個向量q1（xq1，yq1，zq1），q2（xq2，yq2，zq2），其叉積的長度為：

θ是q1逆時針旋轉(zhuǎn)到q2的角度。式（7）可以看作是q1和q2圍成的三角形q1Oq2面積的2倍。因此，任意多邊形的每個頂點可以表示為該平面原點處的向量（參見下頁圖5）并讓最后一個頂點和第一個頂點重合（q1，q2，…qn，q1），按順序?qū)⑾噜徬蛄坎娉?，則多邊形的面積為：

當 π＜θ≤ 2π 時，ε= 1 ；當 0 ≤θ≤ π 時，ε= 2。

dCB的計算原理類似，根據(jù)相似三角形和中點定理推導(dǎo)出三角形的形心公式（8）。計算圖5中每個三角形qi Oqi+1的形心及面積矩，則整個多邊形的形心為可由式（9）求出。

圖5 求多邊形的面積和形心

面元法和切片法有各自的優(yōu)勢：面元法需要先將船體模型劃分成網(wǎng)格模型，但網(wǎng)格模型可以和后續(xù)的有限元結(jié)構(gòu)分析以及水動力分析共享，減少后續(xù)工作量；切片法可以直接使用型值表或型線圖的數(shù)據(jù)進行穩(wěn)性計算，減少前期工作量。當切片數(shù)量或網(wǎng)格數(shù)量足夠多時，兩種方法都能準確地反映船體形狀，因此都是可靠穩(wěn)定的方法。

5 求解非線性方程組

式（3）可以寫成向量形式如下：

f是一個連續(xù)非線性方程組，求解式（10）相當于求解連續(xù)函數(shù)l無約束下的極小值問題：

解決式（11）的基本思想是給出一個初始值x0，然后迭代求使：

直到發(fā)現(xiàn)極值或者找不出新的xk為止。式中：下標k為迭代次數(shù)；pk是單位方向向量；τ是步長。

求無約束函數(shù)極值問題主要有兩類方法：線搜索法和置信域法。線搜索法（比如梯度下降法，牛頓法等）基本思想是先找到搜索方向pk滿足式（12），然后再求出合適的步長τ；而置信域法是先確定一個搜索半徑Δ，然后在Δ內(nèi)同時找到滿足式（12）的pk和τ。如果找不到，則縮小Δ直到找到滿足式（12）的pk和τ為止。

找出pk和τ的方法有多種，本文使用的置信域狗腿法具體原理［9］如下。

根據(jù)泰勒定理將目標函數(shù)l在xk附近近似為二次函數(shù)，如式（13）所示：

當Hk為正定矩陣時，式（14）的解如圖6所示，可分為以下3種情況。

圖6 狗腿算法求p

式中：pB為無約束條件下，式（14）的全局解。

式中：pU為無約束條件下，式（14）在-gk方向的解，

式中：τ為二次方程的解。

這樣做的好處是在搜索半徑內(nèi)找不到式（14）的全局解時，盡可能利用二次導(dǎo)數(shù)信息Hk，在P U和P B的連線上找式（14）的近似解p，參見圖6（c）。而當搜索半徑相對較小時，式（13）可去除二次項，簡化為線性近似，在梯度的負方向求式（14）的近似解p，參見圖6（b）。

當Hk為非正定矩陣時，可以將Hk對角化，并對特征值取絕對值，將Hk轉(zhuǎn)為正定矩陣。無法直接求出Hk時，可用有限差分法法求近似Hk。

搜索半徑Δ的選取通常依據(jù)目標函數(shù)l和二次函數(shù)m的擬合比率r：

圖7 置信域狗腿算法偽代碼

6 收斂結(jié)果

為測試上述置信域狗腿法用于求解浮態(tài)方程的實際收斂效果，用MATLAB編寫船舶穩(wěn)性程序分別以牛頓法（N）和置信域狗腿法（T）計算兩艘實船在不同工況下的浮態(tài)（水線面高度WLz，橫傾角α和縱傾角β），并比較最后目標方程的冗余值‖f‖2和迭代次數(shù)。船舶的排水量和重心坐標為已知，同時給出迭代的初始值WLz0、α0和β0。

首先以圖3和圖4所示Heimdal小型工程船為例。該船主尺度為： 總長36 m、寬8.2 m、高3.2 m。

從下頁表1中的計算結(jié)果可以看出：工況1的兩種方法計算結(jié)果相同，船體處于接近正浮的狀態(tài)，收斂的精度和迭代次數(shù)也基本相當；工況2的排水量增加，重心向船首右舷移動并升至2 m。若使用和工況1相同的初始值時，兩種方法的計算結(jié)果都很理想（橫傾23°，縱傾4°），但使用第3組初始值（0.3，0，0）時，牛頓法會收斂到另一個平衡位置（橫傾-88°，縱傾7°）。雖然這個解也滿足平衡方程，但與實際不符。究其原因是牛頓法并沒有限制步長，在某一步迭代中跨過了最近的解，并收斂到另一個遠處的解；置信域法由于步長限制在搜索半徑內(nèi)，因此計算結(jié)果正常。使用最后一組初始值時，牛頓法在迭代至第4步時失效，原因是該步算出的雅可比矩陣不可逆，因此無法算出這一步的解；置信域法由于保證在搜索半徑內(nèi)找到解，參見式（14）-式（17），因此計算結(jié)果正常。

表1 Heimdal小型工程船浮態(tài)計算結(jié)果

下面再以某大型集裝箱船為例。該船主尺度為：總長288 m、寬32.2 m、高21.6 m。從表2所示計算結(jié)果可以看出：工況1時，初始值為（10，0，0），兩種方法計算結(jié)果相同，但置信域法的迭代次數(shù)較多。這是由于當二次函數(shù)與目標函數(shù)擬合情況較好時，置信域狗腿算法會逐漸增大搜索半徑（如圖6增大1倍）；牛頓法由于沒有限制步長，因此收斂更迅速，但同時也存在風險，如表1中的第3組例子。工況2時，重心向船尾右舷移動，初始值不變，但牛頓法再一次求解失敗，原因和表1中的第4組例子一樣，雅可比矩陣不可逆。

表2 某80 000 t集裝箱船浮態(tài)計算結(jié)果

需要注意的是：本文采用的牛頓法沒有限制步長是為保留其特點并突出和置信域法的區(qū)別，實際中也可以使用回溯法改進牛頓法達到控制步長的目的，但仍然存在雅可比矩陣可能是奇異矩陣的問題。

7 結(jié) 論

從實例的收斂結(jié)果來看，牛頓法較為依賴初始值的選擇，因為雅可比矩陣可能產(chǎn)生不可逆的問題，導(dǎo)致求解中斷；而相比于牛頓法，置信域法不僅更穩(wěn)定、對初始值的依賴程度低，且能確保收斂，因此非常適合用于工業(yè)計算及穩(wěn)性軟件開發(fā)。

此外，求解浮態(tài)平衡方程的原理也可以推廣到其他穩(wěn)性相關(guān)的問題。比如計算復(fù)原力臂，只需將式（3）中的第3個方程去掉，已知橫傾角α求解船體在Z軸和X軸方向的平衡方程，然后計算浮心和重心的Y軸坐標差即可，這種方法同時考慮了自由縱傾的影響。當艙室破損時，可以把進水看成附加質(zhì)量，按照第4節(jié)的方法計算艙室進水的質(zhì)量和重心，疊加到船體本身的質(zhì)量、重心，再求解平衡方程即可。再比如計算海洋平臺的變軸橫傾穩(wěn)定力臂，可以修改第2節(jié)的線性變換矩陣，使船體沿著指定方位角橫傾即可。