#### 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.

#### 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 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 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. 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 and for the link and rotor, respectively.

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: 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.

##### 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 be a vector of real numbers; then is its ^{th} entry, its transpose, is the mean value of the elements of , and is the 2-norm. Let denote an estimate of and let be the estimation error. Given a function , let be the signum function defined such that if , if , and if . If is a vector function, the signum vector function . Given a square real matrix , let indicate that is positive definite; that is, for any nonzero column vector of real numbers. Let map a vector of elements to a diagonal matrix with the ^{th} element of the vector on its ^{th} diagonal entry and zero everywhere else. Similarly, let map the diagonal elements of an matrix to vector of elements with the ^{th} diagonal element of the matrix on the ^{th} element of the vector.

#### 2. Mathematical Model

The Flexible-Joint Robot (FJR) manipulator is considered as an open kinematic chain having rigid bodies; the base and the links are interconnected by revolute joints undergoing deflection and actuated by 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 motor is mounted on link and moves link ; see Figure 2 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)The angular velocity of the rotors is due only to their own spinning

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 coordinates 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, is the symmetric inertia matrix, is the Coriolis and centripetal matrix, is the gravity vector, and is the vector of joint torques which couple the link and rotor subsystems. In the rotor equation, is the diagonal matrix of rotor inertias, is friction acting on the rotor coordinate, is the diagonal matrix of torque constants, and is the torque-generating (quadrature) current obtained from the phase currents via Park’s Transformation. The drive-gain has been calibrated *a* priori with special tests; see, for example, [21].

The dynamic model of the 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 and the rotation matrix

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). 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,with being the identity matrix, and the vector of center of mass positions , is the mass of link , and the symmetric inertia tensor is

#### 3. Dynamics Parametrization

The 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 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.

For the estimation problem to have a unique solution, the parameters must be linearly independent. The set of linearly independent parameters is called *base parameters*. The number of base parameters is [23]where is the number of revolute joints, is the number of prismatic joints, and if joint 1 is revolute; otherwise ; and if the rotation axis of joint 1 is parallel to the direction of the gravitational acceleration; otherwise . 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 . Therefore, a number of inertial parameters are grouped into a fewer number of equivalent parameters using the regrouping relations [24]:

This results in the set of base parameters in Table 3.

The joint torque is expressed aswith the vector of base parameters

##### 3.1. Including the Rotor Dynamics

Friction torques are considered as a sum of estimates and error terms. From experience and empirical observations, the error is assumed to contain Coulomb and viscous friction contributions; that is, . The nonlinear estimates, as presented in [19], describe the friction torques in terms of the angular velocities, load torques, and temperatures. Rotor inertias 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 iswhere 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 Figure 4. 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 (, , and ) are sent back from the UR5e to the Control Box. All the data are then logged by the PC.

The identification procedure is illustrated schematically in Figure 5.

Data is sampled at times , where is the sampling period, and the sampling frequency .

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

The measured trajectory data and are filtered by a 4^{th}-order Butterworth filter in both the forward and reverse directions to eliminate lag of the filtered trajectories 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, . Angular velocities and accelerations , , , 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 samples are ordered in the measurement vector and observation matrix for each joint individually; that is,with being the ^{th} row of the regressor evaluated in the ^{th} 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 4^{th}-order Butterworth filter in both the forward and reverse directions having a cut-off frequency of . The downsampling factor is [10]; that is, every 20^{th} sample is used for parameter estimation. The filtering and downsampling of and produce estimates and , respectively.

##### 4.3. Torque Computation

The filtered and downsampled data are ordered joint-wise in the measurement vector and observation matrix as

The base parameters are estimated by solving the WLS problem:where each weight in is equal to the reciprocal of the estimated standard deviation of the error.with being the number of base parameters related to link . 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:where is the degree of excitation and is the condition number of . The trajectory in Figure 6 is used for yielding a condition number for the regressor . 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.

##### 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,

##### 4.6. Results

Values of the identified base parameters are shown in Table 4. 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.

The parameters are generally well estimated with small relative standard deviations, which demonstrates the effectiveness of the identification procedure. The values of parameters and 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. is a product moment of inertia (off-diagonal element in the inertia tensor) and is therefore likely to be less important in the dynamics. 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 Table 5. 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.

#### 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].

#### Data Availability

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

#### Conflicts of Interest

The authors declare that they have no conflicts of interest reported in this paper.

#### Acknowledgments

This work was supported by the company Universal Robots A/S, Odense, Denmark, and Innovation Fund Denmark (Ref. no. 7038-00058B).