Robustness, Calibration, and Safe Control

Michael BrenndoerferAugust 7, 202658 min read

Part of World Models Handbook

Stress testing, calibrated uncertainty, risk-sensitive planning, and runtime assurance for safe control with imperfect world models.

Choose your expertise level to adjust how many terms are explained. Beginners see more tooltips, experts see fewer to maintain reading flow. Hover over underlined terms for instant definitions.

Article links

Make inline references clickable

Robustness, Calibration, and Safe Control

A hypothetical warehouse robot learns a world model on a dry floor. Your planner trusts narrow prediction intervals when choosing a maneuver. If the floor becomes wet, the relationship between action and motion can change without the model's confidence changing immediately. You may then choose an action that was reasonable under the old dynamics but inappropriate under the new ones. This is a motivating scenario, not a reported robot experiment; the executable example below uses a synthetic oscillator not a wheel-torque model.

Accuracy on the training regime does not establish accuracy after this change. What failed is the assumption that the old transition law and uncertainty estimates still apply. This chapter examines how to test those assumptions and how to limit reliance on them.

The preceding chapter, Failure Modes and Model Exploitation, concerns how world models break. Here we develop an engineering response through four connected pieces:

  • Distribution shift and stress testing: how to evaluate specified changes before deployment and monitor observed mismatches.
  • Calibrated uncertainty and abstention: how to make the model's own confidence a usable decision signal, and how to decline to act when it is not.
  • Risk-sensitive and constrained planning: how to plan against a distribution of futures rather than a single imagined one.
  • Fallback policies and runtime assurance: how to separate a performance controller from a fallback, and what must be proved before calling that fallback safe.

These topics interact. Stress testing identifies specified changes worth evaluating; calibration measures whether predictive probabilities agree with outcomes in an evaluated population. Neither is a universal shift detector. Risk-sensitive planning changes how sampled futures are scored, while runtime assurance checks whether the performance controller may continue operating under independently justified assumptions. Each link needs evidence: a calibration estimate does not prove future coverage, a risk objective does not enforce a hard constraint, and a fallback is not certified merely because its test trajectories stay in bounds.

The central asymmetry

A world model is a predictor. A controller is a decision rule. Good one-step prediction does not establish safe closed-loop behavior. A fallback can reduce dependence on a weak model, but safety requires verified assumptions about the plant, disturbances, state estimates, timing, switching logic, and fallback's operating region. The toy fallback in this chapter has no such certification.

Throughout, we use a one-dimensional positioning task with a damped oscillator as the plant. The printed tables come from executable code, so you can change the assumptions and observe the consequences. This modest example exposes the distinction between prediction accuracy, uncertainty estimates, and closed-loop outcomes. Larger learned systems introduce additional representation and estimation questions; scale alone does not establish that confidence tracks competence.

Distribution Shift and Stress Testing

Distribution shift is the umbrella term for any mismatch between the distribution the world model was trained on and the distribution it is asked to act in. It is worth being precise about what can shift, because different kinds of shift demand different responses, and because a model's own uncertainty carries different amounts of information about each type.

  • Input shift. The distribution of states and actions changes. In pure covariate shift, the conditional transition law given those inputs stays fixed. Observed inputs permit label-free distribution tests, but detection power depends on the change, representation, sample size, and test; visibility does not guarantee detection.
  • Dynamics shift (concept shift). The conditional transition law changes, potentially at inputs that initially look familiar. The floor got wet; the payload got heavier; the actuator got slower. The model can be confidently wrong without obvious input novelty, and changed dynamics may alter the states visited later. This is the failure mode that motivates the entire chapter.
  • Observation shift. The sensor or observation pathway changes: a new camera, changed calibration, lighting, or noise profile. An inaccurate state estimate can undermine an otherwise adequate transition model. Observation diagnostics and transition residuals can both signal a problem, but residuals alone need not identify its cause.
  • Task shift. Goals, rewards, or constraints change. A change from docking to transit may alter the objective and also change visited states and actions. A new operating regime alone is input shift unless the task definition changes too.
  • Hidden-parameter shift. A latent parameter value or its distribution changes, abruptly or gradually. Within a parameterized plant family, drifting friction is one mechanism of dynamics shift; it is not a definition covering every structural change.

The distinction affects what ensemble disagreement can reveal. The members may agree on a wrong transition when the conditional dynamics change. Conversely, disagreement can increase on novel inputs, but that response is not guaranteed: shared model bias can keep all members confident outside the data. We measure these responses in a particular synthetic example rather than assuming either is a reliable detector.

Why stress testing is not the same as held-out validation

Held-out validation measures prediction on its chosen test population. A same-generator split is useful for checking generalization within that generator, but may miss deployment conditions absent from it. Safety-relevant samples can include both rare in-distribution cases and out-of-distribution conditions. A large average RMSE improvement can coexist with serious boundary errors; report conditional and closed-loop outcomes rather than assuming the aggregate captures them.

Stress testing changes the population deliberately: under which plausible conditions does performance degrade, and which consequences follow? It can falsify a claim while also producing measurements. A useful output includes quantitative outcomes and an inventory of tested mechanisms, failures, and untested conditions. Failure to find a problem is evidence about the search performed, not proof of safety or evidence that the search was necessarily inadequate.

Good stress tests are built from a shift taxonomy for the specific domain, in which each category probes a distinct mechanism by which reality can diverge from training:

  • Parameter stress. Perturb parameters over an explicitly justified domain, for example reducing friction, increasing mass, or reducing actuator gain. Search for adverse values within that domain; physical constraints may make it bounded.
  • Input-distribution stress. Test high velocities, large commands, or near-boundary positions with sparse training coverage. Ensemble disagreement may rise there, but sparse coverage does not guarantee high estimated uncertainty.
  • Structural stress. Add an effect outside the model class, such as an unrepresented delay or coupling. More data cannot remove that representational limitation while the class remains fixed. This tests a different limitation from parameter error.
  • Temporal stress. Change the control rate, introduce jitter, or hold commands longer. Whether errors remain tolerable depends on the plant, controller, and timing. Evaluate closed-loop consequences rather than assuming a universal relation between rate difference and degradation.

When a plausible stress condition breaks the model, options include broadening data, changing the model class, restricting the operating domain, redesigning the controller or hardware, or adding detection and fallback. These approaches can complement one another but none universally handles rare or unbounded disturbances. A fallback needs its own operating assumptions. Part VI: Learning World Models at Scale covers learning and adaptation; here we evaluate diagnostic and control responses.

What we are building

Let us set up the running example. The plant is a mass on a spring, integrated with a semi-implicit Euler step at Δt=0.1 \Delta t = 0.1\,s. The state at step tt is the pair (pt,vt)(p_t, v_t). The update advances velocity first and then uses that new velocity to advance position.

at+1=−k pt−c vt−f sign(vt)+g utvt+1=vt+Δt at+1+ϵt,ϵt∼N(0,σp2)pt+1=pt+Δt vt+1\begin{aligned} a_{t+1} &= -k\,p_t - c\,v_t - f\,\mathrm{sign}(v_t) + g\,u_t \\ v_{t+1} &= v_t + \Delta t\, a_{t+1} + \epsilon_t, \qquad \epsilon_t \sim \mathcal{N}(0, \sigma_p^2) \\ p_{t+1} &= p_t + \Delta t\, v_{t+1} \end{aligned}

where:

  • ptp_t: position at step tt, the coordinate the safety envelope constrains
  • vtv_t: velocity at step tt
  • at+1a_{t+1}: the acceleration computed at step tt and used to advance the velocity
  • utu_t: the scalar control action applied at step tt; it is the book's action ata_t under a local notation chosen to distinguish it from acceleration
  • Δt=0.1 \Delta t = 0.1\,s: integration timestep
  • kk: stiffness, the restoring force per unit displacement; cc: viscous damping, the drag proportional to velocity; ff: dry friction, a constant drag opposing motion; gg: actuator gain, how much acceleration one unit of command produces
  • ϵt∼N(0,σp2)\epsilon_t \sim \mathcal{N}(0, \sigma_p^2): process noise on velocity, the part of the dynamics no deterministic model can capture

The first equation computes acceleration; the second advances velocity and adds its innovation; the third uses that new velocity to advance position. Here at+1a_{t+1} denotes acceleration, not the action: the scalar action is utu_t. We observe the full simulated state st=(pt,vt)s_t=(p_t,v_t), so no observation model or belief filter is implemented. Mass is normalized to one and action units are abstract; this oscillator is not a wheel-torque model. The simulator permits dry friction, but every experiment here sets f=0f=0. Semi-implicit Euler advances velocity before position; its stability depends on the timestep and dynamics, not on the method name alone.

In[3]:
Code
import numpy as np
from scipy.stats import norm

DT = 0.1  # simulation step, seconds
U_VALID = 0.35  # action range covered by the training data
U_FALLBACK = 1.20  # stronger actuator authority used by the heuristic fallback
P_LIMIT = 0.60  # tested position envelope, not a proved invariant
SIGMA_P = 0.06  # process-noise standard deviation on velocity


def simulate(params, s0, actions, rng=None, sigma_p=0.0):
    """Semi-implicit Euler rollout of a damped oscillator with dry friction."""
    k, c, gain, fric = params
    p, v = float(s0[0]), float(s0[1])
    traj = np.empty((len(actions) + 1, 2))
    traj[0] = (p, v)
    for t, u in enumerate(actions):
        acc = -k * p - c * v - fric * np.sign(v) + gain * u
        v = v + DT * acc
        if sigma_p > 0.0:
            v = v + rng.normal(0.0, sigma_p)
        p = p + DT * v
        traj[t + 1] = (p, v)
    return traj

Before fitting anything, it is worth seeing the two regimes we keep returning to. The cell below rolls one outward initial condition and one constant braking command through the nominal plant and through a slippery floor whose damping has collapsed, leaving every other parameter unchanged.

In[4]:
Code
# Two parameter tuples for a single illustrative rollout. The nominal tuple sits
# at the centre of the training ranges; the slippery tuple keeps stiffness and
# gain nominal and collapses damping toward zero. Damping is the only thing that
# changes, so any difference in behaviour is due to a hidden-parameter shift.
NOMINAL_PARAMS = (1.00, 0.90, 1.00, 0.0)
# A fixed low-damping regime for an illustrative, noise-free comparison.
SLIPPERY_PARAMS = (1.00, 0.05, 1.00, 0.0)

TRAJ_RNG = np.random.default_rng(20240518)  # fixed seed, rollout is noise-free
TRAJ_STEPS = 60
traj_s0 = np.array([0.30, 0.60])  # outward position and velocity
traj_actions = np.full(TRAJ_STEPS, -0.35)  # constant braking command

traj_nominal = simulate(NOMINAL_PARAMS, traj_s0, traj_actions, TRAJ_RNG, 0.0)
traj_slippery = simulate(SLIPPERY_PARAMS, traj_s0, traj_actions, TRAJ_RNG, 0.0)
traj_time = np.arange(TRAJ_STEPS + 1) * DT
Out[5]:
Visualization
Position over time for nominal and low-damping oscillators with a safety limit line.
Noise-free position rollouts under nominal and low damping with the same initial state and constant negative command. Lower damping changes the oscillation and excursions. Holding this command indefinitely is not a stop policy; this illustration does not certify either trajectory.

The exploration policy samples abstract commands in [−0.35,0.35][-0.35,0.35] from moderate initial conditions. Each episode draws a hidden parameter tuple. With only current state and command as inputs, parameter variation can produce different conditional outcomes at the same input. Even fixed parameters leave Gaussian process noise in this experiment. The regression estimates a conditional mean; its residual spread includes unresolved parameter variation, process noise, and fitting error rather than uncertainty from estimation alone.

In[6]:
Code
K_RANGE = (0.80, 1.20)
C_RANGE = (0.70, 1.10)
GAIN_RANGE = (0.90, 1.10)

P_SCALE = 0.30  # initial-position scale
V_SCALE = 0.75  # initial-velocity scale
EP_LEN = 40  # steps per episode


def sample_nominal(rng):
    return (
        rng.uniform(*K_RANGE),
        rng.uniform(*C_RANGE),
        rng.uniform(*GAIN_RANGE),
        0.0,
    )


def collect_split(
    rng,
    n_episodes,
    params_fn,
    p_scale=P_SCALE,
    v_scale=V_SCALE,
    u_scale=U_VALID,
    ep_len=EP_LEN,
):
    """Roll out episodes and return stacked (s, a) -> delta-s transitions."""
    Xs, Ys = [], []
    for _ in range(n_episodes):
        s0 = rng.uniform(-1.0, 1.0, 2) * np.array([p_scale, v_scale])
        actions = rng.uniform(-u_scale, u_scale, ep_len)
        traj = simulate(params_fn(rng), s0, actions, rng, SIGMA_P)
        Xs.append(np.column_stack([traj[:-1], actions]))
        Ys.append(traj[1:] - traj[:-1])
    return np.vstack(Xs), np.vstack(Ys)


rng = np.random.default_rng(20240517)
X_train, Y_train = collect_split(rng, 100, sample_nominal)
X_calib, Y_calib = collect_split(rng, 25, sample_nominal)
X_recal, Y_recal = collect_split(rng, 25, sample_nominal)
X_test, Y_test = collect_split(rng, 25, sample_nominal)

The quadratic feature map below produces the 10-dimensional feature vector

ϕ(p,v,u)=[1,  p,  v,  u,  p2,  v2,  u2,  pv,  pu,  vu]⊤,\phi(p, v, u) = \left[1,\; p,\; v,\; u,\; p^2,\; v^2,\; u^2,\; pv,\; pu,\; vu\right]^\top,

where:

  • 11: the constant term, allowing an affine offset in the prediction
  • p,v,up, v, u: the raw state and action coordinates
  • p2,v2,u2p^2, v^2, u^2: the squared terms, which let the model curve the prediction within the input range
  • pv,pu,vupv, pu, vu: the pairwise cross terms, which let the effect of one input depend on the value of another
  • ϕ\phi has 1010 entries in total, so it maps each input triple to a 1010-dimensional feature vector.

The quadratic feature map can represent each fixed-parameter plant's linear conditional mean exactly because dry friction is zero here. The difficulty is omitted information: parameters vary across episodes, while the predictor receives only current state and action, not an episode parameter estimate or history. This produces unresolved variation for this representation. A history-conditioned identifier could reduce some of it; it is not necessarily irreducible for every world model. The fitted residual variance also contains finite-fit errors and possible bias, not just process noise. We use that approximate spread for the diagnostics below.

In[7]:
Code
def features(X):
    """Quadratic features of the (position, velocity, action) triple."""
    p, v, u = X[:, 0], X[:, 1], X[:, 2]
    return np.column_stack(
        [np.ones_like(p), p, v, u, p * p, v * v, u * u, p * v, p * u, v * u]
    )


def fit_ensemble(Xtr, Ytr, Xcal, Ycal, n_members, rng):
    """Bootstrap-bagged ridge-free least squares, with a residual-variance head."""
    Phi_tr = features(Xtr)
    Phi_cal = features(Xcal)
    members = []
    for _ in range(n_members):
        idx = rng.integers(0, len(Xtr), len(Xtr))
        W, *_ = np.linalg.lstsq(Phi_tr[idx], Ytr[idx], rcond=None)
        resid = Ycal - Phi_cal @ W
        members.append((W, resid.var(axis=0)))
    return members


RNG_MEMBERS = np.random.default_rng(11)
# Nine independent transition-bootstrap fits. Their spread is only a proxy
# for fit uncertainty, not a calibrated posterior over physical parameters.
members = fit_ensemble(X_train, Y_train, X_calib, Y_calib, 8, RNG_MEMBERS)
extra_member = fit_ensemble(
    X_train, Y_train, X_calib, Y_calib, 1, np.random.default_rng(999)
)[0]
members.append(extra_member)
ALEATORIC = np.mean([var for _, var in members], axis=0)


def predict(members, X):
    """Return (mean, total variance, epistemic variance) for a batch of inputs."""
    Phi = features(X)
    stack = np.stack([Phi @ W for W, _ in members])  # (M, N, 2)
    mean = stack.mean(axis=0)
    epistemic = stack.var(axis=0)
    residual_floor = np.mean(
        [residual_var for _, residual_var in members], axis=0
    )
    return mean, epistemic + residual_floor, epistemic

Two variances are being combined, and the split matters:

  • The epistemic proxy is the spread of the member predictions, Varm[μm(x)]=1M∑m=1M(μm(x)−μˉ(x))2\mathrm{Var}_m[\mu_m(x)] = \frac{1}{M}\sum_{m=1}^{M}(\mu_m(x)-\bar\mu(x))^2, where μˉ(x)\bar\mu(x) is their mean and MM their count. It measures disagreement among these fits; it may grow on novel inputs but need not reveal shared bias.
  • The held-out residual-variance floor is 1M∑m=1MVar(y−μm(x))\frac{1}{M}\sum_{m=1}^{M}\mathrm{Var}(y-\mu_m(x)). It includes process noise, unresolved information for this predictor, and varying fitting errors. It is not a pure estimate of irreducible aleatoric uncertainty; a history-conditioned model may resolve some hidden-parameter variation. Centered variance also does not measure a constant residual bias.

Their sum is the code's approximate predictive variance, not a proven decomposition of true uncertainty. We keep disagreement separate to inspect its response to input novelty. The following output compares realized error with the claimed spread; calibration and decision utility require their own checks.

The next output compares one-step RMSE and mean predicted standard deviation by component. Similar values are a useful scale diagnostic, not proof of calibration; the coverage experiment below is a separate check. Complete trajectories are assigned to different data splits, but the transition bootstrap within each training split ignores temporal dependence and should not be interpreted as a statistically calibrated posterior.

In[8]:
Code
mu, var, epi = predict(members, X_test)
rmse = np.sqrt(np.mean((Y_test - mu) ** 2, axis=0))
mean_pred_std = np.sqrt(var).mean(axis=0)
Out[9]:
Console
one-step RMSE  (delta-p, delta-v): 0.00614, 0.06144
mean predictive std (delta-p, delta-v): [0.00604842 0.06048423]

Building the stress suite

We build two stress sets. The input-shift generator increases initial position and velocity scales and draws abstract actions from [−0.70,0.70][-0.70,0.70], twice the training range, while retaining nominal plant parameter sampling. The dynamics-shift generator retains training-scale initial conditions and commands but draws damping from [0.00,0.05][0.00,0.05] instead of [0.70,1.10][0.70,1.10]. These are distinct interventions on the generator. Later states can change under either intervention, so a dynamics change does not imply that every encountered input remains in-distribution. The measured response of the ensemble, not input visibility alone, tells us what its diagnostic does here.

In[10]:
Code
def collect_input_shift(rng_, n_episodes=25):
    """Expanded initial-state and action ranges with nominal plant sampling."""
    Xs, Ys = [], []
    for _ in range(n_episodes):
        s0 = rng_.uniform(-1.0, 1.0, 2) * np.array([0.95, 1.60])
        actions = rng_.uniform(-0.70, 0.70, EP_LEN)
        traj = simulate(sample_nominal(rng_), s0, actions, rng_, SIGMA_P)
        Xs.append(np.column_stack([traj[:-1], actions]))
        Ys.append(traj[1:] - traj[:-1])
    return np.vstack(Xs), np.vstack(Ys)


def sample_light_damping(rng_):
    """The floor got slippery: damping collapses, everything else is nominal."""
    return (
        rng_.uniform(*K_RANGE),
        rng_.uniform(0.00, 0.05),
        rng_.uniform(0.90, 1.10),
        0.0,
    )


X_in, Y_in = collect_input_shift(rng)
X_dyn, Y_dyn = collect_split(rng, 25, sample_light_damping)


def evaluate_regime(X, Y):
    mu_r, var_r, epi_r = predict(members, X)
    residual = Y - mu_r
    return dict(
        rmse=float(np.sqrt(np.mean(residual**2))),
        pred_sd=float(np.sqrt(var_r).mean()),
        epi_sd=float(np.sqrt(np.maximum(epi_r, 0.0)).mean()),
        rms_z=float(np.sqrt(np.mean(residual**2 / var_r))),
    )


stress_rows = [
    ("in-distribution", evaluate_regime(X_test, Y_test)),
    ("input shift", evaluate_regime(X_in, Y_in)),
    ("dynamics shift", evaluate_regime(X_dyn, Y_dyn)),
]

We record error, predicted spread, standardized residual scale, and ensemble disagreement separately. Here RMSz=N−1∑i(yi−μ(xi))2/σ2(xi)\mathrm{RMS}_z=\sqrt{N^{-1}\sum_i (y_i-\mu(x_i))^2/\sigma^2(x_i)}. Correct conditional means and finite positive conditional variances imply E[RMSz2]=1\mathbb E[\mathrm{RMS}_z^2]=1. Values above one indicate larger squared errors than the predicted scale accounts for; values below one indicate the reverse. Bias and tail shape can affect this diagnostic, so neither result alone establishes interval calibration or useful conservatism. Disagreement may respond to novelty, but need not do so.

Out[11]:
Console
regime                 RMSE   pred sd  epistemic   RMS z
in-distribution      0.0437    0.0333     0.0012    1.02
input shift          0.0451    0.0337     0.0042    1.03
dynamics shift       0.0459    0.0333     0.0021    1.06

The table reports RMSE, predicted standard deviation, epistemic standard deviation, and RMS z for each regime.

For these seeds, compare the epistemic column with the total error and RMS z. The dynamics perturbation changes the plant, yet its aggregate one-step effect is small relative to the injected process noise. This experiment therefore does not demonstrate catastrophic one-step error or a dramatic calibration collapse. It shows why a modest change in aggregate error is not evidence that closed-loop behavior is unchanged. Although the shifted episodes start at training-scale initial states and actions, their later states need not have the same distribution.

Ensemble disagreement alone cannot establish whether the plant still obeys its old dynamics. A monitor can also use observed prediction residuals, independently justified disturbance bounds, or physical constraints. Each signal has blind spots, so stress tests should evaluate detection delay and false alarms, not merely the model's average uncertainty.

Out[12]:
Visualization
Grouped bar chart of prediction error and predicted standard deviation for three evaluation regimes.
Grouped bars compare aggregate one-step RMSE (root mean square across samples and coordinates) with the arithmetic mean of predicted coordinate standard deviations. These different aggregations and coordinate scales can produce different bar heights even for well-scaled component uncertainties; use the per-coordinate output and coverage curves for calibration assessment.

The two coordinates have different scales, so an aggregate comparison can hide component-specific behavior. Inspect the earlier per-coordinate output and the reliability curves as well. This plot is a summary of the specified tests, not a general validation of the uncertainty estimator.

Calibrated Uncertainty and Abstention

A model is calibrated if its stated probabilities match observed frequencies. Concretely: among all the predictions where the model says "there is a 90% chance the true value lies in this interval," roughly 90% should contain the truth. Calibration is a property of the whole predictive distribution, not of any single prediction, so it can only be measured in aggregate. This is a subtle but important point: you cannot ask whether a single interval is calibrated, any more than you can ask whether a single coin flip is fair. You ask it about the set of predictions that share the same confidence level, and you compare the claimed frequency to the realized one.

Calibration vs. sharpness

Coverage and width answer different questions. An interval of (−∞,∞)(-\infty,\infty) has coverage one, not calibrated 90% coverage. A useful interval should be narrow while achieving its stated coverage on the relevant population. Report both width and coverage; also inspect operationally important subgroups rather than relying only on an aggregate average.

In Bayesian Filtering and Belief States, a posterior represents state uncertainty; Stochasticity and Uncertainty distinguishes uncertainty sources. Here we additionally approximate predictions by N(μ(x),σ2(x))\mathcal N(\mu(x),\sigma^2(x)). Producing a mean and variance does not imply Gaussian shape. Quantile-based coverage can be tested without finite moments; standardized residuals using mean and standard deviation require finite moments and a positive scale.

Measuring coverage

For a Gaussian predictive distribution, the natural diagnostic is the standardized residual. If the model is honest, each observation yiy_i should look like a draw from N(μ(xi),σ2(xi))\mathcal{N}(\mu(x_i), \sigma^2(x_i)), so dividing the raw error by the claimed standard deviation puts every prediction on a common scale:

zi=yi−μ(xi)σ(xi).z_i = \frac{y_i - \mu(x_i)}{\sigma(x_i)}.

where:

  • yiy_i: the observed target for the ii-th sample
  • μ(xi)\mu(x_i): the model's predicted mean at input xix_i
  • σ(xi)\sigma(x_i): the model's predicted standard deviation at input xix_i
  • ziz_i: the standardized residual, measuring how many predicted standard deviations the observation deviates from the mean

For the Gaussian location-scale intervals used here, calibration at all central levels requires the absolute standardized residuals to have the corresponding Gaussian quantiles. Calibration alone does not imply that signed residuals are standard normal. The empirical coverage at nominal level α\alpha is:

c^(α)=1N∑i=1N1 ⁣[ ∣zi∣≤Φ−1 ⁣(1+α2)]≈α.\hat{c}(\alpha) = \frac{1}{N} \sum_{i=1}^{N} \mathbf{1}\!\left[\,|z_i| \le \Phi^{-1}\!\left(\tfrac{1+\alpha}{2}\right) \right] \approx \alpha .

where:

  • NN: the number of evaluated samples
  • α∈(0,1)\alpha \in (0,1): the nominal coverage level, the fraction of probability mass the interval is supposed to contain
  • Φ−1\Phi^{-1}: the standard normal quantile function, so Φ−1 ⁣(1+α2)\Phi^{-1}\!\left(\tfrac{1+\alpha}{2}\right) is the half-width of the central interval containing fraction α\alpha of a standard normal
  • 1[⋅]\mathbf{1}[\cdot]: the indicator, equal to 11 when the inequality holds and 00 otherwise
  • c^(α)\hat{c}(\alpha): the empirical fraction of samples whose standardized residual falls inside that central interval

Plotting c^(α)\hat{c}(\alpha) against α\alpha gives a reliability diagram. The diagonal represents matching marginal coverage at the displayed levels; curves below it deliver less coverage than claimed, and curves above it deliver more. Sampling uncertainty remains, and the diagram does not establish conditional coverage at every state. Overcoverage may accompany unnecessarily wide intervals; undercoverage matters when decisions rely on margins that are too narrow. Neither direction alone measures controller safety.

In[13]:
Code
LEVELS = np.array([0.50, 0.68, 0.80, 0.90, 0.95])


def z_scores(X, Y, scale=1.0):
    mu_z, var_z, _ = predict(members, X)
    return ((Y - mu_z) / np.sqrt(var_z * scale**2)).ravel()


def coverage_curve(z, levels=LEVELS):
    return np.array(
        [(np.abs(z) <= norm.ppf(0.5 + lv / 2.0)).mean() for lv in levels]
    )


def fit_scale(X, Y):
    """One global multiplier that makes the average squared z-score equal to one."""
    return float(np.sqrt(np.mean(z_scores(X, Y) ** 2)))


global_scale = fit_scale(X_recal, Y_recal)

z_id = z_scores(X_test, Y_test)
z_dyn = z_scores(X_dyn, Y_dyn)
cover_id = coverage_curve(z_id)
cover_id_scaled = coverage_curve(z_scores(X_test, Y_test, global_scale))
cover_dyn = coverage_curve(z_dyn)
oracle_scale = float(np.sqrt(np.mean(z_dyn**2)))
cover_dyn_scaled = coverage_curve(z_scores(X_dyn, Y_dyn, oracle_scale))
Out[14]:
Console
global scale fitted on in-distribution calibration data: 1.000
scale that would be needed on the shifted regime:        1.065

  nominal   nominal cov   shifted cov
     0.50         0.492         0.493
     0.68         0.667         0.671
     0.80         0.789         0.775
     0.90         0.891         0.862
     0.95         0.947         0.929

The executed global scale is 1.000 and the shifted-set diagnostic is 1.065. At nominal 90% coverage, the measured frequencies are 89.1% on nominal test transitions and 86.2% on shifted transitions. This is a modest empirical coverage loss in this experiment, not evidence of an arbitrary catastrophic collapse.

Read the two coverage columns at the same nominal level. They are measured frequencies for the specified data splits, not guarantees. The shifted-data scale is an oracle diagnostic fitted and evaluated on that same shifted set; it is not a deployable recalibration result or a held-out validation estimate.

The fitted multiplier adjusts average squared standardized residuals, not every conditional probability or tail. It may help on some shifted populations, but coverage does not automatically transfer. Our calibration targets are correlated transitions from independent episodes, so the curve describes their empirical mixture; it is not a finite-sample confidence guarantee for an adaptive controller.

The code pools both coordinate residuals. In this simulator, Δp=Δt (vt+Δv)\Delta p=\Delta t\,(v_t+\Delta v), not Δt Δv\Delta t\,\Delta v. The feature map includes current velocity, so the fitted predictions satisfy Δp^=Δt (vt+Δv^)\widehat{\Delta p}=\Delta t\,(v_t+\widehat{\Delta v}), up to numerical roundoff. Subtracting these predictions from their targets leaves position residuals equal to Δt\Delta t times velocity residuals; the standardized coordinate residuals therefore coincide up to numerical roundoff. In a general multi-output model, pooled coverage can hide opposing coordinate errors; report coordinate-specific diagnostics as well.

Split conformal prediction provides a different route to marginal coverage: fit a model on separate training data, compute nonconformity scores on nn calibration examples, and use the appropriate order statistic, typically rank ⌈(n+1)(1−δ)⌉\lceil(n+1)(1-\delta)\rceil, for target coverage 1−δ1-\delta. If the rank exceeds nn, an infinite threshold is the conservative convention. The finite-sample marginal result requires exchangeability of the calibration and test examples, conditional on the independently fitted predictor. Consecutive transitions from the same episode and adaptive control generally do not satisfy that assumption. Marginal coverage also does not guarantee coverage at every state or for an entire trajectory. This chapter implements Gaussian residual rescaling, not conformal intervals; see Angelopoulos and Bates for the distinction and conditions.

Out[15]:
Visualization
Nominal empirical coverage compared with nominal levels.
Nominal held-out transition coverage at five Gaussian interval levels, with and without a global scale fitted on a separate calibration split. Matching aggregate frequencies do not establish conditional or trajectory coverage.
Out[16]:
Visualization
Shifted empirical coverage with a diagnostic oracle rescaling.
Coverage under the damping shift with raw uncertainty and a shifted-data oracle scaling. The oracle uses the same outcomes for fitting and evaluation, so it is diagnostic rather than independently validated.

Abstention: turning uncertainty into a decision

Calibration describes predictions on a population. Abstention is a decision rule built on top of an uncertainty score. Here we rank predictions by predictive variance and evaluate only a retained fraction. An informative ranking may reduce retained-set error; it does not entitle a model to act safely. An actual abstention system must also specify what happens to the rejected cases and evaluate that alternative.

Possible designs include a planner that defers high-disagreement proposals, a manipulator that requests confirmation, or a controller that selects a predefined fallback. These are examples of design choices, not claims about every deployed system. A percentile trigger asks whether a score is unusual relative to a reference population; an absolute trigger asks whether a specified tolerance is exceeded. Either can be appropriate when its operating population and consequences have been validated.

One way to evaluate this ranking is the risk-coverage curve. Sweep the retained fraction ρ\rho, compute retained-set RMSE, and plot it against ρ\rho. The area summarizes this particular loss and population, so report the loss definition and coverage range. An uncertainty score independent of errors gives approximately unchanged expected error under random retention; finite samples fluctuate. Useful rankings often give lower error at lower coverage, but empirical curves need not be monotone.

In[17]:
Code
def risk_coverage(X, Y, fractions):
    """RMSE of the retained predictions as a function of the retained fraction."""
    mu_rc, var_rc, _ = predict(members, X)
    per_point_err = np.sqrt(np.mean((Y - mu_rc) ** 2, axis=1))
    per_point_sd = np.sqrt(var_rc.mean(axis=1))
    order = np.argsort(per_point_sd)  # most confident first
    curve = []
    for frac in fractions:
        keep = order[: max(1, int(round(frac * len(order))))]
        curve.append(float(np.sqrt(np.mean(per_point_err[keep] ** 2))))
    return np.array(curve)


FRACTIONS = np.linspace(0.15, 1.0, 18)
curve_id = risk_coverage(X_test, Y_test, FRACTIONS)
curve_in = risk_coverage(X_in, Y_in, FRACTIONS)
curve_dyn = risk_coverage(X_dyn, Y_dyn, FRACTIONS)
Out[18]:
Console
retain 0.25 | in-dist 0.0421 | input shift 0.0425 | dynamics shift 0.0424
retain 0.50 | in-dist 0.0430 | input shift 0.0424 | dynamics shift 0.0427
retain 1.00 | in-dist 0.0437 | input shift 0.0451 | dynamics shift 0.0459

At retained fraction 0.25, executed aggregate RMSEs are 0.0421, 0.0425, and 0.0424 for nominal, input-shifted, and damping-shifted data. At full coverage they are 0.0437, 0.0451, and 0.0459. The input-shift curve has 0.0424 at fraction 0.50, slightly below its 0.25 value: the empirical curve is not strictly monotone. These modest differences describe the chosen score and populations, not a universal benefit of abstention.

Three questions guide the reading. Does retaining fewer predictions reduce RMSE in the nominal population? Is that relationship preserved under input shift? Does it persist under dynamics shift? A favorable result in one regime cannot answer the other two. The score combines a residual variance floor with ensemble disagreement, while realized errors also contain unpredictable process noise. That noise can dominate errors and make a useful epistemic ranking difficult to see in aggregate RMSE.

Ranking and calibration are separate properties. A single positive rescaling preserves score rankings while potentially changing interval coverage. A percentile threshold targets a rejection fraction on its reference population, subject to finite-sample quantile and tie conventions, but a rolling window can normalize a sustained harmful shift into the new normal. An absolute threshold can encode a physically meaningful tolerance but needs validation under changed conditions. Neither form automatically survives distribution shift. Choose thresholds using retained-case risk, rejected-case outcomes, and the cost of deferral, not coverage alone.

Out[19]:
Visualization
Line chart of prediction error against retained fraction for three evaluation regimes.
Risk-coverage curves for the nominal, input-shifted, and damping-shifted test populations. Each point reports RMSE after retaining the lowest-variance fraction. Non-monotonic or small differences indicate that this uncertainty ranking is not a reliable selector of lower aggregate error in every regime.

Risk-Sensitive and Constrained Planning

Everything so far has been about the model. Now we change the objective. Minimizing expected cost is risk-neutral with respect to that cost: a certain cost of 10 and an equally likely cost of 0 or 20 have the same expectation. A point prediction does not force risk neutrality; a controller can still use conservative margins, robust constraints, or nonlinear costs. A predictive distribution permits explicit tail objectives, whose suitability depends on the application's losses and requirements.

Two families of alternatives recur in practice:

  • Risk-sensitive objectives use a tail-sensitive functional. For integrable costs, CVaR at level α\alpha averages the worst (1−α)(1-\alpha) probability mass, using fractional boundary mass when necessary. It approaches expected cost as α\alpha decreases to zero; near one it emphasizes increasingly extreme upper-tail costs. Here α=0.75\alpha=0.75 selects the worst quarter.
  • Constrained objectives separate feasibility requirements from performance. A pointwise constraint limits predicted states; a chance constraint such as Pr⁡[∣p∣>pmax⁡]≤δ\Pr[|p|>p_{\max}]\le\delta imposes a hard bound on violation probability, not a promise that every realized state stays inside. The objective can itself be expected cost or a risk-sensitive functional.

The expected-cost objective summarizes average loss; CVaR summarizes a chosen upper tail. These answer different questions. Expected cost can penalize rare failures heavily if their assigned loss is large, and a hard constraint can reject them without CVaR. Conversely, minimizing CVaR can still permit a safety violation if it improves the chosen tail objective. Specify safety requirements separately from performance preferences.

Where the distribution comes from

A risk-sensitive planner needs a distribution of outcomes for a candidate action sequence, not a single rollout. Several sources are available, and they are not equivalent:

  • Ensemble members as samples. Run the candidate sequence through each member. This uses variation between fitted models as an approximation to epistemic uncertainty, i.e., uncertainty about the model fit. Its quality depends on both coverage of plausible dynamics and the weights assigned to members; bootstrap ensembles do not by themselves guarantee calibrated tail probabilities. See Part VIII: Decision-Centric Research Lineages for the lineage.
  • Residual-noise propagation. Add sampled process noise to each member's predicted transition. In this toy plant, position and velocity noise innovations satisfy ϵt(p)=Δt ϵt(v)\epsilon_t^{(p)}=\Delta t\,\epsilon_t^{(v)}, so they must not be sampled independently. This relation concerns the noise contributions, not the full state increments. Combining model fits with this noise defines an approximate rollout distribution; one-step calibration does not validate its multi-step joint distribution.
  • Scenario evaluation. Evaluate explicit parameter and disturbance choices. A finite scenario set is inspectable, but checking it does not prove safety for unsampled conditions. Scenario-optimization guarantees require additional sampling, optimization, and support assumptions; merely naming a chosen schedule worst case supplies none.

The code combines alternative fitted models with residual-noise draws. It does not explicitly add total predictive variance on top of ensemble variation. Nevertheless, held-out residual variance can itself include fitting error, so this construction is not a decomposition into independent uncertainty sources. Its joint rollout spread requires validation; over- or under-dispersion and missed model families can affect planning.

Conditional value at risk

For a random cost CC with E∣C∣<∞\mathbb E|C|<\infty and α∈(0,1)\alpha\in(0,1), CVaR has the threshold representation

CVaRα(C)=min⁡τ{τ+11−αE[(C−τ)+]}.\mathrm{CVaR}_\alpha(C) = \min_{\tau} \left\{ \tau + \frac{1}{1-\alpha} \mathbb{E}\left[(C - \tau)^+ \right] \right\}.

where:

  • CC: the random episode cost
  • α∈(0,1)\alpha \in (0,1): the confidence level, where 1−α1-\alpha is the fraction of the distribution treated as the tail
  • τ\tau: a scalar threshold; an α\alpha-quantile is a minimizer, with possible non-uniqueness for discrete distributions
  • (C−τ)+=max⁡(C−τ,0)(C - \tau)^+ = \max(C - \tau, 0): the positive part, counting only how far each outcome exceeds the threshold

The expression includes the threshold itself plus scaled expected excess. It is convex in τ\tau but generally non-smooth. Pointwise convex loss in controller parameters together with a parameter-independent sampling law is sufficient for convexity of CVaR in those parameters, not necessary in every instance. Nonlinear learned rollouts provide no such guarantee automatically. See Rockafellar and Uryasev's general-loss treatment.

For equal-weight samples, sort costs in descending order, D1≥⋯≥DSD_1\ge\cdots\ge D_S. Let m=(1−α)Sm=(1-\alpha)S, q=⌊m⌋q=\lfloor m\rfloor, and r=m−qr=m-q. The exact empirical-distribution CVaR includes fractional boundary mass:

CVaR^α(C)=∑k=1qDk+rDq+1m.\widehat{\mathrm{CVaR}}_\alpha(C)=\frac{\sum_{k=1}^{q}D_k+rD_{q+1}}{m}.

When mm is an integer, average the largest mm values and omit the boundary term entirely. Otherwise the next value contributes fraction rr. This also avoids indexing past the array when floating-point arithmetic rounds mm to SS at a tiny positive α\alpha. Averaging the largest ⌈m⌉\lceil m\rceil values is generally an approximation. At S=18S=18 and α=0.75\alpha=0.75, four costs contribute fully and the fifth contributes half. Such a small effective tail sample is noisy; evaluate selected controllers independently of the search samples.

The running environment already includes Gaussian process noise. We now add a soft envelope penalty to stage cost when ∣p∣|p| exceeds a chosen threshold. This can amplify expensive excursions, but it is not a control barrier function or a hard constraint. Quadratic costs can also have substantial tails, and neither this penalty nor CVaR guarantees different actions or improved realized tail costs. The experiment measures those outcomes rather than assuming them.

The distinction between averaging across all outcomes and averaging across the worst outcomes is easiest to see directly on a cost distribution. The cell below draws a right-skewed set of synthetic episode costs and computes both summaries, using the same sample estimator the planner optimizes.

In[20]:
Code
# A synthetic right-skewed cost distribution, used only to make the CVaR idea
# concrete. The seed is fixed so the figure is reproducible.
CVAR_RNG = np.random.default_rng(31337)
COST_SAMPLES = CVAR_RNG.lognormal(mean=0.30, sigma=0.60, size=4000)

CVAR_ALPHA = 0.75
cost_expectation = float(COST_SAMPLES.mean())
cvar_k = max(1, int(np.ceil((1.0 - CVAR_ALPHA) * COST_SAMPLES.size)))
cost_cvar = float(np.sort(COST_SAMPLES)[-cvar_k:].mean())
cost_tail_cut = float(np.quantile(COST_SAMPLES, CVAR_ALPHA))
Out[21]:
Visualization
Histogram of synthetic costs with expectation and CVaR lines marked.
Histogram of a right-skewed synthetic distribution of episode costs, with the expectation and the conditional value at risk marked as vertical lines. The expectation lands near the bulk of the distribution, while the CVaR averages only the worst quarter and sits well into the tail that a safety objective cares about, making concrete how a tail-sensitive objective differs from an average one.
In[22]:
Code
SOFT_LIMIT = 0.50
BARRIER_WEIGHT = 5.0
W_P, W_V, W_U = 1.0, 0.05, 0.01

W_STACK = np.stack([W for W, _ in members])  # (M, F, 2)
M_MODELS, F_FEAT = W_STACK.shape[0], W_STACK.shape[1]


def plan(
    state,
    rng_,
    mode="mean",
    horizon=8,
    n_cand=128,
    n_draws=2,
    cvar_alpha=0.75,
    u_max=U_VALID,
    return_sequence=False,
):
    """Random-shooting MPC. Return the first action, or the selected sequence."""
    if mode not in ("mean", "cvar"):
        raise ValueError("mode must be mean or cvar")
    for name, value in (
        ("horizon", horizon),
        ("n_cand", n_cand),
        ("n_draws", n_draws),
    ):
        if (
            isinstance(value, (bool, np.bool_))
            or not isinstance(value, (int, np.integer))
            or value <= 0
        ):
            raise ValueError(f"{name} must be a positive integer")
    if not 0.0 < cvar_alpha < 1.0:
        raise ValueError(
            "CVaR confidence level must be strictly between zero and one"
        )
    S = M_MODELS * n_draws
    W_rows = W_STACK[np.arange(S) % M_MODELS]
    A = rng_.uniform(-u_max, u_max, size=(horizon, n_cand))
    s = np.tile(np.asarray(state, dtype=float), (S, n_cand, 1))
    costs = np.zeros((S, n_cand))
    for t in range(horizon):
        u = np.broadcast_to(A[t], (S, n_cand))
        Phi = features(
            np.column_stack([s[:, :, 0].ravel(), s[:, :, 1].ravel(), u.ravel()])
        )
        delta = np.einsum(
            "snf,sfj->snj", Phi.reshape(S, n_cand, F_FEAT), W_rows
        )
        # A velocity innovation also shifts position by DT times that innovation.
        noise_v = rng_.normal(0.0, np.sqrt(ALEATORIC[1]), size=(S, n_cand))
        delta = delta + np.stack([DT * noise_v, noise_v], axis=-1)
        s = s + delta
        over = np.maximum(0.0, np.abs(s[:, :, 0]) - SOFT_LIMIT)
        costs += (
            W_P * s[:, :, 0] ** 2
            + W_V * s[:, :, 1] ** 2
            + W_U * A[t] ** 2
            + BARRIER_WEIGHT * over**2
        )
    if mode == "mean":
        score = costs.mean(axis=0)
    else:
        tail_mass = (1.0 - cvar_alpha) * S
        whole = int(np.floor(tail_mass))
        fraction = tail_mass - whole
        descending = np.sort(costs, axis=0)[::-1]
        tail_sum = descending[:whole].sum(axis=0)
        if fraction > 0.0 and whole < S:
            tail_sum += fraction * descending[whole]
        score = tail_sum / tail_mass
    selected = A[:, int(np.argmin(score))]
    return selected.copy() if return_sequence else float(selected[0])

At the same state and RNG state, both modes generate the same candidate actions, member assignments, and linked noise draws. Only their cost aggregation changes: mean versus empirical CVaR. Once actions differ, closed-loop states and therefore imagined futures differ too. Using the same seeds pairs the search randomness; it does not make later state-conditioned predictions identical.

Running the comparison

We pair environment noise streams and search seeds across controllers. Different actions still lead to different states. Common random numbers reduce the variance of a difference when paired outcomes have positive covariance: Var(A−B)=Var(A)+Var(B)−2Cov(A,B)\mathrm{Var}(A-B)=\mathrm{Var}(A)+\mathrm{Var}(B)-2\mathrm{Cov}(A,B). They can increase it when covariance is negative. This experiment reports a paired finite comparison, not a measured reduction in required episode count; see Glasserman and Yao's conditions for common random numbers.

In[23]:
Code
N_EPISODES = 30
N_STEPS = 32

env_sampler = np.random.default_rng(4242)
EPISODES = []
for _ in range(N_EPISODES):
    params = sample_nominal(env_sampler)
    p0 = env_sampler.uniform(0.20, 0.40) * env_sampler.choice([-1.0, 1.0])
    v0 = np.sign(p0) * env_sampler.uniform(0.20, 0.50)
    EPISODES.append((params, (p0, v0)))


def run_regulation(mode, params, s0, steps, env_rng, plan_rng):
    p, v = float(s0[0]), float(s0[1])
    total = 0.0
    for _ in range(steps):
        u = plan(np.array([p, v]), plan_rng, mode=mode)
        p, v = simulate(params, (p, v), [u], env_rng, SIGMA_P)[1]
        over = max(0.0, abs(p) - SOFT_LIMIT)
        total += (
            W_P * p * p
            + W_V * v * v
            + W_U * u * u
            + BARRIER_WEIGHT * over * over
        )
    return total


def evaluate_planner(mode, plan_seed):
    env_rng = np.random.default_rng(9999)
    plan_rng = np.random.default_rng(plan_seed)
    return np.array(
        [
            run_regulation(mode, prm, s0, N_STEPS, env_rng, plan_rng)
            for prm, s0 in EPISODES
        ]
    )


costs_mean = evaluate_planner("mean", 101)
costs_cvar = evaluate_planner("cvar", 101)
Out[24]:
Console
planner                mean cost    median   90th pct       max
expected-cost MPC          1.600     1.670      2.372     2.593
CVaR MPC                   1.673     1.727      2.463     2.853

The table gives mean, median, 90th percentile, and maximum episode cost for each planner across thirty paired episodes. These are empirical cost summaries, not population tail guarantees or safety outcomes. Reporting several summaries makes tradeoffs visible that a mean alone would hide.

The executed CVaR controller has mean cost 1.673, 90th-percentile cost 2.463, and maximum 2.853, compared with 1.600, 2.372, and 2.593 for expected-cost MPC. It does not improve the realized upper tail on these thirty episodes. Model mismatch, finite candidate search, and tail-sampling error are possible explanations, not isolated causal findings. Tune α\alpha, sample counts, and cost weights on separate validation populations, then evaluate across independent seeds. Optimizing an imagined tail is not the same as improving the actual one.

Out[25]:
Visualization
Grouped bars comparing empirical cost summaries for two controllers.
Mean, median, and 90th-percentile costs across thirty paired episodes. CVaR has higher realized costs on these summaries; optimizing a modeled tail does not guarantee improving actual outcomes.
Out[26]:
Visualization
Episode-cost histograms comparing expected-cost and CVaR controllers.
Episode-cost histograms for thirty paired initial conditions and process-noise streams. This sample shows no CVaR tail improvement, and one seed design is not a population guarantee.

From soft penalties to explicit constraints

A finite soft penalty permits a tradeoff: a violating candidate can win if reductions in other cost terms outweigh its extra penalty. Our zero-target position cost does not make ∣p∣=0.7|p|=0.7 intrinsically better than ∣p∣=0.55|p|=0.55. Nor does this example establish that no penalty choice could enforce feasibility in any problem: some bounded or finite problems admit exact penalties. The quadratic penalty here supplies no proof of pointwise safety. If an envelope is a requirement, encode and justify it separately rather than treating a tuned weight as a certificate.

The alternatives each have costs:

  • Hard constraints on a nominal rollout enforce ∣p∣≤pmax⁡|p|\le p_{\max} at modeled steps. Feasibility establishes that modeled property, not automatically the real plant's safety. Connecting the two requires justified error and disturbance assumptions.
  • Chance constraints enforce Pr⁡[∣p∣>pmax⁡]≤δ\Pr[|p|>p_{\max}]\le\delta under a specified probability law. Computing or estimating that probability can be difficult. Finite samples can under- or overestimate it, and rare violations can be absent from the sample. A bound on violation probability is a different requirement from all-disturbance invariance.
  • Runtime assurance. Check a proposal against an independently justified safety condition and select a verified fallback when necessary. Guarantees require a sound decision module, an appropriate initial condition, and a fallback with proven properties under the declared plant and disturbance assumptions. Our next experiment illustrates the switching architecture without satisfying those proof obligations.

The code above implements a soft penalty and sampled CVaR, not constrained optimization. Zero observed violations among eighteen scenarios does not imply zero violation probability. A predictive safety filter needs justified uncertainty bounds and recursive feasibility; see Wabersich and Zeilinger. For continuous control-affine dynamics s˙=f(s)+g(s)a\dot s=f(s)+g(s)a, a continuously differentiable barrier hh describes h(s)≥0h(s)\ge0 and constrains actions through h˙(s,a)+γ(h(s))≥0\dot h(s,a)+\gamma(h(s))\ge0, with an appropriate extended class-K function γ\gamma (for example γ(r)=λr\gamma(r)=\lambda r, λ>0\lambda>0). Forward invariance additionally needs regularity, feasible actions, correct dynamics or robust error bounds, an admissible initial state, and satisfaction over time. Sampling and delay require additional treatment. Our quadratic penalty is not this condition; see Ames et al..

Fallback Policies and Runtime Assurance

The architecture we now describe has a history in safety-critical control under various names; a standard formulation is the Simplex pattern. Split the controller into two parts:

  • A performance controller, which may be arbitrarily complex, learned, and unverified. It is the thing that does the task well.
  • A safety controller (the fallback), which is simple enough to be analyzed, tested exhaustively, or verified. It is the thing that keeps the system inside the safe set.

A monitor decides which controller is in charge. A formal safety result needs more than a simple controller and a threshold: the initial condition must admit a safe continuation, the decision module must be correct, and switching must preserve that continuation under the declared disturbances and timing. A permanently safe command sequence and correct decision module are explicit assumptions in the Black-Box Simplex formulation. This relocates the verification burden rather than removing it. A finite simulation of our proportional-derivative fallback is evidence about tested episodes, not certification.

The sections that follow name the important parameters: the declared parameter bounds (k,c,g,f)(k, c, g, f), the innovation threshold ZthresholdZ_{\text{threshold}}, the fallback tail length, and the CVaR level together with the number of draws and candidates.

The recoverable set

Given a fallback controller πfb\pi_{\text{fb}}, we want to know which initial states are guaranteed to stay safe under that controller. The recoverable set R\mathcal{R} answers this question: it is the set of states from which πfb\pi_{\text{fb}} can keep the system inside the safe set forever:

R={x: ∀t≥0, pt∈Psafe when x0=x and ut=πfb(xt)}.\mathcal{R} = \left\{ x : \ \forall t \ge 0, \ p_t \in \mathcal{P}_{\text{safe}} \ \text{when} \ x_0 = x \ \text{and} \ u_{t} = \pi_{\text{fb}}(x_t) \right\}.

where:

  • xx: the initial state (p0,v0)(p_0, v_0)
  • pt,vtp_t, v_t: position and velocity under the fallback rollout
  • Psafe\mathcal{P}_{\text{safe}}: the safe position set (in our example, {p:∣p∣≤Plimit}\{p : |p| \le P_{\text{limit}}\})
  • πfb\pi_{\text{fb}}: the fallback (safety) controller

For a specified plant, disturbance model, and fallback, the monitor should preserve membership in R\mathcal{R}. Under uncertainty the definition must quantify over all allowed disturbances and parameter trajectories as well. Outside this particular fallback's recoverable set, its guarantee is unavailable; another controller may still succeed. Checking a finite horizon or a finite grid does not establish the infinite-time property in the definition. The set is useful precisely because it states which continuation and assumptions must be proved.

For our bounded position envelope and nondegenerate Gaussian velocity innovation, an all-draw robust recoverable set is empty: from any finite state and applied command, an innovation large enough to leave the envelope in the next step has positive probability. A nonempty deterministic invariant argument would require a different, explicitly bounded disturbance model. A probabilistic safety criterion is another option, but is not the all-disturbance definition above.

A monitor that does not use the learned model

The monitor should not assume that accurate training predictions imply a valid safety model under shift. Learned models can participate in a certificate when their error is bounded under defensible assumptions; otherwise a model-only safety check can share the planner's blind spots. Our small damping shift changed empirical coverage modestly, not catastrophically, but that result still does not establish future safety. The toy screen below uses a separate declared scenario to illustrate this separation of responsibilities.

We can posit a separate parameter scenario to illustrate the architecture. Its settings are not independently verified bounds on a physical plant, so calling the scenario conservative does not establish containment. In our example, we use

k≥0.80,c≥0.0,g≥0.90,f≥0,k \ge 0.80, \qquad c \ge 0.0, \qquad g \ge 0.90, \qquad f \ge 0,

where:

  • kk: stiffness, the restoring force per unit displacement
  • cc: viscous damping, the drag proportional to velocity
  • gg: actuator gain, how much acceleration one unit of command produces
  • ff: dry friction, a constant drag opposing motion

These lower limits alone do not define a finite parameter box or prove a worst-case trajectory. Our code fixes stiffness and friction and selects damping/gain endpoints greedily along one trajectory. Its choices are illustrative engineering settings, not verified physical bounds. Greedy outward acceleration need not maximize the eventual peak because it changes subsequent velocity and position. The resulting value is a scenario peak, not an enclosure of all reachable states.

For a bounded uncertainty set B\mathcal B, a disturbance set D\mathcal D, and known initial state ss, a formal endpoint reachable set would be

WH(s)={pH: s0=s, st+1=fbt(st,πfb(st),dt), bt∈B, dt∈D},\mathcal W_H(s)=\{p_H:\ s_0=s,\ s_{t+1}=f_{b_t}(s_t,\pi_{\mathrm{fb}}(s_t),d_t),\ b_t\in\mathcal B,\ d_t\in\mathcal D\},

where:

  • WH(x)\mathcal{W}_H(x): the worst-case forward reachable set of positions after HH steps from state xx
  • pHp_H: the position at the end of the horizon
  • pt,vtp_t, v_t: position and velocity at step tt of the rollout
  • (k,c,g,f)(k, c, g, f): stiffness, damping, gain, and friction parameter tuple
  • B\mathcal{B}: the box of declared parameter bounds
  • ut=πfb(pt,vt)u_t = \pi_{\text{fb}}(p_t, v_t): the fallback action at step tt

Here fbtf_{b_t} is the discrete plant transition with parameter tuple btb_t and disturbance dtd_t; D\mathcal D states the permitted disturbances. Safety over the horizon needs every intermediate reachable set, not only WH\mathcal W_H, to lie inside the safe envelope. State-estimation error requires an initial set rather than a single ss. Infinite-time assurance additionally needs a terminal invariant continuation, and a physical plant requires inter-sample and latency analysis. Gaussian process noise in our simulation is unbounded, so no finite deterministic disturbance envelope covers every draw. A bounded-noise proof would describe a different explicit disturbance model, or a probabilistic result with a stated failure level. None is proved by the scenario screen below.

A second monitor: prediction consistency

A complementary signal comes from innovation monitoring: compare predicted and observed next states after an action, normalized by predictive standard deviation. This is a reactive diagnostic, not a reachable-set certificate. The toy combines it with a finite scenario screen; a certified architecture would replace that screen with a justified enclosure and safe continuation test.

zt=1d∑j=1d(st+1(j)−s^t+1(j))2σt+12 (j).z_t = \sqrt{\frac{1}{d}\sum_{j=1}^{d} \frac{\left( s_{t+1}^{(j)} - \hat{s}_{t+1}^{(j)} \right)^2}{\sigma_{t+1}^{2\,(j)}}} .

where:

  • dd: the state dimension (here d=2d = 2, for position and velocity)
  • st+1(j)s_{t+1}^{(j)}: the jj-th component of the observed next state
  • s^t+1(j)\hat{s}_{t+1}^{(j)}: the jj-th component of the model's predicted mean next state
  • σt+12 (j)\sigma_{t+1}^{2\,(j)}: the model's predicted variance for the jj-th component
  • ztz_t: the root-mean-square standardized innovation; under correct componentwise zero means and variances, E[zt2]=1\mathbb E[z_t^2]=1, not generally E[zt]=1\mathbb E[z_t]=1

The threshold ZthresholdZ_{\text{threshold}} trades sensitivity against false alarms. Changed dynamics may increase the statistic, but benign noise can also trigger it and small changes may remain undetected. This root statistic is not a chi-square innovation test: the position and velocity innovations are correlated by construction, and we do not whiten their full covariance. The code uses a single-transition threshold and a lower reset threshold, not a smoothed window, formal sequential detector, or calibrated false-alarm probability. The earlier filtering material provides context for residual diagnostics, not a guarantee for this heuristic.

Two details matter. The statistic needs an observed transition, so it cannot prevent harm during that first transition. Its scale also depends on the predictive variance, so overestimated variance can hide errors and underestimated variance can create false alarms. There is no universal one- or two-step detection bound here. Logging innovations, threshold crossings, and selected actions is necessary to distinguish detection from scenario-triggered handover.

The two monitors are complementary:

  • The consistency monitor is inexpensive and reactive. It needs at least one observed transition and may miss a shift indefinitely.
  • A certified reachability monitor can reject a proposal before execution, provided its enclosure and timing assumptions hold. Our finite scenario screen has this prospective timing but not that certification.

A design may combine these signals, but combining heuristics does not produce a proof. A sound safety condition must remain enforced even when residuals look benign; an innovation test should not override it merely to improve performance. The two signals address different evidence, and their switching logic needs explicit analysis. In the experiment the shield checks the candidate sequence followed by a fallback tail in one deterministic scenario and also responds to innovation threshold crossings.

In[27]:
Code
# Illustrative parameters for one adverse deterministic schedule, not an
# exhaustive reachable-set calculation or a validated physical envelope.
K_MIN, C_MIN, GAIN_MIN, GAIN_MAX, FRIC_MIN = 0.80, 0.00, 0.90, 1.10, 0.0
C_MAX = 1.10
Z_THRESHOLD = 2.5
FALLBACK_TAIL = 30


def fallback_action(state, kp=3.5, kd=3.0, u_max=U_FALLBACK):
    """Proportional-derivative regulator. No learned component, full authority."""
    p, v = float(state[0]), float(state[1])
    return float(np.clip(-kp * p - kd * v, -u_max, u_max))


def scenario_peak(state, action_fn, steps):
    """Peak along one adverse deterministic parameter schedule, not a bound."""
    p, v = float(state[0]), float(state[1])
    peak = abs(p)
    for t in range(steps):
        u = float(action_fn(t, (p, v)))
        c = C_MIN if p * v > 0.0 else C_MAX  # least damping when escaping
        gain = GAIN_MAX if np.sign(p) * u > 0.0 else GAIN_MIN  # weakest brake
        acc = -K_MIN * p - c * v - FRIC_MIN * np.sign(v) + gain * u
        v = v + DT * acc
        p = p + DT * v
        peak = max(peak, abs(p))
    return peak


def monitor_ok(state, planned_actions, horizon):
    """Finite-horizon scenario screen; omits process noise and is not certified."""

    def policy(t, s):
        return planned_actions[t] if t < horizon else fallback_action(s)

    return scenario_peak(state, policy, horizon + FALLBACK_TAIL) <= P_LIMIT

The declared scenario includes zero damping, but it omits process noise and searches neither all parameter paths nor all disturbance paths. It can trigger unnecessary fallback or accept an unsafe plan. The finite tail is thirty steps, not an invariant-set proof. The shield remains useful as an executable architecture example precisely when these limitations are visible; calling it conservative is not evidence that its peak upper-bounds the actual plant.

Running the shield

We run expected-cost MPC, CVaR MPC, and shielded MPC in nominal and zero-damping environments. Episodes start near the envelope with outward velocity, creating a useful braking stress test. Braking is not necessary in every nominal episode: restoring force and damping can already keep some noise-free trajectories inside. Comparing these controllers measures their combined behavior; it does not isolate manual braking from luck or prove which component caused an improvement.

In[28]:
Code
def run_controller(mode, params, s0, n_steps, env_rng, plan_rng, horizon=8):
    p, v = float(s0[0]), float(s0[1])
    prev_state, prev_action = np.array([p, v]), 0.0
    distrust = False
    violations, shield_steps, cost = 0, 0, 0.0
    trace = np.empty((n_steps + 1, 2))
    shield_flag = np.zeros(n_steps + 1, dtype=bool)
    trace[0] = (p, v)
    for t in range(n_steps):
        state = np.array([p, v])
        if mode == "shield" and t > 0:
            mu_i, var_i, _ = predict(
                members, np.array([[prev_state[0], prev_state[1], prev_action]])
            )
            z = float(
                np.sqrt(np.mean((state - prev_state - mu_i[0]) ** 2 / var_i[0]))
            )
            if z > Z_THRESHOLD:
                distrust = True
            elif (
                z < 1.0
                and scenario_peak(
                    state, lambda tt, s: fallback_action(s), FALLBACK_TAIL
                )
                <= P_LIMIT
            ):
                distrust = False
        planned = plan(
            state,
            plan_rng,
            mode=("cvar" if mode == "cvar" else "mean"),
            horizon=horizon,
            return_sequence=True,
        )
        use_fallback = mode == "shield" and (
            distrust or not monitor_ok(state, planned, horizon)
        )
        if use_fallback:
            u = fallback_action(state)
            shield_steps += 1
        else:
            u = float(planned[0])
        prev_state, prev_action = state, u
        p, v = simulate(params, (p, v), [u], env_rng, SIGMA_P)[1]
        over = max(0.0, abs(p) - SOFT_LIMIT)
        cost += (
            W_P * p * p + W_V * v * v + W_U * u * u + BARRIER_WEIGHT * over**2
        )
        if abs(p) > P_LIMIT:
            violations += 1
        trace[t + 1] = (p, v)
        shield_flag[t] = use_fallback
    return dict(
        cost=cost,
        violations=violations,
        shield_steps=shield_steps,
        trace=trace,
        shield=shield_flag,
    )


NOMINAL_ENV = (1.00, 0.90, 1.00, 0.0)
SHIFTED_ENV = (1.00, 0.00, 1.00, 0.0)
RT_STEPS = 45
RT_MODES = ["mean", "cvar", "shield"]

rt_rng = np.random.default_rng(77)
EPISODES_RT = []
for _ in range(6):
    p0 = rt_rng.uniform(0.25, 0.32) * rt_rng.choice([-1.0, 1.0])
    v0 = np.sign(p0) * rt_rng.uniform(0.78, 0.88)
    EPISODES_RT.append((p0, v0))

rt_rows, rt_traces = [], {}
for env_name, params in (
    ("nominal", NOMINAL_ENV),
    ("damping shift", SHIFTED_ENV),
):
    for mode in RT_MODES:
        plan_rng = np.random.default_rng(1000)
        for i, (p0, v0) in enumerate(EPISODES_RT):
            env_rng = np.random.default_rng(500 + i)
            res = run_controller(
                mode, params, (p0, v0), RT_STEPS, env_rng, plan_rng
            )
            rt_rows.append(
                dict(
                    env=env_name,
                    mode=mode,
                    cost=res["cost"],
                    viol=res["violations"],
                    peak=float(np.abs(res["trace"][:, 0]).max()),
                    shield=res["shield_steps"],
                )
            )
            if i == 0:
                rt_traces[(env_name, mode)] = res
Out[29]:
Console
environment    controller violations  peak |p|  mean cost  shield steps
nominal        mean               7     0.650       4.38           0.0
nominal        cvar               7     0.667       4.53           0.0
nominal        shield             0     0.541       2.99           2.0
damping shift  mean             103     0.947      18.65           0.0
damping shift  cvar             115     0.972      20.42           0.0
damping shift  shield             2     0.604       5.52           6.3

The printed rows pair each environment with each controller and report the total violations, the worst peak position, the mean episode cost, and the mean number of shield steps.

Three columns carry the message. The violations column counts the number of simulated steps in which ∣p∣|p| exceeded the limit. The peak column is how far the worst episode got. The shield steps column counts how often the fallback was in charge. Together they describe the answer to the central question of the chapter: what did the system lose by trusting the learned model, and what did it pay to stop trusting it?

The executed nominal rows have 7, 7, and 0 violating steps for expected-cost, CVaR, and shielded MPC; shifted rows have 103, 115, and 2. The shifted shield's maximum position magnitude is 0.604, above the 0.600 limit. It reduces violations here but does not eliminate them. The shield also has stronger actuator authority than the learned planners, so the reduction cannot be attributed solely to monitoring. Matched seeds pair process noise, not later states. These finite results establish neither invariant safety nor a detection delay.

The shield steps column measures fallback usage. Frequent nominal handovers can reflect the aggressive initial conditions, scenario screen, or innovation false alarms; they do not prove a threshold is too tight. No handover under a shift is not necessarily an error if the proposed actions remain safe. Tune performance only after defining and validating the safety condition, and keep separate counts for each trigger to understand its contribution. The present table is an empirical tradeoff summary, not a soundness test.

Out[30]:
Visualization
Line chart of position over time for two controllers with a horizontal safety limit.
First damping-shifted episode for expected-cost MPC and its heuristic shield. Dashed lines mark both sides of the envelope; vertical shading marks fallback action intervals. The fallback has stronger actuator authority, so this comparison does not isolate the value of monitoring or prove safety.

Vertical shading marks action intervals during which the fallback was selected, including prospective scenario rejections and reactive innovation triggers. This first shifted episode illustrates the combined controller, not a certified recoverable set or a universal detection time. To isolate mechanisms, an additional ablation would compare each trigger separately and equalize actuator authority. Safety analysis must also account for untested initial states, longer horizons, process-noise tails, estimation errors, and switching delays.

Out[31]:
Visualization
Grouped bar chart of simulated envelope-violation counts by controller and environment.
Counts of simulated time steps outside the position envelope for three controllers and two environments. Counts aggregate six seeded episodes per controller. They measure these rollouts, not all possible trajectories or a formal safety guarantee.

Limitations and Impact

The techniques in this chapter are useful, but each has limitations, and a clear account of those limitations matters more than promotion. You may be worse off if you adopt these tools without understanding their limits, because the tools can create false confidence.

Stress testing depends on the mechanisms and ranges considered. Unanticipated couplings, sensor degradation, and rare conditions can remain outside the suite. Revisit it when incidents, near-misses, operating conditions, or requirements reveal missing cases. A calendar year without growth is not by itself evidence that the suite is inadequate; relevance and coverage matter more than its size.

Calibration is a property of specified forecasts and a population. Deployment may sample a population, but adaptive closed-loop data can be dependent and differ from the calibration population. Marginal coverage can hide subgroup failures, including safety-relevant states that need not be the least densely trained region. Broad aggregate diagrams provide some information, not a substitute for conditional and trajectory-level evaluation; see Angelopoulos and Bates.

Ensemble disagreement is not a general-purpose dynamics-shift detector. The tested input shift increased disagreement more than the damping shift, while aggregate one-step error changed little. This is evidence about one fitted ensemble and two populations, not a universal detector comparison. Resampling old data cannot directly observe a newly changed transition law, although changed trajectories may later create unfamiliar inputs that increase disagreement. Different ensemble constructions or shift sizes can behave differently; test the signal rather than treating it as a certificate.

Risk-sensitive planning depends on its rollout distribution. CVaR over a mis-specified model can optimize the wrong tail. With M=9M=9 fits and two draws per fit, S=18S=18 samples represent each candidate; at α=0.75\alpha=0.75 the empirical upper-tail mass is 4.54.5 samples. Four contribute fully and the fifth half. These samples share fitted members and are not eighteen independent estimates of model uncertainty. Increase draws and evaluate across independent training/evaluation seeds to assess sensitivity. Moderate risk levels reduce tail-estimation demands but do not establish safety.

Runtime assurance makes proof obligations explicit. Empirically tuned fallback behavior and convenient parameter ranges do not constitute a safety proof. An independent safety argument can avoid verifying the learned performance model, but still needs justified plant and disturbance assumptions, initial admissibility, sound decision logic, fallback continuation, and timing. That burden may be difficult rather than small or local. Simplex research illustrates aircraft and robot applications; earlier work also discusses pacemaker models. These are research examples, not a claim of universal certified deployment.

Model errors can affect both performance and safety; their consequences depend on the plant and control architecture. The experiments here distinguish empirical prediction diagnostics, tail-sensitive objectives, and fallback switching. They do not show that consequences are bounded in every deployment. Moving from this prototype to a safety argument requires explicit uncertainty/disturbance bounds or probability levels, validation of the fallback and decision module, and timing and state-estimation analysis. The value of the separation is that these obligations can be stated independently of a leaderboard prediction score.

In Part XI: Evaluation and Understanding, prediction was an evaluation axis. This chapter adds closed-loop consequences under stated stresses. These rankings need not agree. Deployment readiness depends on actual requirements and validated operating conditions, not on a prediction score alone or a universal ranking of survival.

Summary

  • Input shift changes the input distribution; dynamics shift changes the transition law and may also change subsequently visited inputs. Observation and task changes can co-occur with either, while hidden parameters provide one mechanism of dynamics change.
  • Stress testing combines falsification and measurement. Specify plausible interventions, measure their consequences, and record what remains untested. Passing a finite suite does not establish safety outside it.
  • Ensemble disagreement is a diagnostic, not a detector guarantee. Our input-shift population produced a larger disagreement increase than the damping-shift population. Aggregate one-step error alone concealed some of the population difference.
  • Calibration means coverage matches its specified claim. Reliability diagrams measure empirical marginal coverage at displayed levels; report interval width too. Global rescaling may help a new population but does not automatically transfer. Conformal guarantees need exchangeability and do not imply conditional or trajectory safety.
  • Abstention needs an informative ranking and a defined alternative. Risk-coverage curves evaluate retained predictions, not the safety of rejected cases. Neither percentile nor absolute thresholds universally survive shift.
  • For integrable costs, CVaR emphasizes a selected upper tail and approaches the mean as its level decreases to zero. Near one it emphasizes extremes. Model mismatch and finite tail samples can prevent realized improvements; a soft penalty is not a hard safety constraint.
  • Runtime assurance splits the controller. A performance controller does the task; a certified fallback maintains the invariant; a monitor decides who is in charge and must trigger before the state leaves the recoverable set.
  • Separate diagnosis from certification. Innovation checks are reactive diagnostics. A prospective reachable-set check carries a guarantee only with sound enclosures, disturbance and timing assumptions, and a safe continuation. Our scenario-based shield demonstrates switching, not that proof.

The next chapter, Security, Ethics, Privacy, and Governance, extends the scope to adversarial, societal, and institutional concerns. The distinction developed here remains useful: a diagnostic can inform a decision without proving that the decision satisfies an engineering or governance requirement.

Quiz

Ready to test your understanding? Take this quick quiz to reinforce what you've learned about robustness, calibration, and safe control.

Robustness, Calibration, and Safe Control

Question 1 of 80 of 8 completed
What distinguishes dynamics shift (concept shift) from pure input shift?

Comments

No comments yet. Be the first to share your thoughts!

Reference

Citation details

Cite or share this article.

BIBTEXAcademic
@misc{brenndoerfer2026robustnesscalibration, author = {Michael Brenndoerfer}, title = {Robustness, Calibration, and Safe Control}, year = {2026}, url = {https://mbrenndoerfer.com/writing/robustness-calibration-safe-control-world-models}, organization = {mbrenndoerfer.com}, note = {Accessed: 2026-10-11} }
APAAcademic
Michael Brenndoerfer (2026). Robustness, Calibration, and Safe Control. Retrieved from https://mbrenndoerfer.com/writing/robustness-calibration-safe-control-world-models
MLAAcademic
Michael Brenndoerfer. "Robustness, Calibration, and Safe Control." 2026. Web. October 11, 2026. <https://mbrenndoerfer.com/writing/robustness-calibration-safe-control-world-models>.
CHICAGOAcademic
Michael Brenndoerfer. "Robustness, Calibration, and Safe Control." Accessed October 11, 2026. https://mbrenndoerfer.com/writing/robustness-calibration-safe-control-world-models.
HARVARDAcademic
Michael Brenndoerfer (2026) 'Robustness, Calibration, and Safe Control'. Available at: https://mbrenndoerfer.com/writing/robustness-calibration-safe-control-world-models (Accessed: October 11, 2026).
SimpleBasic
Michael Brenndoerfer (2026). Robustness, Calibration, and Safe Control. https://mbrenndoerfer.com/writing/robustness-calibration-safe-control-world-models

About the author

Continue with the full handbook

This chapter is part of World Models Handbook. Use the handbook page to browse the complete table of contents and continue reading in sequence.

Explore World Models Handbook
Newsletter

Stay up to date

Get articles, book updates, and news delivered to your inbox.

No spam, unsubscribe anytime.

or

Join the community

Sign in to remove popups, track your reading progress, and join the discussion.