Towards Model-Free Tool Dynamic Identification and Calibration Using Multi-Layer Neural Network

In robot control with physical interaction, like robot-assisted surgery and bilateral teleoperation, the availability of reliable interaction force information has proved to be capable of increasing the control precision and of dealing with the surrounding complex environments. Usually, force sensors are mounted between the end effector of the robot manipulator and the tool for measuring the interaction forces on the tooltip. In this case, the force acquired from the force sensor includes not only the interaction force but also the gravity force of the tool. Hence the tool dynamic identification is required for accurate dynamic simulation and model-based control. Although model-based techniques have already been widely used in traditional robotic arms control, their accuracy is limited due to the lack of specific dynamic models. This work proposes a model-free technique for dynamic identification using multi-layer neural networks (MNN). It utilizes two types of MNN architectures based on both feed-forward networks (FF-MNN) and cascade-forward networks (CF-MNN) to model the tool dynamics. Compared with the model-based technique, i.e., curve fitting (CF), the accuracy of the tool identification is improved. After the identification and calibration, a further demonstration of bilateral teleoperation is presented using a serial robot (LWR4+, KUKA, Germany) and a haptic manipulator (SIGMA 7, Force Dimension, Switzerland). Results demonstrate the promising performance of the model-free tool identification technique using MNN, improving the results provided by model-based methods.


Introduction
Robot control with physical interaction has been widespread and draws a lot of research interests in the past decades [1]. The availability of reliable interaction force information proved to be capable of increasing the control precision and of dealing with the surrounding complex environments, for example, in the context of robot-assisted surgery and bilateral teleoperation [2]. In practical applications, since the force sensor is used for measuring the interaction forces on the tooltip, the force sensor is usually mounted between the end-effector of the robot and its tool. However, the force acquired from the force sensor includes not only the interaction forces but also the gravity force of the tool. Hence the tool dynamic identification is required for accurate dynamic simulation and model-based control. Especially for bilateral teleoperation control, which provides haptic feedback for the surgeon in robot-assisted minimally invasive surgery (RA-MIS) [3,4], achieving accurate

Kinematic Model of the Serial Robot
The kinematic model of an anthropomorphic seven degrees-of-freedom (DoFs) robotic arm (LWR4+, KUKA, Germany) is shown in Figure 1. The corresponding Denavit-Hartenberg (D-H) parameters [26] are listed in Table 1 [27]. Based on the D-H parameters, a homogeneous transformation matrix of two consecutive link frames of the serial robot arm, i − 1 to i can be defined as follows [28]: where the transformation matrix i i−1 T is a composition of rotations and translations to move a frame coincident to frame i − 1 until it coincides with frame i [29]. Parameters of link i − 1 include the link twist angle α i−1 , link length a i−1 , and link offset d i , whereas parameters of link i include joint variable θ i . Rot and Tr are (4×4) matrices of rotational and translational transformations along an axis, respectively [28]. Therefore, the rotation angle can be obtained from its forward kinematic function.
With the D-H matrix at hand, the transformation matrix from the robot base frame to its end effector frame can be computed using joint angles. The robot tool pose can be obtained by multiplying the link transformation matrix as follows [28]: where i i+1 T is the transformation matrix as shown in (1).

Methodology
As shown in Figure 2, in order to transmit the interaction force on the tooltip to the robot, the force sensor is mounted between the end effector and the tool. The interaction force between the tooltip and the environment is measured. It is obvious that the output of the force sensor includes not only the interaction force but also force generated by load of the tool gravity. Hence, it is essential to develop techniques for tool dynamic identification to eliminate the disturbance from the tool weight and transmit the accurate interaction force to the robot control system.

Tool Dynamics Identification
Due to the influence of the gravitational force of the robot tool, the output of force sensor F S ∈ R 3 comprises the gravitational force F ToolGravity ∈ R 3 and the interaction force F Interaction ∈ R 3 between the tool and the environment. The formula can be written as: With the motion of the robot, the components of the gravity force F ToolGravity on F S vary according to the pose of the robot tool. As shown in Figure 3, the force exerted on the force sensor varies with the tool directions θ 1 and θ 2 . Hence, the output of the force sensor is perturbed and it cannot represent the interaction force accurately. It is essential to identify the tool dynamics and eliminate the influence of the gravity for accurate force sensing. 3.1.1. Model-Based Tool Gravity Identification Using Curve Fitting According to Figure 3, the orientation matrix of the tool changes with the robot arm placement, and the impact of the weight on the force sensor measurements also varies with θ 1 and θ 2 . Thus, it is essential to consider the gravity identification depending on the orientation of the tool.
Considering the influence of the tool gravity force on the force sensor output, the mathematical model of the force sensor output (Figure 3) can be defined as following by using the Euler angles: where F x , F y and F z are the outputs of the force sensor. The unknown parameters are the mass m, and the constant coefficients a, b, and c. g represents gravity which is 9.8 m/s 2 . The angles θ 1 and θ 2 are the orientation angles of the tool pose. d is the deviation error angle around the z-axis from the tool installation.
If there is no interaction force on the tool, the output of the force sensor represents the effects of the tool gravity, as follows: With the acquired data including tool pose and the output of the force sensor without interaction force, the parameters listed in Equation (4) can be obtained with the CF technique. The detailed results of CF can be found in our previous works [30].

Model-Free Tool Gravity Identification Using Mnn
However, it is difficult to identify an accurate mathematical model due to manufacturing and assembly variances. Moreover, there should be a deviation error on θ 1 because the tool cannot be installed in a straight way. Hence it is difficult to project the gravity force of the tool on the force sensor with the model proposed in (4). The precise mapping relation between the gravity and the rotation angles of the tool direction is, therefore, too complex to model. Hence, we propose a novel model-free based method to map the relation between the gravity force and the rotation angles as follows: where f is a nonlinear unknown function and the corresponding Euler angles θ x , θ y , θ z are the multiple inputs of the function.
In the past decades, the ANN approach became the most popular method for modeling linear and nonlinear regression problems [31]. Although the capability of ANN models extends to establishing any complex function and nonlinear relationship between a certain set of inputs and outputs with multiple dimensions [32], some limitations remain and need to be solved, such as over-fitting, under-fitting and time-consuming.
To enhance the predictive performance of the built model, such as high-accuracy, stability and high-speed, this work adopts two types of MNN methods to establish a regression model, namely FF-MNN and CF-MNN. Figure 4 shows the structure of MNN mapping the 3D inputs (the degree of angles [θ x , θ y , θ z ]) and 3D outputs (the forces [F x , F y , F z ]). In contrast to the FF-MNN model which connects all neurons in each hidden layer for fitting the nonlinear function, the CF-MNN [33,34] model can acquire the information from all of the previous layers by connecting the outputs to the latter networks (shown in the blue lines).
Hence, the MNN-based nonlinear model for mapping the time-varying multiple inputs x t = [θ x , θ y , θ z ] t to the multi-outputs y t = [F x , F y , F z ] t can be defined as follows: where the whole parameter set Θ accounts for all of the nonlinear weights matrix ω i,j k , i, j k ∈ R + and bias b j k in each layer where k is the order of layer. By substituting the related parameters into FF-MNN model, the nonlinear function can be written as: Φ k is the activation function. In this article, we chose the Broyden-Fletcher-Goldfarb-Shanno (BFGS) quasi-newton function [35]. ω o and b o are the parameters of output layer. γ j k and b j k are the outputs of jth neuron. Similarly, the CF-MNN model can be expressed as follows: It uses not only the output of the former layer but reserves also all of the previous information for obtaining the final results. The MNN regression model aims to search the optimal parameter set Θ by computing the minimum least squares between the predicted resultŷ and the real value y as follows: There are three common evaluation indices to measure the performance of the built MNN models, namely mean square error (MSE), root mean square error (RMSE) and Pearson correlation coefficient ρ defined in Equation (11). .
The time t can be regarded as the number of observations. µŷ and σŷ are the average and standard deviation ofŷ, while µ y and σ y are the same values of y. The best score for the Pearson correlation coefficient ρ is 1 while for the other errors it is 0. The MNN model aims to predict the force close to the measured value.

Force Sensor Calibration
The force sensor is a particularly significant source of feedback in robotic applications to measure forces along x-, yand z-axes at the end effector of the robot increasing sensitivity of the surgeon [36]. To achieve the best possible transparency, the force sensor should be calibrated in the system where it will be used. The SVD of a matrix is a linear algebra tool that has been successfully applied to a wide variety of domains [37]. In this work, the SVD method [38] is adopted to figure out the transformation (calibration) matrix f e T between the reference frames of the slave's end effector and the force sensor, as depicted in Figure 5. Figure 5 demonstrates the input-output of our calibration method, where F R ∈ R 1×3 , F S ∈ R 1×3 , are the robot and sensor forces, respectively. f e T ∈ R 4×4 is the calibration matrix.
where F h ∈ R 1×3 is the final haptic force on the end effector of the robot frame. After the identification and calibration, the output of the system is in the end effector frame, which is able to achieve accurate force sensing for the robot control. The overview of the procedure of the tool dynamic identification and calibration using MNN is shown in Figure 6.

System Description
A brief description of the robot system developed in this project is shown in Figure 7. A redundant robot (LWR4+, KUKA, Augsburg, Germany) served as the slave manipulator torque-controlled through fast research interface (FRI), which could provide direct low-level real-time access to the robot controller (KRC) at rates of up to 1 kHz [39]. The software system was developed with OROCOS (Open Robotic Control Software, http://www.orocos.org/) with a real-time Xenomai-patched Linux kernel and ROS (Robot Operating System, http://www.ros.org/) kinetic in Ubuntu [40]. To guarantee control frequency, force sensor, ROS node and OROCOS torque controller were executed on separate computers with UDP communication [41] between each other: the control loop and the sensing ROS node was executed on separate computers, as shown in Figure 8. The system consists of: • a seven DoFs LightWeight robotic arm (LWR4+, KUKA, Augsburg, Germany) as slave device. • a six-axis force sensor (M8128C6, SRI, Nanning, China) [42] that has the purpose of measuring interaction force between the surgical tool-tip and the environment.  Overview of the developed software system. The "LWR4+ Impedance Controller" is adopted to allow hands-on control to move the surgical robot arm by hands. The "force sensor" is developed using Robot Operating System (ROS) and it communicate with the Open Robotic Control Software (OROCOS) by ROS topic with a frequency of 500 Hz.

Tool Dynamic Identification
Firstly, hands-on control was activated to allow the user to move the robot arm without touching the robot tool and the force sensor. In this way, two groups of data were collected for estimation and validation. All collected signals had 74,195 samples which are divided into a training dataset (42,148 samples) for building the regression models and a testing dataset (32,047 samples) for evaluating the performance of the built models.

Model-Based Tool Dynamic Identification Using Cf
Firstly, we introduced CF for tool dynamic identification. As mentioned before, dynamic identification was implemented with respect to current end effector orientation. By utilizing the CF technique, the constant unknown parameter m, and the constant coefficients a, b and c can be obtained with the first group of sampled data (42,148 samples). Then, the obtained parameters were placed in the mathematical model to predict the force on the force sensor, which is expressed as follows:      F x,estimated = −0.3434 * g * sin (θ 1 ) * cos (θ 2 + 1.401)+0.6 F y,estimated = −0.3434 * g * sin (θ 1 ) * sin (θ 2 + 1.401) + 1.1 F z,estimated = 0.3434 * g * cos (θ 1 ) + 2.0 (13) After the tool dynamic identification with CF, we validated the obtained mathematical model with the testing data (32,047 samples). The RMSE on the x-axis was calculated as 1.696, while it was 2.931 on the y-axis, and 1.057 was obtained on the z-axis. The overall RMSE of the testing data of the prediction error of the norm of the force is 5.684.

Model-Free Tool Dynamic Identification Using Mnn
The model-based tool identification method has shown a big error due to the inaccurate mathematical model. To solve this problem, we adopted MNN to model the tool dynamic without using a mathematical model. To implement a nonlinear regression model for enhanced accuracy, fast computation and strong stability, the experiments were designed to compare and discuss the performance among eight single and multiple layer NN models with different numbers of nodes. As shown in Table 2, the first four MNN models (M1-M4) have two hidden layers with different numbers of nodes. For example, M1 had 30 and 15 neurons in the first and second hidden layers, respectively. Hence, we adopted the notation [30,15] to represent the nodes. The models M1 and M2 use CF networks, while M3 and M4 adopt FF networks. To enhance the validity of the comparison results, four single hidden layer NN models (M5-M8) based on CF and FF networks, namely CF-single layer neural network (SNN) and FF-SNN, are chosen in the experiments. The four aspects for the evaluation the performance of the eight models comprise the training time, the online predictive time, the regression errors (i.e., MSE and RMSE), and the correlation coefficient ρ.
To avoid the phenomenon of over-fitting and under-fitting and to find the most stable NN, each of the models was tested 30 times.
The first experiment was designed for proving the reconstruction ability of the adopted models. Figure 9 shows the results of comparing MSE, RMSE and correlation coefficient ρ among models M1 to M8 on the training dataset. The top three rows show the MSE, RMSE and ρ of each channel (x, y and z). The fourth row displays the sum of errors (MSE and RMSE) and the average of coefficients ρ. By observing the eight pictures on the left two columns, M1 and M3 obtain the lowest errors for modeling the tool dynamic, while M1 had a lower error on the data collected from channel x than model M3. In addition, M1 was the most robust model due to its smaller standard deviation in respect to the other models. Similarly, by comparing the results of correlation coefficient ρ (in the third column), M1 proved to be the best reconstruction model for mapping the Euler angles to the force on the training dataset.
The errors and coefficient only illustrate the ability of reconstruction accuracy and model stability. However, the time-consuming of training a NN model is an other significant issue which is needed to discuss. Table 3 Figure 10 shows the measured and predicted force (by the M1 model) of the three channels. By observing the trend and difference between the measured curves and the predicted curves, the M1 model almost completely reconstructed the shape of the original curves.
After acquiring the regression models, their performance has to be validated on the new collected (testing) dataset. Similarly to the above, the errors (MSE and RMSE), correlation coefficient ρ and testing time were the important indices for measuring the eight models. Figure 11 shows the comparison results of MSE, RMSE and ρ among M1 to M8. Contrary to the results obtained on the training dataset, the M1 and M3 got the worst errors and ρ values with respect to the other models. On the other hand, M2 and M4 are proven to be the best models for predicting force on the testing dataset. Because both of them acquired lower errors and higher ρ values than the other models.
The M1 and M3 approach obtain the worst errors with respect to the other models, but the difference between M1, M3 and other models is not too much. For example, the overall MSE of M1 model is about 0.015, while the best results computed by M8 is 0.018. The M1 and M3 had overfitting problems. The M2 and M4 models acquired better performance on the testing dataset because they are fit the outputs on that dataset. However, when they were adopted on the training dataset, they were underfitting. Finally, the M1 and M3 models were not suitable for predicting on the testing dataset.  The average testing time was a significant index to evaluate detection speed of the built model. Table 4 displays the average and standard deviation of the testing time. Although the M2 model can acquire the highest regression accuracy, it spends 0.0230 s to predict a result which is slower than the other methods. However, this predictive speed was sufficient for the real experiment.  Figure 12 shows the force predicted by the M2 model for each channel. By contrasting the difference between the measured and predicted force, the M2 model can track the real force.

Force Sensor Calibration
Although the interaction force is accurate after tool dynamic identification, the robot does not know the placement of the force sensor, especially for the orientation on its z axis. In this way, the x and y axes were not aligned with the robot coordination frame. Hence a calibration procedure was required to transfer the interaction force into the robot coordinates. As mentioned above, after identifying tool dynamic component by MNN, we performed the force sensor calibration by collecting more data using hands to touch the force sensor and then applying the SVD to solve the transformation matrix. It should be noticed that the influences from the gravity of the force sensor have been eliminated with the tool gravity identification using neural networks together. For the force on the robot end-effector, we used the software package provided by the KUKA Fri interface [40], and we assumed that there were no disturbances from the robot dynamic effects on it. As it is shown in Figure 13, the user applied hand force only on the tool. In this way, the force sensor and the robot were under the same external force. We collect the force from robot [39] and the force from the force sensor. To facilitate the availability of the force from the force sensor for the robot, the transformation matrix between them was calculated using an SVD based method. The results of the calibration are depicted in Figure 14.

Conclusions and Future Work
This paper presents a novel model-free based tool dynamic identification using MNN and force sensor calibration for accurate force sensing. Firstly, tool gravity force was identified by CF and MNN methods. The results showed that a model-free based tool dynamic identification using MNN is more accurate than a model-based identification using CF. Afterwards, force sensor calibration was implemented using SVD. Results showed that the calibrated force was able to represent the force from physical interaction.
Furthermore, the first CF-MNN model (M1) is proven as the best method for reconstructing the training dataset, while the first FF-MNN model (M2) is the best one to predict force on the testing dataset. It is observed that MNN approximation is more accurate than CF in the estimation of the gravity force, thus enhancing the accuracy of the tool dynamic identification and calibration for force sensing. It has been well known that Deep Neural Networks, such as convolutional neural networks (CNNs) and recurrent neural networks (RNN) [43], are capable of learning and modeling complex system [44] with high efficiency. Hence, further works will try to introduce DNN to model the gravity forces and compare their performances with the proposed methodology in this paper. The number of nodes is chosen through trial-and-error in this paper. We put the further analysis of how to select the neurons number as future works. An additional force sensor will also be placed to measure and validate the interaction force after processing. Future works will consider the robot dynamic effects on the force of its end-effector and achieve higher accuracy for the calibration procedure. Future work will also consider more challenging problems (e.g., dead-zone and time-delay) in our robot control framework. The system stability [45][46][47] and tracking accuracy [48,49] might not be guaranteed under these situations, which are a precondition for safety in robot control.

Conflicts of Interest:
The authors declare no conflict of interest. The founders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

Abbreviations
The following abbreviations are used in this manuscript: