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

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. ,is paper presents a new method for calibrating the dynamic model of collaborative robots. ,e 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. ,e 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. ,e 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.


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. e 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. e robot is disassembled to isolate each link. e mass can be evaluated directly. e first moments can be obtained by evaluating the counterbalanced points of each link. e 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). e 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. e 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. e 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. ese methods are based on a closed-loop output error (CLOE) method using both the direct and inverse dynamic models of the robot. e 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. e 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. e 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. e 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 strainwave 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.
irdly, the proposed rotor dynamics compensation is introduced to reduce modeling errors. e novelty lies in the rotor dynamics compensation in which two 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. e 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 (crossvalidation principle). e general methodology of the dynamics calibration in this work is illustrated in Figure 1. e two-encoder setup is illustrated by the Robot outputting two angular position variables q and θ for the link and rotor, respectively. e 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. e rest of the paper is organized as follows: Section 2 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.

Notations.
e notation used in the paper is mostly standard. Let R be the set of real numbers, N be the set of nonnegative integers, and N + be the set of positive integers. Let x ∈ R n be a vector of n real numbers; then x i is its i th entry, x T 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 ⟶ R, let sgn: R ⟶ −1, 0, 1 { } be the signum function defined such that sgn(g) � −1 if g < 0, sgn(g) � 0 if g � 0, and sgn(g) � 1 if g > 0. If g: G ⟶ R n is a vector function, the signum vector function sgn(g) � [sgn(g 1 ) · · · sgn(g n )] T . Given a square real matrix A ∈ R n×n , let A ≻ 0 indicate that A is positive definite; that is, x T A x > 0 for any nonzero column vector x of n real numbers. Let diag: R n ⟶ R n×n map a vector of n elements to a diagonal matrix with the i th element of the vector on its i th diagonal entry and zero everywhere else. Similarly, let diag − 1 : R n×n ⟶ R n map the diagonal elements of an n × n matrix to vector of n elements with the i th diagonal element of the matrix on the i th element of the vector. Assumption (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) e angular velocity of the rotors is due only to their own spinning is 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.

Mathematical Model
To uniquely characterize the manipulator configuration, we choose the generalized coordinates (q θ) ∈ R 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, where, in the link equation, M(q) ≻ 0 ∈ R N×N is the symmetric inertia matrix, C(q, _ q) ∈ R N×N is the Coriolis and centripetal matrix, g(q) ∈ R N is the gravity vector, and τ J ∈ R N is the vector of joint torques which couple the link and rotor subsystems. In the rotor equation, B ≻ 0 ∈ R N×N is the diagonal matrix of rotor inertias, f ∈ R N is friction acting on the rotor coordinate, K τ ≻ 0 ∈ R N×N is the diagonal matrix of torque constants, and i ∈ R N is the torque-generating (quadrature) current obtained from the phase currents via Park's Transformation. e drive-gain K τ has been calibrated a priori with special tests; see, for example, [21]. e dynamic model of the N 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 igure 1: Schematic representation of the methodology of this work showing the interconnection between the identification and validation procedures. Figure 2: Kinematic arrangement of motors and links for the FJR manipulator model. Note that θ i is already scaled by the reduction ratio.
To 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). e 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, with E 3 being the 3 × 3 identity matrix, and the vector of center of mass positions P C,i � P C,i,x P C,i,y P C,i,z T , m i is the mass of link i, and the symmetric inertia tensor is

Dynamics Parametrization
e expressions for the torques obtained from the RNE algorithm can be expressed linearly in the inertial parameters: where the inertial parameters are For a specific robot manipulator, not all 10 N 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. e inertial parameters of a robot can therefore be classified into three groups: fully identifiable, identifiable in linear combinations only, and unidentifiable. is 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.
For the estimation problem to have a unique solution, the parameters must be linearly independent. e set of linearly independent parameters is called base parameters. e number of base parameters b m is [23] b where n r is the number of revolute joints, n p is the number of prismatic joints, and σ 1 � 1 if joint 1 is revolute; otherwise σ 1 � 0; and n g0 � 1 if the rotation axis of joint 1 is parallel to the direction of the gravitational acceleration; otherwise n g0 � 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 b m � 36. erefore, a number of inertial parameters are grouped into a fewer number of equivalent parameters using the regrouping relations [24]: Mathematical Problems in Engineering is results in the set of base parameters in Table 3. e joint torque is expressed as with the vector of base parameters 3.1. Including the Rotor Dynamics. Combining (2) and (10) yields Friction torques are considered as a sum of estimates and error terms. From experience and empirical observations, the error f ≜ f − f is assumed to contain Coulomb and viscous friction contributions; that is, e 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. e system formulated in terms of base parameters and augmented with the rotor dynamics and friction discrepancy is 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.

Identification
is section presents the experimental setup, identification procedure, and results. e experiment is carried out using the setup shown in Figure 4. e system consists of the UR5e collaborative robot manipulator, Teach Pendant, Control Box, and PC. e estimation trajectory is generated using the Teach Pendant and is sent to the Control Box.
e 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. e identification procedure is illustrated schematically in Figure 5. Table 2: e 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. Data is sampled at times t(k) � kT S , k � 1, 2, . . . , M, where T S � 1 ms is the sampling period, and the sampling frequency f S � 1 kHz.

Joint Position, Velocity, and Acceleration Estimation.
e measured trajectory data q and θ are filtered by a 4 thorder 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, 5f dyn � 50 Hz. Angular velocities and accelerations _ q, € q, _ θ, and € θ, respectively, are obtained through a central difference procedure. e combination of the two-pass Butterworth filter and central difference is referred to as the band-pass filtering process.

Parallel Filtering and Downsampling.
e 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 samples k � 1, . . . , M are ordered in the measurement vector y i and observation matrix W i for each joint i � 1, . . . , N individually; that is, with Y B,i (q k , _ q k , € q k ) being the i th row of the regressor evaluated in the k th sample of the filtered trajectory. e parallel filtering of the measurement vector and observation matrix for each joint is conducted by passing the signals through a 4 th -order Butterworth filter in both the forward and reverse directions having a cut-off frequency of 2f dyn � 20 Hz. e downsampling factor is 0.8f S /(4f dyn ) � 20 [10]; that is, every 20 th sample is used for parameter estimation. e filtering and downsampling of y i and W i produce estimates y i and W i , respectively.

Torque Computation.
e filtered and downsampled data are ordered joint-wise in the measurement vector and observation matrix as e base parameters are estimated by solving the WLS problem: where each weight in G is equal to the reciprocal of the estimated standard deviation of the error.
with b m,i being the number of base parameters related to link i. Such weighting operation normalizes the error terms in (13).

Trajectory.
e 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: where α is the degree of excitation and β/α is the condition number of W T W. e trajectory in Figure 6 is used for yielding a condition number for the regressor cond(W) � 66. e 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.

Model Quality Metric.
e model quality is evaluated by the sum of each joint's mean squared error normalized by its average torque magnitude; that is,  Mathematical Problems in Engineering effectiveness of the identification procedure. e values of parameters YZ 4 and ZZR 4 are, however, subject to relative standard deviations of 19% and 20%, respectively. is is an indication of either (1) suboptimality of the chosen trajectory or (2) the parameter being of no big value to the torque computation. YZ 4 is a product moment of inertia (off-diagonal element in the inertia tensor) and is therefore likely to be less important in the dynamics. ZZR 4 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.

Validation
e purpose of the dynamic model calibration is to improve the torque estimation accuracy for arbitrary trajectories.
us, we evaluate the accuracy of our calibrated model on three trajectories different from the one used for parameter estimation. e 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 Table 5. e results show a relative improvement of the calibrated dynamic model compared to the dynamic model

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. is paper presented a new method for the dynamics parametrization and calibration of collaborative industrial robot manipulators. e 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 strainwave transmission, to effectively compensate the rotor inertial torques and nonlinear dynamic friction that were identified in our previous works. e effectiveness of the method was demonstrated by the application to the Universal Robots UR5e collaborative robot manipulator. e 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.
e 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) e 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) e 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) e optimization problem could be constrained to enforce the positive-definiteness of the inertia matrix using either Sylvester's theorem or Cholesky decomposition. is 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].

Data Availability
Data are not available due to confidentiality and third-party rights.