External force estimation for industrial robots using configuration optimization

External force estimation for industrial robots can be applied to the scenes such as human–robot interaction and robot machining. The model-based methods have gained the attention of many researchers because they only need motor signals. However, the performance of the above-mentioned methods depends on the accuracy of robot inverse dynamic model (IDM). Traditional methods focus on improving the accuracy of IDM to obtain a better estimation performance. However, suppressing the negative impact of modelling error can be achieved without reducing the modelling error. Thus, this article proposes an external force estimation method that uses the semi-parametric friction model and applies configuration optimization based on Jacobian condition number (JCN) to reduce the modelling error and its negative impact. First, the IDM is identified. Second, the semi-parametric friction model refines the traditional model in the torque observer to improve the estimation accuracy. Third, the JCN is used as an evaluating indicator to optimize the robot configurations. Finally, several simulations and experiments on a 6-DoF industrial robot demonstrate the validity of the proposed method. This method can enhance the performance of online force estimation by up to 48.28%. In addition, the application of the proposed method is verified on the laboratory-developed machining platform.


Introduction
In recent years, industrial robots have been widely applied to non-interacting operations like palletizing, welding and spraying [1]. However, on the report of the International Federation of Robotics, only some industrial robots perform the contact applications, such as grasping, assembling and machining, where industrial robots must interact with the environment and sense external force [2]. To achieve a compliant and safe movement when interacting with the environment, the force-sensing capabilities are required for industrial robots [3][4][5][6]. However, most industrial robots don't possess joint torque sensors. The Force/Torque sensors are often used to acquire external force information. Nevertheless, a costly sensor will increase system cost and decrease the entire system stiffness because of the potential unstable mechanical integration of the Force/Torque sensor [7]. Therefore, the sensorless techniques have attracted more researchers, which will be of great use in the physical human-robot interaction, robot assembly and robot machining fields.
In the industrial scenes, the contact between the industrial robot and the environment often occurs on the robot end-effector, such as robot machining and end-effector hand guiding. The contacts which occur in other locations of the industrial robot are usually unexpected. Therefore, estimating the external force acting on the end-effector is of great significance for industrial applications. Zhang et al. [8] established a sensorless hand guiding scheme for industrial robots by estimating the external force. Linderoth et al. [9] proposed an external force perception method for the assembly operation. Katsumata et al. [10] proposed a fast robot identification method based on optimal exciting motion and applied the method to the external force estimated in contact painting tasks. A virtual force sensor was implemented on the industrial robot to estimate the polishing force when the robot is machining [11].
The model-free force estimation methods impose restrictions such as imprecise accuracy [12] and necessary prior experience [13]. Various force observers have been designed to estimate the external force based on the robot inverse dynamic model (IDM). The extended state observer [14], the extended sliding mode observer [15], and the high-order finite time observer [16] were transferred to handle the application of contact force estimation. Recently, some works in [17][18][19][20] have made lots of achievements in sensorless interaction force estimation using the Kalman filters. Roveda et al. [21] proposed the estimation algorithms based on the extended Kalman filter that estimated not only the external force but also the environment stiffness, which are very helpful for the application of robot force control. An external force estimation method [22] was presented to estimate any external torque, which was even not applied at the end-effector. This method has strong applicability. The sensorless techniques [14][15][16][17][18][19][20][21][22] enriched the design methods of interaction force estimation and obtained the accurate estimation of timevarying complex external force. The algorithm complexity and the difficulty in parameter tuning should be considered in the real-time force estimation. The generalized momentum observer (GMO) [23] took less computing cost to estimate the external force in joint space without joint acceleration so that the differential noise was not introduced into the estimation result. Furthermore, it is worth mentioning that the observer gain of GMO is easy to tune.
The estimation performance of the above methods depends on the modelling error. Because of the nonlinear effects [24] including backlash, elasticity, or friction, the factors that cause the IDM inaccuracy cannot be ignored. Among these factors, friction is the most dominant one. The robot joint friction is impacted by various factors, such as joint loading, motion direction and other nonlinear factors [25][26][27]. The different joint friction models were designed to improve the IDM accuracy. As the extension work to [26], Roveda et al. [28] compensated friction by optimizing its modelling parameters using the Bayesian optimization. This method minimized the required number of the parameter tuning and contributed to the robot control performance. The internal joint temperature variation was considered by adding temperature dependency in the joint friction model [29]. Madsen et al. [30] described the dynamic characteristics of the robot joint friction by using the refined generalized Maxwell model. The existing models can only be utilized under specified situations and with specific accuracy criteria. To some extent, using the complex friction model can increase the accuracy of friction modelling, but it makes the identification of dynamic parameters more difficult. Reducing the influence of other factors on the friction model improved the identification accuracy of the friction model. Thus, a local IDM identification method was proposed to improve the joint torques prediction accuracy [31]. This approach solved the interference caused by the interaction between the joints' friction and the robot configurations while sacrificing the working space of industrial robots. A simple friction model can reduce the online computing cost, which is conducive to the control algorithm with the real-time requirements. The classical Coulomb-Viscosity model is usually assumed to be the main components of the robot joint friction. In [32][33][34][35][36], the classical friction model was used in the related robotic researches such as control, identification, and estimation.
The neural networks [37][38][39] have achieved some success in robot applications due to their advantages in data-driven reasoning. Especially, the neural networks can be applied to the error compensation. Zuo et al. [40] proposed a sensorless external force detection method, in which the back propagation neural network (BPNN) is used to compensate for the joint torque error of the humanoid robot. Tu et al. [41] proposed a data modelling method that established the relationship between friction, joint velocity and load, by utilizing the BPNN. This method could approximate the nonlinear functions with arbitrary precision. However, a BPNN-based model of fitting the residual nonlinear part of complex friction to enhance the identified friction model has not been reported.
For machining tasks, such as milling, drilling, grinding, and polishing, five degrees of freedom are required, the task path and the normal at each surface point are given. However, due to the large workspace of the robot, the midpoint of the task path and the normal of the end-effector at the start surface point relative to the base frame can be adjusted and optimized. The robot can achieve better configurations when performing the machining task by adjusting the relative relationship between the task path and the robot base frame. Guo et al. [42] presented a posture optimization method to strengthen the stiffness of the industrial robot in drilling processing. Mousavi et al. [43] improved the robot's dynamic behaviour during the milling process and enhanced the robotic machining stability by controlling the robot configuration. Lin et al. [44] optimized the spindle configuration of the robot end-effector to minimize the end-effector deformation.
As indicated in [45][46][47][48][49], 6-DOF industrial robots have a redundant degree of freedom when completing the above machining tasks. Even if the relative relationship between the task path and the robot base frame can't be adjusted, the robot configuration optimization can still be achieved. Due to the existence of a redundant degree of freedom, industrial robots have a redundant angle. For curved surface machining with more constraints, the robot configurations can be optimized by choosing the optimized redundant angle.
The robot Jacobian matrix is often used to characterize the manoeuvrability at a certain point. For instance, the Cartesian stiffness matrix of the robot depends on the Jacobian matrix [50]. The Jacobian condition number (JCN) is one of the important characteristics of the Jacobian matrix and can be used as an evaluation indicator in robot configuration optimization. Zargarbashi et al. [51] chose the JCN as a dexterity index of the conventional industrial robot. A small JCN led to a smoother joint-rate time history, which contributed to the stability of robot processing. Zhang et al. [52] used the JCN to calculate the kinetostatic conditioning index which could be used to quantitatively assess the kinematics performance of the industrial robot. However, the relationship between the JCN and the force estimation error has not yet been investigated. And the improvement of the force estimation accuracy has not been reported through configuration optimization too.
In this article, an external force estimation method, that uses the semi-parametric friction model and applies configuration optimization, is proposed to refine the force estimation performance. The main contributions are listed as follows: (1) The proposed configuration optimization builds the relationship between the JCN and the external force estimation error. The configuration optimization reduces the JCN, which effectively improves the accuracy of force estimation by suppressing the impact generated by the modelling error and noise. The rest of the article is arranged as follows. Section 2 introduces the preliminaries. Section 3 proposes the external force estimation method for specific tasks using the GMO. Section 4 provides the simulation and experimental results. Subsequently, the applications in polishing and grinding are demonstrated in Section 5. At last, Section 6 concludes the article and provides further ideas for future works.

The robot inverse dynamic model and identification
When the robot moves in free space and does not contact the external environment, the joint torques vector τ d can be calculated by the n joints robot IDM which is a function of joint positions and their derivatives: (1) where M(q) ∈ n×n is the inertia matrix of robot links, C(q,q) ∈ n×n is the Coriolis/centrifugal matrix, G(q) ∈ n is the gravity vector, f ∈ n is the joint friction vector, and q,q,q ∈ n are the joint position, velocity, and acceleration vectors respectively. For the ith joint, the joint friction torque is usually described by the traditional friction model which is only related to velocity: where Fc i is the Coulomb's friction, Fv i is the viscous friction coefficient. Many friction models can be selected, and it is worth studying to choose a suitable friction model to improve the accuracy of IDM. The dynamic parameter identification method can be found in the previous works [53,54]. The joint torques can be predicted when the minimal parameter set is identified. It can be assumed approximately that the dynamic parameters of the robot except for the friction are accurately identified. The IDM expressed in Equation (1) can be solved by the method proposed in [55], though the dynamic parameter identification obtains a linear IDM with the minimal parameter set.

The generalized momentum observer
To meet the requirements of real-time estimation, the GMO is chosen due to its advantages in computing cost and parameter tuning. As mentioned in [23], the GMO is constructed based on the robot IDM.
Considering the requirements for industrial scenes, this article focuses exclusively on the external force acting on the robot end-effector and overlooks the external force acting on other locations of the industrial robot. When an external force is exerted to the end-effector, the joint torques of the robot are expressed as: where τ d is the torques generated by the robot's own motion and the gravity field, τ J is the total joint torques which can be obtained by the motor current, τ J consists of two parts, one of which is used to offset the joint torques τ ext generated by the external force, and the other one is used to generate τ d which drives the link to move. From Equation (3), when the current information is obtained and the IDM calculation is completed, the joint torques generated by the external force can be expressed through: Substituting Equation (1) into Equation (4) yields: From Equation (5), it is obvious that the signal-tonoise ratio of joint acceleration will make the calculation of τ ext unreliable. The GMO method is used to address the above problem, as presented in Figure 1.
The GMO method utilizes the difference between the model generalized momentum and the measured generalized momentum to construct a torque observer. The joint torques τ ext can be observed as: where K 0 is the observer gain matrix. The measured generalized momentum is defined as: Derivative on both sides of Equation (7), the following formula can be obtained: Using Equation (4) and integrating, the model generalized momentum can be constructed as: q is a skew-symmetric matrix, and the symmetry of M(q) implies: Substituting Equation (10) into Equation (9) yields:

The proposed external force estimation method
The proposed method contains two parts: the GMO with a semi-parametric friction model and taskoriented configuration optimization based on JCN, as shown in Figure 2. After the dynamic parameter identification, the GMO is constructed based on the robot IDM. The semi-parametric friction model enhances the identified friction model and is added to the GMO to improve the performance of the external force estimation. The configuration optimization based on the performance index using the JCN suppresses the impact of the modelling error on the external force estimation. The robot moves along the optimized task path and estimates the external force in the movement process.

The semi-parametric friction model
Compared with works in [37][38][39][40], the semi-parametric friction model refines the identified friction model by the BPNN-based model, which is used to learn and fit the residual nonlinear parts in the complex friction, as shown in Figure 3. The semi-parametric friction model is composed of different models when the joint velocity belongs to different intervals. Furthermore, the dataset generation method is introduced, which is used to train the BPNN in the semi-parametric friction model. When the robot joint velocity is not near-zero, the joint friction is described by the traditional friction model which is previously identified. The BPNN model is used to establish the mapping relationship between joint friction torque and joint velocity when the joint velocity is near-zero. The input is the joint velocity. The output is the joint friction torque. An input layer (1 output neuron), a hidden layer (10 hidden neurons), and an output layer (1 input neuron) constitute the BPNN model. The hyperbolic tangent sigmoid transfer function and the linear transfer function are selected for the activation function of the hidden layer and the output layer respectively. The number of hidden layers and the number of neurons in the hidden layer represent the complexity of the network. The above suitable sizes of the hidden layer and the neuron are chosen  through multiple training trials to prevent the overfitting phenomenon. Since the semi-parametric friction model is a local supplement to the identified friction model, it is dispensable to use the neural network with complex structure to learn the friction characteristics when the joint velocity is near-zero. BPNN has its advantages in the complexity of the neural network structure and the real-time requirements of the robot controller. Then, the traditional friction model and the BPNN model are merged into a semi-parametric friction model, expressed as: where τ fi (q i ) denotes the friction torque of ith joint described by the semi-parametric friction model, BPNN(q i ) denotes the ith joint friction calculated by the BPNN model at near-zero velocity, andq crit denotes the critical joint velocity. For each joint, an independent BPNN model is trained and applied to predict the joint friction torque. The generation of dataset is an important step to training the BPNN model, which needs the friction data and the velocity data of the robot joint when the joint velocity is near-zero. The joint friction torque can be obtained through characteristics of the IDM and the design of the robot configuration and kinematic state. The motor of each joint of industrial robots is equipped with encoders, so the joint velocity can be directly read. Next, the method of obtaining joint friction torque from the total joint torques is introduced.
As mentioned in [55], the (u, v) element of the matrix C(q,q) is defined as: The terms M(q) uv denotes the (u, v) element of the matrix M(q). From Equation (13), it can be derived that the Coriolis/centrifugal matrix satisfies: For simplicity, it can be assumed that the joint friction torques are only related to the magnitude and direction of joint velocity: Assume two robot configurations, q 1 , q 2 meet the following conditions: Substituting Equation (16) into Equation (1), it can be obtained that: Substituting Equation (15) into Equation (17), it can be found that: According to Equation (18), it can be obtained that: From Equation (19), the joint friction torque at nearzero velocity is separated from the measured joint torque. The joint friction torque corresponding to the joint velocity can be obtained.

Task-oriented configuration optimization
Based on the JCN, the task-oriented configuration optimization is applied to refine the force estimation performance. It is assumed that only the tool centre point (TCP) of the robot end-effector is in contact with the environment. The transformation relation [56] of external force in Cartesian space and joint space can be expressed as follows: (20) where J(q) T is the transpose of the robot Jacobian matrix, F ext is the external force applied to the endeffector in Cartesian space.
According to the properties of condition number [57] and Equation (20), the error of external force estimation δF ext satisfies: Moreover, the JCN is defined as cond , δτ ext is defined as δτ ext = τ ext − τ ext , and denotes the estimation error of the joint torque τ ext in Equation (11), δF ext denotes the estimation error of the external force. From Equation (11),τ ext depends on the IDM accuracy. Therefore, it is assumed that δτ ext is generated by robot dynamic modelling error, and bounded and less than the maximum value of the modelling error. When the robot IDM is established and identified, the ratio of δτ ext to τ ext remains basically the same. The 2-norm is chosen as the Jacobian matrix norm. Obviously, when the values of modelling error in the different configurations are the almost same, a smaller JCN will lead to a smaller error in force estimation. Even if the modelling error is not reduced, the force estimation error can be suppressed by reducing the JCN. From Equation (20), the external force applied to the TCP can be reconstructed as: (22) Considering a common machining scene where the robot performs a polishing task along a linear path. For this scene, the direction of the path can be regulated and the end-effector pose in the base frame can be modified too. By changing the path and the end-effector pose, the robot can obtain an average minimum JCN on this path, so as to minimize the force estimation error on this path. The configuration optimization problem can be transformed into a constrained optimization problem which will be explained as follows.
The path is determined by the direction vector n and the midpoint coordinate m in the base frame which can be expressed as: where n denotes the direction vector, m denotes the midpoint coordinates of the path. Suppose the length of the path is a, the path can be divided into k parts of equal length by evenly inserting k + 1 interpolation points. The coordinates of the ith interpolation point on the path can be described as: where i = 0, 1, 2, . . . , k. Assuming that the end-effector pose does not change during the polishing process, the end-effector pose at the ith interpolation point can be described as: with R ∈ 3×3 being the rotation matrix representing the initial end-effector pose, which is relative to the base frame. The rotation matrix R can be described in the axis-angle representation as θ r = θ r x , r y , r z T , θ is the size of the rotation angle, and r indicates the direction of the rotation axis. From Equations (24) and (25), the homogeneous transformation matrix of the end-effector corresponding to each interpolation point on the path can be expressed as: According to Equation (26), the corresponding robot configuration q i can be obtained through the robot inverse kinematics. The average JCN on this path can be obtained and taken as the objective function. Thus, the configuration optimization problem can be expressed as follows: where q i ∈ S denotes that each interpolation point on the path belongs to the robot's workspace, [m x , m y , m z ] T ∈ V denotes that the position of the midpoint is constrained. This constraint is related to the mutual positional relationship between the workpiece and the robot, and θ min ≤ θ ≤ θ max is the range of the rotation angle. Other constraints indicate the normalization conditions of the unit vector. The constrained optimization problem (27) can be solved by the PSO-GA algorithm. The solving process of task-oriented configuration optimization is shown in Figure 4. When the relative relationship between the task path and the robot base frame can be adjusted, (27) can be solved. If the above requirements are not met, the configuration optimization can still be achieved by the redundant angle optimization. The final output is the optimized configuration.

Simulations and experiments
In this section, several simulations are completed to verify the validity of the proposed method. Based on the simulation results, the experiments are applied to the STS-SRE4-600 industrial robot. As shown in Figure  5, the kinematic model of the STS-SRE4-600 is established by the standard DH convention. The parameters of the standard DH convention are shown in Table 1.
The DH parameters are calibrated by using the Leica  Link

Experiment of the IDM identification
The force estimation observer is designed based on the robot IDM. The precise dynamic parameters of the robot are essential to the force estimation results. The identification methods introduced in works [53,54] are used to obtain the dynamic parameters. The friction model used in the identification experiment is the classical Coulomb-Viscous friction model. The excitation trajectory is designed based on the optimized fifth-order Fourier series. The robot executes the excitation trajectory, then collects the joint kinematics state and current signal to complete the IDM identification.
A new trajectory is used for the cross-validation of the identification results. The five-order Butterworth low-pass filter is used for the raw data collected in the dynamic parameter identification experiment. The filter has a cut-off frequency of 7.5 Hz, which is 15 times the highest frequency of the excitation trajectory. The designed filter can ensure the accuracy of the collected data while removing the noise. The IDM identification result is shown in Figure 6. The prediction error is assessed by the root mean square error (RMSE). The prediction errors for each joint torque are 3.2926, 3.5804, 2.2578, 0.5272, 0.7846, and 0.3602 N·m, which are all within the allowable ranges. The similarity between the joint torques calculated by the IDM and the measured joint torques based on the current signal is more than 90%. The remaining errors are difficult to reduce only using the dynamic parameter identification.

Simulations of the JCN distribution
The simulations verify the influence of the JCN in the external force estimation. To illustrate the effectiveness of the configuration optimization, the robot's workspace and the distribution of the JCN corresponding to the robot configuration are obtained. The JCN is affected not just by the position of the end-effector, but also by its orientation. When the endeffector is in the same position, the JCN value will be different due to the different orientations of the robot end-effector. The different configurations contain different positions and orientations. Half a million robot configurations are randomly generated by Monte Carlo method to obtain the robot's workspace and the JCN corresponding to each configuration. Figure 7 shows the distribution of the JCN in the robot's workspace. The colourbar is used to represent the numerical value of the JCN corresponding to each configuration.
The average value of the JCN corresponding to each configuration is 1052.9. The minimum value of the JCN is 9.1675. It can be seen that the JCNs corresponding to 96.55% of the robot configurations are greater than 11. There is still much room for improvement by optimizing the robot configuration to reduce the JCN.

Simulations of the configuration optimization for linear task paths
In Section 3.2, the task-oriented configuration optimization needs to obtain the optimized task path and the optimized end-effector pose by solving (27). For simplicity, the length of the task path is set as one metre in simulation. The base frame coincides with the robot base frame (frame 0) in Figure 5. The following direction vector, midpoint coordinate and rotation matrix are all described in the base frame.
The initial direction vector of the path is set as n in = [0, 1, 0] T . The initial midpoint coordinates of the path is set as m in = [0.5, 0, 0.45] T . Referring to the commonly used end-effector pose [47] in robot machining, the initial end-effector pose is set as: Several local optimal solutions will be obtained by solving (27). It is necessary to select the appropriate solution or adjust the local optimal solution according to the actual task scene and the machining processing. The optimized direction vector is solved For the determined polishing task, when the endeffector pose changes, the flange connecting the endeffector also needs to be replaced.
The following two simulations, as shown in Figure  8, are used to verify the influence of task-oriented configuration optimization on the force estimation. It is supposed that a constant force is exerted at the endeffector when the robot moves back and forth along the task path. The constant force is set as F ext = [ −3N, 0N, 25N, 0.5N · m, 1.5N · m, 1N · m] T . In the first simulation, the robot performs the task with the initial configurations along the initial path. In the second simulation, the robot completes the task with the optimized configurations along the optimized path. The blue line in Figure 8 shows the JCN value at each point on the task path. In the simulations, it is assumed that the model uncertainty and the noise present a random noise distribution. The mean value of random noise is 0. The amplitude of random noise is 15% of the torques in the joint space generated by the preset constant force. At sample point 0, the robot begins to move at the start point. At sample point 1500, the robot moves to the end point of the task path. At sample point 3000, the robot returns to the start point of the task path. The robot trajectory is planned by uniform acceleration and deceleration. The estimation results are shown in Figure 9. The comparison of the estimation error is shown in Figure 10.
After optimization, the average JCN of the task path has dropped from 20.2378 to 10.8713, with a 46.28% decrease. The maximum JCN in the path has dropped from 59.5668 to 17.5814, with a 70.48% decrease. The estimation accuracy in Y-direction is significantly affected by the JCN than that in other directions. Figure 9 illustrates that when the robot moves along the task path, the larger the JCN at the point on the task path, the larger the external force estimation at that point, especially in Y-direction. From Figure 10, the maximum estimation error has dropped from 24.86 N to 7.30 N after the optimization, with a 70.63% decrease. However, the estimated external torques in X-and Y-directions are insensitive to the configuration optimization.

Simulations of the configuration optimization for task paths on the surface
When the relative relationship between the task path and the robot base frame can't be adjusted, the robot configuration optimization can still be achieved. Due to the existence of a redundant angle, the robot configurations can be optimized by choosing the optimized redundant angle. The redundant angle optimization is based on the orthogonal decomposition method as introduced in [49].
The comparisons of the JCN, robot end-effector poses, and robot joint positions before and after the configuration optimization are shown in Figures 11-13, respectively. From Figure 11, the JCN at each point on the task path is reduced. From Figure 12, the robot configuration optimization is accomplished when the task path and the normal of the end-effector are provided. The optimized robot trajectory maintains the same position and the same normal as the initial trajectory at each sample point. The above configuration optimization can be applied to the curved surface machining.

Experiments of force estimation with the task-oriented configuration optimization
Based on Section 4.3, the same settings, such as the end-effector pose and the task path, are used for experimental verification of the configuration optimization effect. The joint torque prediction error generated by the IDM comes from the model uncertainties and the noise. Thus, the control variates method is used to verify the influence of configuration optimization on the performance of external force estimation.
The design of the experiment should satisfy the joint torque prediction error generated by the IDM is consistent when the robot performs external force estimation before and after configuration optimization. The requirement can be achieved by multiple round trip movements of the robot along the same path. The experimental setup is shown in Figure 14, which contains 5 components: an industrial robot, an ATI Gamma Force/Torque sensor, one OPT polishing tool serving as the end-effector, and two flanges.      The following experiments are shown in Figure 15. First let the robot move back and forth along the initial path with the initial pose, and record the joint torques in the first cycle. The robot doesn't have contact with the external environment or the experimenter in this process. The joint torques are used to maintain the robot's own motion. When the robot's movement enters the second and subsequent cycles, the torques recorded in the first cycle are regarded as the torques τ d in Equation (4). The torques τ J in Equation (4) can be obtained by the real-time motor current information. From Equations (4) and (20), the external force applied to the end-effector can be estimated. Similarly, when the robot moves back and forth along the optimized path with the optimized pose, the experiment can be performed again. In this way, the joint torque prediction error remains consistent in the two experiments.
In the first set of experiments, no external force is exerted to the end-effector. The estimation error is directly obtained when the robot moves along the task path. At 0 s, the robot starts to move at a constant speed. At 200 s, the robot moves to the end point of the task path and begins to move in reverse. At 400 s, the robot returns to the start point of the task path. From Section

4.3, the estimation accuracy in Y-direction is significantly affected by the JCN than that in other directions.
In the subsequent display of experimental results, the estimation results and the measurement results in Ydirection are used to evaluate the effect of the configuration optimization. The force estimation results before and after the configuration optimization are shown in Figure 16. Figure 16 shows that the maximum value of the estimation error in Y-direction has dropped from 49.3864 N to 14.9596 N, with a decrease of 69.71% after the configuration optimization. The RMSE of the estimated external force has dropped from 3.2520 N to 2.0449 N, with a decrease of 37.12%. The experiment results are consistent with the simulation results in Figure 9.
Repeat the above experiments, the difference is that the external force is applied to the end-effector through the human hand. In the subsequent experiments, the reference value of the external force is measured by the ATI F/T sensor. To remove the high-frequency measurement noise, the cutoff frequency for low-pass filtering is chosen as 18 Hz. The force estimation results before and after the configuration optimization are shown in Figure 17. Figure 17(a) shows that when the JCN is large, the result of external force estimation is unreliable under the same modelling error and noise. In 0-20 s, 190-210 s, 390-400 s, the robot is near the start and end of the task path, the JCN corresponding to the robot configuration is near the maximum, and the measured force is almost submerged in the noise of the estimation. Figure  17(b) shows that the estimation error at the start and end of the task path has been reduced after the configuration optimization. This suggests that when the robot is in a good configuration, the force estimation result will be better. After removing the abnormal data points, the RMSE of the estimated external force has dropped from 4.0661 N to 3.2274 N, with a decrease of 20.93%. The maximum estimation error has also dropped from 46.8357 N to 14.4221 N, with a decrease of 66.39%. Finally, the results suggest that even if the modelling error and noise are not reduced, the taskoriented configuration optimization can suppress the force estimation error.

Experiments of force estimation with the semi-parametric friction model
The GMO with a semi-parametric friction model is used for the real-time external force estimation. In the robot controller, the GMO algorithm is executed every millisecond. The external force is applied to the endeffector through the human hand. The observer gain matrix is set as a diagonal matrix with diagonal elements to 20. The result of the semi-parametric friction model is shown in Figure 18. It is assumed that the actual friction torque can be obtained from the total joint torque by the proposed separation method. Through the cross validation, the semi-parametric friction model is closer to the separated friction when the joint velocity is near-zero.
The estimation results using different friction models are shown in Figure 19. It is observed that the GMO with a semi-parametric friction model has a better estimation performance than the GMO with the traditional friction model when the joint velocity is near-zero. The   RMSE of the external force estimation in Z-direction has dropped from 8.1579 N to 5.4073 N, with a decrease of 33.72%. When the time is around 50 s, the estimation error presents peak value. From the zoom of Figure  19(a), the maximum error of external force estimation has also dropped from 49.7432 N to 20.4279 N, with a decrease of 58.93%.

Experiments of real-time external force estimation
The proposed external force estimation method that uses the semi-parametric friction model and applies the configuration optimization is validated in the following experiments. In the first experiment, the GMO with the traditional friction model is used to estimate the real-time external force when the robot performs the task with the initial configurations. In the second experiment, the GMO with the semi-parametric friction model is used to estimate the real-time external force when the robot performs the task with the optimized configurations. In the above experiments, the robot moves uniformly from the start point of the task path to the end point in 100 s. The performance of the external force estimation is evaluated by the RMSE and the correlation coefficient ρ. Figures 20 and 21 show the comparison of the estimation results in X-, Y-, and Z-directions.
Through the proposed method, the RMSE of the external force estimation in X-, Y-, and Z-directions

Discussion for the force estimation experiments
To select more appropriate observer gains, improve the rapidity and accuracy of force estimation, and reduce the interference of noise, the bandwidth of GMO is analysed. As mentioned in [23], the GMO is essentially equivalent to a first-order low-pass filter. Take the derivative on both sides of Equation (11), and then carry out Laplace transform, yield: By adjusting the different values of K 0 , the amplitude-frequency curves of the GMO are shown in Figure  22. The bandwidth analysis shows that choosing a larger observer gain K 0 will obtain a wider bandwidth for the estimated force. However, larger observer gains usually bring stronger observation noise. The GMO estimation results with different observer gains are shown in Figure 23. The result of force estimation is related not only to the observer gain but also to the type of external force. The impact of distinct external force types (step, ramp, and sine) and observer gains on external force estimate performance are explored as follows. The change rate of step-type external force can be considered as 0 after the step is completed, which can be well estimated without a high feedback gain. The change rate of the ramp-type force is a constant value. The estimation result will have a delay with a small observer gain. The change rate of sine-type force is a changing trigonometric function. There will be hysteresis and clipping in the estimation result when the observer gain is not inappropriate.
The proposed method, combined with specific machining tasks, applies configuration optimization to a new field, which can improve the force estimation accuracy without improving the IDM accuracy. Unlike previous works [50][51][52], a new paradigm to enhance the estimation accuracy is present which contributes to the index design for the robot configuration optimization. The proposed method assumes the force location, and can directly reconstruct the external force acting on the end-effector, which is of great significance for the control and feedback in robotic machining. Considering   the typical robotic machining scenes [46], the external torque has little influence on the machining, this study mainly focuses on the external force estimation.
Compared with the previous force estimation methods [18][19][20][21][22], the limitations of this study are that only the external force information can be estimated, not the environmental information, and it also encounters difficulties when multiple forces are acting on different positions of the robot.

Application of the proposed method in robotic machining
In this section, the proposed method is used to estimate the tangential force and normal force of the end-effector in machining. The impedance control algorithm is arranged on the robot controller, which makes the robot and the workpiece keep a desired force contact state on the contact surface. This can ensure that the contact force between the robot and the workpiece remains relatively stable, which is conducive to the stability of the implementation of the force estimation algorithm.

Real-time external force estimation in robotic polishing
For the scene of the robot polishing aluminium plate, the tangential force of the polishing head is estimated by the proposed method. Before the experiments, the polishing is planned in the simulation environment. The robot moved back and forth along the straight path to polish the aluminium plate at a constant speed. The motion simulation and experiment of robotic polishing are shown in Figure 24.
The normal contact force between the end-effector and the aluminium plate is set as 40 N, the polishing path length is 0.27 m. The end-effector carries out three experiments at different velocities. In the polishing experiments, the time for the end-effector to move from the start point to the end point of the polishing path is 15 s, 30 s and 60 s respectively. The tangential force is estimated and measured respectively. The experimental results of tangential force estimation are shown in Figure 25.
In the three experiments of the normal force estimation, the correlation coefficients of the tangential force estimation and the measured results are 0.9534, 0.9538 and 0.9415 respectively, and the RMSE of the external force estimation are 3.0845 N, 3.5468 N and 4.7082 N respectively.

Real-time external force estimation in robotic grinding
For the scene of the robot grinding wind turbine blade, the normal force estimation experiments of the grinding head are conducted. Before experiments, the grinding is planned in the simulation environment. The robot moves back and forth along the desired path to grind the wind turbine blade at a constant speed. The motion simulation and experiment of robotic grinding are shown in Figure 26.
The normal contact force between the end-effector and the wind turbine blade is set as 40 + 20sinϕN, where ϕ is related to the position of the end-effector. In the process of the end-effector moving from the beginning to the end of the grinding trajectory, the desired normal force changes in sine form for two cycles. There are two reasons for setting the desired normal force as a variable force. On the one hand, it is difficult to keep the normal force constant due to the influence of the curved surface in the actual grinding. On the other hand, setting the desired normal force as variable force can better evaluate the performance of the proposed method.
The grinding path length is 0.27 m. The robot carries out three experiments at different velocities. In the grinding experiments, the time for the robot to move from the start point to the end point of the grinding path is 15 s, 30 s and 60 s respectively. The normal force is estimated and measured respectively. The experimental results of external force estimation are shown in Figure 27.

Discussion for the force estimation applications of robotic machining
The above experimental results show that the force estimation accuracy is improved when the speed of the end-effector increases. When the robot moves faster, the modelling torques that drive the links to move increase compared with the unmodeled dynamic interference. The force estimation accuracy in curved surface grinding is lower than that in polishing. One possible explanation is that the joint velocities and joint accelerations of the robot in curved surface grinding are not as smooth as that in polishing, which may introduce disturbances and noises, resulting in the decline of the estimation accuracy. Furthermore, the grinding processing is different from the polishing processing, which may cause some stronger interference to estimation results. The low-pass filter is used to process the force signals collected by the force sensor, and eliminate the partial interference caused by the different machining processes.
In practical applications, the common limitation is that the dynamic parameters of the end-effector need to be previously given or identified. Changing the endeffector will impact the estimation results.

Conclusion
In this article, an external force estimation method that uses the semi-parametric friction model and adopts configuration optimization is proposed to improve the estimation accuracy. A semi-parametric friction model is introduced into the GMO to enhance the estimation performance by up to 33.72%. The dataset generation method for the BPNN model in the semi-parametric friction model is then given. The configuration optimization is applied to the external force estimation for the first time. The robot can obtain an average minimum JCN on the task path, which can minimize the external force estimation error on this path. In the realtime external force estimation, the RMSE of the proposed method in X-, Y-, and Z-directions are 6.5627, 7.6607, and 4.3901 N, respectively. Compared with the initial external force estimation, the estimation accuracy has been greatly improved with 48.28%, 46.53%, and 6.30%. For robotic polishing and grinding applications, the best correlation coefficient and RMSE of the proposed method in tangential force estimation are 0.9534 and 3.0845 N, respectively, and the best correlation coefficient and RMSE of the proposed method in normal force estimation are 0.7901 and 12.9954 N respectively. The application test shows that the proposed method is feasible for robot machining. The limitation of this work is that when the force applied to the robot changes abruptly, the force estimation results cannot be quickly achieved and effectively converge to the stable value, which may cause the steady-state error of the real-time control. In addition, the semi-parametric friction model is designed based on the classical friction model for real-time performance, which would reduce the generalization of the neural network and increase the training time. Future work will focus on further improving the estimation performance of the proposed method by enhancing the friction model and the observer structure. More accurate friction models will be used to characterize the main components of joint friction with inspiration from [26,28], reducing the model training time and improving model generalization ability.

Disclosure statement
No potential conflict of interest was reported by the author(s ).