1. Introduction
Accurate payload dynamics benefit robot applications [Reference Hirzinger, Sporer, Alin Albu-Schaffer, Krenn, Pascucci and Schedl1]. For instance, accurate payload dynamics can provide strong priors for robot planning control [Reference Mavrakis and Stolkin2, Reference Marko Bjelonic, de Viragh, Dhionis Sako, Jenelten and Hutter3], safe interaction [Reference Cherubini, Passama, Crosnier, Lasnier and Fraisse4], and manual guidance [Reference Duque, Prieto and Hoyos5, Reference Canfield, Owens and Zuccaro6].
Generally, the payload parameters need to be identified experimentally. In addition, aiming at some applications in industrial environments such as picking different payload frequently [Reference Bai, Lian, Liu, Wang and Liu7], the identifications have to be provided online [Reference Ding, Zhou, Li and Rong8]. In this area, Abiko et al. [Reference Abiko and Yoshida9] proposed two online identification methods for the payload parameters. One is the vibration motion method, which is to construct the identification model based on the velocity deviation of the base. However, this method needs to measure both linear velocity and angular velocity of the base. The other is the reaction force method, which is to construct the identification model based on the six-dimensional force deviation at the end of the robot. It should be noted that this method requires the installation of an additional force/torque sensor. Kubus et al. [Reference Kubus, Kroger and Wahl10] developed an online method based on the recursive total least squares (RTLS) method [Reference Feng, Zhang, Chang and Zheng11]. Different from the recursive least squares (RLS) [Reference Sánchez, Torres-Torriti and Cheein12, Reference Ma, Wu, Wang, Yi and Liang13], this method can take into account measurement errors to obtain more robust estimation results. In order to cope with the influence of payloads on the balance control of quadruped robots, Tournois et al. [Reference Tournois, Focchi, Prete, Orsolino, Caldwell and Semini14] proposed two iterative-based online identification algorithms. However, neither two algorithms consider the friction model. Farsoni et al. [Reference Farsoni, Landi, Ferraguti, Secchi and Bonfè15, Reference Farsoni, Landi, Ferraguti, Secchi and Bonfe16] proposed an identification method based on Kalman filter [Reference Liu, Wang, Zha, Guo, Jiang and Sun17] and RTLS. Firstly, the linear acceleration, angular velocity, and angular acceleration of the payload are estimated by the Kalman filter based on multi-rate quaternion [Reference Armesto, Tornero and Vincze18]. Then, the inertial parameters of the payload are identified online and in real time by the RTLS. Later, this team proposed an online payload identification method based on collision-free path planning and motion decoupling [Reference Farsoni, Ferraguti and Bonfè19]. When generating recognition trajectory, it regards the operator as an obstacle and uses a mature obstacle avoidance algorithm [Reference Khan, Li and Luo20] to avoid the operator. Additionally, it requires information from visual sensors. Renner et al. [Reference Renner, Wind and Sawodny21] proposed an optimization method based on ordinary least squares method for online payload identification. This method constructs the optimization objective function by minimizing the torque errors [Reference Xu, Fan, Chen, Ng, Ang, Fang, Zhu and Zhao22] and uses the RLS to solve the payload parameters online.
Another way to obtain payload parameters is by the dynamic identification method. In order to obtain good dynamic parameter identification, the design of a suitable excitation trajectory [Reference Chang, An and Ma23, Reference Verdel, Bastide, Vignais, Bruneau and Berret24] is essential to fulfill the continuous excitation conditions for the observation matrix. A classic trajectory is the Fourier series excitation trajectory designed by Swevers et al. [Reference Swevers, Ganseman, Tukel, de Schutter and Van Brussel25]. Optimal trajectory parameters are often obtained by optimizing the condition number of the observation matrix [Reference Park26]. Recently, some studies have demonstrated that incentive effects can be improved by optimizing the sub-matrix condition number [Reference Venture, Ayusawa and Nakamura27, Reference Bonnet, Fraisse, Crosnier, Gautier, González and Venture28]. Additionally, considering measurement noise covariance during optimization and conducting the process iteratively can better capture the dynamic characteristics of the robot [Reference Han, Wu, Liu and Xiong29]. Furthermore, Scalera et al. [Reference Scalera, Nainer, Giusti and Gasparetto30] proposed an interval arithmetic method to deal with dynamics uncertainties in alternative to the online identification of dynamic parameters. Swevers et al. [Reference Swevers, Verdonck, Naumer, Pieters and Biber31] used a maximum likelihood method to estimate the payload parameters based on periodic excitation. The limitation of this method is that it introduces recognition errors as it ignores the dynamics of the last step. To solve this problem, a classical method was proposed in refs. [Reference Swevers, Verdonck and De Schutter32, Reference Khalil, Gautier and Lemoine33]. It distinguished the contribution of the robot itself and the payload to the joint torques. The robot link parameters without load were considered as the prior information for the online payload identification. However, the estimation results in this method are biased. Later, this method was summarized in ref. [Reference Hollerbach, Khalil and Gautier34]. Salah et al. [Reference Salah, Sandoval, Ghiss, Laribi and Zeghloul35] reproduced the online method in refs. [Reference Swevers, Verdonck and De Schutter32, Reference Khalil, Gautier and Lemoine33] for medical applications. Afterward, Gaz et al. [Reference Gaz and De Luca36] proposed a static method to identify the payload parameters. However, this method cannot obtain the inertia tensor of the payload because it lacks the dynamic excitation. Rodriguez et al. [Reference Rodriguez, Ibanez and Battle37] proposed an algebraic method for online payload identification. This method assumes that the total mass of both the link and the payload is concentrated at the end of the link. Then, it transforms the dynamic model to the frequency domain for calculation and identification. Unfortunately, this method is only suitable for a single-link system. In order to achieve accurate motion control, Hu et al. [Reference Hu, Li, Chen and Yao38] proposed an adaptive control algorithm based on generalized momentum for online payload estimation. However, because it sets physical constraints and eliminates many unnecessary dynamic parameters, the accuracy of the torque reconstruction for the robot is insufficient. Thereby, it can increase the cumulative errors of the payload estimation. Recently, Xu et al. [Reference Xu, Fan, Fang, Zhu and Zhao39] have proposed an accurate offline payload identification method based on double weighting technique, but this method cannot yet be applied to online payload identification because of its computation burden. The methods discussed above are all applicable to serial robots. The dynamic modeling, evaluation, and analysis about parallel robots can be found in refs. [Reference Wu, Yu, Gao and Wang40–Reference Wu, Ye, Yu and Huang42].
The authors presented a table as below to compare the proposed work with similar research and state-of-the-art studies identified in the literature review. As shown in the Table I, it can be seen that the payload mass estimation error of the proposed method is the smallest one with respect to other alternative approaches. Khalil et al. [Reference Khalil, Gautier and Lemoine33] and others [Reference Hu, Li, Chen and Yao38] used the simple linear friction model. The friction model was ignored or not used in refs. [Reference Gaz and De Luca36], [Reference Farsoni, Landi, Ferraguti, Secchi and Bonfè15], [Reference Sánchez, Torres-Torriti and Cheein12], and [Reference Kurdas, Hamad, Vorndamme, Mansfeld, Abdolshah and Haddadin43] because the joint torque sensors or F/T sensors were used in these literature. To the best of the authors’ knowledge, a nonlinear friction model (2) is applied to online payload estimation for the first time in this article. The nonlinear friction model can effectively improve the accuracy of robot dynamic parameter identification, thereby improving the accuracy of payload identification (especially for the mass). The methods in refs. [Reference Farsoni, Landi, Ferraguti, Secchi and Bonfè15] and [Reference Sánchez, Torres-Torriti and Cheein12] need external sensors (such as IMU, F/T sensor) for payload estimation. However, other methods only use the proprioceptive sensors to identify payload parameters, which can reduce the cost. Based on the convergence time data presented in the literature, the proposed method has the least convergence times which are about 6.59 s.
After the comparisons, the originality, novelty, and the advantage and disadvantage with respect to alternative approaches are summarized here. The originality and novelty of this article is to derive a symbolic relationship between the parameter difference and the payload parameters and then apply it to the online payload estimation. In addition, a nonlinear friction model is used for online payload estimation for the first time. The main advantage of the proposed method is to achieve higher online estimation accuracy of the payload mass
$m_L$
. In order to increase the online estimation accuracy, on the one hand, a nonlinear friction model (2) was used instead of the linear one. On the other hand, the commanded joint trajectory signals were introduced as the computation input instead of the actual ones. The main disadvantage of the proposed method is that it cannot identify the accurate payload inertia tensor enough online
$[I_{xx,L}, I_{xy,L}, I_{xz,L}, I_{yy,L}, I_{yz,L}, I_{zz,L}]$
, which is consistent with the disadvantage of these alternatively online approaches by using proprioceptive sensors only.
Finally, based on the comparisons, the main contributions of this work are summarized as follows:
1) An online payload identification method without using external sensors is proposed, which is based on the symbolic relationship between robot dynamic parameter differences and payload parameters.
2) Compared with other methods, the mass estimation accuracy of a payload can be significantly improved by using the proposed method.
3) A nonlinear friction model (2) is applied to online payload estimation for the first time to increase the identification accuracy in the proposed method.
Comparisons between the proposed method and other methods.
The remainder of this article is organized as follows. Section 2 presents the preliminaries. Section 3 introduces the proposed method. Section 4 provides the experimental results. Finally, Section 5 concludes this article.
2. Preliminaries
The inverse dynamics of the robot with
$n$
rigid links by considering the motor dynamics can be written as follows:
\begin{equation} \boldsymbol{\tau }=(\boldsymbol{M}(\boldsymbol{q})+\boldsymbol{I}_{a})\ddot{\boldsymbol{q}}+\boldsymbol{C}(\boldsymbol{q},\dot{\boldsymbol{q}})\dot{\boldsymbol{q}}+\boldsymbol{g}(\boldsymbol{q})+\boldsymbol{\tau }_{f}, \end{equation}
where
$\boldsymbol{q},\dot{\boldsymbol{q}},\ddot{\boldsymbol{q}}\in \Re ^{n \times 1}$
are joint position, joint velocity, and joint acceleration of robot, respectively. The
$n$
is the degrees of freedom of the robot.
$\boldsymbol{M}(\boldsymbol{q}) \in \Re ^{n \times n}$
and
$\boldsymbol{I}_{a} \in \Re ^{n \times n}$
are the inertia matrix and the equivalent moment of rotor, respectively. The gravity vector is
$\boldsymbol{g}(\boldsymbol{q}) \in \Re ^{n \times 1}$
. The centrifugal and Coriolis matrix is
$\boldsymbol{C}(\boldsymbol{q},\dot{\boldsymbol{q}}) \in \Re ^{n \times n}$
. The joint torque vector is
$\boldsymbol{\tau } \in \Re ^{n \times 1}$
. The friction torque vector is
$\boldsymbol{\tau }_{f} \in \Re ^{n \times 1}$
.
Compared with the classic linear friction model
$\tau _{f,j}=F_{c,j}\textrm{sign}(\dot{q}_j)+F_{v,j} \dot{q}_j + B_j$
[Reference Zhang and Zhang44], the nonlinear friction model is more conducive to improving the identification accuracy of the robot dynamics [Reference Han, Wu, Liu and Xiong29, Reference Farhat, Mata, Page and Valero45]. It is called a linear friction model because the model can be transformed into a linear product of pure velocity terms and friction parameters, such as
$\tau _{f,j}=\begin{bmatrix} \textrm{sign}(\dot{q}_j) & \dot{q}_j & 1 \end{bmatrix} \cdot \begin{bmatrix} F_{c,j} & F_{v,j} & B_j \end{bmatrix} ^{\mathrm{T}}$
. However, for a nonlinear friction model, the nonlinear friction parameter will be included in the velocity term. For example, a nonlinear friction model was used in this article as
\begin{equation} \begin{aligned} \tau _{f,j}&=(F_{c,j}+F_{v,j}\lvert \dot{q}_j \rvert ^{\alpha _{j}})\textrm{sign}(\dot{q}_j)+B_j\\[5pt] &= \begin{bmatrix} \textrm{sign}(\dot{q}_j) & \lvert \dot{q}_j \rvert ^{\alpha _{j}}\textrm{sign}(\dot{q}_j) & 1 \end{bmatrix} \cdot \begin{bmatrix} F_{c,j} & F_{v,j} & B_j \end{bmatrix} ^{\mathrm{T}}, \end{aligned} \end{equation}
where
$\tau _{f,j}$
is the friction torque of joint
$j$
in
$\boldsymbol{\tau }_{f}$
.
$\alpha _{j}$
is the nonlinear friction parameter of joint
$j$
in
$\boldsymbol{\alpha }$
.
$F_{c,j}$
is the Coulomb friction parameter of joint
$j$
.
$F_{v,j}$
is the viscous friction parameter of joint
$j$
.
$B_j$
is the friction torque offset of joint
$j$
.
$F_{c,j}$
,
$F_{v,j}$
, and
$B_j$
are all the linear friction parameters.
Then (1) can be rewritten by considering (2) as follows:
\begin{equation} \boldsymbol{\tau }=\boldsymbol{W}(\boldsymbol{q},\dot{\boldsymbol{q}},\ddot{\boldsymbol{q}},\boldsymbol\alpha )\boldsymbol\varphi, \end{equation}
where
$\boldsymbol{\varphi }\in \Re ^{14n \times 1}$
is the completed robot dynamic parameters. For
$j$
-th joint,
$\boldsymbol{\varphi }^j$
is equal to
$[\boldsymbol{l}_j,m_{j}\boldsymbol{r}_j,m_j,I_{a,j},F_{c,j},F_{v,j},B_j]^{\mathrm{T}}$
.
$\boldsymbol{l}_{j}$
is the inertia tensor vector which is equal to
$[I_{xx,j},I_{xy,j},I_{xz,j},I_{yy,j},I_{yz,j},I_{zz,j}]$
.
$I_{xx,j}, I_{yy,j},$
and
$I_{zz,j}$
are the moment of inertia of link
$j$
.
$I_{xy,j}, I_{xz,j},$
and
$I_{yz,j}$
are the products of inertia of link
$j$
. It can form the inertia tensor matrix
${{\boldsymbol{I}}_{j}}\in{{\Re }^{3\times 3}}$
with respect to frame
$j$
.
${{m}_{j}}{{\boldsymbol{r}}_{j}}$
is the first moment vector which is equal to
$\left [{{m}_{j}}{{r}_{x,j}},{{m}_{j}}{{r}_{y,j}},{{m}_{j}}{{r}_{z,j}} \right ]$
. The link mass is
$m_j$
.
${r}_{x,j}$
,
${r}_{y,j}$
, and
${r}_{z,j}$
are the position of link center of mass (CoM) on
$X$
axis,
$Y$
axis, and
$Z$
axis, respectively.
According to Gautier and Khalil and others [Reference Gautier and Khalil46–Reference Abedinifar, Ertugrul and Arguz49],
$\boldsymbol{\varphi }$
has some unidentifiable parameters, and they can be identified by the combined form. After deleting or regrouping these parameters, (3) can be modified as follows:
\begin{equation} \boldsymbol{\tau }=\boldsymbol{Y}(\boldsymbol{q},\dot{\boldsymbol{q}},\ddot{\boldsymbol{q}},\boldsymbol\alpha )\boldsymbol\pi, \end{equation}
where
$\boldsymbol{\pi }\in \Re ^{p \times 1}$
is the base parameter set and
$p$
is the number of base parameters.
$\boldsymbol{W}(\boldsymbol{q},\dot{\boldsymbol{q}},\ddot{\boldsymbol{q}},\boldsymbol\alpha )\in \Re ^{n \times 14n}$
is the rank-deficient regressor because of the linearly correlated columns.
$\boldsymbol{Y}(\boldsymbol{q},\dot{\boldsymbol{q}},\ddot{\boldsymbol{q}},\boldsymbol\alpha )\in \Re ^{n \times p}$
is with full rank after eliminating them. In order to simplify the expression, the regressor can be represented by the first letter. For example,
$\boldsymbol{W}(\boldsymbol{q},\dot{\boldsymbol{q}},\ddot{\boldsymbol{q}},\boldsymbol\alpha )$
and
$\boldsymbol{Y}(\boldsymbol{q},\dot{\boldsymbol{q}},\ddot{\boldsymbol{q}},\boldsymbol\alpha )$
can be represented as
$\boldsymbol{W}$
and
$\boldsymbol{Y}$
, respectively.
3. The proposed algorithm
3.1. Algorithm details
Assuming that the payload was rigidly connected to the robot flange, the dynamic equation of robot with payload can then be written by
\begin{equation} \boldsymbol{\tau }_b=\boldsymbol{Y}_{b}\boldsymbol{\pi }_b =\boldsymbol{Y}_b\boldsymbol{\pi }_a+\boldsymbol{W}_L\boldsymbol{\varphi }_L+\boldsymbol{W}_f\Delta \boldsymbol{\varphi }_f, \end{equation}
where subscript
$a$
means the case without payload and subscript
$b$
means the case with payload. Therefore,
$\boldsymbol{\pi }_{a}$
and
$\boldsymbol{\pi }_{b}$
mean the base parameters of the robot without payload and that with payload, respectively.
$\boldsymbol{Y}_b$
is the regressor matrix corresponding to
$\boldsymbol{\pi }_{b}$
.
$\boldsymbol{\tau }_b$
is the measured joint torques of the robot with payload.
$\boldsymbol{\varphi }_L\in \Re ^{10 \times 1}$
is the payload parameters, that is,
$\boldsymbol{\varphi }_L=[I_{xx,L},I_{xy,L},I_{xz,L},I_{yy,L},I_{yz,L},I_{zz,L},m_Lr_{x,L},m_Lr_{y,L},m_Lr_{z,L},m_L]^{\mathrm{T}}$
.
$m_L$
is the payload mass.
${r}_{x,L}$
,
${r}_{y,L}$
, and
${r}_{z,L}$
are the position of payload CoM on
$X$
axis,
$Y$
axis, and
$Z$
axis, respectively.
$I_{xx,L}, I_{yy,L},$
and
$I_{zz,L}$
are the moment of inertia of payload.
$I_{xy,L}, I_{xz,L},$
and
$I_{yz,L}$
are the products of inertia of payload. The corresponding payload regressor matrix is
$\boldsymbol{W}_L\in \Re ^{n \times 10}$
.
$\Delta \boldsymbol{\varphi }_f\in \Re ^{3n \times 1}$
is the linear friction variations due to the payload dynamics, and the corresponding regressor matrix is
$\boldsymbol{W}_f\in \Re ^{n \times 3n}$
.
The difference between the base parameters with payload and the ones without payload can be calculated by
\begin{equation} \boldsymbol{\varepsilon }=\boldsymbol{\pi }_{b}-\boldsymbol{\pi }_{a}, \end{equation}
where
$\boldsymbol{\varepsilon }\in \Re ^{p \times 1}$
is the parameter difference.
The symbolic expressions of
$\boldsymbol{\pi }_{a}$
and
$\boldsymbol{\pi }_{b}$
have been derived in the appendix. Then, the parameter difference
$\boldsymbol{\varepsilon }$
can be computed by (6) and the nonzero elements
$\boldsymbol{\varepsilon }_{nz}$
in
$\boldsymbol{\varepsilon }$
can be then extracted. As shown in the appendix, it can be found that there is a linear relationship between
$\boldsymbol{\varphi }_L$
and
$\boldsymbol{\varepsilon }_{nz}$
. Therefore, once the numerical results of
$\boldsymbol{\varepsilon }_{nz}$
were solved, the
$\boldsymbol{\varphi }_L$
could be then solved according to this relationship. In order to compute
$\boldsymbol{\varepsilon }$
and
$\boldsymbol{\varepsilon }_{nz}$
,
$\boldsymbol{\pi }_{a}$
and
$\boldsymbol{\pi }_{b}$
need to be identified. In fact,
$\boldsymbol{\pi }_{a}$
can be identified offline in advance by the least-squares method. Then, the remaining problem is how to solve
$\boldsymbol{\pi }_{b}$
online. In this paper, the RLS method [Reference Kubus, Kroger and Wahl10] is used to solve
$\boldsymbol{\pi }_{b}$
online. The specific implementation process of the proposed algorithm is shown below.
Firstly, the regressor matrix
$\boldsymbol{Y}_b$
can be computed by the commanded joint data at each sampling moment, namely
\begin{equation} \boldsymbol{Y}_{b}=\boldsymbol{Y}_{b}(\boldsymbol{q}_d,\dot{\boldsymbol{q}}_d,\ddot{\boldsymbol{q}}_d,\boldsymbol{\alpha }_a), \end{equation}
where
$\boldsymbol{q}_d,\dot{\boldsymbol{q}}_d,\ddot{\boldsymbol{q}}_d\in \Re ^{n \times 1}$
is the commanded joint position, commanded joint velocity, and commanded joint acceleration, respectively.
$\boldsymbol{\alpha }_a$
is the nonlinear friction parameters of the robot without payload, which was identified offline in advance with
$\boldsymbol{\pi }_{a}$
together.
Then, the gain matrix
$\boldsymbol{K}^{\ast }\in \Re ^{p \times n}$
at each sampling moment can be computed by
\begin{equation} \boldsymbol{K}^{\ast }=\boldsymbol{P}^{\ast }\boldsymbol{Y}_{b}^{\mathrm{T}}(\xi \textbf{1}+\boldsymbol{Y}_{b}\boldsymbol{P}^{\ast }\boldsymbol{Y}_{b}^{\mathrm{T}})^{-1}, \end{equation}
where
$\boldsymbol{P}^{\ast }\in \Re ^{p \times p}$
is the covariance matrix.
$\xi$
is the forgetting factor which is a positive scalar.
$\textbf{1}$
is a properly sized identity matrix.
The torque residual vector
$\boldsymbol{e}\in \Re ^{n \times 1}$
at each sampling moment can be computed by
\begin{equation} \boldsymbol{e}=\boldsymbol{\tau }_b-\boldsymbol{Y}_{b}\boldsymbol{\pi }_b. \end{equation}
Afterward, the base parameters
$\boldsymbol{\pi }_b$
can be updated by
\begin{equation} \boldsymbol{\pi }_b=\boldsymbol{\pi }_b+\boldsymbol{K}^{\ast }\boldsymbol{e}, \end{equation}
where the initial value of
$\boldsymbol{\pi }_b$
is
$\boldsymbol{\pi }_a$
.
The covariance matrix
$\boldsymbol{P}^{\ast }$
can be updated by
\begin{equation} \boldsymbol{P}^{\ast }=\boldsymbol{P}^{\ast }-\boldsymbol{K}^{\ast }\boldsymbol{Y}_{b}\boldsymbol{P}^{\ast }, \end{equation}
where the initial value of
$\boldsymbol{P}^{\ast }$
is an identity matrix.
The final solution of the iterative process is the base parameters of the robot with payload, namely,
$\boldsymbol{\pi }_b$
. Since the base parameters of the robot without payload
$\boldsymbol{\pi }_a$
can be obtained in advance through offline identification, the RLS iteration is mainly to solve
$\boldsymbol{\pi }_b$
. Once
$\boldsymbol{\pi }_b$
converges and is solved, the above iterative calculation process can end. Then, the payload parameters can be estimated by the derived symbolic relationship between the base parameter difference and payload parameters in the appendix. The numerical results of
$\boldsymbol{\varepsilon }$
can be solved by (6), and the nonzero elements
$\boldsymbol{\varepsilon }_{nz}$
can be extracted from
$\boldsymbol{\varepsilon }$
. According to the symbolic expressions of
$\boldsymbol{\varepsilon }_{nz}$
in the appendix, the following six payload parameters can be directly calculated by
\begin{equation} \begin{aligned} I_{xy,L}&=\varepsilon _{nz,17},\ I_{xz,L}=\varepsilon _{nz,18},\ I_{yz,L}=\varepsilon _{nz,19}, \\ I_{zz,L}&=\varepsilon _{nz,20},\ m_Lr_{x,L}=\varepsilon _{nz,21},\ m_Lr_{y,L}=\varepsilon _{nz,22}. \end{aligned} \end{equation}
where
$\varepsilon _{nz,17}$
,
$\varepsilon _{nz,18}$
,
$\varepsilon _{nz,19}$
,
$\varepsilon _{nz,20}$
,
$\varepsilon _{nz,21}$
, and
$\varepsilon _{nz,22}$
are the 17-
$th$
, 18-
$th$
, 19-
$th$
, 20-
$th$
, 21-
$th$
, and 22-
$th$
elements in
$\boldsymbol{\varepsilon }_{nz}$
, respectively. As explained above,
$\boldsymbol{\varepsilon }_{nz}$
is the non-zero element vector of
$\boldsymbol{\varepsilon }$
and its symbolic expressions have been shown in the appendix.
Then, the payload mass
$m_L$
can be calculated by the following relationship:
\begin{equation} \begin{bmatrix} a_3, & a_4, & d_5 \end{bmatrix} \cdot m_L= \begin{bmatrix} \varepsilon _{nz,5}, & \varepsilon _{nz,9}, & \varepsilon _{nz,12} \end{bmatrix}, \end{equation}
where
$a_3$
and
$a_4$
are the link lengths of link 2 and link 3, respectively.
$d_5$
is the joint distance between joint coordinate 4 and joint coordinate 5. They are all the kinematic parameters of the UR10 robot.
$\varepsilon _{nz,5}$
,
$\varepsilon _{nz,9}$
, and
$\varepsilon _{nz,12}$
are the 5-
$th$
, 9-
$th$
, and 12-
$th$
elements in
$\boldsymbol{\varepsilon }_{nz}$
, respectively.
When
$m_L$
was solved, the
$m_Lr_{z,L}$
can be solved by
\begin{equation} m_Lr_{z,L}=-d_6m_L-\varepsilon _{nz,15}, \end{equation}
where
$d_6$
is the joint distance between joint coordinate 5 and joint coordinate 6.
$\varepsilon _{nz,15}$
is 15-
$th$
element in
$\boldsymbol{\varepsilon }_{nz}$
.
Once
$m_Lr_{x,L}$
,
$m_Lr_{y,L}$
,
$m_Lr_{z,L}$
, and
$m_L$
were solved, the CoM (
$r_{x,L}$
,
$r_{y,L}$
,
$r_{z,L}$
) can be obtained by the ratios of them, such as
$r_{x,L}=\frac{m_Lr_{x,L}}{m_L}$
.
After that,
$I_{yy,L}$
can be solved by the linear equation as follows:
\begin{equation} \begin{bmatrix} 1, & 1 \end{bmatrix} \cdot I_{yy,L}= \begin{bmatrix} \varepsilon _{nz,13}-2d_6m_Lr_{z,L}-d_6^2m_L, & \varepsilon _{nz,14}-2d_6m_Lr_{z,L}-d_6^2m_L \end{bmatrix}, \end{equation}
where
$\varepsilon _{nz,13}$
and
$\varepsilon _{nz,14}$
are the 13-
$th$
and 14-
$th$
elements in
$\boldsymbol{\varepsilon }_{nz}$
, respectively.
Finally,
$I_{xx,L}$
can be solved by
\begin{equation} I_{xx,L}=\varepsilon _{nz,16}+I_{yy,L}, \end{equation}
where
$\varepsilon _{nz,16}$
is the 16-
$th$
element in
$\boldsymbol{\varepsilon }_{nz}$
.
3.2. The overall algorithm
In order to explain the proposed algorithm better, the pseudo-code of the proposed algorithm was shown below. The input of the algorithm include two parts. One part is the base parameters
$\boldsymbol{\pi }_a$
and the nonlinear friction parameters
$\boldsymbol{\alpha }_a$
of the robot without payload, which have been identified offline in advance. The other is the measured joint torques
$\boldsymbol{\tau }_b$
, the commanded joint position
$\boldsymbol{q}_d$
, the commanded joint velocity
$\dot{\boldsymbol{q}}_d$
, and the commanded joint acceleration
$\ddot{\boldsymbol{q}}_d$
of the robot with payload at each sampling moment. The output of the algorithm is the inertial parameters of the payload
$\boldsymbol{\varphi }_L$
. The proposed algorithm mainly includes three steps. The first step is initialization. The forgetting factor
$\xi$
, the covariance matrix
$\boldsymbol{P}^{\ast }$
, and the base parameters
$\boldsymbol{\pi }_b$
of the robot with payload should be initialized. The initial value of forgetting factor
$\xi$
is set to 100. The initial value of covariance matrix
$\boldsymbol{P}^{\ast }$
is set to 10 times
$p \times p$
identity matrix. The initial value of base parameters with payload
$\boldsymbol{\pi }_b$
is set to the one without payload
$\boldsymbol{\pi }_a$
. The second step is the RLS iteration to solve
$\boldsymbol{\pi }_b$
. When the robot tracks an excitation trajectory,
$\boldsymbol{\pi }_b$
can be computed online in real time by the RLS method (namely, the while loop) until it converges. If
$\boldsymbol{\pi }_b$
has not converged, the regression matrix
$\boldsymbol{Y}_b$
can be first computed according to (7). Then the gain matrix
$\boldsymbol{K}^{\ast }$
can be computed by (8). Afterward, the joint torque residual vector
$\boldsymbol{e}$
can be computed by (9). Finally, the base parameters
$\boldsymbol{\pi }_b$
and covariance matrix
$\boldsymbol{P}^{\ast }$
can be updated through (10) and (11), respectively. Once
$\boldsymbol{\pi }_b$
has reached the convergence condition, the iteration process ends. The convergence condition is that the consecutive times for
$\Delta \pi _{b,max} \lt 0.005$
are more than 300.
$\Delta{\pi _{b,max}}$
is the maximum absolute difference of the estimated inertial base parameters such as
$\pi _{b,31}\sim \pi _{b,36}$
at adjacent sampling moments. At the
$k$
-th sampling,
$\Delta{\pi _{b,max}}_k$
can be computed by
$\Delta{\pi _{b,max}}_k=\textrm{max}(\Delta{\pi _{b,31}}_k,\Delta{\pi _{b,32}}_k,\Delta{\pi _{b,33}}_k,\Delta{\pi _{b,34}}_k,\Delta{\pi _{b,35}}_k,\Delta{\pi _{b,36}}_k)$
.
$\rm max(\cdot )$
is a function to find the maximum value among the input. For example,
$\Delta{\pi _{b,31}}_k$
can be computed by
$\Delta{\pi _{b,31}}_k = \lvert{\pi _{b,31}}_k -{\pi _{b,31}}_{k-1} \rvert$
, where
${\pi _{b,31}}_k$
represents the estimated
$\pi _{b,31}$
at
$k$
-th sampling. The third step is to solve the payload parameters online by using the symbolic relationship between the parameter difference and the payload parameters. Once
$\boldsymbol{\pi }_b$
was solved, the numerical results of parameter difference
$\boldsymbol{\varepsilon }$
can be solved by (6). According to the derivation in the appendix, it can be found that there exists a linear relationship between parameter difference and the inertial parameters of the payload. Based on this algebraic relationship, the inertial parameters of the load can be quickly solved. After extracting the nonzero elements
$\boldsymbol{\varepsilon }_{nz}$
from the parameter difference, the payload parameters
$\boldsymbol{\varphi }_L$
can be calculated by (12)–(16). For example, according to the symbolic expressions of the nonzero elements of parameter difference in the appendix,
$I_{xy,L}, I_{xz,L}, I_{yz,L}, I_{zz,L}, m_Lr_{x,L},$
and
$m_Lr_{y,L}$
can be first solved by (12). The payload mass
$m_L$
can be then solved by (13). Afterward,
$m_Lr_{z,L}$
can be solved by (14). After that,
$I_{yy,L}$
can be solved by (15). Finally,
$I_{xx,L}$
can be solved by (16).
The proposed online payload identification algorithm
Actual acceleration and command acceleration.
The proposed algorithm considers the nonlinear friction coefficient
$\boldsymbol{\alpha }_a$
when calculating the regression matrix. In fact, the nonlinear friction model (2) is beneficial to improve the identification accuracy of
$\boldsymbol{\pi }_a$
[Reference Han, Wu, Liu and Xiong29, Reference Xu, Fan, Fang, Zhu and Zhao39]. In addition, the proposed algorithm uses commanded joint trajectory data (
$\boldsymbol{q}_d,\dot{\boldsymbol{q}}_d,\ddot{\boldsymbol{q}}_d$
) as computational input rather than the actual joint trajectory data. This is because the commanded joint acceleration signal has low noise. But the actual joint acceleration signals are obtained by performing two differentials on the actual joint position, which will introduce too much noise as shown in Fig. 1 and reduce the identification accuracy of the payload. In fact, the commanded trajectory signals can replace the actual ones as the calculation input because the commercial robots usually have a good trajectory tracking accuracy. In this article, the command acceleration signals can be acquired from the UR10 robot’s controller. The joint clearances and compliance of links might affect the accuracy of the proposed method when applying the proposed method to other robot platforms with significant joint flexibility (e.g., robots equipped with both harmonic drives and joint torque sensors) because they have effects on the accuracy of dynamic modeling and identification. However, in this article, due to the high stiffness of joints and links of UR10 robot, the effects of joint clearance and link compliance can be usually ignored in dynamic modeling [Reference Gaz, Magrini and De Luca50]. For example, in this article, the identified dynamics can be used to predict the measured torques accurately as shown in Fig. 2, and the root mean square errors (RMSE) are shown in Table II.
Joint torque prediction using identified dynamics of UR10 robot.
4. Experiments
4.1. Experiment platform
The platform used in this paper is UR10 robot. The kinematic frames of the robot is shown in Fig. 3 and the modified DH (MDH) parameters are listed in Table III. The UR10 robot is controlled by a C++ and Python interface
$ur\_rtde$
[Reference Lindvig51]. It can receive data and control the robot by using the real-time data exchange in the robot controller. In this article, the authors use the Python programming language to code the robot. The host computer that controls the robot is a laptop with a Linux operating system. Its CPU is i5-8250U-1.6 Ghz and the memory is 16 GB.
4.2. Data acquisition and data processing
The excitation trajectory for identification was based on Fourier series [Reference Swevers, Ganseman, Tukel, de Schutter and Van Brussel25]. For the
$j$
-th joint, the excitation trajectory equation of joint position
$q_{j}$
is as follows:
\begin{equation} q_{j}(t)=q_{j,0}+\sum _{l=1}^5 [\frac{a_{j,l}}{\omega _{f}l}\textrm{sin}(\omega _{f}lt)-\frac{b_{j,l}}{\omega _{f}l}\textrm{cos}(\omega _{f}lt)] \end{equation}
where both
$a_{j,l}$
and
$b_{j,l}$
are the excitation trajectory coefficients to be optimized, which are shown in Table IV.
$\frac{\omega _{f}}{2\pi }$
is the trajectory fundamental frequency, which is set to 0.1 Hz. The trajectory period is 10 s in this article.
$q_{j,0}$
is the offset constant of joint position, which is set to
$[0,-\pi/2,0-\pi/2,0,0]$
.
RMSE of torque errors between predicted torques and measured torques.
The kinematic frames of UR10 robot (
$n=6$
and
$p=58$
).
As mentioned in the introduction, incentive effects can be improved by optimizing the sub-matrix condition number of regressor [Reference Venture, Ayusawa and Nakamura27, Reference Bonnet, Fraisse, Crosnier, Gautier, González and Venture28]. Additionally, considering measurement noise covariance during optimization and conducting the process iteratively can better capture the dynamic characteristics of the robot [Reference Han, Wu, Liu and Xiong29]. Therefore, the objective function is to minimize the condition numbers of the regressor and its subregressors in this article.
\begin{equation} \begin{aligned} \min _{\boldsymbol{\gamma }} \quad \Vert 10 \cdot \textrm{cond}(\overline{\boldsymbol{Y}}^{\ast }(\boldsymbol{\gamma }))&+5 \cdot \textrm{cond}(\overline{\boldsymbol{Y}}_{i}^{\ast }(\boldsymbol{\gamma })) +\textrm{cond}(\overline{\boldsymbol{Y}}_{g}^{\ast }(\boldsymbol{\gamma }))+\textrm{cond}(\overline{\boldsymbol{Y}}_{f}^{\ast }(\boldsymbol{\gamma })) \Vert _2\\ \textrm{s.t.} \quad & \begin{bmatrix} \sum _{l=1}^5 \frac{-b_{j,l}}{\omega _{f}l}, & \sum _{l=1}^5 a_{j,l}, & \sum _{l=1}^5 \omega _{f}lb_{j,l} \end{bmatrix} =\textbf{0}\\ & \sum _{l=1}^5 \frac{\sqrt{a_{j,l}^2+b_{j,l}^2}}{\omega _{f}l} \le q_{j,max}\\ & \sum _{l=1}^5 \sqrt{a_{j,l}^2+b_{j,l}^2} \le \dot{q}_{j,max}\\ & \sum _{l=1}^5 \omega _{f}l\sqrt{a_{j,l}^2+b_{j,l}^2} \le \ddot{q}_{j,max} \end{aligned} \end{equation}
where
$\boldsymbol{\gamma }$
is the optimization variable.
$\overline{\boldsymbol{Y}}^{\ast }$
is the weighted regressor stacked by
$N$
joint trajectory data with
$\overline{\boldsymbol{Y}}^{\ast }=\boldsymbol{\Lambda }_0^{-\frac{1}{2}}\overline{\boldsymbol{Y}}$
.
$\overline{\boldsymbol{Y}}$
is the stacked version of
$\boldsymbol{Y}$
.
$\boldsymbol{\Lambda }_0$
is the initial value of initial covariance matrix, which can be computed by
$\boldsymbol{\Lambda }_0=\frac{\overline{\boldsymbol{R}}_0\cdot \overline{\boldsymbol{R}}_0^{\mathrm{T}}}{N-p}$
.
$\overline{\boldsymbol{R}}_0$
is the initial torque residual which can be computed by
$\overline{\boldsymbol{R}}_0=\overline{\boldsymbol{\tau }}-\overline{\boldsymbol{Y}}\boldsymbol{\pi }_{ols}$
.
$\overline{\boldsymbol{\tau }}$
is the stacked version of
$\boldsymbol{\tau }$
.
$\boldsymbol{\pi }_{ols}$
is the ordinary least squares solution of
$\boldsymbol{\pi }$
.
$\overline{\boldsymbol{Y}}_{i}^{\ast }$
,
$\overline{\boldsymbol{Y}}_{g}^{\ast }$
, and
$\overline{\boldsymbol{Y}}_{f}^{\ast }$
are the inertia tensor sub-regressor, first moment vector sub-regressor, and linear friction sub-regressor of
$\overline{\boldsymbol{Y}}^{\ast }$
, respectively.
$q_{j,max}$
,
$\dot{q}_{j,max}$
, and
$\ddot{q}_{j,max}$
are the limit of joint position, limit of joint velocity, and limit of joint acceleration, respectively.
$\rm cond(\cdot )$
is a function to solve the condition number of the input.
UR10 MDH parameters.
Trajectory parameters for online payload identification trajectories.
One can refer to refs. [Reference Swevers, Ganseman, Tukel, de Schutter and Van Brussel25, Reference Han, Wu, Liu and Xiong29] for more details about the excitation trajectory optimization. The optimized excitation trajectory with joint position, joint velocity, joint acceleration, and joint torque are shown in Fig. 4, which can give the condition number of the stacked regressor matrix 79.67. As described in ref. [Reference Hollerbach, Khalil and Gautier34], the dynamic characteristics of a robot can be well excited if the condition number of regressor is less than 100. Therefore, the optimized excitation trajectory in this article is well conditioned and adequate for exciting the whole dynamics of the robot. The rate of the real-time interface can be set to 125 Hz for CB-Series UR robot. For the online payload identification, no filtering is used for the sampled joint trajectory data.
Because the UR10 robot has no joint torque sensors, the joint torques can be calculated by
$\boldsymbol{\tau }=\boldsymbol{K}\boldsymbol{i}$
, in which
$\boldsymbol{K}$
is the joint drive gains and
$\boldsymbol{i}$
is the motor currents. The
$\boldsymbol{K}$
of UR10 robot has been identified by the method in [Reference Xu, Fan, Fang, Zhu and Zhao52]. One can refer to our previous work [Reference Xu, Fan, Fang, Zhu and Zhao52] for more details because the identification process of
$\boldsymbol{K}$
is not the main contribution in this paper. The identified results of joint drive gains are shown in Table V. In fact, the joint drive gains
$\boldsymbol{K}$
refer to the product of the motor torque constant
$\boldsymbol{k}$
, gear ratio
$\boldsymbol{\eta }$
, and current amplification factor
$\boldsymbol{\psi }$
, namely,
${\boldsymbol{K}} ={\boldsymbol{k}} \cdot{\boldsymbol{\eta }} \cdot{\boldsymbol{\psi }}$
[Reference Gautier and Briot53]. The identified joint drive gains directly establish the relationship between joint torques and motor currents; therefore, the identified results in Table V include the reduction ratio. For the offline identification of
$\boldsymbol{\pi }_a$
, the Butterworth filter was used for the estimated joint torques, and the frequency filtering was used for estimating the joint velocities and joint accelerations. The Butterworth filter order is set to 4, and the cutoff frequency is set to 40 Hz.
The optimized excitation trajectory for online payload identification.
4.3. Identification results
The
$\boldsymbol{\pi }_a$
and the
$\boldsymbol{\alpha }_a$
of UR10 robot without payload have been identified offline in ref. [Reference Xu, Fan, Fang, Zhu and Zhao39]. Since the identification process of both
$\boldsymbol{\pi }_a$
and
$\boldsymbol{\alpha }_a$
is not the main contribution in this paper, only the identification results of them are shown in Table VI, and the identification process is omitted. One can refer to our previous work [Reference Xu, Fan, Fang, Zhu and Zhao39] for more identification details. When the
$\boldsymbol{\pi }_a$
and
$\boldsymbol{\alpha }_a$
are ready, the online payload identification can start. Firstly, we make the UR10 robot track the designed trajectory with a concentric payload as shown in Fig. 5. The true values of the payload parameters are measured by the CAD software. Then, the proposed algorithm starts running at the same time. Since the identification process is performed online, no filtering is used to improve calculation efficiency. In addition, since the commanded trajectory signal has less noise, the commanded positions, the commanded velocities, and the commanded accelerations are used to compute
$\boldsymbol{Y}_b$
. After online iteration calculation by RLS method, the
$\boldsymbol{\pi }_b$
converged. In order to save space, the convergence process of part of
$\boldsymbol{\pi }_b$
is shown in Fig. 6. From the figure, it can be seen that the base parameters
$\pi _{b,31}\sim \pi _{b,36}$
converged at about 6.59 s because the consecutive times for
$\Delta \pi _{b,max} \lt 0.005$
have exceeded 300 at that moment. Once
$\boldsymbol{\pi }_b$
converged, the payload parameters can be estimated by using the symbolic relationship between the parameter difference and the payload parameters, namely, (12)–(16). Since the calculation of (12)–(16) is simple and very fast, therefore, the computation times for the online payload identification by using the proposed method are about 6.59 s.
Joint drive gains of UR10 robot.
$\boldsymbol{\pi }_a$
and
$\boldsymbol{\alpha }_a$
of UR10 robot by offline identification1.
1 All the units are SI unit.
Concentric payload for online identification.
Once
$\boldsymbol{\pi }_b$
is solved online, the numerical results of
$\boldsymbol{\varepsilon }$
can be computed by (6). Then, the nonzero elements
$\boldsymbol{\varepsilon }_{nz}$
can be extracted from
$\boldsymbol{\varepsilon }$
according to the appendix. Afterward, the payload parameters can be calculated online by (12)–(16). The identified results of
$\boldsymbol{\varphi }_L$
by the proposed algorithm are shown in Table VII. In addition, the identified results of
$\boldsymbol{\varphi }_L$
by another two classical algorithms are also shown in Table VII as comparison. The algorithm in ref. [Reference Khalil, Gautier and Lemoine33] is named as Algorithm 1 and the other one in ref. [Reference Gaz and De Luca36] is named as Algorithm 2. The principle of Algorithm 1 is to use the changes between the joint torques before and after loading the robot on the same trajectory. It distinguished the contribution of the robot itself and the payload to the joint torques. The robot link parameters without load were considered as the prior information for the online payload identification. However, the estimation results in this method are biased. The principle of Algorithm 2 is based on the relationship between robot static coefficients and inertial parameters of payload. However, this method cannot obtain the inertia tensor of the payload because it lacks the dynamic excitation.
Convergence process of the part of
$\boldsymbol{\pi }_b$
during online payload identification.
According to these data, it can be found that the identification error of the proposed algorithm is 0.04% for the payload mass, and the ones of Algorithm 1 and Algorithm 2 are 4.11% and 0.92%, respectively. The proposed algorithm can improve the online identification accuracy of the payload mass by about 103 times and 23 times compared with Algorithm 1 and Algorithm 2, respectively. Therefore, the identification accuracy of the proposed algorithm is higher than that of Algorithm 1 and Algorithm 2 for the payload mass. For the CoM in
$Y$
axis of the payload, the proposed algorithm has the highest identification accuracy. For the CoM in
$Z$
axis of the payload, the identification error of the proposed algorithm (74.32%) is smaller than Algorithm 1 (94.93%) and higher than Algorithm 2 (51.69%). It can be found that all the three algorithms cannot accurately identify the CoM of the payload online. This is because the online payload identification uses a small amount of data to improve the computation efficiency. Therefore, it is difficult to obtain the same identification accuracy of CoM with that of the offline payload identification. For the inertia tensor of payload, both the proposed algorithm and Algorithm 1 are not accurate enough. Usually, it is not easy to accurately identify the inertia tensor of the payload neither online payload identification nor offline payload identification, especially when the payload mass is lighter than the robot mass [Reference Mavrakis and Stolkin2, Reference Xu, Fan, Fang, Zhu and Zhao39].
Online identification results of the payload by three methods.
1
The units of mass, CoM, and inertia tensor are
$\textrm{kg}$
,
$\textrm{m}$
, and
$\textrm{kg}\cdot \textrm{m}^2$
, respectively.
2
The true values were measured by the CAD software, which are relative to the frame
$x_6y_6z_6$
in Fig. 3.
According to the experimental results, it can be found that the proposed algorithm has the highest online identification accuracy for the payload mass. The CoM online identification accuracy of the proposed algorithm is better than Algorithm 1. Although the CoM online identification accuracy of Algorithm 2 is higher than the proposed algorithm, the Algorithm 2 cannot identify the inertia tensor of the payload online. Consistent with existing algorithms, both the proposed algorithm and Algorithm 1 have deficiencies in the online identification accuracy of payload inertia tensor.
As mentioned above, the nonlinear friction model is beneficial to improve the identification accuracy of
$\boldsymbol{\pi }_a$
. And then
$\boldsymbol{\pi }_a$
affects the online identification accuracy of payload. As shown in Fig. 7, the nonlinear friction model fits the estimated friction torque better than the linear one in both unloaded and loaded conditions. The used linear friction model (green circle) is
$\tau _{f,j}=F_{c,j}\textrm{sign}(\dot{q}_j)+F_{v,j} \dot{q}_j + B_j$
[Reference Zhang and Zhang44]. The parameters of this model for joint 2 is
$F_{c,2}=15.07$
,
$F_{v,2}=10.73$
, and
$B_2=8.77\times 10^{-3}$
at the case without payload. The parameters of this model for joint 2 is
$F_{c,2}=15.87$
,
$F_{v,2}=10.77$
, and
$B_2=6.47\times 10^{-8}$
at the case with payload. The used nonlinear friction model (red dot) is
$\tau _{f,j}=(F_{c,j}+F_{v,j}\lvert \dot{q}_j \rvert ^{\alpha _{j}})\textrm{sign}(\dot{q}_j)+B_j$
, namely, (2) [Reference Han, Wu, Liu and Xiong29, Reference Farhat, Mata, Page and Valero45]. The parameters of this model for joint 2 is
$\alpha _2=0.2723$
,
$F_{c,2}=-0.1240$
,
$F_{v,2}=26.4096$
, and
$B_2=0.0582$
at the case without payload. The parameters of this model for joint 2 is
$\alpha _2=0.2603$
,
$F_{c,2}=-0.5500$
,
$F_{v,2}=27.7899$
, and
$B_2=-0.1836$
at the case with payload. Especially when the joint velocity is close to zero, the linear friction model will introduce more fitting errors. As shown in Table VIII, it can be found that the RMSE of the nonlinear friction model is smaller than that of the linear one. Therefore, it benefits to improving the accuracy of online payload identification by using a nonlinear friction model. In addition, it is worth explaining that
$\boldsymbol{\alpha }_a$
was used as the input for online payload identification. On the one hand, it is difficult to identify the nonlinear friction coefficient
$\boldsymbol{\alpha }_b$
with payload in each iteration online. On the other hand, the payload has little effect on the nonlinear friction coefficients except for joint 3 as shown in Table IX. In order to compute the estimated friction torques, the identified friction parameters in
$\boldsymbol{\pi }$
are firstly set to zero, then the torques related to the inertial component in
$\boldsymbol{\pi }$
can be subtracted from the measured torques [Reference Li, Wei, Liu, He, Liu, Zhang and Li54]. For example, the estimated friction torques
$\boldsymbol{\hat{\tau }_f} = \boldsymbol{\tau } - \boldsymbol{Y} \cdot \left [{\boldsymbol{\hat{\pi }}_i}^{\mathrm{T}}, \textbf{0}^{\mathrm{T}} \right ]^{\mathrm{T}}$
, in which
$\boldsymbol{\hat{\tau }_f}$
is the estimated friction torques,
$\boldsymbol{\tau }$
is the measured joint torques,
$\boldsymbol{\hat{\pi }}_i$
is the identified inertial parameter component in
$\boldsymbol{\pi }$
.
Comparison between linear friction and nonlinear friction in joint 2 of UR10 robot without and with payload.
The temperature range of the robot joints at which the experiments are performed is shown in Fig. 8. The authors make the robot to execute the online identification excitation trajectory 10 times and collect the dataset from the robot controller. There was an approximate 10 s interval between each movement. Typically, before the online payload identification experiment, the authors first let the robot run 3 times to ensure the robot’s status and the code being functioning properly. Therefore, the 10 datasets could fully cover the temperature range of the experiment at that time. From the data in the graph, it can be observed that the mean temperature range of all joints is between 25.5
$^{\circ }$
C and 30
$^{\circ }$
C.
4.4. An application of online payload identification
In this subsection, an application experiment of online payload identification is shown. Firstly, one uses the proposed algorithm to identify the inertial parameters of a new payload online. Then, the identified inertia parameters of the payload are compensated to the robot controller. Afterward, the teaching mode of UR10 robot can be activated. In this mode, the user can drag the robot to move freely, and the guidance forces can be recorded by a force/torque sensor. The force/torque sensor used in this work is the ATI MINI45-E. The sensor has a diameter of 45 mm, a height of 15.7 mm, and a mass of 0.0917 kg. The maximum force range is
$\pm$
580 N in X-axis and Y-axis,
$\pm$
1160 N in Z-axis. The maximum torque range is
$\pm$
20 Nm in X-axis, Y-axis, and Z-axis. The force resolution of the sensor is 1/4 N. The torque resolution of the sensor is 1/188 Nm in X-axis and Y-axis and 1/376 Nm in Z-axis. Since the force/torque sensor has a non-detachable connecting wire, it is inconvenient to identify it together with the payload. Therefore, the payload is first identified online to prevent the connecting wire from being torn off when the robot executes the excitation trajectory. Then, the force/torque sensor and the guidance handle are fixed to the end of the payload. After that, the manual guidance experiment can be carried out. In this application experiment, the designed payload and the assembly scheme are shown in Fig. 9. The mass of both new payload and the bolts is 3.223 kg after measurement. By using the proposed algorithm, the inertial parameters of the payload are identified online. The identified mass is 3.237 kg and the error is 0.43%. The identified CoM of the payload is [−0.024, −0.006, 0.045] m. The identified inertia tensor of the payload is [0.006, −0.012, −0.069, −0.106, 0.031, −0.014]
$\textrm{kg}\cdot \textrm{m}^2$
. Since the UR10 robot does not provide a compensation interface for the inertia tensor of payload, the user can only compensate the mass and the CoM of the payload. In fact, the mass of the payload fixed at the end of the robot is usually much smaller than the mass of the robot itself. On the one hand, this will lead to inaccurate estimation of the inertia tensor of the payload. On the other hand, compared with the inertia tensor of the robot itself, the inertia tensor of the payload has little contribution on the dynamics of the robot. It can be speculated that UR robot manufacturer may not provide the compensation interface for payload inertia tensor due to the above reasons.
RMSE1 of different friction models without and with payload.
1
The unit of RMSE is
$\textrm{Nm}$
.
The nonlinear friction coefficients
$\boldsymbol{\alpha }$
of UR10 robot without and with payload.
The mean temperature of each joint of UR10 robot during excitation motions.
New payload for the application experiment. (a) Payload. (b) Connection diagram.
The rectangular guidance path for the UR10 robot.
Measured TCP forces/torques and motor currents of UR10 robot without and with payload compensation during the rectangular guidance path. (a) Measured joint currents without payload compensation. (b) Measured forces/torques without payload compensation. (c) Measured joint currents with payload compensation. (d) Measured forces/torques with payload compensation.
Measured TCP force in manual guidance application experiments along rectangular path.
After completing the online identification of the payload, the force/torque sensor and the guidance handle are fixed to the end of the payload. Then, the user can drag the robot to move along the rectangle path as shown in Fig. 10. Figure 10(b) depicts the position of the robot’s Tool Center Point (TCP) on the
$XY$
plane of the experimental platform as shown in Fig. 10(a) when the author manually dragged the robot end effector. The TCP position of robot is sampled from the robot controller. At the same time, the guidance forces/torques and motor currents of robot are recorded as shown in Fig. 11. The mean and variance of the guidance forces without/with payload compensation are plotted in Fig. 12. From the figure, it can be seen that, compared with the ones without payload compensation, the mean of guidance forces with payload compensation can be reduced by 12.99%, 25.75%, and 44.24% in
$X$
-,
$Y$
-, and
$Z$
-axis, respectively. Similarly, after compensating for the payload, the variance of guidance forces can be reduced by 51.84%, 56.83%, and 64.58% in
$X$
-,
$Y$
-, and
$Z$
-axis, respectively. The variance of force in Fig. 12(b) is defined by
${\sigma ^2} = \frac{1}{M}\sum \limits _{i = 1}^M{{{\left ({{x_i} - \mu } \right )}^2}}$
, in which
$M$
is the data number,
$x_i$
is the measured force vector, and
$\mu$
is the mean value of measured force vector. It can be found that the guidance forces (especially in
$Z$
-axis) are significantly reduced with payload compensation. It is because additional forces are needed to balance the payload gravity if the payload parameters are not compensated. However, the robot controller can increase currents to counteract the payload gravity once the payload parameters are compensated. Therefore, the user will feel that the robot becomes light during the guidance process. The experimental results show that the proposed algorithm can identify accurate payload parameters online (especially the payload mass), and the guidance comfort can be improved by compensating the identified payload parameters.
Based on our experience, it is not easy to drag the UR10 robot in teaching mode to make its end-effector track the rectangle in Fig. 10(b). This difficulty is unrelated to compensating for the dynamic of the load. One possible reason is that the UR10 robot lacks joint torque sensors, which makes it less smooth and easy to drag its end effector to track the rectangle in Fig. 10(b). In fact, whether compensating for load dynamics only affects the amount of force applied by the user during the dragging process. For instance, when load dynamics are not compensated for, more force is required for dragging, whereas when load dynamics are compensated for, less force is needed, as indicated by the data in Fig. 12. It is worth noting that the force/torque sensor used in this work was solely employed for quantitative analysis of dragging forces and was not utilized for any algorithm design and integration. In the teaching mode of the UR10 robot, the force/torque sensor was merely employed to record the forces applied by the user while dragging the robot’s end effector.
5. Conclusions
In this paper, a symbolic relationship between the parameter difference and the payload parameters was derived. Based on the symbolic relationship, an online method is proposed to identify the payload parameters. In order to increase the online estimation accuracy, on the one hand, a nonlinear friction model (2) was used instead of the linear one. On the other hand, the commanded joint trajectory signals were introduced as the computation input instead of the actual ones. Compared with other methods, the online estimation accuracy of the payload mass
$m_L$
can be significantly improved by the proposed method. For example, the accuracy of payload mass estimated by the proposed method can be improved by 103 times and 23 times with respect to Algorithm 1 [Reference Khalil, Gautier and Lemoine33] and Algorithm 2 [Reference Gaz and De Luca36], respectively. According to the experimental results, the computation times by using the proposed method are about 6.59 s. When the identified dynamics of the payload were compensated, the manual guidance forces from the user can be significantly reduced and the guidance comfort can be improved. In this paper, the mean of guidance forces in
$Z$
axis was reduced by 44.24% compared with the case without payload dynamics compensation.
How to accurately identify the inertia tensor of a payload online by using proprioceptive sensors only is still a challenge that needs to be solved. One potential avenue for the future research is to reduce the computational complexity of the online identification model. More joint data can be included for online calculation, and then complex online filtering processing can be performed on the joint data. In the future, the authors will continue to develop new online payload estimation methods and further verify the new methods and the proposed one in this article on a domestic robot.
Acknowledgments
The authors would like to thank the editors and reviewers for their professional comments on the improvement of this article.
Author contributions
Conceptualization, T.X. and H.T.; methodology, T.X. and J.C.; formal analysis, T.X. and H.T.; validation, T.X. and D.S.; data curation, T.X. and Q.F.; writing–original draft preparation, T.X.; writing–review and editing, T.X. and J.F.; supervision, J.Z.; project administration, J.F. and J.Z.
Financial support
This work was supported by the National Key R&D Program of China under Grant No. 2022YFB4701100, 2022YFB4700300, the National Natural Science Foundation of China under Grant No. 92048301, and the Pre-research Task of State Key Laboratory of Robotics and Systems (HIT) under Grant No. SKLRS202411B.
Competing interests
The authors have no relevant financial or non-financial interests to disclose.
Appendix
The symbolic expressions of
$\boldsymbol{\pi }_a$
of UR10 robot without payload are shown as follows:
\begin{equation*} \pi _{a,1}=I_{a,1}+I_{zz,1}+\sum _{j=2}^4 I_{yy,j}+2d_{4}m_{4}r_{z,4}+a_{3}^2m_3+(a_{3}^2+a_{4}^2+d_{4}^2)\sum _{j=4}^6 m_j, \end{equation*}
\begin{equation*} \pi _{a,2}=I_{xx,2}-I_{yy,2}-a_{3}^2\sum _{j=3}^6 m_j,\ \pi _{a,3}=I_{xy,2}, \end{equation*}
\begin{equation*} \pi _{a,4}=I_{xz,2}-a_{3}[m_{3}r_{z,3}+m_{4}r_{z,4}+d_{4}\sum _{j=4}^6 m_j],\ \pi _{a,5}=I_{yz,2}, \end{equation*}
\begin{equation*} \pi _{a,6}=I_{a,2}+I_{zz,2}+a_{3}^2\sum _{j=3}^6 m_j,\ \pi _{a,7}=m_{2}r_{x,2}+a_{3}\sum _{j=3}^6 m_j,\ \pi _{a,8}=m_{2}r_{y,2}, \end{equation*}
\begin{equation*} \pi _{a,9}=I_{xx,3}-I_{yy,3}-a_{4}^2\sum _{j=4}^6 m_j,\ \pi _{a,10}=I_{xy,3}, \end{equation*}
\begin{equation*} \pi _{a,11}=I_{xz,3}-a_{4}m_{4}r_{z,4}-a_{4}d_{4}\sum _{j=4}^6 m_j,\ \pi _{a,12}=I_{yz,3}, \end{equation*}
\begin{equation*} \pi _{a,13}=I_{zz,3}+a_{4}^2\sum _{j=4}^6 m_j,\ \pi _{a,14}=m_{3}r_{x,3}+a_{4}\sum _{j=4}^6 m_j,\ \pi _{a,15}=m_{3}r_{y,3}, \end{equation*}
\begin{equation*} \pi _{a,16}=I_{xx,4}-I_{yy,4}+I_{yy,5}+2d_{5}m_{5}r_{z,5}+d_{5}^2\sum _{j=5}^6 m_j,\ \pi _{a,17}=I_{xy,4}, \end{equation*}
\begin{equation*} \pi _{a,18}=I_{xz,4},\ \pi _{a,19}=I_{yz,4},\ \pi _{a,20}=I_{zz,4}+I_{yy,5}+2d_{5}m_{5}r_{z,5}+d_{5}^2\sum _{j=5}^6 m_j, \end{equation*}
\begin{equation*} \pi _{a,21}=m_{4}r_{x,4},\ \pi _{a,22}=m_{4}r_{y,4}+m_{5}r_{z,5}+d_{5}\sum _{j=5}^6 m_j, \end{equation*}
\begin{equation*} \pi _{a,23}=I_{xx,5}-I_{yy,5}+I_{yy,6}+2d_{6}m_{6}r_{z,6}+d_{6}^2m_6,\ \pi _{a,24}=I_{xy,5}, \end{equation*}
\begin{equation*} \pi _{a,25}=I_{xz,5},\ \pi _{a,26}=I_{yz,5},\ \pi _{a,27}=I_{zz,5}+I_{yy,6}+2d_{6}m_{6}r_{z,6}+d_{6}^2m_{6}, \end{equation*}
\begin{equation*} \pi _{a,28}=m_{5}r_{x,5},\ \pi _{a,29}=m_{5}r_{y,5}-m_{6}r_{z,6}-d_{6}m_{6},\ \pi _{a,30}=I_{xx,6}-I_{yy,6}, \end{equation*}
\begin{equation*} \pi _{a,31}=I_{xy,6},\ \pi _{a,32}=I_{xz,6},\ \pi _{a,33}=I_{yz,6},\ \pi _{a,34}=I_{zz,6},\ \pi _{a,35}=m_{6}r_{x,6}, \end{equation*}
\begin{equation*} \pi _{a,36}=m_{6}r_{y,6},\ \pi _{a,37}=I_{a,3},\ \pi _{a,38}=I_{a,4},\ \pi _{a,39}=I_{a,5},\ \pi _{a,40}=I_{a,6}, \end{equation*}
\begin{equation*} \pi _{a,41}=F_{c,1},\ \pi _{a,42}=F_{c,2}, \ \pi _{a,43}=F_{c,3},\ \pi _{a,44}=F_{c,4},\ \pi _{a,45}=F_{c,5}, \end{equation*}
\begin{equation*} \pi _{a,46}=F_{c,6},\ \pi _{a,47}=F_{v,1},\ \pi _{a,48}=F_{v,2},\ \pi _{a,49}=F_{v,3},\ \pi _{a,50}=F_{v,4}, \end{equation*}
\begin{equation*} \pi _{a,51}=F_{v,5},\ \pi _{a,52}=F_{v,6},\ \pi _{a,53}=B_{1},\ \pi _{a,54}=B_{2},\ \pi _{a,55}=B_{3}, \end{equation*}
\begin{equation*} \pi _{a,56}=B_{4},\ \pi _{a,57}=B_{5},\ \pi _{a,58}=B_{6}. \end{equation*}
Assuming that the payload is rigidly connected to robot flange, the payload parameters are then regarded as a part of the inertia parameters of the link
$n$
. Therefore, the inertial parameters of link
$n$
with considering the payload can be expressed as follows:
\begin{equation*} m'_{\!\!n}=m_{n}+m_L,\; m_{n}\boldsymbol{r}'_{\!\!n}=m_{n}\boldsymbol{r}_{n}+m_L\boldsymbol{r}_{L}, \; \boldsymbol{l}'_{\!\!n}=\boldsymbol{l}_{n}+\boldsymbol{l}_{L}, \end{equation*}
where
$\boldsymbol{l}_{n},m_{n}\boldsymbol{r}_{n}$
, and
$m_{n}$
are the inertial parameters of link
$n$
without payload.
$\boldsymbol{l}_{L},m_{L}\boldsymbol{r}_{L}$
, and
$m_{L}$
are the payload parameters.
$\boldsymbol{l}'_{\!\!n},m_{n}\boldsymbol{r}'_{\!\!n}$
, and
$m'_{\!\!n}$
are the inertial parameters of link
$n$
with payload. These parameters are all expressed relative to the coordinate
$x_ny_nz_n$
.
Therefore, when the payload is fixed at the end of the UR10 robot,
$m_6$
can be substituted by
$m'_{\!\!6}=m_6+m_L$
,
$m_6r_{x,6}$
can be substituted by
$m_6r'_{\!\!x,6}=m_6r_{x,6}+m_Lr_{x,L}$
,
$m_6r_{y,6}$
can be substituted by
$m_6r'_{\!\!y,6}=m_6r_{y,6}+m_Lr_{y,L}$
,
$m_6r_{z,6}$
can be substituted by
$m_6r'_{\!\!z,6}=m_6r_{z,6}+m_Lr_{z,L}$
,
$I_{xx,6}$
can be substituted by
$I'_{\!\!xx,6}=I_{xx,6}+I_{xx,L}$
,
$I_{xy,6}$
can be substituted by
$I'_{\!\!xy,6}=I_{xy,6}+I_{xy,L}$
,
$I_{xz,6}$
can be substituted by
$I'_{\!\!xz,6}=I_{xz,6}+I_{xz,L}$
,
$I_{yy,6}$
can be substituted by
$I'_{\!\!yy,6}=I_{yy,6}+I_{yy,L}$
,
$I_{yz,6}$
can be substituted by
$I'_{\!\!yz,6}=I_{yz,6}+I_{yz,L}$
, and
$I_{zz,6}$
can be substituted by
$I'_{\!\!zz,6}=I_{zz,6}+I_{zz,L}$
. Moreover, the linear friction variations caused by the payload need to be also considered. Then, the symbolic expressions of
$\boldsymbol{\pi }_b$
of UR10 robot with payload can be obtained. In order to save space, the equations of
$\boldsymbol{\pi }_b$
are omitted. According to the symbolic expressions of both
$\boldsymbol{\pi }_a$
and
$\boldsymbol{\pi }_b$
, the difference between them can be easily computed by (6). Then, the nonzero elements
$\boldsymbol{\varepsilon }_{nz}$
in
$\boldsymbol{\varepsilon }$
can be extracted as follows:
\begin{equation*} \varepsilon _{nz,1}=(d_4^2+a_4^2+a_3^2)m_L,\ \varepsilon _{nz,2}=-a_3^2m_L,\ \varepsilon _{nz,3}=-a_3d_4m_L,\ \varepsilon _{nz,4}=a_3^2m_L, \end{equation*}
\begin{equation*} \varepsilon _{nz,5}=a_3m_L,\ \varepsilon _{nz,6}=-a_4^2m_L,\ \varepsilon _{nz,7}=-a_4d_4m_L,\ \varepsilon _{nz,8}=a_4^2m_L, \end{equation*}
\begin{equation*} \varepsilon _{nz,9}=a_4m_L,\ \varepsilon _{nz,10}=d_5^2m_L,\ \varepsilon _{nz,11}=d_5^2m_L,\ \varepsilon _{nz,12}=d_5m_L, \end{equation*}
\begin{equation*} \varepsilon _{nz,13}=I_{yy,L}+2d_6m_Lr_{z,L}+d_6^2m_L,\ \varepsilon _{nz,14}=I_{yy,L}+2d_6m_Lr_{z,L}+d_6^2m_L, \end{equation*}
\begin{equation*} \varepsilon _{nz,15}=-m_Lr_{z,L}-d_6m_L,\ \varepsilon _{nz,16}=I_{xx,L}-I_{yy,L},\ \varepsilon _{nz,17}=I_{xy,L}, \end{equation*}
\begin{equation*} \varepsilon _{nz,18}=I_{xz,L},\ \varepsilon _{nz,19}=I_{yz,L},\ \varepsilon _{nz,20}=I_{zz,L}, \ \varepsilon _{nz,21}=m_Lr_{x,L},\ \varepsilon _{nz,22}=m_Lr_{y,L}, \end{equation*}
\begin{equation*} \varepsilon _{nz,23}=\Delta F_{c,1},\ \varepsilon _{nz,24}=\Delta F_{c,2},\ \varepsilon _{nz,25}=\Delta F_{c,3},\ \varepsilon _{nz,26}=\Delta F_{c,4},\ \varepsilon _{nz,27}=\Delta F_{c,5}, \end{equation*}
\begin{equation*} \varepsilon _{nz,28}=\Delta F_{c,6},\ \varepsilon _{nz,29}=\Delta F_{v,1},\ \varepsilon _{nz,30}=\Delta F_{v,2},\ \varepsilon _{nz,31}=\Delta F_{v,3},\ \varepsilon _{nz,32}=\Delta F_{v,4}, \end{equation*}
\begin{equation*} \varepsilon _{nz,33}=\Delta F_{v,5},\ \varepsilon _{nz,34}=\Delta F_{v,6},\ \varepsilon _{nz,35}=\Delta B_{1},\ \varepsilon _{nz,36}=\Delta B_{2},\ \varepsilon _{nz,37}=\Delta B_{3}, \end{equation*}
\begin{equation*} \varepsilon _{nz,38}=\Delta B_{4},\ \varepsilon _{nz,39}=\Delta B_{5},\ \varepsilon _{nz,40}=\Delta B_{6}. \end{equation*}
In above,
$a_j$
and
$d_j$
are the MDH parameters in Table III. From the above derivation, it can be seen that
$\boldsymbol{\varepsilon }_{nz}$
is made up of
$\varepsilon _{1}$
,
$\varepsilon _{2}$
,
$\varepsilon _{4}$
,
$\varepsilon _{6}$
,
$\varepsilon _{7}$
,
$\varepsilon _{9}$
,
$\varepsilon _{11}$
,
$\varepsilon _{13}$
,
$\varepsilon _{14}$
,
$\varepsilon _{16}$
,
$\varepsilon _{20}$
,
$\varepsilon _{22}$
,
$\varepsilon _{23}$
,
$\varepsilon _{27}$
,
$\varepsilon _{29}\sim \varepsilon _{36}$
, and
$\varepsilon _{41}\sim \varepsilon _{58}$
.
