CN103092077B - 动力定位系统的推力分配方法 - Google Patents

Abstract

一种动力定位系统的推力分配方法，采用以下步骤：一：将推进器分为数个小组，并固定组内推进器方位角的相对角度，但不固定绝对角度；二：由控制指令的低频部分使用优化算法求出一个缓慢变化的推进器方位角，再使用优化算法求出推进器的推力；具体步骤为：将k时刻的控制指令τ_c，k输入到由一阶惯性环节构成的低通滤波器中，然后，再进一步离散化得低频部分τ_k+1；由控制指令τ_c，k采用序列二次规划法计算出各推进器的推力大小f_k和推进器方位角α_k；再应用二次规划法求出推进器的推力f_k；并在所设置的目标函数和部分约束条件的公式中，加入推力折减系数p。本发明可以将推进器的推力和推进器的方位角分开计算，避免了由于线性化造成的误差及奇异结构，降低了推进器的能量消耗。

Description

动力定位系统的推力分配方法

技术领域

本发明涉及动力定位系统，尤其涉及一种动力定位系统的推力分配方法。

背景技术

动力定位系统(DynamicPositioningSystem)是利用自身的推进器系统自动保持船舶的位置和艏向的计算机控制的闭环系统，主要由测量系统、控制系统和推进器系统组成。其功能是：不借助于锚泊系统，而能自动检测出船舶的实际位置与目标位置的偏差，结合风、浪、流等外界扰动力的影响计算出能使船舶恢复到目标位置所需推力和扭矩的大小，并对各推进器进行推力分配，进而使各推进器产生相应的推力和扭矩，使船舶尽可能地保持在海平面上要求的位置上。其中，推进器系统是动力定位系统的一个重要组成部分，作为动力定位系统的执行机构，推进器用于产生力和力矩，用以抗衡作用于船上的干扰力和干扰力矩。推进器一般来说是指螺旋桨，在动力定位系统中除一般的主推外，应用比较多的是槽道螺旋桨、全方位螺旋桨和吊舱螺旋桨。

如图1所示，为深水铺管起重船的推进器布置，由于推进器的数量很多，已经达到了7个，且均为全回转推进器，所以，其存在着无数多个不同的推力和方向的组合，且均满足特定的水平力和船艏摇力矩。这时，需要建立推力分配算法，将推力分配到各个推进器上，不仅使各个推进器产生的推力满足推力指令的要求，还要使推进器消耗的能量最小，以防止推进器过度磨损等。

设τ∈R³为控制器输出的推力指令，包括所需要的纵荡方向的力、横荡方向的力和船艏摇力矩。f∈Rⁿ为各推进器发出的推力，n为推进器数量。则需要满足τ＝B(α)f，其中

$B\left(\mathrm{\α}\right)=\left(\begin{array}{cccc}{\mathrm{cos\α}}_1& {\mathrm{cos\α}}_2& \mathrm{...}& {\mathrm{cos\α}}_n\\ {\mathrm{sin\α}}_1& {\mathrm{sin\α}}_2& \mathrm{...}& {\mathrm{sin\α}}_n\\ -y_1{\mathrm{cos\α}}_1+x_1{\mathrm{sin\α}}_1& -y_2{\mathrm{cos\α}}_2+x_2{\mathrm{sin\α}}_2& \mathrm{...}& -y_n{\mathrm{cos\α}}_n+x_n{\mathrm{sin\α}}_n\end{array}\right)$

B(α)为一关于推进器方位角和推进器坐标位置的矩阵，α＝[α₁，α₂，...，α_n]，α_i为第i个推进器的方位角，(x_i，y_i)为第i个推进器的坐标。

推进器分组法就是将几个推进器分成相应的组，每个组负责发出特定方向的推力。这种方法的优点是计算量小、简单可靠。

上述方法中的目标函数与约束条件为：

$\begin{array}{cc}\min & J_{QP}(\mathrm{\Δ}\mathrm{\α},\mathrm{\Δ}f,s)=\underset{i=1}{\overset{m}{\Σ}}\left(\frac{{\mathrm{dW}}_i}{{\mathrm{df}}_i}\right(f_{0,i}){\mathrm{\Δf}}_i+\frac{d^2{W}_i}{{{\mathrm{df}}_i}^2}(f_{0j}\left){\mathrm{\Δf}}_i^2\right)+s^{T}Qs\end{array}$

s.t.s+B(α₀)Δf＝τ-B(α₀)f₀

f_min-f₀≤Δf≤f_max-f₀

Δf_min≤Δf≤Δf_max

其中，W_i为：第i个推进器功率；f₀为：上一采样时刻计算出的推进器的推力大小；Δf为：推进器推力变化量，f_min、f_max分别为：对应推进器最小推力、最大推力；s为：松弛变量。这一采样时刻所计算出的各推进器的推力大小就为：f＝f₀+Δf。虽然，这种方法考虑了对推进器最大推力的限制，推力最大变化速度的限制，但是，由于各推进器的方位角都是固定的，所以，不能充分利用推进器的性能，不能使推进器消耗的能量最小。

而采用优化算法的推力分配方法的优点是可以最小化推进系统的能耗，包括：燃油能耗、推进器的磨损、推力误差等，而且，还考虑到推进器的推力极限、最大推力变化速率、方位角最大旋转速率、推进器推力的禁区、奇异结构等因素。

上述方法中的目标函数与约束条件的公式为：

$\begin{array}{c}\begin{array}{cc}\min & J_{QP}(\mathrm{\Δ}\mathrm{\α},\mathrm{\Δ}f,s)=\underset{i=1}{\overset{m}{\Σ}}\left(\frac{{\mathrm{dW}}_i}{{\mathrm{df}}_i}\right(f_{0,i}){\mathrm{\Δf}}_i+\frac{d^2{W}_i}{{{\mathrm{df}}_i}^2}(f_{0,i}\left){{\mathrm{\Δf}}_i}^2\right)+s^{T}Qs+{\mathrm{\Δ\α}}^T\mathrm{\Ω}\mathrm{\Δ}\mathrm{\α}\end{array}\\ +\frac{d}{d\mathrm{\α}}{\left(\frac{\mathrm{\δ}}{\mathrm{\ϵ}+\det \left(B\right(\mathrm{\α}\left)B^T\right(\mathrm{\α}\left)\right)}\right)}_{\mathrm{\α}={\mathrm{\α}}_0}\mathrm{\Δ}\mathrm{\α}\end{array}$

$\begin{array}{cc}s.t.& s+B\left({\mathrm{\α}}_0\right)\mathrm{\Δ}f+\frac{\∂}{\∂\mathrm{\α}}\left(B\right(\mathrm{\α}\left)f\right){|}_{\underset{f=f_0}{\mathrm{\α}={\mathrm{\α}}_0}}\·\mathrm{\Δ}\mathrm{\α}=\mathrm{\τ}-B\left({\mathrm{\α}}_0\right)f_0\end{array}$

f_min-f₀≤Δf≤f_max-f₀

Δf_min≤Δf≤Δf_max

α_min-α₀≤Δα≤α_max-α₀

Δα_min≤Δα≤Δα_max

其中，W_i为：第i个推进器功率，f₀为：上一采样时刻计算出的各推进器的推力大小；α₀为：上一采样时刻计算出的各推力器的方位角；Δf、Δα分别为：推进器推力变化量、方位角变化量；f_min、f_max分别为：对应推进器最小推力、最大推力；α_min、α_max分别为推进器最小方位角、最大方位角，s为松弛变量；δ、ε分别为常量。这一采样时刻所计算出的各推进器的推力大小就为f＝f₀+Δf，而各推进器的方位角就为：α₁＝α₀+Δα。虽然，这种方法考虑了对推进器最大推力的限制，推力最大变化速度的限制，方位角的范围，方位角变化的最大速度的限制，并且也避免了奇异结构的出现。这种方法的另外一大优点就是基本不用考虑推进器的推力可行域为非凸集的情况。但是，上述通过采用序列二次规划方法求解上述约束条件下的最优化问题，由于加入了奇异值避免项

因而，导致计算量较大，而且，在算法中往往会做线性化处理，即： $s+B\left({\mathrm{\α}}_0\right)\mathrm{\Δ}f+\frac{\∂}{\∂\mathrm{\α}}\left(B\right(\mathrm{\α}\left)f\right){|}_{\underset{f=f_0}{\mathrm{\α}={\mathrm{\α}}_0}}\·\mathrm{\Δ}\mathrm{\α}=\mathrm{\τ}-B\left({\mathrm{\α}}_0\right)f_0,$ 这也就引入了误差。

由于上述推力分配算法都存在着一定的缺陷，影响了它们的应用，所以需要对这些方法加以改进，以创建一种适合动力定位系统推力分配的算法。

发明内容

本发明的主要目的在于克服现有技术存在的上述缺点，而提供一种动力定位系统的推力分配方法，其可以将推进器的推力和推进器的方位角分开计算，避免了由于线性化造成的误差；而且，考虑了对推进器最大推力、推力最大变化速度的限制、方位角的范围及方位角变化的最大速度的限制；同时，也避免了奇异结构的出现；本发明另一目的是基本不用考虑推进器的推力可行域为非凸集的情况，而且，由于加入了推力折减系数p，不仅计算量小，而且，更容易使推进器的方位角转到一个合适的位置；同时，避免了奇异结构，降低了推进器的能量消耗。

本发明的目的是由以下技术方案实现的：

一种动力定位系统的推力分配方法，其特征在于：采用以下步骤：

第一步：在进行推力分配之前，先将推进器分为数个小组，并固定组内推进器方位角的相对角度，但不固定它们的绝对角度；

第二步：由控制指令的低频部分使用优化算法求出一个缓慢变化的推进器方位角，在已知方位角后，再使用优化算法求出推进器的推力；

其具体步骤如下：

(1)首先，将k时刻的控制指令τ_c，k输入到一个由一阶惯性环节构成的低通滤波器中，然后，再进一步离散化可得缓慢变化的控制指令的低频部分τ_k+1；

(2)由控制指令τ_c，k采用序列二次规划法计算出各推进器的推力大小f_k和推进器方位角α_k；

(3)得到推进器方位角α_k后，再应用二次规划法求出推进器的推力f_k；其中，在计算推进器方位角α_k和推力f_k时，在所设置的目标函数和部分约束条件的公式中，加入推力折减系数p。

所述低通滤波器的形式为：其中，为k时刻控制指令的一阶导数；τ_k为：低通滤波后的k时刻的控制指令，Λ为：一阶惯性环节低通滤波矩阵，且元素均为正的对角矩阵，τ_c，k为第k个时刻的控制指令；τ_k+1为：在第k个时刻预测的第k+1个时刻的控制指令的低频部分；为了便于实时计算，避免求解微分方程，利用基本的求导公式并令其中，h为：采样时间，τ_k+1为：低通滤波后k+1时刻的控制指令，τ_k为低通滤波后k时刻的控制指令；而则与 等价。

所述推进器推力f_k和方位角α_k的计算是分开进行的。

所述目标函数和部分约束条件的具体公式为：

$\begin{array}{cc}\min & J_{QP}(\mathrm{\Δ}\mathrm{\α},\mathrm{\Δ}f,s)=\underset{i=1}{\overset{m}{\Σ}}\left(\frac{{\mathrm{dW}}_i}{{\mathrm{df}}_i}\right(f_{0,i}){\mathrm{\Δf}}_i+\frac{d^2{W}_i}{{\mathrm{df}}_i^2}(f_{0,i}\left){\mathrm{\Δf}}_i^2\right)+s^{T}Qs+{\mathrm{\Δ\α}}^T\mathrm{\Ω}\mathrm{\Δ}\mathrm{\α}---\left(1\right)\end{array}$

s.t.f_min/p-f₀≤Δf≤f_max/p-f₀

min{J＝f^TWf+s^TQs}(2)

s.t.Bf＝τ+s

上述公式：f_min/p-f₀≤Δf≤f_max/p-f₀中，f₀为：前一时刻推进器推力，Δf为当前时刻与前一时刻力的变化量，f_min、f_max分别对应推进器推力最小值、最大值；J_QP(Δα，Δf，s)为：关于Δα、Δf和s的目标函数；

Δα为：m维当前时刻与前一时刻的方位角变化向量；

Δf为：m维当前时刻与前一时刻的推力变化向量；

s为：3维松弛变量；

W_i为：关于第i个推进器功率；

f_i为：第i个推进器当前时刻推力；

f_0，i为第i个推进器前一时刻推力；

Δf_i为：第i个推进器当前时刻与前一时刻推力变化量；

Q为：3阶对角矩阵；

Ω为：m阶对角矩阵；

P为：推力折减系数；

J为：推力分配目标函数；

τ为：控制器输出的推力指令；

W：为推进器个数；

$B=\left(\begin{array}{cccc}{\mathrm{cos\α}}_1& {\mathrm{cos\α}}_2& ...& {\mathrm{cos\α}}_m\\ {\mathrm{sin\α}}_1& {\mathrm{sin\α}}_2& ...& {\mathrm{sin\α}}_m\\ -y_1{\mathrm{cos\α}}_1+x_1{\mathrm{sin\α}}_1& -y_2{\mathrm{cos\α}}_2+x_2{\mathrm{sin\α}}_2& ...& -y_m{\mathrm{cos\α}}_m+x_m{\mathrm{sin\α}}_m\end{array}\right)$

α_i为第i个推进器的方位角，(x_i，y_i)为第i个推进器的坐标，坐标系原点为船舶重心。

本发明的有益效果：

1.将推进器的推力和推进器的方位角分开计算，避免了由于线性化造成的误差；

2.考虑了对推进器最大推力、推力最大变化速度的限制、方位角的范围及方位角变化的最大速度的限制；

3.避免了奇异结构的出现；

4.基本不用考虑推进器的推力可行域为非凸集的情况，而且，由于加入了推力折减系数p，不仅计算量小，而且，更容易使推进器的方位角转到一个合适的位置；

5.避免了奇异结构，降低了推进器的能量消耗。

附图说明：

图1为现有推进器布置示意图。

图2为本发明推力分配算法流程图。

图中主要标号说明：

1.推进器、2.推进器、3.推进器、4.推进器、5.推进器、6.推进器、7.推进器。

具体实施方式

如图1，2所示，本发明采用步骤如下：

第一步：在进行推力分配之前，先将推进器1-7进行分组，即：推进器1和2为一组，推进器3和4为一组，推进器5和6为一组，推进器7为一组，只固定组内推进器方位角的相对角度，而不固定它们的绝对角度；

第二步：由控制指令的低频部分使用优化算法求出一个缓慢变化的推进器方位角，在已知方位角后，再使用优化算法求出推进器的推力；

其具体步骤如下：

(1)首先，将k时刻的控制指令τ_c，k输入到一个由一阶惯性环节构成的低

通滤波器中，低通滤波器的形式为：其中，为k时刻控制指令的一阶导数；τ_k为：低通滤波后的k时刻的控制指令，Λ为：一阶惯性环节低通滤波矩阵，且元素均为正的对角矩阵，τ_c，k为第k个时刻的控制指令，然后，再进一步离散化可得缓慢变化的控制指令的低频部分τ_k+1；τ_k+1为：在第k个时刻预测的第k+1个时刻的控制指令的低频部分。为了便于实时计算，避免求解微分方程，利用基本的求导公式并令其中，h为：采样时间，τ_k+1为：低通滤波后k+1时刻的控制指令，τ_k为低通滤波后k时刻的控制指令；而则与 等价。

(2)由控制指令τ_c，k采用序列二次规划法(二次规划法为常用的求解带有非线性约束的最优化方法)，计算出各推进器的推力大小f_k和推进器方位角α_k；其中，α_k为第k个时刻的各推进器的方位角，f_k为第k个时刻的各推进器发出的推力。

(3)得到推进器方位角α_k后，再应用二次规划法求出推进器的推力f_k。

其中，在计算推进器方位角α_k和推力f_k时，在其所设置的目标函数和部分约束条件中，加入避免奇异结构的推力折减系数p，其具体公式为：

$\begin{array}{cc}\min & J_{QP}(\mathrm{\Δ}\mathrm{\α},\mathrm{\Δ}f,s)=\underset{i=1}{\overset{m}{\Σ}}\left(\frac{{\mathrm{dW}}_i}{{\mathrm{df}}_i}\right(f_{0,i}){\mathrm{\Δf}}_i+\frac{d^2{W}_i}{{\mathrm{df}}_i^2}(f_{0,i}\left){\mathrm{\Δf}}_i^2\right)+s^{T}Qs+{\mathrm{\Δ\α}}^T\mathrm{\Ω}\mathrm{\Δ}\mathrm{\α}\end{array}---\left(1\right)$

s.t.f_min/p-f₀≤Δf≤f_max/p-f₀

min{J＝f^TWf+s^TQs}(2)

s.t.Bf＝τ+s

在构建计算方位角α_k时的目标函数和约束条件中，新加入了避免奇异结构的推力折减系数p，去掉了传统方法中的一项：这样，可以避免复杂的计算，大大减少了计算量，而且，比传统方法更容易使推进器的方位角转到一个合适的位置，从而，降低了能量的消耗。另外，本算法在实时计算推进器推力f_k和方位角α_k时，是分开进行的，这样处理可以大大简化约束条件的非线性方程，使矩阵B(α)变为常数矩阵，避免了传统方法中将约束条件线性化时引入的误差。

由于采用优化算法计算推进器的方位角α_k时，会引入非线性项，这种新算法比传统算法巧妙之处在于，将推进器的推力和推进器的方位角α_k分开计算，避免了由于线性化造成的误差。而且，考虑了对推进器最大推力的限制，推力最大变化速度的限制，方位角α_k的范围，方位角变化的最大速度的限制，也避免了奇异结构的出现。另外，该算法基本不用考虑推进器的推力f_k可行域为非凸集的情况。该算法避免奇异结构的方法是加入推力折减系数p，折减系数p主要体现在采用序列二次规划法求推进器方位角α_k时，力的变化量的约束条件：f_min/p-f₀≤Δf≤f_max/p-f₀中，其中，f₀为：前一时刻推进器推力，Δf为当前时刻与前一时刻力的变化量，f_min、f_max分别对应推进器推力最小值、最大值。J_QP(Δα，Δf，s)为：关于Δα、Δf和s的目标函数；

Δα为：m维当前时刻与前一时刻的方位角变化向量；

Δf为：m维当前时刻与前一时刻的推力变化向量；

s为：3维松弛变量；

W_i为：关于第i个推进器功率；

f_i为：第i个推进器当前时刻推力；

f_0，i为第i个推进器前一时刻推力；

Δf_i为：第i个推进器当前时刻与前一时刻推力变化量；

Q为：3阶对角矩阵；

Ω为：m阶对角矩阵；

P为：推力折减系数；

J为：推力分配目标函数；

τ为：控制器；

W：为推进器个数；

$B=\left(\begin{array}{cccc}{\mathrm{cos\α}}_1& {\mathrm{cos\α}}_2& ...& {\mathrm{cos\α}}_m\\ {\mathrm{sin\α}}_1& {\mathrm{sin\α}}_2& ...& {\mathrm{sin\α}}_m\\ -y_1{\mathrm{cos\α}}_1+x_1{\mathrm{sin\α}}_1& -y_2{\mathrm{cos\α}}_2+x_2{\mathrm{sin\α}}_2& ...& -y_m{\mathrm{cos\α}}_m+x_m{\mathrm{sin\α}}_m\end{array}\right)$

α_i为第i个推进器的方位角，(x_i，y_i)为第i个推进器的坐标，坐标系原点为船舶重心。

这样的处理较其他方法计算量小，而且，更容易使推进器的方位角转到一个合适的位置，不但避免了奇异结构，而且，可以降低推进器的能量消耗。

以上所述，仅是本发明的较佳实施例而已，并非对本发明作任何形式上的限制，凡是依据本发明的技术实质对以上实施例所作的任何简单修改、等同变化与修饰，均仍属于本发明技术方案的范围内。

Claims (4)

1.一种动力定位系统的推力分配方法，其特征在于：采用以下步骤：

第一步：在进行推力分配之前，先将推进器分为数个小组，并固定组内推进器方位角的相对角度，但不固定它们的绝对角度；

第二步：由控制指令的低频部分使用优化算法求出一个缓慢变化的推进器方位角，在已知方位角后，再使用优化算法求出推进器的推力；

其具体步骤如下：

(1)首先，将k时刻的控制指令τ_c，k输入到一个由一阶惯性环节构成的低通滤波器中，然后，再进一步离散化可得缓慢变化的控制指令的低频部分τ_k+1；

(2)由控制指令τ_c，k采用序列二次规划法计算出各推进器的推力大小f_k和推进器方位角α_k；

(3)得到推进器方位角α_k后，再应用二次规划法求出推进器的推力f_k；其中，在计算推进器方位角α_k和推力f_k时，在所设置的目标函数和部分约束条件的公式中，加入推力折减系数p。

2.根据权利要求1所述的动力定位系统的推力分配方法，其特征在于：所述低通滤波器的形式为：其中，为低通滤波后的k时刻控制指令的一阶导数；τ_k为：低通滤波后的k时刻的控制指令，Λ为：一阶惯性环节低通滤波矩阵，且元素均为正的对角矩阵，τ_c，k为第k个时刻的控制指令；τ_k+1为：在第k个时刻预测的第k+1个时刻的控制指令的低频部分；为了便于实时计算，避免求解微分方程，利用基本的求导公式并令 ${\mathrm{\τ}}_{k+1}={\mathrm{\τ}}_k+h*{\stackrel{\·}{\mathrm{\τ}}}_k,$ 其中，h为：采样时间，而 ${\stackrel{\·}{\mathrm{\τ}}}_k=\mathrm{\Λ}({\mathrm{\τ}}_{c,k}-{\mathrm{\τ}}_k)$ 则与 ${\stackrel{\·}{\mathrm{\τ}}}_k=({\mathrm{\τ}}_{k+1}-$ ${\mathrm{\τ}}_k)/h$ 等价。

3.根据权利要求1所述的动力定位系统的推力分配方法，其特征在于：所述推进器推力f_k和推进器方位角α_k的计算是分开进行的。

4.根据权利要求1所述的动力定位系统的推力分配方法，其特征在于：所述目标函数和部分约束条件的具体公式为：

$\begin{array}{cc}\min & J_{QP}(\mathrm{\Δ}\mathrm{\α},\mathrm{\Δ}f,s)=\underset{i=1}{\overset{m}{\Σ}}\left(\frac{{\mathrm{dW}}_i}{{\mathrm{df}}_i}\right(f_{0,i}){\mathrm{\Δf}}_i+\frac{d^2{W}_i}{{\mathrm{df}}_i^2}(f_{0,i}\left){\mathrm{\Δf}}_i^2\right)+s^{T}Qs+{\mathrm{\Δ\α}}^T\mathrm{\Ω}\mathrm{\Δ}\mathrm{\α}---\left(1\right)\end{array}$

s.t.f_min/p-f₀≤Δf≤f_max/p-f₀

min{J＝f^TWf+s^TQs}(2)

s.t.Bf＝τ+s

上述公式：f_min/p-f₀≤Δf≤f_max/p-f₀中，f₀为：前一时刻推进器推力，Δf为m维当前时刻与前一时刻的推力变化向量，f_min、f_max分别对应推进器推力最小值、最大值；J_QP(Δα，Δf，s)为：关于Δα、Δf和s的目标函数；

Δα为：m维当前时刻与前一时刻的方位角变化向量；

Δf为：m维当前时刻与前一时刻的推力变化向量；

s为：3维松弛变量；

W_i为：关于第i个推进器功率；

f_i为：第i个推进器当前时刻推力；

f_0，i为第i个推进器前一时刻推力；

Δf_i为：第i个推进器当前时刻与前一时刻推力变化量；

Q为：3阶对角矩阵；

Ω为：m阶对角矩阵；

P为：推力折减系数；

J为：推力分配目标函数；

τ为：控制器输出的推力指令；

W：为推进器个数；

$B=\left(\begin{array}{cccc}{\mathrm{cos\α}}_1& {\mathrm{cos\α}}_2& ...& {\mathrm{cos\α}}_m\\ {\mathrm{sin\α}}_1& {\mathrm{sin\α}}_2& ...& {\mathrm{sin\α}}_m\\ -y_1{\mathrm{cos\α}}_1+x_1{\mathrm{sin\α}}_1& -y_2{\mathrm{cos\α}}_2+x_2{\mathrm{sin\α}}_2& ...& -y_m{\mathrm{cos\α}}_m+x_m{\mathrm{sin\α}}_m\end{array}\right)$

α_i为第i个推进器的方位角，i＝1-m，(x_i，y_i)为第i个推进器的坐标，坐标系原点为船舶重心。