# Dynamics Parametrization and Calibration of Flexible-Joint Collaborative Industrial Robot Manipulators

Authors: Emil Madsen; Simon Aagaard Timm; Norbert Andras Ujfalusi; Oluf Skov Rosenlund; David Brandt; Xuping Zhang
Journal: Mathematical Problems in Engineering
 (2020)

---

## Abstract

Many collaborative robots use strain-wave-type transmissions due to their desirable characteristics of high torque capacity and low weight. However, their inherent complex and nonlinear behavior introduces significant errors and uncertainties in the robot dynamics calibration, resulting in decreased performance for motion and force control tasks and lead-through programming applications. This paper presents a new method for calibrating the dynamic model of collaborative robots. The method combines the known inverse dynamics identification model with the weighted least squares (IDIM-WLS) method for rigid robot dynamics with complex nonlinear expressions for the rotor-side dynamics to obtain increased calibration accuracy by reducing the modeling errors. The method relies on two angular position measurements per robot joint, one at each side of the strain-wave transmission, to effectively compensate the rotor inertial torques and nonlinear dynamic friction that were identified in our previous works. The calibrated dynamic model is cross-validated and its accuracy is compared to a model with parameters obtained from a CAD model. Relative improvements are in the range of 16.5% to 28.5% depending on the trajectory.


---

## Body

## 1. Introduction

For collaborative industrial robots, it is of crucial importance to acquire accurate predictions of the torques required in order to realize the desired motion or force control task and to ensure a consistently good performance of lead-through programming applications. Being able to accurately predict the torques required to complete the intended task will (1) improve the control performance by being able to react to disturbances before they cause deviations from the reference and (2) improve the robot safety system by being able to more accurately identify external disturbances such as human interference. Accurate torque estimates can be obtained through knowledge about the dynamic properties of the robot. Accurate torque estimates will also improve any possible online estimation procedures such as online estimation of the payload mass and inertia properties [1], friction, and/or wear.The dynamic model of the robot relates the robot motion to the joint torques and it depends on a set of dynamic parameters being the mass, the first moments, and the mass moments of inertia of each link of the robot. Multiple procedures exist for estimating the dynamic parameters of robot manipulators:(1)
Physical Experiments. The robot is disassembled to isolate each link. The mass can be evaluated directly. The first moments can be obtained by evaluating the counterbalanced points of each link. The diagonal elements of the inertia tensor can be evaluated by pendular motions. Such methods are tedious and are not preferred because they require a lot of manual operations to disassemble the robot and carry out the experiments. Furthermore, experiments need to be redone if hardware changes are made to the robot.(2)
Computer Aided Design (CAD). The dynamic parameters of each link are found using their nominal geometric and material characteristics. In the design phase, such investigation can be used in the performance analysis to further improve the design. However, the accuracy of the parameter estimates is reduced because the CAD parts are never identical to the real parts due to the production tolerances.(3)
System Identification. The input/output behavior is analyzed on some planned motion. Parameters are estimated by minimizing the difference between the measured output (possibly the current supplied to the electric actuator) and its mathematical model evaluated in the input (possibly the angular positions of the robot joints). Such procedures are preferred because they generally lead to the most accurate results while offering flexibility in the case of robot hardware changes.For system identification methods, the most common strategy is a combined Inverse Dynamics Identification Model and Least Squares (IDIM-LS) method. For such method, the accuracy of the parameter estimates is generally affected by measurement noise and modeling errors.The issue of measurement noise is often addressed by generating so-called exciting trajectories and/or filtering the noisy measurements [2]. Other identification techniques have also been suggested such as the Extended Kalman Filter (EKF) [3, 4], algorithms based on Linear Matrix Inequality (LMI) tools [5], maximum likelihood (ML) approaches [6], the Set Membership Uncertainty [7], and Huber’s estimator [8]. However, based on the experimental results, these approaches do not improve the IDIM-LS, and they were not validated on 6-degrees-of-freedom (DOF) industrial robots. To eliminate the need for tuning the bandpass filters that are applied to the trajectory data, [9, 10] used the Instrumental Variable (IV) technique, and [11] proposed the Direct and Inverse Dynamic Identification Models (DIDIM) technique. These methods are based on a closed-loop output error (CLOE) method using both the direct and inverse dynamic models of the robot. The direct dynamic model is used to obtain model-based estimates of the position, velocity, and acceleration signals in contrast to the bandpass filtering often coupled with the IDIM-LS method. In [12], the DIDIM and CLOE methods were compared to the IDIM-LS method and it was found that if the IDIM-LS method is coupled with well-tuned bandpass filtering, the DIDIM and CLOE methods do not offer any improvements to the IDIM-LS method. Other methods include identifying the dynamics of a robotics system using neural networks [13].Modeling errors will generally lead to a bias of the parameter estimates and it is an issue yet unsolved in the system identification for industrial robots. Modeling errors arises mainly from neglecting the complex and nonlinear joint dynamics effects resulting in significant deterministic structural errors that cannot be accounted for by random variables. Such nonlinear joint dynamics come, for instance, due to the use of strain-wave type transmissions such as the Harmonic Drive™ which are often used in collaborative robots due to their desirable characteristics of high torque capacity and low weight.The works on the identification of dynamic parameters for collaborative robots are limited. In [14], the essential parameters were identified for the KUKA LWR 4+ collaborative robot assuming a three-parameter friction model. In [1, 15], dynamics parameter identification was performed using the KUKA LWR 4+ collaborative robot with friction neglected. The works on the KUKA LWR 4+ collaborative robot exploited the joint torque sensor located on the output side of the transmission; thus the joint dynamics do not affect the measurements. Such sensor hardware is, however, expensive and is rarely found in industrial robots. In [16], the dynamic parameters for the 7 DOF Franka Emika Panda robot were identified with a constrained optimization procedure to ensure the physical consistency of the parameters. In [17], the parameters for the 2×7 DOF ABB IRB 14000 (YuMi) collaborative robot were identified. The fact that very simple models of the friction are employed with Coulomb and linear viscous friction and that joints are assumed rigid are common in the mentioned works. Such assumptions on the joint dynamics characteristics for strain-wave transmissions are serious simplifications of the real dynamic characteristics.To address the mentioned limitations of the prior art, we propose a new method for estimating the dynamic parameters of collaborative robot manipulators considering the flexible joint dynamics effects. Firstly, the dynamics of the Universal Robots UR5e collaborative robot manipulator are developed in closed form using the modified Denavit-Hartenberg convention and the Recursive Newton–Euler Algorithm. Secondly, the dynamic equations are linearly parametrized and the dimension of the parameter space is reduced to a minimum. Thirdly, the proposed rotor dynamics compensation is introduced to reduce modeling errors. The novelty lies in the rotor dynamics compensation in whichtwo built-in rotary encoders are utilized per joint, one at each side of the transmission element, to effectively compensate the complex nonlinear joint dynamics effects of the Universal Robots UR5e robot manipulator identified prior to this work [18,19]. Any unmodeled friction is handled by augmenting the set of dynamic parameters with Coulomb and viscous friction coefficients for each joint. The parameters are then estimated by a WLS procedure with the weighting equal to the inverse of the estimated covariance matrix. Lastly, the calibrated dynamic model is validated on new trajectories that were not used for the estimation (cross-validation principle). The general methodology of the dynamics calibration in this work is illustrated in Figure 1. The two-encoder setup is illustrated by the Robot outputting two angular position variables q and θ for the link and rotor, respectively.Figure 1
Schematic representation of the methodology of this work showing the interconnection between the identification and validation procedures.The distinguished contributions of this work include the following: (1) a linear parametrization describing the dynamics of the UR5e collaborative robot manipulator has been developed, (2) the complex nonlinear dynamic friction characteristics and rotor inertia have been considered, (3) the minimal set of base parameters that describe the dynamic behavior of the UR5e robot has been accurately estimated, and (4) the performance of the calibrated dynamic model has been validated.The rest of the paper is organized as follows: Section2 describes the mathematical model of the Flexible-Joint Robot (FJR) manipulator. In Section 3, the linear parametrization of the dynamics is described. In Section 4, the identification procedure is described and the results of the identification are presented. In Section 5, the calibrated dynamic model is validated and compared to a model obtained with parameters from a CAD model. Section 6 concludes the work and presents the challenges for our ongoing research.


### 1.1. Notations

The notation used in the paper is mostly standard. Letℝ be the set of real numbers, ℕ be the set of nonnegative integers, and ℕ+ be the set of positive integers. Let x∈ℝn be a vector of n real numbers; then xi is its ith entry, xT its transpose, x¯ is the mean value of the elements of x, and x is the 2-norm. Let x^ denote an estimate of x and let x˜≜x−x^ be the estimation error. Given a function g:G⟶ℝ, let sgn:ℝ⟶−1,0,1 be the signum function defined such that sgng=−1 if g<0, sgng=0 if g=0, and sgng=1 if g>0. If g:G⟶ℝn is a vector function, the signum vector function sgng=sgng1⋯sgngnT. Given a square real matrix A∈ℝn×n, let A≻0 indicate that A is positive definite; that is, xTAx>0 for any nonzero column vector x of n real numbers. Let diag:ℝn⟶ℝn×n map a vector of n elements to a diagonal matrix with the ith element of the vector on its ith diagonal entry and zero everywhere else. Similarly, let diag−1:ℝn×n⟶ℝn map the diagonal elements of an n×n matrix to vector of n elements with the ith diagonal element of the matrix on the ith element of the vector.


## 1.1. Notations

The notation used in the paper is mostly standard. Letℝ be the set of real numbers, ℕ be the set of nonnegative integers, and ℕ+ be the set of positive integers. Let x∈ℝn be a vector of n real numbers; then xi is its ith entry, xT its transpose, x¯ is the mean value of the elements of x, and x is the 2-norm. Let x^ denote an estimate of x and let x˜≜x−x^ be the estimation error. Given a function g:G⟶ℝ, let sgn:ℝ⟶−1,0,1 be the signum function defined such that sgng=−1 if g<0, sgng=0 if g=0, and sgng=1 if g>0. If g:G⟶ℝn is a vector function, the signum vector function sgng=sgng1⋯sgngnT. Given a square real matrix A∈ℝn×n, let A≻0 indicate that A is positive definite; that is, xTAx>0 for any nonzero column vector x of n real numbers. Let diag:ℝn⟶ℝn×n map a vector of n elements to a diagonal matrix with the ith element of the vector on its ith diagonal entry and zero everywhere else. Similarly, let diag−1:ℝn×n⟶ℝn map the diagonal elements of an n×n matrix to vector of n elements with the ith diagonal element of the matrix on the ith element of the vector.


## 2. Mathematical Model

The Flexible-Joint Robot (FJR) manipulator is considered as an open kinematic chain havingN+1 rigid bodies; the base and the N links are interconnected by N revolute joints undergoing deflection and actuated by N electrical actuators. To derive the dynamics of the robot manipulator, the following standard assumptions are made:(i)
The rotors are uniform bodies having their center of mass on the axis of rotation(ii)
Each motori=1,…,N is mounted on link i−1 and moves link i; see Figure 2Assumption (i) is a basic requirement for long life of an electrical drive and implies that the robot dynamics become independent of the angular position of the rotors. For the UR5e manipulator, we take advantage of the presence of large reduction ratios and simply assume the following.(iii)
The angular velocity of the rotors is due only to their own spinningFigure 2
Kinematic arrangement of motors and links for the FJR manipulator model. Note thatθi is already scaled by the reduction ratio.This simplifying assumption was proposed by [20] and is equivalent to neglecting energy contributions due to the inertial couplings between the rotors and the links. It also implies that Coriolis and centripetal terms will be independent of the rotors’ angular velocity.To uniquely characterize the manipulator configuration, we choose the generalized coordinatesqθ∈ℝ2N being, respectively, the positions of the links and rotors reflected through the gear ratios; that is, the rotor positions are seen in the link space. Given assumptions (i)–(iii), the link and rotor dynamics become, respectively,(1)Mqq¨+Cq,q˙q˙+gq=τJ,(2)Bθ¨+f+τJ=Kτi,where, in the link equation, Mq≻0∈ℝN×N is the symmetric inertia matrix, Cq,q˙∈ℝN×N is the Coriolis and centripetal matrix, gq∈ℝN is the gravity vector, and τJ∈ℝN is the vector of joint torques which couple the link and rotor subsystems. In the rotor equation, B≻0∈ℝN×N is the diagonal matrix of rotor inertias, f∈ℝN is friction acting on the rotor coordinate, Kτ≻0∈ℝN×N is the diagonal matrix of torque constants, and i∈ℝN is the torque-generating (quadrature) current obtained from the phase currents via Park’s Transformation. The drive-gain Kτ has been calibrated a priori with special tests; see, for example, [21].The dynamic model of theN link robot manipulator is obtained in closed form using the Denavit-Hartenberg (DH) convention [22] with coordinate systems placed as illustrated in Figure 3 to represent the UR5e manipulator and with the parameters in Table 1 and the Recursive Newton–Euler Algorithm (RNEA) [22]. In the RNEA, the position vector Pii+1=ai−1−disinαi−1dicosαi−1T and the rotation matrix(3)Rii+1=cosqi−sinqi0sinqicosαi−1cosqicosαi−1−sinαi−1sinqisinαi−1cosqisinαi−1cosαi−1.Figure 3
Coordinate systems used for describing the kinematics of the UR5e robot using the Denavit–Hartenberg convention.Table 1
Denavit–Hartenberg parameters of the UR5e robot.
iaiαidiqi10π/2d1q12a200q23a300q340π/2d4q450−π/2d5q5600d6q6To allow a parametrization of the dynamics which is linear in the inertial parameters, the inertia tensor for each link is defined relative to the center of rotation (CoR). The RNE algorithm needs the inertia tensor defined relative to the center of mass (CoM), so the parallel axis theorem (Steiner’s law) is used for translation; that is,(4)ICoM,i=ICoR,i−miPC,iTPC,iE3−PC,iPC,iT,with E3 being the 3×3 identity matrix, and the vector of center of mass positions PC,i=PC,i,xPC,i,yPC,i,zT, mi is the mass of link i, and the symmetric inertia tensor is(5)ICoR,i=XXiXYiXZiXYiYYiYZiXZiYZiZZi.


## 3. Dynamics Parametrization

The expressions for the torques obtained from the RNE algorithm can be expressed linearly in the inertial parameters:(6)τJ=Yq,q˙,q¨γ,where the inertial parameters are(7)γ=γ1T⋯γiT⋯γNTT,γi=XXiXYiXZiYYiYZiZZimXimYimZimiT.For a specific robot manipulator, not all10N inertial parameters can be identified. Not all the inertial parameters have an effect on the dynamic model, while others have an effect only in linear combinations. The inertial parameters of a robot can therefore be classified into three groups: fully identifiable, identifiable in linear combinations only, and unidentifiable. This is due to the kinematic arrangement of the joints as well as the orientation of the manipulator’s base with respect to gravity. Table 2 shows the 49 inertial parameters that appear in the mathematical model.Table 2
The 49 inertial parameters that appear in the dynamic model of the UR5e manipulator when mounted such that the rotation axis of the base joint is oriented parallel to the gravitational acceleration.
iXXiXYiXZiYYiYZiZZimXimYimZimi1–––––ZZ1––––2XX2XY2XZ2YY2YZ2ZZ2m2PC,2,xm2PC,2,y––3XX3XY3XZ3YY3YZ3ZZ3m3PC,3,xm3PC,3,ym3PC,3,zm34XX4XY4XZ4YY4YZ4ZZ4m4PC,4,xm4PC,4,ym4PC,4,zm45XX5XY5XZ5YY5YZ5ZZ5m5PC,5,xm5PC,5,ym5PC,5,zm56XX6XY6XZ6YY6YZ6ZZ6m6PC,6,xm6PC,6,ym6PC,6,zm6For the estimation problem to have a unique solution, the parameters must be linearly independent. The set of linearly independent parameters is calledbase parameters. The number of base parameters bm is [23](8)bm≤7nr+4np−3−σ¯1−2ng0,where nr is the number of revolute joints, np is the number of prismatic joints, and σ¯1=1 if joint 1 is revolute; otherwise σ¯1=0; and ng0=1 if the rotation axis of joint 1 is parallel to the direction of the gravitational acceleration; otherwise ng0=0. For a robot manipulator with the kinematic arrangement in Table 1 and the base joint oriented with its rotation axis parallel to the direction of the gravitational acceleration, the number of base parameters bm=36. Therefore, a number of inertial parameters are grouped into a fewer number of equivalent parameters using the regrouping relations [24]:(9)XXRi=XXi−YYi,XXRi−1=XXi−1+YYi+2dimZi+di2mi,XYRi−1=XYi−1+ai−1sinαimZi+ai−1disinαimi,XZRi−1=XZi−1−ai−1cosαimZi−ai−1dicosαimi,YYRi−1=YYi−1+cos2αiYYi+2dicos2αimZi+αi−12+di2cos2αimi,YZRi−1=YZi−1+cosαisinαiYYi+2dicosαisinαimZi+di2cosαisinαimi,ZZRi−1=ZZi−1+sin2αiYYi+2disin2αimZi+ai−12+di2sin2αimi,mXRi−1=mXi−1+ai−1mi,mYRi−1=mYi−1−sinαimZi−disinαimi,mZRi−1=mZi−1+cosαimZi+dicosαimi,mRi−1=mi−1+mi.This results in the set of base parameters in Table3.Table 3
The 36 base parameters for the UR5e manipulator when mounted such that the rotation axis of the base joint is oriented parallel to the gravitational acceleration.
iXXiXYiXZiYZiZZimXimYi1––––ZZR1––2XXR2XY2XZR2YZ2ZZR2mXR2mY23XXR3XY3XZ3YZ3ZZR3mXR3mY34XXR4XY4XZ4YZ4ZZR4mX4mYR45XXR5XY5XZ5YZ5ZZR5mX5mYR56XXR6XY6XZ6YZ6ZZ6mX6mY6The joint torque is expressed as(10)τJ=Ybq,q˙,q¨γb,with the vector of base parameters(11)γb=ZZR1XXR2⋯mY6T∈ℝbm.


### 3.1. Including the Rotor Dynamics

Combining (2) and (10) yields(12)Kτi−Bθ¨−f=Ybq,q˙,q¨γb.Friction torques are considered as a sum of estimates and error terms. From experience and empirical observations, the errorf˜≜f−f^ is assumed to contain Coulomb and viscous friction contributions; that is, f˜=F˜Csgnθ˙+F˜Vθ˙. The nonlinear estimates, as presented in [19], describe the friction torques in terms of the angular velocities, load torques, and temperatures. Rotor inertias B are considered to be known; however, they could be easily estimated by augmenting the regressor with the angular acceleration of the rotors. The system formulated in terms of base parameters and augmented with the rotor dynamics and friction discrepancy is(13)τJ=YBq,q˙,q¨γB+ρ,(14)τJ=Kτi−Bθ¨−f^,(15)YBq,q˙,q¨=Ybq,q˙,q¨diagsgnθ˙diagθ˙,(16)γB=γbTdiag−1F˜CTdiag−1F˜VTT,where the noise ρ due to model errors and measurement noise is assumed to have zero mean, be serially uncorrelated, and be heteroskedastic, that is, having a diagonal covariance matrix.


## 3.1. Including the Rotor Dynamics

Combining (2) and (10) yields(12)Kτi−Bθ¨−f=Ybq,q˙,q¨γb.Friction torques are considered as a sum of estimates and error terms. From experience and empirical observations, the errorf˜≜f−f^ is assumed to contain Coulomb and viscous friction contributions; that is, f˜=F˜Csgnθ˙+F˜Vθ˙. The nonlinear estimates, as presented in [19], describe the friction torques in terms of the angular velocities, load torques, and temperatures. Rotor inertias B are considered to be known; however, they could be easily estimated by augmenting the regressor with the angular acceleration of the rotors. The system formulated in terms of base parameters and augmented with the rotor dynamics and friction discrepancy is(13)τJ=YBq,q˙,q¨γB+ρ,(14)τJ=Kτi−Bθ¨−f^,(15)YBq,q˙,q¨=Ybq,q˙,q¨diagsgnθ˙diagθ˙,(16)γB=γbTdiag−1F˜CTdiag−1F˜VTT,where the noise ρ due to model errors and measurement noise is assumed to have zero mean, be serially uncorrelated, and be heteroskedastic, that is, having a diagonal covariance matrix.


## 4. Identification

This section presents the experimental setup, identification procedure, and results. The experiment is carried out using the setup shown in Figure4. The system consists of the UR5e collaborative robot manipulator, Teach Pendant, Control Box, and PC. The estimation trajectory is generated using the Teach Pendant and is sent to the Control Box. The Control Box generates torque commands and sends them to the UR5e, and the measurements of actual values (q, θ, and i) are sent back from the UR5e to the Control Box. All the data are then logged by the PC.Figure 4
Experimental setup used for the dynamics identification.The identification procedure is illustrated schematically in Figure5.Figure 5
Schematic representation of the Inverse Dynamics Identification Model with Weighted Least Squares (IDIM-WLS) estimation androtor dynamics compensation procedure. The estimation trajectory sends desired angular positions qd,q˙d,q¨d to the controller, which generates torque commands Kτi to the robot. The angular positions of the links and rotors, q and θ, respectively, are measured and then filtered to generate smoothed estimates q^ and time-derivatives q^˙, q^¨, θ^˙, and θ^¨. The rotor quantities are passed to the Rotor dynamics compensation, and the link quantities are passed to the Linearized dynamic model. The Rotor dynamics compensation augments the measured current with the rotor dynamics based on (15) to estimate joint torques τJ. The Linearized dynamic model generates regressors W¯i for each joint i. The joint torque estimates and regressors are filtered and downsampled to generate, respectively, the measurement vector y and observation matrix W. The measurement vector and observation matrix are passed to the OLS Estimation and WLS estimation procedures. The OLS procedure produces the weighting matrix G used in the WLS procedure. Finally, the WLS estimation provides the dynamic parameters estimates γ^B.Data is sampled at timestk=kTS,k=1,2,…,M, where TS=1ms is the sampling period, and the sampling frequency fS=1kHz.


### 4.1. Joint Position, Velocity, and Acceleration Estimation

The measured trajectory dataq and θ are filtered by a 4th-order Butterworth filter in both the forward and reverse directions to eliminate lag of the filtered trajectories q^ and θ^. To keep the useful signal of the robot dynamics in the filter bandwidth, the cutoff frequency of the filter is chosen to be 5 times the frequency of the robot dynamics; that is, 5fdyn=50Hz. Angular velocities and accelerations q^˙, q^¨, θ^˙, and θ^¨, respectively, are obtained through a central difference procedure. The combination of the two-pass Butterworth filter and central difference is referred to as the band-pass filtering process.


### 4.2. Parallel Filtering and Downsampling

The sampling frequency is much higher than the frequencies of interest in the dynamics, so to reduce the required computational resources the data is parallel-filtered and then decimated/downsampled. Firstly, the samplesk=1,…,M are ordered in the measurement vector yi and observation matrix Wi for each joint i=1,…,N individually; that is,(17)y¯i=τJ,i,1⋯τJ,i,k⋯τJ,i,MT,W¯i=Wi,1T⋯Wi,kT⋯Wi,MTT,Wi,k=YB,iq^k,q^˙k,q^¨k,with YB,iq^k,q^˙k,q^¨k being the ith row of the regressor evaluated in the kth sample of the filtered trajectory. The parallel filtering of the measurement vector and observation matrix for each joint is conducted by passing the signals through a 4th-order Butterworth filter in both the forward and reverse directions having a cut-off frequency of 2fdyn=20Hz. The downsampling factor is 0.8fS/4fdyn=20 [10]; that is, every 20th sample is used for parameter estimation. The filtering and downsampling of y¯i and W¯i produce estimates yi and Wi, respectively.


### 4.3. Torque Computation

The filtered and downsampled data are ordered joint-wise in the measurement vector and observation matrix as(18)y=y1T⋯yiT⋯yNTT∈ℝN⋅M,W=W1T⋯WiT⋯WNTT∈ℝN⋅M×bm.The base parameters are estimated by solving the WLS problem:(19)γ^B=argminγBWTGy−WγB2=WTGW−1WTGy,where each weight in G is equal to the reciprocal of the estimated standard deviation of the error.(20)G=diagS,S=S1…Si…SN,Si=1σ^i,1…1σ^i,j…1σ^i,bm,i,j=1,…,bm,i,σ^i,j2=τi−YB,iγ^B,i2M−bm,i,with bm,i being the number of base parameters related to link i. Such weighting operation normalizes the error terms in (13).


### 4.4. Trajectory

The trajectory used for parameter estimation should allow complete identification of the system; that is, for positive constantsα and β, it should satisfy some persistently exciting condition:(21)βE≥∫0TWTWdt≥αE,where α is the degree of excitation and β/α is the condition number of WTW. The trajectory in Figure 6 is used for yielding a condition number for the regressor condW=66. The trajectory is 31 seconds long; hence, with a sampling frequency of 1000 Hz and a downsampling factor of 20, the total number of samples is 1550. Another approach to the trajectory design is to optimize the condition number of the regressor with respect to the trajectory subject to kinematic and dynamic constraints, for example, position, velocity, acceleration, and current.Figure 6
Trajectory used for the parameter estimation yielding a low condition number of the regressor matrixcondW=66.


### 4.5. Model Quality Metric

The model quality is evaluated by the sum of each joint’s mean squared error normalized by its average torque magnitude; that is,(22)NMSE=∑j=1N∑k=1MτJ,i,k−τ^J,i,k2τJ,i,k.


### 4.6. Results

Values of the identified base parameters are shown in Table4. The effectiveness of the method is demonstrated by considering the accuracy of the dynamic model with the optimized parameters compared to our baseline, a model with parameters obtained through CAD software. The model accuracy improves 81.4% from a NMSE of 506.8 Nm to that of 94.5 Nm.Table 4
Estimated values and relative standard deviation of the 36 base parameters obtained by solving the WLS problem with rotor dynamics compensation.
Base param.Value%σγ^BZZR12.19810.099XXR2−1.65520.130XY2−0.04051.975XZR20.37340.210YZ20.12781.367ZZR2−4.17740.198mXR2−0.01320.001mY2−0.73360.141XXR3−0.15630.166XY30.13020.225XZ3−0.17140.249YZ30.02120.433ZZR3−1.87742.908mXR30.04280.001mY3−0.04280.030XXR4−0.04540.609XY40.01440.422XZ40.01660.767YZ40.000518.894ZZR4−0.000520.082mX4−0.01120.043mYR4−0.18660.001XXR50.02390.450XY5−0.03910.055XZ50.01940.080YZ50.00650.286ZZR50.03150.076mX5−0.00130.137mYR50.04420.005XXR60.01870.132XY60.00590.107XZ60.01870.054YZ6−0.00023.825ZZ60.02090.157mX60.00360.040mY6−0.00430.041The parameters are generally well estimated with small relative standard deviations, which demonstrates the effectiveness of the identification procedure. The values of parametersYZ4 and ZZR4 are, however, subject to relative standard deviations of 19% and 20%, respectively. This is an indication of either (1) suboptimality of the chosen trajectory or (2) the parameter being of no big value to the torque computation. YZ4 is a product moment of inertia (off-diagonal element in the inertia tensor) and is therefore likely to be less important in the dynamics. ZZR4 is a mass moment of inertia (diagonal element in the inertia tensor) and the reduced accuracy in its estimation is likely due to insufficient excitation by the chosen trajectory.


## 4.1. Joint Position, Velocity, and Acceleration Estimation

The measured trajectory dataq and θ are filtered by a 4th-order Butterworth filter in both the forward and reverse directions to eliminate lag of the filtered trajectories q^ and θ^. To keep the useful signal of the robot dynamics in the filter bandwidth, the cutoff frequency of the filter is chosen to be 5 times the frequency of the robot dynamics; that is, 5fdyn=50Hz. Angular velocities and accelerations q^˙, q^¨, θ^˙, and θ^¨, respectively, are obtained through a central difference procedure. The combination of the two-pass Butterworth filter and central difference is referred to as the band-pass filtering process.


## 4.2. Parallel Filtering and Downsampling

The sampling frequency is much higher than the frequencies of interest in the dynamics, so to reduce the required computational resources the data is parallel-filtered and then decimated/downsampled. Firstly, the samplesk=1,…,M are ordered in the measurement vector yi and observation matrix Wi for each joint i=1,…,N individually; that is,(17)y¯i=τJ,i,1⋯τJ,i,k⋯τJ,i,MT,W¯i=Wi,1T⋯Wi,kT⋯Wi,MTT,Wi,k=YB,iq^k,q^˙k,q^¨k,with YB,iq^k,q^˙k,q^¨k being the ith row of the regressor evaluated in the kth sample of the filtered trajectory. The parallel filtering of the measurement vector and observation matrix for each joint is conducted by passing the signals through a 4th-order Butterworth filter in both the forward and reverse directions having a cut-off frequency of 2fdyn=20Hz. The downsampling factor is 0.8fS/4fdyn=20 [10]; that is, every 20th sample is used for parameter estimation. The filtering and downsampling of y¯i and W¯i produce estimates yi and Wi, respectively.


## 4.3. Torque Computation

The filtered and downsampled data are ordered joint-wise in the measurement vector and observation matrix as(18)y=y1T⋯yiT⋯yNTT∈ℝN⋅M,W=W1T⋯WiT⋯WNTT∈ℝN⋅M×bm.The base parameters are estimated by solving the WLS problem:(19)γ^B=argminγBWTGy−WγB2=WTGW−1WTGy,where each weight in G is equal to the reciprocal of the estimated standard deviation of the error.(20)G=diagS,S=S1…Si…SN,Si=1σ^i,1…1σ^i,j…1σ^i,bm,i,j=1,…,bm,i,σ^i,j2=τi−YB,iγ^B,i2M−bm,i,with bm,i being the number of base parameters related to link i. Such weighting operation normalizes the error terms in (13).


## 4.4. Trajectory

The trajectory used for parameter estimation should allow complete identification of the system; that is, for positive constantsα and β, it should satisfy some persistently exciting condition:(21)βE≥∫0TWTWdt≥αE,where α is the degree of excitation and β/α is the condition number of WTW. The trajectory in Figure 6 is used for yielding a condition number for the regressor condW=66. The trajectory is 31 seconds long; hence, with a sampling frequency of 1000 Hz and a downsampling factor of 20, the total number of samples is 1550. Another approach to the trajectory design is to optimize the condition number of the regressor with respect to the trajectory subject to kinematic and dynamic constraints, for example, position, velocity, acceleration, and current.Figure 6
Trajectory used for the parameter estimation yielding a low condition number of the regressor matrixcondW=66.


## 4.5. Model Quality Metric

The model quality is evaluated by the sum of each joint’s mean squared error normalized by its average torque magnitude; that is,(22)NMSE=∑j=1N∑k=1MτJ,i,k−τ^J,i,k2τJ,i,k.


## 4.6. Results

Values of the identified base parameters are shown in Table4. The effectiveness of the method is demonstrated by considering the accuracy of the dynamic model with the optimized parameters compared to our baseline, a model with parameters obtained through CAD software. The model accuracy improves 81.4% from a NMSE of 506.8 Nm to that of 94.5 Nm.Table 4
Estimated values and relative standard deviation of the 36 base parameters obtained by solving the WLS problem with rotor dynamics compensation.
Base param.Value%σγ^BZZR12.19810.099XXR2−1.65520.130XY2−0.04051.975XZR20.37340.210YZ20.12781.367ZZR2−4.17740.198mXR2−0.01320.001mY2−0.73360.141XXR3−0.15630.166XY30.13020.225XZ3−0.17140.249YZ30.02120.433ZZR3−1.87742.908mXR30.04280.001mY3−0.04280.030XXR4−0.04540.609XY40.01440.422XZ40.01660.767YZ40.000518.894ZZR4−0.000520.082mX4−0.01120.043mYR4−0.18660.001XXR50.02390.450XY5−0.03910.055XZ50.01940.080YZ50.00650.286ZZR50.03150.076mX5−0.00130.137mYR50.04420.005XXR60.01870.132XY60.00590.107XZ60.01870.054YZ6−0.00023.825ZZ60.02090.157mX60.00360.040mY6−0.00430.041The parameters are generally well estimated with small relative standard deviations, which demonstrates the effectiveness of the identification procedure. The values of parametersYZ4 and ZZR4 are, however, subject to relative standard deviations of 19% and 20%, respectively. This is an indication of either (1) suboptimality of the chosen trajectory or (2) the parameter being of no big value to the torque computation. YZ4 is a product moment of inertia (off-diagonal element in the inertia tensor) and is therefore likely to be less important in the dynamics. ZZR4 is a mass moment of inertia (diagonal element in the inertia tensor) and the reduced accuracy in its estimation is likely due to insufficient excitation by the chosen trajectory.


## 5. Validation

The purpose of the dynamic model calibration is to improve the torque estimation accuracy for arbitrary trajectories. Thus, we evaluate the accuracy of our calibrated model on three trajectories different from the one used for parameter estimation. The measured joint torques are compared to the torques output by the calibrated dynamic model as well as relative to our CAD model baseline. Improvements in NMSE and improvements relative to our baseline are shown in Table5. The results show a relative improvement of the calibrated dynamic model compared to the dynamic model with CAD parameters of 16.5%–28.5% depending on the trajectory. Time-series torque data for each of the joints of the UR5e dynamic model with calibrated parameters are shown in Figures 7–12 for trajectory no. 1. Time-series torque data for each of the joints of the UR5e dynamic model with CAD model parameters are shown in Figures 13–18 for trajectory no. 1. The reduction in torque prediction error (NMSE) of 16.5%–28.5% of the calibrated dynamic model compared to the dynamic model with CAD model parameters together with the time-series torque data in Figures 7–18 validates the effectiveness of the calibration procedure.Table 5
Normalized Mean Squared Error (NMSE) of the dynamic models; CAD model with parameters obtained from a CAD model of the robot and calibrated with parameters estimated through the WLS procedure with rotor dynamics compensation. Improvements in NMSE and relative to CAD model for three different trajectories.
Traj.NMSE CAD modelNMSE calibratedRelative improvement (%)1325.9233.028.521245.71040.416.53334.0250.924.9Figure 7
Results of measurement and estimate using the calibrated parameters on the validation trajectory.Figure 8
Results of measurement and estimate for joint no. 2 with the calibrated parameters.Figure 9
Results of measurement and estimate for joint no. 3 with the calibrated parameters.Figure 10
Results of measurement and estimate for joint no. 4 with the calibrated parameters.Figure 11
Results of measurement and estimate for joint no. 5 with the calibrated parameters.Figure 12
Results of measurement and estimate for joint no. 6 with the calibrated parameters.Figure 13
Results of measurement and estimate for joint no. 1 with the CAD model parameters.Figure 14
Results of measurement and estimate for joint no. 2 with the CAD model parameters.Figure 15
Results of measurement and estimate for joint no. 3 with the CAD model parameters.Figure 16
Results of measurement and estimate for joint no. 4 with the CAD model parameters.Figure 17
Results of measurement and estimate for joint no. 5 with the CAD model parameters.Figure 18
Results of measurement and estimate for joint no. 6 with the CAD model parameters.


## 6. Conclusion

Collaborative industrial robots often utilize strain-wave type transmissions due to their desirable characteristics of high torque capacity and low weight. However, their inherent complex nonlinear behavior introduces significant errors and uncertainties in the robot dynamics calibration, resulting in decreased performance for motion and force control tasks and lead-through programming applications.This paper presented a new method for the dynamics parametrization and calibration of collaborative industrial robot manipulators. The method combines the IDIM-WLS method for rigid robot dynamics with complex nonlinear expressions for the rotor-side dynamics to obtain increased calibration accuracy. Two angular position measurements per robot joint are utilized, one at each side of the strain-wave transmission, to effectively compensate the rotor inertial torques and nonlinear dynamic friction that were identified in our previous works.The effectiveness of the method was demonstrated by the application to the Universal Robots UR5e collaborative robot manipulator. The results were very accurate estimates of the dynamic parameters. Relative improvement of 16.5% to 28.5% compared to a CAD model baseline was experienced.The distinguished contributions of this work can be summarized as follows: (1) a linear parametrization describing the dynamics of the UR5e collaborative robot manipulator has been developed, (2) the complex nonlinear dynamic friction characteristics and rotor inertia have been considered, (3) the minimal set of base parameters that describe the dynamic behavior of the UR5e robot has been accurately estimated, and (4) the performance of the calibrated dynamic model has been validated.Our ongoing work that we are going to challenge consists of the following:(1)
The number of identifiable parameters can be increased by two (from 36 to 38) if the robot is mounted with the base joint axis of rotation not being parallel to the direction of the gravitational acceleration.(2)
The trajectory used for parameter optimization could possibly be optimized through the use of some optimization procedures. Such trajectory optimization procedures were discussed generally in [25] and applied in [26] on a KUKA 361 IR industrial robot.(3)
The optimization problem could be constrained to enforce the positive-definiteness of the inertia matrix using either Sylvester’s theorem or Cholesky decomposition. This will ensure invertibility of the inertia matrix, which is useful for model-based control design. From Sylvester’s theorem, it is possible to find conditions for the parameters [27], whereas, for Cholesky decomposition, a tolerance is defined and it takes the noise and measurement error into account [28].



---

*Source: 8709870-2020-09-17.xml*