0 引　言

高速多体船具有横向稳定性好、装载量大、操控性强等优点，发展前景广阔^[1]。然而在恶劣海况下，因自身细长的侧体使得纵向恢复力矩较小，高速多体船难以抵消海浪产生的干扰力矩，造成高速多体船垂向稳定性较差，升沉和纵摇运动幅度过大，引起首部砰击、晕船等现象，严重影响适航性^[2]。在船体上安装T型翼和压浪板两类附体可实现协同减摇，为进一步发挥附体的减摇能力，需要设计合理的主动控制策略自动调节附体的攻角，这是目前多体船领域的主要研究内容之一^[3]。

高速多体船的减摇控制存在一定挑战。基于惯性测量的主动式减摇控制中需要引入升沉速度和纵摇角速度，提高减摇效果，但是惯导平台不能直接测量这2个量；同时，高速多体船的升沉和纵摇运动存在强耦合，获得的水动力学系数具有较强不确定性^[4]，目前国内外这方面的研究较少。在升沉速度和纵摇角速度估计方面，赵希人等^[5]将有色海浪扰动白噪化，把白噪声作为系统输入，建立扩张状态的多体船动力学模型，采用卡尔曼滤波器估计状态。但是，该方法需要辨识成型滤波器的参数，计算复杂。在减摇控制方面，CRUZ等^[6]提出工程易于实现的比例微分减摇控制律，但是该方法需花费大量时间离线调试参数，不能适应海况的变化。为提高鲁棒性，ARANDA等^[7]在考虑多体船的水动力学系数最大变化范围下，在频域内设计高阶的定量反馈鲁棒控制器，但控制性能比较保守。IBRAHIM等^[8]考虑优化减摇性能进而提出预测控制减摇方法，降低了鲁棒控制的保守性，减摇效果明显得到提高。但是，预测控制采用序列二次规划进行求解控制序列，在线计算量大，难以实时应用。模型预静态规划基于终端偏差修正的求解方式，在实时性方面有着明显优势^[9]，这为多体船的预测控制减摇提供了可行的解决途径。

考虑多体船的信号估计和减摇控制问题，设计平滑变结构滤波器来估计升沉速度、纵摇角速度，对于多体船的水动力学系数不确定性、以及有色海浪噪声，具有较好的鲁棒性；基于模型预测静态规划提出预测控制减摇，采用迭代更新方式求解控制输入，降低预测控制在线计算量，提高减摇效果。

1 高速多体船垂向运动模型

主动式T型翼和压浪板分别在船首底部和船尾，如图1所示。

基于刚体动量定理和动量矩定理建立高速多体船的非线性数学模型为^[4]

$ \left\{ \begin{array}{l} m(\dot{w}-qu) = 2F_{\text{T - foil}} + F_{\text{flap}} + F_{H} + F_{b}\cos\theta +\\ \qquad\qquad\quad mg\cos\theta+F_{\text{wave}}，\\ I\dot{q}=-2(x_{\text{T - foil}}-x_{\text{G}})F_{\text{T - foil}}-(x_{\text{flap}}-x_{G})\times F_{\text{flap}}-\\ \qquad (x_{H}-x_{G})F_{H}-(x_{b}-x_{G})F_{b}\cos\theta+M_{\text{wave}}。\end{array} \right. $

(1)

式中：$ m $为船质量；$ \theta $为俯仰角；$ I $为船体关于y轴的转动惯量；$ \dot w $为航行$ z $方向的加速度；$ q $为纵摇角速度；$ u $为航速；$ {F_{{\text{T - foil}}}} $和$ {F_{{\text{flap}}}} $分别为减摇附体T型翼和压浪板产生的作用力；$ {F_{H}} $为船体产生的升力；$ {F_{b}} $为船的浮力；$ {x_{G}} $、$ {x_{H}} $、$ {x_{b}} $分别为多体船的重心、升力作用点、浮力作用点到船中心的距离；$ {x_{{\text{T - foil}}}} $、$ {x_{{\text{flap}}}} $分别为T型翼和压浪板作用点到船体中心的距离；${F_{{\text{wave}}}}$和$ {M_{{\text{wave}}}} $分别为随机海浪产生的扰动力和力矩；${F_{{\text{T - foil}}}}$、${M_{{\text{T - foil}}}}$和$ {F_{{\text{flap}}}} $、$ {M_{{\text{flap}}}} $分别为T型翼和压浪板所产生的恢复力和力矩。

高速多体船减摇控制目标是在随机海浪扰动下，自动调节T型翼和压浪板的攻角，控制升沉和纵摇运动幅度尽可能小。

2 平滑变结构滤波器

在基于惯性测量的主动式减摇控制中，无法直接测量升沉速度和纵摇角速度，需要对其进行估计。高速多体船垂向运动模型存在升沉和纵摇运动强耦合、水动力学系数不确定性特点，并且海浪干扰具有非高斯分布特性，因此常规的扩展卡尔曼滤波精度下降。平滑变结构滤波器基于滑模控制的变结构增益策略，通过系统的误差和噪声信息设定状态估计的边界，可提高对模型不确定和非高斯噪声的鲁棒性。本节采用平滑变结构滤波器估计状态^[10]，步骤如下：

对高速多体船的垂向运动模型式（1）在平衡点处$ ({{\boldsymbol{\dot x}}_0},{{\boldsymbol{x}}_0},{{\boldsymbol{u}}_0}) $进行线性化得到

$\left\{ \begin{split} & \dot {\boldsymbol x} = {\boldsymbol{Ax}} + {\boldsymbol{Bu}} + {\boldsymbol{B}_w}{\boldsymbol d}，\\ & {\boldsymbol y} = {\boldsymbol{Cx}} 。\end{split} \right. $

(2)

式中：$ {\boldsymbol x}=[\dot x_h, \dot x_p, x_h, x_p] $分别为升沉速度、纵摇角速度、升沉位移和纵摇角，设计平滑变结构滤波器的滤波增益为：

$ {{\boldsymbol{K}}_k} = {{\boldsymbol{C}}^ + }{\text{diag}}\left[ {\left( {\left| {{{\boldsymbol{e}}_{k|k - 1}}} \right| + \gamma \left| {{{\boldsymbol{e}}_{k - 1}}} \right|} \right) \circ S} \right]{\text{diag}}{\left( {{{\boldsymbol{e}}_{k|k - 1}}} \right)^{ - 1}}，$

(3)

$ \begin{split} & S=\left[\begin{array}{c}\text{sat}\left({e}_{k|k-1,1}\cdot {\psi }_{11}^{-1}\right)\\ \text{sat}\left({e}_{k|k-1,2}\cdot {\psi }_{22}^{-1}\right)\\ \vdots\\ \text{sat}\left({e}_{k|k-1,{n}_{e}}\cdot {\psi }_{{n}_{e}{n}_{e}}^{-1}\right)\end{array}\right],\\ & \text{sat}\left({e}_{k|k-1,i}\cdot {\psi }_{ii}^{-1}\right)=\left\{\begin{array}{ll} 1,& {e}_{k|k-1,i}\cdot {\psi }_{ii}^{-1}\geqslant 1，\\ {e}_{k|k-1,i}\cdot {\psi }_{ii}^{-1}, &1 > {e}_{k|k-1,i}\cdot {\psi }_{ii}^{-1} > -1\\ -1,& {e}_{k|k-1,i}\cdot {\psi }_{ii}^{-1}\leqslant -1。\end{array} \right. \end{split} $

$ {{\boldsymbol{\hat y}}_{k|k - 1}} = {\boldsymbol{C}}{{\boldsymbol{\hat x}}_{k|k - 1}} \text{，} {{\boldsymbol{\hat y}}_{k|k - 1}} = {\boldsymbol{C}}{{\boldsymbol{\hat x}}_{k|k - 1}} \text{，} {{\boldsymbol{e}}_{k - 1}} = {{\boldsymbol{y}}_{k - 1}} - {{\boldsymbol{\hat y}}_{k - 1}}。$

式中：$ \gamma $为SVSF的收敛速率；$ {e_{k|k - 1,i}} $和$ {\psi _{ii}} $分别为先验测量误差$ {{\boldsymbol{e}}_{k|k - 1}} $的第i个元素和平滑边界层宽度矩阵$ {\boldsymbol{\bar \psi }} $的第i个对角线元素，$ {n_{e}} $为$ {{\boldsymbol{e}}_{k|k - 1}} $的维数；diag(v) 为以向量 v的各元素作为主对角线元素，构造对角矩阵；“$ \circ $”为 Hadamard 乘积。

平滑边界层宽度矩阵$ {\boldsymbol{\bar \psi }} $表示SVSF状态估计的不确定性程度，增大$ {\boldsymbol{\bar \psi }} $可以减小估计状态的抖振，状态更为光滑，但$ {\boldsymbol{\bar \psi }} $过大会导致状态估计精度的下降。根据噪声信息和误差信息，提出自适应平滑边界层宽度矩阵$ {{\boldsymbol{\bar \psi }}_k} $如下:

$ \begin{split} {{\boldsymbol{\bar \psi }}_k} = \,& \left\{ {{\left[ {{\text{diag}}\left( {\left| {{{\boldsymbol{e}}_{k|k - 1}}} \right| + \gamma \left| {{{\boldsymbol{e}}_{k - 1}}} \right|} \right)} \right]}^{ - 1}}{\boldsymbol{C}}{{\boldsymbol{P}}_{k|k - 1}}{{\boldsymbol{C}}^{\text{T}}}\right.\cdot\\ & \left.{{\left( {{\boldsymbol{C}}{{\boldsymbol{P}}_{k|k - 1}}{{\boldsymbol{C}}^{\text{T}}} + {{\boldsymbol{R}}_{{v}}}} \right)}^{ - 1}} \right\}^{ - 1} 。\end{split} $

(4)

式中，$ {{\boldsymbol{R}}_{v}} $为量测噪声协方差矩阵，则k时刻的状态估计值$ {{\boldsymbol{\hat x}}_k} $为：

$ {{\boldsymbol{\hat x}}_k} = {{\boldsymbol{\hat x}}_{k|k - 1}} + {{\boldsymbol{K}}_k}{{\boldsymbol{e}}_{k|k - 1}} 。$

(5)

协方差矩阵$ {{\boldsymbol{P}}_k} $为：

$ {{\boldsymbol{P}}_k} = \left[ {{\boldsymbol{I}} - {{\boldsymbol{K}}_k}{\boldsymbol{C}}} \right]{{\boldsymbol{P}}_{k|k - 1}}{\left[ {{\boldsymbol{I}} - {{\boldsymbol{K}}_k}{\boldsymbol{C}}} \right]^{\text{T}}} + {{\boldsymbol{K}}_k}{\boldsymbol{RK}}_k^{\text{T}} 。$

(6)

3 基于模型静态规划的减纵摇控制

在获得升沉位移、升沉速度、纵摇角度、纵摇角速度后，需要设计减摇控制器。减摇控制目的是保持升沉和纵摇运动幅度尽可能小，具有较强的抑制海浪扰动能力，因此对控制系统的鲁棒性和控制性能提出了较高要求。预测控制具有多目标优化特点，并且其滚动优化和模型预测策略能有效抑制模型不确定性及外界扰动，提出具有终端不等式约束的有限时域预测控制减纵摇方法，即下面的优化问题：

$ \begin{gathered} {\text{ }}\mathop {\min }\limits_u J(k) \\ {\mathrm{s}}.{\mathrm{t}}.\left\{ \begin{gathered} \hat x(k + j + 1) = {A_d}\hat x(k + j) + {B_d}u(k + j),{\text{ }}j = 1,{\text{ }}2,{\text{ }}...,N，\\ \hat x{(k + N)^{\mathrm{T}}}P\hat x(k + N) \leqslant {\sigma _{T}}，\\ \end{gathered} \right. \\ \end{gathered} $

(7)

$ \begin{split} J(k) = \;& \frac{1}{2}\sum\limits_{j = 1}^N {{{(y(k + j) - {y^*})}^{\text{T}}}{\boldsymbol{Q}}} (y(k + j) - {y^*}) + \\ & \;\frac{1}{2}\sum\limits_{j = 1}^N {{{(u(k + j))}^{\text{T}}}{\boldsymbol{R}}(u(k + j))}。\end{split} $

式中：$ {\boldsymbol{Q}} \gt 0 $和$ {\boldsymbol{R}} \geqslant 0 $分别为状态和控制加权矩阵；$ N $为控制时域和预测时域；$ {\sigma _{T}} $为设计的终端约束界限常数，保证有限时域优化的闭环稳定性。

传统的预测控制采用序列二次规划求解控制序列^[8]，在线计算量大，难以应用到高速多体船上。为了减少式 （7）在线计算量，采用模型预测静态规划求解，该方法以一组猜想的预测控制序列开始，基于误差修正和最优控制问题得到控制输入的最优迭代增量，直至满足终端状态满足约束^[11]。将多体船标称模型作为预测模型，引入迭代步数$ i $，得到迭代预测模型

$ \left\{ \begin{gathered} {\boldsymbol{x}}_{k+j|k}^{i+1} = {\boldsymbol{x}}_{k+j|k}^i + {\text{d}}{\boldsymbol{x}}_{k+j|k}^i，\\ {\boldsymbol{y}}_{k+j|k}^{i+1} = {\boldsymbol{y}}_{k+j|k}^i + {\text{d}}{\boldsymbol{y}}_{k+j|k}^i，\\ {\boldsymbol{u}}_{k+j|k}^{i+1} = {\boldsymbol{u}}_{k+j|k}^i + {\text{d}}{\boldsymbol{u}}_{k+j|k}^i。\\ \end{gathered} \right. $

(8)

式中：$ {\boldsymbol{x}}_{k|k}^i = {{\boldsymbol{\hat x}}_{{t_k}}} $为k时刻的估计状态；$ j = 0,1,...,N - 1 $；$ {\boldsymbol{x}}_{k + j|k}^i $、$ {\boldsymbol{y}}_{k + j|k}^i $、$ {\boldsymbol{u}}_{k + j|k}^i $分别为在第i次迭代中对k+j时刻的预测状态、输出、输入；$ {\text{d}}{\boldsymbol{x}}_{k + j|k}^i $、$ {\text{d}}{\boldsymbol{y}}_{k + j|k}^i $、$ {\text{d}}{\boldsymbol{u}}_{k + j|k}^i $分别为第i次迭代向第i+1次迭代更新时的预测状态、输出、输入迭代增量。由于$ {\boldsymbol{x}}_{k|k}^i $是$ {t_k} $时刻的状态估计值，不会随迭代变化，所以$ {\text{d}}{\boldsymbol{x}}_{k|k}^i = 0 $，可推得

$ {\text{d}}{\boldsymbol{y}}_{k + j|k}^i{\text{ = }}\sum\limits_{l = 1}^j {{\boldsymbol{CA}}_{d}^{j - l}{{\boldsymbol{B}}_{d}}{\text{d}}{\boldsymbol{u}}_{k + l - 1|k}^i}。$

(9)

第i次迭代的优化目标函数为：

$ \begin{split} {J^i} =\,& \frac{1}{2}\sum\limits_{j = 1}^N {{{\left({\boldsymbol{y}}_{k + j|k}^i + \sum\limits_{l = 1}^j {{\boldsymbol{CA}}_{d}^{j - l}{{\boldsymbol{B}}_{d}}{\text{d}}{\boldsymbol{u}}_{k + l - 1|k}^i} \right)}^{\text{T}}}{\boldsymbol{Q}}}\times\\ & \left({\boldsymbol{y}}_{k + j|k}^i + \sum\limits_{j = 1}^k {{\boldsymbol{CA}}_{d}^{j - l}{{\boldsymbol{B}}_{d}}{\text{d}}{\boldsymbol{u}}_{k + l - 1|k}^i} \right)+\\ &\frac{1}{2}\sum\limits_{j = 1}^N {{{({\text{d}}{\boldsymbol{u}}_{k + j - 1|k}^i)}^{\text{T}}}{\boldsymbol{R}}({\text{d}}{\boldsymbol{u}}_{k + j - 1|k}^i)} 。\end{split} $

(10)

根据最优控制的必要条件，有

$ \frac{{{\text{d}}{J^i}}}{{{\text{d}}({\text{d}}{\boldsymbol{u}}_{k|k,N}^i)}} = 0 。$

(11)

得到第i次迭代增量

$ {\text{d}}{\boldsymbol{u}}_{k|k,N}^{i*} = - {({{\boldsymbol{M}}^{\text{T}}}{{\boldsymbol{Q}}_N}{\boldsymbol{M}} + {{\boldsymbol{R}}_N})^{ - 1}}{{\boldsymbol{M}}^{\text{T}}}{{\boldsymbol{Q}}_N}{\boldsymbol{y}}_{k + 1|k,N}^i 。$

(12)

式中：

$ {\text{d}}{\boldsymbol{u}}_{k|k,N}^i = \left[ {\begin{array}{*{20}{c}} {{\text{d}}{\boldsymbol{u}}_{k|k}^i} \\ {{\text{d}}{\boldsymbol{u}}_{k + 1|k}^i} \\ \vdots \\ {{\text{d}}{\boldsymbol{u}}_{k + N - 1|k}^i} \end{array}} \right] {\boldsymbol{y}}_{k + 1|k,N}^i = \left[ {\begin{array}{*{20}{c}} {{\boldsymbol{y}}_{k + 1|k}^i} \\ {{\boldsymbol{y}}_{k + 2|k}^i} \\ \vdots \\ {{\boldsymbol{y}}_{k + N|k}^i} \end{array}} \right]，$

$ {\boldsymbol{M}} = \left[ {\begin{array}{*{20}{c}} {{\boldsymbol{C}}{{\boldsymbol{B}}_{\text{d}}}}&{\boldsymbol{0}}& \cdots &{\boldsymbol{0}} \\ {{\boldsymbol{C}}{{\boldsymbol{A}}_{\text{d}}}{{\boldsymbol{B}}_{\text{d}}}}&{{\boldsymbol{C}}{{\boldsymbol{B}}_{\text{d}}}}& \cdots &{\boldsymbol{0}} \\ \vdots & \vdots & \ddots & \vdots \\ {{\boldsymbol{CA}}_{\text{d}}^{N - 1}{{\boldsymbol{B}}_{\text{d}}}}&{{\boldsymbol{CA}}_{\text{d}}^{N - 2}{{\boldsymbol{B}}_{\text{d}}}}& \cdots &{{\boldsymbol{C}}{{\boldsymbol{B}}_{\text{d}}}} \end{array}} \right]，$

$ {{\boldsymbol{Q}}_N} = {\text{diag}}\left[ {{\boldsymbol{Q}},{\boldsymbol{Q}}, \ldots ,{\boldsymbol{Q}}} \right] {{\boldsymbol{R}}_N} = {\text{diag}}\left[ {{\boldsymbol{R}},{\boldsymbol{R}}, \ldots ,{\boldsymbol{R}}} \right]。$

第i+1次迭代的最优控制输入序列：

$ {\boldsymbol{u}}_{k|k,N}^{i + 1}{\text{ = }}{\boldsymbol{u}}_{k|k,N}^i + {\text{d}}{\boldsymbol{u}}_{k|k,N}^{i*}，$

(13)

式中，$ {\boldsymbol{u}}_{k|k,N}^i = {\left[ {\begin{array}{*{20}{c}} {{\boldsymbol{u}}_{k|k,N}^{i{\text{T}}}}&{{\boldsymbol{u}}_{k|k,N}^{i{\text{T}}}}& \cdots &{{\boldsymbol{u}}_{k|k,N}^{i{\text{T}}}} \end{array}} \right]^{\text{T}}} $。根据第i+1次迭代的控制序列和预测模型获得终端状态$ {\boldsymbol{x}}_{k + N|k}^{i + 1} $，当终端约束项满足下式时，停止迭代。

$ {\boldsymbol{x}}_{k + N|k}^{{l_p}{\text{T}}}{\boldsymbol{Qx}}_{k + N|k}^{{l_p}} \leqslant {\sigma _{T}} 。$

(14)

最终，得到最优控制输入序列如下：

$ {\boldsymbol{u}}_{k|k,N}^* = {\boldsymbol{u}}_{k|k,N}^{{l_p}}。$

(15)

模型预测静态规划算法在计算中多次使用了小量近似，需要若干次迭代才能满足跟踪要求，一般在5次以内能收敛。

4 仿真研究

以某型三体船为例进行减摇控制仿真^[2]，航速为40 kn，遭遇频率为$ 3.25{\text{ rad/s}} $。将4级海况下的海浪扰动力和力矩作用于三体船上，进行下面仿真研究：

1）平滑变结构滤波器状态估计：量测噪声方差按一级精度的传感器，取R_v=diag[20.3×10⁻⁴ 2.26×10⁻⁶]，根据海浪扰动样本系统噪声方差取Q_w=diag[0.0312 1.1084]×10⁷。采用平滑变结构滤波器的状态估计见图2～图3。可以看出，平滑变结构滤波器能准确估计升沉速度、纵摇角速度，克服了传统的一阶差分加低通滤波带来的时间延迟问题^[4]。

2）基于模型预测静态规划的减摇性能：预测控制的终端约束$ {\sigma _{T}} = 2 $，预测步长为4，加权矩阵分别选为$ {\boldsymbol{Q}} = {\text{diag}}\left[ {\begin{array}{*{20}{c}} 8&8&8&8 \end{array}} \right] $，$ {\boldsymbol{R}} = {\text{diag}}\left[ {\begin{array}{*{20}{c}} 1&1 \end{array}} \right] $。仿真结果如图4～图5所示，在模型预测静态规划控制下升沉幅度比无控减少60%左右，纵摇幅度比无控减少65%左右。模型预测静态规划的平均计算时间为0.00705 s，文献[8]的传统预测控制平均计算时间为0.0467 s，本文所提的方法满足减摇控制的实时性要求。

3）不同预测控制策略的减摇性能对比：将本文所提的模型预测静态规划与文献[8]的传统预测控制进行对比，两者的加权矩阵，预测步长选为相同进行仿真，如图6 ～ 图7所示，可以看出基于模型预测静态规划的高速多体船减摇控制效果明显优于传统的预测控制。

5 结　语

为解决高速多体船的垂向运动幅度过大问题，提出一种基于模型预测静态规划的减纵摇方法。结果表明：

1）采用自适应边界层宽度矩阵的平滑变结构滤波器估计升沉速度、纵摇角速度，对多体船模型不确定性、不符合高斯分布的海浪噪声，具有较好的鲁棒性。

2）基于模型预测静态规划的减纵摇控制策略，通过迭代更新获得预测控制的最优控制量，降低了计算复杂度。