In modern industrial robotics, the demand for precise force sensing has become increasingly critical, particularly in applications such as welding, assembly, and grinding, where real-time monitoring of tool forces is essential. The integration of a six-axis force sensor at the robot wrist enables the measurement of three orthogonal forces and three orthogonal moments, providing a foundation for force feedback in teaching control, trajectory adjustment, and machining processes. However, challenges such as tool wear, environmental changes, and dynamic disturbances during high-speed operations often degrade the accuracy of static calibration methods. To address these issues, we propose an online calibration algorithm based on maximum likelihood estimation, which dynamically computes system parameters including tool center of gravity, installation angles, sensor zero offsets, and load gravity, even under high-speed motion conditions with significant vibrational noise.
Our approach focuses on leveraging real-time data collected during robot motion, such as during tool changes or idle high-speed movements, to continuously update calibration parameters. This method accounts for inertial forces and gravitational effects, ensuring robust performance in practical industrial settings. By employing probabilistic estimation techniques, we aim to mitigate the limitations of traditional least squares methods, which are susceptible to noise in dynamic environments. The core of our work involves formulating the mechanical relationships in the system, implementing optimization strategies, and validating the method through extensive experiments across varying speeds.
The six-axis force sensor plays a pivotal role in this framework, as it captures force and torque data alongside robot trajectory information. We emphasize the importance of accurate sensor calibration to maintain force perception integrity. In the following sections, we detail the mathematical model, experimental setup, and analysis, highlighting how maximum likelihood estimation enhances consistency and reliability compared to conventional approaches.
To begin, we define the key variables involved in the calibration process. The objective is to determine the center of gravity coordinates in the sensor frame, denoted as \( sP = [p_x, p_y, p_z]^T \), the tool mass \( m \) and weight \( G = mg \), the force zero offsets \( sF_0 = [f_{x0}, f_{y0}, f_{z0}]^T \), the moment zero offsets \( sM_0 = [m_{x0}, m_{y0}, m_{z0}]^T \), and the robot installation orientation represented by the rotation matrix \( g_cR \), which depends on angles \( \theta \) and \( \phi \). This matrix transforms coordinates from the robot base frame to the gravity frame and is given by:
$$ g_cR = \begin{bmatrix} \cos\phi & \sin\phi \sin\theta & \cos\phi \sin\theta \\ 0 & \cos\phi & -\sin\phi \\ -\sin\theta & \sin\phi \cos\theta & \cos\phi \cos\theta \end{bmatrix} $$
Known variables include the current robot pose \( c_f R_i \), obtained from the robot’s flange coordinates, and the force and moment readings \( F_i = [f_{xi}, f_{yi}, f_{zi}]^T \) and \( M_i = [m_{xi}, m_{yi}, m_{zi}]^T \) from the six-axis force sensor at each pose \( i \). The sensor installation orientation \( f_s R \) is derived from sensor and connector dimensions. The force and moment measurements are adjusted by subtracting their zero offsets: \( sF_i = F_i – sF_0 \) and \( sM_i = M_i – sM_0 \).
The fundamental force relationships are established by considering the gravity frame, where only the gravitational force acts along the vertical direction, resulting in zero moment about the origin. Thus, in the gravity frame, \( gF = [0, 0, -G]^T \) and \( gM = [0, 0, 0]^T \). Using robot statics principles, the forces and moments in the sensor frame are transformed as follows:
$$ \begin{bmatrix} sF \\ sM \end{bmatrix} = \begin{bmatrix} s_g R & 0 \\ S(s,o_g P) s_g R & s_g R \end{bmatrix} \begin{bmatrix} gF \\ gM \end{bmatrix} $$
Here, \( S(s,o_g P) \) is the skew-symmetric matrix derived from the center of gravity vector:
$$ S(s,o_g P) = \begin{bmatrix} 0 & -p_z & p_y \\ p_z & 0 & -p_x \\ -p_y & p_x & 0 \end{bmatrix} $$
The composite rotation matrix from the sensor frame to the gravity frame is \( g_s R = g_c R \cdot c_f R \cdot f_s R \). This leads to the moment equation in the sensor frame:
$$ sM = S(s,o_g P) \cdot s_g R \cdot gF = S(s,o_g P) \cdot sF $$
Expanding this, we obtain a linear system that facilitates the calculation of the center of gravity and intermediate variables. For instance, the moments can be expressed as:
$$ \begin{bmatrix} M_x \\ M_y \\ M_z \end{bmatrix} = \begin{bmatrix} 0 & f_z & -f_y \\ -f_z & 0 & f_x \\ f_y & -f_x & 0 \end{bmatrix} \begin{bmatrix} k_x \\ k_y \\ k_z \\ p_x \\ p_y \\ p_z \end{bmatrix} $$
where \( k_x, k_y, k_z \) are intermediate terms incorporating zero offsets and gravity components. Similarly, the force relationship is derived as:
$$ F_i = [I \quad c_s R^T_i] [f_{x0} \quad f_{y0} \quad f_{z0} \quad L^T]^T $$
with \( L = [G \cos\phi \sin\theta, -G \sin\phi, -G \cos\phi \cos\theta]^T \). This allows for the computation of force zero offsets, followed by moment zero offsets, gravity, and angles using:
$$ G = \sqrt{L_x^2 + L_y^2 + L_z^2} $$
$$ \phi = \arcsin(-L_y / G), \quad \theta = \arctan(-L_x / L_z) $$
To enhance accuracy, we incorporate inertial compensation by accounting for acceleration at the tool’s center of gravity. The acceleration components \( a_{px,i}, a_{py,i}, a_{pz,i} \) are estimated from the trajectory data, modifying the force equation to:
$$ F_i = [I \quad C_i] [f_{x0} \quad f_{y0} \quad f_{z0} \quad m]^T $$
where
$$ C_i = c_s R^T_i \begin{bmatrix} \cos\phi \sin\theta \\ -\sin\phi \\ -\cos\phi \cos\theta \end{bmatrix} g – \begin{bmatrix} a_{px,i} \\ a_{py,i} \\ a_{pz,i} \end{bmatrix} $$
For data processing, we compare ordinary least squares (LS), covariance-weighted least squares (CWLS), and multivariate maximum likelihood estimation (MLE). While LS is suitable for static conditions with minimal noise, MLE proves more robust in dynamic environments by modeling the probability distribution of errors, thus improving parameter consistency across varying operational speeds.
Our experimental platform consists of a six-axis force sensor, specifically the M4344D2 model, mounted on a Yaskawa MS165 robot. The sensor has a rated accuracy of 1% and a resolution of 0.5 units, with a force range of up to 3,200 N. Data acquisition is synchronized with robot motion at a 12 ms sampling interval. The end-effector includes a floating grinder and an pneumatic grinder, with total load weights of approximately 30.739 N and 19.806 N for different test groups. The robot executes a predefined trajectory involving rotations around its axes at speeds ranging from 10 mm/s to 1000 mm/s, with three repetitions per speed to ensure statistical reliability.
Initial calibration of the six-axis force sensor is performed using standard weights to verify static performance. The results, summarized in the table below, demonstrate high precision and resolution, which is crucial for reliable force control in robotic applications.
| Initial Reading (N) | Added Weight (N) | Final Reading (N) | Measured Force (N) | Error (N) |
|---|---|---|---|---|
| -37.0202 | 0.1961 | -37.1994 | -0.1792 | 0.0169 |
| -37.0328 | 0.3922 | -37.4228 | -0.3900 | 0.0022 |
| -37.0097 | 0.5883 | -37.5619 | -0.5522 | 0.0361 |
| -36.9448 | 0.7844 | -37.7134 | -0.7686 | 0.0158 |
| -36.8860 | 0.9805 | -37.8349 | -0.9489 | 0.0316 |
In the dynamic experiments, we first applied the least squares method to compute parameters such as the center of gravity, zero offsets, and installation angles. However, as shown in the following table, the results exhibited significant variability across different speeds due to vibrational noise and inertial effects, underscoring the limitations of LS in high-motion scenarios.
| Speed (mm/s) | px (mm) | py (mm) | pz (mm) | fx0 (N) | fy0 (N) | fz0 (N) | mx0 (N·m) | my0 (N·m) | mz0 (N·m) | G (N) | φ (°) | θ (°) |
|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 10 | -0.14 | -1.20 | 42.36 | -359.25 | -178.27 | 153.91 | -4.966 | 5.751 | 6.404 | 28.27 | 0.29 | 0.48 |
| 100 | -0.33 | -1.21 | 42.28 | -359.02 | -178.24 | 153.16 | -4.974 | 5.770 | 6.355 | 28.22 | -0.51 | 0.37 |
| 200 | -0.06 | -1.23 | 41.93 | -358.47 | -178.18 | 153.19 | -4.988 | 5.793 | 6.464 | 28.27 | -1.91 | 0.04 |
| 300 | -0.50 | -0.82 | 41.71 | -357.74 | -178.18 | 153.08 | -4.784 | 5.836 | 6.176 | 28.24 | -3.56 | -0.15 |
| 400 | -1.96 | 0.58 | 40.52 | -356.91 | -178.09 | 154.18 | -4.097 | 5.904 | 5.171 | 28.27 | -5.03 | -1.24 |
| 500 | -1.30 | 0.75 | 40.33 | -356.48 | -177.83 | 154.33 | -4.029 | 5.909 | 5.343 | 28.38 | -5.80 | -1.75 |
| 800 | -0.69 | -0.23 | 41.13 | -355.94 | -177.08 | 154.00 | -4.549 | 5.918 | 5.921 | 27.85 | -6.84 | -3.61 |
| 1000 | -0.44 | -0.43 | 41.21 | -356.03 | -177.10 | 154.06 | -4.644 | 5.923 | 6.780 | 27.95 | -7.16 | -2.72 |
To improve consistency, we applied maximum likelihood estimation to compute the center of gravity. The table below compares the results from LS, CWLS, and MLE, showing that MLE significantly reduces standard deviation and range, with the average standard deviation for the center of gravity decreasing from 0.67 mm to 0.23 mm.
| Speed (mm/s) | px (mm) LS | px (mm) CWLS | px (mm) MLE | py (mm) LS | py (mm) CWLS | py (mm) MLE | pz (mm) LS | pz (mm) CWLS | pz (mm) MLE |
|---|---|---|---|---|---|---|---|---|---|
| 10 | -0.138 | 0.061 | 0.220 | -1.200 | -2.281 | -2.146 | 42.357 | 41.744 | 41.659 |
| 100 | -0.333 | -0.145 | 0.124 | -1.214 | -2.382 | -2.054 | 42.279 | 41.617 | 41.443 |
| 200 | -0.059 | 0.187 | 0.400 | -1.233 | -2.673 | -2.236 | 41.929 | 41.199 | 41.093 |
| 300 | -0.504 | -0.141 | 0.115 | -0.823 | -2.694 | -2.050 | 41.709 | 41.048 | 40.999 |
| 400 | -1.963 | -1.155 | -0.330 | 0.582 | -2.390 | -2.121 | 40.822 | 40.549 | 40.721 |
| 500 | -1.295 | -0.482 | 0.235 | 0.750 | -2.057 | -1.795 | 40.332 | 40.256 | 40.660 |
| 800 | -0.685 | 0.045 | -0.068 | -0.231 | -2.779 | -2.093 | 41.130 | 40.727 | 41.238 |
| 1000 | -0.438 | 0.114 | 0.379 | -0.436 | -2.404 | -1.982 | 41.207 | 40.799 | 40.918 |
For the installation angles, we performed a grid search over a range of \([-5^\circ, 5^\circ]\) with a step size of 0.001°, using the log-likelihood function from MLE. The optimal angles were found to be \(\theta = -0.284^\circ\) and \(\phi = 0.022^\circ\), which were then fixed in subsequent calculations to enhance stability.
With the angles determined, we computed the force zero offsets, moment zero offsets, and gravity using MLE with inertial compensation. The consistency of these parameters across speeds is summarized in the following table, demonstrating that MLE reduces the standard deviation of force zeros from 0.73 N to 0.27 N and moment zeros from 0.29 N·m to 0.05 N·m on average.
| Parameter | LS Range | LS Std Dev | MLE Range | MLE Std Dev |
|---|---|---|---|---|
| fx0 (N) | 3.3140 | 1.2435 | 0.5320 | 0.1886 |
| fy0 (N) | 1.1978 | 0.4681 | 0.5091 | 0.1859 |
| fz0 (N) | 1.2475 | 0.4762 | 1.3416 | 0.4274 |
| G (N) | 0.5258 | 0.1692 | 0.9299 | 0.2566 |
| mx0 (N·m) | 0.9596 | 0.3598 | 0.0161 | 0.0056 |
| my0 (N·m) | 0.1727 | 0.0672 | 0.0347 | 0.0144 |
| mz0 (N·m) | 1.2935 | 0.4564 | 0.2945 | 0.0841 |
Additionally, we validated the weight calculation by comparing two load groups. The difference in computed weights should match the known difference of 10.933 N. As shown below, the MLE results exhibit smaller errors and higher consistency, confirming the method’s accuracy in dynamic conditions.
| Speed (mm/s) | Group A G (N) | Group B G (N) | Weight Difference (N) | Error (N) |
|---|---|---|---|---|
| 10 | 27.9693 | 17.2368 | 10.7325 | -0.2001 |
| 100 | 28.3080 | 17.5991 | 10.7089 | -0.2237 |
| 200 | 28.4478 | 17.7309 | 10.7169 | -0.2157 |
| 300 | 28.5490 | 17.4502 | 11.0988 | 0.1662 |
| 400 | 28.4602 | 17.4211 | 11.0391 | 0.1065 |
| 500 | 28.8992 | 17.4015 | 11.4977 | 0.5651 |
| 800 | 28.5880 | 17.1616 | 11.4264 | 0.4938 |
| 1000 | 28.6881 | 17.3370 | 11.3511 | 0.4185 |
To illustrate the impact of inertial compensation, we analyzed force data along the z-axis at 1 m/s. Without compensation, the force zero estimate was 154.29 N, with an average error of 2.49 N and a maximum error of 5.70 N compared to a smoothed reference. With MLE and inertial compensation, the zero estimate improved to 153.45 N, reducing the average error to 1.10 N and maximum error to 3.42 N. This demonstrates that accounting for acceleration effects significantly enhances force perception accuracy.
In conclusion, our proposed maximum likelihood estimation method for online calibration of six-axis force sensors effectively addresses the challenges posed by high-speed robot motion. By dynamically updating parameters and incorporating inertial compensation, we achieve superior consistency and accuracy compared to traditional least squares approaches. The six-axis force sensor is integral to this process, enabling real-time force monitoring in industrial applications. Future work will explore adaptive filtering techniques and integration with other sensor modalities to further improve robustness in complex environments.