-
高速多体船采用流线型支柱对排水体积和主体部分进行连接，具有航行阻力小、机动性和横向稳定性好等优势[1-2]，成为船舶制造领域广泛关注的对象。然而，由于多体船具有细长的侧体，在高速航行时比单体船更容易产生剧烈的垂向运动，使得升沉和纵摇运动幅度过大，易引起失速、艏部砰击、船体结构损坏、乘员晕船等现象，对航行安全和船上设备的运行造成不利的影响[3-4]。因此，应采取有效措施降低升沉和纵摇运动幅度。在船体上安装T型翼和压浪板2类附体是减摇的重要手段之一。T型翼是一种有效的减摇附体，安装在船艏附近，其在航行时会有效增加船舶的纵向阻尼，进而减少升沉和纵摇的运动幅度，改善耐波性[5-6]，但会导致航行阻力加大；而加装压浪板可以改善船舶运动姿态，达到减阻的目的。
将主动控制引入减摇附体控制系统中，随着多体船航态的改变自动调节T型翼和压浪板的攻角，可以显著增加附体的恢复力和力矩，提高其减摇效果[7]。但是，多体船的升沉和纵摇运动存在强耦合[8]，并且T型翼和压浪板有严格的位置饱和约束[1]。T型翼以水平方向为基准在
${\rm{ [ - 1}}{{\rm{5}}^ {\circ} }\sim {\rm{1}}{{\rm{5}}^ {\circ} }{\rm{]}}$ 范围内转动，压浪板以水平方向为基准在$[{0^ {\circ} }\sim {\rm{1}}{{\rm{5}}^ {\circ} }{\rm{] }}$ 范围内转动。减摇附体若长期处于饱和状态，控制性能会下降，甚至会导致动态失速，造成多体船倾覆。与此同时，附体减摇控制信号是影响减摇性能的基础，研究表明，采用纵摇角速度信号控制策略对纵摇角的减摇效果更为明显，采用升沉速度信号控制策略对升沉运动的减摇效果更明显[7]。实际上，传感器只能测量输出升沉和纵摇量，并且受到复杂的海浪有色干扰影响，因此需要估计升沉速度和纵摇角速度。目前，多体船减摇控制和状态估计相关文献较少。针对多体船的升沉/纵摇耦合模型，文献[8]设计了比例微分（PD）减摇控制律，基于传递函数设计升沉/纵摇解耦矩阵，大幅限制了升沉和纵摇运动幅度，但是解耦矩阵鲁棒性较差，需要精确获得多体船水动力学系数。为了提高鲁棒性，ARANDA等建立了从T型翼和压浪板到升沉、纵摇运动的传递函数，设计了高阶定量反馈纵向减摇控制器；但该方法设计的控制器，其控制性能比较保守[9]。为了易于实现减摇控制工程，JAVAD等提出了一种非线性最小时间控制方法，当模型运动方向改变时控制T型翼和压浪板的攻角转到反向最大值，但是采用这种方法减摇附体一直处于饱和的正负约束边界，长期使用对减摇附体的寿命有严重影响[10]。
预测控制是提高多体船减摇性能和解决附体约束的有效途径之一[11-12]，其原因是：1) 预测模型可以有效处理升沉和纵摇的多变量耦合，不需要人为的解耦；2) 预测控制通过约束优化有效处理减摇附体的约束，避免附体的长期饱和。本文拟提出考虑附体输入约束的多体船减纵摇控制方法：首先建立多体船的垂向控制模型，分析升沉和纵摇运动的耦合性；其次基于成型滤波器理论将海浪有色干扰白化处理，构建扩展的多体船状态估计模型，采用自适应卡尔曼在线估计多体船的升沉速度和纵摇速度；最后提出考虑具有反馈校正的有约束的预测控制算法，定义系统实际状态值与预测值之间的误差，获得带有时变误差校正的预测模型，通过误差反馈校正提高减摇控制鲁棒性，通过仿真验证算法的有效性。
-
由于多体船的船体关于纵向中截面对称，故其横向运动和垂向运动无耦合。T水翼和压浪板利用翼面产生的恢复力和力矩来抵消波浪的力和力矩，从而减小升沉和纵摇幅度。在海浪扰动作用下，升沉和纵摇耦合运动模型为[1]：
$$\begin{split} & (m + {a_{33}}){\ddot x_3} + {b_{33}}{\dot x_3} + {c_{33}}{x_3} + {a_{35}}{\ddot x_5} +\\& {\kern 1pt} {\kern 1pt} {b_{35}}{\dot x_5} + {c_{35}}{x_5} = {F_{{\rm{T - foil}}}} + {F_{{\rm{flap}}}} + {F_{{\rm{wave}}}} \end{split}$$ (1) $$\begin{split} & ({{{I}}_{55}} + {a_{55}}){\ddot x_5} + {b_{55}}{\dot x_5} + {c_{55}}{x_5} + {a_{53}}{\ddot x_3} + \\&{b_{53}}{\dot x_3}{\kern 1pt} + {c_{53}}{x_3} = {M_{{\rm{T - foil}}}} + {M_{{\rm{flap}}}} + {M_{{\rm{wave}}}} \end{split}$$ (2) 式中：
$m$ 为多体船的质量；${{{I}}_{55}}$ 为多体船关于$y$ 轴的转动惯量；${a_{33}},{a_{55}}$ 为多体船的附加质量和附加转动惯量；${b_{33}},{b_{55}}$ 为系统的阻尼系数；${c_{33}},{c_{55}}$ 为系统的恢复力系数；${a_{35}},{a_{53}},{b_{35}},{b_{53}},{c_{35}},{c_{53}}$ 为力与力矩的耦合项系数；${x_3},{x_5}$ 分别为升沉位移和纵摇角；${\dot x_3},{\dot x_5}$ 分别为升沉速度和纵摇角速度；${\ddot x_3},{\ddot x_5}$ 分别为升沉加速度和纵摇角加速度；${F_{{\rm{T - foil}}}},{M_{{\rm{T - foil}}}}$ 分别为T型水翼升力和升力矩；${F_{{\rm{flap}}}},{M_{{\rm{flap}}}}$ 分别为压浪板提供的力和力矩；${F_{{\rm{wave}}}},{M_{{\rm{wave}}}}$ 分别为海浪干扰力和力矩。从式（1）和式（2）可以看出，多体船的纵摇和升沉运动具有相互耦合特点，即
${a_{35}}{\ddot x_5} + {\kern 1pt} {\kern 1pt} {b_{35}}{\dot x_5} + {c_{35}}{x_5}$ 和${a_{53}}{\ddot x_3} + {b_{53}}{\dot x_3}{\kern 1pt} + {c_{53}}{x_3}$ 是互相耦合的运动项。为便于控制器设计，进一步将升沉和纵摇的耦合运动学方程转化为下面的状态空间形式：$${\dot{{x}}} = {{{A}}_1}{{x}} + {{{B}}_1}{{u}} + {{{B}}_{{\rm{wave}}}}{{w}}$$ (3) 其中，
$$ \begin{array}{l} {{x}} = \left[ {\begin{array}{*{20}{c}} {{{\dot x}_3}} \\ {{{\dot x}_5}} \\ {{x_3}} \\ {{x_5}} \end{array}} \right],{{u}} = \left[ \begin{array}{l} {\alpha _1} \\ {\alpha _2} \end{array} \right],{{w}} = \left[ {\begin{array}{*{20}{c}} {{F_{{\rm{wave}}}}} \\ {{M_{{\rm{wave}}}}} \end{array}} \right],\;\;{{{A}}_1} = \left[ {\begin{array}{*{20}{c}} { - {{{A}}_0^{ - 1}}{{{B}}_0}}&{ - {{{A}}_0^{ - 1}}{{{C}}_0}} \\ {{{{I}}_{2 \times 2}}}&{{{{0}}_{2 \times 2}}} \end{array}} \right] , \;{{{B}}_{{\rm{wave}}}} = \left[ {\begin{array}{*{20}{c}} { - {{{A}}_0^{ - 1}}} \\ {{{{0}}_{2 \times 2}}} \end{array}} \right]\\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;{{{B}}_1} = \left[ {\begin{array}{*{20}{c}} { - {{{A}}_0^{ - 1}}} \\ {{{{0}}_{2 \times 2}}} \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {\dfrac{1}{2}{{{C}}_{{\rm{L1}}}}\rho {V^2}S}&{\dfrac{1}{2}\rho {{{A}}_{{\rm{T - foil}}}}\dfrac{{\partial {{{C}}_{{\rm{L1}}}}}}{{\partial \alpha }}{V^2}} \\ {\dfrac{1}{2}{{{C}}_{{\rm{L}}1}}\rho {V^2}S \times {l_{{\rm{flap}}}}}&{\dfrac{1}{2}\rho {{{A}}_{{\rm{T - foil}}}}\dfrac{{\partial {{{C}}_{{\rm{L1}}}}}}{{\partial \alpha }}{V^2} \times {l_{{\rm{T - foil}}}}} \end{array}} \right]\\ \;\;\;\;\;\;\;\;\;\;\;\;\;\; {{{A}}_0} = \left[ {\begin{array}{*{20}{c}} {m + {a_{33}}}&{{a_{35}}}\\ {{a_{53}}}&{{I_{55}} + {a_{55}}} \end{array}} \right],{{{B}}_0} = \left[ {\begin{array}{*{20}{c}} {{b_{33}}}&{{b_{35}}}\\ {{b_{53}}}&{{b_{55}}} \end{array}} \right],{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {{{C}}_0} = \left[ {\begin{array}{*{20}{c}} {{c_{33}}}&{{c_{35}}}\\ {{c_{53}}}&{{c_{55}}} \end{array}} \right] \end{array} $$ 式中：
$\rho $ 为海水密度；${{{A}}_{{\rm{T - foil}}}}$ 为T型水翼面积；${{{C}}_{\rm{L}}}$ 为水翼的升力系数；$V$ 为流体相对水翼的速度；${{{C}}_{{\rm{L1}}}}$ 为压浪板升力系数；$S$ 为压浪板的有效面积；${\alpha _1}$ 为压浪板攻角；${\alpha _2}$ 为T型翼攻角；${l_{{\rm{flap}}}},{l_{{\rm{T - foil}}}}$ 分别为压浪板和T型翼的力臂；${{{0}}_{2 \times 2}}$ 为零矩阵；${{{I}}_{2 \times 2}}$ 为单位矩阵。 -
高速多体船的减摇控制目标是控制升沉/纵摇运动，降低升沉/纵摇运动幅度，并抑制参数不确定性和海浪扰动。在多体船减摇控制中，需要获取升沉位移、升沉速度、纵摇、纵摇角速度4个状态参数，而传感器只输出升沉位移和纵摇角2个状态量，需要在线估计升沉速度和纵摇角速度。实际中，窄带随机海浪对船舶的扰动力和力矩是一种平稳随机过程，为有色噪声，若采用传统卡尔曼滤波方法进行数据滤波处理，会导致状态估计失真或者滤波发散。
我们基于有理谱理论建立成型滤波器，对多体船垂向控制系统模型进行扩展，以使系统输入为白噪声，再应用卡尔曼滤波对系统进行最优估计。设
${S_1}(\omega )$ 为有色海浪噪声的功率谱，${S_0}(\omega )$ 为白噪声功率谱，则有：$${S_1}(\omega ) = G(\omega ){G^*}(\omega ){S_0}(\omega )$$ (4) 式中：
$G(\omega )$ 为成型滤波器的传递函数；${G^*}(\omega )$ 为其共轭形式。对于波浪力和力矩这样的窄带谱，经验上可采用下面形式[13]：$$G({{s}}) = \frac{{{b_1}{{s}}}}{{{{{s}}^2} + {a_1}{{s}} + {a_2}}}$$ (5) 式中：
${a_1},{a_2},{b_1}$ 参数可以采用最小二乘估计。分别构造海浪扰动力和力矩的成型滤波器：$${G_3}\left( {{s}} \right) = \frac{{{b_{31}}{{s}}}}{{{{{s}}^2} + {a_{31}}{{s}} + {a_{32}}}}$$ (6) $${G_5}\left( {{s}} \right) = \frac{{{b_{51}}{{s}}}}{{{{{s}}^2} + {a_{51}}{{s}} + {a_{52}}}}$$ (7) 式中：
${{s}}$ 为复变量；${G_3}$ 为海浪扰动力的成形滤波器传递函数；${G_5}$ 为海浪扰动力矩的成形滤波器传递函数；${a_{31}},{a_{32}},{b_{31}}$ 为传递函数${G_3}$ 的模型参数；${a_{51}},{a_{52}},{b_{51}}$ 为传递函数${G_{\rm{5}}}$ 的模型参数。将式（6）和式（7）分别转化为状态空间形式：
$$\left\{ \begin{array}{l} {{\dot{{m}}}_{\rm{3}}} = {{{A}}_{\rm{3}}}{{{m}}_{\rm{3}}} + {{{B}}_{\rm{3}}}{{{w}}_{\rm{3}}} \\ {{{y}}_{\rm{3}}} = {{{C}}_{\rm{3}}}{{{m}}_{\rm{3}}} \end{array} \right.$$ (8) $$\left\{ \begin{array}{l} {{\dot {{m}}}_{\rm{5}}} = {{{A}}_{\rm{5}}}{{{m}}_{\rm{5}}} + {{{B}}_{\rm{5}}}{{{w}}_{\rm{5}}} \\ {{{y}}_{\rm{5}}} = {{{C}}_{\rm{5}}}{{{m}}_{\rm{5}}} \end{array} \right.$$ (9) 式中：
${{{m}}_{\rm{3}}}$ ，${{{y}}_{\rm{3}}}$ ，${{{w}}_{\rm{3}}}$ ，${{{A}}_{\rm{3}}}$ ，${{{B}}_{\rm{3}}}$ ，${{{C}}_{\rm{3}}}$ 分别为海浪扰动力的状态变量、输出变量、输入白噪声、状态转移矩阵、控制矩阵、输出矩阵；${{{m}}_5}$ ，${{{y}}_5}$ ，${{{w}}_5}$ ，${{{A}}_5}$ ，${{{B}}_5}$ ，${{{C}}_5}$ 分别为海浪扰动力矩的状态变量、输出变量、输入白噪声、状态转移矩阵、控制矩阵、输出矩阵。将以上2个状态空间合并得到：
$$ \begin{array}{c} \left\{ \begin{aligned} & {{\dot{{M}}}_1} = {{{A}}_{\rm{f}}}{{{M}}_1} + {{{B}}_{\rm{f}}}{{W}} \\& {{{Y}}_{\rm{f}}} = {{{C}}_{\rm{f}}}{{{M}}_1} \end{aligned}\right.\\ {{{M}}_1^{\rm{T}}} = [ {{{{m}}_3}}\;\;{{{{m}}_5}} ],{{{W}}^{\rm{T}}} = [ {{{{w}}_3}}\;\;{{{{w}}_5}} ],{{Y}}_{\rm{f}}^{\rm{T}} = [ {{{{y}}_3}}\;\;{{{{y}}_5}} ]\\ {{{A}}_{\rm{f}}} \!=\! diag[ {{{{A}}_3}}\;\;{{{{A}}_5}} ],{{{B}}_{\rm{f}}} \!=\! diag[ {{{{B}}_3}}\;\;{{{{B}}_5}} ],{{{C}}_{\rm{f}}} \!=\! diag[ {{{{C}}_3}}\;\;{{{{C}}_5}} ] \end{array}$$ (10) 式中，
${{W}}$ 为系统的输入白噪声。构造扩展状态空间方程：
$$\begin{array}{c} \left[ {\begin{array}{*{20}{c}} {{\dot{{x}}}} \\ {{{{\dot{{M}}}}_1}} \end{array}} \right]{{ = }}\left[ {\begin{array}{*{20}{c}} {{{{A}}_1}}&{{{{B}}_1}{{{C}}_{\rm{f}}}} \\ {{0}}&{{{{A}}_{\rm{f}}}} \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {{x}} \\ {{{{M}}_1}} \end{array}} \right]{{ + }}\left[ {\begin{array}{*{20}{c}} {{{{B}}_1}} \\ {{0}} \end{array}} \right]{{u}} + \left[ {\begin{array}{*{20}{c}} {{0}} \\ {{{{B}}_{\rm{f}}}} \end{array}} \right]{{W}} \\ {{Y }}= \left[ {\begin{array}{*{20}{c}} {{C}}&{{0}} \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {{x}} \\ {{{{M}}_1}} \end{array}} \right] +{{ V}}\\[-15pt] \end{array} $$ (11) 式中，
${{V}}$ 为2维量测噪声，一般情况下可以认为是白噪声。对式(11)进行离散化处理，得到离散形式扩展状态空间方程：
$$ \begin{split} & {{\bar {{x}}}_{k + 1}} = {{\varPhi}} {{\bar {{x}}}_k} + {{\varGamma}} {{{u}}_k} + {{{\varGamma}} _{\rm{f}}}{{{W}}_k} \\&\qquad {{{Y}}_k} = \bar {{C}}{{\bar {{x}}}_k} + {{{V}}_k} \end{split} $$ (12) 式中：
${\bar {{x}}_k}$ 为扩展状态变量；${{{u}}_k}$ 为控制输入；${{{Y}}_k}$ 为输出向量；${{{W}}_k}$ 为输入白噪声${{W}}$ 的采样值；${{{V}}_k}$ 为量测噪声${{V}}$ 的采样值；${{\varPhi}} $ 为状态转移矩阵；${{\varGamma}} $ 为控制矩阵；${{{\varGamma}} _{\rm{f}}}$ 为干扰矩阵；$\bar {{C}}{\rm{ = }}[ {{C}}\;\;{{0}} ]$ 为输出矩阵；$k$ 为采样时刻，$ k=1,2,\cdots ,{\rm{n}}$ 。采用扩展后的离散状态空间方程式（12），可得到以下卡尔曼滤波递推式：
$$ \begin{split} & {\widehat{\overline{{x}}}}_{k\mid k-1}={{\varPhi}} {\widehat{\overline{{x}}}}_{k-1}+{{\varGamma}} {{{u}}}_{k-1} \\ &{{P}}_{k\mid k-1}={{A}}_{1}{{P}}_{k-1}{{A}}_{1}{}^{\rm{T}}+{{{\varGamma}} }_{\rm{f}}{{Q}}_{\rm{w}}{{{\varGamma}} }_{\rm{f}}^{\rm{T}}\\ &{{K}}_{k}={{P}}_{k\mid k-1}{\overline{{C}}}^{\rm{T}}{({\overline{{C}}}{{P}}_{k\mid k-1}{{\overline{{C}}}}^{\rm{T}}+{{Q}}_{\rm{v}})}^{-1}\\ &{\widehat{\overline{{x}}}}_{k}={\widehat{\overline{{x}}}}_{k\mid k-1}+{{K}}_{k}({{Y}}_{k}-{\overline{{C}}}{\widehat{{x}}}_{k\mid k-1})\\ &{{P}}_{k}=({{I}}-{{K}}_{k}{\overline{{C}}}){{P}}_{k\mid k-1} \end{split}$$ (13) 估计的海浪扰动为：
$$\left[ {\begin{array}{*{20}{c}} {{F_{{\rm{wave}}}}} \\ {{{{M}}_{{\rm{wave}}}}} \end{array}} \right]{\rm{ = }}[ {{0}}\;\;{{C}} ]{{\hat{\bar{{x}}}}_k}$$ (14) 式中，
${{{\hat{\bar{{x}}}}}_{k\mid k - 1}}$ 为$k$ 时刻的预测状态；${{{\hat{\bar{{x}}}}}_{k - 1}}$ 为$k{\rm{ - 1}}$ 时刻的状态估计值；${{{\hat{\bar{{x}}}}}_k}$ 为$k$ 时刻的状态估计值；${{{P}}_{k\mid k - 1}}$ 为${{{\hat{\bar{{x}}}}}_{k\mid k - 1}}$ 对应的协方差；${{{P}}_{k - 1}}$ 为${{{\hat{\bar{{x}}}}}_{k - 1}}$ 对应的协方差；${{{K}}_k}$ 为卡尔曼增益；${{{Q}}_{\rm{w}}}$ 为系统噪声方差；${{{Q}}_{\rm{v}}}$ 为量测噪声方差。由于海浪噪声及量测噪声统计特性难以事先确定，采用时变噪声统计的自适应Sage-Husa卡尔曼滤波算法估计状态。为实时在线估计
${{{Q}}_{\rm{w}}}$ 和${{{Q}}_{\rm{v}}}$ ，防止因噪声不确定性造成滤波发散，引入遗忘因子$b$ 限制滤波器的记忆长度，增加新测量值在当前估计的作用，公式为[14]：$${{{Q}}_{\rm{w}}}(k) \!=\! [1 \!-\! {\textit{z}}(k \!-\! 1)]{{{Q}}_{\rm{w}}}(k \!-\! 1)\! +\! {\textit{z}}(k \!-\! 1)[{{{K}}_k}{{{V}}_k}{{{V}}_k}^{\rm{T}}{{K}}_k^{\rm{T}} \!+\! {{{P}}_k}]$$ (15) $${{{Q}}_{\rm{v}}}(k) \!=\! [1 \!-\! {\textit{z}}(k \!-\! 1)]{{{Q}}_{\rm{v}}}(k \!-\! 1)\! +\! {\textit{z}}(k \!-\! 1)[{{{V}}_k}{{{V}}_k}^{\rm{T}} \!-\! {\bar{{C}}}{{{P}}_{k|k \!-\! 1}}{{\bar{{C}}}^{\rm{T}}}]$$ (16) 式中；
${\textit{z}}\left( {k - 1} \right) = \left( {1 - b} \right)/\left( {1 - {b^k}} \right)$ ，为加权系数；遗忘因子$b$ 的取值范围一般为$0.95 \leqslant b \leqslant 0.99$ ；${{{V}}_k} = {{{Y}}_k} - {\bar{{C}}}{{{\hat{\bar{{x}}}}}_k}$ ，为新息序列。 -
高速多体船减纵摇是通过调节压浪板的攻角
${\alpha _1}$ 和T型翼的攻角${\alpha _2}$ ，使多体船的升沉位移和纵摇角达到零值附近。由于压浪板和T型翼附体存在严格的输入攻角约束，采用预测控制解决输入约束。采用一阶前向欧拉法将连续状态方程式(3)转为化为离散状态方程：
$$ \begin{split} & \;\;\;\;\;\;\;\;\;\;\left\{ \begin{array}{l} {{\hat x}}(k + 1) = {{A}}{\hat{{x}}}(k) + {{Bu}}(k) + {\hat{{D}}}(k) \\ \hat{{y}}(k) = {{C}}{\hat{{x}}}(k) \end{array} \right.\\& {\hat{{D}}}(k) = T{{{B}}_{{\rm{wave}}}}[ {{0}}\;\;{{C}} ]{{\hat{\bar{{x}}}}_k} , {{C}} = \left[ {\begin{array}{*{20}{c}} 0&0&1&0 \\ 0&0&0&1 \end{array}} \right] \end{split}$$ (17) 式中：
$T$ 为采样周期；${{B}}$ 为控制矩阵；$\hat {{y}}(k)$ 为$k$ 时刻的输出预测值；${\hat{{D}}}(k)$ 为$k$ 时刻随机海浪干扰的估计值。定义预测状态变量$\hat {{x}}\left( {k + N - 1|k} \right)$ ，则式（17）的预测模型为：$$\left\{ \begin{aligned} & {\hat{{x}}}\left( {k + N|k} \right) = {{A}}{\hat{{x}}}\left( {k + N - 1|k} \right) +\\&\quad {{Bu}}\left( {k + N - 1|k} \right) + {\hat{{D}}}\left( {k + N} \right) \\& {\hat{{y}}}\left( {k + N|k} \right) = {{C}}{\hat{{x}}}\left( {k + N - 1|k} \right) \end{aligned}\right.$$ (18) 式中：
${\hat{{D}}}\left( {k + N} \right)$ 为在$k + N$ 时刻的随机海浪干扰估计值，一般来说，将来时刻的海浪随机干扰无法获得，假设${\hat{{D}}}(k) = {\hat{{D}}}(k + 1) =\cdots = {\hat{{D}}}(k + {N_{\rm{p}}})$ ；${\hat{{x}}}\left( {k + N|k} \right)$ 为在$k$ 时刻预测的$k + N$ 时刻$x$ 的值，$N = 1,2, \cdots ,{\rm{n}}$ 。减纵摇控制目标是减小高速多体船的升沉位移和纵摇角度，采用二次型的目标函数来进行控制，公式为：$$\begin{split} & J\left( {x,u} \right) = \sum\limits_{i = 1}^{{N_{\rm{p}}}} {{{\hat{{y}}}}{{\left( {k + i|k} \right)}^{\rm{T}}}{{Q}}{\hat{{y}}}\left( {k + i|k} \right) + }\\&\qquad \sum\limits_{i = 1}^{{N_{\rm{c}}}} {{{u}}{{\left( {k + i|k} \right)}^{\rm{T}}}{{Ru}}\left( {k + i|k} \right)} \end{split} $$ (19) 式中：
${N_{\rm{p}}}$ 为预测时域；${N_{\rm{c}}}$ 为控制时域；${{Q}}$ 和${{R}}$ 分别为控制性能和控制输入的加权项。对减摇附体控制输入
${{u}}$ 进行约束，即：$${{{u}}_{\min }} \leqslant {{u}} \leqslant {{{u}}_{\max }}$$ (20) 式中：
${{{u}}_{\min }}$ 和${{{u}}_{\max }}$ 分别为控制量${{u}}$ 的下界和上界。高速多体船的预测控制减摇问题是以式(19)为目标函数、以式(20)为输入约束条件的在线优化，即：$$\begin{array}{l} \mathop {{{{u}}^*}\left( {k + i|k} \right)}\limits_{i = 0,1, \cdots ,{N_{\rm{p}}}} = \min J\left( {{{x}},{{u}}} \right) \\ {\rm{s}}{\rm{.t}}{\rm{. \;\;(17) (20)}} \end{array} $$ (21) 由于式(18)考虑系统存在海浪扰动，干扰估计值
${\hat{{D}}}\left( {k + N} \right)$ 存在于预测状态变量中，即干扰被引入到系统的预测模型之中，代入式（19）的目标函数中参与优化，保证了控制器输出在设计的性能指标下仍具有较高控制性能和强抗干扰能力。 -
受到高速多体船的模型参数偏差、海浪时变干扰等因素影响，多体船的模型预测状态和实际状态两者不可避免存在误差，导致系统的控制精度和鲁棒性下降。为了减小模型参数偏差、海浪干扰等不确定因素对减摇性能的影响，本文提出误差校正方法。定义
$k$ 时刻的输出${{y}}\left( k \right)$ 与预测值${{C}}{\hat{{x}}}\left( {k|k - 1} \right)$ 之间的误差${{e}}\left( k \right)$ ：$${{e}}\left( k \right) = {{y}}\left( k \right) - {{C}}{\hat{{x}}}\left( {k|k - 1} \right)$$ (22) 利用误差
${{e}}\left( k \right)$ 校正预测模型(18)在$k + 1$ 时刻的状态预测量，对${{e}}\left( k \right)$ 取校正矩阵${{H}} = diag[ 1 1\;\; \cdots \;\;1 ]$ ，来修正每一步预测值，校正后为：$${{\hat{{x}}}}\left( {k + 1|k} \right) = {{A}}{\hat{{x}}}\left( {k|k} \right) + {{Bu}}\left( {k|k} \right) + {\hat{{D}}}\left( k \right) + {{He}}\left( k \right)$$ (23) 在预测控制的反馈校正中，对于较慢变化的干扰补偿一般采用
${{e}}\left( {k + i|k} \right) = \cdots = {{e}}\left( k \right)$ ，但是多体船的海浪随机干扰变化较快，这种近似常值补偿处理会带来一定的误差，采用线性化误差补偿方式，即$${{e}}\left( {k + i|k} \right) = {{e}}\left( k \right) + i[{{e}}\left( k \right) - {{e}}\left( {k - 1} \right)]$$ (24) 用误差
${{e}}\left( k \right)$ 校正预测模型(18)的状态方程在$k + N$ 时刻($N = 1,2,3, \cdots $ )的状态预测量${\hat{{x}}}\left( {k + N|k} \right)$ ，得到：$$ \begin{split} & \qquad\qquad\qquad{{\hat{{x}}}}\left( {k + N|k} \right) = {{{A}}^N}{{\hat{{x}}}}\left( {k|k} \right) + \\ & \sum\limits_{i = 1}^N {{{{A}}^{N - i}}[ {{{Bu}}\left( {k + i - 1} \right) + {\hat{{D}}}\left( {k + i - 1} \right) + {{He}}\left( {k + i - 1|k} \right)} ]} \end{split} $$ (25) $$ \begin{split} & \qquad\qquad\qquad{\hat{{y}}}\left( {k + N|k} \right) = {{C}}{{{A}}^N}{\hat{{x}}}\left( {k|k} \right) + \\ & {{C}}\sum\limits_{i = 1}^N {{{{A}}^{N - i}}[ {{{Bu}}\left( {k + i - 1} \right) + {\hat{{D}}}\left( {k + i - 1} \right) + {{He}}\left( {k + i - 1|k} \right)} ]} \end{split} $$ (26) 定义预测时域
${N_{\rm{p}}}$ 内的输出向量${{{Y}}_{\rm{1}}}(k)$ ，控制时域${N_{\rm{c}}}$ 内的控制量${{U}}(k)$ ，则预测时域内的预测输出量可用式(27)计算：$${{{Y}}_{\rm{1}}}(k) = {{\varPsi}} {\hat{{x}}}(k) + {{\varTheta}}{{U}}(k) + {{\varOmega }}{{{\varPhi}}_{\rm{1}}}(k) + {{\varOmega}} {{E}}(k)$$ (27) 其中，
$$ \begin{array}{c} {{{Y}}_{\rm{1}}}(t) = \left[ {\begin{array}{*{20}{c}} {{\hat{{y}}}\left( {k + 1|k} \right)} \\ {{\hat{{y}}}\left( {k + 2|k} \right)} \\ \vdots \\ {{\hat{{y}}}\left( {k + {N_{\rm{p}}}|k} \right)} \end{array}} \right] , \\ {{U}}(k) = \left[ {\begin{array}{*{20}{c}} {{{u}}\left( {k|k} \right)} \\ {{{u}}\left( {k + 1|k} \right)} \\ \vdots \\ {{{u}}\left( {k + {N_{\rm{c}}} - 1|k} \right)} \end{array}} \right] , {{\varPsi }} = \left[ {\begin{array}{*{20}{c}} {{{CA}}} \\ {{{C}}{{{A}}^2}} \\ \vdots \\ {{{C}}{{{A}}^{{N_{\rm{p}}}}}} \end{array}} \right],\\ {{\varTheta }} = \left[ {\begin{array}{*{20}{c}} {{{CB}}}&{{{{0}}_{4 \times 2}}}& \cdots &{{{{0}}_{4 \times 2}}} \\ {{{CAB}}}&{{{CB}}}& \cdots &{{{{0}}_{4 \times 2}}} \\ \vdots & \vdots &{}& \vdots \\ {{{C}}{{{A}}^{{N_{\rm{p}}} - 1}}{{B}}}&{{{C}}{{{A}}^{{N_{\rm{p}}} - 2}}{{B}}}& \cdots &{{{C}}{{{A}}^{{N_{\rm{p}}} - {N_{\rm{c}}}}}{{B}}} \end{array}} \right] , \\{{\varOmega }} = \left[ {\begin{array}{*{20}{c}} {{C}}&{{{{0}}_{4 \times 4}}}& \cdots &{{{{0}}_{4 \times 4}}} \\ {{{CA}}}&{{C}}& \cdots &{{{{0}}_{4 \times 4}}} \\ \vdots & \vdots &{}& \vdots \\ {{{C}}{{{A}}^{{N_{\rm{p}}} - 1}}}&{{{C}}{{{A}}^{{N_{\rm{p}}} - 2}}}& \cdots &{{C}} \end{array}} \right],\\ {{E}}(k) = \left[ {\begin{array}{*{20}{c}} {{{He}}(k)} \\ {{{He}}(k + 1)} \\ \vdots \\ {{{He}}(k + {N_{\rm{p}}})} \end{array}} \right] , {{{\varPhi}}_{\rm{1}}}(k) = \left[ {\begin{array}{*{20}{c}} {{\hat{{D}}}(k)} \\ {{\hat{{D}}}(k + 1)} \\ \vdots \\ {{\hat{{D}}}(k + {N_{\rm{p}}})} \end{array}} \right] \end{array} $$ 利用式(27)将二次型的目标函数式(19)转化成向量形式：
$$J\left( {{{x}},{{u}}} \right) = {{{Y}}_{\rm{1}}}{(k)^{\rm{T}}}{\bar{{Q}}}{{{Y}}_{\rm{1}}}(k) + {{U}}{(k)^{\rm{T}}}{\bar{{RU}}}(k)$$ (28) 其中
${\bar{{Q}}} = diag[ {{Q}}\;\;{{Q}}\;\;\cdots \;\;{{Q}} ]$ ，${\bar{{R}}} = diag[ {{R}}\;\;{{R}}\;\;\cdots \;\;{{R}} ]$ 。进一步转化得到：
$$ \begin{split} & \qquad J\left( {{{x}},{{u}}} \right) = {{U}}{(k)^{\rm{T}}}({{{\varTheta }}^{\rm{T}}}{\bar{{Q}}}{{\varTheta }} + {\bar{{R}}}){{U}}(k) + {{{M}}^{\rm{T}}}{\bar{{Q}}}{{M}} +\\ & 2{{{M}}^{\rm{T}}}{\bar{{Q}}}{{\varTheta}} {{U}}(k) = \frac{1}{2}{\left[ {\begin{array}{*{20}{c}} {{{U}}(k)} \\ {{0}} \end{array}} \right]^{\rm{T}}}{{K}}\left[ {\begin{array}{*{20}{c}} {{{U}}(k)} \\ {{0}} \end{array}} \right] + f\left[ {\begin{array}{*{20}{c}} {{{U}}(k)} \\ {{0}} \end{array}} \right] + {{\sigma }} \end{split} $$ (29) 其中
$$ \begin{array}{c} {{M}} = {{\varPsi}} {\hat{{x}}}(k) + {{\varOmega}} {{E}}(k) + {{{{\varGamma}} }}{{{\varPhi}}_{\rm{1}}}(k) , \\ {{K}} = \left[ {\begin{array}{*{20}{c}} {2({{{\varTheta }}^{\rm{T}}}{\bar{{Q}}}{{\varTheta }} + {\bar{{R}}})}&{{0}} \\ {{0}}&{2{{M}}} \end{array}} \right] ,\\ {{f}} = [ {2{{{M}}^{\rm{T}}}{\bar{{Q}}}{{\Theta }}}\;\;{{0}} ] , {{\sigma}} = {{{M}}^{\rm{T}}}{\bar{{Q}}}{{M}} \end{array} $$ 将约束条件式(20)调整为以
${{U}}(k)$ 为优化变量，得到：$$\left[ {\begin{array}{*{20}{c}} {{{{I}}_u}} \\ { - {{{I}}_u}} \end{array}} \right]{{U}}(k) \leqslant \left[ {\begin{array}{*{20}{c}} {{{{U}}_{\max }}} \\ {{{{U}}_{\min }}} \end{array}} \right]$$ (30) 式中：
${{{I}}_u}$ 为${N_{\rm{c}}} \times {N_{\rm{c}}}$ 维度的单位矩阵；${{{U}}_{\max }} = {{{u}}_{\max }}{{{T}}_u}$ ；${{{U}}_{\min }} = {{{u}}_{\min }}{{{T}}_u}$ ；${{{T}}_u}$ 为${N_{\rm{c}}}$ 维度的单位列向量。预测控制最终转化成以式(29)为目标函数、式(30)为约束条件，并以
${{U}}(k)$ 为优化变量的二次规划问题。在每个控制周期完成对式(30)的求解之后，得到控制时域内的一系列控制输入量，即：$${{{U}}^*}(k) = {[{{{u}}^*}(k|k),{{{u}}^*}(k + 1|k), \cdots ,{{{u}}^*}(k + {N_{\rm{c}}} - 1|k)]^{\rm{T}}}$$ (31) 将该控制序列中的第1个元素
${{{u}}^*}(k|k)$ 作为实际的控制输入量，即${{u}}(k) = {{{u}}^*}(k|k)$ ，系统执行这一控制量直到$k + 1$ 时刻，并重复进行这一过程实现滚动优化。 -
根据文献[3]的多体船模型来验证设计的有限时间减摇方法有效性。多体船在高速航行时受到的海况等级为四级海况，海浪采用P-M谱进行仿真。P-M谱为一种重力波谱，其数学表达式相对简单，而且仅与海面上方的风速有关，方便计算，因此得到广泛应用。
$${S_\zeta }(\omega ) = \frac{{8.1 \times {{10}^{ - 3}}{g^2}}}{{{\omega ^5}}}\exp \left[ - 0.74{\left(\frac{g}{{{v_\zeta }\omega }}\right)^4}\right]$$ (32) 式中：
${v_\zeta }$ 为海面以上高度为19.5 m处的平均风速，m/s；g为重力加速度，${\rm{m/}}{{\rm{s}}^{\rm{2}}}$ ；${S_\zeta }(\omega )$ 为功率密度，${{\rm{m}}^2} \cdot {\rm{s}}$ ；$\omega $ 为海浪的角频率，${\rm{rad/s}}$ 。根据切片法，通过数据拟合和叠加，求取不同频率点下海浪作用于多体船的干扰力和干扰力矩。多体船以30 kn航速迎浪航行，仿真过程如下：1）高速多体船的减摇系统量测噪声是白噪声，方差阵初值取为
${R_0} \!=\! diag[20.{\rm{ }}3 \times {10^{ - 4}}\;\;2.{\rm{ }}26 \times {10^{ - 6}}]$ ；经过扩展后的系统噪声是白噪声，方差阵初值取为1。基于海浪干扰力和力矩的成型滤波器，构造扩展卡尔曼滤波器进行滤波，对升沉速度、纵摇角速度进行准确估计，仿真结果如图1和图2所示。
图 1 升沉速度估计
Figure 1. Heave velocity estimation
图 2 纵摇角速度估计
Figure 2. Pitch angular velocity estimation
2）针对带T型翼和压浪板的多体船，在无控制器与加入反馈校正情况下进行仿真，得到的升沉位移和纵摇角仿真结果如图3和4所示。由图可见，在加入反馈校正后，多体船的升沉量减少约40%，纵摇角减少约50%，T型翼和压浪板的攻角如图5和图6所示，满足输入约束，同时避免了附件长时间处于饱和状态。
图 3 多体船升沉位移
Figure 3. Heave of multi-hull ship
图 4 多体船纵摇角
Figure 4. Pitch of multi-hull ship
图 5 T型翼的攻角
Figure 5. Attack angle of T-foil
图 6 压浪板的攻角
Figure 6. Attack angle of flap
3）将反馈校正后的预测控制仿真结果与普通预测控制仿真结果比较，结果如图7和图8所示。由图可见，升沉位移和纵摇角进一步减少，减摇效果得到了有效提高。不同控制策略的标准差比较结果如表1所示，加入反馈校正的标准差比不加入控制器和普通预测控制情况下的更小，验证了算法的有效性。
表 1 不同控制策略的标准差
Table 1. Standard deviation of different control strategies
无控制 普通预测控制 加入反馈校正 升沉位移标准差（米） 0.281 8 0.183 2 0.158 3 纵摇角标准差（度） 0.808 8 0.466 3 0.379 8 
图 7 多体船升沉位移
Figure 7. Heave of multi-hull ship
图 8 多体船纵摇角
Figure 8. Pitch of multi-hull ship
-
针对高速多体船在行驶过程中纵摇和升沉运动幅度过大的问题，提出了一种考虑附体输入约束的高速多体船预测减摇方法。建立了高速多体船的垂向耦合运动模型，基于成型滤波器理论将有色海浪干扰进行白噪化建模，设计了自适应扩展卡尔曼滤波器在线估计多体船状态和海浪扰动。提出了考虑减摇附体约束的预测控制减摇策略，将线性变化反馈校正误差和估计的海浪干扰综合到预测模型中，较好地避免了附体的长期饱和问题，通过滚动优化求解提高预测控制减摇性能和鲁棒性。仿真结果表明，考虑反馈校正的预测减摇控制，有效减少了升沉和纵摇运动幅度，减摇附体输入满足约束。
Predictive control anti-pitching for high-speed multi-hull ship with appendages constraints
-
摘要:
目的 针对高速多体船的纵摇和升沉运动幅度过大，以及减摇附体有着严格输入约束问题，提出一种基于扩展卡尔曼状态估计的预测控制减纵摇方法。 方法 建立由T型翼和压浪板作为减摇附体的多体船垂向控制模型，分析升沉/纵摇运动的耦合性以及模型的不确定性。为了获取减摇控制信号，将海浪有色干扰白化处理，构建扩展的多体船状态估计模型，采用自适应卡尔曼滤波器在线估计多体船的升沉速度和纵摇角速度。在此基础上，提出有输入约束的预测控制减摇算法，定义系统实际状态值与预测值之间的误差，获得带有线性变化误差校正的预测状态模型；通过误差反馈校正提高减摇控制鲁棒性，并将减摇控制问题转化为有输入约束的二次规划问题，通过数值解法实现预测控制滚动优化求解。 结果 仿真结果表明，在考虑反馈校正的预测控制作用下，船体升沉位移减少约40%，纵摇角减少约50%。 结论 采用反馈校正的预测控制提高了减摇系统的控制精度和鲁棒性，在实际工程中具有重要意义。 Abstract:Objectives Considering the large motions of pitch and heave together with the strict input constraints of installed anti-pitching appendages, a predictive control method is proposed for vertical stabilization based on Kalman filtering. Methods The vertical control model of high-speed multihull has been established with T-foil and flap served as anti-pitching appendages, and the motion couplings of heave and pitch has been analyzed. In order to obtain the anti-pitching control signal, the wave-induced colored noise has been whitened, and adaptive extended Kalman filter has been adopted for online estimates of heave velocity and pitch angular velocity. On this basis, the predictive control has been proposed for vertical stabilization with input constraint. Defining the error between the actual state and predicted state, the predictive control model with linear varying error correction can then be obtained. Error feedback correction is used to improve the robustness of anti-pitching control, and the problem of anti-pitching control is transformed into a quadratic programming (QP) problem with input constraints, and the rolling optimization solution of predictive control is realized through numerical solution. Results The simulation results show that under the effect of predictive control considering feedback correction, the hull heave is reduced by about 40%, and the pitch angle is reduced by about 50%. Conclusions The predictive anti-pitching control with feedback correction can improve the control accuracy and robustness of the system, which is of great significance in practical engineering applications. -
表 1 不同控制策略的标准差
Table 1. Standard deviation of different control strategies
无控制 普通预测控制 加入反馈校正 升沉位移标准差（米） 0.281 8 0.183 2 0.158 3 纵摇角标准差（度） 0.808 8 0.466 3 0.379 8 
下载: 导出CSV
-
[1] FOSSEN T I. Handbook of Marine Craft Hydrodynamics and Motion Control[M]. Hoboken: Wiley, 2011. [2] ESTEBAN S, GIRON-SIERRA J M, DE ANDRES-TORO B, et al. Fast ships models for seakeeping improvement studies using flaps and T-Foil[J]. Mathematical and Computer Modelling, 2005, 41(1): 1–24. doi: 10.1016/j.mcm.2004.09.002 [3] GIRON-SIERRA J M, RECAS J, ESTEBAN S. Iterative method based on CFD data for the assessment of seakeeping control effects, considering amplitude and rate saturation[J]. International Journal of Robust and Nonlinear Control, 2011, 21(13): 1562–1573. doi: 10.1002/rnc.1653 [4] FANG C C, CHAN H S. An investigation on the vertical motion sickness characteristics of a high-speed catamaran ferry[J]. Ocean Engineering, 2007, 34(14–15): 1909–1917. doi: 10.1016/j.oceaneng.2007.04.001 [5] 原新, 张欣. 三体船纵向减摇附体设计及减摇效果分析[J]. 武汉理工大学学报(交通科学与工程版), 2017, 41(4): 554–558. YUAN X, ZHANG X. The design of longitudinal damping appendage and the effect on trimaran[J]. Journal of Wuhan University of Technology (Transportation Science & Engineering), 2017, 41(4): 554–558 (in Chinese). [6] 郑义, 董文才. 高速轻型穿浪双体船纵向运动改善措施研究[J]. 中国舰船研究, 2012, 7(2): 14–19. doi: 10.3969/j.issn.1673-3185.2012.02.003 ZHENG Y, DONG W C. Improvement of longitudinal motion performance of high speed light wave-piercing catamaran by hydrofoils[J]. Chinese Journal of Ship Research, 2012, 7(2): 14–19 (in Chinese). doi: 10.3969/j.issn.1673-3185.2012.02.003 [7] 孙一方, 宗智, 姜宜辰. 船舶在波浪上纵向运动与控制研究综述[J]. 中国舰船研究, 2020, 15(1): 1–12, 47. doi: 10.19693/j.issn.1673-3185.01751 SUN Y F, ZONG Z, JIANG Y C. Review of longitudinal motion and controls of ships on waves[J]. Chinese Journal of Ship Research, 2020, 15(1): 1–12, 47 (in Chinese). doi: 10.19693/j.issn.1673-3185.01751 [8] DE LA CRUZ J, ARANDA J, GIRON-SIERRA J M, et al. Improving the comfort of a fast ferry[J]. IEEE Control Systems Magazine, 2004, 24(2): 47–60. doi: 10.1109/MCS.2004.1275431 [9] ARANDA J, DE LA CRUZ J, DÍAZA J M. Design of a multivariable robust controller to decrease the motion sickness incidence in fast ferries[J]. Control Engineering Practice, 2005, 13(8): 985–999. doi: 10.1016/j.conengprac.2004.11.003 [10] JAVAD A, JASON L, MICHAEL R D, et al. An experimental investigation of ride control algorithms for high-speed catamarans Part 1: reduction of ship motions[J]. Journal of Ship Research, 2017, 61(1): 35–49. doi: 10.5957/jsr.2017.61.1.35 [11] ZHANG J, SUN T R, LIU Z L. Robust model predictive control for path-following of underactuated surface vessels with roll constraints[J]. Ocean Engineering, 2017, 143: 125–132. doi: 10.1016/j.oceaneng.2017.07.057 [12] KUCUKDEMIRALA I B, CAKICIB F, YAZICIC H. A model predictive vertical motion control of a passenger ship[J]. Ocean Engineering, 2019, 186: 106100. doi: 10.1016/j.oceaneng.2019.06.005 [13] 赵希人, 唐慧妍, 彭秀艳. 船舶横向运动多变量随机控制研究[J]. 船舶力学, 2004, 8(5): 35–41. doi: 10.3969/j.issn.1007-7294.2004.05.005 ZHAO X R, TANG H Y, PENG X Y. Multivariant stochastic control research on ship lateral movement[J]. Journal of Ship Mechanics, 2004, 8(5): 35–41 (in Chinese). doi: 10.3969/j.issn.1007-7294.2004.05.005 [14] 王福军, 丁小燕, 王前, 等. 自适应强跟踪Sage-Husa卡尔曼滤波器载波环设计[J]. 电光与控制, 2019, 26(10): 12–16. doi: 10.3969/j.issn.1671-637X.2019.10.003 WANG F J, DING X Y, WANG Q, et al. Design of carrier tracking loop based on adaptive strong tracking Sage-Husa Kalman filter[J]. Electronics Optics & Control, 2019, 26(10): 12–16 (in Chinese). doi: 10.3969/j.issn.1671-637X.2019.10.003 -
点击查看大图
图(8) / 表 (1)
计量
- 文章访问数: 5
- HTML全文浏览量: 2
- PDF下载量: 0
- 被引次数: 0
