ABB-IRB2600机器人动力学参数辨识方法研究
Research on Dynamic Parameter Identification Method for ABB-IRB2600 Robot
摘要: 首先通过Lagrange法构建机器人的动力学模型,经过线性变换和动力学参数重组,获取最小惯性参数集。其次,选择傅里叶级数作为激励轨迹,结合MATLAB优化工具箱对傅里叶级数系数进行优化,从而获得最佳的激励轨迹。最后将最优的激励轨迹导入到ADAMS中进行仿真,得出相关数据,并根据最小二乘法计算出待辨识动力学参数估计值,并与理论值进行对比,计算绝对误差。结果表明:仿真实验能够获得较为理想的效果,验证参数辨识方法的正确性。
Abstract: First of all, the dynamic model of the robot is built according to the Lagrange method, and the min-imum set of unbiased parameters is obtained after linear transformation and dynamic recombina-tion of parameters. Secondly, the Fourier series is chosen as the exciting trajectory, and the Fourier series coefficient is optimized in conjunction with the MATLAB optimization tools to achieve the best exciting trajectory. Finally, the best gripping trajectory is imported into ADAMS for modeling, the relevant data is obtained, and the calculated value of the known kinetic parameter is calculated by the method of least squared, and the total error is calculated by comparison with the theoretical value. The results show that by correcting the parameter identification method and verifying the parameter identification method, a more ideal effect can be achieved.
文章引用:周勇宇, 黄广荣. ABB-IRB2600机器人动力学参数辨识方法研究[J]. 建模与仿真, 2023, 12(2): 1428-1440. https://doi.org/10.12677/MOS.2023.122133

1. 引言

随着工业技术的进步,工业上要求机器人的运行轨迹的精度越来越高,执行速度也越来越快 [1]。实现机器人高精度高速度的运行轨迹的重要一个步骤就是建立精确的动力学模型 [2]。机器人的几何尺寸可以通过直接测量获得,但是对于一些动力学参数,无法通过测量直接获得 [3]。因此,对工业机器人动力学参数进行辨识是必不可少的一步。

测量机器人动力学参数方法有三种,分别是解体测量法,CAD法和整体辨识法 [4]。由于需要考虑机器人实际运行中的装配关系和材料属性,解体测量法和CAD法往往导致测量结果的不精确,而整体辨识法可以考虑机器人的实际运行情况 [4] [5],具有更高的辨识精度。因此本文采用整体辨识法进行动力学参数辨识。

机器人动力学参数辨识的步骤一般包括推导动力学模型及线性化处理、激励轨迹的选取和优化、数据采集与处理、参数估计和误差对比 [6]。国内外学者对机器人动力学参数辨识已经做了一些系统性的研究。Gautier [7] [8] 为了获得精确的动力学参数,采用卡尔曼滤波处理数据并根据最小二乘算法对机器人的动力学参数进行估计。Atkeson [9] 使用五次多项式作为机器人的激励轨迹并采用最小二乘辨识算法进行参数估计,获得了不错的辨识结果。Swevers [10] 第一次提出使用周期性傅里叶级数作为机器人的激励轨迹,取代了多项式形式的激励轨迹,并迅速被广泛应用。丁亚东等 [11] 将参数辨识分成多个步骤,多次辨识得到动力学参数,虽然可以准确地辨识出动力学参数,但是辨识过程过于复杂。刘勇 [12] 根据智能算法,利用多目标粒子群优化算法对机器人的关节轨迹进行优化,并取得了良好的效果。吴文祥 [13] 等采用最大似然估计方法对机器人的动力学参数进行估计,也取得了较好的精度效果。

本文根据Lagrange方程建立ABB-IRB2600型机器人前三轴的动力学方程,用傅里叶级数作为激励轨迹并对激励轨迹系数进行优化 [14] [15]。将激励轨迹导入到ADAMS动力学仿真平台进行仿真实验并得到了相关数据,最后根据最小二乘法计算出机器人的动力学参数。

2. 动力学模型的建立

Newton-Euler法和Lagrange法是推导机器人的动力学方程两种常用方法。其中Newton-Euler法是基于各杆件的内力关系通过正向递推和逆向递推建立机器人动力学模型 [16] [17]。Lagrange法是从能量的角度出发,计算机器人关节运动的动势能,从而推出动力学模型 [16]。

2.1. 动力学模型推导

ABB-IRB2600机器人是典型的六轴关节型机器人,一般用于物料搬运、机器上下料、点焊、弧焊、切割、装配、测试、检测、涂胶、研磨和抛光等方面。如图1所示是ABB-IRB2600机器人。

Figure 1. ABB-IRB2600 Robot

图1. ABB-IRB2600机器人

本文用Lagrange方程建立机器人前三关节的动力学的模型,主要对前三关节进行参数辨识,建立动力学方程如下:

τ = M ( q ) q ¨ + C ( q , q ˙ ) q ˙ + G ( q ) (1)

其中, τ R 3 × 1 表示关节力矩, q R 3 × 1 表示关节角度, q ˙ R 3 × 1 表示关节角速度, q ¨ R 3 × 1 表示关节角加速度, M ( q ) R 3 × 3 表示惯性矩阵, C ( q , q ˙ ) R 3 × 1 表示科里奥利和离心力, G ( q ) R 3 × 1 表示重力项。

由于机器人关节摩擦力是一个比较复杂的非线性模型,本文采用库伦摩擦+粘性摩擦的形式来对机器人关节摩擦进行建模 [18],形式如下所示:

τ f = f v q ˙ + f c sgn ( q ˙ ) (2)

其中, τ f 是关节摩擦力矩, f v 是关节粘性摩擦系数, f c 是关节库伦摩擦系数,sgn是符号函数。

因此,机器人动力学方程可以表示为:

τ = M ( q ) q ¨ + C ( q , q ˙ ) q ˙ + G ( q ) + τ f (3)

根据文献 [19],机器人的动能与势能为机器人惯性参数的线性形式,故 M ( q ) C ( q , q ˙ ) G ( q ) τ f 满足如下线性关系:

τ = M ( q ) q ¨ + C ( q , q ˙ ) q ˙ + G ( q ) + τ f = W ( q , q ˙ , q ¨ ) φ (4)

其中, W ( q , q ˙ , q ¨ ) 表示机器人观察矩阵, φ 表示机器人待辨识的动力学参数,且机器人单个臂杆i的动力学参数可表示如下形式:

φ i = [ X X i , Y Y i , Z Z i , X Y i , X Z i , Y Z i , M X i , M Y i , M Z i , M i , f c i , f v i ] (5)

其中, X X i , Y Y i , Z Z i , X Y i , X Z i , Y Z i 表示连杆i的惯性张量, M X i , M Y i , M Z i 表示连杆i的一阶质量矩, M i 是连杆i的质心。

2.2. 动力学参数重组

由于观察矩阵 W ( q , q ˙ , q ¨ ) 是一个非满秩矩阵,因此必须对其进行处理获得一组最小惯性参数,这组最小惯性参数一方面要保证数量最小,另一方面要保证每个参数对机器人的动力学特性有影响。根据文献 [19] 的方法,通过机器人动力学参数重组可以计算出这组最小惯性参数,具体方法可以参考文献 [19]。经过重组后共有21个待辨识的动力学参数,其动力学模型可以表示为如下形式:

τ = W b ( q , q ˙ , q ¨ ) φ b (6)

其中, W b ( q , q ˙ , q ¨ ) 为观测矩阵, φ b 为机器人动力学最小惯性参数也是待辨识的动力学参数。

式中: φ b = [ φ b 1 φ b 2 φ b 21 ] T 且, φ b 1 = I 1 z z + I 2 y y + I 3 y y + a 3 2 m 3 + a 2 2 ( m 2 + m 3 ) φ b 2 = I 2 x x I 2 y y a 3 2 m 3 φ b 3 = I 2 x y φ b 4 = I 2 x z + a 3 m 3 r 3 z φ b 5 = I 2 y z φ b 6 = I 2 z z + a 3 2 m 3 φ b 7 = m 2 r 2 x + a 3 m 3 φ b 8 = m 2 r 2 y φ b 9 = I 3 x x I 3 y y φ b 10 = I 3 x y φ b 11 = I 3 x z φ b 12 = I 3 y z φ b 13 = I 3 z z φ b 14 = m 3 r 3 x φ b 15 = m 3 r 3 y φ b 16 = f c 1 φ b 17 = f v 1 φ b 18 = f c 2 φ b 19 = f v 2 φ b 20 = f c 3 φ b 21 = f v 3

3. 激励轨迹优化设计

激励轨迹的设计是机器人动力学参数辨识的重要步骤,对机器人动力学参数辨识的精度尤其重要。激励轨迹的设计可分为激励轨迹的选择和激励轨迹系数的优化两个部分。本文选用傅里叶级数作为机器人的激励轨迹。

3.1. 傅里叶级数激励轨迹

傅里叶级数激励轨迹模型如式(7)所示:

q i ( t ) = l = 1 N [ a l , i ω f l sin ( w f l t ) b l , i w f l cos ( w f l t ) ] + q i 0 q ˙ i ( t ) = l = 1 N [ a l , i cos ( w f l t ) + b l , i sin ( w f l t ) ] q ¨ i ( t ) = w f l = 1 N [ b l , i l cos ( w f l t ) a l , i l sin ( w f l t ) ] (7)

式中:i表示机器人的关节数;

l表示采用的傅里叶级数的阶数;

N表示傅里叶级数的谐波项数。

选取N = 5,因此机器人每个关节的傅里叶级数有 2 N + 1 = 11 个待优化的参数, ω f 表示基频, f = 0.1 Hz 表示频率, ω f = 2 π f = 0.2 π a l , i b l , i 表示傅里叶级数系数, q i 0 表示偏移量, a l , i b l , i q i 0 表示待优化的激励轨迹系数。

3.2. 激励轨迹系数优化

关于机器人激励轨迹的优化有多种优化标准,其中使用最多的优化标准是条件数法。因此本文采用条件数法,在优化过程中,条件数越小,机器人在各关节约束范围内的动力学特性就发挥的更充分,自身误差带来的干扰就越小,参数辨识结果就更精确。

3.2.1. 目标函数

定义条件数作为目标函数:

c o n d ( W ) = σ max ( W ) σ min ( W ) (8)

式中: σ max ( W ) 表示观测矩阵W奇异值最大值;

σ min ( W ) 表示观测矩阵W奇异值最小值;

c o n d ( W ) 的表达式为 f ( a l , i , b l , i , q i 0 )

矩阵W是每个采样点对应的所有观测矩阵的集合阵:其矩阵大小为 ( N × 3 ) × N b = 30000 × 21 ,其中 N = 10000 为采样点数,3为机器人前三关节个数, N b = 21 为待辨识动力学参数集维数。

3.2.2. 约束条件

定义各关节的约束条件为:

| q i ( t ) | q i , max | q ˙ i ( t ) | q ˙ i , max | q ¨ i ( t ) | q ¨ i , max q i ( t 0 ) = q i ( t f ) = 0 , i , f q ˙ i ( t 0 ) = q ˙ i ( t f ) = 0 q ¨ i ( t 0 ) = q ¨ i ( t f ) = 0 (9)

表1表示机器人各关节具体约束大小:

Table 1. Robot joint limits

表1. 机器人各关节约束范围

根据傅里叶级数激励轨迹的数学模型,将式(7)代入到式(9)并化简后得优化参数 a l , i , b l , i , q i 0 的非线性不等式约束关系 [20] :

| q i ( t ) | = | l = 1 N [ a l , i ω f l sin ( w f l t ) b l , i w f l cos ( w f l t ) ] + q i 0 | l = 1 N 1 w f l a l , i 2 + b l , i 2 + | q i 0 | q i , max | q ˙ i ( t ) | = | l = 1 N [ a l , i cos ( w f l t ) + b l , i sin ( w f l t ) ] | l = 1 N a l , i 2 + b l , i 2 q ˙ i , max | q ¨ i ( t ) | = | w f l = 1 N [ b l , i l cos ( w f l t ) a l , i l sin ( w f l t ) ] | w f l = 1 N l a l , i 2 + b l , i 2 q ¨ i , max (10)

除此之外,设置关节角度、关节角速度和关节角加速度初始值为0, a l , i , b l , i , q i 0 的非线性等式约束条件为:

q i ( t 0 ) = q i ( t f ) = l = 1 N b l , i w f l q i 0 = 0 q ˙ i ( t 0 ) = q ˙ i ( t f ) = l = 1 N a l , i = 0 q ¨ i ( t 0 ) = q ¨ i ( t f ) = l = 1 N w f l b l , i = 0 (11)

3.2.3. Fmincon函数优化

本文采用MATLAB优化工具箱中的Fmincon函数对机器人的激励轨迹系数进行优化。下式表示Fmincon函数的详细说明:

min f ( x ) = { c ( x ) 0 c e q ( x ) = 0 A x b A e q x = b e q l b x u b (12)

式中: c ( x ) 0 代表非线性不等式约束;

c e q ( x ) = 0 代表非线性等式约束;

A x b 代表线性不等式约束;

A e q x = b e q 为线性等式约束;

l b x u b 为待优化参数的上下限。

Fmincon函数的基本形式如下:

[ x , f v a l ] = f m i n c o n ( ' f u n ' , x 0 , A , b , A e q , b e q , l b , u b , ' n o n l c o n ' , o p t i o n s ) (13)

其中,fun为目标函数,后面的nonlcon为非线性约束(包括等式和不等式);x0表示决策变量的初始值,可以随机取一组符合约束条件的数据值;A,b,Aeq,beq分别表示线性不等式约束和等式约束。关于Fmincon函数的详情可见MATLAB帮助文档。

将观测矩阵的条件数作为目标函数,根据各关节得到约束条件,用Fmincon函数对机器人的激励轨迹系数进行优化,优化后的最小条件数大小为cond(W) = 51.5,最优的激励轨迹系数如表2所示:

Table 2. Optimization results of excitation trajectories of each joint

表2. 各关节激励轨迹优化结果

机器人最优激励轨迹曲线如图2~4所示,由图中可知优化后的各关节激励轨迹都在约束的范围之内,并且优化后的激励轨迹所对应的观测矩阵条件数最小,可以充分地激发机器人的动力学特性,降低观测矩阵自身带来的计算误差,提高动力学参数的辨识精度。

Figure 2. Joint angle trajectory

图2. 关节角度轨迹

Figure 3. Joint angular velocity trajectories

图3. 关节角速度轨迹

Figure 4. Joint angular acceleration trajectories

图4. 关节角加速度轨迹

4. 采样数据处理与参数估计

4.1. 采样数据处理

将最优激励轨迹输入到ADAMS仿真软件中,使机器人在ADAMS中运行最优激励轨迹,并测量各关节的角度和力矩数据,多次重复运行激励轨迹,并求取多次测量数据的均值以便于后续的计算。将测量得到的关节力矩和关节角度进行滤波处理和关节角度数据的差分运算,最终得到仿真过程中的力矩,角度,角速度及角加速度数据。

4.2. 最小二乘估计算法

当采集完数据后需要进行机器人的参数估计,现已有多种算法可进行参数估计,本文用最小二乘法进行动力学参数的估计。

对于如下模型,如果定义:

h ( k ) = [ y ( k 1 ) , y ( k 2 ) , , y ( k n ) , u ( k 1 ) , u ( k 2 ) , , u ( k n ) ] (14)

θ = [ a 1 , a 2 , , a n , b 1 , b 2 , , b n ] T (15)

z ( k ) = h ( k ) θ + v ( k ) (16)

在式中, θ 代表待估计的参数。

若令 k = 1 , 2 , , m ,则有:

Z m = [ z ( 1 ) z ( 2 ) z ( m ) ] (17)

H m = [ h ( 1 ) h ( 2 ) h ( m ) ] = [ y ( 0 ) y ( 1 n ) u ( 0 ) u ( 1 n ) y ( 1 ) y ( 2 n ) u ( 1 ) u ( 2 n ) y ( m 1 ) y ( m n ) u ( m 1 ) u ( m n ) ] (18)

θ = [ a 1 a n b 1 b n ] T (19)

V m = [ v ( 1 ) v ( 2 ) v ( m ) ] T (20)

故:

Z m = H m θ + V m (21)

最小二乘算法的原理就是找到一个 θ 的估计值 θ ^ ,使得各次测量的 Z i ( i = 1 , , m ) 与由估计 θ ^ 所确定的 Z ^ i = H i θ ^ 之差的平方和最小,即:

J ( θ ^ ) = ( Z m H m θ ^ ) T ( Z m H m θ ^ ) = min (22)

J θ | θ = θ ^ = 2 H m T ( Z m H m θ ^ ) = 0 (23)

H m T H m θ ^ = H m T Z m (24)

如果 H m 的行数大于等于列数,即 m 2 n ,并且 H m T H m 满秩,即 rank ( H m T H m ) = 2 n ,则 ( H m T H m ) 1 存在。因此 θ 的最小二乘估计为:

θ ^ = ( H m T H m ) 1 H m T Z m (25)

5. 仿真实验与结果分析

本文通过ADAMS动力学分析软件对机器人进行仿真实验,并借助MATLAB对数据进行处理,以验证上述理论的正确性。仿真过程中用到的机器人模型与实际机器人实体一致,实验步骤如下:

步骤一:动力学参数辨识实验仿真。将优化后的激励轨迹导入到ADAMS软件中,如图5所示,并定义机器人各杆件理论惯性参数等属性,使机器人运行最优激励轨迹,测量关节力矩和角度数据。

步骤二:数据处理。对测量得到的关节力矩和关节角度进行滤波处理和角度的差分运算,计算出关节角速度和角加速度数据。

步骤三:动力学参数辨识。将测量和计算得到数据代入机器人动力学方程,并根据最小二乘法计算机器人动力学参数。

步骤四:验证与分析。将机器人理论惯性参数值与本文辨识算法计算的惯性参数值作比较并计算绝对误差,以验证参数辨识的准确性。

Figure 5. The motion path under the robot excitation trajectory

图5. 机器人激励轨迹下的运动路径

Figure 6. Joint 1 after filtering measures the angle

图6. 滤波后关节1的测量角度

Figure 7. Joint 2 after filtering measures the angle

图7. 滤波后关节2的测量角度

Figure 8. Joint 3 after filtering measures the angle

图8. 滤波后关节3的测量角度

Figure 9. Joint 1 Torque

图9. 滤波后关节1力矩

Figure 10. Joint 2 Torque

图10. 滤波后关节2力矩

Figure 11. Joint 3 Torque

图11. 滤波后关节3力矩

图6~8分别是滤波后得到的各关节角度数据,图9~11分别是滤波后得到的各关节力矩数据。最后将处理后的数据代入到最小二乘辨识算法中进行参数估计,最终得到机器人的动力学参数值如表3所示:

Table 3. Dynamic inertial parameter identification results

表3. 动力学惯性参数辨识结果

表3的第4列所示是机器人的理论惯性参数和辨识算法计算得到的惯性参数的绝对误差,其绝对误差较小,由此可以证明本文的辨识算法的正确性。

6. 总结

1) 由Lagrange法推导工业机器人前三关节的动力学方程,处理方程得到待辨识的动力学参数。

2) 选用傅里叶级数作为机器人激励轨迹,并借助MATLAB的优化工具箱对机器人的激励轨迹进行优化,最终得到一组最优激励轨迹。

3) 将优化后的最优激励轨迹导入到动力学仿真软件中,使机器人运行最优激励轨迹,测量关节角度和关节力矩,并进行数据处理,将处理得到的数据代入机器人动力学方程,并根据最小二乘法计算机器人动力学参数。

4) 将计算值与理论值进行对比。结果表明:理论值与辨识值绝对误差较小,证明本文辨识方法的正确性和有效性,为机器人实际动力学参数辨识提供了参考依据。

参考文献

[1] 泮求亮. 基于免疫算法和神经网络的机械臂动力学参数辨识与验证[D]: [硕士学位论文]. 杭州: 浙江大学, 2021.
[2] 吉阳珍. 6R工业机器人运动学及动力学关键问题研究[D]: [硕士学位论文]. 成都: 四川大学, 2021.
[3] 中国人工智能学会. 中国人工智能学会第9届全国学术年会论文集[M]. 北京: 北京邮电大学出版社, 2001.
[4] 严浩, 白瑞林, 吉峰. 一种改进的SCARA机器人动力学参数辨识方法[J]. 中国机械工程, 2017, 28(22): 2707-2713.
[5] Benimeli, F., Mata, V. and Valero, F. (2006) A Comparison between Direct and Indirect Dynamic Parameter Identification Methods in Industrial Ro-bots. Robotica, 24, 579-590.
https://doi.org/10.1017/S0263574706002645
[6] 杨永雄. 固晶机摆臂的动力学建模与运动控制系统设计[D]: [硕士学位论文]. 广州: 广东工业大学, 2022.
[7] Gautier, M. and Poignet, P. (2001) Extended Kal-man Filtering and Weighted Least Squares Dynamic Identification of Robot. Control Engineering Practice, 9, 1361-1372.
https://doi.org/10.1016/S0967-0661(01)00105-8
[8] Gautier, M. and Khalil, W. (1988) A Direct Determination of Min-imum Inertial Parameters of Robots. IEEE International Conference on Robotics and Automation, Philadelphia, 24-29 April 1988, 1682-1687.
[9] Atkeson, C.G., An, C.H. and Hollerbach, J.M. (1986) Estimation of Inertial Parameters of Manipulator Loads and Links. The International of Robotic Research, 5, 101-119.
https://doi.org/10.1177/027836498600500306
[10] Swevers, J., Ganseman, C., Tukel, D.B., et al. (1997) Optimal Robot Excitation and Identification. IEEE Transactions on Robotics and Automation, 13, 730-740.
https://doi.org/10.1109/70.631234
[11] 丁亚东, 陈柏, 吴洪涛. 一种工业机器人动力学参数的辨识方法[J]. 华南理工大学学报, 2015, 43(3): 49-56.
[12] 刘勇, 庆轩, 陈钢等. 基于多目标粒子群优化算法的自由漂浮空间机器人负载最大化轨迹优化[J]. 机器人, 2014, 36(4): 402-410.
[13] 吴文祥. 多自由度串联机器人关节摩擦分析与低速高精度运动控制[D]: [博士学位论文]. 杭州: 浙江大学, 2013.
[14] 冯利民, 俞经虎, 王延玉, 刘佳怡. 基于IRLS算法的机器人动力学参数辨识[J]. 现代制造工程, 2022, 23(4): 37-46.
[15] 方骏玮, 王斌锐, 谢胜龙, 任海军. 工业机器人动力学参数分步辨识与零力控制示教[J]. 机械设计与制造, 2021, 5(1): 240-244.
[16] 田莉莉. 串联机器人动力学特性及结构优化设计研究[D]: [硕士学位论文]. 济南: 山东大学, 2020.
[17] 赵京, 李玉龙. 冗余度机器人动力学容错性能[J]. 北京工业大学学报, 2011, 37(10): 1441. 1445.
[18] Hur, S.M., Kim, S.K., Oh, Y., et al. (2012) Joint Torque Ervo of a High Friction Robot Manipulator Based on Time-Delay Control with Feed-Forward Friction Compensation. 2012 IEEE RO-MAN: The 21st IEEE International Symposium on Robot and Human Interactive Communication, Paris, 9-13 September 2012, 37-42.
https://doi.org/10.1109/ROMAN.2012.6343728
[19] 霍伟. 机器人动力学与控制[M]. 北京: 高等教育出版社, 2005: 89-121.
[20] Jin, J. and Gans, N. (2015) Parameter Identification for Industrial Robots with a Fast and Robust Trajectory De-sign. Robotics and Computer-Integrated Manufacturing, 31, 21-29.
https://doi.org/10.1016/j.rcim.2014.06.004