Use of a Force-Torque Sensor for Self-Calibration of a 6-DOF Medical Robot

The aim of this paper is to improve the position accuracy of a six degree of freedom medical robot. The improvement in accuracy is achieved without the use of any external measurement device. Instead, this work presents a novel calibration approach based on using an embedded force-torque sensor to identify the robot’s kinematic parameters and thereby enhance the positioning accuracy. A simulation study demonstrated that our calibration approach is effective, whether or not any measurement noise is present: the position error is improved, inside the robot target workspace, from 12 mm to 0.320 mm, for the maximum values, and from 9 mm to 0.2771 mm, for the mean errors.


Introduction
Medical robots show a promising future in various health issues in the most recent decades. With recent developments in sensors and control theory, medical robots provide many inspiring solutions in the fields of: diagnosis, surgery, orthopedics, rehabilitation, prosthetics and exoskeletons, etc. [1,2]. The force-torque (wrench) sensor is an essential component of these medical robot applications. MedRUE [3], OTELO [4] and Hippocrate [5] robot systems were developed for the ultrasound scanning of vascular diseases. In all three robot systems, force-torque sensors are employed to maintain proper contact with the patient's body during the examination. The Black Falcon system, which is a fundamental study for many other surgical robot systems, allows the surgeon to feel the interaction with tissue and thereafter improve the surgical performance. The Da Vinci system is a popular surgery robot widely used in hospitals [6]. The adoption of force-torque sensors in Da Vinci system has been studied in depth [7]. The Robodoc assistant system is a medical orthopedic robot for use during total knee replacement, and it achieves results that are comparable to technician performance [8]. The force-torque sensors are used in Robodoc for both control and safety reasons.
Medical robots have also contributed to the field of rehabilitation. The InMotion ARM, which is based on the MIT-Manus project, is an interactive robotic system for upper-extremity rehabilitation therapy [9,10], and the robotic stepper is a device, developed by the National Aeronautics and Space Administration (NASA), to help patients with lower-extremity rehabilitation [11]. Force-torque sensors are employed in these rehabilitation robots to measure the strength and the capability of the patient. Medical robots have also been developed as substitutes for malfunctioning parts of the human body. The I-limb ultrasound system and the ReWalk system are exoskeleton robots for hand prosthetics and leg prosthetics, respectively [12,13]. Force-torque sensors are used in the prosthetic and exoskeleton robots to control the joints and to evaluate the power of the limb movements.
Medical robots often need to be accurate, not just repeatable, which means that they must be calibrated. Most robot calibration approaches are based on minimizing the pose residual, which

Robot Description and the Main Reference Frames
The MedRUE robot ( Figure 1a) is a medical robot dedicated for vascular ultrasound examination. MedRUE is a six degrees of freedom (6-DOF) hybrid serial-parallel robot. It is composed of two five-bar mechanisms (Figure 1b), which are symmetrically assembled. These mechanisms are considered to be perfectly parallel to each other and perpendicular to the robot base. The robot base is fixed to a linear guide actuated through a servomotor SM 1 ; the corresponding joint variable is denoted by q 1 . Five other servomotors are also used in order to actuate the robot revolute joints: SM 2 and SM 3 , for the left five bar mechanism, and SM 4 and SM 5 for the right side (Figure 1a). The sixth servomotor-attached Sensors 2016, 16, 798 3 of 19 to the link having L 14 as its length-is used to rotate the probe about the x axis of the last reference frame (F 6 ), which is defined as follows: the origin of F 6 is located midway between G 1 and G 2 ; its x axis (x 6 ) is defined to pass through G 1 and G 2 , and z 6 is pointing toward the probe center.
Sensors 2016, 16,798 3 of 20 base. The robot base is fixed to a linear guide actuated through a servomotor SM1; the corresponding joint variable is denoted by q1. Five other servomotors are also used in order to actuate the robot revolute joints: SM2 and SM3, for the left five bar mechanism, and SM4 and SM5 for the right side ( Figure 1a). The sixth servomotor-attached to the link having 14 L as its length-is used to rotate the probe about the x axis of the last reference frame (F6), which is defined as follows: the origin of F6 is located midway between G1 and G2; its x axis (x6) is defined to pass through G1 and G2, and z6 is pointing toward the probe center. Each five-bar mechanism i (i = 1, 2) has five links: the distance di between the anchor points of the two proximal links, and the four mobile links having Lij (j = 1, 2, 3, 4) as lengths. The five links connect five revolute joints (Ai, Bi, Ci, Di, Ei), among which only two (Ai and Ci) are actuated through servomotors SMi×2 and SMi×2+1: the corresponding two angles are denoted qi×2 and qi×2+1, respectively. A total of five angles of active joints are considered (q2, …, q6). Each five-bar mechanism i (i = 1, 2) has five links: the distance d i between the anchor points of the two proximal links, and the four mobile links having L ij (j = 1, 2, 3, 4) as lengths. The five links connect five revolute joints (A i , B i , C i , D i , E i ), among which only two (A i and C i ) are actuated through servomotors SM iˆ2 and SM iˆ2+1 : the corresponding two angles are denoted q iˆ2 and q iˆ2+1 , respectively. A total of five angles of active joints are considered (q 2 , . . . , q 6 ).
Finally, joints E 1 and E 2 are linked through the probe support (Figure 1a), which has a universal joint at each extremity G i . The x coordinate of G i with respect to the base frame is denoted by d 4i . In our calibration process, nine reference frames are considered: The reference frame of the robot base, located on the robot base at O 0 . As shown in Figure 1a, the x axis (x 0 ) is aligned with the axis of the linear guide, and z 0 is normal to the plane defined by the platform of the robot base. The translation T world 0 " " x 0 y 0 z 0 ı T and the orientation (α 0 , β 0 , γ 0 ), described in XYZ fixed Euler angles, of F 0 with respect to F world , are expected to be identified by the calibration process.

‚
F world : The world reference frame (Figure 1a), associated with the robot work-cell. It has approximately the same orientation as F 0 .
‚ F 1 to F 6 : The reference frames associated with the joints (F 1 , F 2 , . . . , F 6 ). These frames are not shown.
‚ F tool : The tool reference frame associated with the robot probe ( Figure 2). The origin of F tool is described to be the center of the end-effector (i.e., the probe), and its orientation is considered to be the same as that of F 6 . Knowing that the end-effector orientation is not used in our calibration process, therefore, only the translation T 6 tool " " x t y t z t ı T of F tool with respect to F 6 is considered. Finally, joints E1 and E2 are linked through the probe support (Figure 1a), which has a universal joint at each extremity Gi. The x coordinate of Gi with respect to the base frame is denoted by 4i d .
In our calibration process, nine reference frames are considered: The reference frame of the robot base, located on the robot base at O0. As shown in Figure 1a, the x axis (x0) is aligned with the axis of the linear guide, and z0 is normal to the plane defined by the platform of the robot base. The translation T world x y z  T and the orientation (α0, β0, γ0), described in XYZ fixed Euler angles, of F0 with respect to Fworld, are expected to be identified by the calibration process.  Fworld: The world reference frame (Figure 1a), associated with the robot work-cell. It has approximately the same orientation as F0.  F1 to F6: The reference frames associated with the joints (F1, F2, …, F6). These frames are not shown.  Ftool: The tool reference frame associated with the robot probe ( Figure 2). The origin of Ftool is described to be the center of the end-effector (i.e., the probe), and its orientation is considered to be the same as that of F6. Knowing that the end-effector orientation is not used in our calibration process, therefore, only the translation  

Position Equations
Given the vector ψ = [q1, q2, …, q6] T of the active joint variables, the end-effector's pose with respect to the world frame is represented by homogeneous matrices as follows:

Position Equations
Given the vector ψ = [q 1 , q 2 , . . . , q 6 ] T of the active joint variables, the end-effector's pose with respect to the world frame is represented by homogeneous matrices as follows: where A b a denotes the homogeneous matrix representing a frame a with respect to a frame b, and can be represented in the matrix form of The rotation matrix is often represented by R(α,β,γ) = R x (γ)R y (β)R z (α), and the translation matrix by T(x,y,z) = T x (x)T y (y)T z (z), where R u (φ) and T u (d) are the rotation/translation operator along the u axis with value φ/d. Since robot joints are only included in A 0 wrist , then A world base , A wrist sensor and A sensor tool in Equation (1) are constant matrices, which can be defined directly by parameter sets [x 0 , y 0 , z 0 , α 0 , β 0 , γ 0 ], [x S , y S , z S , α S , β S , γ S ] and [x T , y T , z T , α T , β T , γ T ]. The following paragraphs present the calculation of A 0 wrist . The coordinates of B i and D i are expressed with respect to the local frame F i on ith five-bar mechanism as follows: where L i1 and L i3 are the lengths of the four swinging links as shown in Figure 1, and δq i is the offset of ith active joint. The vector r O i A i is calculated as follows: and The coordinates of E i with respect to a frame F i are obtained as follows: where The coordinates of E i w.r.t. F base are obtained by a transformation matrix A 0 i as follows: where and θ i , which is the angle between r A i C i and the normal of the x 0 y 0 plane, is calculated as follows: The orientation of F wrist w.r.t. F base is obtained from the corresponding rotation matrix R base wrist pα, β, γq, where α, β and γ are the fixed XYZ Euler angles. The rotation angle along x 0 is directly obtained as γ " q DE 1`q 6`δ q 6 . According to the design of the tool part shown in Figure 2, Then α and β are obtained as: where, Ψ " u y sinγ´u z cosγ. The translation of F wrist w.r.t. F base is calculated as follows: Finally, the pose of F wrist w.r.t. F base is expressed as follows:

Force and Torque Equations
The gravity frame F gravity is assigned at the tool part's center of gravity, as shown in Figure 3. When the robot is not in contact with its environment, the gravity force f G in frame F gravity is the cause of the force and torque on the force sensor. The forward kinematic solution to obtain the force and torque in force sensors is: where F is a 6ˆ1 wrench vector composed of force and torque, and A is a 6ˆ6 transformation matrix between screws. The orientation of F gravity is in alignment with F world , rather than fixed relative to the tool part. The wrench of the gravity force of the tool part w.r.t. the F gravity is with m Tool is the mass of the tool part and g is the gravitational constant. Since gravity is a pure force, τ The transformation matrix A writst gravity can be expressed as Since F gravity keeps the same orientation with F world , then R wrist gravity " R wrist world "`R world base R base wrist˘T .
p wrist gravity " " x G y G z G ı T is the translation offset of the origin of F gravity w.r.t. F wrist . p wrist gravityˆi s the vector product operation, and it is equal to Similar to Equation (19), A sensor wrist is characterized by parameters describing F sensor w.r.t. F wrist : The sensor reference frame F sensor is defined in the information given by the sensor manufacturer. During assembly, its orientation w.r.t. F wrist is expressed by Euler-XYZ angles: Similar to Equation (20), p wrist sensorˆi s the vector product operation of p wrist sensor " with mTool is the mass of the tool part and g is the gravitational constant. Since gravity is a pure force, The transformation matrix  writst gravity can be expressed as Similar to Equation (19), sensor wrist  is characterized by parameters describing Fsensor w.r.t. Fwrist:

Parameters Used During Calibration
In our robot calibration process, the following parameters were considered:

‚
The lengths of the ten links of the two five-bar mechanisms: L 11 , L 12 , L 13 , L 14 , L 21 , L 22 , L 23 and L 24 .

‚
The y and z coordinates of the anchor points of the two proximal links of the five-bar mechanisms: The offsets of the six active joints: δq 1 , δq 2 , δq 3 , δq 4 , δq 5 , δq 6 .

‚
The position of the tool frame with respect to the wrist frame: x T , y T , z T , α T , β T , γ T .

‚
The parameters to describe the sensor frame w.r.t. the wrist frame: x S , y S , z S , α S , β S , γ S . ‚ The parameters to describe the offset of gravity frame w.r.t. the wrist frame: x G , y G , z G .

‚
The mass of the tool part: m Tool .
Of the 46 parameters that we considerd, a total of 17 parameters are non-identifiable, which means that we need to reduce the number of identifiable parameters to 29.

Calibration Process
Our calibration process is explained in detail in Sections 4.1-4.3. Its main steps are presented in what follows:

1.
Develop the calibration model: the forward kinematics, presented in Section 2.2.

2.
Create a pool Ω of 40,000 configurations uniformly distributed inside the whole robot workspace. Create a set Ω t of 336 configurations uniformly distributed inside the target workspace (see Section 4.3). We note that the configurations of the set Ω t are different from these of Ω.

3.
Select 100 configurations to be used in the identification process. These configurations are chosen through an observability analysis, as explained in Section 4.1.

4.
Take the force and torque measurements, for all robot configurations (Ω t and Ω). Measurements are done by using the robot force-torque sensor. We note that in this paper all measurements are generated by simulation, as explained in Sections 4.3 and 5.

5.
Identify the robot parameter values by using the calibration configurations selected in step 3; the identification approach is presented in details in Section 4.2. 6.
Evaluate the accuracy after calibration, as explained in Section 5.

Selection of Calibration Configurations
After creating the calibration model, and generating a pool Ω of 40,000 configurations uniformly distributed inside the whole robot workspace, a set of 100 calibration configurations is selected among Ω. This is done by using an approach commonly called observability analysis. This analysis is used to obtain the optimal set of the calibration configurations, and is based on the singular value decomposition (SVD) of the identification Jacobian matrix J. The matrix J is composed of the derivatives of the end-effector force and torque vector (Equation (17)), with respect to all of the robot independent parameters. The Jacobian matrix is also used in the linearization of the force and torque equations (Equation (17)), around the calibration configurations (i.e., Taylor approximation). This linearization allows identifying the parameter values, as explained in Section 4.2. The nominal values of the robot's independent parameters are represented by the vector p nom . The matrix J is calculated as follows, for i = 1 . . . n.
where J i is the 6ˆm Jacobian matrix at the ith calibration configuration, n is the number of calibration configurations, and m is the number of considered parameters (not all of which are necessarily identifiable). In our case, n = 100 and m = 49. J i is given by: The matrix J is also used to find the non-identifiable parameters. The rank r J of the Jacobian J represents the number of identifiable parameters. If r J < m, then m´r J parameters are non-identifiable; the corresponding columns should be removed from J. This procedure is carried out using an algorithm that is based on the approach proposed in [21]. The stop criteria is when m becomes equal to r J . The algorithm proceeds as follows: 1.
Remove all zero columns from J. The corresponding parameters have no impact on the calibration model.

2.
Calculate the condition number, c J , of J. The condition number is used to evaluate how good is the matrix J for the parameter identification. With a bad condition number (high value), the solutions are unstable with respect to small changes in measurement errors. Therefore, to have a robust identification system, the condition number should be as small as possible.

3.
Remove, one at a time, the column related to each parameter from J, and calculate both the new rank and condition number (rJ and cJ ) for the new Jacobian matrix J * . The column that, if eliminated, results in the maximum reduction of the condition number and gives the same rank (rJ = r J ), is definitively removed (i.e., the corresponding parameter will be not subject to the identification process).

4.
Replace J with J*, and repeat the process from step (2).
As a result, of the 49 parameters considered in our calibration model, 22 are non-identifiable and are indicated by the symbol '*' in Table 1.
The fact that some parameters are non-identifiable is mainly due to redundancy, or the fact that they have no impact on the force and torque equations that represent the calibration model.
The parameter identification is achieved by minimizing the residual of the end-effector force and torque, which are measured by the robot's force sensor ( Figure 3); no external measurement device is required. Further, only the gravity effect of the end-effector is used to apply forces to the robot's end-effector. Therefore, to change the applied force on the end-effector, it is necessary to change its orientation.
In our identification process, 100 calibration configurations are selected from among the 40,000 configurations. The remaining 39,900 configurations are used for validation purpose. The 40,000 configurations are uniformly distributed on three layers, within the entire robot workspace. Note that several orientations are generated for each position. The calibration configurations are selected using an observability analysis, which allows us to identify the most appropriate configurations and thus identify the most effective parameters. This analysis is based on using the first observability index, denoted by O 1 and calculated by using the singular value of the Jacobian identification matrix (i.e., the sensitivity matrix). The procedure of selecting the calibration configurations is based on the DETMAX algorithm, which was initially proposed in [24].
According to [25,26], the index O 1 seems to be the most appropriate index for the kinematic calibration. This was also confirmed by our simulation, through a comparison of the five observability indices that were presented in the literature and thoroughly detailed in [26]. The convergence of O 1 is represented in Figure 4, and is calculated as follows: where n is the number of calibration configurations, σ 1 . . . σ m are the singular values of the Jacobian identification matrix for the m = 29 identifiable parameters. In our identification process, 100 calibration configurations are selected from among the 40,000 configurations. The remaining 39,900 configurations are used for validation purpose. The 40,000 configurations are uniformly distributed on three layers, within the entire robot workspace. Note that several orientations are generated for each position. The calibration configurations are selected using an observability analysis, which allows us to identify the most appropriate configurations and thus identify the most effective parameters. This analysis is based on using the first observability index, denoted by O1 and calculated by using the singular value of the Jacobian identification matrix (i.e., the sensitivity matrix). The procedure of selecting the calibration configurations is based on the DETMAX algorithm, which was initially proposed in [24].
According to [25,26], the index O1 seems to be the most appropriate index for the kinematic calibration. This was also confirmed by our simulation, through a comparison of the five observability indices that were presented in the literature and thoroughly detailed in [26]. The convergence of O1 is represented in Figure 4, and is calculated as follows: where n is the number of calibration configurations, σ1 … σm are the singular values of the Jacobian identification matrix for the m = 29 identifiable parameters.

Parameter Identification Process
The parameter identification process is based on using the Jacobian matrix J, which relates the force and torque errors to the 29 unknown parameter values. The matrix J is built by the linearization of the forward kinematics model (Equation (17)) around each calibration configuration.

Parameter Identification Process
The parameter identification process is based on using the Jacobian matrix J, which relates the force and torque errors to the 29 unknown parameter values. The matrix J is built by the linearization of the forward kinematics model (Equation (17)) around each calibration configuration. The parameter values are identified by means of an iterative algorithm, in which the parameters' vector is initialized by p nom , and is updated at each iteration (i.e., replaced by the vector p identified of the identified values). The matrix J is also iteratively updated, since its calculation is based on p.
The identification algorithm is presented in Figure 5 and has the following steps: A system of linear equations is formed by the measured force and torque errors, the unknown robot's parameter errors, and the Jacobian matrix J. In order to maintain acceptable variance of each parameter (i.e., proper convergence in the linear system), parameter scaling is implemented, by using the column scaling approach proposed in [24]. The scaled matrix obtained is denoted by J scal , and it is used to identify the robot's scaled parameter errors (∆ scal ), as follows:  (17) The parameter errors, which represent the difference between the real values and the nominal values of the parameters, are denoted by ∆, and calculated as follows: where D j , (j = 1, 2, . . . , m) are the scaling coefficients, defined as follows: is the number of calibration configurations, and J ij is the element of the Jacobian matrix located at the ith row and the jth column.
Finally, the vector of the identified parameter values is To converge towards a solution for the unknown parameter values, an iterative Newton-based procedure was used. After p identified has been calculated, the p vector is replaced by the last p identified vector obtained, and the estimation process is restarted from step (a).
Steps (a) to (d) are repeated until reaching a convergence criterion, which is the root mean square error (RMSE) between two successive iterations. The RMSE is evaluated between the vector of the latest identified parameters and the previous one. The convergence criterion was set to 10´1 6 , and the system converged towards a solution after five iterations.

Validation after Calibration
After calibration, the accuracy is validated using 336 configurations that are uniformly distributed inside the Cartesian target workspace. The target workspace is intended to correspond to the area where the patient's leg will be located ( Figure 6b). Also, the accuracy after calibration is evaluated by using the 39,900 configurations (denoted by Ωw), which are the remaining configurations among the initial set Ω composed of 40,000 configurations (Figure 6a) uniformly distributed within the whole robot workspace: 100 calibration configurations are selected from Ω, through the observability analysis, to be used in the parameter identification process, and 39,900 configurations are used in the validation after calibration.
After achieving the parameter identification process, the identified parameters are used in the robot kinematics, instead of the nominal parameter values. The next step consists to evaluate the robot accuracy for all validation configurations (Ωw and Ωt), by using the following algorithm:

Loop 1
For each validation set, Ωw and Ωt:

Loop 2
For each validation configuration:

Validation after Calibration
After calibration, the accuracy is validated using 336 configurations that are uniformly distributed inside the Cartesian target workspace. The target workspace is intended to correspond to the area where the patient's leg will be located ( Figure 6b). Also, the accuracy after calibration is evaluated by using the 39,900 configurations (denoted by Ω w ), which are the remaining configurations among the initial set Ω composed of 40,000 configurations (Figure 6a) uniformly distributed within the whole robot workspace: 100 calibration configurations are selected from Ω, through the observability analysis, to be used in the parameter identification process, and 39,900 configurations are used in the validation after calibration.
After achieving the parameter identification process, the identified parameters are used in the robot kinematics, instead of the nominal parameter values. The next step consists to evaluate the robot accuracy for all validation configurations (Ω w and Ω t ), by using the following algorithm:

Loop 1
For each validation set, Ω w and Ω t :

Loop 2
For each validation configuration: Calculate the composed error ( a Ex 2`E y 2`E z 2 ) by using results obtained in step (c).

End Loop 2
Calculate the mean, the maximum and the standard deviation of all composed errors obtained in Loop 1.

End Loop 1
The force and torque validation is achieved by using the same algorithm as for position. The only difference is using Equation (17)  Ex Ey Ez   ) by using results obtained in step (c).

End Loop 2
Calculate the mean, the maximum and the standard deviation of all composed errors obtained in Loop 1.

End Loop 1
The force and torque validation is achieved by using the same algorithm as for position. The only difference is using Equation (17)

Simulation Study
The efficiency of our calibration process is evaluated through a simulation. This process also evaluates the sensitivity of our identification process to the measurement noise, and verifies the effectiveness of the observability analysis for choosing the calibration configurations. Finally, the

Simulation Study
The efficiency of our calibration process is evaluated through a simulation. This process also evaluates the sensitivity of our identification process to the measurement noise, and verifies the effectiveness of the observability analysis for choosing the calibration configurations. Finally, the calibration results (i.e., position accuracy) that were obtained from each of the five observability indices are compared, and the index that gives the best accuracy is used in the actual calibration.
For simulation purposes, the actual parameters' values are simulated by introducing randomly-generated errors of˘2 mm for the distances, and˘1˝for the angles. The differences between the nominal and the actual parameter errors simulate the behavior of a robot with poor accuracy. By using the calibration process, the identified parameters will be as close as possible to their actual values, despite the presence of the measurement errors. The measurement errors that were used in our simulation are˘1 N and˘0.2 N¨m, for the force and torque, respectively. These errors are an exaggeration of the accuracy of the robot's force-torque sensor (a Mini 40 from ATI), the details of which are provided by its manufacturer, through a calibration certificate. From this information, force accuracy according to x, y, and z axes was˘0.25 N,˘0.2 N, and˘0.45 N, respectively. The torque accuracy was˘0.0125 N.m,˘0.0125 N.m, and˘0.02 N.m, for x, y, and z axes.
The measurement errors are generated according to a normal distribution, for each axis (i.e., errors for F x are generated within˘1 N, and similarly for F y and F z ). The data acquisition is simulated by generating 100 measurements (i.e., force and torque errors) for each calibration configuration of the robot. As it is known that the number of identifiable parameters is 29, the number of calibration configurations that are used in the identification process is 100, in order to over-constrain the calibration model.
The measured wrench vector (composed of force and torque) is simulated for each calibration configuration, by substituting the corresponding active joint angles and the actual parameter values (Table 1) in Equation (17). A vector of measurement noise is then added to the obtained wrench vector.
Once all force-torque vectors are generated, the robot parameters are identified as explained in Section 4.2. The identified parameter values are presented in Table 1.
Once the parameters have been identified, the calibration process is validated. This validation is carried out, as explained in Section 4.3, on two levels: the robot's position accuracy is assessed in both the whole robot workspace (39,900 configurations) and by using the 336 positions that are uniformly distributed within the target workspace. The force, torque and position errors are summarized in Table 2 and Table 3, and it shows that the position accuracy was improved from 8.9135 mm before calibration to 286 µm, after calibration, inside the target workspace. The wrench errors (Table 4) were highly improved (better than the position accuracy improvement) because in our identification process, only the residuals of force and torque were minimized in the objective function Equation (26). The distribution of the robot's position errors (before and after calibration) is presented in Figure 7 and Figure 8, which represent the number of occurrences (frequency) of robot xyz composed position errors that lie within the ranges of error, presented on the horizontal axis.
A deep analysis of the position errors, for the whole workspace, shows that 0.33% of the evaluated positions have an accuracy lower than 0.2024 mm (mean´2ˆSTD), 94.37% are within the range mean 2ˆSTD, and only 5.3% of positions present the poorer accuracy, which is higher than 0.3696 (mean + 2ˆSTD). The same analysis was achieved for the target workspace, and it shows that 93.0952% of positions have an accuracy within 0.2457 mm and 0.3085 mm (mean˘2ˆSTD), and only 6.9048% of positions have errors higher than 0.3085 mm (mean + 2ˆSTD).
For illustrative purpose, Table 4 shows the accuracy obtained by using each observability index, separately, in the calibration process. Results confirm that O 1 is the most appropriate index for calibrating our robot, since it gives the smallest position (and force/torque) errors, after calibration. Moreover, deeper statistical analyses were carried out on the results obtained by the five indices. First we verified whether the data distributions are Gaussian or not, and then used parametric or non-parametric tests accordingly. Therefore a Kolmogorov-Smirnov test (not shown) was achieved, and it showed that all observability indices provide Gaussian distribution. Based on these results, we decided to use parametric analyses (i.e., ANOVA analysis, and t-test).
For illustrative purpose, Table 4 shows the accuracy obtained by using each observability index, separately, in the calibration process. Results confirm that O1 is the most appropriate index for calibrating our robot, since it gives the smallest position (and force/torque) errors, after calibration. Moreover, deeper statistical analyses were carried out on the results obtained by the five indices. First we verified whether the data distributions are Gaussian or not, and then used parametric or non-parametric tests accordingly. Therefore a Kolmogorov-Smirnov test (not shown) was achieved, and it showed that all observability indices provide Gaussian distribution. Based on these results, we decided to use parametric analyses (i.e., ANOVA analysis, and t-test).    Initially, an ANOVA analysis, with a probability threshold α = 0.05, is used to confirm the objectivity of comparing the five indices (i.e., confirm that there is actually differences between the use of the five indices). Results show that the F value is significantly higher than the F criteria, which leads to reject the null hypothesis, and therefore conclude that the comparison of results (position accuracy) obtained by using the five indices is meaningful (i.e., results are different, and some indices are better than others).
The ANOVA does not tell where the difference between indices lies. Therefore, an additional test is considered (t-test). The t-test is used to compare each pair of indices. However, the position errors obtained by using O 3 and O 4 were clearly poorer than results obtained by the other indices (O 1 , O 2 , and O 5 ). Thus, only O 1 , O 2 , and O 5 are considered in the t-Test. The results of this test are shown in Table 5, and they show that the position errors obtained by these three indices are quite different, since the t Stat value is significantly lower than -t_Critical_two-tail, for all pairs of comparisons. Furthermore, to take into account multiple comparison effects, a post-hoc correction is included (Bonferroni correction). As summarized in Table 5, results show that all the t-tests were statistically significant (i.e., there are significant differences between the performances of the observability indices). Based on the aboves tests, we conclude that the fives observability indices allow different results, regarding the robot accuracy after calibration. The analysis of the mean and maximum errors (Table 4) shows that the index O 1 leads to the best robot accuracy: it gives not only the smallest mean errors, but also the smallest maximum errors. Also, O 1 has a small standard deviation, which means that position errors are closely distributed around the mean value.
We recall that in our simulation the used measurement noise was˘1 N. This range of error is an exaggerated error of our wrench sensor (Mini 40, from API). For illustrative propose, we achieved other simulations by considering lower measurement noise (Table 6). Results demonstrate that the accuracy after calibration is much better in case of low measurement errors. However, the impact of these errors can be reduced by: -Using continuous tracking approach: taking several measurements for each robot calibration configuration (i.e., the same applied force), and then averaging the collected data. Most sensors, and data collection card allow a frequency upper than 100 Hz. -Calibrating the wrench sensor only in a limited range, in which the sensor will be actually used. This will reduce the measurement uncertainty.

Conclusions
We have presented a self-kinematic calibration approach using a force-torque sensor. With this approach, the position accuracy of a 6-DOF medical robot was significantly improved. The robustness of our calibration model, regarding measurement noise, was confirmed through a simulation study, which also allowed us to conduct an observability analysis in order to identify the most appropriate calibration configurations. The simulation demonstrated that the robot's position errors were reduced from 12 mm to 0.320 mm for the maximum values, and from 9 mm to 0.277 mm for the mean errors. The calibration method presented in this paper will be tested experimentally in further work.