In MuJoCo a contact is a soft constraint: it produces force only by being violated a little, and how much force a given violation produces is set by parameters you choose. Three ideas carry the whole model:

Every contact number earlier in the course follows from these: the 0.367 mm a resting ball sinks (Lesson 0.1), why that does not depend on its mass, why a joint stop sinks further under a heavier arm (Lesson 5.1), and why inverse dynamics with contact is well defined (Lesson 7.2).

Visual intuition

The lab drops a ball on a plane with a frictionless contact whose parameters you set, and prints the resting penetration it measures beside the one the equations predict.

```lab simlab {“dock”: true, “title”: “Contact softness, predicted and measured”, “height”: 220, “xml”: “<body name='ball' pos='0 0 0.15'></body>”, “camera”: {“azimuth”: -90, “elevation”: 3, “distance”: 0.6, “target”: [0, 0, 0.06]}, “controls”: [ {“label”: “time constant”, “unit”: “s”, “apply”: “model.geom_solref[0] = value; model.geom_solref[2] = value”, “min”: 0.004, “max”: 0.05, “step”: 0.001, “value”: 0.02, “digits”: 3, “reset”: true}, {“label”: “damping ratio”, “apply”: “model.geom_solref[1] = value; model.geom_solref[3] = value”, “min”: 0.1, “max”: 2, “step”: 0.05, “value”: 1, “digits”: 2, “reset”: true}, {“label”: “impedance d0”, “apply”: “model.geom_solimp[0] = value; model.geom_solimp[5] = value”, “min”: 0.5, “max”: 0.99, “step”: 0.01, “value”: 0.9, “digits”: 2, “reset”: true}, {“label”: “impedance dw”, “apply”: “model.geom_solimp[1] = value; model.geom_solimp[6] = value”, “min”: 0.5, “max”: 0.99, “step”: 0.01, “value”: 0.95, “digits”: 2, “reset”: true}, {“label”: “ball mass”, “unit”: “kg”, “apply”: “model.body_mass[1] = value; mj.mj_setConst(model, data)”, “min”: 0.1, “max”: 10, “step”: 0.1, “value”: 1, “digits”: 1, “reset”: true}, {“type”: “button”, “label”: “Drop”, “run”: “mj.mj_resetData(model, data)”}], “readouts”: [ {“label”: “penetration, measured (mm)”, “expr”: “1000 * (0.05 - data.qpos[2])”, “digits”: 4}, {“label”: “penetration at rest, predicted (mm)”, “expr”: “(() => { const s = model.geom_solimp, tc = model.geom_solref[0], dr = model.geom_solref[1], d0 = s[0], dw = s[1], w = s[2], mid = s[3], p = s[4], g = -model.opt.gravity[2]; let r = 1e-4; for (let i = 0; i < 200; i++) { const x = Math.min(r / w, 1); const y = x >= 1 ? 1 : (x <= mid ? x ** p / mid ** (p - 1) : 1 - (1 - x) ** p / (1 - mid) ** (p - 1)); const d = d0 + y * (dw - d0); const k = d / (dw * dw * tc * tc * dr * dr); r = (1 - d) * g / (d * k); } return 1000 * r; })()”, “digits”: 4}, {“label”: “normal force (N)”, “expr”: “lib.normalForce()”, “digits”: 3}], “plots”: [{“ylabel”: “height (mm)”, “window”: 2, “traces”: [{“label”: “ball surface above plane (mm)”, “expr”: “1000 * (data.qpos[2] - 0.05)”}]}], “note”: “Both geoms get the same parameters, so no mixing takes place. The prediction solves the resting balance d k r = (1 − d) g for r, with d = d(r) from solimp; it assumes the ball is at rest, so compare it once the bouncing has stopped.”}


## Mathematics

### Resting penetration

At rest, $a_1 = 0$ and $v = 0$, so the constraint equation reduces to $d\,k\,r = (1-d)\,a_0$. With the time-constant format,

$$|r| = \frac{(1-d)\,a_0\,d_w^2\,\text{timeconst}^2\,\text{dampratio}^2}{d^2}, \qquad d = d(r),$$

an implicit equation because the impedance depends on the violation; a few fixed-point iterations solve it. For the default parameters it gives $|r|$ = 0.3672 mm, with $d$ = 0.9135, exactly what MuJoCo measures (output 1).

> [!derivation] Why the penetration does not depend on mass, but a joint stop's does
> For a body resting on a plane, the unconstrained acceleration along the normal is $a_0 = g$ whatever the mass, so $r$ is the same for 0.1, 1 and 10 kg (output 1). For a joint stop, $a_0$ is the joint's acceleration if the stop were removed: the gravity torque divided by the effective inertia about that joint. Making one link heavier raises both, but not in proportion, which is the +26% of Lesson 5.1's quiz. The parameters set a *reference motion*, not a stiffness in newtons per metre.

### From parameters to force

The regularizer is $R_{ii} = \frac{1-d_i}{d_i}\hat A_{ii}$, with $\hat A$ a diagonal approximation of $A = JM^{-1}J^\top$ computed once at `qpos0`. When the unconstrained minimizer of the solver's problem is admissible, the force is

$$f = (A + R)^{-1}(a_\text{ref} - a_0), \qquad a_1 = A(A+R)^{-1}a_\text{ref} + R(A+R)^{-1}a_0 ,$$

so the constrained acceleration interpolates between the reference and the unconstrained one: hard as $R \to 0$ ($d \to 1$), absent as $R \to \infty$ ($d \to 0$). Output (2) reads every quantity for the resting ball out of MuJoCo (`efc_KBIP`, `efc_aref`, `efc_R`, `efc_force`) and recomputes each from these formulas. Two details show up: `efc_KBIP` stores $k$ without the factor $d(r)$, which MuJoCo applies when it forms $a_\text{ref} = -b v - d\,k\,r$; and for this free ball $\hat A = A = 1/m$, so the regularizer matches its definition exactly.

### The optimization MuJoCo solves

All constraint forces come from one convex problem, a generalization of Gauss's principle of least constraint:

$$\min_{x,\,y}\ \|x - M^{-1}(\tau - c)\|_M^2 + \|y - a_\text{ref}\|^2_{R^{-1}} \quad \text{s.t.}\ \ J_\mathcal{E}x - y_\mathcal{E} = 0,\ \ J_\mathcal{F}x - y_\mathcal{F} = 0,\ \ J_\mathcal{C}x - y_\mathcal{C} \in \mathcal{K}^*,$$

where $x$ is the joint acceleration, $y$ a slack in constraint space, $\mathcal{K}^*$ the dual of the friction cones, and the norm on $y$ becomes a Huber function for friction loss. The first term keeps the motion close to what it would be without constraints, measured in the inertia metric; the second pulls each constraint towards its reference acceleration, as strongly as $R^{-1}$ says. Because the problem is convex, its solution is unique and the solvers (Newton by default, also CG and PGS) find it reliably; because $R > 0$, it is never the hard, possibly ill-posed rigid-contact problem.

### Bounce

The damping ratio is the same in the effective constraint-space oscillator as in `solref` (MuJoCo absorbs one factor of $d$ into $k$ but not into $b$ for that reason). Lower ratios bounce: output (3) shows rebound heights of 0, 7%, 33% and 57% of the drop for ratios 1, 0.5, 0.2 and 0.1, and a fully elastic bounce with zero damping in the direct format. The linear-oscillator formula for restitution underestimates these because a contact releases as soon as its force would turn negative, before the spring has returned its energy symmetrically (inference from the unilateral constraint).

## Implementation

```io
INPUT: a ball on a plane with frictionless contacts
PROCESS: predict and measure resting penetration across parameters, masses and mixed parameters; read MuJoCo's constraint quantities and recompute them; measure bounce against the damping ratio; probe refsafe
OUTPUT: printed tables

```python file=examples/l9_2_soft_contact.py “"”Lesson 9.2: the soft contact model: solref, solimp, the regularizer and the solver.

INPUT a ball on a plane, written below (frictionless contacts where friction would only add rows to the comparison) PROCESS (1) resting penetration predicted from the model’s equations and measured, for several time constants, a constant impedance, the direct solref format and three masses; (2) MuJoCo’s own constraint quantities for the resting contact: impedance d, stiffness k, damping b, reference acceleration, regularizer R, and the force (A + R)^-1 (aref - a0) against efc_force; (3) bounce: rebound height against the damping ratio; (4) refsafe: a time constant below twice the timestep is raised to 2h OUTPUT printed tables

Core concept

Visual intuition

Debugging

Exercise

Challenge

Research connection

Next