An Adaptive External Torque Estimation Algorithm for Collision Detection in Robotic Arms

1. Introduction

In physical human–robot interaction (pHRI) [1,2], quantifying the severity of collisions on the human body provides critical guidance for higher-level systems in managing post-collision events. Haddadin et al. [3,4] systematically investigated this problem through extensive experiments, quantifying the severity of impacts on different parts of the human body. Their results demonstrated that relative velocity, effective mass, and the elastic modulus of the contact surface all influence collision severity, with the relative velocity at impact identified as the dominant factor. In contrast, the contribution of the robot’s reflected mass exhibits a saturation effect [5,6]. Building on this work, Haddadin later introduced the concept of a “collision event pipeline” [7], a generic stage-based framework that decomposes the collision process into multiple phases to facilitate systematic analysis.

A generalized momentum–based collision detection method was first introduced by Haddadin [3] and experimentally validated in various pHRI scenarios. This approach is simple to implement, requires no external sensors, and has demonstrated broad applicability. At its core, the first-order momentum observer (FOMO) estimates external torques by computing momentum residuals, without relying on joint angular accelerations or the inverse of the inertia matrix, thereby substantially reducing sensitivity to measurement noise [8]. Building on this framework, the second-order momentum observer (SOMO) [9] exhibits stronger robustness to velocity noise, albeit at the expense of slower response and without a complete theoretical proof. An improved second-order observer [10,11] enhances dynamic performance by incorporating a feedforward term, which shortens response time but compromises disturbance rejection relative to the standard SOMO. More recently, Liu [12] proposed a higher-order generalized momentum observer that achieves stronger disturbance rejection by adopting a higher cutoff frequency, though at the cost of increased computational burden inherent to higher-order systems.

Kalman filtering has shown strong potential for suppressing noise in external torque estimation. Wahrburg et al. [13] first applied Kalman filtering to robotic manipulators, formulating a state-space model based on generalized momentum that effectively attenuated noise but exhibited limited performance in handling nonlinearities such as friction. Moreover, due to the inherent properties of Kalman filtering, inaccuracies in covariance matrices led to amplified estimation errors. To address these limitations, Yigit et al. [14] employed the extended Kalman filter (EKF), which improved robustness against nonlinear factors including friction, though estimation accuracy remained suboptimal. Wang et al. [15] further advanced this approach by integrating encoder and accelerometer measurements within an EKF framework to estimate link-side states, achieving substantial improvements in estimation precision and dynamic responsiveness. However, the reliance on accelerometers increased both system cost and complexity. More recently, Dong et al. [14] proposed a sliding-window variational adaptive Kalman filtering method grounded in an error-probability model. By adaptively updating the external force parameters of the ECM within the process model, this approach significantly reduced estimation error without requiring additional sensors.

In recent years, extensive efforts have been devoted to adaptive Kalman filtering algorithms and external torque estimation [16,17,18] in robotic systems. A defining feature of adaptive Kalman filters is their ability to provide real-time estimates of noise covariance. For example, Ref. [19] employed perturbation analysis of the discrete-time algebraic Riccati equation (DARE) to examine the sensitivity of Kalman filters to deviations in noise covariance estimates, thereby assessing estimation accuracy. Building on this line of work, Zhenwei Li et al. [20] developed a variational Bayesian maximum correntropy cubature Kalman filter with adaptive kernel bandwidth, where the bandwidth was updated online using Mahalanobis distance. This approach substantially reduced the influence of non-Gaussian measurement noise on the estimation results. However, the method suffers from high computational cost and sensitivity to parameter tuning, which constrains its applicability in broader scenarios.

In this work, we propose a Kalman filtering–based algorithm for estimating interaction torques in robotic manipulators under conditions of uncertain or time-varying process noise covariance. The method employs variational inference to estimate unknown process noise covariances, decoupling noise parameters from manipulator configurations and payload characteristics. At each measurement cycle, the posterior probability density function of the noise covariance is iteratively updated within a closed-form variational Bayesian framework, yielding an approximate Gaussian posterior of the process noise. This enables optimal external torque estimation through Kalman filtering. The approach is experimentally validated on a Franka manipulator and benchmarked against several state-of-the-art methods for external torque estimation in real robotic systems.

2. Materials and Methods

2.1. System Dynamics of Robotic Manipulators

When a robotic manipulator interacts with the external environment through force exchange, and assuming rigid joints and rigid links, the system dynamics can be described by the following equation of motion [21]:

$$\mathit{M}\left(q\right)\ddot{q}+\mathit{C}\left(q,\dot{q}\right)\dot{q}+\mathit{g}\left(q\right)=\tau ={\tau }_m+{\tau }_e,$$

(1)

where the position and velocity of the linkage are denoted by $q,\dot{q}\in {\mathbb{R}}^n$ where n is the number of joints, respectively, and the two variables are the states of the robot. $\mathit{M}\left(q\right)\in {\mathbb{R}}^{n\times n}$denotes the symmetric positive definite inertia matrix, $\dot{\mathit{C}(q,q)}\in {\mathbb{R}}^{n\times n}$denotes the Coriolis matrix satisfying passivity and $\mathit{g}\left(q\right)\in {\mathbb{R}}^n$ is the gravity torque vector. The torque generated by the motor $\tau \in {\mathbb{R}}^n$ and the external torque ${\tau }_e$ are the inputs to the system.

From the reformulation of the manipulator dynamics under the momentum concept, we obtain:

$$\dot{\mathit{p}}={\tau }_m+{\mathit{C}}^T\left(\mathit{q},\dot{\mathit{q}}\right)\dot{\mathit{q}}-\mathit{g}\left(\mathit{q}\right)+{\tau }_e,$$

(2)

The standard KF provides the optimal linear estimate of system states under the assumption of an accurate model, even in the presence of noisy measurements. To estimate the interactive torque of a robotic arm using the KF, a model of the robotic arm must be formulated. In this context, ${\tau }_a={\tau }_m+{\mathit{C}}^T(\dot{q},q)-\mathit{g}\left(q\right)$ and a noise term $w_p\sim N(0,{\mathit{Q}}_p)$. The time derivative of momentum can be interpreted as the resultant of all forces acting on the system at that moment. Consequently, the dynamic equation of the robotic arm can be rewritten as:

$$\dot{\mathit{p}}=\overline{\tau }+{\mathit{J}}^T\mathit{f}+{\omega }_p,$$

(3)

Following the design principle of the classical disturbance observer, the driving torque is modeled as a constant, whereas the derivative of the external torque is treated as noise-driven. Set $\dot{\mathit{f}}={\omega }_f$, ${\omega }_f\sim \mathrm{N}(0,{\mathit{Q}}_{c,f})$, and the system state variables are defined as $\mathit{x}=\left[\begin{array}{c}\mathit{p}\\ \mathit{f}\end{array}\right]$, can thus be formulated in the following state-space form:

$$\left[\begin{array}{c}\dot{\mathit{p}}\\ \dot{\mathit{f}}\end{array}\right]=\left[\begin{array}{cc}{0}_{n\times n}& -{\mathit{J}}^T\\ 0_{n_{ext}\times n}& 0_{n_{ext}\times n_{ext}}\end{array}\right]\left[\begin{array}{c}\mathit{p}\\ \mathit{f}\end{array}\right]+\left[\begin{array}{c}{\mathit{I}}_n\\ 0_{n_{ext}\times n}\end{array}\right]\overline{\tau }+\omega$$

(4)

$${\mathit{p}}_{meas}=\left[\begin{array}{cc}{\mathit{I}}_n& 0_{n\times n_{ext}}\end{array}\right]\left[\begin{array}{c}\mathit{p}\\ \mathit{f}\end{array}\right]+\upsilon$$

(5)

Here, ${\mathit{p}}_{meas}$ denotes the measured generalized momentum, $\upsilon \sim \mathrm{N}(0,{\mathit{R}}_c)$ represents the system measurement noise, and by augmenting the system state variable $\mathit{x}$, it follows that the process noise $\omega$ is assumed to follow the normal distribution:

$$\omega \sim \mathrm{N}\left(0,\left[\begin{array}{cc}{\mathit{Q}}_{c,p}& 0_{n\times n_{ext}}\\ 0_{n_{ext}\times n}& {\mathit{Q}}_{c,f}\end{array}\right]\right)$$

(6)

Here, ${\mathit{Q}}_{c,p}$ and ${\mathit{Q}}_{c,f}$ denote components of the system process noise covariance matrix, while ${\mathit{R}}_c$ represents the system measurement noise covariance matrix.

$$\mathit{x}=\left[\begin{array}{cc}\mathit{p}& \mathit{J}{\mathit{\tau }}_e\end{array}\right]$$

(7)

$${\mathit{A}}_C=\left[\begin{array}{cc}{0}_{n\times n}& -{\mathit{J}}^T\\ 0_{n_{ext}\times n}& 0_{n_{ext}\times n_{ext}}\end{array}\right]$$

(8)

$${\mathit{B}}_c=\left[\begin{array}{c}{\mathit{I}}_n\\ 0_{n_{ext}\times n}\end{array}\right]$$

(9)

$$\mathit{u}=\overline{\tau }$$

(10)

$$\mathit{y}=p_{meas}$$

(11)

$${\mathit{C}}_c=\left[\begin{array}{cc}{\mathit{I}}_n& 0_{n\times n_{ext}}\end{array}\right]$$

(12)

$${\mathit{Q}}_c=\left[\begin{array}{cc}{\mathit{Q}}_{c,p}& 0_{n\times n_{ext}}\\ 0_{n_{ext}\times n}& {\mathit{Q}}_{c,f}\end{array}\right]$$

(13)

Accordingly, the state-space representation can be rewritten as:

$$\dot{\mathit{x}}={\mathit{A}}_c\mathit{x}+{\mathit{B}}_c\mathit{u}+\omega$$

(14)

$$\mathit{y}={\mathit{C}}_c\mathit{x}+\upsilon$$

(15)

where, $\omega \sim \mathrm{N}(0,{\mathit{Q}}_c)$, $\upsilon \sim \mathrm{N}(0,{\mathit{R}}_c)$, upon discretization, the equation can be written as:

$${\mathit{x}}_{k+1}=\mathit{A}{\mathit{x}}_k+\mathit{B}{\mathit{u}}_k+{\omega }_k$$

(16)

$${\mathit{y}}_k=\mathit{C}{\mathit{x}}_k+{\upsilon }_k$$

(17)

Discretization of the above equation yields:

$$\left[\begin{array}{cc}\mathit{A}& \mathit{B}\\ 0& \mathit{I}\end{array}\right]=exp\left(\left[\begin{array}{cc}{\mathit{A}}_c& {\mathit{B}}_c\\ 0& 0\end{array}\right]T_s\right)$$

(18)

$$\left[\begin{array}{cc}{\mathit{M}}_{11,k}& {\mathit{M}}_{12,k}\\ 0& {\mathit{M}}_{22,k}\end{array}\right]=exp(\mathit{H}·T_s)$$

(19)

$${\mathit{Q}}_k={\mathit{M}}_{12,k}{\left({\mathit{M}}_{11,k}\right)}^T$$

(20)

Meanwhile, the parameters of the system state equation can be expressed as: $\mathit{C}={\mathit{C}}_c$, ${\mathit{R}}_k={\displaystyle \frac{1}{T_s}}{\mathit{R}}_c$, $\mathit{H}=\left[\begin{array}{cc}\mathit{A}& {\mathit{Q}}_c\\ 0& -{\mathit{A}}^T\end{array}\right]$, and $T_s$ is sampling time.

2.2. Effects of Noise Modeling Inaccuracies on External Torque Estimation

The optimality of the Kalman filter fundamentally depends on accurate knowledge of process and measurement noise covariances. However, in robotic manipulators, process noise varies with configuration and joint velocity, making fixed covariance assumptions unrealistic. Such mismatches reduce estimation accuracy, leading to false collision detections or erroneous external torque estimates that may compromise system safety. Prior work [14,22] analyzed these effects by comparing filter-calculated (FMSE), true (TMSE), and ideal (IMSE) mean squared errors, providing a quantitative view of estimation error under covariance mismatch.

To analyze the performance of the filter under conditions of process noise covariance mismatch, the process noise covariance matrix ${\mathit{Q}}_{k,k-1}$ is defined as follows:

$${\mathit{Q}}_{k,k-1}^u={\mathit{Q}}_{k,k-1}+\Delta {\mathit{Q}}_{k,k-1}$$

(21)

Here, ${\mathit{Q}}_{k,k-1}^u$ denotes the process noise covariance matrix employed in the filtering equations, ${\mathit{Q}}_{k,k-1}$ represents the true process noise covariance matrix of the system, and $\Delta {\mathit{Q}}_{k,k-1}$ denotes the deviation between them. When the practically applied process noise covariance matrix ${\mathit{Q}}_{k,k-1}^u$ is substituted into the Kalman filter framework, the recursive form of the filter can be expressed as:

$${\widehat{\mathit{x}}}_{k|k-1}^f={\mathit{A}}_{k,k-1}{\widehat{\mathit{x}}}_{k-1|k-1}^f$$

(22)

$${\mathit{P}}_{k|k-1}^f={\mathit{A}}_{k,k-1}{\mathit{P}}_{k-1|k-1}^f{\mathit{A}}_{k,k-1}^{\prime }+{\mathit{Q}}_{k,k-1}^u$$

(23)

$${\mathit{K}}_k^f={\mathit{P}}_{k|k-1}^f{\mathit{C}}_k^{\prime }{({\mathit{C}}_k{\mathit{P}}_{k|k-1}^f{\mathit{C}}_k^{\prime }+{\mathit{R}}_k)}^{-1}$$

(24)

$${\widehat{\mathit{x}}}_{k|k}^f={\widehat{\mathit{x}}}_{k|k-1}^{\phantom{,}}+{\mathit{K}}_k^f({\mathit{z}}_k-{\mathit{C}}_k{\widehat{\mathit{x}}}_{k|k-1}^{\phantom{,}})$$

(25)

$${\mathit{P}}_{k|k}^f=(\mathit{I}-{\mathit{K}}_k^f{\mathit{C}}_k){\mathit{P}}_{k|k-1}^f{(\mathit{I}-{\mathit{K}}_k^f{\mathit{C}}_k)}^{\prime }+{\mathit{K}}_k^f{\mathit{R}}_k{\left({\mathit{K}}_k^f\right)}^{\prime }$$

(26)

Here, the superscript f denotes the output term computed by the filter. According to the recursive equations, at time step k, the covariance matrix ${\widehat{\mathit{x}}}_{k|k}^f$ must be calculated under the current state estimate ${\mathit{P}}_{k|k}^f$. However, due to the substitution of ${\mathit{Q}}_{k,k-1}^u$, the state estimates across two consecutive time steps are given by:

$${\tilde{\mathit{x}}}_{k|k-1}^f=\mathit{x}-{\widehat{\mathit{x}}}_{k|k-1}^f={\mathit{A}}_{k,k-1}{\tilde{\mathit{x}}}_{k-1|k-1}^f+{\mathit{w}}_{k,k-1}$$

(27)

$${\tilde{\mathit{x}}}_{k|k}^f=\mathit{x}-{\widehat{\mathit{x}}}_{k|k}^f=(\mathit{I}-{\mathit{K}}_k^f{\mathit{C}}_k){\mathit{A}}_{k,k-1}{\tilde{\mathit{x}}}_{k-1|k-1}^f+(\mathit{I}-{\mathit{K}}_k^f{\mathit{C}}_k){\mathit{w}}_{k,k-1}-{\mathit{K}}_k^f{\mathit{v}}_k$$

(28)

In the presence of covariance mismatch $\Delta {\mathit{Q}}_{k,k-1}$, the covariance matrix ${\mathit{P}}_{k|k}^{\phantom{,}}$ under the ideal mean squared error (IMSE) fails to capture the true estimation error of the Kalman filter, thereby distorting the update of the Kalman gain $\mathit{K}$.

2.3. Variational Inference–Based Adaptive External Torque Estimation Algorithm

The complexity of robotic systems makes process noise difficult to model, which limits the accuracy of Kalman filter–based external torque estimation. To address time-varying and non-Gaussian noise, we employ variational inference to approximate latent variable distributions from historical system data, enabling online estimation of the noise covariance. This adaptive framework allows dynamic parameter adjustment in response to noise variations, thereby improving the accuracy of external torque estimation for robotic systems.

Building on the feasibility of Kalman filtering for external torque estimation in robotic manipulators, the modeling error of ${\mathit{Q}}_k$ ultimately leads to inaccuracies in torque estimation. The problem addressed in this section therefore reduces to deriving an approximate distribution of ${\mathit{Q}}_k$ in order to improve the accuracy of external torque estimation.

Using variational inference, this subsection provides a detailed derivation of the implementation of variational Kalman filtering within robotic dynamics, with the estimation procedure encapsulated into an algorithmic module. It is assumed that the system noise ${\mathit{Q}}_k$ is governed by an independent dynamic model, decoupled from the system state variable ${\mathit{x}}_k$:

$$p({\mathit{x}}_k,{\mathit{Q}}_k|{\mathit{x}}_{k-1},{\mathit{Q}}_{k-1})=p\left({\mathit{x}}_k\right|{\mathit{x}}_{k-1})p\left({\mathit{Q}}_k\right|{\mathit{Q}}_{k-1}).$$

(29)

According to the Chapman–Kolmogorov relation, the prior of the system state ${\mathit{x}}_k$ jointly with the process noise covariance ${\mathit{Q}}_k$ is given by:

$$p\left({\mathit{x}}_k,{\mathit{Q}}_k|{\mathit{y}}_{1:k-1}\right)=\int p\left({\mathit{x}}_k\right|{\mathit{x}}_{k-1},{\mathit{Q}}_k)p\left({\mathit{Q}}_k\right|{\mathit{Q}}_{k-1})\times p({\mathit{x}}_{k-1},{\mathit{Q}}_{k-1}|{\mathit{y}}_{1:k-1})\mathrm{d}{\mathit{x}}_{k-1}\mathrm{d}{\mathit{Q}}_{k-1},$$

(30)

with the incorporation of the subsequent measurement, Bayes’ theorem yields the posterior distribution:

$$p\left({\mathit{x}}_k,{\mathit{Q}}_k|{\mathit{y}}_{1:k}\right)\propto p\left({\mathit{y}}_k\right|{\mathit{x}}_k,{\mathit{Q}}_k)p({\mathit{x}}_k,{\mathit{Q}}_k|{\mathit{y}}_{1:k-1})$$

(31)

Variational inference approximates an intractable posterior with a tractable distribution, reformulating inference as the optimization of a closed-form variational distribution. Recalling the prediction step of the classical Kalman filter:

$$\begin{array}{c}{\mathit{x}}_{k\mid k-1}={\mathit{A}}_k{\mathit{x}}_{k-1\mid k-1}\hfill \\ {\mathit{P}}_{k\mid k-1}={\mathit{A}}_k{\mathit{P}}_{k-1\mid k-1}{\mathit{A}}_k^T+{\mathit{Q}}_k\hfill \end{array}$$

(32)

In this formulation, ${\mathit{A}}_k$ denotes the state transition matrix. The predicted covariance matrix depends not only on the process noise covariance matrix ${\mathit{Q}}_k$, but also on ${\mathit{A}}_k{\mathit{P}}_{k-1\mid k-1}{\mathit{A}}_k^T$. To fully characterize the predicted covariance, we introduce a latent variable $\mathit{\theta }={\mathit{A}}_k{\mathit{x}}_{k-1\mid k-1}$, leading to the representation $\Sigma ={\mathit{A}}_k{\mathit{P}}_{k-1\mid k-1}{\mathit{A}}_k^T$.

The process noise covariance matrix ${\mathit{Q}}_k$ is assumed to have a prior governed by an Inverse-Wishart (IW) distribution.According to mean-field theory, one obtains:

$$p\left({\mathit{x}}_k,{\mathit{\theta }}_k,{\mathit{Q}}_k|{\mathit{y}}_{1:k}\right)\approx q({\mathit{x}}_k,{\mathit{\theta }}_k,{\mathit{Q}}_k)\approx q\left({\mathit{x}}_k\right)q\left({\mathit{\theta }}_k\right)q\left({\mathit{Q}}_k\right)$$

(33)

For brevity, $\mathit{z}$ is used in place of $({\mathit{x}}_k,{\mathit{\theta }}_k,{\mathit{Q}}_k)$, yielding:

$$p\left(\mathit{z}|{\mathit{y}}_{1:k}\right)\approx q\left(\mathit{z}\right)\approx q\left({\mathit{z}}_1\right)q\left({\mathit{z}}_2\right)q\left({\mathit{z}}_3\right)$$

(34)

Formally, the divergence between $q\left(\mathit{z}\right)$ and $p\left(\mathit{z}|{\mathit{y}}_{1:k}\right)$ is defined as:

$$\mathrm{KL}\left(q_z\left(\mathit{z}\right)\left|\right|p\left(\mathit{z}|{\mathit{y}}_{1:k}\right)\right)=-\int q_z\left(\mathit{z}\right)log{\displaystyle \frac{p\left(\mathit{z}|{\mathit{y}}_{1:k}\right)}{q_z\left(\mathit{z}\right)}}d\mathit{z}$$

(35)

By the evidence lower bound (ELBO) property, we obtain:

$$L\left(q_z\left(\mathit{z}\right)\right)=\int q_z\left(\mathit{z}\right)log{\displaystyle \frac{p\left(\mathit{z},{\mathit{y}}_k|{\mathit{y}}_{1:k-1}\right)}{q_z\left(\mathit{z}\right)}}d\mathit{z}$$

(36)

$$logq\left({\mathit{x}}_k\right)={\mathrm{E}}_{q({\mathit{\theta }}_k,{\mathit{Q}}_k)}[logp\left({\mathit{x}}_k,{\mathit{\theta }}_k,{\mathit{Q}}_k,{\mathit{y}}_k|{\mathit{y}}_{1:k-1}\right)]$$

(37)

$$logq\left({\mathit{\theta }}_k\right)={\mathrm{E}}_{q({\mathit{x}}_k,{\mathit{Q}}_k)}[logp\left({\mathit{x}}_k,{\mathit{\theta }}_k,{\mathit{Q}}_k,{\mathit{y}}_k|{\mathit{y}}_{1:k-1}\right)]$$

(38)

$$logq\left({\mathit{Q}}_k\right)={\mathrm{E}}_{q({\mathit{x}}_k,{\mathit{\theta }}_k)}[logp\left({\mathit{x}}_k,{\mathit{\theta }}_k,{\mathit{Q}}_k,{\mathit{y}}_k|{\mathit{y}}_{1:k-1}\right)]$$

(39)

Expanding $logp\left({\mathit{x}}_k,{\mathit{\theta }}_k,{\mathit{Q}}_k,{\mathit{y}}_k|{\mathit{y}}_{1:k-1}\right)$ yields:

$$\begin{array}{c}logp\left({\mathit{x}}_k,{\mathit{\theta }}_k,{\mathit{Q}}_k,{\mathit{y}}_k|{\mathit{y}}_{1:k-1}\right)=logp\left({\mathit{x}}_k|{\mathit{y}}_{1:k-1},{\mathit{\theta }}_k,{\mathit{Q}}_k\right)+logp\left({\mathit{\theta }}_k|{\mathit{y}}_{1:k-1}\right)+\hfill \\ logp\left({\mathit{Q}}_k|{\mathit{y}}_{1:k-1}\right)+logp\left({\mathit{y}}_k|{\mathit{x}}_k\right)\hfill \\ =logN\left({\mathit{x}}_k|{\mathit{\theta }}_k,{\mathit{Q}}_k\right)+logN\left({\mathit{\theta }}_k|{\mathit{x}}_{k|k-1},{\Sigma }_k\right)+\hfill \\ logIW\left({\mathit{Q}}_k|{\mathit{u}}_{k|k-1},{\mathit{U}}_{k|k-1}\right)+logN\left({\mathit{y}}_k\right|{\mathit{H}}_k{\mathit{x}}_k,{\mathit{R}}_k)\hfill \end{array}$$

(40)

It is straightforward to derive the approximate posterior distributions of $q\left({\mathit{x}}_k\right)$, $q\left({\mathit{\theta }}_k\right)$, and $q\left({\mathit{Q}}_k\right)$.

$$q\left({\mathit{x}}_k\right)=\mathrm{N}\left\{{\mathit{x}}_k|\mathrm{E}\left({\mathit{\theta }}_k\right),{\left[\mathrm{E}\left({{\mathit{Q}}_k}^{-1}\right)\right]}^{-1}\right\}\mathrm{N}\left({\mathit{y}}_k|{\mathit{H}}_k{x}_k,{\mathit{R}}_k\right)=\mathrm{N}\left({\mathit{x}}_k|{\mathit{x}}_{k|k},{\mathit{P}}_{k|k}\right)$$

(41)

$$q\left({\mathit{\theta }}_k\right)=\mathrm{N}\left\{\mathrm{E}\left({\mathit{x}}_k\right)|{\mathit{\theta }}_k,{\mathit{A}}_k\right\}\mathrm{N}\left({\mathit{\theta }}_k|{\mathit{x}}_{k|k-1},{\Sigma }_k\right)=\mathrm{N}\left({\mathit{\theta }}_k|{\mathit{\theta }}_{k|k},{\mathit{P}}_{\theta ,k|k}\right)$$

(42)

$$logq\left({\mathit{Q}}_k\right)=-{\displaystyle \frac{1}{2}}\left({\mathit{u}}_{k|k-1}+{\mathit{n}}_x+2\right)log\left|{\mathit{Q}}_k\right|-{\displaystyle \frac{1}{2}}tr\left[\left({\mathit{U}}_{k|k-1}+{\mathit{B}}_k\right){{\mathit{Q}}_k}^{-1}\right]+{{\mathit{C}}_3}^{\prime }$$

(43)

where ${{\mathit{C}}_3}^{\prime }$ is an indeterminate constant, and set ${\mathit{B}}_k=\mathit{E}\left[\left({\mathit{x}}_k-{\mathit{\theta }}_k\right){\left({\mathit{x}}_k-{\mathit{\theta }}_k\right)}^T\right]$, we obtain:

$$\begin{array}{c}{\mathit{B}}_k=\left[\mathrm{E}\left({\mathit{x}}_k\right)-\mathrm{E}\left({\mathit{\theta }}_k\right)\right]{\left[\mathrm{E}\left({\mathit{x}}_k\right)-\mathrm{E}\left({\mathit{\theta }}_k\right)\right]}^T+\mathrm{cov}\left({\mathit{x}}_k{{\mathit{x}}_k}^T\right)+\mathrm{cov}\left({\mathit{\theta }}_k{{\mathit{\theta }}_k}^T\right)\hfill \\ =\left[\mathrm{E}\left({\mathit{x}}_k\right)-\mathrm{E}\left({\mathit{\theta }}_k\right)\right]{\left[\mathrm{E}\left({\mathit{x}}_k\right)-\mathrm{E}\left({\mathit{\theta }}_k\right)\right]}^T+{\mathit{P}}_{k|k}+{\mathit{P}}_{\theta ,k|k}\hfill \end{array}$$

(44)

Under the stated assumption, $q\left({\mathit{Q}}_k\right)\propto IW\left({\mathit{Q}}_k\right|{\mathit{u}}_{k|k},{\mathit{U}}_{k|k})$, ${\mathit{u}}_{k|k}={\mathit{u}}_{k|k-1}+1$, ${\mathit{U}}_{k|k}={\mathit{U}}_{k|k-1}+{\mathit{B}}_k$:

$$\mathrm{E}\left({{\mathit{Q}}_k}^{-1}\right)={\mathit{u}}_{k|k}{{\mathit{U}}_{k|k}}^{-1}$$

(45)

accordingly, we have:

$${\mathit{A}}_k={\displaystyle \frac{{\mathit{U}}_{k|k}}{{\mathit{u}}_{k|k}}}$$

(46)

The variational Kalman-based external torque estimation method employs a closed-loop iterative scheme to estimate the process noise covariance, which is then incorporated into the momentum-based Kalman filter torque observer. This formulation preserves the linear optimality of the Kalman filter even under covariance mismatch.

3. Results

To evaluate the performance of the proposed torque estimation algorithm, experiments were conducted on the Franka robotic arm, chosen as the platform for collision detection and external torque estimation. Developed by Franka Emika in Germany, Franka integrates high precision, flexibility, and safety, and is widely employed in research and industry. Each joint is equipped with a high-resolution torque sensor (<0.05 N), making it particularly suitable for experimental validation.

Figure 1. Robotic arm experimental platform.

Figure 1. Robotic arm experimental platform.

The experimental protocol comprised three scenarios: torque estimation under fixed-configuration, no-load conditions to assess noise suppression; rigid-box impacts at the end-effector to evaluate collision detection; and transient human-hand interactions to examine external torque estimation under compliant contact.

Figure 2. (a)Transient hand-induced impact and push on the robotic arm. (b) Corresponding sensor data recorded during the interaction.

Figure 2. (a)Transient hand-induced impact and push on the robotic arm. (b) Corresponding sensor data recorded during the interaction.

External torque of the Franka arm under collision-free conditions was estimated using ESO, KFMO, SH-AKFMO, and VBKFMO-Q.

Figure 3 indicates that, even in the absence of external forces, all considered momentum observers exhibit non-negligible estimation bias arising from modeling inaccuracies. During external torque estimation, abrupt transients may be erroneously interpreted as collisions. However, the observations presented in Figure 4 reveal that the absolute estimation errors remain low across the different observers, minimizing the likelihood of false collision detection. This pattern is further substantiated by the results summarized in Table 1.

Figure 5 and Figure 6 depict the observer outputs and corresponding joint-wise external torque estimation errors during collisions with a rigid object. Notably, for joints 1 through 4, ESO and KF exhibit substantial overshoot and discernible phase lag, whereas joints 5 through 7 show comparable performance across all four observers. This discrepancy arises because the collision occurs near the end-effector, with proximal joints (5–7) experiencing minimal perturbation, resulting in estimation behavior akin to collision-free conditions.

Evaluation based on root-mean-square error indicates that VBKFMO-Q achieves marginally superior external torque estimation relative to SH-AKFMO, while ESO and KF demonstrate considerably lower accuracy.

Figure 7 depicts joint torque estimation during a brief ( 0.5 s) interaction with a compliant object (human hand) at the robotic arm’s end-effector, showing the external torque observations from the different observers. The corresponding joint-wise estimation error curves are presented in Figure 8, and Table 2 summarizes the root-mean-square errors for each observer throughout this interaction. Notably, the red dashed line represents the proposed VBKFMO-Q observer, which exhibits superior rapid response and estimation accuracy under soft contact conditions, closely approximating the actual external torque. Its smooth output profile further underscores its suitability for practical engineering applications.

Table 2. Joint torque estimation RMS errors during high-stiffness impacts.

Algorithm	ESO	KFMO	SH-AKFMO	VBKFMO-Q
RMSE(J1)	1.3122	1.3756	0.4347	0.331
RMSE(J2)	2.9709	2.7695	1.4114	1.3468
RMSE(J3)	2.4096	2.1156	1.0478	1.1456
RMSE(J4)	0.8097	0.8055	0.2118	0.1577
RMSE(J5)	0.429	0.4187	0.4582	0.4479
RMSE(J6)	0.1144	0.099	0.1173	0.1154
RMSE(J7)	0.1181	0.0917	0.1124	0.1143

Table 2. Joint torque estimation RMS errors during high-stiffness impacts.

Algorithm	ESO	KFMO	SH-AKFMO	VBKFMO-Q
RMSE(J1)	1.3122	1.3756	0.4347	0.331
RMSE(J2)	2.9709	2.7695	1.4114	1.3468
RMSE(J3)	2.4096	2.1156	1.0478	1.1456
RMSE(J4)	0.8097	0.8055	0.2118	0.1577
RMSE(J5)	0.429	0.4187	0.4582	0.4479
RMSE(J6)	0.1144	0.099	0.1173	0.1154
RMSE(J7)	0.1181	0.0917	0.1124	0.1143

Table 3. Joint torque estimation RMS errors during compliant contact (human hand).

Algorithm	ESO	KFMO	SH-AKFMO	VBKFMO-Q
RMSE(J1)	1.0266	0.8956	1.0259	0.286
RMSE(J2)	0.8367	0.6369	0.8364	0.4485
RMSE(J3)	1.3626	1.2528	1.3612	1.0145
RMSE(J4)	1.0637	0.9181	1.0622	0.4996
RMSE(J5)	0.1237	0.0791	0.1244	0.042
RMSE(J6)	0.1017	0.0876	0.1439	0.0968
RMSE(J7)	0.0715	0.0661	0.0737	0.0762

Table 3. Joint torque estimation RMS errors during compliant contact (human hand).

Algorithm	ESO	KFMO	SH-AKFMO	VBKFMO-Q
RMSE(J1)	1.0266	0.8956	1.0259	0.286
RMSE(J2)	0.8367	0.6369	0.8364	0.4485
RMSE(J3)	1.3626	1.2528	1.3612	1.0145
RMSE(J4)	1.0637	0.9181	1.0622	0.4996
RMSE(J5)	0.1237	0.0791	0.1244	0.042
RMSE(J6)	0.1017	0.0876	0.1439	0.0968
RMSE(J7)	0.0715	0.0661	0.0737	0.0762

Root-mean-square error analysis demonstrates that, in line with the preceding two experimental scenarios, the proposed algorithm consistently yields lower RMS errors, with torque trajectories exhibiting both smoothness and rapid dynamic response. A consistent observation across all three experiments is that, given end-effector contact, torque estimations and corresponding RMS errors progressively smooth with increasing joint indices. This phenomenon reflects that, under sufficiently accurate modeling, joints distal to the collision site are subject to attenuated perturbations, thereby providing a rigorous basis for subsequent collision isolation and identification.

4. Conclusions

In this work, an adaptive external torque estimation algorithm for robotic manipulators was proposed, capable of automatically adjusting the system process noise covariance matrix. First, commonly used generalized momentum observers were reviewed and derived based on the robot dynamics model to address human-robot safety requirements. Then, an external torque estimation method grounded in generalized momentum was described, and the effects of covariance mismatch on Kalman filter estimation were analyzed. To handle practical situations where the system process noise is unknown and difficult to model, a method was developed to adaptively adjust the process noise covariance matrix using historical data, which was then incorporated into the Kalman filter recursive equations to improve estimation accuracy. The algorithm was validated experimentally on the Franka manipulator under three scenarios: static with no external torque, collision with a rigid object, and brief contact with a compliant object. The results demonstrate that the proposed method provides fast and accurate external torque estimates without additional sensors, offering a reliable basis for safe operation in human-robot collaborative tasks.