A physically consistent stiffness formulation for contact-rich manipulation

Abstract

Ensuring symmetric stiffness in impedance-controlled robots is crucial for physically meaningful and stable interaction in contact-rich manipulation. Conventional approaches neglect the change of basis vectors in curved spaces, leading to an asymmetric joint-space stiffness matrix that violates passivity and conservation principles. In this work, we derive a physically consistent, symmetric joint-space stiffness formulation directly from the task-space stiffness matrix by explicitly incorporating Christoffel symbols. This correction resolves long-standing inconsistencies in stiffness modeling, ensuring energy conservation and stability. We validate our approach experimentally on a robotic system, demonstrating that omitting these correction terms results in significant asymmetric stiffness errors. Our findings bridge theoretical insights with practical control applications, offering a robust framework for stable and interpretable robotic interactions.

Introduction

In the realm of robotics, the concept of “controller design in the physical domain” (Sharon et al., 1989, 1991) and the associated methodology of “control by interconnection” (Stramigioli, 2001; Van der Schaft, 2016) emphasize that robot controllers should be more than mere signal processors. Instead, they should have a physical interpretation (Lachner, 2022), which is especially important for robots that physically interact with the environment (Dietrich and Hogan, 2022; Hogan, 1988, 2022). This paper delves into impedance control (Hogan, 1984) during physical interaction, specifically focusing on the symmetry of the stiffness matrix and its role in ensuring passive physical equivalent robot controllers.

Passivity is a fundamental property for ensuring coupled stability when interacting with arbitrary passive objects (Colgate and Hogan, 1988). Stability in robotics can be achieved by monitoring and controlling the energy supplied by the controller (Colgate and Hogan, 1987, 1988; Stramigioli, 2015). In impedance-controlled robots, this monitoring is particularly straightforward, as energy is stored in virtual springs (potential energy) and transferred into kinetic energy during movement (Lachner et al., 2021).

During physical interaction, stiffness plays a crucial role, as it defines how energy is stored and exchanged between the robot and its environment. Task-space stiffness determines interaction forces due to contact, which is especially important at low frequencies (e.g., steady-state). Ensuring symmetry of the stiffness matrix in impedance-controlled robots is crucial for passive physical equivalent control, as an asymmetric stiffness matrix implies non-conservative force fields that contradict fundamental physical principles and violate passivity.

Many prior studies have addressed task-space stiffness asymmetry theoretically (Chen and Kao, 2000; Ciblak and Lipkin, 1994; Howard et al., 1998; Žefran et al., 1999; Žefran and Kumar, 1996; Zefran and Kumar, 1997, 1998), and empirically (Lee et al., 2025; Woolfrey et al., 2024) in the context of robotic manipulation tasks. This work builds upon these insights to present a physically consistent resolution that guarantees symmetric stiffness.

In this paper, we demonstrate how to derive a symmetric stiffness matrix in joint-space coordinates. By explicitly accounting for the robot’s kinematic changes under external forces, we identify two asymmetric components whose sum results in a symmetric stiffness. To our knowledge, this is the first work to derive a physically consistent symmetric kinematic stiffness, incorporating correction terms via Christoffel symbols. Unlike previous studies concluding with an asymmetric stiffness matrix, our approach ensures energy conservation and can be effectively used for robot control during physical interactions. To validate our approach, we implement our method on a real robot and provide open-source code for users to apply our findings.

Methods

In this paper, we employ tensor notation and the Einstein summation convention to explicitly consider basis vectors (Appendix A). All notations related to twists ξ ∈ se (3) and wrenches F ∈ se^⋆(3) can be found in Appendix B.

Affine connection and Christoffel symbols

An affine connection provides a systematic way to compare vectors at different points on a manifold, defining how basis vectors evolve along coordinate directions and enabling differentiation in curved spaces. The Christoffel symbols ${\mathrm{\Gamma }}_{ij}^k$ are the connection coefficients associated with an affine connection, capturing how basis vectors change as they move along the manifold.

A connection can be symmetric or asymmetric, depending on whether it has torsion. In a symmetric connection, the Christoffel symbols satisfy ${\mathrm{\Gamma }}_{ij}^k={\mathrm{\Gamma }}_{ji}^k$ , meaning the basis change is independent of the order of differentiation. In contrast, an asymmetric connection introduces torsion, leading to Christoffel symbols that do not necessarily exhibit this symmetry.

To perform differentiation on curved manifolds, it is essential to account for how the basis vectors e _i evolve along coordinate directions ξ ^j . The vectors e _i represent basis vector fields of the tangent space, which describe unit directions at different points on the manifold. In contrast, the coordinate-induced basis vectors ∂/∂ξ ^j arise from the coordinate chart and describe differentiation along coordinate lines. The Christoffel symbols quantify how the basis vectors change due to curvature, ensuring that the differentiation aligns with the underlying geometry (Figure 1).

Figure 1.

Chart map of a toroidal manifold patch (highlighted in purple). The coordinate-induced basis vectors ∂/∂ξ¹ and ∂/∂ξ² at different locations illustrate the rotation of the basis due to the manifold’s curvature. This necessitates correction terms (Christoffel symbols) to properly account for the basis change when computing derivatives.

The dependence of basis vectors e _i on the coordinate components ξ ^j leads to changes that can be expressed as

$\frac{\partial e_i}{\partial {\xi }^j}={\mathrm{\Gamma }}_{ij}^m{e}_m.$

(1)

Equation (1) shows that the Christoffel symbols ${\mathrm{\Gamma }}_{ij}^m$ describe the rate at which the basis vector e _i changes as we move along the coordinate ξ ^j . This change is expressed in terms of the basis vector e _m , effectively measuring how the basis itself changes due to the curvature of the space.

In a Cartesian space with zero curvature, the Christoffel symbols vanish, that is, ${\mathrm{\Gamma }}_{ij}^m=0$ . This is the case for pure translations, which are free vectors—shifting the origin does not alter their intrinsic properties, such as direction and magnitude.

However, in curved spaces, the Christoffel symbols become essential, as the basis vectors themselves change with position. A key example is spatial rotations, which belong to the Lie group SO(3), a curved and nonlinear manifold. Unlike Cartesian vectors, rotations do not follow standard vector operations: vector addition is not defined, composition is non-commutative (R₁R₂ ≠ R₂R₁). These same principles extend to directional derivatives, where the order of differentiation matters due to the underlying geometry.

This leads to a fundamental property of the Christoffel symbols, which follows from the non-commutativity of directional derivatives:

$\frac{\partial e_i}{\partial {\xi }^j}-\frac{\partial e_j}{\partial {\xi }^i}=C_{ij}^m{e}_m.$

Here, the structure coefficients $C_{ij}^m$ quantify the non-commutativity of the coordinate-induced basis vectors. In the context of robotic control, these coefficients correspond to the Lie algebra se (3), whose explicit form is provided in Appendix C.

This directly implies that the Christoffel symbols satisfy:

${\mathrm{\Gamma }}_{ij}^m-{\mathrm{\Gamma }}_{ji}^m=C_{ij}^m.$

(2)

Thus, the presence of non-zero structure coefficients $C_{ij}^m$ introduces an asymmetry in the Christoffel symbols, reflecting the geometric structure of the coordinate basis. This asymmetry influences how differential changes in position and orientation affect the local basis vectors, which in turn influences how twists and wrenches are transformed under coordinate changes in task-space.

Since stiffness is defined as the rate of change of a wrench with respect to spatial displacement, it is inherently tied to how wrenches transform under coordinate variations. As we will see in the next section, the Christoffel symbols appear explicitly in the stiffness formulation.

In this paper, we adopt the Kinematic Connection, originally introduced by Žefran et al. (1999) for body-fixed coordinates (also referred to as the left-invariant connection). A detailed derivation of this connection is presented in Appendix D.

Task-space stiffness

The task-space stiffness matrix K _ij characterizes how wrenches F _i change in response to spatial displacements ξ ^j . It is defined as

$K_{ij}=\frac{\partial F_i}{\partial {\xi }^j}.$

(3)

The wrench components F _i can be defined from the gradient of the potential energy function $\mathcal{U}$ (Figure 2(b)):

$F_i=\frac{\partial \mathcal{U}}{\partial {\xi }^j}$

(4)

Figure 2.

(a) Quadratic potential function $\mathcal{U}({\xi }^1,{\xi }^2)$ with coordinate-induced basis vectors ∂/∂ξ¹ and ∂/∂ξ² at different locations. The global minimum is marked by the black point. (b) At the green point, the presence of an external force (negative gradient $\partial /\partial {\xi }^k(\mathcal{U})$ ) induces a basis shift and rotation. This basis change, due to the curvature of the potential field, requires correction terms (Christoffel symbols).

In coordinate-free notation, the wrench can be described as the basis vector field e _i , acting as a directional derivative:

$F_i=e_i(\mathcal{U})$

(5)

Here, the notation $e_i(\mathcal{U})$ denotes the application of the directional derivative e _i to the scalar field $\mathcal{U}$ , rather than a product.

Substituting the expression for F _i from Equations (5) into (3), we obtain:

$K_{ij}=\frac{\partial }{\partial {\xi }^j}\left(e_i(\mathcal{U})\right).$

(6)

Applying the Leibniz rule for differentiation, we expand the derivative as

$K_{ij}=\left(\frac{\partial e_i}{\partial {\xi }^j}\right)(\mathcal{U})+e_i\left(\frac{\partial \mathcal{U}}{\partial {\xi }^j}\right).$

(7)

The first term, $\left(\partial e_i/\partial {\xi }^j\right)(\mathcal{U})$ , accounts for the variation of the basis vector field e _i with respect to coordinate ξ ^j . The second term, $e_i\left(\partial \mathcal{U}/\partial {\xi }^j\right)$ , represents the directional derivative of the gradient of $\mathcal{U}$ along e _i .

This formulation explicitly considers changes in the basis vectors, ensuring a physically meaningful representation of stiffness in curved spaces.

This step is crucial in the derivation, as the term ∂e _i /∂ξ ^j explicitly captures how the basis vectors e _i must be adjusted when shifted along the curved coordinate components ξ ^j . An illustrative example is shown in Figure 2: an external force prevents the system from reaching the global minimum of the potential energy function. Instead, due to the curvature of the potential energy landscape, the local coordinate-induced basis is rotated relative to the global minimum. This correction is systematically encoded by the Christoffel symbols, as shown in Section 2. While conventional matrix formulations neglect changes in the basis, tensor notation explicitly accounts for them, ensuring a complete representation of the system’s behavior.

From the definition of the Christoffel Symbols (Equation (1)), we know that

$\frac{\partial e_i}{\partial {\xi }^j}={\mathrm{\Gamma }}_{ij}^m\frac{\partial }{\partial {\xi }^m}.$

(8)

Using the coordinate-induced basis representation e _i = ∂/∂ξ ⁱ , we rewrite Equation (7) as

$K_{ij}=\frac{{\partial }^2\mathcal{U}}{\partial {\xi }^i\partial {\xi }^j}+{\mathrm{\Gamma }}_{ij}^m\frac{\partial \mathcal{U}}{\partial {\xi }^m}.$

(9)

Since the wrench components are defined as $F_m=\partial \mathcal{U}/\partial {\xi }^m$ , we arrive at the final expression:

$K_{ij}=\frac{{\partial }^2\mathcal{U}}{\partial {\xi }^i\partial {\xi }^j}+{\mathrm{\Gamma }}_{ij}^m{F}_m.$

(10)

As can be seen, the task-space stiffness matrix consists of two components: the second partial derivatives of the potential energy function and a correction terms, involving the external wrench F _m , which arises due to coordinates basis change.

Problem statement

Conventional matrix formulations overlook the change of basis vectors, leading to the conclusion that the task-space stiffness matrix is solely given by the second partial derivatives of the potential energy function. This formulation inherently yields a symmetric stiffness matrix.

By explicitly accounting for basis changes, Equation (10) introduces a correction term ${\mathrm{\Gamma }}_{ij}^m{F}_m$ which depends on the Christoffel symbols. This means that the symmetry of the task-space stiffness K _ij is directly influenced by the symmetry properties of the Christoffel symbols themselves.

The asymmetry of the task-space stiffness matrix in robotic control has been extensively reported in the literature (Chen and Kao, 2000; Ciblak and Lipkin, 1994; Howard et al., 1998; Lee et al., 2025; Woolfrey et al., 2024; Žefran et al., 1999; Žefran and Kumar, 1996; Zefran and Kumar, 1997, 1998).

An asymmetric stiffness is not usable for passive physical equivalent controllers, as it is equivalent to implementing a perpetual motion machine of the first kind, as will be shown in the following.

The task-space stiffness K _ij can be decomposed into symmetric and anti-symmetric (skew-symmetric) components:

$K_{ij}=\underset{\text{Symmetric\hspace{0.17em}part}{K}^S}{\underbrace{\frac{1}{2}(K_{ij}+K_{ji})}}+\underset{\text{Antisymmetric\hspace{0.17em}part}{K}^A}{\underbrace{\frac{1}{2}(K_{ij}-K_{ji})}}.$

(11)

The antisymmetric component evokes a force orthogonal to displacement. The (incremental) mechanical work done along a displacement trajectory is given by δW = F _i dξ ⁱ . When the stiffness matrix includes a skew-symmetric component $K_{ij}^A$ , the work over a closed trajectory C of non-zero size becomes

$\mathrm{\Delta }W=\frac{1}{2}{\oint }_C{\xi }^i{K}_{ij}^A\mathrm{d}{\xi }^j.$

(12)

This integral is generally non-zero when $K_{ij}^A\ne 0$ . As a result, the system exhibits a net energy gain or loss per cycle, leading to energy accumulation or dissipation. Traversing this orbit multiple times would implement a perpetual motion machine of the first kind. No physical spring can do this.

A robot-implemented virtual spring may be able to implement this behavior. However, this has a clear disadvantage: violation of passivity. The virtual spring becomes non-conservative when $K_{ij}^A\ne 0$ , as indicated by ∮ _C F _i dξ ⁱ ≠ 0.

In the next section, we derive a physically meaningful, symmetric stiffness matrix in joint-space coordinates, which is obtained from the task-space stiffness matrix. By explicitly incorporating Christoffel symbols, we demonstrate how to recover a physically consistent stiffness formulation that accounts for basis vector changes, adheres to conservation principles, and ensures passive control in contact-rich manipulation.

Results

In this section, Greek letters denote joint-space coordinate components, while Latin indices represent task-space coordinate components, consistent with the notation in Section 2.

Theory

A joint-space stiffness matrix ${\mathcal{K}}_{\alpha \beta }$ is the gradient of joint torque τ _β with respect to joint displacement q ^α :

${\mathcal{K}}_{\alpha \beta }=\frac{\partial {\tau }_{\beta }}{\partial q^{\alpha }}.$

(13)

We can map the task-space wrench F _k to torques by the Jacobian map $b_{\beta }^k$ :

${\tau }_{\beta }=b_{\beta }^k{F}_k.$

(14)

Substituting Equations (14) into (13) yields

${\mathcal{K}}_{\alpha \beta }=\frac{\partial }{\partial q^{\alpha }}(b_{\beta }^k{F}_k).$

(15)

Applying the Leibniz rule for differentiation, this expands to

${\mathcal{K}}_{\alpha \beta }=\frac{\partial b_{\beta }^k}{\partial q^{\alpha }}{F}_k+b_{\beta }^k\frac{\partial F_k}{\partial q^{\alpha }}.$

(16)

The term ∂F _k /∂q ^α can be expanded using the chain rule:

$\frac{\partial F_k}{\partial q^{\alpha }}=\frac{\partial F_k}{\partial {\xi }^l}\frac{\partial {\xi }^l}{\partial q^{\alpha }}=\frac{\partial F_k}{\partial {\xi }^l}{a}_{\alpha }^l,$

(17)

where $a_{\alpha }^l$ represents the Jacobian mapping between the coordinate systems. Substituting Equations (17) into (16) yields:

$K_{\alpha \beta }=\frac{\partial b_{\beta }^k}{\partial q^{\alpha }}{F}_k+b_{\beta }^k{a}_{\alpha }^l\frac{\partial F_k}{\partial {\xi }^l}.$

(18)

The wrench component F _k is the partial derivative of the potential function $\mathcal{U}$ with respect to coordinate ξ ^k , representing the gradient in the k-th direction. In coordinate-free notation, the basis vector e _k of the tangent space acts as a directional derivative, expressed as $e_k(\mathcal{U})$ :

$F_k=\frac{\partial \mathcal{U}}{\partial {\xi }^k}=e_k(\mathcal{U}).$

(19)

Substituting the expression for F _k from Equations (19) into (18), we obtain:

$K_{\alpha \beta }=\frac{\partial b_{\beta }^k}{\partial q^{\alpha }}{F}_k+b_{\beta }^k{a}_{\alpha }^l\frac{\partial }{\partial {\xi }^l}\left(e_k(\mathcal{U})\right),$

(20)

where e _k is a vector field and $\mathcal{U}$ is a function. The differentiation of their product follows the Leibniz rule and we expand the derivative as

$K_{\alpha \beta }=\frac{\partial b_{\beta }^k}{\partial q^{\alpha }}{F}_k+b_{\beta }^k{a}_{\alpha }^l\left(\frac{\partial e_k}{\partial {\xi }^l}(\mathcal{U})+e_k\frac{\partial \mathcal{U}}{\partial {\xi }^l}\right).$

(21)

Equation (21) represents a key step in the derivation, as it explicitly accounts for the correction of basis vectors e _k along the curved coordinate components ξ ^l (Figure 2). This aspect is overlooked in conventional matrix formulations (Hogan, 1985; Mussa-Ivaldi et al., 1985; Mussa-Ivaldi and Hogan, 1991).

For a basis vectors, we have e _k = ∂/∂ξ ^k . Additionally, from Equation (1), we know that the basis derivatives satisfy

$\frac{\partial e_k}{\partial {\xi }^l}={\mathrm{\Gamma }}_{kl}^m\frac{\partial }{\partial {\xi }^m},$

(22)

which allows us to rewrite Equation (21) as

$K_{\alpha \beta }=\frac{\partial b_{\beta }^k}{\partial q^{\alpha }}{F}_k+b_{\beta }^k{a}_{\alpha }^l\left({\mathrm{\Gamma }}_{kl}^m\frac{\partial \mathcal{U}}{\partial {\xi }^m}+\frac{{\partial }^2\mathcal{U}}{\partial {\xi }^k\partial {\xi }^l}\right).$

(23)

Since $\partial \mathcal{U}/\partial {\xi }^m$ corresponds to the wrench component F _m and ${\partial }^2\mathcal{U}/\partial {\xi }^k\partial {\xi }^l$ represents the task-space stiffness matrix K _kl , the basis transformation yield three distinct contributions to the joint-space stiffness:

$K_{\alpha \beta }=\underset{\begin{array}{c}\text{Kinematic}\\ \text{stiffness}\end{array}}{\underbrace{\frac{\partial b_{\beta }^k}{\partial q^{\alpha }}{F}_k}}+\underset{\begin{array}{c}\text{Connection}-\text{induced}\\ \text{correction\hspace{0.17em}term}\end{array}}{\underbrace{b_{\beta }^k{a}_{\alpha }^l{\mathrm{\Gamma }}_{kl}^m{F}_m}}+\underset{\begin{array}{c}\text{Pullback\hspace{0.17em}of}\\ \text{task}-\text{space\hspace{0.17em}stiffness}\end{array}}{\underbrace{b_{\beta }^k{a}_{\alpha }^l{K}_{kl}}}$

(24)

This is the central finding of our paper. The term $(b_{\beta }^k{a}_{\alpha }^l{K}_{kl})$ is inherently symmetric, as it derives from the Hessian matrix, which possesses this property. In contrast, the kinematic stiffness $(\partial b_{\beta }^k/\partial q^{\alpha }{F}_k)$ can be asymmetric when rotations (i.e., moment components in F _k ) are involved. An example can be found in Appendix F. Crucially, we establish that their combination with the Christoffel term,

$\frac{\partial b_{\beta }^k}{\partial q^{\alpha }}{F}_k+b_{\beta }^k{a}_{\alpha }^l{\mathrm{\Gamma }}_{kl}^m{F}_m,$

(25)

is symmetric since we provide a connection that is consistent with the robot kinematics (see example in Appendix G).

This result is significant because it demonstrates that the inclusion of Christoffel symbols ensures symmetry in the joint-space stiffness matrix, a property essential for stable contact-rich manipulation. While Žefran et al. (1999) showed that the task-space stiffness in general asymmetric, our approach corrects this, ensuring a physically equivalent control framework that adheres to fundamental principles like energy conservation.

While this derivation is theoretical, the next section presents experimental validation, highlighting its direct impact on real-world robotic interaction.

Practical experiments

The experiments were conducted on a KUKA LBR iiwa with seven DOFs, using KUKA’s Fast Robot Interface (FRI) for torque control. Built-in gravity and Coriolis/centrifugal compensation remained active throughout. Kinematic and dynamic matrices were computed via the Exp[licit]-FRI interface. ¹ The external wrenches were acquired with an ATI Gamma force/torque transducer, mounted on the robot flange. To follow standard robot control conventions, we present matrix notation in this subsection.

We implemented a joint-space impedance controller, where the control torque $\mathit{\tau }\in {\mathbb{R}}^n$ was determined by joint-space stiffness $\mathcal{K}\in {\mathbb{R}}^{n\times n}$ , acting on the deviation between the zero-force configuration q ₀ and the joint configuration q , and joint-space damping $\mathcal{B}\in {\mathbb{R}}^{n\times n}$ , acting on the joint velocity $\dot{\mathit{q}}\in {\mathbb{R}}^n$ :

$\mathit{\tau }=\mathcal{K}({\mathit{q}}_0-\mathit{q})-\mathcal{B}\dot{\mathit{q}}.$

(26)

The joint-space damping matrix $\mathcal{B}$ was designed based on $\mathcal{K}$ and the robot’s inertia matrix $\mathcal{M}(\mathit{q})\in {\mathbb{R}}^{7\times 7}$ (Albu-Schaffer et al., 2003). The joint-space stiffness was computed using Equation (24), incorporating the Jacobian matrix $\mathit{J}(\mathit{q})\in {\mathbb{R}}^{6\times 7}$ , task-space stiffness $\mathit{K}\in {\mathbb{R}}^{6\times 6}$ , and correction terms from Christoffel symbols and external wrenches $(\mathbf{\Gamma }{\mathit{F}}_{\text{ext}})\in {\mathbb{R}}^{6\times 6}$ (Appendix E). In matrix form, this can be written as

$\mathcal{K}=\frac{\partial \mathit{J}{(\mathit{q})}^T}{\partial \mathit{q}}{\mathit{F}}_{\text{ext}}+\mathit{J}{(\mathit{q})}^T\left[\mathit{K}+\mathbf{\Gamma }{\mathit{F}}_{\text{ext}}\right]\mathit{J}(\mathit{q}).$

(27)

Two experiments were conducted: (1) the robot moving through a parkour course (Section 3.2.1), and (2) the robot wiping a bowl (section 3.2.2). Both experiments were performed using the proposed approach (with the correction term ΓF_ext) and the standard approach (without the correction term ΓF_ext), as convention in conventional matrix formulations (Hogan, 1985; Lee et al., 2025; Mussa-Ivaldi et al., 1985; Mussa-Ivaldi and Hogan, 1991; Woolfrey et al., 2024).

The robot controller’s sample time was set to 5 ms for both experiments. Two threads were implemented using the Boost C++ Library with the mutex option. The main thread computed the controller based on Equation (27), including all kinematic and dynamic matrices, and acquired the external wrenches. The second thread calculated the joint-space stiffness $\mathcal{K}$ (Equation (27)). The kinematic stiffness term ∂ J ( q ) ^T /∂ q was derived symbolically using Exp[licit]-Matlab ² (Lachner et al., 2024) and integrated into the real-time control code.

The mean computation time of the first thread across both experiments was 855 μs. The second thread required 60 μs without the correction terms and 80 μs with the correction terms.

The full code is publicly available in the paper’s GitHub repository. ³

Parkour

The first experiment was to move the robot end-effector through a parkour course (Figure 3), while exerting significant external wrenches, measured by force/torque transducer attached to the robot end-effector.

Figure 3.

Robot movement, visualized as numbered overlaid robot configurations during the parkour experiments.

The zero-force trajectory q ₀ was generated using imitation learning (Nah et al., 2024) and pre-processed using Exp[licit]-Matlab. The code can be found on the paper’s GitHub repository.

As can be seen in Figure 4, the robot trajectory for both approaches were nearly identical.

Figure 4.

Recorded joint configurations during the parkour experiments. Solid lines indicate trials with correction terms (Christoffel symbols), while dashed lines represent trials without. The nearly overlapping trajectories suggest minimal differences between the two approaches.

To assess the practical relevance of both approaches, we computed the average RMSE of the error term ( q ₀ − q ) across all joints. The RMSE was 2.40° for the standard approach and 2.45° for the approach that included the Christoffel symbols. These results indicate that both approaches produced very similar behavior, and that the proposed method does not negatively impact the robot’s trajectory.

The measured external forces and moments during both experimental trials can be seen in Figure 5.

Figure 5.

Measured external forces and moments from the ATI force-torque sensor during the parkour experiments, expressed in the robot base frame. No filtering was applied to the signals. The solid curve represents one representative trial with the proposed correction terms (Christoffel symbols), and the dashed curve corresponds to a trial without them. The numbered shaded regions indicate distinct obstacle interactions, as referenced in Figure 3.

The measured external forces ranged from −22 N to 20 N, and the measured external moments ranged from −5 Nm to 5 Nm. These forces and moments influenced the symmetry of the joint-space stiffness matrix. To quantify the asymmetry, the stiffness matrix was decomposed into symmetric and anti-symmetric components (Equation (11)).

Figure 6 shows the maximum singular values of both components over time. The results confirm that the correction term ensured the symmetry of the joint-space stiffness matrix (cf. zero anti-symmetric component in Figure 6). In contrast, omitting the correction term led to an anti-symmetric part with singular values reaching approximately 8.4 Nm/rad around 12 s. To put this into perspective, this corresponds to a joint torque error of 8.4 Nm for just a 1 rad (about 57°) deviation in joint position, comparable to holding a 10 kg weight at the end of a 8.6 cm lever arm.

Figure 6.

Maximum singular values of the symmetric (left) and anti-symmetric (right) parts of the joint-space stiffness matrix during the parkour experiments. The solid line represents the trial with correction terms (Christoffel symbols), and the dashed line without.

Wiping a bowl

The second robot experiment was to wipe a bowl in a semi-circular motion (Figure 7). We selected this experiment to showcase the effect of substantial external moments.

Figure 7.

Robot movement, visualized as overlaid robot configurations and directions during the bowl-wiping experiments.

Again, the correction terms didn’t have a major influence on the robot trajectory, as can be seen in Figure 8.

Figure 8.

Recorded joint configurations during the wiping experiments. Solid lines indicate trials with correction terms (Christoffel symbols), while dashed lines represent trials without. The nearly overlapping trajectories suggest minimal differences between the two approaches.

Compared to the Parkour experiment, the average RMSE of the error term ( q ₀ − q ) across all joints was higher but remained comparable for both approaches, with 8.11° for the standard approach and 8.16° for the approach including the Christoffel symbols.

The recorded external forces and moments during both experimental trials can be seen in Figure 9.

Figure 9.

Measured external forces and moments recorded by the ATI force-torque sensor during the wiping experiments, expressed in the robot base frame. No filtering was applied to the signals. The solid curve shows one representative trial using the proposed correction terms (Christoffel symbols), while the dashed curve corresponds to a trial without correction terms.

The applied moments at the end-effector, ranging between −11 Nm and 19 Nm, affected the symmetry of the joint-space stiffness matrix. To quantify the asymmetry, again, the stiffness matrix was decomposed into symmetric and anti-symmetric components (Equation (11)).

Figure 10 shows the maximum singular values of both components over time.

Figure 10.

Maximum singular values of the symmetric (left) and anti-symmetric (right) parts of the joint-space stiffness matrix during the wiping experiments. The solid line represents the trial with correction terms (Christoffel symbols), and the dashed line without.

The results confirm that the correction term ensured the symmetry of the joint-space stiffness matrix (cf. zero anti-symmetric component in Figure 10). In contrast, omitting the correction term led to an anti-symmetric part with singular values reaching approximately 18.5 Nm/rad around 1s.

This stiffness error is significant. For comparison, the torque required to tighten a standard M10 bolt with a hand wrench is about 20 Nm. The observed stiffness error corresponds approximately to that effort, highlighting the practical impact of neglecting the correction terms.

Discussion

This work demonstrates that ensuring the symmetry of the stiffness matrix is essential for passive physically equivalent control in contact-rich manipulation. Conventional matrix formulations neglect the change of basis vectors in curved spaces, leading to an asymmetric kinematic stiffness matrix in joint-space coordinates. By explicitly incorporating Christoffel symbols, we derived a symmetric stiffness formulation that adheres to conservation principles, ensuring passivity and stability.

The proposed approach resolves long-standing inconsistencies in task-space stiffness modeling, where prior work acknowledged asymmetries but lacked a physically consistent correction. Our experimental validation on a real robotic system confirms that including correction terms significantly impacts the system’s response, ensuring that the joint-space stiffness remains symmetric under external wrenches.

Although the robot trajectories were nearly identical for both approaches, the correction term reshapes the robot’s impedance behavior in response to external wrenches. This led to slight differences in the resulting interaction wrenches. For instance, in the parkour experiment (Figure 5), the approach without correction term exhibited higher external moments, such as m_x around 9 s and m_z around 12 s. In contrast, in the bowl-wiping experiment (Figure 9), adding the correction term changed the stiffness response, producing higher f_x, f_z, and torque values around 1 second.

In our experiments, we observed slight oscillatory behavior in the measured external forces and moments during contact. Further cross-correlation analysis revealed that these fluctuations were caused by stick–slip: forces led end-effector velocity with a lag of 0.5 s. While friction may affect the magnitude of the Christoffel correction, the symmetry of the joint-space stiffness matrix is preserved by construction in the proposed approach.

The proposed approach restores symmetry in the joint-space stiffness matrix by compensating for curvature effects that arise from configuration-dependent kinematics. This can be viewed as a counterpart to the Conservative Congruence Transformation (CCT) framework (Chen, 2003; Chen and Kao, 2002), which attributes asymmetry in Cartesian stiffness to the use of non-coordinate twist bases under SE(3). While CCT preserves energy consistency in task space, our formulation ensures geometric consistency in joint space. Both approaches could be unified under the principle of virtual work, providing a theoretical bridge between joint- and task-space formulations.

Recent learning-based approaches have focused on generating impedance behavior by learning symmetric positive-definite (SPD) stiffness matrices. These methods guarantee symmetry by construction, often by representing stiffness on Riemannian manifolds or projecting estimated matrices onto the SPD space (Abu-Dakka et al., 2018; Beik-Mohammadi et al., 2023; Jaquier et al., 2021; Sun and Figueroa, 2024). While these approaches draw on concepts from differential geometry and yield physically plausible stiffness estimates, they typically do not account for configuration-dependent kinematic effects, such as changes in the Jacobian and the resulting kinematic stiffness.

In contrast, our formulation using Christoffel symbols maintains stiffness symmetry not only by construction but also under coordinate transformations, ensuring physical consistency during contact-rich manipulation. This formulation could be integrated with learning-based methods to enable task-specific stiffness profiles that remain valid across varying kinematics.

By bridging theoretical foundations with practical implementation, our work provides a robust framework for impedance control that respects fundamental physical principles, paving the way for stable robotic interactions.

Conclusion

This paper presents a passive physically equivalent approach for contact-rich manipulation, resolving a long-standing issue of asymmetric stiffness in impedance-controlled robots. By explicitly accounting for the geometric properties of robot kinematics, we derived a symmetric stiffness formulation in joint-space coordinates, incorporating correction terms derived from Christoffel symbols. Unlike previous formulations that resulted in inherently asymmetric stiffness matrices, our method guarantees compliance with fundamental physical principles, including passivity and energy conservation.

We demonstrated that the presence of asymmetric stiffness terms leads to non-conservative force fields, violating passivity and introducing artificial energy sources into the system. Through theoretical analysis and experimental validation, we showed that including Christoffel symbol-based corrections ensures a symmetric stiffness matrix, preventing unphysical energy accumulation and enabling stable interaction forces in real-world robotic tasks.

Our experimental results on a KUKA LBR iiwa robot confirmed that neglecting these correction terms can lead to substantial stiffness errors, potentially compromising the stability of impedance-controlled robots in contact-rich environments. By integrating the proposed corrections, we restored the expected physical behavior, yielding stable and interpretable robotic interactions.