, 陈志雄2, 刘春3, 薛俊41. 清华大学 机械工程系, 北京 100084;
2. 电子科技大学 机械与电气工程学院, 成都 611173;
3. 成都飞机工业(集团)有限责任公司, 成都 610092;
4. 中国航空制造技术研究院, 北京 100025
收稿日期：2020-10-29
基金项目：国家重点研发计划（2017YFB1301702）
作者简介：关立文(1972-), 男, 研究员。E-mail: guanlw@tsinghua.edu.cn
摘要：机器人自动钻铆技术突破了飞机装配质量和灵活性瓶颈，在装配效率和成本上有巨大优势，已成为飞机数字化、柔性化装配发展的必然选择。目前机器人钻铆系统大多以6自由度工业串联机器人为平台，此类机器人具有结构弱刚性的固有缺点，在钻铆作业过程中承受较大的压紧力和制孔力等载荷时，机器人末端和执行器会发生一定的变形甚至颤振，严重影响钻铆精度和质量。该文以钻铆机器人为对象建立机器人静刚度模型，设计关节刚度辨识实验获取了机器人关节刚度值，进一步结合刚度性能评价指标分析了机器人工作空间范围内刚度性能分布特征，根据实际作业工况，采用粒子群算法对机器人进行位姿优化，从而增强了机器人刚度性能，对于保证系统作业稳定性和作业质量具有重要意义。
关键词：工业机器人静刚度关节刚度辨识刚度性能粒子群算法位姿优化
Static stiffness modeling for optimizing drilling and riveting robots
GUAN Liwen1
, CHEN Zhixiong2, LIU Chun3, XUE Jun41. School of Mechanical Engineering, Tsinghua University, Beijing 100084, China;
2. School of Mechanical and Electrical Engineering, University of Electronic Science and Technology of China, Chengdu 611173, China;
3. Chengdu Aircraft Industrial (Group) Co., Ltd., Chengdu 610092, China;
4. AVIC Manufacturing Technology Institute, Beijing 100025, China
Abstract: Automatic robotic drilling and riveting greatly improve aircraft assembly quality and flexibility while reducing costs, so automatic drilling and riveting are being widely used for aircraft assembly. However, almost all robotic drilling and riveting systems are based on six degrees of freedom serial robots. However, these robots have an inherent weakness due to their poor rigidity. Large pressing and drilling forces during drilling and riveting can lead to deformation and even flutter of the robot end actuator, which seriously affect the drilling and riveting accuracy and quality. A static stiffness model of a drilling and riveting robot is developed with joint stiffness measurements to predict robot joint stiffnesses. A stiffness evaluation index is then used to characterize the robotic arm stiffness in the working space. The stiffness index can be used to optimize the robot position and posture for a specific operation using a particle swarm algorithm to improve the stiffness. This improves the system stability and quality.
Key words: industrial robotstatic stiffnessjoint stiffness identificationstiffness performanceparticle swarm algorithmpose optimization
工业机器人以其低成本、高效率、高柔性的优势广泛应用于机械制造、航空航天、海洋探测等多个领域。随着现代飞机对性能和寿命，以及航空制造过程对飞机制造的高质量、高效率、低成本的要求越来越高，机器人自动钻铆技术以其独有的优势成为飞机数字化、柔性化装配的必然选择[1]。然而，自动钻铆系统普遍采用的6自由度串联机器人存在机构弱刚性问题，很大程度影响着机器人钻铆加工精度和加工质量[2]。为此，通过分析机器人的静刚度和优化机器人系统作业时的位姿来增强机器人刚度性能，对于改善系统加工质量、保证系统作业稳定性具有重要意义。
在机器人刚度建模方面，Chen等提出守恒合同变换法(conservation contract transformation，CCT)，建立了6自由度旋转关节机器人6×6的Descartes刚度矩阵[3]。Abele等完成了机器人铣削系统工业机器人的刚度建模及参数辨识，分析了机器人末端滑移变形规律，并对其加工性能作出评价[4]。曲巍葳等分析了KUKA KR360的奇异性和灵巧性，通过实验法获得机器人关节刚度，并利用遗传算法对机器人刚度性能进行了优化[5-6]。布音研究了机器人制孔系统中机器人运动学性能、刚度建模以及刚度辨识，获取机器人关节刚度参数，分析了机器人作业典型位姿下刚度特性对制孔质量的影响规律[7]。
机器人刚度是衡量机器人负载能力的重要指标，也是影响机器人钻铆加工质量的关键参数。本文以KUKA KR500L340型钻铆机器人为研究对象，建立了静刚度模型；然后利用机器人系统、API激光跟踪仪和Kistler测力仪搭建实验平台，完成了机器人关节刚度辨识实验，获取了机器人关节刚度值；在此基础上，结合刚度性能评价指标，采用Monte Carlo法仿真得到机器人工作空间内的刚度性能云图，为在机器人作业时选择较优刚度性能区域给出参考；进一步结合机器人钻铆系统作业时的具体约束条件，采用粒子群算法优化了机器人的位姿，通过对优化前后机器人末端承受相同载荷时的变形量进行对比，证明了通过优化机器人位姿可实现增强其刚度性能的目的，对改善机器人系统作业质量和稳定性具有重要意义。
1 机器人静刚度建模KUKA KR500L340型钻铆机器人为6自由度旋转关节机器人。本文首先采用改进Denavit-Hartenberg(DH)参数法[8-9]在机器人各关节处建立坐标系，进而得到其运动学模型并列出DH参数。该机器人具体结构如图 1所示，DH参数见表 1，运动学模型如图 2所示。
 |
| 图 1 (网络版彩图)KUKA KR500L340型机器人 |
| 图选项 |
表 1 KUKA KR500L340机器人DH参数表
| 关节i | ai-1/m | αi-1/(°) | di/m | θi/(°) |
| 1 | 0 | 0 | 1.045 | 0 |
| 2 | 0.5 | 90 | 0 | 90 |
| 3 | 1.3 | 0 | 0 | 0 |
| 4 | 0.055 | 90 | 1.525 | 0 |
| 5 | 0 | -90 | 0 | 0 |
| 6 | 0 | 90 | 0.29 | 180 |
表选项
 |
| 图 2 机器人运动学模型 |
| 图选项 |
根据机器人运动学模型和DH参数推导出机器人运动学正解，采用微分变换法构造机器人运动学Jacobi矩阵(机器人关节速度与末端操作速度的映射关系式[10])。由于重载机器人连杆是由铸铁和高强度钢等材料构成，关节处由电机、减速器等传动机构组成，在承受载荷时，机器人关节处的变形远大于连杆变形，因此机器人静刚度建模时假定连杆为刚性[11]，且机器人末端变形全部由关节变形引起。根据Jacobi矩阵和Hooker定律有关公式，建立机器人静刚度模型为
| $\boldsymbol{F}=\boldsymbol{J}^{-\mathrm{T}} \boldsymbol{K}_{\theta} \boldsymbol{J}^{-1} \boldsymbol{X} . $ | (1) |
机器人运动学Jacobi矩阵随位姿变化而变化，机器人末端变形不仅与外力有关，还与其位姿有关。
2 机器人关节刚度辨识获取机器人关节刚度值的方法有2种：1) 通过关节处的结构、材料和尺寸等参数进行计算[12]；2) 设计关节刚度辨识实验[13]，通过处理实验数据得到。由于第1种方法过于理想化，关节处具体结构等参数很难完全掌握，因此本文选择通过关节刚度辨识实验的方法获取机器人关节刚度值。
2.1 关节刚度辨识原理将静刚度模型(式(1))变形处理得到
| $\boldsymbol{X}=\boldsymbol{J} \boldsymbol{K}_{\theta}^{-1} \boldsymbol{J}^{\mathrm{T}} \boldsymbol{F} . $ | (2) |
| $\left[\begin{array}{c}\boldsymbol{U}_{1} \\\boldsymbol{U}_{2} \\\boldsymbol{U}_{3} \\\boldsymbol{U}_{\mathrm{R} 1} \\\boldsymbol{U}_{\mathrm{R} 2} \\\boldsymbol{U}_{\mathrm{R} 3}\end{array}\right]=\left[\begin{array}{ccc}\boldsymbol{J}_{11} \sum\limits_{i=1}^{6} \boldsymbol{J}_{i 1} \bf{F}_{i} & \cdots & \boldsymbol{J}_{16} \sum\limits_{i=1}^{6} \boldsymbol{J}_{i 6} \bf{F}_{i} \\\vdots & & \vdots \\\boldsymbol{J}_{61} \sum\limits_{i=1}^{6} \boldsymbol{J}_{i 1} \boldsymbol{F}_{i} & \cdots & \boldsymbol{J}_{66} \sum\limits_{i=1}^{6} \boldsymbol{J}_{i 6} \boldsymbol{F}_{i}\end{array}\right]\left[\begin{array}{c}\bf{K}_{\theta_{1}}^{-1} \\\vdots \\K_{\theta_{6}}^{-1}\end{array}\right] . $ | (3) |
进一步处理得到
| $\boldsymbol{X}=\boldsymbol{A} \boldsymbol{c} . $ | (4) |
通过关节刚度辨识实验测量机器人在不同位姿下的力与变形数据，即可求解得到c，进而得到关节刚度矩阵Kθ。
2.2 关节刚度辨识实验机器人关节刚度辨识实验分3步：
步骤1?? 标定机器人基坐标系与激光跟踪仪坐标系的转换关系；
步骤2?? 标定测力仪坐标系与激光跟踪仪坐标系的转换关系；
步骤3?? 测量机器人末端的力与变形。
图 3所示为机器人关节刚度辨识实验现场。保持标定完成后的机器人系统、激光跟踪仪和测力仪等设备的位置不变，将激光跟踪仪靶镜安装在末端执行器靠近机器人末端法兰处。实验过程中，首先调节机器人各关节角来调整位姿使末端执行器靠近工件表面，记录此时的给关节角与末端位姿参数；然后控制压脚机构使其贴于工件表面并施加压紧力，记录测力仪与激光跟踪仪的数据；上述步骤重复6次，最终测得6组位姿，每组位姿可以得到16个测量点的力与变形数据，由于扭转变形量UR1、UR2、UR3无法测得，因此取U1、U2、U3，总计可以得到由288个方程构成的超定方程组。
 |
| 图 3 (网络版彩图)关节刚度辨识实验现场 |
| 图选项 |
2.3 辨识结果首先求解超定方程组。采用最小二乘法可以解得c，进一步将其变换成矩阵形式求逆，得到机器人关节刚度矩阵Kθ，最终得到的KUKA KR500L340机器人的关节刚度矩阵为
| $\begin{aligned}&\boldsymbol{K}_{\theta}=\operatorname{diag}\left(7.65 \times 10^{6}, 5.76 \times 10^{6}, 2.43 \times 10^{7},\right. \\&\left.8.49 \times 10^{6}, 3.43 \times 10^{6}, 3.88 \times 10^{6}\right) \mathrm{N} \cdot \mathrm{m} / \mathrm{rad} .\end{aligned} $ | (5) |
 |
| 图 4 机器人末端理论预测变形量与实际变形量对比 |
| 图选项 |
3 机器人工作空间刚度性能分析结合所辨识出的关节刚度矩阵和机器人处于不同位姿时的Jacobi矩阵，根据式(6)可以得到机器人在不同位姿下的末端操作刚度矩阵Kx，
| $\boldsymbol{K}_{x}=\boldsymbol{J}^{-\mathrm{T}} \boldsymbol{K}_{\theta} \boldsymbol{J}^{-1} . $ | (6) |
| $\left[\begin{array}{l}\boldsymbol{f} \\\boldsymbol{m}\end{array}\right]=\left[\begin{array}{ll}\boldsymbol{k}_{f d} & \boldsymbol{k}_{f \hat{o}} \\\boldsymbol{k}_{m d} & \boldsymbol{k}_{m \delta}\end{array}\right]\left[\begin{array}{l}\boldsymbol{d} \\\boldsymbol{\delta}\end{array}\right] . $ | (7) |
为重点研究力与平移变形关系，令δ=0，由式(7)可以得到
| $\boldsymbol{f}=\boldsymbol{k}_{f d} \cdot \boldsymbol{d} . . $ | (8) |
| $\left|\boldsymbol{f}\right|^{2}=\boldsymbol{f}^{\mathrm{T}} \boldsymbol{f}=\boldsymbol{d}^{\mathrm{T}} \boldsymbol{k}_{f d}^{\mathrm{T}} \boldsymbol{k}_{f d} \boldsymbol{d}=1 \text {. } $ | (9) |
同理，将式(7)进行变形处理得到
| $\left[\begin{array}{l}\boldsymbol{d} \\\boldsymbol{\delta}\end{array}\right]=\left[\begin{array}{cc}\boldsymbol{c}_{f d} & \boldsymbol{c}_{m d} \\\boldsymbol{c}_{f \hat{o}} & \boldsymbol{c}_{m \delta}\end{array}\right]\left[\begin{array}{l}\boldsymbol{f} \\\boldsymbol{m}\end{array}\right] . $ | (10) |
可见，力-线位移柔度矩阵cfd也对应着一个椭球体，可以称之为力椭球，表示使机器人末端产生单位平移变形所需的外力。椭球的短半轴对应着矩阵cfd的最大特征值λmax(cfd)，表示引起机器人末端单位平移变形所需的最小外力，可以将其作为刚度性能评价指标。λmax(cfd)越大，说明使机器人末端产生单位平移变形所需的外力越大，机器人的刚度性能越好。
3.2 机器人工作空间刚度性能分析本文利用Monte Carlo法和机器人DH参数仿真分析机器人工作空间[17]，根据机器人工作空间内每一点位姿和操作刚度矩阵，得到机器人工作空间内的刚度性能云图，图 5和6所示为分别采用λmin(kfd)和λmax(cfd)为评价指标得到的该机器人的工作空间刚度性能云图和XOZ平面投影。
 |
| 图 5 (网络版彩图)以λmin(kfd)为评价指标的机器人工作空间刚度性能云图及XOZ平面投影图 |
| 图选项 |
 |
| 图 6 (网络版彩图)以σmax(cfd)为评价指标的机器人工作空间刚度性能云图及XOZ平面投影图 |
| 图选项 |
当以λmin(kfd)为评价指标时，深蓝色区域代表刚度性能较好区域，橙黄色区域代表刚度性能较差区域，而以λmax(cfd)为评价指标时，橙黄色区域代表刚度性能较好区域，深蓝色代表刚度性能较差区域。由图 5和6的分析结果可知，两种评价指标下刚度性能较好和较差区域完全一致。
4 基于刚度性能的位姿优化4.1 优化问题数学模型机器人钻铆系统作业要求刀具刀尖点与目标加工点重合，刀具进给方向与目标点法线方向一致，此时约束了机器人5个自由度，还存在绕目标点法线方向旋转自由度没有约束。此时，对于6自由度机器人，满足加工条件的位姿有无穷多个，为此以刚度性能最大化为目标对机器人作业位姿进行优化。
优化数学模型包括目标函数、优化变量和约束条件。选择力-线位移刚度矩阵的最小特征值λmin(kfd)作为刚度性能评价指标和优化目标函数；优化变量为机器人关节角θ，对应着机器人位姿；约束条件为机器人各个关节角限位和加工时位姿约束，即刀尖点与目标点重合、刀具进给方向与目标点法线方向一致。最终得到优化问题的数学模型为
| $\left\{ \begin{array}{l}\min \;\;\;\;\;\;\;\;\;\;\;\;\;\;f(\theta )\\{\rm{s}}{\rm{. t}}\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;H(\theta ) = 0,\\\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;{\theta _{i\min }} \leqslant {\theta _i} \leqslant {\theta _{{i \mathop{\rm max}\nolimits} }}\quad (i = 1,2, \cdots ,6).\end{array} \right. $ | (11) |
4.3 优化结果选取一组机器人常用作业位姿，根据优化模型和算法，利用MATLAB求解得到优化后机器人位姿。表 2为位姿优化前后机器人各关节角和刚度性能指标对比，机器人位姿状态对比见图 7。
表 2 机器人位姿优化前后关节角及刚度性能指标对比
| 位姿 | θ1/(°) | θ2/(°) | θ3/(°) | θ4/(°) | θ5/(°) | θ6/(°) | 刚度性能指标k/106 |
| 优化前 | 0.001 | -78.824 | 108.327 | 6.14 | 61.734 | 180.058 | 2.583 |
| 优化后 | 12.822 | -67.635 | 102.931 | -45.049 | 75.828 | 213.499 | 2.265 |
表选项
 |
| 图 7 (网络版彩图)机器人优化前后位姿对比 |
| 图选项 |
4.4 优化结果验证假设在优化前后两组位姿下机器人末端执行器上刀尖点处，沿进给反方向施加400~1 000 N的外力，结合两组位姿下Jacobi矩阵和关节刚度矩阵，计算得到优化前后两组位姿下机器人末端理论变形量，结果对比如图 8所示。可见，优化后的变形量明显降低。
 |
| 图 8 机器人位姿优化前后末端受力变形对比 |
| 图选项 |
5 结论结合机器人Jacobi矩阵和Hooker定律建立了某型钻铆机器人静刚度模型；设计并完成了机器人关节刚度辨识实验，获取了机器人关节刚度值，验证了所辨识出的关节刚度值的准确性；以力-线位移柔度矩阵的最小特征值和力-线位移柔度矩阵的最大特征值两种刚度性能评价指标，结合Monte Carlo法得到机器人工作空间内刚度性能云图。以机器人刚度性能为目标，采用粒子群算法优化了机器人位姿。优化前后机器人末端受力变形情况分析结果表明，机器人位姿优化可提高机器人作业刚度，对保证机器人系统作业质量和稳定性都具有重要意义。
参考文献
| [1] | 杜兆才, 姚艳彬, 王健. 机器人钻铆系统研究现状及发展趋势[J]. 航空制造技术, 2015(4): 26-31. DU Z C, YAO Y B, WANG J. Research status and development trends of robot drilling and riveting system[J]. Aeronautical Manufacturing Technology, 2015(4): 26-31. (in Chinese) |
| [2] | LI C J, QU R X, LI S Q, et al. Robotic three-dot force feedback to suppress surface contact slipping in robot drilling[J]. Applied Mechanics and Materials, 2013, 404: 650-656. DOI:10.4028/www.scientific.net/AMM.404.650 |
| [3] | CHEN S F, KAO I. Geometrical approach to the conservative congruence transformation (CCT) for robotic stiffness control [C]//Proceedings 2002 IEEE International Conference on Robotic and Automation. Washington DC, USA, 2002: 544-549. |
| [4] | ABELE E, WEIGOLD M, ROTHENBüCHER S. Modeling and identification of an industrial robot for machining applications[J]. CIRP Annals, 2007, 56(1): 387-390. DOI:10.1016/j.cirp.2007.05.090 |
| [5] | 曲巍崴, 侯鹏辉, 杨根军, 等. 机器人加工系统刚度性能优化研究[J]. 航空学报, 2013, 34(12): 2823-2832. QU W W, HOU P H, YANG G J, et al. Research on the stiffness performance for robot machining systems[J]. Acta Aeronautica et Astronautica Sinica, 2013, 34(12): 2823-2832. (in Chinese) |
| [6] | 侯鹏辉. 机器人加工系统刚度性能优化研究[D]. 杭州: 浙江大学, 2013. HOU P H. Study on the stiffness performance optimization for robot machining system [D]. Hangzhou: Zhejiang University, 2013. (in Chinese) |
| [7] | 布音. 工业机器人精密制孔系统刚度特性研究[D]. 南京: 南京航空航天大学, 2017. BU Y. Analysis of stiffness properities for robotic precise drilling system [D]. Nanjing: Nanjing University of Aeronautics and Astronautics, 2017. (in Chinese) |
| [8] | HARTENBERG R S, DENAVIT J. Kinematic synthesis of linkages[M]. New York, USA: McGraw-Hill, 1964. |
| [9] | KHALIL W, DOMBRE E. Compliant motion control[M]. Amsterdam, Netherlands: Elsevier, 2004: 377-393. |
| [10] | 蔡自兴. 机器人学[M]. 北京: 清华大学出版社, 2000. CAI Z X. Robotics[M]. Beijing: Tsinghua University Press, 2000. (in Chinese) |
| [11] | ALICI G, SHIRINZADEH B. Enhanced stiffness modeling, identification and characterization for robot manipulators[J]. IEEE Transactions on Robotics, 2005, 21(4): 554-564. DOI:10.1109/TRO.2004.842347 |
| [12] | 陈首彦. 机器人的切削加工模型及其控制方法研究[D]. 广州: 华南理工大学, 2017. CHEN S Y. Modeling and control research on robotic machining process [D]. Guangzhou: South China University of Technology, 2017. (in Chinese) |
| [13] | 陈玉山. 6R型工业机器人关节刚度辨识与实验研究[D]. 武汉: 华中科技大学, 2011. CHEN Y S. Joint stiffness identification of 6R industrial robot and experimental verification [D]. Wuhan: Huazhong University of Science and Technology, 2011. (in Chinese) |
| [14] | CAMPOLO D. Cartesian stiffness for wrist joints: Analysis on the Lie group of 3D rotations and geometric approximation for experimental evaluation[J]. Computer Methods in Biomechanics and Biomedical Engineering, 2013, 16(9): 975-986. DOI:10.1080/10255842.2011.646392 |
| [15] | 张永贵, 刘晨荣, 刘鹏. 6R工业机器人刚度分析[J]. 机械设计与制造, 2015(2): 257-260. ZHANG Y G, LIU C R, LIU P. 6R industrial robot stiffness analysis[J]. Machinery Design & Manufacture, 2015(2): 257-260. (in Chinese) |
| [16] | 汪博文. 多机械臂协同加工系统静刚度建模与优化研究[D]. 上海: 上海大学, 2018. WANG B W. The research on stiffness modeling and optimization of multiple coordinated robots system [D]. Shanghai: Shanghai University, 2018. (in Chinese) |
| [17] | 蔡蒂, 谢存禧, 张铁, 等. 基于蒙特卡洛法的喷涂机器人工作空间分析及仿真[J]. 机械设计与制造, 2009(3): 161-162. CAI D, XIE C X, ZHANG T, et al. Study on workspace analysis and simulation of 6-DOF painting robot based on Monte-Carlo method[J]. Machinery Design & Manufacture, 2009(3): 161-162. (in Chinese) |
| [18] | 方峻. 粒子群算法及其应用研究[D]. 成都: 电子科技大学, 2006. FANG J. Research on particle swarm optimization and its application [D]. Chengdu: University of Electronic Science and Technology of China, 2006. (in Chinese) |
