1. Introduction
The inverse dynamics model of industrial robots plays a key role in improving robot control accuracy and operation efficiency. These inverse dynamics models allow robots to perform complex tasks more accurately and effectively by predicting and calculating the forces and torques required. In high-performance manufacturing, automated assembly, and fine operation, the application of inverse dynamics models is particularly important. However, these inverse dynamics models still face challenges when dealing with complex variables and uncertainties in the actual industrial environment. For example, environmental noise, joint friction, dynamic load change, and uncertain factors such as flexibility and joint clearance of the robotic arm may lead to inaccurate inverse dynamics model predictions, affecting industrial robot performance. Considering these problems and urgent needs, academia and industry are actively researching improving the accuracy and computational efficiency of inverse dynamics models to better meet the needs of modern industry for high-precision and high-efficiency industrial robots.
The uncertainty residuals of the inverse dynamics model originated from influencing factors such as the greater influence of the motor by temperature, the flexible connection of the harmonic reducer, and the greater influence of temperature on nonlinear friction. Therefore, Gao et al. [Reference Gao, Yuan, Han, Wang and Wang1] proposed that an accurate and reasonable friction model is important for studying the inverse dynamics mode. Iskandar et al. [Reference Iskandar and Wolf2] proposed a collaborative joint friction model of industrial robots that comprehensively considered the effects of speed, temperature, and load torque. Callar [Reference Çallar and Böttger3] proposed a physically excited friction model with a parametric description of a nonlinear dependence of temperature and velocity as well as the dependence on external loads.
The friction model alone is insufficient; the model still has many gaps, which cause uncertain residuals. Fitting models such as neural networks can provide efficient data-driven optimization with a small deviation. Romeres et al. [Reference Romeres, Zorzi, Camoriano, Traversaro and Chiuso4] developed a learning method for an inverse dynamics model of a long short-term memory (LSTM) network based on time complexity. Seeger et al. [Reference Seeger5] discussed an online robot inverse dynamics modelling algorithm. Dalla et al. [Reference Libera, Giacomuzzo, Carli, Nikovski and Romeres6] used the Gaussian process (GP) framework to model the inverse dynamics model. Michael et al. [Reference Lutter, Listmann and Peters7] proposed a method for forward dynamics estimation from data-driven inverse dynamics learning, which derived the physical dynamics equations of rigid body dynamics (RBD) by learning the inverse dynamics model and estimating each dynamics component from it. However, these purely data-driven inverse kinetics models usually have low data efficiency, cannot guarantee the stability of the extrapolation method, and can achieve good performance only near the training data, which prevents them from being used in safety-critical applications in semistructured environments.
Hybrid models inspired by physics-inspired neural networks, such as deep Lagrangian networks (DeLaNs) [Reference Greydanus, Dzamba and Yosinski8] or Hamiltonian networks [Reference Gupta, Menda, Manchester and Kochenderfer9], ensure physically reasonable dynamic models and thus save energy. These physics-inspired neural networks incorporate the system dynamics structure into the neural network architecture, are limited to modelling conservative forces, and capture all nonconservative effects using an additional neural network. The residual hybrid model combines the rigid body model with an additional residual model to learn the residuals of the RBD model [Reference Meier, Kappler, Ratliff and Schaal10]. A GP or a neural network (NN) can be used for the residual model. Nguyen-Tuong and Eters [Reference Nguyen-Tuong and Peters11] used the GP model to learn the residuals of the inverse dynamics model. Several recent works have combined standard multilayer perceptron NNs with rigid body models [Reference Wang and Liu12–Reference Sun, Shao, Qu, Guan and Tan15]. Similar methods have been applied to physics-inspired NNs to capture nonconservative forces. Our previous work successfully applied different recurrent NNs to the black-box inverse dynamics model [Reference Shaj, Becker, Buchler, Pandya, Niels van Duijkeren, Hanheide and Neumann16]. The LSTM algorithm was applied to train the residual mixture model, and the results were better than those of GP.
Industrial robot is a complex nonlinear system with multiple inputs and multiple outputs, with time-varying, strongly coupled, nonlinear, and other complex dynamics, thus making the high-performance control of the robot very complex and difficult [Reference Spong17]. There are many control methods for robotic arms nowadays; for example, An et al. proposed a sparrow search algorithm [Reference An, Chen, Sun, Yin, Song, Zhuang, Chang and Liu18] based on the Levy flight operator for self-correcting robot controller parameters to improve the accuracy of motion trajectory. Chen et al. proposed a target tracking control system based on a bionic closed-loop central pattern generator to achieve accurate target tracking [Reference Chen, Xu, Yang, Hu, Cheng, Zhu, Zhang, Shi and Chai19]. Wu et al. formulated a systematic approach for the elastic dynamics analysis of tandem robot manipulators [Reference Wu20]. Chen et al. proposed a segmented dynamic modelling approach to construct propulsion and lift models for the beaver-like tail of a robot. Combining fluid dynamics and material mechanics, a theory of sectorial flexible dynamics was developed [Reference Chen, Xu, Wang, Tu, Hu, Chen, Xu, Chai, Zhang and Shi21]. Model-based predictive control (MPC) algorithms [Reference Tatjewski22] provide better results with less accurate models avoiding the need for dynamic modelling of disturbances. However, traditional Proportion Integration Differentiation (PID) control as well as MPC control can hardly meet the demand of high speed and high-precision control in modern industry, and feedforward control based on the inverse dynamics model is an effective way to improve the control performance of robots. However, this control method has the problem of model imprecision. The actual industrial robot system has parameter uncertainty, non-parametric uncertainty, external environment interference, etc. It is difficult to establish an accurate inverse dynamics model of the robot, and the inaccuracy of the model will lead to a decrease in the control performance. In this paper, under the assumption that the robot is an ideal rigid body and elastic deformation is not taken into account, we carry out a research to address the above problems and propose a hybrid inverse dynamics model compensation method for industrial robots oriented to the feedforward control method of robots. However, only the LSTM algorithm is still insufficient for the compensation accuracy of the inverse dynamics of industrial robots, and there is an urgent need to propose new algorithms to further improve the compensation accuracy. Therefore, this paper proposes a combination of bootstrap aggregating (bagging) algorithm and LSTM network and at the same time introduces the linear layer as the network optimization layer. The BLL residual prediction algorithm (BLL algorithm for short) is constructed to improve the accuracy of the residual model of industrial robots.
First, the rigid body inverse dynamics model and the hybrid inverse dynamics model based on the compensation of the residual prediction algorithm are introduced. Then, the composition and principle of the BLL algorithm are introduced. Finally, experiments are conducted based on the public dataset of the Franka Emika Panda robot to validate the average residual value of each joint of the hybrid inverse dynamics model of the robot decreased from 0.5651 N · m to 0.1096 N · m, which improved the calculation accuracy of the robot inverse dynamics model. Compared with the LSTM and GP algorithms in this field, the proposed BLL algorithm improves the mean squared error (MSE), root mean square residual error (RMSE), mean absolute error (MAE), R2, and the reduction rate of the mean value of the residuals after compensation, validating the effectiveness of the proposed algorithm. Figure 1 shows the flow chart of the model compensation technique.
2. Robot inverse dynamics model
2.1. Body inverse dynamics model
The primary purpose of industrial robot dynamics is to achieve real-time control and an accurate dynamic model; accurate dynamic parameters are the key. Robot inverse dynamics modelling methods include Lagrangian methods, recursive Newton–Euler methods, symbolic computational methods, and methods based on a mixture of symbolic and numerical methods. Among these, the Newton–Euler [Reference Zhang, Xu, Yao and Zhao23] method has a relatively fast computational speed, especially when using a recursive form, which is well suited for robot control systems and is computationally efficient [Reference Yilmaz, Wu, Kazanzides and Tumerdem24], so the Newton–Euler method is used here, with the following expression:
$\mathrm{q},\dot{\mathrm{q}},\ddot{\mathrm{q}}$ are the vectors of the joint position, velocity, and acceleration, respectively, and $n$ represents robot degrees of freedom (DOF). ${\unicode[Arial]{x03C4}}$ is the joint torque vector, $\mathrm{M}(\mathrm{q})\in \mathrm{R}^{\mathrm{n}\times \mathrm{n}}$ is the inertia matrix, $\mathrm{H}(\mathrm{q},\dot{\mathrm{q}})\in \mathrm{R}^{\mathrm{n}}$ is the vector of the centrifugal force and the Coriolis force, $\mathrm{G}(\mathrm{q})\in \mathrm{R}^{\mathrm{n}}$ is the gravitational torque or force, and $\xi \in R^{n}$ is the offset torque, which represents the dynamic uncertainty influencing factor. The nonlinear friction factor is one of the factors affecting the high-precision movement of a robot [Reference Wang and Wang25]. The dynamic identification model of industrial robots usually uses the viscosity and Coulomb friction model to represent the effect of friction on robot dynamics:
where $f_{v},f_{c}$ are the viscosity and the Coulomb friction coefficient, respectively, and ${\unicode[Arial]{x03F5}} (\mathrm{q},\dot{\mathrm{q}},\ddot{\mathrm{q}})\in \mathrm{R}^{\mathrm{n}}$ represents the factors not modelled in the physical kinetics model. Predicting the unmodelled part is the focus of this paper in Sections 3 and 4 and is described in detail.
2.2. Hybrid inverse dynamics model compensated based on the BLL residual prediction algorithm
Based on Formula (2.a), this paper proposes a hybrid architecture that couples the RBD model, friction model, and residual error prediction algorithm. The BLL algorithm models the residual errors of the RBD model and accounts for the uncaptured effects, such as the flexibility of the robot. The architecture proposed in this paper is composed of the following three parts:
$f_{\mathrm{RBD}}$ is a partial rigid body model using the motion equations in Eq. (2.a). $f_{\mathrm{F}} f_{BLL}$ is a residual prediction model that captures some observable effects, such as joint and link flexibility and more complex friction effects.
This paper proposes a hybrid inverse dynamics model framework, as shown in Fig. 2. First, the Newton–Euler algorithm is used for inverse dynamics joint torque calculations. Then, the values of the friction model are calculated. Finally, the calculated joint torques are used as the BLL network inputs to train the BLL network model. The network-predicted values are compared with those of the RBD model. The friction models are combined to obtain a hybrid inverse dynamics model. In this paper, the main objective is to predict the residual error.
3. Principle of the BLL residual prediction algorithm
The BLL residual prediction algorithm proposed in this paper comprises three layers: a sample data sampling layer based on ensemble learning bagging, a network LSTM prediction layer, and a linear optimization fitting layer to predict the residual joint torque values (Fig. 3).
Bootstrap aggregating (bagging) was first proposed by Leo Breiman in 1996 [Reference Breiman26]. The bagging algorithm can be combined with other classification and regression algorithms to improve the accuracy and stability and reduce the variance in the results simultaneously to avoid overfitting.
LSTM NNs are a special type of recurrent neural network (RNN) [Reference Graves27]. During the training of the original RNN, with prolonged training time and an increase in the number of network layers, the problem of gradient explosion or gradient disappearance is prone to occur, and the RNN cannot process longer sequence data and thus cannot obtain long-distance data. LSTM was proposed to solve this problem.
The linear layer, also called the fully connected layer, is a basic network layer structure in deep learning [Reference Abdi28]. A linear layer usually contains an input matrix, a weight matrix, and a bias vector.
In this paper, the bagging algorithm is first used to sample n groups, and an independent LSTM network is used for training. At this time, the learning rate of the fully connected layer is set to 0, and no training is performed. The trained n groups of LSTM network models are validated using the same validation set, and n sets of validation results are obtained. These results still have a residual error $\epsilon$ from the actual value. To further optimize the residuals, the n-group LSTM validation set results are averaged and combined with the true values $y_{m}$ of the original validation set to form a new training set, which is used to train the subsequent fully connected layer to optimize the residuals. The learning rate of the fully connected layer is changed to 0.01. Then, a unified whole is validated with the test dataset after the training is complete. The validation is performed using the same validation set to validate the n-group model and perform residual optimization.
Bagging generates different training subsets based on random sampling, such that some data in each subset are not sampled. In this way, more data changes can be introduced, which can increase the difference between models and improve the diversity of the ensemble model. In addition, as a deep learning network, LSTM has many parameters that require training and is prone to overfitting the training data. Through bagging, since each learner is trained on a different subset, overfitting can be reduced, and the generalization ability of the model can be improved. Finally, the linear layer is trained and optimized based on the LSTM prediction result to reduce the residual error further.
In this paper, three algorithms, namely, bagging, LSTM, and linear, are coupled, and the mathematical mapping formula corresponding to the BLL algorithm is created according to the coupling mechanism. The relevant mathematical symbols are shown in Table I.
Bagging-based random sampling:
Next, the LSTM algorithm makes independent predictions. Figure 4 shows the principle of LSTM.
The underlying logic of LSTM is to selectively retain and forget historical information [Reference Liu, Li, Hao, Yang, Hu, Xue and Wang30] based on various gating weight matrices. The specific formulas are as follows:
Forget gate:
Input gate:
Update vector:
Cell status update:
Output gate:
Cell state output:
Output value:
After n groups of LSTM networks are independently trained, the average value of the n models is taken at each moment after validation is performed using test_data:
Assuming there are K time points in the validation set, the average vector predicted by LSTM is:
The validation results of the test_data and actual result are used as the input of the linear layer to train a NN with a total of three layers, including input, output, and linear hidden layers. The linear weight matrix is:
The linear hidden layer is:
The output result of the full connection is:
4. Residual model validation
4.1. Dataset preprocessing
The dynamic parameters and public datasets of the Franka Emika Panda robot used in this paper are from the “Dynamic Identification of the Franka Emika Panda Robot With Retrieval of Feasible Parameters Using Penalty-Based Optimization” provided by Claudio Gaz et al. [Reference Gaz, Cognetti, Oliva, Robuffo Giordano and De Luca31].
In this paper, the angular acceleration is calculated by derivation based on the time, angular displacement, and angular velocity of the seven joints in the dataset. The inverse Newton–Euler algorithm [Reference Lutter32] was used to calculate the torque variation in each joint of the rigid body inverse dynamics model (2.a).
In this paper, the data collected at 12.9 S–13.2 S are selected as the results. The input layer of the inverse Newton–Euler algorithm has seven channels: the joint angle, joint velocity, and joint acceleration of the seven DOF of the Franka robot; the output is seven channels, which are the calculated torques of the seven joints of the Franka robot. The calculation results of the inverse Newton–Euler algorithm are compared with the actual acquisition results, as shown in Fig. 5.
The residual data of subtracting the calculated and actual torque are shown in Fig. 6.
The above figure shows a large gap exists between the calculated and actual results. The model training process and results are described in detail in Section 4.2.
4.2. Analysis of the model training set results
The hardware conditions of this model were an Intel i7-13700H CPU, 16 GB of memory, an NVIDIA RTX4060 graphics card with 8 GB of video memory, an operating system of Windows 11, and environmental configurations of Anaconda Python 3.9, PyTorch-Cuda 11.8, and PyTorch 2.1.0.
The input layer of the BLL algorithm is the torque residual of a single joint, which is the difference between the calculated torque and the actual torque of the RBD of the robot; the output layer of the BLL is a single channel, which is the predicted value of the torque residual of a single joint. The figure below shows the decline curve of the loss function when the BLL algorithm trains seven joint data points alone.
The fitness value function decreases as the number of iterations increases. After 800 iterative training tests, most converge within 130 iterations. Therefore, in this paper, the number of epochs is adjusted to 130. The training loss is shown in Fig. 7. The training data are presented in Table II.
Table II shows that the training results of the decrease rate of the fitness function of the model training were between 94.4% and 96.4%. Additionally, the model was able to fit the data points better. Next, the BLL parameter settings for training are introduced (Table III).
In Table III, the dimensions of the BLL hidden layer include the original LSTM and fully connected layers. Therefore, the hidden dimension needs to be divided into LSTM_hidden_dim and linear_hidden_dim1-linear_hidden_dim3. The effect of a batch size of approximately 32 is the best. The following is the evaluation of the test set and the model comparison.
4.3. Results analysis on the model test set
To verify the effectiveness and improvement effect of the proposed BLL network, 290 verification points were arbitrarily selected within the range of 25,000 datasets to verify the training set. The comparison results of the actual residuals with the BLL-predicted residuals and the comparison between the actual residuals and model-compensated residuals are shown in Figs. 9 and 10, respectively.
The BLL network proposed in this paper uses four indicators: the MSE, the RMSE, the MAE, and the R2 score. We conducted a test set analysis [Reference Spong, Hutchinson and Vidyasagar33]. The definitions and descriptions of these four indicators are given as follows:
where $y_{i}$ represents the true value of the dataset, $\hat{y}_{i}$ represents the predicted value of the dataset, $\overline{y}_{i}$ represents the mean value of the predicted value of the dataset, and n represents the number of datasets and the variance in the data. R2 represents the determination coefficient of the model. The best score is 1.0, indicating that the model predicts the true value perfectly. It can also be negative because the model can vary arbitrarily; that is, there is no mapping fit between the predicted and real data. Table IV shows the MSE, RMSE, MAR, and R2 values.
The determination coefficient of R2 of the predicted values of the specific residuals is shown in Fig. 8. R2 is approximately 0.84–0.98. The predicted values (blue dots) are closely distributed around the true value (red line), with R2 close to 1, indicating that the predicted values have a high correlation with the actual value and a high fitting accuracy.
During the BLL network training, the training and iterations were performed in a single dimension, so the characteristics of each dimension were inconsistent. In addition, the accuracy of the robot end was affected by nonlinear factors such as the accuracy of the dataset and the environmental conditions, which caused some differences in the effect of residual compensation of the end-accuracy of the robots, such as joint 6.
Table IV shows that the MSE, RMSE, and MAE of the predicted values of the joint robot residuals are all close to zero. Therefore, the proposed machine learning model is adaptable to the prediction of the residuals of the inverse dynamics model.
As Figs. 9 and 10 show, before compensation, the absolute values of the Joint1 residuals were approximately evenly distributed at [4.0374, 0.0006]; the absolute values of the Joint2 residuals were distributed at [5.4627, 0.0121]; the absolute values of the Joint3 residuals were distributed at [4.3725, 0.0066]; the absolute values of the Joint4 residuals were distributed at [2.9014, 0.0109]; the absolute values of the Joint5 residuals were distributed at [0.7850, 0.0039]; the absolute values of the Joint6 residuals were distributed at [0.6646, 0.0014]; and the absolute values of the Joint7 residuals were evenly distributed at [0.4305, 0.0088]. Using the BLL algorithm-based residual compensation method proposed in this paper, the residual errors of the seven robot joints are distributed at approximately 0, the fluctuations near $\pm 0.2N\cdot m$ , and the fluctuation range is very small, indicating that the compensated accuracy has high stability and can improve the operating accuracy of the robot.
The results of the static statistical analysis of the robot inverse dynamics model before and after residual compensation are shown in Table V. Compared with the precompensation, the model residual of the seven robot joints was reduced [62.76–92.72%] to good effect.
The robot dataset is based on Franka Emika Panda, a highly collaborative robot. The torque residual range is much smaller than that of traditional heavy-duty industrial robots. Therefore, using LSTM or GP [Reference Rezaei-Shoshtari, Meger and Sharf34] for feature extraction, model training, and optimization is more difficult. In this paper, the BLL algorithm is proposed to predict and compensate for the residuals of the Franka Emika Panda robot’s inverse dynamics model, and a joint residual mapping model and a hybrid inverse dynamics model of the industrial robot are established. Figure 11 shows the seven joints reproducing the actual torques after rearranging the errors to compensate for them in the order in which they were collected. Note that the range and accuracy of the vertical axis are different for each figure. Therefore, the seven joints cannot be compared. The results show that the reproduction results of the 7 joints were better, and the residual compensation model was effective.
4.4. Results analysis on the model test set
To validate the test results, the method in this paper was compared with the LSTM algorithm [Reference Zheng, Gong, Huo-Ming, Zhi-Cheng, Wen-Lin, Ji-min, Jian and Xing35] and the GP algorithm [Reference Rueckert, Nakatenus, Tosatto and Peters36], which are commonly used in previous studies. The results are shown in Tables VI and VII. After offline compensation, the algorithm proposed in this paper outperforms the LSTM and GP algorithms in terms of the MSE, RMSE, MAE, and R2. Compared with those of the LSTM, the mean values of the seven joints increase with respect to the MSE by 19.23%, RMSE by 12.68%, MAE by 21.69%, R2 by 5.84%, and the average residual reduction rate by 21.48%. Compared with those of the GP, the average increase in MSE is 90.89%, t in RMSE is 70.63%, in MAE is 65.23%, in R2 is 42.37%, and the residual reduction rate is 65.23%. Therefore, compared with the LSTM and GP methods, the RMSE, MAE, R2, and reduction rate of the average residual in the robot inverse dynamics model based on the BLL method proposed in this paper increased.
5. Conclusions
In order to improve the residual prediction accuracy of inverse dynamics models for industrial robots, this paper combines the bootstrap aggregating (bagging) algorithm with LSTM network, introduces linear layer as the network optimization layer, and proposes a method based on the BLL residual prediction algorithm, which is introduced as a network optimization layer. Additionally, hybrid inverse dynamics model compensation method for robots based on BLL residual prediction algorithm is proposed, and the framework of BLL residual prediction algorithm is given. The BLL residual prediction algorithm can achieve accurate prediction of residuals, and the algorithm can accurately estimate the residual value under the current value, which provides a prerequisite for the establishment of a hybrid inverse dynamics model.
To verify the effectiveness of the BLL residual model-based compensation method, we compared it with the rigid body inverse dynamics model. The torque residual of each robot joint decreased from 0.5651 N · m to 0.1096 N.m on average for the hybrid inverse dynamics model based on the BLL residual prediction algorithm. The calculation accuracy of the robot inverse dynamics model is improved, which is beneficial for aiding the robot to perform manipulation tasks more accurately. The BLL algorithm proposed in this paper was compared and analysed with the LSTM and GP algorithms in this field, and the BLL algorithm proposed in this paper achieved improvements in MSE, RMSE, MAE, R2, and the reduction rate of average residuals after compensation, which validated the effectiveness of the proposed algorithm in this paper. The study laid the foundation for accurately modelling the inverse dynamics of industrial robots.
Author contribution
Conceptualization, Y.T. and S.C.; methodology, Y.T.; software, H.L. and S.C.; validation, J.W.; formal analysis, S.C.; investigation, H.W.; resources, H.L.; data curation, S.C.; writing – original draft preparation, S.C.; writing – review and editing, S.C.; visualization, S.C.; supervision, J.W.; project administration, Y.T.; funding acquisition, Y.T. All authors have read and agreed to the published version of the manuscript.
Financial support
This research was funded by the Ministry of Industry and Information Technology of the People’s Republic of China, National Key Research and Development Plan “Intelligent Robot” Project No. 2022YFB4700400.
Competing interests
The authors declare no competing interests.