3.1. Recursive Newton–Euler formulation derivation

Owing to the effects of inertia, centrifugal force, Coriolis force, gravity, and external force vectors in joint space, the associated inverse dynamics for an n-DOF serial robot can be described by the second-order differential equation:

(1) \begin{align}{{\boldsymbol{\tau }}_{d}}=\textbf{M}(\textbf{q})\ddot{\textbf{q}}+\textbf{V}\left ( \textbf{q},\dot{\textbf{q}} \right )\dot{\textbf{q}}+\textbf{G}(\textbf{q})+{{\textbf{J}}^{\textit{T}}}(\textbf{q}){{\textbf{F}}_{ext}} \end{align}

where $\textbf{M}(\textbf{q})$ is the inertia matrix, Coriolis and centripetal effects are captured by matrix $\textbf{V}( \textbf{q},\dot{\textbf{q}})$ , $\textbf{G}(\textbf{q})$ is the gravity vector, ${\textbf{J}}(\textbf{q})$ signifies the corresponding Jacobian matrix, and ${\textbf{F}}_{ext}$ is the external wrench exerted by the end-effector.

Due to Newton–Euler-based formulations for dynamics, using Lie theory require only four equations, resulting in a more concise dynamic modeling process, which involves a two-step recursion, that is, forward recursion and backward recursion [Reference Park, Bobrow and Ploen27, Reference Yang, Yang and Wang44]. As shown in Figure 2, taking the $i$ -th joint as an example, the corresponding recursion process can be summarized as follows.

Figure 2. Variables relationship of adjacent links.

3.1.1. Forward recursion

The six-dimensional generalized velocity $\mathbf{V}_{i}$ and acceleration $\mathbf{A}_{i}$ of $i$ -th link expressed in the frame $\{i\}$ can be computed by the ones $\mathbf{V}_{i-1}$ and $\mathbf{A}_{i-1}$ of $(i-1)$ -th link, and the velocity $\dot{\textbf{q}}_i$ and acceleration $\ddot{\textbf{q}}_i$ of $i$ -th joint:

(2) \begin{align} &{\mathbf{V}_{i}}={{Ad}_{\mathbf{g}_{i-1,i}^{-1}}}\left ({\mathbf{V}_{i-1}} \right )+{{\boldsymbol{\xi } }_{i}}{{{\dot{\textbf{q}}}}_{i}} \nonumber \\ &{\mathbf{A}_{i}}={{Ad}_{\mathbf{g}_{i-1,i}^{-1}}}\left ({\mathbf{A}_{i-1}} \right )-{{ad}_{{{Ad}_{\mathbf{g}_{i-1,i}^{-1}}}\left ({\mathbf{V}_{i-1}} \right )}}\left ({{{\boldsymbol{\xi } }_{i}}{{{\dot{\textbf{q}}}}_{i}}} \right )+{{\boldsymbol{\xi } }_{i}}{{{{\ddot{\textbf{q}}}}}_{i}} \end{align}

where $\boldsymbol{\xi }_{i}$ is the twist coordinate of the $i$ -th joint in the link frame.

3.1.2. Backward recursion

On the basis of the kinematic variables in Eq. (2), the generalized force $\mathbf{F}_i$ can be calculated by the generalized force $\mathbf{F}_{i+1}$ of the $(i+1)$ -th link on the $i$ -th link, as following:

(3) \begin{align} &{\mathbf{F}_i}=Ad_{\mathbf{g}_{i,i+1}^{-1}}^{\textit{T}}\left ({\mathbf{F}_{i+1}} \right )+{\mathbf{E}_{i}}{\mathbf{A}_{i}}-ad_{{\mathbf{V}_{i}}}^{\textit{T}}\left ({\mathbf{E}_{i}}{\mathbf{V}_{i}} \right ) \nonumber \\ &{{\boldsymbol{\tau }}_{di}}={\boldsymbol{\xi }}_{i}^{\textit{T}}{\mathbf{F}_i} \end{align}

where ${\mathbf{E}_{i}}=\left [\begin{matrix}{\mathbf{I}_{i}}-{{m}_{i}}{\hat{\mathbf{r}}}_{i}^{2} &{{m}_{i}}{{\hat{\mathbf{r}}}_{i}} \\ -{{m}_{i}}{{\hat{\mathbf{r}}}_{i}} &{{m}_{i}}{{\textbf{I}}_{3}} \\ \end{matrix} \right ]$ is the generalized inertial matrix with respect to $i$ -th link.

On the basis of the above analysis, the Newton–Euler dynamic model of the robot based on Lie theory can be expressed as (refer to [Reference Fu, Pan, Spyrakos-Papastavridis, Lin, Zhou, Chen and Dai6] for the detailed derivation):

(4) \begin{align}{{\boldsymbol{\tau }}_{d}}&={{\boldsymbol{\xi } }^{\textit{T}}}{{\textbf{B}}^{\textit{T}}}{\textbf{E}}{\textbf{B}}{\boldsymbol{\xi }}{\ddot{\textbf{q}}}-{{\boldsymbol{\xi } }^{\textit{T}}}{{\textbf{B}}^{\textit{T}}}({\textbf{E}}{\textbf{B}}{{\boldsymbol{{\boldsymbol{\Lambda } }} }_{{\boldsymbol{\xi }}{\dot{\textbf{q}}}}}\boldsymbol{\Gamma } +{{{\boldsymbol{\Lambda }}}_{\textbf{V}}}{\textbf{E}}){\textbf{B}}{\boldsymbol{\xi }}{\dot{\textbf{q}}} +{{\boldsymbol{\xi } }^{\textit{T}}}{{\textbf{B}}^{\textit{T}}}{\textbf{E}}{\textbf{B}}{{\boldsymbol{\Gamma } }_{{{\textbf{V}}_{0}}}}{{\textbf{A}}_{0}}+{{\boldsymbol{\xi } }^{\textit{T}}}{{\textbf{B}}^{\textit{T}}}\boldsymbol{\Gamma } _{{{\textbf{F}}_{0}}}^{\textit{T}}{{\textbf{F}}_{ext}} \end{align}

where

(5) \begin{align} & \textbf{M}(\textbf{q})={{\boldsymbol{\xi } }^{\textit{T}}}{{\textbf{B}}^{\textit{T}}}{\textbf{E}}{\textbf{B}}{\boldsymbol{\xi }} \nonumber \\ & \textbf{V}\left ( \textbf{q},\dot{\textbf{q}} \right )=-{{\boldsymbol{\xi } }^{\textit{T}}}{{\textbf{B}}^{\textit{T}}}({\textbf{E}}{\textbf{B}}{{{\boldsymbol{\Lambda } } }_{{\boldsymbol{\xi }}{\dot{\textbf{q}}}}}\boldsymbol{\Gamma } +{{{\boldsymbol{\Lambda }}}_{\textbf{V}}}{\textbf{E}}){\textbf{B}}{\boldsymbol{\xi }} \nonumber \\ & \textbf{G}(\textbf{q})={{\boldsymbol{\xi } }^{\textit{T}}}{{\textbf{B}}^{\textit{T}}}{\textbf{E}}{\textbf{B}}{{\boldsymbol{\Gamma } }_{{{\textbf{V}}_{0}}}}{{\textbf{A}}_{0}} \nonumber \\ &{\textbf{J}}(\textbf{q})={{\boldsymbol{\Gamma } }_{{{\textbf{F}}_{0}}}}{\textbf{B}}{\boldsymbol{\xi }} \end{align}

3.2. Nonlinear friction modeling

The joint friction, which is relatively complex and presents typical nonlinear characteristics, should be considered for accurate dynamic modeling in practical applications. This is because the friction torque decreases with increasing velocity after overcoming static friction. After the improvement of the nonlinear Tustin friction model [Reference Tustin35], the nonlinear friction model is constructed as:

(6) \begin{align}{\boldsymbol{\tau }_{f}} ={\boldsymbol{f}}_{v}{\dot{\textbf{q}}} +{\boldsymbol{f}}_{c}sgn\left ({{{\dot{\textbf{q}}}}}\right )\left (1-e^{-\left (sgn\left ({{{\dot{\textbf{q}}}}}\right ){\dot{\textbf{q}}}/{\dot{\textbf{q}}_e}\right )^{\delta _s}}\right ) +{\boldsymbol{f}}_{s}sgn\left ({{{\dot{\textbf{q}}}}}\right )e^{-\left (sgn\left ({{{\dot{\textbf{q}}}}}\right ){\dot{\textbf{q}}}/{\dot{\textbf{q}}_e}\right )^{\delta _s}} \end{align}

where ${\boldsymbol{f}}_v$ is viscous friction, ${\boldsymbol{f}}_c$ represents the Coulomb friction, ${\boldsymbol{f}}_s$ is static or stiction friction, $\dot{\textbf{q}}_e$ is Stribeck velocity, $\delta _s$ is the exponent parameter of the Stribeck nonlinearity, and $sgn\left ({{{\dot{\textbf{q}}}}} \right )$ is sign function of $\dot{\textbf{q}}$ .

In the Gaussian parametrization as used in [Reference Tustin35], the exponent parameter is $\delta _s=2$ . Moreover, the parameter $\frac{1}{{\dot{\textbf{q}}_e}}$ that corresponds to ${\boldsymbol{f}}_{s}$ and ${\boldsymbol{f}}_{c}$ terms is replaced by two empirical parameters, namely $\boldsymbol{f}_{s1}$ and $\boldsymbol{f}_{c1}$ [Reference Dong, Xu, Zhou, Zhu and Yu5]. Thus, the nonlinear friction model at the commutation and low velocities can be modified through introducing a third-order polynomial as follows:

(7) \begin{align}{\boldsymbol{\tau }_{f}} ={\boldsymbol{f}}_{v_1}{\dot{\textbf{q}}} +{\boldsymbol{f}}_{v_2}{sgn\left ({{{\dot{\textbf{q}}}}}\right )}{\dot{\textbf{q}}^2} +{\boldsymbol{f}}_{v_3}{\dot{\textbf{q}}^3} +{\boldsymbol{f}}_{s}sgn\left ({{{\dot{\textbf{q}}}}}\right )e^{-\left ({\dot{\textbf{q}}}/\boldsymbol{f}_{s1}\right )^2} +{\boldsymbol{f}}_{c}sgn\left ({{{\dot{\textbf{q}}}}}\right )\left (1-e^{-\left ({\dot{\textbf{q}}}/\boldsymbol{f}_{c1}\right )^2}\right ) \end{align}

A significant dependence on the temperature and external load has been observed for friction torque. In previous studies (see ref. [Reference Gao, Yuan, Han, Wang and Wang7, Reference He, Wang, Bao, Yuan, Du, Ma and Wan14]), we find that friction increases exponentially with temperature, while load torque has a positive correlation with friction, resulting in a Stribeck effect at low speeds. The corresponding joint temperature reaches a relatively dynamic equilibrium state when each joint was run continuously for a period of time at a set trajectory. Thus, the nonlinear friction model (7) can be further extended as follows:

(8) \begin{align}{\boldsymbol{\tau }_{f}} & = {\boldsymbol{f}}_{v_1}{\dot{\textbf{q}}} +{\boldsymbol{f}}_{v_2}{sgn\left ({{{\dot{\textbf{q}}}}}\right )}{\dot{\textbf{q}}^2} +{\boldsymbol{f}}_{v_3}{\dot{\textbf{q}}^3} +{\boldsymbol{f}}_{c}sgn\left ({{{\dot{\textbf{q}}}}}\right )\left (1-e^{-\left ({\dot{\textbf{q}}}/\boldsymbol{f}_{c1}\right )^2}\right ) +{\boldsymbol{f}}_{s}sgn\left ({{{\dot{\textbf{q}}}}}\right )e^{-\left ({\dot{\textbf{q}}}/\boldsymbol{f}_{s1}\right )^2} \nonumber \\ &\quad +\left ({{\boldsymbol{f}}_{{{\boldsymbol{\tau }}_{ext}}}}\vert{{\boldsymbol{\tau }}_{ext}}\vert +{{\boldsymbol{f}}_{{{\boldsymbol{\tau }}_{ext1}}}}{e}^{-\boldsymbol{f}_{{{\boldsymbol{\tau }}_{ext2}}}\vert{{\boldsymbol{\tau }}_{ext}}\vert /{\dot{\textbf{q}}}}\right ) sgn\left ({{{\dot{\textbf{q}}}}} \right ) \end{align}

where ${\boldsymbol{\tau }}_{ext}$ represents the load torque vector, including gravity torque of the robot itself, inertial torque, coupling torque, and external load torque; ${\boldsymbol{f}}_{c{{\boldsymbol{\tau }}_{ext}}}$ , ${\boldsymbol{f}}_{s{{\boldsymbol{\tau }}_{ext}}}$ , and ${{\dot{\textbf{q}}}}_{{{\tau }_{0}}}$ are three parameters of the load-dependent friction term.

3.3. Lie-theory-based dynamics generation

Combining the Newton–Euler dynamics (1) based on the Lie theory with an improved nonlinear friction model $\left (8\right )$ , the dynamic model for a serial robot consisting of $n$ rigid links interconnected via joints can be generated as:

(9) \begin{align}{{\boldsymbol{\tau }}_{Joint}} &={{\boldsymbol{\tau }}_{d}}+{{\boldsymbol{\tau }}_{f}}=\textbf{M}(\textbf{q})\ddot{\textbf{q}}+\textbf{V}\left ( \textbf{q},\dot{\textbf{q}} \right )\dot{\textbf{q}}+\textbf{G}(\textbf{q})+{{\textbf{J}}^{\textit{T}}}(\textbf{q}){{\textbf{F}}_{ext}}+{{\boldsymbol{\tau }}_{\boldsymbol{f}}} \nonumber \\ &={{\boldsymbol{\xi } }^{\textit{T}}}{{\textbf{B}}^{\textit{T}}}{\textbf{E}}{\textbf{B}}{\boldsymbol{\xi }}{\ddot{\textbf{q}}}-{{\boldsymbol{\xi } }^{\textit{T}}}{{\textbf{B}}^{\textit{T}}}({\textbf{E}}{\textbf{B}}{{{\boldsymbol{\Lambda } } }_{{\boldsymbol{\xi }}{\dot{\textbf{q}}}}}\boldsymbol{\Gamma } +{{{\boldsymbol{\Lambda }}}_{\textbf{V}}}{\textbf{E}}){\textbf{B}}{\boldsymbol{\xi }}{\dot{\textbf{q}}}+{{\boldsymbol{\xi } }^{\textit{T}}}{{\textbf{B}}^{\textit{T}}}{\textbf{E}}{\textbf{B}}{{\boldsymbol{\Gamma } }_{{{\textbf{V}}_{0}}}}{{\textbf{A}}_{0}}+{{\boldsymbol{\xi } }^{\textit{T}}}{{\textbf{B}}^{\textit{T}}}\boldsymbol{\Gamma } _{{{\textbf{F}}_{0}}}^{\textit{T}}{{\textbf{F}}_{ext}}\nonumber \\ & \quad + {\boldsymbol{f}}_{v_1}{\dot{\textbf{q}}} +{\boldsymbol{f}}_{v_2}{sgn\left ({{{\dot{\textbf{q}}}}}\right )}{\dot{\textbf{q}}^2} +{\boldsymbol{f}}_{v_3}{\dot{\textbf{q}}^3} +{\boldsymbol{f}}_{c}sgn\left ({{{\dot{\textbf{q}}}}}\right )\left (1-e^{-\left ({\dot{\textbf{q}}}/\boldsymbol{f}_{c1}\right )^2}\right ) +{\boldsymbol{f}}_{s}sgn\left ({{{\dot{\textbf{q}}}}}\right )e^{-\left ({\dot{\textbf{q}}}/\boldsymbol{f}_{s1}\right )^2}\nonumber \\ & \quad + \left ({{\boldsymbol{f}}_{{{\boldsymbol{\tau }}_{ext}}}}\vert{{\boldsymbol{\tau }}_{ext}}\vert +{{\boldsymbol{f}}_{{{\boldsymbol{\tau }}_{ext1}}}}{{e}^{-\boldsymbol{f}_{{{\boldsymbol{\tau }}_{ext2}}}\vert{{\boldsymbol{\tau }}_{ext}}\vert /{\dot{\textbf{q}}}}}\right ) sgn\left ({{{\dot{\textbf{q}}}}} \right ) \end{align}

where ${\boldsymbol{\tau }}_{Joint}=\left [{\boldsymbol{\tau }}_{Joint1} \cdots{\boldsymbol{\tau }}_{Jointi} \cdots{\boldsymbol{\tau }}_{Jointn}\right ]^{\textit{T}} \in{\Re }^n$ , ${\boldsymbol{\tau }}_{d}=\left [{\boldsymbol{\tau }}_{d1} \cdots{\boldsymbol{\tau }}_{di} \cdots{\boldsymbol{\tau }}_{dn}\right ]^{\textit{T}} \in{\Re }^n$ , and ${\boldsymbol{\tau }}_{f}=\left [{\boldsymbol{\tau }}_{f1} \cdots{\boldsymbol{\tau }}_{fi} \cdots{\boldsymbol{\tau }}_{fn}\right ]^{\textit{T}} \in{\Re }^n$ .

4.1. Dynamic parameter identification

Traditionally, a lot of the related studies employ a linear form to friction modeling in the context of dynamic identification for a serial robot. In general, the three coefficients are adopted to capture the viscous friciton and Coulomb friciton [Reference Han, Wu, Liu and Xiong13, Reference Huang, Ke, Zhang and Ota15], incorporating the Lie-theory-based Newton–Euler dynamics, the dynamic model of a serial robot is written as follows:

(10) \begin{align}{{\boldsymbol{\tau }}_{Joint}}={{\boldsymbol{\xi } }^{\textit{T}}}{{\textbf{B}}^{\textit{T}}}({\textbf{E}}{\textbf{B}}{\boldsymbol{\xi }}{\ddot{\textbf{q}}}-\left ({\textbf{E}}{\textbf{B}}{{{\boldsymbol{\Lambda } } }_{{\boldsymbol{\xi }}{\dot{\textbf{q}}}}}\boldsymbol{\Gamma }+{{{\boldsymbol{\Lambda }}}_{\textbf{V}}}{\textbf{E}}\right ){\textbf{B}}{\boldsymbol{\xi }}{\dot{\textbf{q}}}+{\textbf{E}}{\textbf{B}}{{\boldsymbol{\Gamma } }_{{{\textbf{V}}_{0}}}}{{\textbf{A}}_{0}} )\nonumber \\ +{{\boldsymbol{\xi } }^{\textit{T}}}{{\textbf{B}}^{\textit{T}}}\boldsymbol{\Gamma } _{{{\textbf{F}}_{0}}}^{\textit{T}}{{\textbf{F}}_{ext}}+{\boldsymbol{f}}_{v}{\dot{\textbf{q}}} +{\boldsymbol{f}}_{c}sgn\left ({{{\dot{\textbf{q}}}}}\right )+{\boldsymbol{f}}_{b} \end{align}

where $\boldsymbol{f}_c$ represents Coulomb friction coefficient vector, $\boldsymbol{f}_{v}$ denotes viscous friction coefficient vector, and $\boldsymbol{f}_b$ is the corresponding friction offset vector.

Due to the proposed friction model being strongly nonlinear, the dynamics (9) cannot be expressed in the parameter-linearized form. Simultaneously, the friction model in Eq. (10) is inaccurate for capturing the main characteristics of viscous friction in the high-velocity range. Hence, it is worth mentioning that one assumption is made to get the accurate inertia parameters, which is operating the robot with a constant temperature and zero external payload, so as to eliminate the impacts from external load. On the basis of the above assumption, the proposed dynamic model (9) can be simplified as follows:

(11) \begin{align}{{\boldsymbol{\tau }}_{Joint}} =\, & {{\boldsymbol{\xi } }^{\textit{T}}}{{\textbf{B}}^{\textit{T}}}({\textbf{E}}{\textbf{B}}{\boldsymbol{\xi }}{\ddot{\textbf{q}}}-\left ({\textbf{E}}{\textbf{B}}{{{\boldsymbol{\Lambda } } }_{{\boldsymbol{\xi }}{\dot{\textbf{q}}}}}\boldsymbol{\Gamma }+{{{\boldsymbol{\Lambda }}}_{\textbf{V}}}{\textbf{E}}\right ){\textbf{B}}{\boldsymbol{\xi }}{\dot{\textbf{q}}}+{\textbf{E}}{\textbf{B}}{{\boldsymbol{\Gamma } }_{{{\textbf{V}}_{0}}}}{{\textbf{A}}_{0}} ) \nonumber \\ & + {\boldsymbol{f}}_{v_1}{\dot{\textbf{q}}} +{\boldsymbol{f}}_{v_2}{sgn\left ({{{\dot{\textbf{q}}}}}\right )}{\dot{\textbf{q}}^2} +{\boldsymbol{f}}_{v_3}{\dot{\textbf{q}}^3}+{\boldsymbol{f}}_{c}sgn\left ({{{\dot{\textbf{q}}}}}\right )+{\boldsymbol{f}}_{b} \end{align}

where $\boldsymbol{f}_{v_1}$ , $\boldsymbol{f}_{v_2}$ , and $\boldsymbol{f}_{v_3}$ denote viscous friction coefficient vectors.

Using the Kronecker product and the Sylvester equation (for a more detailed derivation, please refer to ref. [Reference Fu, Pan, Spyrakos-Papastavridis, Lin, Zhou, Chen and Dai6]), Eq. (11) can be written as a linear equation form:

(12) \begin{align} \mathbf{Y}&=\textbf{J}\left (\textbf{q},{\dot{\textbf{q}}},{\ddot{\textbf{q}}} \right )\cdot{{\left [ \begin{matrix}{{\mathbf{X}}_{Iner}} &{{\mathbf{X}}_{Fric}} \\ \end{matrix} \right ]}^{\textit{T}}} \nonumber \\ & =\textbf{J}\left (\textbf{q},{\dot{\textbf{q}}},{\ddot{\textbf{q}}} \right )\cdot{\mathbf{X}} \end{align}

where $\mathbf{X}\in{\Re }^{15n \times 1}$ is the standard dynamic parameter vector.

However, it is impossible to identify all dynamic parameters in Eq. (12). When the null columns and linear correlations of the matrix $\textbf{J}\left (\textbf{q},{\dot{\textbf{q}}},{\ddot{\textbf{q}}} \right )$ are removed, Eq. (12) is rewritten as:

(13) \begin{align} \mathbf{Y}&=\overset{\scriptscriptstyle \frown }{\textbf{J}}\left (\textbf{q},{\dot{\textbf{q}}},{\ddot{\textbf{q}}} \right )\cdot \overset{\scriptscriptstyle \frown }{\mathbf{X}} \end{align}

where $\overset{\scriptscriptstyle \frown }{\mathbf{X}}$ is the base dynamic parameters.

Thus, the dynamic parameters can be obtained by CVX [Reference Han, Wu, Liu and Xiong13], as follows:

(14) \begin{align} \arg \min &\left ({{\left \|{\mathbf{Y}}- \overset{\scriptscriptstyle \frown }{\textbf{J}}\cdot \overset{\scriptscriptstyle \frown }{\mathbf{X}} \right \|}^{2}} \right ) \end{align}

4.2. Nonlinear friction parameter estimation

To identify the friction coefficients with velocity and load dependencies, we analyze the components of the joint torque, thereby calculating the friction torque separately. To this end, we command the joint in velocity mode by constructing a motion trajectory in the joint space [Reference Dong, Xu, Zhou, Zhu and Yu5]. When the joint rotates forward and backward at the same velocity, the friction torque is symmetrical about the zero point of the coordinate, indicating opposite torque with equal absolute values. Thus, we can get the friction torque by subtracting the corresponding joint torque of the robot, as follows:

(15) \begin{align}{{\boldsymbol{\tau }}_{Joint}}\left (t\right )-{{\boldsymbol{\tau }}_{Joint}}\left (T-t\right )&=\textbf{M}\left ( \textbf{q}\left (t\right ) \right )\ddot{\textbf{q}}\left (t\right )+\textbf{V}\left ( \textbf{q}\left (t\right ),\dot{\textbf{q}}\left (t\right ) \right )\dot{\textbf{q}}\left (t\right )+\textbf{G}\left ( \textbf{q}\left (t\right ) \right )\nonumber \\ &+{{\textbf{J}}^{\textit{T}}}\left ( \textbf{q}\left (t\right ) \right ){{\textbf{F}}_{ext}}+{{\boldsymbol{\tau }}_{\boldsymbol{f}}}\left (\dot{\textbf{q}}\left (t\right ),{\boldsymbol{\tau }}_{ext}\right )-\left (\right .\textbf{M}\left ( \textbf{q}\left (T-t\right ) \right )\ddot{\textbf{q}}\left (T-t\right )\nonumber \\ &+\textbf{V}\left ( \textbf{q}\left (T-t\right ),\dot{\textbf{q}}\left (T-t\right ) \right )\dot{\textbf{q}}\left (T-t\right )+\textbf{G}\left ( \textbf{q}\left (T-t\right ) \right )\nonumber \\ &+{{\textbf{J}}^{\textit{T}}}\left ( \textbf{q}\left (T-t\right ) \right ){{\textbf{F}}_{ext}}+{{\boldsymbol{\tau }}_{\boldsymbol{f}}}\left (\dot{\textbf{q}}\left (T-t\right ),{\boldsymbol{\tau }}_{ext}\right )\left .\right )\nonumber \\ &=2{{\boldsymbol{\tau }}_{\boldsymbol{f}}}\left (\dot{\textbf{q}}\left (t\right ),{\boldsymbol{\tau }}_{ext}\right ) \end{align}

Subsequently, the proposed friction model (8) is utilized to fit the friction curves obtained by Eq. (13), and the friction coefficients can be also identified offline by the Levenberg–Marquardt nonlinear least square (L-MNLS) method.

4.3. Accurate dynamic model

This model is to enhance the versatility of the dynamic model and ensure the accuracy of the estimated torques. On the basis of the parameter identification, we can give an accurate dynamic model in the following manner:

(16) \begin{align}{\boldsymbol{\tau }}_{Joint}=\left \{ \begin{matrix} \textbf{M}(\textbf{q})\ddot{\textbf{q}}+\textbf{V}\left ( \textbf{q},\dot{\textbf{q}} \right )\dot{\textbf{q}}+\textbf{G}(\textbf{q})+{\boldsymbol{\tau }}_{f_{sim}} & \textbf{F}_{ext} = 0 \\ \textbf{M}(\textbf{q})\ddot{\textbf{q}}+\textbf{V}\left ( \textbf{q},\dot{\textbf{q}} \right )\dot{\textbf{q}}+\textbf{G}(\textbf{q})+{{\textbf{J}}^{\textit{T}}}(\textbf{q}){{\textbf{F}}_{ext}}+{\boldsymbol{\tau }}_{f_{com}} & \textbf{F}_{ext} \neq 0 \\ \end{matrix} \right . \end{align}

where ${\boldsymbol{\tau }}_{f_{sim}}$ represents the simplified joint friction model and ${\boldsymbol{\tau }}_{f_{com}}$ is the comprehensive friction model.

Letting

(17) \begin{align} {\boldsymbol{\tau }}_{f_{sim}}&={\boldsymbol{f}}_{v_1}{\dot{\textbf{q}}} +{\boldsymbol{f}}_{v_2}{sgn\left ({{{\dot{\textbf{q}}}}}\right )}{\dot{\textbf{q}}^2} +{\boldsymbol{f}}_{v_3}{\dot{\textbf{q}}^3}+{\boldsymbol{f}}_{c}sgn\left ({{{\dot{\textbf{q}}}}}\right )+{\boldsymbol{f}}_{b} \end{align}

(18) \begin{align}{\boldsymbol{\tau }}_{f_{com}} =& {\boldsymbol{f}}_{v_1}{\dot{\textbf{q}}} +{\boldsymbol{f}}_{v_2}{sgn\left ({{{\dot{\textbf{q}}}}}\right )}{\dot{\textbf{q}}^2} +{\boldsymbol{f}}_{v_3}{\dot{\textbf{q}}^3}+{\boldsymbol{f}}_{c}sgn\left ({{{\dot{\textbf{q}}}}}\right )\left (1-e^{-\left ({\dot{\textbf{q}}}/\boldsymbol{f}_{c1}\right )^2}\right ) \nonumber \\ &+{\boldsymbol{f}}_{s}sgn\left ({{{\dot{\textbf{q}}}}}\right )e^{-\left ({\dot{\textbf{q}}}/\boldsymbol{f}_{s1}\right )^2}+\left ({{\boldsymbol{f}}_{{{\boldsymbol{\tau }}_{ext}}}}\vert{{\boldsymbol{\tau }}_{ext}}\vert +{{\boldsymbol{f}}_{{{\boldsymbol{\tau }}_{ext1}}}}{{e}^{-\boldsymbol{f}_{{{\boldsymbol{\tau }}_{ext2}}}\vert{{\boldsymbol{\tau }}_{ext}}\vert /{\dot{\textbf{q}}}}}\right ) sgn\left ({{{\dot{\textbf{q}}}}} \right ) \end{align}

5.1. Excitation trajectory

In order to discover excitation trajectories that offer enhanced signal-to-noise ratio and analytical derivatives, one can define the FS-based periodic trajectory for the $i$ -th joint as a function of the time $t$ , and it is made up of a finite sum of $N$ harmonic sine and cosine functions in the following manner:

(19) \begin{align}{{q}_{i}}\left ( t \right )={{q}_{i,0}}+\sum \limits _{k=1}^{N}{\left ({\frac{{a}_{i,k}}{k{{\omega }_{fi}}}}\sin \left ( k{{\omega }_{fi}}t \right )-{\frac{{b}_{i,k}}{k{{\omega }_{fi}}}}\cos \left ( k{{\omega }_{fi}}t \right ) \right )} \end{align}

where $\omega _{fi}$ is the fundamental frequency, $q_{i,0}$ is the initial position for the $i$ -th joint, and $a_{i,k}$ and $b_{i,k}$ are the Fourier coefficients, which represent the amplitudes of the cosine and sine functions, respectively. Note that we can obtain the velocity and acceleration of each joint through differentiating Eq. (12).

5.2. Proposed optimization criterion

In the literature, several ways are applied to optimize excitation trajectories, including minimizing the Cond $({\cdot})$ [Reference Bonnet, Fraisse, Crosnier, Gautier, González and Venture3, Reference Calafiore, Indri and Bona4, Reference Fu, Pan, Spyrakos-Papastavridis, Lin, Zhou, Chen and Dai6, Reference Liu, Ma, Chen, Cao, Fu and qi Li22, Reference Presse and Gautier28, Reference Wu, Zhu, Wang and Liu42], $\frac{1}{\prod{({\cdot})}}$ [Reference Dong, Xu, Zhou, Zhu and Yu5, Reference Jin and Gans18], and Log{det $({\cdot})$ } [Reference Calafiore, Indri and Bona4, Reference Swevers, Ganseman, Tukel, De Schutter and Van Brussel33, Reference Vantilt, Aertbeliën, De Groote and De Schutter36]. Considering the benefits of three different optimization criteria [Reference Luo, Bao, Du, Hu, Liu and Yuan24], we proposed a new optimization criterion, where add $\text{Log}\left \{\text{det}({\cdot})\right \}$ as the optimization index feedback when minimizing $\text{Cond}({\cdot})$ . In the meantime, the $\frac{1}{\prod{({\cdot})}}$ is introduced to quickly find a feasible initial guess in the shortest duration time as follows:

(20) \begin{align} &\arg \min \left (\frac{1}{\prod \nolimits _{g=1}^{p}{{\overset{\scriptscriptstyle \frown }{\textbf{J}}_{g}^{s}}}} \right ) \qquad \qquad \xrightarrow [{q}_{i,0}]{{a}_{i,k},{b}_{i,k}} \qquad &\arg \min \left (\textrm{Cond}\left ( \overset{\scriptscriptstyle \frown }{\textbf{J}}\right ) \right ) \leq \Delta _{1} \nonumber \\ & s.t.\left \{ \begin{matrix}{{\textbf{q}}_{\min }}\le{\textbf{q}}\left ( \Upsilon \right )\le{{\textbf{q}}_{\max }} \\[3pt] \left |{\dot{\textbf{q}}}\left ( \Upsilon \right ) \right |\le{{{\dot{\textbf{q}}}}_{\max }} \\[3pt] \left |{\ddot{\textbf{q}}}\left ( \Upsilon \right ) \right |\le{{{{\ddot{\textbf{q}}}}}_{\max }} \\[3pt] \left |{{\textbf{V}}_c} \right | \le{{{\textbf{V}}_c}_{\max }}\\[3pt] {\textbf{q}}\left ({{\Upsilon }_{s}} \right )={\textbf{q}}\left ({{\Upsilon }_{e}} \right )=0 \\[3pt] {\dot{\textbf{q}}}\left ({{\Upsilon }_{s}} \right )={\dot{\textbf{q}}}\left ({{\Upsilon }_{e}} \right )=0 \\[3pt] {\ddot{\textbf{q}}}\left ({{\Upsilon }_{s}} \right )={\ddot{\textbf{q}}}\left ({{\Upsilon }_{e}} \right )=0 \\ \end{matrix} \right . & s.t.\left \{ \begin{matrix} \text{Log}\left \{\text{det}\left (\textbf{J}_{log}\right )\right \} \leq \Delta _{2} \\[3pt] {{\textbf{q}}_{\min }}\le{\textbf{q}}\left ( \Upsilon \right )\le{{\textbf{q}}_{\max }} \\[3pt] \left |{\dot{\textbf{q}}}\left ( \Upsilon \right ) \right |\le{{{\dot{\textbf{q}}}}_{\max }} \\[3pt] \left |{\ddot{\textbf{q}}}\left ( \Upsilon \right ) \right |\le{{{{\ddot{\textbf{q}}}}}_{\max }} \\[3pt] \left |{{\textbf{V}}_c} \right | \le{{{\textbf{V}}_c}_{\max }}\\[3pt] {\textbf{q}}\left ({{\Upsilon }_{s}} \right )={\textbf{q}}\left ({{\Upsilon }_{e}} \right )=0 \\[3pt] {\dot{\textbf{q}}}\left ({{\Upsilon }_{s}} \right )={\dot{\textbf{q}}}\left ({{\Upsilon }_{e}} \right )=0 \\[3pt] {\ddot{\textbf{q}}}\left ({{\Upsilon }_{s}} \right )={\ddot{\textbf{q}}}\left ({{\Upsilon }_{e}} \right )=0 \\ \end{matrix} \right . \end{align}

where $\textbf{J}_{log}=\overset{\scriptscriptstyle \frown }{\textbf{J}}^{T}\overset{\scriptscriptstyle \frown }{\textbf{J}}$ , ${\textbf{q}}_{min}$ , ${\textbf{q}}_{max}$ , ${\dot{\textbf{q}}}_{max}$ , ${\ddot{\textbf{q}}}_{max}$ , and ${{\boldsymbol{\tau }}_m}_{max}$ denote the minimum and maximum joint positions, maximum joint velocities, accelerations, and torque, respectively. ${{\textbf{V}}_c}_{\max }$ is the maximum Cartesian velocity. Unexpected behavior will result at the start and end points if ${\dot{\textbf{q}}}\left ({{\Upsilon }_{s}} \right ) \neq{\dot{\textbf{q}}}\left ({{\Upsilon }_{e}} \right )$ and ${\ddot{\textbf{q}}}\left ({{\Upsilon }_{s}} \right ) \neq{\ddot{\textbf{q}}}\left ({{\Upsilon }_{e}} \right )$ , for example, at the start and end of the motion, the robot experiences sudden changes in velocity and acceleration, leading to vibrations.

Hence, with respect to the problem of determining an optimal trajectory, the optimization procedure can be synopsized in Table I.

6.1. Experimental setup

A $\text{Franka Emika}$ robot system with seven revolute joints is used as the experiment setup. Figure 3 shows the initial pose of the robot with an external load of $3.035$ $kg$ . The position variable of joint $i$ is defined as $q_{i}$ , the positions of joints $4$ and $6$ in the present relationship are acquired from the corresponding measured ones using the relations, $q_{4}=q_{4Panda}-\pi /2$ and $q_{6}=q_{6Panda}+\pi /2$ , and the temperature sensors are respectively installed on the surface of the corresponding joint to measure the joint temperature in real time. The aims of experiments are mainly twofold: 1) the accurate parameter identification for dynamic model with the nonlinear friction and 2) the performance between the overall required time and the quality for the excitation trajectory optimization. As seen in Figure 3, the kinematic model of the robot is built by utilizing the spatial representation of its geometric twists, where the kinematic parameters are provided in Table II. The experiments are conducted through the Franka Control Interface (FCI) via the libfranka interface. For the experiments and verification, we utilize MATLAB 2019a to run the complete computer programs. The optimization has been done over Inter(R) Core(TM) i9-12900H CPU @ $2.5$ GHZ, $16$ GB RAM.

Table I. Optimal trajectory design process.

Table II. Kinematic parameters of the robot.

Figure 3. Experimental setup of Franka Emika robot system.

6.2. Data acquisition and processing

As described in Section 5, we choose the excitation trajectories with a period of $20$ s. After that, the robot joint is commanded to move in a designed trajectory using velocity mode; the measurement time is $100$ s, which equals five periods of the excitation trajectory. The motor-side encoder can track sine and cosine curves for each robot joint, so positions and velocities can be gained. Additionally, the force/torque sensors can provide the joint torques. Since the excitation trajectory is periodic, we first improve the signal-to-noise ratio by averaging data. Afterward, the measured joint positions and velocities undergo additional high-frequency noise reduction through the application of a zero-phase low-pass filter. Next, the joint accelerations can be computed via the central difference algorithm and then smoothed through RLOESS. At the same time, we can also filter the averaged joint torques by using RLOESS to remove the high-frequency noise and its ripples.

6.3. Excitation trajectory optimization

To make a comprehensive comparison between the proposed optimization criterion and others in terms of robust, fast, and accurate identification, we analyze the computation time of the trajectory optimization process and evaluate the quality of optimization results.

According to Section 5, five groups of optimization excitation trajectories were generated for the evaluation of various optimization criteria, in which the condition number is typically optimized to be below 100 for the seven DoFs serial robot [Reference Liu, Ma, Chen, Cao, Fu and qi Li22], that is, $\Delta _{1}=100$ . Simultaneously, the single criterion $\frac{1}{\prod{({\cdot})}}$ is optimized to obtain the matched $\text{Log}\left \{\text{det}({\cdot})\right \}$ of the corresponding optimization value as the goal weight $\Delta _{2}$ of the proposed optimization criterion. Furthermore, we record the computational time during the optimization process and the matched $\textrm{Cond}({\cdot})$ and $\text{Log}\left \{\text{det}({\cdot})\right \}$ for each minimum value of the objective function, as listed in Table III. The result shows that the average execution time for $C_1$ is approximately $29.95\%$ of $C_2$ , and the matched $C_4$ for $C_1$ is better than that for the value of $C_2$ when only $C_2$ is used as the optimization criterion, although the matched $C_2$ after optimization of $C_1$ was close to the corresponding value of $C_2$ when only $C_2$ is used as the optimization criterion. Additionally, the matched $C_4$ is close to that when using only $C_3$ as the optimization criterion, but the matched $C_2$ is better than that of the latter.

Table III. Comparison of different optimization criteria.

1. 1. $C_1$ , $C_2$ , $C_3$ , and $C_4$ are the proposed optimization criterion, $\text{Cond}({\cdot})$ , $\frac{1}{\prod{({\cdot})}}$ , and $\text{Log}\left \{\text{det}({\cdot})\right \}$ , respectively.

2. 2. To perform optimization of $C_1$ , we set $\Delta _1 = 100$ and $\Delta _2 = 265$ , which are derived from the objective function values of $C_2$ and $C_3$ , respectively.

Reference [Reference Jia, Zhang, Li, Gao, Zang and Zhao17] has shown that reducing the values of $C_2$ and $C_4$ can effectively improve the anti-interference capability of the identification results against noises and the accuracy of the estimated parameters, respectively. Hence, compared to other optimization criteria, the optimal trajectory selected through the proposed criterion can achieve shorter overall trajectory optimization time while maintaining parameter identification accuracy and robustness. In a word, a better compromise approach can effectively balance the quality of identification results and the total required time for the optimization process, thereby offering significant advantages in achieving fast, accurate, and robust identification.

6.4. Identification experiment results analysis

6.4.1. Identification experiment of the dynamic parameters in a no-load scenario

The main objective of this section is to get the accurate dynamic parameters by the experiment. Furthermore, we can also give a fair and intuitive comparison by the prediction torque errors in the static framework.

First, each joint of the robot runs continuously for a while at the trajectory with the minimal condition number, and the corresponding joint temperature reaches a relatively dynamic equilibrium state, in which the surface temperature of the joints is $37.8$ – $39.1$ $^{\circ }$ C, $41.4$ – $42.9$ $^{\circ }$ C, $41.0$ – $42.7$ $^{\circ }$ C, $42.0$ – $43.5$ $^{\circ }$ C, $38.7$ – $40.0$ $^{\circ }$ C, $39.3$ – $40.4$ $^{\circ }$ C, and $37.5$ – $38.1$ $^{\circ }$ C under the environment temperature of $24.5$ – $25.8$ $^{\circ }$ C, respectively. Subsequently, according to the filtered joint positions, velocities, and accelerations, we can acquire the corresponding regressor matrix $\overset{\scriptscriptstyle \frown }{\textbf{J}}\left (\textbf{q},{\dot{\textbf{q}}},{\ddot{\textbf{q}}} \right )$ . Combining with the measured joint torques after filtering, the final dynamic parameter results including the nonlinear friction can be estimated via Eq. (14), as listed in Table IV. We can evaluate the effectiveness of the identified results through the two evaluation criteria of root-mean-square (RMS) error $\varepsilon _{RMS}$ and correlation coefficient $\sigma$ [Reference Dong, Xu, Zhou, Zhu and Yu5].

Table IV. The identification values of dynamic parameters.

Figure 4. Comparison between the measured and estimated joint torques for the parameter identification.

Thus, Figure 4 compares the estimated joint torques, which are computed by the identified dynamic parameters, with the measured joint torques. The calculated RMS errors of each joint, which quantify the difference between the measured and estimated torques, and the correlation coefficients of the seven joints are shown in Table V. We can notice that the torque estimations are almost perfectly superimposing the measured ones for joints $2\dots 4$ and $6$ , and the estimated errors of other joints show some discrepancies. When the joint velocity is close to zero, the errors are larger due to the friction modeling errors at low speed: this behavior is typical of joint friction. Simultaneously, the proposed dynamic model with a nonlinear friciton exhibits smaller residuals for a comparison with others.

6.5. Verification experiment results analysis

6.5.1. Verification experiment of accuracy dynamic model without the external load

To verify the accuracy of the dynamic models in the absence of an external load, the robot was operated under a new excitation trajectory, which is different from the excitation trajectory used for parameter identification. Using the identified dynamic parameters, the joint torques along the same exciting trajectory are calculated via Eq. (10), Eq. (16), and the conventional model [Reference Fu, Pan, Spyrakos-Papastavridis, Lin, Zhou, Chen and Dai6]. The measured and the predicted torques are compared, and the results of this comparison are reported in Figure 6 and Table VI, showing good result. In Table VI, we can observe that the predicted torques calculated by the dynamic model (16) exhibit smaller residuals in all the joints. That is, our dynamic model has better friction characteristics than the model from ref. [Reference Fu, Pan, Spyrakos-Papastavridis, Lin, Zhou, Chen and Dai6]. The simplified friction model can accurately describe the robot dynamics in normal motion without external load, despite the coupling of joint frictions in velocity–load-dependent items due to its ability to capture the main contribution of friction.

Table V. RMS errors and correlation coefficients for the identification experiment.

Table VI. RMS errors and correlation coefficients for the verification experiment.

Figure 5. Dependency of the total friction torque on load and velocity.

6.5.2. Verification experiment of nonlinear friction model with different loads

In order to verify the accuracy and reliability of the nonlinear friction model (8), three different external payloads, including $1.015$ $kg$ , $2.020$ $kg$ , and $2.525$ $kg$ , are attached to the end-effector of the robot, respectively. Simultaneously, three kinds of friction models of refs. [Reference Gao, Yuan, Han, Wang and Wang7, Reference He, Wang, Bao, Yuan, Du, Ma and Wan14, Reference Madsen, Rosenlund, Brandt and Zhang25] are used for comparison. Taking joint $3$ as an example, we command a single robot joint with different angular velocities via the constructed motion trajectory and the other joints remain stationary. Subsequently, we can obtain the actual joint friction torques according to the actual joint velocities $\dot{\textbf{q}}$ and link side torques ${\boldsymbol{\tau }}_{Joint}$ , and the estimated friction torques are computed by four kinds of friction models. Note that all the parameters in different friction models can be obtained under the identification methodology. The comparisons between the actual and estimated friction values are depicted in Figure 7, which shows that the proposed nonlinear friction model can better characterize the friction and has also been successfully validated in the other six joints.

Figure 6. Comparison between the measured and estimated torques for the parameter verification.

Figure 7. RMS errors between the actual and estimated friction values for different models.