2.1. Inverse dynamics

The dynamics of a n-DOF serial robot can be expressed as:

(1) \begin{equation} M(q)\ddot{q} + C(q,\dot{q})\dot{q} + G(q) +{F_{fric}} = \tau +{\tau _{ext}} \end{equation}

where $q \in \mathbb{R}^n$ , $\dot{q} \in \mathbb{R}^n$ , $\ddot{q} \in \mathbb{R}^n$ stand for the joint position, velocity and acceleration, respectively. $M(q) \in \mathbb{R}^{n \times n}$ , $C(q,\dot{q}) \in \mathbb{R}^{n \times n}$ and $G(q) \in \mathbb{R}^{n}$ represent the inertia matrix, the Coriolis/centrifugal matrix and the gravity vector, respectively. ${F_{fric}} \in \mathbb{R}^{n}$ and $\tau \in \mathbb{R}^{n}$ are the friction torque vector and control torque vector of the joint. ${\tau _{ext}} \in \mathbb{R}^{n}$ is the external torque vector due to the physical contact with the environment.

In the case of robot dynamic identification with electric current measured by drive actuation, the joint control torque can be calculated by:

(2) \begin{equation} \tau = KI \end{equation}

where $I\in \mathbb{R}^{n}$ is the measured joint electric current and $K \in \mathbb{R}^{n \times n}$ is a diagonal matrix consisting of the torque constant of motor.

2.2. Friction modeling

It is assumed that the friction torques are uncoupled among the joints and the frictional model form is not changed due to assembly of joint modules for collaborative robots [Reference Madsen, Rosenlund, Brandt and Zhang1]. And frictional torque $F_{fric}$ is frequently described as a combination of three components, which characterize the effects of joint angular velocity, temperature and load torque, respectively [Reference Madsen, Rosenlund, Brandt and Zhang1]. Specifically, the friction torque of the i-th joint can be calculated by using a composite model

(3) \begin{equation}{F_{fric,i}} ={f_v}({\dot{q_i}}) +{f_T}({\dot{q_i}},{T_i}) +{f_\tau }(q,{\dot{q_i}}), |\dot{q_i}| \geq \delta \dot{q_i} \end{equation}

where $T_i$ represents the joint temperature. ${f_v}({\dot{q_i}})$ , ${f_T}({\dot{q_i}},{T_i})$ and ${f_\tau }(q,{\dot{q_i}})$ stand for the dependence of friction torque on velocity, temperature and load torque, respectively. $\delta \dot{q_i}$ is the minimum resolution of velocity data. Three groups of experiments were conducted on collaborative robots with easy-assemblable joints to investigate these terms.

Figure 1.

The velocity and temperature dependence of friction.

The first group of experiments were carried out in each single joint with zero load and at a constant temperature to investigate the velocity-dependent friction model. During these experiments, each joint was placed horizontally and driven with a sinusoidally changing speed. The velocity, temperature and driving moment information of the joints were collected during these experiments as shown in Fig. 1(a). Fig. 1(a) indicates that the friction torque increases exponentially with velocity but is discontinuous at near-zero velocities. At different temperatures, the static friction torque almost keeps unchanged, but the viscous friction torque varies significantly as in [Reference Gao, Yuan, Han, Wang and Wang24, Reference Zeyu, Hongxing, Ziyi and Gang25]. Thus, the components, ${f_v}({\dot{q_i}})$ can be described as:

(4) \begin{equation}{f_v}({{\dot{q}}_i}) = [{f_{c,i}} +{f_{v,i}}{\|{{{\dot{q}}_i}}\|^{{\alpha _i}}}]sign({{\dot{q}}_i}) +{f_{{{vo,i}}}} \end{equation}

where $f_{c,i} \gt 0$ represents the Coulomb friction coefficient. $f_{v,i} \gt 0$ and $\alpha _i \gt 0$ are two coefficients that describe the viscous friction behavior. $f_{{{vo,i}}}$ stands for constant offset.

The second groups of experiments were carried out in each single joint with zero load in an ordinary working environment to investigate the relationship between friction and temperature. A temperature sensor was fixed on the harmonic reducer to measure the joint temperature and each joint was also driven with the same trajectory as the first experiment. In this case, robot joints naturally heat up over 20 $^{\circ }$ C. Fig. 1(b) presents the experimental results, from which it can be seen that temperature-dependent friction is a function of speed and temperature. Therefore, the components ${f_T}({\dot{q_i}},{T_i})$ of friction torque can be expressed as in [Reference Madsen, Rosenlund, Brandt and Zhang1]:

(5) \begin{equation}{f_T}({{\dot{q}}_i},{T_i}) = sign({{\dot{q}}_i})\sqrt{\|{{\dot{q}}_i}\|} ({\beta _{1,i}} +{\beta _{2,i}}{T_i} +{\beta _{3,i}}{T_i}^{ - 3}) \end{equation}

where $\beta _{1,i}$ , $\beta _{2,i}$ and $\beta _{3,i}$ are temperature coefficients.

The third groups of experiments were introduced to model the relationship between friction and load. During these experiments, each joint was placed horizontally, subjected to a known load and moving at a constant speed in both forward and backward directions, and the joint temperature was kept close to ambient temperature. Fig. 2 indicates the experimental results, from which it can be seen that friction torque has a quadratic-form dependence on the payload torque, and its component ${f_\tau }(q,{\dot{q_i}})$ is expressed as:

(6) \begin{equation}{f_\tau }(q,{\dot{q_i}}) ={f_{L,i}}sign({\dot{q_i}}){\tau _{l,i}}{(q)^2} +{f_{{{lo,i}}}} \end{equation}

where ${\tau _{l,i}}{(q)}$ stands for the load torque of the i-th joint. $f_{L,i}$ represents a coefficient that describes the effect of joint payload on friction. $f_{{{lo,i}}}$ stands for constant offset and it can be combined with $f_{{{vo,i}}}$ in Eq. (4) to define a unified zero-speed offset as $f_{{{o,i}}}$ .

Figure 2.

Frictional experiments under different loads with the curves fitted by using Eq. (6).

2.3. Iterative methodology for dynamic identification

When $\tau _{ext} = 0$ , one compensates the components, $f_T$ and $f_\tau$ , which are calculated by Eqs. (5) and (6) – with $\alpha _i$ being regarded as a known constant. Then, Eq. (1) becomes

(7) \begin{equation} M(q)\ddot{q} + C(q,\dot{q})\dot{q} + G(q) +{f_v}(\dot{q}) = \tau ' \end{equation}

and Eq. (1) can be further linearized [Reference Khalil and Kleinfinger30] as:

(8) \begin{equation} \tau ' = \Phi (q,\dot{q},\ddot{q}){\theta _{linear}} \end{equation}

where $\tau '{\rm{ = }}\tau -{\hat f_T}(\dot{q},T) -{\hat f_\tau }(q,\dot{q})$ . $\Phi \in \mathbb{R}^{p \times 13n}$ is a regression matrix, $p$ is the sample times and ${\theta _{linear}} \in \mathbb{R}^{13n}$ is a parameter vector defined as:

(9) \begin{equation} \begin{aligned}{\theta _{linear}} &={[{\theta _D}\;{\theta _F}]^T}\\{\theta _D} &={[m\;{S_x}\;{S_y}\;{S_z}\;{I_{xx}}\;{I_{yy}}\;{I_{zz}}\;{I_{xy}}\;{I_{xz}}\;{I_{yz}}]}\\{\theta _F} &={[{f_c}\;{f_v}\;{f_o}]} \end{aligned} \end{equation}

where $\theta _{linear}$ consists of dynamic inertia parameter vector $\theta _D$ and frictional parameter vector $\theta _F$ . $m$ is the link’s mass. $S$ and $I$ stand for the first and second-order inertia moments, respectively.

One can perform singular value decomposition or QR decomposition operation on the regression matrix to get the full-rank matrix $\Phi _{min}$ and the minimum parameter set ${\theta _{min}} = \{{\theta _{D,\min }},{\theta _{F,\min }}\}$ . $\theta _{D,\min }$ and $\theta _{F,\min }$ are two re-organized parameter sets of $\theta _D$ and $\theta _F$ , respectively. The symbolic expressions of the minimum parameter set are listed in Appendix.

By using the minimum parameter set, the dynamic equation can be recast as:

(10) \begin{equation} \tau ' ={\Phi _{min}}(q,\dot{q},\ddot{q}){\theta _{min}} \end{equation}

For the purpose of reducing the interference of outliers on calculation, the covariance matrix of errors $\Omega$ as defined in [Reference Han, Wu, Liu and Xiong11] is utilized to normalize the regression matrix and respond vector in the first loop (NB: The first loop is the innermost loop denoted as $L_1$ ). Then, the iterative weight vector $W \in{\mathbb{R}^{mn \times 1}}$ is updated with a T-class hard weight function to gradually remove outliers in the second loop (i.e., $L_2$ ). Thus, Eq. (10) becomes

(11) \begin{equation}{\tau ^{\#} } = \Phi _{min}^{\#} (q,\dot{q},\ddot{q}){\theta _{min}} \end{equation}

and

(12) \begin{equation} \begin{aligned}{\tau ^{\#} } &= W.*({\Omega ^{ - \frac{1}{2}}} \cdot \tau ')\\ \Phi _{min}^{\#} &={W_e}.*({\Omega ^{ - \frac{1}{2}}} \cdot{\Phi _{min}}) \end{aligned} \end{equation}

where $W_e \in{\mathbb{R}^{mn \times p}}$ is an extension matrix consisting of p-column vectors $W$ .

When the solution of Eq. (11) is considered as the result of a linear regression process, the solution of $\hat{\theta }$ can be estimated by using LS to minimize the squared residuals. When the physical feasibility of the solution of Eq. (11) is considered, the constraints of parameters can be written in the form of LMI. Then, a convex optimization problem is reached, and it can be solved by using the SDP technique [Reference Sousa and Cortesão12].

In this case, the normalized error is calculated by

(13) \begin{equation}{\mathcal{R}^{\#} } ={\tau ^{\#} } - \Phi _{min}^{\#} (q,\dot{q},\ddot{q}){\hat \theta _{min}} \end{equation}

Then, the covariance matrix is updated by

(14) \begin{equation} \Omega = \frac{{{\Omega ^{\frac{1}{2}}} \cdot{E^{\#} } \cdot{E^{\# T}} \cdot{\Omega ^{\frac{1}{2}}}}}{{m - p}} \end{equation}

where ${E^{\#} } \in{\mathbb{R}^{n \times m}}$ is a reshaped matrix of $R^{\#}$ .

2.4. Improved iterative methodology for friction identification

In the above section, the components of $F_{fric}$ are considered as a priori knowledge used in Eq. (8). For purpose of reducing the interference of this priori on calculation and enhancing the model’s accuracy, the friction-model parameters need to be calculated and updated during the third loop (i.e., $L_3$ ) of I-IRLS.

When the solution $\hat \theta _{min}$ is estimated, the friction torque $F_{fric}$ can be calculated as:

(15) \begin{equation}{F_{fric}}({\theta _{Fric}}) = \tau -{\Phi _{min}}(q,\dot{q},\ddot{q}) \cdot{[{\hat \theta _{D,min}},{\rm{0}}]^T} \end{equation}

where ${\theta _{Fric}} = [{f_c}\;{f_v}\;\alpha \;{f_{\rm{{o}}}}\;{f_L}]$ are the friction-model parameters described in section 2.2. The calculation of $\hat \theta _{Fric}$ is feasible if a solution exists in its physically constrained space. This can be represented as an optimization problem as:

(16) \begin{equation} \begin{aligned} \mathop{{\rm{min\;}}}\limits _{{\theta _{Fric}}} &\|{F_{fric}} -{{\hat F}_{fric}}\|\\{\rm{s.t.\ }}{f_{c,\max }} &\ge{\rm{ }}{f_c} \gt 0\\{f_{v,\max }} &\ge{f_v} \gt 0\\{\alpha _{\max }} &\ge \alpha \gt 0\\{f_{o,\max }} &\ge{f_{o}} \ge{f_{o,\min }}\\{f_{L,\max }} &\ge{f_L} \gt 0 \end{aligned} \end{equation}

where $f_{c,\max }$ , $f_{v,\max }$ , $\alpha _{\max }$ , $f_{o,\max }$ , $f_{o,\min }$ and $f_{L,\max }$ are the boundary values for the parameters’ feasible domain.

The optimization problem in Eq. (16) can be solved by the multi-start interior-point method to get the global optimal solution. In experiments, the duration of the excitation trajectory is too short to reflect the effect of temperature on friction, so the result of the single-joint experiment mentioned in section 2.2 can be adopted to determine the component ${f_T}(\dot{q},T)$ .

Subsequently, the residual error of the dynamic model in Eq. (1) is calculated by

(17) \begin{equation} \begin{aligned} R &= \tau - \hat \tau \\ &= \tau -{{\hat F}_{fric}} -{\Phi _{min}}(q,\dot{q},\ddot{q}) \cdot{[{{\hat \theta }_{D,min}},{\rm{0}}]^T} \end{aligned} \end{equation}

The load torque vector ${\tau _l}(q)$ , which is used to calculate ${f_\tau }(q,\dot{q})$ in Eq. (6), can be replaced by the approximated gravity vector [Reference Gao, Yuan, Han, Wang and Wang24]. By importing the identification result, the load torque vector ${\tau _l}(q)$ can be updated by using the gravity vector during the outermost fourth loop $L_4$ :

(18) \begin{equation}{\hat \tau _l}(q) ={\Phi _{min}}(q,0,0) \cdot{[{\hat \theta _{D,min}}{\rm{ }}, 0]^T} \end{equation}

The I-IRLS friction model can be linearized by considering some model parameters as priori knowledge and estimating other parameters together with the dynamic parameters set. When the friction model cannot be linearized, the I-IRLS can be designed as a decoupled form (abbr. as ID-IRLS), which decouples the dynamic model identification from the friction-model identification simply by redefining $\tau '$ and $\theta _{linear}$ as:

(19) \begin{equation} \begin{aligned} \tau ' &= \tau -{f_v}(\dot{q}) -{f_T}(\dot{q},T) -{f_\tau }(q,\dot{q})\\{\theta _{linear}} &= [m\;{S_x}\;{S_y}\;{S_z}\;{I_{xx}}\;{I_{yy}}\;{I_{zz}}\;{I_{xy}}\;{I_{xz}}\;{I_{yz}}] \end{aligned} \end{equation}

2.5. Excitation trajectory

The excitation trajectory is necessary for identification. It needs to sufficiently excite the robotic system to obtain a well-conditioned regression matrix so that parameters are insensitive to the noise and thus can be fully estimated. When a fifth-order finite Fourier series [Reference Swevers, Ganseman, Tükel, Schutter and a. Brussel31] is applied, the position trajectory of the i-th joint can be written as:

(20) \begin{equation}{q_i}(t) ={q_{i0}} + \sum \nolimits _{k = 1}^{N = 5}{(a_k^i\sin (k{w_f}t)} + b_k^i\cos (k{w_f}t)) \end{equation}

where $w_f$ is a fundamental pulse, $q_{i0}$ is an offset, and $a_k^i$ and $b_k^i$ are the amplitudes of harmonic sine and cosine functions, respectively.

For the purpose of calculating the coefficients in Eq. (20), the regression matrix is divided into two sub-matrices according to the significance of parameters, where $\Phi _1$ and $\Phi _2$ are related to the first-order and second-order moments, respectively. Then, an optimization problem can be formulated in the form of a weighted function as follows:

(21) \begin{equation} \begin{aligned} \mathop{{\rm{min\;}}}\limits _{{q_{i0}},a_{k = 1 \ldots,5}^i,b_{k = 1 \ldots,5}^i} &W_1 cond(\Phi _1) + W_2 cond(\Phi _2)\\{\rm{s.t.\ }}{q_{\min }} &\le{q} \le{q_{\max }}\\{{\dot{q}}_{\min }} &\le{{\dot{q}}} \le{{\dot{q}}_{\max }}\\{{\ddot{q}}_{\min }} &\le{{\ddot{q}}} \le{{\ddot{q}}_{\max }}\\{\tau _{{\min }}} &\le \Gamma{(q,\dot{q},\ddot{q})} \le{\tau _{\max }}\\{r_{\min }} &\le r_{{\rm{ee}}} (q_{\rm{end}}) \le{r_{\max }}\\{h_{\min }} &\le h_{{\rm{ee}}} (q_{\rm{end}}) \le{h_{\max }} \end{aligned} \end{equation}

where $W_1$ and $W_2$ are the weight coefficients of the condition numbers of matrices $\Phi _1$ and $\Phi _2$ , respectively. $q_{min}$ , $q_{max}$ , $\dot{q}_{min}$ , $\dot{q}_{max}$ , $\ddot{q}_{min}$ , $\ddot{q}_{max}$ , $\tau _{min}$ and $\tau _{max}$ stand for the ranges of angle, velocity, acceleration and driving torque of the joints. $r_{min}$ , $r_{max}$ , $h_{min}$ and $h_{max}$ represent the ranges of radius and height of the robotic end effector.

Algorithm 1.

I-IRLS

2.6. Overall framework

Algorithm 1 shows how the model parameters are estimated. In $L_1$ , dynamic parameters are estimated by using the normalized observation matrix $\Phi _{min}^{\#}$ and the normalized observation response $\tau ^{\#}$ in Eqs. (11) and (12). In $L_2$ , the weight matrix $W$ is updated by the normalized residuals $R^{\#}$ calculated from Eq. (13) so that the outliers in the measured data are eliminated.

In the present work, IRLS is improved as I-IRLS by modifying $L_3$ and adding $L_4$ . In $L_3$ , a comprehensive model that considers the effects of velocity, temperature and load torque is developed to describe the friction torque of joints, and its parameters are identified based on a globally optimal algorithm with constraints by using Eqs. (15) and (16). In this way, the interference of a priori knowledge on the identification results can be reduced. Since the friction components and the coefficient $\alpha$ are considered as a priori knowledge for calculating the observed response $\tau '$ and the observation matrix $\Phi _{min}$ , they should be re-computed from Eqs. (8)-(10) after the new solutions $\hat{\theta }_F$ are obtained.

In $L_4$ , the joint load torque ${\hat \tau _l}(q)$ is updated. Since the ${\hat \tau _l}(q)$ vector can be approximated by the gravity vector $G(q)$ in Eq. (1), its accuracy is also influenced by the estimation of dynamics parameters set $\theta _D$ and friction torque $F_{fric}$ . It has to be re-computed after the new solutions of $\hat{\theta }_D$ and $\hat{\theta }_F$ are obtained from Eq. (18). Therefore, the components $f_T$ and $f_\tau$ are calculated and compensated in $L_3$ and $L_4$ , and this effectively avoids their influence on the linearization of dynamics in $L_2$ .