Skip to content

MuJoCo Dynamics – Single DOF Trace Validation

Overview

Internally, MuJoCo uses

M(q) q̈ = qfrc_actuator + qfrc_passive + qfrc_constraint − qfrc_bias

Each term is reconstructed independently and compared with MuJoCo's internal solver output.

1. System Equation

MuJoCo implements rigid-body dynamics in the form:

M(q) q̈ + C(q, q̇) = τ
Symbol Meaning
M(q) Mass matrix
C(q,q̇) Coriolis + centrifugal + gravity
τ All external generalized forces

In MuJoCo, τ is not just the actuator term — it is the sum of all non-inertial forces:
motors, springs/dampers, contact forces, joint limits, and equality constraints.

Rewritten in solver form:

M(q) q̈ = τ_actuator + τ_passive + τ_constraint − C(q,q̇)
Term Meaning
M(q) Mass matrix
τ_actuator Actuator torques
τ_passive Passive joint forces
τ_constraint Contact and constraint forces
τ_bias Coriolis + centrifugal + gravity

MuJoCo defines: τ_bias := C(q,q̇) — exactly the value returned by mj_rne() with q̈ = 0.

qfrc ≡ τ stands for “q-space forces” — forces in the space of generalized coordinates q.

2. RNEA – Recursive Newton-Euler Algorithm

RNEA computes inverse dynamics according to Featherstone (Newton/Euler, O(n)):

Question: Which joint torques τ are required to generate a motion (q, q̇, q̈)?

Algorithm:

  1. Forward pass (base → end effectors):
  2. propagate angular velocities
  3. propagate accelerations

  4. accumulate forces

  5. Backward pass (end effectors → base):

  6. accumulate forces and torques

If q̈ = 0 is set:

mujoco.mj_rne(model, data, 0, result)

Then: τ_bias = C(q, q̇) — that is, all velocity- and gravity-dependent forces (gravity, Coriolis, centrifugal).

3. Control Input – Actuator

For a position actuator:

τ = kp · (ctrl − q) − kv · q̇

MuJoCo stores this internally in affine form:

τ = gain · ctrl + b0 + b1·q + b2·q̇

MuJoCo rewrites this algebraically as:

τ = kp·ctrl − kp·q − kv·q̇
  = kp·ctrl + b1·q + b2·q̇    with b1 = −kp, b2 = −kv

The actuator force is then transformed into generalized coordinates:

qfrc_actuator = J_act^T · τ

Transmission Jacobian J_act: - Describes actuator length as a function of DOF motion: J_act ∈ R^(n_act × n_dof) - l̇ = J_act · q̇ - J_act^T transforms force from actuator space into DOF space

For joint actuators (<position joint="...">): l = q, therefore J_act = 1
→ qfrc_actuator[dof] = τ_actuator

4. Passive Forces

Computed in mj_passive():

τ_passive = − stiffness · (q − q_ref)   # Hooke: q > q_ref → Feder zieht zurück
          − damping · q̇                 # viscous damping, opposite to velocity
          − frictionloss · sign(q̇)      # Coulomb friction, always braking

Explanation:

  • Stiffness: If the joint is displaced, a restoring force acts on it
  • Damping: Viscous damping proportional to velocity
  • Coulomb friction: q̇ > 0 → negative torque, q̇ < 0 → positive torque

If frictionloss ≤ 0, MuJoCo ignores the Coulomb term.
All terms are dissipative — they reduce energy.

Parameters come from the XML: <joint damping="..." stiffness="..." frictionloss="..."/>

5. Constraint Forces

Constraints (contacts, joint limits, equality constraints) generate Lagrange multipliers λ.

MuJoCo solves the KKT system:

[ M   Jᵀ ] [ q̈ ]   [ τ ]
[ J   0  ] [ λ  ] = [ b ]
Symbol Meaning
J Constraint Jacobian
λ Lagrange multipliers (reaction forces)
b Constraint velocity / stabilization

The generalized constraint torque:

τ_constraint = Jᵀ · λ

qfrc_constraint contains all constraint reactions:

  • contact forces (normal + friction)
  • joint limits
  • equality constraints (welds, tendons, closed kinematic chains)
  • friction constraints at contacts

qfrc_constraint is therefore not “just contact forces,” but the projection of all constraint reaction forces into joint space.

Constraint Jacobian in Detail

A constraint is an equation: ϕ(q) = 0

The constraint Jacobian is its derivative:

J(q) = ∂ϕ(q) / ∂q

For a small DOF delta δq: δϕ = J(q) · δq

  • Each row → one constraint
  • Each column → one DOF
  • J describes how DOF motion affects constraint violation
  • Jᵀ distributes the constraint force back to the DOFs: τ_constraint = Jᵀ · λ

Every constraint that prevents motion must generate a reaction force — this is projected into generalized coordinates as τ_constraint = Jᵀ λ.

5.1 Contact Forces

Each contact generates multiple constraint equations depending on the friction model:

  • elliptic friction → 3 constraints
  • pyramidal friction → 4–5 constraints

These generate corresponding λ entries in data.efc_force.

Cartesian contact force:

mujoco.mj_contactForce(model, data, contact_id, force_array)

Generalized contribution: 

τ_contact = J_contact^T · λ_contact

5.2 Joint Limits

If a joint reaches its limit (φ(q) ≥ 0), a constraint force λ_limit is generated:

τ_limit = J_limit^T · λ_limit

It is included in τ_constraint.

5.3 Equality Constraints

Enforce φ(q) = 0. Examples: weld constraints, loop closures, tendon-length constraints.
Contribution: τ_equality = J_equality^T · λ_equality

6. Total Torque & Acceleration

τ_total = qfrc_actuator + qfrc_passive + qfrc_constraint − qfrc_bias

q̈ = M⁻¹(q) · τ_total

All generalized forces are summed → equation of motion M(q)q̈ = τ_total → the solver solves for q̈.

Validation:

  1. Extract mass matrix: mj_fullM
  2. Solve the linear system with NumPy
  3. Compare against data.qacc

Expected numerical deviation: ~1e−2 to 1e−3

7. Integration

With the implicitfast-Integrator:

v_{t+1} = v_t + a_{t+1} · dt
q_{t+1} = q_t + v_{t+1} · dt

This is implicit integration, not explicit Euler.
Therefore, small deviations arise when comparing against v_t + a_t·dt.

8. Indexing Rule

Array Index
qvel, qacc DOF index
qpos qpos index (jnt_qposadr)

DOF index ≠ qpos index for free joints or quaternions.

9. Expected Numerical Tolerances

Term Typical deviation
Actuator ~1e−7
Bias (RNEA) ~1e−3
Constraint solver-dependent
Acceleration ~1e−2
Position ~1e−4