An Integrated Compensation Method for the Force Disturbance of a Six-Axis Force Sensor in Complex Manufacturing Scenarios

As one of the key components for active compliance control and human–robot collaboration, a six-axis force sensor is often used for a robot to obtain contact forces. However, a significant problem is the distortion between the contact forces and the data conveyed by the six-axis force sensor because of its zero drift, system error, and gravity of robot end-effector. To eliminate the above disturbances, an integrated compensation method is proposed, which uses a deep learning network and the least squares method to realize the zero-point prediction and tool load identification, respectively. After that, the proposed method can automatically complete compensation for the six-axis force sensor in complex manufacturing scenarios. Additionally, the experimental results demonstrate that the proposed method can provide effective and robust compensation for force disturbance and achieve high measurement accuracy.


Introduction
Nowadays, robots have become indispensable equipment in industrial manufacturing systems [1]. With the urgent demand for flexible manufacturing and fast development of sensor technology, higher requirements are put forward for the intelligence and selfadaptability of robots [2]. Two prominent issues are robot-environment interaction [3] and human-robot collaboration [4,5]. The development of a multi-axis force sensor, especially a six-axis force sensor, makes it possible to obtain force information in robot operation environments. Additionally, six-axis force sensors have been applied to robotic systems that perform contact tasks, such as grinding robots, surgical robots, etc. [6,7].
In order to detect six-axis forces and moments in space, a cylindrical six-axis force sensor [8,9] is often used, which can measure three-axis orthogonal forces and moments relative to the sensor frame, i.e., the force components Fx, Fy, Fz along the X, Y, Z axes of the sensor and the moments Mx, My, Mz around the X, Y, Z axes of the sensor. When external forces are applied to the sensor, the internal elastomers are deformed and the corresponding strain signal according to its electrical characteristics are output [10][11][12]. Typically, a six-axis force sensor is mounted between the end flange of an industrial robot and the end-effector to perceive the external forces [13,14]. In this case, the sensor generally measures a combination of forces including the contact forces, the gravity force acting on the tool, self-gravity, system error, and other force disturbances. The tool gravity and self-gravity continuously affect measurement output as the robot pose changes. The selfgravity refers to the gravitational force acting on the six-axis force sensor. To accurately obtain the contact forces applied to the six-axis force sensor, it is necessary to compensate for force disturbance.
For several decades, researchers have been conducting wide studies to compensate for force disturbances of a six-axis force sensor at the end of a robot. Shetty B R et al. [15][16][17]. showed that the output of the sensor consisted of load gravity and interaction force, and where f bias is the drift of the sensor, f sys is the system error, f o_gravity is the effect caused by the self-gravity of the sensor, f t_gravity is the effect caused by the gravity force acting on the tool, f inertia is the force output due to robot vibration and inertia, f contact is the pure contact force applied to the end-effector. In addition to the pure contact force, f bias , f sys , and f o_gravity are parameters related to the sensor itself, f t_gravity relies on the selected execution tool. Since the robot is in low speed, the effect of f inertia is not considered in this paper. Let the effects caused by the sensor itself be attributed to f zero . Therefore, Equation (1) is rewritten as: where f zero = f bias + f o_gravity + f sys is the output of the sensor without any external force, which is defined as the zero-point. The goal of the method is to calculate the zero-point of the sensor for different poses of the robot and the effect of the tool gravity on the output. Based on the above calculation results, the effect of force interference factors on the contact force perception is eliminated.

Integrated Compensation Method
In view of the effects of drift and system error, as well as the interference of the force sensor self-gravity and end-effector gravity, we need to compensate for the actual output of the force sensor in order to accurately perceive the contact forces. Therefore, we separate the force sensor and the end-effector, use deep learning to obtain the mapping relationship between the attitude and the zero-point of the force sensor, and use the least squares method to estimate the tool load. Then, the compensation of the force sensor is completed, and the influence of force interference is eliminated.

Zero-Point Estimation of a Six-Axis Force Sensor Based on Deep Learning
The defined zero-point includes the drift, system error, and the influence of its selfgravity. Generally, for simplicity, the drift is often regarded as a constant, and a linear function is established based on the relationship between the self-gravity and the posture. However, in real-life application scenarios of robots, the established linear relationship does not fully reflect the relationship between the posture and output as well as the effects of drift and system error due to the installation and positioning of force sensors, etc. Deep learning is a machine learning technique that has made significant progress in the past decade. With sufficient training data and suitable network structures, neural networks (NNs) can approximate arbitrary nonlinear functions, and this powerful fitting ability has been widely applied to natural language processing and image recognition [26][27][28][29]. Studies have shown that neural networks have a better effect on error compensation [30,31]. Given its capability described above, NN is used in this section to estimate the zero-point of the force sensor.
Thus, it is necessary to establish a mapping between the pose of the force sensor and the zero-point in the absence of any load on the force sensor, i.e., where θ is the model parameter.
The pose of a rigid body can generally be described by rotation vector, rotation matrix, quaternion or Euler angles, etc. To reduce the computational complexity, the literature [15,16,22,23] chose a rotation matrix to derive the relationship between the effect of the self-gravity and the rotation matrix in the ideal case. According to the existing relationship, this paper designs an NN as shown in Figure 1 to estimate the zero-point.
The input to the NN is the nine elements of the rotation matrix R ij (i = 1, 2, 3; j = 1,2,3), and the general robot control system can directly obtain the data related to the pose and thus the rotation matrix can be calculated. Considering the symmetry of the force sensor, the theoretical output of Mz should be 0 in the non-load condition. Additionally, its output is small and fluctuates irregularly during the experiment. To reduce its unreliable effect on model training, the output data of the NN are the five outputs Fx, Fy, Fz, Mx, My of the force sensor, and the fluctuating output of Mz caused by incidental factors is not estimated. Theoretically, the accuracy of approximation can be improved by increasing the number of neurons by a moderate amount [32]. However, too many neurons require higher computation power, which inevitably affects the response time of the sensor in practice. Therefore, with a combination of fitting accuracy and sensor response speed, we chose a neural network with two hidden layers based on the experiments, the first with 45 nodes and the second with 90, as the regressor to fit the mapping function (p; θ). Currently, ReLU(x) = max(0, x) is often recommended as the activation function of NN, but in the force sensor application scenario, this activation function may cause some neurons to never be activated, resulting in the parameters not being updated. Therefore, this paper chooses the derivative of this function Leaky_ReLU(x) = max(αx, x) as the activation function, where α = 0.01. The activation function of the output layer uses the constant function. The parameter θ in the model is learned by minimizing the following loss function: , and the general robot control system can directly obtain the data related to the pose and thus the rotation matrix can be calculated. Considering the symmetry of the force sensor, the theoretical output of Mz should be 0 in the non-load condition. Additionally, its output is small and fluctuates irregularly during the experiment. To reduce its unreliable effect on model training, the output data of the NN are the five outputs , , , , Fx Fy Fz Mx My of the force sensor, and the fluctuating output of Mz caused by incidental factors is not estimated. Theoretically, the accuracy of approximation can be improved by increasing the number of neurons by a moderate amount [32]. However, too many neurons require higher computation power, which inevitably affects the response time of the sensor in practice. Therefore, with a combination of fitting accuracy and sensor response speed, we chose a neural network with two hidden layers based on the experiments, the first with 45 nodes and the second with 90, as the regressor to fit the mapping function ( ; ) is often recommended as the activation function of NN, but in the force sensor application scenario, this activation function may cause some neurons to never be activated, resulting in the parameters not being updated. Therefore, this paper chooses the derivative of this function Leaky_ReLU( ) max( , ) x

Tool Load Identification Based on Least Squares Method
Generally, different end-effectors are attached to the end of the force sensor in order to achieve various tasks. As the pose of the robot changes, the pose of the tool changes accordingly, and its influence on the sensor changes as well. In order to achieve a quick compensation of the gravitational influence of the end-effector after its replacement, this section uses the least squares method to estimate the tool load.
Let the robot base frame be { } B and the world frame be { } W , and assume that w Z is reversed with respect to gravity. Since the error in the robot installation will make { } B Figure 1. The neural network for zero-point estimation of a six-axis force sensor. A [0] is the input layer; A [1] is the first hidden layer; A [2] is the second hidden layer; A [3] is the output layer.

Tool Load Identification Based on Least Squares Method
Generally, different end-effectors are attached to the end of the force sensor in order to achieve various tasks. As the pose of the robot changes, the pose of the tool changes accordingly, and its influence on the sensor changes as well. In order to achieve a quick compensation of the gravitational influence of the end-effector after its replacement, this section uses the least squares method to estimate the tool load.
Let the robot base frame be {B} and the world frame be {W}, and assume that Z w is reversed with respect to gravity. Since the error in the robot installation will make {B} and {W} not completely coincide, assuming that {B} can be obtained by rotating {W} around the X w by angle α and around the Y w by angle β, then, the rotation matrix of the robot base frame relative to the world frame can be written as: The frame of the force sensor is {S}, and the relationship between the frames is shown in Figure 2.
The frame of the force sensor is { } S , and the relationship between the frames i shown in Figure 2. Assuming that the gravity force acting on the tool is G, . Through force transformation, the tool gravity in the ith attitude can be converted from the world frame to the force sensor frame: which can be equated to: , and B S i R is th rotation matrix of the sensor frame at the ith pose with respect to the robot base frame. The gravity force acting on the tool in the force sensor frame is shown in Figure 3. Assuming that the gravity force acting on the tool is G, then it is expressed in {W} as W G = 0 0 −G T . Through force transformation, the tool gravity in the ith attitude can be converted from the world frame to the force sensor frame: which can be equated to: where S G i = [G xi G yi G zi ] T , B G = G cos α sin β −G sin α −G cos α cos β T , and B S R i is the rotation matrix of the sensor frame at the ith pose with respect to the robot base frame.
The gravity force acting on the tool in the force sensor frame is shown in Figure 3. The frame of the six-axis force sensor is a spatial Cartesian rectangular frame. As suming the coordinate of the tool's centre of gravity in the sensor frame as ( , , ) x y z , ac cording to the relationship between forces and moments, it is obtained that: The frame of the six-axis force sensor is a spatial Cartesian rectangular frame. Assuming the coordinate of the tool's centre of gravity in the sensor frame as (x, y, z), according to the relationship between forces and moments, it is obtained that: Stacking all forces S G i as S G, moments S M i as S M, rotation matrix B S R T i as R, and matrix A i as A in n robot poses, we have: Therefore, the gravity force acting on the tool and centre of gravity can be estimated as: where A * is the matrix calculated from the optimal value of the tool gravity. Therefore, according to the least squares method, the optimal estimate of the load gravity and centre of gravity coordinates is:

Compensation of Force Disturbance
According to the proposed method, the zero-point f zero (p) = (p; θ) of the force sensor is predicted by using the NN model when the robot is in an arbitrary attitude p. In addition, the optimal value of the tool gravity in the robot base frame B G * and the tool centre of gravity in the force sensor frame r * can be quickly estimated in a small number of poses. Further, when the robot is operated to an arbitrary attitude p, the effect of the tool gravity on the sensor output is calculated as follows: According to Equation (2), the pure contact force applied to the end-effector is: Finally, the compensation of force disturbance can be completed according to Equation (16). See Appendix A for details of the pseudo-code of the proposed integrated compensation method.

Experimental Results
In order to verify the feasibility of the method, this paper designs the force interference compensation experiment of the six-axis force sensor at the end of the robot. The flow chart of the experiment is shown in Figure 4. The robot used in the experiment is the UR10 collaborative robot. The pose of the six-axis force sensor can be provided from the robot control system by positioning the installation of the six-axis force sensor and setting the position of the robot's TCP (Tool Centre Point). The force sensor is a contact force sensor in a tandem force sensor [33] independently developed by our laboratory, using self-developed acquisition box and software for data collection. The technical parameters of the force sensor are shown in Table 1. A total of 200 sets of data are collected in each posture to ensure the stability and accuracy of the data, and eliminate the deviation caused by the robot vibration and inertia. Finally, the median value is taken as the response value in that posture. pensation method.

Experimental Results
In order to verify the feasibility of the method, this paper designs the force interference compensation experiment of the six-axis force sensor at the end of the robot. The flow chart of the experiment is shown in Figure 4. The robot used in the experiment is the UR10 collaborative robot. The pose of the six-axis force sensor can be provided from the robot control system by positioning the installation of the six-axis force sensor and setting the position of the robot's TCP (Tool Centre Point). The force sensor is a contact force sensor in a tandem force sensor [33] independently developed by our laboratory, using self-developed acquisition box and software for data collection. The technical parameters of the force sensor are shown in Table 1. A total of 200 sets of data are collected in each posture to ensure the stability and accuracy of the data, and eliminate the deviation caused by the robot vibration and inertia. Finally, the median value is taken as the response value in that posture.

Zero-Point Estimation of a Six-Axis Force Sensor Based on Deep Learning
It is necessary to collect sufficient output of the force sensor in different poses to learn the parameter θ, however, there is no smooth continuous trajectory that allows the robot to traverse all possible poses in the workspace. Xiong [34] et al. pointed out that the first three joints of the robot mainly determine the location of the wrist, and the latter three joints mainly determine the posture of the wrist. Considering the time required to collect the data so that the robot can traverse as many poses as possible while moving at a small angle each time, the following experimental scheme is designed. Keeping the first three joints of the robot fixed and considering the interference problem of the force sensor during the motion and the common working range, the data acquisition range of the fourth, fifth and sixth joints of the robot is restricted to [−135 respectively, as shown in Figure 5. The step size of the three joints during data collection is 10 • , 10 • , 15 • , so a total of 15 * 29 * 24 = 10, 440 sets of data are collected. Each set of collected data includes the robot's pose and the output of the force sensor in the current pose. After the original data collection, the order of the data is disordered randomly, and 9966 sets of data are selected as the training set, while the remaining 474 sets are used as the test set. respectively, as shown in Figure 5. The step size of the three joints during data collection is 10 ,10 ,15 o o o , so a total of 15 29 24 10, 440 * * = sets of data are collected. Each set of collected data includes the robot's pose and the output of the force sensor in the current pose. After the original data collection, the order of the data is disordered randomly, and 9966 sets of data are selected as the training set, while the remaining 474 sets are used as the test set. The experimental environment is as follows: CPU is Intel(R) Core (TM) i7-9750H CPU@2.60 GHz, graphics card is GTX1660Ti, operating system is Windows 10, and Pytorch deep learning framework is used. The complete training data set of 9966 is used in each iteration, and the optimizer uses Adam with 100,000 iterations. The learning rate settings for the iterations are shown in Table 2. To prevent overfitting, L2 regularization is added to the loss function. Figure 6 shows the change of loss function for the 50,000th-100,000th iteration during the training. The experimental environment is as follows: CPU is Intel(R) Core (TM) i7-9750H CPU@2.60 GHz, graphics card is GTX1660Ti, operating system is Windows 10, and Pytorch deep learning framework is used. The complete training data set of 9966 is used in each iteration, and the optimizer uses Adam with 100,000 iterations. The learning rate settings for the iterations are shown in Table 2. To prevent overfitting, L2 regularization is added to the loss function. Figure 6 shows the change of loss function for the 50,000th-100,000th iteration during the training.

Number of Iterations
Learning Rate 0th-10,000th 0.005 After training, the model is tested using the test set data, and the results are shown in Figure 7, where the unit of forces is Newton (N) and the unit of moments is Newton-cm (Ncm). To illustrate the distribution of the zero-point estimation errors, a kernel density distribution of the errors was made, as shown in Figure 8. For comparison, a zero-point estimation method of the force sensor based on a dynamical model (DM) is conducted, and the results are shown in Figure 9. The mean (µ error ), standard deviation (σ error ), and maximum value of the absolute value of the error (max error ) obtained from the data are shown in Table 3.      Compared to Figure 9, the mean value of the error data on the test set in Figure 7 is close to 0. As shown in Table 3, max error of the test dataset of the NN method is much less than that of the DM method, except for Mx, but the errors are extremely similar. In particular, the error σ of the NN method is much smaller than that of the DM method. The data shows that the error data distribution of the NN method is more concentrated and uniformly disordered on both sides of 0, as can also be illustrated by Figure 8.
To verify that the NN method outperforms the DM method, we reduce the size of the training set for experiments. The experimental results are presented in Appendix B. The bar graphs of 3 error error μ σ + are shown in Figure 10 for training data of 1000, 3000, 5000, and 7000 sets, respectively. Figure 10 shows that the DM method still performs much worse than the NN method even with the same training samples. Furthermore, the result shows that the NN model generalizes well and the errors obtained after fitting are caused by noise, so the deep learning-based sensor zero-point estimation method is feasible.  Compared to Figure 9, the mean value of the error data on the test set in Figure 7 is close to 0. As shown in Table 3, max error of the test dataset of the NN method is much less than that of the DM method, except for Mx, but the errors are extremely similar. In particular, the σ error of the NN method is much smaller than that of the DM method. The data shows that the error data distribution of the NN method is more concentrated and uniformly disordered on both sides of 0, as can also be illustrated by Figure 8.
To verify that the NN method outperforms the DM method, we reduce the size of the training set for experiments. The experimental results are presented in Appendix B. The bar graphs of µ error + 3 σ error are shown in Figure 10 for training data of 1000, 3000, 5000, and 7000 sets, respectively. Figure 10 shows that the DM method still performs much worse than the NN method even with the same training samples. Furthermore, the result shows that the NN model generalizes well and the errors obtained after fitting are caused by noise, so the deep learning-based sensor zero-point estimation method is feasible. Sensors 2021, 21, x FOR PEER REVIEW 12 of 18

Tool Load Identification by Least Squares Method
In order to verify the feasibility of the least squares-based tool load identification algorithm, a vacuum chuck is installed at the end of the six-axis force sensor to perform the experiment. The experimental site is shown in Figure 11. The robot is controlled by the host computer to adjust to 15 attitudes (rotation vectors) to collect data, respectively. The attitude data are shown in Appendix C.

Tool Load Identification by Least Squares Method
In order to verify the feasibility of the least squares-based tool load identification algorithm, a vacuum chuck is installed at the end of the six-axis force sensor to perform the experiment. The experimental site is shown in Figure 11. The robot is controlled by the host computer to adjust to 15 attitudes (rotation vectors) to collect data, respectively. The attitude data are shown in Appendix C.

Tool Load Identification by Least Squares Method
In order to verify the feasibility of the least squares-based tool load identification a gorithm, a vacuum chuck is installed at the end of the six-axis force sensor to perform th experiment. The experimental site is shown in Figure 11. The robot is controlled by th host computer to adjust to 15 attitudes (rotation vectors) to collect data, respectively. Th attitude data are shown in Appendix C.   as expected. In reality, the mass of the vacuum chuck is 571 g, the gravity is 5.592 N. The gravity obtained from the identification is | B G * | = 5.117 N, the error of the identification result is less than 8.5%. However, because the true value of the position of the centre of gravity cannot be measured accurately, and it is only used as the intermediate value in the compensation process, the error of the obtained result is not analysed. Therefore, we used Equations (15) and (16) in Section 3.3 to calculate the compensation results to analyse the feasibility of the method.

Compensation Results
When the tool at the end of the six-axis force sensor is not in contact with any object, the sensor does not receive any external force except for gravity, that is f contact = 0. According to Formula (16), the error after the compensation is: According to the proposed method, the robot is adjusted to a series of random postures in the work area and collect force sensor data to verify the feasibility of the algorithm. The robot is controlled to move to 215 different random poses in the work area, and the corresponding robot pose data and force sensor output data are collected. The error between the response values f of the force sensor and the predicted values of the algorithm f zero + f t_gravity , i.e., the compensation error of the disturbance force, is calculated and the box plot of the error is shown in Figure 12. The figure shows that the error distribution of the test data is relatively concentrated, and the percentage of abnormal data is less than 0.9%. However, the error values of Mx and My are large relative to the forces, which is due to the fact that they are influenced by both the forces and the load centre of gravity, resulting in an unavoidable accumulation of errors.
Sensors 2021, 21, x FOR PEER REVIEW 13 of the centre of gravity cannot be measured accurately, and it is only used as the intermedia value in the compensation process, the error of the obtained result is not analysed. Ther fore, we used Equations (15) and (16) in Section 3.3 to calculate the compensation resul to analyse the feasibility of the method.

Compensation Results
When the tool at the end of the six-axis force sensor is not in contact with any objec the sensor does not receive any external force except for gravity, that is . A cording to Formula (16), the error after the compensation is: According to the proposed method, the robot is adjusted to a series of random po tures in the work area and collect force sensor data to verify the feasibility of the algorithm The robot is controlled to move to 215 different random poses in the work area, and th corresponding robot pose data and force sensor output data are collected. The error b tween the response values f of the force sensor and the predicted values of the alg rithm _ zero t gravity + f f , i.e., the compensation error of the disturbance force, is calculate and the box plot of the error is shown in Figure 12. The figure shows that the error distr bution of the test data is relatively concentrated, and the percentage of abnormal data less than 0.9%. However, the error values of Mx and My are large relative to the force which is due to the fact that they are influenced by both the forces and the load centre gravity, resulting in an unavoidable accumulation of errors.  Table 4 shows the maximum absolute error (MAX) and the mean absolute err (MAE) of the compensation, where Bias+LSM is the compensation algorithm proposed the literature [22], Double-LSM is the method proposed in the literature [23], which co siders the sensor and the tool load as a whole, and NN + LSM is the proposed metho Table 4 shows that, compared with Bias + LSM and Double-LSM, the proposed method this paper can effectively reduce the MAX and the MAE of , , , , Fx Fy Fz Mx My after th compensation. The compensation error level of Mz is comparable to that of the Doubl LSM, which is due to the fact that the zero-point estimation in this algorithm does n  Table 4 shows the maximum absolute error (MAX) and the mean absolute error (MAE) of the compensation, where Bias+LSM is the compensation algorithm proposed in the literature [22], Double-LSM is the method proposed in the literature [23], which considers the sensor and the tool load as a whole, and NN + LSM is the proposed method. Table 4 shows that, compared with Bias + LSM and Double-LSM, the proposed method in this paper can effectively reduce the MAX and the MAE of Fx, Fy, Fz, Mx, My after the compensation. The compensation error level of Mz is comparable to that of the Double-LSM, which is due to the fact that the zero-point estimation in this algorithm does not compensate for Mz and only compensates for it in the load gravity identification stage. However, in general, the present algorithm can accomplish the compensation of robot end force disturbance with higher accuracy. According to the technical parameters of the force sensor provided in Table 1, the errors after the compensation is completed by this method are 0.83%F.S., 1.39%F.S., 0.71%F.S., 0.45%F.S., 0.84%F.S., and 0.19%F.S., respectively.
To quantify the degree of error more accurately, the disturbance force compensation degree factor φ i is defined as: where n i is the output of the ith axis in the test. Therefore, this method reduces the effects of force sensor self-gravity, drift, system error, and tool gravity on the sensor output by 89.33%, 83.83%, 83.83%, 92.91%, 82.89%, and 76.48% in each dimension of the output, respectively. Through analysis, we conclude that the factors leading to incomplete compensation are: robot control error, the influence of the collection cable during the collection process, the calibration error of the six-axis force sensor, and the random error during collection. However, the compensated force error is still controlled within 1.50 N (1.5%F.S.) and the moment error is controlled within 4.3 Ncm (0.86%F.S.) by this algorithm.
Compared to the literature [24], the present algorithm is able to perform the compensation task efficiently while meeting the accuracy requirements. Although the test was conducted for the particular robot and end-effector, the method is general.
In conclusion, compared to existing methods, the proposed method can significantly reduce the influence of the end of the robot interference force, and can meet the demand of active compliance control.

Conclusions
In this paper, we proposed an integrated method to compensate for the external forces applied to the six-axis force sensor. Considering the interference of drift, system error, and gravity, the proposed method used a deep learning model to predict the zero-point of the sensor and identified the tool load by using the least squares method. In the experiment, we designed and trained an NN model to predict the zero-point of the sensor, based on a 9966-item dataset obtained in the non-load condition. Moreover, we collected 15 groups of data when the robot arm was in different poses, including robot poses and sensor output. Such data were used to identify the tool load based on the least squares method. Finally, the force compensation values for the sensor were calculated by integrating the zero-point prediction and tool load identification obtained above. The experiment results show that the proposed method can perform more accurate and effective force compensation in different conditions compared with existing methods. This method can automatically complete force compensation of the six-axis force sensor in complex manufacturing scenarios, which can significantly improve the intelligence of the robot.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix B
Results for training data of different sizes. NN stands for Neural Network method and DM stands for Dynamic Model method. where l error = µ error + 3 σ error .

Appendix C
The robot postures data required in the experiment of tool load identification based on least squares method.