Backward sequential approach for dynamic parameter identification of robot manipulators

Abstract

This article proposes a sequential optimization approach to efficiently identify non-minimal dynamic parameters of robot manipulators, possibly having large degrees of freedom. A back-substitution-based parameter identification from the last link to inward links is enabled due to the block upper triangular form of inherent regressor matrix. Starting with the dynamic model using the non-minimal parameters, we derive a generic compact formulation for the linear regression equation. We then establish a sequential optimization procedure taking into account physical feasibility of parameters. Numerical case examples demonstrate the validity of the proposed approach.

Introduction

Dynamic parameter identification of robot manipulators is an indispensable step for sophisticated motion control and better utilization of the robot systems. ¹ For dynamic parameter identification, we need to have three fundamental elements: (i) first an efficient form of dynamic equation that can be transformed into a linear regression equation in terms of the parameters, ^2,3 (ii) a thoughtfully designed experimental trajectory that fully excites the dynamic nature of the robot system, ⁴ and (iii) a proper numerical method to solve for a feasible set of dynamic parameters by using input torques and motion data. ^{5 –7} Our study is closely related to efficiently constructing and solving the regression equation while maintaining physical aptness or feasibility of the solution.

Among abundant researches on parameter identification, ⁸ some noticeable works relevant to our study are listed as follows. Armstrong et al. ⁹ proposed a minimal set of dynamic parameters of the PUMA560 industrial robot from the symbolic form of whole dynamic equation, even if the minimal set showed a nonlinear combination of the link inertial parameters. Gautier and Khalil ¹⁰ proposed a systematic way to select the set of minimally independent parameters (the base parameters) by the nullspace decomposition of the regressor matrix and later Yoshida and Khalil ¹¹ showed that in some cases, the solution of the minimal set of parameters could violate the positive definiteness of the inertial matrix, leading to a physically incorrect result. Mata et al. ¹² proposed a quadratic programming method to obtain a feasible set of parameters that best approximates the unconstrained least square solution; and they showed its validity by testing on a 6-degree-of-freedom (DoF) industrial robot manipulator. Ayusawat and Nakamura ¹³ tried to identify a set of minimal parameters for a large-scaled robot systems by solving an optimization problem fulfilling both the least square error and the feasibility condition approximated by convex cones. Venture and Gautier ¹⁴ proposed a singular value decomposition–based calibration method for human inertial parameters, so as to satisfy physical consistency with a priori anatomical data. Gautier et al. proposed an observer-like closed-loop dynamic parameter identification method ¹⁵ which does not require the measurements (or estimates) of joint velocity and acceleration. Sousa and Cortesão ⁷ formulated the identification problem as a semi-definite programming (SDP) with linear inequality constraints that satisfy feasible conditions. Their algorithm was verified on three-link and seven-link robot manipulators. More recently, Joukov et al. ¹⁶ used an extended Kalman filter method for parameter identification of robot manipulators after imposing physical feasibility conditions into Sigmoid functions.

These previous approaches, however, would not be adequate nor efficient for general cases that possibly possess a large number of links. Although Ayusawat and Nakamura ¹³ tried to find the inertial parameters of a large DoF system, the algorithm basically required precise bounding box information on each link, which is not possible. Reorganizing all parameters into a set of minimal parameters, as was done for Programmable Universal Machine for Assembly (PUMA) robot ⁹ and KUKA robot, ¹⁷ must be even farther from possibility when the number of links becomes large. So, in this article, we propose an efficient backward sequential identification approach that divides the whole identification problem into several smaller ones, by utilizing the upper triangular nature of the regression matrix. ¹⁸ The regressor associated with non-minimal parameters is constructed from our compact analytic form of dynamics. The method is general and scalable to robot manipulators of any DoFs. The computation time is linearly increasing with the number of DoFs, at least, in theory.

The organization of the rest of this article is as follows. “Dynamic equations of robot manipulators” section addresses a succinct form of dynamics suitable for constructing regressor form; “Sequential procedure for parameter identification” section introduces the sequential identification of dynamic parameters in detail; “Case examples” section illustrates some case examples to show the validity of the proposed method; and finally, “Concluding remarks” section makes the conclusion of our study.

Dynamic equations of robot manipulators

Revisit: Newton–Euler equation

Consider an intermediate link i of an n-DoF robot manipulator as shown in Figure 1, where $f_{i,\text{eq}}$ and $\tau { }_{i,\text{eq}}$ are the equivalent external force and moment exerted at the center of mass (O_i) of the body, respectively. The linear momentum of link i is

$\nu { }_{v,i}=m_i{V}_i=m_i{J}_{vi} \dot{q}$

1

where m_i and V _i , respectively, denote the mass and the linear velocity at O_i, and $J_{vi}\in {ℝ}^{3\hspace{0.17em}\times \hspace{0.17em}n}$ is the kinematic map from joint velocity $\dot{q}$ to V _i . Similarly, the angular momentum of link i is

$\nu { }_{w,i}={\mathbb{I}}_i \mathit{\omega }{ }_i=R_i {\mathbb{I}}_i^i R_i^T{J}_{wi} \dot{q}$

2

where ω _i is the angular velocity of link i; ${\mathbb{I}}_i\in {ℝ}^{3\hspace{0.17em}\times \hspace{0.17em}3}$ and ${\mathbb{I}}_i^i\in {ℝ}^{3\hspace{0.17em}\times \hspace{0.17em}3}$ denote the inertia tensors of link i at O_i w.r.t. the reference frame, {0}, and w.r.t. the body frame, {i}, respectively; $R_i\in {ℝ}^{3\hspace{0.17em}\times \hspace{0.17em}3}$ is the rotation matrix from frame {0} to frame {i}; and $J_{wi}\in {ℝ}^{3\hspace{0.17em}\times \hspace{0.17em}n}$ is the kinematic map from joint velocity $\dot{q}$ to ω _i . Note that ${\mathbb{I}}_i^i$ is a constant matrix, whereas ${\mathbb{I}}_i$ is varying with the configuration of the robot.

Figure 1.

Free-body diagram of link i of a robot manipulator.

Now, taking time derivatives on equations (1) and (2) yields

$\left[\begin{array}{c} \dot{\nu }{ }_{vi}\\ \dot{\nu }{ }_{wi}\end{array}\right]=-\left[\begin{array}{c}{G}_{vi}\\ 0\end{array}\right]+\left[\begin{array}{c}{f}_{i,c}\\ \tau { }_{i,c}\end{array}\right]+\left[\begin{array}{c}{f}_{i,\text{eq}}\\ {\tau }_{i,\text{eq}}\end{array}\right]$

3

where $G_{vi}:=m_i{g}_r$ denotes the gravitational force with g _r being the gravitational acceleration vector, and f _i,c and $\tau { }_{i,c}$ , respectively, are the equivalent constraint force and moment acting on O_i. We can combine the above relation over the entire system as

$\dot{\nu } =-G+F_c+F$

4

where

$\begin{array}{ll} \nu \hfill & =[ \nu { }_{v,1}^T \nu { }_{v,2}^T\dots \nu { }_{v,n}^T| \nu { }_{w,1}^T \nu { }_{w,2}^T\dots \nu { }_{w,n}^T{]}^T\hfill \\ G\hfill & =[G_{\mathit{v},1}^{\mathit{T}}{G}_{\mathit{v},2}^{\mathit{T}}\dots G_{\mathit{v},\mathit{n}}^{\mathit{T}}|{\mathbf{0}}^{\mathit{T}}{\mathbf{0}}^{\mathit{T}}\cdots {\mathbf{0}}^{\mathit{T}}{]}^{\mathit{T}}\hfill \\ F_{\mathit{c}}\hfill & =[f_{\mathit{c},1}^{\mathit{T}}{f}_{\mathit{c},2}^{\mathit{T}}\dots f_{\mathit{c},\mathit{n}}^{\mathit{T}}| \tau { }_{\mathit{c},1}^{\mathit{T}} \tau { }_{\mathit{c},2}^{\mathit{T}}\dots \tau { }_{\mathit{c},\mathit{n}}^{\mathit{T}}{]}^{\mathit{T}}\hfill \\ F\hfill & =[f_{1,\text{eq}}^{\mathit{T}}{f}_{2,\text{eq}}^{\mathit{T}}\dots f_{\mathit{n},\text{eq}}^{\mathit{T}}| \tau { }_{1,\text{eq}}^{\mathit{T}} \tau { }_{2,\text{eq}}^{\mathit{T}}\dots \tau { }_{\mathit{n},\text{eq}}^{\mathit{T}}{]}^{\mathit{T}}\hfill \end{array}$

The momentum, on the other hand, can be cast into the following form

$\nu =\mathcal{M}(q)J(q) \dot{q}$

5

where

$\begin{array}{ll}\mathcal{M}(q)\hfill & =\hspace{0.17em}\text{blockdiag} [m_1{I}_{3\hspace{0.17em}\times \hspace{0.17em}3}\hspace{0.17em}\dots \hspace{0.17em}{m}_n{I}_{3\times 3}|{\mathbb{I}}_1\hspace{0.17em}\dots \hspace{0.17em}{\mathbb{I}}_n]\hfill \\ J(q)\hfill & ={\left[J_{v1}^T{J}_{v2}^T\hspace{0.17em}\dots \hspace{0.17em}{J}_{vn}^T|J_{w1}^T{J}_{w2}^T\hspace{0.17em}\dots \hspace{0.17em}{J}_{wn}^T\right]}^T\hfill \end{array}$

with $I_{3\hspace{0.17em}\times \hspace{0.17em}3}$ being the 3 × 3 identity matrix. Plugging equation (5) into equation (4) and then pre-multiplying $J^T(q)$ yield

$\begin{array}{r}{J}^T\mathcal{M}J\ddot{\mathbf{q}}\hspace{0.17em}+\hspace{0.17em}\left(J^T\dot{\mathcal{M}}J\hspace{0.17em}+\hspace{0.17em}{J}^T\mathcal{M}\dot{J}\right) \dot{q}\hspace{0.17em}=-J^{T}G\hspace{0.17em}+\hspace{0.17em}{J}^T{F}_c\hspace{0.17em}+\hspace{0.17em}{J}^{T}F\hspace{0.17em}\hspace{0.17em}\end{array}$

6

where we omitted the dependency of q in J and $\mathcal{M}$ for simplicity. Since the virtual work principle gives

$\tau =J^{T}F \text{and} J^T{F}_c=0$

where $\tau \in {ℝ}^n$ is the joint torque, equation (6) becomes

$M(q)\ddot{\mathbf{q}}+C(q, \hspace{0.17em}\dot{q} ) \dot{q} +g(q)= \tau$

7

where

$M=J^T\mathcal{M}J,C=J^T\dot{\mathcal{M}}J+J^T\mathcal{M}\dot{J}, \text{and} g=J^{T}G.$

If written in a summation form, equation (7) finally becomes

$\begin{array}{l}{\displaystyle \sum _{i=1}^n}\left(m_i{J}_{vi}^T{J}_{vi}+J_{wi}^T{R}_i{\mathbb{I}}_i^i{R}_i^T{J}_{wi}\right)\ddot{q}+{\displaystyle \sum _{i=1}^n}(m_i{J}_{vi}^T{\dot{J}}_{vi}\\ +J_{wi}^T{R}_i{\mathbb{I}}_i^i{R}_i^T{\dot{J}}_{wi}+\\ J_{wi}^{T}S( \mathit{\omega }{ }_i)R_i{\mathbb{I}}_i{R}_i^T{J}_{wi}) \dot{q} \\ +{\displaystyle \sum _{i=1}^n}{J}_{vi}^T{G}_{vi}= \tau \end{array}$

8

where $S( {\mathit{\omega }}_i):=( {\mathit{\omega }}_i\times )$ is a matrix operator of vector cross product.

In spite of its elegance, equation (8) is not useful in its current form to be cast into a parameter regression form, because J _vi ’s (and their derivatives) are not computable without prior knowledge of O_i, the center of mass of link i, which is, however, unknown before identification. In the next subsection, equation (8) shall be modified by translating the point of interest from O_i to a known kinematic position to yield an effective dynamic equation giving us a convenient regression model.

Regressor form from translated output coordinate

Suppose ${\overline{V}}_i$ be the linear velocity at E_i, the known distal end of link i, and also suppose r _i be the distance vector from E_i to O_i (refer to Figure 2). Then, it is true that

$V_i={\overline{V}}_i+{\mathit{\omega }}_i\times r_i=\left({\overline{J}}_{vi}-S(r_i)J_{wi}\right) \dot{q}$

9

where ${\overline{J}}_{vi}$ is defined such that ${\overline{V}}_i={\overline{J}}_{vi} \dot{q}$ . Thus, we obtain a Jacobian relation as

$J_{vi}={\overline{J}}_{vi}-S(r_i)J_{wi}.$

10

Figure 2.

Velocity of link i of a robot manipulator in translated coordinate.

Plugging the above two relations into equation (4) and pre-multiplying J^T yield

$\begin{array}{l}{\displaystyle \sum _{i\hspace{0.17em}=\hspace{0.17em}1}^n}{m}_i{\left({\overline{J}}_{vi}-S(r_i)J_{\mathit{\omega }i}\right)}^T\left({\dot{\overline{V}}}_i+{\dot{\mathit{\omega }}}_i\times r_i+{\mathit{\omega }}_i\times {\dot{r}}_i\right)+\\ {\displaystyle \sum _{i\hspace{0.17em}=\hspace{0.17em}1}^n}{J}_{wi}^T\left({\mathbb{I}}_i {\dot{\mathit{\omega }}}_i+{\dot{\mathbb{I}}}_i {\mathit{\omega }}_i\right)+\\ {\displaystyle \sum _{i\hspace{0.17em}=\hspace{0.17em}1}^n}{\left({\overline{J}}_{vi}^T-S(r_i)J_{\mathit{\omega }i}\right)}^T(m_i{g}_r)=\tau .\end{array}$

11

After expanding and doing mathematical manipulation, the above equation reduces to

$\begin{array}{l}{\displaystyle \sum _{i=1}^n}{m}_i{\overline{J}}_{vi}^T{\dot{\overline{V}}}_i+\\ {\displaystyle \sum _{i=1}^n}{m}_i\left({\overline{J}}_{vi}^{T}S({\dot{\omega }}_i)+{\overline{J}}_{vi}^{T}S({\omega }_i)S({\omega }_i)-J_{wi}^{T}S({\dot{\overline{V}}}_i)\right)r_i+\\ {\displaystyle \sum _{i=1}^n}{J}_{wi}^T\left(R_i{\hspace{0.17em}}^t{\mathbb{I}}_i^i{R}_i^T{\dot{\omega }}_i+S({\omega }_i)R_i{\hspace{0.17em}}^t{\mathbb{I}}_i^i{R}_i^T{\omega }_i\right)+\\ {\displaystyle \sum _{i=1}^n}{({\overline{J}}_{vi}^T-S(r_i)J_{\omega i})}^T(m_i{g}_r)=\tau \end{array}$

12

where ${\hspace{0.17em}}^t{\mathbb{I}}_i^i:={\mathbb{I}}_i^i-m_{i}S(r_i^i)S(r_i^i)$ means the invariant inertia tensor computed at the translated coordinate, E_i (this is called as Steiner’s theorem ¹⁹ ). Also, note that $r_i^i$ means r _i represented in the body frame {i}.

To derive a linear regression equation from equation (12), we rearrange the parameters of ${\hspace{0.17em}}^t{\mathbb{I}}_i^i$ , by following the way suggested by Lu and Meng, ²⁰ as

$\begin{array}{l}\left[\begin{array}{ccc}{\hspace{0.17em}}^t{I}_{xx}^i& {\hspace{0.17em}}^t{I}_{xy}^i& {\hspace{0.17em}}^t{I}_{xz}^i\\ {\hspace{0.17em}}^t{I}_{xy}^i& {\hspace{0.17em}}^t{I}_{yy}^i& {\hspace{0.17em}}^t{I}_{yz}^i\\ {\hspace{0.17em}}^t{I}_{xz}^i& {\hspace{0.17em}}^t{I}_{yz}^i& {\hspace{0.17em}}^t{I}_{zz}^i\end{array}\right]\left[\begin{array}{c}{d}_1\\ d_2\\ d_3\end{array}\right]=\\ \underset{B(d)}{\underbrace{\left[\begin{array}{cccccc}{d}_1& 0& 0& -d_2& -d_3& 0\\ 0& d_2& 0& -d_1& 0& -d_3\\ 0& 0& d_3& 0& -d_1& -d_2\end{array}\right]}}\left[\begin{array}{c}{\hspace{0.17em}}^t{I}_{xx}^i\\ {\hspace{0.17em}}^t{I}_{yy}^i\\ {\hspace{0.17em}}^t{I}_{zz}^i\\ {\hspace{0.17em}}^t{I}_{xy}^i\\ {\hspace{0.17em}}^t{I}_{xz}^i\\ {\hspace{0.17em}}^t{I}_{yz}^i\end{array}\right]\\ \leftrightarrow {}^t{\mathbb{I}}_i^{i}d=B(d){\text{γ}}_i\end{array}$

where $B(d)$ means a matrix operator for any vector $d\in {ℝ}^3$ and ${\gamma }_i\in {ℝ}^6$ is the vector of the parameters of the symmetric matrix ${\hspace{0.17em}}^t{\mathbb{I}}_i^i$ . Applying the above transformation to equation (12) yields

$\begin{array}{l}{\displaystyle \sum _{i=1}^n}{m}_i{\overline{J}}_{vi}^T\left({\dot{\overline{V}}}_i+g_r\right)+{\displaystyle \sum _{i=1}^n}{m}_i({\overline{J}}_{vi}^{T}S({\dot{\omega }}_i)+{\overline{J}}_{vi}^{T}S({\omega }_i)S({\omega }_i)\\ -J_{wi}^{T}S({\dot{\overline{V}}}_i+g_r))R_i{r}_i^i{\displaystyle \sum _{i=1}^n}{J}_{wi}^T\left(R_{i}B(\dot{\omega }{}_i^i)+S({\omega }_i)R_{i}B({\omega }_i^i)\right){\text{γ}}_i=\tau \end{array}$

13

where ${\dot{\mathit{\omega }}}_i^i=R_i^T{\dot{\mathit{\omega }}}_i$ is the angular acceleration with respect to frame {i}. All the terms in the above equation are expressed linearly in the inertial parameters: m_i, $m_i{r}_i^i$ , and γ _i .

To add the effect of the frictional forces in the joints, the external input τ can be modified to

$\tau ={\tau }_a- {\tau }_f$

14

where τ_a denotes the pure active joint torque and τ_f is the friction torque in the joints which is commonly modeled as

${\tau }_f= \mathit{\alpha } \text{sgn}(\dot{q})+\beta \dot{q} ={\displaystyle \sum _{i\hspace{0.17em}=\hspace{0.17em}1}^n}{H}_i{\mu }_i$

15

where $\mathit{\alpha }=\text{diag} [{\alpha }_1,\hspace{0.17em}\dots \hspace{0.17em},{\alpha }_n]$ and $\mathit{\beta }=\text{diag}[{\beta }_1,\hspace{0.17em}\dots \hspace{0.17em},{\beta }_n]$ are the matrices of coefficients of Coulomb’s friction and viscous friction, respectively; $\text{sgn}(\dot{q})$ is the component-wise vector sign function; $H_i\in {ℝ}^{n\hspace{0.17em}\times \hspace{0.17em}2}$ is a matrix whose j th row is $[\text{sgn}({\dot{q}}_i)1]$ for i = j, but is zero for j ≠ i; and ${\mathit{\mu }}_i=[{\alpha }_i\hspace{1em}{\beta }_i{]}^T$ is the vector of friction coefficients related to the i th joint.

Ultimately, by combining equations (13) to (15), the regression equation in the full dynamic parameters is constructed as

${\displaystyle \sum }_{i=1}^n{\Gamma }_i{Y}_i={\displaystyle \sum }_{i=1}^n[a_{i,1} a_{i,2}{a}_{i,3}{a}_{i,4}]\left[\begin{array}{c}{m}_i\\ m_i{r}_i^i\\ {\gamma }_i\\ {\mu }_i\end{array}\right]={\tau }_a$

16

where $a_{i,1}:={\overline{J}}_{vi}^T\left({\dot{\overline{V}}}_i+g_r\right)$ , $a_{i,2}:=\left({\overline{J}}_{vi}^{T}S({\dot{\mathit{\omega }}}_i)+{\overline{J}}_{vi}^{T}S({\mathit{\omega }}_i)\times S({\mathit{\omega }}_i)-J_{wi}^{T}S({\dot{\overline{V}}}_i+g_r)\right)R_i$ , $a_{i,3}:=J_{wi}^T\left(R_{i}B({\dot{\mathit{\omega }}}_i^i)+S({\mathit{\omega }}_i)R_{i}B({\mathit{\omega }}_i^i)\right)$ , and $a_{i,4}:=H_i$ ; also, ${\mathbf{\Gamma }}_i\in {ℝ}^{n\hspace{0.17em}\times \hspace{0.17em}12}$ and $Y_i\in {ℝ}^{12}$ , respectively, denote the regressor matrix and the dynamic parameter vector associated with link i. The regression equation (16) can be cast into a single big-sized matrix equation as

$\mathbf{\Gamma }Y={\tau }_a$

17

where $\mathbf{\Gamma }:=[{\mathbf{\Gamma }}_1{\mathbf{\Gamma }}_2\dots {\mathbf{\Gamma }}_n]\in {ℝ}^{n\hspace{0.17em}\times \hspace{0.17em}12n}$ and $Y:=[Y_1^T{Y}_2^T\dots Y_n^T{]}^T\in {ℝ}^{12n}$

It is important to note that all the rows below the i th row are null in Γ _i by the construction of ${\overline{J}}_{vi}$ and $J_{\mathit{\omega }i}$ . It means that inertial effect of link i is applied only to its inner joints, that is, joint 1 through joint i (see, for example, the works of Ha et al. ¹⁸ and Gautier ²¹ ). Ultimately, Γ becomes a block upper triangular matrix as shown below

$\mathbf{\Gamma } =\left[\begin{array}{ccccc}*& *& \cdots & *& *\\ \mathbf{0}& *& \cdots & *& *\\ ⋮& ⋮& ⋮& ⋮& ⋮\\ \mathbf{0}& \mathbf{0}& \cdots & *& *\\ \mathbf{0}& \mathbf{0}& \cdots & \mathbf{0}& *\end{array}\right].$

Sequential procedure for parameter identification

To identify the dynamic parameters, input torque and motion data of each joint have to be collected from a set of test motions, and using these, the dynamic parameters are sought such that the sum of squared errors between the actual input torque and the calculated torque by equation (17) is to be minimum. The necessary minimization is formulated as a least square problem

$\underset{\widehat{Y}}{{\displaystyle \text{min}}}J=\frac{1}{2}{\displaystyle \sum _{k=1}^N}{e}^T(k)e(k)$

18

where $e(k)\hspace{0.17em}=\hspace{0.17em}\tau { }_a(k)- \mathbf{\Gamma } (k)\widehat{Y}$ is the torque estimation error at t = kT, $\widehat{Y}$ is the estimate of Y, and N means the total number of measurements. (Hereafter, variables with $(\hspace{0.17em}\widehat{\hspace{0.17em}\cdot \hspace{0.17em}}\hspace{0.17em})$ would mean the estimates of the original variables.) By the block upper triangular nature of Γ(k), e_j(k), the j th component of e(k), can be written as

$e_j(k)={\tau }_{aj}(k)-{\displaystyle \sum _{i=j}^n} \text{row} ({\mathbf{\Gamma }}_i(k),j){\widehat{Y}}_i$

where ${\tau }_{aj}(k)$ denotes the jth component of $\tau { }_a(k)$ and $\text{row} ({\mathbf{\Gamma }}_i(k),j)$ is an operator extracting the jth row vector of ${\mathbf{\Gamma }}_i(k)$ . Rearranging the estimation errors and input torques for entire measurements as

$E=[E_1^T{E}_2^T\dots E_n^T{]}^T,T=[T_1^T{T}_2^T\dots T_n^T{]}^T$

where $E_j:=[e_j(1)e_j(2)\dots e_j(N{)]}^T$ and $T_j:=[{\tau }_{aj}(1){\tau }_{aj}(2)\dots {\tau }_{aj}(N{)]}^T$ , we establish

$E=T-\widetilde{ \mathbf{\Gamma } }\widehat{Y}=T-\left[\begin{array}{cccc}{\tilde{\mathbf{\Gamma }}}_{11}& {\tilde{\mathbf{\Gamma }}}_{12}& \cdots & {\tilde{\mathbf{\Gamma }}}_{1n}\\ 0& {\tilde{\mathbf{\Gamma }}}_{22}& \cdots & {\widetilde{ \mathbf{\Gamma }}}_{2n}\\ ⋮& ⋮& \cdots & ⋮\\ 0& 0& \cdots & {\widetilde{ \mathbf{\Gamma }}}_{nn}\end{array}\right]\widehat{Y}$

19

where $\widetilde{ \mathbf{\Gamma } }\hspace{0.17em}\in {ℝ}^{nN\hspace{0.17em}\times \hspace{0.17em}12n}$ is again a block upper triangular matrix, whose (i, j) th block component, ${\widetilde{\mathbf{\Gamma } }}_{ij}\in {ℝ}^{N\hspace{0.17em}\times \hspace{0.17em}12}$ , is made by stacking up $\text{row} ({\mathbf{\Gamma }}_j(k),i),k=1,2,\cdots ,N$ . The cost function in equation (18) becomes equivalent to $J=1/2E^{T}E$ .

Note that some of the dynamic parameters or some of their linear combinations are destined to be unobservable through the dynamic response due to the inherent structure of any given robot manipulator. Such redundancy of the dynamic parameters makes $\widetilde{ \mathbf{\Gamma } }$ in equation (19) rank-deficient, meaning that infinitely many $\widehat{Y}$ ’s can possibly be the solution of equation (18). Among those possible solutions, to obtain a physically meaningful acceptable $\widehat{Y}$ , we need to impose constraints on $\widehat{Y}$ such as the positivity of inertia and mass or even some bounds of parameters if we have guesses on them, as many previous studies did. By taking these constraints into account, the minimization problem can be reformulated as an SDP

$\begin{array}{ll}\hspace{0.17em}\hfill & \underset{\widehat{Y}}{{\displaystyle \text{min}\hspace{1em}}}u\hfill \\ \text{s.t.} \hfill & ||E{||}^2\le u,D(\widehat{Y})\succcurlyeq 0\hfill \end{array}$

20

where u is a positive optimization variable, and $D(\widehat{Y})$ is a set of all feasibility conditions to meet the physical conditions (refer to the works of Sousa and Cortesão ⁷ and Boyd and Vandenberghe ²² for introduction to SDP in the parameter identification).

However, a drawback associated with equation (20) is that it would take too much time to solve the problem if the numbers of DoFs (n) and measurements (N) are large. The computation time tends to grow geometrically as they increase. Moreover, if no solution is found after a long computation, then we have to resolve the SDP by relaxing conditions of some constraints. After all, the parameter identification of a large DoF robot manipulator via SDP (equation (20)) would not be completed within a reasonable time duration in practice.

Hence, it seems more reasonable to divide the big-sized SDP in equation (20) into smaller SDPs by taking advantage of the block upper triangular form of $\widetilde{ \mathbf{\Gamma } }$ and reduce the overall computation time. That is, first ${\widehat{Y}}_n$ is independently obtained by taking into account the last N equations of equation (19), then is ${\widehat{Y}}_{n\hspace{0.17em}-\hspace{0.17em}1}$ after combining the previously determined ${\widehat{Y}}_n$ , and so on. By repeatedly solving n times of small SDPs, $\widehat{Y}$ can be completely determined. The idea is similar to the back-substitution to solve a set of linear equations after Gaussian elimination. ²³ Advantage is that the time taken for parameter identification grows only linearly with the number of DoFs in theory, although some variation of computation time may occur since the convergence time reaching to solution is not deterministic in the SDP.

In detail, the i th SDP formulation under this sequential approach can be written as

$\begin{array}{l}\underset{{\widehat{Y}}_i}{{\displaystyle \text{min}\hspace{1em}}}{\lambda }_i\hfill \\ \text{s.t.} \hfill & \left[\begin{array}{cc}I& {\widetilde{ \mathbf{\Gamma } }}_{ii}{\widehat{Y}}_i-(T_i{)}^{\prime }\\ {\widehat{Y}}_i^T{\widetilde{ \mathbf{\Gamma } }}_{ii}^T-(T_i^T{)}^{\prime }& {\lambda }_i\end{array}\right]\succcurlyeq 0\text{,} D_i({\widehat{Y}}_i)\succcurlyeq 0\hfill \end{array}$

21

where ${\lambda }_i>0$ is a positive optimization variable and $(T_i{)}^{\prime }:=T_i\hspace{0.17em}-\hspace{0.17em}{\displaystyle {\sum }_{j\hspace{0.17em}=\hspace{0.17em}i\hspace{0.17em}+\hspace{0.17em}1}^n}{\widetilde{ \mathbf{\Gamma } }}_{ij}{\widehat{Y}}_j$ . The first linear matrix inequality is equivalent to $||{\widetilde{ \mathbf{\Gamma } }}_{ii}{\widehat{Y}}_i-(T_i{)}^{\prime }{||}^2\le {\lambda }_i$ , and $D_i({\widehat{Y}}_i)\succcurlyeq 0$ is the set of feasibility conditions associated with ${\widehat{Y}}_i$ . In equation (21), $(T_i{)}^{\prime }$ is a given term since ${\widehat{Y}}_j\text{,}\hspace{1em}j\hspace{0.17em}=\hspace{0.17em}i\hspace{0.17em}+\hspace{0.17em}1,i\hspace{0.17em}+\hspace{0.17em}2,\dots ,n$ are implicitly assumed to be known during prior steps. Eventually solving equation (21) by n times will complete the identification of dynamics parameters.

Case examples

To demonstrate how effective the proposed modeling and identification approach, we carried out simulation study by using two example systems, one from the conventional model of industrial manipulators and the other from a specially designed 9-DoF robot manipulator, as shown in Figures 3 and 5, respectively. The simulation was run under the MATLAB software (MathWorks, Inc., Natick, MA, USA) on a dedicated PC (CPU: Intel Core i5) with no other tasks running simultaneously.

Figure 3.

A PUMA industrial robot manipulator and its joint trajectory. PUMA: Programmable Universal Machine for Assembly. (a) An industrial robot manipulator; (b) Joint trajectory for data collection.

Example 1: Industrial manipulator

A 6-DoF PUMA-type industrial robot manipulator was used in the simulation, where only the first 3-DoF were included in the dynamic model after neglecting the DoFs in the wrist which were not significant to the overall dynamics. The dynamic simulation was done by using the summation form of Newton–Euler equation in equation (8). A friction model in equation (15) was inserted in the dynamic simulation to mimic the frictional effect. Gaussian random noise and high-frequency sinusoidal disturbance (A sin(100t)) were added to the input torque to enhance the realism of simulation.

Then, the regressor $\widetilde{ \mathbf{\Gamma }}$ in equation (19) was constructed by stacking up N = 2000 samples from the set of data for 5 s, resulting in a block upper triangular form. The next step was to prepare $D_i({\widehat{Y}}_i)$ to solve equation (21). For the condition that $D_i({\widehat{Y}}_i)\succcurlyeq 0,$ the following inequalities were concatenated in series

$\left[\begin{array}{cc}{\hspace{0.17em}}^t {\widehat{\mathbb{I}}}_i^i& S^T({\widehat{m}}_i{\widehat{\hspace{0.17em}r}}_i^i)\\ S({\widehat{m}}_i{\widehat{\hspace{0.17em}r}}_i^i)& {\widehat{m}}_{i}I\end{array}\right]\succcurlyeq 0 \text{and} {\underset{\_}{\eta }}_{ij}\le {\widehat{Y}}_{ij}\le {\overline{\eta }}_{ij}$

22

where ${\underset{\_}{\eta }}_{ij}$ and ${\overline{\eta }}_{ij}$ denote the lower and upper guessed bounds of ${\widehat{Y}}_{ij}$ . The matrix inequality in equation (22) implies the simultaneous positivity of ${\widehat{\mathbb{I}}}_i^i$ and ${\widehat{m}}_i$ from the Schur complement. ²⁴

With ${\underset{\_}{\eta }}_{ij}$ and ${\overline{\eta }}_{ij}$ being the values deviated by ±50% from their true values, we solved equation (21) from i = 3 to i = 1 using YALMIP foundation ²⁵ with SeDuMi SDP solver ²⁶ under MATLAB. As summarized in Table 1, the results of the non-minimal parameter identification were fairly accurate to the true parameters; some noticeable estimation errors could be observed in the parameters of the first link, but it was due to the non-observability of a few parameters in the regression equation and did not affect the overall dynamics. The computation time was about 51 min, although it slightly varied depending upon cases. For comparison, we solved the same problem using equation (20) where all parameters and constraints were simultaneously taken into account for. The identification results from the sequential and the simultaneous approaches were almost the same as shown in Table 1, but the simultaneous approach took nearly 10 h for completion, which was almost 11 times longer than that by using the proposed sequential approach.

Table 1.

Comparison of identified parameters of 3-DoF robot.^a

		m	rx	ry	rz	Ixx	Iyy	Izz	Ixy	Ixz	Iyz	α	β
Joint 1	True	10	−0.1	−0.05	−0.05	0.5	0.6	0.55	0.2	0.3	0.25	3	1.2
	Estimated 1	12.3302	−0.0999	−0.0500	−0.0500	0.5669	0.5336	0.6270	0.1959	0.2818	0.2403	3.0194	1.1673
	Estimated 2	12.3662	−0.0999	−0.0500	−0.0500	0.5679	0.6210	0.6275	0.1965	0.2811	0.2416	3.0111	1.1963
Joint 2	True	8	−0.15	−0.05	−0.05	0.3	0.2	0.36	0.14	0.24	0.14	2.5	0.8
	Estimated 1	8.0690	−0.1556	−0.0482	−0.0503	0.3400	0.2359	0.3633	0.1435	0.2296	0.1643	2.4546	0.8162
	Estimated 2	8.0443	−0.1502	−0.0421	−0.0500	0.2871	0.1981	0.3603	0.1478	0.2392	0.1514	2.4731	0.8087
Joint 3	True	6	−0.18	−0.01	0.02	0.15	0.15	0.30	0.05	0.12	0.06	2	0.6
	Estimated 1	6.0138	−0.1935	−0.0102	0.0172	0.1736	0.1554	0.3325	0.0618	0.1191	0.0641	1.9765	0.5658
	Estimated 2	5.9870	−0.1809	−0.0100	0.0253	0.1476	0.1504	0.2921	0.0517	0.1221	0.0604	2.0292	0.5479

DoF: degree-of-freedom.

^aTrue: true values; estimated 1: the proposed sequential approach; and estimated 2: the simultaneous approach.

The input torques and the estimated torques from the proposed sequential approach were compared in Figure 4. As shown, the estimated torques faithfully followed the actual input torques, even with noises and disturbance added in the input torques. This verifies that the sequential SDP was solved correctly.

Figure 4.

Comparison of the input torque and the reconstructed torque after optimization for an industrial robot manipulator. (a) τ1; (b) τ2; (c) τ3.

Example 2: 9-DoF robot manipulator

As the second example, a 9-DoF robot manipulator, as shown in Figure 5, was employed to verify the validity of the proposed approach in high-DoF systems. Since every link possesses 12 dynamic parameters, total 108 parameters had to be identified in this case. Like in the previous example, the recursive form of dynamics in equation (8) was reused for simulation. The Gaussian noises and sinusoidal disturbances were also included during the simulation. After collecting input torques and motion data, we formed the block upper triangular matrix $\widetilde{ \mathbf{\Gamma } }$ with N = 1000 samples. The optimization problem in equation (21) was then solved with feasibility constraints (equation (22)) in the backward direction, that is, from i = 9 to i = 1.

Figure 5.

A 9-DoF robot manipulator and its joint trajectory. DoF: degree-of-freedom. (a) A 9-DoF robot manipulator, (b) Joint trajectory for data collection.

Table 2 shows the numerical values of the true and the identified parameters by the proposed sequential method. The identified parameters all satisfied the feasibility conditions specified by inequalities. The computation time taken by the proposed method was approximately 20 min, a level that was slightly higher than our linearity expectation from the previous 3-DoF case. Overall, most of the parameters were identified close to their true values, but the identification error became larger, when compared to the previous 3-DoF example. This was because some parameters of the 9-DoF robot manipulator were too small or less involved in the motion of the given joint trajectory to be identified exactly; put in different respect, the trajectory seemed rather monotonous and its time duration was short to fully excite the dynamics of the whole nine links. With more data from different joint motions, we would achieve better identification result at the cost of increased computational burden.

Table 2.

Comparison of identified parameters of 9-DoF robot.^a

		m	rx	ry	rz	Ixx	Iyy	Izz	Ixy	Ixz	Iyz	α	β
Joint 1	True	12.9000	−0.0011	−0.0812	−0.0058	0.2327	0.2188	0.1820	−0.0019	0.0001	0.0083	3	1.2
	Estimated 1	15.8876	−0.0011	−0.0811	−0.0058	0.2573	0.1912	0.2012	−0.0019	0.0001	0.0083	3.9962	0.9665
	Estimated 2	15.9449	−0.0011	−0.0811	−0.0058	0.2578	0.1893	0.2016	−0.0019	0.0001	0.0083	2.7622	1.2387
Joint 2	True	14.0503	−0.2356	0.0007	−0.1975	0.0830	0.4836	0.4884	−0.0020	0.0243	0.0001	3	1.2
	Estimated 1	15.9405	−0.1942	0.0011	−0.2962	0.0768	0.6873	0.4195	−0.0030	0.0365	0.0002	4.5000	0.6203
	Estimated 2	14.6928	−0.2431	0.0004	−0.1901	0.0861	0.4851	0.4721	−0.0010	0.0232	0.0002	2.7262	1.3655
Joint 3	True	2.9731	−0.0158	0.0404	−0.0177	0.0131	0.0089	0.0135	0.0020	0.0006	−0.0020	3	1.2
	Estimated 1	2.4866	−0.0237	0.0506	−0.0189	0.0066	0.0045	0.0203	0.0010	0.0003	−0.0010	2.4530	1.2283
	Estimated 2	3.1513	−0.0237	0.0202	−0.0266	0.0197	0.0045	0.0068	0.0030	0.0003	−0.0010	3.2119	1.1004
Joint 4	True	1.6650	0.0166	−0.0305	−0.0348	0.0071	0.0074	0.0047	0.0006	−0.0011	0.0017	2.5	0.8
	Estimated 1	2.4974	0.0250	−0.0152	−0.0174	0.0036	0.0037	0.0070	0.0003	−0.0005	0.0009	2.4859	1.2000
	Estimated 2	1.9855	0.0250	−0.0237	−0.0248	0.0036	0.0111	0.0070	0.0009	−0.0005	0.0026	2.7793	0.7553
Joint 5	True	2.6874	0.0187	−0.0577	−0.0167	0.0155	0.0068	0.0154	0.0022	−0.0007	0.0024	2.5	0.8
	Estimated 1	3.3437	0.0194	−0.0300	−0.0800	0.0103	0.0034	0.0137	0.0011	−0.0003	0.0036	2.7557	0.9955
	Estimated 2	3.3586	0.0203	−0.0865	−0.0148	0.0093	0.0034	0.0104	0.0033	−0.0003	0.0036	2.6040	0.7507
Joint 6	True	1.2656	0.0170	−0.0271	−0.0287	0.0039	0.0040	0.0030	0.0005	−0.0007	0.0010	2.5	0.8
	Estimated 1	0.6328	0.0085	−0.0135	−0.0431	0.0036	0.0042	0.0045	0.0002	−0.0004	0.0005	1.2500	0.9250
	Estimated 2	0.6328	0.0085	−0.0151	−0.0144	0.0059	0.0020	0.0028	0.0007	−0.0004	0.0015	2.5803	0.8189
Joint 7	True	2.2301	0.0000	0.0662	−0.0136	0.0137	0.0042	0.0130	−0.0000	−0.0000	−0.0018	2	0.6
	Estimated 1	2.6900	0.0000	0.0617	−0.0163	0.0144	0.0038	0.0154	−0.0000	−0.0000	−0.0027	1.9556	0.6300
	Estimated 2	1.7765	0.0000	0.0773	−0.0068	0.0069	0.0063	0.0130	−0.0000	−0.0000	−0.0009	1.9556	0.6300
Joint 8	True	1.6912	0.0412	−0.0213	−0.0198	0.0040	0.0067	0.0065	0.0017	−0.0011	0.0006	2	0.6
	Estimated 1	1.6857	0.0337	−0.0229	−0.0248	0.0038	0.0053	0.0086	0.0017	−0.0016	0.0003	1.9982	0.5991
	Estimated 2	2.5244	0.0396	−0.0253	−0.0099	0.0059	0.0035	0.0066	0.0008	−0.0017	0.0003	1.6126	0.4580
Joint 9	True	0.5510	−0.0180	0.0028	0.0434	0.0061	0.0064	0.0090	−0.0000	−0.0000	0.0000	2	0.6
	Estimated 1	0.6267	−0.0212	0.0028	0.0475	0.0073	0.0077	0.0087	−0.0000	−0.0000	0.0000	2.0013	0.6061
	Estimated 2	0.6256	−0.0197	0.0028	0.0414	0.0039	0.0071	0.0071	−0.0000	−0.0000	0.0000	2.1847	0.5681

DoF: degree-of-freedom.

^aTrue: true values; estimated 1: the proposed sequential approach; and estimated 2: the simultaneous approach.

For the purpose of comparison, when we tried to solve equation (20) by taking all parameters into account simultaneously, it took more than 4 days to obtain the solution. This shows that the simultaneous SDP tends to tremendously increase the computation time as the number of DoFs of the robot manipulator becomes larger, whereas the computation time of the proposed sequential method increases just reasonably. The numerical values from this simultaneous SDP, called as “estimated 2” in Table 2, looked close to those from the proposed sequential SDP.

The profiles of the actual input torques and the estimated torques by the proposed approach, as shown in Figure 6, were very close. These results supported the validity and effectiveness of the proposed approach.

Figure 6.

Comparison of the input torques and the estimated torques for a 9-DoF robot manipulator. DoF: degree-of-freedom (a) τ1; (b) τ2; (c) τ3; (d) τ4; (e) τ5; (f) τ6; (g) τ7; (h) τ8; (i) τ9.

Concluding remarks

In this article, a sequential computational algorithm was proposed to solve the optimization problem for parameter identification of robot manipulators that may have a large number of links. We showed how to construct the regressor equation from a succinctly arranged form of robot dynamic equation. The associated regressor matrix was bounded to be of an upper triangular matrix, which allowed us to solve the linear optimization by the sequential backward manner. We combined inequality conditions, while formulating the optimization, as needed so that the identified parameters have valid physical meaning. Two numerical examples showed the validity and effectiveness of the proposed algorithm. Since the proposed algorithm is general and scalable, we believe that a computer code for the algorithm, once programmed, can be reused to any kind of serial-type robot manipulator.