一种多约束条件下的dp船推力优化分配仿真方法
技术领域
1.本发明涉及海洋技术领域，尤其涉及一种多约束条件下的dp船推力优化分配仿真方法。
背景技术：2.动力定位技术可以通过推进器保持位置和艏向于固定点或设定航迹上，并能够保证控制的精度，灵活有效，尤其是在锚泊受限的深海区，能够展现其强大的优势。对于动力定位船舶，推力分配模块是不可或缺的十分关键的模块。推力分配问题的本质是解决一个优化问题，如何从众多的推力组合方式中，选取最优的结果。
3.通常是选取某些技术指标构成优化问题的目标函数，根据实际物理系统的输出特性确定约束条件，然后通过相应的优化算法求解该问题获得最优解。推力分配方法的设计往往需要以目标船舶推进系统所配备推进器的特性为出发点，针对船舶所处的海洋环境、操作工况、推进器约束发生改变时分别设计相适应的推力分配算法。
4.在dp船推力优化分配仿真计算的过程中，以下问题；(1)常规的推力分配一般采用直接分配法和伪逆法，前者是通过求解多元方程来求出推力的大小和方向、后者是通过求逆矩阵来求出推力的大小和方向。虽然可以求解出推力，但是不能考虑或不能同时考虑各种约束条件；(2)由于两相近推进器之间存在桨-桨干扰问题，推进器的尾流会影响其相邻推进器产生的推力，因此会影响动力定位的效果。
技术实现要素：5.本发明要解决的技术问题是如何提供一种多约束条件下的dp船推力优化分配仿真方法，实现有推进约束的情况下的推力优化分配，保证系统在不同海洋环境下考虑推进器输出推力大小和回转可行区域的情况下能准确地输出推力，完成不同情况下的作业任务，为后续动力定位系统模拟器的设计奠定基础。
6.为了解决上述技术问题，本发明提供了一种多约束条件下的dp船推力优化分配仿真方法，通过以下步骤完成：
7.(1)建立动力定位船舶的运动坐标系，包括北东坐标系和随船坐标系；
8.(2)建立船舶的运动状态在两坐标系的转换关系；
9.(3)建立船舶刚体动力学模型；
10.(4)建立船舶水动力学模型；
11.(5)根据动力定位船推进器的布置对其推进系统建立数学模型；
12.(6)确认各推进器的安装位置；
13.(7)对该船的合推力进行建模；
14.(8)建立推进器的干扰模型，确定推力区间和回转禁止域；
15.(9)建立推力优化分配的数学模型；
16.(10)利用起作用集法求解得出推力分配结果，完成对船舶运动的一次控制仿真。
17.进一步的，在步骤2中，船舶运动状态下两坐标的转换方法为：
[0018][0019][0020]
式中：为船舶在北东坐标系下的速度向量；
[0021]
j(ψ)为两坐标系之间的旋转矩阵。
[0022][0023]
υ＝[u,v,r]
t
为船舶在随体坐标系下的速度向量；
[0024]
η＝[n,e,ψ]
t
是船舶在北东坐标系下的位置及姿态向量。
[0025]
进一步的，在步骤(3)中，在建立船舶刚体运动学时设船舶重心坐标为(xg，0)，根据牛顿-欧拉方程，可以得到船舶在随体坐标系下的刚体动力学方程：
[0026][0027]
式中：m
rb
代表船体惯性矩阵：
[0028][0029]crb
代表船体的科里奥利—向心力矩阵：
[0030][0031]
m代表船舶刚体的质量；
[0032]iz
代表船体绕obzb轴的转动惯量；
[0033]
υ＝[u,v,r]
t
代表船舶控制点处的速度向量；
[0034]
τ
rb
＝[x,y,n]
t
代表作用在控制点处的合外力，x代表纵向力，y代表横向力，n代表转艏力矩。
[0035]
进一步的，在步骤(4)中，建立船舶水动力模型时，做了出三个假设：(1)研究的船舶具有一般性；(2)将船舶视为刚体；(3)船舶运动区域为理想状态的零流速水域；
[0036]
基于上述三个假设，建立如下关系式：
[0037][0038]
式中：ma为附加质量矩阵：
[0039]
[0040]dl
为线性水阻尼系数矩阵：
[0041][0042]dn
(υ)为非线性水阻尼矩阵：
[0043][0044]
τh＝[xhyhnh]
t
：为作用在船中心点上的三自由度水动力。
[0045]
因此综合船舶的刚体动力学和水动力学可以得到船舶综合动力学方程：
[0046][0047]
式中：m＝m
rb
+ma：为系统惯量矩阵；
[0048]
c(υ)＝c
rb
(υ)+ca(υ)：为船舶科里奥利-向心力矩阵；
[0049]
d(υ)＝d
l
+dn(υ)：为船舶的水阻尼矩阵；
[0050]
τ：为推进系统对船舶的作用合力。
[0051]
进一步的，在步骤(8)中，推进器的桨-桨干扰模型是通过建立起螺旋桨的三维模型，并在cfd中进行桨-桨干扰的仿真计算，从而得到桨-桨干扰的作用域，进而得到每个螺旋桨的回转可行域；
[0052]
推力区间的设置根据推力减额因数的大小来确定，推力减额因数的计算的公式为：
[0053][0054]
式中：t
φ
为推力减额因数，推力减额因数的定义是用后方螺旋桨产生的推力与前方螺旋桨产生的推力之间的比值，其中t是φ＝0时的推力减额因数。
[0055]
进一步的，在步骤(9)中，推力优化分配的数学模型为
[0056][0057]
s.t.:s＝τ-b(α)f(a)
[0058]fmin
≤f≤f
max
(b)
[0059]
δf
min
≤f-f0≤δf
max
(c)
[0060]
α
min
≤α≤α
max
(d)
[0061]
δα
min
≤α-α0≤δα
max
(e)
[0062]-∞≤s≤∞(f)
[0063]
目标函数中第一项为降低功率消耗，第二项为降低分配误差，第三项为控制全回转推进器回转速度，最后一项为避免全回转推进器进入奇点位置；
[0064]
在约束条件中，条件(a)为等式约束，条件(b)为推力极值限制，条件(c)为推力响应限制，条件(d)为全回转推进器回转可行域，条件(e)为全回转推进器回转速度限制；
[0065]
式中：τ＝[x
t y
t n
t
]
t
为期望推力；f为推进器的推力；s为推力分配的误差；α为全回转推进器回转角；w、q、ω为权值矩阵；ρ、ε为常数；b(α)为推进器布置矩阵。
[0066]
本发明的技术效果在于：
[0067]
(1)本发明根据动力定位船推进器需要考虑的推进约束条件，建立了推进器的干扰模型，并对推进器回转可行域进行了分析和计算，得到了不影响实际输出推力的推进器回转可行域；
[0068]
(2)对求解方式进行了改进，采用起作用集法对推力优化分配问题进行求解，最终结果充分考虑了功率最优、降低分配误差、回转速度及推进器响应能力等推力优化分配的约束条件；
[0069]
(3)由于起作用集法是将上一时刻的最优解作为下次计算的初始值，这样可以降低寻找最优解的迭代次数，提高计算速度，提高了系统的响应速度，因此推力分配效果优于其他方法。
附图说明
[0070]
图1北东坐标系和随船坐标系。
[0071]
图2某动力定位船推进器的布置图。
[0072]
图3各推进器的合力(矩)。
[0073]
图4全回转推进器间干扰模型。
[0074]
图5全回转推进器推力可行域。
具体实施方式
[0075]
下面结合附图和具体实施例对本发明作进一步说明，以使本领域的技术人员可以更好的理解本发明并能予以实施，但所举实施例不作为对本发明的限定。
[0076]
本发明实施例提供了一种多约束条件下的dp船推力优化分配仿真方法，结合附图1-5所示，每一步具体实施如下：
[0077]
(1)建立某动力定位船舶的运动坐标系：北东坐标系和随船坐标系，如图1所示，北东坐标系为在地球表面局部区域的船舶运动依赖的坐标系，而随船坐标系则主要用来描述船舶运动速度、角速度等物理量，如图1所示。
[0078]
(2)将船舶的运动状态在两坐标系下进行转换，方法如下：
[0079][0080][0081]
式中：为船舶在北东坐标系下的速度向量；
[0082]
j(ψ)为两坐标系之间的旋转矩阵。
[0083][0084]
υ＝[u,v,r]
t
为船舶在随体坐标系下的速度向量；
[0085]
η＝[n,e,ψ]
t
是船舶在北东坐标系下的位置及姿态向量；
[0086]
(3)建立船舶刚体动力学模型
[0087]
设船舶重心坐标为(xg，0)，根据牛顿-欧拉方程，可以得到船舶在随体坐标系下的刚体动力学方程：
[0088][0089]
式中：m
rb
代表船体惯性矩阵：
[0090][0091]crb
代表船体的科里奥利—向心力矩阵：
[0092][0093]
m代表船舶刚体的质量；
[0094]iz
代表船体绕obzb轴的转动惯量；
[0095]
υ＝[u,v,r]
t
代表船舶控制点处的速度向量；
[0096]
τ
rb
＝[x,y,n]
t
代表作用在控制点处的合外力，x代表纵向力，y代表横向力，n代表转艏力矩。
[0097]
(4)建立船舶水动力学模型
[0098]
考虑船舶的动力学问题时要用到牛顿定律，所以必须在惯性参考系下进行，船舶的运动数学模型是通过cfd数值模拟进行计算，得到水动力学参数，以保证数学模型的精确性和可信度。
[0099]
根据动力定位系统研究理论，船舶的水动力学模型主要考虑附加质量惯性力、潜在阻尼力以及由于重力和浮力产生的阿基米德回复力这三部分，通常是通过试验测量得到的相关数据，且试验时选取船舶的中心点作为控制点。在海洋中航行的船舶，其主要受力就是水动力，所以，无论是研究动力定位系统还是其他相关的船舶运动理论，研究船舶水动力都是前提。
[0100]
由于水动力的形成机制非常复杂，不仅包含线性因子同时也包含非线性因子，因此，本发明提出以下三种建模假设：研究的船舶具有一般性；将船舶视为刚体；船舶运动区域为理想状态的零流速水域。
[0101][0102]
式中：ma为附加质量矩阵：
[0103][0104]dl
为线性水阻尼系数矩阵：
[0105][0106]dn
(υ)为非线性水阻尼矩阵：
[0107][0108]
τh＝[x
h y
h nh]
t
：为作用在船中心点上的三自由度水动力。
[0109]
因此综合船舶的刚体动力学和水动力学可以得到船舶综合动力学方程：
[0110][0111]
式中：m＝m
rb
+ma：为系统惯量矩阵；
[0112]
c(υ)＝c
rb
(υ)+ca(υ)：为船舶科里奥利-向心力矩阵；
[0113]
d(υ)＝d
l
+dn(υ)：为船舶的水阻尼矩阵；
[0114]
τ：为推进系统对船舶的作用合力。
[0115]
(5)根据动力定位船推进器的布置对其推进系统建立数学模型
[0116]
根据某动力定位船推进器的布置对其推进系统建立数学模型，推进系统数学模型是通过机理分析法，通过船型数据给出的推进器的坐标，建立起各推进器的相对位置的数学模型。由于全回转推进器输出的推力是矢量，需要将各个推力分解到沿着船体坐标系的两坐标轴方向上，将同一方向上的所有分力代数求和、将各分力与其对应的力臂距离做积再对各积求和，便得到推进系统作用于船体上的合力和合力矩，图2为其推进器布置图。
[0117]
(6)经查阅相关资料，获取到各推进器的安装位置如下表，设1～7号推进器的在随船坐标系下的坐标分别为(x1,y1)，(x2,y2)，(x3,y3)，(x4,y4)，(x5,y5)，(x6,y6)，(x7,y7)，具体如表1所示
[0118]
表1各推进器安装位置
[0119]
[0120]
(7)对该船的合推力建模，船体坐标系下对船产生的合力(矩)如图3所示；
[0121]
则7台推进器在船体坐标系中心点处产生的，纵向推力x
t
，横向推力y
t
，转艏力矩n
t
公式如下：
[0122][0123]
将公式(12)转化为矢量的形式为：
[0124]
τ＝b(α)t(13)
[0125]
式中：
[0126]
τ＝[x
t
,y
t
,n
t
]
t
为作用船体上旋转中心上的三自由度合力；
[0127]
t＝[t1,t2,t3,t4,t5,t6,t7]
t
为推进器输出的推力矢量；
[0128]
b(α)为推进器布置矩阵：
[0129][0130]
式中：c代表cos()，s代表sin()。
[0131]
(8)确定推力区间和回转禁止域
[0132]
推力分配模型需要对推进器输出推力大小和回转方位角有一定的限制。推力大小是根据推进器的推进能力限定，而方位角的限制主要是针对于全回转推进器，避免全回转推进器排出的尾流影响其后面的推进器的推力产生，使被干扰的推进器在原有指令的控制下输出的推力下降，影响推进器实际输出推力。
[0133]
许多拖曳水池都针对螺旋桨之间的干扰问题展开过试验研究，虽然使用的螺旋桨类型不尽相同，但结果一致性很高。结果表明，螺旋桨间干扰造成的推力损失最高可达70％，两螺旋桨间距越小，推力损失越严重。主要原因是上游螺旋桨的尾流，改变了下游螺旋桨的入流速度。由螺旋桨进速比公式
[0134][0135]
可知当螺旋桨转速n及直径d一定时，进速va增大会导致进速系数j增大。由螺旋桨敞水特性、推力系数及转矩系数公式
[0136][0137][0138]
可知，进速系数j变大，推力系数k
t
和转矩系数kq减小。当螺旋桨转速一定时，推力和转矩都会随之下降；总之，上游螺旋桨的尾流速度以及方位角直接影响下游螺旋桨的水动力性能。此外，下游螺旋桨的运转也会对上游螺旋桨的水动力性能产生影响。如下游螺旋
桨加快上游螺旋桨尾流速度，使得上游螺旋桨尾流中心压力下降，形成类似“形状阻力”的水动力损失。当然，推力损失与否要视具体情况而定,不同的工作状态下，结果可能完全不同，甚至会出现改善水动力性能的情况。
[0139]
推进器的桨-桨干扰模型是通过建立起螺旋桨的三维模型，并在cfd中进行桨-桨干扰的仿真计算，从而得到桨-桨干扰的作用域，进而得到每个螺旋桨的回转可行域。通常螺旋桨排出的尾流以一定的角度从推进器向外扩展，一般为8～10度。为避免推进器间的干扰，需要为全回转推进器设置推力禁区。推力禁区的设置根据推力减额因数的大小来确定。推力减额因数的计算的公式：
[0140][0141]
式中：t
φ
为推力减额因数，推力减额因数的定义是用后方螺旋桨产生的推力与前方螺旋桨产生的推力之间的比值。其中t是φ＝0时的推力减额因数。全回转推进器之间的干扰模型如图4所示。
[0142]
结合实际情况，该动力定位船7个推进器的推力限制和回转禁区的设置如下：
[0143]
表2推进器推力限制和回转禁区
[0144][0145]
根据上表，容易得出7个全回转推进器的推力可行域，如图5所示
[0146]
(9)建立推力优化分配的数学模型
[0147][0148]
该推力优化分配数学模型考虑到了功率最优，降低分配误差，考虑了推进器的回
转禁止域，回转速度以及推进器推力响应能力，以及避免全回转推进器之间的奇异性。目标函数中第一项为降低功率消耗，第二项为降低分配误差，第三项为控制全回转推进器回转速度，最后一项为避免全回转推进器进入奇点位置；在约束条件中，条件(a)为等式约束，条件(b)为推力极值限制，条件(c)为推力响应限制，条件(d)为全回转推进器回转可行域，条件(e)为全回转推进器回转速度限制。
[0149]
式中：τ＝[x
t y
t n
t
]
t
为期望推力；f为推进器的推力；s为推力分配的误差；α为全回转推进器回转角；w、q、ω为权值矩阵；ρ、ε为常数；b(α)为推进器布置矩阵。
[0150]
(10)利用起作用集法，对于一般正定二次规划式(17)，先求出问题(17)式的一个可行点x
(k)
，计算点x
(k)
处有起作用约束指标集jk。
[0151][0152]
求解相应的等式约束的二次规划问题(18)式。
[0153][0154]
设其最优解为乘子向量为λk。令即认为pk是从点x
(k)
出发至的方向，如求出了pk也就找到了也就找到了因此求解(18)式可以化成求解：
[0155][0156]
设的起作用约束指标集为j
k+1
，则根据与x
(k)
之间不同关系来调整jk，使目标函数值不断减少。按照这种思路继续，得到j
*
，从而得到(17)式的最优解x
*
。
[0157]
至此，起作用集法求解结束，得出推力分配结果，完成对船舶运动的一次控制。本发明中建立的船舶运动数学模型和得到的推力优化分配结果，与实船动力定位能力分析报告进行对比，仿真结果与实船动力定位能力分析报告基本一致，具有较高的精度和可信度。
[0158]
采用起作用集法将上一时刻的最优解作为下次计算的初始值，以此降低寻找最优解的迭代次数，提高了计算速度和系统的响应速度，并成功应用于海洋船舶的动力定位仿真系统的推力分配实时求解中，并有效避免了推力饱和和回转禁止域。
[0159]
由于使用直接分配法和伪逆法解决推力分配问题不能考虑到同时存在的多种约束条件，而起作用集法正适合于有不等式约束的二次规划问题的求解，因此在求解多约束条件下的推力优化分配仿真问题是恰到好处的。
[0160]
以上所述实施例仅是为充分说明本发明而所举的较佳的实施例，本发明的保护范围不限于此。本技术领域的技术人员在本发明基础上所作的等同替代或变换，均在本发明的保护范围之内。本发明的保护范围以权利要求书为准。
技术特征：1.一种多约束条件下的dp船推力优化分配仿真方法，通过以下步骤完成，其特征在于：(1)建立动力定位船舶的运动坐标系，包括北东坐标系和随船坐标系；(2)建立船舶的运动状态在两坐标系的转换关系；(3)建立船舶刚体动力学模型；(4)建立船舶水动力学模型；(5)根据动力定位船推进器的布置对其推进系统建立数学模型；(6)确认各推进器的安装位置；(7)对该船的合推力进行建模；(8)建立推进器的干扰模型，确定推力区间和回转禁止域；(9)建立推力优化分配的数学模型；(10)利用起作用集法求解得出推力分配结果，完成对船舶运动的一次控制仿真。2.根据权利要求1所述的多约束条件下的dp船推力优化分配仿真方法，其特征在于：在步骤2中，船舶运动状态下两坐标的转换方法为：船舶运动状态下两坐标的转换方法为：式中：为船舶在北东坐标系下的速度向量；j(ψ)为两坐标系之间的旋转矩阵。υ＝[u,v,r]
t
为船舶在随体坐标系下的速度向量；η＝[n,e,ψ]
t
是船舶在北东坐标系下的位置及姿态向量。3.根据权利要求1所述的多约束条件下的dp船推力优化分配仿真方法，其特征在于：在步骤(3)中，在建立船舶刚体运动学时设船舶重心坐标为(x
g
，0)，根据牛顿-欧拉方程，可以得到船舶在随体坐标系下的刚体动力学方程：式中：m
rb
代表船体惯性矩阵：c
rb
代表船体的科里奥利—向心力矩阵：m代表船舶刚体的质量；
i
z
代表船体绕o
b
z
b
轴的转动惯量；υ＝[u,v,r]
t
代表船舶控制点处的速度向量；τ
rb
＝[x,y,n]
t
代表作用在控制点处的合外力，x代表纵向力，y代表横向力，n代表转艏力矩。4.根据权利要求1所述的多约束条件下的dp船推力优化分配仿真方法，其特征在于：在步骤(4)中，建立船舶水动力模型时，做了出三个假设：(1)研究的船舶具有一般性；(2)将船舶视为刚体；(3)船舶运动区域为理想状态的零流速水域；基于上述三个假设，建立如下关系式：式中：m
a
为附加质量矩阵：d
l
为线性水阻尼系数矩阵：d
n
(υ)为非线性水阻尼矩阵：τ
h
＝[x
h y
h n
h
]
t
：为作用在船中心点上的三自由度水动力。因此综合船舶的刚体动力学和水动力学可以得到船舶综合动力学方程：式中：m＝m
rb
+m
a
：为系统惯量矩阵；c(υ)＝c
rb
(υ)+c
a
(υ)：为船舶科里奥利-向心力矩阵；d(υ)＝d
l
+d
n
(υ)：为船舶的水阻尼矩阵；τ：为推进系统对船舶的作用合力。5.根据权利要求1所述的多约束条件下的dp船推力优化分配仿真方法，其特征在于：在步骤(8)中，推进器的桨-桨干扰模型是通过建立起螺旋桨的三维模型，并在cfd中进行桨-桨干扰的仿真计算，从而得到桨-桨干扰的作用域，进而得到每个螺旋桨的回转可行域；推力区间的设置根据推力减额因数的大小来确定，推力减额因数的计算的公式为：
式中：t
φ
为推力减额因数，推力减额因数的定义是用后方螺旋桨产生的推力与前方螺旋桨产生的推力之间的比值，其中t是φ＝0时的推力减额因数。6.根据权利要求1所述的多约束条件下的dp船推力优化分配仿真方法，其特征在于：在步骤(9)中，推力优化分配的数学模型为s.t.:s＝τ-b(α)f(a)f
min
≤f≤f
max
(b)δf
min
≤f-f0≤δf
max
(c)α
min
≤α≤α
max
(d)δα
min
≤α-α0≤δα
max
(e)-∞≤s≤∞(f)目标函数中第一项为降低功率消耗，第二项为降低分配误差，第三项为控制全回转推进器回转速度，最后一项为避免全回转推进器进入奇点位置；在约束条件中，条件(a)为等式约束，条件(b)为推力极值限制，条件(c)为推力响应限制，条件(d)为全回转推进器回转可行域，条件(e)为全回转推进器回转速度限制；式中：τ＝[x
t
y
t
n
t
]
t
为期望推力；f为推进器的推力；s为推力分配的误差；α为全回转推进器回转角；w、q、ω为权值矩阵；ρ、ε为常数；b(α)为推进器布置矩阵。
技术总结本发明公开了一种多约束条件下的DP船推力优化分配仿真方法，建立了船舶的运动数学模型、推进器间干扰模型以及推进系统数学模型，在能保持船舶目标运动或姿态的同时，综合考虑推进器总功率、推进器的磨损和推进器的回转可行区域等优化条件，实现有推进约束的情况下的推力优化分配，保证系统在不同海洋环境下考虑推进器输出推力大小和回转可行区域的情况下能准确地输出推力，完成不同情况下的作业任务，为后续动力定位系统模拟器的设计奠定基础。础。础。
技术研发人员：刘培林 田立锋 王屹 赵福臣 李新飞 李斌 喻龙 张广磊 昝英飞
受保护的技术使用者：哈尔滨工程大学
技术研发日：2022.03.18
技术公布日：2022/7/5
