移动装弹机械臂的逆运动学多种群灰狼算法求解方法

Solving method on inverse kinematics of mobile loading missile manipulator by multi-population grey wolf optimization algorithm

摘要

针对移动装弹机械臂逆运动学求解性能需要提升的问题，提出了一种基于多种群灰狼算法的逆运动学求解方法。首先，将机械臂的逆运动学问题转化成了一个等效最优问题，并依据优化目标建立了适应度函数。其次，基于灰狼算法拓展了灰狼种群数量，引入了粒子群算法及最优个体反向指导的位置更新方式，并且设置了随机重组阈值淘汰机制。然后，应用多种群灰狼算法循环迭代求逆，当适应度函数值趋近于0即可得到逆解。最后，通过与其他算法求解仿真对比表明：本文所提出的逆运动学求解方法具有更好的收敛性、求解精度及重复精度。

Abstract

Aiming at the problem that the inverse kinematics performance of the mobile loading manipulator needs to be improved， an inverse kinematics solution method based on the multi-population gray wolf algorithm is proposed. Firstly， the inverse kinematics problem of the mobile loading manipulator is transformed into an equivalent optimal problem， and the fitness function is established according to the optimization objective. Secondly， based on the gray wolf algorithm， the grey wolf population is expanded， the particle swarm algorithm and the position updating method of the optimal individual inverse guidance are introduced， and the random reorganization threshold elimination mechanism is set. Then， the algorithm is applied to iteratively invert so that the fitness function approaches 0 to obtain the inverse solution. Finally， the simulation comparison with other solving methods shows that the proposed method has better convergence， solution accuracy and repeatability.

Graphical abstract

关键词

控制理论与控制工程 / 移动装弹机械臂 / 逆运动学求解 / 多种群灰狼算法 / 适应度函数

Key words

control science and control engineering / mobile loading missile manipulator / kinematic inversion / multi-population grey wolf optimization algorithm / fitness function

0 引 言

移动装弹机械臂自动化程度高，转载炮弹时平稳性和快速性好，受到越来越多的学者和军事专家关注。在影响其打击性能的诸多要素中，逆解精度决定了路径规划和位置控制研究的可靠性，更影响着军事打击的精度与重复精度。

构造简单的机械臂可以采用解析法或数值法进行逆运动学的求解^［1］。相比之下，移动装弹机械臂的冗余和多关节类型特性，使它的逆运动学问题变成了一个冗余、多变量、强耦合的非线性问题^［2］。此时若仍采用传统求解方式，整个过程会烦琐且不易操作。而智能优化算法的高速发展为该问题的公共难点——正运动建模及逆运动求解^［3］的解决提供了新思路。文献［4］采用粒子群算法（Particle swarm optimization algorithm， PSO）求解了6自由度库卡KR 6 R900机械臂逆运动学问题；文献［5，6］在此基础上，分别引入边界限制和自适应权重调整策略改进了PSO算法；文献［7］通过增添激励函数改进算法，解决了传统BP神经网络的收敛问题，并将其用于求解6自由度机器人的逆运动学；文献［8］为解决典型冗余模型（人体上肢）的逆解问题，提出了搜索空间改进的差分进化自适应控制参数算法，在上肢康复机器人的具体应用中有重大意义；文献［9］改进了人工蜂群算法（Artificial bee colony algorithm， ABC），提高了算法性能，更好地解决一种5自由度机械臂的逆运动学问题。但是，受限于问题本身的复杂性和算法性能，搜寻逆运动学最优解时不可避免地会陷入局部最优。同时，这些改进优化算法对求解能力的重复表现效果并未进行分析，不能满足移动装弹机械臂的实际应用中的要求。

本文以移动装弹机械臂为研究对象，为解决提高其逆运动学求解精度和重复精度，提出了一种多种群灰狼优化算法（Multi-population grey wolf optimization algorithm， MGWO）。MGWO算法通过增加灰狼种群，保证了种群多样性，改善了算法求解后期全局搜索能力减退的问题。同时，为平衡整体与局部的关系，提高局部勘探精度，引入粒子群位置更新机制贴合“狼王”反向指导最差个体过程，并设置随机阈值淘汰最差种群，进一步提高收敛精度。算法经过测试函数的性能验证后，进行移动装弹机械臂的求解仿真，结果表明：与其他算法相比，MGWO算法求解移动装弹机械臂逆运动学问题的精度更高且稳定性更强，验证了本文方法的有效性。

1 正运动学模型

1.1　移动装弹机械臂结构分析

移动装弹机械臂属于军事应用的特种机械臂，是基座安装在车上的8自由度串联机械臂（见图1）：关节1和8为回转关节，关节2、3和7为旋转关节，关节4~6为伸缩关节，多类型的关节使得其工作空间更大。每个关节由液压系统驱动，机械强度十分可靠，由于油缸内置传感器，所以系统可精确反馈关节状态。

1.2　移动装弹机械臂工作流程

在执行自动装弹任务时，移动装弹机械臂固定在2号车平台上，装载运输炮弹的1号车、2号车与放置炮筒和炮弹残骸的具体区域分布，如图2所示，工作流程如下（将炮弹残骸取出并放置于空地处的过程与下述流程相似）：

（1）1号车与2号车平行停靠，通过调整机械臂旋转关节的大小，控制机械臂从原始折叠位置点（即机械臂在运输过程中的折叠储存状态）运行到工作初始位置点（即机械臂在可以开始工作前的初始状态）。

（2）利用测量设备检测两车及炮筒的空间位置，再通过空间位置参数进行装载炮弹过程的轨迹规划。之后，关节1旋转一定角度，使末端执行器转至1号车的炮弹上方，伸展臂架使末端执行器下降，下降至一定高度后，末端执行器夹手张开，臂架前伸，夹手到达合适位置后夹紧炮弹，臂架携带末端执行器及炮弹上升，达到一定高度。

（3）自动控制2号车的机械臂将1号车的炮弹起吊并填装至空地处的炮筒中；然后，臂架回转至炮筒上方；最后，臂架前伸或后缩选择最优位置填充炮弹至炮筒中。

1.3　正运动学建模

针对移动装弹机械臂系统，以其连杆方向为z轴，以纸面平面内与z轴垂直的方向为x轴（两个回转关节为了直观凸显空间关系将坐标轴进行了旋转），再利用右手定则确定y轴，根据上述准则在机械臂各关节处建立空间直角坐标系，如图3所示。

基于齐次变换矩阵法，可以得到相应的位姿变换矩阵如下（规定逆时针为正）：

{{}_{}{}^{\mathrm{0}}}_{\mathrm{1}}=\left[\begin{array}{cccc}\mathrm{c}\mathrm{o}\mathrm{s}{\theta }_{\mathrm{1}}& -\mathrm{s}\mathrm{i}\mathrm{n}{\theta }_{\mathrm{1}}& \mathrm{0}& \mathrm{0}\\ \mathrm{s}\mathrm{i}\mathrm{n}{\theta }_{\mathrm{1}}& \mathrm{c}\mathrm{o}\mathrm{s}{\theta }_{\mathrm{1}}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{1}& l_{\mathrm{1}}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{1}\end{array}\right]（1）

{{}_{}{}^{\mathrm{1}}}_{\mathrm{2}}=\left[\begin{array}{cccc}\mathrm{c}\mathrm{o}\mathrm{s}{\theta }_{\mathrm{2}}& -\mathrm{s}\mathrm{i}\mathrm{n}{\theta }_{\mathrm{2}}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& -\mathrm{1}& \mathrm{0}\\ \mathrm{s}\mathrm{i}\mathrm{n}{\theta }_{\mathrm{2}}& \mathrm{c}\mathrm{o}\mathrm{s}{\theta }_{\mathrm{2}}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{1}\end{array}\right]（2）

{{}_{}{}^{\mathrm{2}}}_{\mathrm{3}}=\left[\begin{array}{cccc}\mathrm{c}\mathrm{o}\mathrm{s}{\theta }_{\mathrm{3}}& -\mathrm{s}\mathrm{i}\mathrm{n}{\theta }_{\mathrm{3}}& \mathrm{0}& l_{\mathrm{2}}\\ \mathrm{s}\mathrm{i}\mathrm{n}{\theta }_{\mathrm{3}}& \mathrm{c}\mathrm{o}\mathrm{s}{\theta }_{\mathrm{3}}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{1}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{1}\end{array}\right]（3）

{{}_{}{}^{\mathrm{3}}}_{\mathrm{4}}=\left[\begin{array}{cccc}\mathrm{1}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& -\mathrm{1}& -l_{\mathrm{3}}-d_{\mathrm{4}}\\ \mathrm{0}& \mathrm{1}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{1}\end{array}\right]（4）

{{}_{}{}^{\mathrm{4}}}_{\mathrm{5}}=\left[\begin{array}{cccc}\mathrm{1}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{1}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{1}& d_{\mathrm{5}}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{1}\end{array}\right]（5）

{{}_{}{}^{\mathrm{5}}}_{\mathrm{6}}=\left[\begin{array}{cccc}\mathrm{1}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{1}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{1}& d_{\mathrm{6}}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{1}\end{array}\right]（6）

{{}_{}{}^{\mathrm{6}}}_{\mathrm{7}}=\left[\begin{array}{cccc}\mathrm{c}\mathrm{o}\mathrm{s}{\theta }_{\mathrm{7}}& \mathrm{0}& \mathrm{s}\mathrm{i}\mathrm{n}{\theta }_{\mathrm{7}}& \mathrm{0}\\ \mathrm{0}& \mathrm{1}& \mathrm{0}& \mathrm{0}\\ -\mathrm{s}\mathrm{i}\mathrm{n}{\theta }_{\mathrm{7}}& \mathrm{0}& \mathrm{c}\mathrm{o}\mathrm{s}{\theta }_{\mathrm{7}}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{1}\end{array}\right]（7）

{{}_{}{}^{\mathrm{7}}}_{\mathrm{8}}=\left[\begin{array}{cccc}\mathrm{c}\mathrm{o}\mathrm{s}{\theta }_{\mathrm{8}}& \mathrm{s}\mathrm{i}\mathrm{n}{\theta }_{\mathrm{8}}& \mathrm{0}& \mathrm{0}\\ -\mathrm{s}\mathrm{i}\mathrm{n}{\theta }_{\mathrm{8}}& \mathrm{c}\mathrm{o}\mathrm{s}{\theta }_{\mathrm{8}}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{1}& l_{\mathrm{4}}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{1}\end{array}\right]（8）

式中：{\theta }_{\mathrm{1}}\sim {\theta }_{\mathrm{8}}为8个关节参数；{}^{i-\mathrm{1}}{\mathit{T}}_i为坐标系\left\{i-\mathrm{1}\right\}到坐标系\left\{i\right\}的齐次变换矩阵。

将以上各相邻关节坐标系之间的齐次变换矩阵依次相乘即可得到移动装弹机械臂的正运动学模型：

{}^{\mathrm{0}}{\mathit{T}}_{\mathrm{8}}{=}^{\mathrm{0}}{{\mathit{T}}_{\mathrm{1}}}^{\mathrm{1}}{{\mathit{T}}_{\mathrm{2}}}^{\mathrm{2}}{{\mathit{T}}_{\mathrm{3}}}^{\mathrm{3}}{{\mathit{T}}_{\mathrm{4}}}^{\mathrm{4}}{{\mathit{T}}_{\mathrm{5}}}^{\mathrm{5}}{{\mathit{T}}_{\mathrm{6}}}^{\mathrm{6}}{{\mathit{T}}_{\mathrm{7}}}^{\mathrm{7}}{\mathit{T}}_{\mathrm{8}}（9）

机械臂的末端位姿\mathit{T}为：

\mathit{T}=\left[\begin{array}{cccc}{n}_x& o_x& a_x& p_x\\ n_y& o_y& a_y& p_y\\ n_z& o_z& a_z& p_z\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{1}\end{array}\right]=\left[\begin{array}{cc}{\mathit{R}}_{\mathrm{3}\times \mathrm{3}}& {\mathit{P}}_{\mathrm{3}\times \mathrm{1}}\\ {\mathrm{0}}_{\mathrm{1}\times \mathrm{3}}& \mathrm{1}\end{array}\right]={{}_{}{}^{\mathrm{0}}\mathit{T}}_{\mathrm{8}}（10）

式中：\mathit{T}为齐次变换矩阵；{\mathit{R}}_{\mathrm{3}\times \mathrm{3}}为移动装弹机械臂末端的姿态矩阵，其3个列向量{\left(n_x,n_y,n_z\right)}^{\mathrm{{\rm T}}}、{\left(o_x,o_y,o_z\right)}^{\mathrm{{\rm T}}}、{\left(a_x,a_y,a_z\right)}^{\mathrm{{\rm T}}}依次为法向量、滑动向量和接近向量，表示末端姿态的坐标系中3个轴对应的单位方向向量在参考坐标系中的投影；p_x\mathrm{、}{p}_y\mathrm{、}{p}_z为坐标原点在参考系中的位置；{\mathit{P}}_{\mathrm{3}\times \mathrm{1}}={\left(p_x,p_y,p_z\right)}^{\mathrm{{\rm T}}}为移动装弹机械臂末端位置矢量。

1.4　工作空间分析

分析移动装弹机械臂的工作空间前，需要明确各个关节参数的约束条件。由于军用车辆在行驶过程中需要对车体所载物体进行遮蔽，所以运输过程中的机械臂为折叠状态，如图4（a）所示。为了保证机械臂在伸展及工作过程中不与车体平台所在的水平平面发生碰撞，此时仅考虑连杆1~3的位置即可。故而对关节角2和3通过具体的动态位置分析，找到图4（b）和（c）中实线描述的临界状态，并结合表1中连杆1~3的长度参数和三角几何关系，对角度做了限制，图中的临界状态满足：\mathrm{c}\mathrm{o}\mathrm{s}{\theta }_{\mathrm{2}}=l_{\mathrm{1}}/l_{\mathrm{2}}，\mathrm{s}\mathrm{i}\mathrm{n}{\theta }_{\mathrm{3}}=l_{\mathrm{1}}/l_{\mathrm{3}}，代入连杆长度后可以得到：{\theta }_{\mathrm{2}}\ge \mathrm{45.89}°，{\theta }_{\mathrm{3}}\le \mathrm{35.45}°。

假设机械臂工作时末端执行器不会低于车载平台所在的平面，由上述分析可以得到各关节角参数的定义区间：

\begin{array}{cc}{\theta }_{\mathrm{1}}\in \left.\left(-{\mathrm{180}}^{\circ },{\mathrm{180}}^{\circ }\right.\right],& {\theta }_{\mathrm{2}}\in \left[-\mathrm{45}.{\mathrm{89}}^{\circ },\mathrm{45}.{\mathrm{89}}^{\circ }\right],\end{array}

\begin{array}{cc}{\theta }_{\mathrm{3}}\in \left[-\mathrm{35}.{\mathrm{45}}^{\circ },\mathrm{35}.{\mathrm{45}}^{\circ }\right],& d_{\mathrm{4}}\in \left[\mathrm{0,0.5}\right],\end{array}

\begin{array}{cc}{d}_{\mathrm{5}}\in \left[\mathrm{0,0.5}\right],& d_{\mathrm{6}}\in \left[\mathrm{0,0.5}\right],\end{array}

\begin{array}{cc}{\theta }_{\mathrm{7}}\in \left.\left(-{\mathrm{180}}^{\circ },{\mathrm{180}}^{\circ }\right.\right],& {\theta }_{\mathrm{8}}\in \left.\left(-{\mathrm{180}}^{\circ },{\mathrm{180}}^{\circ }\right.\right].\end{array}

通过各关节参数的定义区间，按几何关系即可得机械臂的工作范围，但求解过程烦琐且难以保证结果的正确性。而借助计算机数据处理能力的蒙特卡洛法是基于大数定律的简单数值计算方法，已被广泛应用于太空维修机械臂^［10］、拟人机械臂^［11］等机器人的工作空间求解过程中。为了应用蒙特卡洛法求取移动装弹机械臂的工作空间，从式（10）中提取出末端在空间坐标系位置矢量如下：

\mathit{P}={\left[\begin{array}{ccc}{p}_x& p_y& p_z\end{array}\right]}^{\mathrm{{\rm T}}}（11）

在定义空间内生成关节变量的随机值：

q_i=q_{i\mathrm{m}\mathrm{i}\mathrm{n}}+\mathrm{r}\mathrm{a}\mathrm{n}\mathrm{d}\times \left(q_{i\mathrm{m}\mathrm{a}\mathrm{x}}-q_{i\mathrm{m}\mathrm{i}\mathrm{n}}\right),i=\mathrm{1,2},\cdots ,\mathrm{8}（12）

式中：q_i为第i个关节值；q_{i\mathrm{m}\mathrm{a}\mathrm{x}}和q_{i\mathrm{m}\mathrm{i}\mathrm{n}}分别为第i个关节参数的上下界；rand函数生成\left[\mathrm{0},\mathrm{1}\right]随机数。

将式（12）代入式（11）中，即可得到一个机械臂可达的位置坐标。重复N=\mathrm{10}\mathrm{000}次该随机过程，即可获得接近实际的移动装弹机械臂工作空间，如图5所示，其中原点处是机械臂底座固定点位。

2 优化目标函数

2.1　逆运动学问题分析

为了应用群体智能算法求解移动装弹机械臂的运动学逆解，首先需将逆运动学问题转换成等效的最优问题，即以末端执行器位姿误差最小为优化目标建立适应度函数，再通过算法迭代获得最优解。

2.2　适应度函数

设移动装弹机械臂末端期望位姿{\mathit{T}}^{\mathrm{*}}为：

{\mathit{T}}^{\mathrm{*}}=\left[\begin{array}{cccc}{n}_x^{\mathrm{*}}& o_x^{\mathrm{*}}& a_x^{\mathrm{*}}& p_x^{\mathrm{*}}\\ n_y^{\mathrm{*}}& o_y^{\mathrm{*}}& a_y^{\mathrm{*}}& p_y^{\mathrm{*}}\\ n_z^{\mathrm{*}}& o_z^{\mathrm{*}}& a_z^{\mathrm{*}}& p_z^{\mathrm{*}}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{1}\end{array}\right]=\left[\begin{array}{cc}{\mathit{R}}_{\mathrm{3}\times \mathrm{3}}^{\mathrm{*}}& {\mathit{P}}^{\mathrm{*}}\\ {\mathrm{0}}_{\mathrm{1}\times \mathrm{3}}& \mathrm{1}\end{array}\right]（13）

将算法搜寻到的候选解，即一组关节角参数x=\left[{\theta }_{\mathrm{1}},{\theta }_{\mathrm{2}},{\theta }_{\mathrm{3}},d_{\mathrm{4}},d_{\mathrm{5}},d_{\mathrm{6}},{\theta }_{\mathrm{7}},{\theta }_{\mathrm{8}}\right]代入式（9），即可得到末端的实际位姿\mathit{T}。理想的候选解可使{\mathit{T}}^{\mathrm{*}}=\mathit{T}，此时的候选解就是机械臂的逆运动学解。为了判断每次迭代后所产生候选解的优劣性，设计适应度函数如下：

fit\left(X^{\mathrm{*}}\right)=\mathrm{m}\mathrm{i}\mathrm{n}(\mu (\left|p_x^{\mathrm{*}}-p_x\right|+\left|p_y^{\mathrm{*}}-p_y\right|+\left|p_z^{\mathrm{*}}-p_z\right|)+(\left|n_x^{\mathrm{*}}-n_x\right|+\left|n_y^{\mathrm{*}}-n_y\right|+

\left|n_z^{\mathrm{*}}-n_z\right|+\left|o_x^{\mathrm{*}}-o_x\right|+\left|o_y^{\mathrm{*}}-o_y\right|+

\left.\left|o_z^{\mathrm{*}}-o_z\right|+\left|a_x^{\mathrm{*}}-a_x\right|+\left|a_y^{\mathrm{*}}-a_y\right|+\left|a_z^{\mathrm{*}}-a_z\right|\right)（14）

式中：\mu为调节姿态误差和位置误差之间的调节因子。

以期望位姿与实际位姿的误差最小为优化目标，即式（14）的值越小，得到的逆运动学解越精确。实际上，姿态除了用旋转矩阵\mathit{R}表示，还可通过只需要3个变量的欧拉角来表示。故而，在应用Matlab进行机械臂末端位姿求解过程中，将齐次矩阵\mathit{T}转化为欧拉角（\varsigma \mathrm{、}\sigma \mathrm{、}\tau）以减少变量，从而优化目标函数可以简化为：

fit\left(X^{\mathrm{*}}\right)=\mathrm{m}\mathrm{i}\mathrm{n}\left(\mu \left(\left|p_x^{\mathrm{*}}-p_x\right|+\left|p_y^{\mathrm{*}}-p_y\right|\right.\right.+

\left.\left|p_z^{\mathrm{*}}-p_z\right|\right)+\left.\left(\left|{\varsigma }_{}^{\mathrm{*}}-\varsigma \right|+\left|{\sigma }_{}^{\mathrm{*}}-\sigma \right|+\left|{\tau }_{}^{\mathrm{*}}-\tau \right|\right)\right)（15）

式中：\varsigma \mathrm{、}\sigma \mathrm{、}\tau为Z-Y-Z形式的欧拉角；\left|{\varsigma }_{}^{\mathrm{*}}-\varsigma \right|+\left|{\sigma }_{}^{\mathrm{*}}-\sigma \right|+\left|{\tau }_{}^{\mathrm{*}}-\tau \right|为姿态误差；\left|p_x^{\mathrm{*}}-p_x\right|+\left|p_y^{\mathrm{*}}-p_y\right|+\left|p_z^{\mathrm{*}}-p_z\right|为位置误差。

3 算法原理

3.1　多种群灰狼算法

GWO算法是一种生物启发式群智能算法^［12］，灰狼个体在解空间内，以灰狼群体捕猎的数学模型为依据，通过不断更新自身位置找到全局的最优解。该算法具有结构简单、需要调节的参数少、容易实现等优点，且存在自适应调整的收敛因子以及信息反馈机制，能够在局部寻优与全局搜索之间实现平衡，因此，在对最优问题的求解精度和收敛速度方面都有良好的性能^［13］。

由于GWO算法为随机生成初始种群，所以缺乏种群的多样性；且其搜索机制会导致后期收敛速度慢，易陷入局部最优。针对上述问题，在初始化过程中增加灰狼的种群个数，并增设迭代过程中个体的位置更新标准，同时引入PSO算法提升性能，改进框架原理如图6所示，具体改进如下：

（1）通过提升种群多样性，改善GWO算法在处理高维度复杂问题时易陷入局部最优的问题。

在完成灰狼群体的初始化以后，将整个种群分裂为m个子种群，每个子种群按照基本灰狼算法中的个体位置更迭方式进行猎物的搜寻。每个子种群在每次搜寻完成之后都需要将该次的搜寻情况互通，从所有子种群的头狼{\alpha }_i\left(i=\mathrm{1},\cdots ,m\right)中产生整体种群里最优秀的头狼，将其记为狼王I_{\alpha }：

I_{\alpha }=\mathrm{m}\mathrm{i}\mathrm{n}\{{\alpha }_{\mathrm{1}},{\alpha }_{\mathrm{2}},\cdots ,{\alpha }_m\}（16）

式中：，{\alpha }_i\left(i=\mathrm{1},\cdots ,m\right)为每个子种群中适应度最好的头狼；n为整个种群中的灰狼个体数目。

（2）采用狼王I_{\alpha }反向指导子种群中适应度最差灰狼的策略，提高GWO算法后期的收敛速度。

在灰狼种群完成了多种群的裂变后，考虑的核心内容是各个子种群的信息交流问题。多种群的引入能够提高群体的多样性，在初期寻优过程中会有明显提升，但是同时每个种群的个体数目会变少，如果不做处理，后期的收敛速度会下降。

具体地，由于PSO算法有操作模型简单，参数调节量少且高效的优势，又因其思路与狼王反向指导过程有着极大的贴合性——PSO算法中每个粒子的移动只跟历史最佳粒子与当前最佳粒子有关，故而将PSO算法的思想应用于具体的狼王引导中。

PSO算法中将第i粒子的位置记为一个D维的向量{\mathit{X}}_i，其移动速度记为V_i：

\begin{array}{cc}{\mathit{X}}_i=\left(x_{i\mathrm{1}},x_{i\mathrm{2}},\dots ,x_{iD}\right),& i=\mathrm{1,2},\dots ,n\end{array}（17）

\begin{array}{cc}{V}_i=\left(v_{i\mathrm{1}},v_{i\mathrm{2}},\dots ,v_{iD}\right),& i=\mathrm{1,2},\dots ,n\end{array}（18）

粒子的更新规则为：

v_{id}=w\times v_{id}+c_{\mathrm{1}}{r}_{\mathrm{1}}\left(p_{id}-x_{id}\right)+c_{\mathrm{2}}{r}_{\mathrm{2}}\left(p_{gd}-x_{id}\right)（19）

x_{id}=x_{id}+v_{id}（20）

式中：p_{id}为个体的最优解；p_{gd}为种群的历史最优解。

狼王适应度值作为整体种群的最优解，替代式中：p_{gd}子种群中\alpha狼的个体适应度值作为个体最优解，替代其中p_{id}，然后基于如下的反向指导规则，对每个子种群中适应度最低的个体进行位置更新：

v_d^k=w\times v_d^k+c_{\mathrm{1}}{r}_{\mathrm{1}}\left(I_{\alpha }-x_d^k\right)+c_{\mathrm{2}}{r}_{\mathrm{2}}\left(I_k-x_d^k\right),

k=\mathrm{1,2},\dots ,m（21）

x_d^k=x_d^k+v_d^k（22）

式中：v_d^k为第k个种群适应度值最差的狼的速度；I_{\alpha }为狼王的适应度值；I_k为第k个种群中\alpha狼的个体适应度值。

（3） 设置随机重组阈值r，有效避免GWO算法陷入局部最优。

r=\gamma G_{\mathrm{t}}（23）

式中：\gamma为\left(\mathrm{0,1}\right)的比例系数；G_{\mathrm{t}}为总进化代数。

在寻优进程中，对每次迭代后筛选出的狼王适应度值进行标记操作，并判断下次迭代后狼王适应度是否发生变化：若变化，则标记清零并重新计数；若未发生变化，且标记次数大于设置的随机重组阈值r时，随机重组适应度最差的子种群。

综上，MGWO算法通过对比\alpha狼的适应度值来竞争狼王，没有破坏其他的头狼，增加了种群的多样性。同时，用狼王指导的方式更新子种群的最劣个体，保证了算法的寻优速度。最后，通过设置随机阈值的方式改善了算法易陷入局部最优的问题。

3.2　基于MGWO算法的逆运动学求解

MGWO算法求解移动装弹机械臂逆运动学解的流程如图7所示，改进后的新增步骤通过色彩填充区分，具体求解步骤如下：

Step 1　设定机械臂的期望位姿，在机械臂工作空间内随机取一组关节角参数：

x^{\mathrm{*}}=\left[{{\theta }_{\mathrm{1}}}^{\mathrm{*}},{{\theta }_{\mathrm{2}}}^{\mathrm{*}},{{\theta }_{\mathrm{3}}}^{\mathrm{*}},{d_{\mathrm{4}}}^{\mathrm{*}},{d_{\mathrm{5}}}^{\mathrm{*}},{d_{\mathrm{6}}}^{\mathrm{*}},{{\theta }_{\mathrm{7}}}^{\mathrm{*}},{{\theta }_{\mathrm{8}}}^{\mathrm{*}}\right]（24）

将式（24）代入式（9）得到{\mathit{T}}^{\mathrm{*}}，即设定的期望位姿，此时适应度函数f\left(x^{\mathrm{*}}\right)=\mathrm{0}，求取过程中需要的具体连杆参数见表1。

Step 2　设置GWO算法中的基本系数A\mathrm{、}C\mathrm{、}a，以及改进算法中新增的阈值系数\gamma和阈值计数值r等。设定迭代次数t为0，灰狼种群的初始位置也即移动装弹机械臂各关节角的初始化：X_{\mathrm{0}}=\left\{x_{\mathrm{1}},x_{\mathrm{2}},\cdots ,x_n\right\}，x_i=[{\theta }_{\mathrm{1}},{\theta }_{\mathrm{2}},{\theta }_{\mathrm{3}},d_{\mathrm{4}},d_{\mathrm{5}},d_{\mathrm{6}},{\theta }_{\mathrm{7}},{\theta }_{\mathrm{8}}],i=\mathrm{1,2},\cdots ,n。

Step 3　设定子种群个数m，将整个灰狼种群（n匹灰狼）分裂为m个子种群，以增加灰狼种群的多样性从而提升算法后期的全局搜索能力。

Step 4　将各个子种群中灰狼的空间位置代入适应度函数中计算适应度值，计算结果比较后可以得到每个子种群最优的3匹灰狼：\alpha狼、\beta狼和\delta狼。最终，可以得到m组领导层灰狼，把m个子种群中的\alpha狼作为初始化关节角的最优候选解集合\left\{{\alpha }_{\mathrm{1}},{\alpha }_{\mathrm{2}},\cdots {\alpha }_m\right\}。同时，标记好每个子种群中最差的灰狼个体，作为此时最差的候选解。

Step 5　各子种群直接进行交流，竞选出整个群体中适应度值最好的狼王，根据式（16），从最优候选解集合\left\{{\alpha }_{\mathrm{1}},{\alpha }_{\mathrm{2}},\cdots ,{\alpha }_m\right\}中选出的最优候选解I_{\alpha }，并记录阈值计数值r——狼王的适应度值不发生变化的次数。

Step 6　用狼王根据式（21）（22）指导最差的灰狼移动。再通过判断Step 4中的计数值与随机阈值的关系，依据改进规则对适应度最差的子种群进行重组，重组后阈值计数值r清零。

Step 7　应用基本GWO算法中的个体位置移动公式分别对各个子种群的个体位置进行更新。

Step 8　判断迭代次数t是否达到初始化设置的值，如果达到了执行下一步；反之返回Step 4。

Step 9　输出移动装弹机械臂逆运动学的最优解，即最后一次迭代后的关节角（Best position）、该关节角所对应的适应度函数值（Best cost）。

4 仿真分析

4.1　基准函数测试与分析

为了验证MGWO算法的性能，选择单峰值函数F₁和F₂测试算法的求解精度与速度；多峰值函数F₃和F₄测试算法的全局寻优能力，如表2所示。设置算法的最大迭代次数为200，种群规模为30。分别应用MGWO算法和标准GWO算法对4个基准测试函数进行50次单独的极值求解，得到的极值搜索迭代过程如图8所示，相关统计结果如表3所示。从图8可以看出，与标准的GWO算法相比，MGWO算法在迭代次数相同时达到了更小的适应度值，且达到同一个适应度值时迭代次数更少。

由表3可知，基准测试函数F₁~F₃的收敛效果改进显著。函数F₁和F₂最小值求解的最差适应度值比改进前最优适应度值降低了3个数量级。应用函数F₄时，平均适应度反映出MGWO算法在50次求解中发挥更稳定，整体求解效果更好。由此可见，MGWO算法的收敛速度更快，收敛精度更高，全局搜索能力也得到了提升。

4.2　逆运动学求解仿真结果

为了验证MGWO算法对移动装弹机械臂逆运动学问题求解的有效性，将其与GWO^［12］、VAGWO^［14］和PGWO进行了对比仿真分析。设置算法的最大迭代次数为200，种群规模为30，适应度函数中的\lambda =\mathrm{1}。

给定移动装弹机械臂的一个期望位姿x^{\mathrm{*}}=\left[\mathrm{\pi }/\mathrm{4},\mathrm{\pi }/\mathrm{6},\mathrm{\pi }/\mathrm{4,0.8,1},\mathrm{1},\mathrm{\pi }/\mathrm{6},-\mathrm{\pi }/\mathrm{4}\right]，根据3.2节中的步骤，用4种算法分别进行50次重复独立求解，得到适应度值平均收敛曲线和重复定位适应度函数曲线，分别如图9和图10所示，性能指标的统计见表4。

根据表4中的数据，结合图9~11可知，MGWO算法在求解的搜索迭代结束时适应度值最小。VAGWO和PGWO的求解精度高时有优势，但整体上表现不稳定。MGWO相比GWO的求解精度有提升，虽然单次最优适应度值非最佳，但最差和平均适应度值均最小，且方差远小于其他算法，说明其求解该问题的稳定性很好。

以上仅对一个位姿进行了定位，算法求解能力的体现不具有普遍性。因此，考虑到机械臂的整个工作空间，通过1.4节工作空间分析中的各关节参数范围及随机点的生成方式，得到50个随机工作点，再应用4种算法分别对这些点进行求解，仿真得到随机定位点的适应度函数曲线，如图11所示。MGWO算法的求解精度优于GWO和PGWO算法，且求解稳定性好。VAGWO算法在某些位姿的求解中精度最佳，但对其他点的求解精度又过差，不能满足实际的定位需求。综上，MGWO算法求解的效果更能满足机械臂的实际工程应用需求。

求解的精度与重复精度的提升，不可避免地牺牲了时间成本，而考虑到移动装弹机械臂的实时性能在实际作业时也十分重要，因此，结合其填弹过程多为重复性操作的情况，采用离线求解、在线操作的方式，规避时间上的性能损失。

5 结束语

本文针对移动装弹机械臂逆运动学求解问题，融合了丰富灰狼种群策略、狼王反向指导策略和跟踪狼王策略，提出了MGWO算法。首先通过结构分析建立正运动学模型；应用MGWO算法，以机械臂末端位姿误差最小为优化目标进行迭代寻优，完成了逆运动学求解。仿真结果中重复求解的平均适应度值达到了\mathrm{0.012}\mathrm{2}，方差仅为\mathrm{8.036}\mathrm{2}\times {\mathrm{10}}^{-\mathrm{5}}，表明该算法的寻优能力得到了提升，求解精度更高且稳定性更好，可以有效解决移动装弹机械臂的逆运动学问题。针对本文研究的不足之处——算法的实时性能与求解精度的平衡问题，以及求解精度的进一步提升，考虑结合VAGWO求解部分实例的高精度与MGWO算法的强稳定性，得到更优异的求解性能与时间性能。

参考文献

[1] 霍希建, 刘伊威, 姜力, 具有关节限位的7R仿人机械臂逆运动学优化[J]. 吉林大学学报: 工学版, 2016, 46(1): 213-220.

[2] Huo Xi-jian, Liu Yi-wei, Jiang Li, et al. Inverse kinematics optimization of a 7R humanoid robot arm with joint limits [J]. Journal of Jilin University(Engineering and Technology Edition ), 2016, 46(1): 213-220.

[3] 石建平, 刘鹏, 陈冬云. 基于改进粒子群优化算法的冗余机械臂逆运动学求解[J]. 机械传动, 2021, 45(2): 69-75.

[4] Shi Jian-ping, Liu Peng, Chen Dong-yun. Inverse kinematics of redundant manipulator based on improved particle swarm optimization algorithm[J]. Journal of Mechanical Transmission, 2021, 45(2): 69-75.

[5] 冷舒, 吴克, 居鹤华. 机械臂运动学建模及解算方法综述[J]. 宇航学报, 2019, 40(11): 1262-1273.

[6] Leng Shu, Wu Ke, Ju He-hua, Kinematics modeling and calculating method of mechanical arm review [J]. Journal of Aerospace, 2019, 40(11): 1262-1273.

[7] Alkayyali M, Tutunji T A. PSO-based algorithm for inverse kinematics solution of robotic arm manipulators[C]∥20th International Conference on Research and Education in Mechatronics (REM),IEEE, Wels, Austria, 2019: 1-6.

[8] Ghosh A, Singh O, Ray A K. Inverse kinematic solution of a 7 DOF robotic manipulator using boundary restricted particle swarm optimization[J]. IFAC-PapersOnLine, 2022, 55(1): 101-105.

[9] Liu Y, Xi J, Bai H, et al. A general robot inverse kinematics solution method based on improved PSO algorithm[J]. IEEE Access, 2021, 9: 32341-32350.

[10] Gao R. Inverse kinematics solution of robotics based on neural network algorithms[J]. Journal of Ambient Intelligence and Humanized Computing, 2020, 11(12): 6199-6209.

[11] Nguyen T, Bui T, Pham H. Using proposed optimization algorithm for solving inverse kinematics of human upper limb applying in rehabilitation robotic[J]. Artificial Intelligence Review, 2022, 55(1): 679-705.

[12] Ahmed E S, Elhosseini M A, Haikal A Y. A new ABC variant for solving inverse kinematics problem in 5 DOF robot arm[J]. Applied Soft Computing, 2018, 73: 24-38.

[13] 杨凯, 黄晋英. 一种8自由度空间机械臂运动学及工作空间分析[J]. 机械传动, 2021, 45(3): 147-152.

[14] Yang Kai, Huang Jin-ying. Kinematics and workspace analysis of a 8-DOF space manipulator[J]. Mechanical Transmission, 2021, 45(3): 147-152.

[15] 张清松, 段帅臣, 夏热. 基于Matlab的拟人机械臂工作空间分析及仿真[J]. 机械传动, 2020, 44(12): 99-105.

[16] Zhang Qing-song, Duan Shuai-chen, Xia Re. Workspace analysis and simulation of anthropomorphic robotic arm based on Matlab[J]. Mechanical Transmission, 2020, 44(12): 99-105.

[17] Mirjalili S, Mirjalili S M, Lewis A. Grey wolf optimizer[J]. Advances in Engineering Software, 2014, 69: 46-61.

[18] 张晓凤, 王秀英. 灰狼优化算法研究综述[J]. 计算机科学, 2019, 46(3): 30-38.

[19] Zhang Xiao-feng, Wang Xiu-ying, Review of grey wolf optimization algorithm[J]. Computer Science, 2019, 46(3): 30-38.

[20] Rezaei F, Safavi H R, Abd Elaziz M, et al. An enhanced grey wolf optimizer with a velocity-aided global search mechanism[J]. Mathematics, 2022, 10(3):No.351.

基金资助

吉林省教育厅科学研究重点项目(JJKH20221007KJ)

吉林省科技发展计划项目(20220101120JC)