智能船舶的路径跟踪是指在不考虑时间约束下，使船舶跟踪预先规划的期望路径，但多源时变环境干扰与执行装置的物理约束，使船舶难以准确跟踪期望路径；因此，如何在存在时变不确定干扰下，提高跟踪精度并考虑执行装置的约束与降低能耗的目标，是本研究控制策略的出发点．为实现主动抑制扰动，韩京清教授最早提出了自抗扰控制(active disturbance rejection control，ADRC)．文献[1]将模型不确定性、时变干扰和虚拟控制量作为总扰动，设计了有限时间自抗扰控制器．文献[2]基于模糊可变船长的视线(line of sight，LOS)引导律与自抗扰航向控制器实现了在真实环境下的航迹控制．文献[3]采用基于径向基函数神经网络的非线性ADRC技术，实现了对欠驱动船舶的航迹控制．模型预测控制(model predictive control，MPC)因其能显性考虑约束条件与优化指标、可利用干扰的任何可用信息及预测模型的建立方式灵活等优点在船舶的路径跟踪控制中得到应用．文献[4]考虑了船舶的运动状态和执行装置饱和约束．文献[5]使用最小二乘矩阵处理控制输入的幅值约束．同时，为加强系统对时变干扰的抑制，文献[6]将鲁棒模型预测控制应用于船舶的路径跟踪控制中，然而干扰有界的假设将导致过于保守的结果与算法的不可行[7]．为增强控制器的鲁棒性，文献[7]采用基于干扰补偿的MPC策略实现了船舶的航向控制．文献[8]提出了一种结合扰动观测器的鲁棒非线性模型预测控制方法．文献[9]通过径向基函数神经网络对误差数据进行训练，实现对外界干扰的逼近和补偿．文献[10]利用带有扰动观测器的三阶模型预测控制在实船实验中取得了理想的控制效果．基于以上研究，为实现主动抑制扰动与最优化目标的达成，同时为提高模型预测控制器对时变干扰的鲁棒性，本研究设计了一种自抗扰模型预测控制算法．采用在线引入系统已知模型的MESO对系统的总未知干扰和状态进行实时估计；基于状态和总未知干扰的估计值构建线性预测模型；建立以跟踪精度和能量消耗为指标的目标函数，在控制输入的幅值和变化量幅值约束下进行优化问题的求解；并对观测误差造成的模型失配进行鲁棒补偿．最后对所提出的算法进行稳定性证明和仿真实验验证．1 船舶路径跟踪建模及问题描述1.1　船舶路径跟踪数学建模水面船通常只考虑纵向、横向和艏向三个自由度的运动，由固定坐标系O0-X0Y0Z0和随船坐标系O-XYZ进行描述，如图1所示，图中：x为北向位置；y为东向位置；ψ为艏向角；Vc为流速；Vw为风速；βc为流向角；βwind为风向角；βwave为浪向角．10.13245/j.hust.239269.F001图1固定坐标系与随船坐标系船舶航行受风、浪、流的影响，建立船舶运动学模型和动力学模型η˙=R(ψ)υ; (1)υ˙=M-1τ+M-1τE-M-1CRB(υ)υ-M-1CA(υr)υr-M-1Dυr, (2)式中：η=[x,y,ψ]T为船舶位置；υ=[u,v,r]T为船舶速度；R(ψ)为坐标转换矩阵；M为系统惯性矩阵；D为线性阻尼矩阵；CRB(υ)为刚体科里奥利向心力矩阵；CA(υr)为附加科里奥利向心力矩阵；υr=[ur,vr,rr]T为相对流速矩阵．各模型矩阵可表示为：M=m1000m2m30m4m5；D=d1000d2d30d4d5；CRB(υ)=00-c100-c2c1c20；CAυr=00-c300-c4c3c40．设τ为控制力矩阵，欠驱动船舶可表示为τ=[τX,0,τN]T；τX和τN分别为纵向控制力和转艏控制力矩．υr=υ-υc，υc为流速矩阵，可由流速Vc和流向角βc描述．τE=τwind+τwave，τwind和τwave分别表示风载荷和波浪载荷，分别由风速Vw、风向角βwind和平均波高Hwave、浪向角βwave描述[11]．动力学模型(2)多变量耦合、非线性程度高，为降低所设计的观测器和控制器的结构复杂度与所设计算法的计算量，对其进行线性化处理．定义f(υ)=-M-1(C(υ)υ+Dυ)；fE(υ,υc,τE)=M-1(τE+Dυc+CA(υ)υc+CA(υc)υ-CA(υc)υc);B0=M-1，式中C(υ)=CRB(υ)+CA(υr)，为科里奥利向心力矩阵．将式(2)在期望速度矩阵υd处泰勒展开可得到线性化的动力学模型，有υ˙=f(υd)+f'(υd)(υ-υd)+d+B0τ，(3)式中：d=ο(f(υ))+fE(υ,υc,τE)为总未知干扰，ο(f(υ))为f(υ)的泰勒展开高阶项；f'(υd)为f(υ)的偏导在υd处的值．可令h(υ,υd)=f(υd)+f'(υd)(υ-υd)．使用LOS引导律将船舶对路径的三维跟踪转化为对期望艏向角ψd和期望速度υd的跟踪[12]．定义控制状态量为x=[u,v,r,ψ]T，由式(1)和(3)建立线性化的船舶运动状态空间方程，有x˙=F(xd)+F'(xd)(x-xd)+Dm+Bmτ，(4)式中：F'(xd)=diag[f'(υd),∂/∂t]；F(xd)=[f(υd)T,rd]T；xd=[ud,vd,rd,ψd]T；Bm=diag[B0，0]；Dm=[dT,0]T为状态空间方程的干扰项．定义状态跟踪误差为e=x-xd，不考虑系统干扰项的参考轨迹方程为x˙d=F(xd)+Bmτd．(5)将式(4)与式(5)做差，可得跟踪误差方程e˙=Ame+Bmu+Dm，(6)式中：Am=F'(xd)；u=τ-τd．1.2　船舶路径跟踪控制问题描述设计观测器对状态量x和干扰项Dm进行估计．控制器的设计须考虑：满足执行装置的物理约束；实现船舶对期望路径的准确跟踪，即e收敛至0；对能耗进行优化．自抗扰模型预测控制算法原理如图2所示．在给定的期望速度υd和期望艏向角ψd下，由微分跟踪器(TD)得到参考路径点(xd,τd)；MESO根据可测的船舶位置信息η估计得到的状态估计值x^和总未知干扰的估计值d^；针对MESO的观测误差ε，采用鲁棒补偿(robust compensation，RC)的方式求得补偿控制输入u'；将(xd,τd)，u'，x^与d^输入至鲁棒补偿的模型预测控制算法(RC-MPC)中，在约束条件和优化指标下求得最优控制输入u*；将u'和u*之和作为总控制输入u作用于船舶，推动船舶运动．10.13245/j.hust.239269.F002图2自抗扰模型预测控制算法原理图2 MESO设计与分析假设1 存在一正常数Cd，对于任意时间t，满足||d||≤Cd．假设2 系统的总未知干扰d对其所有自变量是连续可微的；存在一正常数Ch，对于任意时间t，使得d关于时间的导数||d˙ ||≤Ch[13]．假设3 h(υ,υd)关于υ满足全局李普希茨条件，即存在一个正常数Cf，使其满足||h(υ1,υd)-h(υ2,υd)||≤Cf||υ1-υ2||．η是可测的，为降低观测器的带宽和提高扰动的估计精度，设计MESO为η^˙=R(ψ)υ^+K1(η-η^);υ^˙=h(υ^,υd)+d^+B0τ+RT(ψ)K2(η-η^);d^˙=RT(ψ)K3(η-η^), (7)式中：(⋅^˙)为(⋅˙)的估计值；Ki(i=1,2,3)为增益矩阵；h(υ^,υd)=f(υd)+f'(υd)(υ^-υd)为模型已知信息．令观测误差ε=[η-η^,υ-υ^,d-d^]，将式(1)、(3)与式(7)做差，得观测误差动态方程为ε˙=A1ε+B1(h(υ,υd)-h(υ^,υd))+B2d˙，(8)式中：A1=-K1R(ψ)03×3-RT(ψ)K203×3I3×3-RT(ψ)K303×303×3;B1=03×3I3×303×3；B2=03×303×3I3×3．定理1 若假设1~假设3成立，则存在一个连续的李雅普诺夫(Lyapunov)函数V1=εTΓε(Γ为正定矩阵)和一个正常数C1，使得观测误差ε一致最终有界，最终界为C1，且增益矩阵的取值可使MESO闭环系统具有全局一致渐近稳定性．证明 A1是赫尔维茨的，必存在一个Γ和一个正常数C2，满足A1Γ+ΓA1T=-C2I．令λmax和λmin分别为Γ的最大和最小特征值，可得V1正定有界，即λmin||ε||2≤V1≤λmax||ε||2．对V1求导可得V˙1≤(-λmin+2Cf||Γ||)||ε||2+2||Γ|| ||ε||B2Ch．(9)由文献[14]可知，∃κ0，使得V˙1≤-ω1V1+ω2，(10)式中：ω1=(λmin-2Cf||Γ||-κ)/λmax，且ω10；ω2=||B2Ch||2||Γ||2/κ．对式(10)求解，两边开方可得ε(t)≤λmaxλmin||ε(0)||2e-ω1t+ω2[1-exp(-ω1t)]ω1λmin1/2, (11)ε是一致最终有界的，且存在一正常数C1=ω2/(ω1λmin)，使||ε(t)||≤C1．若观测器增益矩阵K1，K2，K3的取值满足Γω2ω1λmin||ε(0)||2, (12)则V˙10，即V˙1负定有界．||ε||→∞，V1→∞，定理1证毕．3 RC-MPC算法的设计3.1　观测误差鲁棒补偿控制算法设Am=F'(xd)，Bm=diag[B0,0]，Dm=[dT,0]T．若不考虑观测误差，令离散步长为T，由式(6)可得含有干扰的线性仿射离散时间系统e(k+1k)=Ake(kk)+Dk+Bku(k), (13)式中：Ak=I+TAm；Bk=TBm；Dk=TDm；(∙) (k+1k)表示基于k时刻的状态预测的k+1时刻的变量(∙)的值．若考虑观测误差，可将式(13)改写为e(k+1k)=e*(k+1k)+E^k+Bku'(k), (14)式中：e^=x^-xd；E^k为未知的观测误差项，由补偿控制输入u'对其进行补偿，使满足E^k+Bku'(k)≤0．(15)根据定理1，基于离线数据，可采用如下鲁棒补偿算法求解u'，即ℏ(C1)+Bku'(k)=0, (16)式中ℏ(∙)为E^k和C1的维度转换函数．3.2　鲁棒补偿的模型预测控制算法在船舶航行过程中，干扰无法被准确预测，若预测时域内均采用干扰仿射反馈控制，则将使约束被过度收缩，算法保守性增加[15]．基于以上现实及出于简化算法的目的，可做如下合理假设．假设4 xd和τd均为缓变的．假设5 预测时域内除k+1外的其他时刻系统退化为无观测误差与干扰项的名义系统，即e(k+i+1k)=A^ke(k+ik)+Bku(k+i), (17)此时，u(k+i)=u*(k+ik)，e(k+i)=e*(k+ik)．定义预测时域和控制时域均为N，由式(14)和(17)，可递推得到预测时域内各时刻的预测状态跟踪误差e(k+ik)．可用推力的平方表示功率[16]，为优化跟踪精度和船舶能耗，可设置目标函数J(k)=∑i=1N-1||e(k+ik)||Q2+∑i=0N-1||u(k+i)||R2+||e(k+Nk)||P2, (18)式中：Q和R分别为跟踪误差项与控制输入项的权重矩阵；||e(k+Nk)||P2为终端惩罚函数；P为终端跟踪误差项的权重矩阵．同时对状态和控制进行约束，将造成对终端状态的过强惩罚，从而导致可行域减小与优化问题不可行．为实现系统的递归可行性与优先保证预测终端的状态约束，可仅考虑e(k+Nk)[17]，即设置状态约束为e(k+Nk)∈Ω, (19)式中终端不变集Ω为原点的一个领域．船舶执行装置的运动范围和速率受限，须对控制力的幅值和变化量幅值进行约束．记控制力单位时间内变化量为Δτ(k)=τ(k+1)-τ(k)，因此系统须满足的控制约束为τmin≤τ(k+i-1)≤τmax；(20)Δτmin≤Δτ(k+i-1)≤Δτmax，(21)式中：τmin和τmax分别为推进装置能输出的控制力的最小值和最大值，Δτmin和Δτmax为控制力单位时间内变化量的最小值和最大值．这样，RC-MPC算法可转化为须考虑不等式约束，即式(19)、(20)和(21)的二次型最优值问题min  Ju*．(22)求解式(22)可得u*(k+i-1k)．由u(k)=u*(kk)+u'(k), (23)可得该时刻作用于船舶的总控制输入u(k)．4 级联系统稳定性分析自抗扰模型预测控制系统可视为级联系统，即∑1:e˙=g1(t,e,ε);∑2:ε˙=g2(t,ε), (24)式中：RC-MPC闭环系统∑1与MESO闭环系统∑2分别为该级联系统的驱动子系统与驱动系统；u*和u'分别为e和ε的反馈控制律．假设6 ∑1对各变量是连续可微的，且是局部李普希茨的；存在一个𝒦类函数ρ，使得∀||e||≥ρ(||ε||)0．假设7 RC-MPC闭环系统在k时刻存在最优解u⌢(k+i-1k)．定理2 在定理1和满足假设1~7的前提下，若e作为状态变量，ε作为输入变量，∑1具有输入-状态稳定性，则级联系统式(24)的原点是全局一致渐近稳定的[18]．证明 定义Jmin=minJ(e)为RC-MPC的最优值函数，Jf=J(e)为一可行值函数．选用V2=Jmin+V1作为候选李雅普诺夫函数．因(A^k,  Bk)可稳，必存在状态反馈u=Ke使得A^k+BkK渐近稳定，使得e最终趋于稳定[19]．函数具有连续性，由假设7可知，必存在Ω:={eeTPe≤α}使系统及其轨迹满足状态约束与控制约束，其中，α0，P0是P=Q-P+(A^k+BkK)ΤP(A^k+BkK)+KΤRK的解．由文献[20]可知：存在一个常数α1∈(0,α)，使Ω1:={eeTPe≤α1}∈Ω满足控制增量约束．在k+1时刻，选择一可行控制输入序列：[u⌢(k+1k),u⌢(k+2k),⋯,u⌢(k+N-1k), Ke⌢(k+Nk)]T;Jf(k+1)=Jmin(k)-||e⌢(k+1k)||Q2-||u⌢(kk)||R2 .则最优值函数满足Jmin(k+1)≤Jf(k+1)≤Jmin(k)．(25)由式(25)、假设6和定理1可得∀||e||≥ρ(||ε||)，使V˙2≤0．另定义𝒦∞函数：α1(e)=Jmin，α2(e)=Jmin+λmax||ρ-1(||e||)||2，使V2满足α1(||e||)≤V2≤α2(||e||)，可证系统∑1具有输入-状态稳定性．定理2证毕．5 仿真实验5.1　仿真参数设置以一艘652 t欠驱动船舶为研究对象，船舶模型矩阵系数为：m1=1.52×106；m2=1.47×106；m3=3.99×106；m4=3.51×106；m5=1.38×107；d1=1.84×104；d2=2.30×104；d3=-2.55×105；d4=3.80×105；d5=1.37×107；c1=6.52×105×(4.6r+v)；c2=-6.52×105×u；c3=105×(8.16vr+ 9.88rr)；c4=-105×8.66ur．多源时变环境载荷的参数为：Vc=0.2 m/s，βc=-π/9 rad；Vw=0.2 m/s；βwind=π/9 rad；Hwave=0.05 m；βwave=-π/2 rad．船舶的初始位置η0=[0 m，0 m，1.5 rad]T，初始速度υ0=[0 m/s，0 m/s，1.5 rad/s]T．路径点设置为[0 m，0 m]，[240 m，1 200 m]，[2 160 m，1 200 m]和[2 160 m，2 640 m]．期望路径的生成方式及引导律见文献[21-22]．其中，转弯圆弧半径均为500 m．船舶期望纵向速度ud=3 m/s．执行装置τmax=[424 kN,0,424 kN∙s]；τmin=-[424 kN,0,424 kN∙s]；Δτmax=[14 kN,0,40 kN∙s]；Δτmin=-[14 kN,0,40 kN∙s]．在以上仿真条件下，进行自抗扰模型预测控制算法与自抗扰控制算法的仿真对比实验．其中，自抗扰控制算法中线性误差反馈控制律为τXk=Kp(ud-u^k)-KdΔu^k;τNk=kp(ψd-ψ^k)-kdΔψ^k,式中：u^k为uk的估计值；ψ^k为ψk的估计值；Δu^k=u^k-u^k-1；Δψ^k=ψ^k-ψ^k-1；Kp，Kd，kp和kd为控制参数．MESO的增益矩阵均为：K1=diag[2,2,2]；K2=diag[1.3,1.3,1.3]；K3=diag[0.25,0.25，0.25]．在自抗扰模型预测控制算法中：N=4；Q=50000001.5×1010；R=20000001.4．自抗扰控制中的参数为：Kp=2×105；Kd=0；kp=1×106；kd=2×105．5.2　仿真结果分析图3为航行路径对比图．由图3可知两种方法能较准确地跟踪期望路径．图4为状态估计结果，图5为时变干扰估计结果．由图4和图5可知：MESO能在存在时变环境干扰下较准确地对状态和总未知干扰进行估计，并为控制算法提供可靠的输入．10.13245/j.hust.239269.F003图3航行路径对比图10.13245/j.hust.239269.F004图4状态估计结果10.13245/j.hust.239269.F005图5时变干扰估计结果图6为跟踪误差对比图．由图6可知：自抗扰模型预测控制的跟踪精度明显优于自抗扰控制，特别是当路径曲率发生变化时，该算法具有预测能力，能预先考虑未来变化的期望值，从而提高跟踪精度．10.13245/j.hust.239269.F006图6跟踪误差对比图1—自抗扰控制；2—自抗扰模型预测控制(下同)．图7为控制力对比图．由图7可知：二者均满足控制力幅值限制，但自抗扰模型预测控制的控制力变化更为平缓．图8为单位时间内控制力变化量对比图，图8(a)纵轴为单位时间内纵向控制力变化量，图8(b)纵轴为单位时间内转艏控制力矩变化量．由图8可知：自抗扰模型预测控制单位时间内控制力变化量均在其幅值上下限内，而自抗扰控制由于未对控制力变化率进行约束，因此控制力变化量超出幅值限制，将使得执行装置无法响应．10.13245/j.hust.239269.F007图7控制力对比图10.13245/j.hust.239269.F008图8单位时间内控制力变化量对比图根据文献[16]中的方式计算两种控制方法所消耗功率的数值，自抗扰模型预测控制功率消耗量优于自抗扰控制，这是目标函数式(18)对航行能耗进行显式优化的结果．6 结语本研究采用MESO对总未知干扰和状态进行估计，提高了观测精度，并将路径跟踪控制系统转化为含有干扰的线性仿射系统，降低了控制器的设计复杂度．为避免鲁棒模型预测控制造成的结果保守，本研究提出了一种RC-MPC算法，离线获得鲁棒补偿项，在保证系统跟踪精度与能耗最优性的同时，提高了系统对干扰的鲁棒性．该自抗扰模型预测控制算法实现了对船舶航行过程中未知干扰的估计和补偿，提高了船舶路径跟踪控制的精度，并减小了执行装置的负担和磨损，降低了航行能耗，有利于提高船舶航行的安全性和经济性．