Core concept

Regression (Lesson 18.1) needs the dynamics written linearly in the parameters and needs velocities and accelerations differentiated from noisy positions. Identification by simulation needs neither. Pick parameters, simulate the recorded commands from the recorded initial state, compare the simulated sensor readings with the measured ones, and adjust the parameters to shrink the difference. Any parameter the simulator exposes can be fitted this way: contact softness, actuator delay, friction loss. The price is a simulation per evaluation, and a cost that is only as smooth as the simulation. This lesson fits the 18.1 pendulum with MuJoCo’s own toolbox, adds an actuator delay the model does not have, and then tries a contact parameter, where the method’s weaknesses show.

[!established] mujoco.sysid MuJoCo 3.14.0 ships a system identification toolbox, mujoco.sysid, installed with the extra pip install "mujoco[sysid]==3.14.0" (the course’s pip install -e ".[sysid]"). Without the extra, import mujoco.sysid still succeeds but the module exports nothing; check for mujoco.sysid.optimize. A Parameter carries bounds and a modifier that writes its value into an MjSpec. ModelSequences pairs the spec with measured initial states, commands and sensor readings, and build_residual_fn turns them into a function that recompiles the model, runs mujoco.rollout and returns the sensor errors. optimize minimizes those errors by box-bounded nonlinear least squares (mujoco.minimize.least_squares, a regularized Gauss-Newton method with a finite-difference Jacobian), and calculate_intervals turns the Jacobian at the optimum into 95% intervals.

Visual intuition

There is no browser lab for this lesson: the toolbox runs only in Python. The 18.1 lab is the regression version of part (1); compare its estimates with the table below. The picture to keep is a cost surface. For the pendulum the surface is a smooth bowl, so the optimizer reaches the bottom in a few dozen simulations from any reasonable start. For a bouncing ball the same surface is jagged, and the optimizer stops wherever its local slope says it should.

Mathematics

Simulation error

For parameters $\theta$, simulate the recorded commands $u_{0:N}$ from the recorded initial state $x_0$ and read the predicted sensors $\hat y_k(\theta)$. Minimize

Each evaluation is one simulation, and a finite-difference Jacobian costs one simulation per parameter, so a fit takes tens to hundreds of simulations. The errors $\hat y_k - y_k$ are not linear in $\theta$, so there is no closed form and no guarantee that a minimum the optimizer finds is the global one.

(1) The pendulum: no derivatives, no friction model to get wrong

On the 18.1 record (same hidden values, same multisine, same noise), starting from 1.0 kg, 0.01 and 0.01:

The correctly stamped fit recovers every parameter within its interval, and on the held-out excitation its predicted angle differs from the measurement by 0.50 mrad RMS, which is the measurement noise (0.5 mrad) and what the true parameters also give. The 18.1 regression on the same record gave 2.57 mrad. Simulation error has no derivative noise and uses MuJoCo’s own friction-loss constraint, so neither of the two biases of 18.1 can occur.

This is an inverse crime: the data come from the model class being fitted, so exactness is expected. What it shows is that the method adds no bias of its own; the errors of 18.1 came from differentiation and from the friction model, not from the noise. Real data afford neither luxury, which is what parts (2) and (3) are about.

[!established] Stamp samples the way the simulator does For step $k$, mujoco.rollout returns the state after the step, at time $(k+1)h$, and the sensor values computed during the step, which MuJoCo evaluates from the state at time $kh$ before integrating: after two steps from rest the pendulum’s qpos is 3.504e-05 rad while its angle sensor reads 1.168e-05, and mj_forward brings the sensor to 3.504e-05. The toolbox times the sensor readings by the state’s time. A measured angle taken at $kh$ and stamped $kh$ is therefore compared with a prediction one step older. Measured: that 2 ms offset moves the friction loss by 9% and the damping by 3.5%, with intervals 10 to 20 times narrower than the errors. Stamping the measurements one step later removes it.

(2) An actuator delay the model does not have

The measured pendulum now applies its torque 10 ms late. MuJoCo 3.14.0 models this per actuator: delay (s), a history buffer of nsample past commands, and interp, the interpolation between stored commands (interp="zoh" holds the last sample, "linear" interpolates; in MjSpec the same choice is the integer 0 or 1). Fitted without a delay, and then with the delay as a fourth parameter:

Without the delay the fit absorbs the lag into the other parameters: damping 18% low and friction loss 48% high, with intervals of 1.5% or less. A delayed torque $\tau(t-d) \approx \tau(t) - d\,\dot\tau(t)$ adds a term the model can only mimic with other velocity-dependent terms (inference). The estimates are confidently wrong, and the held-out error, about 23 times the noise, is the only sign.

The second row shows why a delay is awkward to fit. Driving the actuator with a ramp and measuring how far its force lags behind gives, with linear interpolation and a 2 ms step:

Any delay between zero and one step acts as exactly one step, so the cost is flat there. The optimizer, started at 5 ms with a lower bound of 0, ended at the bound, where the Jacobian column of the delay is zero and its interval infinite. With the lower bound at one step, it finds 9.99 ms and the other three parameters to within 0.15%. With interp="zoh" the lag is rounded up to whole steps (1, 1, 2 and 2 steps for 1 to 4 ms), so the cost is a staircase in the delay and a finite-difference gradient is zero almost everywhere; fit delays with interpolation on.

(3) Contact: a jagged cost

A ball dropped from 0.5 m bounces with a contact softness solref = (time constant 0.015 s, damping ratio 0.25) to identify from 1.2 s of height measured with 0.5 mm noise. The height error along each parameter, with the other at its true value:

The truth is a sharp minimum at the noise level, and along the time constant the cost jumps up and down between neighbouring values. A plausible cause (inference): each bounce ends on a whole step, and a small change in softness moves the step on which the ball leaves the floor, which shifts every later bounce. A gradient taken by finite differences sees only the local wiggle. Fits from five starting points ended at five different places:

None is near the truth, whose error is 0.50 mm. The first stopped at the 4 ms lower bound, which is there because MuJoCo clamps a time constant below two timesteps: at a damping ratio of 0.25 a time constant of 1 ms and one of 4 ms give the same highest point after the drop, so the cost is flat below 4 ms. That row also shows where the contact stops being physical. The ball, dropped from 500 mm, rises to 4001 mm at 1 and 4 ms, to 1324 mm at 5 ms, and no higher than its drop from 6 ms up: a stiff, lightly damped contact this close to the step size adds energy. The fourth fit stopped in that region, with the ball bouncing to 666 mm. A 40 × 40 grid over both parameters on log scales did better, (0.0127, 0.2687) with a 1.28 mm error, and still missed the basin of the truth.

Fit what the data constrain

The first two apex heights of the measured bounces, 152.5 and 53.8 mm, are a smooth function of the parameters. Fitting them alone, by a scan and a one-dimensional refinement of the damping ratio with the time constant fixed, gives 0.2473, 0.2514 and 0.2467 for time constants fixed at 0.01, 0.02 and 0.03 s: within 1.3% of the true 0.25 whichever value is fixed. Bounce heights measure the restitution, which these numbers show depends mainly on the damping ratio; the time constant is barely constrained by them. On a held-out drop from 0.3 m:

The apex fits predict the held-out drop at least as well as the best raw fit, with a 1-D search instead of a 2-D one, and the remaining error depends on the time constant they could not identify. Pinning that down takes data that depend on it, such as contact duration or penetration, which needs a faster sensor or a force measurement (inference).

[!recommendation] A checklist for identification by simulation Stamp measurements by the simulator’s convention, then check the fit on noise-free simulated data before using real data. Bound every parameter to its physical and numerical range: below one step for a delay and below two steps for a contact time constant the cost is flat, and some ranges add energy. For contact, plot the cost along each parameter before trusting a gradient method; if it is jagged, fit smooth features of the motion, or search on a grid and then refine. Intervals from the Jacobian measure noise only: they were at most 1.5% wide around estimates 18% and 48% wrong. Validate on held-out motion and compare the error with the noise level.

Implementation

INPUT: pendulum.xml with hidden parameters and the 18.1 record; the same pendulum with a 10 ms actuator delay; a bouncing ball with hidden solref
PROCESS: mujoco.sysid fits with two time stampings; delay as a parameter with two lower bounds; the contact cost along each parameter, five starts, a grid, an apex-height fit; validation on held-out motion
OUTPUT: printed tables

The backend prints a line per iteration, and optimize(verbose=False) does not silence it, because the backend’s own verbose flag cannot be passed through optimize’s keyword arguments; the script captures standard output and reads the number of simulations from it. In 3.14.0 a joint’s damping in MjSpec is a polynomial (three coefficients, the linear one first), so it is written as [b, 0, 0].

```python file=examples/l18_2_sysid_simulation.py “"”Lesson 18.2: identification by simulation, with MuJoCo’s sysid toolbox (Level 18).

INPUT pendulum.xml with hidden mass 1.3 kg, damping 0.08 N m s/rad and friction loss 0.05 N m, driven and measured as in Lesson 18.1 (multisine up to 2.5 N m, 10 s, 0.5 mrad angle noise); the same pendulum with a 10 ms actuator delay that the model does not have; a ball dropped on a floor with hidden contact softness solref = (0.015 s, 0.25) PROCESS (1) mujoco.sysid simulates the record and minimizes the angle error over (mass, damping, friction loss): with each sample stamped at the time it was measured, then stamped the way mujoco.rollout stamps its sensors; next to the 18.1 regression on the same record; (2) the delayed pendulum, fitted without a delay and with the delay as a fourth parameter (lower bound 0, then one timestep); how MuJoCo applies a delay shorter than one step; validation on a held-out excitation; (3) contact: the height error over 1.2 s of bounces along each solref parameter, fits from five starting points, a 40 x 40 grid; then a fit of the bounce apex heights alone, with the time constant fixed; validation on a drop from 0.3 m OUTPUT printed tables

Run: python examples/l18_2_sysid_simulation.py (about 2 minutes) “””

requires: sysid

import contextlib import io import os import re

import mujoco import numpy as np from mujoco import sysid from scipy.optimize import minimize_scalar from scipy.signal import find_peaks, savgol_filter

from mjcourse import model_path

FAST = os.environ.get(“MJC_FAST”) == “1” TRUE = {“mass”: 1.3, “damping”: 0.08, “frictionloss”: 0.05} GUESS = {“mass”: 1.0, “damping”: 0.0, “frictionloss”: 0.0} L, I_C, G, NOISE, H = 0.5, 1e-4, 9.81, 0.0005, 0.002

def pendulum_spec(p: dict, delay: float | None = None) -> mujoco.MjSpec: “"”pendulum.xml with parameters p; with delay, the actuator keeps a 20-sample history (linear).””” spec = mujoco.MjSpec.from_file(str(model_path(“pendulum”))) spec.body(“pole”).mass = p[“mass”] spec.joint(“hinge”).damping = [p[“damping”], 0.0, 0.0] # polynomial damping: linear term first spec.joint(“hinge”).frictionloss = p[“frictionloss”] if delay is not None: act = spec.actuator(“torque”) act.nsample, act.interp, act.delay = 20, 1, delay return spec

def multisine(t: np.ndarray, amplitude: float, seed: int) -> np.ndarray: rng = np.random.default_rng(seed) freqs, phases = rng.uniform(0.2, 3.0, 5), rng.uniform(0, 2 * np.pi, 5) return amplitude * sum(np.sin(2 * np.pi * f * t + p) for f, p in zip(freqs, phases)) / 5

def simulate(spec: mujoco.MjSpec, torque: np.ndarray) -> np.ndarray: “"”Angle at t = k h, measured before step k applies torque[k] (as in Lesson 18.1).””” model = spec.compile() data = mujoco.MjData(model) q = np.empty(len(torque)) for k, tau in enumerate(torque): q[k] = data.qpos[0] data.ctrl[0] = tau mujoco.mj_step(model, data) return q

def setter(path: str): “"”Modifier writing a scalar parameter into the spec (sysid calls it before each compile).””” def apply(spec, p): v = float(p.value[0]) if path == “mass”: spec.body(“pole”).mass = v elif path == “damping”: spec.joint(“hinge”).damping = [v, 0.0, 0.0] elif path == “frictionloss”: spec.joint(“hinge”).frictionloss = v else: spec.actuator(“torque”).delay = v return apply

def run_optimize(params, residual_fn) -> tuple[sysid.ParameterDict, object, int]: “"”sysid.optimize, with the backend’s per-iteration log captured; returns the residual evaluations.””” log = io.StringIO() with contextlib.redirect_stdout(log): opt, result = sysid.optimize(initial_params=params, residual_fn=residual_fn, optimizer=”mujoco”, verbose=False) found = re.search(r”Residual evals: (\d+)”, log.getvalue()) return opt, result, int(found.group(1)) if found else -1

def fit_pendulum(t, torque, q, stamp_shift: float, delay_bounds=None) -> tuple[dict, np.ndarray, int]: “"”Fit (mass, damping, friction loss[, delay]) to a measured angle; returns estimate, 95% half-widths, evals.””” spec = pendulum_spec(GUESS, None if delay_bounds is None else 0.0) model = spec.compile() data = mujoco.MjData(model) x0 = sysid.create_initial_state(model, data.qpos, data.qvel, data.act) sequences = sysid.ModelSequences( “pendulum”, spec, “record”, x0, sysid.TimeSeries(t, torque[:, None]), sysid.TimeSeries.from_names(t + stamp_shift, q[:, None], model, names=[“angle”]), allow_missing_sensors=True) # the angular-velocity sensor is not measured params = sysid.ParameterDict() params.add(sysid.Parameter(“mass”, 1.0, 0.2, 5.0, modifier=setter(“mass”))) params.add(sysid.Parameter(“damping”, 0.01, 0.0, 1.0, modifier=setter(“damping”))) params.add(sysid.Parameter(“frictionloss”, 0.01, 0.0, 0.5, modifier=setter(“frictionloss”))) if delay_bounds is not None: params.add(sysid.Parameter(“delay”, 0.005, delay_bounds[0], delay_bounds[1], modifier=setter(“delay”))) residual_fn = sysid.build_residual_fn(models_sequences=[sequences]) opt, result, evals = run_optimize(params, residual_fn) residuals, _, _ = residual_fn(result.x, opt) _, half_widths = sysid.calculate_intervals(residuals, result.jac) return {k: float(opt[k].value[0]) for k in opt.keys()}, half_widths, evals

def regression_sg(q, torque) -> np.ndarray: “"”The Lesson 18.1 estimate: Savitzky-Golay derivatives and least squares on the inverse dynamics.””” qf, qd, qdd = (savgol_filter(q, 51, 3, deriv=k, delta=H) for k in range(3)) sl = slice(30, len(q) - 30) x = np.c_[L * L * qdd + G * L * np.sin(qf), qd, np.sign(qd)][sl] return np.linalg.lstsq(x, (torque - I_C * qdd)[sl], rcond=None)[0]

def row(label: str, est: dict, hw=None) -> str: cells = [f”{est[k]:.5f}” + (f” +- {w:.5f}” if hw is not None else “”) for k, w in zip(est, hw if hw is not None else [None] * len(est))] return f” {label:<46}” + ““.join(f”{c:>21}” for c in cells)

def rms_mrad(a, b) -> float: return 1000 * float(np.sqrt(np.mean((a - b) ** 2)))

BALL = “””

””” BALL_TRUE = (0.015, 0.25)

def ball_spec(tc: float, dr: float, z0: float = 0.5) -> mujoco.MjSpec: spec = mujoco.MjSpec.from_string(BALL) for g in (“floor”, “ball”): spec.geom(g).solref = [tc, dr] spec.body(“ball”).pos = [0, 0, z0] return spec

def drop(tc: float, dr: float, z0: float = 0.5, seconds: float = 1.2) -> np.ndarray: “"”Height of the ball’s centre (m) at t = k h, before step k.””” model = ball_spec(tc, dr, z0).compile() data = mujoco.MjData(model) z = np.empty(round(seconds / H)) for k in range(len(z)): z[k] = z0 + data.qpos[0] mujoco.mj_step(model, data) return z

def apexes(z: np.ndarray, count: int) -> np.ndarray: “"”Heights of the first count bounce apexes of a smoothed height record (fewer if it stops bouncing).””” smooth = savgol_filter(z, 15, 2) peaks, _ = find_peaks(smooth, prominence=0.005) return smooth[peaks[:count]]

def ball_fit(start, z_meas, t) -> tuple[float, float]: “"”sysid fit of (time constant, damping ratio) to the raw height record from a starting point.””” spec = ball_spec(*start) model = spec.compile() data = mujoco.MjData(model) x0 = sysid.create_initial_state(model, data.qpos, data.qvel, data.act) sequences = sysid.ModelSequences( “ball”, spec, “drop”, x0, sysid.TimeSeries(t, np.zeros((len(t), 0))), sysid.TimeSeries.from_names(t + H, (z_meas - 0.5)[:, None], model, names=[“height”]))

def solref(i):
    def apply(s, p):
        for g in ("floor", "ball"):
            ref = list(s.geom(g).solref)
            ref[i] = float(p.value[0])
            s.geom(g).solref = ref
    return apply
params = sysid.ParameterDict()
params.add(sysid.Parameter("timeconst", start[0], 0.004, 0.1, modifier=solref(0)))   # MuJoCo clamps below 2 h
params.add(sysid.Parameter("dampratio", start[1], 0.01, 2.0, modifier=solref(1)))
opt, _, _ = run_optimize(params, sysid.build_residual_fn(models_sequences=[sequences]))
return float(opt["timeconst"].value[0]), float(opt["dampratio"].value[0])

if name == “main”: t = np.arange(5000) * H torque = multisine(t, 2.5, seed=1) q = simulate(pendulum_spec(TRUE), torque) + np.random.default_rng(0).normal(0.0, NOISE, len(t)) held_torque = multisine(t, 2.5, seed=7) q_held = simulate(pendulum_spec(TRUE), held_torque) + np.random.default_rng(2).normal(0.0, NOISE, len(t)) print(f”MuJoCo {mujoco.version}; pendulum record of 10 s at {1000 * H:.0f} ms, angle noise {1000 * NOISE} mrad”) print(f” {‘’:<46}” + ““.join(f”{k:>21}” for k in (“mass (kg)”, “damping”, “friction loss (N m)”))) print(row(“truth”, TRUE)) model = pendulum_spec(TRUE).compile() data = mujoco.MjData(model) data.ctrl[0] = 1.0 for _ in range(2): mujoco.mj_step(model, data) stepped = data.sensordata[0] mujoco.mj_forward(model, data) print(f” after two steps: time {data.time:.3f} s, qpos {data.qpos[0]:.3e}, angle sensor {stepped:.3e} “ f”(after mj_forward: {data.sensordata[0]:.3e})”) print(“(1) simulation error minimized by mujoco.sysid, starting from 1.0, 0.01, 0.01”) for label, shift in ((“samples stamped at their measurement time”, 0.0), (“samples stamped one step later, as rollout”, H)): est, hw, evals = fit_pendulum(t, torque, q, shift) print(row(label, est, hw) + f” ({evals} simulations)”) if shift: best = est print(row(“Lesson 18.1 regression (Savitzky-Golay)”, dict(zip(TRUE, regression_sg(q, torque))))) print(“ held-out excitation, angle RMS error: sysid estimate “ f”{rms_mrad(simulate(pendulum_spec(best), held_torque), q_held):.2f} mrad, truth “ f”{rms_mrad(simulate(pendulum_spec(TRUE), held_torque), q_held):.2f} mrad (the noise floor)”)

print("(2) the measured pendulum has a 10 ms actuator delay")
for interp, name in ((1, "linear interpolation"), (0, "the last sample held")):
    lags = []
    for delay in (0.0, 0.001, 0.002, 0.003, 0.004):
        spec = pendulum_spec(TRUE, delay)
        spec.actuator("torque").interp = interp
        model = spec.compile()
        data = mujoco.MjData(model)
        for k in range(3):
            data.ctrl[0] = k + 1.0                                     # a ramp: the command at step k is k + 1
            mujoco.mj_step(model, data)
        lags.append(f"{1000 * delay:.0f} ms -> {3.0 - data.actuator_force[0]:.1f}")
    print(f"  delay with {name}, steps by which the force lags a ramp: " + ", ".join(lags))
q_delay = simulate(pendulum_spec(TRUE, 0.010), torque) + np.random.default_rng(0).normal(0.0, NOISE, len(t))
q_delay_held = simulate(pendulum_spec(TRUE, 0.010), held_torque) + np.random.default_rng(2).normal(0.0, NOISE, len(t))
print(f"  {'':<46}" + "".join(f"{k:>21}" for k in ("mass (kg)", "damping", "friction loss (N m)", "delay (s)")))
fits = {"model without a delay": fit_pendulum(t, torque, q_delay, H),
        "delay as a parameter, bounds 0 to 30 ms": fit_pendulum(t, torque, q_delay, H, (0.0, 0.03)),
        "delay as a parameter, bounds 2 to 30 ms": fit_pendulum(t, torque, q_delay, H, (H, 0.03))}
for label, (est, hw, evals) in fits.items():
    print(row(label, est, hw) + f"   ({evals} simulations)")
for label, (est, _, _) in fits.items():
    spec = pendulum_spec(est, est.get("delay"))
    print(f"  held-out excitation, {label:<40} angle RMS error {rms_mrad(simulate(spec, held_torque), q_delay_held):7.2f} mrad")

print("(3) contact: a ball dropped from 0.5 m, height measured with 0.5 mm noise for 1.2 s; truth solref (0.015, 0.25)")
z_meas = drop(*BALL_TRUE) + np.random.default_rng(0).normal(0.0, 0.0005, 600)
t_ball = np.arange(600) * H
rms_mm = lambda tc, dr: 1000 * float(np.sqrt(np.mean((drop(tc, dr) - z_meas) ** 2)))
print("  highest point (mm) after the drop from 500 mm, damping ratio 0.25, time constant 1, 4, 5, 6, 8, 10, 15 ms: "
      + " ".join(f"{1000 * drop(tc, 0.25).max():.0f}" for tc in (0.001, 0.004, 0.005, 0.006, 0.008, 0.010, 0.015)))
tcs, drs = np.linspace(0.010, 0.020, 11), np.linspace(0.20, 0.30, 11)
print("  height RMS error (mm) along the time constant 0.010-0.020 s at the true damping ratio: "
      + " ".join(f"{rms_mm(a, 0.25):.1f}" for a in tcs))
print("  height RMS error (mm) along the damping ratio 0.20-0.30 at the true time constant:   "
      + " ".join(f"{rms_mm(0.015, b):.1f}" for b in drs))
starts = [(0.02, 1.0), (0.02, 0.5)] if FAST else [(0.02, 1.0), (0.02, 0.5), (0.02, 0.1), (0.005, 0.25), (0.04, 0.05)]
results = []
for start in starts:
    tc, dr = ball_fit(start, z_meas, t_ball)
    results.append((rms_mm(tc, dr), tc, dr))
    print(f"  sysid fit from {start}: time constant {tc:.4f} s, damping ratio {dr:.4f}, height RMS error "
          f"{results[-1][0]:.2f} mm, highest point {1000 * drop(tc, dr).max():.0f} mm")
n = 10 if FAST else 40
grid = [(rms_mm(a, b), a, b) for a in np.geomspace(0.004, 0.1, n) for b in np.geomspace(0.02, 2.0, n)]
g_best = min(grid)
print(f"  best of a {n} x {n} log grid: time constant {g_best[1]:.4f} s, damping ratio {g_best[2]:.4f}, "
      f"height RMS error {g_best[0]:.2f} mm; at the truth {rms_mm(*BALL_TRUE):.2f} mm")
a_meas = apexes(z_meas, 2)
print(f"  bounce apex heights measured: {', '.join(f'{1000 * a:.1f}' for a in a_meas)} mm")
z_held = drop(*BALL_TRUE, z0=0.3) + np.random.default_rng(1).normal(0.0, 0.0005, 600)
candidates = {"best raw-trajectory fit": min(results + [g_best])[1:]}
for tc in (0.01, 0.02, 0.03):
    def cost(dr, tc=tc):
        a = apexes(drop(tc, dr), 2)
        return float(np.sum((a - a_meas) ** 2)) if len(a) == 2 else 1.0     # no second bounce: far off
    scan = np.geomspace(0.02, 1.5, 60)
    i = int(np.argmin([cost(d) for d in scan]))
    dr = minimize_scalar(cost, bounds=(scan[max(i - 1, 0)], scan[min(i + 1, 59)]), method="bounded",
                         options={"xatol": 1e-5}).x
    candidates[f"apex fit, time constant fixed at {tc}"] = (tc, float(dr))
print("  held-out drop from 0.3 m (apex heights measured: "
      + ", ".join(f"{1000 * a:.1f}" for a in apexes(z_held, 2)) + " mm):")
for label, (tc, dr) in candidates.items():
    z = drop(tc, dr, z0=0.3)
    print(f"    {label:<40} ({tc:.4f}, {dr:.4f}): apexes "
          + ", ".join(f"{1000 * a:.1f}" for a in apexes(z, 2)) + f" mm, height RMS error {1000 * np.sqrt(np.mean((z - z_held) ** 2)):.2f} mm") ``` Output:

MuJoCo 3.14.0; pendulum record of 10 s at 2 ms, angle noise 0.5 mrad
                                                            mass (kg)              damping  friction loss (N m)
  truth                                                       1.30000              0.08000              0.05000
  after two steps: time 0.004 s, qpos 3.504e-05, angle sensor 1.168e-05 (after mj_forward: 3.504e-05)
(1) simulation error minimized by mujoco.sysid, starting from 1.0, 0.01, 0.01
  samples stamped at their measurement time        1.30449 +- 0.00021   0.08276 +- 0.00027   0.04536 +- 0.00030   (39 simulations)
  samples stamped one step later, as rollout       1.30001 +- 0.00008   0.07995 +- 0.00010   0.05004 +- 0.00011   (35 simulations)
  Lesson 18.1 regression (Savitzky-Golay)                     1.32655              0.08676              0.04166
  held-out excitation, angle RMS error: sysid estimate 0.50 mrad, truth 0.50 mrad (the noise floor)
(2) the measured pendulum has a 10 ms actuator delay
  delay with linear interpolation, steps by which the force lags a ramp: 0 ms -> 0.0, 1 ms -> 1.0, 2 ms -> 1.0, 3 ms -> 1.5, 4 ms -> 2.0
  delay with the last sample held, steps by which the force lags a ramp: 0 ms -> 0.0, 1 ms -> 1.0, 2 ms -> 1.0, 3 ms -> 2.0, 4 ms -> 2.0
                                                            mass (kg)              damping  friction loss (N m)            delay (s)
  model without a delay                            1.27683 +- 0.00091   0.06585 +- 0.00102   0.07384 +- 0.00110   (178 simulations)
  delay as a parameter, bounds 0 to 30 ms          1.28112 +- 0.00082   0.06753 +- 0.00086   0.07050 +- 0.00094       0.00000 +- inf   (176 simulations)
  delay as a parameter, bounds 2 to 30 ms          1.29999 +- 0.00009   0.07993 +- 0.00010   0.05007 +- 0.00012   0.00999 +- 0.00002   (44 simulations)
  held-out excitation, model without a delay                    angle RMS error   11.32 mrad
  held-out excitation, delay as a parameter, bounds 0 to 30 ms  angle RMS error    9.16 mrad
  held-out excitation, delay as a parameter, bounds 2 to 30 ms  angle RMS error    0.50 mrad
(3) contact: a ball dropped from 0.5 m, height measured with 0.5 mm noise for 1.2 s; truth solref (0.015, 0.25)
  highest point (mm) after the drop from 500 mm, damping ratio 0.25, time constant 1, 4, 5, 6, 8, 10, 15 ms: 4001 4001 1324 500 500 500 500
  height RMS error (mm) along the time constant 0.010-0.020 s at the true damping ratio: 8.2 18.0 3.6 4.4 6.1 0.5 6.9 4.1 4.3 2.9 3.2
  height RMS error (mm) along the damping ratio 0.20-0.30 at the true time constant:   23.8 23.5 21.6 14.6 9.3 0.5 3.2 4.5 7.3 10.6 14.0
  sysid fit from (0.02, 1.0): time constant 0.0040 s, damping ratio 0.6190, height RMS error 13.72 mm, highest point 500 mm
  sysid fit from (0.02, 0.5): time constant 0.0103 s, damping ratio 0.2122, height RMS error 29.14 mm, highest point 500 mm
  sysid fit from (0.02, 0.1): time constant 0.0193 s, damping ratio 0.1035, height RMS error 88.23 mm, highest point 500 mm
  sysid fit from (0.005, 0.25): time constant 0.0055 s, damping ratio 0.2569, height RMS error 333.12 mm, highest point 666 mm
  sysid fit from (0.04, 0.05): time constant 0.0585 s, damping ratio 0.0231, height RMS error 104.10 mm, highest point 500 mm
  best of a 40 x 40 log grid: time constant 0.0127 s, damping ratio 0.2687, height RMS error 1.28 mm; at the truth 0.50 mm
  bounce apex heights measured: 152.5, 53.8 mm
  held-out drop from 0.3 m (apex heights measured: 82.8, 45.2 mm):
    best raw-trajectory fit                  (0.0127, 0.2687): apexes 72.0, 39.9 mm, height RMS error 5.75 mm
    apex fit, time constant fixed at 0.01    (0.0100, 0.2473): apexes 77.8, 37.7 mm, height RMS error 4.06 mm
    apex fit, time constant fixed at 0.02    (0.0200, 0.2514): apexes 84.4, 40.8 mm, height RMS error 2.18 mm
    apex fit, time constant fixed at 0.03    (0.0300, 0.2467): apexes 88.1, 44.5 mm, height RMS error 5.72 mm

Debugging

mujoco.sysid has no attribute optimize. The extras are missing: pip install "mujoco[sysid]==3.14.0".

The fit is good but every estimate is off by a few percent with tight intervals. A timing convention: compare a fit on simulated data stamped both ways.

A parameter ends at a bound with an infinite interval. The cost is flat there: a delay shorter than one step, a held (interp="zoh") delay, a time constant under two steps. Move the bound or change how the parameter acts.

Different starts give different answers. The cost has many local minima. Plot it along each parameter, search coarsely first, or fit features that vary smoothly.

The identified contact bounces higher than it fell. The optimizer found a numerically unstable region. Bound the parameters away from it and check energy on validation.

Exercise

Identify the pendulum from the first 2 s of the record and from the full record with 2 mrad of angle noise. For each, report the estimates, the intervals and the held-out error, and compare with the 18.1 regression under the same conditions.

Challenge

Make the ball’s raw-trajectory fit work. Split the record at each apex, start each segment from the measured height and a filtered velocity, and fit (time constant, damping ratio) to the sum of the segment errors (multiple shooting). Plot the cost along each parameter and report whether the five starts now converge, and to what.

Research connection

[!research] Identification for sim-to-real Tan et al. (arXiv:1804.10332) improved their quadruped simulator with system identification, an actuator model and simulated latency, and also randomized the physical environment during training. Hwangbo et al. (Science Robotics 4(26), 2019, eaau5872; arXiv:1901.08652) learned an “actuator net” on the physical robot by self-supervised learning and used it in the simulation loop for each of ANYmal’s 12 joints. Read both for which parameters they identified, from what data, and how they validated, and compare with this lesson’s findings (inference): unmodeled dynamics such as a delay hide inside identified parameters with tight intervals, and contact parameters fitted to trajectories meet a jagged cost. Lesson 19.1 measures what a delay does to a controller.

{"id": "18.2-check", "title": "Knowledge check", "questions": [
  {"kind": "mcq", "q": "Why does simulation error recover the pendulum's damping and friction loss when the 18.1 regression did not?",
   "options": ["It uses more data", "It needs no differentiated accelerations and uses MuJoCo's own friction-loss constraint", "Its optimizer is better", "It ignores noise"],
   "answer": 1,
   "explain": "<p>The 18.1 biases came from derivatives and from modeling friction loss as f sign(qd). On data generated by the same model class (an inverse crime), simulation error recovers the truth.</p>"},
  {"kind": "mcq", "q": "A fit on data stamped at the measurement time is off by 9% in friction loss with a 0.6% interval. What is the likeliest cause?",
   "options": ["Too much noise", "A one-step timing offset between measured samples and the simulator's sensor stamps", "Too few iterations", "The friction is nonlinear"],
   "answer": 1,
   "explain": "<p>mujoco.rollout pairs the state after step k with sensors computed before it; stamp measurements one step later.</p>"},
  {"kind": "numeric", "q": "With linear interpolation and a 2 ms timestep, how many steps does the force lag behind the command for a 1 ms delay?",
   "answer": 1, "tol": 0.01,
   "explain": "<p>Any delay in (0, h] acts as one full step, so the cost is flat there; bound a fitted delay below by one step.</p>"},
  {"kind": "mcq", "q": "Bounce apex heights determine the damping ratio to within 1.3% but not the contact time constant. Why?",
   "options": ["The time constant has no effect", "Apex heights measure restitution, which here depends mainly on the damping ratio", "The optimizer failed", "The noise is too large"],
   "answer": 1,
   "explain": "<p>Fixing the time constant at 0.01, 0.02 or 0.03 s gives damping ratios 0.2473, 0.2514, 0.2467; the time constant needs data sensitive to contact duration.</p>"}
]}

Next

Level 19 puts the pieces together for a real robot: what a delay, a mismatched actuator and noise do to a controller tuned in simulation, and how to design a transfer experiment that fails informatively.