Parameter Estimation for Robotic Manipulator Systems ^†

Abstract

In this paper, a novel methodology for estimating the parameters of robotic manipulator systems is proposed. It can be seen that, for the purpose of parameter estimation, the input torque to each joint motor is designed as a linear combination of sinusoids. After the transient responses of joint angles exponentially converge to zero, the steady states of joint angle outputs can be extracted. Since the steady states of joint angles are the equivalent finite Fourier series, the coefficients of the steady state components of joint angles can be further extracted in a fundamental period. With the amazing finding that the steady states contain all dynamic information of manipulator systems, all unknown parameters of the system model can be accurately estimated with the extracted coefficients in finite frequency bands. The simulation results for a two-link manipulator are carried out to illustrate the effectiveness and robustness against measurement noise of the proposed method.

1. Introduction

In modern manufacturing industry, robot manipulator systems are widely adopted in many aspects such as assembling, welding, painting, etc. [1,2,3]. In order to ensure the stability and reliable performance of robot manipulators, many model-based controllers have been proposed [4,5,6,7,8,9,10,11]. The stability of the control performance highly depends on the accuracy of the system dynamic model. To establish the accurate dynamic model, the knowledge of system dynamic parameters is required [12]. However, internal parameter perturbation may occur during long-time continuous operation in many practical applications, which may degrade the reliability of the system operations. Thus, the accurate estimation of system dynamic parameters is of great significance.

Various methodologies have been developed to solve the parameter estimation problem for manipulator systems [13,14,15,16,17,18,19]. One commonly adopted method for identifying the unknown manipulator parameters is the recursive least square method (RLS) [20,21]. The common characteristics of the existing RLS methods are that a cost function needs to be designed; the optimization process will start at a random initial position in the parameter space and proceed along with the surface of that cost function in a series manner [22]. Although many RLS-based approaches have been proposed to complete variable parameter estimation tasks, these methods still have several limitations. For instance, most RLS methods in the time domain have always suffered from low convergence speed [23]; the robustness of the RLS is weak against external disturbances due to its low sensitivity to the new measurements in long-time recursive iteration [24]. Moreover, the local minima problem will happen if the start point has not been selected properly, which may also lead to the degeneration of the estimation performance [25,26].

The Kalman filter (KF) and its variants including extended the Kalman filter (EKF) and unscented the Kalman filter (UKF) have also been proved to be effective of estimating unknown dynamic parameters of manipulator systems [27,28,29,30,31,32]. The KF algorithms employ the stochastic system model, and the model uncertainties and measurement noises are all assumed to be Gaussian. With properly designed KF gains, the posteriori error covariance is minimized, and the optimal estimates of model parameters are obtained in the form of a weighted combination of the priori estimates and measurements. However, the implementation of KF algorithms requires a nominal system model and full system states information [33]. The estimation accuracy will also be degraded due to the unknown distribution of disturbance [34]. In manipulator systems, if the joint angular velocity and angular acceleration are not available, it will be quite hard to formulate the KF algorithm for either estimating system model parameters or system states [35,36]. Furthermore, the robustness of the Kalman filter cannot be guaranteed as well [37,38].

In this paper, a novel parameter estimation scheme is developed for manipulator systems in the frequency domain. To start with, the manipulator system dynamic model has been re-constructed using the linearization technique. It is shown that the information from the other joint can be decoupled and only the joint angular position measurement is required for further parameter estimation. Considering the fact that the transient response of the joint angular position and angular velocity are all damped and will exponentially converge to zero with fast convergence speed, the angular position and angular velocity are dominated by their steady state after sufficient long time. The steady state of angular position can then be extracted from the measurement. In this study, the manipulator system is excited by a linear combination of sinusoidal components. The steady state of the output can also be equivalent to a finite Fourie series. It is shown that, the coefficients of the steady state joint position measurement are all polynomials about the system dynamic parameters and system’s dynamic information has been fully embedded in these coefficients within one fundamental period. Then, the coefficients of the steady state of the angular position can be extracted from the angular position measurement using joint position measurement in one fundamental period. Furthermore, utilizing the relationship among the coefficients of the steady state sinusoidal components of joint position measurement and the system’s dynamic parameters in the frequency domain, all the unknown parameters including the mass of the joint, the length of the joint, etc. can then be estimated based on extracted coefficients.

The advantages of the proposed method can be concluded as: (i) only the steady state part of the measurement in one fundamental period is required to extract the coefficients of the steady state sinusoidal components. (ii) All unknown parameters within the system model can be accurately estimated based on extracted coefficients while the estimation convergence speed is faster than that of traditional recursive estimation schemes and the global optimization mechanism is guaranteed. (iii) The averaging process of extracting the coefficients of steady state sinusoidal components contributes to compensate the effects of the measurement noises. The proposed estimation scheme is robustness against measurements.

Compared to existing parameter estimation approaches, in the proposed method, a nominal model and prior information are not required, only the steady state of the joint position measurement in one fundamental period is needed; the coefficients of the steady state sinusoidal components can all be extracted, based on which all unknown parameters can be estimated in finite frequency bands. In addition, due to the orthogonal properties of the trigonometric base in the sinusoidal components, the computational cost of matrix inversion of the high dimension can be greatly reduced; therefore, the proposed estimation method exhibits high convergence speed, since no recursive process is involved. Moreover, the proposed estimation scheme will show a global optimization mechanism.

The rest of this paper is organized as follows: in Section 2, the dynamic model of a manipulator is reconstructed using a linearization technique. The system’s output response in frequency domain is analyzed. In Section 3, the input torque is designed as a linear combination of sinusoids and the coefficients of steady state sinusoidal components of the measurement are derived and extracted. The unknown parameters are then estimated based on these coefficients. In Section 4, the simulation results are carried out to verify the theoretical analysis and good estimation performance of the proposed method; the robustness of the proposed method is also discussed. In Section 5, the conclusion is briefly summarized.

2. Problem Formulation

Consider a two-link manipulator; the angular position and angular velocity are denoted as $q={\left[q_1,q_2\right]}^T$ and $\dot{q}={\left[{\dot{q}}_1,{\dot{q}}_2\right]}^T,$ respectively. The dynamic model of this 2 degrees of freedom (DOF) manipulator can be derived as:

$$D\left(q\right)\ddot{q}+C\left(q,\dot{q}\right)\dot{q}+G\left(q\right)=\tau$$

(1)

where $\tau$ is the vector of torques applied at two joints and

$$D\left(q\right)=\left[\begin{array}{cc}{d}_{11}& d_{12}\\ d_{21}& d_{22} \end{array}\right]$$

(2)

$$\{\begin{array}{l}{d}_{11}=\left(m_1+m_2\right)r_1^2+m_2{r}_2^2+2m_2{r}_1{r}_2\cos \left(q_2\right)+J_1\hfill \\ d_{12}=d_{21}=m_2{r}_2^2+m_2{r}_1{r}_2\cos \left(q_2\right)\hfill \\ d_{22}=m_2{r}_2^2+J_2\hfill \end{array}$$

(3)

$$C\left(q,\dot{q}\right)\dot{q}=\left[\begin{array}{c}{\beta }_{12}\left(q_2\right){\dot{q}}_1^2+2{\beta }_{12}\left(q_2\right){\dot{q}}_1{\dot{q}}_2\\ -{\beta }_{12}\left(q_2\right){\dot{q}}_2^2\end{array}\right]$$

(4)

$${\beta }_{12}\left(q_2\right)=m_2{r}_1{r}_2\sin \left(q_2\right)$$

(5)

$$G\left(q\right)=\left[\begin{array}{c}{\gamma }_1\left(q_1,q_2\right)g\\ {\gamma }_2\left(q_1,q_2\right)g\end{array}\right]$$

(6)

$${\gamma }_1\left(q_1,q_2\right)=-\left(\left(m_1+m_2\right)r_1\cos \left(q_2\right)+m_2{r}_2\cos \left(q_1+q_2\right)\right)$$

(7)

where $m_1$ and $m_2$ are the masses of link 1 and link 2. $l_1$ and $l_2$ are the length of link 1 and link 2, respectively. Due to the manufacturing error, the actual values of these dynamic quantities are different from their nominal values. In addition, the parameter perturbation during the operation also degenerates the reliability of the model. Thus, it is of great significance to develop an online parameter estimation scheme to estimate real time values of all robotic mechanical quantities. The system’s model (1) can be further re-written as:

$$\ddot{q}=\alpha \left(q,\dot{q}\right)+\beta \left(q\right)\tau =f\left(q,\tau \right)$$

(8)

where

$$\alpha \left(q,\dot{q}\right)=D^{-1}\left(q\right)\left(-C\left(q,\dot{q}\right)\dot{q}-G\left(q\right)\right)$$

(9)

$$\left(q\right)=D^{-1}\left(q\right)$$

(10)

In order to implement the proposed parameter identification algorithm, the Taylor series expansion method has been applied to system (1) for the approximation of a linear system, considering the fact that the manipulator system can be stabilized, and the manipulator joint can remain at an equilibrium position $q_0=\left[q_{10},q_{20}\right]$. Meanwhile, the system’s joint angular position and velocity are all assumed to change slowly. Thus, the system can be linearized regarding to the equilibrium position $q_0$ as:

$${\ddot{q}}_0+\delta \ddot{q}=f\left(q_0+\delta q,{\tau }_0+\delta \tau \right)\approx f\left(q_0,{\tau }_0\right)+\frac{\partial f}{\partial q}{|}_{q_0,{\tau }_0}\delta q+\frac{\partial f}{\partial \tau }{|}_{q_0,{\tau }_0}\delta \tau$$

(11)

where ${\ddot{q}}_0=f\left(q_0,{\tau }_0\right)$; thus,

$$\delta \ddot{q}=\frac{\partial f}{\partial q}{|}_{q_0,{\tau }_0}\delta q+\frac{\partial f}{\partial \tau }{|}_{q_0,{\tau }_0}\delta \tau$$

(12)

The linearized system model (12) in the time domain can be further re-expressed in state-space form as:

$$\left[\begin{array}{c}\delta {\ddot{q}}_1\\ \delta {\ddot{q}}_2\end{array}\right]=\left[\begin{array}{cc}{a}_{11}& a_{12}\\ a_{21}& a_{22}\end{array}\right]\left[\begin{array}{c}\delta q_1\\ \delta q_2\end{array}\right]+\left[\begin{array}{cc}{b}_{11}& b_{12}\\ b_{21}& b_{22}\end{array}\right]\left[\begin{array}{c}\delta {\tau }_1\\ \delta {\tau }_2\end{array}\right]$$

(13)

where

$$\{\begin{array}{l}{a}_{11}=\frac{g\ast \left(2m_1\sin \left(q_{10}\right)\right)+m_2\sin \left(q_{10}\right)-m_2\sin \left(q_{10}+2q_{20}\right)}{l_1\left(2m_1+m_2-m_2\cos \left(2q_{20}\right)\right)}\hfill \\ a_{12}=-\frac{g{m}_2\left(\sin \left(q_{10}+2q_{20}\right)+\sin \left(q_{10}\right)\right)}{l_1\left(2m_1+m_2-m_2\cos \left(2q_{20}\right)\right)}\hfill \\ a_{21}=\frac{g\left(\begin{array}{c}{l}_2{m}_2{\cos }^2\left(q_{20}\right)\sin \left(q_{10}\right)-l_2{m}_2\sin \left(q_{10}\right)-\\ l_2{m}_1\sin \left(q_{10}\right)+l_1{m}_1\cos \left(q_{10}\right)\sin \left(q_{20}\right)+\\ l_2{m}_2\cos \left(q_{10}\right)\cos \left(q_{20}\right)\sin \left(q_{20}\right)\end{array}\right)}{l_1{l}_2\left(m_1+m_2-m_2{\cos }^2\left(q_{20}\right)\right)}\hfill \\ a_{22}=\frac{l_1^2{m}_1+l_1^2{m}_2+2l_1{l}_2{m}_2\cos \left(q_{20}\right)}{l_1^2{l}_2^2{m}_2\left(m_1+m_2-m_2{\cos }^2\left(q_{20}\right)\right)}\hfill \end{array}$$

(14)

$$\{\begin{array}{l}\begin{array}{c}{b}_{11}=\frac{1}{l_1^2\left(m_1+m_2-m_2{\cos }^2\left(q_{20}\right)\right)}\\ b_{12}=-\frac{l_2+l_1\cos \left(q_{20}\right)}{l_1^2{l}_2\left(m_1+m_2-m_2{\cos }^2\left(q_{20}\right)\right)}\end{array}\hfill \\ \begin{array}{c}{b}_{21}=b_{12}\\ b_{22}=\frac{l_1^2{m}_1+l_1^2{m}_2+l_2^2{m}_2+2l_1{l}_2{m}_2\cos \left(q_{20}\right)}{l_1^2{l}_2^2{m}_2\left(m_1+m_2-m_2{\cos }^2\left(q_{20}\right)\right)}\end{array}\hfill \end{array}$$

(15)

The transfer functions of system (13) can then be formulated as:

$$\begin{array}{c}{Q}_1\left(s\right)=\frac{b_{11}s+a_{12}{b}_{21}-a_{22}{b}_{11}}{s^2-\left(a_{11}+a_{22}\right)s+a_{11}{a}_{22}-a_{12}{a}_{21}}{U}_1\left(s\right)+\frac{b_{12}s+a_{12}{b}_{22}-a_{22}{b}_{12}}{s^2-\left(a_{11}+a_{22}\right)s+a_{11}{a}_{22}-a_{12}{a}_{21}}{U}_2\left(s\right)\\ =\frac{\left(b_{11}s+p_1\right)U_1\left(s\right)}{{\left(s-{\sigma }_0\right)}^2+{\omega }_d^2}+\frac{\left(b_{12}s+p_2\right)U_2\left(s\right)}{[{\left(s-{\sigma }_0\right)}^2+{\omega }_d^2]}\end{array}$$

(16)

$$\begin{array}{c}{Q}_2\left(s\right)=\frac{b_{21}s+a_{21}{b}_{11}-a_{11}{b}_{21}}{s^2-\left(a_{11}+a_{22}\right)s+a_{11}{a}_{22}-a_{12}{a}_{21}}{U}_1\left(s\right)+\frac{b_{22}s+a_{21}{b}_{12}-a_{11}{b}_{22}}{s^2-\left(a_{11}+a_{22}\right)s+a_{11}{a}_{22}-a_{12}{a}_{21}}{U}_2\left(s\right)\\ =\frac{\left(b_{21}s+p_3\right)U_1\left(s\right)}{{\left(s-{\sigma }_0\right)}^2+{\omega }_d^2}+\frac{\left(b_{22}s+p_4\right)U_2\left(s\right)}{[{\left(s-{\sigma }_0\right)}^2+{\omega }_d^2]}\end{array}$$

(17)

where ${\sigma }_0=\frac{a_{11}+a_{22}}{2}$, ${\omega }_d^2=\left(a_{11}{a}_{22}-a_{12}{a}_{21}\right)-{\frac{a_{11}+a_{22}}{2}}^2={\omega }_0^2-{\sigma }_0^2$. $p_1=a_{12}{b}_{21}-a_{22}{b}_{11}$, $p_2=a_{12}{b}_{22}-a_{22}{b}_{12}$, $p_3=a_{21}{b}_{11}-a_{11}{b}_{21}$, $p_4=a_{21}{b}_{12}-a_{11}{b}_{22}$ In this case, $|{\sigma }_0|$ is the damping factor, ${\omega }_0$ is the undamped frequency, and ${\omega }_d$ is the damping frequency. In this section, for convenience, it is assumed that:

$${\omega }_0^2-{\sigma }_0^2\ge 0$$

(18)

In systems (16) and (17), parameters $a_{11}$− $a_{22}$, $b_{11 }$ − $b_{22}$ are all polynomials with respect to the dynamic quantities. In order to estimate the dynamic quantities of manipulator systems, system parameters need to be obtained first. In the next section, the parameter estimation algorithm is proposed to extract the parameters from the steady state of the joint position output.

3. Main Results

Regarding the system’s transfer function (16) and (17), the input signal has been designed as a linear combination of sinusoidal components as:

$$u_k\left(t\right)={{\displaystyle \sum }}_{i=1}^n{d}_{ki}\cos \left({\omega }_{ki}t\right),k=1,2$$

(19)

The corresponding transfer function of the above input signal is given by:

$$U_k\left(s\right)={{\displaystyle \sum }}_{i=1}^n\frac{d_{ki}s}{s^2+{\omega }_{ki}^2}, k=1,2$$

(20)

Substituting (20) into (16) and (17) yields:

$$Q_1\left(s\right)={{\displaystyle \sum }}_{i=1}^n(\frac{\left(b_{11}s+p_1\right)d_{1i}s}{[{\left(s-{\sigma }_0\right)}^2+{\omega }_d^2](s^2+{\omega }_{1i}^2)}+\frac{\left(b_{12}s+p_2\right)d_{2i}s}{[{\left(s-{\sigma }_0\right)}^2+{\omega }_d^2](s^2+{\omega }_{2i}^2)})$$

(21)

$$Q_2\left(s\right)={{\displaystyle \sum }}_{i=1}^n(\frac{\left(b_{21}s+p_3\right)d_{1i}s}{[{\left(s-{\sigma }_0\right)}^2+{\omega }_d^2](s^2+{\omega }_{1i}^2)}+\frac{\left(b_{22}s+p_4\right)d_{2i}s}{[{\left(s-{\sigma }_0\right)}^2+{\omega }_d^2](s^2+{\omega }_{2i}^2)})$$

(22)

3.1. Estimation of Steady State Coefficients

Using the partial fraction expansion method, the $i_{\mathrm{th}}$ component of (21) can be expressed as:

$$Q_{1i}\left(s\right)=\frac{A_{11i}\left(s-{\sigma }_0\right)+B_{11i}{\omega }_d}{{\left(s-{\sigma }_0\right)}^2+{\omega }_d^2}+\frac{C_{11i}s+D_{11i}{\omega }_{1i}}{s^2+{\omega }_{1i}^2}+\frac{A_{12i}\left(s-{\sigma }_0\right)+B_{12i}{\omega }_d}{{\left(s-{\sigma }_0\right)}^2+{\omega }_d^2}+\frac{C_{12i}s+D_{12i}{\omega }_{2i}}{s^2+{\omega }_{2i}^2}$$

(23)

Based on (21) and (23), the following relationships can be obtained:

$$\{\begin{array}{l}\begin{array}{l}{A}_{11i}+C_{11i}=0\hfill \\ -A_{11i}{\sigma }_0+B_{11i}{\omega }_d-2C_{11i}{\sigma }_0+D_{11i}{\omega }_{1i}=d_{1i}{b}_{11} \hfill \end{array}\hfill \\ \begin{array}{l}{A}_{11i}{\omega }_{1i}^2+C_{11i}{\omega }_0^2-2D_{11i}{\omega }_{ii}{\sigma }_0=d_{1i}{p}_1\hfill \\ -A_{11i}{\sigma }_0{\omega }_{11i}^2+B_{11i}{\omega }_d{\omega }_{11i}^2+D_{11i}{\omega }_{11i}{\omega }_0^2=0\hfill \end{array}\hfill \end{array}$$

(24)

By solving the equation set (24), $A_{11i}$, $B_{11i}$, $C_{11i}$ and $D_{11i}$ can be calculated as:

$$\{\begin{array}{l}{A}_{11i}=\frac{2{\omega }_{11i}{\sigma }_0{d}_{1i}{b}_{11}+d_{1i}{p}_1\left({\omega }_{1i}^2-{\omega }_0^2\right)}{4{\omega }_{1i}^2{\sigma }_0^2+{\left({\omega }_{1i}^2-{\omega }_0^2\right)}^2}\hfill \\ B_{11i}=\frac{(d_{1i}{b}_{21})}{{\omega }_d}+\frac{{\sigma }_0}{{\omega }_d}\frac{-2{\omega }_{11i}{\sigma }_0{d}_{1i}{b}_{11}-d_{1i}{p}_1\left({\omega }_{1i}^2-{\omega }_0^2\right)}{4{\omega }_{1i}^2{\sigma }_0^2+{\left({\omega }_{1i}^2-{\omega }_0^2\right)}^2}\hfill \\ +\frac{{\omega }_{1i}}{{\omega }_d}\frac{2{\omega }_{1i}{\sigma }_0{d}_{1i}{p}_1-d_{1i}{b}_{11}{\omega }_{1i}\left({\omega }_{1i}^2-{\omega }_0^2\right)}{4{\omega }_{1i}^2{\sigma }_0^2+{\left({\omega }_{1i}^2-{\omega }_0^2\right)}^2}\hfill \\ C_{11i}=\frac{-2{\omega }_{11i}{\sigma }_0{d}_{1i}{b}_{11}-d_{1i}{p}_1\left({\omega }_{1i}^2-{\omega }_0^2\right)}{4{\omega }_{1i}^2{\sigma }_0^2+{\left({\omega }_{1i}^2-{\omega }_0^2\right)}^2}\hfill \\ D_{11i}=-\frac{2{\omega }_{1i}{\sigma }_0{d}_{1i}{p}_1-d_{1i}{b}_{11}{\omega }_{1i}\left({\omega }_{1i}^2-{\omega }_0^2\right)}{4{\omega }_{1i}^2{\sigma }_0^2+{\left({\omega }_{1i}^2-{\omega }_0^2\right)}^2}\hfill \end{array}$$

(25)

Similarly, $A_{21i}$, $B_{21i}$, $C_{21i}$ and $D_{21i}$ can also be obtained as:

$$\{\begin{array}{l}{A}_{21i}=\frac{2{\omega }_{21i}{\sigma }_0{d}_{2i}{b}_{12}+d_{2i}{p}_2\left({\omega }_{2i}^2-{\omega }_0^2\right)}{4{\omega }_{2i}^2{\sigma }_0^2+{\left({\omega }_{2i}^2-{\omega }_0^2\right)}^2}\hfill \\ B_{21i}=\frac{(d_{2i}{b}_{22})}{{\omega }_d}+\frac{{\sigma }_0}{{\omega }_d}\frac{-2{\omega }_{21i}{\sigma }_0{d}_{2i}{b}_{12}-d_{2i}{p}_2\left({\omega }_{2i}^2-{\omega }_0^2\right)}{4{\omega }_{2i}^2{\sigma }_0^2+{\left({\omega }_{2i}^2-{\omega }_0^2\right)}^2}\hfill \\ +\frac{{\omega }_{2i}}{{\omega }_d}\frac{2{\omega }_{2i}{\sigma }_0{d}_{2i}{p}_2-d_{2i}{b}_{12}{\omega }_{2i}\left({\omega }_{2i}^2-{\omega }_0^2\right)}{4{\omega }_{2i}^2{\sigma }_0^2+{\left({\omega }_{2i}^2-{\omega }_0^2\right)}^2}\hfill \\ C_{21i}=\frac{-2{\omega }_{21i}{\sigma }_0{d}_{2i}{b}_{12}-d_{2i}{p}_2\left({\omega }_{2i}^2-{\omega }_0^2\right)}{4{\omega }_{2i}^2{\sigma }_0^2+{\left({\omega }_{2i}^2-{\omega }_0^2\right)}^2}\hfill \\ D_{21i}=-\frac{2{\omega }_{2i}{\sigma }_0{d}_{2i}{p}_2-d_{2i}{b}_{12}{\omega }_{2i}\left({\omega }_{2i}^2-{\omega }_0^2\right)}{4{\omega }_{2i}^2{\sigma }_0^2+{\left({\omega }_{2i}^2-{\omega }_0^2\right)}^2}\hfill \end{array}$$

(26)

Therefore, the $i_{th}$ component of the output joint angle (23) can be re-expressed in the time domain as:

$$q_{1i}\left(t\right)=q_{1ti}\left(t\right)+q_{1si}\left(t\right)\#\left(27\right)$$

(27)

where

$$\begin{array}{c}{q}_{1ti}\left(t\right)=A_{11i}{e}^{{\sigma }_{0}t}\cos \left({\omega }_{1i}t\right)+B_{11i}{e}^{{\sigma }_{0}t}\sin \left({\omega }_{1i}t\right)+A_{21i}{e}^{{\sigma }_{0}t}cos\left({\omega }_{2i}t\right)+B_{21i}{e}^{{\sigma }_{0}t}sin\left({\omega }_{2i}t\right)\end{array}$$

(28)

$$q_{1is}\left(t\right)=C_{11i}\cos \left({\omega }_{1i}t\right)+D_{11i}\sin \left({\omega }_{1i}t\right)+C_{12i}cos\left({\omega }_{2i}t\right)+D_{12i}\sin \left({\omega }_{2i}t\right)$$

(29)

where $q_{1ti}\left(t\right)$ is the transient component of the $i_{\mathrm{th}}$ term of joint position output while $q_{1is}\left(t\right)$ is the steady state component. As time $t$ moves towards infinity, with the excitation input signal (19) and the transient response component $q_{1ti}\left(t\right)$ will exponentially decrease to $0$. If the time constant $1/|{\sigma }_0|$ is small enough, the joint position output will be dominated by the steady state component. The steady state output can then be formulated as:

$$q_{1s}\left(t\right)={\displaystyle \sum }_{i=1}^n{C}_{11i}\cos \left({\omega }_{1i}t\right)+D_{11i}\sin \left({\omega }_{1i}t\right)+C_{12i}cos\left({\omega }_{2i}t\right)+D_{12i}\sin \left({\omega }_{2i}t\right)$$

(30)

To extract the coefficients $C_{11i}$, $D_{11i}$, $C_{12i}$ and $D_{12i}$ (30) can be expanded into the matrix from as:

$$q_{1s}\left(t\right)=[\cos \left({\omega }_{11}\right),\sin \left({\omega }_{11}\right),\cos \left({\omega }_{21}\right),\sin \left({\omega }_{21}\right),\cdots \cos \left({\omega }_{1n}\right),\sin \left({\omega }_{1n}\right),\cos \left({\omega }_{2n}\right),\sin \left({\omega }_{2n}\right)]\left[\begin{array}{c}\begin{array}{c}{C}_{11}\\ D_{11}\end{array}\\ ⋮\\ \begin{array}{c}{C}_{2n}\\ D_{2n}\end{array}\end{array}\right]$$

(31)

Taking $m$ samples of each function in (29) gives:

$${\mathit{Q}}_1=[q_{1s}\left(t-\left(m-1\right)\Delta \right),q_{1s}\left(t-\left(m-2\right)\Delta \right),\cdots {q_{1s}\left(t-\Delta \right),q_{1s}\left(t\right)]}^T$$

(32)

where $\Delta$ is the sampling period. Substituting (28) into (29) yields:

$$\begin{gathered} {\mathit{Q}}_{1 }={\mathit{A}}_1{\mathit{P}}_1=[\begin{array}{cc}\begin{array}{cc}\begin{array}{c}\begin{array}{c}\cos \left({\omega }_{11}\left(t-\left(m-1\right)\Delta \right)\right)\\ \cos \left({\omega }_{11}\left(t-\left(m-2\right)\Delta \right)\right)\end{array}\\ \begin{array}{c}⋮\\ \cos \left({\omega }_{11}t\right)\end{array}\end{array}& \begin{array}{c}\begin{array}{c}\sin \left({\omega }_{11}\left(t-\left(m-1\right)\Delta \right)\right)\\ \sin \left({\omega }_{11}\left(t-\left(m-2\right)\Delta \right)\right)\end{array}\\ \begin{array}{c}⋮\\ \sin \left({\omega }_{11}t\right)\end{array}\end{array}\end{array}& \begin{array}{c}\begin{array}{c}\cdots \\ \cdots \end{array}\\ \begin{array}{c}⋮\\ \cdots \end{array}\end{array}\end{array} \\ \begin{array}{cc}\begin{array}{c}\begin{array}{c}\cos \left({\omega }_{2n}\left(t-\left(m-1\right)\Delta \right)\right)\\ \cos \left({\omega }_{2n}\left(t-\left(m-2\right)\Delta \right)\right)\end{array}\\ \begin{array}{c}⋮\\ \cos \left({\omega }_{2n}t\right)\end{array}\end{array}& \begin{array}{c}\begin{array}{c}\sin \left({\omega }_{2n}\left(t-\left(m-1\right)\Delta \right)\right)\\ \sin \left({\omega }_{2n}\left(t-\left(m-2\right)\Delta \right)\right)\end{array}\\ \begin{array}{c}⋮\\ \sin \left({\omega }_{2n}t\right)\end{array}\end{array}\end{array}]\left[\begin{array}{c}\begin{array}{c}{C}_{11}\\ D_{11}\end{array}\\ ⋮\\ \begin{array}{c}{C}_{2n}\\ D_{2n}\end{array}\end{array}\right] \end{gathered}$$

(33)

Using the least square optimization method, the coefficients vector ${\mathit{P}}_1$ can be estimated as:

$${\widehat{\mathit{P}}}_1={\left({\mathit{A}}_1^T{\mathit{A}}_1\right)}^{-1}{\mathit{A}}_1^T{\mathit{Q}}_1$$

(34)

It is noteworthy that, in the design of the input signal (19), each sinusoidal component satisfies the following orthogonal properties:

$${{\displaystyle \int }}_0^T\cos \left({\omega }_{ki}\right)\cos \left({\omega }_{kj}\right)dt=0, k=1,2, i\ne j$$

(35)

$${{\displaystyle \int }}_0^T\sin \left({\omega }_{ki}\right)\sin \left({\omega }_{kj}\right)dt=0, k=1,2, i\ne j$$

(36)

$${{\displaystyle \int }}_0^T\sin \left({\omega }_{ki}\right)\cos \left({\omega }_{kj}\right)dt=0, k=1,2, i\ne j$$

(37)

In addition,

$${{\displaystyle \int }}_0^T{\sin }^2\left({\omega }_{ki}\right)dt=\frac{T}{2}, k=1,2$$

(38)

$${{\displaystyle \int }}_0^T{\cos }^2\left({\omega }_{ki}\right)dt=\frac{T}{2}, k=1,2$$

(39)

In practice, in order to obtain the accurate estimation of the coefficients $C_{11i}$, $C_{12i}$, $D_{11i}$ and $D_{12i}$, for $i=1,2,\cdots ,n$, the output measurement needs to be sampled within a fundamental period $T$. Furthermore, the sampling period $\Delta$ is selected to be small enough, (35)–(39) can then be approximated as follows:

$$\underset{0}{\overset{T}{{\displaystyle \int }}}\cos \left({\omega }_{ki}\right)\cos \left({\omega }_{kj}\right)dt\approx {\displaystyle \sum }_{n=1}^m\cos \left({\omega }_{ki}-\left(n-1\right)\Delta \right)\cos \left({\omega }_{kj}-\left(n-1\right)\Delta \right)\Delta \approx 0, k=1,2, i\ne j$$

(40)

$$\underset{0}{\overset{T}{{\displaystyle \int }}}\sin \left({\omega }_{ki}\right)\sin \left({\omega }_{kj}\right)dt\approx {\displaystyle \sum }_{n=1}^m\sin \left({\omega }_{ki}-\left(n-1\right)\Delta \right)\sin \left({\omega }_{kj}-\left(n-1\right)\Delta \right)\Delta \approx 0, k=1,2, i\ne j$$

(41)

$$\underset{0}{\overset{T}{{\displaystyle \int }}}\sin \left({\omega }_{ki}\right)\cos \left({\omega }_{kj}\right)dt\approx {\displaystyle \sum }_{n=1}^m\sin \left({\omega }_{ki}-\left(n-1\right)\Delta \right)\cos \left({\omega }_{kj}-\left(n-1\right)\Delta \right)\Delta \approx 0, k=1,2, i\ne j$$

(42)

$$\underset{0}{\overset{T}{{\displaystyle \int }}}{\sin }^2\left({\omega }_{ki}\right)dt\approx {\displaystyle \sum }_{n=1}^m{\sin }^2\left({\omega }_{ki}-\left(n-1\right)\Delta \right)\Delta \approx \frac{T}{2}, k=1,2$$

(43)

$$\underset{0}{\overset{T}{{\displaystyle \int }}}{\cos }^2\left({\omega }_{ki}\right)dt\approx {\displaystyle \sum }_{n=1}^m{\cos }^2\left({\omega }_{ki}-\left(n-1\right)\Delta \right)\Delta \approx \frac{T}{2}, k=1,2$$

(44)

Based on (40)–(44), the term ${\left(A_1^T{A}_1\right)}^{-1}$ in (34) can then be approximate as:

$${\left({\mathit{A}}_1^T{\mathit{A}}_1\right)}^{-1}\approx {\left[\begin{array}{ccc}\frac{T}{2\Delta }& 0& 0\\ ⋮& \ddots & ⋮\\ 0& \dots & \frac{T}{2\Delta }\end{array}\right]}^{-1}=\frac{2\Delta }{T}\left[\begin{array}{ccc}1& 0& 0\\ ⋮& \ddots & ⋮\\ 0& \dots & 1\end{array}\right]$$

(45)

Furthermore,

$${\mathit{A}}_1^T{\mathit{Q}}_1\approx \frac{T}{2\Delta }{\mathit{P}}_1$$

(46)

Thus,

$${\widehat{\mathit{P}}}_1={\left({\mathit{A}}_1^T{\mathit{A}}_1\right)}^{-1}{\mathit{A}}_1^T{\mathit{Q}}_1\approx \frac{2\Delta }{T}{\mathit{A}}_1^T{\mathit{Q}}_1\approx {\mathit{P}}_1$$

(47)

The steady state coefficients $C_{11i}$, $D_{11i}$, $C_{12i}$ and $D_{12i}$. can then be numerically estimated according to (47) as:

$$\{\begin{array}{l}{\widehat{C}}_{11i}=\frac{2}{T}{\displaystyle {\displaystyle \sum }_{k=1}^m}{q}_{1s}\left(t-\left(k-1\right)\Delta \right)\cos \left({\omega }_{1i}\left(t-\left(k-1\right)\Delta \right)\right)\Delta \hfill \\ {\widehat{D}}_{11i}=\frac{2}{T}{\displaystyle {\displaystyle \sum }_{k=1}^m}{q}_{1s}\left(t-\left(k-1\right)\Delta \right)\sin \left({\omega }_{1i}\left(t-\left(k-1\right)\Delta \right)\right)\Delta \hfill \\ {\widehat{C}}_{12i}=\frac{2}{T}{\displaystyle {\displaystyle \sum }_{k=1}^m}{q}_{1s}\left(t-\left(k-1\right)\Delta \right)\cos \left({\omega }_{2i}\left(t-\left(k-1\right)\Delta \right)\right)\Delta \hfill \\ {\widehat{D}}_{12i}=\frac{2}{T}{\displaystyle {\displaystyle \sum }_{k=1}^m}{q}_{1s}\left(t-\left(k-1\right)\Delta \right)\sin \left({\omega }_{2i}\left(t-\left(k-1\right)\Delta \right)\right)\Delta \hfill \end{array}$$

(48)

It is shown from (45) that the computational cost of the matrix inverse calculation has been greatly decreased after the numerical approximation of the sinusoidal components based on their orthogonal properties. Moreover, the averaging process in (48) also contributes to reduce the effects of the measurement noises and, as a result, guarantee the robustness of the proposed algorithm against measurement noises.

Since the accurate estimation of the coefficient vector ${\mathit{P}}_1$of the steady state of the joint position measurement can be obtained from (48), the identification algorithm for estimating the unknown dynamic parameters can be further developed.

3.2. Estimation of Unknown Dynamic Parameters

As discussed in (27)–(30), the joint angular position output is dominated by its steady state components after sufficient long time; thus:

$$q_{1i}\left(t\right)\approx q_{1si}\left(t\right)=C_{11i}\cos \left({\omega }_{1i}t\right)+D_{11i}\sin \left({\omega }_{1i}t\right)+C_{12i}cos\left({\omega }_{2i}t\right)+D_{12i}\sin \left({\omega }_{2i}t\right)$$

(49)

Based on (21) and (23), the following relations can be derived:

$$\begin{array}{c}{C}_{11i}{\omega }_{1i}^2=C_{11i}{\omega }_0^2-2D_{11i}{\omega }_{1i}{\sigma }_0-d_{1i}{P}_1\end{array}$$

(50)

$$D_{1i}{\omega }_{1i}^2=D_{1i}{\omega }_0^2+2C_{1i}{\omega }_{1i}{\sigma }_0+d_{1i}{\omega }_{1i}{b}_{11}$$

(51)

Since $C_{11i}$ and $D_{11i}$ can be accurately estimated via (45), (50) and (51) can then be re-written as:

$${\widehat{C}}_{11i}{\omega }_{1i}^2={\widehat{C}}_{11i}{\widehat{\omega }}_0^2-2{\widehat{D}}_{11i}{\omega }_{1i}{\widehat{\sigma }}_0-d_{1i}{\widehat{P}}_1$$

(52)

$${\widehat{D}}_{11i}{\omega }_{1i}^2={\widehat{D}}_{11i}{\widehat{\omega }}_0^2+2{\widehat{C}}_{11i}{\omega }_{1i}{\widehat{\sigma }}_0+d_{1i}{\omega }_{1i}{\widehat{b}}_{11}$$

(53)

To estimate unknown parameters, sufficient number of coefficients $C_{11}$ and $D_{11}$ need to be obtained via (47). Then, the following data equation in matrix form can be derived:

$${\mathit{Q}}_2={\mathit{A}}_2{\mathit{P}}_2$$

(54)

where

$${\mathit{Q}}_2={[\begin{array}{ccc}\begin{array}{cc}{\widehat{C}}_{111}{\omega }_{11}^2& {\widehat{C}}_{112}{\omega }_{12}^2\end{array}& \begin{array}{cc}\dots & {\widehat{D}}_{111}{\omega }_{11}^2\end{array}& \begin{array}{cc}\dots & {\widehat{C}}_{11n}{\omega }_{1n}^2\end{array}\end{array}]}^T$$

(55)

$${\mathit{A}}_2=\left[\begin{array}{cc}\begin{array}{cc}\begin{array}{c}\begin{array}{c}{\widehat{C}}_{111}\\ {\widehat{C}}_{112}\end{array}\\ \begin{array}{c}⋮\\ \begin{array}{c}{\widehat{D}}_{111}\\ \begin{array}{c}⋮\\ {\widehat{D}}_{11n}\end{array}\end{array}\end{array}\end{array}& \begin{array}{c}\begin{array}{c}-2{\widehat{D}}_{112}{\omega }_{12}\\ -2{\widehat{D}}_{112}{\omega }_{12}\end{array}\\ \begin{array}{c}⋮\\ \begin{array}{c}2{\widehat{C}}_{111}{\omega }_{11}\\ \begin{array}{c}⋮\\ 2{\widehat{C}}_{11n}{\omega }_{1n}\end{array}\end{array}\end{array}\end{array}\end{array}& \begin{array}{cc}\begin{array}{c}\begin{array}{c}-d_{11}\\ -d_{12}\end{array}\\ \begin{array}{c}⋮\\ \begin{array}{c}0\\ \begin{array}{c}⋮\\ 0\end{array}\end{array}\end{array}\end{array}& \begin{array}{c}\begin{array}{c}0\\ 0\end{array}\\ \begin{array}{c}⋮\\ \begin{array}{c}{d}_{11}{\omega }_{11}\\ \begin{array}{c}⋮\\ d_{1n}{\omega }_{1n}\end{array}\end{array}\end{array}\end{array}\end{array}\end{array}\right]$$

(56)

$${\mathit{P}}_2={\left[\begin{array}{cc}\begin{array}{cc}{\widehat{\omega }}_0^2& {\widehat{\sigma }}_0\end{array}& \begin{array}{cc}{\widehat{P}}_1& {\widehat{b}}_{11}\end{array}\end{array}\right]}^T$$

(57)

Using the least square method, the parameter vector ${\mathit{P}}_2$ can be estimated as:

$${\mathit{P}}_2={\left({\mathit{A}}_2^T{\mathit{A}}_2\right)}^{-1}{\mathit{A}}_2^T{\mathit{Q}}_2$$

(58)

where all elements in ${\mathit{P}}_2$ are equations with respect to the manipulator’s mechanical quantities including the mass of the link, the length of the link, etc.

Similarly, the $i_{\mathrm{th}}$ component of (22) can also be expanded using the partial fraction expansion method as:

$$Q_{2i}\left(s\right)=\frac{A_{21i}\left(s-{\sigma }_0\right)+B_{21i}{\omega }_d}{{\left(s-{\sigma }_0\right)}^2+{\omega }_d^2}+\frac{C_{21i}s+D_{21i}{\omega }_{1i}}{s^2+{\omega }_{1i}^2}+\frac{A_{22i}\left(s-{\sigma }_0\right)+B_{22i}{\omega }_d}{{\left(s-{\sigma }_0\right)}^2+{\omega }_d^2}+\frac{C_{22i}s+D_{22i}{\omega }_{2i}}{s^2+{\omega }_{2i}^2}$$

(59)

Applying the same process, $m$ samples of the second joint’s angular position measurement are sampled within one fundamental period; sample period $\Delta$ is selected to be sufficiently small. The coefficients of the steady state component of (59) can be numerically estimated as:

$$\{\begin{array}{c}\begin{array}{c}{\widehat{C}}_{21i}=\frac{2}{T}{\displaystyle {\displaystyle \sum }_{k=1}^m}{q}_{2s}\left(t-\left(k-1\right)\Delta \right)\cos \left({\omega }_{1i}\left(t-\left(k-1\right)\Delta \right)\right)\Delta \\ {\widehat{D}}_{21i}=\frac{2}{T}{\displaystyle {\displaystyle \sum }_{k=1}^m}{q}_{2s}\left(t-\left(k-1\right)\Delta \right)\sin \left({\omega }_{1i}\left(t-\left(k-1\right)\Delta \right)\right)\Delta \end{array}\\ \begin{array}{c}{\widehat{C}}_{22i}=\frac{2}{T}{\displaystyle {\displaystyle \sum }_{k=1}^m}{q}_{2s}\left(t-\left(k-1\right)\Delta \right)\cos \left({\omega }_{2i}\left(t-\left(k-1\right)\Delta \right)\right)\Delta \\ {\widehat{D}}_{22i}=\frac{2}{T}{\displaystyle {\displaystyle \sum }_{k=1}^m}{q}_{2s}\left(t-\left(k-1\right)\Delta \right)\sin \left({\omega }_{2i}\left(t-\left(k-1\right)\Delta \right)\right)\Delta \end{array}\end{array}$$

(60)

Based on (22) and (59), the following relationships can be obtained:

$$C_{21i}{\omega }_{2i}^2=C_{21i}{\omega }_0^2-2D_{21i}{\omega }_{2i}{\sigma }_0-d_{2i}{P}_3$$

(61)

$$D_{2i}{\omega }_{2i}^2=D_{2i}{\omega }_0^2+2C_{2i}{\omega }_{2i}{\sigma }_0+d_{2i}{\omega }_{2i}{b}_{21}$$

(62)

Since $C_{21i}$ and $D_{21i}$ can be accurately estimated via (50), (51) and (60) can then be re-written as:

$${\widehat{C}}_{21i}{\omega }_{2i}^2={\widehat{C}}_{21i}{\widehat{\omega }}_0^2-2{\widehat{D}}_{21i}{\omega }_{2i}{\widehat{\sigma }}_0-d_{2i}{\widehat{\mathit{P}}}_3$$

(63)

$${\widehat{D}}_{21i}{\omega }_{2i}^2={\widehat{D}}_{21i}{\widehat{\omega }}_0^2+2{\widehat{C}}_{21i}{\omega }_{2i}{\widehat{\sigma }}_0+d_{2i}{\omega }_{1i}{\widehat{b}}_{21}$$

(64)

To estimate unknown parameters, a sufficient number of coefficients $C_{21}$ and $D_{21}$ needs to be obtained via (60). Then, the following data equation in matrix form can be derived:

$${\mathit{Q}}_3={\mathit{A}}_3{\mathit{P}}_3$$

(65)

where

$${\mathit{Q}}_3={\left[\begin{array}{ccc}\begin{array}{cc}{\widehat{C}}_{211}{\omega }_{21}^2& {\widehat{C}}_{212}{\omega }_{22}^2\end{array}& \begin{array}{cc}\dots & {\widehat{D}}_{211}{\omega }_{21}^2\end{array}& \begin{array}{cc}\dots & {\widehat{C}}_{21n}{\omega }_{2n}^2\end{array}\end{array}\right]}^T$$

(66)

$${\mathit{A}}_3=\left[\begin{array}{cc}\begin{array}{cc}\begin{array}{c}\begin{array}{c}{\widehat{C}}_{211}\\ {\widehat{C}}_{212}\end{array}\\ \begin{array}{c}⋮\\ \begin{array}{c}{\widehat{D}}_{211}\\ \begin{array}{c}⋮\\ {\widehat{D}}_{21n}\end{array}\end{array}\end{array}\end{array}& \begin{array}{c}\begin{array}{c}-2{\widehat{D}}_{212}{\omega }_{22}\\ -2{\widehat{D}}_{212}{\omega }_{22}\end{array}\\ \begin{array}{c}⋮\\ \begin{array}{c}2{\widehat{C}}_{211}{\omega }_{21}\\ \begin{array}{c}⋮\\ 2{\widehat{C}}_{21n}{\omega }_{2n}\end{array}\end{array}\end{array}\end{array}\end{array}& \begin{array}{cc}\begin{array}{c}\begin{array}{c}-d_{21}\\ -d_{22}\end{array}\\ \begin{array}{c}⋮\\ \begin{array}{c}0\\ \begin{array}{c}⋮\\ 0\end{array}\end{array}\end{array}\end{array}& \begin{array}{c}\begin{array}{c}0\\ 0\end{array}\\ \begin{array}{c}⋮\\ \begin{array}{c}{d}_{21}{\omega }_{21}\\ \begin{array}{c}⋮\\ d_{2n}{\omega }_{2n}\end{array}\end{array}\end{array}\end{array}\end{array}\end{array}\right]$$

(67)

$${\mathit{P}}_3={\left[\begin{array}{cc}\begin{array}{cc}{\widehat{\omega }}_0^2& {\widehat{\sigma }}_0\end{array}& \begin{array}{cc}{\widehat{p}}_3& {\widehat{b}}_{21}\end{array}\end{array}\right]}^T$$

(68)

Using the least square method, the parameter vector ${\mathit{P}}_3$ can be estimated as:

$${\mathit{P}}_3={\left({\mathit{A}}_3^T{\mathit{A}}_3\right)}^{-1}{\mathit{A}}_3^T{\mathit{Q}}_3$$

(69)

Once ${\mathit{P}}_2$ and ${\mathit{P}}_3$ can be estimated via (58) and (69), respectively, all the unknown system dynamic parameters as well as the dynamic quantities can then be derived.

In the proposed method, the parameter estimator for each joint is separately designed. Since the dynamic information is embedded within the steady state of the joint angle measurement, only the joint angle measurement in one fundamental period is required to estimate all dynamic parameters. Meanwhile, the orthogonality of the trigonometric base function contributes to reducing the computation load of the high dimension matrix. As a result, the proposed method can be extended to high DOF manipulators.

4. Simulation Results

In this section, simulation with a 2-DOF manipulator is performed to verify the feasibility of the proposed parameter estimation algorithm. The nominal values of robotic manipulator parameters are given in the

Table 1. Parameter of the robotic manipulator.

Parameters	Values
$\mathrm{Mass} \mathrm{of} \mathrm{link} 1 {\mathrm{m}}_1$$\left(\mathrm{kg}\right)$	5
$\mathrm{Mass} \mathrm{of} \mathrm{link} 2 {\mathrm{m}}_2$$\left(\mathrm{kg}\right)$	1.5
$\mathrm{Inertial} \mathrm{of} \mathrm{link} 1 j_1$$(\mathrm{kg} ·{\mathrm{m}}^2)$	5
$\mathrm{Inertial} \mathrm{of} \mathrm{link} 2 j_2$$(\mathrm{kg} ·{\mathrm{m}}^2)$	5
$\mathrm{Length} \mathrm{of} \mathrm{link} 1 l_1$$\left(\mathrm{m}\right)$	1
$\mathrm{Length} \mathrm{of} \mathrm{link} 2 l_2$$\left(\mathrm{m}\right)$	0.8

[39]:

In this study, the masses of two links are supposed to be unknown. The input torque is designed to be a linear combination of 10 sinusoidal components as shown in

Table 2. Sinusoidal components of the input torque signal ${\tau }_1$.

${\mathit{i}}_{\mathbf{th}}$ Component	Frequency	Amplitude
1	$3{\omega }_0$	0.2
2	$5{\omega }_0$	1.5
3	$7{\omega }_0$	0.3
4	$9{\omega }_0$	0.5
5	$11{\omega }_0$	0.7
6	$13{\omega }_0$	0.8
7	$17{\omega }_0$	1.0
8	$21{\omega }_0$	1.2
9	$27{\omega }_0$	1.4
10	$42{\omega }_0$	2.0

[39]. The fundamental frequency is selected to be $f_0=0.01 \mathrm{Hz}$ and ${\omega }_0=2\pi f_0$.

The system input torque with 10 sinusoidal components is shown in Figure 1 [39]. The above torque has been applied to both joints at the same time.

The comparison between the theoretical value and calculated value of coefficients $C_{11}$ $C_{21}$, $D_{11}$ and $D_{21}$ is illustrated in Figure 2 and Figure 3 [39], respectively:

It can be seen from Figure 2 and Figure 3 that, with the fundamental frequency $f_0=0.01 \mathrm{Hz}$ as the sampling interval and 10 consecutive components, the coefficients of the steady state components of the joint angular position can be accurately estimated based on (48) and (60).

In practical implementations, the measurement of the angular position is vulnerable to being affected by the sensor noise. In order to illustrate the robustness of the proposed identification algorithm against the measurement noises, a white noise $d=0.1\ast rand()$ is manually added to system output to simulate the actual measurement. The estimation performance of two unknown link masses can be seen in Figure 4 and Figure 5 [39], respectively.

It can be seen from Figure 4 and Figure 5 that the unknown parameters derived based on (58) and (65) can well approximate their actual values. The effect of measurement noises has been reasonably compensated by the averaging process during the calculation of the steady state coefficients in (48) and (60).

To illustrate the accuracy of the estimation, the RMSEs (root-mean square error) for two joints are calculated:

$$rms{e}_{m_1 }=1.0568$$

(70)

$$rms{e}_{m_2 }=1.5420$$

(71)

5. Conclusions

In this paper, a novel mythology to estimate the unknown parameters of the robotic manipulator system in frequency domain has been proposed. It is seen that the joint angular position measurement will be dominated by its steady state component after sufficient long time. By designing the input signal to be a linear combination of sinusoidal components, the steady state component of the measurement is equivalent to a finite Fourie series. It is shown that all system dynamic information is embedded in the coefficients of the steady state component within one fundamental period. Thus, the coefficients of the steady state sinusoidal component of the joint angular position output are first extracted. By utilizing the relationship among these coefficients and the system’s dynamic parameters in the frequency domain, all unknown dynamic parameters and mechanical quantities can be estimated accurately. It is expected that this new kind of new estimation scheme can be widely applied in various industrial applications where the systems suffer from uncertainties and disturbances.