基于改进SHADE 算法的船舶电力系统推力分配

朱昱雯，肖健梅，王锡淮

上海海事大学物流工程学院，上海201306

0 引 言

随着航运业的不断发展，人们对海洋资源的依赖程度越来越高，传统的锚泊定位系统已经无法满足船舶的深海作业需求。目前，很多船舶都配置了动力定位系统，可以精准定位船舶，从而控制船舶在期望方位上航行。作为动力定位系统的重要执行模块，推力分配模块可以根据动力定位控制器给出的期望推力和力矩指令，对船舶推进器的推力和力矩进行合理分配，同时抵消环境扰动，从而使船舶在期望的艏向和位置上航行。此外，通过合理分配各个推进器的推力和方向角，可以解决推进器过度旋转导致的磨损问题和推力损失引起的船舶电力系统功率波动问题。

在保持期望艏向和位置的前提下，推力分配模块还需保证推进器产生的功率是最小推力组合，因此，推力分配问题可以视为具有约束条件的单目标优化问题。目前，大多采用序列二次规划方法来解决推力分配问题，但该方法较为繁琐且不利于编程实现，同时受约束条件的影响较大，故其实际可应用的范围非常有限［1］。徐云成等［1］利用了序列二次规划来解决推力分配问题，结果表明：推进器推力和方向角的变化较为平缓，且误差较小；但该研究未结合推力分配模块与动力定位控制器进行实时仿真，即仿真系统没有根据上一时刻船舶的位置和航速来改变推力力矩，而是仅对已知的推力力矩指令进行推力分配测试，故其应用范围较窄。而智能算法可以很好地解决推力分配问题，其过程也比序列二次规划方法更简明。目前，适用于推力分配问题的智能算法主要有遗传算法、人工鱼群算法和粒子群算法等［2］。Zhao 等［3］将遗传算法应用于推力分配问题，结果表明，遗传算法可以响应动力定位控制器根据上一时刻船舶位置而产生的期望指令，但其推进器推力和方向角的变化存在震荡，所以增加了推进器的磨损。刘明等［4］结合伪逆法与混沌粒子群算法来解决推力分配问题，结果表明，推进器推力和方向角的前期变化存在波动，但后期变化较为平缓。虽然智能算法在处理速度、简便程度、可实施性方面均优于序列二次规划方法，但在算法后期易陷入局部最优问题，从而导致期望推力与实际推力的力矩指令之间存在一定误差。

因此，本文拟采用基于成功历史记录的自适应参数差分进化（success-history based parameter adaptation for differential evolution，SHADE）改进算法来解决船舶的推力分配问题。差分进化算法在前期收敛速度较快，后期因种群多样性减少而过早停滞，传统SHADE 算法已针对这一问题进行了改进。然而，对于存在环境扰动的推力分配问题，SHADE 算法的误差较大且响应速度较慢，故本文拟基于传统SHADE 算法进行改进：根据当前迭代次数，在动态范围内选择较优的个体进行变异操作，以提高后期的收敛速度；根据算法的迭代程度，对用以产生交叉率CR 的正态分布方差和用以动态调整位置参数的权重c 进行自适应优化，从而提高算法的稳定性。最后，通过与传统SHADE算法进行对比，证明改进SHADE 算法的可行性和有效性。

1 动力定位数学模型

典型的动力定位控制系统结构如图1 所示，本文将采用PID 控制器作为动力定位系统的控制器，用以控制船舶纵荡、横荡和艏摇3 个自由度的运动。首先，PID 控制器通过对比实际位置η、实际航速v 和参考位置ηd，从而产生期望推力力矩指令τd；然后，基于改进SHADE 算法，即可得到推进器总功率最小时各个推进器的推力f 和方位角α 的最优组合；最后，将该组合分配给各个推进器以产生实际合力τ ，根据动力学模型估算船舶的当前位置和速度信息并反馈至动力定位控制器，从而实时控制船舶定位。

图1 控制系统结构图Fig.1 Structure diagram of control system

PID 控制器的控制律为［5］

式中：RT(ψ)为旋转矩阵的转置，其中ψ 为船舶的艏摇；Kp，Ki，Kd均为PID 控制器的参数；s 为拉氏变换中的复变量。

1.1 船舶数学模型

船舶的动力定位系统通常仅考虑纵荡、横荡和艏摇3 个自由度，这类船舶在北东坐标系NED下的位置和船体坐标系xyz 下的速度关系为［6］

1.2 推力分配数学模型

1.2.1 推进器数学模型

其中

1.2.2 推力分配目标函数

推力分配问题不仅需考虑每个推进器产生的推力和方向角是否满足动力定位控制器的期望推力和力矩指令，还需考虑推进器推力和方向角的变化率。在每次更新期望推力和力矩指令之后，推力分配模块均需快速响应每个推进器需产生的推力和方向角，但推进器在短时间内的推力增加量有限。同时，全回转推进器的过度旋转将导致推进器磨损和功耗增加。此外，当多个推进器同时工作时，相邻推进器将互相影响，一个推进器螺旋桨产生的尾流会进入相邻推进器的螺旋桨工作区域，对相邻推进器产生水动力干扰，从而导致推进器的推力损失，因此需限制全回转推进器的方向角，从而避免推力损失并减缓推进器方向角的变化速率［8］。

推力分配的目标函数为

约束条件为

2 改进SHADE 算法

SHADE 算法是一种利用前期成功试验向量的控制参数来调整当前迭代控制参数的改进差分进化算法，其计算思路为：首先，利用存储空间记录多组可生成较优控制参数的均值参数；然后，利用这些均值参数来生成每次迭代过程的控制参数，并在每次迭代结束之后利用生成成功试验向量的控制参数来调整存储空间的均值参数。这种自适应调整控制参数的方法不仅加快了差分进化算法的收敛速度，还提高了鲁棒性。

SHADE 算法在解决无约束问题时的表现较好，但处理复杂问题时将过于依赖随机生成的初始种群。例如，对于推力分配问题，如果动力定位控制器产生的期望指令过小，则SHADE 算法的鲁棒性将有所下降，且易出现过早停滞的现象，从而导致偏差较大。针对这些问题，本文将对传统SHADE 算法进行改进。

2.1 算法原理

2.1.1 变异操作

式中：g 为当前迭代次数，对应第g 代种群；G 为最大迭代次数。

基于式（9），SHADE 算法的后期迭代过程可以减少在非精英区域搜索的范围，从而加快后期的收敛速度并避免陷入局部最优。本文将参数p的最小值设定为round(0.025×NP) ，当种群规模NP=200 时，参数p 随迭代次数的变化曲线如图2 所示。

图2 当G=500，NP=200 时，参数p 随迭代次数的变化曲线Fig.2 The curve of parameter p with the number of iterations at G=500，NP=200

2.1.2 交叉操作

由于推力分配优化问题含有2 个变量：推进器推力f 和方向角α，因此需对2 个变量均进行交叉操作。本文将采用二项式交叉方式对目标向量和变异向量进行交叉操作，具体如下：

通过对每个变量向量中的至少一维分量进行交叉操作，即可增加推进器推力和方向角的组合方式，用以保持种群的多样性，从而提高算法的前期收敛速度。

2.1.3 选择操作

分别计算试验向量Ui，g和目标向量Xi，g的适应度值，从而选取可以产生更优适应度值的个体进入下一代，其选择操作如下：

2.1.4 参数自适应调整

由于差分进化算法依赖于控制参数（例如，缩放因子F 和交叉率CR）的设置，并且动力定位控制器每一次更新期望推力和力矩指令之后的工况均有所不同，因此，如果人为设置固定的控制参数，则将导致船舶航行过程中某一次期望推力及力矩指令与实际推力及力矩指令之间的误差过大。虽然在SHADE 算法中已采用自适应的方式对控制参数进行了调整，但用于生成控制参数分布的某些参数仍依赖于人为设定，因此，为了避免该误差导致船舶偏离航线，本文将对这部分人为设定参数进行自适应调整。

本文将分别采用柯西分布和正态分布随机生成每个个体的缩放因子Fi，g和交叉率，即

式中：Fi，g∈［0.2，0.8］；CRi，g∈［0，1］；randc 为柯西分布函数；randn 为正态分布函数；μF为柯西分布的位置参数；μCR为正态分布的均值；σ2CR为正态分布的方差。

如果Fi，g和CRi，g超出阈值范围，则将采用截断的方式对其进行限制［11］。如果f（Ui，g）＜f（Xi，g），即表示生成了成功试验向量Ui，g，其对应的Fi，g和CRi，g将被分别存储至空间SF和SCR中，用以在每一代结束后动态调整μF和μCR。同时，调整之后的μF和μCR将存储至空间W 中，用以在每一个目标向量变异之前从W 中随机选取μF和μCR，其中空间W 的大小为

将空间W 中μF和μCR的初始值分别设为0.1和0.5，并在算法初始化时设定一个常数k，用以控制μF和μCR在空间W 中的存储位置；同时，利用μF和μCR更新空间W 中第k 个位置的数值［12］，当k ＞Wsize时，令k=1，重新开始计数。

针对柯西分布中的位置参数μF，本文基于SHADE 算法对其更新方式增加了一个权重参数c，用以控制μF的更新速率。虽然空间W 的存在可以预防μF更新过快所导致的算法陷入局部最优的问题，但空间W 的大小并不容易控制，不同的期望推力指令需要不同大小的空间W，才能达到期望范围内的偏差。当船舶需要稳定在某一点时，过大的空间W 将导致算法迭代缓慢；但在船舶航行过程中，过小的空间W 将容易导致算法陷入局部最优。因此，本文将通过新增的权重参数c来进行平衡，从而减小空间W 对船舶航行过程中各种工况下推力分配的影响。位置参数μF的动态调整方法具体如下：

对于控制μF更新的权重参数c，本文将根据迭代次数进行动态调整，即在算法前期更新速度较慢以保证算法的稳定性，而在后期则将提高更新速度以满足精度需求。权重c 的更新方式为

正态分布虽与柯西分布类似，但其存在均值和方差，故可较好地反映交叉率CR 的变化情况，因此，生成成功试验向量的交叉率CR 可以引导均值和方差变化，从而自适应调整交叉率CR。但是，如果在均值μCR的动态调整过程中增加了权重c，反而将降低最优均值μCR的寻找速度。因此，本文将在正态分布的参数调整过程取消权重c，则均值μCR和方差σ2CR的动态调整方式如下［13］：

2.2 改进SHADE 算法的流程及步骤

在基本差分进化算法最优值的基础上，改进SHADE 算法增加了存储历史参数的空间和利用历史参数调整控制参数的步骤，并在变异操作中增加了随机选择的精英个体。利用改进SHADE算法计算最优推进器推力和方向角的流程如图3所示。

图3 改进SHADE 算法程序框图Fig.3 Block diagram of improved SHADE algorithm program

采用改进SHADE 算法对推力分配问题进行优化的具体步骤如下：

1）根据动力定位控制器给出的期望推力和力矩指令τd以及推进器在前一时刻的推力和方向角，在相应的可行解空间中产生种群规模为NP的初始种群（包含每个推进器的推力和方向角）。根据式（6）计算松弛变量s0，然后计算每个个体的初始适应度值并进行排序。

2）进行参数初始化，将空间SF和SCR清零。

3）根据式（9）计算精英个体数p。从空间W中随机选取μF和μCR，根据式（12）和式（13）计算Fi，g和CRi，g。

4）根据式（8）的变异规则，选择3 个不同个体，对当前目标向量按式（8）进行变异操作，生成变异个体Vi，g。

5）根据式（10）进行交叉操作，生成试验向量Ui，g，并采用截断的方 式对Ui，g的数值 进行边界处理。

6）根据式（11）进行选择操作，生成下一代的目标向量Xi，g+1，并将失败的目标向量Xi，g存入集合A 中，将生成成功试验向量的Fi和CRi分别存储至空间SF和SCR。

7）判断第g 代种群的个体是否全部完成变异、交叉和选择操作，如果已完成，执行步骤8；否则，返回步骤3。

8）判断是否满足终止条件，即迭代次数为最大迭代次数或s0达到最小要求误差。如果满足，输出实际推力和力矩指令；否则，按式（15）～（21）更新空间W 中的μF，μCR，σ2CR，并返回步骤2。

3 仿真分析

为了验证改进的SHADE 算法解决推力分配问题的可行性，本文将在相同的仿真环境和参数设置下，分别应用改进SHADE 算法和传统SHADE 算法对推进器的推力进行优化分配。本文将采用文献［14］中带有动力定位系统的船舶进行仿真验证，其中整个系统采用Simulink 进行仿真，而数值优化部分则采用Matlab Function 模块对算法进行编程实现。该船长76.2 m，宽18.8 m，质量为4.591×106kg，配置了4 个推进器（2 个全回转推进器和2 个隧道推进器），本文将每个推进器产生的最大推力设为船舶干重的1/60。推进器的具体布置如图4 所示，其相关参数如表1 所示。

图4 推进器的详细布置图Fig.4 Specific layout position of thruster

表1 推进器参数Table 1 Parameters of thruster

设定船舶从初始位置η0=［0，0，0］T航行至期望位置ηd=［60，40，30］T，本文将采用PID 控制器实现动力定位模式仿真，该控制器所产生的期望推力和力矩指令如图5 所示。

在仿真过程中，本文为船舶施加了一个环境扰动b，该环境扰动由白噪声经过二阶近似海浪模型而得，如图6 所示。

图5 动力定位控制器产生的期望推力和力矩指令Fig.5 Expected thrust and torque instruction generated by dynamic positioning controller

图6 环境扰动图Fig.6 Environmental disturbance

由图5 和图6 可见，经仿真验证，改进SHADE算法可以对推力进行合理分配，保证船舶由初始位置航行至期望位置，并在期望位置保持稳定的环境扰动抑制能力。其航行路线和动力定位响应如图7 和图8 所示。由图8 可以看出，当船舶到达期望位置后，在环境扰动下基本可以稳定在原位，并且在3 个自由度上几乎无波动。

图9 所示为改进SHADE 算法和传统SHADE算法的作用效果对比图，在相同的仿真环境下，虽然2 种算法求解推力分配问题时的船舶动力定位响应无明显区别，但传统SHADE 算法所计算的各个推进器推力均较大，这将导致功率消耗和方向角变化也较大。

图7 船舶的航行路线Fig.7 Navigation route of ship

图8 船舶动力定位响应Fig.8 Dynamic positioning response of ship

由图9（a）可以看出，基于改进SHADE 算法的推进器推力的变化曲线较为平缓，没有大幅度振荡的情况，并且未出现2 个相邻推进器所产生的推力相互影响而导致的推力损失问题。由图9（b）可以看出，当船舶到达期望位置后，较小的推力即可使船舶在环境扰动下基本稳定在原位，其中隧道推进器的推力值在-0.1～0.1 kN 范围内，全回转推进器则为-0.5～0.5 kN。由图9（c）中可以看出，基于传统SHADE 算法的推进器推力的变化波动较大，并且在船舶航行过程中出现了相邻推进器推力方向相反所导致的推力损失问题。由图9（d）可以看出，当船舶到达期望位置后，推力变化曲线的振荡较大，且每个推进器均需产生较大的推力才能使船舶在环境扰动下基本稳定在原位。

图9 基于2 种算法的推进器推力变化对比图Fig.9 Thrust change comparison diagram of thruster based on two algorithms

图10 所示为基于2 种算法的全回转推进器方向角的变化对比图。由图10 可以看出，基于改进SHADE 算法的推进器方向角的变化曲线较为平缓，没有出现突变情况，并且当船舶到达期望位置后，全回转推进器的方向角趋于稳定。而基于传统SHADE 算法的推进器方向角在时刻变化，当船舶到达期望位置后也无法抑制环境扰动的影响，故其角度无法稳定。由此可见，改进SHADE 算法可以满足船舶的基本航行要求和推力、方向角变化的约束要求，从而降低全回转推进器的磨损。

图11 所示为基于2 种算法的期望推力与实际推力力矩指令的误差对比图。由图11 可以看出，由动力定位控制器产生的期望推力力矩指令与改进SHADE 算法得到的实际推力力矩指令之间的误差非常小，其中：仿真前期的纵荡和横荡误差在±5×10-5kN 范围内，艏摇力矩误差在±3×10-5kN·m 范围内；仿真后期 的纵荡和横荡误 差 在±1×10-7kN 范 围 内，艏 摇 力 矩 误 差在±1×10-8kN·m 范围内。然而，传统SHADE 算法的误差则较大，且其在船舶到达期望位置后的误差高于航行阶段。因此，利用改进SHADE 算法来解决推力分配问题时，可以满足较高的精度要求，并可将整个过程的误差控制在较小范围内。

图10 基于2 种算法的全回转推进器方向角变化对比图Fig.10 Direction change comparison diagram of full-rotating propulsion based on two algorithms

图11 基于2 种算法的期望推力与实际推力力矩指令的误差对比图Fig.11 Error comparison diagram of torque instruction of expected thrust and actual thrust based on two algorithms

图12 所示为基于2 种算法的推进器总功率变化对比图。由图12 可以看出，基于改进SHADE算法的推进器总功率变化较为平缓，未出现大幅振荡，其在环境扰动下可以使船舶稳定在期望位置，同时保证较小的功率消耗。而基于传统SHADE 算法的推进器方向角和推力变化均存在振荡现象，无法在较短时间内达到最优，故其推进器总功率的振幅和波动较大。因此，基于改进SHADE 算法对船舶进行推力分配，可以较好地抑制环境扰动对船舶航行的影响。

图12 基于2 种算法的推进器总功率变化对比图Fig.2 Total power change comparison diagram of thruster based on two algorithms

4 结 语

本文提出了一种运用改进SHADE 算法解决船舶推力分配问题的方法，仿真验证对比结果表明，在整个航行过程中，推进器的推力和方向角变化较为平缓，可以降低推进器磨损和相互作用影响；推进系统的功率消耗变化几乎无波动，期望推力和实际推力的力矩指令误差可以控制在较小的范围内，其中纵荡和横荡误差在±5×10-5kN 范围内，艏摇力矩误差在±3×10-5kN·m 范围内。

由此可见，该算法可以较好地解决船舶推力分配问题，且在环境扰动下可以保持航向稳定；当船舶到达期望位置后，较小的推力即可使船舶基本稳定在原位。此外，该成果可为如何利用推力分配维持船舶电力系统的稳定性研究奠定基础。