@@ -10,6 +10,7 @@
 #include <momentum/character/skeleton_state.h>
 #include <momentum/character_sequence_solver/acceleration_sequence_error_function.h>
 #include <momentum/character_sequence_solver/jerk_sequence_error_function.h>
+#include <momentum/character_sequence_solver/joint_to_joint_sequence_error_function.h>
 #include <momentum/character_sequence_solver/model_parameters_sequence_error_function.h>
 #include <momentum/character_sequence_solver/sequence_solver.h>
 #include <momentum/character_sequence_solver/state_sequence_error_function.h>
@@ -21,6 +22,7 @@
 #include <pybind11/stl.h>
 #include <Eigen/Core>
 
+#include <pymomentum/python_utility/eigen_quaternion.h>
 #include <pymomentum/solver2/solver2_sequence_error_functions.h>
 #include <pymomentum/solver2/solver2_utility.h>
 
@@ -272,6 +274,92 @@ the gradient contribution is set to zero.
           &mm::VelocityMagnitudeSequenceErrorFunction::getTargetSpeeds,
           "Returns the per-joint target speeds.");
 
+  py::class_<
+      mm::JointToJointSequenceErrorFunction,
+      mm::SequenceErrorFunction,
+      std::shared_ptr<mm::JointToJointSequenceErrorFunction>>(
+      m, "JointToJointSequenceErrorFunction")
+      .def(
+          py::init<>([](const mm::Character& character,
+                        float weight,
+                        float positionWeight,
+                        float orientationWeight) {
+            validateWeight(weight, "weight");
+            validateWeight(positionWeight, "position_weight");
+            validateWeight(orientationWeight, "orientation_weight");
+
+            auto result = std::make_shared<mm::JointToJointSequenceErrorFunction>(character);
+            result->setWeight(weight);
+            result->setPositionWeight(positionWeight);
+            result->setOrientationWeight(orientationWeight);
+
+            return result;
+          }),
+          R"(A sequence error function that penalizes frame-to-frame changes in the relative
+transform between two joints.
+
+For each constraint with source joint A and reference joint B, this function computes
+the relative transform at two consecutive frames and penalizes any change between them:
+
+- Position error: ``relPos_{i+1} - relPos_i``
+- Orientation error: ``relRot_{i+1} - relRot_i``
+
+This keeps the relative transform (e.g. a hand's grip on an object) stable from frame
+to frame.
+
+:param character: The character to use.
+:param weight: The weight of the error function. Defaults to 1.0.
+:param position_weight: The weight of the position error. Defaults to 1.0.
+:param orientation_weight: The weight of the orientation error. Defaults to 1.0.)",
+          py::arg("character"),
+          py::kw_only(),
+          py::arg("weight") = 1.0f,
+          py::arg("position_weight") = 1.0f,
+          py::arg("orientation_weight") = 1.0f)
+      .def(
+          "add_constraint",
+          [](mm::JointToJointSequenceErrorFunction& self,
+             size_t sourceJoint,
+             size_t referenceJoint,
+             float constraintWeight,
+             const std::string& name) {
+            self.addConstraint(sourceJoint, referenceJoint, constraintWeight, name);
+          },
+          R"(Add a constraint between two joints.
+
+:param source_joint: Index of the source joint (the joint whose relative pose is tracked).
+:param reference_joint: Index of the reference joint (the frame of reference).
+:param weight: Weight for this constraint. Defaults to 1.0.
+:param name: Name for debugging. Defaults to empty string.)",
+          py::arg("source_joint"),
+          py::arg("reference_joint"),
+          py::kw_only(),
+          py::arg("weight") = 1.0f,
+          py::arg("name") = "")
+      .def(
+          "clear_constraints",
+          &mm::JointToJointSequenceErrorFunction::clearConstraints,
+          "Clears all constraints from the error function.")
+      .def_property_readonly(
+          "num_constraints",
+          &mm::JointToJointSequenceErrorFunction::getNumConstraints,
+          "Returns the number of constraints.")
+      .def_property(
+          "weight",
+          &mm::JointToJointSequenceErrorFunction::getWeight,
+          &mm::JointToJointSequenceErrorFunction::setWeight,
+          "The weight applied to the error function.")
+      .def_property(
+          "position_weight",
+          &mm::JointToJointSequenceErrorFunction::getPositionWeight,
+          &mm::JointToJointSequenceErrorFunction::setPositionWeight,
+          "The weight of the position error.")
+      .def_property(
+          "orientation_weight",
+          &mm::JointToJointSequenceErrorFunction::getOrientationWeight,
+          &mm::JointToJointSequenceErrorFunction::setOrientationWeight,
+          "The weight of the orientation error.");
+
   py::class_<
       mm::StateSequenceErrorFunction,
       mm::SequenceErrorFunction,
 
@@ -3253,3 +3253,85 @@ def test_transform_pose_single(self) -> None:
         np.testing.assert_allclose(
             skel_out[:, :3], skel_init[:, :3] + transform[:3], atol=1e-2
         )
+
+    def test_joint_to_joint_sequence_error_function(self) -> None:
+        """Test JointToJointSequenceErrorFunction creation and basic behavior."""
+        character = pym_test_utils.create_test_character(num_joints=4)
+        n_params = character.parameter_transform.size
+
+        # Create error function
+        ef = pym_solver2.JointToJointSequenceErrorFunction(
+            character,
+            weight=1.0,
+            position_weight=2.0,
+            orientation_weight=0.5,
+        )
+        self.assertEqual(ef.weight, 1.0)
+        self.assertEqual(ef.position_weight, 2.0)
+        self.assertEqual(ef.orientation_weight, 0.5)
+
+        # Add constraints
+        ef.add_constraint(source_joint=2, reference_joint=0)
+        ef.add_constraint(
+            source_joint=1,
+            reference_joint=3,
+            weight=0.5,
+            name="test_constraint",
+        )
+        self.assertEqual(ef.num_constraints, 2)
+
+        # Same parameters on both frames: zero error
+        model_params_same = np.zeros((2, n_params), dtype=np.float32)
+        error = ef.get_error(model_params_same)
+        self.assertAlmostEqual(error, 0.0, places=6)
+
+        # Different parameters on both frames: non-zero error
+        rng = np.random.RandomState(42)
+        model_params_diff = rng.randn(2, n_params).astype(np.float32) * 0.1
+        error = ef.get_error(model_params_diff)
+        self.assertGreater(error, 0.0)
+
+        # Clear constraints
+        ef.clear_constraints()
+        self.assertEqual(ef.num_constraints, 0)
+
+        # After clearing, error should be zero
+        error = ef.get_error(model_params_diff)
+        self.assertAlmostEqual(error, 0.0, places=6)
+
+    def test_joint_to_joint_sequence_error_function_solve_reduces_error(
+        self,
+    ) -> None:
+        """Test solve_sequence() can reduce joint-to-joint sequence error."""
+        character = pym_test_utils.create_test_character(num_joints=4)
+        n_params = character.parameter_transform.size
+        n_frames = 2
+
+        rng = np.random.RandomState(7)
+        model_params_init = np.zeros((n_frames, n_params), dtype=np.float32)
+        model_params_init[1] = (rng.randn(n_params) * 0.25).astype(np.float32)
+
+        rel_error = pym_solver2.JointToJointSequenceErrorFunction(
+            character,
+            weight=10.0,
+            position_weight=10.0,
+            orientation_weight=10.0,
+        )
+        rel_error.add_constraint(source_joint=2, reference_joint=0)
+        initial_error = rel_error.get_error(model_params_init)
+        self.assertGreater(initial_error, 1e-4)
+
+        solver_function = pym_solver2.SequenceSolverFunction(character, n_frames)
+        solver_function.set_enabled_parameters(np.ones(n_params, dtype=bool))
+        solver_function.add_sequence_error_function(0, rel_error)
+
+        solver_options = pym_solver2.SequenceSolverOptions()
+        solver_options.max_iterations = 25
+        solver_options.regularization = 1e-6
+
+        model_params_final = pym_solver2.solve_sequence(
+            solver_function, model_params_init, solver_options
+        )
+
+        final_error = rel_error.get_error(model_params_final)
+        self.assertLess(final_error, initial_error * 1e-3)