喷水推进器利用喷水推进泵喷出高速水流的反作用力推动船舶前进，通过调节泵转速、方向舵舵角和倒车斗收放角可以实现推力大小和方向无级可调，具有抗空泡能力强、附体阻力小、保护性能好、噪声低、操纵性好等优点，近年来已被广泛应用于高性能船舶和各类军用舰艇[1]．相比于传统推进机构，采用喷水推进器可以在一定程度上提高动力定位系统的控制效果和定位能力．因此，开展基于喷水推进器的动力定位推力分配方法研究具有重要意义．Borrett等[2]构建了喷水推进器的推力模型，并开展了喷水推进器的控制系统及推力分配相关研究．高双[3]考虑转向机构和倒航机构对推力的影响，完善了喷水推进器的推力模型研究．Fossen等[4]将动力定位中的推力分配问题分为线性无约束问题、线性约束问题和非线性约束问题三大类．Johansen等[5]将功耗、推进器磨损和操纵性等作为综合优化目标，采用序列二次规划(sequential quadratic programming，SQP)方法进行求解．金超等[6]结合罚函数法和滤子法，对优化的目标函数进行了改进．李立国等[7]以某型起重船为对象，设计并进行了SQP法和遗传算法仿真对比实验．谢笑颖等[8]以某型双喷水推进船为对象，采用直接求解算法进行了数值仿真验证．对于动力定位船舶中的推力分配问题，传统的直接分配算法如推进器分组分配[9]，具有原理简单、实时性好、结果可行等优点而被广泛应用，然而该方法无法考虑复杂的非线性约束和功耗磨损目标函数，分配结果往往不能满足操作者的运动控制要求．对于非线性约束下的目标函数最优化问题，SQP是目前公认的最有效求解方法之一，但求解复杂的非线性约束优化问题时，迭代初值的选取不当或可行域的非凸都将导致算法收敛速度较慢且求解精度较低，甚至陷入局部最优，无法得到满足要求的解[10]．为了解决喷水推进器模型中复杂非线性约束下的推力分配问题，本研究提出一种区域二次规划算法，从优化初值和优化动态可行域两个方面提升算法的准确性和实时性．最后，以一艘装配喷水推进器的动力定位船舶为研究对象进行仿真实验，验证算法的有效性．1 推力分配问题描述1.1　问题描述推力分配实际上是求解一个最优化问题，输入是控制器输出的控制力和力矩指令，输出为各个推进器的推力与方位角指令，映射关系为以功率最小化等为优化目标，同时满足力的平衡、推进器推力变化率、角度变化率及推进器之间的水动力干扰等约束条件的优化过程．系统的工作原理如图1所示，控制器通过实际状态与期望状态的偏差计算出船舶运动所需控制力指令τ，推力分配单元将τ转化为各推进器的推力T和方位角α，各推进器的输出为喷泵转速n、方向舵舵角δ和倒车斗收放角β．10.13245/j.hust.249545.F001图1动力定位系统原理1.2　推力分配数学模型设动力定位船舶推进系统有p个推进器，随船坐标系为xoy，其原点通过船舶旋转中心，推进器i(i=1,2,...,p)的安装位置为(Lxi,Lyi)，如图2所示，图中：Ti为第i个推进器的推力；Txi和Tyi分别为该推力在纵向和横向上的分量．喷水推进器的推力可行域与文献[2]相同．10.13245/j.hust.249545.F002图2动力定位船舶示意图参考文献[11-15]，建立推力分配数学模型min    f(u,s,α)=uTWu+sTQs+ηε+det(BBT)+(α-α0)TG(α-α0);s.t.    equ   τ=Bu,       inequ   Ti∈Ω   (i=1,2,...,p), (1)式中：u为各推进器推力；W,Q和G为权值矩阵；s为松弛变量；ε为任意小的正常数；η为正的权衡因子；α为当前时刻推进器方位角；α0为上一时刻推进器方位角；uTWu为推进器的总功耗；sTQs为推进器分配输出与控制力指令之间的误差；η/(ε+det(BBT))为奇异惩罚项；(α-α0)TG(α-α0)为推进器磨损；τ=[X,Y,N]T为控制器输出的力和力矩；Ti为推进器i的推力矢量；B为推进器配置矩阵，B=cosα1sinα1lx1sinα1-ly1cosα1cosα2sinα2lx2sinα2-ly2cosα2⋮cosαpsinαplxpcosβp-lypcosαpT. (2)其中，lxi和lyi分别为第i个推进器在随船坐标系下的横纵坐标值；αi为第i个推进器的推力方位角．根据图2所示，可行域Ω可表示为Ω=Tx'≤g1Ty';Tx'≥g2Ty';Tx'≤g3Ty';Tx'≥g4Ty', (3)式中：g1~g4为图2所示四条边界曲线的对应函数；Tx'=(Tx/Tmax)×100%，Ty'=(Ty/Tmax)×100%分别为相应的归一化纵向和横向分量，Tx和Ty分别为喷水推进器推力在纵向和横向的分量，Tmax为正车时纵向最大推力值．1.3　喷水推进器模型考虑喷水推进器的转速指令调制时，通常采用等功率线形式的推进性能曲线来指导喷水推进器的操纵[16]．喷泵功率(Pe)与叶轮转速(n)、推进器推力(T)与航速(vs)的关系式如下Pe=C(n/1 000)3ηp；(4)T=ρwQ(vj-αsvs)=ρwQ(kv-as)vs，(5)式中：C为推进器功率系数；ηp为轴系效率；ρw为海水密度；Q为水泵流量；kv为喷速比；as为进流影响系数．为简化倒车斗收放角的影响，假设倒车斗只有正车(完全收起，倒斗角为0°)和倒车(完全放下，倒斗角为90°)两种工况[17]，方向舵舵角与方位角对应．选择Volvo Penta D4-150A系列喷泵发动机，最大转速为3 400 r/min．2 区域二次规划算法提出的区域二次规划算法拟从两方面提升算法的求解精度：一是对初值的优化，即通过预求解的方式得到一个满足约束条件的局部可行解，再将其作为初值代入二次规划求解算法中进行优化求解；二是采用区域分配策略，将喷水推进器的瞬时反向特性视做是动态可行域的空间不完全对称变化，通过区域分配的方式确定优化求解的动态可行域范围，同时保证其凸性．2.1　初值预求解根据上述推力分配数学模型，待分配变量为Ti=(αi,Ti)，为便于计算，引入扩展推力Ti=(Txi,Tyi)，Txi=Ticosαi，Tyi=Tisinαi．扩展型配置矩阵Be可表示为Be=1010⋯10010101-ly1lx1-ly2lx2-lyplxp. (6)采用伪逆法进行预求解，则有u0=Be+τ;Be+=We-1BeT(BeWe-1BeT)-1, (7)式中：u0为伪逆法求解结果；Be+分别为推进器扩展型配置矩阵的伪逆矩阵；We为扩展型权值矩阵．伪逆法求解过程不会影响算法的收敛速度，且求解的结果是满足等式约束和最小二乘目标函数的局部可行解[18]，将其作为迭代求解算法的初值能够加快算法的收敛速度，同时避免局部极值出现．2.2　区域分配策略进行喷水推进动力定位的动态推力分配时，还须考虑推进器的推力变化率ΔTlim、转角变化率Δαlim及正倒车状态变化时纵向推力的瞬时反向特性等，此时动态可行域为如图3所示非凸域，而传统SQP方法在求解该型非凸域内的优化问题时，会出现求解精度下降、收敛速度缓慢甚至陷入局部最优的情况．10.13245/j.hust.249545.F003图3喷水推进器动态可行域为了解决上述问题，采取区域分配的方式，首先将喷水推进器的推力可行域Ω根据正倒车状态划分为Ω1和Ω2两个子区域Ω1={Tx≥0   (Tx,Ty)∈Ω};Ω2={Tx0   (Tx,Ty)∈Ω}, (8)对图3中的动态可行域进行近似线性凸化处理，Θ'=sin(αtmax)Tx-cos(αtmax)Tx+2ωcos(αtmax)ΔTt≥0;       -sin(αtmin)Tx-cos(αtmin)Tx+2ωcos(αtmin)ΔTt≥0;cos(αt-1)Tx+sin(αt-1)Ty-Ttmin≥0;      -cos(αtmax)Tx-sin(αtmax)Ty+cos(Δαt)(Ttmax)≥0;       -cos(αtmin)Tx-sin(αtmin)Ty+cos(-Δαt)(Ttmax)≥0, (9)式中：Θ'为凸化处理后的动态可行域；αtmax=αt-1+Δαt；αtmin=αt-1-Δαt；Ttmax=Tt-1+ΔTt；Ttmin=Tt-1-ΔTt；αt-1和Tt-1分别为上一时刻的推进器方位角和推力值；t为相应的采样时刻；ω为任意小的正常数；Δαt和ΔTt为方位角和推力值增量，Δαt=Δαlim      (Tx≥Te),(90-αt-1+Δαlim)/2     (TxTe);ΔTt=ΔTlim        (Tx≥0),ΔTlim/2    (Tx0), (10)其中，Te为动态可行域变化的纵向推力阈值，Te=maxbtanΔαlimk-tanΔαlim,Ttmaxtan(Δαlim)；(11)k和b为推力可行域边界线性拟合系数．结合式(9)，设计如下区域分配逻辑Φ'=Θ'⋂Ω1(Tx0≥0);Θ'⋂Ω2(Tx00), (12)式中：Φ'为区域分配策略处理后的凸化求解域；Tx0为预求解的结果．2.3　组合偏置优化求解将式(7)计算出的结果作为迭代初值，结合式(12)所得的凸化求解域，基于前述推力分配模型进行二次规划优化求解，可以得到满足目标要求和推进器动态特性的最优分配结果，但在低速工况下，推力过小会导致喷水推进器工作效率和推力精度降低．为了解决上述问题，采用组合偏置的策略来避免推进器工作在推力较小的低效率区，同时可减小方位角的变化幅度，提高船舶的操纵性．偏置量计算公式如下[19]Δ'k=μFΔσ，式中：Δ'为偏置量；μ为偏置量系数；FΔ为偏置组内推进器推力上限的最小值；σ为自适应偏置因子，σ(t+1)=1(F(t+1)FΔ≤χ1;σ(k)χ1(F(t+1)FΔ≤χ2;0 χ2(F(t+1)FΔ,其中，F(t+1)为t+1时刻的推进器推力；χ1和χ2为阈值，0≤χ1≤χ2≤1．2.4　算法流程算法流程如图4所示．10.13245/j.hust.249545.F004图4区域二次规划法流程图3 仿真实验3.1　仿真条件为了进一步验证所设计方法的有效性和优越性，基于推力分配模型，以一艘652 t喷水推进动力定位船舶为对象，将所设计的推力分配单元集成至动力定位系统进行在线仿真实验验证．该船配备四台喷水推进器，船舶模型参数：船长69.9 m，质量652 t，惯性矩1.38×107 kg∙m2，重心纵向位置xg=4.6 m，水动力导数分别为Xu˙=-8.66×105，Xu=-1.84×104，Yv˙=-8.16×105，Yr˙=-9.88×105，Yv=-2.30×104，Yr=2.55×105，Nv˙=-5.08×105，Nr˙=-2.27×105，Nv=-3.80×105，Nr=-1.37×107．推进器参数如表1所示，表中：nmax为推进器转速上限；δmax为推进器方向舵舵角上限；βlim为推进器倒车斗收放角上限．环境参数：风速为2 m/s，风向角为20°，平均波浪幅值为0.1 m，浪向角为-90°，流速为0.2 m/s，流向角为-20°，环境力的计算均采用经验公式法[20]．10.13245/j.hust.249545.T001表1推进器参数设置编号nmax/(r∙min-1)δmax/(°)βlim/(°)位置/m13 4003090/0(-27.6，-2.5)23 4003090/0(-27.6，-1.05)33 4003090/0(-27.6，2.5)43 4003090/0(-27.6，1.05)仿真实验中的相关参数设置如下：最大迭代步数K=500，W=I4×4，We=I8×8，Q=diag[1×103,1×103,1×103]，G=diag[1×103,1×103，1×103]，奇异项小量ε=0.01，权衡因子η=1 000，推力可行域边界线性拟合系数k=9.4，b=2.94，近似线性项小量ω=1×10-5．仿真工况为低速循迹工况，期望航速为6 kn，期望路径为的光滑参数化曲线路径(0≤θ≤110)，表达式如下：xd=300sin(0.05θ)+15θ;yd=10θ.船舶初始状态为ς0=[30 m,0 m,1.53 rad]，v0=[0 m/s,0 m/s,0 rad/s]．航行过程中受突变载荷影响，纵向推力出现反向突变．在相同的仿真条件下，将区域二次规划法(RQP法)与SQP法[7]对比．3.2　仿真结果分析图5和图6为两种方法下轨迹和速度对比，图中u，v和r分别为船舶的纵向、横向和转艏速度．两种方法下船舶都能从初始位置运动至期望轨迹，然后按照期望轨迹航行，纵向速度最终也都能稳定在6 kn；相较于SQP法，RQP法下的船舶以更短的过渡航行轨迹进入期望路径上，且船舶各向自由度上的速度变化平滑．在完成动力定位循迹任务的前提下，RQP法更具有效性和优越性．10.13245/j.hust.249545.F005图5两种方法下的轨迹对比10.13245/j.hust.249545.F006图6两种方法下的速度对比图7为两种方法的分配结果对比，图中：τu，τv和τr分别为船舶纵向力、横向力和转艏力矩的分配结果．SQP法无法分配受超调和突变载荷影响下突变的控制力，分配结果产生明显的抖振现象；RQP法考虑了推进器动态特性，可以成功分配超调和突变载荷影响下突变的控制力．具有更高的求解精度．10.13245/j.hust.249545.F007图7两种方法的分配结果对比图8为两种方法下的推进器转速对比．在初始阶段和突变载荷影响阶段，SQP法下推进器转速频繁变化，且变化幅度较大，而RQP法下推进器转速变化平滑，避免了转速突变对动力系统的影响．10.13245/j.hust.249545.F008图8两种方法下的推进器转速对比图9为两种方法下的推进器倒车斗收放角对比．在初始阶段和突变载荷影响阶段，SQP法下推进器倒车斗正负向频繁切换，而RQP倒车斗状态无须频繁切换，较大程度降低了倒车斗的磨损．图10为两种方法下的方向舵舵角对比．在初始阶段和突变载荷影响阶段，SQP法下推进器方向舵舵角频繁变化，而RQP法下推进器方向舵舵角指令变化平滑，较大程度地降低了喷口转向机构的磨损．10.13245/j.hust.249545.F009图9两种方法下的推进器倒车斗收放角对比10.13245/j.hust.249545.F010图10两种方法下的方向舵舵角对比表2为两种方法的精度和速度量化对比．表中：MAE取三自由度上的最大绝对误差值，单次耗时取平均耗时．相较于SQP法，RQP法的平均绝对误差(MAE)降低至近乎为0 kN，单次耗时降低了约38.1%，求解精度和收敛速度得到较大的提升．10.13245/j.hust.249545.T002表2两种方法的精度和速度量化对比算法MAE/kN单次耗时/msSQP2.72×10227.6RQP8.52×10-1417.14 结论a．RQP法通过预求解和区域分配策略来处理喷水推进器的非线性约束并得到局部可行的迭代初值，比SQP法的求解精度更高，收敛速度更快，仿真实验也证明了该方法在分配突变控制力指令时的有效性和优越性．b．低速循迹工况下，采用组合偏置策略可有效解决推进器推力指令过小会导致推进器陷入低效率区的问题，并且能减少方位角的变化幅度．