Michael Brenndoerfer | Data & AI, Private Equity, Technology

Michael Brenndoerfer

Sim-to-Real Deployment and Continual Correction › Four Candidate Models

Michael BrenndoerferAugust 10, 202612 min read

World Models HandbookAll chapters

Part 12 of 12Reliable World Models

Chapter 71Sim-to-Real Deployment and Continual Correction › Four Candidate Models

Reading section 5 of 13 · Four Candidate Models · 12 minSection 5 of 13 · 12 min

Section 5 of 13 · 12 min read

Sim-to-Real Deployment and Continual Correction › Four Candidate ModelsOne section · Sections in this chapter5 / 13 ↓

Sections in this chapter13

Reading tools Definitions & inline links

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

We compare four predictors. Each receives the same randomized training data where applicable, and each is used to drive the same controller, so differences in closed-loop cost come from the predictor and not from the planner. Holding the controller fixed is what makes the comparison a test of the predictor; if the controller changed with the model, we could not tell which change produced the difference.

Nominal sim model: the hand-written simulator, with $k=1.0$ , $c=0.4$ , $b=1.0$ , and clipping at $u_{\max}$ . It is exactly the model you would ship if you trusted your simulator.
Domain-randomized, parameter-blind: a single predictor fitted to the pooled transitions from all 24 plants, with no context input.
Domain-randomized, parameter-conditioned: the same pooled data, but the predictor also receives the plant parameters. At deployment we feed it parameters identified from the adaptation episode.
Nominal plus residual: the nominal simulator plus an explicitly shaped residual estimator fitted to the adaptation episode only.

The predictor class is deliberately modest: a quadratic feature map over $(s_t, u_t)$ appended with the ones vector, followed by a linear solve. Concretely, the feature map $\phi(s_t, u_t) \in \mathbb{R}^{10}$ is:

\phi(s_t, u_t) = \bigl[\,p_t,\; v_t,\; u_t,\; p_t v_t,\; p_t u_t,\; v_t u_t,\; p_t^2,\; v_t^2,\; u_t^2,\; 1\,\bigr]^\top

The feature map transforms the raw state and action into a fixed set of basis functions, allowing the linear predictor to capture quadratic nonlinearities. The feature map is a fixed, human-chosen transformation; the learning happens in the linear weights applied to it, so the expressive power of the predictor is exactly the span of these ten basis functions. Each element of the feature vector is defined below.

where:

$p_t$ : the position component of the state
$v_t$ : the velocity component of the state
$u_t$ : the issued scalar command
the bilinear terms $p_t v_t, p_t u_t, v_t u_t$ : couplings between state components ( $p_t v_t$ ), and between state and action ( $p_t u_t, v_t u_t$ )
the quadratic terms $p_t^2, v_t^2, u_t^2$ : curvature in each variable
$1$ : the constant bias term

The base predictor is $\hat{s}_{t+1}=W^\top\phi(s_t,u_t)$ with $W\in\mathbb{R}^{10\times2}$ : ten learned weights for each state component. The conditioned extension appends $k,c,b,kp,cv,bu$ , giving sixteen features and a $16\times2$ weight matrix. These interactions represent the per-plant semi-implicit transition exactly within the tested command range. The small explicit class removes neural optimization complexity, but does not isolate training-distribution quality from context, feature capacity or adaptation information. A neural predictor could implement the same batched interface and evaluation workflow; its fit and numerical comparisons would need fresh validation.

In[10]:

Code

def build_pairs(data, context=False):
    """Return (S, U, S_next), or (S, U, C, S_next) when context=True."""
    S, U, S_next, C = [], [], [], []
    for d in data:
        S.append(d["s"][:-1])
        U.append(d["u"])
        S_next.append(d["s"][1:])
        if context:
            p = d["params"]
            C.append(np.tile([p["k"], p["c"], p["b"]], (len(d["u"]), 1)))
    S, U, S_next = np.vstack(S), np.concatenate(U), np.vstack(S_next)
    if context:
        return S, U, np.vstack(C), S_next
    return S, U, S_next


def feat_base(S, U):
    p, v = S[:, 0], S[:, 1]
    u = np.asarray(U, dtype=float).reshape(-1)
    return np.column_stack(
        [p, v, u, p * v, p * u, v * u, p * p, v * v, u * u, np.ones_like(p)]
    )


def feat_cond(S, U, C):
    base = feat_base(S, U)
    p, v = S[:, 0], S[:, 1]
    u = np.asarray(U, dtype=float).reshape(-1)
    k, c, b = C[:, 0], C[:, 1], C[:, 2]
    return np.column_stack([base, k, c, b, k * p, c * v, b * u])


def fit_linear(X, Y, ridge=1e-8):
    A = X.T @ X + ridge * np.eye(X.shape[1])
    return np.linalg.solve(A, X.T @ Y)


S_tr, U_tr, S_next_tr = build_pairs(train_data)
W_blind = fit_linear(feat_base(S_tr, U_tr), S_next_tr)

S_trc, U_trc, C_tr, S_next_trc = build_pairs(train_data, context=True)
W_cond = fit_linear(feat_cond(S_trc, U_trc, C_tr), S_next_trc)

The predictors themselves expose the same batched interface predict_batch(S, U), so the planner uses the same code for every candidate. In this noiseless linear plant, the identified parameter-conditioned model and the target-fitted residual model recover essentially the same dynamics. Their predictions differ only at numerical fitting precision, and their tracking costs agree to the displayed precision. These are two fitting routes, not two empirically distinct controllers in this experiment; the comparison cannot establish that one route is better.

In[11]:

Code

def nominal_predict_batch(S, U):
    """The simulator: knows the nominal dynamics and its own actuator limit."""
    p, v = S[:, 0], S[:, 1]
    u = np.clip(np.asarray(U, dtype=float).reshape(-1), -U_MAX, U_MAX)
    acc = -NOMINAL["k"] * p - NOMINAL["c"] * v + NOMINAL["b"] * u
    v_next = v + DT * acc
    return np.column_stack([p + DT * v_next, v_next])


def blind_predict_batch(S, U):
    """Pooled fit: one compromise predictor for the whole randomized family."""
    return feat_base(S, U) @ W_blind


def make_cond_batch(ctx):
    """Context-conditioned fit; at deployment the context must be identified."""
    ctx = np.asarray(ctx, dtype=float).reshape(1, -1)

    def predict_batch(S, U):
        C = np.tile(ctx, (len(S), 1))
        return feat_cond(S, U, C) @ W_cond

    return predict_batch


def make_residual_batch(W_res):
    """Nominal dynamics plus an explicitly shaped additive residual."""

    def predict_batch(S, U):
        return nominal_predict_batch(S, U) + feat_base(S, U) @ W_res

    return predict_batch

The deployment plant and a bounded-command adaptation episodeLink Copied

The deployment plant sits inside the randomized support but is not the nominal plant: $k = 1.25$ , $c = 0.62$ , $b = 0.82$ . We collect 140 steps of smooth excitation from it. That episode stands for a bounded, supervised data-collection run, not for autonomous exploration. The distinction matters: the experiment assumes we were permitted to excite the plant gently for a fixed, short interval under supervision, and it assumes the excitation was chosen to cover the directions we care about rather than to maximize information at any cost.

From that episode we do two things. First, system identification: invert the Euler step to get an acceleration target $\hat{a}_t = (v_{t+1} - v_t)/\Delta t$ , then solve the following linear system. The system identification step recovers the plant parameters $(k, c, b)$ by fitting a linear model to acceleration data. The relationship is:

\hat{a}_t = -k\,p_t - c\,v_t + b\,u_t^{\text{app}}

where:

$\hat{a}_t$ : the acceleration target recovered by finite-differencing the observed velocity
$p_t, v_t$ : the state components at time $t$
$u_t^{\text{app}}$ : the applied (post-clip) command at time $t$
$(k, c, b)$ : the unknown plant parameters we want to estimate

The system is solved via least squares over the adaptation episode to recover the estimates $(\hat{k}, \hat{c}, \hat{b})$ . The inversion works because the semi-implicit Euler step is affine in the parameters: given $(p_t, v_t, v_{t+1}, u_t^{\text{app}})$ , the parameter vector $(k, c, b)$ appears linearly on the right-hand side. A linear least-squares problem in three unknowns is about as cheap as identification gets, which is the point: when the parametric part of the gap is identifiable and you have the excitation to see it, closing it costs almost nothing. Second, residual fitting: compute the residual target $r_t = s_{t+1} - f_{\text{nom}}(s_t, u_t)$ using the nominal model, and fit a shaped estimator to it.

The residual target must use the same nominal transition, timestep and input/target alignment as deployment. A different integrator or a one-step pairing error can introduce a learnable discrepancy that is numerical rather than physical. A residual estimator may fit that artifact with small training error, yet introduce error when applied to correctly aligned deployment inputs. Check the convention before interpreting a fitted correction as evidence of changed plant dynamics.

In[12]:

Code

# Deployment plant (inside randomized support) and a plant outside it.
DEPLOY = {"k": 1.25, "c": 0.62, "b": 0.82}
SHIFT = {"k": 1.85, "c": 1.05, "b": 0.70}

# Bounded adaptation episode from the deployment plant.
adapt_u = exploration_actions(rng, 140)
adapt_states, _ = rollout(DEPLOY, adapt_u, s0=np.array([0.15, -0.05]), rng=rng)
S_a, U_a, S_next_a = build_pairs(
    [{"params": DEPLOY, "s": adapt_states, "u": adapt_u}]
)

# System identification from the same episode.
acc_target = (S_next_a[:, 1] - S_a[:, 1]) / DT
X_sys = np.column_stack([-S_a[:, 0], -S_a[:, 1], np.clip(U_a, -U_MAX, U_MAX)])
theta_hat = np.linalg.lstsq(X_sys, acc_target, rcond=None)[0]

# Residual target, aligned to the nominal model and the same (s, u) pairing.
res_target = S_next_a - nominal_predict_batch(S_a, U_a)
W_res = fit_linear(feat_base(S_a, U_a), res_target)

# Held-out validation episode from the same plant, used only for gating.
val_u = exploration_actions(rng, 120)
val_states, _ = rollout(DEPLOY, val_u, s0=np.array([-0.2, 0.1]), rng=rng)

Out[13]:

Console

Identification from the bounded adaptation episode
  k_hat = 1.250    true k = 1.250
  c_hat = 0.620    true c = 0.620
  b_hat = 0.820    true b = 0.820
  nominal simulator: k = 1.00, c = 0.40, b = 1.00

Adaptation episode length : 140 steps
Validation episode length : 120 steps

Even a simple least-squares identification recovers the deployment parameters to a useful precision from 140 steps of excitation. That is the parametric part of the gap, and it is cheap to close when you have excitation covering the directions you care about. The gap between what a cheap least-squares solve recovers here and what a wide-domain-randomized problem demands sets up the comparison in the next section. Unique noiseless recovery requires the complete design matrix with columns [-p, -v, applied u] to have rank three. Position and velocity must not be collinear, and the command column must add independent information. Varying commands are a useful excitation design, but a constant nonzero command can suffice when the transient states supply the other independent directions. Zero command makes the gain unidentifiable in this regression.

We check both the generating family and the finite sample box against the evaluation plants. The figure places plant parameters beside training draws. The parameter-blind predictor does not receive these parameters, so membership in either box does not classify its state/action predictions as interpolation or extrapolation. The conditioned predictor does receive parameter context, with coverage still depending on joint feature inputs.

In[14]:

Code

# Declared generating family, separate from observed sample extrema.
distribution_bounds = {"k": (0.6, 1.6), "c": (0.2, 0.8), "b": (0.7, 1.3)}
sample_bounds = {
    "k": (float(ks.min()), float(ks.max())),
    "c": (float(cs.min()), float(cs.max())),
    "b": (float(bs.min()), float(bs.max())),
}


def inside_parameter_box(plant, bounds):
    """Coordinate-wise interval membership, not a performance or density test."""
    return all(
        bounds[name][0] <= plant[name] <= bounds[name][1] for name in bounds
    )


deploy_inside_family = inside_parameter_box(DEPLOY, distribution_bounds)
shift_inside_family = inside_parameter_box(SHIFT, distribution_bounds)
deploy_inside_sample_box = inside_parameter_box(DEPLOY, sample_bounds)
shift_inside_sample_box = inside_parameter_box(SHIFT, sample_bounds)
assert deploy_inside_family and deploy_inside_sample_box
assert not shift_inside_family and not shift_inside_sample_box
# A valid family member can lie outside the observed twenty-four-draw box.
assert inside_parameter_box(
    {"k": 0.65, "c": 0.5, "b": 1.0}, distribution_bounds
)
assert not inside_parameter_box({"k": 0.65, "c": 0.5, "b": 1.0}, sample_bounds)

Out[15]:

Visualization

Scatter of training plant parameters with deployment and shifted plants marked. — Twenty-four training draws and two evaluation plants in the (k, c) plane. The outline marks the declared generating-family bounds, not the smaller extrema of the finite sample. The deployment plant is inside that family; the shifted plant is outside its stiffness and damping ranges. This parameter check does not determine whether a parameter-blind predictor interpolates in its state/action inputs. Gain is checked separately in code.

Sample-based planning with a shared actuator envelopeLink Copied

For closed-loop evaluation we use a short-horizon random-shooting model predictive controller. It samples candidate action sequences within the actuator bounds, rolls the predictor forward, and applies the first action of the lowest-cost sequence. The same planner and candidate-action draws are used for each predictor within a regime. Different predicted states can still produce different selected actions and subsequent trajectories; this is a comparison of each predictor coupled to that fixed planner.

At each control tick, the current reference position is held constant across all eight predicted steps. This is not lookahead tracking of the future reference sequence. The stage cost adds squared position error to an effort penalty, $0.02\,u^2$ ; with position in metres and command in metres per second squared, the coefficient has units of seconds to the fourth power so both terms have units of square metres. The chosen weight, horizon and 64 sampled sequences are experimental settings, not optimized controller parameters.

In[16]:

Code

def plan_action(
    predict_batch, s, p_ref, horizon=8, n_samples=64, u_max=U_MAX, rng=None
):
    """Random-shooting MPC over a fixed actuator envelope."""
    cand = rng.uniform(-u_max, u_max, size=(n_samples, horizon))
    S = np.tile(np.asarray(s, dtype=float), (n_samples, 1))
    cost = np.zeros(n_samples)
    for h in range(horizon):
        U = cand[:, h]
        S = predict_batch(S, U)
        cost += (S[:, 0] - p_ref) ** 2 + 0.02 * U**2
    return float(cand[int(np.argmin(cost)), 0])


def close_loop(
    predict_batch, params, p_ref_seq, delay_steps=0, u_max=U_MAX, seed=7
):
    """Closed-loop tracking. The controller sees a delayed observation of the true state."""
    plan_rng = np.random.default_rng(seed)
    n = len(p_ref_seq)
    states = np.zeros((n + 1, 2))
    applied = np.zeros(n)
    for t in range(n):
        obs = states[max(0, t - delay_steps)].copy()
        u = plan_action(
            predict_batch, obs, p_ref_seq[t], u_max=u_max, rng=plan_rng
        )
        states[t + 1], applied[t] = step_plant(
            states[t], u, params, u_max=u_max
        )
    return states, applied

Notice what delay_steps does: the plant evolves from the true states[t], while the controller receives states[max(0, t - delay_steps)]. At zero delay, the controller has exact current-state observations. This is an idealized information condition, not a mathematical lower bound on tracking cost. The finite-horizon sampled planner is not globally optimal, and a change in delay can change its chosen actions in either direction. The zero-delay nominal-model run below is a reference for this specific controller, not an achievable-performance guarantee or a hardware baseline.