External force estimation of fracture reduction robot based on force residual method

Robot-assisted fracture reduction surgery should be highly safe and accurate. It is necessary to accurately determine and control the reduction force of the reduction robot. In this article, a fracture reduction robot is designed with the configuration of the 6-universal-prismatic-universal (6-UPU). The vector method is used to analyze the kinematics of the reduction robot. The Jacobian matrix of the reduction robot and the second-order influence coefficient matrix of the acceleration are obtained. The Lagrangian method is adopted to analyze the dynamic of the reduction robot. The muscle contraction force of the femoral shaft fracture is analyzed based on the Hill model to determine the fracture reduction force. According to changing of the reduction force during the fracture reduction operation, a kind of external force estimation method based on force residual analysis is proposed. When there is no external force sensor at the end of the reduction robot, this method can be used to detect the reduction force in real time. According to the displacement, velocity, and output torque of each branch chain, the reduction force during the reduction operation of the reduction robot can be estimated. A simulation system is conducted and the simulation results show that the fracture reduction force can be estimated and accurately tracked in real time, which is of great significance for the safe operation of the fracture reduction robot.


Introduction
The traditional fracture reduction surgery relies mainly on the manual operation of the surgeon to achieve reduction. Since the reduction surgery effect depends on the surgeon experience, the quality of the surgery is unstable and the surgical effect is unsatisfactory. 1 With the development of robotics and image navigation technology, the robot used for fracture reduction surgery has gradually become a hot topic. Compared with series robots, parallel robots have the advantages of high rigidity, high precision, and stable movement, which can be widely used as machine tools, vibration test platforms, precision operations, and medical robots. 2 During the robot-assisted fracture reduction surgery, the robot needs to provide the reduction force in order to balance the muscle force caused by the muscle force. 1 The reduction force needs to be estimated and limited in real time to avoid excessive stretching force to pull the muscle and nerve tissue, which may cause secondary injury to the patient. On the other hand, the collision is likely to occur among the peripheral tissue of distal and proximal of fracture during the reduction process, so the reduction force needs to be detected to avoid serious collisions. If the external force is detected by the force sensors fixed on the robot, they are not only expensive but also easy to be damaged and difficult to be fixed. Meanwhile, there are many noises in the measured information and so on. [3][4][5] To overcome the shortcomings above, the external force estimation method with sensorless has attracted widespread attention. At present, many scholars adopt the motor dynamics model to design the disturbance observer to estimate the external force. This method is simple in structure, and the external forces can be estimated by obtaining motor current and motor position information. An additional force sensor needn't to be installed on the robot system. [6][7][8][9] For the series manipulator, the contact force and contact position can be estimated simultaneously by using the nonlinear optimization method, and the algorithm is experimentally verified in the robot arm. 10 To achieve the direction adjustment at the force action point while obtaining the desired velocity, the residual force method is used to estimate and locate the contact force. A generalized hybrid force/velocity control method is adopted at the contact point. 11 For the shortcomings of easily generating noise using the disturbance observer, a disturbance observer based on Kalman filter is designed to estimate the external force. This method has the advantages of wide bandwidth, no sensor noise, and so on. 12,13 According to the robot dynamics model, the generalized momentum observer method is adopted to estimate the external force. This method can estimate the external force only by measuring the joint displacement and velocity information, and it does not need to calculate the acceleration and inertia matrix. 14-18 The contact particle filter method is used to detect and locate the external contact of the robot without an external sensor. Even if the external torque measurement is corrupted by noise, it can return a reasonable estimate value. 19 The Kalman filter is adopted to extend the generalized momentum observer, and the contact force and moment at the center point of the end tool are estimated online, which can increase robustness with respect to disturbances such as uncertainty in joint friction. 20 The neural network is trained as robot dynamics with or without external contact by the Levenberg-Marquardt algorithm, and the robot sensor is used to detect the collision quickly and efficiently. 21 The H 1 filter is used to estimate the external force of the manipulator. The dynamic model of the manipulator does not contain the acceleration signal, and the external force is regarded as the state variable of the state space equation. 22 To ensure the safety of interaction between humans and robots, a torque control strategy and the residual error of the tracking error are adopted to estimate the external force, which can avoid all unintended contact, and a suitable collision strategy is considered to reduce the damage caused by the collision. 23 Various sensorless external force estimation methods above are currently mainly used for series robots. There are few relevant literatures on the external force estimation of parallel robots.
For robot-assisted fracture reduction surgery, considering the operating environmental conditions, it is difficult to install the force sensor at the end of the reduction robot. In order to detect the force information in real time during reduction surgery operation for improving the force control and the surgery safety, the external force estimation method of the robot-assisted femoral fracture reduction is studied. The kinematics of the reduction robot is analyzed by the vector method, and the Jacobian matrix of the reduction robot and the second-order influence coefficient matrix of the acceleration are obtained. The Lagrangian method is used to analyze the dynamics.
Based on the Hill model, the muscle contraction force during the fracture region is analyzed to determine the fracture reduction force. A force residual observer for the robot-assisted fracture reduction system is constructed to detect the reduction force at the distal end of the fracture during the reduction operation. Simulation experiments are carried out to verify the effectiveness of the proposed method for estimating the fracture reduction force by the force residual observer. Simulation results show that the force residual observer only needs to obtain the displacement, velocity, and output torque of each electric cylinder of the parallel robot, then the fracture reduction force can be estimated and the change of the reduction force can be accurately tracked.

Reduction robot
The robot-assisted fracture reduction surgery system with 6-UPU configuration is shown in Figure 1. It consists of a distal gripper, a proximal gripper, a reduction robot, a lifting platform, a movable base, and an optical surgical navigation system. The proximal gripper is mounted on the operating bed, and the proximal end of the femur is fixed by the bone needle. The distal gripper is mainly used to fix the distal end of the femur. The distal gripper is mounted on the reduction robot and connects the distal end of the femur to the reduction robot by the bone needle. The reduction robot can translate and rotate along the x, y, and z axes, with six degrees of freedom to meet the motion requirements of fracture reduction. The movable base is mainly used to install the lifting platform and the robot controller, and the reduction robot is mounted on the lifting platform. The height of the robot can be adjusted manually by the lifting platform to meet the needs of robotic surgical operations. The optical surgical navigation system is mainly used for the positioning of the proximal and distal ends of the fracture during the fracture reduction.
The structural characteristics of the reduction robot with six-degree-of-freedom parallel robot are shown in Figure 2(a). It is mainly composed of six electric cylinders, an upper platform, and a lower platform. The upper platform is connected to the lower platform by the electric cylinders. The end of each electric cylinder is equipped with two universal hinges.
The structure parameters of the parallel robot are shown in Figure 2(b). The height of the bottom surface of the lower hinge center is h A , and the distance from the center plane of the upper hinge to the top of the upper platform is h B . During the movement of the robot, the center plane of the hinge is unchanged with respect to the position and orientation of the upper platform and the lower platform, respectively, so h A and h B are constant. Therefore, the bottom center plane of the lower hinge and the bottom surface of the lower platform can be converted by h A , and the center plane of the upper hinge and the top surface of the upper platform are converted by h B . A 0 A i has a length of R, and A'A 1 has an angle with the coordinate axis X A . The length of B 0 B i is r, and the angle between B'B 1 and the coordinate axis X B is q. L i is the electric cylinder length vector.
Using the method of homogeneous coordinate transformation, setting the position coordinate of the origin of the coordinate system of the upper platformfBg relative to the coordinate system of the lower platform fAg is ½ x y z T , and the orientation angle of fBg with respect to fAg is ½ a b g T The coordinate system fBg coincides with the coordinate system fAg, first rotates the a angle around the Z B axis, then rotates the b angle around the temporary axis Y 0 B , and then rotates the g angle around the temporary axis X 00 B , and finally the fBg coordinate system origin shift (x, y, z), the homogeneous coordinate transformation matrix is where cq ¼ cos q and sq ¼ sin q.
In the space coordinate system, the coordinates of any point on the moving coordinate system fBg can be  expressed in the fixed coordinate system fAg by the homogeneous coordinate transformation matrix T

Kinematic analysis
The generalized velocity of the moving platform of the parallel robot is expressed as Euler angle ½ a b g T represents the orientation of the moving platform, but the Euler angle's first and second derivatives with respect to time are not the angular velocity and angular acceleration of the moving platform.
The angular velocity v of the moving platform in the fixed coordinate system fAg and the Euler angle ½ a b g T satisfy the following relationship The forward speed of the parallel robot is where _ L i is the electric cylinder speed, and J is the firstorder influence coefficient matrix of the moving platform velocity to the input speed of the electric cylinder, which is the Jacobian matrix of the robot.
The expression of the positive acceleration solution of the parallel robot is 24 where € L i is the electric cylinder acceleration, H q o is the second-order influence coefficient matrix solution of acceleration, and H q

Dynamic analysis
The Lagrangian method is used to analyze the dynamics of the reduction robot. The reduction robot is divided into two subsystems: the moving platform and the six electric cylinders. The kinetic energy and potential energy of two subsystems are analyzed, respectively. The Lagrangian formulation of the robot can be established. Finally, the partial derivative of the Lagrangian formulation is derived, and the dynamic equation of the robot can be obtained. 25,26 The kinetic energy of the moving platform in the generalized coordinate space is where m M is the mass of the moving platform; and I x , I y , and I z are the moment of inertia of the moving platform about the x, y, and z axes, respectively. The parameters of the electric cylinder are shown in Figure 3. m 1 and m 2 are the mass of the cylinder and the piston, respectively. l i is the unit vector along the A i B i ! direction. l gi is the centroid position of the electric cylinder, and v gi is the centroid speed of the electric cylinder.
The kinetic energy of a single electric cylinder is wherel , and The kinetic energy of six electric cylinders is where . Therefore, equation (9) can be expressed as Then, the inertia matrix of the electric cylinder can be obtained as follows In the operating space, the Lagrangian dynamic equation can be expressed as where q, _ q , and € q are the displacement, velocity, and acceleration of the moving platform of the parallel robot in the operating space, respectively. F is the generalized force. M ðq Þ is the inertia matrix. C ðq ; _ q Þ is Coriolis force and centripetal force matrix. G ðq Þ is the gravity matrix vector.
The Coriolis force/central force coefficient matrix C ðq ; _ q Þ of the robot can be calculated 27 by The Coriolis force/central force coefficient matrix C M of the moving platform is calculated according to equation (13) and is expressed as follows where the parameters are detailed in the literature. 26 As the rotating angle of the piston with respect to the cylinder is relatively small, then the Coriolis force/centripetal force of the piston is very small. It is usually not considered, C L ¼ 0.
The gravity of the moving platform is: The gravity of the electric cylinder is where @L i @q is ith row of the inverse matrix of the Jacobian. The related parameters of the model are measured in the software, as shown in Table 1.
Generally, it is difficult to obtain the parameters of the moving platform. The displacement and velocity of the electric cylinder can be easily obtained from the encoder. The driving force can be obtained from the force sensing of the electric cylinder. The robot dynamics is expressed in the joint space as follows where L, _ L, and € L are the displacement, velocity, and acceleration of the electric cylinder of the parallel robot in the joint space, respectively. D ðq Þ, H ðq ; _ q Þ, and P ðq Þ represent the inertia matrix, the Coriolis force and centripetal force matrix, and the gravity matrix vector in the joint space, respectively, and t ¼ J T F is the driving force.
Analysis on fracture reduction force As there are tensor fascia lata and semitendinosus and semimembranosus in the fracture area, these muscles are still in a continuous state after fracture. The force generated by the muscle stretching is the reduction force, which is the most important load during the fracture reduction and has a great influence on the accuracy and effect of the reduction. 28 Compared with other complex muscle force models, the Hill model is a phenomenon model, which is more reasonable for evaluating of overall muscle force. 29 AV Hill proposed a skeletal muscle three-element model to describe muscle force and assumed that the muscle force F ct and the tendon force F se are instantaneously equal. 30 Considering that the role of skin is very limited, so the effect of skin is ignored. The Hill model is shown in Figure 4. The mathematical expression of the Hill model is where F me is the active contractile element of skeletal muscle. F pe is the passive elastic element of muscle connective tissue. F se is the series elastic element of tendons. aðtÞ is the pinnate angle between the muscle fiber and the tendon, and a 0 is the pinnate angle of the optimal muscle fiber length. L is the muscle length factor where L 0 is the optimal muscle fiber length. L c is the muscle fiber length after stretching. Before the femoral fracture is restored, the relaxation agent is injected into the fracture area in order to prevent the nerve impulses being transmitted to the muscles. At this time, the muscles have only passive stretching force. 31 The expression of muscle force during reduction is where L 0i , F 0i , and a 0i are constants; i ¼1, 2, . . . n, n is the number of muscles. F 0i is the maximum isometric contraction force. The fracture reduction robot is shown in Figure 5. The acting position and direction of the reduction force is shown, which is expressed as   As the reduction force is the largest along the direction of the femoral shaft axis, so the reduction force can be simplified as F c ¼ ½ F cx 0 0 0 0 0 T .
The maximum reduction force of the femoral shaft fracture is 396 N. Setting the feather angle between muscle fibers and tendons, a(t) is 0. 28 The optimal muscle fiber length L 0 is 90 mm. 29 Combined with the path planning of this subject, the maximum displacement along the x-axis during the reduction process is 25.94 mm, and the maximum length of muscle fibers after muscle stretching is 115.94 mm.
The reduction force of the robot changing with the xaxis distance can be calculated by equation (23), as shown in Figure 6. The reduction force gradually increases with the increase of the x-axis displacement. When the x-axis displacement is 25.94 mm, the reduction force reaches maximum value of 396 N.

External force estimation
The external force estimation based on the force residual is widely used in the series robot, which can realize the external force estimation and collision detection during the movement of the manipulator. 14,32,33 At present, for the fracture reduction robot, a force sensor is generally fixed at the end of the reduction robot to obtain the reduction force during the fracture reduction. 34,35 Due to the particularity of the fracture reduction surgical environment, installing a force sensor will limit the clamping of the distal segment of the broken bone. Therefore, the external force estimation method based on force residual is proposed for the fracture reduction robot, especially for the robot with the parallel configuration. The proposed method can effectively avoid the installation of sensors at the end of the parallel robot.
When the reduction robot collides with the environment during the fracture reduction, the generalized momentum of the robot is generally changed greatly. 36,37 The force residual method based on the generalized momentum is proposed in this article. The fracture reduction force is estimated in real time by detecting the displacement, velocity, and output torque of the electric cylinder during the robot reduction operation. The acceleration of the electric cylinder need not to be detected. The force residual observer is equivalent to using a virtual force sensor between the surgical robot and the fraction region.
Considering the reduction force, the dynamics equation is modified as where t c is the equivalent moment of the reduction force in the joint space. The reduction force F c acts on the distal end of the fracture and its equivalent moment t c can be transformed by the Jacobian matrix of the robot as follows where J c ðq Þ is the Jacobian matrix of the robot.
Defining the generalized momentum of the reduction robot is as Doing the partial deviation of equation (26) _ D ðq Þ À 2H ðq ; _ q Þ 2 R nÂn is an anti-symmetric matrix, and D ðq Þ is a symmetric and positive definite matrix. 26 Therefore, From equations (24), (27), and (28), the following equation can be obtained Since equation (29) shows that _ p contains the external moment component t c , the observer can be designed with generalized momentum to detect whether the robot has an abnormal collision.
The force residual observer is defined based on the momentum deviation to estimate the external force value, and the force residual signal r is defined as where K 1 > 0, K 2 > 0 is the gain matrix._ p is the momentum estimate value._ p ¼ H ðq ; _ q Þ T _ L þ t À P ðq Þ þ r . Equation (30) can be rewritten as Figure 6. The relationship between axial reduction force and muscle stretch length.
When only the equation above is used as the observer, a large delay and oscillation are generally generated. 33 So constructing an adjustment function is as follows Taking f c as the feedforward adjustment of equation (31), the effective residual observer is where K 3 > 0 is the gain matrix. The structure diagram of the force residual observer is designed, as shown in Figure 7.
Find the second-order partial derivative for equation (33) and bring the result into equation (29) For ith electric cylinder component of the parallel robot, the Laplace transform is performed on equation (34), and the transfer function of the observer system is From equation above, it can be seen that the force residual observer is a second-order system, whose output is the force residual observation value r, and r varies with the input value t c . The system can be stable by changing the observer gains K 1 , K 2 , and K 3 .
It can be concluded from equation (36) that the reduction robot can estimate the reduction force of the fraction region during the robot reduction operation by using the force residual observation value r.
Simulation MATLAB/Simulink 2018b software is adopted to simulate the robot-assisted femoral fracture reduction, and the force residual observer is designed to estimate in real time the reduction force during the robot reduction operation. Simulation system with force residual observer is shown in Figure 8. The desired path of the reduction robot is preoperative planned. Then, the displacement of each electric cylinder is calculated by the kinematics inverse solution based on the desired path. The proportion integration differentiation (PID) controller block controls the displacement of the linear electric cylinder, and the muscular force block is used to add an external reduction force to the marked point P c . The parallel robot block is a parallel robot simulation model obtained by introducing the 3-D model established by SolidWorks into MATLAB/Simulink, then the actual displacement, speed, and driving force of the electric cylinder can be obtained. Then, the force residual observer block estimates the reduction force of the marked point on the distal segment in real time and tracks the changing of the reduction force during the fracture reduction.
The reduction path is planned according to the fraction model by the A* algorithm, and the planned reduction path is indicated by " " in Figure 9. Selecting a feature point on the distal fracture as the marking point P c . During the simulation of the robot reduction operating, the moving path of the marking point will be recorded in real time. At the end of the reduction operation, the coordinates of the marking point P c are (25.94, 32.96, and 2.72 mm).

Simulation without reduction force
When the reduction force is not considered, the robotassisted fracture reduction is simulated. The actual movement trajectory of the marked point P c is obtained, which is indicated by " ," as shown in Figure 9. The error between the planned reduction path and the actual path of the marking point P c is shown in Figure 10. The maximum position error is about 1 mm, which is along the y-axis direction at the beginning of the robot reduction. At the end of the reduction operation, the position and orientation error is [À0.2 mm, 0.3 mm, 0.01 mm, 0 , 0.04 , 0.1 ] T . Simulation results show that the robot can reset with the PID controller according to the planned reduction path, and the closed-loop control of the robotassisted femoral fracture is feasible.
When the reduction force is not considered, the robot is only affected by its own inertial force, Coriolis force, and centrifugal force and gravity. Simulation results of the driving force of each cylinder is shown in Figure 11, which is mainly used to balance the robot gravity.
As shown in Figure 12, during the robot fracture reduction operating, due to the low speed of the robot, the inertial force, Coriolis force, and centrifugal force of the robot are much smaller than itself gravity and can be ignored.
The force residual observer is used to estimate the external force during the robot reduction operating, and the external force obtained includes the gravity and inertial force. In order to detect the external reduction force more clearly, the force residual observer needs to be initialized    before the robot reduction. Thus, the reduction force observer output value is approximately 0 N. When no reduction force is applied, the external force value outputted by the force residual observer is approximately 0 N, as shown in Figure 13.

Simulation with reduction force
Applying a reduction force to the fracture reduction robot, the simulation is performed. The position and direction of the reduction force are shown in Figure 5, and the variation is shown in Figure 6. The position and orientation error of the marker point P c is shown in Figure 14, which gradually increases with the increase in reduction force. The maximum error is about À1.3 mm, which occurs in the x-axis direction. Then, the magnitude of the reduction force remains unchanged. The reduction movement gradually tends to be stable with adjusting the PID controller. At the end of the reduction operation, the position and orientation error is approximately [À0.2 mm, 0.3 mm, 0.01 mm, 0 , 0 , 0 ] T . The results show that the robot could make the fracture distal segment move according to the planned reduction path by adopting the PID controller.
According to the simulation results, the driving force of each electric cylinder is shown in Figure 15. The driving    force of the electric cylinder gradually increases with increasing the fracture reduction force. At the end of the reduction operation, the maximum driving force of the electric cylinder reaches À411 N, which can be used as a reference for the selection of the electric cylinders.
The driving force of the electric cylinder is calculated by the Jacobian matrix and expressed in the generalized coordinate space. As shown in Figure 16 at the end point of the fracture reduction, the generalized driving force of the electric cylinder is [À394 N, À5 N, À106 N, À3 NÁm, À19 NÁm, À11 NÁm] T during the reduction process, and the driving force of the robot cylinder is reduction. The robot's gravity and fracture reduction force are balanced. The torque generated in Figure 16 is due to the fact that the center of gravity of the reduction robot gravity and the direction of the reduction force do not coincide with the direction of the origin of the coordinate of the reduction robot.
The fracture reduction force estimated by the force residual method for robot-assisted fracture reduction system is shown in Figure 17(a) and (b). At the reduction end point, the force-reduction observer estimates the restoring force to be [392 N, 0, 0, 0, 13 NÁm, 5 NÁm] T . The force residual observer along the x-axis direction measured as 392 N.
From Figure 6, the maximum of the reduction force of the femur is 396 N, and the error is about 4 N (1.01%). Comparing Figures 6 and 17, the trend of reduction force is consistent. The force residual method can effectively estimate the reduction force. As shown in Figure 17(b), there is external moment of the simulation results. This is because the direction of the reduction force does not pass the coordinate origin of the robot.

Conclusion
In order to improve the safety of robot-assisted fracture reduction operation, a kind of external force estimation method based on the force residual analysis is proposed. By this method, the reduction robot can detect the external force in real time without installing an external sensor. The external force estimation only needs to measure the displacement and velocity information of the robotic cylinder to detect the reduction force at the distal end of the fracture without detecting and calculating the acceleration. In order to verify the validity and correctness of the external force estimation method, the simulation of the robot-assisted fracture reduction system is conducted in the MATLAB/Simulink software. The external force estimation can detect the reduction force in real time and accurately track its changing during the reduction operating. At the end of the reduction operation, the error of the external force estimate value is 4 N. The proposed method has certain universality, which can be applied to estimate the external force not only the femoral fracture reduction but also the various fracture reduction, such as pelvis and tibia.
In future, the force/position hybrid control algorithm of the reduction robot will be studied based on the fracture reduction force detected in real time by the force residual analysis, which will improve the safety of the robotassisted fracture reduction.

Declaration of conflicting interests
The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.