动力定位(dynamic positioning，DP)船舶是一种重要的海洋工程船舶，可应用于管道铺设、援潜救生等多种工程作业中^[1]。船舶在定位过程中受外界扰动的影响，选择适当的滤波方法，将其与控制方法相结合用于控制器设计有着重要的意义。船舶DP系统的滤波问题，通常采用扩展卡尔曼滤波^[2](extended Kalman filtering，EKF)方法来解决，该方法采用线性化模型进行滤波，但线性化可能带来误差过大和计算困难等问题。文献^[3]采用粒子滤波(particle filtering，PF)方法估计船舶的状态，该方法无需对模型进行线性化，但存在样本退化和对量测模型要求较高等问题。文献^[4]中给出了一种基于UKF的滚动时域估计(moving horizon estimation，MHE)方法，该方法通过无迹卡尔曼滤波(unscented Kalman filtering，UKF)计算到达代价，将状态估计问题转化为有限时域优化问题，利用在线滚动优化使得估计问题得到动态满足^[5]。该方法采用确定性采样方法，避免了样本退化，并且具有无需线性化的特点，能够避免线性化带来的误差过大和计算困难等问题。

船舶的动力学特性具有明显的非线性特征，为了解决船舶DP控制问题，科研人员将非线性预测控制方法引入到船舶DP控制系统中，例如基于支持向量机的模型预测控制^[6]和切换解析模型预测控制^[7](switch analytic model predictive control,SAMPC)等。SAMPC具有模型预测控制的特点，并吸取了反馈线性化的优点，但是运用该方法设计的控制器需要在不同状态间切换，可能引起系统振荡。张国银等^[8]在前人的工作基础上提出了基于相关度的非切换解析模型预测控制算法(non-switch analytic model predictive control,NSAMPC)，避免了控制器在不同状态间切换。王元慧等^[9]采用该算法设计动力定位非线性控制器，取得了较好的控制效果。

鉴于基于UKF的MHE与NSAMPC的特点，本文将二者相结合，提出一种船舶DP非线性控制器设计方法。对承受外界扰动的某供应船，采用本文所提出的方法设计DP非线性控制器，并进行仿真验证。

1 动力定位船舶运动模型

对于水面DP船舶，只考虑纵荡、横荡和艏摇这三个水平面上的自由度的运动。假设船舶运动模型的惯性矩阵M₁和阻尼矩阵D为已知，并且为定常矩阵，船舶做低速运动，可以忽略二阶项，得到三自由度船舶动力学模型^[9]:

$\left\{ \begin{matrix} \dot{\eta }=R\left( \psi \right)\nu \\ {{M}_{1}}\dot{\nu }+D\nu =\tau \\ \end{matrix} \right.$

(1)

式中：η=[x,y,ψ]^T为北东坐标系下船舶的位置和艏向角度向量；ν=[u,v,r]^T为随船坐标系下船舶的纵荡、横荡和艏摇角速度向量；R(ψ)为北东坐标系和随船坐标之间的转换矩阵；τ=[τ_X,τ_Y,τ_N]^T代表控制向量，由推进器在纵荡、横荡和艏摇三个方向上产生的力和力矩组成；M₁是惯性矩阵，满足正定要求M₁=M₁^T ＞0；D是水动力阻尼矩阵，严格正定。本文中所有关于时间t的函数在不产生歧义的情况下简写为函数名，例如，x(t)简写为x。矩阵定义：

$R\left( \psi \right)\left[ \begin{matrix} cos\psi & -sin\psi & 0 \\ sin\psi & cos\psi & 0 \\ 0 & 0 & 1 \\ \end{matrix}= \right]$

(2)

$D=\left[ \begin{matrix} {{d}_{11}} & 0 & 0 \\ 0 & {{d}_{22}} & {{d}_{23}} \\ 0 & {{d}_{32}} & {{d}_{33}} \\ \end{matrix} \right]$

(3)

${{M}_{1}}=\left[ \begin{matrix} {{m}_{11}} & 0 & 0 \\ 0 & {{m}_{22}} & {{m}_{23}} \\ 0 & {{m}_{32}} & {{m}_{33}} \\ \end{matrix} \right]$

(4)

为了便于控制器设计，根据式(2)~(4)，将式(1)转换为微分方程组形式: x·=ucos ψ－vsin ψy·=usin ψ+vcos ψψ·=r

$\left\{ \begin{matrix} \dot{x}=ucos\psi -vsin\psi \\ \dot{y}=usin\psi +vcos\psi \\ \dot{\psi }=r \\ \end{matrix} \right.$

(5)

$\left\{ \begin{matrix} {{y}_{1}}=x \\ {{y}_{2}}=y \\ {{y}_{3}}=\psi \\ \end{matrix} \right.$

(6)

$\left\{ \begin{matrix} \dot{u}=-\frac{{{d}_{11}}}{{{m}_{11}}}u+\frac{{{\tau }_{X}}}{{{m}_{11}}} \\ \dot{v}=-\frac{{{m}_{23}}}{{{m}_{22}}}\dot{r}-\frac{{{d}_{22}}}{{{m}_{22}}}v-\frac{{{d}_{23}}}{{{m}_{22}}}r+\frac{{{\tau }_{Y}}}{{{m}_{22}}} \\ \dot{r}=-\frac{{{m}_{32}}}{{{m}_{33}}}\dot{v}-\frac{{{d}_{32}}}{{{m}_{33}}}v-\frac{{{d}_{33}}}{{{m}_{33}}}r+\frac{{{\tau }_{N}}}{{{m}_{33}}} \\ \end{matrix} \right.$

(7)

将系统写成如下的规范形式：

$\left\{ \begin{matrix} \dot{x}=f\left( x \right)+g\left( x \right)u \\ y=h\left( x \right) \\ \end{matrix} \right.$

(8)

式中：x∈Rⁿ为系统的状态向量；u∈R^l为系统的控制输入；y∈R^m为系统输出；本文取n=6，l=3，m=3。

$x={{[{{x}_{1}},{{x}_{2}},{{x}_{3}},{{x}_{4}},{{x}_{5}},{{x}_{6}}]}^{T}}={{\left[ x,y,\psi ,u,v,r \right]}^{T}}$

(9)

$f\left( x \right)=\left[ \begin{matrix} {{x}_{4}}cos{{x}_{3}}-{{x}_{5}}sin{{x}_{3}} \\ {{x}_{4}}sin{{x}_{3}}+{{x}_{5}}cos{{x}_{3}} \\ {{x}_{6}} \\ {{a}_{11}}{{x}_{4}} \\ {{c}_{11}}{{{\dot{x}}}_{6}}+{{a}_{21}}{{x}_{5}}+{{a}_{22}}{{x}_{6}} \\ {{c}_{22}}{{{\dot{x}}}_{5}}+{{a}_{31}}{{x}_{5}}+{{a}_{32}}{{x}_{6}} \\ \end{matrix} \right]$

(10)

$g\left( x \right)=[{{g}_{1}}\left( x \right),{{g}_{2}}\left( x \right),{{g}_{3}}\left( x \right)]={{\left[ \begin{matrix} 0 & 0 & 0 & {{b}_{11}} & 0 & 0 \\ 0 & 0 & 0 & 0 & {{b}_{22}} & 0 \\ 0 & 0 & 0 & 0 & 0 & {{b}_{33}} \\ \end{matrix} \right]}^{T}}$

(11)

$u={{[{{u}_{1}},{{u}_{2}},{{u}_{3}}]}^{T}}={{[{{\tau }_{X}},{{\tau }_{Y}},{{\tau }_{N}}]}^{T}}$

(12)

$y={{[{{y}_{1}},{{y}_{2}},{{y}_{3}}]}^{T}}$

(13)

$h\left( x \right)={{[{{h}_{1}}\left( x \right),{{h}_{2}}\left( x \right),{{h}_{3}}\left( x \right)]}^{T}}={{\left[ x,y,\psi \right]}^{T}}$

(14)

式中：$\begin{align} & {{a}_{11}}=-\frac{{{d}_{11}}}{{{m}_{11}}},{{a}_{21}}=-\frac{{{d}_{22}}}{{{m}_{22}}},{{a}_{22}}=-\frac{{{d}_{23}}}{{{m}_{22}}},-{{a}_{31}}=-\frac{{{d}_{32}}}{{{m}_{33}}},{{a}_{32}}=-\frac{{{d}_{33}}}{{{m}_{33}}}, \\ & {{b}_{11}}=\frac{1}{{{m}_{11}}},{{b}_{22}}=\frac{1}{{{m}_{22}}},{{b}_{33}}=\frac{1}{{{m}_{33}}},{{c}_{11}}=-\frac{{{m}_{23}}}{{{m}_{22}}},{{c}_{22}}=-\frac{{{m}_{32}}}{{{m}_{33}}} \\ \end{align}$。

考虑系统扰动和测量噪声对船舶的影响:

$\left\{ \begin{matrix} \dot{x}=f\left( x \right)+g\left( x \right)u+\omega \\ z=h\left( x \right)+\upsilon \\ \end{matrix} \right.$

(15)

式中：z=[z₁,z₂,z₃]^T为测量输出，ω=[ω₁,ω₂,ω₃]^T为零均值系统扰动向量，υ为零均值测量噪声向量。

2 动力定位控制器设计

在外界环境中存在着各种扰动，选择适当的滤波方法，能够减少不必要的操作、节能减排和延长设备使用寿命。选用适当的滤波方法设计滤波器，是DP控制系统设计中所需解决的重要问题。基于UKF的MHE，具有无需线性化的特点，可以有效避免线性化产生的误差过大和计算困难等问题。UKF是一种非线性滤波方法，以无迹变换为基础，采用卡尔曼滤波框架和确定性采样方式进行参数估计。UKF是对非线性函数的概率密度分布进行近似，即该算法具有实现难度不会随着系统模型复杂度的增加而增加的特点^[10]。

船舶的动力学特性具有强耦合、非线性、大时滞和大惯性等特点，采用非线性模型预测控制方法，能够较好解决船舶DP控制问题。本文将基于UKF的MHE与NSAMPC相结合，用于解决随机扰动下的船舶DP非线性控制问题，应用基于UKF的MHE结合差分法所得的离散模型设计非线性滤波器；依照船舶动力学模型，运用NSAMPC方法设计控制器，将船舶当前周期状态的估计值作为下一预测周期状态的初值，用于解决随机扰动下的船舶DP非线性控制问题。NSAMPC是一种基于优化的控制策略，符合模型预测控制的基本原理，包括三个部分：预测模型、滚动优化和反馈校正。预测模型只需要其预测功能，不苛求其结构形式，为建模带来了便利；滚动优化与反馈校正相结合不但避免了全局优化带来的计算量大的问题，而且还顾及了扰动的影响能够及时加以修正，具有较好的鲁棒性。模型预测控制可以克服过程的不确定性、非线性和关联性，并能够处理被控变量和控制变量中的约束。

由控制系统、测量系统和推进系统组成的船舶DP系统结构图如图 1所示。

结合带有随机扰动的船舶动力学模型设计控制器，为了计算简便，在不影响控制结果的前提下，令m₂₃=0,m₃₂=0。

2.1 滤波器设计

运用差分法将式(15)离散化，然后应用基于UKF的MHE结合离散模型设计非线性滤波器。

对式(8)进行差分变换:

$\left\{ \begin{matrix} x\left( k+1 \right)=x\left( k \right)+\delta f\left( x\left( k \right) \right)+ \\ \delta g\left( x\left( k \right) \right)u\left( k \right)+\delta \omega \left( k \right) \\ z\left( k \right)=h\left( x\left( k \right) \right)+\upsilon \left( k \right) \\ \end{matrix} \right.$

(16)

式中：x(k)=[x₁(k),x₂(k),x₃(k),x₄(k),x₅(k),x₆(k)]^T为系统的离散状态向量；δ为差分步长；系统的控制输入为u(k)=[u₁(k),u₂(k),u₃(k)]^T；系统测量输出为z(k)=[z₁(k),z₂(k),z₃(k)]^T；ω(k)为系统扰动向量，假定均值为零；υ(k)为测量噪声向量，假定均值为零，k=0,1,2,…，为采样时刻。

$f\left( x\left( k \right) \right)=\left[ \begin{matrix} {{x}_{4}}\left( k \right)cos{{x}_{3}}\left( k \right)-{{x}_{5}}\left( k \right)sin{{x}_{3}}\left( k \right) \\ {{x}_{4}}\left( k \right)sin{{x}_{3}}\left( k \right)+{{x}_{5}}\left( k \right)cos{{x}_{3}}\left( k \right) \\ {{x}_{6}}\left( k \right) \\ {{a}_{11}}{{x}_{4}}\left( k \right) \\ {{a}_{21}}{{x}_{5}}\left( k \right)+{{a}_{22}}{{x}_{6}}\left( k \right) \\ {{a}_{31}}{{x}_{5}}\left( k \right)+{{a}_{32}}{{x}_{6}}\left( k \right) \\ \end{matrix} \right]$

(17)

$\begin{align} & g\left( x\left( k \right) \right)=[{{g}_{1}}\left( x\left( k \right) \right),{{g}_{2}}\left( x\left( k \right) \right),{{g}_{3}}\left( x\left( k \right) \right)]= \\ & {{\left[ \begin{matrix} 0 & 0 & 0 & {{b}_{11}} & 0 & 0 \\ 0 & 0 & 0 & 0 & {{b}_{22}} & 0 \\ 0 & 0 & 0 & 0 & 0 & {{b}_{33}} \\ \end{matrix} \right]}^{T}} \\ \end{align}$

(18)

$\omega ={{[{{\omega }_{1}}\left( k \right),{{\omega }_{2}}\left( k \right),{{\omega }_{3}}\left( k \right)]}^{T}}$

(19)

将式(16)规范为式(20)的形式:

$\left\{ \begin{matrix} x\left( k \right)=f\left( x\left( k-1 \right),u\left( k-1 \right),\omega \left( k-1 \right) \right) \\ z\left( k \right)=h\left( x\left( k \right) \right)+\upsilon \left( k \right) \\ \end{matrix} \right.$

(20)

针对式(20)描述的系统，滚动时域估计的目标函数如式(21)所示:

${{\Phi }_{T}}\left( \left\{ \omega \left( k \right) \right\},x\left( T-N \right) \right)=\sum\limits_{k=T-N}^{T-1}{\|\upsilon }\left( k \right){{\|}^{2}}_{{{R}^{-1}}}+\|\omega \left( k \right){{\|}^{2}}_{{{Q}^{-1}}}+{{\Theta }_{T-N}}\left( x\left( T-N \right) \right)$

(21)

满足式(22):

$\left\{ \begin{matrix} \hat{x}\left( k+1 \right)=f\left( \hat{x}\left( k \right),u\left( k \right),\omega \left( k \right) \right),k\in \left[ T-N,T-1 \right] \\ \hat{x}\left( T-N \right)=x\left( T-N \right) \\ \end{matrix} \right.$

(22)

和如下约束：

$\left\{ \begin{matrix} \hat{x}\left( k \right)\in X,k\in \left[ T-N,T \right] \\ \upsilon \left( k \right)\in V;\omega \left( k \right)\in W,k\in \left[ T-N,T-1 \right] \\ \end{matrix} \right.$

式中：Θ_T－N(x(T－N))为到达代价函数；N为估计时域长度；R为噪声向量υ(k)的协方差矩阵，Q为扰动向量ω(k)的协方差矩阵；${\hat{x}}$(T－N)为x(T－N)的估计；Χ、V和W为对应变量的约束。

选取到达代价函数如式(23)所示：

${{\Theta }_{T-N}}\left( x\left( T-N \right) \right)=\|x\left( T-N \right)-x-\left( T-N \right){{\|}^{2}}{{_{P}}^{-1}}+{{\varphi }^{*}}_{T-N}$

(23)

式中:φ_T－N^*为T-N时刻的最优到达代价；x－(T－N)为起始状态x(T－N)的先验估计；P为加权系数矩阵，表示对该先验估计的信心。

下面运用UKF方法计算到达代价。选取Sigma点，Sigma点矩阵X(k)定义为

$X\left( k \right)=\left[ m\left( k \right),\ldots ,m\left( k \right) \right]+\sqrt{c}\left[ 0,\sqrt{P\left( k \right)},-\sqrt{P\left( k \right)} \right]$

(24)

式中：m(k)为Sigma点状态均值，c=n+λ，λ=α²(n+κ)－n；α为设定参数；n为状态维度；κ为调节参数；P(k)为Sigma点状态协方差矩阵；X(k)∈R^n×d为2n+1个Sigma点的矩阵，d=2n+1。

与Sigma点对应的权重表示如下：

${{W}^{m}}_{0}=\lambda /\left( n+\lambda \right)$

(25)

${{W}^{c}}_{0}=\lambda /\left( n+\lambda \right)+(1-{{\alpha }^{2}}+\beta )$

(26)

${{W}^{m}}_{i}=1/\left[ 2\left( n+\lambda \right) \right],i=1,2,\ldots ,2n$

(27)

${{W}^{c}}_{i}=1/\left[ 2\left( n+\lambda \right) \right],i=1,2,\ldots ,2n$

(28)

将选定的Sigma点代入式(20)进行非线性变换，变换后的Sigma点集中第i个点对应的状态和测量输出分别为x_i(k)和x_i(k)。下面计算状态均值和协方差${\hat{x}}$(k|k－1)和p_x(k|k－1)，测量输出的均值和协方差${\hat{x}}$(k)和p_z(k)，组合协方差p_xz(k)。

$\begin{align} & \hat{x}\left( k|k-1 \right)={{\sum\limits_{i=0}^{2n}{{{W}^{m}}}}_{i}}{{x}_{i}}\left( k \right),\hat{z}\left( k \right)={{\sum\limits_{i=0}^{2n}{{{W}^{m}}}}_{i}}{{z}_{i}}\left( k \right){{p}_{x}}\left( k|k-1 \right)=\sum\limits_{i=0}^{2n}{{{W}^{c}}_{i}} \\ & [{{x}_{i}}\left( k \right)-\hat{x}\left( k|k-1 \right)]\times {{[{{x}_{i}}\left( k \right)-\hat{x}\left( k|k-1 \right)]}^{T}}{{p}_{z}}\left( k \right)=\sum\limits_{i=0}^{2n}{{{W}^{c}}_{i}}[{{z}_{i}}\left( k \right)-\hat{z}\left( k \right)] \\ & {{[{{z}_{i}}\left( k \right)-\hat{z}\left( k \right)]}^{T}}{{p}_{xz}}\left( k \right)=\sum\limits_{i=0}^{2n}{{{W}^{c}}_{i}}[{{x}_{i}}\left( k \right)-\hat{x}\left( k|k-1 \right)]{{[{{z}_{i}}\left( k \right)-\hat{z}\left( k \right)]}^{T}} \\ \end{align}$

(29)

通过滤波器增益K(k)计算协方差p(k)，并代入式(23)得:

$K\left( k \right)={{p}_{xz}}\left( k \right){{({{p}_{z}}\left( k \right))}^{-1}}$

(30)

$p\left( k \right)={{p}_{x}}\left( k|k-1 \right)-K\left( k \right){{p}_{z}}\left( k \right){{K}^{T}}\left( k \right)$

(31)

${{\Theta }_{k}}\left( x\left( k \right) \right)=\|x\left( k \right)-\bar{x}\left( k \right){{\|}^{2}}{{_{P(k)}}^{-1}}+{{\varphi }^{*}}_{k}$

(32)

求解目标函数(21)取得最小值时的最优解，记为(x^*(T－N),{ω^*(k)}_T－N^T－1)，则T时刻的状态估计值${\hat{x}}$(T)可表示为

$\hat{x}\left( k+1 \right)=f\left( \hat{x}\left( k \right),u\left( k \right),{{\omega }^{*}}\left( k \right) \right),k\in \left[ T-N,T-1 \right]\hat{x}\left( T-N \right)={{x}^{*}}\left( T-N \right)$

(33)

在下一时刻，将新测量值z(T)放入数据序列，去掉最旧的那个数据，保持N个测量数据不变。根据滚动优化原理，令T:=T+1，重新求解目标函数的最优解。本文采用前一时刻的估计值作为当前时刻对起始状态的先验估计:

$\bar{x}\left( T-N|T \right)=\hat{x}\left( T-N|T-1 \right)$

(34)

$\hat{x}\left( T-N|T-1 \right)=f({{x}^{*}}\left( T-N-1|T-1 \right),u\left( T-N-1|T-1 \right),{{\omega }^{*}}\left( T-N-1|T-1 \right))$

(35)

式(34)和式(35)联立得到式(36):

$\bar{x}\left( T-N|T \right)=f({{x}^{*}}\left( T-N-1|T-1 \right),u\left( T-N-1|T-1 \right),{{\omega }^{*}}\left( T-N-1|T-1 \right))$

(36)

在不引起歧义的情况下，省去“|T－1”和“T”，得到先验估计如式(37)所示:

$\bar{x}\left( T-N \right)=f({{x}^{*}}\left( T-N-1 \right),u\left( T-N-1 \right),{{\omega }^{*}}\left( T-N-1 \right))$

(37)

2.2 非切换解析模型预测控制器设计

本小节参照文献^[9]，结合系统相关度概念，应用NSAMPC方法，结合船舶动力学模型，设计船舶DP非线性控制器。

假设1 系统输出(x,y,ψ)和期望输出(x_d,y_d,ψ_d)连续，可做足够次数的微分运算。

考虑到预测控制只关心每个预测周期中控制量的初值，并且在数字控制技术及其工程实践中，控制量可近似为分段函数，因此做出假设2。

假设2 控制量在滚动时域[t,t+T_p]内为常数，即

$\hat{u}\left( t+\tau \right)=u\left( t \right),\tau \in [0,{{T}_{p}}]$

(38)

模型预测控制设计需要给出具体的性能指标。本文采用的在滚动时域内的性能指标函数:

$\begin{align} & J=\frac{1}{2}{{\mu }_{1}}{{\left[ \hat{y}\left( t+{{T}_{p}} \right)-{{{\hat{y}}}_{d}}\left( t+{{T}_{p}} \right) \right]}^{T}}[\hat{y}\left( t+{{T}_{p}} \right)-{{{\hat{y}}}_{d}}(t+{{T}_{p}})]+12{{\int }^{{{T}_{p}}}}_{0} \\ & [{{\mu }_{2}}{{(\hat{y}\left( t+\tau \right)-{{{\hat{y}}}_{d}}\left( t+\tau \right))}^{T}}\cdot (\hat{y}\left( t+\tau \right)-{{{\hat{y}}}_{d}}\left( t+\tau \right))+{{\mu }_{3}}{{{\hat{u}}}^{T}}\left( t+\tau \right)\hat{u}\left( t+\tau \right)]d\tau \\ \end{align}$

(39)

式中：T_p表示预测周期，μ₁、μ₂、μ₃是非负的(通常情况下，μ₂是正数，μ₁、μ₃可以为零)，分别反映了输出终端约束、跟踪误差和控制量所占的权重；${\hat{y}}$(t+T_p)和${\hat{y}}$_d(t+T_p)分别为系统输出和期望输出的预测值，带有上标“＾”的变量，代表预测变量，期望输出的预测值是设定值。

模型预测控制可表述在一个滚动时域[t,t+T_p]内，用${\hat{x}}$(t+τ)表示系统状态，用${\hat{u}}$(t+τ)表示控制输入，那么系统在滚动时域内的动态方程可以表示为

$\left\{ \begin{matrix} \hat{x}\prime \left( t+\tau \right)=f\left( \hat{x}\left( t+\tau \right) \right)+g\left( \hat{x}\left( t+\tau \right) \right)\hat{u}\left( t+\tau \right) \\ \hat{y}\left( t+\tau \right)=h\left( \hat{x}\left( t+\tau \right) \right),\tau \in [0,{{T}_{p}}] \\ \end{matrix} \right.$

(40)

其初始状态为系统的当前状态，即

$\hat{x}\left( t+\tau \right){{|}_{\tau =0}}=x\left( t \right)$

(41)

根据式(40)和(41)的约束，可以预测[t,t+T_p]时间段内系统的输出。最小化性能指标J的问题就是寻找[t,t+T_p]时间段内的最优控制输入${\hat{u}}$^*(t+τ)。则非线性模型预测控制可以描述为

$\underset{{\hat{u}}}{\mathop{min}}\,J$

(42)

式(42)可以无需计算τ∈[0,T_p]内所有的控制输入，而只计算其初始值。与其他预测控制方法类似，实际控制量u(t)取最优控制输入${\hat{u}}$^*(t+τ)的初始值，即

$u\left( t \right)={{\hat{u}}^{*}}\left( t+\tau \right){{|}_{\tau =0}}$

(43)

在NSAMPC中，系统控制输入总是取使性能指标J最小的控制输入，并且只关心其初始值。当滚动时域时，性能指标J逐渐减小，同时系统输出逐渐接近期望值。

为了得到预测控制律，将滚动时域[t,t+T_p]内系统的输出以及性能指标进行适当阶次的泰勒级数展开。

系统输出${\hat{y}}$(t+τ)以及期望输出$\hat{y}$_d(t+τ)可以表示为

$\left\{ \begin{matrix} \hat{y}\left( t+\tau \right)={{\tau }^{T}}\left( \tau \right)\hat{J}\left( t \right) \\ {{{\hat{y}}}_{d}}\left( t+\tau \right)={{\tau }^{T}}\left( \tau \right){{{\hat{y}}}_{d}}\left( t \right) \\ \end{matrix} \right.$

(44)

式中：

$\tau \left( \tau \right)=\left[ \begin{matrix} {{\tau }_{1}} & {} & {} \\ {} & {{\tau }_{2}} & {} \\ {} & {} & {{\tau }_{3}} \\ \end{matrix} \right]$

(45)

${{\tau }_{1}}={{\tau }_{2}}={{\tau }_{3}}={{[1\tau \ldots \frac{{{\tau }^{n}}}{N!}]}^{T}}$

(46)

$\begin{align} & \hat{y}\left( t \right)=[{{{\hat{y}}}_{1}}\left( t \right){{{\hat{y}}}_{1}}\prime \left( t \right)\ldots {{{\hat{y}}}_{1}}^{(N)}\left( t \right){{{\hat{y}}}_{2}}\left( t \right){{{\hat{y}}}_{2}}\prime \left( t \right)\ldots {{{\hat{y}}}_{2}}^{(N)}\left( t \right){{{\hat{y}}}_{3}}\left( t \right){{{\hat{y}}}_{3}}\prime \left( t \right)\ldots \\ & {{{\hat{y}}}_{3}}(N)\left( t \right){{]}^{T}}{{{\hat{y}}}_{d}}\left( t \right)=[{{{\hat{y}}}_{1d}}\left( t \right){{{\hat{y}}}_{1d}}\prime \left( t \right)\ldots {{{\hat{y}}}_{1d}}^{(N)}\left( t \right){{{\hat{y}}}_{2d}}\left( t \right){{{\hat{y}}}_{2d}}\prime \left( t \right)\ldots {{{\hat{y}}}_{2d}}^{(N)} \\ & \left( t \right){{{\hat{y}}}_{3d}}n\left( t \right){{{\hat{y}}}_{3d}}\prime \left( t \right)\ldots {{{\hat{y}}}_{3d}}^{(N)}\left( t \right){{]}^{T}}n \\ \end{align}$

(47)

其中，${\hat{y}}$_j⁽ⁱ⁾(t)和${\hat{y}}$⁽ⁱ⁾_jd(t),(i=0,1,…,N;j=1,2,3)分别代表$\hat{y}$_j(t)和$\hat{y}$_jd(t)的i阶时间导数；N是泰勒级数展开阶次；期望输出$\hat{y}$_d(t+τ)=y_d(t+τ),$hat{y}$_d(t)=y_d(t)。

在假设2下，${\hat{u}}$的各阶导数均为零。因此式(39)所示的性能指标J的N阶泰勒级数展开可以近似为

$J\approx \frac{1}{2}{{[\hat{y}\left( t \right)-{{\hat{y}}_{d}}\left( t \right)]}^{T}}M[\hat{y}\left( t \right)-{{\hat{y}}_{d}}\left( t \right)]+\frac{1}{2}{{\mu }_{3}}{{T}_{p}}{{\hat{u}}^{T}}\hat{u}$

(48)

式中：

$M={{\mu }_{1}}\tau \left( \tau \right){{\tau }^{T}}\left( \tau \right){{|}_{\tau ={{T}_{p}}}}+{{\mu }_{2}}{{\int }^{{{T}_{p}}}}_{0}\tau \left( \tau \right){{\tau }^{T}}\left( \tau \right)d\tau $

(49)

因此，式(42)等价于：

${{\left( \frac{\partial \hat{y}\left( t \right)}{\partial \hat{u}} \right)}^{T}}M[\hat{y}\left( t \right)-{{\hat{y}}_{d}}\left( t \right)]+{{\mu }_{3}}{{T}_{p}}\hat{u}=0$

(50)

${{\left( \frac{\partial \hat{y}\left( t \right)}{\partial \hat{u}} \right)}^{T}}M\left( \frac{\partial \hat{y}\left( t \right)}{\partial \hat{u}} \right)+{{(\frac{{{\partial }^{2}}\hat{y}\left( t \right)}{\partial {{{\hat{u}}}^{2}}})}^{T}}\cdot M[\hat{y}\left( t \right)-{{\hat{y}}_{d}}\left( t \right)]+{{\mu }_{3}}{{T}_{p}}>0$

(51)

应用泰勒级数展开来分析$\hat{y}$_j(t)(j=1,2,3)的各阶导数，由假设2可知${\hat{u}}$的导数为零，即${\hat{u}}$⁽ⁱ⁾=0(i≥1)，从而可以去掉$\hat{y}$_j(t)的导数中的${\hat{u}}$的导数项，限于篇幅只给出与本文相关的内容：

$\left\{ \begin{matrix} {{{\hat{y}}}_{1}}\left( t \right)={{h}_{1}}=q{{1}_{0,0}},{{{\hat{y}}}_{2}}\left( t \right)={{h}_{2}}=q{{2}_{0,0}} \\ {{{\hat{y}}}_{3}}\left( t \right)={{h}_{3}}=q{{3}_{0,0}} \\ \end{matrix} \right.$

(52)

$\left\{ \begin{matrix} \hat{y}{{\prime }_{1}}\left( t \right)={{L}_{f}}{{h}_{1}}=q{{1}_{1,0}},\hat{y}{{\prime }_{2}}\left( t \right)={{L}_{f}}{{h}_{2}}=q{{2}_{1,0}}, \\ \hat{y}{{\prime }_{3}}\left( t \right)={{L}_{f}}{{h}_{3}}=q{{3}_{1,0}} \\ \end{matrix} \right.$

(53)

$\begin{align} & {{{\hat{y}}}_{1}}^{(2)}\left( t \right)={{L}^{2}}_{f}{{h}_{1}}+[{{L}_{g}}_{1}{{L}_{f}}{{h}_{1}},{{L}_{g}}_{2}{{L}_{f}}{{h}_{1}},{{L}_{g}}_{3}{{L}_{f}}{{h}_{1}}]\hat{u}=q{{1}_{2,0}}+q{{1}_{2,1}}\hat{u}, \\ & {{{\hat{y}}}_{2}}^{(2)}\left( t \right)={{L}^{2}}_{f}{{h}_{2}}+[{{L}_{g}}_{1}{{L}_{f}}{{h}_{2}},{{L}_{g}}_{2}{{L}_{f}}{{h}_{2}},{{L}_{g}}_{3}{{L}_{f}}{{h}_{2}}]\hat{u}=q{{2}_{2,0}}+q{{2}_{2,1}}\hat{u}, \\ & {{{\hat{y}}}_{3}}^{(2)}\left( t \right)={{L}^{2}}_{f}{{h}_{3}}+[{{L}_{g}}_{1}{{L}_{f}}{{h}_{3}},{{L}_{g}}_{2}{{L}_{f}}{{h}_{3}},{{L}_{g}}_{3}{{L}_{f}}{{h}_{3}}]\hat{u}=q{{3}_{2,0}}+q{{3}_{2,1}}\hat{u} \\ \end{align}$

(54)

$\begin{align} & {{{\hat{y}}}_{1}}^{(3)}\left( t \right)=q{{1}_{3,0}}+q{{1}_{3,1}}\hat{u}+q{{1}_{3,2}}{{{\hat{u}}}^{2}}{{{\hat{y}}}_{2}}^{(3)}\left( t \right)=q{{2}_{3,0}}+q{{2}_{3,1}}\hat{u}+q{{2}_{3,2}}{{{\hat{u}}}^{2}}{{{\hat{y}}}_{3}}^{(3)} \\ & \left( t \right)=q{{3}_{3,0}}+q{{3}_{3,1}}\hat{u}+q{{3}_{3,2}}{{{\hat{u}}}^{2}}\ldots \ldots \\ \end{align}$

(55)

式中：$\hat{u}=\left[ \begin{matrix} {{{\hat{u}}}_{1}} \\ {{{\hat{u}}}_{2}} \\ {{{\hat{u}}}_{3}} \\ \end{matrix} \right],{{\hat{u}}^{2}}={{\left[ {{{\hat{u}}}^{2}}_{1},{{{\hat{u}}}_{1}}{{{\hat{u}}}_{2}},{{{\hat{u}}}^{2}}_{2},{{{\hat{u}}}_{2}}{{{\hat{u}}}_{3}},{{{\hat{u}}}^{2}}_{3},{{{\hat{u}}}_{1}}{{{\hat{u}}}_{3}} \right]}^{T}}$。

讨论一下确定和不确定相关度的问题。从式(52)~(55)可以导出对于所有状态x满足L_g3L_fh₃=b₃₃≠0，因此ρ₃₃是确定相关度。ρ₁₁,ρ₁₂,ρ₁₃,ρ₂₁,ρ₂₂,ρ₂₃,ρ₃₁,ρ₃₂在某些状态下为零，则ρ₁₁,ρ₁₂,ρ₁₃,ρ₂₁,ρ₂₂,ρ₂₃,ρ₃₁,ρ₃₂是不确定相关度。其中，ρ_ij(1≤i≤l,1≤j≤m)表示控制输入u_i到系统输出y_j的相关度。

本文选取N=3，L=2。将$\hat{y}$_j(t)各阶导数中${\hat{u}}$ⁱ项的系数表示为

${{q}_{.,i}}\left( x \right)={{[q{{1}_{.,i}}\left( x \right),q{{2}_{.,i}}\left( x \right),q{{3}_{.,i}}\left( x \right)]}^{T}}$

(56)

式中：$\begin{align} & q{{k}_{.,i}}\left( x \right)=\left[ \begin{matrix} q{{k}_{0,i}}\left( x \right) \\ q{{k}_{1,i}}\left( x \right) \\ \vdots \\ q{{k}_{N,i}}\left( x \right) \\ \end{matrix} \right],k=1,2,3;i=0,1\text{ }\ldots ,L; \\ & \left\{ \begin{matrix} q{{k}_{0,i}}\left( x \right)=0,0＜i\le L \\ q{{k}_{j,i}}\left( x \right)=0,1\le j\le i \\ \end{matrix} \right.,j=1,2,,\ldots ,N;i=0,1,\ldots ,L \\ \end{align}$。

因此向量${\hat{y}}$(t)重写为

$\hat{y}\left( t \right)=[{{q}_{.,0}}\left( x \right),{{q}_{.,1}}\left( x \right),\ldots ,{{q}_{.,L}}\left( x \right)]\begin{matrix} 1 \\ {\hat{u}} \\ \vdots \\ {{{\hat{u}}}^{L}} \\ \end{matrix}=Q\left( x \right)\hat{U}$

(57)

简单起见，忽略${\hat{u}}$的二次方项及更高次项，即

$\hat{U}\approx {{\left[ 1\hat{u}0 \right]}^{T}}$

(58)

将式(57)和(58)代入式(50)，得到：

${{q}_{.,1}}{{\left( x \right)}^{T}}M\left( {{q}_{.,0}}\left( x \right)-{{y}_{d}} \right)+({{q}_{.,1}}{{\left( x \right)}^{T}}M{{q}_{.,1}}\left( x \right)+{{M}_{3}})\hat{u}=0$

(59)

式中M₃=diag(μ₃T_p,μ₃T_p,μ₃T_p)。

由式(59)，得到NSAMPC的解：

$\hat{u}\left( t+\tau \right)=-\frac{{{q}_{.,1}}{{\left( x \right)}^{T}}M({{q}_{.,0}}\left( x \right)-{{y}_{d}})}{{{q}_{.,1}}{{\left( x \right)}^{T}}M{{q}_{.,1}}\left( x \right)+{{M}_{3}}},\tau \in [0,{{T}_{p}}]$

(60)

取其初值，即为NSAMPC的船舶DP非线性控制律：

$\hat{u}\left( t \right)=-\frac{{{q}_{.,1}}{{\left( x \right)}^{T}}M({{q}_{.,0}}\left( x \right)-{{y}_{d}})}{{{q}_{.,1}}{{\left( x \right)}^{T}}M{{q}_{.,1}}\left( x \right)+{{M}_{3}}}$

(61)

鉴于q_.,1(x)^TMq_.,1(x)+M₃是正定的，具体证明参见文献^[8]，控制律(61)解决了由奇异点带来的不确定相关度问题；而且该控制律对所有状态都是连续的，避免了控制器在不同状态间切换产生的系统振荡。

3 供应船动力定位仿真

本文以某DP供应船为对象进行仿真，验证文中提出的控制方法。主要参数：船长76.2 m；船舶净重4 200 t；主发动机功率3 533 kW；纵向推力[-1 000 kN,1 000 kN]；横向推力[-300 kN,300 kN]；艏摇力矩[-7 620 kN·m,7 620 kN·m]；无因次量m₁₁=1.127 4，m₂₂=1.890 2，m₃₃=0.127 8，d₁₁=0.035 8，d₂₂=0.118 3，d₂₃=－0.012 4，d₃₂=－0.004 1，d₃₃=0.030 8。差分步长δ=0.1 s，估计时域长度N=10，系统扰动的均值为0，协方差为diag(10^－2,10^－2,10^－2,10^－2,10^－2,10^－2)，测量噪声的均值为0,协方差为diag(1,1,10^－1)；α=0.5，β=2，κ=－2；预测周期T_p取1.4 s，μ₁=1.0，μ₂=0.005，μ₃=0，仿真时长取500 s；船舶的起始位置为(0 m,0 m,0°)，期望位置为(50 m,50 m,0°)。

根据设定的参数进行仿真，系统稳态输出测量值未滤波和滤波后的均值和方差如表 1所示。

由表 1可以看出，在随机扰动下系统稳态输出测量值的均值与设定值相差不大，但存在一定的方差；系统稳态输出测量值滤波后与原值相比，均值基本相当，方差明显减小。

采用本文设计控制器的DP系统仿真结果如图 2和图 3所示。

从图 2(a)和图 2(b)可以看出，船舶在100 s内到达指定位置，北向位置和东向位置的超调量较小，调节时间较短，并且曲线较为平滑；从图 2(c)可以看出船舶艏向角度超调量不大，但是所需调节时间较长；从图 3可以看出，控制力和力矩满足约束，其曲线连续，并且波动较小。仿真结果表明，本文设计的船舶DP控制器，能够较好地解决在随机扰动下的船舶DP控制问题，使船舶较快到达并保持在指定位置。

4 结论

本文提出了一种船舶DP控制器设计方法，应用基于UKF的MHE方法获取船舶运动状态的估计值，将该估计值作为NSAMPC控制的反馈值，通过滚动优化得到船舶运动的控制力和力矩。仿真验证了本文所设计的控制器的有效性，并得出如下结论：

1) 应用基于UKF的MHE方法对船舶运动状态进行估计，能够使滤波后信号的方差比原信号的方差明显减小。

2)将基于UKF的MHE与NSAMPC相结合设计的船舶DP控制器的控制力和力矩曲线是连续的，并且波动较小，能够防止由控制器的控制力和力矩大幅波动导致的推进器的机械磨损。

下一步的工作将考虑船舶动力定位系统中存在的参数不确定问题，对船舶动力定位鲁棒控制进行研究。