基于监督与切换及群智能的推力分配算法研究

刘 明，华 亮

（1.南通大学 电气工程学院，江苏 南通 226019；2.南通大学杏林学院，江苏 南通 226236）

0 引 言

随着海洋产业的快速兴起，动力定位系统已经成为许多海洋工程船舶不可或缺的系统，特别是深海作业类船舶。动力定位船舶一般采用过驱系统，即将系统需求的合力指令合理地分配给各个推进器，故动力定位控制系统推力分配问题主要是约束优化问题，即在满足一定合力需求及推进器本身物理条件约束的同时，使功率消耗、机械磨损最小化，避免奇异性等。

目前国内外已有大量科研人员、学者对推力分配算法进行了深入研究，Winchers 等[1]用推进器分组法来解决推力分配的问题；摩根[2]最早提出了分配逻辑法来解决推力分配问题；序列二次规划法[3-4]是一种常见的推力分配策略，Johansen 等[5]提出了一种考虑避免奇异性的序列二次规划推力分配算法；Fossen[6]通过对推力分配算法的深入研究指出，伪逆法[7-9]和直接分配法是两种实时性比较强的推力分配算法，但必须解决存在的饱和问题；李新飞[10]把伪逆法和控制量归一化方法相结合解决了直接利用伪逆法存在的饱和问题；针对推力分配过程中尾流造成的推力损失问题，杨世知[11]和吴显法等[12]提出设置推力分配禁区方法进行解决；施小成等[13]、徐海洋等[14]提出了组合偏置推力分配算法。这些方法及其改进、衍生的方法在一定的条件下很好解决了推力分配问题，然而在综合考虑推力分配约束、优化条件时尚存在一些困难和问题。近年来，随着新兴的群智能算法的不断发展及在工程领域的应用，由于该类算法不依赖于对象问题的特征函数和解的形式，能有效解决多约束、非线性、多维度的复杂优化问题，故该类方法为推力分配问题的解决寻到了一条新的出路，目前已有部分学者如乔东生等[15]对该类算法在推力分配问题上的应用进行了初步研究，其本身具备的较大优化能力使得其在推力分配问题的处理上具有较好的应用前景。

1 船舶动力定位推力分配数学模型

1.1 推力分配问题约束条件和数学模型

对于动力定位船舶，推力优化分配目标是使控制器输出期望的力矩和控制力，除此之外推力分配还要考虑能量损耗、机械磨损、船舶操纵性等问题。推力分配问题约束条件可以分为期望控制力和力矩的等式约束条件、推进器本身物理限制的约束条件。

考虑所述的推力分配目标和约束条件，建立推力分配问题的数学模型，如下所示：

式（1）中，f∈为各个推进器产生的推力；α∈为各推进器的方位角；B（α）为推进器的位置构造矩阵；S为松弛变量，保证推力分配问题始终有可行解。式（2）中，（lxn,lyn）为第n个推进器在船体坐标系中的位置，αn为第n个推进器的方位角。

1.2“能量”最优推力分配问题的广义逆解

动力定位推力分配单元将推进系统需求的合力τc转化为各个推进器产生的力的大小和方位角，满足如下关系式：

推力分配等式约束是非线性方程，“能量”最优推力分配通常描述为如下形式的非线性推力分配优化问题：

为避免复杂的推力方程求解，Sørdalen 引入了扩展推力向量[16]概念，即将推进器的推力沿船体坐标系下X方向和Y方向分解成纵荡和横荡两个相互正交的力：

则式（3）可表示为

式中，

此时推力分配问题转化为线性等式约束最小化问题：

该线性等式约束最小化问题的广义逆解可以由拉格朗日乘数法求取：

式中，λ为拉格朗日算子。根据拉格朗日极小值存在条件得

式中，w为各推力的权重系数，且为正定对角矩阵。

2“双向穷尽”混群算法

2.1 算法概述

群体智能优化算法是指以模仿大自然中生物群体（而非生物个体）所能表现出来的分工、协同合作行为机制为目标的仿生智能优化算法，常用的群体智能优化算法有粒子群算法[17]、人工蜂群算法[18]、细菌觅食算法等。为了提高算法的寻优能力，一方面可以从算法自身机制上进行改进，另一方面可以借鉴其它算法的优点进行相互融合。本文提出的群体智能算法是融合粒子群算法、人工鱼群算法、果蝇算法、细菌觅食算法部分思想，以粒子群算法结构为主体，充分挖掘、发挥个体的寻优能力，保证群体的全局和局部寻优能力。一种群体智能优化算法的优劣体现在能否具备可靠的快速收敛性，而算法的快速收敛性在执行上主要体现为算法的全局和局部寻优能力，二者是相辅相成又互为矛盾的统一体。为了提高算法的全局寻优能力，本文算法引入了混沌算子、飞蝇算子和浓度算子；为了提高算法的局部寻优能力，本文算法一方面采取双向寻优、穷尽求取种群进化方式，另一方面每次迭代选取一定比例的最优个体进行局部寻优。

2.2 算法流程

算法执行的具体步骤描述如下：

（1）初始化种群，其中个体采用混沌算子进行相应初始化。

研究表明参数搜索空间和个体初始化分散性对寻优结果具有一定的影响。初始化种群，如果个体分布较均匀，分散性较好，则能为全局搜索多样性奠定基础。

混沌是非常普遍的一种非线性现象，具有随机性、遍历性等特性。运用混沌算子进行个体初始化比随机初始化更具优越性，能够避免算法陷入局部极值点，提高算法全局寻优能力。

（2）计算个体适应度值，计算全局最优值。

（3）采用双向寻优、穷尽求取法进行种群进化、个体更新。

设个体X={X1,X2,…,XN} ，每个个体最优解记做Pibest，全局最优解记做Pgbest。种群进化步骤如下：

①每个个体根据式（14）进行更新：

计算个体适应值并判别是否得到优化，是就执行步骤③，否则执行步骤②。

②每个个体根据式（15）进行更新：

计算个体适应值并判别是否得到优化，是就执行步骤④，否则执行步骤⑤。

③每个个体根据式（16）进行更新：

计算个体适应值并判别是否得到优化，是就根据式（17）更新c，执行步骤③，否则执行步骤⑦。

④每个个体根据式（18）进行更新：

计算个体适应值并判别是否得到优化，是就根据式（19）更新c，执行步骤④，否则执行步骤⑦。

⑤判断尝试次数是否超过设定值m，是则执行步骤⑥，否则执行步骤①。

⑥尝试m次后个体适应度值仍然不能得到改进且不是全局最优，则执行随机算子，如是全局最优则保留当前个体相应值。

⑦个体进化结束，更新全局最优值。

（4）选取p1个最优个体根据式（20）进行局部优化：

为了提高局部优化能力，式中v1的取值不能太大，每个个体根据式（20）以父代为基础衍生V1 个子代，计算个体适应度值，更新个体值和全局最优值。

（5）选取p2个最差个体根据式（21）进行局部优化：

为使个体跳出当前不利位置，式中v2的取值不能太小，每个个体根据式（18）以父代为基础衍生V2个子代，计算个体适应度值，更新个体值和全局最优值。

（6）运用浓度算子进行种群更新。

为了避免在迭代初期大量个体过于集中在某个区域，造成全局寻优能力减弱，在算法执行初期引入浓度算子（包括距离算子和适应度值算子）。

设第i个个体的值为Xi，适应度值为Yi，第j个个体的值为Xj，适应度值为Yj，则以第i个个体为中心，与第j个个体之间的距离定义为D，且D=‖Xi-Xj‖，同样适应度值之差定义为E，且

距离浓度算子定义为

式中，Dset为设定的常数。

适应度值算子定义为

式中，Eset为设定的常数。

则该算子的执行步骤如下：

①判断当前总迭代次数是否小于iter1，是则执行步骤③，初始化i=1，否则执行步骤②。

②判断当前总迭代次数是否小于iter2，是则执行④，初始化i=1，否则执行步骤⑤。

③运用距离浓度算子进行种群更新，根据式（22）和式（24）计算、判别第i个个体是否需要个体更新。

判别dcu（i）是否大于，是则根据式（25）进行第i个个体更新，否则第i个个体保持不变。

式中，Xpg为距离浓度低且适应值高的个体。

计算第i个个体的适应度值，并进行全局最优值更新，判别当前第i个个体是否为种群最后一个个体，如是则执行步骤⑤，否则i=i+1，执行步骤③。

④运用距离浓度算子和适应度值浓度算子进行种群更新，根据式（22）、式（23）和式（24）计算、判别第i个个体是否需要个体更新。

判别dfu(i)是否大于，是则根据式（26）进行第i个个体更新，否则第i个个体保持不变。

计算第i个个体的适应度值，并进行全局最优值更新，判别当前第i个个体是否为种群最后一个个体，如是则执行步骤⑤，否则i=i+1，执行步骤④。

⑤算法结束。

（7）判别迭代次数是否达到设定的最大值，是则结束，否则执行步骤（3）。

3 带监督与切换的推力分配

3.1 推力分配中监督与切换概述

对于动力定位船舶，一方面根据当前推力需求及权衡一些性能指标来确定推力分配方案，另一方面随着外界环境的改变又对推力分配有着不同的要求，针对特定的需求，利用监督与切换[19]，选出当前最合适的推力分配方案。

对于带监督与切换的推力分配器在进行推力分配时有两个问题应该考虑：首先是如何确定何种推力分配方案为最优方案；其次是当前推力分配事件至最优推力分配事件是否属于切换事件及切换事件能否发生。对于第一个问题，可以根据动力定位船舶当前测量的信息及推力分配数学模型中的部分或全部性能指标来决定哪种推力分配方案为当前最优的任务分配方案。对于第二个问题，一般采用停留时间切换技术和滞后切换技术，一方面避免不必要的切换，另一方面避免频繁的切换带来推进器的抖动和机械结构磨损的增加，甚至导致系统的不稳定、失控。

3.2 带监督与切换的推力优化分配

3.2.1 扇区约束问题

为了避免各个推进器之间的相互干扰，往往对各个推进器设置推力禁止产生的区域，即推力禁区，图1 为Cybership Ⅲ禁区设置示意图，图中阴影扇形区为对应推进器的推力禁区，由于采用群智能优化算法进行寻优，不需要考虑子扇区的凸与不凸问题，推力扇区组合一共被分为了4 个子扇区组合，具体如表1所示。

表1 扇区组合Tab.1 Sector division

图1 Cybership Ⅲ禁区设置示意图Fig.1 Schematic diagram of Cybership III forbidden zone setting

3.2.2 子扇区组合应用逻辑

推力禁区的存在使得相应的推力分配器的旋转角度不能连续变化，多个推进器的扇区被分成若干个子扇区的组合，在推力分配及其优化时有时需要进行扇区切换，即推力分配优化问题要解决子扇区组合应用逻辑问题。该问题中监督与切换功能的主要任务是在确保τc-B（α）f满足一定要求的情况下，降低推进器的总功耗，减少推进器的机械磨损，避免奇异现象。

带监督与切换的“双向穷尽”混群优化算法在推力分配优化中的应用主要通过如下几个步骤进行迭代求取：

步骤1：初始化。初始化参数包括：推进器推力大小及对应方位角、合力及合力矩大小、推力系数矩阵、每个推进器的推力极限值、每个推进器的推力最大变化率、每个推进器的方位角最大变化率、位置构造矩阵、功率与推力关系参数、式（6）齐次线性方程组的解等等，其中运用式（6）的齐次线性方程组的解进行初始化是为了能够确保获取合理的推力分配解，“双向穷尽”混群优化算法并没有选取推力及方位角作为个体，而是选取式（6）的齐次线性方程组的解系数c1、c2、…、ck为个体，这样以式（6）的齐次线性方程组的解为通解，以式（13）的最优“能量”为特解便可保证该算法能够收敛并获得合理的推力分配解，而且该个体选择方法有助于降低个体维度（例如本文实例个体维度可由6 维度降为2 维度）提高算法的实时性，对固定的船舶结构该解的解结构值是固定的，而特解值随需求合力及力矩大小改变而变化，并采取增量模式求解进一步提高算法实时收敛性。

步骤2：基于前一个时刻的推力分配解（fP，αP），运用“双向穷尽”混群优化算法（代价函数包括能耗、力与力矩偏差和奇异结构）搜索所有子扇区，最终得到最优推力分配的解（fbest，αbest）和最优推力分配的子扇区Ibest。

然而，如何运用“双向穷尽”混群优化算法获得最优推力分配的解和最优推力分配的子扇区呢？根据前面推力分配数学模型公式可知推力器由于自身的物理条件限制，当前时刻到下一个时刻Δt时间最大旋转角度为ΔαmaxΔt，即如果考虑角度变化率条件约束，则运用“双向穷尽”混群优化算法进行搜索只能在当前子扇区进行部分区域的寻优，不能对所有子扇区进行寻优，每组子扇区降低总的能耗效果的好坏就不能被立刻察觉到。此外，如果当前子扇区没有满足τc-B（α）f一定需求的解，运用“双向穷尽”混群优化算法搜索时，由于推力禁区的设置，使得推力器不能获得可靠的推力分配解，可能最终导致系统的不稳定。故运用“双向穷尽”混群优化算法获取最优推力分配的解和最优推力分配的子扇区时，其代价函数的设置不用考虑角度变化率约束条件，就能对整个子扇区进行搜索寻优。

步骤3：依据停留时间切换技术和滞后切换技术，综合考虑各种因素，确定一个最适合当前推力分配情况的扇区组合。

虽然步骤2 运用“双向穷尽”混群优化算法获取最优推力分配的解和最优推力分配的子扇区，但没有考虑角度变化率这个条件约束使得所获取的最优推力分配的子扇区未必是最适合当前推力分配情况的扇区组合。本文通过扇区切换逻辑来获取最适合当前推力分配情况的子扇区组合，该计算逻辑可以用式（27）进行表示。Ibest是全局寻优获得的最优子扇区组合；In是当前最优子扇区组合；IP是指In之前被选定的最优子扇区组合；t是指IP切换至In时到现在时刻为止的时间；tdts是指停留时间切换的最小切换时间；th是为了防止瞬态干扰导致的切换所设置的最大返回切换时间，其取值不大于tdts；JPlim是滞后切换时需要限制的代价函数变化最小值；S（·）是松弛变量累加和；Sdts是松弛变量累加和变化最大值，取值为百分数。

步骤4：根据选择的最适合子扇区组合实际情况，运用相应的优化算法，最终确定每个推进器的方位角和推力大小。

本文根据子扇区组合是否切换以及当前子扇区情况分为三种情况进行相应推进器的方位角和推力大小的求解：

（1）当In=inow且Ibest=inow时，inow为当前所处的扇区组合，运用“双向穷尽”混群优化算法（代价函数包括能耗、力与力矩偏差、方位角变化率和奇异结构）求解，根据所求最优角度解再结合全局最优角度解确定最终推进器的方位角，最后根据所求方位角利用双向穷尽”混群优化算法（代价函数包括能耗和力与力矩偏差）求得最终推力大小的解。

（2）当In=inow且Ibest≠inow时，或者In≠inow且inow不在禁区时，直接运用“双向穷尽”混群优化算法（代价函数包括能耗、力与力矩偏差、方位角变化率和奇异结构）进行推进器的方位角和推力大小的求解。

（3）其它情况时，首先确定推进器方位角的大小，然后在方位角已知情况下运用“双向穷尽”混群优化算法（代价函数包括能耗和力与力矩偏差）进行推进器推力大小的求解。

步骤5：输出推力分配结果，并判别此次推力分配是否是最后一次推力分配求解，是则结束，否则执行步骤2。

4 仿真实例分析

本文运用上述基于“双向穷尽”混群优化算法以及监督与切换机制来对动力定位船舶推力分配问题进行优化处理，并采用挪威科技大学仿真船模Cybership III 为对象进行仿真分析,推进器的部分参数如表2 所示，表中Kpt为功率与推力关系参数，此外两个吊舱式全回转推进器在1 s 时间内最大旋转角度为10°。

表2 推进器的部分参数Tab.2 Some parameters of thrusters

仿真实例一：

该实例的x轴推力需求在前100个采样周期始终为1 N，后100个采样周期除了在第121个采样周期出现突变力-2 N，其它都为2 N，y轴推力和系统力矩都为0，仿真结果如图2～5所示。

图2 各个推进器的推力大小Fig.2 Thrust of each thruster

由已知条件分析可知，如果不考虑奇异结构问题，该推力分配的理想结果为船艏推力分配器不作用，船尾两个推进器推力大小相等各占需求推力一半，且推进器方位角一致并为0°。由图2可知船艏推力分配器推力几乎为0，船尾两个推进器推力大小相等并接近需求推力一半，达到了很好的推力分配预期效果。由图3可以看出船尾两个推进器的方位角并不是保持一致为0°，但是几乎保持正负对称相等，这主要是因为目标函数考虑了奇异结构问题，该项权值的大小影响推进器的方位角大小。由图4和图5可以看出期望的推力和力矩几乎与实际的推力和力矩重合，它们的偏差几乎为0，推力分配效果良好。此外在第121个采样周期出现突变力情况下，由图2和图4可以看出系统能迅速反应，由于推进器本身物理约束条件限制，此时难以满足设定力和力矩存在一定偏差，但是推力迅速降为0，将偏差降为最小，且减少能量损耗，在第122个采样周期又能迅速达到稳定状态，可以避免干扰带来的影响。结合4 幅仿真图可以看出，当系统的最优扇区处于当前扇区时，所提算法能很好地在当前扇区进行推力分配，即使出现干扰使得干扰最优扇区不在当前扇区时系统也能很好地处理该问题，防止误操作的产生。

图3 各个推进器的方位角大小Fig.3 Azimuth of each thruster

图4 期望/实际的推力和力矩Fig.4 Expected/actual thrusts and moments

图5 期望/实际的推力和力矩的偏差Fig.5 Errors between expected and actual thrusts and moments

仿真实例二：

该实例的x轴推力需求在前100 个采样周期始终为1 N，后100 个采样周期始终为-1 N，y轴推力和系统力矩都为0，仿真结果如图6～9所示。

图6 各个推进器的推力大小Fig.6 Thrust of each thruster

由已知条件分析可知，该仿真实例需要船尾两个推进器从一个扇区组合切换到另一个扇区组合，而且都要经过禁区。由图6 和图7 可以看出，前100 个采样周期和仿真实例一运行结果一样，由于考虑奇异结构问题，船艏推力分配器不作用，船尾两个推进器推力大小几乎相等各占需求推力一半，推进器方位角几乎保持正负对称相等，在101 个采样周期时，期望推力发生巨变，此时最优扇区组合随即发生改变，图6 中船尾两个推进器推力大小先是迅速减少至0，然后增加并保持一段时间后再迅速增加，然后又平滑下降直至稳定，这是由于采用了停留时间切换技术。从图7 可以看出，船尾两个推进器的方位角先是以最快速度通过禁区，达到新的最优扇区组合后方位角相对平稳达到新的平衡态。由图8和图9可以看出除了切换过渡过程中外，其它时间段的期望的推力和力矩与实际的推力和力矩几乎都是重合的，它们的偏差几乎为0，用了大约10 个采样周期实现了大区域的切换，可见达到了预期的目的，推力分配效果理想。由图6 和图8 可以看出在最优扇区切换过程中，由于推进器本身物理约束条件限制，推进器推力变化出现了突变，力和力矩的设定值和实际值存在一定偏差，而运用本文算法可以使得这些不利因素降低。结合4 幅仿真图可以看出，当系统的最优扇区不处于当前扇区而需要切换时，所提算法能很好并快速地进行切换并保证推力分配最优化。

图7 各个推进器的方位角大小Fig.7 Azimuth of each thruster

图8 期望/实际的推力和力矩Fig.8 Expected/actual thrusts and moments

图9 期望/实际的推力和力矩的偏差Fig.9 Errors between expected and actual thrusts and moments

仿真实例三：

前面两个实例对所提推力分配算法进行了“静态”测试，达到了预期的效果。为了更好说明该算法的可行性，本实例将对所提推力分配算法进行“动态”测试。该实例推力需求在前100个采样周期：x轴1 N 左右随机产生，y轴始终为0.5 N，力矩为0；后100 个采样周期：x轴-1 N 左右随机产生，y轴在0.5 N左右随机产生，力矩为0。仿真结果如图10～13所示。

图10 各个推进器的推力大小Fig.10 Thrust of each thruster

由图10可以看出在前100个采用周期，由于y轴存在固定期望推力，故船艏推进器作用，船尾两个推进器随着x轴的变化能同步改变，后100 个采用周期，由于y轴期望推力发生了随机改变，此时船尾两个推进器推力不再同步改变，而且船艏推进器不再是一个定值，但所有推进器的推力都没有发生失控现象，且变化幅度除了切换时间时都比较小。由图11 可以看出，整个过程中船尾两个推进器的方位角变化都很平稳，说明该算法能尽量减少推进器的机械磨损。由图12 和图13 可以看出，尽管给定的推力发生了随机改变，但是期望的推力和力矩几乎与实际的推力和力矩除了切换过渡过程中其它时间段几乎都是重合的，它们的偏差几乎为0，同样用了大约10 个采样周期实现了大区域的切换，达到了预期的目的，推力分配效果理想。结合4幅仿真图可以看出，所提算法能实现推力分配的动态最优化。

图11 各个推进器的方位角大小Fig.11 Azimuth of each thruster

图12 期望/实际的推力和力矩Fig.12 Expected/actual thrusts and moments

图13 期望/实际的推力和力矩的偏差Fig.13 Errors between expected and actual thrusts and moments

通过以上实例分析可知，本文所提出的推力分配策略在处理推进器推力极限、推力变化率极限、方位角变化率极限、机械磨损、禁区限制，奇异性等约束优化问题上相对二次规划法和伪逆法等推力分配算法体现出了较大的优势和灵活性，然而这些优势的获取是以相对牺牲实时性为前提的，针对实例3，通过20次仿真运行，结果显示单位采样周期最大运行时间为0.5633 s，最小运行时间为0.4977 s，平均运行时间为0.5153 s，可见满足一定的实时性要求，然而实时性与二次规划法和伪逆法等推力分配算法相比仍然处于劣势，这也是由本文算法本身结构特点决定的。虽然实时性上相比较而言处于劣势，然而可以获得更加灵活、更好的其它性能，使得该算法具备更好的应用前景，而且随着量子计算技术的成熟，该劣势会逐渐降低。

5 结 语

本文研究的基于“双向穷尽”混群优化算法以及监督与切换机制是一种新型的动力定位船舶推力分配求解方法，其中“双向穷尽”混群优化算法融合了多种群智能算法的优势，使得该算法在大范围搜索空间内局部及全局寻优能力和收敛性更强，把该算法和监督与切换机制相结合使船舶推力分配的全局最优解更容易获得。本文船模Cybership III 实例“静态”和“动态”分析也表明该算法在推力分配过程中能取得较好效果，具有良好的应用前景和推广价值。