随着海洋资源的开采，适用于海洋作业的船舶定位方式成为海洋界研究的焦点，常用的定位方式主要有锚泊定位和动力定位两种。但是锚泊定位随着定位水深的增加，定位精度低，成本高，安全性降低；动力定位系统则前期投入成本高，不符合节能减排的要求。针对这一问题，学者们提出了锚泊辅助动力定位系统，该定位系统不仅结合了锚泊定位和动力定位的优点，也保证了作业的安全，提高了全天候作业能力。

锚泊定位一般在水深300~1 000 m作业，动力定位则在水深1 000 m以上的海域作业。在锚泊辅助动力定位系统中，控制方法主要分为PID控制^[1]、切换控制^[2]和锚泊线集成控制^[3]。有学者将智能控制方法应用于混合系统的研究中，如：雷正玲等^[4]将基于扩张状态观测器的自抗扰控制技术应用于锚泊辅助动力定位系统的研究；Kwon等^[5]采用神经网络自适应控制方法对锚泊辅助动力定位系统控制器进行了设计。

本文将广义预测控制算法(general predictive control，GPC)应用于锚泊辅助动力定位系统的控制方法。采用隐式算法求解控制律，避免了求解丢番图方程，减少了计算量，并针对目标船进行了仿真验证。

1 船舶运动模型

在整个的锚泊辅助动力定位系统数学模型中，通常用到3种坐标系，分别是固定坐标系、船体坐标系和平行坐标系。对锚泊系统进行分析时，有时还会需要海底坐标系。建立船舶运动坐标系如图 1所示。对于定位船舶，一般情况下只考虑水平面三自由度运动，即纵荡、横荡和艏摇^[6-7]。

建立船舶的三自度运动学模型：

$ \mathit{\boldsymbol{\dot \eta }} = \mathit{\boldsymbol{R}}\left( \psi \right)\mathit{\boldsymbol{\nu }} $

(1)

其中，$\mathit{\boldsymbol{R}}\left(\psi \right)=\left[\begin{matrix} \text{cos}\psi &-\text{sin}\psi & 0 \\ \text{sin}\psi & \text{cos}\psi & 0 \\ 0 & 0 & 1 \\ \end{matrix} \right]$

船舶在海平面内的运动主要分为低频运动和波频运动两部分。低频运动由慢变环境力产生；波频运动由高频扰动力产生。

锚泊系统作用下的船舶低频动力学方程可表示为

$ \begin{array}{*{20}{c}} {\mathit{\boldsymbol{M\dot v}} + \mathit{\boldsymbol{C}}\left( \mathit{\boldsymbol{v}} \right)\mathit{\boldsymbol{v}} + \mathit{\boldsymbol{D}}\left( \mathit{\boldsymbol{v}} \right)\mathit{\boldsymbol{v}} + \mathit{\boldsymbol{g}}\left( \mathit{\boldsymbol{\eta }} \right) = }\\ {{\mathit{\boldsymbol{\tau }}_{{\rm{env}}}} + {\mathit{\boldsymbol{\tau }}_{{\rm{thrust}}}} + {\mathit{\boldsymbol{\tau }}_{{\rm{moor}}}}} \end{array} $

(2)

式中：v=[u v w]^T是速度向量；M是质量矩阵；C(v)是科里奥利向心力矩阵；D(v)是阻尼矩阵；g(η)是静态回复力及力矩；τ_env是海洋环境干扰力；τ_thrust是推进系统产生的力及力矩；τ_moor是锚链的总的张力。

2 隐式广义预测控制算法描述 2.1 CARIMA模型

广义预测控制采用的是受控自回归积分滑动平均(controlled auto-regressive integrated moving average，CARIMA)模型^[8]：

$ \begin{array}{l} A\left( {{z^{ - 1}}} \right)y\left( k \right) = B\left( {{z^{ - 1}}} \right)u\left( {k - 1} \right) + \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;C\left( {{z^{ - 1}}} \right)\omega \left( k \right)/\Delta \end{array} $

(3)

式中：y(k)和u(k)分别是系统的输出与输入；z^-1表示后移算子；其中A(z^-1)、B(z^-1)和C(z^-1)分别为：

$ A\left( {{z^{ - 1}}} \right) = 1 + {a_1}{z^{ - 1}} + {a_1}{z^{ - 2}} + \cdots + {a_{{n_a}}}{z^{ - {n_a}}} $

$ B\left( {{z^{ - 1}}} \right) = {b_0} + {b_1}{z^{ - 1}} + {b_1}{z^{ - 2}} + \cdots + {b_{{n_b}}}{z^{ - {n_b}}} $

$ C\left( {{z^{ - 1}}} \right) = 1 + {a_1}{z^{ - 1}} + {a_1}{z^{ - 2}} + \cdots + {a_{{n_a}}}{z^{ - {n_c}}} $

式中n_a为模型输入阶次；n_b和n_c为输出阶次。Δ=1-z^-1；ω(k)表示均值为零的白噪声扰动。

引入下式的丢番图(Diophantine)方程来预测超前时刻的未来变量：

$ 1 = {E_j}\left( {{z^{ - 1}}} \right)A\left( {{z^{ - 1}}} \right)\Delta + {z^{ - j}}{F_j}\left( {{z^{ - 1}}} \right) $

(4)

其中，$\left\{ \begin{align} & {{E}_{j}}({{z}^{-1}})={{e}_{0}}+{{e}_{1}}{{z}^{-1}}+\ldots \ldots +{{e}_{j-1}}{{z}^{-j+1}} \\ & {{F}_{j}}({{z}^{-1}})={{f}_{0}}^{j}+{{f}_{1}}^{j}{{z}^{-1}}+\ldots \ldots +{{f}_{na}}^{j}{{z}^{-na}} \\ \end{align} \right.$

根据式(3)、(4)，可以得到系统的预测输出为

$ \mathit{\boldsymbol{\hat y}} = \mathit{\boldsymbol{G}}\Delta \mathit{\boldsymbol{u}} + \mathit{\boldsymbol{f}} $

(5)

式中：

$ \mathit{\boldsymbol{\hat y}} = {\left[ {\hat y\left( {k + 1} \right),\hat y\left( {k + 2} \right), \cdots ,\hat y\left( {k + n} \right)} \right]^{\rm{T}}} $

$ \Delta \mathit{\boldsymbol{u}} = {\left[ {\Delta u\left( k \right),\Delta u\left( {k + 1} \right), \cdots ,\Delta u\left( {k + n - 1} \right)} \right]^{\rm{T}}} $

$ \mathit{\boldsymbol{f}} = {\left[ {f\left( {k + 1} \right),f\left( {k + 2} \right), \cdots ,f\left( {k + n} \right)} \right]^{\rm{T}}} $

$ {\mathit{\boldsymbol{G}}_{N \times M}} = \left[ {\begin{array}{*{20}{c}} {{g_0}}&{}&{}&{}\\ {{g_1}}&{{g_0}}&{}&{}\\ \vdots &{}&{}&{}\\ {{g_{M - 1}}}&{{g_{M - 2}}}& \cdots &{{g_0}}\\ \vdots &{}&{}&{}\\ {{g_{N - 1}}}&{{g_{N - 2}}}& \cdots &{{g_{N - M}}} \end{array}} \right] $

2.2 隐式自校正算法

预测控制的最终目标就是使被控对象的输出y(k+j)跟踪上期望输出y_r(k+j)。为了使被控对象的实际输出y(k)能够稳定的到达期望值，常用一阶滤波方程来实现平滑过渡，并用柔化因子α来调节柔化程度^[9-14]。

系统的性能指标：

$ \begin{array}{*{20}{c}} {J = \sum\limits_{j = 1}^N {{{\left[ {\hat y\left( {k + j} \right) - {y_r}\left( {k + j} \right)} \right]}^2}} + }\\ {\sum\limits_{j = 1}^M {{\lambda _j}{{\left[ {\Delta u\left( {k + j - 1} \right)} \right]}^2}} } \end{array} $

(6)

令$\bar{v}\left(k+j \right)={{y}_{r}}\left(k+j \right)$和$\frac{\partial J}{\partial \Delta \mathit{\boldsymbol{u}}}=0$，求解得到最优控制律：

$ \Delta \mathit{\boldsymbol{u}} = {\left( {{\mathit{\boldsymbol{G}}^{\rm{T}}}\mathit{\boldsymbol{G}} + \lambda \mathit{\boldsymbol{I}}} \right)^{ - 1}}{\mathit{\boldsymbol{G}}^{\rm{T}}}\left( {\mathit{\boldsymbol{\bar v}} - \mathit{\boldsymbol{f}}} \right) $

(7)

由式(7)可以看出，欲求Δu由G和f控制。隐式自校正方法就是利用输入输出数据，根据预测直接求得矩阵G和f。

根据式(5)可以得到n个并列预测器：

$ \begin{array}{*{20}{c}} {y\left( {k + 1} \right) = {g_0}\Delta u\left( k \right) + f\left( {k + 1} \right) + {E_1}\xi \left( {k + 1} \right)}\\ {y\left( {k + 2} \right) = {g_1}\Delta u\left( k \right) + {g_0}\Delta u\left( {k + 1} \right) + }\\ {f\left( {k + 2} \right) + {E_2}\xi \left( {k + 2} \right)}\\ \cdots \\ {y\left( {k + n} \right) = {g_{n - 1}}\Delta u\left( k \right) + {g_0}\Delta u\left( {k + n - 1} \right) + }\\ {f\left( {k + n} \right) + {E_n}\xi \left( {k + n} \right)} \end{array} $

(8)

由式(8)可以看出，只需要对上式的最后一个方程辨识即可得到矩阵G。

将上式的最后一个方程写为:

$ y\left( {k + n} \right) = \mathit{\boldsymbol{\omega }}\left( k \right)\mathit{\boldsymbol{\theta }}\left( k \right) + {E_n}\xi \left( {k + n} \right) $

(9)

其中：

$ \mathit{\boldsymbol{\omega }}\left( k \right) = \left[ {\Delta u\left( k \right),\Delta u\left( {k + 1} \right), \cdots ,\Delta u\left( {k + n - 1} \right),1} \right] $

$ \mathit{\boldsymbol{\theta }}\left( k \right) = {\left[ {\begin{array}{*{20}{c}} {{g_{n - 1}},}&{{g_{n - 2}},}&{ \cdots ,}&{{g_0},}&{f\left( {k + n} \right)} \end{array}} \right]^{\rm{T}}} $

所以输出预测值为:

$ y\left( {k + n\left| k \right.} \right) = \mathit{\boldsymbol{\omega }}\left( k \right)\mathit{\boldsymbol{\theta }}\left( k \right) $

(10)

根据递推最小二乘法，即可得到矩阵G和向量f中的元素。

3 锚泊辅助动力定位系统隐式GPC设计 3.1 混合定位系统GPC预测模型建立

本文研究的锚泊辅助动力定位系统采用CARIMA模型作为预测模型。欲得到船舶预测模型须进行如下操作：

1) 对${\mathit{\boldsymbol{\dot{v}}}}$线性离散化。

锚泊定位时，船舶速度很小，受速度影响的矩阵项可以忽略不计，简化的数学模型^[15]：

$ \mathit{\boldsymbol{M\dot v}} + \mathit{\boldsymbol{Dv}} = {\mathit{\boldsymbol{\tau }}_{{\rm{thrust}}}} + {\mathit{\boldsymbol{\tau }}_{{\rm{ext}}}} $

(11)

其中，${\mathit{\boldsymbol{\tau }}_{{\rm{ext}}}} = {\mathit{\boldsymbol{\tau }}_{{\rm{wind}}}} + {\mathit{\boldsymbol{\tau }}_{{\rm{current}}}} + {\rm{ }}{\mathit{\boldsymbol{\tau }}_{{\rm{wave}}}} + {\mathit{\boldsymbol{\tau }}_{{\rm{moor}}}}$

将式(11)变形可得：

$ \mathit{\boldsymbol{\dot v}} = - {\mathit{\boldsymbol{M}}^{ - 1}}\mathit{\boldsymbol{Dv}} + {\mathit{\boldsymbol{M}}^{ - 1}}\left( {{\mathit{\boldsymbol{\tau }}_{{\rm{thrust}}}} + {\mathit{\boldsymbol{\tau }}_{{\rm{ext}}}}} \right) $

(12)

令$\mathit{\boldsymbol{A}} =-{\mathit{\boldsymbol{M}}^{-1}}\mathit{\boldsymbol{D}}, \mathit{\boldsymbol{B}} = {\mathit{\boldsymbol{M}}^{-1}}, \mathit{\boldsymbol{U}} = {\mathit{\boldsymbol{\tau }}_{{\rm{thr}}}} + {\mathit{\boldsymbol{\tau }}_{{\rm{ext}}}}$，则式(12)可以写成

$ \mathit{\boldsymbol{\dot v}} = - \mathit{\boldsymbol{Av}} + \mathit{\boldsymbol{BU}} $

(13)

对${\mathit{\boldsymbol{\dot v}}}$进行线性离散化并代入式(14)，取采样时间为T：

$ \mathit{\boldsymbol{v}}\left( {k + 1} \right) = \left( {\mathit{\boldsymbol{I}} + T\mathit{\boldsymbol{A}}} \right)\mathit{\boldsymbol{v}}\left( k \right) + T\mathit{\boldsymbol{BU}}\left( k \right) $

(14)

令G=I+TA, H=TB则

$ \mathit{\boldsymbol{v}}\left( {k + 1} \right) = \mathit{\boldsymbol{Gv}}\left( k \right) + \mathit{\boldsymbol{HU}}\left( k \right) $

(15)

2) 求v(k)的表达式。

船舶运动过程中，需要控制的是船舶位置和艏向η=[x, y, ψ]^T，而不是速度向量v=[u, v, r]^T，且在锚泊辅助动力定位系统中，船舶艏向看作是固定值^[16-19]，所以可以对运动方程式(1)变形得到:

$ \mathit{\boldsymbol{R}}{\left( \psi \right)^{\rm{T}}}\mathit{\boldsymbol{\dot \eta }} = \mathit{\boldsymbol{\nu }} $

(16)

令Y=R(ψ)η，并式(17)两边对时间积分, 取积分时间为一个采样周期T，则有：

$ \mathit{\boldsymbol{Y}}\left( t \right) = \int\limits_k^{k + T} {\mathit{\boldsymbol{v}}\left( t \right){\rm{d}}t} $

(17)

整理得：

$ \frac{{\mathit{\boldsymbol{Y}}\left( {k + 1} \right) - \mathit{\boldsymbol{Y}}\left( k \right)}}{T} = \frac{{\mathit{\boldsymbol{v}}\left( {k + 1} \right) + \mathit{\boldsymbol{v}}\left( k \right)}}{2} $

(18)

将式(15)代入到式(18)中可得：

$ \begin{array}{l} \mathit{\boldsymbol{v}}\left( k \right) = {\left( {\mathit{\boldsymbol{G}} + \mathit{\boldsymbol{I}}} \right)^{ - 1}}\left( {\frac{2}{T}\mathit{\boldsymbol{Y}}\left( {k + 1} \right) - \mathit{\boldsymbol{Y}}\left( k \right)} \right) - \\ \;\;\;\;\;\;\;\;\;\;\left. {\mathit{\boldsymbol{HU}}\left( k \right)} \right) \end{array} $

(19)

3) 求Y(k)的表达式。

将式(19)往后递推两个采样周期：

$ \begin{array}{l} \mathit{\boldsymbol{v}}\left( {k - 1} \right) = {\left( {\mathit{\boldsymbol{G}} + \mathit{\boldsymbol{I}}} \right)^{ - 1}}\left( {\frac{2}{T}\mathit{\boldsymbol{Y}}\left( k \right) - \mathit{\boldsymbol{Y}}\left( {k - 1} \right)} \right) - \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\left. {\mathit{\boldsymbol{HU}}\left( {k - 1} \right)} \right) \end{array} $

(20)

同理可得v(k-2)的表达式，联立式(15)消掉v(k)可得：

$ \begin{array}{l} \mathit{\boldsymbol{Y}}\left( k \right) - \left( {\mathit{\boldsymbol{I}} + \left( {\mathit{\boldsymbol{G}} + \mathit{\boldsymbol{I}}} \right)\mathit{\boldsymbol{G}}{{\left( {\mathit{\boldsymbol{G}} + \mathit{\boldsymbol{I}}} \right)}^{ - 1}}} \right)\mathit{\boldsymbol{Y}}\left( {k - 1} \right) + \\ \;\;\;\;\;\;\;\;\;\;\left( {\mathit{\boldsymbol{G}} + \mathit{\boldsymbol{I}}} \right)\mathit{\boldsymbol{G}}{\left( {\mathit{\boldsymbol{G}} + \mathit{\boldsymbol{I}}} \right)^{ - 1}}\mathit{\boldsymbol{Y}}\left( {k - 2} \right) = \\ \;\;\;\;\;\;\;\;\;\;\frac{T}{2}\mathit{\boldsymbol{HU}}\left( {k - 1} \right) + \frac{T}{2}\left( {\mathit{\boldsymbol{G}} + \mathit{\boldsymbol{I}}} \right)\mathit{\boldsymbol{H}} - \\ \;\;\;\;\;\;\;\;\;\;\left. {\left( {\mathit{\boldsymbol{G}} + \mathit{\boldsymbol{I}}} \right)\mathit{\boldsymbol{G}}{{\left( {\mathit{\boldsymbol{G}} + \mathit{\boldsymbol{I}}} \right)}^{ - 1}}\mathit{\boldsymbol{H}}} \right)\mathit{\boldsymbol{U}}\left( {k - 2} \right) \end{array} $

(21)

4) 求船舶CARIMA预测模型与式(3)相对比，可以得到：

$ \begin{array}{*{20}{c}} {\left( {\mathit{\boldsymbol{I}} + {\mathit{\boldsymbol{A}}_1}{\mathit{\boldsymbol{z}}^{ - 1}} + {\mathit{\boldsymbol{A}}_2}{\mathit{\boldsymbol{z}}^{ - 2}}} \right)\mathit{\boldsymbol{Y}}\left( k \right) = }\\ {\left( {{\mathit{\boldsymbol{B}}_0} + {\mathit{\boldsymbol{B}}_1}{\mathit{\boldsymbol{z}}^{ - 1}}} \right)\mathit{\boldsymbol{\tau }}\left( {k - 1} \right) + \mathit{\boldsymbol{C}}\left( {{\mathit{\boldsymbol{z}}^{ - 1}}} \right){\mathit{\boldsymbol{\tau }}_{{\rm{ext}}}}\left( {k - 1} \right)} \end{array} $

(22)

其中，n_a=2, n_b=1且

$ {\mathit{\boldsymbol{A}}_1} = - \left( {\mathit{\boldsymbol{I}} + \left( {\mathit{\boldsymbol{G}} + \mathit{\boldsymbol{I}}} \right)\mathit{\boldsymbol{G}}{{\left( {\mathit{\boldsymbol{G}} + \mathit{\boldsymbol{I}}} \right)}^{ - 1}}} \right), $

$ {\mathit{\boldsymbol{A}}_2} = \left( {\mathit{\boldsymbol{G}} + \mathit{\boldsymbol{I}}} \right)\mathit{\boldsymbol{G}}{\left( {\mathit{\boldsymbol{G}} + \mathit{\boldsymbol{I}}} \right)^{ - 1}}, $

$ {\mathit{\boldsymbol{B}}_0} = \frac{T}{2}\mathit{\boldsymbol{H}}, $

$ {\mathit{\boldsymbol{B}}_1} = \frac{T}{2}\left( {\left( {\mathit{\boldsymbol{G}} + \mathit{\boldsymbol{I}}} \right)\mathit{\boldsymbol{H}} - \left( {\mathit{\boldsymbol{G}} + \mathit{\boldsymbol{I}}} \right)\mathit{\boldsymbol{G}}{{\left( {\mathit{\boldsymbol{G}} + \mathit{\boldsymbol{I}}} \right)}^{ - 1}}\mathit{\boldsymbol{H}}} \right), $

$ \mathit{\boldsymbol{A}}\left( {{\mathit{\boldsymbol{z}}^{ - 1}}} \right) = \mathit{\boldsymbol{I}} + {\mathit{\boldsymbol{A}}_1}{\mathit{\boldsymbol{z}}^{ - 1}} + {\mathit{\boldsymbol{A}}_2}{\mathit{\boldsymbol{z}}^{ - 2}}, $

$ \mathit{\boldsymbol{B}}\left( {{\mathit{\boldsymbol{z}}^{ - 1}}} \right) = {\mathit{\boldsymbol{B}}_0} + {\mathit{\boldsymbol{B}}_1}{\mathit{\boldsymbol{z}}^{ - 1}},\mathit{\boldsymbol{C}}\left( {{\mathit{\boldsymbol{z}}^{ - 1}}} \right) = \mathit{\boldsymbol{B}}\left( {{\mathit{\boldsymbol{z}}^{ - 1}}} \right) $

将式(22)写成下式即为船舶CARIMA模型：

$ \begin{array}{l} \mathit{\boldsymbol{\tilde A}}\left( {{\mathit{\boldsymbol{z}}^{ - 1}}} \right)\mathit{\boldsymbol{Y}}\left( k \right) = \mathit{\boldsymbol{\tilde B}}\left( {{\mathit{\boldsymbol{z}}^{ - 1}}} \right)\mathit{\boldsymbol{\tau }}\left( {k - 1} \right) + \mathit{\boldsymbol{\zeta }}\left( k \right)\\ \mathit{\boldsymbol{\zeta }}\left( k \right) = \mathit{\boldsymbol{C}}\left( {{\mathit{\boldsymbol{z}}^{ - 1}}} \right){\mathit{\boldsymbol{\tau }}_{{\rm{ext}}}}\left( {k - 1} \right) + \mathit{\boldsymbol{\omega }}\left( k \right) \end{array} $

(23)

式中:ζ(k)表示随机扰动; ω(k)表示零均值高斯白噪声，且

$ \mathit{\boldsymbol{\tilde A}}\left( {{\mathit{\boldsymbol{z}}^{ - 1}}} \right) = \left[ {\begin{array}{*{20}{c}} {{{\mathit{\boldsymbol{\tilde A}}}_1}\left( {{\mathit{\boldsymbol{z}}^{ - 1}}} \right)}&0&0\\ 0&{{{\mathit{\boldsymbol{\tilde A}}}_2}\left( {{\mathit{\boldsymbol{z}}^{ - 1}}} \right)}&0\\ 0&0&{{{\mathit{\boldsymbol{\tilde A}}}_3}\left( {{\mathit{\boldsymbol{z}}^{ - 1}}} \right)} \end{array}} \right] $

$ \mathit{\boldsymbol{\tilde B}}\left( {{\mathit{\boldsymbol{z}}^{ - 1}}} \right) = \left[ {\begin{array}{*{20}{c}} {{{\tilde B}_1}\left( {{z^{ - 1}}} \right)}&0&0\\ 0&{{{\tilde B}_2}\left( {{z^{ - 1}}} \right)}&0\\ 0&0&{{{\tilde B}_3}\left( {{z^{ - 1}}} \right)} \end{array}} \right] $

$ {{\mathit{\boldsymbol{\tilde A}}}_i}\left( {{\mathit{\boldsymbol{z}}^{ - 1}}} \right) = 1 + {a_{i1}}{z^{ - 1}} + \cdots + {a_{i{n_a}}}{z^{ - na}}\left( {i = 1,2,3} \right) $

$ {{\mathit{\boldsymbol{\tilde B}}}_i}\left( {{\mathit{\boldsymbol{z}}^{ - 1}}} \right) = {b_{i0}} + {b_{i1}}{z^{ - 1}} + \cdots + {b_{i{n_b}}}{z^{ - nb}} $

3.2 混合定位系统隐式自校正方法设计

采用隐式算法求解预测控制器的设计步骤：

1) 预测输出的求解。

欲得到船舶隐式广义预测控制的预测输出，需引入丢番图方程^[20]，联立Diophantine方程和式(23)可得预测输出矢量为:

$ \begin{array}{l} \mathit{\boldsymbol{Y}}\left( {k + j} \right) = {\mathit{\boldsymbol{E}}_j}\left( {{\mathit{\boldsymbol{z}}^{ - 1}}} \right)\mathit{\boldsymbol{B}}\left( {{\mathit{\boldsymbol{z}}^{ - 1}}} \right)\Delta \tau \left( {k + j - 1} \right) + \\ \;\;\;\;\;\;\;\;\;\;\;\;\;{\mathit{\boldsymbol{F}}_j}\left( {{\mathit{\boldsymbol{z}}^{ - 1}}} \right)\mathit{\boldsymbol{Y}}\left( k \right) + {\mathit{\boldsymbol{E}}_j}\left( {{z^{ - 1}}} \right)\Delta \zeta \left( {k + j} \right) \end{array} $

(24)

由于k时刻以后的系统扰动是不确定的，所以最优预测输出为:

$ \begin{array}{l} \mathit{\boldsymbol{\hat Y}}\left( {k + j} \right) = {\mathit{\boldsymbol{E}}_j}\left( {{\mathit{\boldsymbol{z}}^{ - 1}}} \right)\mathit{\boldsymbol{B}}\left( {{\mathit{\boldsymbol{z}}^{ - 1}}} \right)\Delta \mathit{\boldsymbol{\tau }}\left( {k + j - 1} \right) + \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;{\mathit{\boldsymbol{F}}_j}\left( {{\mathit{\boldsymbol{z}}^{ - 1}}} \right)\mathit{\boldsymbol{Y}}\left( k \right) \end{array} $

(25)

也可写为:

$ \begin{array}{l} \mathit{\boldsymbol{\hat Y}}\left( {k + j} \right) = {\mathit{\boldsymbol{G}}_j}\left( {{\mathit{\boldsymbol{z}}^{ - 1}}} \right)\Delta \mathit{\boldsymbol{\tau }}\left( {k + j - 1} \right) + \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;{\mathit{\boldsymbol{F}}_j}\left( {{\mathit{\boldsymbol{z}}^{ - 1}}} \right)\mathit{\boldsymbol{Y}}\left( k \right) \end{array} $

(26)

其中，因为锚泊辅助动力定位系统是多变量系统，所以E_j(z^-1)和F_j(z^-1)为：

$ {\mathit{\boldsymbol{E}}_j}\left( {{\mathit{\boldsymbol{z}}^{ - 1}}} \right) = \left[ {\begin{array}{*{20}{c}} {{E_{1j}}\left( {{z^{ - 1}}} \right)}&0&0\\ 0&{{E_{2j}}\left( {{z^{ - 1}}} \right)}&0\\ 0&0&{{E_{3j}}\left( {{z^{ - 1}}} \right)} \end{array}} \right] $

$ {\mathit{\boldsymbol{F}}_j}\left( {{\mathit{\boldsymbol{z}}^{ - 1}}} \right) = \left[ {\begin{array}{*{20}{c}} {{F_{1j}}\left( {{z^{ - 1}}} \right)}&0&0\\ 0&{{F_{2j}}\left( {{z^{ - 1}}} \right)}&0\\ 0&0&{{F_{3j}}\left( {{z^{ - 1}}} \right)} \end{array}} \right] $

各对角线上的多项式元素为：

$ \begin{array}{*{20}{c}} {{E_{ij}}\left( {{z^{ - 1}}} \right) = {e_{i,0}} + {e_{i,1}}{z^{ - 1}} + \cdots + {e_{i,j - 1}}{z^{ - j + 1}},}\\ {i = 1,2,3;} \end{array} $

$ \begin{array}{*{20}{c}} {{F_{ij}}\left( {{z^{ - 1}}} \right) = {f_{ij,0}} + {f_{ij,1}}{z^{ - 1}} + \cdots + {f_{ij,na}}{z^{ - na}},}\\ {j = 1,2, \cdots ,N;} \end{array} $

$ \begin{array}{*{20}{c}} {{\mathit{\boldsymbol{G}}_j}\left( {{\mathit{\boldsymbol{z}}^{ - 1}}} \right) = {\mathit{\boldsymbol{E}}_j}\left( {{\mathit{\boldsymbol{z}}^{ - 1}}} \right)\mathit{\boldsymbol{B}}\left( {{\mathit{\boldsymbol{z}}^{ - 1}}} \right) = }\\ {\left[ {\begin{array}{*{20}{c}} {{G_{1j}}\left( {{z^{ - 1}}} \right)}&0&0\\ 0&{{G_{2j}}\left( {{z^{ - 1}}} \right)}&0\\ 0&0&{{G_{3j}}\left( {{z^{ - 1}}} \right)} \end{array}} \right]} \end{array} $

对角线上的多项式元素为:

$ {G_{ij}}\left( {{z^{ - 1}}} \right) = {E_{ij}}{B_i} = {g_{i1}} + {g_{i2}}{z^{ - 1}} + \cdots $

根据式(26)可知，在确定好未来控制策略的前提下，由预测模型便可得到船舶未来时刻的位置和艏向。

2) 最优控制律的求解。

设Y(k)=[Y₁(k) Y₂(k) Y₃(k)]^T为船舶在未来j个时刻船舶位置和艏向参考值组成的向量，其中：

$ {\mathit{\boldsymbol{Y}}_i}\left( k \right) = \left[ {\begin{array}{*{20}{c}} {{y_i}\left( {k + 1} \right)}&{{y_i}\left( {k + 2} \right)}& \cdots &{{y_i}\left( {k + j} \right)} \end{array}} \right]_{1 \times j}^{\rm{T}} $

其中，Y_i(k)为未来j个时刻船舶的期望位置和艏向组成的向量，其中：

$ {y_i}\left( {k + l} \right) = \alpha _i^l{y_{si}}\left( k \right) + \left( {1 - \alpha _i^l} \right){y_{ri}},l = 1,2, \cdots ,j $

式中：α_i为船舶期望位置和艏向的柔化因子，y_ri为船舶位置和艏向的期望值，y_si(k)为k时刻船舶位置和艏向的实际输出值。

取k时刻的优化性能指标为:

$ \begin{array}{l} \mathit{\boldsymbol{J}}\left( k \right) = {\left( {{\mathit{\boldsymbol{Y}}_r}\left( k \right) - \mathit{\boldsymbol{\dot Y}}\left( k \right)} \right)^{\rm{T}}}\mathit{\boldsymbol{Q}}\left( {{\mathit{\boldsymbol{Y}}_r}\left( k \right) - \mathit{\boldsymbol{\dot Y}}\left( k \right)} \right) + \\ \;\;\;\;\;\;\;\;\;\;\;\Delta {\mathit{\boldsymbol{\tau }}^{\rm{T}}}\left( k \right)\mathit{\boldsymbol{R}}\Delta \mathit{\boldsymbol{\tau }}\left( k \right) \end{array} $

(27)

即$\mathit{\boldsymbol{J}} = {({\mathit{\boldsymbol{Y}}_r}-\mathit{\boldsymbol{\hat Y}})^{\rm{T}}}\mathit{\boldsymbol{Q}}({\mathit{\boldsymbol{Y}}_r}-\mathit{\boldsymbol{\hat Y}}) + \Delta {\mathit{\boldsymbol{\tau }}^{\rm{T}}}\mathit{\boldsymbol{R}}\Delta \mathit{\boldsymbol{\tau }}$

加权系数矩阵Q、R分别表示为：

$ \begin{array}{*{20}{c}} {\mathit{\boldsymbol{Q}} = {\rm{diag}}\left( {{Q_1},{Q_2}, \cdots ,{Q_N}} \right),{\mathit{\boldsymbol{Q}}_m} = {\rm{diag}}\left( {\begin{array}{*{20}{c}} {{q_1}}&{{q_2}}&{{q_3}} \end{array}} \right),}\\ {\left( {m = 1,2, \cdots ,N} \right)} \end{array} $

$ \begin{array}{*{20}{c}} {\mathit{\boldsymbol{R}} = {\rm{diag}}\left( {{r_1},{r_2}, \cdots ,{r_M}} \right),{\mathit{\boldsymbol{r}}_l} = {\rm{diag}}\left( {\begin{array}{*{20}{c}} {{r_1}}&{{r_2}}&{{r_3}} \end{array}} \right),}\\ {\left( {l = 1,2, \cdots ,M} \right)} \end{array} $

令$\frac{{\partial \mathit{\boldsymbol{J}}}}{{\partial \Delta \mathit{\boldsymbol{\tau }}}} = 0$可得最优控制律为:

$ \Delta \mathit{\boldsymbol{\tau }}\left( k \right) = {\left( {{\mathit{\boldsymbol{G}}^{\rm{T}}}\mathit{\boldsymbol{QG}} + \mathit{\boldsymbol{R}}} \right)^{ - 1}}{\mathit{\boldsymbol{G}}^{\rm{T}}}\mathit{\boldsymbol{Q}}\left( {{\mathit{\boldsymbol{Y}}_r} - \mathit{\boldsymbol{f}}\left( k \right)} \right) $

(28)

则当前时刻的最优控制量为:

$ \mathit{\boldsymbol{\tau }}\left( k \right) = \mathit{\boldsymbol{\tau }}\left( {k - 1} \right) + \mathit{\boldsymbol{d}}\Delta \mathit{\boldsymbol{\tau }}\left( k \right) $

(29)

其中, d=[1, 0, …, 0]_1×M。

根据2.2节讲述的隐式广式预测控制算法，利用递推最小二乘法^[21]，求出矩阵G，即可求出最优控制律。

4 船舶隐式GPC仿真结果

设定船舶期望位置和艏向为[0  0  30°]，控制算法的相关参数设计如下：采样周期T=0.4 s，预测步长N=120，控制步长M=5；遗忘因子γ₁=γ₂=γ₃=0.98，柔化系数α₁=α₂=α₃=0.3辨识参数P₁=P₂=P₃=10⁶，θ₁=θ₂=θ₃=10^-5；加权系数r₁=1×10^-4，r₂=1.3×10^-5，r₃=1.9×10^-6；q₁=q₂=q₃=1。分别采用隐式GPC算法、GPC算法和工程中常用的PID控制方法对此系统进行了设计，三、四级海况仿真结果对比如图 2~5所示。

三级海况中风(10.3 m/s，0°)，浪(1.5 m，0°)，初始北向位置：-5 m，初始东向位置：-5 m，艏向30°保持不变。图 2(a)~(c)分别表示的是船舶北向、东向位置和艏向的输出响应结果，可以看出采用隐式GPC算法的位置和艏向输出响应要优于采用PID算法。隐式GPC与传统GPC的位置和艏向响应曲线大致一样。图 3(a)~(c)分别表示的是船舶纵向、横向控制力和艏向控制力矩输出。从图中可以看出采用隐式GPC算法的控制力和力矩输出振幅比采用PID算法小。并且采用隐式GPC算法的控制力和力矩输出振幅也比采用传统GPC算法小。

四级海况中，风(13.5 m/s，0°)，浪(2.0 m，0°)，船舶的初始北向位置：-5 m，初始东向位置：-5 m，艏向30°保持不变。下面的图 4、图 5分别表示的是船舶的北向、东向位置和艏向的输出响应结果。图 4、图 5分别表示的是船舶纵向、横向控制力和艏向控制力矩输出。从图中可以看出，四级海况时采用的GPC算法也能较好的对船舶进行定位。

综上所述，采用的隐式广义预测控制算法能很好地跟踪船舶锚泊辅助动力定位系统的位置和艏向，并且稳定性表现也不错。

5 结论

1) 根据算法的基本原理得到了混合定位系统控制器的最优控制律。

2) 仿真通过与传统的PID算法和传统的GPC算法相对验证了所采用的隐式GPC算法能有效的改善控制系统的实时性问题，从而提高了混合定位系统的性能。