2.1. Physical dynamic model

Use the Lagrangian approach to establish the robot dynamic model.

(1) \begin{equation}M\!\left(q\right)\ddot{q}+C\!\left(q,\dot{q}\right)\dot{q}+G\!\left(q\right)=\tau +\tau _{e}\end{equation}

(2) \begin{equation}\tau =K_{t}I\end{equation}

where $M(q)$ , $C(q,\dot{q})$ , and $G(q)$ represent the inertia matrix, the Coriolis matrix, and the gravity matrix, respectively, $q$ , $\dot{q}$ , and $\ddot{q}$ represent the joint position, joint velocity, and joint acceleration, respectively, $\tau _{e}$ is the external torque applied at the end of the joint, $\tau$ is the drive torque of joint module which is proportion of drive current according to Eq. (2), and $K_{t}$ is the torque coefficient.

For a lightweight robot, the whole body can be regarded as an integration of each joint module. Hence, the single joint module’s model and its control are the foundation of the whole robot’s performance. For the joint module, a test platform is established as shown in Figure 3 which includes the motor, reducer, coupling, torque sensor, and load which can simulate different external work conditions.

Figure 3. Schematic diagram of a single joint system.

Considering joint friction, the dynamic model in Figure 3 is represented by Eq. (3). Since the motor rotates around a fixed horizontal axis, both the Coriolis matrix $C(q,\dot{q})$ and the gravity matrix $G(q)$ in Eq. (1) are zero. The inertia matrix $M(q)$ has only one term, $I_{zz}$ , which represents the system’s moment of inertia about the z-axis.

(3) \begin{equation}\tau +\tau _{e}=I_{zz}\ddot{q}+\tau _{f}\end{equation}

(4) \begin{equation}\tau _{f}=f_{c}\mathrm{sgn}\!\left(\dot{q}\right)+f_{v}\dot{q}\end{equation}

where $\tau _{f}$ is the friction torque, described by the classical Coulomb-viscous friction model, as shown in Eq. (4), and $f_{c}$ and $f_{v}$ represent the coefficient of Coulomb friction and viscous friction, respectively.

However, some issues should be mentioned. The Coulomb-Viscous friction model is only suitable for describing simple frictional behaviors. However, in practical applications, friction forces are often more complex. The Coulomb-Viscous model has limitations when addressing these complexities in unstable or nonlinear friction phenomena. For instance, the friction force exhibits significant variation across different velocity ranges, the Stribeck effect [Reference Liu, Wei, Ni, Lu, Zhang and Li20]. It’s difficult for the traditional Coulomb friction model to accurately capture this nonlinear characteristic, especially in low speeds or static states, where friction forces may exhibit dramatic changes.

2.2. Parameter identification of physical dynamic model

The single joint module dynamic model which is described as Eq. (3) is linear so that its dynamic parameters can be identified using the least squares method (LSM) directly. The core idea of LSM is to minimize the sum of squared residuals between the actual observed values and the model estimates to seek optimal parameters fitting model. To promote the efficiency of identification, the joint mechanism must be operated without a load, in other words, $\tau _{e}$ is zero. The dynamic model can be further represented in Eq. (5) according to Eq. (3) and Eq. (4).

(5) \begin{equation}\tau =Y\left(\dot{q},\ddot{q}\right)P\end{equation}

where $Y(\dot{q},\ddot{q})$ is the measurement matrix, specifically $[\ddot{q},\mathrm{sgn}(\dot{q}),\dot{q}]$ , and $P$ is the parameter identification matrix, specifically $[{I_{zz}},{f_{c}},{f_{v}}]^{T}$ .

Then use Eq. (6) to identify parameter identification matrix.

(6) \begin{equation}P=\left[Y^{T}Y\right]^{-1}Y^{T}\tau\end{equation}

The parameter identification procedure uses the robot’s trajectory experiments to form the datasets. Hence, the excitation trajectory design would influence the accuracy of robot dynamic parameter identification [Reference Swevers, Ganseman, Tukel, d. Schutter and Brussel21]. The sampling datasets furtherly form the coefficient matrix $Y(\dot{q},\ddot{q})$ . While the condition number of coefficient matrix is large or ill-conditioned condition number, the dynamic model would be numerical sensitive, which means tiny input’s changes (robot state variables’ change) would affect the output’s accuracy (torque’s accuracy in Eq. (5). Therefore, the excitation trajectory should be designed with its coefficient matrix having a small condition number, specifically a well-conditioned condition number to reduce errors and uncertainties caused. Additionally, the excitation trajectory should be periodic and smooth for repeatable robot motion so that the Fourier series [Reference Park22], composed of sine and cosine functions with the periodic feature and continuous higher-order derivatives is designed to meet the above requirements. In this article, a third-order Fourier series, as shown in Eq. (7), is chosen as the excitation trajectory for dynamic parameters identification.

(7) \begin{align} q\!\left(t\right)=a+\sum _{n=1}^{3}\left[\frac{b_{n}}{\omega _{f}n}\sin\! \left(\omega _{f}nt\right)-\frac{c_{n}}{\omega _{f}n}\cos\! \left(\omega _{f}nt\right)\right] \end{align}

where $t$ is time, $a$ is a constant, $b_{n}$ and $c_{n}$ are the Fourier coefficients, $\omega _{f}$ is the angular frequency.

To ensure that the joint mechanism adheres to its hardware’s physical motion limitations while executing the trajectory, and robot’s smooth start or stop, the boundary condition of the excitation trajectory is shown in Eq. (8).

(8) \begin{align} \left\{\begin{array}{l} \dot{q}\left(t\right)\leq \dot{q}_{max}\\[3pt] q\!\left(t_{0}\right)=q\!\left(t_{f}\right)=0\\[3pt] \dot{q}\!\left(t_{0}\right)=\dot{q}\!\left(t_{f}\right)=0\\[3pt] \ddot{q}\!\left(t_{0}\right)=\ddot{q}\!\left(t_{f}\right)=0 \end{array}\right. \end{align}

where $\dot{q}_{max}$ is the upper limit of the joint angular velocity, $t_{0}$ is the initial time of the trajectory, and $t_{f}$ is the final time of the trajectory.

3.1. Data processing of physic-free model

To improve the data quantity of our external torque estimation BPNN model, data preprocessing is necessary to filter some error values such as outers and normalize data to promote the training performance. The specific filtering method is described in Section 4.2. Besides, the data preprocessing is performed using linear normalization as shown in Eq. (9) and its denormalization is shown in Eq. (10).

(9) \begin{equation}y^{\mathrm{*}}=\frac{y-y_{min}}{y_{max}-y_{min}}\end{equation}

where $y$ is the original data, $y_{max}$ is the maximum value of $y$ , $y_{min}$ is the minimum value of $y$ , and $y^{\mathrm{*}}$ is the preprocessed data.

(10) \begin{equation}y=\left(y_{max}-y_{min}\right)\mathrm{*}y^{\mathrm{*}}+y_{min}\end{equation}

3.2. BPNN structure in physic-free model

The BP neural network of physic-free model is made of an input layer, several hidden layers, and an output layer as Figure 5 shows. The output layer and input layer’s neurons are determined by the dimension of input data (three dimensions which are the velocity, acceleration, and current) and output data (one dimension which is the external torque) respectively.

Figure 5. BP neural network of Physic-Free Model.

Figure 6. The principle of neuron’s mechanism.

In the hidden layers, every neuron has its inputs and outputs. For a neuron, each of the inputs is summarized together by the linear combination and furtherly functioned by a nonlinear transformation to output its result as Eq. (11) shows. The principle of the $i^{th}$ neuron in the $m^{th}$ layer is designed in Figure 6.

(11) \begin{equation}x_{m}^{j}=f_{m}\left(\sum _{i=1}^{N_{m-1}}\omega _{m}^{i,j}x_{m-1}^{i}\right),j\in \left[1,2,\ldots, N_{m}\right]\end{equation}

where $x_{m}^{j}$ is the output of the ${neuron}_{m}^{j}$ , $f_{m}$ is the activation function of $m^{th}$ layer, which is always a nonlinear function. $\omega _{m}^{i,j}$ is the weight between ${neuron}_{m}^{j}$ and ${neuron}_{m-1}^{i}$ , and all weights of the neural network are shown in Eq. (12).

(12) \begin{equation}\omega =\left[\omega _{1},\cdots \omega _{m},\cdots, \omega _{M+1}\right],\omega _{m}=\left[\omega _{m}^{1,1},\cdots, \omega _{m}^{i,j},\cdots, \omega _{m}^{N_{m-1},N_{m}}\right]^{T}\end{equation}

While the neuron principle is determined, the number of hidden layers and the number of neurons in each layer need to be designed. The more complex the hidden layers, the stronger the network’s fitting ability, but excessively complex hidden layers may lead to overfitting, reducing the generalization ability of the network [Reference Karsoliya23]. The design of the hidden layers mostly relies on empirical knowledge and can refer to Eq. (13).

(13) \begin{equation}N_{m}=\sqrt{N_{in}+N_{out}}+d\end{equation}

where $N_{m}$ , $N_{in}$ , and $N_{out}$ represent the number of neurons in the hidden layers, input layer, and output layer, respectively, $d$ is a constant between 0 and 10. After determining the structure of a network, training is continuously needed to obtain the optimal parameters of the network.

3.3. BPNN training in physic-free model

To verify the model in Section 3.2, 85% of data sets with a size of $N_{tr}$ are used for network training, and the other part with a size of $N_{val}$ is used for verifying the network performance, which are the training set and the validation set, respectively. The training process of BPNN includes two parts: forward propagation and backward propagation, and the aim of the training is to optimize the initial weights until the optimal solution is solved continuously. In our model, all of data are processed before which are $[{I^{\mathrm{*}}},{\dot{q}^{\mathrm{*}}},{\ddot{q}^{\mathrm{*}}}]^{T}$ to ${\tau _{p}}^{\mathrm{*}}$ . For the forward propagation, each set of data inputs into the network, passes through the hidden layer in order, and finally goes through the output layer to obtain the estimated result ${\tau _{p}}^{\mathrm{*}}$ according to Eq. (11) and Eq. (14).

(14) \begin{equation}{\tau _{p}}^{\mathrm{*}}=f_{M+1}\left(\sum _{i=1}^{N_{M}}\omega _{M+1}^{i,1}x_{M}^{i}\right)\end{equation}

where ${\tau _{p}}^{\mathrm{*}}$ is the output of the output layer neuron. For the convenience of subsequent tasks, the mathematical operation $\Phi$ is determined for reshaping $\omega$ into a column vector $\Omega$ of length $n$ in Eq. (15).

(15) \begin{equation}\Omega =\Phi \left(\omega \right)=[{{\omega _{1}}^{T}},\cdots {{\omega _{m}}^{T}},\cdots, {{\omega _{M+1}}^{T}}]^{T}\end{equation}

The actual $\tau$ is processed by normalization to nominal $\tau ^{\mathrm{*}}$ in this model. Furtherly, an error model is designed in Eq. (16). To drive model training, the initial weights $\Omega ^{0}$ are randomly generated within the range of 0 to 1, and the loss function based on output error is designed in Eq. (17).

(16) \begin{equation}e_{r}\left(\Omega \right)=\tau ^{\mathrm{*}}\left(r\right)-{\tau _{p}}^{\mathrm{*}}\left(\Omega, x_{0}^{1}\left(r\right),x_{0}^{2}\left(r\right),x_{0}^{3}\left(r\right)\right)\end{equation}

(17) \begin{equation}L\left(\Omega \right)=\frac{1}{N_{tr}}\sum _{r=1}^{N_{tr}}\left[e_{r}\left(\Omega \right)\right]^{2}\end{equation}

where $e_{r}(\Omega )$ is the error between the actual value $\tau ^{\mathrm{*}}$ and the estimated value $\tau _{p}^{\mathrm{*}}$ of the $r^{th}$ sample, $L(\Omega )$ is the loss function, which uses the mean square error.

Using Levenberg-Marquardt algorithm [Reference Hagan and Menhaj24] for backpropagation to optimize loss problem, the optimal solution $\Omega ^{opt}$ can be calculated as Eq. (18) shows. There are some different optimization algorithms such as gradient descent, Newton-Euler, and LM which are shown in Eq. (19), Eq. (20), and Eq. (22) with the weights updating condition $L(\Omega ^{k+1})\lt L(\Omega ^{k})$ .

(18) \begin{equation}\Omega ^{opt}={argmin}\!\left\{L\left(\Omega \right)\right\}\end{equation}

(19) \begin{equation}\Omega ^{k+1}=\Omega ^{k}-pg\left(\Omega ^{k}\right)\end{equation}

(20) \begin{equation}\Omega ^{k+1}=\Omega ^{k}-H(\Omega ^{k})^{-1}g\left(\Omega ^{k}\right)\end{equation}

(21) \begin{equation}H\left(\Omega ^{k}\right)=\left[\begin{array}{c@{\quad}c@{\quad}c@{\quad}c} \dfrac{\partial ^{2}L\left(\Omega ^{k}\right)}{\partial (\Omega _{1}^{k})^{2}} & \dfrac{\partial ^{2}L\left(\Omega ^{k}\right)}{\partial \Omega _{1}^{k}\partial \Omega _{2}^{k}} & \cdots & \dfrac{\partial ^{2}L\left(\Omega ^{k}\right)}{\partial \Omega _{1}^{k}\partial \Omega _{n}^{k}}\\[3pt] \dfrac{\partial ^{2}L\left(\Omega ^{k}\right)}{\partial \Omega _{2}^{k}\partial \Omega _{1}^{k}} & \dfrac{\partial ^{2}L\left(\Omega ^{k}\right)}{\partial (\Omega _{2}^{k})^{2}} & \cdots & \dfrac{\partial ^{2}L\left(\Omega ^{k}\right)}{\partial \Omega _{2}^{k}\partial \Omega _{n}^{k}}\\[3pt] \vdots & \vdots & \ddots & \vdots \\[3pt] \dfrac{\partial ^{2}L\left(\Omega ^{k}\right)}{\partial \Omega _{n}^{k}\partial \Omega _{1}^{k}} & \dfrac{\partial ^{2}L\left(\Omega ^{k}\right)}{\partial \Omega _{n}^{k}\partial \Omega _{2}^{k}} & \cdots & \dfrac{\partial ^{2}L\left(\Omega ^{k}\right)}{\partial (\Omega _{n}^{k})^{2}} \end{array}\right]\end{equation}

where $\Omega ^{k+1}$ is the specific value of the parameters in the neural network after the iteration $(k+1)^{th}$ , $g(\Omega ^{k})$ is the search direction of $L(\Omega ^{k})$ at point $\Omega ^{k}$ , $p$ is the step size, $H(\Omega ^{k})$ is the Hessian matrix of $L(\Omega ^{k})$ at point $\Omega ^{k}$ .

Reasonable step size $p$ is vital for nonlinear problem convergence. Generally, a big $p$ could hinder convergence and vice versa [Reference Hagan, Demuth and Beale25]. For backpropagation, there are Newton’s method and the LM algorithm. Newton’s method improves convergence speed by using the Hessian matrix of the loss function to determine the step size. However, the calculation of the Hessian matrix requires computing the second derivative of the loss function, which incurs substantial computational costs. However, the Hessian matrix might not be positive definite, which could cause the algorithm’s instability. Therefore, Newton’s method is not typically used as the neural network backpropagation algorithm. The LM algorithm combines the advantages of gradient descent and Newton’s method, improving convergence speed and ensuring stability. Its core idea is to adjust the step size by approximating the Hessian matrix, which makes the LM algorithm have better search performance and quickly converge to the optimal solution. Therefore, this paper chooses the LM algorithm for neural network backpropagation as Eq. (22) to Eq. (25) show.

(22) \begin{equation}\Omega ^{k+1}=\Omega ^{k}-\left[J^{T}\left(\Omega ^{k}\right)J\left(\Omega ^{k}\right)+\mu U\right]^{-1}\left[J^{T}\left(\Omega ^{k}\right)e\left(\Omega ^{k}\right)\right]\end{equation}

(23) \begin{equation}J\left(\Omega ^{k}\right)=\left[\begin{array}{c@{\quad}c@{\quad}c@{\quad}c} \dfrac{\partial e_{1}\left(\Omega ^{k}\right)}{\partial \Omega _{1}^{k}} & \dfrac{\partial e_{1}\left(\Omega ^{k}\right)}{\partial \Omega _{2}^{k}} & \cdots & \dfrac{\partial e_{1}\left(\Omega ^{k}\right)}{\partial \Omega _{n}^{k}}\\[3pt] \dfrac{\partial e_{2}\left(\Omega ^{k}\right)}{\partial \Omega _{1}^{k}} & \dfrac{\partial e_{2}\left(\Omega ^{k}\right)}{\partial \Omega _{2}^{k}} & \cdots & \dfrac{\partial e_{2}\left(\Omega ^{k}\right)}{\partial \Omega _{n}^{k}}\\[3pt] \vdots & \vdots & \ddots & \vdots \\ \dfrac{\partial e_{{N_{tr}}}\left(\Omega ^{k}\right)}{\partial \Omega _{1}^{k}} & \dfrac{\partial e_{{N_{tr}}}\left(\Omega ^{k}\right)}{\partial \Omega _{2}^{k}} & \cdots & \dfrac{\partial e_{{N_{tr}}}\left(\Omega ^{k}\right)}{\partial \Omega _{n}^{k}} \end{array}\right]\end{equation}

(24) \begin{equation}\mu =\mu \cdot \beta\end{equation}

(25) \begin{equation}\mu =\frac{\mu }{\beta }\end{equation}

where $J(\Omega ^{k})$ is the Jacobian matrix, $e(\Omega ^{k})$ is an $n$ -dimensional column vector composed of whole $e_{r}(\Omega ^{k})$ , $\mu$ (MU) is a damping factor, which is divided by the constant $\beta$ when the performance of the neural network increases, otherwise $\mu$ will be multiplied by $\beta$ , $U$ is an identity matrix, adding $\mu U$ is to ensure that the matrix is positively definite.

Neural networks use forward propagation and backward propagation to adjust their parameters. In practical training, it is essential to establish stopping criteria by consideration of iterations, efficiency, synchronicity, overfitting, and so on. In engineering applications, there are five criteria as follows.

1. 1. $k\geq k_{{target}}$

2. 2. $L(\Omega ^{k})\leq L_{{target}}$

3. 3. $| g(\Omega ^{k})| \leq g_{{target}}$

4. 4. $\xi \geq \xi _{{target}}$

5. 5. $\mu \geq \mu _{{target}}$

where $k$ is the iteration, setting the maximum iteration $k_{{target}}$ to prevent training from continuing indefinitely. The neural network training approaches the optimum when $L(\Omega ^{k})$ achieves sufficiently small, the gradient magnitude $| g(\Omega ^{k})|$ achieves sufficiently small, or $\mu$ achieves sufficiently large. Additionally, $\xi$ (VAL CHECK) is the number of consecutive increases in validation set error, which is used for preventing overfitting for the neural network. When $\xi$ continues to rise $\xi _{{target}}$ times, it indicates a decline in the performance of the neural network, so the training should be stopped.

4.1. Experimental platform

A joint’s external torque testing platform is furnished as shown in Figure 7 to compare the physics-based model and physic-free model, which mainly consists of three parts: joint module, torque sensor, and powder brake. The joint module consists of a motor, harmonic reducer, and hardware to drive the actuator motion and monitor state variables, and the powder brake supplies a stable payload to simulate actual external torque, and the torque sensor records the value of different real-time external torque. Specifically, the experimental platform consists of six devices and its information as shown in Table I, in which there is a joint module (integrates a harmonic reducer, frameless torque motor, incremental encoder, and servo driver), Beckhoff controller and EL3064 (drives the actuator and receives the external torque information from the sensor), NCTE torque sensors, the powder brake and its tension controller (supplies the stable external torque) and power supply.

Table I. Characteristics of experimental platform.

Figure 7. Experiment platform.

4.2. Experimental principle

Figure 8 represents the experimental principle: the PC sends control commands following the excitation trajectory to the Beckhoff controller to drive the joint module motion. Then powder brake will work to generate different external torques by changing its control voltage signal. Meanwhile, the corresponding external torque value can be measured by the torque sensor, which is the analog voltage signal that can be converted into a digital signal through EL3064 and then transmitted into an actual physical quantity by Eq. (26). The joint data $I,\dot{q},\ddot{q}$ in Beckhoff controller and external torque $\tau _{e}$ are collected at a sampling frequency of $1000\,Hz$ into the PC simultaneously. Subsequently, design a second-order low-pass Butterworth filter [Reference Proakis26] with a cutoff frequency of $10\,Hz$ , whose frequency response characteristics in the z-domain can be represented by the transfer function in Eq. (27), where the coefficients are obtained using the MATLAB “butter” function [27]. Then, apply this transfer function to filter those data both in forwards and backwards procedure. The filtered data can be further used for experimental result analysis.

(26) \begin{equation}\tau _{e}=\frac{\left(\frac{2.5}{2^{13}}\mathrm{*}V_{out,d}-2.5\right)\mathrm{*}1000}{26.38}\end{equation}

(27) \begin{equation}H\left(z\right)=\frac{9.4469\times 10^{-4}+0.0019z^{-1}+9.4469\times 10^{-4}z^{-2}}{1-1.9112z^{-1}+0.9150z^{-2}}\end{equation}

Table II. Excitation trajectory parameters.

5.1. Experimental results of parameter identification

Before comparing the physics-based model and the physic-free model of torque estimation, the physical model method as well as physics-based model should be tested to show its performance. Hence, excitation trajectories are designed to identify the parameters by the MOB, and the coefficients and position figure of the identification trajectory and the validation trajectory are shown in Table II and Figure 9, respectively, which are with well-conditioned condition numbers and joint velocity limits $\dot{q}_{max}$ of $\frac{\pi }{3}/s$ and $\frac{\pi }{6}/s$ according to the Eq. (7) and Eq. (8).

Table III. Identified dynamic parameters of external torque estimation based on physical model.

Figure 8. Experimental schematic diagram.

Table IV. The error evaluation of different trajectories by identified dynamic model.

Figure 9. Excitation trajectory with the parameters in Table II.

Using the identification trajectory, and collecting data sets, the identified dynamic parameters can be calculated by Eq. (6) and the results are shown in Table III.

Then these identified dynamic parameters can construct the model and use the identification trajectory and validation trajectory to test its accuracy with the Eq. (2). The results are shown in Figure 10. From Figure 10, it’s obvious that the estimation accuracy is better in the identification trajectory compared to the validation trajectory. The reason is that the parameter identification is based on the identification trajectory. Moreover, the fitting procedure shows that drive torque estimation by the identified dynamic model has a bad generalization ability, proving that some lightweight robots’ drag-and-drop teaching sometimes behaves badly when the operators drag different joints with different forces.

Figure 10. Joint driving torque estimation in different excitation trajectory. (a) Identification trajectory and (b) validation trajectory.

Furthermore, the errors in different excitation trajectories are compared in Figure 11. The results show that the huge impulse errors always exist around the change of velocity direction, in which the velocity is around zero, as well as around the orange point zone in Figure 11. Table IV shows different evaluation way with the RSME and R-squared, which reveals that the identified dynamic model in identification trajectory has better performance with smaller RSME value and larger R-squared value.

Table V. The parameters of external torque estimation based on physic-free model.

Table VI. The error evaluation by MOB model and BPNN model with the same identification trajectory.

Figure 11. Estimation errors of different trajectory by identified dynamic model.

5.2. Comparison of external torque estimation with two model

To improve the model generalization ability by the BPNN method, it is crucial to establish a diverse training set, collected by many different scenarios and variations. In our experiments, there are six control groups with different voltages ( $1V$ , $2V$ , $3V$ , $4V$ , $5V$ , $6V$ ) of powder brake, corresponding to different external torque for the joint module from small to large. Then select the identification trajectory in Figure 9 as the excitation trajectory and use the control groups with $1V$ , $2V$ , $3V 5V$ , and $6V$ as training sets with $4V$ as a testing set to build up the BPNN model. Through the training, the model parameters are shown in Table V.

Compare physics-based model (MOB model) and physic-free model (BPNN model) with the identification trajectory with some payload ( $4V$ of powder brake) to verify the accuracy. In this MOB model, the proportion coefficient $k_{p}$ is 10. The estimated torque by different model is shown in Figure 12, in which the BPNN model has better fitting ability compared to the MOB model.

Figure 12. Estimated torque by MOB model and BPNN model with same identification trajectory.

Furthermore, the error of different model is shown in Figure 13, which also shows that the huge impulse errors are around the change of velocity direction (orange zone in Figure 13). Table VI shows different evaluation way with the RSME and R-squared, which reveals that the BPNN model has better performance with smaller RSME value and larger R-squared value.

Figure 13. Estimated torque error by MOB model and BPNN model with same identification trajectory.