Feed-forward control of elastic-joint industrial robot based on hybrid inverse dynamic model

Abstract

A novel feedforward control method of elastic-joint robot based on hybrid inverse dynamic model is proposed in this paper. The hybrid inverse dynamic model consists of analytical model and data-driven model. Firstly, the inverse dynamic analytical model of elastic-joint robot is established based on Lie group and Lie algebra, which improves the efficiency of modeling and calculation. Then, by coupling the data-driven model with the analytical model, a feed-forward control method based on hybrid inverse dynamics model is proposed. This method can overcome the influence of the inaccuracy of the analytical inverse dynamic model on the control performance, and effectively improve the control accuracy of the robot. The data-driven model is used to compensate for the parameter uncertainties and non-parameter uncertainties of the analytical dynamic model. Finally, the proposed control method is proved to be stable and the multi-domain integrated system model of industrial robot is developed to verify the performance of the control scheme by simulation. The simulation results show that the proposed control method has higher control accuracy than the traditional torque feed-forward control method.

Introduction

Industrial robots are widely used in modern industrial manufacture, such as paint spraying, welding, grinding ¹ milling, ² drilling, ³ and so on. In recent years, with the rise of intelligent manufacturing technology and the close combination of the artificial intelligence technology and industrial robots, the application of robots has been expanding to all walks of life. Since the rapid development of modern enterprises requires higher production efficiency, product quality and faster product iteration speed, the industrial robots have entered the development stage of high speed, high precision, heavy load, lightweight, and intelligent.

The traditional control method of joint industrial robot is PID ⁴ control algorithm which is not based on model. Each joint is regarded as an independent second-order linear system, and other coupled joint dynamics is regarded as external interference. This control algorithm is simple and easy to design. In some cases (such as low speed and stable working conditions) it can meet the accuracy requirements of the control system. However, the industrial robot is a complex and multi-domain coupling nonlinear system, the traditional PID control is usually hard to meet the high-speed and high-precision control requirements. With the development of industrial robot control technology, many control algorithms have been proposed such as fuzzy control, ⁵ optimal control, ⁶ dynamic surface control, ⁷ sliding mode control,^8,9 fraction order control, ¹⁰ robust neural control, ¹¹ torque feedforward control, ¹² and so on. Among them, torque feedforward control is the most intuitive nonlinear control method. It is an effective control method which relies on the robot inverse dynamics, so sometimes it is also called inverse dynamic control. Compared with the model-free control methods, the inverse dynamic control usually can achieve better control performance, such as high tracking accuracy, low feedback gain coefficient, high adaptability to compliance control, and low energy consumption.

Robot inverse dynamic control technology firstly needs to establish inverse dynamic model to predict the feedforward torque. In most researches of the robot inverse dynamic control, the robot is regarded as rigid body. In effect, the synchronous conveyor with significant flexible characteristics and RV harmonic reducer are commonly used in the joints of industrial robot. The joint flexibility effect is significant, which may lead to mechanical resonance and affect the stability and tracking accuracy of the servo control system. Compared with the rigid-joint robot, the elastic-joint robot presents serious problems such as complexity, high nonlinearity, and large amount of calculation. The research results show that the difficulty of the inverse dynamic modeling of elastic-joint robot lies in the fact that the high order of robot joint angel and joint driving torque under the given trajectory are needed. ¹³ To deal with the elastic joints in the inverse dynamic modeling, Hopler and Thümmel ¹³ propose a symbolic method base on Lagrange theory to deduce the high-order dynamic equations. However, the method has inherent drawbacks: depending on the chosen coordinates and complex computation for symbolic operation in high-order derivation. Guarino Lo Bianco ¹⁴ think that Newton-Euler method is more efficient to derive the high-order kinematics and dynamics, and a recursive method based on Newton-Euler was proposed to compute the derivatives of actuator torques. More recently, Buondonno and De Luca^15,16 proposed a method to analyze the inverse dynamics of the elastic-joint robot. But in the study, it is too complex to derive the high-order kinematics and dynamics because of using the Denavit-Hartenberg (D-H) technique. In recent years, as an analysis tool of modern geometric mathematics, the application of Lie group and Lie algebra in robot kinematics and dynamics has attracted more and more attention.^17–19 Park et al. ¹⁷ firstly introduced Lie group to formulate the dynamic model of rigid-joint robot. Based on the study of the Park, Yang et al. ¹⁹ proposed a recursive Newton–Euler method based on Lie group to deal with elastic joints in the inverse dynamic analysis, but the studies were not applied to the robot inverse dynamic control.

The existence of flexible robot joints also makes the high-performance control of robot very complex and difficult. In Cai et al., ²⁰ a kind of position control method based on double speed closed loops for elastic-joint robot is proposed to reduce turning error and nonlinear dead zone. This control method is pure feedback control so that the improvement of control performance is limited. To achieve excellent control performance, it is an effective method to add inverse dynamics feedforward on the basis of feedback control, which needs accurate inverse dynamic model to predict torque. Robot dynamic model is usually derived from expert domain knowledge. This modeling method is called analytical modeling, also known as knowledge-driven, physics-based, or mechanical modeling. In Kim and Croft, ²¹ and Moberg and Hanssen ²² the elastic-joint robot torque feedforward control is based on the fully analytical dynamic model. But the accurate analytical model is difficult to be established. In order to obtain the accurate dynamic model, parameter identification method is widely used to modify the nominal dynamic model. In Swevers et al., ²³ a series of excitation trajectories are designed for the robot, and then different optimization algorithm are used to fit the parameters in the nominal dynamic model. In Schaal et al., ²⁴ the local weight mapping regression method is used to estimate the dynamic model of the robot through simulation and experiment, and the estimated dynamic model is used in the feed-forward control. However, industrial robot is a complex multi-domain coupling nonlinear system, which has parameter uncertainty and non-parameter uncertainty. Only by analytical modeling and parameters identification, the accuracy of robot dynamic model is limited. In Zeinali and Notash, ²⁵ an uncertainty estimator to compensate for the un-modeled dynamics, external disturbances, and time-varying parameters was used in the sliding control of the robot manipulator. The robot inverse dynamic model can also be established by machine learning techniques, which is called the data-driven modeling method. Machine learning is one of today’s most rapidly growing technical fields ²⁶ ; there has been growing interest in applying machine learning to draw insights gained from the data in engineering. ²⁷ In Nguyen-Tuong et al., ²⁸ the inverse dynamics model of robot is established by nonparametric regression method. The fully data-based models lack of interpretability and it may suffer from inaccuracies as well when the working condition varies. So in Guo et al., ²⁹ and Mei et al. ³⁰ the hybrid analytical and data-driven modeling method was used to establish the dynamical model of the feed drive system. Reinhart et al.^31,32 have used the hybrid analytical and data-driven modeling method to establish the kinematics model for soft robots and the inverse dynamic model for rigid robots. But the author’s research mainly focused on rigid robot and the inverse dynamic model is not applied to the control.

From the above analysis, it can be seen that the inverse dynamic control of elastic-joint robot mainly faces two problems: one is the derivation calculation of high-order kinematics and dynamics is complex and tedious, the other is that the accurate dynamic model is hard to be established. To deal with these problems, an inverse dynamic control approach for elastic-joint robot based on hybrid dynamic model is proposed in this paper. Firstly, for simplicity and efficiency, the analytical inverse dynamic model of elastic-joint robot is established based on Lie group and Lie algebra. Then, a radial basis function (RBF) neural network model is combined with the analytical inverse dynamic model and the inverse dynamic control based on the hybrid model is proposed. The control law is proved to be stable by mathematical. At last, a multi-domain integrated system model of industrial robot is established and the proposed control method is validated by simulation. The simulation results show that the proposed control method is superior to the traditional inverse dynamics feedforward control.

The rest of the paper is organized as follows. In Sect.2, some basic background of Lie group and Lie algebra is introduced. In Sect.3, the elastic-joint industrial robot inverse dynamic model is established based on Lie group. In Sect.4, the inverse dynamic control method based on hybrid analytical inverse dynamic and data-driven model is proposed and the mathematical proof is given. In Sect.5, the multi-domain integrated system model of industrial robot is established and the proposed control method is validated by simulation. Finally, the conclusions are summarized in Sect.6.

Basic background of Lie group and Lie algebra

In the analysis of robot kinematics and dynamics, the problem of position and orientation transformation is often faced. Based on the traditional D-H technique, it is hard to calculate the high-order derivation of the robot kinematics and dynamics. In effect, the transformation space of all possible rigid body motions is an element of Lie group. In this section, some basic background of Lie group and Lie algebra is introduced. The detail of related theories can be referenced to the work of Selig. ³³

Previous knowledge of Lie group and Lie algebra

Lie group is not only a group, but also a smooth manifold. The tangent space of the unit element of a Lie group is called the Lie algebra of this Lie group. All the ration matrix of rigid body is a subgroup of Lie groups:

$$SO(3)=\{\mathbf{R}\in {ℝ}^{3\times 3}|\mathbf{R}{\mathbf{R}}^T=1,\det (\mathbf{R})=1\}$$

(1)

The position and orientation of rigid body in 3D space can be described by the transformation matrix

$$\mathbf{G}=[\begin{array}{l}\mathbf{R}\\ \mathbf{0}\end{array}\begin{array}{l}\mathbf{r}\\ \mathbf{1}\end{array}]$$

(2)

where $\mathbf{R}\in \mathrm{SO}(3)$ and $\mathbf{r}\in {\mathbb{R}}^3$ . All of this transformation matrix is also a subgroup of Lie groups and it is a Special Euclidean Group:

$$\mathrm{SE}(3)=\{\mathbf{G}=\left[\begin{array}{cc}\mathbf{R}\hfill & \mathbf{r}\\ \multicolumn{1}{c}{0}& 1\end{array}\right]\in {\mathbb{R}}^{4\times 4}|\mathbf{R}\in \mathrm{SO}(3),\mathbf{r}\in {\mathbb{R}}^3\}$$

(3)

Every Lie group has its corresponding Lie algebra, which describes the local properties of Lie groups. The corresponding Lie algebra of $\mathrm{SO}(3)$ and $\mathrm{SE}(3)$ are denoted with $\mathrm{so}(3)$ and $\mathrm{se}(3)$ :

$$\hat{\mathbf{\omega }}=[\mathbf{\omega }]=[\begin{array}{c}0\hfill \\ \multicolumn{1}{c}{{\omega }_z}\\ \multicolumn{1}{c}{-{\omega }_y}\end{array}\begin{array}{c}-{\omega }_z{\omega }_y\hfill \\ \multicolumn{1}{c}{0-{\omega }_x}\\ \multicolumn{1}{c}{{\omega }_{x}0}\end{array}]\in \mathrm{so}(3)$$

(4)

$$\hat{\mathbf{\xi }}=[\begin{array}{c}\hat{\mathbf{\omega }}\hfill \\ \multicolumn{1}{c}{0}\end{array}\begin{array}{c}\mathbf{v}\hfill \\ \multicolumn{1}{c}{0}\end{array}]\in \mathrm{se}(3)$$

(5)

The transformation from Lie algebra to Lie group can be described by using the following exponential mapping:

$$\mathbf{R}=e^{\hat{\mathbf{\omega }}}={\mathbf{I}}_3+\frac{\sin (\Vert \mathbf{\omega }\Vert )}{\Vert \mathbf{\omega }\Vert }\hat{\mathbf{\omega }}+\frac{1-\cos (\Vert \mathbf{\omega }\Vert )}{{\Vert \mathbf{\omega }\Vert }^2}{\hat{\mathbf{\omega }}}^2$$

(6)

$$\mathbf{G}=e^{\hat{\mathbf{\xi }}}=\left[\begin{array}{c}{e}^{\hat{\mathbf{\omega }}}\frac{({\mathbf{I}}_3-e^{\hat{\mathbf{\omega }}})(\mathbf{\omega }\times \mathbf{v})+\mathbf{\omega }{\mathbf{\omega }}^T\mathbf{v}}{\Vert \mathbf{\omega }\Vert }\hfill \\ \multicolumn{1}{c}{0^{\mathbf{T}}1}\end{array}\right]$$

(7)

An element of a Lie group can also be identified with a linear mapping between its Lie algebra via the adjoint representation. Suppose $\mathbf{G}\in \mathrm{SE}(3)$ , its adjoint map acting on an element $\hat{\mathbf{\xi }}\in \mathrm{se}(3)$ is defined as:

$$A{d}_{\mathbf{G}}(\mathbf{\xi })=(\mathbf{G}\hat{\mathbf{\xi }}{\mathbf{G}}^{-1}{)}^{\vee }$$

(8)

Similarly, an element of a Lie algebra can also be identified with a linear mapping between its Lie algebra via the Lie bracket. If ${\hat{\mathbf{\xi }}}_{\mathbf{a}}$ and ${\hat{\mathbf{\xi }}}_{\mathbf{b}}$ are elements of $\mathrm{se}(3)$ , its adjoint representation is:

$$a{d}_{{\mathbf{\xi }}_{\mathbf{a}}}({\mathbf{\xi }}_{\mathbf{b}})=({\hat{\mathbf{\xi }}}_{\mathbf{a}}{\hat{\mathbf{\xi }}}_{\mathbf{b}}-{\hat{\mathbf{\xi }}}_{\mathbf{b}}{\hat{\mathbf{\xi }}}_{\mathbf{a}}{)}^{\vee }$$

(9)

These two mappings can be calculated in matrix forms:

$$A{d}_{\mathbf{G}}=\left[\begin{array}{c}\mathbf{R}0\hfill \\ \multicolumn{1}{c}{\hat{\mathbf{r}}\mathbf{R}\mathbf{R}}\end{array}\right]$$

(10)

$$a{d}_{\mathbf{\xi }}=\left[\begin{array}{c}\hat{\mathbf{\omega }}0\hfill \\ \multicolumn{1}{c}{\hat{\mathbf{v}}\hat{\mathbf{\omega }}}\end{array}\right]$$

(11)

In the same way, the corresponding dual operator $A{d_{\mathbf{G}}}^{*}$ and $a{d_{\mathbf{\xi }}}^{*}$ are as follows:

$${\mathrm{Ad}}_{\mathbf{G}}^{*}=(A{d}_{\mathbf{G}}{)}^T=\left[\begin{array}{c}{\mathbf{R}}^T{\mathbf{R}}^T{\hat{\mathbf{r}}}^T\hfill \\ \multicolumn{1}{c}{0{\mathbf{R}}^T}\end{array}\right]$$

(12)

$${\mathrm{ad}}_{\mathbf{\xi }}^{*}=(\mathrm{a}{d}_{\mathbf{\xi }}{)}^T=\left[\begin{array}{c}-\hat{\mathbf{\omega }}-\hat{\mathbf{v}}\hfill \\ \multicolumn{1}{c}{0-\hat{\mathbf{\omega }}}\end{array}\right]$$

(13)

Rigid-robot dynamics using Lie group

For series industrial robot, the joint between adjacent rigid links is revolve joint. Each revolve joint can be associated with a twist and the twist coordinate expressed in the reference frame is defined as follows:

$$\mathbf{\xi }=\left[\begin{array}{c}\mathbf{\omega }\hfill \\ \multicolumn{1}{c}{\mathbf{v}}\end{array}\right]=\left[\begin{array}{c}\mathbf{\omega }\hfill \\ \multicolumn{1}{c}{\mathbf{p}\times \mathbf{\omega }}\end{array}\right]$$

(14)

where $\mathbf{\omega }$ is the unit vector along the joint axis and $\mathbf{p}$ is a position vector locating at any point on the joint axis. As previously described, the twist $\mathbf{\xi }$ can also be described in a matrix form by the operator “∧,” as shown in equation (5).

The Euclidean transformation that describes the position and orientation of rigid body is an element of Lie group $\mathrm{SE}(3)$ . If ${\mathbf{\xi }}_{\mathbf{i}}$ and $q_i$ represents the twist and rotation angel associated with the $i^{\mathrm{th}}$ joint respectively, the transformation matrix from link i to link $i-1$ is:

$${\mathbf{G}}_{\mathbf{i}-1,\mathbf{i}}={\mathbf{M}}_{\mathbf{i}}{e}^{{\hat{\mathbf{\xi }}}_{\mathbf{i}}{q}_i}\in \mathrm{SE}(3)$$

(15)

where $M_i$ is the initial transformation matrix ( ${\theta }_i$  = 0). Thus, the frame fixed at link i is related to that at the base by the product:

$${\mathbf{G}}_{\mathbf{i}}={\mathbf{G}}_{0,1}{\mathbf{G}}_{1,2}...{\mathbf{G}}_{\mathbf{i}-1,\mathbf{i}}$$

(16)

Then by forward recursion (for $i=1,2...n$ ), the six dimensional generalized velocity, acceleration can be calculated as follows:

$${\mathbf{V}}_i={\mathbf{G}}_{\mathbf{i}}^{-1}\frac{d({\mathbf{G}}_{\mathbf{i}})}{\mathrm{dt}}=A{d}_{{\mathbf{G}}_{\mathbf{i}-1,\mathbf{i}}^{-1}}({\mathbf{V}}_{i-1})+{\mathbf{\xi }}_{\mathbf{i}}{\stackrel{·}{q}}_i$$

(17)

$${\mathbf{A}}_{\mathbf{i}}=\frac{d({\mathbf{V}}_i)}{\mathrm{dt}}=A{d}_{{\mathbf{G}}_{\mathbf{i}-1,\mathbf{i}}^{-1}}({\mathbf{A}}_{\mathbf{i}-1})-a{d}_{{\mathbf{\xi }}_{\mathbf{i}}\stackrel{·}{\theta }}({\mathbf{V}}_i)+{\mathbf{\xi }}_{\mathbf{i}}{\stackrel{··}{q}}_i$$

(18)

According to the theory of Newtonian rigid body dynamics, the resultant force and torque acting on the rigid body in the local coordinate frame $\{i\}$ is:

$$\left[\begin{array}{c}{\mathbf{T}}_{\mathbf{i}}\hfill \\ \multicolumn{1}{c}{{\mathbf{f}}_{\mathbf{i}}}\end{array}\right]=\left[\begin{array}{c}{\mathbf{I}}_{\mathbf{i}}{\stackrel{·}{\mathbf{\omega }}}_{\mathbf{i}}+{\mathbf{\omega }}_{\mathbf{i}}\times {\mathbf{I}}_{\mathbf{i}}{\mathbf{\omega }}_{\mathbf{i}}+{\mathbf{r}}_{\mathbf{i}}\times m_i{\mathbf{a}}_{\mathbf{ci}}\hfill \\ \multicolumn{1}{c}{m_i{\mathbf{a}}_{\mathbf{ci}}}\end{array}\right]$$

(19)

Based on the adjoint transformation property of Lie algebra, the above equation can also be simplified and expressed as follows:

$$\left[\begin{array}{c}{\mathbf{T}}_{\mathbf{i}}\hfill \\ \multicolumn{1}{c}{{\mathbf{f}}_{\mathbf{i}}}\end{array}\right]={\mathbf{E}}_{\mathbf{i}}{\mathbf{A}}_{\mathbf{i}}-{\mathrm{ad}}_{{\mathbf{V}}_i}^{*}({\mathbf{E}}_{\mathbf{i}}{\mathbf{V}}_{\mathbf{i}})$$

(20)

where ${\mathbf{E}}_{\mathbf{i}}$ is the mass matrix of link i,

$${\mathbf{E}}_{\mathbf{i}}=\left[\begin{array}{c}{\mathbf{I}}_{\mathbf{i}}-m_i{\hat{\mathbf{r}}}_{\mathbf{i}}^2{m}_i{\hat{\mathbf{r}}}_i\hfill \\ \multicolumn{1}{c}{m_i{\hat{\mathbf{r}}}_i{m}_i{\mathbf{I}}_3}\end{array}\right]$$

(21)

Furthermore, according to the interaction principle of force, the resultant force and torque can also be expressed:

$$\left[\begin{array}{c}{\mathbf{T}}_{\mathbf{i}}\hfill \\ \multicolumn{1}{c}{{\mathbf{f}}_{\mathbf{i}}}\end{array}\right]={\mathbf{E}}_{\mathbf{i}}{\mathbf{A}}_{\mathbf{i}}-{\mathrm{ad}}_{{\mathbf{V}}_i}^{*}({\mathbf{E}}_{\mathbf{i}}{\mathbf{V}}_{\mathbf{i}})={\mathbf{F}}_{\mathbf{i}}-{\mathrm{Ad}}_{{\mathbf{G}}_{\mathbf{i},\mathbf{i}+1}^{-1}}^{*}({\mathbf{F}}_{\mathbf{i}+1})$$

(22)

Therefore, by backward recursion (for $i=n,n-1...1$ ), the generalized forces can be calculated as follows:

$${\mathbf{F}}_{\mathbf{i}}={\mathrm{Ad}}_{{\mathbf{G}}_{\mathbf{i},\mathbf{i}+1}^{-1}}^{*}({\mathbf{F}}_{\mathbf{i}+1})+{\mathbf{E}}_{\mathbf{i}}{\mathbf{A}}_{\mathbf{i}}-{\mathrm{ad}}_{{\mathbf{V}}_i}^{*}({\mathbf{E}}_{\mathbf{i}}{\mathbf{V}}_{\mathbf{i}})$$

(23)

And the driving torque of each joint is:

$${\mathbf{\tau }}_{\mathbf{i}}={\mathbf{\xi }}_i^T{\mathbf{F}}_i$$

(24)

Based on the derivation of the above recurrence equation, the matrix form of the robot rigid body dynamics model can be represented as follows:

$$\mathbf{\tau }=\mathbf{M}(q)\stackrel{··}{q}+\mathbf{C}(q,\stackrel{·}{q})+\mathbf{G}(q)+{\mathbf{\tau }}_{\mathbf{L}}$$

(25)

where $\mathbf{M}(\mathbf{q})\in {\mathbb{R}}^{n\times n}$ is the inertia matrix of robot, $\mathbf{C}(\mathbf{q},\stackrel{·}{\mathbf{q}})\in {\mathbb{R}}^{n\times 1}$ is centrifugal force and Coriolis force matrix, $\mathbf{G}(\mathbf{q})\in {\mathbb{R}}^{n\times 1}$ is the gravity vector and the ${\mathbf{\tau }}_{\mathbf{L}}$ is the load torque vector.

Elastic-joint robot analytical inverse dynamic model based on Lie group

When the joint flexibility is considered, the derivation of the robot inverse dynamic model is becoming extremely complex. By traditional D-H parameters method, it is hard to deal with this kind of problem. So in this section, the analytical inverse dynamic model of elastic-joint robot is derived using Lie group and Lie algebra. The related researches can also be seen in Yang et al. ¹⁹

Elastic joint dynamic model

The dynamic model of elastic joint based of linear torsion spring is first proposed by Spong. ³⁴ In the study, the elastic joint of the robot is simplified as a linear torsion spring, and the torsion deformation is assumed to be within the range of elastic deformation. Based on the Spong’s research, the elastic-joint industrial robot is presented schematically in Figure 1.

Figure 1.

Schematic diagram of single elastic joint of robot.

For n-DOF industrial robot, its analytical dynamic model can be simplified as the elastic joint and the rigid link are connected in series in turn. The output torque of all the motors are represent as the vector ${\mathbf{\tau }}_{\mathbf{m}}$ , and the rotation angle is ${\mathbf{q}}_{\mathbf{m}}$ . $\mathbf{N}\in {\mathbb{R}}^{n\times n}$ is the diagonal matrix of transmission ratio and $\mathbf{K}\in {\mathbb{R}}^{n\times n}$ is the diagonal matrix of spring stiffness. In the link side, the rotation angle is $\mathbf{q}$ , the driving torque of the joint is $\mathbf{\tau }$ . According to the principle of torque balance, the elastic joint equation can be derived:

$${\mathbf{\tau }}_{\mathbf{m}}={\mathbf{J}}_{\mathbf{m}}{\stackrel{··}{\mathbf{q}}}_{\mathbf{m}}+{\mathbf{N}}^{-1}{\mathbf{\tau }}^{\text{'}}$$

(26a)

$$\mathbf{\tau }\text{'}=\mathbf{K}({\mathbf{N}}^{-1}{\mathbf{q}}_{\mathbf{m}}-\mathbf{q})=\mathbf{K}({\mathbf{q}}_0-\mathbf{q})$$

(26b)

$$\mathbf{K}({\mathbf{q}}_0-\mathbf{q})-\mathbf{\tau }=0$$

(26c)

In the robot control, the given command signal is usually the trajectory of end-effectors. By inverse kinematic and inverse dynamic calculating, only the rotation angel $\mathbf{q}$ and driving torque $\mathbf{\tau }$ can be obtained, so the above equation should be further deducted. Combined with the equation (25), the inverse dynamic equation of elastic-joint robot can be calculated as follows,^15,16 where ${\mathbf{f}}_0$

$$\mathbf{u}={\mathbf{\tau }}_{\mathbf{m}}={\mathbf{J}}_{\mathbf{m}}\mathbf{N}[{\mathbf{K}}^{-1}\stackrel{··}{\mathbf{\tau }}+\stackrel{··}{\mathbf{q}}]+{\mathbf{N}}^{-1}\mathbf{\tau }+{\mathbf{f}}_0$$

(27a)

$$\mathbf{\tau }=\mathbf{M}(q)\stackrel{··}{q}+\mathbf{C}(q,\stackrel{·}{q})+\mathbf{G}(q)+{\mathbf{\tau }}_{\mathbf{L}}$$

(27b)

represents the nonlinear uncertainties.

High-order kinematics and dynamics

From the equation of above, it can be seen that the key problem of the inverse dynamics calculation of the elastic-joint robot is how to efficiently calculate the second derivative of the rigid robot driving torque. In this paper, the high-order kinematics and dynamics is derived based on Lie group and Lie algebra.

With the given Cartesian trajectory $\mathbf{X}$ of end-effector, the joint trajectory $\mathbf{q}$ can be obtained by inverse kinematic calculation, and the joint velocity $\stackrel{·}{\mathbf{q}}$ can be derived by the Jacobian matrix $\mathbf{J}(\mathbf{q})$ :

$$\stackrel{·}{\mathbf{q}=}{\mathbf{J}}^{-1}(\mathbf{q})\stackrel{·}{\mathbf{X}}$$

(28)

Differentiating the above equation, the joint acceleration $\stackrel{··}{\mathbf{q}}$ , jerk $\stackrel{···}{\mathbf{q}}$ , and snap $\stackrel{···}{\mathbf{q}}$ can be calculated.

The six dimensional generalized jerk and snap can be derived by differentiating the equation (18). Combined with the equations (17) and (18), through forward recursion (for i = 1, 2, …n), we can obtain:

$${\mathbf{V}}_i={\mathbf{G}}_{\mathbf{i}}^{-1}\frac{d({\mathbf{G}}_{\mathbf{i}})}{\mathrm{dt}}=A{d}_{{\mathbf{G}}_{\mathbf{i}-1,\mathbf{i}}^{-1}}({\mathbf{V}}_{i-1})+{\mathbf{\xi }}_{\mathbf{i}}{\stackrel{···}{\mathbf{q}}}_i$$

(29a)

$${\mathbf{A}}_{\mathbf{i}}=\frac{d({\mathbf{V}}_i)}{\mathrm{dt}}=A{d}_{{\mathbf{G}}_{\mathbf{i}-1,\mathbf{i}}^{-1}}({\mathbf{A}}_{\mathbf{i}-1})-a{d}_{{\mathbf{\xi }}_{\mathbf{i}}{\stackrel{·}{q}}_i}({\mathbf{V}}_i)+{\mathbf{\xi }}_{\mathbf{i}}{\stackrel{···}{\mathbf{q}}}_i$$

(29b)

$$\begin{array}{c}\multicolumn{1}{c}{\mathbf{J}{\mathbf{e}}_{\mathbf{i}}=\frac{d({\mathbf{A}}_{\mathbf{i}})}{\mathrm{dt}}=A{d}_{{\mathbf{G}}_{\mathbf{i}-1,\mathbf{i}}^{-1}}(\mathbf{J}{\mathbf{e}}_{\mathbf{i}-1})-a{d}_{{\mathbf{\xi }}_{\mathbf{i}}{\stackrel{·}{q}}_i}(A{d}_{{\mathbf{G}}_{\mathbf{i}-1,\mathbf{i}}^{-1}}({\mathbf{A}}_{\mathbf{i}-1}))}\\ \multicolumn{1}{c}{-a{d}_{{\mathbf{\xi }}_{\mathbf{i}}{\stackrel{··}{q}}_i}({\mathbf{V}}_i)-a{d}_{{\mathbf{\xi }}_{\mathbf{i}}{\stackrel{·}{q}}_i}({\mathbf{A}}_{\mathbf{i}})+{\mathbf{\xi }}_{\mathbf{i}}{\stackrel{···}{\mathbf{q}}}_i}\end{array}$$

(29c)

$$\begin{array}{c}{\mathbf{S}}_{\mathbf{i}}=\frac{d(\mathbf{J}{\mathbf{e}}_{\mathbf{i}})}{\mathrm{dt}}=A{d}_{{\mathbf{G}}_{\mathbf{i}-1,\mathbf{i}}^{-1}}({\mathbf{S}}_{\mathbf{i}-1})-2a{d}_{{\mathbf{\xi }}_{\mathbf{i}}{\stackrel{·}{q}}_i}(A{d}_{{\mathbf{G}}_{\mathbf{i}-1,\mathbf{i}}^{-1}}(\mathbf{J}{\mathbf{e}}_{\mathbf{i}-1}))\hfill \\ \multicolumn{1}{c}{-(a{d}_{{\mathbf{\xi }}_{\mathbf{i}}{\stackrel{··}{q}}_i}-{\mathrm{ad}}_{{\mathbf{\xi }}_{\mathbf{i}}{\stackrel{·}{q}}_i}^2)(A{d}_{{\mathbf{G}}_{\mathbf{i}-1,\mathbf{i}}^{-1}}({\mathbf{A}}_{\mathbf{i}-1}))-a{d}_{{\mathbf{\xi }}_{\mathbf{i}}{\stackrel{···}{\mathbf{q}}}_i}({\mathbf{V}}_{\mathbf{i}})}\\ -2a{d}_{{\mathbf{\xi }}_{\mathbf{i}}{\stackrel{··}{q}}_i}({\mathbf{A}}_{\mathbf{i}})-a{d}_{{\mathbf{\xi }}_{\mathbf{i}}{\stackrel{·}{q}}_i}(\mathbf{J}{\mathbf{e}}_{\mathbf{i}})+{\mathbf{\xi }}_{\mathbf{i}}{q}_i^{(4)}\hfill \end{array}$$

(29d)

Similarly, by differentiating the equation (23), the first derivative and second derivative of generalized forces $\stackrel{·}{\mathbf{F}}$ and $\stackrel{··}{\mathbf{F}}$ can be derived. Through backward recursion (for i = n, n–1,…1),

$${\mathbf{F}}_{\mathbf{i}}={\mathrm{Ad}}_{{\mathbf{G}}_{\mathbf{i},\mathbf{i}+1}^{-1}}^{*}({\mathbf{F}}_{\mathbf{i}+1})+{\mathbf{E}}_{\mathbf{i}}{\mathbf{A}}_{\mathbf{i}}-{\mathrm{ad}}_{{\mathbf{V}}_i}^{*}({\mathbf{E}}_{\mathbf{i}}{\mathbf{V}}_{\mathbf{i}})$$

(30a)

$$\begin{array}{c}\multicolumn{1}{c}{{\stackrel{·}{\mathbf{F}}}_{\mathbf{i}}={\mathrm{Ad}}_{{\mathbf{G}}_{\mathbf{i},\mathbf{i}+1}^{-1}}^{*}({\stackrel{·}{\mathbf{F}}}_{\mathbf{i}+1})-{\mathrm{Ad}}_{{\mathbf{G}}_{\mathbf{i},\mathbf{i}+1}^{-1}}^{*}({\mathrm{ad}}_{{\mathbf{\xi }}_{\mathbf{i}}{\stackrel{·}{q}}_i}^{*}({\mathbf{F}}_{\mathbf{i}+1}))+{\mathbf{E}}_{\mathbf{i}}\mathbf{J}{\mathbf{e}}_{\mathbf{i}}}\\ \multicolumn{1}{c}{-{\mathrm{ad}}_{{\mathbf{A}}_{\mathbf{i}}}^{*}({\mathbf{E}}_{\mathbf{i}}{\mathbf{V}}_{\mathbf{i}})-{\mathrm{ad}}_{{\mathbf{V}}_{\mathbf{i}}}^{*}({\mathbf{E}}_{\mathbf{i}}{\mathbf{A}}_{\mathbf{i}})}\end{array}$$

(30b)

$$\begin{array}{c}\multicolumn{1}{c}{{\stackrel{··}{\mathbf{F}}}_{\mathbf{i}}={\mathrm{Ad}}_{{\mathbf{G}}_{\mathbf{i}-1,\mathbf{i}}^{-1}}^{*}(({\mathrm{ad}}_{{\mathbf{\xi }}_{\mathbf{i}}{\stackrel{·}{q}}_i}^{*}{\mathrm{ad}}_{{\mathbf{\xi }}_{\mathbf{i}}{\stackrel{·}{q}}_i}^{*}-{\mathrm{ad}}_{{\mathbf{\xi }}_{\mathbf{i}}{\stackrel{··}{q}}_i}^{*})({\mathbf{F}}_{\mathbf{i}+1}))-2{\mathrm{Ad}}_{{\mathbf{G}}_{\mathbf{i}-1,\mathbf{i}}^{-1}}^{*}}\\ \multicolumn{1}{c}{({\mathrm{ad}}_{{\mathbf{\xi }}_{\mathbf{i}}{\stackrel{·}{q}}_i}^{*}({\stackrel{·}{\mathbf{F}}}_{\mathbf{i}+1}))+{\mathrm{Ad}}_{{\mathbf{G}}_{\mathbf{i},\mathbf{i}+1}^{-1}}^{*}({\stackrel{··}{\mathbf{F}}}_{\mathbf{i}+1})+{\mathbf{E}}_{\mathbf{i}}{\mathbf{S}}_{\mathbf{i}}-}\\ {\mathrm{ad}}_{{\mathbf{J}}_{\mathbf{i}}}^{*}({\mathbf{E}}_{\mathbf{i}}{\mathbf{V}}_{\mathbf{i}})-2{\mathrm{ad}}_{{\mathbf{A}}_{\mathbf{i}}}^{*}({\mathbf{E}}_{\mathbf{i}}{\mathbf{A}}_{\mathbf{i}})-{\mathrm{ad}}_{{\mathbf{V}}_{\mathbf{i}}}^{*}({\mathbf{E}}_{\mathbf{i}}\mathbf{J}{\mathbf{e}}_{\mathbf{i}})\hfill \end{array}$$

(30c)

Then the $\stackrel{·}{\mathbf{\tau }}$ and $\stackrel{··}{\mathbf{\tau }}$ can be obtained,

$${\stackrel{·}{\mathbf{\tau }}}_{\mathbf{i}}={\mathbf{\xi }}_i^T{\stackrel{·}{\mathbf{F}}}_i$$

(31a)

$${\stackrel{··}{\mathbf{\tau }}}_{\mathbf{i}}={\mathbf{\xi }}_i^T{\stackrel{··}{\mathbf{F}}}_i$$

(31b)

Combined with the equation (27), the elastic-joint industrial robot inverse dynamic calculation can be summarized as a flowchart as shown in Figure 2.

Figure 2.

Elastic-joint robot inverse dynamic calculation flowchart.

The inverse dynamics modeling method of elastic-joint robot based on Lie group and Lie algebra can simplify the modeling process and improve the computing efficiency compared with the traditional D-H modeling method. Therefor this modeling method has more advantages in inverse dynamics feedforward control. In order to analyze and prove this, the computational complexity of two different modeling methods is analyzed, and the simulation models of two different modeling methods are established to compare the simulation efficiency. For the sake of simplicity, the input of the two different models is given the same joint angular displacement of sinusoidal signal, and the simulation time is given as 5 s. The comparison of complexity analysis and calculation efficiency is shown in Table 1.

Table 1.

Comparison of computation complexity and efficiency of two different modeling methods.

Modeling method	The complexityof calculating ${\mathbf{\tau }}_m$	Calculating time of ${\mathbf{\tau }}_m$
Modeling method based on Lie groupand Lie algebra	$\begin{array}{c}\multicolumn{1}{c}{491n-15(\times )}\\ \multicolumn{1}{c}{390n-10(+)}\end{array}$	0.286 s
Traditional modelingmethod based onD-H parameters	$\begin{array}{c}\multicolumn{1}{c}{639n-192(\times )}\\ \multicolumn{1}{c}{553n-182(+)}\end{array}$	0.37 s

It can be seen from Table 1 that the modeling method based on Lie group and Lie algebra is more efficient to calculate the inverse dynamics of flexible joint robot. Compared with the traditional D-H parameter method, the efficiency is increased by 22.7%

Inverse dynamic control based on hybrid model

Principle of hybrid modeling

The principle of this hybrid modeling method is shown in Figure 3. It can be seen from the figure that the hybrid model is composed of analytical model, data-driven model, and the coupling relationship between them. If $f_k({\mathit{x}}_{\mathrm{in}},{\mathit{p}}_{\mathbf{k}})$ is the analytical model and $f_d({\mathit{x}}_{\mathrm{in}},{\mathit{p}}_d)$ is the data-driven model, the hybrid model can be expressed as follows:

$$\mathit{y}=f_k({\mathit{x}}_{\mathrm{in}},{\mathit{p}}_{\mathit{k}})\oplus f_d({\mathit{x}}_{\mathrm{in}},{\mathit{p}}_d)$$

(32)

Figure 3.

Schematic diagram of the hybrid modeling method.

In the formula, ${\mathit{p}}_{\mathit{k}}$ and ${\mathit{p}}_{\mathit{d}}$ are the parameter vectors of analytical model and data-driven model respectively, ${\mathit{x}}_{\mathrm{in}}$ and ${\mathit{y}}_{\mathrm{out}}$ are the input vector and output vector respectively. The symbol “⊕” represents a coupling operation between the analytical model and data-driven model. The parameters of analytical model are usually obtained by parameter identification, while the parameters of data-driven model are generally obtained by training, including on-line training and off-line training. The coupling relationship between analytical model and data-driven model is very complex, which usually needs to be reasonably designed according to specific problems. Series coupling and parallel coupling are two relatively simple and common coupling relations, which can be expressed by the following expressions respectively

$${\mathit{y}}_{\mathrm{out}}=f_k({\mathit{x}}_{\mathrm{in}},f_d({\mathit{x}}_{\mathrm{in}},{\mathit{p}}_d),{\mathit{p}}_{\mathit{k}})$$

(33)

$${\mathit{y}}_{\mathrm{out}}=f_k({\mathit{x}}_{\mathrm{in}},{\mathit{p}}_{\mathit{k}})+f_d({\mathit{x}}_{\mathrm{in}},{\mathit{p}}_d)$$

(34)

Design of the controller

As mentioned in the introduction, the industrial robot is a complex nonlinear system, the pure analytical inverse dynamic model can suffer from inaccuracy because of parameter uncertainties and non-parameter uncertainties. The inaccuracy inverse dynamic model is not conductive to high speed and high precision control. So in this paper, the hybrid analytical and data-driven inverse dynamic model is developed and the inverse dynamic control method based on the hybrid model is proposed. The analytical model is derived in the Sect.3, the data-driven mode is based on RBF neural network and the coupling relationship between the analytical model and the data-driven model is parallel coupling. The weights of RBF neural network are adjusted adaptively and updated online. The proposed control method based on the hybrid inverse dynamic model can be seen in Figure 4.

Figure 4.

Inverse dynamic control based on hybrid model.

The RBF neural network has the universal approximation property, ³⁵ so the data-driven model is established using RBF. If $\mathbf{x}$ is the input of the network, j is the $j^{\mathrm{th}}$ node of the network hidden layer, $\mathbf{h}=[{\mathrm{h}}_j{]}^{\mathrm{T}}$ is the Gauss function of the network and

$$h_j=\exp (\frac{||\mathbf{x}-{\mathbf{c}}_j|{|}^2}{2b_j^2})$$

(35)

then the output of RBF can be expressed as equation in which $\mathbf{W}$ is the weights of the network.

$$f(\mathbf{x})={\mathbf{W}}^{\mathrm{T}}\mathbf{h}(\mathbf{x})$$

(36)

In the Figure 4, the input of the RBF neural network is a vector composed of joint angular displacement and its derivatives, and joint angular control error and its derivatives.

Stablity analysis

The proposed control law can be expressed as follows,

$$\mathbf{u}={\mathbf{u}}_{\mathbf{ff}}+{\mathbf{u}}_{\mathbf{fb}}$$

(37)

where ${\mathbf{u}}_{\mathbf{ff}}$ is feedforward control and ${\mathbf{u}}_{\mathbf{fb}}$ is feedback control. ${\mathbf{u}}_{\mathbf{ff}}$ is computed by the hybrid inverse dynamic model,

$${\mathbf{u}}_{\mathbf{ff}}={\mathbf{J}}_{\mathbf{m}}\mathbf{N}{\mathbf{K}}^{-1}{\stackrel{··}{\mathbf{\tau }}}_{\mathbf{d}}+{\mathbf{J}}_{\mathbf{m}}\mathbf{N}{\stackrel{··}{\mathbf{q}}}_{\mathbf{d}}+{\mathbf{N}}^{-1}{\mathbf{\tau }}_{\mathbf{d}}+{\mathbf{J}}_{\mathbf{m}}{\mathbf{N}}^2{\mathbf{K}}^{-1}\mathbf{M}\hat{\mathbf{f}}$$

(38)

By applying equation (34) to equation (27a), the following equation can be obtained,

$${\mathbf{J}}_{\mathbf{m}}\mathbf{N}{\mathbf{K}}^{-1}\mathbf{\Delta }\stackrel{··}{\tau }+{\mathbf{J}}_{\mathbf{m}}\mathbf{N}\stackrel{··}{e}+{\mathbf{N}}^{-1}\mathbf{\Delta }\tau +{\mathbf{J}}_{\mathbf{m}}{\mathbf{N}}^{\mathbf{2}}{\mathbf{K}}^{-1}\mathbf{M}\widehat{f}-{\mathbf{f}}_{\mathbf{0}}+{\mathbf{u}}_{\mathbf{f}b}=0$$

(39)

where $\mathbf{e}={\mathbf{q}}_{\mathbf{d}}-\mathbf{q}$ and $\mathbf{\Delta }\mathbf{\tau }={\mathit{\tau }}_{\mathit{d}}-\mathbf{\tau }$ .

By equation (27b), the joint driving torque error is calculated as

$$\mathbf{\Delta }\mathbf{\tau }=\mathit{M}\stackrel{··}{\mathit{e}}+\mathbf{\Delta }{\mathit{f}}_1$$

(40)

Its second derivations is

$$\stackrel{··}{\mathbf{\Delta }}\mathbf{\tau }=\mathit{M}{\mathit{e}}^{(4)}+2\stackrel{·}{\mathit{M}}\stackrel{···}{\mathit{e}}+\stackrel{··}{\mathit{M}}\stackrel{··}{\mathit{e}}+\mathbf{\Delta }{\stackrel{··}{\mathit{f}}}_1$$

(41)

As the fourth derivative of the link error is much larger than the second or third derivatives, the terms, $2\stackrel{·}{\mathit{M}}\stackrel{···}{\mathit{e}}$ and $\stackrel{··}{\mathbf{M}}\stackrel{··}{\mathit{e}}$ are negligible compared with $\mathit{M}{\mathit{e}}^{(4)}$ for a typical motion. ²⁰ Therefore, the equation (38) is written as

$$\stackrel{··}{\mathbf{\Delta }}\mathbf{\tau }=\mathit{M}{\mathit{e}}^{(4)}+\mathbf{\Delta }{\mathit{f}}_2$$

(42)

By applying equations (37) and (39) to equation (36), the following equation can be obtained,

$$G{e}^{(4)}+\left(J_m{N}^2+M\right)\stackrel{··}{e}+\mathbf{G}\left(\widehat{f}-f\right)+\mathbf{N}{\mathbf{u}}_{fb}=0$$

(43)

where $\mathbf{G}={\mathit{J}}_{\mathit{m}}{\mathit{N}}^2{\mathit{K}}^{-1}\mathit{M}$ and $\mathit{f}={\mathit{G}}^{-1}\mathbf{N}{\mathit{f}}_0-\mathbf{\Delta }{\mathit{f}}_2-{\mathit{G}}^{-1}\mathbf{\Delta }{\mathit{f}}_1$

The feedback control law can be referred to Moberg and Hanssen, ²²

$$\mathbf{N}{\mathbf{u}}_{fb}=\mathbf{G}\left(K_3\stackrel{⃛}{e}+K_2\stackrel{··}{e}+K_1\stackrel{·}{e}+K_0\mathbf{e}\right)-\left(J_m{N}^2+M\right)\stackrel{··}{e}$$

(44)

where ${\mathit{K}}_{\mathit{i}}=\mathrm{diag}({\mathrm{K}}_{i1},{\mathrm{K}}_{i2},\dots {\mathrm{K}}_{\mathrm{in}})$ is a constant diagonal matrix. Then the error dynamics can be obtained:

$${\mathit{e}}^{(4)}+{\mathit{K}}_3\stackrel{···}{\mathit{e}}+{\mathit{K}}_2\stackrel{··}{\mathit{e}}+{\mathit{K}}_1\stackrel{·}{\mathit{e}}+{\mathit{K}}_0\mathbf{e}={\mathit{U}}_{\mathit{f}}=\mathit{f}-\hat{\mathit{f}}$$

(45)

If the $\mathbf{W}$ is the ideal weight of RBF neural network and $\hat{\mathbf{W}}$ is the estimation weight of the neural network, then

$${\mathit{U}}_{\mathit{f}}={\mathit{W}}^{\mathit{T}}\mathit{h}\left(\mathit{x}\right)-{\hat{\mathit{W}}}^{\mathit{T}}\mathit{h}\left(\mathit{x}\right)={\tilde{\mathit{W}}}^{\mathit{T}}\mathit{h}(\mathit{x})$$

(46)

The servo error $\mathit{E}$ is defined,

$$\mathit{E}=\stackrel{··}{\mathit{e}}+{\mathbf{\Psi }}_2\stackrel{··}{\mathit{e}}+{\mathbf{\Psi }}_1\stackrel{·}{\mathit{e}}+{\mathbf{\Psi }}_0\mathit{e}$$

(47)

In equation (44), ${\mathbf{\Psi }}_2=\mathrm{diag}({\mathrm{\Psi }}_{21},{\mathrm{\Psi }}_{22},\dots {\mathrm{\Psi }}_{2n})$ with ${\mathrm{\Psi }}_{2i}>0$ , ${\mathbf{\Psi }}_1=\mathrm{diag}({\mathrm{\Psi }}_{11},{\mathrm{\Psi }}_{12},\dots {\mathrm{\Psi }}_{1n})$ with ${\mathrm{\Psi }}_{1i}>0$ , and ${\mathbf{\Psi }}_0=\mathrm{diag}({\mathrm{\Psi }}_{01},{\mathrm{\Psi }}_{02},\dots {\mathrm{\Psi }}_{0n})$ with ${\mathrm{\Psi }}_{0i}>0$ . The value $\mathit{E}$ can be computed from sensor readings. Then

$$\frac{E_j}{U_{\mathrm{fj}}}=\frac{s^3+{\mathrm{\Psi }}_{2j}{s}^2+{\mathrm{\Psi }}_{1j}s+{\mathrm{\Psi }}_{0j}}{s^4+K_{3j}{s}^3+K_{2j}{s}^2+K_{1j}s+K_0}$$

(48)

The ${\mathrm{\Psi }}_{2j}$ , ${\mathrm{\Psi }}_{1j}$ , and ${\mathrm{\Psi }}_{0j}$ are chosen so that the above transfer function is strictly positive real (SPR). Then by the positive real lemma, we are assured of the existence of the positive definite matrices ${\mathbf{P}}_{\mathit{j}}\in {\mathbb{R}}^{4\times 4}$ and − ${\mathbf{Q}}_{\mathit{j}}\in {\mathbb{R}}^{4\times 4}$ , such that

$${\mathit{A}}_{\mathit{j}}^{\mathit{T}}{\mathbf{P}}_{\mathit{j}}+{\mathbf{P}}_{\mathit{j}}\mathit{A}=-{\mathbf{Q}}_{\mathit{j}}$$

(49a)

$${\mathbf{P}}_{\mathit{j}}{\mathbf{B}}_{\mathit{j}}={\mathit{C}}_{\mathit{j}}^{\mathit{T}}$$

(49b)

where the matrices ${\mathbf{A}}_{\mathit{j}}$ , ${\mathbf{B}}_{\mathit{j}}$ , and ${\mathbf{C}}_{\mathit{j}}$ are the matrices of a minimal state-space realization of the transfer function (45).

$$\stackrel{.}{{\mathit{x}}_{\mathit{j}}}={\mathbf{A}}_{\mathit{j}}{\mathit{x}}_{\mathit{j}}+{\mathbf{B}}_{\mathit{j}}{U}_{\mathrm{fj}}$$

(50a)

$${\mathrm{E}}_j={\mathit{C}}_{\mathit{j}}^{\mathit{T}}{\mathit{x}}_{\mathit{j}}$$

(50b)

The state vector ${\mathit{x}}_{\mathit{j}}$ is a linear transformation of $[{\mathbf{e}}_{\mathbf{j}},{\stackrel{·}{\mathbf{e}}}_{\mathbf{j}},{\stackrel{··}{\mathbf{e}}}_{\mathbf{j}},{\stackrel{···}{\mathbf{e}}}_{\mathbf{j}}{]}^T$ .

Then the error dynamic equation of the entire robot system in state-space form is given by

$$\stackrel{·}{\mathit{X}}=\mathit{A}\mathbf{X}+\mathit{B}{\mathit{U}}_{\mathit{f}}$$

(51a)

$$\mathit{E}=\mathit{CX}$$

(51b)

where $\mathit{A}$ , $\mathit{B}$ , and $\mathbf{C}$ are all block diagonal with ${\mathbf{A}}_{\mathit{j}}$ , ${\mathbf{B}}_{\mathit{j}}$ , and ${\mathbf{C}}_{\mathit{j}}$ on the diagonals, respectively. $\mathit{X}=[{\mathit{x}}_1,{\mathit{x}}_2,\dots {\mathit{x}}_{\mathit{n}}{]}^{\mathit{T}}$ , $\mathbf{P}=\mathrm{diag}({\mathit{P}}_1,{\mathit{P}}_2,\dots {\mathit{P}}_{\mathit{n}})$ , and $\mathbf{Q}=\mathrm{diag}({\mathit{Q}}_1,{\mathit{Q}}_2,\dots {\mathit{Q}}_{\mathit{n}})$ .It is obvious that

$${\mathit{A}}^{\mathit{T}}\mathit{P}+\mathit{PA}=-\mathit{Q}$$

(52a)

$$\mathit{PB}={\mathit{C}}^{\mathit{T}}$$

(52b)

Now use Lyapunov theory to analyze the stability. The following Lyapunov function is defined

$$\mathit{v}(\mathit{X},\tilde{\mathit{W}})={\mathit{X}}^{\mathit{T}}\mathit{PX}+\mathrm{tr}({\tilde{\mathit{W}}}^T{\mathbf{\Gamma }}^{-1}\tilde{\mathit{W}})$$

(53)

where $\mathbf{\Gamma }=\mathrm{diag}({\mathrm{\lambda }}_1,{\mathrm{\lambda }}_2,\dots {\mathrm{\lambda }}_r)$ and ${\mathrm{\lambda }}_i>0$ . Differentiate the Lyapunov function and make use of equation (49), we obtain

$$\stackrel{·}{\mathit{v}}(\mathit{X},\tilde{\mathit{W}})=-{\mathit{X}}^{\mathit{T}}\mathit{QX}+2\mathrm{tr}({\tilde{\mathit{W}}}^T(\mathbf{h}\left(\mathbf{x}\right){\mathit{E}}^{\mathit{T}}+{\mathbf{\Gamma }}^{-1}\stackrel{.}{\stackrel{~}{\mathit{W}}}))$$

(54)

Therefore, if we chose

$$\stackrel{.}{\hat{\mathit{W}}}=\mathbf{\Gamma }\mathbf{h}(\mathbf{x}){\mathit{E}}^{\mathit{T}}$$

(55)

we have

$$\stackrel{·}{\mathit{v}}(\mathit{X},\tilde{\mathit{W}})=-{\mathit{X}}^{\mathit{T}}\mathit{QX}$$

(56)

which is non-positive because $\mathit{Q}$ is positive definite. Then the control system is proved to be stable. That is to say, if the weight of RBF neural network satisfies the equation (52), the proposed control law can guarantee the stability of the closed-loop system.

Simulation and verification

Multi-domain integrated robot system model

In order to validate the performance of the proposed control method, a multi–domain industrial robot system model is established. Different from traditional model, the model established in this paper integrated mechanical, control and electrical subsystems. It is much closer to the real robot system. The system model is mainly composed of mechanical part, electrical subsystem, control subsystem and trajectory planning algorithm. The model topology can be seen in the Figure 5.

Figure 5.

Topology of the multi-domain integrated robot system model.

The robot is the typical 6-DOF industrial robot and the main parameters are show in the following Table 2.

Table 2.

Main parameters of industrial robot system.

Parameter	Value
Mass of links/(kg)	M = [8.949, 13, 9.174, 4.413, 3.089, 0.6635]
Transmission ratio of each joint	N = [81, 100, 80, 76.5, 81, 50]
Moment of inertia of each motor rotor/( $\mathrm{Kg}*c{m}^2$ )	${\mathit{J}}_{\mathit{m}}=$  [1.3, 1.3, 0.44, 0.27, 0.27, 0.052]
Torsional stiffness of each joint/( $\mathrm{N}.\mathrm{m}/\mathrm{rad}$ )	K = [43, 20, 58, 500, 800, 1000]

The initial orientation ( $\mathbf{\theta }=0$ ) of the robot mechanical part, main geometric size and the link frame configuration can be clearly seen in the Figure 6.

Figure 6.

Mainly geometric size and frame configuration.

Simulation and analysis

In order to analyze the proposed control algorithm, two types of typical trajectories are given for simulation: one is the line trajectory and another is circle trajectory. Fifteen segments S-type velocity planning is used for motion planning. The trajectory and motion planning can be seen in the Figure 7.

Figure 7.

Trajectory and motion planning.

Because of the accurate robot dynamic model is hardly to be established, the control performance of the traditional inverse dynamic feedforward control method is limited. In order to demonstrated that the control method based on hybrid inverse dynamic model proposed in this paper can avoid the decrease of control accuracy caused by the inaccuracy of the model to some extent, three different simulations are performed. In the first case, the inverse dynamic model is completely accurate, that is the inverse dynamic control is based on accurate model (IDC-AM). In the second case, it is assumed that the estimated dynamic parameters (joint stiffness, link mass, and link inertia tensor) are 90% of the real value. That is to say the inverse dynamic control is based on inaccurate model (IDC-IM). In the three case, the inaccurate model is combined with the RBF neural network model, the inverse dynamic control is based on the hybrid model (IDC-HM). The simulation results in these three cases are compared.

This research is mainly focused on trajectory tracking precision. When the straight line command trajectory is given, the simulation results can be seen in Figure 8. For conveniently reading the figure, the straight line contour error is magnified by 40 times.

Figure 8.

Simulation results of straight line trajectory tracking performance.

It can be seen directly from the Figure 8 that trajectory tracking precision has decreased a lot when the inverse dynamic model is inaccurate. When the inaccurate model is combined with the RBF neural network model, the trajectory tracking precision has improved. Similarly, when the circle command trajectory is given, the same conclusion can be found, as can be seen in the Figure 9 (The contour error is magnified by 40 times). The compared simulation results shows that the inverse dynamic control based on hybrid model can improve the control performance.

Figure 9.

Simulation results of circle trajectory tracking performance.

The curve of contour error respect to time in the three cases can be seen in Figures 10 and 11. The control performance can be compared more clearly from the figures. In order to evaluate the contour error quantitatively, the maximum absolute error (MAE) and root mean square error (RSME) are evaluated. The details can be seen in Table 3. It can be seen in the table that the MAE and RSME of IDC-HM contour error are both smaller than that of IDC-IM whether the command trajectory is straight line or circle. The simulation results shows that the elastic-joint industrial robot inverse dynamic control based on hybrid model proposed in this paper can overcome the bad influence on the control performance caused by the inaccuracy of the dynamic model.

Figure 10.

Comparison of the straight line contour error.

Figure 11.

Comparison of the circle contour error.

Table 3.

Comparison of contour error of different control method.

Type of trajectory	Control method	MAE	RMSE
Straight linetrajectory	IDC-AM	0.18 mm	0.14 mm
	IDC-IM	1.15 mm	0.81 mm
	IDC-HM	0.54 mm	0.34 mm
Circletrajectory	IDC-AM	0.22 mm	0.14 mm
	IDC-IM	0.79 mm	0.51 mm
	IDC-HM	0.57 mm	0.20 mm

The influence of load on the robot control performance is also analyzed. Assume that under the condition of circular trajectory planning in Figure 7, a 200 N resistance is applied at the end of the robot along the tangent direction of the arc trajectory at the time of 0.6 s, and the direction of the force is always opposite to the direction of the speed, as shown in Figure 12. Using the traditional inverse dynamics feedforward control method, this sudden increase of load is difficult to be considered in the inverse dynamics model. Therefore, the uncertainty of robot end-effector load will inevitably affect the performance of traditional inverse dynamics feedforward control method. If the feedforward control method based on hybrid inverse dynamic model proposed in this paper is adopted, the influence of load uncertainty on control accuracy will be weakened. In the case of sudden change load applied at the end-effector of the robot, the control performance of different control schemes is compared through the simulation analysis based on the multi domain system synthesis model of the robot, and the superiority of the control method proposed in this paper is verified.

Figure 12.

The given load force curve of robot end-effector.

As shown in the Figure 13, the simulation results of trajectory tracking performance under two different control schemes are compared, that is, the traditional inverse dynamics feedforward control with inaccurate model and the proposed feedforward control based on hybrid inverse dynamic model. It is obvious from the figure that the trajectory tracking accuracy of the proposed feedforward control method based on hybrid inverse dynamic model is higher.

Figure 13.

Position tracking error comparison of different control schemes under load disturbance.

Industrial robot is a nonlinear complex system, it is difficult to get accurate dynamic model based on mechanism modeling, so the inverse dynamic control based on pure analytical model is hard to achieve high control performance. The simulation results verified this conclusion and from the simulation results, if can be seen that if the analytical mode is combined with the data-driven model, the inverse dynamic control based on the hybrid model can improve the trajectory tracking performance. The advantage of the proposed control method is demonstrated by the simulation results.

Conclusions

This paper proposed an inverse dynamic feedforward control method of industrial robot with elastic joints based on hybrid model. The hybrid model is composed of analytical inverse dynamic model and data-driven model. Considering the elastic joint of the industrial robot, the inverse dynamic analysis is extremely complex, so in this paper the analytical model is derived based on the theory of Lie group and Lie algebra. Because of the complexity of the real robot system, the pure analytical model can hardly achieve high accuracy. The data-driven model is used to compensate the parameter uncertainties and non-parameter uncertainties. In order to validate the propose control method, the multi-domain integrated robot system model is established and by simulation, the control performance is compared with the traditional inverse dynamic control. The simulation results demonstrated the advantage of the proposed control method. The control method proposed in this paper can overcome the influence of the uncertainty of the inverse dynamic model on the control performance to a certain extent, improve the control performance of the flexible joint robot, and has important application value for realizing the high-precision control of the robot.