由已有研究成果可知，抑制柔性关节残余振动的方法主要有3种：结构优化、主动控制、输入整形[1-3]。其中，结构优化方法虽然可以抑制机器人的振动，但是由于本文中的物料装填机器人负载质量较大且空间受限，低阶固有频率很难有较大的提升。主动控制方法是振动抑制方法研究的一个热点，最优控制[4]、滑模控制[5]、鲁棒控制[6]等方法已被研究证实可应用于振动控制，但主动控制方法中，多需要增加额外的传感器，并对控制率进行设计，这些都会额外增加成本。区别于传统的输入指令整形方法，许多****通过对关节轨迹进行优化，即优化轨迹指令的方式进行抑振。Rew和Kim[7]通过设计非对称S型曲线，对机器人的残余振动进行抑制。Federico等[8]以机器人振动幅度最小为优化目标，研究了机器人在线轨迹优化的方法。Wisnu和Kazuhiko[9]将三次样条曲线与浮动插值点技术相结合，以减小末端振动，并优化关节运动时间。Guo等[10]以关节应变能最小为优化目标，使用遗传算法对三自由度平面机器人的关节轨迹进行优化。吴明月等[11]通过优化振幅方程中的S型曲线参数，使得硅片传输机器人的低阶振动振幅达到最小值。解则晓等[12]以关节电机近似输出能量函数为目标，通过优化Lame曲线参数，对Delta并联机器人的残余振动进行控制。
本文以设计的物料装填机器人为研究对象，建立了机器人的刚柔耦合动力学模型，并推导得到特定轨迹下的动态响应计算方法。以关节残余弹性势能最小为优化目标，使用最大最小蚂蚁系统(Max-Min Ant System，MMAS)[13]优化插值点时间，间接实现关节轨迹优化，从而达到振动抑制的目的。
1 物料装填机器人总体方案 本文中，输送的圆柱物料为规格一致的生产原料，物料高度约1 085 mm，直径约130 mm，质量约35 kg。由图 1可知，装填过程中，圆柱物料需要在空间中3个自由度上移动，且在1个自由度上转动，同时，考虑到物料加工过程中，加工设备会有一个回撤动作，因此需要物料装填机器人具有一个冗余自由度以躲避物料加工设备。因此，将物料装填机器人设计成图 1(a)所示的“RPP-RR”构型。各个关节依次称作整体转动关节(关节1)、平移关节(关节2)、举升关节(关节3)、调姿关节(关节4)、避障关节(关节5)和末端抓持器，关节处采用谐波减速器以提高传动精度和承载能力。其中，关节1、关节2和关节3用于调整圆柱物料的平面位置，关节2和关节3采用移动关节后，可以使得结构更加紧凑，适合在载具中使用。关节4用于调整圆柱物料姿态，使物料与加工设备进料口轴线一致。关节5提供冗余自由度，用于避开加工设备的回撤动作。
 |
| 图 1 物料装填机器人总体方案 Fig. 1 General scheme of material loading robot |
| 图选项 |
物料装填过程如下：首先解除避障状态，关节5转动90°至抓持位置，同时各关节协调运动，使末端抓持器对准圆柱物料，该过程如图 1(a)所示。随后，末端抓持器抓稳物料，关节1、关节2、关节3、关节4协调运动，将物料送至进料口，并使物料轴线与进料口轴线一致，该过程如图 1(b)所示。末端抓持器松开物料，推送机构将物料推入加工设备，该过程如图 1(c)所示。物料被送入加工设备后，加工设备会有一个快速的回撤动作，如图 1(d)所示，因此，在将物料推送入加工设备后，关节5回归零位，以避免物料加工设备与机器人发生碰撞，即机器人进入避障状态。加工作业完成后，各关节回位，准备进行下一次装填作业。
当圆柱物料到达图 1(c)所示的位置时，垂直于物料输送方向的振动，会对圆柱物料轴线位置产生较大影响；而物料输送方向的振动，对圆柱物料轴线位置的影响相对较小。因此，着重考虑抑制垂直于物料输送方向的振动。图 1(c)中，实线箭头方向为物料最终的输送方向，虚线箭头方向为需要抑制振动的方向，本文中对虚线箭头正反方向的振动均加以考虑。
2 建模与动态特性分析 2.1 动力学建模 考虑到拉格朗日方法具有物理含义明确、形式简洁的优点，本文中使用拉格朗日方法建立物理状态机器人的动力学模型。在动力学建模过程中，将关节等效为均质杆，关节质量等效为集中质量，关节刚度等效为无质量线性弹簧，系统阻尼和关节阻尼简化为小阻尼。经过上述假设后，物料装填机器人可由图 2所示的简化模型表示。
 |
| 图 2 物料装填机器人简化模型 Fig. 2 Simplified model of material loading robot |
| 图选项 |
简化模型中各物理量参数定义及数值见表 1。其中，各关节的等效刚度为仿真结果。
表 1 物料装填机器人物理参数定义 Table 1 Physical parameter definition of material loading robot
| 变量 | 物理意义 | 数值 |
| θ1/rad | 关节1关节角 | |
| d2/m | 关节2关节位移 | |
| d3/m | 关节3关节位移 | |
| θ4/rad | 关节4关节角 | |
| θ5/rad | 关节5关节角 | |
| m1/kg | 关节1质量 | 90.8 |
| m2/kg | 关节2质量 | 64.1 |
| m3/kg | 关节3质量 | 13.4 |
| m4/kg | 关节4质量 | 26.2 |
| m5/kg | 关节5和末端抓持器等效质量 | 46.4 |
| k1/((N·m)·rad-1) | 关节1等效线性刚度 | 15 000 |
| k2/(N·m-1) | 关节2等效线性刚度 | 15 000 |
| k3/(N·m-1) | 关节3等效线性刚度 | 15 000 |
| k4/((N·m)· rad-1) | 关节4等效线性刚度 | 10 000 |
| k5/((N·m)· rad-1) | 关节5等效线性刚度 | 15 000 |
| I1/(kg·m2) | 关节1在坐标系O1x1y1z1中相对于质心的惯量 |  |
| I2/(kg·m2) | 关节2在坐标系O2x2y2z2中相对于质心的惯量 |  |
| I3/(kg·m2) | 关节3在坐标系O3x3y3z3中相对于质心的惯量 |  |
| I4/(kg·m2) | 关节4在坐标系O4x4y4z4中相对于质心的惯量 |  |
| I5/(kg·m2) | 关节5在坐标系O5x5y5z5中相对于质心的惯量 |  |
表选项
基于拉格朗日方程对物料装填机器人进行动力学建模，其动力学模型可表示为
 | (1) |
式中：K为刚度矩阵;
 | (2) |
 |
 |
 |
式中：Iyy2为惯量矩阵I2的yy分量；r1y为关节1的质心矢量r1在y方向的分量。
刚度矩阵K可写成K=KT+KG的形式，其中KT为弹性势能推导得到的刚度矩阵，KG为重力势能推导得到的刚度矩阵。KT和KG表达式如下：
 | (4) |
式中：
 |
 |
其中：g为重力加速度，本文中取g=9.8 m/s2。
阻尼矩阵采用比例阻尼模型[14]进行简化。比例阻尼模型如下:
 | (5) |
式中：α和β分别为质量矩阵M(q)和刚度矩阵K的比例系数。
2.2 动态特性分析 将式(1)写成式(6)所示的形式进行动态特性分析：
 | (6) |
根据模态分析理论[14]，得到系统的固有频率及模态振型:
 | (7) |
式中：ωn为固有角频率；fn为固有频率。
采用表 1所示的系统参数进行仿真，得到物料装填机器人前3阶固有频率特性，如图 3所示。在仿真过程中，设定物料装填机器人初始时刻关节参数向量q0=[0 rad, -0.028 m, -0.321 m, 0 rad, 0 rad]T，举升结束时刻关节参数向量q1=[0 rad, 0.372 m, -0.821 m, π/2 rad, 0 rad]T。举升时间为3 s，且举升过程中，各关节匀速运动。
 |
| 图 3 物料装填机器人固有频率特性 Fig. 3 Characteristics of natural frequencies of material loading robot |
| 图选项 |
仿真结果表明，物料装填机器人的第1阶和第2阶固有频率位于10~30 Hz的低频段，故主要针对1阶和2阶固有频率下的振动进行抑制；物料装填机器人系统的固有频率受关节变量的影响。
3 末端动态响应分析 本节给出举升过程的末端动态响应分析方法。举升过程需要关节2、关节3和关节4协调动作，所以在举升过程中，关节1和关节5保持不动。考虑到物料装填机器人工作在约束空间中，所以预先设置举升过程中的几个插值节点。各个插值节点处关节2、关节3和关节4的相对位移如表 2所示。
表 2 轨迹插值点 Table 2 Trajectory interpolation points
| 插值点序号 | 关节2相对位移/mm | 关节3相对位移/mm | 关节4相对位移/(°) |
| 0 | 0 | 0 | 0 |
| 1 | 195 | 0 | 0 |
| 2 | 305 | 200 | 50 |
| 3 | 405 | 386 | 90 |
表选项
采用五次多项式规划插值点之间的关节轨迹。五次多项式曲线的位移、速度、加速度为
 | (8a) |
 | (8b) |
 | (8c) |
各个插值点处的约束关系为
 | (9) |
式中：i=0, 1, 2, 3；Xi为机器人关节轨迹在插值点i的位置；Vi为机器人关节驱动在插值点i处的速度；Ai为机器人关节驱动在插值点i处的加速度。
根据式(9)的约束条件可以发现，对每个五次多项式，当插值时间t确定时，五次多项式系数均可唯一确定，即可以通过优化插值时间t，优化五次多项式。
关节按五次多项式曲线运动时，关节驱动力F(t)可按照牛顿-欧拉迭代法[15]求解得到，以关节驱动力F(t)作为输入力向量，计算物料装填机器人的动态响应。注意在计算过程中，认为举升动作结束后，关节驱动力仅用于克服重力。
设u为系统的模态矩阵，引入坐标变换方程:
 | (10) |
并注意模态向量的正交性，可得到模态运动方程：
 | (11) |
式中：模态质量矩阵Mr、模态阻尼矩阵Cr、模态刚度矩阵Kr、广义力Qr表达式分别为
 | (12a) |
 | (12b) |
 | (12c) |
 | (12d) |
式(11)表示5个独立方程，且具有与单自由度系统相同的形式，采用杜哈梅积分进行求解，第i个模态坐标的响应为
 | (13) |
式中：f(τ)为第i个模态坐标对应的激励函数。
将式(13)所示的响应方程代入式(10)所示的转换方程中，可得到原系统坐标下的动力学响应方程。进一步地，根据物料装填机器人正向运动学方程[15]，可计算得到的机器人末端误差的动态变化。
4 轨迹优化 采用第3节中的末端动态响应分析方法，可以得到特定轨迹下，物料装填机器人末端的动态响应。因为机器人采用了串联结构，所以需要同时优化关节2、关节3和关节4的轨迹，以抑制机器人末端的残余振动。MMAS是应用较为成功的多变量优化算法之一，目前已有深入的研究[13]。MMAS与传统的蚁群算法不同，体现在以下几个方面：①只有迄今最优游历的精英蚂蚁或迭代最优游历的蚂蚁，才被允许释放信息素。②信息素的取值范围和初值均受到限制，以使算法在迭代初期具有更多的探索性且防止出现停滞现象。③出现停滞或局部最优解时，重新初始化边上的信息素。因此，与传统的蚁群算法相比，MMAS更适合于求解多变量优化问题。
本节中，使用MMAS对关节轨迹进行优化。优化变量取各个插值节点之间的运动时间，优化目标为装填动作结束时，物料装填机器人关节具有的残余弹性势能最小。
优化过程中，各个插值点间运动时间的取值范围如表 3所示。其中，插值时间t的最小值根据关节位移量和关节最高运动速度求得。表中：t1表示插值点0与插值点1之间的运动时间，t2表示插值点1与插值点2之间的运动时间，t3表示插值点2与插值点3之间的运动时间。
表 3 插值时间取值范围 Table 3 Value range of interpolation time?
| 运动时间 | 最大值 | 最小值 |
| t1 | 2 | 0.5 |
| t2 | 2 | 0.5 |
| t3 | 2 | 0.5 |
表选项
仿真得到插值时间与残余振动导致的末端误差幅值的关系，如图 4所示。图中：x方向为搬运动作结束时，坐标系O4x4y4z4中的x4轴方向，z方向为搬运动作结束时，坐标系O4x4y4z4中的z4轴方向。
 |
| 图 4 插值时间与末端误差幅值的关系 Fig. 4 Relation between interpolation time and amplitude of end-effector error |
| 图选项 |
为便于后续优化，定义Ai为装填动作结束时第i个关节的振幅，定义如式(14)所示的残余弹性势能做为优化指标：
 | (14) |
MMAS的信息素更新公式为
 | (15) |
式中：ρ为信息素挥发率。
 | (16a) |
 | (16b) |
式中：Q为常数；ω为权重参数；Le与Lib分别为精英蚂蚁与迭代最优蚂蚁所完成的游历长度。应用MMAS算法优化机器人关节轨迹的流程，如图 5所示。
 |
| 图 5 MMAS流程 Fig. 5 Flowchart of MMAS |
| 图选项 |
优化中，蚁群数ma=20，初始信息素浓度τ0=1.2，信息素挥发率ρ=0.1，信息素浓度影响因子α1=1，启发式信息值影响因子β1=2.3。常数Q=1，权重参数ω=50。优化前后参数的对比如表 4所示。
表 4 优化前后参数的对比 Table 4 Comparison of parameters before and after optimization
| 参数 | 优化前 | 优化后 |
| 运动时间t1/s | 2.04 | 1.77 |
| 运动时间t2/s | 1.96 | 1.23 |
| 运动时间t3/s | 1.02 | 1.17 |
| 时间和t/s | 5.02 | 4.17 |
| 残余弹性势能Ek/J | 0.878 | 0.576 |
| x方向误差幅值/mm | 5.771 | 3.487 |
| z方向误差幅值/mm | 3.28 | 2.238 |
表选项
优化过程中，平均适配度的迭代变化如图 6所示。优化前后，物料装填机器人末端的误差变化如图 7所示。
 |
| 图 6 平均适配度变化 Fig. 6 Variation of average fitness |
| 图选项 |
 |
| 图 7 物料装填机器人末端误差变化 Fig. 7 Variation of material loading robot end-effector error |
| 图选项 |
由表 4和图 7可知，优化后，关节残余弹性势能减小约34.4%，末端x方向误差幅值减小约39.6%，z方向误差幅值减小约31.8%。仿真结果验证了本文得到的优化轨迹在机器人振动抑制方面的有效性。
5 结论 1) 本文提出一种“RPP-RR”构型的新型圆柱物料装填机器人，该机器人应用2个平动关节实现物料在竖直平面内的搬运，具有结构紧凑、适于狭小空间的特点。
2) 以新设计的物料装填机器人为研究对象，对机器人末端振动的抑制方法进行研究。以减小关节残余弹性势能为目标，通过使用最大最小蚁群算法优化插值运动时间完成关节轨迹优化，从而抑制残余振动。
3) 建立封闭形式的物料装填机器人动力学模型，并推导得到物料装填机器人沿特定轨迹运动时的末端动态响应方程。在此基础上，对优化结果进行仿真验证。仿真结果表明，优化后，末端误差在x方向和z方向分别降低了39.6%和31.8%，机器人末端振动振幅降低约37.6%。
参考文献
| [1] | KANG B S, PARK G J, ARORA J. A review of optimization of structure subjected to transient loads[J]. Structural and Multidisciplinary Optimization, 2006, 31(2): 81-95. DOI:10.1007/s00158-005-0575-4 |
| [2] | BRIDGES M M, DAWSON D M, ABDALLAH C T. Control of rigid-link, flexible-joint robots:A survey of backstepping approaches[J]. Journal of Robotic Systems, 1995, 12(3): 199-216. DOI:10.1002/rob.4620120305 |
| [3] | SINGHOSE W. Command shaping for flexible systems:A review of the first 50 years[J]. International Journal of Precision Engineering and Manufacturing, 2009, 10(4): 153-168. DOI:10.1007/s12541-009-0084-2 |
| [4] | 易江, 朱秋国, 吴俊, 等. 基于最优控制的仿人机器人行走振动抑制[J]. 机器人, 2018, 40(2): 129-135. YI J, ZHU Q G, WU J, et al. Walking vibration suppression for humanoid robot based on optimal control[J]. Robot, 2018, 40(2): 129-135. (in Chinese) |
| [5] | ZHANG L J, CHEN L. Adaptive time varying sliding mode control and vibration suppression for space manipulator with fully flexible arms and uncertain parameters[J]. Revista de la Facultad de Ingenieria, 2016, 31(9): 232-245. |
| [6] | 陈志勇, 张婷婷, 郭益深. 弹性基和弹性关节空间机器人的自适应鲁棒抗扰控制及振动抑制[J]. 自动化学报, 2018, 44(7): 1271-1281. CHEN Z Y, ZHANG T T, GUO Y S. Adaptive robust anti-interference control and vibration suppression for an elastic-base elastic-joint space robot[J]. Acta Automatica Sinica, 2018, 44(7): 1271-1281. (in Chinese) |
| [7] | REW K H, KIM K S.Using asymmetric S-curve profile for fast and vibrationless motin[C]//IEEE International Conference on Robotics and Automation.Piscataway, NJ: IEEE Press, 2011: 1-4. |
| [8] | FEDERICO C, CORRADO G, MARCO L. Minimum-jerk online planning by a mathematical programming approach[J]. Engineering Optimization, 2014, 46(6): 763-783. DOI:10.1080/0305215X.2013.806916 |
| [9] | WISNU A, KAZUHIKO T.An implementation of cubic spline optimization with floating via points and vibration suppression for trajectory planning of industrial robot arm[C]//SICE Annual Conference, 2013: 288-293. |
| [10] | GUO C Q, GAO H B, NI F L, et al.A vibration suppression method for flexible joints manipulator based on trajectory optimization[C]//Proceedings of 2016 IEEE International Conference on Mechatronics and Automation.Piscataway, NJ: IEEE Press, 2016: 338-343. |
| [11] | 吴明月, 刘延杰, 蔡鹤皋. 基于最优S型曲线轨迹的硅片传输机器人末端振动抑制方法[J]. 机器人, 2014, 36(4): 446-454. WU M Y, LIU Y J, CAI H G. Method of vibration suppression of the end-effector of wafer handling robot based on optimal S-curve profile[J]. Robot, 2014, 36(4): 446-454. (in Chinese) |
| [12] | 解则晓, 商大伟, 任凭. 基于Lame曲线的Delta并联机器人拾放操作轨迹的优化与试验验证[J]. 机械工程学报, 2015, 51(1): 52-59. XIE Z X, SHANG D W, REN P. Optimization and experimental verification of pick-and-place trajectory for a Delta parallel robot based on Lame curves[J]. Journal of Mechanical Engineering, 2015, 51(1): 52-59. (in Chinese) |
| [13] | 孙增圻, 邓志东, 张再兴. 智能控制理论与技术[M]. 2版. 北京: 清华大学出版社, 2011: 387-398. SUN Z Q, DENG Z D, ZHANG Z X. Intelligent control theory and technology[M]. 2nd ed. Beijing: Tsinghua University Press, 2011: 387-398. (in Chinese) |
| [14] | 鲍文博, 白泉, 陆海燕. 振动力学基础与MATLAB应用[J]. 北京:清华大学出版社, 2015, 157-170. BAO W B, BAI Q, LU H Y. Vibration mechanics foundation and application of MATLAB[J]. Beijing:Tsinghua University Press, 2015, 157-170. (in Chinese) |
| [15] | 熊有伦. 机器人技术基础[M]. 武汉: 华中科技大学出版社, 1996: 72-92. XIONG Y L. Robot technology foundation[M]. Wuhan: Huazhong University of Science and Technology Press, 1996: 72-92. (in Chinese) |
