Rigid-Body Dynamics (Newton-Euler & Lagrangian) — Robotics Reference

1. At a glance

Robot dynamics is the study of the relationship between joint torques τ (the only thing a motor can produce) and joint accelerations q̈ (the thing the world will respond with). For a serial-chain manipulator with n joints, every algorithm — control, simulation, force estimation, MPC, RL world models — eventually reduces to one matrix equation, the manipulator equation:

$M(q)\ddot{q}+C(q,\dot{q})\dot{q}+g(q)=\tau +J^{\top }(q)F_{\text{ext}}$

• M(q) — joint-space mass / inertia matrix (n × n, symmetric, positive-definite)
• C(q, q̇)·q̇ — Coriolis + centripetal vector (n × 1)
• g(q) — joint-space gravity vector (n × 1)
• τ — applied joint torques (motor commands minus drivetrain losses)
• J(q) — geometric Jacobian (6 × n)
• F_ext — external wrench applied at the end-effector (6 × 1)

The same equation appears in industrial-arm torque control, cobot gravity compensation, humanoid centroidal-momentum control, drone attitude dynamics, and RL simulators — only the dimension and the contact constraints change.

Where it sits in the design stack: kinematics gives you q ↔ pose; dynamics gives you τ ↔ q̈ and is the prerequisite for any model-based controller, any physics simulator, and any model-based RL pipeline.

Common variants:

• Fixed-base serial (industrial arm, cobot) — root link is bolted down; q ∈ ℝⁿ.
• Floating-base (humanoid, quadruped, drone) — base has 6 unactuated DOF (SE(3)); q ∈ ℝ⁶⁺ⁿ.
• Constrained / closed-chain (parallel manipulator, robot in contact) — add Lagrange multipliers λ for J_c·q̈ + J̇_c·q̇ = 0.
• Hybrid dynamics — some joints torque-controlled, others position-controlled; mixed inverse/forward solve.

First ask before applying: Do you need inverse dynamics (τ from q̈, e.g. computed-torque control, MPC, ID-based feedforward — O(n) via RNEA) or forward dynamics (q̈ from τ, e.g. simulation, predictive rollouts — O(n) via ABA)? The two questions take different algorithms; running the wrong one is the most common newcomer mistake.

2. First principles

Rigid body in SE(3)

A rigid body has 6 DOF: 3 translation + 3 rotation. Its state is the pair (p, R) ∈ SE(3) where p ∈ ℝ³ is the position of the body-fixed origin and R ∈ SO(3) is the orientation (rotation matrix). Its velocity is the spatial twist V = (v, ω) ∈ ℝ⁶ — linear velocity v of the body origin and angular velocity ω. Its inertia is the pair (m, I_com) — scalar mass m (kg) and the 3 × 3 inertia tensor I_com about the centre of mass (kg·m²).

The inertia tensor at an arbitrary frame origin O is related to the COM tensor through the parallel-axis (Steiner) theorem:

$I_O=I_{\text{com}}+m(\Vert {\mathbf{r}}_{O\to \text{com}}{\Vert }^2{I}_3-{\mathbf{r}}_{O\to \text{com}}{\mathbf{r}}_{O\to \text{com}}^{\top })$

URDF / SDF / MJCF inertia tags must specify the frame the tensor is expressed in — a frequent source of silent 10–30 % errors when authors confuse link-origin and COM frames.

A rigid body has 10 inertial parameters: mass m, three COM coordinates (c_x, c_y, c_z), and six independent entries of the symmetric I_com (I_xx, I_yy, I_zz, I_xy, I_xz, I_yz). For an n-link manipulator, the standard parameter set has 10n entries; the number of identifiable (base) parameters is lower (typically 6–8 per link for a serial 6R arm; Khalil & Dombre 2002).

Newton’s second law for a rigid body

Translation of the centre of mass:

$\sum \mathbf{F}=m{\ddot{\mathbf{p}}}_{\text{com}}$

Rotation about the centre of mass (Euler’s equation, body frame):

$\sum {\mathbf{\tau }}_{\text{com}}=I_{\text{com}}\mathbf{\alpha }+\mathbf{\omega }\times (I_{\text{com}}\mathbf{\omega })$

The cross-product term ω × (I·ω) is the gyroscopic torque — it vanishes for ω = 0 or for a spherically symmetric inertia, and it dominates fast-spinning bodies like reaction wheels and helicopter rotors.

Spatial-vector form (Featherstone)

Stacking translation and rotation into a single 6-vector eliminates duplicate bookkeeping. Define the spatial twist V = [ω; v] and spatial wrench F = [τ; f]. The 6 × 6 spatial inertia is

$\mathcal{I}=\left[\begin{array}{cc}{I}_{\text{com}}+m[\mathbf{c}{]}_{\times }^{\top }[\mathbf{c}{]}_{\times }& m[\mathbf{c}{]}_{\times }\\ m[\mathbf{c}{]}_{\times }^{\top }& m{I}_3\end{array}\right]$

where [c]_× is the 3 × 3 skew-symmetric matrix of the COM-offset vector c. Newton-Euler in spatial form becomes the deceptively compact

$\mathcal{F}=\mathcal{I}\dot{\mathcal{V}}+\mathcal{V}{\times }^{\ast }(\mathcal{I}\mathcal{V})$

The cross-product operator ×* is the spatial-force cross product (Featherstone 2008 §2.9). This single equation, applied recursively through a kinematic chain, is what RNEA, CRBA, and ABA all expand into code.

Newton-Euler recursive algorithm (RNEA)

For a serial chain of n rigid bodies connected by joints, Luh, Walker & Paul (1980) showed that ΣF = m·a applied link-by-link can be turned into a two-pass recursion:

• Forward pass (base → tip): propagate position, velocity, and acceleration outward, computing link spatial velocities V_i and accelerations A_i from joint values (q_i, q̇_i, q̈_i).
• Backward pass (tip → base): propagate forces and moments inward, computing the wrench on each link and projecting it onto the joint axis to extract the joint torque τ_i.

Cost: O(n) per call. This is the workhorse of inverse dynamics — Pinocchio, RBDL, Drake, MuJoCo all use a Featherstone-style RNEA at their core.

The forward pass propagates spatial velocity V_i = V_{p(i)} ⊕ S_i q̇_i and spatial acceleration A_i = A_{p(i)} + S_i q̈_i + V_i × S_i q̇_i, where p(i) is the parent of link i and S_i is the 6 × 1 joint motion subspace (a unit vector along the joint axis for a 1-DOF revolute or prismatic joint). The backward pass computes the net wrench F_i = 𝓘_i A_i + V_i ×* (𝓘_i V_i) − F_i^ext, then accumulates child wrenches up the chain and extracts τ_i = S_i^⊤ F_i.

Lagrangian formulation

The energy-based alternative defines the Lagrangian L = T − V where

$T=\frac{1}{2}{\dot{q}}^{\top }M(q)\dot{q},V=V_{\text{grav}}(q)+V_{\text{elastic}}(q)$

The Euler-Lagrange equations,

$\frac{d}{dt}\left(\frac{\partial L}{\partial {\dot{q}}_i}\right)-\frac{\partial L}{\partial q_i}={\tau }_i$

produce the manipulator equation directly, with M symmetric by construction and the Coriolis matrix derivable from Christoffel symbols of the first kind:

$C_{ij}(q,\dot{q})={\sum }_{k=1}^n\frac{1}{2}\left(\frac{\partial M_{ij}}{\partial q_k}+\frac{\partial M_{ik}}{\partial q_j}-\frac{\partial M_{jk}}{\partial q_i}\right){\dot{q}}_k$

Equivalent to Newton-Euler for any rigid-body system; chosen by hand-derivation because it eliminates internal constraint forces automatically.

The two formulations are interchangeable: Newton-Euler is preferred for algorithmic implementation (the recursion structure is natural), Lagrangian is preferred for hand derivation (kinetic and potential energies are scalar and additive) and for systems with non-trivial holonomic constraints (which can be folded into a reduced generalised-coordinate set).

Hamilton’s principle and port-Hamiltonian view

The Lagrangian itself is the time integral of L = T − V — the action functional. Hamilton’s principle, δ∫L dt = 0, gives the same equations of motion from a variational standpoint, and underpins the discrete-mechanics integrators (variational integrators, Marsden & West 2001) that preserve symplectic structure exactly in simulation. The Hamiltonian H = q̇^⊤ p − L = T + V (the total energy) leads to port-Hamiltonian formulations used in passivity-based control and energy-shaping designs (van der Schaft, Ortega).

Featherstone’s spatial vector algebra

Featherstone (1987, expanded 2008) unified twists, wrenches, transforms, and inertias into 6-D spatial vectors, yielding three canonical algorithms that every modern dynamics library implements:

• RNEA — Recursive Newton-Euler — inverse dynamics in O(n).
• CRBA — Composite Rigid-Body Algorithm — joint-space mass matrix M(q) in O(n²).
• ABA — Articulated-Body Algorithm — forward dynamics q̈ from τ in O(n), without forming M⁻¹.

ABA is the reason real-time forward simulation of a 30-DOF humanoid is feasible on commodity hardware; it avoids the O(n³) cost of forming and inverting M. ABA works by propagating an “articulated-body inertia” 𝓘_i^A and bias force p_i^A inward (tip → base), each accounting for the inertia of the subtree distal to link i with all distal joints free to move. The forward pass then propagates accelerations outward (base → tip), now using each joint’s articulated inertia to solve a 1 × 1 (for 1-DOF joints) inertial inverse — never the n × n joint-space inverse.

The Carpentier-Mansard analytical derivatives of RNEA, CRBA, and ABA (Pinocchio 2.0+) make the Jacobians ∂τ/∂q, ∂τ/∂q̇, ∂q̈/∂q, ∂q̈/∂q̇ available in closed form at O(n²) cost — the enabler of efficient DDP / iLQR trajectory optimisation, sensitivity analysis, and second-order optimisation in policy learning.

Passivity and structural properties

For any choice of C(q, q̇) satisfying the Christoffel construction:

• M(q) is symmetric positive-definite for all q.
• Ṁ(q) − 2C(q, q̇) is skew-symmetric: x^⊤ (Ṁ − 2C) x = 0 for all x.
• Total mechanical energy Ḣ = q̇^⊤ τ — the system is passive with respect to the (q̇, τ) port.

Passivity is the linchpin of most stability proofs in robot control — energy-based Lyapunov arguments, impedance control, port-Hamiltonian frameworks all lean on it. If your code’s M is non-symmetric or Ṁ − 2C fails the skew-symmetry test, your dynamics implementation is wrong.

Floating-base extension

For a floating-base robot (humanoid, quadruped, drone), append a 6-DOF unactuated joint between the world frame and the root link. The configuration is then q = (p_root, R_root, q_joints) with nq = 7 + n_joints (7 because the orientation is a unit quaternion) and nv = 6 + n_joints (velocity uses 3 angular + 3 linear components, not a quaternion derivative). The right-hand side of the manipulator equation has zero in the first 6 rows because the base is unactuated; all motion of the base must come through reaction forces (gravity + contact). Configuration integration is not simple Euler: q ← q ⊕ Δt v with ⊕ being the Lie-group exponential on SE(3). Skipping this step is the most common floating-base bug.

3. Practical math

Computing M, C, g from a dynamics library

In practice, almost no one derives M, C, g symbolically beyond n ≤ 3. Standard recipes for an arbitrary n-DOF URDF:

• Gravity vector g(q): one RNEA call at (q, q̇ = 0, q̈ = 0) with gravity enabled.
• Mass matrix M(q): either CRBA (O(n²), preferred) or n RNEA calls at (q, q̇ = 0, q̈ = e_j) with gravity disabled — column j is M·e_j.
• Coriolis vector c(q, q̇) = C(q, q̇)·q̇: one RNEA call at (q, q̇, q̈ = 0) with gravity disabled.
• Coriolis matrix C(q, q̇): rarely needed explicitly; if you do need it, Christoffel construction from M’s partial derivatives.

In Pinocchio:

import pinocchio as pin
model = pin.buildModelFromUrdf("ur5.urdf")
data  = model.createData()
M = pin.crba(model, data, q)              # O(n^2)
g = pin.computeGeneralizedGravity(model, data, q)
c = pin.nonLinearEffects(model, data, q, v) - g  # C v
tau = pin.rnea(model, data, q, v, a)      # full inverse dynamics
a   = pin.aba(model, data, q, v, tau)     # forward dynamics, O(n)

Inverse dynamics (the control workhorse)

Given a desired trajectory (q_d, q̇_d, q̈_d), compute the feedforward torque:

${\tau }_{\text{ff}}=M(q_d){\ddot{q}}_d+C(q_d,{\dot{q}}_d){\dot{q}}_d+g(q_d)$

This is the computed-torque (a.k.a. inverse-dynamics) control law. Add an outer PD on the joint error and the closed-loop becomes a linear, decoupled double integrator per joint — the foundational result enabling high-bandwidth manipulator control. Real-time cost on a 6-DOF arm: ~50 µs single-threaded on a modest ARM Cortex-A53; ~5 µs on a desktop x86-64.

The full computed-torque control law adds an outer PD loop with feedback linearisation:

$\tau =M(q)[{\ddot{q}}_d+K_p(q_d-q)+K_d({\dot{q}}_d-\dot{q})]+C(q,\dot{q})\dot{q}+g(q)$

Closed-loop error dynamics: ë + K_d ė + K_p e = 0, decoupled, per joint. Choose K_p = ω_n², K_d = 2ζω_n with ω_n ≈ 5–20 rad/s and ζ = 1 for critically damped tracking. The catch: full computed-torque requires accurate M, C, g — model error degrades performance gracefully but does not cause instability, while sensor noise on q̇ entering through K_d limits achievable bandwidth.

Forward dynamics (the simulator workhorse)

Given current (q, q̇) and applied τ, compute q̈ = M(q)⁻¹ · [τ − C(q, q̇)·q̇ − g(q)]. Naïvely O(n³) from the matrix inverse; ABA gives the same answer in O(n) without forming M⁻¹. Integrate q̈ → q̇ → q with a stable scheme (semi-implicit Euler for moderate stiffness, Runge-Kutta 4 for accuracy, MuJoCo’s implicit-in-velocity scheme for contact).

Operational-space (Cartesian) inertia

The end-effector behaves like a rigid body of effective inertia

$\Lambda (q)=(J(q)M(q{)}^{-1}{J}^{\top }(q){)}^{-1}$

at the operational point. This 6 × 6 matrix appears in operational-space control (Khatib 1987), in impedance shaping, and in collision-energy bounds. Its condition number blows up near singularities — the same singularities that plague kinematic inverses also plague Λ inversion.

The operational-space dynamics, projected onto end-effector coordinates x = f(q), is

$\Lambda (q)\ddot{x}+\mu (q,\dot{q})\dot{x}+p(q)=F_{\text{ee}}$

with μ = Λ J M⁻¹ C J^⊤ Λ − Λ J̇ J⁺ and p = Λ J M⁻¹ g. The mapping back to joint torques is τ = J^⊤ F_ee + (I − J^⊤ J^{†⊤})·τ_null, where the second term is a null-space torque that does not affect the end-effector — used for posture control and redundancy resolution on 7-DOF arms (Franka Panda, KUKA LBR iiwa).

Centroidal momentum (floating-base shortcut)

For locomotion, the centroidal momentum matrix A_G(q) maps joint velocity to total spatial momentum at the centre of mass: h_G = A_G(q) q̇. Its time derivative ḣ_G equals the sum of external wrenches at the COM (gravity + ground reactions). This 6-dimensional view collapses an 18- to 36-dimensional whole-body problem to a 6-dim “wrench balance” problem and is the basis of centroidal-dynamics MPC (Orin & Goswami 2008, Wensing & Orin 2016). All modern quadruped/biped control stacks use this representation at the planning layer.

Worked example A — 2-link planar arm (Lagrangian)

A canonical textbook problem. Two revolute joints, link lengths ℓ₁, ℓ₂, point masses m₁, m₂ at the end of each link, q = (q₁, q₂), gravity g = 9.81 m/s² acting in −y.

Kinetic energy (after computing the centre-of-mass velocities of each mass):

$T=\frac{1}{2}(m_1+m_2){\ell }_1^2{\dot{q}}_1^2+\frac{1}{2}{m}_2{\ell }_2^2({\dot{q}}_1+{\dot{q}}_2{)}^2+m_2{\ell }_1{\ell }_2\cos q_2{\dot{q}}_1({\dot{q}}_1+{\dot{q}}_2)$

Potential energy:

$V=(m_1+m_2)g{\ell }_1\sin q_1+m_{2}g{\ell }_2\sin (q_1+q_2)$

Apply Euler-Lagrange. The mass matrix:

$M(q)=\left[\begin{array}{cc}(m_1+m_2){\ell }_1^2+m_2{\ell }_2^2+2m_2{\ell }_1{\ell }_2\cos q_2& m_2{\ell }_2^2+m_2{\ell }_1{\ell }_2\cos q_2\\ m_2{\ell }_2^2+m_2{\ell }_1{\ell }_2\cos q_2& m_2{\ell }_2^2\end{array}\right]$

Coriolis matrix (let h = −m₂ ℓ₁ ℓ₂ sin q₂):

$C(q,\dot{q})=\left[\begin{array}{cc}h{\dot{q}}_2& h({\dot{q}}_1+{\dot{q}}_2)\\ -h{\dot{q}}_1& 0\end{array}\right]$

Gravity vector:

$g(q)=\left[\begin{array}{c}(m_1+m_2)g{\ell }_1\cos q_1+m_{2}g{\ell }_2\cos (q_1+q_2)\\ m_{2}g{\ell }_2\cos (q_1+q_2)\end{array}\right]$

Numerical evaluation at ℓ₁ = ℓ₂ = 0.30 m, m₁ = 2.0 kg, m₂ = 1.5 kg, q = (0, π/2), q̇ = (1.0, −0.5) rad/s:

• M = [[0.450, 0.135], [0.135, 0.135]] kg·m²
• C·q̇ = [(−m₂ℓ₁ℓ₂ sin q₂)·(2 q̇₁ q̇₂ + q̇₂²), (m₂ℓ₁ℓ₂ sin q₂)·q̇₁²] = [−(0.135)·(−0.75), (0.135)·1.0] = [0.101, 0.135] N·m
• g = [(3.5)(9.81)(0.30)·cos 0 + (1.5)(9.81)(0.30)·cos(π/2), (1.5)(9.81)(0.30)·cos(π/2)] = [10.30, 0] N·m

For q̈ = 0 (static hold), τ = C·q̇ + g = [10.40, 0.135] N·m. Verify M is symmetric ✓ and positive-definite (det = 0.0425 > 0, trace > 0) ✓.

Worked example B — UR5 trajectory inverse dynamics

A UR5 (6-DOF, ~18 kg arm mass, 5 kg payload, total ~23 kg lumped) executes a 1-s quintic from q₀ = (0, −π/2, π/2, 0, π/2, 0) to q_f = q₀ + (0.5, 0.5, 0.5, 0.5, 0.5, 0.5) rad. Peak |q̇| ≈ 0.94 rad/s (at t = 0.5 s), peak |q̈| ≈ 5.8 rad/s². Use Pinocchio with the official UR5 URDF (ur_description/urdf/ur5.urdf.xacro).

Order-of-magnitude torques (per joint, peak):

Joint	g(q) hold	M·q̈ peak	C·q̇ peak	τ peak
shoulder pan (1)	0.5 N·m	22 N·m	4 N·m	~26 N·m
shoulder lift (2)	95 N·m	35 N·m	6 N·m	~136 N·m
elbow (3)	32 N·m	14 N·m	3 N·m	~49 N·m
wrist 1 (4)	1.5 N·m	1.0 N·m	0.2 N·m	~2.7 N·m
wrist 2 (5)	0.4 N·m	0.4 N·m	0.1 N·m	~0.9 N·m
wrist 3 (6)	0.1 N·m	0.05 N·m	0.02 N·m	~0.2 N·m

Shoulder lift (J2) is gravity-dominated by an order of magnitude — typical of any shoulder joint on a horizontally extended arm. This is why gravity compensation alone captures 70–90% of the steady-state torque demand for most cobot use cases.

Worked example B’ — Dynamic parameter identification for a 6-DOF arm

Given a recorded trajectory {q_k, q̇_k, q̈_k, τ_k} for k = 1, …, N, the inverse-dynamics equation is linear in the inertial parameters (Khalil & Dombre 2002, Atkeson et al. 1986):

${\tau }_k=Y(q_k,{\dot{q}}_k,{\ddot{q}}_k)\pi$

where π ∈ ℝ^{10n} stacks all (m, mc_x, mc_y, mc_z, I_xx, I_yy, I_zz, I_xy, I_xz, I_yz) per link, and Y is the regressor matrix computed by Pinocchio’s computeJointTorqueRegressor. Stacking N samples:

$\underset{T\in {\mathbb{R}}^{nN}}{\underbrace{\left[\begin{array}{c}{\tau }_1\\ ⋮\\ {\tau }_N\end{array}\right]}}=\underset{\mathcal{Y}}{\underbrace{\left[\begin{array}{c}{Y}_1\\ ⋮\\ Y_N\end{array}\right]}}\pi$

Least squares: π̂ = (𝒴^⊤𝒴)⁻¹ 𝒴^⊤ T. Caveats: (1) only base parameters (typically 50–60 for a 6-DOF arm, not the full 60) are identifiable — the rest sit in the null space of 𝒴 and need to be regularised or removed via QR-based base-parameter projection; (2) the trajectory must be exciting (rich in q, q̇, q̈) — typical practice is to optimise a periodic Fourier-series trajectory to minimise cond(𝒴) (Swevers et al. 1997); (3) physical constraints (m > 0, I ≻ 0, COM inside link’s convex hull) need an SDP / NLP formulation to enforce strictly (Sousa & Cortesão 2014).

Practical result for a UR10e identified from a 60-s excitation trajectory at 500 Hz: torque residuals drop from ~15 % of nominal (factory URDF) to ~2 % (identified model). The remaining 2 % is friction, transmission flexibility, and unmodelled cable forces.

Worked example C — Floating-base quadruped (Spot-class)

A 12-leg-DOF quadruped (3 actuated joints × 4 legs) on a 6-DOF floating base has total dimension n = 6 + 12 = 18. Generalized coordinates: q = (p_base, R_base, q_legs). The dynamics:

$\left[\begin{array}{cc}{M}_{bb}& M_{bj}\\ M_{bj}^{\top }& M_{jj}\end{array}\right]\left[\begin{array}{c}{\ddot{q}}_b\\ {\ddot{q}}_j\end{array}\right]+\left[\begin{array}{c}{h}_b\\ h_j\end{array}\right]=\left[\begin{array}{c}0\\ {\tau }_j\end{array}\right]+{\sum }_c\left[\begin{array}{c}{J}_{c,b}^{\top }\\ J_{c,j}^{\top }\end{array}\right]{\lambda }_c$

with h = C·q̇ + g, and λ_c the contact wrench at each foot in contact (typically 3 components per point contact, 6 per surface contact). The first six rows have zero generalised force on the right-hand side because the base is unactuated; this is the unactuated dynamics constraint that forces all locomotion to happen through ground reaction.

A typical 50 Hz MPC over a 0.5-s horizon solves a QP with q̈, λ as decision variables, friction-cone constraints λ_z ≥ 0 and ‖λ_{xy}‖ ≤ μ λ_z, plus dynamics. Solver cost on a Spot-class robot (Pinocchio + OSQP or HPIPM): 5–15 ms per QP on an Intel NUC i7.

Concrete numerical example — a 12 kg quadruped trotting at 1.0 m/s on flat ground with two diagonal feet in stance:

• Base inertia (estimated): I_base ≈ diag(0.07, 0.30, 0.32) kg·m².
• Per-leg inertia at hip ≈ 0.02 kg·m² (mainly thigh + calf).
• Gravity load per stance foot: 12 · 9.81 / 2 ≈ 59 N vertical.
• Friction coefficient μ ≈ 0.6 (rubber-foam on indoor flooring) → max horizontal force ≈ 35 N per foot.
• Peak base linear acceleration during trot: ~3 m/s² → total horizontal force ≈ 36 N — right at the friction limit, which is why trot transitions can slip on slick floors.
• MPC horizon: 10 nodes × 50 ms = 0.5 s; ~180 decision variables, ~150 constraints; QP solve ≤ 5 ms on a Jetson Orin NX.

4. Design heuristics

Joint-torque sizing

Rule of thumb for arm joint sizing:

${\tau }_{\text{rated}}\ge 1.3\times ({\tau }_{\text{gravity,max}}+{\tau }_{\text{dynamic,max}})+{\tau }_{\text{friction}}$

The 1.3 multiplier covers payload variance, model error, and worst-case configuration. For a horizontal-extension cobot (UR5/Panda class):

• τ_gravity (shoulder lift, fully extended, max payload) is usually 70–90% of nameplate τ_rated.
• τ_dynamic at typical cobot speeds (q̇ ≤ 1 rad/s, q̈ ≤ 5 rad/s²) is 10–30% of τ_gravity.
• Industrial high-speed pick-and-place (Delta, SCARA): τ_dynamic can reach 150% of τ_gravity at peak.

Mass-matrix conditioning

cond(M) > 1000 ⇒ numerical conditioning problems in operational-space inverse, in computed-torque feedback, and in QP-based controllers. Common causes:

• Two co-linear joint axes (kinematic singularity propagates to dynamic singularity).
• Very different link inertias on the same chain (e.g., a small wrist link with a heavy payload at the flange).
• Floating-base coordinates expressed in body frame near a gimbal-lock orientation.

Use Damped Least Squares (DLS) when forming Λ = (JM⁻¹J^⊤)⁻¹ near singularities:

${\Lambda }_{\text{DLS}}=(J{M}^{-1}{J}^{\top }+{\lambda }^{2}I{)}^{-1},\lambda \approx 0.01-0.1$

Gravity-compensation control

The simplest model-based controller:

$\tau =g(q)+K_p(q_d-q)+K_d({\dot{q}}_d-\dot{q})$

Removes the largest configuration-dependent torque term, dramatically improving low-speed tracking and zero-impedance “hand-guiding” feel. Standard on Universal Robots e-series, Franka Emika Panda, Kinova Gen3. Adding inertia and Coriolis terms (full computed-torque) gains another 5–10× tracking accuracy at the cost of model-fidelity sensitivity.

Forward-dynamics integration

• Explicit Euler — q̇ ← q̇ + Δt·q̈, q ← q + Δt·q̇. Unstable for stiff systems; usable for soft, slow toys.
• Semi-implicit (symplectic) Euler — q̇ ← q̇ + Δt·q̈, q ← q + Δt·q̇_new. Preserves energy structure; default in MuJoCo, Bullet, PhysX.
• RK4 — accurate for smooth dynamics, expensive (4 RNEA per step). Used in offline simulation and trajectory optimisation.
• Implicit — required for very stiff contact or actuator dynamics; cost is 1 nonlinear solve per step.

Typical Δt: 1 ms for cobot sim, 0.5–1 ms for legged-robot whole-body MPC, 2–4 ms for drone hover.

A useful CFL-like stability criterion for spring-damper joint approximations: Δt ≤ 2 / ω_n where ω_n = √(k/m) is the highest natural frequency in the system. For k = 10⁴ N·m/rad, I = 10⁻³ kg·m², ω_n ≈ 3160 rad/s, so Δt ≤ 0.6 ms. Halve that for safety with explicit integrators.

Inverse-dynamics feedforward without C term

The simplified law τ = M(q) q̈_d + g(q) (omitting Coriolis) is always stable when paired with a PD on (q, q̇) — the omitted C·q̇ term is bounded for bounded q̇ and acts as a disturbance the PD can reject. The full C term contributes most at high q̇; for slow cobots (q̇ ≤ 0.5 rad/s) omitting it costs < 5 % tracking error. Many production controllers ship without it for code-simplicity reasons.

Choosing K_p, K_d in computed-torque control

After feedback linearisation the per-joint error dynamics are ë + K_d ė + K_p e = 0. Pick:

• desired closed-loop natural frequency ω_cl (typically 5–20 rad/s for cobot, 50–200 rad/s for direct-drive legged)
• damping ratio ζ (1.0 critically damped; 0.7 for slight overshoot but faster rise)
• then K_p = ω_cl², K_d = 2 ζ ω_cl

Bandwidth ceiling: ω_cl ≤ ω_joint_resonance / 5. Exceeding it excites unmodelled flexible-joint modes; the symptom is high-frequency chatter that PD cannot remove (because PD assumed rigid dynamics — see vibration-dynamics).

Linearisation around an equilibrium

For control design, linearise the manipulator equation around an equilibrium configuration q_eq with q̇_eq = 0 and τ_eq = g(q_eq):

$M(q_{\text{eq}})\Delta \ddot{q}+\frac{\partial g}{\partial q}{|}_{q_{\text{eq}}}\Delta q=\Delta \tau$

This is a multivariable second-order linear system. Diagonalising via the generalised eigenproblem (∂g/∂q)·v = ω² M·v yields dynamic mode shapes and natural frequencies at the equilibrium. For a horizontal-extended robot, the lowest dynamic eigenfrequency typically falls between 1–5 Hz — the mode the shoulder gravity load “wants” to oscillate at when perturbed. Closed-loop bandwidth above this frequency is straightforward; below it, the open-loop system has inverted-pendulum-like dynamics that PID can stabilise but not without care.

State-space form for LQR / Kalman design:

$\frac{d}{dt}\left[\begin{array}{c}\Delta q\\ \Delta \dot{q}\end{array}\right]=\left[\begin{array}{cc}0& I\\ -M^{-1}\frac{\partial g}{\partial q}& -M^{-1}{C}_{\text{eq}}\end{array}\right]\left[\begin{array}{c}\Delta q\\ \Delta \dot{q}\end{array}\right]+\left[\begin{array}{c}0\\ M^{-1}\end{array}\right]\Delta \tau$

See state-space-lqr for the controller-side treatment.

When to skip dynamics altogether

For a tracking task with low q̇ and low q̈ relative to the joint’s torque envelope, a position-controlled robot driven by an outer kinematic controller (Jacobian-pseudoinverse Cartesian control) is dynamics-free at the user level — the manufacturer’s joint servo already handles M, C, g internally. This describes 90 % of real industrial deployments. Reach for explicit dynamics only when:

• forces / contacts are an explicit task variable (assembly, polishing, surgery)
• accelerations approach the joint’s torque envelope (high-speed pick-and-place)
• the base is unactuated (mobile manipulation, legged, aerial)
• the controller is model-based by design (impedance, MPC, OSC)

Contact and constraints

Add J_c^⊤ λ to the right-hand side. The complementarity condition λ_n · gap_n = 0 with λ_n ≥ 0 and gap_n ≥ 0 makes contact a Linear Complementarity Problem (LCP) — solved by Lemke (Bullet), projected Gauss-Seidel (MuJoCo, ODE), or a relaxed convex QP (Drake’s hydroelastic contact, MuJoCo’s smooth contact). Smooth contact trades physical exactness for differentiability — essential for any gradient-based controller or RL pipeline.

Real-time budget

Dynamics-library cost on a Cortex-A53 @ 1.4 GHz (Raspberry Pi 4 ballpark, Pinocchio benchmarks 2019):

Robot	n	RNEA	CRBA	ABA
UR5	6	12 µs	18 µs	20 µs
Franka Panda	7	15 µs	24 µs	26 µs
Atlas humanoid	36	110 µs	320 µs	180 µs
Spot quadruped	18	60 µs	130 µs	95 µs

On desktop x86-64 (Intel i7-10700, single thread), divide by ~10. GPU-parallel sims (Isaac Lab, Genesis) achieve 10⁵–10⁶ env-steps/s by batching thousands of independent rollouts.

5. Components & sourcing (software ecosystem)

The “components” of a robot-dynamics pipeline are libraries, not hardware. As of 2026-05:

Library	Vendor / License	Strength	Where it fits
Pinocchio	INRIA, BSD-2	Fastest spatial-algebra C++/Python; analytic derivatives	Canonical for control & trajectory optimisation
RBDL	M. Felis (Heidelberg), LGPL	Featherstone reference impl; clearest source	Academic, education
Drake	TRI / MIT, BSD-3	Symbolic + numerical, multibody plant, contact	Verification, simulation-driven design
MuJoCo	DeepMind (was Roboti LLC), Apache-2	Fastest contact-rich sim; differentiable	RL training, manipulation research
Bullet Physics	Open-source, Zlib	Game-physics roots; broad use	ROS Gazebo classic, hobby robotics
PhysX	NVIDIA, BSD-3 (5.x)	GPU-accelerated; production game engine	Isaac Lab backbone
Isaac Lab / Isaac Sim	NVIDIA, NVIDIA OSL	Massively parallel GPU RL	Sim-to-real RL pipelines
Genesis	CMU/SAIL, Apache-2 (2024)	Pure-GPU sim; benchmark-leading throughput	Frontier RL, world models
iDynTree	IIT (iCub team), BSD	Humanoid-focused; URDF-to-DH	Whole-body MPC, iCub
Crocoddyl	LAAS-CNRS, BSD	DDP/iLQR trajectory optimisation on Pinocchio	Legged robot OCS
MATLAB Robotics System Toolbox	MathWorks, commercial	GUI-based, Simulink integration	Industry prototyping
MATLAB Simscape Multibody	MathWorks, commercial	CAD-import multi-body	Mechanism design

Recommended stack for a new project as of 2026:

• Modelling: URDF / MJCF (see robotics-control).
• Dynamics + control: Pinocchio + Crocoddyl.
• Simulation: MuJoCo (for contact-rich); Isaac Lab (for parallel RL).
• Visualisation: meshcat (Pinocchio), MuJoCo viewer.

The MuJoCo Menagerie (https://github.com/google-deepmind/mujoco_menagerie) is the gold-standard open repository of validated MJCF models for 30+ commercial robots (UR5e, UR10e, Franka Panda, Kinova Gen3, Boston Dynamics Spot, Unitree Go2 / H1, Allegro Hand, Shadow Hand, ANYmal, …). Start there before authoring a URDF from CAD.

6. Reference data

Algorithm complexity (Featherstone 2008)

Algorithm	Computes	Complexity	Typical use
RNEA	τ = ID(q, q̇, q̈)	O(n)	Inverse dynamics, computed torque
CRBA	M(q)	O(n²)	Mass matrix for OSC, MPC
ABA	q̈ = FD(q, q̇, τ)	O(n)	Simulation
CRBI	Composite inertia at each link	O(n)	Centroidal momentum
Newton-Euler symbolic	Closed-form expression	O(?)	Hand derivation, n ≤ 3

Joint-space dimension (typical platforms)

Platform	Type	DOF	Floating base?	Total q dim
Universal Robots UR5e	Fixed serial 6R	6	no	6
Franka Emika Panda	Fixed serial 7R	7	no	7
Kinova Gen3	Fixed serial 7R	7	no	7
ABB IRB 1200	Fixed serial 6R	6	no	6
Stäubli TX-90	Fixed serial 6R	6	no	6
ANYbotics ANYmal C	Quadruped	12	yes	6 + 12 = 18
Boston Dynamics Spot	Quadruped	12	yes	18
Unitree Go2	Quadruped	12	yes	18
Unitree H1	Bipedal humanoid	19	yes	25
Agility Cassie	Bipedal	20 (incl. springs)	yes	26
Boston Dynamics Atlas (electric)	Humanoid	28	yes	34
Shadow Dexterous Hand	Hand only	24	no (attached)	24
Quadrotor (DJI F450)	Drone	4 (rotor speeds)	yes (6)	10

UR5 dynamic parameters (Kufieta 2014 ID, used in MuJoCo Menagerie)

Link	Mass (kg)	COM (m, link-frame)	Inertia diag (kg·m²)
Base	4.0	(0, 0, 0)	(0.0044, 0.0044, 0.0072)
Shoulder	3.7	(0, 0, 0.028)	(0.0103, 0.0103, 0.0066)
Upper arm	8.4	(0.21, 0, 0.14)	(0.226, 0.226, 0.015)
Forearm	2.3	(0.20, 0, 0.013)	(0.025, 0.025, 0.004)
Wrist 1	1.2	(0, 0.005, 0)	(0.002, 0.002, 0.002)
Wrist 2	1.2	(0, 0.005, 0)	(0.002, 0.002, 0.002)
Wrist 3	0.19	(0, 0, −0.02)	(1.3e-4, 1.3e-4, 8.6e-5)

Total moving mass ≈ 17 kg; rated payload 5 kg.

Common dynamics-library throughput (Pinocchio 2.6, 2023 benchmarks)

Single-threaded, Intel i7-10700K @ 3.8 GHz, double precision, all timings in microseconds per call.

Robot	n_v	RNEA	CRBA	ABA	dRNEA/dq
2-DOF planar	2	0.6	0.7	1.0	1.2
Kinova Gen3	7	1.9	2.7	2.9	4.5
Franka Panda	7	1.9	2.7	3.0	4.5
KUKA LBR iiwa 14	7	1.9	2.7	3.0	4.5
UR10e	6	1.7	2.4	2.6	4.1
Atlas (Boston Dyn)	36	17	49	28	45
ANYmal-C	18	8.0	19	14	22
HRP-4 humanoid	36	17	50	28	45
Talos humanoid	38	18	53	30	48

ARM Cortex-A72 (Raspberry Pi 4) numbers are roughly 5–8× higher. M1 Apple Silicon comparable to i7-10700K.

Comparison — manipulator dynamic parameters

Robot	Total mass	Payload	Reach	Peak τ (J2 shoulder)	Peak speed (TCP)
UR5e	20.6 kg	5 kg	0.85 m	150 N·m	1.0 m/s
UR10e	33.5 kg	12.5 kg	1.30 m	330 N·m	1.0 m/s
Franka Panda	18 kg	3 kg	0.85 m	87 N·m	2.0 m/s
Kinova Gen3	8.2 kg	4 kg	0.90 m	39 N·m	0.5 m/s
KUKA LBR iiwa 7 R800	23.9 kg	7 kg	0.80 m	176 N·m	2.0 m/s
ABB IRB 1200	54 kg	5 kg	0.90 m	~600 N·m	4.0 m/s

Manipulator Coriolis vs gravity dominance

For an n-DOF arm, a useful diagnostic is the ratio of peak Coriolis torque to peak gravity torque along a representative trajectory:

Class	Typical q̇_max	‖C·q̇‖ / ‖g‖	Implication
Slow cobot (UR, Franka teach mode)	0.5 rad/s	< 5 %	Gravity-comp PD adequate
Cobot at rated speed	1.5 rad/s	10–25 %	Add C term to feedforward
High-speed industrial 6R	4 rad/s	50–100 %	Full computed-torque required
Direct-drive research arm	8 rad/s	100–300 %	Coriolis dominates; OSC mandatory
Quadruped leg swing	15 rad/s	200–500 %	Whole-body QP essential

The ratio scales as q̇²/q̇ = q̇ to leading order (because C ∝ q̇, g constant), so doubling speed roughly doubles the ratio.

Reflected motor inertia (typical cobot joints)

Drivetrain	Gear ratio N	Motor inertia I_m (kg·m²)	Reflected I (kg·m²)
Harmonic drive, large cobot	100	5·10⁻⁵	0.5
Harmonic drive, small cobot	50	2·10⁻⁵	0.05
Cycloidal, industrial 6-axis	80	2·10⁻⁴	1.28
Planetary, legged QDD	6	5·10⁻⁵	1.8·10⁻³
Direct drive (none)	1	1·10⁻³	1·10⁻³
Belt + spur, hobby servo	4	1·10⁻⁵	1.6·10⁻⁴

Real-time inverse-dynamics budgets (production deployments)

Stack	Robot	Loop rate	ID share of budget
UR e-Series internal	UR5e / UR10e	500 Hz	~25 µs (~1.3 %)
Franka Control Interface	Panda	1000 Hz	~30 µs (~3.0 %)
KUKA Sunrise.OS	LBR iiwa 14	1000 Hz	~50 µs (~5.0 %)
Boston Dynamics Spot SDK	Spot quadruped	333 Hz	~150 µs (~5.0 %)
ANYbotics ANYmal control PC	ANYmal C	400 Hz	~120 µs (~4.8 %)
MIT Mini Cheetah open-source	Mini Cheetah	1000 Hz	~80 µs (~8 %)
PX4 attitude inner loop	Quadrotor	4000 Hz	~5 µs (~2 %)

ID is rarely the bottleneck — the QP / OCS solver, the contact pre-processing, and the estimator typically dominate.

Typical contact-model parameters

Model	Stiffness	Damping	Penetration tolerance
Penalty (hard)	10⁶–10⁸ N/m	10²–10³ N·s/m	< 0.1 mm
Penalty (soft)	10⁴–10⁵ N/m	10–10² N·s/m	1–5 mm
Hydroelastic (Drake)	E_modulus ≈ 10⁶ Pa	viscous dissipation	continuous
MuJoCo solref/solimp	time-constant Δt_ref ≈ 0.02 s	softness band ≈ 1 mm	adaptive

Penalty methods are easy but inject stiffness; constraint-based (MuJoCo, Drake) solves a complementarity problem each step but tolerates Δt up to 5 ms with stiff contact.

Recommended Δt by simulation target

Target	Δt (typical)	Integrator
Off-line trajectory verification	1–4 ms	RK4
Real-time MPC rollout	5–20 ms	Semi-implicit Euler
Contact-rich RL training	2–5 ms	MuJoCo implicitfast
High-bandwidth control inner loop	0.5–1 ms	Semi-implicit Euler
Aerial vehicle attitude	2–5 ms	RK4 (low-stiffness)
Legged whole-body QP	0.5–1 ms	Constraint-based contact

Typical joint-space gains (computed-torque)

Robot class	K_p (rad⁻²)	K_d (rad⁻¹·s)	ω_cl (Hz)
Slow industrial 6R (ABB IRB, KUKA KR)	400	40	~3
Cobot (UR, Franka, Kinova)	100	20	~1.5
Direct-drive 7R research arm	2500	100	~8
Quadruped leg joint (QDD)	1600	50	~6
High-bandwidth manipulator wrist	10000	200	~16
Drone attitude inner loop	90000	600	~50

Numbers are starting points; final tuning depends on identified flexibility, sensor noise, and Δt. As a rule, never set bandwidth above 1/5 of the first observed structural resonance.

Floating-base coordinate counts

Robot	nq (config)	nv (velocity)	Notes
Fixed 6R arm	6	6	Standard
Fixed 7R arm (Panda, iiwa)	7	7	Redundant
Quadrotor	7	6	quaternion + 3 pos
Quadruped 12-joint	19	18	7 base + 12 joints
Bipedal 24-joint	31	30	7 base + 24 joints
Humanoid 36-joint	43	42	7 base + 36 joints
Dexterous hand 24 + arm 7	31	31	Fixed base

Stiffness / damping ratios for joint approximations

A geared joint can be modelled as a rotational spring-damper for the purposes of frequency-domain analysis. Typical values for a harmonic-drive cobot joint:

• Stiffness k_j: 5,000–20,000 N·m / rad (UR-class), 30,000–100,000 N·m / rad (industrial robot)
• Damping c_j: 5–50 N·m·s / rad
• First structural resonance: 10–40 Hz (cobot), 30–80 Hz (industrial 6-axis)
• Closed-loop joint-bandwidth ceiling: ~ω_n / 3 (Brogliato 2007 rule of thumb)

See vibration-dynamics for the underlying SDOF / MDOF treatment.

Drivetrain inertia reflection

A motor of rotor inertia I_m driving the load through a gear ratio N (output speed = motor speed / N) reflects an inertia

$I_{\text{reflected}}=N^2{I}_m$

to the joint output. For a typical cobot joint with harmonic drive N = 100 and I_m = 5·10⁻⁵ kg·m², the reflected motor inertia is 0.5 kg·m² — frequently larger than the link inertia itself. Practical consequence: at high gear ratio the dynamics are dominated by motor inertia, the off-diagonal couplings of M shrink, and a single-joint PID becomes adequate (because the system has decoupled itself mechanically). Direct-drive robots (MIT Mini Cheetah, Unitree A1/Go2 with quasi-direct-drive QDD actuators at N ≈ 6) intentionally avoid this regime to retain rich multi-joint coupling for compliant control.

Backdrivability and torque transparency

A drivetrain is backdrivable if an external torque applied at the output produces visible motion at the input. Backdrivability scales roughly as 1 / (N² · η), where η is the gear efficiency. Harmonic drives (N ≈ 100, η ≈ 0.7) are barely backdrivable; quasi-direct-drive planetary (N ≈ 6, η ≈ 0.9) are highly backdrivable. Torque-transparency is the prerequisite for proprioceptive force estimation (no torque sensor needed) — the basis of the MIT Cheetah / ANYmal / Unitree class of robots.

Common dynamic-parameter pitfalls (CAD vs reality)

Source of error	Typical magnitude
Motor + harness mass missing from CAD	10–30 % link mass underestimate
Default uniform-density inertia from CAD	20–50 % I_zz error
Inertia expressed in link-origin frame instead of COM	error scales with m·d²
Tool / end-effector mass missing	100 % at flange link
Cable-loom routing mass	5–15 % distributed
Counterweight or spring missing	model bias, configuration-dependent

Calibrate against real τ from slow-motion identification trajectories — never trust CAD-derived dynamic parameters for anything safety-critical.

7. Failure modes & debugging

Wrong inertia data from CAD

Mass and inertia tags in URDF often start from CAD bulk-density estimates. Real robots have motors, harnesses, brackets — densities are non-uniform. Symptom: feed-forward torques are off by 20–50% at high acceleration; tracking error grows with speed. Fix: dynamic identification — execute exciting trajectories, record (q, q̇, q̈, τ), regress for the 10n Newton-Euler base parameters (or the reduced minimal set; Khalil/Dombre 2002). Pinocchio includes pin.computeJointTorqueRegressor; the OpenSoT / Crocoddyl ecosystem has end-to-end tooling.

Unmodelled friction

The single most common cause of poor low-speed tracking. Add to the model:

${\tau }_{\text{friction}}=F_v\dot{q}+F_c\mathrm{sgn}(\dot{q})+{\tau }_{\text{stribeck}}(\dot{q})$

Typical magnitudes for a harmonic-drive cobot joint: F_v = 0.1–1.0 N·m·s/rad, F_c = 0.5–5 N·m. Coulomb friction alone halves your low-speed tracking error if your robot is geared.

Stick-slip oscillation

Geared joints with high static friction can stick during slow motion and then break free abruptly. Visible as low-frequency (1–5 Hz) sawtooth in position error. Cures (in order): smoother trajectory generation, dither (small high-frequency torque), reduced gear ratio, direct-drive.

Numerical drift in long sims

Floating-base orientations represented as quaternions drift off the unit-norm manifold. Renormalise every step: q ← q / ‖q‖. Rotation matrices drift off SO(3) — re-orthogonalise via SVD or Gram-Schmidt every ~100 steps. Time integration error compounds; never trust a position integrated past ~10⁴ steps without a measurement update.

Stiff contact in simulation

Penalty-based contact (k_contact ≈ 10⁵–10⁷ N/m) imposes the stiffest dynamics in the system. Symptom: explicit integrators explode; required Δt drops to 10–100 µs. Cures:

• Use a constraint-based contact solver (MuJoCo, Drake) that doesn’t introduce spring stiffness.
• Use semi-implicit / implicit integration (default in MuJoCo with integrator="implicitfast").
• Lower contact stiffness with mild “softness” (mujoco solref, solimp) — trades realism for stable simulation.

Mass-matrix singularity

cond(M) blows up near kinematic singularities and at certain end-effector configurations. Operational-space control fails first (because of M⁻¹), then joint-space computed-torque (because Ṁ ≈ 2C breaks). Detect with numpy.linalg.cond(M); fall back to DLS or to a joint-space controller in singular zones.

Coriolis sign / passivity check

The skew-symmetry test x^⊤(Ṁ − 2C)x = 0 is a runtime sanity check for any dynamics implementation. Compute it on a few random (q, q̇, x) triples after any C derivation or library change. A non-zero result means C is inconsistent with M’s time derivative — most often a sign or index error in a Christoffel symbol.

Joint limits in simulation

Hard joint limits enforced by Lagrange multipliers (Drake style) preserve energy and complementarity but are expensive. Penalty-spring limits (MuJoCo default) are cheap but inject stiffness. Practical heuristic: use penalty limits with safety margin δ_lim ≈ 0.05 rad before the hard mechanical stop, soft cubic-spring profile.

Cartesian vs joint coordinate confusion

A frequent bug: applying joint-space M to a Cartesian command, or operational-space Λ to a joint command. Always tag your vectors with their frame. Cartesian dynamics needs Λ = (J M⁻¹ J^⊤)⁻¹, not M projected by J.

Wrong base type

A humanoid simulated with a fixed base falls trivially because the unactuated balance dynamics are missing. A fixed-base arm simulated with a floating base has 6 ill-defined unactuated DOF and integrates to garbage. Check model.nq, model.nv and confirm they match your design intent before running anything.

Energy-drift diagnostic

For a passive system simulated without input torque or external forces, total mechanical energy H = ½ q̇^⊤ M q̇ + V(q) should be conserved. Plot H(t) over a 10-second free-swing rollout. A monotone decrease indicates numerical damping from the integrator (often acceptable, ≤ 0.1 %/s). A monotone increase indicates the integrator is injecting energy — the simulation will eventually blow up; switch to a symplectic integrator (MuJoCo implicit, Drake implicit_euler) or shorten Δt. Erratic, non-monotone change is usually a contact-solver bug (penetration spikes, regularisation chatter).

Quaternion vs rotation-matrix representation

For floating-base orientation, quaternions (4 floats) are compact and singularity-free but need unit-norm renormalisation each integration step. Rotation matrices (9 floats) are easier to compose but drift off SO(3) and need periodic SVD-based re-orthogonalisation. Modified Rodrigues parameters (3 floats) are minimal but singular at 360°. Pinocchio, Drake, and MuJoCo all use quaternions internally; URDF stores Euler angles (RPY) but converts on load. Never store accumulated orientation as Euler angles in code — gimbal lock kills you near ±90° pitch.

Friction-cone violations

In contact-rich simulation, requested tangential contact force can exceed μ · F_normal — physically impossible. Penalty-based solvers silently clip or oscillate; constraint-based solvers (MuJoCo, Drake) treat it as a hard inequality and the QP becomes infeasible. Symptom: the planner asks for an acceleration the contact cone cannot provide, controller falls back to a default, robot “skids”. Diagnostic: log ‖F_tangent‖ / (μ F_normal) per contact; if it exceeds 1, your reference trajectory is dynamically infeasible.

Symptom-cause table

Symptom	Likely cause	Fix
Inverse-dynamics torque too small at high q̈	Inertias under-estimated in URDF	Dynamic identification
Inverse-dynamics torque too small at low q̇	Friction not modelled	Add F_v q̇ + F_c sgn(q̇)
Forward sim explodes after a few seconds	Explicit integrator + stiff contact	Switch to implicit / constrained solver
Quaternion drifts off unit norm	Missing renormalisation step	q ← q / ‖q‖ each step
Position drifts during a “hold”	Gravity vector computed in wrong frame	Re-check world-frame gravity sign
Robot drifts sideways on flat ground (legged)	Friction coefficient too low	Raise μ to 0.6+ or check contact normals
Tracking error grows with speed	Coriolis term sign error or missing	Verify Ṁ − 2C skew-symmetry
OSC controller blows up near a singularity	M⁻¹ ill-conditioned	DLS (λ ≈ 0.05)
Whole-body QP infeasible	Friction-cone or torque-limit conflict	Soften task weights, relax cone
Joint angle integrates past mechanical stop in sim	Penalty limits too soft	Stiffen or use hard-constraint solver
Floating-base robot accelerates up off the ground	Gravity sign wrong, or base mass too low	Recheck world gravity vector and root mass
Reaction-wheel sat overshoots commanded attitude	Gyroscopic ω × Iω term missing	Include in Euler eq for spinning rotors
Inverse-dynamics produces NaN	q̇ or q̈ contains NaN from upstream filter	Sanitise inputs; check estimator
Identified inertias non-physical (negative diag)	No SDP constraints in regression	Use SDP-constrained ID (Sousa & Cortesão)
ZMP / contact wrench outside support polygon	Reference trajectory infeasible	Replan with support-polygon constraint

7b. Implementation pitfalls (library-specific)

URDF / MJCF authoring

• Inertia frame: URDF puts the <inertia> tag inside an <inertial> block with its own <origin> — the inertia tensor must be expressed in that frame, not in the link frame. Most CAD exporters get this wrong by default; check by re-computing m, c, I from the mesh and comparing.
• Joint axis convention: URDF joint axes are expressed in the parent link frame; SDF uses the child link frame; MJCF uses the body frame of the joint’s body. Switching formats without re-aligning axes is an instant source of mirrored-joint bugs.
• Continuous vs revolute joints: A continuous URDF joint has no limits, modelled in dynamics with q ∈ ℝ (no wrap). A revolute joint has hard limits applied as penalty or constraint. Confusing the two breaks limit-aware planners.
• Mimic joints: URDF <mimic> ties two joints by a scalar; Pinocchio supports them, MuJoCo does not (use an equality constraint instead).

Pinocchio specifics

• pin.RNEA and pin.ABA require q of length model.nq and v, a of length model.nv. For a fixed-base 6R arm both equal 6; for a floating-base robot they differ (nq = nv + 1 per quaternion). Mixing them silently produces garbage.
• The frame placement returned by pin.framesForwardKinematics is in the world frame by default; the Jacobian from pin.computeFrameJacobian defaults to LOCAL (body) frame — call with pin.ReferenceFrame.WORLD to align.
• Analytic derivatives (pin.computeRNEADerivatives) are exact only if the URDF is twice-differentiable, which excludes Coulomb friction unless smoothed.

MuJoCo specifics

• The MJCF <inertial> tag’s diaginertia is in the body’s principal-axis frame; full off-diagonal inertia is set via fullinertia and is in the body frame. Conversion from URDF is via urdf2mjcf (mujoco CLI).
• Default integrator is Euler (semi-implicit); switch to implicit or implicitfast for high-stiffness contact.
• The actuator tag couples motor torque to a joint; MuJoCo’s force-vs-position controllers differ subtly from a textbook computed-torque law in their treatment of velocity damping (the damping parameter on the joint adds a passive term independent of the controller).

Drake specifics

• MultibodyPlant<double> is the canonical class; MultibodyPlant<AutoDiffXd> gives gradients via dual numbers. Switching between them is template-only at compile time.
• Drake’s contact model is hydroelastic by default — gives smooth, continuously differentiable contact forces ideal for trajectory optimisation but slightly stiffer to integrate than point-contact + LCP.
• Time stepping with discrete_period > 0 uses a TAMSI solver for contact; with discrete_period == 0 uses continuous-time integration with event detection.

8. Case studies

Universal Robots e-series gravity compensation

UR’s e-series cobots (CB3 + e-Series) use a gravity-compensated PD scheme as the user-facing “freedrive” / “hand-guidance” mode. Internal dynamics model is parameter-identified per arm at the factory and refined per joint via accelerometer and torque-sensor data; the per-arm dh.conf and calibration.conf files include link masses, COMs, and inertia tensors usable in Pinocchio with the official ur_description URDFs (https://github.com/UniversalRobots/Universal_Robots_ROS2_Description). Published gravity-comp tracking error is ≤ 5° at any TCP pose without external calibration; with offline pin.computeJointTorqueRegressor ID refinement, this drops to ~1°. The same model drives the 500 Hz teach-pendant force-mode torque estimator that gives UR its certified 250 N collision-force ISO 10218 compliance.

Boston Dynamics Spot / ANYmal / Unitree quadrupeds — whole-body MPC

The dominant control architecture for modern quadrupeds is Centroidal Momentum Dynamics + Whole-Body MPC:

1. Compute floating-base dynamics (ABA, n = 18) with stance-foot contact constraints.
2. Plan a 0.5-s footstep + base trajectory at 50–100 Hz (MPC over a reduced centroidal model, Di Carlo et al. 2018 “Dynamic Locomotion in the MIT Cheetah 3”).
3. Track the trajectory with a per-step whole-body QP at 200–500 Hz that solves for (q̈, λ_contact, τ) given friction cones, dynamics, and motion tasks.
4. Apply τ via high-bandwidth proprioceptive actuators (~1 kHz torque loop).

Stack used in published research: Pinocchio for dynamics, Crocoddyl or OCS² for trajectory opt, OSQP or HPIPM for the QP, ROS 2 for wiring. The same framework is now common to Spot, ANYmal C, Unitree Go2 / B2, Xiaomi CyberDog, and Ghost Robotics V60.

Franka Emika Panda — model-based torque control

The Franka Panda exposes a 1 kHz torque interface (franka_ros2 / libfranka) directly accessible to user code. The robot’s onboard controller computes M(q), C(q,q̇), g(q) from a factory-identified URDF and exposes them as part of the robot state, so user controllers can simply add their own feedforward and let the firmware handle stability monitoring. Reference open-source stacks (https://github.com/frankaemika/franka_ros2) include a Pinocchio-based gravity-compensated impedance controller as the canonical example. The published reachable torque envelope at the J2 shoulder is 87 N·m continuous / 175 N·m peak; gravity at full horizontal extension with 3 kg payload is ~60 N·m, so dynamic headroom is ~30 N·m continuous — enough for q̈ ≤ 4 rad/s² without violating limits. This headroom calculation is exactly what robot-arm sizing engineers do on every new payload integration.

MIT Cheetah / Mini Cheetah — proprioceptive force control

The MIT Cheetah series (Seok et al. 2015; Katz, Di Carlo & Kim 2019 on Mini Cheetah) intentionally used low-gear-ratio quasi-direct-drive (QDD) actuators (gear ratio 6:1, planetary) with custom high-pole-count BLDC stators. Low N gives the joint backdrivability and transparent force estimation from current sensing — no SEA spring, no in-line torque sensor, just I_motor · K_t = τ_joint − τ_friction. The resulting force-control bandwidth is > 200 Hz, enabling the high-bandwidth impact mitigation that lets a Mini Cheetah land a 1-m drop without breaking gearboxes.

The dynamics implications: with N = 6, reflected motor inertia 36 · 5·10⁻⁵ = 1.8·10⁻³ kg·m² is smaller than typical leg-link inertia, so the cross-coupling terms in M are large and full Newton-Euler dynamics matter (unlike high-N industrial arms where decoupled SISO PID works). All Cheetah-class control stacks therefore use whole-body QP control with Pinocchio-class dynamics under the hood; Cheetah, Cheetah 3, ANYmal C, Unitree A1/Go2, Stanford Doggo, MIT HERMES, and Unitree H1 all sit in this regime.

KUKA LBR iiwa 14 R820 — joint-torque sensors + impedance control

The KUKA LBR iiwa family is the canonical commercial example of a torque-controlled cobot. Each joint has an integrated strain-gauge joint-torque sensor measuring τ_meas directly (resolution ~0.01 N·m, bandwidth ~1 kHz). This obviates the need for proprioceptive torque estimation from motor current and enables high-bandwidth Cartesian impedance control:

$F_{\text{cmd}}=K_x(x_d-x)+D_x({\dot{x}}_d-\dot{x})+\Lambda (q){\ddot{x}}_d+\mu (q,\dot{q})\dot{x}+p(q)$

mapped through τ = J^⊤ F + null-space posture term. The KUKA Sunrise.OS controller computes M, C, g internally at 1 kHz on a per-arm-identified dynamic model, exposes them via the FRI (Fast Robot Interface) protocol, and lets the user inject feedforward terms. Performance envelope: Cartesian stiffness up to 5000 N/m linear, 300 N·m/rad rotational; force resolution at the tool ≈ 0.5 N. Used widely in assembly, hand-guidance teaching, surgical pre-clinical research (Cyberknife, Mako), and learning-from-demonstration research.

Quadrotor attitude dynamics

A quadrotor is a rigid body with 6 DOF whose only actuation is four scalar rotor thrusts, mapped to a 4-D wrench (collective thrust + 3 body torques) by the rotor geometry. The body-frame Euler equations are

$\dot{\mathbf{\omega }}=I_b^{-1}({\mathbf{\tau }}_b-\mathbf{\omega }\times I_b\mathbf{\omega })$

with I_b a diagonal inertia tensor (typical 1 kg drone: I_xx = I_yy ≈ 0.01 kg·m², I_zz ≈ 0.02 kg·m²). The gyroscopic ω × Iω term is small at hover and dominant in aggressive flight; PX4 and ArduPilot include it explicitly in their inner attitude loops. The differential-flatness property (Mellinger & Kumar 2011) maps a smooth desired position trajectory directly to required attitude and thrust, sidestepping inverse dynamics in the classical sense — but the same rigid-body Newton-Euler equations underpin the derivation.

Industrial pick-and-place — high-speed delta robot

Parallel manipulators such as the ABB IRB 360 FlexPicker or Adept Quattro use a closed-chain delta architecture: three or four parallel kinematic arms with a shared end-effector platform. Each arm has 1–2 actuated joints near the base; the remaining DOFs are passive. The end-effector is constrained to translate (and optionally rotate about one axis); typical accelerations are 5–10 g (50–100 m/s²) and cycle rates exceed 200 picks per minute.

Dynamics analysis is non-trivial because of the closed loops: the constrained manipulator equation is

$M(q)\ddot{q}+C(q,\dot{q})\dot{q}+g(q)=B\tau +J_c^{\top }\lambda ,J_c\ddot{q}+{\dot{J}}_c\dot{q}=0$

with B selecting actuated joints and λ the internal-loop multipliers. Solved as a coupled DAE or by Lagrangian elimination of dependent coordinates. The very small moving mass (≈ 0.5 kg platform + 50 g per arm) and short, stiff carbon-fibre arms put the first structural resonance at 80–120 Hz, well above the closed-loop bandwidth, enabling the order-of-magnitude higher accelerations than serial arms achieve.

8b. Case studies — emergent themes

Across all six platforms above, three patterns recur:

1. Whoever owns the dynamics model wins. UR, Franka, KUKA, ANYbotics and Boston Dynamics each publish or sell carefully identified dynamic models because controller performance is bottlenecked by model fidelity, not by control law. Even open-source efforts (mujoco_menagerie, ros_pinocchio, franka_ros2) standardise on factory-validated parameter sets.
2. Backdrivability is a design choice with deep dynamics consequences. High-N harmonic drives → decoupled SISO control suffices, but no transparent force sensing; low-N QDD → rich M coupling, must use full whole-body QPs, but proprioceptive force estimation from motor current is “free”.
3. Floating-base + contact + MPC has converged. Every modern legged-robot stack is some variant of Pinocchio + a QP solver + a centroidal-momentum or whole-body MPC layer. The differences are tuning, not architecture.

MuJoCo Menagerie — reference benchmark library

DeepMind publishes mujoco_menagerie as a curated, version-controlled set of MJCF models (https://github.com/google-deepmind/mujoco_menagerie). Each model includes geometry, mass / inertia, actuator and joint-limit calibration validated against published spec sheets and (where available) against the manufacturer’s reference dynamics. As of 2026-05 it covers 30+ commercial robots and is the de-facto reference for cross-paper benchmarks in RL, model-based control, and sim-to-real research. Every paper that claims “we ran on the Franka Panda” can now run on a fixed, reproducible MJCF — eliminating a decade of one-off URDF parameter inconsistency.

8c. End-to-end dynamics-pipeline checklist

A complete, deployable dynamics pipeline for a new robot — from CAD to controller — has these steps:

1. Author URDF / MJCF from CAD. Confirm joint axes, joint types, parent-child hierarchy. Validate against rendering in meshcat / mujoco viewer.
2. Audit inertias. Re-compute m, c, I from each mesh assuming uniform density; flag any link whose CAD inertia differs by > 30 %. Document mass of motors, harnesses, brackets not present in CAD.
3. Load into Pinocchio (or chosen library). Verify model.nq, model.nv, joint names match design intent. For floating base, ensure root joint is JointModelFreeFlyer (Pinocchio) or freejoint (MuJoCo).
4. Sanity-check g(q). At the home pose, compute g(q) and compare against hand calculation of static gravity torque. Discrepancies > 10 % indicate inertia or COM errors.
5. Sanity-check M(q). Confirm M is symmetric and positive-definite at random q samples. Compute cond(M); flag configurations with cond > 1000.
6. Run skew-symmetry test. For 100 random (q, q̇, x) triples, verify |x^⊤ (Ṁ − 2C) x| < 10⁻⁹ ‖x‖². Failure indicates a C-derivation bug.
7. Dynamic identification if intended for force / torque control. Generate exciting Fourier trajectory, record (q, q̇, q̈, τ) at ≥ 200 Hz for 60 s, regress for base parameters. Validate by predicting τ on a held-out trajectory.
8. Add friction model. Calibrate F_v, F_c per joint via slow constant-velocity sweeps. Document Stribeck regime if low-speed performance matters.
9. Choose integrator and Δt for the simulator. Default to semi-implicit Euler at 1 ms; tighten for stiff contact, loosen for off-line trajectory verification.
10. Define contact model (point contact with friction cone, soft compliance, hydroelastic). Pick friction coefficient μ to match deployment surface.
11. Build controller. Start with gravity-comp PD; layer in computed-torque feedforward only after baseline tracking is verified.
12. Verify on hardware with a low-speed trajectory and conservative gains. Compare logged τ against τ_model; iterate ID until residual < 5 %.

A team that skips any of steps 4–7 will spend the next six months debugging “controller bugs” that are actually model bugs.

Cost-of-fidelity table

Model fidelity	What’s included	Use case
Kinematic only	Forward / inverse kinematics, Jacobian	Position-controlled robot, slow motion
g(q) only	+ gravity vector	Hand-guidance, gravity-compensated PD
Rigid-body, no friction	+ M, C	Computed torque, OSC, basic MPC
+ Coulomb + viscous friction	+ F_v q̇ + F_c sgn(q̇)	High-fidelity tracking, force control
+ Joint flexibility	k_j and c_j per joint	High-bandwidth control, vibration analysis
+ Drivetrain dynamics	motor inertia, gear backlash	Direct-drive, harmonic-drive nonlinearities
+ Contact constraints	J_c, friction cone, λ	Manipulation, locomotion
+ Aerodynamics / hydrodynamics	external wrenches	Drones, AUVs
+ Soft-body / continuum	FEM coupling	Soft robots, cable-driven

Each row adds 5–20× implementation effort and 2–5× compute cost; pick the lowest row that solves your problem.

9. Cross-references

• kinematics-dh — kinematics is the prerequisite; q ↔ pose feeds the J that dynamics uses.
• motors-electric — torque sizing rules in §4 feed motor selection.
• pid-control — gravity-compensated PID is the simplest practical use of g(q).
• state-space-lqr — linearisation of M·q̈ + … around an equilibrium produces the state-space model.
• impedance-control — operational-space Λ = (JM⁻¹J^⊤)⁻¹ is the inertia term in impedance control.
• legged-robotics — floating-base dynamics + contact + MPC.
• trajectory-generation — output trajectories (q_d, q̇_d, q̈_d) feed inverse-dynamics feedforward.
• vibration-dynamics — mass-stiffness coupling, modal analysis, structural resonance.
• state-space-methods — control-theoretic framing of the manipulator equation.
• mpc-control — MPC requires forward dynamics in the prediction model.
• classical-control — PD on top of g(q) is classical-control + model-based feedforward.
• statics-fundamentals — the a = 0 limit of robot dynamics; manipulator statics is just g(q) + J^⊤F.
• mechanics-of-materials — link stiffness sets the flexible-joint terms ignored in rigid-body dynamics.
• robotics-control — URDF / MJCF / SDF inertia tags consumed by Pinocchio / MuJoCo / Drake.
• ros2-robotics-config — ROS 2 launch / parameter wiring around the dynamics pipeline.

10. Citations

1. Featherstone, R. Rigid Body Dynamics Algorithms. Springer, 2008. ISBN 978-0-387-74314-1. The canonical reference; RNEA, CRBA, ABA, articulated-body inertia, all derived in spatial-vector form. Every modern dynamics library is built on this.
2. Featherstone, R. and Orin, D.E. “Robot dynamics: equations and algorithms.” IEEE ICRA, 2000. Earlier, unified preview of the three core algorithms.
3. Luh, J.Y.S., Walker, M.W., and Paul, R.P.C. “On-line computational scheme for mechanical manipulators.” ASME J. Dynamic Systems, Measurement, and Control, 102(2):69–76, 1980. First O(n) recursive Newton-Euler formulation.
4. Khalil, W. and Kleinfinger, J.-F. “A new geometric notation for open and closed-loop robots.” IEEE ICRA, 1986. Modified DH convention that pairs with the dynamics recursions.
5. Khalil, W. and Dombre, E. Modeling, Identification and Control of Robots. Hermès Penton / Butterworth-Heinemann, 2002. ISBN 978-1-903996-13-8. Standard reference for dynamic-parameter identification.
6. Craig, J.J. Introduction to Robotics: Mechanics and Control, 4th ed., Pearson, 2017. ISBN 978-0-13-349796-5. Textbook standard; chapters 6–10 cover dynamics and computed torque.
7. Siciliano, B., Sciavicco, L., Villani, L., and Oriolo, G. Robotics: Modelling, Planning and Control, 2nd ed., Springer, 2009. ISBN 978-1-84628-641-4. The European reference; rigorous Lagrangian derivation, operational space, redundancy resolution.
8. Lynch, K.M. and Park, F.C. Modern Robotics: Mechanics, Planning, and Control, Cambridge University Press, 2017. ISBN 978-1-107-15630-2. Free online at https://hades.mech.northwestern.edu/index.php/Modern_Robotics. Body-frame and screw-theory treatment.
9. Murray, R.M., Li, Z., and Sastry, S.S. A Mathematical Introduction to Robotic Manipulation, CRC Press, 1994. ISBN 978-0-8493-7981-9. Free PDF from Caltech. Lie-group foundations.
10. Spong, M.W., Hutchinson, S., and Vidyasagar, M. Robot Modeling and Control, 2nd ed., Wiley, 2020. ISBN 978-1-119-52404-9. The control-engineer’s robotics textbook; rigorous passivity arguments.
11. Khatib, O. “A unified approach for motion and force control of robot manipulators: The operational space formulation.” IEEE J. Robotics and Automation, RA-3(1):43–53, 1987. Foundational operational-space control.
12. Carpentier, J., Saurel, G., Mansard, N., et al. “The Pinocchio C++ library — A fast and flexible implementation of rigid body dynamics algorithms and their analytical derivatives.” IEEE/SICE SII, 2019. The Pinocchio paper; published benchmarks against RBDL, RobCoGen, KDL.
13. Felis, M.L. “RBDL: an efficient rigid-body-dynamics library using recursive algorithms.” Autonomous Robots, 41(2):495–511, 2017. RBDL design and benchmarks.
14. Tedrake, R. et al. Drake: Model-based design and verification for robotics. https://drake.mit.edu/. Continuously updated software documentation; the Drake “Multibody Dynamics” notes are the clearest single derivation of contact-constrained dynamics.
15. Todorov, E., Erez, T., and Tassa, Y. “MuJoCo: A physics engine for model-based control.” IEEE/RSJ IROS, 2012. The MuJoCo paper; contact dynamics formulation.
16. Wensing, P.M., Wang, A., Seok, S., Otten, D., Lang, J., and Kim, S. “Proprioceptive actuator design in the MIT Cheetah: Impact mitigation and high-bandwidth physical interaction for dynamic legged robots.” IEEE T-RO, 33(3):509–522, 2017. Reference for high-bandwidth-actuator robot dynamics.
17. Di Carlo, J., Wensing, P.M., Katz, B., Bledt, G., and Kim, S. “Dynamic Locomotion in the MIT Cheetah 3 Through Convex Model-Predictive Control.” IEEE/RSJ IROS, 2018. The reference convex-MPC paper for quadrupeds.
18. Mason, M.T. Mechanics of Robotic Manipulation. MIT Press, 2001. ISBN 978-0-262-13396-8. Quasi-static manipulation, contact mechanics, complementarity — the bridge between rigid-body dynamics and manipulation.
19. Hwangbo, J., Lee, J., Dosovitskiy, A., et al. “Learning agile and dynamic motor skills for legged robots.” Science Robotics, 4(26):eaau5872, 2019. End-to-end RL on ANYmal; demonstrates sim-to-real with MuJoCo-class dynamics.
20. ISO 10218-1:2011 + Amd 1:2024. Robots and robotic devices — Safety requirements for industrial robots — Part 1: Robots. Standard governing force/power-limited cobot operation; gravity-comp control is a foundational enabler.
21. Khalil, W. “Dynamic modeling of robots using recursive Newton-Euler techniques.” ICINCO, 2010 (keynote). Concise modern survey of the recursion variants.
22. Atkeson, C.G., An, C.H., and Hollerbach, J.M. “Estimation of inertial parameters of manipulator loads and links.” Int. J. Robotics Research, 5(3):101–119, 1986. Foundational paper on regressor-based dynamic identification.
23. Swevers, J., Ganseman, C., Tükel, D.B., De Schutter, J., and Van Brussel, H. “Optimal robot excitation and identification.” IEEE T-RA, 13(5):730–740, 1997. Fourier-trajectory excitation design for ID.
24. Sousa, C.D. and Cortesão, R. “Physical feasibility of robot base inertial parameter identification: A linear matrix inequality approach.” Int. J. Robotics Research, 33(6):931–944, 2014. SDP-constrained identification ensuring physical inertia.
25. Orin, D.E. and Goswami, A. “Centroidal momentum matrix of a humanoid robot: Structure and properties.” IEEE/RSJ IROS, 2008. Foundational centroidal-dynamics paper.
26. Wensing, P.M. and Orin, D.E. “Improved computation of the humanoid centroidal dynamics and application for whole-body control.” Int. J. Humanoid Robotics, 13(1):1550039, 2016. The standard reference for centroidal-dynamics-based whole-body control.
27. Mellinger, D. and Kumar, V. “Minimum snap trajectory generation and control for quadrotors.” IEEE ICRA, 2011. Differential-flatness control for quadrotors; foundation for modern aggressive flight.
28. Marsden, J.E. and West, M. “Discrete mechanics and variational integrators.” Acta Numerica, 10:357–514, 2001. Theoretical basis for symplectic / variational integrators preserving energy in simulation.
29. Carpentier, J. and Mansard, N. “Analytical derivatives of rigid body dynamics algorithms.” Robotics: Science and Systems, 2018. The Pinocchio analytic-derivative paper.
30. Seok, S., Wang, A., Chuah, M.Y., Otten, D., Lang, J., and Kim, S. “Design principles for energy-efficient legged locomotion and implementation on the MIT Cheetah robot.” IEEE/ASME T-Mechatronics, 20(3):1117–1129, 2015. The proprioceptive-actuator design philosophy paper.

11. Glossary

• ABA — Articulated-Body Algorithm. O(n) forward dynamics.
• CRBA — Composite Rigid-Body Algorithm. O(n²) mass matrix.
• DH parameters — Denavit-Hartenberg parameters for joint kinematics; see kinematics-dh.
• Floating base — A robot whose root link is not fixed to the world (humanoids, quadrupeds, drones).
• g(q) — Joint-space gravity-torque vector.
• Inverse dynamics — Given q, q̇, q̈, compute τ.
• Forward dynamics — Given q, q̇, τ, compute q̈.
• M(q) — Mass / inertia matrix; n × n, symmetric, positive-definite.
• Manipulator equation — M·q̈ + C·q̇ + g = τ + J^⊤ F_ext.
• MPC — Model Predictive Control. Repeated finite-horizon optimisation.
• OSC — Operational-Space Control. Control in end-effector / task coordinates.
• QDD — Quasi-Direct-Drive actuator. Low gear ratio, high backdrivability.
• RNEA — Recursive Newton-Euler Algorithm. O(n) inverse dynamics.
• Spatial vector — 6-D twist or wrench; Featherstone’s unified formalism.
• τ — Joint torque vector.
• URDF / SDF / MJCF — Robot description formats; see robotics-control.
• ZMP — Zero Moment Point. Locomotion stability criterion.

12. Practical follow-ups

• Pick one dynamics library (recommend Pinocchio for control, MuJoCo for sim, Drake for verification) and master its API before mixing them in production.
• Use mujoco_menagerie as the source of truth for any benchmark or paper-reproduction work — eliminates one of the biggest reproducibility hazards in robotics literature.
• For any robot you intend to torque-control, allocate calendar time to dynamic identification — it almost always uncovers a 15–30 % discrepancy in the factory URDF that explains otherwise mysterious controller behaviour.
• Read Featherstone 2008 cover-to-cover at least once if you intend to write any non-trivial dynamics code; the algorithms are deep and the published-paper derivations skip steps that the book fills in.
• When debugging, always check passivity (Ṁ − 2C skew-symmetric) and energy conservation before suspecting controller bugs. These two tests catch 80 % of dynamics-implementation errors.