MuJoCo, and how it becomes robot.zig

What a reduced-coordinate robotics simulator actually does — and the parts of it worth having in zimr, given that zimrphysics.zig already exists.

That is the whole of MuJoCo in one line, and it is not the equation zimrphysics solves. Understanding why is most of the work; the rest of this document is the consequences.

zimrphysics.zig is a Jolt port. It thinks in maximal coordinates: every body carries a full 6-DOF pose in world space, joints are constraints that remove unwanted freedom, and a solver iterates until the constraint violations are small. That is the right model for a game — ragdolls, crates, vehicles, thousands of loose objects.

MuJoCo thinks in generalized coordinates. A robot arm with four hinges has exactly four numbers, not four bodies × seven numbers with constraints holding them together. The joints are not enforced; they are structurally impossible to violate, because the unwanted freedom was never represented. This is the difference between a simulator that drifts and one that cannot.

1. The two worldviews, side by side

Neither is better. They are good at different things, which is exactly why robot.zig should be a sibling of zimrphysics rather than a replacement or a rewrite of it.

2. The model / data split

MuJoCo's single most transferable architectural idea has nothing to do with physics. It splits everything in two:

• mjModel — constant. The robot's structure: tree topology, joint types and axes, body inertias, actuator parameters, solver settings. Compiled once from XML, then never written. One mjModel can be shared by many simulations.
• mjData — everything that changes. Positions, velocities, all intermediates, and every derived quantity. One per simulation instance.

This is why MuJoCo can run thousands of parallel rollouts of the same robot cheaply (the whole premise of MJX, its JAX/GPU sibling), and why its API is a flat list of free functions mj_xxx(m, d) that read m, write d, and hold no state of their own.

3. The kinematic tree, and why nq ≠ nv

Bodies form a tree. Each body carries zero or more joints, and each joint contributes coordinates. There are exactly four joint types:

A body with no joint is welded to its parent — free and correct, no constraint required. This is how you build a rigid multi-geom link.

Position and velocity have different dimensions. Orientations live on the unit sphere in 4D (SO(3)); velocities live in its 3D tangent space. So qpos has length nq, qvel has length nv, and nq > nv in any model containing a ball or free joint. You cannot subtract two qpos vectors to get a velocity — MuJoCo provides mj_differentiatePos for that, and mj_integratePos for the reverse. Every quaternion-aware operation in the codebase exists because of this one fact.

4. Smooth dynamics: five passes over a tree

The "smooth" part of MuJoCo — everything except constraints — is five tree traversals, and it is much smaller than you would guess.

mj_kinematics — forward pass

Walk the tree parent→child, composing each body's local pose with its joint's contribution, producing global position and orientation for every body, geom and site. Normalizes quaternions on the way through.

mj_comPos — center-of-mass frames

MuJoCo does its dynamics in global frames (so collision detection shares them) but translated to each kinematic subtree's center of mass. That translation is purely for floating-point accuracy: inertias expressed about a distant origin lose precision badly. This produces cdof (each DOF's motion axis as a 6D spatial vector) and cinert (each body's 10-parameter inertia).

mj_crb — the mass matrix

Composite Rigid Body. Backward pass accumulating each subtree's total inertia into its parent, then for each DOF i, walk up the DOF ancestor chain writing M[i][j] = cdof_j · (crb_i × cdof_i). The entire algorithm is about sixty lines of C. Two things fall out of it:

• M is sparse, and its sparsity pattern is exactly the ancestor relation in the tree — M[i][j] is nonzero only if DOF i and j lie on a common root path. This is fixed by the model, not by the state.
• Armature (rotor inertia reflected through a gearbox) is simply added to the diagonal. One number per DOF, and it is the single most effective stabilizer for stiff geared robots.

mj_factorM / mj_solveM — sparse L^TDL

Because M is needed as M-1x repeatedly, it is factorized once per step into L^TDL, in an elimination order that preserves the tree sparsity exactly — no fill-in. Then every M-1x is two sparse back-substitutions.

mj_rne — bias forces

Recursive Newton-Euler run with acceleration set to zero, which yields exactly c(q,v): Coriolis, centrifugal and gravitational forces combined. Forward pass propagates velocities, backward pass accumulates forces. RNE is also how inverse dynamics is computed, with the real acceleration instead of zero.

5. Jacobians, the bridge between spaces

J maps joint velocities into constraint (or task) space, and JT maps forces back the other way. This one object serves several jobs that look unrelated until you see them written down:

• a contact is a row of J in the contact normal direction;
• a joint limit is a row of J that is just a unit vector on that DOF;
• an equality constraint is ∂r/∂q for its residual;
• an actuator's moment arm is ∇l(q), the gradient of its transmission length;
• an end-effector task ("move the gripper this way") is the body Jacobian.

Getting Jacobians right once buys all five. In a tree, the Jacobian of a point is assembled by walking that body's ancestor chain and taking each DOF's contribution — the same chain the CRB walk uses.

6. The soft constraint model

This is MuJoCo's real signature, and the thing most worth stealing.

Every other rigid-body engine treats contact as a hard complementarity problem: either the bodies touch and there is force, or they separate and there is none. That is an LCP, it is NP-hard with friction, and it has no unique solution in the general case.

MuJoCo drops strict complementarity and makes every constraint soft, giving a convex optimization problem instead — always solvable, always unique, differentiable. Contacts can push back slightly before touching, and penetrate slightly under load. That is not a bug being tolerated; it is the model.

Each scalar constraint gets three numbers:

• impedance d ∈ (0,1) — how hard the constraint is. d→1 is rigid, d→0 is absent. It interpolates the achieved acceleration between the reference and the unconstrained one.
• stiffness k and damping b — the reference acceleration is a spring-damper pulling the residual to zero. Authored as solref = (timeconst, dampratio), so you specify how fast and how bouncy, not raw gains.

The crucial trick is the Â in the regularizer: R is scaled by an approximation of the constraint-space inertia, precomputed at the reference pose. That is what makes solref mean the same thing on a 1 g finger and a 100 kg torso — the penetration under load becomes independent of effective mass. Without that scaling, every contact would need hand-tuned gains. It is a small idea with an enormous usability payoff.

Constraints are assembled in a fixed order — equality, friction loss, limit, contact — into one system of nc rows, and one solver handles all of them uniformly. Limits are literally frictionless contacts internally.

7. Actuators: three independent pieces

MuJoCo factors an actuator into three parts you configure separately, and the combinations give you everything from a DC motor to a muscle:

1. Transmission — where the force goes. Attach to a joint, a tendon, a site (Cartesian thrust: jets, propellers), a body (adhesion, vacuum grippers), or a slider-crank. Each defines a scalar length l(q); its gradient ∇l is the moment-arm vector that maps scalar force to joint forces.
2. Activation dynamics — optional internal state w with first-order dynamics, making the system third-order. Types: integrator (ẇ = u), filter (ẇ = (u−w)/τ), filterexact (the same, integrated analytically so it cannot diverge when τ < timestep), muscle, and DC-motor/PID variants.
3. Force generation — affine: p = a·(w or u) + b₀ + b₁l + b₂l̇. Choosing gain and bias turns the same machinery into a torque motor (b=0), a position servo (b₁ = −kp), or a velocity servo (b₂ = −kv).

Because force generation is affine, inverse dynamics can recover the controls from the required force with a pseudo-inverse. That property is the reason the whole model is deliberately kept affine.

8. Tendons

A tendon is a scalar length defined over the model — either a fixed linear combination of joint coordinates, or a spatial path through sites, wrapping around spheres and cylinders with proper minimal-length geodesics. It can have its own spring, damper, friction and limits, and actuators can pull on it.

Two joints coupled by a fixed tendon is how you build a differential, a car's steering linkage, or a coupled finger — without introducing a closed kinematic loop. In the car model above, "forward" and "turn" are tendon actuators over the two wheel joints, so two controls drive a differential drive directly.

9. Sensors, and why they are in the engine

Sensors are declared in the model and evaluated at three points in the pipeline — after position, after velocity, after force. That is not a convenience feature: a physically meaningful accelerometer or force/torque reading is only available at the right stage, and computing it afterwards from finite differences would be both wrong and noisy. For a robotics simulator the sensor set is part of the physics, not part of the application.

10. Integrators

Implicit-in-velocity is worth understanding because it is cheap and it buys a lot. Instead of v←v+h·a(v), it solves v←v+h·a(vnew) by one Newton step, giving (M − hD)vnew = … where D = ∂(forces)/∂v. Velocity-dependent forces — damping, fluid drag, velocity servos — stop being an explicit-integration stability hazard. You get to raise damping without lowering the timestep.

11. The pipeline

26 stages, strictly ordered, each consuming the last. The numbers are the point: this is a dependency chain, and the reason MuJoCo's API exposes mj_step1/mj_step2 is so a controller can run between stage 19 and 20 — after everything derived from state, before anything derived from force.

■ port into robot.zig ■ reuse zimrphysics / zimrmath ■ out of scope

12. Making it Zig, and making it zimr

A faithful transliteration of MuJoCo would be a mistake. Most of its complexity is a consequence of two constraints zimr does not have: it is C, and its models arrive as XML at runtime. Both facts force the same thing — everything must be dynamically sized, so everything becomes an index into a flat array with a runtime-built sparsity table. That is where mjModel's 900 lines of parallel arrays come from.

The central bet: the model is comptime

A robot's structure does not change while it runs. Its tree, its joint types, nq, nv, and therefore the entire sparsity pattern of M, are all known before the program starts. Zig can know them at compile time.

const Arm = robot.Model(.{
    .bodies = .{
        .{ .name = "base",     .parent = null,       .inertia = … },
        .{ .name = "upper",    .parent = "base",     .joint = .{ .hinge = .{ .axis = zm.vec3(0,0,1) } } },
        .{ .name = "fore",     .parent = "upper",    .joint = .{ .hinge = .{ .axis = zm.vec3(0,1,0) } } },
        .{ .name = "hand",     .parent = "fore",     .joint = .{ .ball = .{} } },
    },
    .actuators = .{
        .{ .name = "shoulder", .on = .{ .joint = "upper" }, .kind = .{ .position = .{ .kp = 40 } } },
    },
});

// nq, nv and the sparsity of M are comptime constants of the TYPE.
var data: Arm.Data = .init;          // fixed-size arrays. no allocator.
Arm.step(&data, dt);

What this buys, concretely:

• No allocator in the hot path, ever. Data is a plain struct of [nv]f32-shaped fields. It can live on the stack, in an arena, or in a GPU buffer.
• The sparsity tables disappear. MuJoCo carries M_rownnz, M_rowadr, dof_parentid and walks them at runtime. With a comptime tree, the ancestor chains are known and the CRB inner loop can be an inline for over exactly the nonzero entries. The sparse structure stops being data and becomes control flow.
• Wrong models fail at compile time. A parent that does not exist, an actuator on a missing joint, a ball joint given a 1-D gear — all @compileError with a message, not a runtime validation pass.
• It is the zimr house style already. Designated struct literals as the API, Options-structs, no builders. Compare zp.createBody(world, .{ … }): the same shape, one level more comptime.

Where robot.zig meets zimrphysics.zig

This is the hinge that makes the two feel like one library rather than two. MuJoCo contains a whole collision engine — GJK, box-box, convex, height fields, SDFs, broad-phase. zimrphysics already has all of that, device-tested, and it is the single largest chunk of MuJoCo by volume.

So robot.zig should not port collision at all. It should consume it:

// zimrphysics owns: shapes, GJK/SAT, the BVH broad-phase, manifolds.
// robot.zig owns:   the tree, M, c, J, actuators, the constraint solve.
//
// The seam is one function: a contact manifold in, constraint rows out.
const contacts = zp.collideAll(&scene, robot_geom_poses, scratch);
robot.addContactConstraints(&data, contacts);

The same seam works in the other direction later: a robot can be an actor in a zimrphysics world, pushing loose objects around, with the robot solved in generalized coordinates and the debris solved in maximal ones. That is a genuinely good architecture and it is what most robotics stacks converge on anyway.

Likewise zimrmath already provides Vec, Quat, Mat, and the quaternion vocabulary. robot.zig adds only what is missing: spatial (6D) motion and force vectors, and the 10-parameter inertia. Those are genuinely new — and they belong in robot.zig rather than zimrmath, because zimrmath is shader-safe and must stay small.

Naming and idiom

MuJoCo's names are terse C (qfrc_bias, cdof, efc_J). zimr's rules say otherwise. A translation table, applied consistently:

Three ideas that are more natural in zimr than in MuJoCo

1. The control hook is a comptime function, not a callback pointer. MuJoCo's mjcb_control is a global function pointer — a mutable module global, which zimr bans outright. Passing the controller as a comptime parameter to step is faster and legal: it inlines, and there is no global.
2. CPU/GPU duality is already a zimr pattern. MJX exists because batched rollouts on an accelerator are the whole game for modern robot learning. zimr has kompute and a comptime-shaped Data struct that maps directly onto a GPU buffer, and a house tradition of side-by-side CPU|GPU demos. A batched rollout demo is a natural flagship, not a stretch.
3. Determinism is testable headlessly. This whole subsystem is pure CPU maths with no device dependency — the rarest thing in zimr. Every algorithm here can be pinned by a host test, and the strongest tests are the physical invariants: energy conservation under RK4 with no damping, momentum conservation for a free body, and forward-then-inverse dynamics round-tripping to the original force. That last one is a complete self-check of M, c and J at once, and MuJoCo ships it as a diagnostic (stage 25) precisely because it is so sharp.

13. The build plan for robot.zig

Ordered so that every phase ends somewhere demonstrable, and each one is verifiable on the host before the next begins.

Phases 0–4 are the irreducible core and are worth having even if nothing after them is ever built: they are what turns zimr from "a game physics engine" into "a thing that can simulate a robot arm correctly".

What is deliberately not in robot.zig

• Collision detection — zimrphysics owns it. This is the point.
• XML / MJCF / URDF parsing — the model is Zig. A loader can come later, as a generator that emits Zig source.
• Deformables (flex), SDF collision, plugins, muscles — large, specialized, and not on the path to a working robot.
• Rendering — zimr already draws; the robot is just poses.

14. The one risk worth naming up front

The soft constraint solver (phase 6) is the part that looks small in the equations and is not. MuJoCo's engine_core_constraint.c is 3,500 lines, it has four solver algorithms, two friction cone types, and the Â diagonal approximation that makes solref behave sanely across mass scales. Everything before it is textbook and well-conditioned; this is the part with the taste in it.

The mitigation is ordering, and it is already in the plan: phases 0–5 produce a genuinely useful articulated-body simulator with no constraint solver at all — forward and inverse dynamics, Jacobians, actuators, integrators. A robot arm moving under torque control needs none of phase 6. Start the solver with a single algorithm (projected Gauss-Seidel over a pyramidal cone), get joint limits working, and only reach for the elliptic cone and the Newton solver if something demands it.

Sources: MuJoCo doc/computation/index.rst (the authoritative algorithm description), src/engine/engine_core_smooth.c, engine_core_constraint.c, engine_forward.c, include/mujoco/mjtype.h, mjmodel.h, mjdata.h, and the model/car and model/humanoid examples. Read alongside src/zimrphysics.zig's CharacterVirtual header for the house conventions on porting an upstream engine and recording the divergences.