1 引言
高速多体船(包括两体船和三体船)是当前新型船舶的重要发展趋势之一，与单体船相比具有更小的兴波阻力、横向稳定性好、装载量大、抗打击能力强，在军事、航运领域有着广泛的应用前景[1]。然而，高速多体船细长的侧体结构使得纵向恢复力矩变小，造成船舶在高速航行时会产生较大的升沉和纵摇幅度，出现“埋首”运动[2]，严重破坏船体结构。剧烈的升沉和纵摇运动产生较大的垂向加速度，造成乘客晕船、仪器失灵。因此，如何保证高速多体船的垂向运动稳定，减少升沉和纵摇运动幅度，成为该领域研究的重点。目前，在多体船船体上安装主动式T型翼和压浪板两类附体可实现闭环减摇控制，T型翼可以提高船体的附加质量和阻尼，压浪板用于改变船体纵倾[3]。
安装T型翼和压浪板的高速多体船是多输入输出系统，升沉和纵摇运动存在强耦合，并且航速较高导致水动力参数变化剧烈，模型具有快时变性、强不确定性等特点[4-5]。文献[4]中设计解耦的比例微分(PD)控制律实现双体船的垂向稳定控制，但是需要精确辨识水动力学参数。文献[5]中提出扩张观测器和PD结合的解耦控制策略，提高了鲁棒性。文献[6]中设计定量反馈控制律，但是控制性能比较保守。文献[7]中提出有限时间控制器和有限时间观测器复合的控制策略，增强了多体船抑制海浪扰动的能力。上述研究主要从镇定角度来设计解耦或者鲁棒控制器，没有考虑垂向控制性能的优化，而垂向控制目的是实现升沉和纵摇变化幅度和加速度越小越好。
针对升沉和纵摇耦合的多体船垂向运动模型，预测控制不仅能有效处理多变量耦合，并且实现控制性能优化。文献[8]中提出多体船的解析预测控制，避免在线复杂的数值优化。文献[9]中对随机海浪扰动力和力矩进行建模，采用卡尔曼滤波器估计状态，提出显式预测控制。然而，上述预测控制只考虑升沉和纵摇运动幅度的均方差最小，没有考虑两者的暂态响应过程。实际上在海浪扰动激励下升沉和纵摇的输出不仅幅度大，波动变化率也大，造成垂向加速度很大，严重影响高速多体船的耐波性。针对上述问题，提出基于增量预测控制的垂向运动稳定方法，其创新点为：
1) 分析随机海浪扰动的频域和时域特性，基于积分扰动模型在线估计海浪扰动力和力矩；为减少海浪平均扰动力和扰动力矩的影响，将高速多体船舶的垂向运动模型转化为增量模型；
2) 为减少随机海浪扰动下的升沉和纵摇幅度和波动变化率，将升沉和纵摇变化率和幅
度同时引入到优化目标函数中，推导出解析的预测控制律。
具体章节安排如下：首先，建立升沉和纵摇描述的高速多体船垂向运动模型，分析随机海浪的时域特性。其次，将高速多体船的垂向运动模型转化为增量模型，基于滚动优化策略推导预测控制律。最后，通过仿真验证所提算法的有效性。
2 高速多体船的垂向运动模型
高速多体船的船体是相对于中间截面对称的，T型翼和压浪板分别在船艏底部和船艉。在航行中通过调节T型翼和压浪板攻角，产生恢复力和力矩，抵消海浪引起的扰动力和力矩，抑制升沉和纵摇剧烈运动。当在康明斯方程中考虑流体记忆效应时[10]，基于升沉和纵摇两自由度描述的高速多体船垂向运动模型为
(1)
(2)
式中：
分别为升沉速度、纵摇角速度、升沉位移和纵摇角；a33,a55分别表示多体船的附加质量和附加转动惯量；b33,b55为系统的阻尼系数；c33,c55为系统的恢复力系数；a35,a53,b35,b53,c35,c53为力与力矩的耦合项系数；K33为升沉通道中由升沉运动引起的脉冲响应函数或延迟函数；K35为升沉通道中由纵摇运动引起的脉冲响应函数；K55为纵摇通道中由纵摇运动引起的脉冲响应函数；K53为纵摇通道中由升沉运动引起的脉冲响应函数；Fwave和Mwave分别表示随机海浪产生的扰动力和力矩；FT-foil、MT-foil和Fflap、Mflap分别表示T型翼和压浪板所产生的恢复力和力矩。
式中：ρ为海水密度；S1、S2分别为T型翼和压浪板的有效面积；CL1和CL2分别为T型翼和压浪板的升力系数；V表示流体相对于水翼的速度；α1、α2和lT-foil、lflap分别表示T型翼和压浪板的攻角和力臂。
为后续设计方便，将式(1)与式(2)写成下面的状态空间形式
(3)
式中：
且
从式 (3) 可以看出，高速多体船模型存在升沉和纵摇运动的强耦合，以及存在强的海浪扰动。高速多体船垂向运动控制目的是在海浪扰动下，自动调节减摇附体的攻角使得升沉和纵摇运动的幅度和加速度越小越好。
3 海浪扰动力和力矩估计
高速多体船的升沉和纵摇幅度变化剧烈，主要原因是受到海浪扰动力和力矩作用，本节研究四级海浪产生扰动的时域和频域特性，为后续的减摇控制设计提供基础。本研究中采用ITTC单参数谱作为海浪仿真功率谱函数[1,11]，即
(4)
式中:h1/3为海浪有义波高；ω为海浪主导频率。根据功率谱反演和随机信号合成理论，得到作用于两体船的扰动力和扰动力矩，如图1、图2所示。
图1 多体船的海浪扰动力
Fig.1 Disturbance wave force of the multihull
图2 多提船的海浪扰动力矩
Fig.2 Disturbance wave moment of the multihull
可以看出多体船的扰动力和力矩是一种平稳随机过程。通过统计数据知道海浪扰动力和扰动力矩均值分别为0.429 7×103 N，2.028 2×103 N·m，均不为零，并且海浪变化率Δd(k)≠0，这会影响下一节的预测控制模型准确性，需要在线实时估计d(k)，在文献[10]中假设Δd(k)=0是有局限性的。
海浪扰动力和力矩的估计常采用白噪化方法，将白噪声作为系统新的输入，对系统状态进行扩张，采用卡尔曼滤波器估计状态[12]。但是，该方法需要辨识成型滤波器的参数，并且卡尔曼滤波维数较大，计算复杂。这里将海浪扰动力和力矩作为扩张状态，基于积分扰动模型理论近似估计海浪扰动力和力矩[13]。将高速多体船模型方程式(3) 进行离散化可得
(5)
其中,量测噪声v(k)为零均值白噪声，量测方差为
假设海浪扰动的随机噪声模型为
d(k+1)=d(k)+ζk(k)
(6)
其中：ζk(k)为零均值白噪声；噪声方差为Qr。将海浪扰动作为状态的一部分，扩张后的状态为
扩展的随机线性模型如下
xe(k+1)=Aexe(k)+Beu(k)+Bewζ(k)
(7)
其中，
当量测方差按照仪表精度选取为
对式(7) 所示的模型系统采用Kalman滤波器进行估计，估计的四级海浪扰动力和力矩如图3、图4所示，可以看出在变化平滑地方海浪扰动估计较好，在峰峰值附近有估计误差，并且有一定的滞后。因此，在下一节反馈设计中，采用对时间滞后和不确定性有一定抑制能力的预测控制。
图3 估计的海浪扰动力
Fig.3 Estimated disturbance wave force
图4 估计的海浪扰动力矩
Fig.4 Estimated disturbance wave moment
4 减少变化率的增量预测控制
高速多体船的垂向运动模型方程式 (3) 是多变量，并且垂向运动控制目的是实现升沉和纵摇输出和加速度越小越好，本节提出具有优化功能的垂向稳定预测控制[13-14]。由上一节分析知道，海浪扰动力和扰动力矩均值都不为零，需要有积分作用来提高抑制海浪扰动能力。因此，本节采用增量预测模型设计预测控制律。
定义状态增量Δx(k)=x(k)-x(k-1)，控制增量Δu(k)=u(k)-u(k-1)，扰动增量
将式 (5) 的相邻时刻的状态作差，可得
Δx(k+1)=AdΔx(k)+BdΔu(k)+BwdΔd(k)
(8)
由于当前k时刻的扰动d是已知的，但未来时刻k+1，k+2…却是未知并且不可测的，所以假设在k时刻之后的可测扰动近似不变，即Δd(k+i)=0,i=1,2,…,Np。
对于系统的输出y(k)，有
y(k+1)=y(k)+CdAdΔx(k)+
CdBdΔu(k)+CdBwdΔd(k)
(9)
将y(k)扩展为一个新的状态变量，扩展后的状态空间模型为
(10)
式中，
对k时刻，定义Np步预测输出和Nc步控制输入
YNp(k+1)=[y(k+1|k),y(k+2|k),…,y(k+Np|k))]T
ΔU(k)=[Δu(k),Δu(k+1),…,Δu(k+Nc-1)]T
则基于第k步系统状态，未来k+Np步的预测输出为
YNp(k+1)=Fxu(k)+ΘΔU(k)+GΔd(k)
(11)
式中，
高速多体船的垂向运动稳定控制不仅局限于减少纵摇和升沉运动幅度，还需要减少纵摇和升沉的变化率，其原因是过大的纵摇和升沉变化率会引起多体船垂向加速度过大，进而使得晕船率过大。这里将预测时域内的纵摇和升沉变化率也引入优化目标函数，提出了下面的预测控制目标函数
J=(Yref(k+1)-YNp(k+1))TQ(Yref(k+1)-
YNp(k+1))+ΔU(k)TRΔU(k)+
(YNp(k+1)-YNp(k))TK(YNp(k+1)-YNp(k))
(12)
式中：Q=diag(q1,q2,…,qNp×2)为误差加权矩阵；R=diag(r1,r2,…,rNc×2)为控制增量加权矩阵；K=diag(k1,k2,…,k(Np-1)×2)为变化率的加权矩阵；Yref(k+1)为参考轨迹序列。联立式 (11) 与式 (12) ，根据极值必要条件，
时，J取得极小，则ΔU(k)的最优值为
ΔU(k)=(ΘTQΘ+R+ΦTKΦ)-1 [ΘTQ(Yref(k+1)-
Fxu(k)-GΔd(k))-ΦTK(Γxu(k)+ΨΔd(k))]
(13)
式中，
为了提高预测控制的鲁棒性，只取ΔU(k)中的第一个元素作为实际的控制增量，在下一采样时刻的控制增量根据下一采样时刻的系统状态重新计算。令矩阵H=[I2×2 0 … 0]，则当前时刻的MPC控制量
Δu(k)=H(ΘTQΘ+R+ΦTKΦ)-1 [ΘTQ(Yref(k+1)-
Fxu(k)-GΔd(k))-ΦTK(Γxu(k)+ΨΔd(k))]
(14)
为了验证上述闭环增量MPC的稳定性，将增量MPC控制器方程式 (14) 代入系统方程式 (10) 可得
xu(k+1)=BuH(ΘTQΘ+R+ΦTKΦ)-1·
ΘTQYref(k+1)+(Au-BuH(ΘTQΘ+R+ΦTKΦ)-1·
(ΘTQF+ΦTKΓ))xu(k)+
(Bwu-BuH(ΘTQΘ+R+ΦTKΦ)-1·
(ΘTQG+ΦTKΨ))Δd(k)
(15)
闭环系统的稳定性通过矩阵Au-BuH(ΘTQΘ+R+ΦTKΦ)-1(ΘTQF+ΦTKΓ)的特征值来判断。当加权矩阵为Q=diag(6.2×103，4.8×104)，R=diag(3.4×10-3，2.4×10-5)，K=diag(1.2×106，10×107)，计算可得6个特征值分别是0.963 7+0.138 4i，0.963 7-0.138 4i，0.976 7+0.067 5i，0.976 7-0.067 5i，0.995 2+0.000 0i，0.999 9+0.000 0i 均位于单位圆内，因此离散系统是闭环稳定。
5 仿真分析
在仿真中高速多体船以两体船为例，取采样周期T=0.05 s，动力学模型参数为[12]：
为验证增量预测控制减摇的有效性，采用四级海浪进行仿真，仿真条件同文献[9]，进行下面仿真：
1) 分别对多体船进行无控制器和增量预测控制的仿真，升沉位移和纵摇角图5和图6所示。无控时升沉位移均值为0.240 6 m，方差为0.089 1 m2，纵摇角均值为 0.736 3 deg，方差为0.014 6 deg2。在加入增量预测控制后，升沉位移均值为0.097 1 m，方差为0.015 4 m2，纵摇角均值为0.382 0 deg，方差为0.004 2 deg2。通过对比可以看出，在增量预测控制作用下多体船升沉位移减少了59%，纵摇角减少了47%，减摇控制效果明显。
图5 升沉运动
Fig.5 Comparison of heave motion
图6 纵摇运动
Fig.6 Comparison of pitch angle motion
2) 为进一步验证所提增量预测控制性能，与传统预测控制减摇进行比较[10]，结果如图7和图8所示，可以看出所提的预测控制使得升沉位移幅度进一步减少16%，纵摇运动幅度减少约11%，波动变化率明显变缓，增强了多体船抑制海浪扰动的能力。
图7 升沉运动
Fig.7 Comparison of heave motion
图8 纵摇运动
Fig.8 Comparison of pitch angle motion
6 结论
提出基于增量预测控制的高速多体船垂向运动稳定方法，抑制过大的升沉和纵摇运动幅度。仿真结果表明：
1) 采用积分扰动模型在线估计海浪扰动力和力矩，将估计值引入到预测模型中，增强多体船抑制海浪扰动的能力。
2) 基于增量模型设计预测控制减摇方法，使得升沉位移减少59%，纵摇角减少47%。
[1] Fossen T I.Handbook of marine craft hydrodynamics and motion control[M].New Jersey,John Wiley &Sons,Ltd,2011.
[2] Zheng L H,Liu Z L,Yang S,et al.Optimal design and installation of appendages for improving stability performance of a trimaran[J].Ocean Engineering,2021,238:109719.
[3] Jiang Y C,Bai J Y,Liu S J,et al.Experimental investigation of T-foil hybrid control strategy for ship motion reduction in head seas[J].Ocean Engineering,2022,243:110251.
[4] Cruz J,Aranda J,Giron-sierra J M,et al.Improving the comfort of a fast ferry[J].IEEE Control Systems Magazine,2004,24(02):47-60.
[5] Zhang J,Liu Z L,Dai X K,et al.Robust decoupled anti-pitching control of a high-speed multihull[J].Journal of Marine Science and Technology,2021,26:1112-1125.
[6] Aranda J,Cruz J,Diaz J M.Design of a multivariable robust controller to decrease the motion sickness incidence in fast ferries[J].Control Engineering Practice,2005,13(08):985-999.
[7] 张军,刘志林,戴小康,等.一种高速多体船的有限时间减纵摇控制方法[J].武汉理工大学学报(交通科学与工程版),2021,45(03):512-516.
Zhang J,Liu Z L,Dai X K,et al.A finite time anti-pitching control method for a high-speed multihull[J].Journal of Wuhan University of Technology (Transportation Science &Engineering),2021,45(03):512-516.
[8] 张军,温昊,刘志林,等.高速多体船的解析预测控制减纵摇方法[J].兵工学报,2021,42(09):1962-1971.
Zhang J,Wen H,Liu Z L,et al.Analytic prediction control method of high-speed multihull[J].Acta Armamentarii,2021,42(09):1962-1971.
[9] Zhang J,Liu Z,Sun T,et al.Explicit stochastic model predictive control for anti-pitching a high-speed multihull[J].Applied Ocean Research,2022,119:102917.
[10] Ibrahim B K,Fredi C,Hakan Y.A model predictive vertical motion control of a passenger ship[J].Ocean Engineering,2019,186:106100.
[11] 梁利华,王保华,贾鹤鸣.基于扰动观测器的 SWATH 船运动非线性预测控制[J].控制与决策,2014,29(03):431-436.
Liang L H,Wang B H,Jia H M.Nonlinear predictive control of SWATH ship motions with disturbance observer[J].Control and Decision,2014,29(03):431-436.
[12] 赵希人,陈虹丽,艾晓庸,等.舶运动纵向受扰力和受扰力矩估计方法研究[J].哈尔滨工程大学学报,2004,25(02):167-170.
Zhao X R,Chen H L,Ai X Y,et al.Estimation of longitudinal disturbing force and tortuosity of ship movement[J].Journal of Harbin Engineering University,2004,25(02):167-170.
[13] Morari M,Maeder U.Nonlinear offset-free model predictive control[J].Automatica,2012,48(09):2059-2067.
[14] Betti G,Farina M,Scattolini R.A robust MPC algorithm for offset-free tracking of constant reference signals[J].IEEE Transactions on Automatic Control,2013,58(09):2394-2400.