置信域优化法求解船舶浮态方程

张 弢

（巴斯大学 理学院 英国 巴斯）

引 言

船舶的浮态平衡问题是静力学中最基本的问题，也是任何浮体在设计过程中首先考虑的问题之一，中外学者多年来均对该问题进行了深入研究。传统的船舶静力学原理基于等体积假设，即体积微幅倾斜时，浮体旋转轴必定通过当前水线面的漂心位置。赵晓非［1］根据静力学原理建立平衡方程并使用牛顿法迭代求解，但每次迭代需要计算包含水线面惯性矩等参数的雅可比矩阵。桑松［2］在此基础上利用最小功原理将该方法扩展到适用于任何浮体的变轴横倾稳性上，并推导出浮式结构空间稳性臂计算公式。此外，孙承猛［3］使用有限差分法近似计算雅可比矩阵，避免了复杂的公式推导。求方程的根其实是求函数极值的特殊情况，可以转化为数值优化中求函数最小值的问题。因此，马坤［4］基于非线性规划法设定惩罚函数来求目标函数的最小值，并将排水体积和浮心表示为3个自变量——吃水、横倾角和纵倾角的函数，进而避免计算水线面惯性矩等参数。之后，陆丛红［5］使用约束优化的遗传算法来求解浮态方程；金宁［6］在此基础上进一步优化遗传算法，使目标方程收敛更快；张维英［7］则使用差分进化法来求解优化问题。

其实，浮态平衡方程是二阶可导的非线性方程组，使用非约束优化（比如线搜索法，置信域法）就可以求解；遗传算法、差分进化法等却更适合求解约束优化的全局最小值问题，并且由于不需要计算梯度，目标方程可以是非连续不可微的函数，因此被广泛应用在生物和经济学领域（若用于解决船舶浮态平衡问题，则比牛顿法和置信域法的超线性收敛速度要慢很多）。

本文前4节将建立浮态平衡的数学模型并论述两种计算排水量和浮心的通用方法，第5和第6节介绍置信域Dogleg算法（又称“狗腿算法”），并通过MATLAB算例验证该方法用于求解浮态方程的实际效果。

1 坐标系

图1 正浮状态下的坐标系

2 欧拉角和线性变换

刚体和随刚体运动的局部坐标系OX′Y′Z′相对于静止的全局坐标系OXYZ的方位可以用3个欧拉角来表示，但同样的方位可以用不同组合的3个欧拉角表示。因此，船舶工程中通常将绕Z′轴的转动称为首摇（γ），绕Y′轴的转动称为纵摇（β），绕X′轴的转动称为横摇（α），并规定顺序为首摇-纵摇-横摇。在船舶静力学中，首摇并不会改变船体的稳性，因此欧拉角可以省略为2个。在本文中，欧拉角的正负号按照右手螺旋法则定义，即当船首下沉时纵摇为正，右舷下沉时横摇为正。

船体从正浮状态（OX′Y′Z′与OXYZ重合）纵摇，可以看作是3空间中的线性变换，用矩阵P表示：

船体从正浮状态横摇的线性变换用矩阵R表示为：

根据定义的欧拉角顺序（即先绕Y′轴纵摇再绕X′轴横摇），合并后的线性变换表示为矩阵T：

式中：下标G表示基底为的线性变换，下标L表示基底为的线性变换。因Y′轴与Y轴在正浮时重合，所以［P］G=［P］L。从式（1）中可以看出先以Y′轴纵摇再以X′轴横摇等效于先以X轴横摇再以Y轴纵摇。由于此线性变换只发生在3空间内，因此也可以理解为局部坐标不变进行基底变换，求得的矩阵相同。至此，船体倾斜后的型线坐标（x，y，z）可由式（2）求出：

式中：（x0，y0，z0）为船体正浮的型线坐标。

3 浮态平衡方程

船舶静力学的模型通常只涉及到垂荡、横摇、纵摇3个自由度，如图2根据牛顿定律可以表示为由3个未知量以及3个等式组成的方程组参见式（3）：

其中：

水密度ρ、船体排水量m为已知量，水线面的Z轴坐标WLz、纵倾角β和横倾角α是未知量，重心坐标CG是β和α的函数。由此，式（2）便可表示为该式中：（CGx0CGy0CGz0）为船体正浮时的重心坐标；下标x、y和z表示全局参考系OXYZ下的坐标。除特别说明以外，之后所有的坐标均参考全局坐标系OXYZ。

图2 浮态平衡

排水体积V和浮心坐标CB是WLz、β以及α的函数。本文所用方法不基于等体积假设，因此坐标系原点无需定于水线面漂心处，而是定于如图1和图2所示的船尾中心基线处。这样做的好处是：

（1）和船体型值表的参考系一致，不需进行额外的坐标转换。

（2）对船长远大于船宽的普通船舶来说，纵倾对浮心变化的影响远远大于横倾，因此将纵倾的旋转轴设在船尾可以减小浮心的变化速度。用数值优化方法求解时，合理的缩放目标函数可以提高收敛的稳定性。

（3）旋转轴的平移并不会影响平衡后的横倾和纵倾角度，只会影响水线面高度WLz。而此模型中，因水线面不受约束，故可在全局坐标系Z轴上任意平移，而对最终计算结果没有影响。

4 浮心和排水体积

船体模型通过式（2）线性变换后与水线面相交，水线面切割模型（通常由离散的多边形组成，见下页图3、图4）后得到水下的部分，其算法的原理是找到水线面和多边形的交点并舍弃水线以上部分，再按顺序重新连接交点和多边形在水下的顶点。接下来，计算式（3）中的排水体积V和浮心坐标CB主要有面元法［8］和切片法。

4.1 面元法

如下页图3所示，船体表面可划分成很多平面单元，对水下单元的表面积分计算其浮力、浮力矩和体积矩等，然后用斯托克斯定理（或格林定理）将表面积分转换为单元轮廓的封闭曲线积分，见式（4）。由于单元的轮廓为多边形，转换为曲线积分以后，只需对每条边的直线方程积分即可。通过这个原理对所有水下单元的积分求和，即可得到总的排水体积和浮心。例如，当式 （4） 中curl F·n取常数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），计算每个切片水下区域的面积Aw和形心dCB，然后使用式（5）和式（6）进行数值积分（如梯形法）求出排水体积和浮心，其中Lwl为水线长。

文献［4］按照上述面元法的思路，以格林定理来求Aw和dCB，但本文给出一种更简洁的表述方法，并可直接应用于三维平面。在3中有任意两个向量q1（xq1，yq1，zq1），q2（xq2，yq2，zq2），其叉积的长度为：

θ是q1逆时针旋转到q2的角度。式（7）可以看作是q1和q2围成的三角形q1Oq2面积的2倍。因此，任意多边形的每个顶点可以表示为该平面原点处的向量（参见下页图5）并让最后一个顶点和第一个顶点重合（q1，q2，…qn，q1），按顺序将相邻向量叉乘，则多边形的面积为：

当 π＜θ≤ 2π 时，ε= 1 ；当 0 ≤θ≤ π 时，ε= 2。

dCB的计算原理类似，根据相似三角形和中点定理推导出三角形的形心公式（8）。计算图5中每个三角形qi Oqi+1的形心及面积矩，则整个多边形的形心为可由式（9）求出。

图5 求多边形的面积和形心

面元法和切片法有各自的优势：面元法需要先将船体模型划分成网格模型，但网格模型可以和后续的有限元结构分析以及水动力分析共享，减少后续工作量；切片法可以直接使用型值表或型线图的数据进行稳性计算，减少前期工作量。当切片数量或网格数量足够多时，两种方法都能准确地反映船体形状，因此都是可靠稳定的方法。

5 求解非线性方程组

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

f是一个连续非线性方程组，求解式（10）相当于求解连续函数l无约束下的极小值问题：

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

直到发现极值或者找不出新的xk为止。式中：下标k为迭代次数；pk是单位方向向量；τ是步长。

求无约束函数极值问题主要有两类方法：线搜索法和置信域法。线搜索法（比如梯度下降法，牛顿法等）基本思想是先找到搜索方向pk满足式（12），然后再求出合适的步长τ；而置信域法是先确定一个搜索半径Δ，然后在Δ内同时找到满足式（12）的pk和τ。如果找不到，则缩小Δ直到找到满足式（12）的pk和τ为止。

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

根据泰勒定理将目标函数l在xk附近近似为二次函数，如式（13）所示：

当Hk为正定矩阵时，式（14）的解如图6所示，可分为以下3种情况。

图6 狗腿算法求p

式中：pB为无约束条件下，式（14）的全局解。

式中：pU为无约束条件下，式（14）在-gk方向的解，

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

这样做的好处是在搜索半径内找不到式（14）的全局解时，尽可能利用二次导数信息Hk，在P U和P B的连线上找式（14）的近似解p，参见图6（c）。而当搜索半径相对较小时，式（13）可去除二次项，简化为线性近似，在梯度的负方向求式（14）的近似解p，参见图6（b）。

当Hk为非正定矩阵时，可以将Hk对角化，并对特征值取绝对值，将Hk转为正定矩阵。无法直接求出Hk时，可用有限差分法法求近似Hk。

搜索半径Δ的选取通常依据目标函数l和二次函数m的拟合比率r：

图7 置信域狗腿算法伪代码

6 收敛结果

为测试上述置信域狗腿法用于求解浮态方程的实际收敛效果，用MATLAB编写船舶稳性程序分别以牛顿法（N）和置信域狗腿法（T）计算两艘实船在不同工况下的浮态（水线面高度WLz，横倾角α和纵倾角β），并比较最后目标方程的冗余值‖f‖2和迭代次数。船舶的排水量和重心坐标为已知，同时给出迭代的初始值WLz0、α0和β0。

首先以图3和图4所示Heimdal小型工程船为例。该船主尺度为： 总长36 m、宽8.2 m、高3.2 m。

从下页表1中的计算结果可以看出：工况1的两种方法计算结果相同，船体处于接近正浮的状态，收敛的精度和迭代次数也基本相当；工况2的排水量增加，重心向船首右舷移动并升至2 m。若使用和工况1相同的初始值时，两种方法的计算结果都很理想（横倾23°，纵倾4°），但使用第3组初始值（0.3，0，0）时，牛顿法会收敛到另一个平衡位置（横倾-88°，纵倾7°）。虽然这个解也满足平衡方程，但与实际不符。究其原因是牛顿法并没有限制步长，在某一步迭代中跨过了最近的解，并收敛到另一个远处的解；置信域法由于步长限制在搜索半径内，因此计算结果正常。使用最后一组初始值时，牛顿法在迭代至第4步时失效，原因是该步算出的雅可比矩阵不可逆，因此无法算出这一步的解；置信域法由于保证在搜索半径内找到解，参见式（14）-式（17），因此计算结果正常。

表1 Heimdal小型工程船浮态计算结果

下面再以某大型集装箱船为例。该船主尺度为：总长288 m、宽32.2 m、高21.6 m。从表2所示计算结果可以看出：工况1时，初始值为（10，0，0），两种方法计算结果相同，但置信域法的迭代次数较多。这是由于当二次函数与目标函数拟合情况较好时，置信域狗腿算法会逐渐增大搜索半径（如图6增大1倍）；牛顿法由于没有限制步长，因此收敛更迅速，但同时也存在风险，如表1中的第3组例子。工况2时，重心向船尾右舷移动，初始值不变，但牛顿法再一次求解失败，原因和表1中的第4组例子一样，雅可比矩阵不可逆。

表2 某80 000 t集装箱船浮态计算结果

需要注意的是：本文采用的牛顿法没有限制步长是为保留其特点并突出和置信域法的区别，实际中也可以使用回溯法改进牛顿法达到控制步长的目的，但仍然存在雅可比矩阵可能是奇异矩阵的问题。

7 结 论

从实例的收敛结果来看，牛顿法较为依赖初始值的选择，因为雅可比矩阵可能产生不可逆的问题，导致求解中断；而相比于牛顿法，置信域法不仅更稳定、对初始值的依赖程度低，且能确保收敛，因此非常适合用于工业计算及稳性软件开发。

此外，求解浮态平衡方程的原理也可以推广到其他稳性相关的问题。比如计算复原力臂，只需将式（3）中的第3个方程去掉，已知横倾角α求解船体在Z轴和X轴方向的平衡方程，然后计算浮心和重心的Y轴坐标差即可，这种方法同时考虑了自由纵倾的影响。当舱室破损时，可以把进水看成附加质量，按照第4节的方法计算舱室进水的质量和重心，叠加到船体本身的质量、重心，再求解平衡方程即可。再比如计算海洋平台的变轴横倾稳定力臂，可以修改第2节的线性变换矩阵，使船体沿着指定方位角横倾即可。