基于改进粒子群优化算法的船舶纵向运动参数辨识

戴运桃 ， 赵希人 ， 刘利强

（哈尔滨工程大学 a.理学院；b.自动化学院，哈尔滨 150001）

1 引 言

水动力参数是船舶操纵设计的重要参数，但是有些参数难以用实验方法获得，而理论计算值和实验值都与实航实验有较大区别。因此，水动力参数的获得一般都是应用系统辨识技术，从船舶的运动状态观测数据中辨识出来，从而直接建立船舶的水动力参数和运动状态之间的数学模型。水动力参数辨识方法一般采用经典的辨识方法如极大似然辨识方法[1]和预报误差法[2]等。但由于纵向运动状态对阻尼力和力矩等参数不敏感，灵敏度非常低，而常规辨识算法对灵敏度较低的参数是非常难以求解到真值的。

粒子群优化算法（Particle Swarm Optimization，PSO）[4]是一种新的基于群体智能的优化算法，它仿效生物学中进化和遗传的过程，同时考察多个候选解，淘汰劣质解，鼓励发展优质解，逐步提高解群体的质量，从而逼近所研究问题的最优解，为解决优化问题提供了新思路[10]。本文针对常规辨识方法在水动力参数辨识中的局限性，探索了模拟生物智能的粒子群优化算法在船舶水动力参数辨识中的应用。

2 问题描述及数学模型

2.1 状态方程

根据船舶水动力理论[8]，船舶纵向运动的简化方程如公式（1）所示。

式中：z表示船体的升沉，θ表示船体的纵摇，I5表示纵摇惯性矩，Δ表示排水量，FR表示水平舵升力，XR为舵升力中心至船体重心的纵向距离，Fw3表示海浪升沉干扰力，Mw5表示海浪纵摇干扰力矩；a33表示附加质量；a35，a53表示质量矩（kg·m）；a55表示转动惯量（kg·m2）；b33，b35表示阻尼力系数；c33，c35表示复原力系数；b53，b55表示阻尼力矩系数；c53，c55表示复原力矩系数。 a33，a35，a53，a55，b33，b35，b53，b55，c33，c35，c53，c55统称为水动力参数。

令 x1=z，x3=θ，设状态变量为

则得到系统状态方程如下：

其中，A=E-1A*，B=E-1B*，C=E-1C*，

para=[a33，a35，a53，a55，b33，b35，b53，b55，c33，c35，c53，c55]为待辨识参数，W表示海浪扰动力和力矩。

2.2 观测方程

选择状态x1、x3为测量状态，则

式中：

Y为观测向量，V为二维测量噪声，通常情况下可以认为是白噪声，它的方差阵按一级精度的传感器可以取为Qvv=diag[ 20.3e-4 2.25e-6]。

2.3 离散方程

在进行计算机仿真前，首先需要对状态方程和观测方程进行离散化。设Ts为采样周期，在连续的控制对象前有一个零阶保持器，即u（t）=u（k ），kTs＜t＜（k+ 1） Ts。 忽略量化效应，并将连续的控制对象和它前面的零阶保持器一起离散化，从而使整个系统变为纯粹的离散系统。则离散化后状态方程为：

其中，Φ=eATs，

3 粒子群优化算法

粒子群优化（Particle Swarm Optimization，PSO）算法是由Kennedy和Eberhart借鉴鸟类寻找食物的自然现象提出的一类基于种群的随机全局优化技术[4-5]。在鸟群的飞行中，每只鸟在初始状态下处于随机位置，且朝各个方向随机飞行，但随着时间的推移，这些初始处于随机状态的鸟通过相互学习、相互跟踪，自组织地聚集成一个个小的群落，并以相同的速度朝着相同的方向飞行，最终整个群体聚集到同一位置—食物源。

PSO算法在可行解空间中初始化产生一组粒子，每个粒子在搜索空间中以一定的速度飞行，并根据对个体和集体的飞行经验的综合分析来动态调整个体自身的飞行速度，同时向着自己以前经历过的最好位置和其它微粒曾经经历过的最好位置飞行，依靠这种个体间的信息交换来达到整个群体的共同演化。实现复杂空间全局最优解的搜索。与其他进化算法相比，PSO算法的主要特点是：概念简单，实现容易，需要调整的参数相对较少，可以用于解决复杂非线性优化问题。下面以一个最小化问题的求解为例，给出PSO算法的流程。

步骤1：在整个搜索空间中随机地初始化粒子群中各粒子的位置xi和速度vi，并计算每个粒子在当前位置处的适应度值fitnessi=fitness（ xi）。相应的初始化粒子自身的最优解（个体极值pbesti）和整个粒子群的当前最优解（全局极值 gbest），即 pbesti=fitnessi，gbest=min（fitness1,…， fitnessn），gbest越小越好。

步骤2：在每次迭代过程中，粒子通过跟踪其自身的个体极值和全局极值，按公式（6）来进行自身速度和位置的更新。

式中：vi表示粒子的速度向量；xi代表粒子的位置向量；c1、c2为常数，称为学习因子；r1、r2是[0，]1上均匀分布的随机数；w为惯性权重。

步骤3：更新个体极值pbesti和全局极值gbest。

步骤4：如果已经满足停止准则，如gbest达到某个阈值或者已经达到最大迭代次数，算法中止计算；否则跳转到“步骤2”。

4 改进粒子群优化算法

4.1 算法基本思想

许多优化算法在寻优的过程中都会出现早熟问题，粒子群算法也不例外。在粒子群算法中，各粒子是根据自身的最优位置与全体粒子的最优位置来调整搜索方向的，在算法运行的初始阶段，收敛速度比较快，但运行一段时间后，速度开始减慢甚至停滞，当所有粒子的速度几乎为0时，粒子群丧失了进一步进化的能力，可以认为算法执行已经收敛。但在许多情况下，算法并没有收敛到全局极值，甚至连局部极值也未达到。这是因为发生该现象时粒子群高度聚集，严重缺乏多样性，粒子群会长时间或永远跳不出聚集点[9]。

针对这一问题，本文设计了一种改进粒子群优化算法，在算法运行的每一迭代步中计算全局最优粒子的适应值globalBestValue（即全局最小适应值），若在设定的迭代步changeNum内，全局最小适应值globalBestValue一直都不变化，则代表各粒子已经开始高度聚集，算法在当前的位置和速度下已无法找到更优解。此时，利用高斯变异算子变异各粒子的位置，使粒子从新的位置开始继续搜索。随着粒子变异次数的增加，算法不断缩小高斯变异算子的变异范围，从而使粒子在逐渐缩小的搜索空间进行精细寻优，直到搜索到全局最优解。

4.2 高斯变异算子设计

选择高斯密度函数作为各粒子的位置密度分布模型。如（7）式所示。

其中，xmini，d为适应值最小的粒子的位置，σd表示本次变异中的第d个参数的方差。

算法对粒子的位置密度分布函数进行采样，得到各粒子的位置。采样过程中，使用函数p（ xi）所对应的随机数生成器进行，该生成器如公式（8）所示。

其中，randl为均匀分布在[0，1]之间的一个随机变量。

随着算法变异次数的增加，利用公式（9）不断缩小粒子位置的变异范围。

其中，ρ（ 0＜ρ＜ 1）为收缩因子，k为变异次数。 在算法开始运行阶段，方差 σd较大，这样可以使粒子在较大的范围内变异，从而保持粒子的多样性；随着粒子位置不断变异，方差σd越来越小，算法开始在更小范围内进行精细搜索。

5 船舶纵向运动参数辨识算法设计

本文设计的参数辨识算法的思想是将粒子的位置作为待辨识的水动力参数。首先，用给定的水动力参数求得一组理论观测值；然后，在算法运行过程中，利用各粒子的位置（辨识得到的水动力参数）求得实际观测值，并使用适应值评价函数对本组水动力参数作出评价。这一评价即为粒子当前位置的适应度值。各粒子使用个体最优适应度值和全局最优适应度值进行速度更新和位置更新，算法迭代运行，从而引导各粒子向着适应度值高的方向前进。

适应值评价函数是评估算法辨识得到的水动力参数好坏的重要依据。如公式（10）所示，在算法中我们使用理论观测值和实际观测值之间对应数据的误差平方和的均值作为适应值评价函数。函数值越小，代表该组辨识得到的水动力参数越好。

其中，Yreali是根据实船试验得到的水动力参数计算所得的理论观测值，Yi是根据辨识得到的水动力参数计算得到的实际观测值，dataNum是观测次数。算法的流程图及运算步骤如下：

步骤1：初始化，设定学习因子c1、c2和惯性权重w，在允许范围内，对各粒子的初始位置和速度进行随机选取，利用给定的水动力参数和船舶纵向运动方程解算理论观测值；

步骤2：使用每个粒子对应的位置信息（水动力参数）和船舶纵向运动方程解算实际观测值，并使用适应值评价函数得到各粒子当前位置的适应度值；

步骤3：对各粒子比较其当前位置适应度值和个体极值pbest，如果更好，则用当前位置适应度值来更新pbest；

步骤4：用每个粒子的个体极值pbest与全局极值gbest比较，若更好，则更新gbest；

步骤5：判断是否满足终止条件，若满足，则输出全局极值gbest和其所对应的位置；否则，转到“步骤 6”。

步骤6：判断全局最优适应值是否有连续changeNum次未改变，如果是，则按照公式（8）重新生成粒子位置，否则，转入“步骤2”。

6 仿真分析

图1 参数辨识算法流程图Fig.1 Flowchart of parameter identification algorithm

针对上述设计的参数辨识算法，本文使用Visual C++和MATLAB作为工具进行了计算机仿真。由文献[3]可知，c33,c35,c53,c55为常数，因此，在辨识过程中，将c33,c35,c53,c55当作已知参数，然后对a33,a35,a53,a55，b33,b35,b53,b55进行辨识。在仿真实验中，二维测量噪声V取均值为零的高斯白噪声，输入舵角是幅值±10 为伪随机数，理论观测值使用文献[8]中海情3 级航速 18 节，航向分别为 0°，45°，90°，135°，180°的五组水动力参数计算得到。

本次实验采用基本的PSO算法与本文改进的PSO算法进行比较，两算法基本参数设置相同，各参数如下：粒子数60，c1=c2=1.496 2，惯性权重随迭代次数减少，即从初始的0.9线性递减到0.4，各参数的搜索范围为±200%，算法迭代次数为1 000次。在改进算法特有的参数中，收缩因子ρ=0.9，changeNum=20。 针对海情 3 级航速 18 节，航向分别为 0°，45°，90°，135°，180°的五组水动力参数的辨识问题，算法收敛的标准是所有的水动力参数的辨识误差均收敛到±10%以内。

两算法对于每个航向各运行100次，算法收敛速度和计算成功率如下表1所示。表1中各参数解释如下：meanN表示所有达到收敛要求的迭代次数的均值；minN表示达到收敛要求的最小迭代次数；而maxN则是达到收敛要求的最大迭代次数；succeed是100次搜索中达到收敛要求的次数的比例。

表1 PSO与改进PSO算法在水动力参数辨识问题的对比Tab.1 Comparison between PSO and improved PSO in hydrodynamic parameter identification

从表1中我们可以看到，改进的PSO算法在求解船舶纵向运动水动力参数辨识问题时稳定性远高于基本的PSO算法，而且平均迭代次数也小于PSO算法的迭代次数。

表2给出了辨识得到的船舶纵向运动水动力参数仿真结果。其中，水动力参数的理论值是基于切片理论通过船模在水池做大量实验并经测试计算得到，而辨识值是通过改进的PSO算法辨识得到，N表示收敛到次数是±10%以内所需的迭代次数，相对误差的计算是：相对误差=（理论值-辨识值）/理论值×100%。

表2 船舶纵向运动参数辨识结果Tab.2 Parameter identification result of ship vertical motion

图2给出了3级海情下，航速18Kn，航向45°时算法运行过程中各水动力参数的变化曲线。

图2中横坐标表示迭代次数，纵坐标为各水动力参数在不同迭代次数时的辨识值。直线表示各水动力参数的理论值，曲线为各参数在不同迭代次数时的值。从图中可以看到，粒子出现震荡收敛现象，在最初始迭代时震荡较大，而后缓慢收敛。这是因为在初始阶段，变异的粒子的位置是在较大方差下给出的，保证了粒子的多样性，增加了跳出局部最优的概率。在随后的过程中方差逐渐减小，粒子位置缓慢趋于理论值。

图3和图4分别为升沉曲线图和纵摇曲线图。横坐标为观测的次数，纵坐标分别为升沉值（单位为米）和纵摇值（单位为度）。“*”为理论上的观测值，曲线表示通过辨识的参数得到的模型输出值。从图上可以看出，理论观测值与模型输出值基本吻合，说明辨识结果能较准确地描述船舶的动力学特性。

图2 海情3级，航速18Kn，航向45°的船舶纵向运动参数收敛曲线图Fig.2 Parameter identification of the ship vertical motion for 3 sea condition,18Kn of ship speed,45°of course

图3 升沉的辨识曲线Fig.3 Identification of heave

图4 纵摇的辨识曲线Fig.4 Identification of pitch

7 结 论

船舶水动力参数辨识是获取船舶水动力参数的一种重要手段。本文使用改进的粒子群优化算法设计了一种船舶纵向运动水动力参数辨识算法，介绍了参数辨识算法的思想，给出了算法流程和适应值函数，并进行了计算机仿真。仿真结果表明，改进的PSO算法在求解船舶纵向运动水动力参数辨识问题时稳定性好，而且辨识精度高，能够满足实际需要。

[1]丁文镜,罗仁凡等.受控航行体水动参数的极大似然辨识[J].清华大学学报（自然科学版）,1992,2(32):89-95.

[2]丁文镜,罗仁凡等.辨识航行器水动力参数的预报误差法[J].舰船科学技术,1994(5):22-27.

[3]Dong R R,Joseph Katz.On the structure of bow waves on a ship model[J].Journal of Fluid Mechanics,1997,346:77-115.

[4]Kennedy J,Eberhart R.Particle Swarm Optimization[C]//IEEE Int’l Confon Neural Networks.Perth,Australia,1995:1942-1948.

[5]Eberhart R,Kennedy J.A new optimizer using particle swarm theory[C]//In:Proc of the Sixth International Symposium on Micro Machine and Human Science.Nagoya,Japan,1995:39-43.

[6]Shi Y,Erberhart R C.A modified particle swarm optimizer[J].In:Proc of the IEEE CEC,1998:69-73.

[7]Liu Bo,Wang Ling,et al.Improved particle swarm optimization combined with chaos[J].In:Chaos,Solitons and Fractals,2005,25:1261-1271.

[8]陈虹丽等.船舶纵向运动水动力参数模型的建立及分析[J].仪器仪表学报,2006,27(6):1120-1121.

[9]李 莉,李洪奇.基于混合粒子群算法的高维复杂函数求解[J].计算机应用,2007,27(7):1754-1756.

[10]陈纬琪,颜 开等.水下航行体水动力参数智能辨识方法研究[J].船舶力学,2007,11(1):40-46.