Lie-theory-based dynamic model identification of serial robots considering nonlinear friction and optimal excitation trajectory

Abstract

Accurate dynamic model is essential for the model-based control of robotic systems. However, on the one hand, the nonlinearity of the friction is seldom treated in robot dynamics. On the other hand, few of the previous studies reasonably balance the calculation time-consuming and the quality for the excitation trajectory optimization. To address these challenges, this article gives a Lie-theory-based dynamic modeling scheme of multi-degree-of-freedom (DoF) serial robots involving nonlinear friction and excitation trajectory optimization. First, we introduce two coefficients to describe the Stribeck characteristics of Coulomb and static friction and consider the dependency of friction on load torque, so as to propose an improved Stribeck friction model. Whereafter, the improved friction model is simplified in a no-load scenario, a novel nonlinear dynamic model is linearized to capture the features of viscous friction across the entire velocity range. Additionally, a new optimization algorithm of excitation trajectories is presented considering the benefits of three different optimization criteria to design the optimal excitation trajectory. On the basis of the above, we retrieve a feasible dynamic parameter set of serial robots through the hybrid least square algorithm. Finally, our research is supported by simulation and experimental analyses of different combinations on the seven-DoF Franka Emika robot. The results show that the proposed friction has better accuracy performance, and the modified optimization algorithm can reduce the overall time required for the optimization process while maintaining the quality of the identification results.

1. Introduction

Robotic dynamics has broad applications in industrial production scenarios [5, 30]. Model-based torque predictions are crucial in developing serial industrial robots to get superior performance in terms of accuracy and speed [25]. Hence, the accurate dynamic model is important and necessary for serial robots, which are used to design advanced control laws with superior performance and offline programming with task optimization [34]. However, due to many uncertainties in manufacturing and application situations, the dynamic parameters are typically unknown or only partially provided [2]. Experimental identification of dynamic parameters is a relatively low-cost and reliable approach for obtaining an accurate model, which is divided into five parts: modeling, excitation trajectory optimization, data acquisition and processing, parameter identification, and model validation [17, 41].

Over the years, a lot of related studies have been conducted about dynamic models of a robot in the context of parameter identification, which usually include two parts: the inertial dynamics [1, 6, 9, 11, 18, 19, 21, 32, 38, 43] and the joint frictions [7, 8, 10, 12, 14, 16, 20, 26, 31, 37, 39]. In the former, the Lagrangian approach or the recursive Newton–Euler equation can effectively construct the inertial dynamic model. Accurate robot dynamic modeling is challenging for the latter due to inherent joint friction effects. So far, the classical static friction model, which combines Coulomb friction with linear viscous friction, is usually applied in parameter identification. Although the classical static model accounts for the major friction features at medium and high velocities, but it fails to fully capture the behavior during startup and low velocities in practice [5], resulting in the derivation between the estimated friction torques and its real ones because the real friction has strong nonlinearity. That is, the accuracy of dynamic modeling will be inevitably affected. Several schemes have been employed so that the issue mentioned above can be effectively solved. Zhang et al. [45] built the comprehensive dynamics of industrial robots and calibrated the corresponding parameters by a nonlinear method. Han et al. [13] proposed an iterative algorithm of the dynamic model identification for serial robots to estimate the nonlinear frictions in a unified way with linear regression. Wu et al. [40] identified the inertial parameters without the assumed friction model via utilizing the baseplate sensor, and the friction coefficient of the joint can be estimated by fitting the joint friction model after subtracting the inertial torques from the joint torques. Dong et al. [5] modified the Tustin friction model, in which two parameters were used to describe the Stribeck features of static and Coulomb frictions, and the friction and inertial parameters were identified based on the LS with SOS algorithm. Roveda et al. [29] exploited a 6D virtual sensor to quantify the joint friction effects and proposed a Bayesian optimization-based algorithm to enhance the impedance control performance by tuning the model-based friction compensator parameters. Huang et al. [15] proposed an IHLS-BPNN approach of parameter identification for robot dynamics, where the Stribeck friction model was substituted for the Coulomb-viscous friction model. For the methods mentioned above, the velocity-dependent effects on the friction, including the Stribeck effect, are only taken into account. In a word, the accurate nonlinear friction model has not yet been built because the friction coefficients related to the temperature–load-dependent were neglected.

Since an exciting trajectory fundamentally affects the robustness and accuracy of dynamic model identification results, optimizing the excitation trajectory is also significant among experimental identification parts. To simplify the analysis and conditions of the measurements, we typical select finite-Fourier-series (FS)-based periodic trajectory as the excitation trajectory [33, 34]. Generally, two commonly used optimality criteria are employed to determine the optimal trajectory parameters. One criterion focuses on minimizing the condition number of the regression matrix [3, 4, 6, 22, 28, 42], which effectively mitigates the impact of noise errors on parameter identification; another is minimizing the logarithm of the determinant of the Fisher information matrix [4, 23, 33, 36], resulting in parameter estimation with smaller uncertainty bounds. To reduce the computational complexity and optimization time, Jin et al. [18] given a new optimization criterion by employing Hadamard’s inequality. However, the identification results are easily interfered by system noise because the number of matching conditions is greatly increased. Recently, combinatorial optimization based on different criteria has been conducted concerning the mutual restriction of the objective function values for the single optimization criterion. Jia et al. [17] gave an improved optimization method of the excitation trajectories through introducing the feedback index as a constraint target. However, generating feasible initial parameters for multi-degree-of-freedom (DoF) robots is hard due to the stochastic effect of the starting iteration values in the optimization. Besides, the calculation time is greatly enhanced in finding the optimal trajectory parameters.

Figure 1. Schematic of the entire identification procedure.

As an alternative mathematical tool, formulas for Lie-theory-based inertial dynamics analysis proved to be a simpler and more computationally efficient approach compared with one based on the Denavit–Hartenberg law [6, 44]. Motivated by the above discussion, we propose a dynamic modeling scheme for multi-DOF serial robots with/without external loads involving nonlinear friction and excitation trajectory optimization. The significant contributions of this study are as follows:

1. 1) We propose an improved Stribeck friction model, in which the dependency of friction is related to velocity and load when the joint temperature reaches a relatively dynamic equilibrium.

2. 2) A novel Lie-theory-based nonlinear dynamic model is linearized to capture the features of viscous friction across the entire velocity range and introduce an accurate dynamic modeling scheme.

3. 3) We put forward a novel excitation trajectory optimization algorithm based on the advantages of three different optimization criteria.

The remainder of this article is arranged as follows. Section 2 reviews the overview of the proposed dynamic identification. Section 3 establishes a Lie-theory-based accurate dynamic model for robots in joint space, which includes the improved friction model related to the velocity–load dependency. Later, Section 4 describes the whole identification procedure of dynamic parameters. In Section 5, a new optimality criterion is presented, thereby finding the optimization excitation trajectory in the shortest duration. Additionally, Section 6 provides simulation and experiments to validate the performance of improved optimization methods and an accurate dynamic framework. Finally, Section 7 gives the conclusions and future works of the paper.

2. Overview

The overall approach for the dynamic identification is illustrated in Figure 1, as far as the authors know, (1) a classical static friction model is usually adopted to the process of parameter identification, and nonlinearity of the frictions related to the factors including velocity and load dependency is seldom treated, causing the deviations from the real ones. Therefore, it’s crucial to give a accurate dynamic model for serial robots with/without a external payload across the entire velocity range. (2) The excitation trajectories are a basal and significant part of improving identification accuracy, but the most of previous studies ignore its calculation time when improving the optimization quality. Thus, it’s essential to find a reasonable balance between computational efficiency and optimization quality. This study presents a Lie-theory-based dynamic modeling scheme of multi-DoF serial robots involving nonlinear friction and excitation trajectory optimization.

3. Lie-theory-based robot dynamics analysis

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 [27, 44]. 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 [6] 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 [35], 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 [35], 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}$ [5]. 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. [7, 14]), 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. Model identification

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 [13, 15], 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. [6]), 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 [13], 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 [5]. 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. Optimal excitation trajectory design

In parameter identification, appropriate excitation trajectories should be designed and optimized to fully simulate the mechanical properties of a real robot system. To this end, we proposed an improved optimization criterion to find an optimal excitation trajectory of each joints.

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})$ [3, 4, 6, 22, 28, 42], $\frac{1}{\prod{({\cdot})}}$ [5, 18], and Log{det $({\cdot})$ } [4, 33, 36]. Considering the benefits of three different optimization criteria [24], 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. Experiment and verification

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 [22], 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 [17] 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$ [5].

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.4.2. Estimation experiment of the nonlinear friction coefficients

To identify $10$ coefficients in the improved friction model (8) related to the velocity–load-dependent, an external payload of $3.035$ $kg$ is attached to the end-effector of the robot to increase the maximum possible load torque. According to the constructed motion trajectory, we command a single robot joint with different angular velocities in velocity mode, and all other joints remain stationary. Then, we can obtain the actual joint velocities $\dot{\textbf{q}}$ and link side torques ${\boldsymbol{\tau }}_{Joint}$ through the libfranka interface. The joint torques corresponding to points with the same position but opposite velocities are selected, the velocity–load-dependent frictions are obtained according to Eq. (15), and the optimum parameter of the model (8) can be found by the L-MNLS algorithm (using MATLAB function, lsqnonlin). That is, the friction coefficients of each joint can be estimated at different velocities, in which the velocity of each joint is set to $0.1$ , $0.2$ , $0.4$ , $0.5$ , $0.8$ , $1$ , $1.2$ , $1.5$ , $1.6$ , $1.8$ , $2$ , $2.2$ , $2.5$ , $2.8$ , $3$ , $3.5$ , $4$ , $4.5$ , $5$ , $6$ , $8$ , $10$ , $12$ , $15$ , $16$ , $20$ , $25$ , and $30$ $(\text{all in} \;{}^{\circ }/s)$ . The fitting diagram from the identification for the dependency of friction on velocity–load is depicted in Figure 5, the calculated $\textrm{RMS}$ errors of each joint are $0.0105$ , $0.0264$ , $0.0213$ , $0.0099$ , $0.0550$ , $0.0503$ , and $0.0033$ (all in Nm), respectively, and the corresponding correlation coefficients are $0.9437$ , $0.9572$ , $0.9853$ , $0.9703$ , $0.9507$ , $0.9776$ , and $0.9567$ .

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 [6]. 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. [6]. 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. [7, 14, 25] 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.

7. Conclusion

In our work, a Lie-theory-based accurate dynamic modeling scheme is given for multi-DOF serial robots with/without external loads, where we propose the improved Stribeck friction model involving the nonlinear dependence of friction on the velocity–load and introduce a novel linearizable nonlinear dynamic model. On the basis of the interaction between different optimization criteria, we modify the optimization technique for the design of optimal excitation trajectories used in dynamic identification. Finally, several experiments are carried out on the seven DoFs Franka Emika robot, and the experimental results reveal twofold: (1) the proposed dynamics scheme has better fitting accuracy and higher versatility and (2) the optimal excitation trajectory generated via the proposed criterion requires shorter optimization time while ensuring the quality of identification results compared to others, which can provide advantages for fast, robust, and accurate identification.

In the next work direction, the time-varying temperature-dependent friction phenomena will be researched for fine modeling and compensation. Simultaneously, the developed friction will be seamlessly extended to the dynamic friction model and applied to robot dynamics in a unified way. Concurrently, there is a need for further exploration at the robot planning level in conjunction with advanced intelligent control theories.