A joint-space controller decides each joint’s torque from where the joint is and where it should be. The simplest useful one is PD: push towards the target in proportion to the error, and resist velocity.
\[\tau = k_p\,(q^\star - q) - k_d\,\dot q .\]Two numbers, and nearly every question about them has a quantitative answer. How much does it overshoot? The damping ratio says. Why does it stop short of the target under gravity? The stiffness says by how much. How do you remove that error? With a model of gravity, or with integral action, which brings its own failure at the torque limit. And should you write this controller yourself, or use MuJoCo’s servo actuators? It depends on the integrator, for a reason Lesson 5.2 already measured.
This lesson answers each with the pendulum (1 kg at 0.5 m, inertia 0.2501 kg m² about the hinge, a motor limited to ±3 N m), where every prediction is a short formula.
[!established] MuJoCo’s servo actuators
<position kp kv>applies $k_p(u - q) - k_v\dot q$ with $u$ the control,<velocity kv>applies $k_v(u - \dot q)$, and since 3.12<pid kp kv ki imax>adds integral action, with the integral of the position error stored in the actuator’s activation (act) and clamped atimax(anti-windup). Underimplicitfastandimplicit, MuJoCo integrates the velocity-dependent part of these actuators implicitly; underEulerit does not.
```lab simlab {“dock”: true, “title”: “Tuning one joint”, “model”: “pendulum”, “height”: 230, “camera”: {“azimuth”: -90, “elevation”: 5, “distance”: 1.9, “target”: [0, 0, 0.75]}, “trail”: “tip”, “setup”: “ctx.ie = 0”, “controls”: [ {“label”: “target angle”, “unit”: “rad”, “param”: “qd”, “min”: -1.2, “max”: 1.2, “step”: 0.01, “value”: 0.5, “digits”: 2}, {“label”: “kp”, “unit”: “N m/rad”, “param”: “kp”, “min”: 1, “max”: 60, “step”: 0.5, “value”: 20, “digits”: 1}, {“label”: “damping ratio ζ”, “param”: “zeta”, “min”: 0, “max”: 1.5, “step”: 0.05, “value”: 0.7, “digits”: 2}, {“label”: “ki”, “unit”: “N m/(rad s)”, “param”: “ki”, “min”: 0, “max”: 60, “step”: 1, “value”: 0, “digits”: 0}, {“type”: “select”, “label”: “integral clamp (anti-windup)”, “param”: “clamp”, “options”: [[“off”, 0], [“±0.07 rad s”, 0.07]], “value”: 0}, {“type”: “select”, “label”: “gravity”, “set”: “model.opt.gravity[2]”, “options”: [[“on”, -9.81], [“off”, 0]], “value”: -9.81}, {“type”: “select”, “label”: “gravity compensation”, “param”: “gc”, “options”: [[“off”, 0], [“on: m g L sin q”, 1]], “value”: 0}, {“type”: “button”, “label”: “Restart from hanging”, “run”: “mj.mj_resetData(model, data); ctx.ie = 0”}], “controller”: “const h = model.opt.timestep, I = 0.2501, g = -model.opt.gravity[2]; const kd = 2 * ctx.zeta * Math.sqrt(ctx.kp * I); const e = ctx.qd - data.qpos[0]; ctx.ie += e * h; if (ctx.clamp > 0) ctx.ie = Math.min(Math.max(ctx.ie, -ctx.clamp), ctx.clamp); data.ctrl[0] = ctx.kp * e - kd * data.qvel[0] + ctx.ki * ctx.ie + ctx.gc * g * 0.5 * Math.sin(data.qpos[0])”, “readouts”: [ {“label”: “error (mrad)”, “expr”: “1000 * (ctx.qd - data.qpos[0])”, “digits”: 1}, {“label”: “motor torque (N m)”, “expr”: “data.actuator_force[0]”, “digits”: 2}, {“label”: “overshoot predicted, no gravity (%)”, “expr”: “ctx.zeta >= 1 ? 0 : 100 * Math.exp(-Math.PI * ctx.zeta / Math.sqrt(1 - ctx.zeta * ctx.zeta))”, “digits”: 1}, {“label”: “PD steady-state error predicted (mrad)”, “expr”: “1000 * 4.905 * Math.sin(ctx.qd) / (ctx.kp + 4.905 * Math.cos(ctx.qd)) * (-model.opt.gravity[2] / 9.81)”, “digits”: 1}], “plots”: [{“ylabel”: “angle (rad)”, “window”: 6, “traces”: [ {“label”: “angle (rad)”, “expr”: “data.qpos[0]”}, {“label”: “target (rad)”, “expr”: “ctx.qd”, “dash”: true}]}], “note”: “The controller runs before every step: kd is set from the damping ratio as 2ζ√(kp I). The steady-state prediction is for PD alone (no compensation, no integral). The motor saturates at ±3 N m, so large steps or large targets show the limit too.”}
## Mathematics
### PD as a mass-spring-damper
Without gravity, the joint obeys $I\ddot q = k_p(q^\star - q) - k_d\dot q$: a mass on a spring with a damper. Its natural frequency and damping ratio are
$$\omega_n = \sqrt{k_p / I}, \qquad \zeta = \frac{k_d}{2\sqrt{k_p I}},$$
and for $\zeta < 1$ a step response overshoots by $e^{-\pi\zeta/\sqrt{1-\zeta^2}}$ of the step, peaking at $t_p = \pi / (\omega_n\sqrt{1-\zeta^2})$. Output (1) checks both for five damping ratios: the measured overshoot is within 0.2 percentage points and the peak time within 2 ms of the formulas. Choose $k_p$ for speed (bandwidth $\approx \omega_n$) and $k_d$ for the overshoot you accept; $\zeta \approx 0.7$ is the usual compromise.
### The error gravity leaves
With gravity, PD settles where its spring balances the weight: $k_p(q^\star - q) = mgL\sin q$. Linearizing around the target gives
$$e_{ss} \approx \frac{mgL\sin q^\star}{k_p + mgL\cos q^\star},$$
which shrinks with $k_p$ but never vanishes. Output (2) measures 162, 96 and 43 mrad for $k_p$ = 10, 20 and 50 against predictions of 164, 97 and 43.
Two remedies. **Gravity compensation** adds the model's gravity torque, $mgL\sin q$ here or `qfrc_bias` on a general robot (Lesson 7.1): exact when the model is exact, and it removes the error instantly (output 2: zero to the printed precision). **Integral action** adds $k_i\int e\,dt$: it needs no model and removes any constant disturbance, but only gradually, and it has a failure mode of its own.
### Windup
While the motor is saturated, the error stays large and the integral keeps growing. When the joint finally reaches the target, the stored integral pushes it past. Output (3): with $k_i$ = 40 and a 0.55 rad target, the arm overshoots by 87 mrad. **Anti-windup** clamps the integral. The clamp must exceed the integral the hold actually needs ($mgL\sin q^\star / k_i$ = 0.064 rad s here): at 0.07 the overshoot falls to 9 mrad and settling from 1.6 s to 0.5 s; at 0.05 the controller can never supply the holding torque and stops 23 mrad short. MuJoCo's `<pid>` with `imax` = 0.07 behaves almost identically (10.5 mrad), with the integral visible in `act`.
### Computing torques yourself versus a servo actuator
The same PD can be computed in your code and written to `ctrl` for a `motor`, or delegated to a `<position>` actuator. The torque is the same; what differs is how the integrator treats the damping. A torque you compute is a constant for the step, so its damping is explicit and stable only for $k_d < 2I/h$ (Lesson 1.3), whatever the integrator. A `<position>` actuator's `kv` is integrated implicitly under `implicitfast`. Output (4) uses $k_v$ = 300 N m s/rad, above the explicit limit of 250: the Python PD diverges under both integrators, the actuator diverges under `Euler` and is stable under `implicitfast`.
> [!recommendation] Which to use
> Use MuJoCo's servo actuators (with `implicitfast`) when the controller is a fixed PD or PID per joint, especially with stiff or heavily damped gains: they are cheaper and more stable. Compute torques yourself when the controller needs the model (gravity compensation, computed torque, operational space: Lessons 8.2 and 8.3), couples joints, or must match a specific real controller's discrete implementation, and keep its damping below $2I/h$.
## Implementation
```io
INPUT: pendulum.xml
PROCESS: measure PD step responses against the damping-ratio formulas; measure and remove the steady-state error under gravity; show windup and anti-windup, in Python and with MuJoCo's pid actuator; compare computed PD with a position actuator under Euler and implicitfast
OUTPUT: printed tables
```python file=examples/l8_1_joint_space.py “"”Lesson 8.1: joint-space control of one joint: P, PD, PID and MuJoCo’s servo actuators.
INPUT pendulum.xml (1 kg at 0.5 m on a hinge, torque motor limited to +-3 N m)
PROCESS (1) gravity off: PD step responses against the damping-ratio predictions;
(2) gravity on: the steady-state error of PD, and removing it with gravity
compensation or with integral action;
(3) integral windup at the torque limit, a hand-written anti-windup clamp, and
MuJoCo’s
Run: python examples/l8_1_joint_space.py “””
import math
import mujoco import numpy as np
from mjcourse import model_path
M, L, G = 1.0, 0.5, 9.81
I_PIVOT = M * L**2 + 1e-4 # kg m^2, from pendulum.xml’s
def pendulum(gravity: bool = True, torque_limit: float = 3.0): model = mujoco.MjModel.from_xml_path(str(model_path(“pendulum”))) if not gravity: model.opt.gravity[:] = 0 model.actuator_ctrlrange[0] = [-torque_limit, torque_limit] return model, mujoco.MjData(model)
def run(model, data, policy, seconds: float) -> np.ndarray: “"”Step with ctrl = policy(data); return the angle after each step.””” out = np.zeros(round(seconds / model.opt.timestep)) for k in range(len(out)): data.ctrl[0] = policy(data) mujoco.mj_step(model, data) out[k] = data.qpos[0] return out
def step_responses() -> None: kp, target = 20.0, 0.5 omega = math.sqrt(kp / I_PIVOT) print(f” kp {kp} N m/rad, I {I_PIVOT} kg m^2 -> natural frequency {omega:.3f} rad/s; step of {target} rad”) print(f” {‘zeta’:>6}{‘kd’:>8}{‘overshoot predicted’:>21}{‘measured’:>10}{‘peak time predicted’:>21}{‘measured’:>10}”) for zeta in (0.1, 0.3, 0.5, 0.7, 1.0): kd = 2 * zeta * math.sqrt(kp * I_PIVOT) model, data = pendulum(gravity=False, torque_limit=100.0) q = run(model, data, lambda d: kp * (target - d.qpos[0]) - kd * d.qvel[0], 3.0) over = 100 * max(q.max() - target, 0) / target pred = 100 * math.exp(-math.pi * zeta / math.sqrt(1 - zeta2)) if zeta < 1 else 0.0 tpk_pred = f”{math.pi / (omega * math.sqrt(1 - zeta2)):.3f} s” if zeta < 1 else “none” tpk = f”{(np.argmax(q) + 1) * model.opt.timestep:.3f} s” if zeta < 1 else “none” print(f” {zeta:>6.1f}{kd:>8.3f}{pred:>19.1f} %{over:>8.1f} %{tpk_pred:>21}{tpk:>10}”)
def gravity_error() -> None: target = 0.5 tau_g = M * G * L * math.sin(target) for kp in (10.0, 20.0, 50.0): kd = 2 * math.sqrt(kp * I_PIVOT) model, data = pendulum() q = run(model, data, lambda d: kp * (target - d.qpos[0]) - kd * d.qvel[0], 6.0) # Equilibrium: kp (target - q) = m g L sin q; to first order e = m g L sin(target) / (kp + m g L cos(target)). pred = tau_g / (kp + M * G * L * math.cos(target)) print(f” PD, kp {kp:>4.0f}: steady-state error {1000 * (target - q[-1]):6.1f} mrad, “ f”first-order prediction {1000 * pred:6.1f} mrad”) kp = 20.0 kd = 2 * math.sqrt(kp * I_PIVOT) model, data = pendulum() q = run(model, data, lambda d: kp * (target - d.qpos[0]) - kd * d.qvel[0] + M * G * L * math.sin(d.qpos[0]), 6.0) print(f” PD + gravity compensation m g L sin q, kp {kp:.0f}: steady-state error {1000 * (target - q[-1]):.4f} mrad”) integral = [0.0]
def pid(d, ki=20.0):
e = target - d.qpos[0]
integral[0] += e * model.opt.timestep
return kp * e - kd * d.qvel[0] + ki * integral[0]
model, data = pendulum()
q = run(model, data, pid, 6.0)
print(f" PID, kp {kp:.0f}, ki 20: error after 6 s {1000 * (target - q[-1]):.4f} mrad")
def windup() -> None: target, kp, ki = 0.55, 20.0, 40.0 # holding 0.55 rad takes 2.56 N m of the 3 N m available kd = 2 * math.sqrt(kp * I_PIVOT) print(f” holding {target} rad needs an integral of m g L sin(target) / ki = {M * G * L * math.sin(target) / ki:.4f} rad s”) for clamp in (None, 0.05, 0.07): integral = [0.0] model, data = pendulum()
def pid(d):
e = target - d.qpos[0]
integral[0] += e * model.opt.timestep
if clamp is not None:
integral[0] = min(max(integral[0], -clamp), clamp)
return kp * e - kd * d.qvel[0] + ki * integral[0]
q = run(model, data, pid, 8.0)
label = "no anti-windup" if clamp is None else f"integral clamped at +-{clamp} rad s"
settle = np.flatnonzero(np.abs(q - target) > 0.01)
print(f" Python PID, {label:<32}: peak {1000 * (q.max() - target):+6.1f} mrad from target, "
f"inside 10 mrad after {(settle[-1] + 1) * model.opt.timestep if len(settle) else 0:.2f} s"
f"{'' if settle[-1] + 1 < len(q) else ' (never: the clamp is below the integral the hold needs)'}")
xml = model_path("pendulum").read_text().replace(
'<motor name="torque" joint="hinge" gear="1" ctrlrange="-3 3"/>',
f'<pid name="pid" joint="hinge" kp="{kp}" kv="{kd:.6f}" ki="{ki}" imax="0.07" input="pos" forcerange="-3 3"/>')
model = mujoco.MjModel.from_xml_string(xml)
data = mujoco.MjData(model)
q = run(model, data, lambda d: target, 8.0)
settle = np.flatnonzero(np.abs(q - target) > 0.01)
print(f" MuJoCo <pid> actuator, ki {ki:.0f}, imax 0.07{'':<15}: peak {1000 * (q.max() - target):+6.1f} mrad from target, "
f"inside 10 mrad after {(settle[-1] + 1) * model.opt.timestep:.2f} s; act holds the integral: {data.act[0]:.4f}")
def servo_versus_python() -> None:
kp, kv, target = 100.0, 300.0, 0.3 # kv beyond the explicit limit 2 I / h = 250 N m s/rad
print(f” kp {kp}, kv {kv} N m s/rad, h 2 ms: explicit damping is stable only for kv < 2 I / h = {2 * I_PIVOT / 0.002:.0f}”)
for integrator in (“mjINT_EULER”, “mjINT_IMPLICITFAST”):
model, data = pendulum(gravity=False, torque_limit=1e6)
model.opt.integrator = getattr(mujoco.mjtIntegrator, integrator)
q = run(model, data, lambda d: kp * (target - d.qpos[0]) - kv * d.qvel[0], 4.0)
xml = model_path(“pendulum”).read_text().replace(
‘
if name == “main”:
print(“(1) PD step responses, gravity off”)
step_responses()
print(“(2) under gravity: steady-state error, and two ways to remove it (target 0.5 rad)”)
gravity_error()
print(“(3) integral windup at the 3 N m limit (target 0.55 rad, kp 20, ki 40)”)
windup()
print(“(4) the same PD in Python and as a
Output:
```text
(1) PD step responses, gravity off
kp 20.0 N m/rad, I 0.2501 kg m^2 -> natural frequency 8.942 rad/s; step of 0.5 rad
zeta kd overshoot predicted measured peak time predicted measured
0.1 0.447 72.9 % 72.9 % 0.353 s 0.352 s
0.3 1.342 37.2 % 37.1 % 0.368 s 0.366 s
0.5 2.237 16.3 % 16.1 % 0.406 s 0.404 s
0.7 3.131 4.6 % 4.4 % 0.492 s 0.490 s
1.0 4.473 0.0 % 0.0 % none none
(2) under gravity: steady-state error, and two ways to remove it (target 0.5 rad)
PD, kp 10: steady-state error 162.4 mrad, first-order prediction 164.4 mrad
PD, kp 20: steady-state error 96.3 mrad, first-order prediction 96.8 mrad
PD, kp 50: steady-state error 43.3 mrad, first-order prediction 43.3 mrad
PD + gravity compensation m g L sin q, kp 20: steady-state error 0.0000 mrad
PID, kp 20, ki 20: error after 6 s -0.0316 mrad
(3) integral windup at the 3 N m limit (target 0.55 rad, kp 20, ki 40)
holding 0.55 rad needs an integral of m g L sin(target) / ki = 0.0641 rad s
Python PID, no anti-windup : peak +87.2 mrad from target, inside 10 mrad after 1.62 s
Python PID, integral clamped at +-0.05 rad s: peak -22.9 mrad from target, inside 10 mrad after 8.00 s (never: the clamp is below the integral the hold needs)
Python PID, integral clamped at +-0.07 rad s: peak +9.4 mrad from target, inside 10 mrad after 0.52 s
MuJoCo <pid> actuator, ki 40, imax 0.07 : peak +10.5 mrad from target, inside 10 mrad after 0.86 s; act holds the integral: 0.0641
(4) the same PD in Python and as a <position> actuator, gravity off, step to 0.3 rad
kp 100.0, kv 300.0 N m s/rad, h 2 ms: explicit damping is stable only for kv < 2 I / h = 250
euler PD in Python: angle after 4 s -1463 rad (max |q| 1.46e+03); <position> actuator: -81.59 rad (max |q| 1.27e+04)
implicitfast PD in Python: angle after 4 s -1463 rad (max |q| 1.46e+03); <position> actuator: 0.221 rad (max |q| 0.221)
MuJoCo also printed WARNING: Nan, Inf or huge value in QACC for the diverging runs in (4) and reset them; the final angles are what remained.
The joint oscillates forever at a moderate gain. No damping, or $k_d$ set without regard to the inertia. Compute $\zeta$ from the joint’s actual inertia (mj_fullM diagonal plus armature) and set $k_d = 2\zeta\sqrt{k_pI}$.
The joint stops a few degrees short of the target. Gravity against a finite $k_p$ (output 2). Add gravity compensation, or integral action with a clamp.
After a large move, the joint overshoots and then crawls back. Integral windup during saturation (output 3). Clamp the integral above the steady-state value it needs, or stop integrating while the motor is saturated.
Raising $k_d$ makes the joint blow up. Explicit damping beyond $2I/h$. Use a <position> actuator under implicitfast, lower $k_d$, or reduce the timestep (output 4).
For the pendulum held at 0.5 rad by PD with $k_p$ = 20, measure the steady-state error, then double the bob’s mass and measure again. Predict both from the formula above before running.
Give arm7 a <pid> actuator on every joint with implicitfast, choose $k_p$ for a 3 Hz bandwidth and $\zeta$ = 0.7 at the home pose’s diagonal inertias, choose $k_i$ and imax from the steady-state gravity torques, and step every joint by 0.3 rad. Report overshoot and settling per joint, and explain why the joints with the most armature behave most like the single-joint prediction.
[!research] Low-level controllers shape what policies learn Most learned manipulation policies output joint or end-effector targets for a PD loop like this one, not torques. The gains decide how fast and how compliantly targets are reached, so a policy trained with one set of gains can fail with another on the same robot, and reported results depend on gains that papers often omit (inference). Domain randomization over controller gains (Level 17) and identifying the real robot’s servo behaviour (Level 18) both start from the formulas in this lesson.
{"id": "8.1-check", "title": "Knowledge check", "questions": [
{"kind": "numeric", "q": "A joint with inertia 0.1 kg m² has kp = 40 N m/rad. What kd (N m s/rad) gives a damping ratio of 0.7?",
"answer": 2.8, "tol": 0.005, "unit": "N m s/rad",
"explain": "<p>$k_d = 2\\zeta\\sqrt{k_p I} = 1.4 \\times \\sqrt{4}$ = 2.8.</p>"},
{"kind": "numeric", "q": "What overshoot (percent of the step) does a damping ratio of 0.5 produce?",
"answer": 16.3, "tol": 0.1, "unit": "%",
"explain": "<p>$e^{-\\pi \\cdot 0.5/\\sqrt{0.75}}$ = 0.163; output (1) measured 16.1%.</p>"},
{"kind": "mcq", "q": "An integral term with anti-windup clamp 0.05 rad s never brings the pendulum to its target, while 0.07 does. Why?",
"options": ["The clamp is too large", "Holding the target needs an integral of about 0.064 rad s, which the 0.05 clamp forbids", "Integral action cannot remove gravity errors", "The motor limit is too high"],
"answer": 1,
"explain": "<p>At steady state the integral alone must supply the gravity torque: $k_i I = mgL\\sin q^\\star$, so $I$ = 2.56/40 = 0.064 rad s.</p>"},
{"kind": "mcq", "q": "Under implicitfast, a PD with kv beyond 2I/h is stable as a position actuator but unstable when the same torque is computed in Python and written to ctrl. Why?",
"options": ["Python is slower", "MuJoCo integrates the actuator's velocity term implicitly, but treats a torque in ctrl as a constant for the step", "Position actuators have no damping", "implicitfast ignores ctrl"],
"answer": 1,
"explain": "<p>Output (4). The integrator can only treat implicitly the damping it knows about.</p>"}
]}
Lesson 8.2 controls the whole arm along a trajectory and adds the model: gravity compensation, feedforward of the inverse dynamics, and computed-torque control, each measured with a correct model and with a wrong one.