基于遗传算法的半潜式平台动力定位系统动态约束可行域推力分配法

2018-11-02梁海志乔东生李芦钰张纪刚欧进萍

船舶力学 2018年10期

梁海志，乔东生，李芦钰，张纪刚，欧进萍

（1.大连理工大学海岸和近海工程国家重点实验室，辽宁大连 116024；2.青岛理工大学土木工程学院，山东青岛266033）

0 引言

动力定位系统最早在钻井船等大型船舶上得到广泛应用，随着海洋石油开采业向深海进军，动力定位系统以其不受作业水深影响等优势逐渐取代传统锚泊系统被应用到深海海洋平台上。在风浪流联合作用下，动力定位推力系统提供控制力使半潜式平台位移响应保持在安全范围内，考虑到推力器出力方式及运行经济性，仅控制纵荡、横荡和艏摇三个自由度的低频慢漂运动，而推力系统由8个推力器组成以保证系统的冗余度，因此整个系统为过驱动系统。在动力定位模块化设计过程中，如图1所示，首先由控制算法得到三自由度控制力，然后通过分配算法将其分配到底层执行机构推力器上，最终实现定位目标。

动力定位系统的底层推力系统一般采用全回转推力器、舵等，在推力器方位角已知的情况

下，Sørdalen[1]采用伪逆法对推力直接进行分配，分配过程中存在的奇异构型问题，通过修改奇异值来解决。付海军等[2]提出了基于自适应遗传算法的加权伪逆策略，改善推进器角度饱和问题。伪逆法的特点是快捷有效，但是无法给出最优推力分配方案。针对分配过程中的奇异构型问题，徐胜文等[3]提出了一种高效地避免奇异结构的简单方法。序列二次规划法是推力分配策略中采用较为普遍的优化方法，如文献[4-6]等。推力分配过程中，影响动力定位系统性能及其运行经济性的主要因素为推力器间尾流相互干扰造成的推力损失，针对该问题，文献[4-8]提出了对推力器推力方向设置推力禁区的方法，使推力器间的相互影响降到最低，但由此得到的推力可行域为非凸集，无法直接采用序列二次规划法，处理方法较为复杂。

本文针对以上存在的问题，基于文献[5]建立了以推力变化率和转角变化速率为控制参数，以耗能最小为目标的优化模型。在考虑推力器推力损失的基础上通过动态设置推力可行域，利用遗传算法对控制力进行分配，在提高运行经济性的同时也避免了序列二次规划法对可行域要求严格的问题。

图1 动力定位系统工作流程Fig.1 Workflow diagram of DPS

1 推力分配模型

为使动力定位系统能准确执行顶层控制策略得到的控制力，推力器推力和控制力之间需满足如下关系

式中，τ包括纵荡、横荡控制力和艏摇控制力矩，T为推力器推力，T=［T1…Tm］T，m为底层推力系统推力器个数，B=［B（ α1）…B（ αm）］，［lxi,lyi］为第i个推力器在图2所示随体坐标下的布置坐标，αi为第i个推力器转角，Ti为第i个推力器的推力值，如图3所示。

推力分配策略在满足控制力要求的前提下，更应使整个系统燃油消耗最小以提高运行经济性，系统耗能功率为

式中：KQ、KT分别为螺旋桨扭矩系数和推力系数，ρ为海水密度，D为螺旋桨直径。针对推力分配系统，还需进一步考虑推力方向变化速率，由此得到动力定位系统推力分配策略的优化目标

图2 推力器布置示意图Fig.2 Propeller layout diagram

图3 推力器推力及推力方向示意图Fig.3 Propeller’s thrust and thrust direction

约束条件：

（3）式中，T∈Rm×1，Δα∈Rm×1为优化参数；Ω 为推力方向变化权矩阵，Ω∈Rm×m。（3）、（4）式建立了以推力T和推力方向变化Δα为变量的推力分配基本数学模型，但是对于推力器而言，受其自身动力性能限制，推力器很难实时响应由上述模型得到的推力T，若以推力变化ΔT作为优化变量将更满足客观

现实，为此采用ΔT替换推力T对上述模型进行改进，对（2）式系统耗能功率进行泰勒展开，可得

式中：Ti0为t0时刻第i号推力器推力，将上式代入（3）式，则得到以推力变化ΔT和推力方向变化Δα为变量的优化目标方程

同时，约束条件式(4)中的等式约束方程相应的线性化为

推力变化ΔT和转角变化Δα约束条件为

由此得到如图4所示的推力器运行可行域，并可通过α=α0+Δα，T=T0+ΔT得到推力器出力方向和推力，最后还需考虑推力幅值约束条件

式中：Tmin、Tmax为推力器的推力上下限。

图4 单个推力器可行域Fig.4 Schematic diagram of single propeller’s thrust feasible region

综上所述，建立了以（6）式为优化目标方程，u= ［ΔT1×mΔα1×m］T为优化变量，（7）、（8）式和（9）式为约束方程的优化模型。但是对于半潜式平台动力定位系统而言，从图2中可看出，同一浮箱上推力器间安装距离较近，受推力器间尾流的相互影响，势必会造成推力器推力损失问题，若不加以考虑，不仅会增加运行成本，更会导致定位失位问题，本文从约束推力方向入手，最大程度地降低因桨—桨干扰造成的推力损失问题，提高系统可靠性及经济性。

2 推力方向动态约束域设置

2.1 推力损失模型

推力分配策略的首要任务是使各推力器协调作用以满足控制力要求，整个系统中影响定位性能的除推力器自身动力性能外，还有推力器之间的相互干扰造成推力损失，致使推力器实际推力达不到指令要求造成定位失位。

针对推力损失问题，Ekstrom等[9]进行了桨-桨干扰敞水试验，从试验结果中可看出尾流在相当大的范围内会对下游螺旋桨产生影响。Nienhuis[10]给出了平板下桨-桨干扰试验结果。为了能在实际工程中对下游螺旋桨推力减额进行定量估算，Dang等[11]分析了桨间距和推力方向对推力损失的影响，并总结出了桨—桨干扰推力损失经验公式，前后推力器同向布置，桨间距变化时推力损失经验公式为

式中：t为推力减额因数，T0为敞水系柱推力，T为下游螺旋桨所产生的推力，x为两螺旋桨间距离。进一步，固定桨间距，推力器间相对推力方向变化时的推力损失经验公式为

式中：Φ为两推力器之间的相对夹角，t为x/D为定值、Φ=0时的推力减额系数。

2.2 推力可行域设置

在动力定位系统中，若不对推力器推力方向加以限制，总会出现如图5所示的情况，下游推力器完全处于上游推力器的尾流中，由推力损失模式经验公式可知，此时下游推力器推力损失最严重。文献[8]针对推力损失问题总结了常用的3种处理方法：（1）接受推力损失，通过增加推力转速来弥补推力损失；（2）设置推力禁区，如图6（a）阴影部分所示，对推力方向不施加约束，当推力方向在该禁区时设置推力为零，；（3）设置推力禁区，始终保持推力器出力方向处于禁区外。方法（1）最易实现，但会明显增加燃油消耗；方法（2）中，当1号推力器推力为零时（如图6（a）所示），由于整体系统的推力突然减小从而导致短时失位问题，并且若其余推力器已满负荷工作，那么由于无法提供足够的推力亦会导致定位失位；方法（3）中如果分配指令要求1号推力器的推力方向从a快速到达b，那么由于禁止角的存在，将使整个系统耗能无法在最短时间内实现最优。基于以上方法的缺点，本文依据推力损失模型，通过动态设置下游推力器推力可行域来实现推力分配。

图6 推力器推力可行域示意图Fig.6 Schematic diagram of thrust feasible region

针对如图2所示定位系统，本文将相临的两个推力器分为一组，如将1号和2号分为一组，首先，对每组中的各推力器推力方向按图6（a）进行分区，阴影区域为桨-桨干扰最严重范围；其次，若上游推力器推力方向进入阴影区域，启动推力方向约束条件以避免桨-桨干扰，下游推力器推力方向可行域根据主推力器的推力角进行设置，下面针对第1组推力器介绍其可行域的具体设置方法。

如图6（b）所示，当1号推力器位于上游时，其推力、推力角可行域仅需满足推力器物理性能即可，其可行域如图4所示，图中Δα±为转角变化范围，ΔT±为推力变化范围；对于2号推力器，除需满足物理性能外，还应通过设置两推力器之间的夹角以避免受1号推力器尾流影响而造成的推力损失，因此当1号推力器进入阴影区域时，将2号推力器可行域设置为图6（b）中所示，其约束方程为

综合（8）、（9）、（12）式及（13）式，推力器可行域归纳为

（14）式建立了第一组推力器的约束条件，其余组可依此设计。若同组推力器未发生推力干扰问题，则（14b）式约束条件将不被考虑。（14）式根据各个推力器的实时推力方向来决定其可行域，从而实现了动态约束。

3 数值模拟分析

本文以某半潜式平台动力定位系统为研究对象，该平台总长114 m，总宽90 m，总高110 m，最大甲板可变载荷9 000 t，平台作业吃水19.0 m，作业水深为1 500 m。波浪载荷采用Pierson-Moskowitz谱，有义波高为4.0 m，跨零周期为7.9 s，入射方向由45°到60°，利用模糊PID控制算法[15]得到控制力和控制力矩时程，如图7所示，0-1 500 s波浪力入射方向为45°，1 500-3 000 s入射方向为60°。选择W-B4-70导管螺旋桨作为推力器，各个螺旋桨额定功率为4 600 kW，螺旋桨的主要参数采用文献[16]所给值，分别为：系柱状态下推力系数KT0=0.445，扭矩系数KQ0=0.06，直径D=4 m，Tmax=540 kN，Tmin=0 kN，螺旋桨转速变化率为±0.2 s-2，则推力变化率范围为 ΔT=±4.67 kN/s，转角变化率范围为 Δα=±1°/s。本文在对螺旋桨进行布置时选择了错位布置，如图2所示，错位布置的主要目的是为了加大螺旋桨之间的距离，降低桨-桨干扰，表1为螺旋桨布置坐标值。角度权矩阵Ω=10×I8×8。采用Matlab提供的遗传算法工具箱，染色体用双精度向量编码，种群数为20，遗传代数100代；适应度函数为（6）式；采用随机均匀选择算子；选用Scattered法交叉，交叉概率0.2；变异过程采用默认函数。

表1 螺旋桨坐标值Tab.1 Propeller coordinates

图7 控制力、控制力矩时程曲线Fig.7 Time-histories of DP control force

如图6（a）所示，1号和2号推力器中心连线与x轴的夹角δ=55.5°。根据文献[13]可知，在系柱状态下离上游推力器约3倍直径时，尾流速度分布为沿轴线最大的单峰状态，尾流影响范围直径约为2D，基于此本文取θ=36°，即相临推力器间理想的最小出力夹角为18°，此时下游推力器推力减额因数由（11）式计算得到为 0.882 5。1号推力器的禁止角为[37.5°,73.5°]，同理，2号推力器的禁止角为[217.5°,253.5°]，其余三组推力器均依此设置。

图8 推力方向时程曲线Fig.8 Time-histories of thrust directions

图9 推力时程曲线Fig.9 Thrust time-histories

图11 推力变化du时程曲线Fig.11 Time-histories of thrust changes du

图8和图9分别为ψ=18°时的推力角时程曲线和推力时程曲线。图8（a）为第一组推力器的推力角时程曲线，图中斜线所示区域为1号推力器禁止角，1号推力器进入该区时成为主推力器，此时2号推力器在依据1号推力器出力角设置的可行域内出力。在设置次推力器可行域时，充分考虑了推力器的物理性能，即其受转角转动速率的限制，因此在短时间内无法完全实现避免推力器间的推力损失问题，从图8（a）的700-950 s区间，以及图8（b）和8（c）中均可看出。当主推力器在禁止区作用时间较长时，推力器的转角得到充分转动，从图8（a）的1 800 s到2 550 s间可看出两推力器存在明显的夹角并在约1 900 s处达到理想状态，此时下游推力器推力损失降低到了最小。当推力器均未处于推力禁止角时，两者出力角基本一致。图9（a）为第一组推力器推力时程曲线，由于2号推力器出力角受限，此时1号推力器出力较大。从图中可看出，第三组情况与第一组基本相同。图8（b）中，在约850 s以后推力角大于2π，主要原因是推力器出力角范围未受限制，其出力角以2π为周期。为了方便比较，对3号推力器的推力角进行取余，使其推力角落在[-2π,2π]内，修正后可看出，当该组推力器均未处于禁止角时出力方向基本一致。

为了保证动力定位系统的安全性，每个推力器推力需留出20%的安全裕度，从图9中可看出推力均未超过432 kN。图10和图11分别为推力角变化da和推力变化du的时程曲线，从图中可以看出推力变化未超出[-4.67 kN/s 4.67 kN/s]范围，转角变化率亦满足[-1°/s 1°/s]。图10和图11说明该分配方法满足推力器的物理性能。