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

Abstract

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. ¹ 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. ²

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. ¹ 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 –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 –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. ¹⁰ 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. ¹¹ 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. ¹⁹ 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. ²⁰ 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. ²¹ The $H_{\infty }$ 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. ²² 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. ²³ 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.

Figure 1.

Fracture reduction robot system.

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.

Figure 2.

(a) 6-UPU parallel robot. (b) Structure parameters.

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^{\prime }{A}_i$ has a length of R, and A’A ₁ has an angle ϕ with the coordinate axis X_A . The length of $B^{\prime }{B}_i$ is r, and the angle between B’B ₁ and the coordinate axis X_B is θ. 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 platform{B} relative to the coordinate system of the lower platform {A} is ${\left[\begin{array}{ccc}x& y& z\end{array}\right]}^T$ , and the orientation angle of {B} with respect to {A} is ${\left[\begin{array}{ccc}\alpha & \beta & \gamma \end{array}\right]}^T$ The coordinate system {B} coincides with the coordinate system {A}, first rotates the α angle around the Z_B axis, then rotates the β angle around the temporary axis $Y_B^{\text{'}}$ , and then rotates the γ angle around the temporary axis $X_B^{″}$ , and finally the {B} coordinate system origin shift (x, y, z), the homogeneous coordinate transformation matrix is

${}_B{}^A\mathit{T}=\left[\begin{array}{cc}\mathit{R}& \mathit{P}\\ 0& 1\end{array}\right]=\left[\begin{array}{cccc}c\alpha c\beta & c\alpha s\beta s\gamma -s\alpha c\gamma & c\alpha s\beta c\gamma +s\alpha s\gamma & x\\ s\alpha c\beta & s\alpha s\beta s\gamma +c\alpha c\gamma & s\alpha s\beta c\gamma -c\alpha s\gamma & y\\ -s\beta & c\beta s\gamma & c\beta c\gamma & z\\ 0& 0& 0& 1\end{array}\right]$

1

where $c\theta =\text{cos }\theta$ and $s\theta =\text{sin }\theta$ .

In the space coordinate system, the coordinates of any point on the moving coordinate system {B} can be expressed in the fixed coordinate system {A} by the homogeneous coordinate transformation matrix T

${}_{}{}^A\mathit{q}={}_B{}^A\mathit{T}{}_{}{}^B\mathit{q}$

2

Kinematic analysis

The generalized velocity of the moving platform of the parallel robot is expressed as $\dot{\mathit{q}}={\left[\begin{array}{cccccc}\dot{x}& \dot{y}& \dot{z}& {\omega }_x& {\omega }_y& {\omega }_z\end{array}\right]}^T$ . Euler angle ${\left[\begin{array}{ccc}\alpha & \beta & \gamma \end{array}\right]}^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 ω of the moving platform in the fixed coordinate system {A} and the Euler angle ${\left[\begin{array}{ccc}\alpha & \beta & \gamma \end{array}\right]}^T$ satisfy the following relationship

$\mathit{\omega }=\left[\begin{array}{ccc}0& -s\alpha & c\alpha c\beta \\ 0& c\alpha & s\alpha c\beta \\ 1& 0& -s\beta \end{array}\right]\left[\begin{array}{c}\dot{\alpha }\\ \dot{\beta }\\ \dot{\gamma }\end{array}\right]=\mathit{G}{\left[\begin{array}{ccc}\dot{\alpha }& \dot{\beta }& \dot{\gamma }\end{array}\right]}^T$

3

The forward speed of the parallel robot is

$\dot{\mathit{q}}=\mathit{J}{\dot{\mathit{L}}}_{i}i=1,\cdots ,6$

4

where ${\dot{\mathit{L}}}_i$ is the electric cylinder speed, and J is the first-order 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 ²⁴

$\ddot{\mathit{q}}=\mathit{J}{\ddot{\mathit{L}}}_i+{\dot{\mathit{L}}}_i{}^T{\mathit{H}}_o^q{\dot{\mathit{L}}}_{i}i=1,\cdots ,6$

5

where ${\ddot{\mathit{L}}}_i$ is the electric cylinder acceleration, ${\mathit{H}}_o^q$ is the second-order influence coefficient matrix solution of acceleration, and ${\mathit{H}}_o^q=-{\mathit{J}}^T\left(\mathit{J}*{\mathit{H}}_q^o\right)\mathit{J}\in {\mathit{R}}^{6\times 6\times 6}$ is the scalar matrix.

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

${\mathit{K}}_M=\frac{1}{2}{\dot{\mathit{q}}}^T{\mathit{M}}_M\left(\mathit{q}\right)\dot{\mathit{q}}$

6

${\mathit{M}}_{\text{M}}\left(\mathit{q}\right)=[\begin{array}{cccc}{\text{m}}_{\text{M}}& 0& 0& 0\\ 0& {\text{m}}_{\text{M}}& 0& 0\\ 0& 0& {\text{m}}_{\text{M}}& 0\\ 0& 0& 0& {\text{I}}_{\text{x}}{\text{s}}^2\text{β}+{\text{I}}_{\text{y}}{\text{s}}^2{\text{γc}}^2\text{β}+{\text{I}}_{\text{z}}{\text{c}}^2{\text{γc}}^2\text{β}\\ 0& 0& 0& \left({\text{I}}_{\text{y}}-{\text{I}}_{\text{Z}}\right)\text{cγsγcβ}\\ 0& 0& 0& -{\text{I}}_{\text{x}}\text{sβ}\end{array}\begin{array}{cc}0& 0\\ 0& 0\\ 0& 0\\ \left({\text{I}}_{\text{y}}-{\text{I}}_{\text{Z}}\right)\text{cγsγcβ}& -{\text{I}}_{\text{x}}\text{sβ}\\ {\text{I}}_{\text{y}}{\text{c}}^2\text{γ}+{\text{I}}_{\text{z}}{\text{s}}^2\text{γ}& 0\\ 0& {\text{I}}_{\text{x}}\end{array}]$

7

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 ₁ and m ₂ are the mass of the cylinder and the piston, respectively. l _i is the unit vector along the $\overrightarrow{A_i{B}_i}$ direction. $l_{gi}$ is the centroid position of the electric cylinder, and ${\mathit{v}}_{gi}$ is the centroid speed of the electric cylinder.

Figure 3.

Electric cylinder structure diagram.

The kinetic energy of a single electric cylinder is

${\mathit{K}}_{L_i}=\frac{1}{2}\left(m_1+m_2\right)\left[h_i{\mathit{v}}_{Bi}^T{\mathit{v}}_{Bi}-k_i{\mathit{v}}_{Bi}^T{\mathit{l}}_i{\mathit{l}}_i^T{v}_{Bi}\right]$

8

where $\widehat{l}=\frac{l_1{m}_1-l_2{m}_2}{m_1+m_2}$ , $h_i={\left(\frac{\widehat{l}}{L_i}+\frac{m_2}{m_1+m_2}\right)}^2$ , and $k_i=h_i-{\left(\frac{m_2}{m_1+m_2}\right)}^2$ .

The kinetic energy of six electric cylinders is

$$\begin{gathered} {\mathit{K}}_L=\sum _{i=1}^6{\mathit{K}}_{L_i} \\ =\frac{1}{2}\left(m_1+m_2\right)\left[{\left[\begin{array}{c}{\mathit{v}}_{B1}\\ {\mathit{v}}_{B2}\\ ⋮\\ {\mathit{v}}_{B6}\end{array}\right]}^T\mathit{H}\left[\begin{array}{c}{\mathit{v}}_{B1}\\ {\mathit{v}}_{B2}\\ ⋮\\ {\mathit{v}}_{B6}\end{array}\right]-{\left[\begin{array}{c}{\dot{L}}_1\\ {\dot{L}}_2\\ ⋮\\ {\dot{L}}_6\end{array}\right]}^T\mathit{K}\left[\begin{array}{c}{\dot{L}}_1\\ {\dot{L}}_2\\ ⋮\\ {\dot{L}}_6\end{array}\right]\right] \end{gathered}$$

9

where $\mathit{H}=\text{diag}\left\{\begin{array}{ccccccc}{h}_1,& h_1,& h_1,& \cdots ,& h_6,& h_6,& h_6\end{array}\right\}\in {\mathit{R}}^{18\times 18}$ and $\mathit{K}=\text{diag}\left\{\begin{array}{cccccc}{k}_1,& k_2,& k_3,& k_4,& k_5,& k_6\end{array}\right\}\in {\mathit{R}}^{6\times 6}$ .

Therefore, equation (9) can be expressed as

${\mathit{K}}_L=\frac{1}{2}{\dot{\mathit{q}}}^T{\mathit{M}}_L\left(\mathit{q}\right)\dot{\mathit{q}}$

10

Then, the inertia matrix of the electric cylinder can be obtained as follows

${\mathit{M}}_L=\left(m_1+m_2\right)\left[{\mathit{J}}_2^T\left[\mathit{H}-{\mathit{J}}_1^T\mathit{K}{\mathit{J}}_1\right]{\mathit{J}}_2\right]$

11

In the operating space, the Lagrangian dynamic equation can be expressed as

$\mathit{M}\left(\mathit{q}\right)\ddot{\mathit{q}}+\mathit{C}\left(\mathit{q},\dot{\mathit{q}}\right)\dot{\mathit{q}}+\mathit{G}\left(\mathit{q}\right)=\mathit{F}$

12

where q , $\dot{\mathit{q}}$ , and $\ddot{\mathit{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. $\mathit{M}\left(\mathit{q}\right)$ is the inertia matrix. $\mathit{C}\left(\mathit{q},\dot{\mathit{q}}\right)$ is Coriolis force and centripetal force matrix. $\mathit{G}\left(\mathit{q}\right)$ is the gravity matrix vector. $\mathit{M}\left(\mathit{q}\right)=M_{\mathit{M}}+{\mathit{M}}_L$ , $\mathit{C}\left(\mathit{q},\dot{\mathit{q}}\right)=C_M+{\mathit{C}}_L$ , and $\mathit{G}\left(\mathit{q}\right)={\mathit{G}}_M+{\mathit{G}}_L$ .

The Coriolis force/central force coefficient matrix $\mathit{C}\left(\mathit{q},\dot{\mathit{q}}\right)$ of the robot can be calculated ²⁷ by

$\mathit{C}\left(\mathit{q},\dot{\mathit{q}}\right)=\frac{1}{2}\left[\dot{\mathit{M}}\left(\mathit{q}\right)+{\mathit{U}}_m^T-{\mathit{U}}_m\right]$

13

Where ${\mathit{U}}_m=\left({\mathit{I}}_n\otimes {\dot{\mathit{q}}}^T\right) \left[\frac{\partial {\mathit{M}}^T}{\partial \mathit{q}}\right]$ and $\frac{\partial {\mathit{M}}^T}{\partial \mathit{q}}= \left[\begin{array}{cccc}\frac{\partial \mathit{M}}{\partial q_1}& \frac{\partial \mathit{M}}{\partial q_2}& \cdots & \frac{\partial M}{\partial q_n}\end{array}\right]$ .

The Coriolis force/central force coefficient matrix C_M of the moving platform is calculated according to equation (13) and is expressed as follows

${\mathit{C}}_M=\left[\begin{array}{cccccc}0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0\\ 0& 0& 0& c_{44}& c_{45}& c_{46}\\ 0& 0& 0& c_{54}& c_{55}& c_{56}\\ 0& 0& 0& c_{64}& c_{65}& 0\end{array}\right]$

14

where the parameters are detailed in the literature. ²⁶

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, ${\mathit{C}}_L=0$ .

The gravity of the moving platform is:

${\mathit{G}}_M=\frac{\partial {\mathit{P}}_M}{\partial \mathit{q}}={\left[\begin{array}{cccccc}0& 0& m_{M}g& 0& 0& 0\end{array}\right]}^T$

15

The gravity of the electric cylinder is

$\begin{array}{l}{\mathit{G}}_L=\left(m_1+m_2\right)g\sum _{i=1}^6\left[\frac{\widehat{\mathit{l}}}{L_i^2}{\left(\frac{\partial {\mathit{L}}_i}{\partial \mathit{q}}\right)}^T\right] \left(z_P-x_{B_i}^{\text{'}}s\beta +y_{B_i}^{\text{'}}c\beta s\gamma +z_{B_i}^{\text{'}}c\beta c\gamma \right)\\ +\left(m_1+m_2\right)g\sum _{i=1}^6\left[\frac{\widehat{l}}{L_i}+\frac{m_2}{m_1+m_2}\right]\left[\begin{array}{c}0\\ 0\\ 1\\ 0\\ -x_{B_i}^{\text{′}}c\beta -y_{B_i}^{\text{′}}s\beta s\gamma -z_{B_i}^{\text{′}}s\beta s\gamma \\ y_{B_i}^{\text{′}}c\beta c\gamma -z_{B_i}^{\text{′}}c\beta s\gamma \end{array}\right]\end{array}$

16

where $\frac{\partial {\mathit{L}}_i}{\partial \mathit{q}}$ is ith row of the inverse matrix of the Jacobian. $x_{B_i}^{\text{'}} y_{B_i}^{\text{'}} z_{B_i}^{\text{'}}$ indicates the position of the upper hinge in the moving coordinate system.

The three-dimensional (3-D) model of the fracture reduction robot is established in SolidWorks 2018 software. The related parameters of the model are measured in the software, as shown in Table 1.

Table 1.

Dynamics parameters of 6-UPU robot.

Parameter	Variable	Value	Unit
mM	Moving platform weight	5.292	kg
s Ix	Inertia to the x-axis	0.023	kg·m2
Iy	Inertia to the y-axis	0.023	kg·m2
Iz	Inertia to the z-axis	0.045	kg·m2
m 1	Cylinder weight	1.248	kg
m 2	Piston weight	0.221	kg
l 1	Cylinder center of gravity	93	mm
l 2	Piston center of gravity	45	mm
g	Gravity acceleration	9.8	m/s2

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

$\mathit{D}\left(\mathit{q}\right)\ddot{\mathit{L}}+\mathit{H}\left(\mathit{q},\dot{\mathit{q}}\right)\dot{\mathit{L}}+\mathit{P}\left(\mathit{q}\right)=\mathit{\tau }$

17

where L , $\dot{\mathit{L}}$ , and $\ddot{\mathit{L}}$ are the displacement, velocity, and acceleration of the electric cylinder of the parallel robot in the joint space, respectively. $\mathit{D}\left(\mathit{q}\right)$ , $\mathit{H}\left(\mathit{q},\dot{\mathit{q}}\right)$ , and $\mathit{P}\left(\mathit{q}\right)$ represent the inertia matrix, the Coriolis force and centripetal force matrix, and the gravity matrix vector in the joint space, respectively, and $\mathit{\tau }={\mathit{J}}^T\mathit{F}$ is the driving force.

$\mathit{D}\left(\mathit{q}\right)={\mathit{J}}^{\mathit{T}}\mathit{M}\left(\mathit{q}\right)\mathit{J}$

18

$\mathit{H}\left(\mathit{q},\dot{\mathit{q}}\right)={\mathit{J}}^T\mathit{C}\left(\mathit{q},\dot{\mathit{q}}\right)\mathit{J}+{\mathit{J}}^T\mathit{M}\left(\mathit{q}\right)\left({\dot{L}}_i{}^T{\mathit{H}}_o^q\right)$

19

$\mathit{P}\left(\mathit{q}\right)={\mathit{J}}^T\mathit{G}\left(\mathit{q}\right)$

20

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. ²⁸ Compared with other complex muscle force models, the Hill model is a phenomenon model, which is more reasonable for evaluating of overall muscle force. ²⁹ 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. ³⁰ Considering that the role of skin is very limited, so the effect of skin is ignored. The Hill model is shown in Figure 4.

Figure 4.

Hill skeletal muscle force model.

The mathematical expression of the Hill model is

$\{\begin{array}{c}{F}_{\text{ct}}=\left(F_{\text{me}}+F_{\text{pe}}\right)\cdot \text{cos}\left[\alpha \left(t\right)\right]\\ F_{\text{ct}}=F_{\text{se}} \\ \alpha \left((t)={\text{sin}}^{-1}\left(\frac{\text{sin}{\alpha }_0}{L}\right)\right) \end{array}$

21

where $F_{\text{me}}$ is the active contractile element of skeletal muscle. $F_{\text{pe}}$ is the passive elastic element of muscle connective tissue. $F_{\text{se}}$ is the series elastic element of tendons. $\alpha \left(t\right)$ is the pinnate angle between the muscle fiber and the tendon, and ${\alpha }_0$ is the pinnate angle of the optimal muscle fiber length. L is the muscle length factor

$L=\frac{L_c}{L_0}$

22

where L ₀ 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. ³¹ The expression of muscle force during reduction is

$\begin{array}{l}{F}_{\text{c}i}=1.3\cdot F_{0i}\cdot \text{arctan}\left[0.1\cdot {\left(\frac{L_{\text{c}i}}{L_{0i}}-0.22\right)}^{10}\right]\\ \cdot \text{cos}\left[\text{arcsin}\left(\frac{L_{0i}\cdot \text{sin}{\alpha }_{0i}}{L_{\text{c}i}}\right)\right]\end{array}$

23

where $L_{0i}$ , $F_{0i}$ , and ${\alpha }_{\text{0}i}$ 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 ${\mathit{F}}_{\text{c}}={\left[\begin{array}{cccccc}{F}_{\text{cx}}& F_{\text{cy}}& F_{\text{cz}}& M_{\text{cx}}& M_{\text{cy}}& M_{\text{cz}}\end{array}\right]}^T$ ∈ R ^{6 × 1} .

Figure 5.

Reduction force action diagram.

As the reduction force is the largest along the direction of the femoral shaft axis, so the reduction force can be simplified as ${\mathit{F}}_{\text{c}}={\left[\begin{array}{cccccc}{F}_{\text{cx}}& 0& 0& 0& 0& 0\end{array}\right]}^T$ .

The maximum reduction force of the femoral shaft fracture is 396 N. Setting the feather angle between muscle fibers and tendons, α(t) is 0. ²⁸ The optimal muscle fiber length L ₀ is 90 mm. ²⁹ 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 x-axis 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.

Figure 6.

The relationship between axial reduction force and muscle stretch length.

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

$\mathit{D}\left(\mathit{q}\right)\ddot{\mathit{L}}+\mathit{H}\left(\mathit{q},\dot{\mathit{q}}\right)\dot{\mathit{L}}+\mathit{P}\left(\mathit{q}\right)=\mathit{\tau }+{\mathit{\tau }}_{\text{c}}$

24

where ${\mathit{\tau }}_{\text{c}}$ is the equivalent moment of the reduction force in the joint space.

The reduction force ${\mathit{F}}_{\text{c}}$ acts on the distal end of the fracture and its equivalent moment ${\mathit{\tau }}_{\text{c}}$ can be transformed by the Jacobian matrix of the robot as follows

${\mathit{\tau }}_{\text{c}}={\mathit{J}}_{\text{c}}^T\left(\mathit{q}\right){\mathit{F}}_{\text{c}}$

25

where ${\mathit{J}}_{\text{c}}\left(\mathit{q}\right)$ is the Jacobian matrix of the robot.

Defining the generalized momentum of the reduction robot is as

$\mathit{p}=\mathit{D}\left(\mathit{q}\right)\dot{\mathit{L}}$

26

Doing the partial deviation of equation (26)

$\dot{\mathit{p}}=\mathit{D}\left(\mathit{q}\right)\ddot{\mathit{L}}+\dot{\mathit{D}}\left(\mathit{q}\right)\dot{\mathit{L}}$

27

$\dot{\mathit{D}}\left(\mathit{q}\right)-2\mathit{H}\left(\mathit{q},\dot{\mathit{q}}\right)\in R^{n\times n}$ is an anti-symmetric matrix, and $\mathit{D}\left(\mathit{q}\right)$ is a symmetric and positive definite matrix. ²⁶

Therefore,

$\dot{\mathit{D}}\left(\mathit{q}\right)=\mathit{H}\left(\mathit{q},\dot{\mathit{q}}\right)+{\mathit{H}}^T\left(\mathit{q},\dot{\mathit{q}}\right)$

28

From equations (24), (27), and (28), the following equation can be obtained

$\dot{\mathit{p}}={\mathit{H}}^T\left(\mathit{q},\dot{\mathit{q}}\right)\dot{\mathit{q}}+\mathit{\tau }-\mathit{P}\left(\mathit{q}\right)+{\mathit{\tau }}_{\text{c}}$

29

Since equation (29) shows that $\dot{\mathit{p}}$ contains the external moment component ${\mathit{\tau }}_{\text{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

$\mathit{r}\left(t\right)={\mathit{K}}_1{\int }_0^t\left[-{\mathit{K}}_2\mathit{r}+\left(\mathit{p}-\widehat{\mathit{p}}\right)\right]dt$

30

where ${\mathit{K}}_1>0$ , ${\mathit{K}}_2>0$ is the gain matrix. $\widehat{\dot{\mathit{p}}}$ is the momentum estimate value. $\widehat{\dot{\mathit{p}}}=\mathit{H}{\left(\mathit{q},\dot{\mathit{q}}\right)}^T\dot{\mathit{L}}+\mathit{\tau }-\mathit{P}\left(\mathit{q}\right)+\mathit{r}$ . Equation (30) can be rewritten as

$\mathit{r}\left(t\right)={\mathit{K}}_1{\int }_0^t[-{\mathit{K}}_2\mathit{r}+\mathit{p}\left(t\right)-{\int }_0^t\left(\mathit{\tau }+{\mathit{H}}^T\left(\mathit{q},\dot{\mathit{q}}\right)\dot{\mathit{L}}-\mathit{P}\left(\mathit{q}\right)+\mathit{r}\right)dt]dt$

31

When only the equation above is used as the observer, a large delay and oscillation are generally generated. ³³ So constructing an adjustment function is as follows

${\mathit{f}}_c=\mathit{p}\left(\mathit{t}\right)-{\int }_0^{\mathit{t}}\left(\mathit{\tau }+H^T\left(\mathit{q},\dot{\mathit{q}}\right)\dot{\mathit{L}}-\mathit{P}\left(\mathit{q}\right)+\mathit{r}\right)dt$

32

Taking f_c as the feedforward adjustment of equation (31), the effective residual observer is

$\begin{array}{l}\mathit{r}\left(t\right)={\mathit{K}}_1{{\displaystyle \int }}_0^t[-{\mathit{K}}_{2}r+\mathit{p}\left(t\right)-\\ {{\displaystyle \int }}_0^{\mathit{t}}\left(\mathit{\tau }+{\mathit{H}}^T\left(\mathit{q},\dot{\mathit{q}}\right)\dot{\mathit{L}}-\mathit{P}\left(\mathit{q}\right)+\mathit{r}\right)\text{d}t]\text{d}t+{\mathit{K}}_1{\mathit{K}}_3{\mathit{f}}_c\end{array}$

33

where ${\mathit{K}}_3>0$ is the gain matrix.

The structure diagram of the force residual observer is designed, as shown in Figure 7.

Figure 7.

Force residual observer block diagram.

Find the second-order partial derivative for equation (33) and bring the result into equation (29)

$\ddot{\mathit{r}}=-{\mathit{K}}_1{\mathit{K}}_2\dot{\mathit{r}}-{\mathit{K}}_1\mathit{r}+{\mathit{K}}_1{\mathit{\tau }}_{\text{c}}+{\mathit{K}}_1{\mathit{K}}_3{\dot{\mathit{\tau }}}_{\text{c}}$

34

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

${\mathit{H}}_{\text{r}}\left(s\right)=\frac{{\mathit{r}}_i\left(s\right)}{{\mathit{\tau }}_{\text{c},i}\left(s\right)}=\frac{{\mathit{K}}_{3,i}s+1}{\frac{1}{{\mathit{K}}_{1,i}}{s}^2+{\mathit{K}}_{2,i}s+1}i=1,\cdots ,6$

35

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 ${\mathit{\tau }}_{\text{c}}$ . The system can be stable by changing the observer gains K ₁, K ₂, and K ₃.

${\mathit{r}}_i\approx {\mathit{\tau }}_{\text{c},i}={\mathit{J}}_i^T\left(\mathit{q}\right){\mathit{F}}_{\text{c},i}i=1,\cdots ,6$

36

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_{\text{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.

Figure 8.

Simulation system with force residual observer.

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.

Figure 9.

Marking point trajectory.

Selecting a feature point on the distal fracture as the marking point $P_{\text{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_{\text{c}}$ are (25.94, 32.96, and 2.72 mm).

Simulation without reduction force

When the reduction force is not considered, the robot-assisted fracture reduction is simulated. The actual movement trajectory of the marked point $P_{\text{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_{\text{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 robot-assisted femoral fracture is feasible.

Figure 10.

Motion error of marked point.

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.

Figure 11.

Driving force of each electric cylinder.

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.

Figure 12.

Driving force of electric cylinder expressed in operating space.

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.

Figure 13.

External force value by the force residual observer.

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_{\text{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.

Figure 14.

Motion error of distal segment with acting the reduction force.

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.

Figure 15.

Driving force of each electric cylinder with acting the reduction force.

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.

Figure 16.

Driving force of electric cylinder in operating space with acting the reduction force.

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.

Figure 17.

(a) External force by force residual observer. (b) Torque by the force residual observer.

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 robot-assisted fracture reduction.