基于改进蚁群算法的海上风电运维船系泊系统逆运动学求解研究*

郭辰

（华能新能源股份有限公司，北京 100036）

基于改进蚁群算法的海上风电运维船系泊系统逆运动学求解研究*

郭辰

（华能新能源股份有限公司，北京 100036）

本文针对一种球形电动机驱动的系泊系统，提出了基于改进蚁群算法的系泊系统逆运动学求解方法。通过对改进蚁群算法应用于系泊系统逆运动学求解的合理性，及参数设置规律等进行仿真研究，对采用此逆运动学求解方法后，系泊系统的控制效果进行了验证。仿真结果表明，本文提出的逆运动学求解方法能够快速准确地找到系泊系统逆运动学解，从而使系泊系统实现精确的轨迹跟踪控制。

系泊系统；逆运动学求解；改进蚁群算法；全局搜索；局部搜索；球面定位；轴向伸缩

0 引言

截至2012年年底，全球海上风电装机容量已达到5410MW。欧洲海上风电发展较快，其中，英国海上风电装机容量超过2900MW，位居世界第一；我国海上风电装机容量也接近400MW[1]。随着我国风电产业技术的发展和相关支持性政策的逐步出台，我国海上风电场将在未来得到迅速发展。

海上风电场运行维护主要依靠海上运维船只。这种船只在进入目标风电机组数米范围内水域时，需依靠其接近系统（根据船只大小及需要，接近系统可能包括搭乘系统和系泊系统）将船身与风电机组连接起来，以起到固定船身、搭载人员或小型部件、维修工具等目的。

1 一种海上运维船舶系泊系统

海上运维船舶的搭乘系统和系泊系统需要建立船身与风电机组特定构件间的稳定联系，故上述系统需要具备空间定位能力。海上运维船及系泊系统示意图如图1所示。

一种系泊系统的执行机构可采用球形电动机。球形电动机结构简单、体积小、重量轻、损耗小、里能指标高、便于控制，能实现自转、俯仰、偏航等三个自由度的运动，可应用于机器人关节等做空间多自由度运动的精密装置中[2-3]。一种球形电动机如图2、图3所示[4-5]。

图1 海上运维船舶及系泊系统示意图

图2 一种球形电动机结构示意图

图3 转子的三自由度运动

配合机械臂等装置，可在一定角度范围内实现空间定位。用球形电动机作为驱动机构实现的系泊系统示意图如图4所示。球形电动机输出轴与机械臂实现刚性连接；末端为机械锁扣装置（机械爪环），可在接触到风电机组挂靠杆时，通过机械机构实现扣锁。

球形电动机实现俯仰、偏航运动时，可带动机械臂进行固定球面内的空间定位；球形电动机自转时，通过输出轴法兰、螺栓杆、传动螺母等机械结构，机械臂能实现伸缩运动。于是，通过对球形电动机的广义欧拉角进行控制，就能实现机械臂末端机械抓环的空间定位，从而实现船身与风电机组结构的可靠连接。

要实现上述目的，关键在于通过对球形电动机欧拉角的控制，实现精确、快速的空间定位；这就需要对球形电动机的逆运动学进行求解，即根据空间定位的要求，迅速求解得到欧拉角的控制要求。本文结合图4所示系泊系统，提出一种基于改进蚁群算法的逆运动学求解方法，从而实现对系泊系统进行精确快速的控制。

球形电动机的运动学分为正运动学和逆运动学。其中，前者是根据各个自由度的欧拉角[4]变化情况来求解电动机转子输出轴的运动情况，是从欧拉角空间到笛卡尔空间的映射的求解问题；后者则是相反的过程，是从笛卡尔空间到欧拉角空间的映射的求解问题。球形电动机的逆运动学求解问题是其进行运动控制，运动分析，离线编程和轨迹规划等的基础。目前，关于球形电动机逆运动学分析求解的文献还很少。

意大利学者M. Dorigo在1991年首次提出了蚁群算法[6]。蚁群算法是一种全局优化的搜索算法,具有较强的鲁棒性，且易于与其他算法结合，已经成功地应用于旅行商问题(TSP)、资源二次分配等经典优化问题，取得了良好的效果[7-9]。

图4 系泊系统结构示意图

2 系泊系统运动学逆问题

由前所述，系泊系统运动学逆问题，实际上是球形电动机运动学逆问题与轴向伸缩问题的组合，因此，球形电动机的逆运动学问题是难点所在。以下对其进行研究。

球形电动机的定子位置和转子位置分别用静坐标系xyz和动坐标系dqp来定义，转子输出轴与dqp系中的p坐标轴重合。若x轴从xyz系旋转α角到x1y1z1系，y轴从x1y1z1旋转β角到x2y2z2系，z轴从x2y2z2系旋转γ角到dqp系，则所产生的角度θ＝(α,β,r)T称为广义欧拉角[5-10]。角度余弦cos简记为c，角度正弦sin简记为s，则所产生的旋转矩阵A可表示如下：

此旋转矩阵满足如下关系式：

由式（2）可知，永磁球形电动机输出轴在某一时刻t的位置向量和欧拉角向量之间的关系可以表示如下：

其中，F为一个(3×s)的矩阵，其中，s的具体取值取决于初始坐标值式（3）即为球形电动机的正向运动学方程。

一种逆运动学求解策略是建立球形电动机的微分运动关系，根据转子输出轴在xyz坐标系中的速度向量来求解其对应的欧拉角速度向量。对式（3）两边同时求导数，可得：

于是可以求出欧拉角速度向量：

式中，J-1(θ(t))为雅可比矩阵的逆矩阵。

然后根据初始条件对式（5）两边同时求积分，即可得到永磁球形电动机的逆运动学方程。

用上述方法求球形电动机的雅可比矩阵及其逆矩阵，计算比较复杂；且对于不同的初始位置坐标，雅可比矩阵的形式不同。因此采用这种方法进行逆运动学求解比较困难。

3 一种改进蚁群算法及其实现

蚁群算法利用一群人工蚂蚁来模拟真实蚂蚁的行为，通过人工蚂蚁之间的协作来寻找较为优化的解。每只人工蚂蚁代表一个计算单元，在每次迭代中，负责构建解问题的一条路径，计算并存储该路径的目标函数值，同时在所经过的路径上释放一定数量的信息素。在后一步的迭代中，蚂蚁能够检测到前一步迭代中积累的信息素的浓度，并据此选择自己的前进方向；同时，信息素会随着时间的推移逐渐挥发掉。于是，路径的长短及该路径上通过蚂蚁的多少就对残余信息素的浓度产生了影响；同样，残余信息素浓度的大小又指导着后来蚂蚁的行动方向。因此，某条路径上走过的蚂蚁越多，则后来的蚂蚁选择该路径的概率就越大。

由式（3）可知，本文中的逆运动学问题是一个三维函数F(θ)的求解问题，即通过优化算法求得3个欧拉角α，β和γ的数值。由转子输出轴的初始位置坐标和所求得的欧拉角，可以确定其旋转后的位置坐标令目标函数为：

式中(xd, yd, zd)为给定的转子输出轴的位置坐标值。目标函数值越小，代表根据逆运动学解法得到的(xe, ye, ze)与给定的(xd, yd, zd)越接近，即该解法的求解精度越高。

在永磁球形电动机逆运动学求解中，待求解的参数共有3个，分别为α，β和γ。普通蚁群算法在每个参数的取值范围内随机产生N个值作为备选数值点，这样产生的备选数据点就有可能没有覆盖最优解附近的范围，从而导致寻优结果不理想。本文中，这3个参数的取值范围分别设置为[-1, 1]，则逆运动学的解空间为[-1, 1]×[-1, 1]×[-1, 1]。将这3个参数的取值范围分别均匀地离散化为N个值，在本文中，称N为单维离散化率。

3.1 全局搜索

由前所述，逆运动学的解空间中备选数值点的数目为N×N×N。对于其中任意一个参数，将其取值范围内的每个值看作一个元素，则N个元素形成一个集合，设为Ii(i=1,2,3)。定义蚂蚁的数目为m，全部蚂蚁从蚁巢出发寻找食物。每只蚂蚁从集合I1出发，根据集合中每个元素的信息素状态和式（7），独立随机地从每个集合Ii中唯一地选择一个元素；蚂蚁在所有集合中完成元素的选择（即完成全局搜索），然后在所选择元素周围完成局部搜索后，它就到达了食物源，之后调节集合中各个元素的信息素。这一过程反复进行，直至找到最优解。全局搜索的步骤如下：

（1）初始条件：令集合Ii(i=1,2,3)中的元素j(j=1,…,N)的信息素初始值phj(Ii)(0)=C，迭代次数初值Nc=1，设置最大迭代次数Nc_max。

（2）启动所有蚂蚁，每只蚂蚁从集合Ii(i=1,2,3)开始，按照下述规则依次在每个集合中选择一个元素，直到蚂蚁全部选择完毕。

路径选择规则：对于集合Ii(i=1,2,3)，任意一只蚂蚁k(k=1,…,m)，根据下式计算的概率随机地选择它的第j个元素。

（3）当每只蚂蚁在每个集合中都选择一个元素后，计算由各个蚂蚁所选数值作为欧拉角参数时的目标函数值，并记录其中的最小值及其对应的欧拉角参数。设上述蚂蚁觅食过程经历了m个时间单位，对所有集合Ii(i=1,2,3)中的各个元素的信息素按照下式进行调整：

其中，参数ρ(0≤ρ＜1)表示信息素的持久性，则1-ρ表示信息素的消逝程度。表示在本次循环中第k只蚂蚁在集合Ii的第j个元素上留下的信息素，可以用下式来表示：

其中，Q为常数，用来调节信息素的调整速度；Fk是以蚂蚁k选择的三个元素分别作为欧拉角数值时的目标函数值；η为Fk的指数，不同的η值可以得到不同的寻优效果。由式（10）可以看出，目标函数越小，信息素的增量就越大。普通蚁群算法中，η=1；在本文所提出的改进蚁群算法中，η=2。由后面的仿真分析可以看出，η=2时蚁群对欧拉角参数的寻优速度会明显提高。当然，η取值过大可能会造成寻优过程不稳定，并且容易陷入局部极小。

3.2 局部搜索

由上所述，全局搜索是根据解空间中的N×N×N个备选数值对应的目标函数值来确定每次最优的目标函数值对应的数值点。全局搜索的对象是各个离散点，点与点之间的值就会被忽略掉，由此影响解的质量。为了得到更加优化的解，本改进蚁群算法相比普通蚁群算法增加了局部搜索，即让蚂蚁在其所选择的数值点周围的一个小的邻域内进行搜索，并对原数值点进行相应的微量移动，以使得数值点对应的目标函数值更小，由此提高解的质量。

在局部搜索过程中，需要判断蚂蚁进行微量移动的方向和移动的步长。为了确定移动的方向，先根据当前集合Ii(i=1,2,3)中第j(j=1,…,N)个元素对应的信息素来确定该元素上应该有的蚂蚁数目。具体计算公式如下：

要判断蚂蚁的移动方向，还要借助于蚂蚁所在元素左边和右边的所有元素上的应有蚂蚁数目和实际蚂蚁数目。集合Ii(i=1,2,3)中第j(j=1,…,N)个元素左边所有元素上的应有蚂蚁数目如下：

集合Ii(i=1,2,3)中第j(j=1,…,N)个元素右边所有元素上的应有蚂蚁数目如下：

集合Ii(i=1,2,3)中第j(j=1,…,N)个元素左边、右边以及第j个元素上的实际蚂蚁数目Nsl-j(Ii)、 Nsr-j(Ii)和Nsj(Ii)可以通过蚂蚁当前的分布直接得出。蚂蚁的移动方向如表1所示。

表1中所列的7种情况以外的其他情况下，蚂蚁将不移动。由表1可以看出，蚂蚁的移动方向总是向着实际蚂蚁数目少于应有蚂蚁数目的方向。离散空间中的离散点越多，即元素数目越大，各个元素上应有蚂蚁数跟实际蚂蚁数的差别就越小，此时的离散空间也就越接近于连续空间。本文的局部搜索就是要在不改变解空间中元素数目的情况下，让离散解空间最大限度地接近连续空间，从而提高解的质量。定义差别度变量D，对于集合Ii(i=1,2,3)中第j(j=1,…,N)个元素来说，差别度D如下所示：

其中，b由如下规则赋值：

于是，每次微量移动的步长定义如下：

其中，ε为局部搜索步长因子。在本文中，取ε=0.02/m，其中m为蚂蚁数目。

在蚂蚁完成局部搜索后，将新的元素值代替局部搜索前的元素值，重新进行全局搜索。全局搜索和局部搜索交替进行，直到得到最优解，或者达到最大迭代次数为止。事实上，局部搜索的运算量较小，蚁群算法的运算量和运算耗时主要体现在全局搜索阶段。蚂蚁数目m和单维离散化率N较小时，蚁群算法的运算量和运算耗时也较小。

本文提出的改进蚁群算法运行流程图如图2所示。图中，Nc表示迭代次数，Nc_max表示最大迭代次数。

表1 蚂蚁移动方向判断规则

改进蚁群算法实现了在任意初始位置下，从笛卡尔空间到广义欧拉角空间的转换。然而，对于球面求解问题而言，笛卡尔空间三个变量x, y, z之间并非完全独立，而满足球形约束关系，只有其中两个向量是完全独立的；这就是说，（3）式中F在固定球面内的秩为2；广义欧拉角向量也只需两个欧拉角变量即可在特定球面上确定空间位置。结合欧拉角的物理意义可知，自转欧拉角γ在固定球面空间的定位问题中为冗余解；也就是说，在球面定位过程中，自转欧拉角γ的控制指令为零。

由以上的分析可知，系泊系统要想从图1(a)的状态到达图1(b)的状态，需要首先实现球面内的定位，然后实现伸缩运动，即通过自转运动推进此机械臂，使系泊系统末端的机械抓环到达挂靠杆。

以下建立自转运动对应的广义欧拉角与伸缩量之间的关系。设螺栓杆的螺距为，则两者关系如下：

4 仿真研究

图5 改进蚁群算法运行流程图

以下对本文提出的改进蚁群算法进行仿真研究，取Q=200，ρ=0.7。设转子球体的半径为R，仿真中，电动机输出轴的转子球面位置点的初始坐标为(xi,yi,zi)T=(0, 0, R)T，转子旋转后的坐标为(xe, ye, ze)T=(0.433R, 0.500R, 0.750R)T。首先比较在采用普通蚁群算法和改进蚁群算法时目标函数最小值随迭代次数的变化情况，如图6所示。图6中，蚂蚁数目m=40，解空间内的单维离散化率N=20，最大循环次数设置为Nc_max=30。

由上图可以看出，在具有相同起点的情况下，采用本文提出的改进蚁群算法时目标函数最小值的收敛要明显快于普通蚁群算法。事实上，当m和N取其他数值时，也有类似的规律。这是由于改进蚁群算法中η=2，信息素的增量对不同大小的目标函数值更加敏感，每次迭代中使得目标函数值最小的欧拉角数值就会被加速强化，提高其下次被选择的概率，从而提高整个算法的收敛速度。由于局部搜索的作用，改进蚁群算法更容易找到最优解，因此图6中两种算法所得到的目标函数最小值不同。采用改进蚁群算法，到第5步时目标函数最小值为0.0370，到第20步时目标函数最小值为0.0351。这是由于在第3步以后，局部搜索起主要作用，找到的数值点仍然向着最优的方向微量移动。

为了便于理解蚁群算法，以下对图6中采用改进蚁群算法时蚂蚁分布随迭代次数的变化情况进行仿真，如图7所示。其中，图7(a)表示Nc=1，即迭代开始时的蚁群分布情况，图7(b)表示Nc=5时的蚁群分布情况，图7(c)表示Nc=10时的蚁群分布情况，图7(d)表示Nc=15时的蚁群分布情况。

从上图可以看出，在迭代开始时，蚂蚁分散地分布在解空间的元素上。根据前述信息素的初值及各个元素被选择概率的关系，蚁群此时的分布是随机的，选择解空间内任意元素的概率是相等的。随着迭代的进行，蚁群的分布越来越集中，并逐渐汇集于一点，如图7(d)所示。这个点正是蚁群寻优的收敛点，是改进蚁群算法在迭代中找到的最优解。

由图6可以看出，改进蚁群算法在运行至第5步迭代时，目标函数最小值已稳定在一个基值上（以后各步中局部搜索起主要作用），说明目标函数最小值基本已经找到，但是并不说明此时所有蚂蚁都聚集在一点。这从图7(b)中也可以看出。

改进蚁群算法中的各个参数设置会改变仿真的结果。为了进行准确的逆运动学求解，应当选取合适的参数。与普通蚁群算法一样，本文提出的算法中多数参数根据经验设置。在仿真实验中，Q和ρ两个参数的取值对仿真结果的影响不大。对仿真结果影响最大的是解空间的单维离散化率N和人工蚁群中蚂蚁的数目m。当蚁群中蚂蚁数目m=15，最大迭代次数为Nc_max=20时，目标函数最小值和逆运动学求解结果随单维离散化率N的变化情况如图8所示。

由上图可以看出，随着N的增大，目标函数最小值逐渐减小，逆运动学求解结果也趋于稳定。事实上，N越大，解空间的划分就越细密，越有利于蚂蚁在较好的解空间区域中进行更加细密的局部搜索。因此，N越大就越容易搜索到较好的目标函数最小值和全局优化的解。

当单维离散化率N=20，最大迭代次数为Nc_max=20时，目标函数最小值和逆运动学求解结果随蚁群蚂蚁数目m的变化情况如图9所示。

由上图可以看出，蚂蚁数目较小时，目标函数的最小值出现波动，且逆运动学求解的结果也不稳定，这是由于，此时蚁群的正反馈作用较强，全局随机搜索能力较弱，容易陷入局部极小。随着蚂蚁数目m的增大，目标函数最小值开始稳定减小，表明逆运动学求解结果在逐渐变得精确。这是因为，蚂蚁数目较大时，蚁群的全局随机搜索能力就会加强，就越容易搜索到较好的解。

图10表示系泊系统将采用本文提出逆运动学求解方法得到的欧拉角作为控制系统的输入量，并采用PD控制策略时，系泊系统末端机械抓环在笛卡尔空间的运动轨迹。设球形电动机转子中心为圆心，转子中心到机械抓环之间的距离为半径，此半径设为单位1。图中，球面上的红色圆圈为球面定位运动的指令位置，球面以外的红色菱形表示伸缩运动的指令位置。

从图10(a)可以看出，系泊系统完成球面定位后，机械抓环的位置与红色圆圈表示的指令位置重合度较高。从图10(b)可以看出，完成伸缩运动后，机械抓环末端位置与红色菱形表示的指令位置几乎完全重合。图10说明，采用本文提出的逆运动学求解方法后，系泊系统进行系泊操作的定位精度较高。

5 结论

图6 普通蚁群算法和改进蚁群算法中目标函数最小值的变化情况

图7 蚂蚁分布随迭代次数的变化情况

图8 单维离散化率N对目标函数最小值和逆运动学求解结果的影响

本文针对一种球形电动机驱动的系泊系统，提出基于改进蚁群算法的系泊系统逆运动学求解方法。改进蚁群算法能够在解空间中对逆运动学解进行全局与局部搜索寻优，得到的解具有较高的精度，且比普通蚁群算法的收敛速度快。基于改进蚁群算法的逆运动学求解方法，能够实现准确的球面定位和伸缩定位，从而实现系泊系统的轨迹跟踪控制。仿真结果验证了该算法的合理性，对该算法中蚂蚁数目和单维离散化率参数与求解结果之间的影响机制进行了分析研究，为算法的参数设置提供了依据；仿真结果显示，以逆运动学求解得到的欧拉角作为控制系统输入信号，在采用PD控制算法时，系泊系统可以得到良好的轨迹跟踪效果。

图9 蚁群蚂蚁数目m对目标函数最小值和逆运动学求解结果的影响

图10 系泊系统采用本文提出逆运动学求解方法及PD控制策略时的轨迹跟踪效果

[1] 2012年中国风电装机容量统计. 中国可再生能源协会风能专业委员会.

[2] Liang Yan, I-Ming Chen, Chee Kian Lim, Guilin Yang, Wei Lin, Kok-Meng Lee. Design and Analysis of a Permanent Magnet Spherical Actuator[C]. IEEE/ASME Transactions on Mechatronics. 2008, 13(2): 239-248.

[3] Klemens Kahlen, Ingo Voss, Christian Priebe, Rik W. De Doncker. Torque Control of a Spherical Machine with Variable Pole Pitch[J]. IEEE Transactions on Power Electronics, 2004, 19(6): 1628-1634.

[4] 夏长亮, 李洪凤, 宋鹏. 基于Halbach阵列的永磁球形电动机磁场研究[J]. 电工技术学报, 2007, 22(7): 126-130.

Xia Changliang, Li Hongfeng, Song Peng. Magnetic field model of a PM spherical motor[J]. Transactions of China Electrotechnical Society, 2007, 22(7): 126-130 (in Chinese).

[5] Changliang Xia, Chen Guo and Tingna Shi. A Neural Network Identifier and Fuzzy Controller Based Algorithm for Dynamic Dycoupling Control of Permanent Magnet Sperical Motor [J], IEEE Transactions on Industrial Electronics, pp. 361-372, vol. 59(1), 2012.

[6] Marco Dorigo, Luca Maria Gambardella. Ant Colony System: a Cooperative Learning Approach to the Traveling Salesman Problem[J]. IEEE Transactions on Evolutionary Computation, 1997, 1(1): 53-66.

[7] 翟海保, 程浩忠, 吕干云, 陈小良, 马则良.基于模式记忆并行蚁群算法的输电网规划[J]. 中国电机工程学报, 2005, 5, 25(9): 17-22.

[8] 赵庆杞, 黎明, 张化光. 基于蚁群算法的灵敏负荷调度[J]. 中国电机工程学报, 2006, 12, 26(25): 15-21.

[9] 余玲, 刘康, 李开世. 蚁群算法的连续空间算法研究[J]. 机械设计与研究, 2006, 4, 22(2): 6-9.

[10] 黄声华. 三维电动机及其控制系统[M]. 武汉:华中科技大学出版社, 1998.

Study on Problem of Inverse Kinematics of Mooring System of Offshore O&M Vassel Based on Advanced Ant Colony Algorithm

Guo Chen
（Huaneng Renewables Co., Ltd., Beijing 100036, China）

For a Spherical Motor drived MS of ofshore O&M vessel, one Advanced Ant Colony Algorithm (AACA) based IKS strategy is proposed in this paper. Through the axial moment of Spherical Motor, Axial Telescopic (AT) control of MS can be completed and the Trajectory Tracking of MS is realized. The proposed method for problem of Inverse Kinematics is validated, and the parameter seThing rule of AACA is researched. Simulation result indicated that the proposed IKS strategy for MS can solve the problem of IKS accurately and actual trajectory tracking.

Mooring System; Inverse Kinematics Solution; Advanced Ant Colony Algorithm; global search; local search; Spherical Positioning; telescopic positioning

TM614

A

1674-9219（2013）12-0064-08

国家863计划课题(2012AA051706)和(2012AA051703)资助项目。

2013-10-02。

郭辰（1982-），男，工学博士，主要从事风能资源评估及微观选址技术、海上风电场接近技术及系统、非线性系统动力学建模及仿真、智能控制策略及其应用、特种电机及其智能控制等方面的研究。