新闻资讯

联系方式

青岛恒昌机器人科技有限公司
地址：青岛市黄岛区临港工业园陡楼山路367号手机：13806390681
服务热线：0532-86131102
邮箱：qdhengchangkeji@163.com
www.hengchangrobot.com

新闻资讯

当前位置：网站首页 > 新闻资讯

工业机器人动力学参数辨识与自适应控制方法研究

摘要：

研究六自由度工业机器人动力学模型小惯性参数辨识和模型参数自适应PD控制方法。首先分析六自由度机器人动力学模型及其小惯性参数;研究基于位置、速度、加速度约束条件的傅里叶级数型激励轨迹优化方法;依据激励轨迹跟踪实验获取的关节角位置、速度、加速度和力矩数据,研究基于小二乘的小惯性参数估计方法。在此基础上,研究六自由度机器人模型自适应PD控制方法。后,构建了基于Codesys平台的六自由度机器人控制系统,利用SYMORO+推导回归矩阵元素,结合Hadamard不等式,利用MatlabFunction函数生成目标函数并将其代入Matlab fmincon函数计算激励轨迹参数,通过激励轨迹跟踪实验辨识出小惯性参数;通过机器人跟踪激励轨迹和验证轨迹实验,比较实测力矩与基于辨识模型估计力矩的均方根误差验证参数辨识方法有效性;通过期望轨迹跟踪实验验证了自适应PD控制算法可行性。

关键词：

工业机器人; 轨迹跟踪; 参数辨识; 自适应控制;

作者简介： 徐建明(1970—),男,江西吉安人,教授,博士,博士生导师,研究方向为鲁棒控制、运动控制、机器人控制及其应用,E-mail:xujm@zjut.edu.cn。;

收稿日期：2019-09-11

基金： 国家自然科学基金-浙江省自然科学基金联合基金两化融合项目(U1709213); 国家自然科学基金面上项目(61374103);

Research on dynamic parameter identification and adaptive control method of industrial robot

XU Jianming ZHAO Shuai

College of Information Engineering,Zhejiang University of Technology

Abstract：

In this paper, the identification of minimum inertia parameters and the adaptive PD control method of model parameters for a 6-DOF(six-degree-of-freedom) industrial robot dynamics model are studied. Firstly, the dynamic model of 6-DOF robot and its minimum inertia parameters are analyzed. The Fourier series excitation trajectory optimization method is studied based on position, velocity and acceleration constraints. The least square method of minimum inertial parameter estimation is studied according to the joint angular position, velocity, acceleration and moment data obtained from the excitation trajectory tracking experiment. On this basis, the model adaptive PD control method of 6-DOF robot is studied. Finally, a 6-DOF robot control system based on Codesys platform is constructed. The regression matrix elements are deduced by SYMORO+ and combined with Hadamard's inequality. The objective function is generated by MatlabFunction function and substituted into Matlab fmincon function to calculate the excitation trajectory parameters. The minimum inertia parameters are identified by the experiment of tracking the excitation trajectories. Through the experiments of tracking the excitation trajectories and verifying the trajectories, the parameter identification method is validated by analyzing the root mean square error between the measured moment and the estimated moment based on the identification model. The feasibility of the adaptive PD control algorithm is verified by the experiment of tracking the expected trajectory.

Keyword：

industrial robot; trajectory tracking; parameter identification; adaptive control;

Received： 2019-09-11

工业机器人广泛应用于工业生产过程中的搬运、装配、涂装等方面,并且正向着高速、高精度、高智能化的方向发展,加之产业需求的不断扩大与提升,诸如新型的激光焊接、激光切割、多机器人协作等任务对机器人的控制精度提出了更高要求[1,2]1-2]。设计基于机器人动力学模型的控制器是实现高速高精度运动控制的有效途径[3]3],但该类控制器的设计需要已知机器人的动力学参数。一般情况下,机器人的动力学参数难以直接测量,需采用特定的辨识实验获得[4]4]。Atkeson等[5]5]选择关节空间的五次多项式作为激励轨迹,采用小二乘法作为参数估计的方法。Swevers等[6]6]首次提出了采用基于傅里叶级数的激励轨迹进行辨识实验,该方法产生较大影响,在随后的辨识研究中,基本均采用傅里叶级数作为激励轨迹[7]7]。傅里叶级数是指由不同幅值和周期的三角函数构成的非线性函数,该函数具有连续的多阶导数,其二阶导数具有较大幅值,能充分激发系统的动态特性;并且通过合理设置参数,能保证函数起始和结束时刻具有相同值,便于机器人重复执行激励轨迹。但传统傅里叶级数形式的激励轨迹无法保证各关节的角速度和角加速度在轨迹的起始和停止时刻为零,导致机器人在执行激励轨迹时存在抖动。文献[8]在此基础上提出了一种基于改进的傅里叶级数的辨识方法,该方法考虑速度、加速度边界条件,保证激励轨迹起始阶段的稳定性,提高了辨识精度。目前,确定激励轨迹参数有两种应用较多的方法,一种是基于回归矩阵的条件数小化方法[9]9];另一种是基于对数行列式的优化方法[10]10]。参数估计的方法也有很多,其中应用较多的两种方法为小二乘法[11]11]和大似然估计法[8]8]。其他的方法包括扩展卡尔曼滤波法[12]12],加权小二乘法[12]12],非线性小二乘法[13]13]等。在机器人控制方面,Slotine等做了大量深入的研究,提出了直接自适应、间接自适应、复合自适应等控制算法[14,15,16,17]14-17];针对机器人工作环境限制、负载变化等问题,提出将自适应理论结合阻抗控制,并证明了控制系统的稳定性[18]18];焦晓红等[19]19]将自适应结合鲁棒控制,针对误差扰动上界确定和不确定两种情况提出控制律,并通过仿真证明了控制算法的可行性。文献[20]提出了自适应神经网络控制算法。利用径向基神经网络来调整滑模控制增益,消除滑模控制输入抖动并优化系统性能。

笔者以ZCR07S六自由度工业机器人为实验研究对象,基于Codesys平台完成机器人动力学参数辨识实验与自适应控制实验。辨识出一组动力学参数并实现了自适应控制算法。

1 机器人动力学模型与参数辨识

1.1 机器人动力学特性

n连杆刚性机器人动力学方程可由二阶非线性微分方程[21]21]描述为

$Η (q) \ddot{q} + C (q, \dot{q}) \dot{q} + G (q) + τ_{f} = τ$ (1)

式中:q, $\dot{q}$ , $\ddot{q}$ ∈Rn×1分别代表关节角位置、角速度和角加速度向量;τ∈Rn×1为关节驱动力矩;H(q)∈Rn×n为惯性矩阵;C(q, $\dot{q}$ )∈Rn×1为哥氏力和离心力矩阵;G(q)∈Rn×1为重力项;τf∈Rn×1为摩擦力项,本研究中摩擦力模型为 $τ_{f} = f_{v} \dot{q} + f_{c} sgn (\dot{q})$ ,其中fv,fc分别为黏性摩擦系数与库伦摩擦系数,sgn代表符号函数。

机器人系统的动力学结构特性如下:

特性1 可适当定义 $C (q, \dot{q}), \dot{Η} - 2 C$ 为斜对称矩阵,即对任意x∈Rn均有

$x^{Τ} [\dot{Η} - 2 C] x = 0$ (2)

特性2 存在一个依赖于机械臂参数的参数向量a,使得H(q),C(q, $\dot{q}$ ),G(q),τf满足线性关系

$\begin{array}{l} Η (q) \ddot{q} + C (q, \dot{q}) \dot{q} + G (q) + \\ τ_{f} = Y (q, \dot{q}, \ddot{q}) a (3) \end{array}$

式中:Y(q, $\dot{q}$ , $\ddot{q}$ )为n×m的已知关节变量的回归矩阵,它是机器人广义坐标及其各阶导数的已知函数矩阵,a为m×k描述机器人质量特性的未知定常参数向量。即机器人的动力学特性是机器人惯性参数的线性函数。

式(3)中Y的有些列恒为零,有些列间存在常值线性关系,所以不是所有惯性参数都对力矩有影响,通过线性关系重组标准惯性参数[11,22]或数值解析方法[23]获得可辨识的小惯性参数集amin。式(3)可重写为

$\begin{array}{l} Η (q) \ddot{q} + C (q, \dot{q}) \dot{q} + G (q) + τ_{f} = \\ Y (q, \dot{q}, \ddot{q}) a = Y^{'} (q, \dot{q}, \ddot{q}) a_{\min} (4) \end{array}$

式中Y′(q, $\dot{q}$ , $\ddot{q}$ )各列线性无关。六自由度串联机器人小惯性参数如表1所示,惯性参数重组关系为

Izz1r=(m3+m4+m5+m6)(d $_{2}^{2}$ +d $_{3}^{2}$ +r $_{3}^{2}$ )+

Ia1+d $_{2}^{2}$ m2+2mz3r3+Iyy2+Iyy3+Izz1

Ixx2r=-(m3+m4+m5+m6)·

(d $_{3}^{2}$ +r $_{2}^{2}$ -r $_{3}^{2}$ )+Ixx2-Iyy2

Ixz2r=Ixz2-d3(m3+m4+m5+m6+mz3)r3

Izz2r=Ia2+d $_{3}^{2}$ (m3+m4+m5+m6)+Izz2

mx2r=d3(m3+m4+m5+m6)+mx2

Ixx3r=2mz4r4+(m4+m5+m6)r $_{4}^{2}$ +

Ixx3-Iyy3+Iyy4

Izz3r=2mz4r4+(m4+m5+m6)r $_{4}^{2}$ +

Iyy4-Izz3

my3r=my3-mz4-(m4+m5+m6)r4

Ixx4r=Ixx4-Iyy4+Iyy5

Izz4r=Iyy5+Izz4

my4r=Imy4+Imz5

Ixx5r=Ixx5-Iyy5+Iyy6

Izz5r=Iyy6+Izz5

my5r=Imy5-Imz6

Ixx6r=Ixx6-Iyy6

表1 小惯性参数 导出到EXCEL

Table 1 Minimum inertial parameters

连杆1	连杆2	连杆3	连杆4	连杆5	连杆6
0	Ixx2r	Ixx3r	Ixx4r	Ixx5r	Ixx6r
0	Ιxy2	Ιxy3	Ιxy4	Ixy5	Ιxy6
0	Ixz2r	Ixz3r	Ixz4r	Ixz5r	Ιxz6
0	Ιyz2	Ιyz3	Ιyz4	Ιyz5	Ιyz6
Izz1r	Izz2r	Izz3r	Izz4r	Izz5r	Ιzz6
0	mx2r	mx3	mx4	mx5	mx6
0	my2	my3r	my4r	my5r	my6
0	0	Ia3	Ia4	Ia5	Ia6
fv1	fv2	fv3	fv4	fv5	fv6
fc1	fc2	fc3	fc4	fc5	fc6

1.2 动力学参数辨识原理及实验

为实现基于动力学模型的控制,需要对结构已知的模型参数进行辨识,获取机器人的一组小惯性参数值,即式(4)中的amin。机器人的参数辨识是一个复杂的系统工程,包括运动学与动力学建模、激励轨迹的选取、实验设计、数据采集与处理、参数估计、模型验证等环节,其流程图如图1所示。

图1 辨识实验流程图 下载原图

Fig.1 Flow chart of identification experiment

ZCR07S机器人6 个关节均为旋转关节,本体如图2所示,运动学模型的建立采用文献[11]中的方法,连杆坐标系如图3所示,MDH参数如表2所示。

图2 实验机器人本体 下载原图

Fig.2 Experimental robot body

图3 连杆坐标系 下载原图

Fig.3 Connecting rod coordinate system

表2 MDH参数 导出到EXCEL

Table 2 MDH parameters

连杆	ri/mm	αi/(°)	di/mm	θi/(°)
1	0	0	0	θ1
2	0	90	25	θ2+90
3	0	0	455	θ3
4	420	90	35	θ4
5	0	90	0	θ5
6	0	-90	0	θ6

动力学模型由SYMORO+[24]软件推导,以自定义符号表达式形式在robot-name.dyn文件中给出,由于篇幅限制,文中不再给出。

各关节激励轨迹为N次谐波的正弦和余弦函数有限项和,激励轨迹关节角位置表达式为

$\begin{array}{l} q_{i} (t) = \sum_{l = 1}^{Ν} \frac{a_{l}}{ω_{f} l} \sin (ω_{f} l t) - \\ \frac{b_{l}}{ω_{f} l} \cos (ω_{f} l t) + q_{i 0} (5) \end{array}$

关节角速度、角加速度分别为一阶、二阶导数。实验中:基频为0.1 Hz;周期为10 s;N为5;qi0为位置偏移量;al,bl确定了三角函数的幅值。

确定激励轨迹参数al,bl,qi0的问题可描述为非线性约束优化问题:

优化目标为

q*(t)=argmin(J) (6)

约束条件为

式中:J为待确定的目标函数;qmax, $\dot{q}$ max, $\ddot{q}$ max分别为关节限制大位置、速度、加速度。实验所用机器人的关节位置、速度限位如表3所示。

表3 各关节位置、速度限位 导出到EXCEL

Table 3 Position and speed limitation of joints

关节	位置限位/(°)	速度限位/((°)·s-1)
1	±130	150.0
2	-18,+90	112.5
3	-20,+140	150.0
4	±250	204.5
5	±130	225.0
6	±350	360.0

式(7)中后面的3 个约束条件保证了激励轨迹起始和终止位置、速度、加速度为零。根据Hadamard不等式,正定矩阵的行列式小于等于其对角项的乘积,可计算det(YTY)上界,可小化矩阵条件数。目标函数J[25]为

$J = \frac{1}{\prod_{j = 1}^{n} \sum_{i = 1}^{n} Y_{i j}^{2}}$

式中Y为m×n的回归矩阵,Yij为矩阵Y的第i行第j列元素。通过MatlabFunction函数计算目标函数参数,将目标函数与式(7)中的约束条件代入到Matlab fmincon函数计算激励轨迹参数。激励轨迹如图4所示。

图4 激励轨迹 下载原图

Fig.4 Excitation trajectory

辨识实验设计与数据采集通过Codesys平台完成,同时采集关节角位置、速度和力矩,采样频率为250 Hz。由Codesys采集到的力矩经过如下转换关系换算为机械臂端力矩:机械臂端力矩=采集的力矩相对值×电机额定力矩×减速器减速比/1 000。实验所用电机额定力矩和减速器减速比如表4所示。

表4 电机额定力矩和减速器减速比 导出到EXCEL

Table 4 Rated torque of motor and reduction ratio of speed reducer

关节	电机额定力矩/(N·m)	减速器减速比
1	2.390	120
2	2.390	120
3	1.270	100
4	0.318	116
5	0.318	80
6	0.318	50

实验中采集的关节位置轨迹光滑,无需滤波处理,关节速度、加速度通过中心差分算法计算得到,并通过Matlab smooth函数对关节加速度、力矩进行平滑处理。

小惯性参数值的估计采用标准小二乘法,去除掉一些相对标准偏差较大(RSD>30%)的惯性参数后,得到29 个描述6 个关节机器人动力学模型的小惯性参数,如表5所示。

表5 小惯性参数值 导出到EXCEL

Table 5 Minimum inertial parameter values

小惯性参数	$\hat{a}$	RSD/%
Izz1r	39.560 689 257 7	13.15
fv1	31.972 108 039 4	1.26
fc1	16.061 131 863 2	1.25
Ixx2r	-42.671 108 237 0	15.73
Ixz2r	0.347 377 258 1	25.04
Ιyz2	-4.336 088 223 2	28.62
mx2r	9.078 391 556 5	8.02
my2	-0.210 724 339 0	14.43
fv2	-59.361 832 520 0	4.31
fc2	-13.513 298 318 7	1.55
Ixx3r	8.847 127 449 7	26.17
Izz3r	-5.644 584 976 7	17.19
mx3	-0.784 824 398 6	22.06
my3r	-2.284 085 631 8	0.87
Ia3	6.317 742 287 6	24.07
fv3	-11.772 151 253 0	23.69
fc3	-7.879 692 889 6	2.75
Ιyz4	-0.201 137 447 0	24.43
fv4	9.930 327 124 3	3.47
fc4	9.812 399 300 8	2.18
Ixz5r	-0.221 956 985 4	28.06
Ia5	0.039 832 177 0	27.59
fv5	2.717 465 875 1	11.18
fc5	2.461 814 008 5	7.32
Ιxy6	0.087 176 081 2	26.16
Ιyz6	0.138 817 631 4	17.87
Ia6	-0.059 928 086 2	25.28
fv6	0.441 493 161 0	29.04
fc6	0.921 412 730 6	17.76

通过比较实测力矩与由辨识模型所估计的力矩的均方根误差来判断所辨识的模型是否符合精确度要求。图5为重新生成的一组验证轨迹。

图5 验证轨迹 下载原图

Fig.5 Verification trajectory

图6,7分别代表激励轨迹和验证轨迹下的实测力矩与估计力矩。

图6 6 个关节激励轨迹下实测力矩与估计力矩图 下载原图

Fig.6 Measured and estimated torque diagrams under 6 joint excitation trajectories

图7 6 个关节验证轨迹下实测力矩与估计力矩图 下载原图

Fig.7 Measured and estimated torque diagrams under 6 joint verification trajectories

图7 6 个关节验证轨迹下实测力矩与估计力矩图 下载原图

Fig.7 Measured and estimated torque diagrams under 6 joint verification trajectories

表6为激励轨迹和验证轨迹下的实测力矩与估计力矩的均方根误差,误差值较小,辨识的模型可用。

表6 均方根误差 导出到EXCEL

Table 6 RMS error

关节	激励轨迹/(N·m)	验证轨迹/(N·m)
1	4.130 1	7.631 3
2	4.333 2	6.196 1
3	2.370 1	3.413 0
4	2.484 8	2.743 3
5	1.197 2	1.143 3
6	0.470 1	0.475 7

2 自适应控制算法实验

2.1 自适应控制器的实现

选择控制律和自适应律[14]为

$τ = \hat{Η} (q) \ddot{q}_{r} + \hat{C} (q, \dot{q}) \dot{q}_{r} + \hat{G} (q) + \hat{F} (\dot{q}) - Κ_{D} s (8)$

$\dot{\hat{a}} = - Γ^{- 1} Y^{Τ} (q, \dot{q}, \dot{q}_{r}, \ddot{q}_{r}) s$ (9)

式中: $\hat{Η}$ (q), $\hat{C}$ (q, $\dot{q}$ ), $\hat{G}$ (q), $\hat{F}$ ( $\dot{q}$ )分别为惯性矩阵、离心力和哥氏力、重力、摩擦力的估计值,其中 $\hat{f}$ v, $\hat{f}$ c分别为辨识的黏性摩擦系数与库伦摩擦系数;KD,Γ分别为PD控制器增益、自适应增益,均选为正定对角常数矩阵; $\dot{q}$ r, $\ddot{q}$ r分别为参考速度、参考加速度, $\dot{q}_{r} = \dot{q}^{*} - Λ \tilde{q}, \ddot{q}_{r} = \ddot{q}^{*} - Λ \dot{\tilde{q}},$ 其中 $\dot{q}^{*}, \ddot{q}^{*}$ 分别为期望速度、期望加速度;q, $\dot{q}$ 分别为实际位置、实际速度,其中 $\tilde{q} = q - q^{*}, \dot{\tilde{q}} = \dot{q} - \dot{q}^{*}$ ;s为位置误差和速度误差的加权和,其表达式为 $s = \dot{q} - \dot{q}_{r} = \dot{\tilde{q}} + Λ \tilde{q}$ 。由式(4)知:控制律式(8)可写为

$τ = Y (q, \dot{q}, \dot{q}_{r}, \ddot{q}_{r}) \hat{a}_{\min} - Κ_{D} s$ (10)

式中: $\hat{a}$ min为由辨识得到的小惯性参数;回归矩阵Y可直接由符号计算软件SYMORO+得出,其中矩阵元素Yij以自定义形式在文件robot-name.dim中给出,对应于robot-name.dim中的DGiKj。则由得到的小惯性参数 $\hat{a}$ min与回归矩阵Y即可实现自适应控制器。自适应控制器结构如图8所示。

图8 自适应控制器结构图 下载原图

Fig.8 Structure diagram of adaptive controller

2.2 实验装置

研究中使用的六轴工业机器人控制系统采用一主多从的控制模式,主站采用ARM+Linux控制方案,从站采用带EtherCAT接口的高创CDHD系列伺服驱动,主站与从站、从站与从站之间的通讯采用工业以太网EtherCAT技术,软件平台以Codesys作为开发环境,以ST和CFC作为主要编程语言实现机器人动力学控制,机器人控制系统结构如图9所示。

图9 控制系统结构 下载原图

Fig.9 Control system architecture

2.3 实验结果

实验的期望轨迹为

${\begin{cases} q_{1}^{*} = 5 t^{2} - \frac{10}{9} t^{3} 0 \leq t < 3 \\ q_{2}^{*} = 3 t^{2} - \frac{2}{3} t^{3} 0 \leq t < 3 \\ q_{3}^{*} = 5 t^{2} - \frac{10}{9} t^{3} 0 \leq t < 3 \\ q_{4}^{*} = - 10 t^{2} + \frac{20}{9} t^{3} 0 \leq t < 3 \\ q_{5}^{*} = 10 t^{2} - \frac{20}{9} t^{3} 0 \leq t < 3 \\ q_{6}^{*} = 10 t^{2} - \frac{20}{9} t^{3} 0 \leq t < 3 \\ q_{i}^{*} = q_{i}^{*} (3) 3 \leq t, i = 1, 2, \dots, 6 \end{cases}$

程序运行采样周期为4 ms,自适应增益均为1,在实际的调节过程中,根据系统的稳定性与跟踪性能之间的权衡,确定KD=diag[400,350,130,22,9,2],Λ=diag[10,10,10,30,30,25](KP=ΛKD)。实验结果如图10所示,图10中曲线由上至下分别代表期望位置y1,实际位置y2,期望速度y3,实际速度y4,位置误差y5,速度误差y6。跟踪过程中,关节1至关节6大位置误差分别为-0.06°,0.22°,0.34°,-0.29°,0.45°,0.85°。关节1至关节6位置定位误差分别为-0.06°,0.02°,0.04°,-0.11°,-0.11°,0.09°。

图10 关节1至关节6的跟踪结果图 下载原图

Fig.10 Tracking results of joint 1 to joint 6

图10 关节1至关节6的跟踪结果图 下载原图

Fig.10 Tracking results of joint 1 to joint 6

3 结论

采用一种计算效率高的优化准则设计辨识实验的激励轨迹,通过Codesys平台完成辨识实验,去掉相对标准偏差较大的惯性参数后,得到29 个描述机器人动力学特性的小惯性参数,并通过验证轨迹验证了辨识参数的精确性;接着基于所辨识的小惯性参数和机器人动力学模型,将一种自适应控制算法应用于六自由度工业机器人轨迹跟踪控制,通过实验验证了所研究方法的有效性。主要体现为应用创新:第一个创新点是辨识出了ZCR07S工业机器人动力学参数;第二个创新点是将自适应控制算法应用于ZCR07S工业机器人,并取得了良好的控制效果。