System identification estimates a model’s parameters from measured motion. For a rigid-body system the inverse dynamics are linear in the inertial parameters and in viscous and Coulomb friction coefficients, so given positions, velocities, accelerations and torques, the parameters follow from ordinary least squares. The difficulty is everything around that one line of algebra: velocities and accelerations are not measured but differentiated from noisy positions; the friction model in the regression may not be the friction in the system; and the data must actually excite each parameter. This lesson identifies the mass, damping and friction loss of a pendulum whose parameters are hidden from the estimator, and measures what goes wrong at each of those three points.
[!established] The model being identified The pendulum of
pendulum.xmlis a point mass $m$ at $L$ = 0.5 m with a small rotor inertia $I_c = 10^{-4}$ kg m², a joint damping $b$ and a joint friction loss $f$. With torque $\tau$, $(mL^2 + I_c)\ddot q + b\dot q + mgL\sin q + \text{friction loss} = \tau$. MuJoCo implements friction loss as a constraint with a bounded force $f$ (its documentation calls it dry friction), softened like every MuJoCo constraint, so near zero velocity it is not exactly $f\,\text{sign}(\dot q)$.
New hidden pendulum draws a mass, damping and friction loss you cannot see; Record 10 s drives it with a multisine torque and measures its angle with noise; Fit estimates the parameters and validates them on a different excitation, plotting the measured angle against the estimated model’s prediction.
```lab sysid {“dock”: true, “title”: “Identify a pendulum”}
## Mathematics
### A linear regression
Move the known terms to the left and the unknowns to the right:
$$\tau - I_c\ddot q = m\,(L^2\ddot q + gL\sin q) + b\,\dot q + f\,\text{sign}(\dot q) = \phi(q, \dot q, \ddot q)^\top \theta,\qquad \theta = (m, b, f).$$
Stacking 5 000 samples gives $Y = \Phi\theta + \epsilon$ and $\hat\theta = (\Phi^\top\Phi)^{-1}\Phi^\top Y$. The hidden values are $m$ = 1.3 kg, $b$ = 0.08 N m s/rad, $f$ = 0.05 N m; the torque is a sum of five sines (peak 2.5 N m) for 10 s; the angle is measured with Gaussian noise of 0.5 mrad.
### Derivatives decide the answer
| Derivatives | $m$ | $b$ | $f$ |
|---|---|---|---|
| finite differences of the noisy angle | 0.0007 | 0.1006 | 0.0514 |
| Savitzky-Golay (cubic fits over 0.1 s) | 1.3266 | 0.0868 | 0.0417 |
| finite differences of the noise-free angle | 1.2999 | 0.0835 | 0.0460 |
Differencing twice at a 2 ms step multiplies the angle noise by about $\sqrt 6/h^2 \approx 6\times 10^5$ s⁻², so 0.5 mrad of noise becomes roughly 300 rad/s² of noise in $\ddot q$, far larger than the true accelerations. The regressor of $m$ is then mostly noise, and least squares, which assumes the regressors are exact, shrinks its coefficient toward zero: the errors-in-variables bias. The mass estimate is 0.0007 kg. A Savitzky-Golay filter fits a cubic to each 0.1 s window and differentiates the fit, which recovers the mass within 2%.
### Uncertainty that is honest, and bias that is not covered
The least-squares standard errors assume independent residuals. Filtered derivatives make neighbouring residuals strongly correlated (autocorrelation 0.906 at one sample), and a block bootstrap that resamples contiguous stretches of the record gives standard errors about twice as large (0.0028 against 0.0016 kg for the mass, 0.0038 against 0.0017 for the damping). Even the bootstrap interval misses the truth for damping and friction loss by about two standard errors. That is not noise but model error: refitting without the 9% of samples slower than 0.1 rad/s gives 1.3238, 0.0797 and 0.0501, essentially the hidden values. Near zero velocity MuJoCo's friction loss is a soft constraint, not $f\,\text{sign}(\dot q)$, and those samples pulled damping and friction loss toward each other. No standard error accounts for a wrong model; only a check of the residuals and of the assumptions does.
### Excitation
With a tenth of the torque (0.25 N m), the pendulum never moves faster than 0.10 rad/s, and the estimate is $m$ = 0.61, $b$ = 0.75, $f$ = 0.017: wrong by factors of 2 to 9, with standard errors that claim a precision of a few percent. The condition number of the column-normalized regressor is 4 for both excitations, so it did not warn. Every sample is in the regime where the friction model is wrong and the acceleration signal is small next to its noise. Excitation design is part of the identification: the trajectory must move every parameter's regressor well above its noise, and away from regimes the model does not describe.
### Validation
A held-out trajectory, driven by a different multisine and measured with noise, is simulated open loop with each estimate:
| Parameters | Angle RMS error over 10 s | Over the first second |
|---|---|---|
| Savitzky-Golay estimate | 2.57 mrad | 1.12 mrad |
| the same, without slow samples | 2.53 mrad | 1.11 mrad |
| finite-difference estimate | 907 mrad | 949 mrad |
| a nominal guess (1.0 kg, no damping or friction) | 129 mrad | 24.4 mrad |
The two good estimates differ by 8% in damping and 17% in friction loss and predict the held-out motion equally well: this trajectory constrains a combination of the two more than each one. Validation can reject a bad model; agreement on one trajectory does not prove every parameter right.
> [!recommendation] An identification checklist
> Filter before differentiating, and check the filter on noise-free simulated data. Excite every parameter, and check the speeds and accelerations the data actually reach. Report block-bootstrap intervals, then look at the residuals against velocity for model error. Validate on a trajectory not used for fitting, and for parameters that trade off, test on a trajectory that separates them.
## Implementation
```io
INPUT: pendulum.xml with hidden mass, damping and friction loss; a noisy angle record
PROCESS: differentiate by finite differences and by Savitzky-Golay; fit by least squares; estimate uncertainty with a block bootstrap; vary the excitation; validate open loop on a held-out trajectory
OUTPUT: printed tables
```python file=examples/l18_1_sysid_least_squares.py “"”Lesson 18.1: identifying a pendulum’s mass, damping and friction loss by least squares (Level 18 checkpoint).
INPUT pendulum.xml with hidden parameters: mass 1.3 kg, joint damping 0.08 N m s/rad, friction loss 0.05 N m; the angle measured with Gaussian noise (0.5 mrad), the torque known (it is the command) PROCESS the inverse dynamics tau = m (L^2 qdd + g L sin q) + I_c qdd + b qd + f sign(qd) are linear in (m, b, f) given q, qd, qdd. From 10 s of multisine excitation: (1) qd and qdd by finite differences of the noisy angle, and by a Savitzky-Golay filter; least-squares estimates from each; (2) uncertainty: the ordinary least-squares standard errors against a block bootstrap; (3) identifiability: a weak excitation (a tenth of the amplitude), the condition number of the regressor, how much of the record is slow, and a fit that leaves out samples slower than 0.1 rad/s, where the friction model is wrong; (4) validation: a held-out trajectory with a different excitation, simulated open loop with the true and with the estimated parameters OUTPUT printed tables
Run: python examples/l18_1_sysid_least_squares.py “””
import mujoco import numpy as np from scipy.signal import savgol_filter
from mjcourse import model_path
TRUE = {“mass”: 1.3, “damping”: 0.08, “frictionloss”: 0.05} L, I_C, G = 0.5, 1e-4, 9.81 NOISE = 0.0005 # angle noise (rad)
def pendulum(params: dict) -> tuple[mujoco.MjModel, mujoco.MjData]: model = mujoco.MjModel.from_xml_path(str(model_path(“pendulum”))) body, joint = model.body(“pole”).id, model.joint(“hinge”).id model.body_mass[body] = params[“mass”] model.dof_damping[joint] = params[“damping”] model.dof_frictionloss[joint] = params[“frictionloss”] data = mujoco.MjData(model) mujoco.mj_setConst(model, data) # derived constants after editing the mass (Lesson 17.1) mujoco.mj_resetData(model, data) return model, data
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(params: dict, torque: np.ndarray) -> np.ndarray: model, data = pendulum(params) 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 regressor(q, qd, qdd) -> np.ndarray: “"”Columns multiply (m, b, f); the known rotor inertia term goes to the left-hand side.””” return np.c_[L * L * qdd + G * L * np.sin(q), qd, np.sign(qd)]
def derivatives(q: np.ndarray, h: float, method: str) -> tuple[np.ndarray, np.ndarray, np.ndarray]: if method == “finite differences”: qd = np.gradient(q, h) return q, qd, np.gradient(qd, h) window = 51 # 0.1 s at 2 ms, cubic fits return (savgol_filter(q, window, 3), savgol_filter(q, window, 3, deriv=1, delta=h), savgol_filter(q, window, 3, deriv=2, delta=h))
def fit(q, qd, qdd, tau, trim: int = 30, min_speed: float = 0.0) -> tuple[np.ndarray, np.ndarray, np.ndarray]: “"”Least squares on the samples away from the ends (and, optionally, faster than min_speed); returns (estimate, standard errors, residuals).””” sl = np.arange(trim, len(q) - trim) sl = sl[np.abs(qd[sl]) >= min_speed] x, y = regressor(q[sl], qd[sl], qdd[sl]), tau[sl] - I_C * qdd[sl] theta, *_ = np.linalg.lstsq(x, y, rcond=None) resid = y - x @ theta cov = resid.var(ddof=3) * np.linalg.inv(x.T @ x) return theta, np.sqrt(np.diag(cov)), resid
def block_bootstrap(q, qd, qdd, tau, blocks: int = 20, n: int = 500, seed: int = 0) -> np.ndarray: “"”Refit on resampled contiguous blocks of the record: errors that are correlated in time stay together.””” rng = np.random.default_rng(seed) idx = np.array_split(np.arange(30, len(q) - 30), blocks) estimates = [] for _ in range(n): pick = np.concatenate([idx[i] for i in rng.integers(0, blocks, blocks)]) x, y = regressor(q[pick], qd[pick], qdd[pick]), tau[pick] - I_C * qdd[pick] estimates.append(np.linalg.lstsq(x, y, rcond=None)[0]) return np.std(estimates, axis=0)
def row(name: str, theta, se=None) -> str: cells = [f”{v:.4f}” + (f” +- {s:.4f}” if se is not None else “”) for v, s in zip(theta, se if se is not None else theta)] return f” {name:<44}” + ““.join(f”{c:>22}” for c in cells)
if name == “main”: model, _ = pendulum(TRUE) h = model.opt.timestep t = np.arange(round(10.0 / h)) * h rng = np.random.default_rng(0) torque = multisine(t, 2.5, seed=1) q_true = simulate(TRUE, torque) q = q_true + rng.normal(0.0, NOISE, len(q_true)) print(f”MuJoCo {mujoco.version}; 10 s at {1000 * h:.0f} ms, angle noise {1000 * NOISE} mrad, “ f”angle range {q_true.min():.2f} to {q_true.max():.2f} rad, peak speed {np.abs(np.gradient(q_true, h)).max():.2f} rad/s”) print(f” {‘’:<44}” + ““.join(f”{k:>22}” for k in (“mass (kg)”, “damping (N m s/rad)”, “friction loss (N m)”))) print(row(“truth”, list(TRUE.values()))) print(“(1) derivatives and estimates”) for method in (“finite differences”, “Savitzky-Golay”): theta, se, _ = fit(derivatives(q, h, method), torque) print(row(method, theta)) theta, se, resid = fit(derivatives(q_true, h, “finite differences”), torque) print(row(“finite differences of the noise-free angle”, theta)) print(“(2) uncertainty of the Savitzky-Golay estimate”) qs = derivatives(q, h, “Savitzky-Golay”) theta, se, resid = fit(qs, torque) print(row(“ordinary least-squares standard errors”, theta, se)) print(row(“block bootstrap standard errors”, theta, block_bootstrap(qs, torque))) lag1 = np.corrcoef(resid[:-1], resid[1:])[0, 1] print(f” residual autocorrelation at one sample: {lag1:.3f}”) print(“(3) excitation and identifiability”) for label, amplitude in ((“full excitation (2.5 N m)”, 2.5), (“weak excitation (0.25 N m)”, 0.25)): tq = multisine(t, amplitude, seed=1) qq = simulate(TRUE, tq) + np.random.default_rng(1).normal(0.0, NOISE, len(t)) dq = derivatives(qq, h, “Savitzky-Golay”) x = regressor((d[30:-30] for d in dq)) theta, _, _ = fit(dq, tq) print(row(f”{label}, cond {np.linalg.cond(x / np.abs(x).max(axis=0)):.0f}”, theta, block_bootstrap(dq, tq))) speed = np.abs(dq[1][30:-30]) fast = np.sum(speed >= 0.1) refit = (“, “.join(f”{v:.4f}” for v in fit(dq, tq, min_speed=0.1)[0]) if fast >= 100 else f”only {fast} samples, too few to fit”) print(f” peak speed {speed.max():.2f} rad/s; {np.mean(speed < 0.1):.0%} of samples slower than 0.1 rad/s; “ f”fit without them: {refit}”) print(“(4) validation on a held-out trajectory (different multisine), open-loop simulation”) estimate = dict(zip(TRUE, fit(derivatives(q, h, “Savitzky-Golay”), torque)[0])) held_torque = multisine(t, 2.5, seed=7) q_held = simulate(TRUE, held_torque) + np.random.default_rng(2).normal(0.0, NOISE, len(t)) filtered = dict(zip(TRUE, fit(derivatives(q, h, “Savitzky-Golay”), torque, min_speed=0.1)[0])) for label, params in ((“estimated parameters”, estimate), (“estimated without samples slower than 0.1 rad/s”, filtered), (“finite-difference estimate”, dict(zip(TRUE, fit(*derivatives(q, h, “finite differences”), torque)[0]))), (“nominal guess (1.0 kg, no damping, no friction)”, {“mass”: 1.0, “damping”: 0.0, “frictionloss”: 0.0})): q_sim = simulate(params, held_torque) err = q_sim - q_held print(f” {label:<48} angle RMS error {1000 * np.sqrt(np.mean(err2)):7.2f} mrad, “ f”after 1 s {1000 * np.sqrt(np.mean(err[:round(1 / h)]2)):6.2f} mrad”)
Output:
```text
MuJoCo 3.14.0; 10 s at 2 ms, angle noise 0.5 mrad, angle range -0.46 to 0.44 rad, peak speed 1.94 rad/s
mass (kg) damping (N m s/rad) friction loss (N m)
truth 1.3000 0.0800 0.0500
(1) derivatives and estimates
finite differences 0.0007 0.1006 0.0514
Savitzky-Golay 1.3266 0.0868 0.0417
finite differences of the noise-free angle 1.2999 0.0835 0.0460
(2) uncertainty of the Savitzky-Golay estimate
ordinary least-squares standard errors 1.3266 +- 0.0016 0.0868 +- 0.0017 0.0417 +- 0.0016
block bootstrap standard errors 1.3266 +- 0.0028 0.0868 +- 0.0038 0.0417 +- 0.0041
residual autocorrelation at one sample: 0.906
(3) excitation and identifiability
full excitation (2.5 N m), cond 4 1.3247 +- 0.0028 0.0823 +- 0.0034 0.0470 +- 0.0037
peak speed 1.94 rad/s; 9% of samples slower than 0.1 rad/s; fit without them: 1.3238, 0.0797, 0.0501
weak excitation (0.25 N m), cond 4 0.6073 +- 0.0485 0.7471 +- 0.1055 0.0171 +- 0.0038
peak speed 0.10 rad/s; 100% of samples slower than 0.1 rad/s; fit without them: only 0 samples, too few to fit
(4) validation on a held-out trajectory (different multisine), open-loop simulation
estimated parameters angle RMS error 2.57 mrad, after 1 s 1.12 mrad
estimated without samples slower than 0.1 rad/s angle RMS error 2.53 mrad, after 1 s 1.11 mrad
finite-difference estimate angle RMS error 906.74 mrad, after 1 s 949.33 mrad
nominal guess (1.0 kg, no damping, no friction) angle RMS error 129.31 mrad, after 1 s 24.42 mrad
An inertial parameter comes out near zero. Noisy accelerations in its regressor. Filter, or lower the noise, or use a derivative-free formulation (Lesson 18.2).
Intervals are tight and the estimate is still wrong. Correlated residuals (use a block bootstrap) or a wrong model (look at residuals against velocity, position and time).
Damping and friction come out traded against each other. They are hard to separate from motion that rarely changes speed. Excite a wide range of speeds, leave out the regime the friction model does not describe, or fix one from a separate experiment.
The estimate validates well on one trajectory and fails on another. The first did not excite what the second needs. Validate on trajectories that differ from the identification data in speed and amplitude.
Identify the three parameters with 2 mrad of angle noise. Find the shortest Savitzky-Golay window that keeps the mass within 5%, and report what that window costs in bias of the damping estimate.
Add the rotor inertia $I_c$ as a fourth unknown and identify it together with the other three. Show from the regressor why it is hard to separate from the mass with this pendulum, and design an excitation, within the actuator’s 3 N m limit, that separates them better.
[!research] Identification as a sim-to-real tool Inverse-dynamics least squares is the classic first step in identifying robot arms, where the inertial parameters enter linearly through a regressor matrix and friction is added per joint; excitation trajectories are often optimized to make that regressor well conditioned. The lesson’s results carry over directly (inference): filtering decides whether inertial parameters survive, friction models fail near zero velocity, and residual correlation makes naive intervals too narrow. When parameters do not enter linearly, such as contact stiffness or friction coefficients between objects, identification has to simulate, which is Lesson 18.2.
{"id": "18.1-check", "title": "Knowledge check", "questions": [
{"kind": "numeric", "q": "Angle noise of 0.5 mrad, second differences at a 2 ms step. Roughly how large is the noise in the acceleration (rad/s^2)? Use sqrt(6) sigma / h^2.",
"answer": 306, "tol": 5, "unit": "rad/s^2",
"explain": "<p>$\\sqrt 6 \\times 0.0005 / 0.002^2$ = 306 rad/s², far above the true accelerations.</p>"},
{"kind": "mcq", "q": "Why does noise in a regressor bias a least-squares coefficient toward zero?",
"options": ["Least squares assumes the regressors are exact; noise in them inflates their variance without adding signal", "The noise averages out", "The torque is noisy", "The model is nonlinear"],
"answer": 0,
"explain": "<p>Errors-in-variables: the fitted slope is the true one times the signal's share of the regressor's variance.</p>"},
{"kind": "mcq", "q": "The bootstrap interval misses the true damping by two standard errors, and refitting without slow samples hits it. What was wrong?",
"options": ["Too few samples", "The friction model: near zero velocity MuJoCo's friction loss is a soft constraint, not f sign(qd)", "The filter window", "The random seed"],
"answer": 1,
"explain": "<p>Output (3): without the 9% of samples slower than 0.1 rad/s the estimate is 1.3238, 0.0797, 0.0501.</p>"},
{"kind": "mcq", "q": "Two parameter sets differ by 8% in damping and 17% in friction loss but predict a held-out trajectory equally well. What does that mean?",
"options": ["Both are exactly right", "That trajectory constrains a combination of the two more than each one", "The validation is broken", "Friction does not matter"],
"answer": 1,
"explain": "<p>Validation rejects bad models; it does not prove each parameter right.</p>"}
]}
Lesson 18.2 identifies parameters that do not enter the dynamics linearly, by fitting simulations to measurements.