# On Hybrid Inverse Dynamic Modeling for Industrial Robots

Sacha Clavel, Mazen Alamir, Julien Favre-Faure

<figure><img src="image_1.png" /></figure>

## To cite this version:

Sacha Clavel, Mazen Alamir, Julien Favre-Faure. On Hybrid Inverse Dynamic Modeling for Industrial Robots. Mechatronics, 2026, 115 (April 2026), pp.103451. 10.1016/j.mechatronics.2025.103451. hal-05433203

HAL Id: hal-05433203

https://hal.science/hal-05433203v1

Submitted on 6 Jan 2026

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# On Hybrid Inverse Dynamic Modeling for Industrial Robots

S. Clavel* M. Alamir** J. Faure-Favre* S. Blanc*

* Stäubli company, 74210 Faverges, France (e-mail: s.clavel@staubli.com).

** University of Grenoble Alpes, CNRS, Grenoble INP, GIPSA-lab, 38000 Grenoble, France (e-mail: mazen.alamir@grenoble-inp.fr)

**Abstract:** Inverse dynamic modeling plays a crucial role in developing optimal feedforward control laws in robotics. A promising approach to improve its accuracy involves using *hybrid* approaches combining physics-based models with data-driven models. This paper presents an in-depth study of such an approach, based on Gated Recurrent Unit (GRU) neural networks for application on 4-joints and 6-joints Stäubli industrial robots. We start with a comprehensive review of the recent literature. Subsequently, we demonstrate that our model significantly surpasses more traditional black-box models in terms of accuracy and extrapolation. Furthermore, we explore the various factors that influence its accuracy with real data.

**Keywords:** Robots manipulators, Dynamic modeling, Gray-box modeling, Machine learning, Deep learning

## 1. INTRODUCTION

The dynamic modeling of industrial robots has been widely explored in the literature, notably due to its critical importance in achieving compliant, precise, and responsive control. Indeed, the inverse dynamic model is commonly used in robots' control laws as a feedforward mechanism, making its accurate identification crucial for attaining high control performances: the more precise the model, the better the control.

The conventional method for modeling the inverse dynamics of a robot relies on the derivation of *white-box* models, generally based on Eulerian, Newtonian, or Lagrangian physics equations. With such white-box approaches, the model's accuracy depends on both the granularity of the mathematical modeling and the precision of the physical parameters identification. In practice, due to these two factors, the accuracy of white-box models often reaches a limit. To avoid time-consuming modeling and identification efforts, only basic physical effects are modeled and physical parameters are frequently inaccurately identified.

To overcome these limitations, the use of machine learning and even deep learning algorithms is now a well established practice with rather remarkable success. Many data-driven algorithms have already demonstrated convincing results for modeling the dynamics of robots, including Locally Linear Regressors, Support Vector Regressors, Gaussian Process Regressors, or Neural Networks. However, Recurrent Neural Networks are considered highly promising due to their design, which enables them to effectively leverage temporal correlations in data, making them ideal for addressing tasks involving time-series. Popular recurrent architectures such as Long Short-Term Memory (LSTM) networks or Gated Recurrent Unit (GRU) networks have been designed to exploit both short-term and long-term

temporal dependencies, addressing challenges like vanishing or exploding gradients commonly faced when analyzing long time-series. The use of LSTMs or GRUs for dynamic modeling of robots enables reaching accuracy levels that outperform those achieved through traditional methods.

In fact, GRU cells are a newer and lighter variation of the well-known LSTM cells: whereas LSTMs have complex and heavy internal architectures, GRUs are characterized by their simpler design, reducing computational complexity while often achieving comparable accuracy (Chung et al., 2014; Mukhopadhyay et al., 2019). Thus, their streamlined architecture can be beneficial in embedded systems, such as robots, where speed, as well as memory, and computational efficiency are critical. Furthermore, since GRU structure has fewer parameters compared to deeper architectures, it may help limit overfitting. This could be advantageous for learning the dynamics of industrial robots that have highly variable operating conditions, in a context where data are often limited and noisy.

There are two main approaches: black-box modeling and gray-box modeling. First, black-box models aim at learning the whole robots' dynamics, without any prior information, completely replacing the physical equations. In their research, (Liu et al., 2019) studied a black-box model based on an LSTM network and reached high accuracies. The same approach was used by (Rueckert et al., 2017) with real data, who demonstrated its superiority in terms of accuracy and training time compared to a more conventional method based on Gaussian Processes (GPs). Similar conclusions were drawn by (Mukhopadhyay et al., 2019) where the authors compared different black-box models, based on GPs, simple RNNs, LSTMs, GRUs or basic Multi-Layer Perceptrons and showed that the use of LSTMs or GRUs generally leads to better results.
The second approach is the use of gray-box models, also known as *hybrid* models, which exploit the strengths of white-box and black-box models by combining physical equations with data-driven algorithms. A first way to do so is to build a residual model whose goal is to learn the errors of a prior physical model. This approach was used by (Liu et al., 2020) for example, where the authors built a residual LSTM model that led to significantly better accuracy than their previous black-box approach. The same method was used by (Wang et al., 2020), who trained a residual LSTM model with attention mechanisms for a 7-axis robot. Similarly, in (Valencia-Vidal et al., 2023), the authors learned residual bidirectional GRUs for a collaborative robot with elastic joints but did not observe a clear superiority of their hybrid models compared to their black-box version. In their research, (Çallar and Böttger, 2022) tried different models, including residual ones, to learn the dynamics of a 7-axis robot performing specific movements with highly non-linear dynamics. A second way to implement a gray-box model is by using physics-informed deep learning, which consists of integrating physics laws in the model architecture and the learning process. This type of work has been done by (Lutter et al., 2019, 2021) who integrated Lagrangian and Newtonian mechanics in deep learning models, by constraining their architecture and loss function based on physics equations. In their study, (Reuss et al., 2022) combined this method with a residual model based on LSTMs, resulting in a significant improvement compared to their LSTM black-box model. As a matter of fact, the study of hybrid models based on physics and data-driven algorithms has a growing interest in many fields of robotics research such as hydraulic manipulator control (Yao et al., 2024, 2025), soft robot modeling (Zhang et al., 2022), motion planning (Yang et al., 2024).

In the present contribution, we developed a hybrid model based on GRUs, to derive the inverse dynamics of a 4-axis TS2-80 and a 6-axis TX2-90 Stäubli robots. We explored a wide range of hyperparameters to tune its accuracy and compared it with a black-box approach. To give the reader an idea of the comprehensiveness of this paper, we have listed below the different research items that we have covered and Table 1 synthesizes, for each research item, which paper has addressed it. The research items are:

(a) Comparison of hybrid models with black-box models

(b) Impact of the temporal span of input time-series

(c) Impact of the input features

(d) Impact of the model's size

(e) Impact of the training set's size

<table><caption>Table 1. Research items addressed in the literature.</caption><thead><tr><th></th><th>(a)</th><th>(b)</th><th>(c)</th><th>(d)</th><th>(e)</th></tr></thead><tbody><tr><td>(Liu et al., 2019)</td><td>✗</td><td>✗</td><td>✓</td><td>✓</td><td>✗</td></tr><tr><td>(Liu et al., 2020)</td><td>✓</td><td>✗</td><td>✗</td><td>✗</td><td>✗</td></tr><tr><td>(Mukhopadhyay et al., 2019)</td><td>✗</td><td>✗</td><td>✗</td><td>✓</td><td>✓</td></tr><tr><td>(Reuss et al., 2022)</td><td>✓</td><td>✗</td><td>✗</td><td>✗</td><td>✗</td></tr><tr><td>(Rueckert et al., 2017)</td><td>✗</td><td>✗</td><td>✗</td><td>✓</td><td>✓</td></tr><tr><td>(Tao et al., 2024)</td><td>✗</td><td>✗</td><td>✗</td><td>✗</td><td>✗</td></tr><tr><td>(Valencia-Vidal et al., 2023)</td><td>✓</td><td>✓</td><td>✗</td><td>✓</td><td>✓</td></tr><tr><td>(Wang et al., 2020)</td><td>✗</td><td>✓</td><td>✓</td><td>✓</td><td>✗</td></tr><tr><td>(Çallar and Böttger, 2022)</td><td>✓</td><td>✓</td><td>✓</td><td>✓</td><td>✗</td></tr><tr><td><strong>This paper</strong></td><td>✓</td><td>✓</td><td>✓</td><td>✓</td><td>✓</td></tr></tbody></table>

The paper is organized as follows. First, Section 2 recalls the problem formulation and presents the architecture of the proposed model, then some aspects of its behavior are analysed in Section 3. Subsequently, Section 4 highlights the benefits of using a hybrid model towards a black-box model. Finally, Section 5 presents an extensive study of the main hyperparameters that affect the accuracy of the hybrid model before a conclusion is given.

## 2. METHODOLOGY

In this section, we introduce our hybrid model $M_H$, which is composed of a physics-based model $M_P$ and a residual data-driven model $M_θ$. Here, we first formulate the problem $M_H$ should solve, and then describe its architecture.

### 2.1 Problem Formulation

Our method aims at deriving the inverse dynamics of any $n$-joint robot. Such dynamics is classically described as a sum of physics effects:

$$ \tau = \underbrace{I(q)\ddot{q} + C(q, \dot{q}) + g(q) + F(q, \dot{q}, \ddot{q}, \dots)}_{\tau_p} + \tau_e \quad (1) $$

where $\tau \in \mathbb{R}^n$ represents the motor torques, $q, \dot{q}$ and $\ddot{q} \in \mathbb{R}^n$ are the joint positions, velocities and accelerations, $I(q) \in \mathbb{R}^{n \times n}$ stands for the symmetric positive-definite inertia matrix, $C$ denotes the Coriolis and centrifugal effects term, $g$ refers to the gravitational term, and $F(q, \dot{q}, \ddot{q}, \dots)$ concentrates additional *easy-to-model* physics effects (e.g. basic friction, springs, motor shafts inertia or external forces, etc.). The first four terms on the right side of (1) can be gathered in one term $\tau_p$ that represents all the *easy-to-model* torques, that is to say, all the torques that can simply be modeled with a white-box approach. Consequently, the term $\tau_e \in \mathbb{R}^n$ should be a small residual term which only includes really complex and *hard-to-model* physics phenomena (e.g., backslashes, joints' elasticities, complex friction effects, etc.).

To model this sum of effects, we first build the physics-based model $M_P$ to represent the terms gathered in $\tau_p$ with the classic physics equations, so that $M_P(q, \dot{q}, \ddot{q})$ provides $\hat{\tau}_p$ which is an estimate of $\tau_p$. The $P$ of $M_P$ represents all the identified physics parameters such as mass, gravity centers, inertia, or friction coefficients. Then, the data-driven model $M_θ$ is built to learn the residual effects that are not modeled by $M_P$, so that $M_θ(q, \dot{q}, \ddot{q}, \dots)$ provides $\hat{\tau}_e$ which is an estimation of $\tau_e$, or more exactly of $\tau_e := \tau - \hat{\tau}_p$. Here, the $\theta$ of $M_θ$ represents all the trainable parameters of the learning model. Please note that, in our work, the model is intended to be used in a pure feedforward, or open-loop, manner. Consequently, $q, \dot{q}$ and $\ddot{q}$ represent here the desired joint positions, velocities, and accelerations, and $\tau$ is the desired torque required by the joint motors and computed by the robot's PID controller. This detail is of great importance as this approach is underexplored in the literature, as it can be seen in Table 2 where a comparison with related hybrid dynamics studies is shown.

After having built $M_P$ and $M_θ$, the hybrid model $M_H$ is simply defined as $M_H := M_P + M_θ$, which predicts $\hat{\tau}$ which should be an accurate estimation of the torques
<table><caption>Table 2. Comparison to related hybrid dynamics studies made with real robots. The subscript letter *m* denotes measured signals while *d* denotes desired signals. Papers whose application have been labeled "Identification" do not specify explicitly applications but often mention control.</caption><thead><tr><th></th><th>Input features</th><th>Target variable</th><th>Approach</th><th>Application</th></tr></thead><tbody><tr><td>(Liu et al., 2020)</td><td>q<sub>m</sub>, ḋ<sub>m</sub>, ḋ̇<sub>m</sub></td><td>τ<sub>m</sub></td><td>Residual LSTM</td><td>Identification</td></tr><tr><td>(Reuss et al., 2022)</td><td>q<sub>m</sub>, ḋ<sub>m</sub>, ḋ̇<sub>d</sub></td><td>τ<sub>d</sub></td><td>Physics inspired with residual LSTM</td><td>Feedforward control</td></tr><tr><td>(Lutter et al., 2019)</td><td>q<sub>m</sub>, ḋ<sub>m</sub>, ḋ̇<sub>m</sub></td><td>τ<sub>d</sub></td><td>Physics-informed NN</td><td>Feedforward control</td></tr><tr><td>(Valencia-Vidal et al., 2023)</td><td>q<sub>m</sub>, ḋ<sub>m</sub></td><td>τ<sub>d</sub></td><td>Residual Bidirectionnal GRU</td><td>Feedforward control</td></tr><tr><td>(Wang et al., 2020)</td><td>q<sub>m</sub>, ḋ<sub>m</sub>, ḋ̇<sub>m</sub></td><td>τ<sub>m</sub></td><td>Residual LSTMs with Attention</td><td>Identification</td></tr><tr><td>(Çallar and Böttger, 2022)</td><td>q<sub>m</sub>, ḋ<sub>m</sub>, ḋ̇<sub>m</sub>, r<sub>m</sub></td><td>τ<sub>m</sub></td><td>Residual LSTMs or Transformer</td><td>Identification</td></tr><tr><td><strong>This paper</strong></td><td>q<sub>d</sub>, ḋ<sub>d</sub>, ḋ̇<sub>d</sub></td><td>τ<sub>d</sub></td><td>Residual GRU</td><td>Feedforward control and fault detection</td></tr></tbody></table>

$\tau$. This model should reach high accuracy levels, on one hand, since $M_\theta$ has been trained to specifically fit the inaccuracies of $M_P$. On the other hand, it should maintain good generalization and extrapolation capabilities, given that $M_P$ is a white-box model. Indeed, considering that $M_P$ is purely physics-based, its mathematical equations remain valid and invariant for a wide variety of operating conditions regardless of payload, position space, or trajectory profile. This contrasts with data-driven models that are prone to overfitting on specific operating conditions since they learn directly and blindly from data. Thus, the two models complement each other: where $M_\theta$ is generally fragile in terms of generalization and extrapolation, $M_P$ compensates. Conversely, where $M_P$ struggles to capture very complex physical phenomena, $M_\theta$ addresses these challenges.

In this paper, the models' accuracies are evaluated thanks to two dimensionless metrics that are the normalized Mean Absolute Error (nMAE) and the normalized Mean Squared Error (nMSE), defined as follows.

$$ \text{nMAE}(\hat{\tau}, \tau) = \frac{\sum_{i=0}^{T} |\hat{\tau}_i - \tau_i|}{\sum_{i=0}^{T} |\tau_i|} \qquad (2) $$

$$ \text{nMSE}(\hat{\tau}, \tau) = \frac{\sum_{i=0}^{T} (\hat{\tau}_i - \tau_i)^2}{\sum_{i=0}^{T} \tau_i^2} \qquad (3) $$

The next sections provide an in-depth exploration of the architecture of each component of $M_H$.

## 2.2 Physics-based component $M_P$

The general architecture of the physics-based model is described in Fig. 1. In the first layer of $M_P$, we use the Newton-Euler algorithm to compute $\tau_{jnt}$ which represents the joint torques induced by the inertial effects of the motions of each axis. Since each axis of the robot is connected to another, their dynamics have a strong coupling, which can be difficult to model as the number of axes increases. The Recursive Newton-Euler algorithm is a well-known method that allows us to model these complex dynamics in a computationally efficient way (Khalil and Dombre, 2004).

The second layer of $M_P$ involves a simple base change that brings the joint torques $\tau_{jnt}$, along with the joint positions $q$, velocities $\dot{q}$ and accelerations $\ddot{q}$ into the motor reference frame, to obtain the motor torques with motor positions, velocities and accelerations.

<figure><img src="image_3.png" /><figcaption>Fig. 1. General architecture of $M_P$.</figcaption></figure>

Following the transmission layer, friction layers are introduced. First, a component models viscous friction basically as a power of the motor velocities, for each axis $i$:

$$ \mathcal{V}_i(t) = k_i^{(1)} \times \text{sign}(\dot{q}_i(t)) \times |\dot{q}_i(t)|^{k_i^{(2)}} \qquad (4) $$

where $k_i^{(1)}, k_i^{(2)} \in \mathbb{R}$ represent the usual viscous friction parameters of the axis $i$. In parallel, a dry friction component is added which is somewhat more complex. First, if the axis $i$ is moving (i.e. $\dot{q}_i(t) \neq 0$), the dry friction torque is calculated based on a modified Coulomb model which is then filtered using a distance relaxation mechanism, as described in (5), (6) and (7). Conversely, if the axis $i$ is not moving (i.e. $\dot{q}_i(t) = 0$)<fnref n="1" />, the dry friction torque is progressively reduced to 0 following a time relaxation law as described in (8).

$$ \mathcal{D}_i(t)|_{\dot{q} \neq 0} = \tilde{\mathcal{D}}_i(t) + (\mathcal{D}_i(t)(t - \delta t) - \tilde{\mathcal{D}}_i(t)) \times e^{d_i(t)} \qquad (5) $$

$$ \text{with } \tilde{\mathcal{D}}_i(t) = \text{sign}(\dot{q}_i(t)) \times (k_i^{(3)} + k_i^{(4)} |\tau_{jnt_i}(t)|) \qquad (6) $$

$$ \text{and } d_i(t) = \frac{-|q_i(t) - q_i(t - \delta t)|}{k_i^{(5)}} \qquad (7) $$

$$ \mathcal{D}_i(t)|_{\dot{q}=0} = \mathcal{D}_i(t - \delta t) \times (1 - \frac{\delta t}{k_i^{(6)}}) \qquad (8) $$

Here, $\delta t$ denotes the controller sampling time. The parameter $k_i^{(3)}$ is the basic Coulomb parameter and $k_i^{(4)}$ enables us to model the dependence on the applied torque $\tau_{jnt}$. The parameter $k_i^{(5)}$ is the relaxation distance which enables us to model the hysteresis effects observed during small displacements, also known as Dahl friction. Finally,

<footnote n="1">¹ Note that this condition would never hold with measured signals, due to measurement noise. However, we recall that our models rely on desired signals so the nullification of $\dot{q}$ is possible.</footnote>
$k_i^{(6)}$ is the relaxation time. It should be noted that, to avoid cumbersome notation, the variable names $q_i$, $\dot{q}_i$ and $\tau_{\text{int}_i}$ are used in (4), (6) and (7) for simplicity. However, these variables correspond here to the positions, velocities, and torques within the motor reference frame due to the transmission layer and are no longer in the joint reference frame, as shown in Fig. 1.

Ultimately, for each axis $i$, a motor inertia layer models the motor shaft inertia with a basic proportional model:

$$ \mathcal{J}_i(t) = k_i^{(7)} \times \ddot{q}_i(t) \qquad (9) $$

where $k_i^{(7)}$ is the shaft inertia of motor $i$. After this layer, the estimation of the motor torque $\hat{\tau}_p$ is finally obtained by summing all the previously computed torques: the torques due to shaft inertia, viscous and dry frictions, and arm inertia.

The physical parameters $P$ of $M_P$ are either known by construction (for example, the mass, gravity center, and inertia matrix of the robot joints are known thanks to the CAD files) or have been identified using experimental procedures (such as friction parameters) by Stäubli. These identification procedures involve the acquisition of torque and motion data during specific identification trajectories (i.e., at specific speeds, in steady state, free of gravity if necessary, or during small displacements, etc.). Based on these specific acquisition points, the parameters are then identified with Grey-box identification tools. The exact details of these identification procedures are beyond the scope of this paper and will not be described here.

With this architecture, the physics-based model achieves good accuracy. The usual accuracy metrics values reached by $M_P$ can be seen in Table 4. It is important to note that given that we use a hybrid model, achieving high accuracy with $M_P$ is still of great interest. Indeed, as formulated in (1), $\tau = \tau_p + \tau_e$. Therefore, as the accuracy of $M_P$ increases, the residual torques decreases, so $M_\theta$ will only have to model small residual phenomena.

## 2.3 Data-driven component $M_\theta$

The goal of the data-driven model $M_\theta$ is to simulate the residual $\tau_e := \tau - \hat{\tau}_p$, and is mainly based on several layers of gated recurrent neural networks.

Gated recurrent cells, including GRUs and LSTMs cells, have been introduced to tackle the traditional challenges of vanishing and exploding gradients encountered when training deep networks with long time-series data. These networks are composed of neural layer combinations, commonly referred to as gates, which are designed to concentrate on relevant long-term and short-term dependencies while discarding irrelevant ones. The GRU cells are a more recent and simplified variation of the well-known LSTM cells. Their structure is recalled in Fig. 2.

In its final architecture, the model $M_\theta$ receives as input sequences of features $q$, $\dot{q}$, $\ddot{q}$, and $\hat{\tau}_p$, of length 128. These input sequences will be denoted respectively $\tilde{q}$, $\tilde{\dot{q}}$, $\tilde{\ddot{q}}$, and $\tilde{\hat{\tau}}_p$. These time-series first pass through 4 layers of GRU cells, each with a 64-dimensional hidden state, stacked on each other. The hidden state of the final GRU layer is then

<figure><img src="image_4.png" /><figcaption>Fig. 2. Structure of a GRU cell annotated with its elementary operators. The cell is based on 2 gates: the reset gate determines what information to discard from the previous hidden state $h_{t-1}$, and the input gate decides what information from the input $x_t$ to store in the current hidden state $h_t$.</figcaption></figure>

fed to a Fully Connected (FC) layer which predicts $\hat{\tau}_e$, the $n$-dimensional vector that estimates the residual $\tau_e$. This final configuration required several experiments, which will be discussed in Section 5.

To train the model $M_\theta$, large datasets containing time-series of $q$, $\dot{q}$, $\ddot{q}$, $\hat{\tau}_p$ for model input and $\tau_e$ for model output must be gathered from the robots. Two datasets were collected, on a 6-axis TX2-90 and a 4-axis TS2-80 Stäubli robots, which enabled us to train a $M_\theta$ for each type of robot. To ensure the models’ robustness across different robots’ trajectories, these datasets must encompass a wide variety of motion scenarios. To do so, the robots performed random point-to-point movements with randomly chosen points in almost their entire workspace, random velocities, and accelerations. During these random movements, the required time-series of $q$, $\dot{q}$, $\ddot{q}$, and $\tau$ were sampled at 250Hz, then the remaining required signals were post-calculated. First, $\hat{\tau}_p$ was computed with $M_P(q, \dot{q}, \ddot{q})$, then the residual torque $\tau_e$ was finally computed as the difference $\tau - \hat{\tau}_p$. No filter was applied to any of the signals. A total of 200 random trajectories were executed for each robot using this data collection process, which led to a first dataset containing 1,623,000 samples for the TS2-80 and a second of 2,575,000 samples for the TX2-90. These datasets are then normalized with the usual Min-Max normalizer and split into 3 subsets: 80% of the data are kept for the training set, 17% of the data are allocated for the test set, and the remaining 3% form a small validation set.

The models are then trained, offline, to perform a regression task with the MAE loss function. The optimization of the trainable parameters $\theta$ of $M_\theta$ is carried out with the well-known Adam optimizer with an exponentially decreasing learning rate from $10^{-4}$ to $10^{-5}$ and the default parameters $\beta_1 = 0.9$, $\beta_2 = 0.999$ and $\epsilon = 10^{-8}$. The mini-batch gradient descent is performed with batches of size 64 and a maximum of 200 epochs. An early stopping mechanism based on the validation loss is settled to prevent overfitting: the model that reaches the minimal
<table><caption>Table 3. Datasets description.</caption><thead><tr><th></th><th>TX2-90</th><th>TS2-80</th></tr></thead><tbody><tr><td>Movements per trajectory</td><td>30</td><td>24</td></tr><tr><td>Motion type</td><td>Joint space</td><td>Joint space</td></tr><tr><td>Joints limits</td><td>q<sup>(1)</sup>: [-180°, 180°]<br>q<sup>(2)</sup>: [-90°, 90°]<br>q<sup>(3)</sup>: [-60°, 60°]<br>q<sup>(4)</sup>: [-180°, 180°]<br>q<sup>(5)</sup>: [-90°, 90°]<br>q<sup>(6)</sup>: [-180°, 180°]</td><td>q<sup>(1)</sup>: [-180°, 180°]<br>q<sup>(2)</sup>: [-120°, 120°]<br>q<sup>(3)</sup>: [0m, 0.2m]<br>q<sup>(4)</sup>: [-500°, 500°]</td></tr><tr><td>Number of trajectories</td><td>200</td><td>200</td></tr><tr><td>Sampling frequency</td><td>250Hz</td><td>250Hz</td></tr><tr><td>Sampled signals</td><td>q<sup>(i)</sup>, ḋ<sup>(i)</sup>, ḋ̃<sup>(i)</sup>, τ<sup>(i)</sup><br>for i ∈ {1, ..., 6}</td><td>q<sup>(i)</sup>, ḋ<sup>(i)</sup>, ḋ̃<sup>(i)</sup>, τ<sup>(i)</sup><br>for i ∈ {1, ..., 4}</td></tr><tr><td>Post-calculated signals</td><td>τ̂<sub>p</sub><sup>(i)</sup> and τ<sub>e</sub><sup>(i)</sup><br>for i ∈ {1, ..., 6}</td><td>τ̂<sub>p</sub><sup>(i)</sup> and τ<sub>e</sub><sup>(i)</sup><br>for i ∈ {1, ..., 4}</td></tr><tr><td>Total samples</td><td>2.575×10<sup>6</sup></td><td>1.623×10<sup>6</sup></td></tr><tr><td>Duration</td><td>171 min</td><td>108 min</td></tr></tbody></table>

validation loss is considered optimal. Once training is done, test loss and test metrics are computed in order to assess the performance of the obtained model and to compare various hyperparameter configurations, as will be done in Sections 4 and 5.

## 2.4 Resulting Hybrid model

Once $M_P$ and $M_θ$ are built, the resulting hybrid model can be formed as described in Fig. 3. Since $M_θ$ takes time-series as input, a buffer mechanism is added after $M_P$ to create the appropriate sequences of the features $q, \dot{q}, \ddot{q}$, and $\hat{\tau}_p$. Please note that $M_θ$ takes as input the data of all the axes of the robot ($q, \dot{q}, \ddot{q}, \hat{\tau}_p \in \mathbb{R}^n$) to produce the residual torque estimate of all the axes ($\hat{\tau}_e \in \mathbb{R}^n$). The sum of $\hat{\tau}_p$ and $\hat{\tau}_e$ gives $\hat{\tau}$ which should be a robust and precise estimation of $\tau$ defined in (1). We recall here that $q, \dot{q}, \ddot{q}$ are the desired joint's positions, velocities, and accelerations, while $\hat{\tau}$ is an estimation of the desired motor torques.

<figure><img src="image_5.png" /><figcaption>Fig. 3. Final architecture of $M_H$.</figcaption></figure>

With this hybrid model, significantly higher accuracy can be reached compared to the initial white-box model $M_P$, as we can see in Table 4.

## 3. DETAILED ANALYSIS OF $M_H$ BEHAVIOR

This section provides insights and results regarding the behavior of $M_H$. A brief study on the contribution of each prior physical sub-model is first proposed, followed by a visualization of the phenomena not captured by the model. Finally, results from a short feedforward control experiment utilizing $M_H$ are presented.

<table><caption>Table 4. Accuracy metrics of $M_P$ and $M_H$ for the two robots on the 40 test trajectories. The rows correspond to the different joints of the robots. The columns correspond to the metrics computed between $\tau$ and $\hat{\tau}_p$ for $M_P$ and between $\tau$ and $\hat{\tau}$ for $M_H$.</caption><thead><tr><th></th><th></th><th colspan="2">nMAE</th><th colspan="2">nMSE</th></tr><tr><th></th><th></th><th>M<sub>P</sub></th><th>M<sub>H</sub></th><th>M<sub>P</sub></th><th>M<sub>H</sub></th></tr></thead><tbody><tr><td rowspan="6">06-2XST</td><td>1</td><td>0.179 ± 0.009</td><td>0.139 ± 0.008</td><td>0.033 ± 0.003</td><td>0.020 ± 0.003</td></tr><tr><td>2</td><td>0.114 ± 0.007</td><td>0.035 ± 0.002</td><td>0.015 ± 0.001</td><td>0.002 ± 0.000</td></tr><tr><td>3</td><td>0.110 ± 0.013</td><td>0.050 ± 0.004</td><td>0.015 ± 0.003</td><td>0.003 ± 0.000</td></tr><tr><td>4</td><td>0.287 ± 0.082</td><td>0.225 ± 0.067</td><td>0.158 ± 0.100</td><td>0.117 ± 0.082</td></tr><tr><td>5</td><td>0.399 ± 0.030</td><td>0.105 ± 0.009</td><td>0.155 ± 0.019</td><td>0.013 ± 0.002</td></tr><tr><td>6</td><td>0.378 ± 0.038</td><td>0.116 ± 0.007</td><td>0.123 ± 0.021</td><td>0.016 ± 0.002</td></tr><tr><td rowspan="4">08-2SST</td><td>1</td><td>0.211 ± 0.007</td><td>0.055 ± 0.004</td><td>0.044 ± 0.002</td><td>0.003 ± 0.000</td></tr><tr><td>2</td><td>0.183 ± 0.007</td><td>0.064 ± 0.005</td><td>0.040 ± 0.003</td><td>0.004 ± 0.001</td></tr><tr><td>3</td><td>0.193 ± 0.012</td><td>0.078 ± 0.006</td><td>0.042 ± 0.004</td><td>0.007 ± 0.001</td></tr><tr><td>4</td><td>0.231 ± 0.012</td><td>0.188 ± 0.012</td><td>0.067 ± 0.007</td><td>0.050 ± 0.006</td></tr></tbody></table>

### 3.1 Contributions of prior physical models

To enhance the explainability of our hybrid architecture, the impact of each prior physical sub-model has been studied. For this purpose, $M_H$ was tested while progressively activating the sub-models of $M_P$: starting with none activated (Ø), then activating only the Newton-Euler algorithm (NE), adding friction models (NE+F), and finally adding motor shaft inertia models (NE+F+MS). For each configuration of $M_P$, the residual model $M_θ$ was trained accordingly. This experiment was conducted on the TS2-80 dataset. The metrics values reached by each configuration on the test trajectories are gathered in Table 5.

<table><caption>Table 5. Impact of each prior model. Results presented for the nMAE metric only.</caption><thead><tr><th></th><th></th><th colspan="4">M<sub>P</sub></th></tr><tr><th></th><th></th><th>Ø</th><th>NE</th><th>NE+F</th><th>NE+F+MS</th></tr></thead><tbody><tr><td rowspan="4">08-2SST</td><td>1</td><td>1.0 ± 0.0</td><td>0.562 ± 0.033</td><td>0.326 ± 0.005</td><td>0.211 ± 0.007</td></tr><tr><td>2</td><td>1.0 ± 0.0</td><td>0.541 ± 0.044</td><td>0.204 ± 0.009</td><td>0.183 ± 0.007</td></tr><tr><td>3</td><td>1.0 ± 0.0</td><td>0.577 ± 0.046</td><td>0.203 ± 0.014</td><td>0.193 ± 0.012</td></tr><tr><td>4</td><td>1.0 ± 0.0</td><td>0.981 ± 0.024</td><td>0.353 ± 0.025</td><td>0.231 ± 0.012</td></tr><tr><th></th><th></th><th colspan="4">M<sub>H</sub></th></tr><tr><th></th><th></th><th>Ø</th><th>NE</th><th>NE+F</th><th>NE+F+MS</th></tr><tr><td rowspan="4">08-2SST</td><td>1</td><td>0.083 ± 0.005</td><td>0.059 ± 0.005</td><td>0.057 ± 0.004</td><td>0.055 ± 0.004</td></tr><tr><td>2</td><td>0.081 ± 0.005</td><td>0.061 ± 0.005</td><td>0.061 ± 0.006</td><td>0.064 ± 0.005</td></tr><tr><td>3</td><td>0.082 ± 0.007</td><td>0.083 ± 0.007</td><td>0.079 ± 0.007</td><td>0.078 ± 0.006</td></tr><tr><td>4</td><td>0.191 ± 0.012</td><td>0.189 ± 0.012</td><td>0.189 ± 0.012</td><td>0.188 ± 0.012</td></tr></tbody></table>

The results show, for $M_P$, that the Newton-Euler algorithm first captures a significant part of the torques, except on the last joint. Then, the friction and motor shaft inertia models have reasonable contributions, especially on the 4th joint, where they help to reduce the residual to a level comparable to the other joints, revealing that the dynamics of this joint is particularly governed by friction and motor's inertia. In contrast, within the overall hybrid architecture, gradually activating friction and motor shaft inertia models yields only minor accuracy gains. This suggests that $M_θ$ is greatly assisted by the Newton-Euler algorithm, which notably models the complex coupled dynamics of robot joints, but has no difficulty in replacing the terms of friction and motor inertia if they are lacking from the physical model.
## 3.2 Vibratory effects

By looking closely at the predictions of $M_H$, we found that our model was particularly good at simulating low-frequency phenomena, but had difficulty to capture some vibratory effects. For example, this can be seen in the carefully chosen zooms on the predictions of $M_H$ for the TX2-90 dataset, in Fig. 4. First, it is clear that $M_\theta$ succeeds in compensating the inaccuracies of $M_P$. However, we can see that the model struggles to reproduce the vibratory phenomena starting at 191.7s and is finally completely out of touch at 192s, especially for the joints 1 and 6.

These observations are confirmed by looking at the Power Spectral Densities of our model on joint 1 and 6, shown in Fig. 5. We see that $M_H$ manages to model precisely the effects below 10Hz, but misses important parts of the frequency content between 10 and 50Hz. These vibration phenomena could be caused by joint elasticity which the model struggles to capture for long durations, probably because of vanishing gradients.

<figure><img src="image_6.png" /><figcaption>Fig. 4. Zooms on $M_H$ predictions for the TX2-90 robot.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Fig. 5. Power spectral densities for joints 1 and 6.</figcaption></figure>

## 3.3 Feedforward control experiment

To validate the relevance of our model, we used it in a short feedforward control experiment on the TS2-80. A simplified diagram of the initial control loop of Stäubli robots is shown in Fig. 6. It is made up of a classic PID equipped with a feedforward torque unit named Advanced Torque Control (ATC).

<figure><img src="image_8.png" /><figcaption>Fig. 6. Control loop of Stäubli robots.</figcaption></figure>

In our experiment, we ordered the robot to perform exactly the same predefined movements under the same operating conditions, first with the original control loop and then with the ATC unit being replaced by $M_H$. The motions were new motions performed in the joint space, unseen by our model during its training, and lasted 31 seconds. During these motions, we recorded the position tracking error signal $q_e$ at 250Hz. Excerpts of the measured absolute tracking error can be compared in Fig. 7.

<figure><img src="image_9.png" /><figcaption>Fig. 7. Tracking error comparison.</figcaption></figure>

The average values of $|q_e|$ obtained for each joint are presented in Table 6. We can observe a general improvement in the tracking errors in almost all joints, which is an encouraging result that validates the relevance of the model in a feedforward control context. Nevertheless, a slight deterioration can be seen for the third joint. This side effect is still under investigation and could be explained by the fact that this axis is the only translation axis of the robot. The mechanical transmission of this axis is provided by a ball screw system, making its dynamics unique and significantly influenced by friction phenomena, which could hinder the modeling task of $M_H$.

# 4. BLACK-BOX VS HYBRID

We also wanted to validate the relevance of using our hybrid model rather than a black-box model. To do so,
<table><caption>Table 6. Mean absolute tracking errors.</caption><thead><tr><th></th><th></th><th>with ATC</th><th>with M<sub>H</sub></th><th>Unit</th></tr></thead><tbody><tr><td rowspan="4">TS2-80</td><td>1</td><td>0.0098</td><td>0.0058</td><td>rad</td></tr><tr><td>2</td><td>0.035</td><td>0.023</td><td>rad</td></tr><tr><td>3</td><td>0.0036</td><td>0.0041</td><td>mm</td></tr><tr><td>4</td><td>0.0019</td><td>0.0016</td><td>rad</td></tr></tbody></table>

we trained a model with the same architecture as $M_θ$ to predict the whole signal $τ$ instead of the residual $τ - \hat{τ}_p$. This model, called $M_B$, is a black box version of $M_H$ since it uses the same input features $q, \dot{q}$, and $\ddot{q}$ to model the same target variable $τ$. This model was trained with the exact same optimization configuration used for $M_θ$ (i.e. same optimizer, same number of epochs, same early stopping mechanism, etc.). We performed this experiment while increasing the number of trainable parameters of the two models to see if a black-box model with a sufficient complexity would be able to reach a better or at least similar precision than a hybrid model.

The results are shown in Fig. 8, with the displayed metric being the nMAE obtained on the TS2-80 test dataset, averaged across all axes of the robot. The right y-axis shows the relative accuracy deviation for each model compared to the optimal model. This second y-axis will be present in all remaining graphs to facilitate interpretation. The results reveal that the black-box version has significantly lower prediction performances than the hybrid model, even when increasing the model's complexity up to more than a million parameters. The use of black-box models leads to a loss of at least 8% in accuracy regardless of the number of trainable parameters. This confirms the idea that all the prior knowledge about Stäubli robots contained in $M_P$ is of great importance for $M_H$.

<figure><img src="image_10.png" /><figcaption>Fig. 8. Comparison of $M_H$ and $M_B$ on the test dataset.</figcaption></figure>

## 4.1 Generalization to a new motion type

Subsequently, we intended to verify the hypothesis that $M_H$ had better generalization capabilities than $M_B$. To achieve this, a first experiment on the TS2-80 was conducted and consisted in validating the models on an entirely new type of motion. Therefore, the robot executed predefined movements in the linear space, contrasting with the training movements performed in the joint space. These new movements differed significantly in terms of positions, velocities, and accelerations from those used for training. The trajectory lasted approximately 12 seconds and the appropriate data were sampled at 250Hz.

The performance metrics of both models during this trajectory are presented in Table 7. We see a clear superiority

of $M_H$, especially on the first three axes. This suggests that $M_H$ is still more accurate than $M_B$ when it comes to generalize. In future work, this validation experiment will be performed on longer trajectories.

<table><caption>Table 7. Generalization metrics on a movement in the linear space.</caption><thead><tr><th rowspan="2"></th><th rowspan="2"></th><th colspan="2">nMAE</th><th colspan="2">nMSE</th></tr><tr><th>M<sub>B</sub></th><th>M<sub>H</sub></th><th>M<sub>B</sub></th><th>M<sub>H</sub></th></tr></thead><tbody><tr><td rowspan="4">TS2-80</td><td>1</td><td>0.118</td><td>0.098</td><td>0.016</td><td>0.010</td></tr><tr><td>2</td><td>0.190</td><td>0.177</td><td>0.033</td><td>0.028</td></tr><tr><td>3</td><td>0.170</td><td>0.156</td><td>0.034</td><td>0.029</td></tr><tr><td>4</td><td>0.320</td><td>0.309</td><td>0.114</td><td>0.109</td></tr></tbody></table>

## 4.2 Extrapolation on unknown position space

A second and more robust experiment was conducted to validate the superiority of $M_H$ for extrapolation tasks. For this experiment, a third training dataset was collected with the exact same acquisition settings as the initial TS2-80 dataset described in Table 3, except for the joints limits. The limits of the joint positions were deliberately reduced such that $q^{(1)}, q^{(2)}, q^{(4)} \in [-20°, 20°]$, and $q^{(3)} \in [0m, 0.04m]$.

The two models $M_H$ and $M_B$ were first trained with this dataset collected on a limited workspace. The models were then validated on the 40 test trajectories of the initial dataset collected in a wide workspace, which therefore includes a large majority of positions not seen by the models during training. Specifically, only 17 samples out of the $3.12 \times 10^5$ test samples were collected with the robot's joints positions located within the training set limits. This dataset was therefore legitimate for assessing the extrapolation abilities of the models on completely unknown positions.

The performance metrics obtained during this experiment are shown in Table 8. The accuracy gap between $M_H$ and $M_B$ is even clearer in this experiment: a ratio of at least 3 is observed across all metrics. These results confirm the great superiority of $M_H$ for extrapolating and emphasizes the advantage of using a physics-based component.

<table><caption>Table 8. Extrapolation metrics on positions not seen during training.</caption><thead><tr><th rowspan="2"></th><th rowspan="2"></th><th colspan="2">nMAE</th><th colspan="2">nMSE</th></tr><tr><th>M<sub>B</sub></th><th>M<sub>H</sub></th><th>M<sub>B</sub></th><th>M<sub>H</sub></th></tr></thead><tbody><tr><td rowspan="4">TS2-80</td><td>1</td><td>1.090 ± 0.201</td><td>0.316 ± 0.044</td><td>1.214 ± 0.452</td><td>0.084 ± 0.022</td></tr><tr><td>2</td><td>1.060 ± 0.180</td><td>0.341 ± 0.041</td><td>1.189 ± 0.418</td><td>0.103 ± 0.020</td></tr><tr><td>3</td><td>0.973 ± 0.170</td><td>0.262 ± 0.036</td><td>1.047 ± 0.268</td><td>0.076 ± 0.014</td></tr><tr><td>4</td><td>2.353 ± 0.199</td><td>0.772 ± 0.074</td><td>5.439 ± 0.738</td><td>0.585 ± 0.100</td></tr></tbody></table>

# 5. TUNING OF HYPERPARAMETERS

Here we present some of the experiments that we have conducted in order to find the best hyperparameter configuration that minimizes the precision metrics computed on the test data. These experiments were essentially conducted with the data collected on the TS2-80 robot, however we believe the results are transposable to some extent.

## 5.1 Model's size

First, the most crucial hyperparameters of $M_θ$ that must be tuned are those determining the complexity of the
model. This can be done with the parameter `n layer` which corresponds to the number of GRU layers used in the model, and the parameter `hidden dim` which corresponds to the dimension of the GRUs' hidden states. These two parameters drastically influence the number of trainable parameters of $M_θ$. For example, a configuration with `n layer = 4` and `hidden dim = 256` has nearly $1.6 \times 10^6$ trainable parameters, while taking `n layer = 1` with `hidden dim = 16` leads to only $1.2 \times 10^4$ trainable parameters. Therefore, these must be cautiously tuned to achieve good predicting performances while avoiding overfitting.

To do so, several experiments have been conducted by varying these two parameters, Fig. 9 shows some of the results obtained for `n layer` varying from 1 to 4 and `hidden dim` varying from 16 to 128. The results show that, beyond 64 hidden dimensions, adding more layers does not have an impact on the model's accuracy and can even lead to overfitting. Based on Fig. 9, the optimal architecture is composed of 4 layers with 64 hidden dimensions and is therefore the configuration that will be used in the following sections.

<figure><img src="image_11.png" /><figcaption>Fig. 9. Impact of `n layer` and `hidden dim`.</figcaption></figure>

## 5.2 Input features and temporal span

We also wanted to determine the optimal input sequences for $M_θ$. First, the temporal span of the input sequences had to be tuned. This was done by adjusting the `seq len` hyperparameter which represents the length of the input sequences. For example, a `seq len = 125` corresponds to input time-series that have a temporal span of 0.5s, since the data are sampled at 250Hz. A second way to tune this parameter was to adjust the hyper-parameter `x sampling` that corresponds to the subsampling ratio of the input data. For example, time-series of `seq len = 125` with `x sampling = 1` (i.e. no subsampling) correspond to 0.5s, while the same time-series with `x sampling = 2` correspond to 1.0s, etc. Several experiments were conducted varying these two parameters, Fig. 10 shows some results obtained for `seq len` varying from 4 to 256 and `x sampling` varying from 1 to 4. The results show that, once the corresponding temporal span is sufficiently large, these two hyperparameters do not greatly affect the model's accuracy. However, using single point, non-sequential inputs (i.e. `seq len = 1`) would lead to poor precision performances. Finally, the optimal configuration is achieved rapidly with `seq len` set to 128 and `x sampling` set to 2, which corresponds to a temporal span of 1.024s.

<figure><img src="image_12.png" /><figcaption>Fig. 10. Impact of `seq len` and `x sampling`.</figcaption></figure>

We also had to choose which input features were the best to use, among all the collected signals. To do so, we performed experiments using various subsets of input features. Table 9 shows the performance metrics obtained for several subsets. The results show that the use of $τ_p$ has a strong impact on the accuracy of the model. Indeed, once this feature is used, adding more input features leads only to slight improvements in accuracy. This could be explained by the fact that this feature is a combination of the others, since it is computed by $M_P$ with $q, \dot{q}$ and $\ddot{q}$. However, it can be observed that the use of all features $q, \dot{q}, \ddot{q}$ and $τ_p$ seems to be the optimal configuration.

<table><caption>Table 9. Impact of the choice of input features. The metrics are averaged over all robot's axes.</caption><thead><tr><th></th><th>nMAE</th><th>nMSE</th></tr></thead><tbody><tr><td>q̇</td><td>0.116</td><td>2.11 × 10<sup>-3</sup></td></tr><tr><td>q̇̇</td><td>0.111</td><td>1.90 × 10<sup>-3</sup></td></tr><tr><td>q̇̇̇</td><td>0.108</td><td>1.74 × 10<sup>-3</sup></td></tr><tr><td>τ̂<sub>p</sub></td><td>0.102</td><td>1.14 × 10<sup>-3</sup></td></tr><tr><td>q, q̇, q̇̇</td><td>0.100</td><td>1.13 × 10<sup>-3</sup></td></tr><tr><td>q, τ̂<sub>p</sub></td><td>0.099</td><td>1.07 × 10<sup>-3</sup></td></tr><tr><td>q, q̇, q̇̇, τ̂<sub>p</sub></td><td>0.098</td><td>1.04 × 10<sup>-3</sup></td></tr></tbody></table>

## 5.3 Size of the training set

Finally, given the complexity of the chosen model, we wanted to know how much data was necessary to properly train it. This question is of interest, as collecting large training datasets on a robot can take several hours, and Stäubli might need to minimize this data collection time. The amount of data required depends, of course, on the richness of the scenarios realized by the robot. However, we believe that it could be useful for the community to know, for this type of robot and model, an order of magnitude of the time required to collect a sufficient amount of data with the robot performing random trajectories.

Experiments were carried out with various portions of the original TS2-80 dataset of 1,623,000 samples. Fig. 11 shows the performance metrics reached by different models, for different portions starting from 5% (equivalent to 4 minutes of recording) to 100% (83 minutes of recording). The x-axis has been converted in equivalent recording time for easier interpretation.

The results show that good prediction performances are rapidly achieved, with only 50% of our initial dataset, which corresponds to 43 minutes of data recording on the robot. Beyond this point, the precision gain becomes
<figure><img src="image_13.png" /><figcaption>Fig. 11. Impact of the training set's size.</figcaption></figure>

very small. We can also notice that models of varying complexities learn equivalently, regardless of the amount of data used. Therefore, the idea of collecting more data to use more complex models with the aim of achieving better performances would be pointless: it seems that we have reached a limit in terms of accuracy with this hybrid structure.

Naturally, industrial contexts could prohibit hours of data acquisition due to costs caused by downtime or because of technical issues related to the production environment. Strategies to reduce the required amount of operational data could therefore be considered. As an example, one strategy would involve pre-training a base model using an extensive dataset historically acquired by Stäubli. This model would then be fine-tuned with transfer learning techniques, using small datasets quickly collected at the end customer's production site. Incremental learning strategies could also be used to ensure the consistency of the model throughout the robot's lifetime. The base model would be periodically re-trained with the original dataset judiciously updated with small amounts of data frequently collected on the robot during its use on site.

## 6. CONCLUSION AND FUTURE WORK

This paper has proposed an extensive study of hybrid models based on GRU neural networks for accurately capturing the dynamics of Stäubli robots in a complete feedforward manner, which has never been done in the literature, as far as we know. The key hyperparameters to optimize its precision have been presented and analyzed. In addition, valuable insights have been proposed on the explainability and behavior of such hybrid architectures, and their great superiority over black-box architectures has been demonstrated and validated in generalization and extrapolation experiments.

The use of the model for feedforward control seems promising and will be further explored in the future. Future work could also involve the study of transfer learning or incremental learning to address the challenges associated with real world industrial contexts. Finally, utilizing the model for anomaly detection tasks is highly considered.

## REFERENCES

* Chung, J., Gülçehre, Ç., Cho, K., and Bengio, Y. (2014). Empirical evaluation of gated recurrent neural networks on sequence modeling. *CoRR*, abs/1412.3555.

* Khalil, W. and Dombre, E. (2004). *Modeling, identification and control of robots*.

* Liu, N., Li, L., Hao, B., Yang, L., Hu, T., Xue, T., and Wang, S. (2019). Modeling and simulation of robot inverse dynamics using lstm-based deep learning algorithm for smart cities and factories. *IEEE Access*, 7, 173989–173998.

* Liu, N., Li, L., Hao, B., Yang, L., Hu, T., Xue, T., Wang, S., and Shao, X. (2020). Semiparametric deep learning manipulator inverse dynamics modeling method for smart city and industrial applications. *Complexity*, 2020, 1–11.

* Lutter, M., Ritter, C., and Peters, J. (2019). Deep lagrangian networks: Using physics as model prior for deep learning. *CoRR*, abs/1907.04490.

* Lutter, M., Silberbauer, J., Watson, J., and Peters, J. (2021). A differentiable newton-euler algorithm for real-world robotics. *CoRR*, abs/2110.12422.

* Mukhopadhyay, R., Chaki, R., Sutradhar, A., and Chattopadhyay, P. (2019). Model learning for robotic manipulators using recurrent neural networks. In *TENCON 2019 IEEE Region 10 Conference*, 2251–2256.

* Reuss, M., van Duijkeren, N., Krug, R., Becker, P., Shaj, V., and Neumann, G. (2022). End-to-end learning of hybrid inverse dynamics models for precise and compliant impedance control.

* Rueckert, E., Nakatenus, M., Tosatto, S., and Peters, J. (2017). Learning inverse dynamics models in $o(n)$ time with lstm networks.

* Tao, Y., Chen, S., Liu, H., Wan, J., Wei, H., and Wang, T. (2024). Robot hybrid inverse dynamics model compensation method based on the bll residual prediction algorithm. *Robotica*, 1–18.

* Valencia-Vidal, B., Ros, E., Abadia, I., and Luque, N. (2023). Bidirectional recurrent learning of inverse dynamic models for robots with elastic joints: a real-time real-world implementation. *Frontiers in Neurorobotics*, 17.

* Wang, S., Shao, X., Yang, L., and Liu, N. (2020). Deep learning aided dynamic parameter identification of 6-dof robot manipulators. *IEEE Access*.

* Yang, W.T., Liao, J.M., and Lin, P.C. (2024). Physics-data hybrid dynamic model of a multi-axis manipulator for sensorless dexterous manipulation and high-performance motion planning.

* Yao, Z., Liang, X., Wang, S., and Yao, J. (2025). Model-data hybrid driven control of hydraulic euler-lagrange systems. *IEEE/ASME Transactions on Mechatronics*, 30(1), 131–143.

* Yao, Z., Xu, F., Jiang, G.P., and Yao, J. (2024). Data-driven control of hydraulic manipulators by reinforcement learning. *IEEE/ASME Transactions on Mechatronics*, 29(4), 2673–2684.

* Zhang, Z., Jin, Z., and Gu, G.X. (2022). Efficient pneumatic actuation modeling using hybrid physics-based and data-driven framework. *Cell Reports Physical Science*, 3(4), 100842.

* Çallar, T.C. and Böttger, S. (2022). Hybrid learning of time-series inverse dynamics models for locally isotropic robot motion. *IEEE Robotics and Automation Letters*, 8, 1061–1068.