顾海华 ，胡兴柳 ，毛军辉
1.盐城工学院 电气工程学院，江苏 盐城 224002；
2.金陵科技学院 智能科学与控制工程学院，江苏 南京 211169；
3.南京邮电大学 自动化学院、人工智能学院，江苏 南京 210003
随着社会与科学技术的不断发展，工业机器人已经成为生产制造和科学研究领域的热点。研究机器人首先要对机器人进行运动学分析，而运动学分析通常分为正运动学分析与逆运动学分析。机器人的正运动学分析是已知关节驱动力或者力矩，求解机器人末端执行器的位姿，如标准D-H参数法与改进的D-H参数法等；
逆运动学分析为已知末端执行器位姿，求解对应关节产生的位移，如Newton-Euler 法、Lagrange 法、Kane法等。这些方法被广泛运用到众多场景中，目前，在MATLAB 与ADAMS 中进行联合仿真是有效分析机器人运动特性的一种方式［1］，其不但可以保证轨迹规划的精确性，也可测试系统的各项特性［2］。
一般工业机器人均由转动副和移动副组成，6R 机器人作为目前最为广泛使用的工业机器人，其6 个关节均是由旋转副组成。图1 为6 自由度机器人的KUKA KR16 操作臂，可实现3 个空间位置和3个控件姿态的任意运动。
图1 KUKA KR16操作臂Fig. 1 KUKA KR16 manipulator
在针对KUKA KR16 进行运动学分析前，需要建立其数学模型。为便于分析，建立机器人的连杆坐标系，如图2 所示，即以机器人的第i- 1个关节的轴线Zi-1为Z轴，以Zi-1和Zi的公垂线Xi-1为X轴建立坐标系，其中X轴的方向为关节i- 1 指向关节i，而Y轴的Yi-1则按右手法则确定，即Yi-1=Zi-1⋅Xi-1。因此，坐标系中的每个连杆均可通过ai-1、αi-1、di和θi这4个参数进行描述，其中ai-1表示连杆长度，为Zi-1和Zi两关节轴线之间的公垂线长度，αi-1表示连杆扭角，为Zi-1和Zi两关节轴线之间的夹角，di表示相邻两关节之间的线位移，为Xi-1和Xi之间在Zi方向上的长度，θi表示相邻两关节之间的角位移，为Xi-1到Xi绕Zi转过的角度。
图2 连杆坐标系Fig. 2 l ink-pole coordinate system
根据机器人的各项机械特性参数，采用目前使用最为广泛的改进D-H 参数法，建立改进的D-H参数表［3］，如表1所示。
表1 改进的D-H参数表Table 1 Improved D-H parameters table
通过MATLAB中的Link和SerialLink函数建立KUKA KR16机械臂模型［4］。由于采用改进的D-H参数表，MATLAB中Link函数的一般表达式为L= Link([θi,d i,ai-1,αi-1,offset],′modified′)。
根据表1，可以得到Link 函数建立机械臂的命令为：
在MATLAB 中运行以上程序后，即可得到图3所示的机器人模型。
图3 机器人MATLAB模型Fig. 3 Robot MATLAB model
根据表1 和图2、图3 机器人各参数之间的关系［5］，可推导出两个连杆坐标系之间的齐次变换矩阵，如式（1）所示。
将表1 数据带入式（1），可以得出各个关节的变换矩阵，分别为：
根据机器人正运动学原理，机器人的关节变换矩阵可写为：
式中：
根据机器人逆运动学中的反代数法，将机器人各关节变量分离出来，再逐次、逐个地求解出各个关节的关节变量(θ1，θ2，θ3，θ4，θ5，θ6)。具体求解方式如下：
设机器人初始位姿（机器人第1～6 个关节的角度）为θstart=[0，- pi/2， 0， pi， pi， 0]，通过正运动学方程求解，结果如式（10）所示；
通过MATLAB/Robotics Toolbox 工具箱求解，结果如式（11）所示。
比较式（10）、式（11）可以看出，式（10）的二行二列为-1，式（11）的二行二列为1，其余各行各列的数值都相同。二行二列两个数值的不同是因为机器人关节正反转的不同，不影响最终位姿呈现的结果。故正运动学方程推导出的结果和MATLAB 工具箱求解的结果是一致的，证明了DH参数选取和齐次变换矩阵的正确性。
采用MATLAB/Robotics Toolbox 工具箱中的ikine 函数对机器人进行逆运动学计算，得到各关节的初始角度，如式（12）所示。由于运动学逆解存在不唯一解，在比较式（12）与初始位姿θstart=[0，- pi/2， 0， pi， pi， 0]时，可以发现式（12）与初始位姿θstart设定的结果基本一致。
机器人轨迹规划分为关节空间轨迹规划与笛卡尔空间轨迹规划。机器人运动轨迹确定方法很多，为获得更高的平滑性，采用关节空间轨迹规划中的五次多项式插值的方法进行轨迹规划［6］。
设机器人各关节角位移函数为：
式中：s0、s1、s2、s3、s4、s5分别为五次多项式的系数。
对式（13）位移函数求一阶导数，即可得到机器人各关节运动角的速度函数：对位移函数求二阶导数即可得到机器人各关节运动角加速度函数：
为得到唯一的关节运动轨迹曲线，必须对其起始时刻的关节角度、关节角速度、角加速度共计6参数进行约束。
式中：θ0i为第i个关节的初始角度，（°）；
为第i个关节的初始角速度，（°）/s；
为第i个关节的初始角加速度，(°)/s2；
θfi为第i个关节的终止角度，（°）为第i个关节的终止角速度，（°)/s为第i个关节的终止角加速度，（°)/s2；
t0i为第i个关节运动的起始时刻，s；
tfi为第i个关节从初始角度运动至终止角度的时间，s。
将式（16）带入式（13）～式（15），得：
为得到机器人运动轨迹，在MATLAB 仿真前需先确定机器人各关节的初始角度θintial= [0，0，0，0，0，0]、终 止 角 度θtarget= [pi/2，pi/4，-pi/4，pi/3， pi/6，pi]、仿真时间10 s、采样50 次，再根据式（17）及式（13）～式（15），在MATLAB 中仿真即可得出各关节对应的角度、角速度、角加速度的变化曲线，如图4～图6所示。
由图4 可知，图形准确地显示了各关节的起始角度和终止角度；
由图5 可知，图中曲线平滑，无突兀的激增和陡降，说明各关节角速度合适；
由图6可知，图中曲线平稳变化。显然，采用五次多项式插值的方法进行轨迹规划是比较成功的［7］。
图4 各关节角度轨迹规划Fig. 4 Angle trajectory planning of each joint
图5 各关节角速度轨迹规划Fig. 5 Angular velocity trajectory planning of each joint
图6 各关节角加速度轨迹规划Fig. 6 Angular acceleration trajectory planning of each joint
ADAMS 是一款集建模、求解、可视化等技术于一体，在虚拟样机领域被广泛使用的软件。利用ADAMS 软件，用户可以直接建立机械系统的几何模型，或者从其他三维软件（SolidWorks、Pro/E 等）中导入更为精细的模型［8］。本文利用Solid⁃Works 软件建立模型并保存为Parasolid 格式，然后将其导入至ADAMS。需要注意的是，为使得仿真结果更贴近实际，同时还能减小计算量，通常在仿真之前先将细小的、零碎的、不影响仿真结果的部件删除，然后将相关零件通过布尔和合为一体。另外，在模型导入至ADAMS 后，各零部件的材质、属性、质量等参数都会丢失。因此，在模型导入至ADAMS 后，需要为每个零部件重新定义材质，并添加合适的约束。本文将所有材料属性定义为steel，各关节之间均取旋转约束（本文为6R机器人）。
利用MATLAB 软件将五次多项式插值法轨迹规划的数值计算出来，并将其保存为txt 格式，再导入至ADAMS 中；
在ADAMS 中为各约束添加合适的驱动力后，利用ADAMS 中的spline 函数将各关节角度作用至对应的关节即可完成驱动，再利用轨迹工具画出机器人运动轨迹线，如图7～图9 所示。机器人各关节之间的约束关系及驱动关系分别如表2～表3所示。
表2 机器人各关节约束关系Table 2 Constraint relation of each joint of robot
表3 机器人各关节驱动关系Table 3 Driving relationship of each joint of robot
对比图7～图9，可知各曲线的数值大小、弧度等都近似，说明联合仿真比较成功。
图8 机器人各个关节在ADAMS中的角度曲线Fig. 8 Angles of each joint in ADAMS
图9 机器人各关节在MATLAB中的角度曲线Fig. 9 Angles of robot joints in MATLAB
在ADAMS 轨迹规划中，也可使用“一般点驱动”功能对机器人进行轨迹规划［9-10］，操作步骤为：在模型的末端执行器上设置一个标记点，在该标记点上添加“一般点驱动”，在驱动编辑里设置相应的驱动方程。
利用“一般点驱动”功能对常见且较为复杂的圆锥螺旋曲线进行仿真，步骤如下：
（1）确定圆锥螺旋曲线的数学表达式
式中：S为螺旋线导程，mm；
T为运动周期，s；
a为螺旋曲线的螺距，mm；
t为仿真时间，s。
（2）设置相关参数：螺旋线导程S= 30 mm，运动周期T= 2s，螺旋曲线的螺距a= 15 mm，仿真时间t=6 s，步数300步。
（3）选择交互仿真。
仿真结束后在末端执行器上通过轨迹工具即可得到末端执行器的运动轨迹曲线，如图10所示。
图10 机器人末端执行器轨迹Fig. 10 Trajectory diagram of robot end-effector
在ADAMS 中，也可以使用spline 函数进行仿真，方法如下：清除之前设置的点驱动，导入6 个关节的角度spline 函数，给每个关节添加相对应的驱动（驱动函数由spline 函数完成），并设置仿真时间6 s，步数300 步［11］。通过DX、DY、DZ 函数得到末端执行器在X、Y、Z轴方向上的运动曲线，如图11 所示；
通过Ax/Ay/Az Projected Rotation 得到机器人各关节角度变化曲线，如图12 所示；
通过Relative Angular Velocity 得到机器人各关节角加速度变化曲线，如图13所示。
图11 机器人末端执行器运行曲线Fig. 11 Running curve of robot end-effector
图12 机器人各关节角度变化曲线Fig. 12 Angle variation curve of each joint of robot
图13 机器人各关节角速度变化曲线Fig. 13 Angular velocity curve of each joint of robot
将图11～图13 与图10 进行对比，可知仿真结果正确。
在ADAMS 中，通过“一般点驱动”与spline 函数的2次仿真，可以发现仿真目标均为空间形状，仿真后最终所呈现的末端执行器的运行轨迹与机器人运动学分析相符，说明ADAMS 下运动学仿真分析是正确的，验证了从SW 导入至ADAMS的KUKA KR16机器人模型的可靠性。
针对6 自由度工业机器人的运动学与轨迹规划问题，在改进的D-H 参数法基础上，通过矩阵的齐次变换推导出机器人的运动学方程，再通过正逆运动学进行求解；
采用MATLAB/Robotics Toolbox 工具箱对机器人进行正逆运动学求解，验证了运动学方程的正确性，再通过五次多项式插值法进行轨迹规划；
最后在ADAMS 中导入KUKA KR16 机器人的3D 模型，利用ADAMS 与MATLAB 联合仿真技术验证了运动学模型的正确性。