1. Introduction
With the improvement of the ratio from the payload mass to the manipulator mass, the robot payload brings about a much larger contribution to the actuator torques and has an increasing importance in robotics [Reference Hirzinger, Sporer, Albu-Schaffer, Hahnle, Krenn, Pascucci and Schedl1]. For example, the inertial parameters of the payload can provide a strong priori for planning and controlling of robots [Reference Mavrakis and Stolkin2–Reference Li, Wang, Wang, Li, Zhao and Ma4]. In addition, inaccurate payload inertial parameters can lead to approximately 50
$\%$
of robot overloading in path planning [Reference Swevers, Verdonck, Naumer, Pieters and Biber5]. Therefore, the accurate knowledge of the payload inertial parameters becomes increasingly important. In this article, only this case is considered, in which a payload was rigidly connected to the end effector of the robot. For this case, there are mainly two experimental identification approaches. One experimental method is based on the force–torque sensor mounted at the robot wrist [Reference Duan, Liu, Bin, Cui and Dai6]. By collecting the measured forces/torques and the measured joint angles, the inertial parameters of the payload can be solved by some linear equations. In this regard, some classical works can be found in refs. [Reference Atkeson, An and Hollerbach7–Reference Schedlinski and Link9]. However, this approach is not generally applicable because most industrial robots are not equipped with these sensors.
Since the payload is rigidly attached and can be considered as a part of the last robot link, the other experimental method considers the payload identification as a regular dynamic identification problem. In this regard, Swevers et al. [Reference Swevers, Verdonck, Naumer, Pieters and Biber5] proposed a method based on the periodic robot excitation and the maximum likelihood estimation. In this method, the mass and the inertia tensor of the last robot link was negligible. The payload was then assumed as the last link, which would obviously introduce estimation errors. To solve this problem, a classical payload identification model was derived in ref. [Reference Swevers, Verdonck and De Schutter10]. It was based on the observation that the robot links and the payload each contributed separately to the joint torques. The base parameter set of the robot without payload can be estimated in advance. Then, it can be considered as a priori information for the payload identification. However, only the ordinary least squares (OLS) solution is obtained in this method. Then, Khalil et al. presented two impressive methods in ref. [Reference Khalil, Gautier and Lemoine11] to deal with this issue. One is to compute the differences of joint torques with the same trajectories twice. The first trajectory was with payload and the second one was without payload. But this method assumes that the difference of the joint variables is negligible. The other one is to estimate robot link parameters and payload parameters together in one step based on the weighted least squares (WLS) method. Later, these methods were summarized in ref. [Reference Hollerbach, Khalil and Gautier12]. After that, Gaz et al. [Reference Gaz and De Luca13] proposed a different method. It was based on the linear relationship between the differences of the robot static coefficients and the inertial parameters of payload. However, this method can only identify the mass and the center of mass (CoM) of the payload. Recently, some researchers [Reference Bahloul, Tliba and Chitour14,Reference Farsoni and Bonfè15] have also tried to apply the nonlinear optimization techniques with physical constraints and the intelligent algorithms [Reference Li, Gu, Sun, Li and Tang16] to deal with the problem of payload estimation. But the estimation accuracy of the payload parameters in these methods is unsatisfactory. In addition, some online payload estimation methods [Reference Ding, Zhou, Li and Rong17] have been proposed, such as the recursive least squares (RLS) methods [Reference Kubus, Kroger and Wahl18–Reference Sánchez, Torres-Torriti and Cheein20], the Kalman filter methods [Reference Farsoni, Landi, Ferraguti, Secchi and Bonfe21,Reference Farsoni, Landi, Ferraguti, Secchi and Bonfè22], the adaptive algorithms [Reference Hu, Li, Chen and Yao23], and the algebraic methods [Reference Rodriguez, Ibanez and Battle24]. However, the identification accuracy of these online methods is usually worse than that of off-line methods. This is because the sampling data amounts are reduced to improve online computation efficiency. In addition, noise cannot be dealt with better by the online filtering than that by the off-line filtering. Undoubtedly, these excellent works mentioned above provide the good theoretical basis and the valuable reference for the payload identification. However, the data weighting technique, which is to put additional weighting matrices on data points so that the outliers contribute less to the final results and data points after weighting can fit the theoretical model better, has not been applied to the payload identification yet. This technique was originally proposed in ref. [Reference Holland and Welsch25]. Then, it was previously introduced in robot dynamic identification by Janot et al. [Reference Janot, Vandanjon and Gautier26] and Han et al. [Reference Han, Wu, Liu and Xiong27]. In this article, the data weighting technique was explored and used to further improve the identification accuracy of payload parameters.
In this article, the main contribution is to propose an accurate identification framework based on double weighting for inertial parameters of robot payloads. (1) In order to obtain the weighting matrices, a modified dynamic parameter identification method with two loops is firstly proposed. In the loops, the data weighting technique is applied which can reduce the effects of outliers on the identification results. When the iterations finished, the converged covariance matrices and the converged data weighting matrices can be obtained. (2) Based on the identified results of dynamic parameters, a payload identification model based on double weighting is then constructed. The first weighting is to use the covariance matrices to normalize the constructed model. It helps to obtain the WLS solution. The second weighting is to use the data weighting matrices to make the data points fit the theoretical model well. It helps to improve the identification results. (3) Experimental comparisons between our method and another four methods are conducted. The results confirm that our method shows the best performance, especially on improving the identification accuracy of the mass and the CoM of the payload.
The remainder of this article is organized as follows. Section 2 presents the preliminaries. Section 3 introduces the proposed method which includes the identification of both dynamic parameters and payload parameters of robots. Section 4 provides the experimental results and comparative results between the proposed method and another four methods. Finally, Section 5 concludes this article.
2. Preliminaries
By considering the inverse dynamic model on motor sides and mapping it to the link sides, the differential equations of motion, which describe the inverse dynamic model of a serial robot with
$n$
rigid bodies, can be written in the form:
\begin{equation} \boldsymbol{\tau }=(\boldsymbol{{M(q)}}+\boldsymbol{{I}}_{a})\ddot{\boldsymbol{{q}}}+{\boldsymbol{{C}}(\boldsymbol{{q}},\dot{\boldsymbol{{q}}})\dot{\boldsymbol{{q}}}}+\boldsymbol{{g(q)}}+\boldsymbol{\tau }_{f}, \end{equation}
where
$\boldsymbol{{q}},\dot{\boldsymbol{{q}}},\ddot{\boldsymbol{{q}}}\in \Re ^{n \times 1}$
, with which
$n$
is the degrees of freedom (DoF) of the robot, are the vectors of joint positions, joint velocities, and joint accelerations, respectively.
$\boldsymbol{{M(q)}} \in \Re ^{n \times n}$
is the symmetric and positive-definite inertia matrix of the links.
$\boldsymbol{{I}}_{a} \in \Re ^{n \times n}$
is a diagonal constant matrix which includes the equivalent moment of inertia of rotor and transmission system referred to the link sides.
${\boldsymbol{{C}}(\boldsymbol{{q}},\dot{\boldsymbol{{q}}})} \in \Re ^{n \times n}$
is the centrifugal and Coriolis matrix of the links.
$\boldsymbol{{g(q)}} \in \Re ^{n \times 1}$
is the gravity vector of the links.
$\boldsymbol{\tau }_{f} \in \Re ^{n \times 1}$
is the joint friction.
$\boldsymbol{\tau } \in \Re ^{n \times 1}$
is the joint torques.
In order to improve the identification results of dynamic model, a nonlinear friction model [Reference Han, Wu, Liu and Xiong27,Reference Farhat, Mata, Page and Valero28] is used here instead of the common linear friction model [Reference Zhang, Wang, Jing and Tan29,Reference Xu, Fan, Chen, Ng, Ang, Fang, Zhu and Zhao30], that is, for the
$j$
-th joint,
\begin{align} \tau _{f,j} &= \left(F_{c,j}+F_{v,j}|\dot{q}_j|^{\alpha _{j}}\right)\mathrm{sign}(\dot{q}_j)+B_j \nonumber\\[3pt] &= \left [\mathrm{sign}(\dot{q}_j)\ |\dot{q}_j|^{\alpha _{j}}\mathrm{sign}(\dot{q}_j)\ 1\right ] \cdot \left [F_{c,j}\ F_{v,j}\ B_j\right ]^{\mathrm{T}}, \end{align}
where
$F_{c,j}$
,
$F_{v,j}$
, and
$B_j$
are the linear friction parameters.
$\alpha _{j}$
is the
$j$
-th element in the nonlinear friction parameter vector
$\boldsymbol{\alpha }$
. Then, (1) can be written as a linear relation to a set of dynamic parameters 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 dynamic parameter set. For joint
$j$
, the elements are
$\boldsymbol{\varphi }^j=[\boldsymbol{{l}}_j, m_{j}\boldsymbol{{r}}_j, m_j, I_{a,j}, F_{c,j}, F_{v,j}, B_j]^{\mathrm{T}}$
.
$\boldsymbol{{I}}_{j}\in \Re ^{3\times 3}$
is the inertia tensor matrix with respect to frame
$j$
but of which only six elements are considered due to its symmetry, that is,
$\boldsymbol{{l}}_{j} = \left [I_{xx,j}, I_{xy,j}, I_{xz,j}, I_{yy,j}, I_{yz,j}, I_{zz,j}\right ]$
.
$m_{j} \boldsymbol{{r}}_{j} = \left [m_{j} r_{x,j}, m_{j} r_{y,j},m_{j}r_{z,j} \right ]$
is the first moments with respect to frame
$j$
.
$m_j$
is the mass of the link.
$I_{a,j}$
is the
$j$
-th diagonal element of
$\boldsymbol{{I}}_a$
.
According to Gautier and Khalil [Reference Gautier and Khalil31] and Khalil and Bennis [Reference Khalil and Bennis32], some inertial parameters in
$\boldsymbol{\varphi }$
are completely unidentifiable and some others can be only identified in linear combinations. Therefore, (3) can be further simplified 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.
$p$
is the element number in
$\boldsymbol{\pi}\!.\ {\boldsymbol{{Y}}(\boldsymbol{{q}},\dot{\boldsymbol{{q}}},\ddot{\boldsymbol{{q}}},\boldsymbol\alpha )}\in \Re ^{n \times p}$
is the regressor matrix with full rank after removing the linearly correlated columns in
${\boldsymbol{{W}}(\boldsymbol{{q}},\dot{\boldsymbol{{q}}},\ddot{\boldsymbol{{q}}},\boldsymbol\alpha )}\in \Re ^{n \times 14n}$
. In the remainder of this article, all the regressor matrices will be represented by only the first letter to simplify the expression, such as
$\boldsymbol{{Y}}={\boldsymbol{{Y}}(\boldsymbol{{q}},\dot{\boldsymbol{{q}}},\ddot{\boldsymbol{{q}}},\boldsymbol\alpha )}$
.
3. The Proposed Payload Identification Framework
3.1. Dynamic parameter identification
The initial base parameter set of the robot can be firstly computed by the OLS method as follows:
\begin{equation} \boldsymbol{\pi }_\textrm{OLS}=(\overline{\boldsymbol{{Y}}}^{\mathrm{T}}\overline{\boldsymbol{{Y}}})^{-1}\overline{\boldsymbol{{Y}}}^{\mathrm{T}}\overline{\boldsymbol{\tau }}, \end{equation}
where
$\overline{\boldsymbol{{Y}}}\in \Re ^{Nn \times p}$
is the stacked regressor matrix and
$\overline{\boldsymbol{\tau }}\in \Re ^{Nn \times 1}$
is the stacked torque vector by collecting
$N$
samples of joint data.
Then, the torque residuals can be computed by:
\begin{equation} \overline{\boldsymbol{{R}}}=\overline{\boldsymbol{\tau }}-\overline{\boldsymbol{{Y}}}\boldsymbol{\pi }_\textrm{OLS}, \end{equation}
where
$\overline{\boldsymbol{{R}}}\in \Re ^{Nn \times 1}$
is the torque residual vector.
Assuming the noise in the measured torques of different joints to be correlated, the initial covariance matrix of the torque residual vector can be then estimated by:
\begin{equation} \boldsymbol{\Lambda }=\frac{\overline{\boldsymbol{{R}}}\cdot \overline{\boldsymbol{{R}}}^{\mathrm{T}}}{N-p}, \end{equation}
where
$\boldsymbol{\Lambda } \in \Re ^{Nn \times Nn}$
is the covariance matrix of the torque residual vector.
In order to reduce the effects of bad data on the identification results, additional data weighting matrix
$\boldsymbol{{P}}\in \Re ^{Nn \times 1}$
and its extended matrix
$\boldsymbol{{P}}_{e}\in \Re ^{Nn \times p}$
with each column being
$\boldsymbol{{P}}$
can be applied as:
\begin{equation} \left \{\begin{array}{l} \widetilde{\boldsymbol{{Y}}}=\boldsymbol{{P}}_{e}\odot (\boldsymbol{\Lambda }^{-\frac{1}{2}}\overline{\boldsymbol{{Y}}}),\\ \\[-9pt] \widetilde{\boldsymbol{\tau }}=\boldsymbol{{P}}\odot (\boldsymbol{\Lambda }^{-\frac{1}{2}}\overline{\boldsymbol{\tau }}), \end{array}\right. \end{equation}
where
$\odot$
is the elementwise multiplication operator.
Then, the base parameter set can be solved by the WLS method with considering data weighting as follows:
\begin{equation} \boldsymbol{\pi }_\textrm{WLS}=(\widetilde{\boldsymbol{{Y}}}^{\mathrm{T}}\widetilde{\boldsymbol{{Y}}})^{-1}\widetilde{\boldsymbol{{Y}}}^{\mathrm{T}}\widetilde{\boldsymbol{\tau }}. \end{equation}
When the solution of
$\boldsymbol{\pi }_\textrm{WLS}$
is ready, the normalized residual vector can be computed as follows:
\begin{equation} \widetilde{\boldsymbol{{R}}}=\widetilde{\boldsymbol{\tau }}-\widetilde{\boldsymbol{{Y}}}\boldsymbol{\pi }_\textrm{WLS}. \end{equation}
Then, the covariance matrix
$\boldsymbol{\Lambda }$
of the torque error vector can be updated with its initial value computed by (7) as follows:
\begin{equation} \boldsymbol{\Lambda }_\textrm{new}=\frac{\boldsymbol{\Lambda }^{\frac{1}{2}}\widetilde{\boldsymbol{{R}}}\cdot \widetilde{\boldsymbol{{R}}}^{\mathrm{T}}\boldsymbol{\Lambda }^{\frac{1}{2}}}{N-p}. \end{equation}
With the normalized residual vector
$\widetilde{\boldsymbol{{R}}}$
as input, the new data weighting matrix
$\boldsymbol{{P}}_\textrm{new}$
can be updated with initial
$\boldsymbol{{P}}$
being an identity vector as follows:
\begin{equation} \boldsymbol{{P}}_\textrm{new}=\boldsymbol{{P}}\odot \Phi (\widetilde{\boldsymbol{{R}}},\delta ), \end{equation}
where
$\Phi (\cdot )$
is a function. It can return a vector with size of
$Nn\times 1$
, in which the elements will be set to 0 when the corresponding elements in
$\widetilde{\boldsymbol{{R}}}$
exceed to a threshold
$\delta$
, otherwise they will be set to 1.
In order to compute the estimated friction torque, one way is to set the friction parameters to zero after identifying
$\boldsymbol{\pi }_\textrm{WLS}$
. Then, the torques related to the inertial component in
$\boldsymbol{\pi }_\textrm{WLS}$
were subtracted from the measured ones. In addition, the data weighting matrices are applied here to make the estimated friction fit the friction model (2) better, that is,
\begin{equation} \overline{\boldsymbol{{f}}}=\boldsymbol{{P}}\odot \overline{\boldsymbol{\tau }}-(\boldsymbol{{P}}_{e}\odot \overline{\boldsymbol{{Y}}})\cdot \boldsymbol{\pi }_{\textrm{WLS},i}, \end{equation}
where
$\overline{\boldsymbol{{f}}} \in \Re ^{Nn \times 1}$
is the estimated friction torque vector.
$\boldsymbol{\pi }_{\textrm{WLS},i}$
is the identified inertial parameters in
$\boldsymbol{\pi }_\textrm{WLS}$
.
Finally, the nonlinear friction parameter vector
$\boldsymbol{\alpha }$
can be solved by the following nonlinear least squares problem, for
$j$
-th joint, that is,
\begin{equation} \begin{array}{r@{\quad }l} \min \limits _{\boldsymbol{\vartheta }_{j}} & \Vert \boldsymbol{\tau }_{f,j}-\overline{\boldsymbol{{F}}}_{j} \Vert _2\\ \\[-8pt]\mathrm{s.t.} & 0\lt \alpha _{j}\lt 1\\ \\[-8pt]& j=1,2,\cdots,n \end{array}, \end{equation}
where
$\boldsymbol{\vartheta }_{j}=[\alpha _{j},F_{c,j},F_{v,j},B_{j}]$
is the optimization variable. The initial value of
$\alpha _j$
is 1.
$\overline{\boldsymbol{{F}}} \in \Re ^{n \times N}$
is the reshaped version of
$\overline{\boldsymbol{{f}}}$
.
$\overline{\boldsymbol{{F}}}_j$
is the
$j$
-th row of
$\overline{\boldsymbol{{F}}}$
.
Compared with ref. [Reference Han, Wu, Liu and Xiong27], there are mainly three differences in this subsection. Firstly, the method in ref. [Reference Han, Wu, Liu and Xiong27] consists of three loops. The outer loop was used for updating nonlinear friction coefficients, the middle loop was used for updating the data weighting matrix, and the inner loop was used for updating the covariance matrix. In this subsection, however, the original middle loop and inner loop in ref. [Reference Han, Wu, Liu and Xiong27] are merged into one loop. It means that the data weighting matrix and the covariance matrix are updated in one loop. In this way, one can greatly improve the speed of convergence without effecting the accuracy of the identified parameters. In addition, the initial covariance matrix is estimated by (7) here instead of an identity matrix in ref. [Reference Han, Wu, Liu and Xiong27]. Secondly, the actual currents are used for identification here instead of the command currents in ref. [Reference Han, Wu, Liu and Xiong27]. Although the command currents have low noise, the command currents will not change even though a payload is fixed on the robot compared with the case without payload. This is because the command currents are related to the command trajectories. If the command trajectories are the same in both the case with payload and the case without payload, the command currents will be also the same in these two cases. However, the payload identification methods rely on the current variations caused by the payload. Therefore, the identification results of dynamic parameters cannot be applied to the payload identification if command currents were used. Finally, the data weighting is also considered in friction torque estimation formula (13) here, which can reduce the outliers and improve the fitting accuracy of friction model. But data weighting was not considered in the friction torque estimation formula in ref. [Reference Han, Wu, Liu and Xiong27].
3.2. Payload parameter identification
Based on the identified results of dynamic parameters, an accurate payload identification model based on double weighting is then constructed. By considering the condition with payload fixed on the end effector of the robot, two sets of equations can be then grouped. The first one uses a trajectory without payload (
$N_1$
data samples). The second one uses a trajectory with payload (
$N_2$
data samples). In addition, the variations of friction parameters caused by the payload are also considered in the model to obtain better results. Then, in order to obtain the WLS solution, the converged covariance matrices are used to normalize the constructed model, which is called the first weighting as follows:
\begin{equation} \left [\begin{array}{l} \boldsymbol{\Lambda }_a^{-\frac{1}{2}}\overline{\boldsymbol{\tau }}_a \\ \\[-8pt] \boldsymbol{\Lambda }_b^{-\frac{1}{2}}\overline{\boldsymbol{\tau }}_b \end{array}\right ] = \left [\begin{array}{l@{\quad}l@{\quad}l} \boldsymbol{\Lambda }_a^{-\frac{1}{2}}\overline{\boldsymbol{{Y}}}_a & \textbf{0}_{N_1n\times 10} & \textbf{0}_{N_1n\times 3n} \\ \\[-8pt] \boldsymbol{\Lambda }_b^{-\frac{1}{2}}\overline{\boldsymbol{{Y}}}_b & \boldsymbol{\Lambda }_b^{-\frac{1}{2}}\overline{\boldsymbol{{W}}}_{L} & \boldsymbol{\Lambda }_b^{-\frac{1}{2}}\overline{\boldsymbol{{W}}}_{f} \end{array}\right ] \cdot \boldsymbol{\phi }_\textrm{tot}, \end{equation}
where the subscript
$a$
and
$b$
mean the case without payload and the case with payload, respectively.
$\boldsymbol{\Lambda }_a$
and
$\boldsymbol{\Lambda }_b$
are obtained from the above dynamic parameter identification.
$\boldsymbol{\phi }_\textrm{tot}=[\boldsymbol{\pi }^{\mathrm{T}}\ \boldsymbol{\varphi }_{L}^{\mathrm{T}}\ \Delta \boldsymbol{\varphi }_{f}^{\mathrm{T}}]^{\mathrm{T}} \in \Re ^{(p+10+3n) \times 1}$
is the total parameter set.
$\boldsymbol{\pi }$
is the base parameter set in the case without payload.
$\boldsymbol{\varphi }_L \in \Re ^{10 \times 1}$
is the inertial parameters of the payload and
$\overline{\boldsymbol{{W}}}_{L} \in \Re ^{N_2n \times 10}$
is the corresponding payload regressor matrix.
$\Delta \boldsymbol{\varphi }_{f} \in \Re ^{3n \times 1}$
is the variations of linear friction parameters caused by the payload and
$\overline{\boldsymbol{{W}}}_{f} \in \Re ^{N_2n \times 3n}$
is the corresponding friction regressor matrix. The variations of nonlinear friction parameters are considered in
$\overline{\boldsymbol{{Y}}}_a$
and
$\overline{\boldsymbol{{Y}}}_b$
.
In order to reduce the effects of outliers on the results, the converged data weighting matrices are used to further improve the model, which is called the second weighting as follows:
\begin{equation} \left [\begin{array}{l} \boldsymbol{{P}}_a\odot \overline{\boldsymbol{\tau }}_a^{\ast } \\ \\[-8pt] \boldsymbol{{P}}_b\odot \overline{\boldsymbol{\tau }}_b^{\ast } \end{array}\right ] = \left [\begin{array}{l@{\quad}c@{\quad}c} \boldsymbol{{P}}_{e,a}\odot \overline{\boldsymbol{{Y}}}_a^{\ast } & \textbf{0}_{N_1n\times 10} & \textbf{0}_{N_1n\times 3n} \\ \\[-8pt] \boldsymbol{{P}}_{e,b}\odot \overline{\boldsymbol{{Y}}}_b^{\ast } & \boldsymbol{{P}}_{e,b_1}\odot \overline{\boldsymbol{{W}}}_{L}^{\ast } & \boldsymbol{{P}}_{e,b_2}\odot \overline{\boldsymbol{{W}}}_{f}^{\ast } \end{array}\right ] \cdot \boldsymbol{\phi }_\textrm{tot}, \end{equation}
where superscript
$\ast$
denotes the results weighted by covariance matrices. For example,
$\overline{\boldsymbol{\tau }}_a^{\ast }=\boldsymbol{\Lambda }_a^{-\frac{1}{2}}\overline{\boldsymbol{\tau }}_a$
and
$\overline{\boldsymbol{\tau }}_b^{\ast }=\boldsymbol{\Lambda }_b^{-\frac{1}{2}}\overline{\boldsymbol{\tau }}_b$
. The data weighting matrix
$\boldsymbol{{P}}_a$
(
$\boldsymbol{{P}}_b$
) and its extended matrix
$\boldsymbol{{P}}_{e,a}$
(
$\boldsymbol{{P}}_{e,b}$
) are converged and obtained from the above dynamic parameter identification.
$\boldsymbol{{P}}_{e,b_1}\in \Re ^{N_2n \times 10}$
and
$\boldsymbol{{P}}_{e,b_2}\in \Re ^{N_2n \times 3n}$
are the matrices with each column being
$\boldsymbol{{P}}_b$
. With a little abuse of notation, (16) can be concisely expressed as follows:
\begin{equation} \widetilde{\boldsymbol{\tau }}_\textrm{tot}=\widetilde{\boldsymbol{{W}}}_\textrm{tot}\boldsymbol{\phi }_\textrm{tot}. \end{equation}
Then, the total parameter set
$\boldsymbol{\phi }_\textrm{tot}$
can be solved as follows:
\begin{equation} \boldsymbol{\phi }_\textrm{tot}=(\widetilde{\boldsymbol{{W}}}_\textrm{tot}^{\mathrm{T}}\widetilde{\boldsymbol{{W}}}_\textrm{tot})^{-1}\widetilde{\boldsymbol{{W}}}_\textrm{tot}^{\mathrm{T}}\widetilde{\boldsymbol{\tau }}_\textrm{tot}. \end{equation}
Finally, the identified inertia parameters of the payload can be extracted from the
$(p+1)$
-th element to
$(p+10)$
-th element of
$\boldsymbol{\phi }_\textrm{tot}$
as follows:
\begin{equation} \boldsymbol{\varphi }_{L}=\boldsymbol{{S}}\boldsymbol{\phi }_\textrm{tot}, \end{equation}
where
$\boldsymbol{{S}}\in \Re ^{10\times (p+10+3n)}$
is the selection matrix to extract
$\boldsymbol{\varphi }_{L}$
from
$\boldsymbol{\phi }_\textrm{tot}$
.
3.3. The overall approach
The proposed identification framework for inertial parameters of robot payloads is shown in Algorithm 1. In order to obtain the weighting matrices, the dynamic parameter identification needs to be firstly conducted in both the case without payload and the case with payload. For description convenience, a function named
$\mathrm{DynIdent}(\cdot )$
is defined, in which two while loops are included.
The proposed payload identification algorithm
In the inner while loop, the WLS solution
$\boldsymbol{\pi }_\textrm{WLS}$
of base parameters is constantly solved and updated by (9). When the solution of
$\boldsymbol{\pi }_\textrm{WLS}$
is ready, the normalized residual vector
$\widetilde{\boldsymbol{{R}}}$
can be solved by (10). Then, the covariance matrix can be updated by (11). The data weighting matrix
$\boldsymbol{{P}}$
with its extended version
$\boldsymbol{{P}}_e$
can be updated by (12). Initially,
$\boldsymbol{{P}}$
is a vector in which all the elements are 1. As iterations go within this loop,
$\boldsymbol{{P}}$
converges to a constant vector. Only data points that fit the dynamic model and have low measurement noise will reside. Thus, the effects of bad data on the identification results can be reduced. In the outer while loop, the nonlinear friction parameters
$\boldsymbol{\alpha }$
can be solved and updated by (14). Initially,
$\boldsymbol{\alpha }$
is a vector in which all the elements are 1. As iterations go within this loop,
$\boldsymbol{\alpha }$
converges to a constant vector. Then, the friction can get properly modeled, which can improve the identification accuracy of base parameters.
In order to identify the inertial parameters of the payload, two sets of dynamic parameters of the robot have to be identified by running function
$\mathrm{DynIdent}(\cdot )$
twice. One set was identified by the data without payload. The other was identified by the data with payload. Then, the converged covariance matrices
$\boldsymbol{\Lambda }_a$
and
$\boldsymbol{\Lambda }_b$
were used for constructing model (15) as the first weighting. In addition, the variations of both nonlinear friction parameters and linear friction parameters caused by the payload were also included. Afterwards, the converged data weighting matrices
$\boldsymbol{{P}}_a$
(
$\boldsymbol{{P}}_{e,a}$
) and
$\boldsymbol{{P}}_b$
(
$\boldsymbol{{P}}_{e,b}$
) were applied to model (16) as the second weighting. Based on the linear relationship (17), the total parameter set can be solved by (18). Finally, the payload parameters can be obtained by (19).
4. Experiments
4.1. Experiment setup
The experiment platform is UR10 robot by Universal Robots in this article. As shown in Fig. 1, the modified Denavit–Hartenberg (DH) convention is used to build the kinematic frames of the robot. The modified DH parameters are listed in Table I. The robot is controlled by the
$ur\_rtde$
software package [Reference Lindvig33], which is a fast and lightweight interface for controlling and receiving data from the controller of the UR robot.
The modified DH parameters of the UR10 robot.
Modified DH model of the UR10 robot (
$n=6$
).
4.2. Data acquisition and data processing
In order to obtain a well-conditioned regressor matrix, the trajectory was parametrized using Fourier series [Reference Swevers, Ganseman, Tukel, De Schutter and Van Brussel34]. The period of each trajectory is set to 10 s. The condition numbers of regressor
$\boldsymbol{\Lambda }_0^{-\frac{1}{2}}\overline{\boldsymbol{{Y}}}$
and its subregressors are chosen to be the optimization criterion.
$\boldsymbol{\Lambda }_0$
is the initial value computed by (7). In addition, the optimization was performed iteratively so that optimized motions at the present step could complement previous motions [Reference Han, Wu, Liu and Xiong27]. By considering the optimization problem described above, 15 motions are designed and experimented to collect data, in which 10 motions are used for identification as shown in Fig. 2, which give the condition numbers of
$\mathrm{cond}(\boldsymbol{\Lambda }_0^{-\frac{1}{2}}\overline{\boldsymbol{{Y}}})=85.75$
. The remaining five motions are used for model validation. Then, the velocity commands are sent to the robot controller to track the designed trajectories. In order to construct model (15), the robot needs to track the designed trajectories twice, namely, one without payload and the other with payload. Meanwhile, the motor currents and the joint positions are collected with a sampling frequency of 250 Hz on the host computer. Since the computation of regressor matrix requires joint velocities and joint accelerations, the frequency filtering method originally proposed in ref. [Reference Swevers, Verdonck and De Schutter10] is used to compute them.
Assuming a linear relationship between motor currents
$\boldsymbol{{i}}$
and joint torques
$\boldsymbol{\tau }$
, the joint torques can be estimated by the product of joint drive gains
$\boldsymbol{{K}}$
and motor currents, namely
$\boldsymbol{\tau }=\boldsymbol{{K}}\boldsymbol{{i}}$
[Reference Gautier and Briot35,Reference Briot and Gautier36]. The joint drive gains
$\boldsymbol{{K}}$
of UR10 robot have been identified in ref. [Reference Xu, Fan, Fang, Zhu and Zhao37]. It was based on the current-level dynamic identification and considered the effects of the payload on the joint currents. Then, the identification model was constructed by two dynamic equations, in which one is with respect to the case without payload and the other is with respect to the case with payload. Because the identification of joint drive gains is not the main contribution in this paper, the identification process of
$\boldsymbol{{K}}$
is omitted. One can refer to ref. [Reference Xu, Fan, Fang, Zhu and Zhao37] for more details. The identified results are shown in Table II. Since
$\boldsymbol{{K}}$
will amplify the noise in motor currents, the estimated joint torques are then filtered by a low-pass Butterworth filter with an order of 4 and a cutoff frequency of 40 Hz. In addition, the threshold
$\delta$
in
$\Phi (\cdot )$
is set to 2.5 in this paper.
Joint drive gains
$\boldsymbol{{K}}$
of the UR10 robot.
Ten identification trajectories for the UR10 robot.
4.3. Dynamic parameter identification results
By using the method proposed in Section 3.1, the base parameters
$\boldsymbol{\pi }$
and the nonlinear friction parameters
$\boldsymbol{\alpha }$
of UR10 robot were identified. The identification results with standard variances of percentage
$\sigma \%$
have been shown in Table III. Then, in order to validate the model prediction accuracy, the measured torques and the predicted ones for validation trajectories are shown in Fig. 3. The measured torques are estimated by the motor currents without filtering. The predicted torques are computed by the identified dynamic parameters. Figure 3 shows that the predicted torques can fit the measured torques accurately. It reveals that the identified dynamic parameters are accurate by the proposed method.
Identified
$\boldsymbol{\pi }$
and
$\boldsymbol{\alpha }$
of the UR10 robot1.
1 All the units are in SI unit.
$\sigma =\sqrt{\mathrm{diag}((\widetilde{\boldsymbol{{Y}}}^{\mathrm{T}}\widetilde{\boldsymbol{{Y}}})^{-1}})$
.
$N_1=\textrm{10,253},N_2=\textrm{10,241},p=58$
.
In order to compare the method in ref. [Reference Han, Wu, Liu and Xiong27] (named Method 1) with the modified method in Section 3.1 (named Method 2), the root mean square errors (RMSE) of torque residuals of two methods are shown in Table IV. The torque residuals are computed by the differences between predicted torques and measured torques. Notably, Method 2 is the modified dynamic parameter identification method based on two loop iterations proposed in the Section 3.1. It is used to identify the weighting matrices and the nonlinear friction parameters. From Table IV, it can be found that the RMSE of two methods are both small and close. It indicates that the identification accuracy is almost not reduced even though only two loops are used in Method 2. In fact, it is not necessary to use another loop to update the covariance matrix separately. It can be updated in the same loop with the data weighting matrices together. When the initial covariance matrix was estimaed by (7) instead of an identity matrix, it also helps to avoid a separately converged process of the covariance matrix. In addition, the iteration numbers and the iteration time of two methods are shown in Table V. The codes were tested on Windows 10 with processor Intel(R) Core(TM) i5-6500 CPU @ 3.20 GHz and 24-GB RAM. From Table V, it can be found that the iteration numbers and the iteration time are greatly reduced in Method 2 compared with Method 1. Specifically, the iteration numbers were reduced by 70.48
$\%$
and the iteration time was reduced by 59.87
$\%$
. The results show a significant improvement in our method.
RMSE1 of torque residuals for validation trajectories without payload.
1 The unit of RMSE is
$\mathrm{Nm}$
.
Iteration numbers and iteration time of two different methods.
Cross-validation of identified results for validation trajectories without payload (joints 1
$\sim$
6 are from the top to the bottom).
In order to observe the effects of data weighting on the identification results, the identified friction torques versus joint velocities are shown in Fig. 4. From the figure, it can be found that the friction torques estimated without data weighting (blue dots and black circles together) fit the friction model (2) (red dots) badly because of the outliers. However, the ones estimated with data weighting (only blue dots) fit the friction model (2) well. It is because the outliers are eliminated by the data weighting matrices. The results reveal that the data weighting can reduce the effects of outliers and significantly improve the identification results. In addition, the change of elements in data weighting matrix
$\boldsymbol{{P}}$
was shown in Fig. 5. From the figure, it can be found that the sum of
$\boldsymbol{{P}}_j$
decreases as the iterations go. It also means that the outliers are eliminated. Moreover, the sum of
$\boldsymbol{{P}}_j$
for joint 1 and joint 4
$\sim$
joint 6 is smaller than that for joint 2
$\sim$
joint 3. It indicates that more outliers in the four joints are eliminated than that in another two joints. This is consistent with the presentation in Fig. 4. In addition, as shown in Fig. 6, the variances of torque residuals for each joint also decrease. This is because the outliers are eliminated as the iterations go. It means that the variances of torque residuals become smaller and the identification results become more accurate. Similarly, the reduction multiple of variances of torque residuals in the four joints are larger than that in another two joints. To save pages and avoid repeating, only the model validation results without payload are shown.
Identified friction torques versus joint velocities in the case without payload. The blue dots denote the estimated friction torques with data weighting. The blue dots and the black circles together denote the estimated friction torques without data weighting. The black circles denote the outliers. The red dots denote the fitted friction model (2) with data weighting.
The sum of
$\boldsymbol{{P}}_j$
in each iteration of outer loop in the case without payload.
$\boldsymbol{{P}}_j$
denotes the elements for joint
$j$
in
$\boldsymbol{{P}}$
. The outer loop is to update nonlinear friction coefficient
$\boldsymbol{\alpha }$
. The total number of iterations in outer loop is 9. The initial sum of
$\boldsymbol{{P}}_j$
is all 10,253, namely the number of data points.
$\boldsymbol{\Lambda }_j$
in each iteration of outer loop in the case without payload.
$\boldsymbol{\Lambda }_j$
denotes the variance of torque residuals for joint
$j$
.
4.4. Payload parameter identification results
Method 3 is the classical method in ref. [Reference Swevers, Verdonck and De Schutter10]. It is based on the observation that the robot links and the payload each contributed separately to the joint torques. The base parameter set of the robot without payload can be estimated in advance. Then, it can be considered as a priori information for the payload identification. Method 4 is one of the two impressive methods in ref. [Reference Khalil, Gautier and Lemoine11]. It is to compute the differences of joint torques with the same trajectories twice. Method 5 is the other of the two impressive methods in ref. [Reference Khalil, Gautier and Lemoine11]. It is to estimate robot link parameters and payload parameters together in one step based on the WLS method. Method 6 is the method in ref. [Reference Gaz and De Luca13]. It is based on the linear relationship between the differences of the robot static coefficients and the inertial parameters of payload. In particular, a simple static friction model
$\tau _{f,j}=F_{c,j}\mathrm{sign}(\dot{q}_j)+B_j$
is considered in Method 6 instead of originally ignoring joint static friction effects in ref. [Reference Gaz and De Luca13]. This is because the UR10 robot is not equipped with joint torque sensors. The brief descriptions about these four methods can be also found in the Section 1. The proposed method based on double weighting in this paper is named Method 7. The 10 inertial parameters of payload in Fig. 7 estimated by these methods are shown in Table VI. The relative errors of both the mass and the CoM in
$Z$
axis of the payload estimated by these methods are shown in Fig. 8.
Concentric payload fixed on the end effector of the UR10 robot.
Relative errors of
$m_L$
and
$r_{z,L}$
estimated by different methods.
As shown in Fig. 8, the relative errors of both the mass and the CoM in
$Z$
axis for Method 7 are all the smallest ones by comparing Method 7 with another four methods. Specifically, the mass relative error of Method 7 is about 10 times smaller than that of Method 5 and about 43 times smaller than that of Method 4. Moreover, the CoM in
$Z$
axis of the payload estimated by Method 7 is more accurate than that by another four methods. For the inertia tensor of the payload, no one method can match the true values accurately. In fact, it is usually difficult to accurately estimate the inertia tensor of the payload, especially when the payload is lightweight with respect to the robot mass. But the ones estimated by Method 5 and Method 7 are more close to the true values than that by Method 3 and Method 4. In addition,
$I_{xz,L}$
,
$I_{yy,L}$
. and
$I_{zz,L}$
estimated by Method 7 are better than that by Method 5. The final results reveal that our method is better than another four methods. The reasons are discussed below.
Based on the dynamic parameter identification results, the converged covariance matrices and the converged data weighting matrices in two cases are obtained. As shown in Fig. 6, the variances become minimum when the covariance matrices are converged. When the converged covariance matrices are weighted to the payload identification model, the WLS solution can be then obtained. As shown in Fig. 5, the outliers are eliminated when the data weighting matrices are converged. When the converged data weighting matrices are weighted to the payload identification model, the effects of outliers can be reduced and the results can be improved. In addition, the nonlinear friction coefficients in two cases are also considered in the model. But the variations of nonlinear friction parameters
$\boldsymbol{\alpha }$
are not obvious except for
$\alpha _3$
as shown in Table III. Therefore, the influence of the variations of nonlinear friction parameters on payload identification results might be slight. But the nonlinear friction model is useful to improve the dynamic parameter identification results [Reference Han, Wu, Liu and Xiong27]. The variations of linear friction parameters caused by the payload are also included in the model. In fact, as shown in Table III, 11 bold parameters are significantly affected by the payload. These parameters include seven inertial parameters (
$\pi _1$
,
$\pi _2$
,
$\pi _6$
,
$\pi _7$
,
$\pi _9$
,
$\pi _{13}$
, and
$\pi _{14}$
) of the robot and four linear friction parameters (
$\pi _{41}$
,
$\pi _{43}$
,
$\pi _{48}$
, and
$\pi _{49}$
). Therefore, compared with the nonlinear friction parameters, the variations of linear friction parameters have more effects on the payload identification results. The identification results of the payload can be improved by considering the variations. In fact, the double weighting technique and the considered friction variations benefit to improve the accuracy of payload identification in our method.
5. Conclusion
In this paper, an accurate identification method based on double weighting for the inertial parameters of robot payloads is proposed. Firstly, in order to obtain the weighting matrices and nonlinear friction coefficients, a modified dynamic parameter identification method with two loops is proposed. Compared with the original one in ref. [Reference Han, Wu, Liu and Xiong27], the iteration speed in the modified one can be significantly improved without loss of accuracy. Based on the dynamic identification results, a payload identification model based on double weighting is then constructed. The first weighting is to use the converged covariance matrices to normalize the constructed model. It can help to obtain the WLS solution. The second weighting is to use the converged data weighting matrices to make the data points fit the theoretical model well. It can reduce the effects of outliers and help to improve the identification results. In addition, the variations of both nonlinear friction parameters and linear friction parameters caused by the payload are also included to improve the results. Finally, experimental comparisons between our method and another four methods are shown. The results indicate that the proposed framework shows the best performance, especially on improving the identification accuracy of the mass and the CoM of the payload. It is worth noting that the proposed method in this article cannot be applied to the real-time payload identification, since the model (15) can be only constructed after the sampling joint data of the manipulator in both the case without payload and the case with payload have been fully collected. In the future, the authors would try to modify and improve the method proposed in this article for a real-time application.
Acknowledgments
The authors would like to thank the reviewers for their professional comments on the improvement of this article.
Author Contributions
Conceptualization, T.X.; methodology, T.X.; formal analysis, T.X. and Q.F.; validation, T.X. and Y.Z.; 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 Research and Development Plan [grant number 2018YFB1306901]; and the National Natural Science Foundation of China [grant number 92048301].
Conflicts of Interest
The authors have no relevant financial or nonfinancial interests to disclose.
