Consider a robot with state \( \mathbf{x}_t \in \mathbb{R}^n \) (for example, joint angles and velocities) and control input \( \boldsymbol{\tau}_t \in \mathbb{R}^m \). The real discrete-time dynamics are

\[ \mathbf{x}_{t+1}^{\text{real}} = \mathbf{f}_{\text{real}}(\mathbf{x}_t,\boldsymbol{\tau}_t), \]

but in simulation we only have an approximate model \( \mathbf{f}_{\text{phys}} \), typically obtained from rigid-body dynamics with simplified friction, contact models, or parameter errors:

\[ \mathbf{x}_{t+1}^{\text{sim}} = \mathbf{f}_{\text{phys}}(\mathbf{x}_t,\boldsymbol{\tau}_t). \]

The modeling error is the residual mapping

\[ \mathbf{e}(\mathbf{x}_t,\boldsymbol{\tau}_t) = \mathbf{f}_{\text{real}}(\mathbf{x}_t,\boldsymbol{\tau}_t) - \mathbf{f}_{\text{phys}}(\mathbf{x}_t,\boldsymbol{\tau}_t). \]

Residual physics fits a function approximator \( \mathbf{r}_\theta \) to this error and uses

\[ \mathbf{f}_{\text{hyb}}(\mathbf{x}_t,\boldsymbol{\tau}_t) = \mathbf{f}_{\text{phys}}(\mathbf{x}_t,\boldsymbol{\tau}_t) + \mathbf{r}_\theta(\mathbf{x}_t,\boldsymbol{\tau}_t), \]

whereas pure learning discards the physics and directly learns \( \mathbf{f}_{\theta}^{\text{pure}} \approx \mathbf{f}_{\text{real}} \) from data. The choice has strong implications for sample complexity, extrapolation, and sim-to-real robustness.

flowchart TD
  A["Real robot data (x_t, tau_t, x_{t+1})"] --> B["Analytical model f_phys(x_t, tau_t)"]
  B --> C["Learn residual r_theta(x_t, tau_t)"]
  C --> D["Hybrid model f_phys + r_theta"]
  A --> E["Ignore physics"]
  E --> F["Learn f_theta_pure(x_t, tau_t)"]
  D --> G["Sim-to-real deployment"]
  F --> G
        

Intuitively, residual physics says: “Most of the structure is known; learn only what we do not know.” Pure learning says: “Assume nothing; let the network discover everything.” In the rest of this lesson, we formalize this intuition and derive optimization problems and properties for both approaches.

Recall the continuous-time manipulator dynamics for generalized coordinates \( \mathbf{q} \) and velocities \( \dot{\mathbf{q}} \):

\[ \mathbf{M}(\mathbf{q})\ddot{\mathbf{q}} + \mathbf{C}(\mathbf{q},\dot{\mathbf{q}})\dot{\mathbf{q}} + \mathbf{g}(\mathbf{q}) + \mathbf{f}_d(\dot{\mathbf{q}}) = \boldsymbol{\tau} + \mathbf{d}(\mathbf{q},\dot{\mathbf{q}},\boldsymbol{\tau}). \]

Here:

• \( \mathbf{M}(\mathbf{q}) \) is positive definite inertia,
• \( \mathbf{C}(\mathbf{q},\dot{\mathbf{q}})\dot{\mathbf{q}} \) are Coriolis/centrifugal terms,
• \( \mathbf{g}(\mathbf{q}) \) is gravity,
• \( \mathbf{f}_d(\dot{\mathbf{q}}) \) is a simple friction model,
• \( \mathbf{d}(\mathbf{q},\dot{\mathbf{q}},\boldsymbol{\tau}) \) collects unmodeled effects (flexibilities, complex friction, contact, etc.).

In simulation we typically ignore the last term and use a nominal model

\[ \mathbf{M}_0(\mathbf{q})\ddot{\mathbf{q}} + \mathbf{C}_0(\mathbf{q},\dot{\mathbf{q}})\dot{\mathbf{q}} + \mathbf{g}_0(\mathbf{q}) + \mathbf{f}_{d,0}(\dot{\mathbf{q}}) = \boldsymbol{\tau}. \]

After time discretization (e.g. forward Euler with step \( \Delta t \)), we can write the nominal discrete-time dynamics as

\[ \mathbf{x}_{t+1}^{\text{sim}} = \mathbf{f}_{\text{phys}}(\mathbf{x}_t,\boldsymbol{\tau}_t), \quad \mathbf{x}_t = \begin{bmatrix}\mathbf{q}_t \\ \dot{\mathbf{q}}_t\end{bmatrix}. \]

The true discrete-time dynamics are instead

\[ \mathbf{x}_{t+1}^{\text{real}} = \mathbf{f}_{\text{real}}(\mathbf{x}_t,\boldsymbol{\tau}_t) = \mathbf{f}_{\text{phys}}(\mathbf{x}_t,\boldsymbol{\tau}_t) + \mathbf{e}(\mathbf{x}_t,\boldsymbol{\tau}_t), \]

with error term \( \mathbf{e}(\mathbf{x}_t,\boldsymbol{\tau}_t) \) induced by parameter mismatch and neglected physics. Sim-to-real methods must reduce the impact of \( \mathbf{e} \) on the closed-loop behavior.

A residual physics model augments the nominal simulator with a learned correction:

\[ \mathbf{f}_{\theta}^{\text{res}}(\mathbf{x}_t,\boldsymbol{\tau}_t) = \mathbf{f}_{\text{phys}}(\mathbf{x}_t,\boldsymbol{\tau}_t) + \mathbf{r}_\theta(\mathbf{x}_t,\boldsymbol{\tau}_t), \]

where \( \mathbf{r}_\theta \) is usually a neural network with parameters \( \theta \). Given real transition data

\[ \mathcal{D} = \big\{(\mathbf{x}_t,\boldsymbol{\tau}_t, \mathbf{x}_{t+1}^{\text{real}})\big\}_{t=1}^N, \]

we define the empirical loss

\[ L_{\text{res}}(\theta) = \frac{1}{N}\sum_{t=1}^N \left\| \mathbf{x}_{t+1}^{\text{real}} - \mathbf{f}_{\text{phys}}(\mathbf{x}_t,\boldsymbol{\tau}_t) - \mathbf{r}_\theta(\mathbf{x}_t,\boldsymbol{\tau}_t) \right\|^2. \]

In the limit of infinite data and model capacity, with squared loss and i.i.d. sampling from the real transition distribution, we can derive the optimal residual.

Proposition (Optimal residual as conditional mean error).

Assume that:

• the training samples are drawn i.i.d. from a distribution over \( (\mathbf{x}_t,\boldsymbol{\tau}_t) \),
• \( \mathbf{r}_\theta \) can represent any square-integrable function,
• we minimize the population loss

\[ \mathcal{L}_{\text{res}}(\theta) = \mathbb{E}\left[ \left\| \mathbf{x}_{t+1}^{\text{real}} - \mathbf{f}_{\text{phys}}(\mathbf{x}_t,\boldsymbol{\tau}_t) - \mathbf{r}_\theta(\mathbf{x}_t,\boldsymbol{\tau}_t) \right\|^2 \right]. \]

Then any minimizer \( \theta^\star \) satisfies

\[ \mathbf{r}_{\theta^\star}(\mathbf{x},\boldsymbol{\tau}) = \mathbb{E}\!\left[ \mathbf{e}(\mathbf{x},\boldsymbol{\tau}) \,\middle|\, \mathbf{x},\boldsymbol{\tau} \right], \]

i.e., the residual network learns the conditional mean of the modeling error.

Sketch of proof.

Fix \( (\mathbf{x},\boldsymbol{\tau}) \). Define random variable \( \mathbf{Z} = \mathbf{x}_{t+1}^{\text{real}} - \mathbf{f}_{\text{phys}}(\mathbf{x},\boldsymbol{\tau}) \). Then the integrand of \( \mathcal{L}_{\text{res}}(\theta) \) is

\[ \mathbb{E} \big[ \|\mathbf{Z} - \mathbf{r}_\theta(\mathbf{x},\boldsymbol{\tau})\|^2 \,\big|\, \mathbf{x},\boldsymbol{\tau} \big]. \]

The unique minimizer of \( \mathbb{E}[\|\mathbf{Z}-\mathbf{u}\|^2] \) over \( \mathbf{u} \) is the conditional mean \( \mathbb{E}[\mathbf{Z}] \). Substituting back \( \mathbf{Z} = \mathbf{f}_{\text{real}} - \mathbf{f}_{\text{phys}} = \mathbf{e} \) yields the claim.

Thus, residual physics acts as a statistical estimator of the unmodeled dynamics on top of a structured analytical prior.

Pure learning discards the simulator and directly parameterizes \( \mathbf{f}_{\theta}^{\text{pure}} \), for example with a deep neural network:

\[ \mathbf{f}_{\theta}^{\text{pure}} : (\mathbf{x}_t,\boldsymbol{\tau}_t) \mapsto \hat{\mathbf{x}}_{t+1}. \]

The typical loss is

\[ L_{\text{pure}}(\theta) = \frac{1}{N}\sum_{t=1}^N \left\| \mathbf{x}_{t+1}^{\text{real}} - \mathbf{f}_{\theta}^{\text{pure}}(\mathbf{x}_t,\boldsymbol{\tau}_t) \right\|^2. \]

Proposition (Optimal pure model as conditional mean).

Under the same assumptions as in the previous section, any minimizer \( \theta^\star \) of the population version of \( L_{\text{pure}} \) satisfies

\[ \mathbf{f}_{\theta^\star}^{\text{pure}}(\mathbf{x},\boldsymbol{\tau}) = \mathbb{E}\!\left[ \mathbf{x}_{t+1}^{\text{real}} \,\middle|\, \mathbf{x},\boldsymbol{\tau} \right]. \]

Bias and variance comparison.

Decompose the real dynamics as

\[ \mathbf{f}_{\text{real}} = \mathbf{f}_{\text{phys}} + \mathbf{e}. \]

If the physical model already explains most of the structure, then \( \mathbf{e} \) has smaller magnitude and lower complexity than \( \mathbf{f}_{\text{real}} \) itself. A residual model only has to learn \( \mathbf{e} \), which can significantly reduce:

• Approximation error (simpler function to approximate),
• Estimation error (fewer effective degrees of freedom to fit),
• Extrapolation error (the physical prior handles out-of-distribution regions).

In contrast, pure learning must discover conservation laws, constraints, and contact behavior from scratch, which typically requires much more data and careful regularization. However, if the analytical model is severely misspecified (e.g. wrong topology, wrong contact structure), residual physics may inherit its biases, while pure learning has the potential to escape them given enough data.

For residual physics in sim-to-real, a common objective combines next-state prediction accuracy with regularization that encourages the residual to remain small:

\[ L_{\text{total}}(\theta) = L_{\text{res}}(\theta) + \lambda \, \mathbb{E}\big[ \|\mathbf{r}_\theta(\mathbf{x}_t,\boldsymbol{\tau}_t)\|^2 \big], \quad \lambda > 0. \]

This penalizes overly aggressive corrections, steering the optimizer to trust the analytical model wherever consistent with the data.

Multi-step error accumulation.

Let \( \mathbf{x}_{t+k}^{\text{real}} \) and \( \mathbf{x}_{t+k}^{\theta} \) be the states obtained by rolling out the real dynamics and the learned model under the same control sequence. Assume that both maps are Lipschitz with constant \( L \) in \( \mathbf{x} \) and that the one-step prediction error of the learned model is uniformly bounded:

\[ \left\| \mathbf{f}_{\theta}(\mathbf{x},\boldsymbol{\tau}) - \mathbf{f}_{\text{real}}(\mathbf{x},\boldsymbol{\tau}) \right\| \leq \delta \quad \text{for all } (\mathbf{x},\boldsymbol{\tau}). \]

Then by induction one can show the bound

\[ \big\| \mathbf{x}_{t+k}^{\text{real}} - \mathbf{x}_{t+k}^{\theta} \big\| \leq \delta \, \frac{L^k - 1}{L - 1}, \]

so reducing \( \delta \) even moderately can have a large effect on long-horizon prediction error. Residual physics aims to reduce \( \delta \) by leveraging analytical structure.

flowchart TD
  D1["Collect real transitions (x_t, tau_t, x_{t+1})"] --> S1["Evaluate f_phys(x_t, tau_t)"]
  S1 --> R1["Compute labels: x_{t+1} - f_phys(x_t, tau_t)"]
  R1 --> N1["Train residual network r_theta"]
  N1 --> H1["Hybrid model f_phys + r_theta"]
  H1 --> P1["Use in planning / control in simulation"]
  P1 --> R2["Deploy policy on real robot"]
        

In sim-to-real, the residual can be trained directly on real-world data (on-policy or off-policy) while keeping the bulk of planning and policy learning in simulation using \( \mathbf{f}_{\text{phys}} \).

We illustrate a simple Python implementation of residual dynamics. We assume access to a rigid-body dynamics library (e.g., Pinocchio, RBDL-like bindings, or a simulator) providing simulate_phys(x, tau, dt). We build a small neural network to predict the residual and train it using PyTorch-style code.

import numpy as np
import torch
import torch.nn as nn
import torch.optim as optim

# Example dimensions for a 7-DOF manipulator (q, qdot)
n_q = 7
n_x = 2 * n_q
n_u = n_q
dt = 0.01

class ResidualNet(nn.Module):
    def __init__(self, state_dim, act_dim, hidden_dim=128):
        super().__init__()
        self.net = nn.Sequential(
            nn.Linear(state_dim + act_dim, hidden_dim),
            nn.Tanh(),
            nn.Linear(hidden_dim, hidden_dim),
            nn.Tanh(),
            nn.Linear(hidden_dim, state_dim),
        )

    def forward(self, x, u):
        # x: (batch, n_x), u: (batch, n_u)
        inp = torch.cat([x, u], dim=-1)
        return self.net(inp)

def simulate_phys(x_np, u_np, dt):
    """
    Placeholder for a call to a robotics physics engine.
    x_np: np.array with shape (n_x,)
    u_np: np.array with shape (n_u,)
    Returns next state x_next_np.
    """
    # Here we just do identity plus small drift as a stand-in.
    # In practice, call your simulator or analytic dynamics integrator.
    return x_np + dt * np.concatenate([x_np[n_q:], np.zeros_like(x_np[n_q:])])

# Build dataset from real rollouts (collected on the real robot).
# real_data is assumed to be an array of triples (x_t, u_t, x_tp1_real).
real_data = []  # fill from experiments
# Example: real_data.append((x_t_np, u_t_np, x_tp1_real_np))

model = ResidualNet(n_x, n_u)
optimizer = optim.Adam(model.parameters(), lr=1e-3)
lambda_reg = 1e-3

def train_epoch(real_data, batch_size=64):
    np.random.shuffle(real_data)
    total_loss = 0.0
    for i in range(0, len(real_data), batch_size):
        batch = real_data[i:i + batch_size]
        x_list, u_list, x_next_real_list = zip(*batch)
        x = torch.tensor(np.stack(x_list), dtype=torch.float32)
        u = torch.tensor(np.stack(u_list), dtype=torch.float32)
        x_next_real = torch.tensor(np.stack(x_next_real_list), dtype=torch.float32)

        # Nominal physics prediction
        x_next_phys_list = []
        for j in range(x.shape[0]):
            x_np = x[j].detach().numpy()
            u_np = u[j].detach().numpy()
            x_next_phys_np = simulate_phys(x_np, u_np, dt)
            x_next_phys_list.append(x_next_phys_np)
        x_next_phys = torch.tensor(np.stack(x_next_phys_list), dtype=torch.float32)

        # Residual prediction
        r = model(x, u)
        x_next_hyb = x_next_phys + r

        loss_pred = ((x_next_hyb - x_next_real) ** 2).mean()
        loss_reg = (r ** 2).mean()
        loss = loss_pred + lambda_reg * loss_reg

        optimizer.zero_grad()
        loss.backward()
        optimizer.step()
        total_loss += loss.item() * x.shape[0]

    return total_loss / len(real_data)

# After training, use f_hyb(x, u) = simulate_phys(x, u, dt) + model(x, u)
      

In practice, simulate_phys calls a differentiable or black-box simulator derived from your analytical dynamics. The residual network corrects for the mismatch observed on the real robot.

In C++, we often use libraries like Eigen and Pinocchio or RBDL for dynamics. Below is a skeleton for integrating a residual model (here, assumed to be a small multilayer perceptron) with an existing physics step. The network weights could be loaded from disk after offline training.

#include <Eigen/Dense>

// Example dimensions
constexpr int NQ = 7;
constexpr int NX = 2 * NQ;
constexpr int NU = NQ;

// Toy residual MLP with one hidden layer
struct ResidualMLP {
  Eigen::Matrix<double, NX + NU, NX> W1;
  Eigen::Matrix<double, NX, 1> b1;
  Eigen::Matrix<double, NX, NX> W2;
  Eigen::Matrix<double, NX, 1> b2;

  Eigen::Matrix<double, NX, 1> forward(
      const Eigen::Matrix<double, NX, 1> &x,
      const Eigen::Matrix<double, NU, 1> &u) const {
    Eigen::Matrix<double, NX + NU, 1> in;
    in << x, u;
    Eigen::Matrix<double, NX, 1> h = (W1.transpose() * in + b1).array().tanh();
    Eigen::Matrix<double, NX, 1> out = W2.transpose() * h + b2;
    return out;
  }
};

// Placeholder for analytical discrete-time dynamics
Eigen::Matrix<double, NX, 1> f_phys(
    const Eigen::Matrix<double, NX, 1> &x,
    const Eigen::Matrix<double, NU, 1> &tau,
    double dt) {
  Eigen::Matrix<double, NX, 1> x_next;
  // Here we just do a dummy integrator; in practice call your rigid-body dynamics
  x_next.head(NQ) = x.head(NQ) + dt * x.tail(NQ);
  x_next.tail(NQ).setZero();
  return x_next;
}

Eigen::Matrix<double, NX, 1> f_hybrid(
    const Eigen::Matrix<double, NX, 1> &x,
    const Eigen::Matrix<double, NU, 1> &tau,
    double dt,
    const ResidualMLP &residual) {
  Eigen::Matrix<double, NX, 1> x_phys = f_phys(x, tau, dt);
  Eigen::Matrix<double, NX, 1> r = residual.forward(x, tau);
  return x_phys + r;
}

int main() {
  ResidualMLP res;
  // TODO: load weights into res.W1, res.b1, res.W2, res.b2

  Eigen::Matrix<double, NX, 1> x;
  Eigen::Matrix<double, NU, 1> tau;
  x.setZero();
  tau.setZero();

  double dt = 0.01;
  Eigen::Matrix<double, NX, 1> x_next = f_hybrid(x, tau, dt, res);

  // Use x_next in your controller or simulator
  return 0;
}
      

Training of the residual network parameters can be performed in Python or C++ (e.g., with a deep learning library), then exported and loaded as matrices and vectors in C++ for runtime evaluation.

In Java, robotics libraries are less standardized, but one can implement the residual around any available dynamics code or even directly around a simulator. Below is a small example using arrays only (no external libraries) for clarity.

public class ResidualModel {
    private final int nX;
    private final int nU;
    private final int hidden;

    // Simple fully-connected residual net with tanh activation
    private double[][] W1; // (hidden, nX + nU)
    private double[] b1;   // (hidden)
    private double[][] W2; // (nX, hidden)
    private double[] b2;   // (nX)

    public ResidualModel(int nX, int nU, int hidden) {
        this.nX = nX;
        this.nU = nU;
        this.hidden = hidden;
        W1 = new double[hidden][nX + nU];
        b1 = new double[hidden];
        W2 = new double[nX][hidden];
        b2 = new double[nX];
        // TODO: initialize or load from file
    }

    public double[] forward(double[] x, double[] u) {
        double[] in = new double[nX + nU];
        System.arraycopy(x, 0, in, 0, nX);
        System.arraycopy(u, 0, in, nX, nU);

        double[] h = new double[hidden];
        for (int i = 0; i != hidden; ++i) {
            double sum = b1[i];
            for (int j = 0; j != nX + nU; ++j) {
                sum += W1[i][j] * in[j];
            }
            h[i] = Math.tanh(sum);
        }

        double[] out = new double[nX];
        for (int i = 0; i != nX; ++i) {
            double sum = b2[i];
            for (int j = 0; j != hidden; ++j) {
                sum += W2[i][j] * h[j];
            }
            out[i] = sum;
        }
        return out;
    }

    // Placeholder analytic dynamics; replace with your model
    public double[] fPhys(double[] x, double[] tau, double dt) {
        double[] xNext = new double[nX];
        int nQ = nX / 2;
        for (int i = 0; i != nQ; ++i) {
            xNext[i] = x[i] + dt * x[i + nQ]; // q_{t+1}
        }
        for (int i = 0; i != nQ; ++i) {
            xNext[i + nQ] = 0.0; // very simple
        }
        return xNext;
    }

    public double[] fHybrid(double[] x, double[] tau, double dt) {
        double[] xPhys = fPhys(x, tau, dt);
        double[] r = forward(x, tau);
        double[] xNext = new double[nX];
        for (int i = 0; i != nX; ++i) {
            xNext[i] = xPhys[i] + r[i];
        }
        return xNext;
    }
}
      

The Java side is primarily for runtime evaluation. Training can be done offline and weights serialized into the Java arrays W1, W2, b1, and b2.

9.1 MATLAB Script with Simulink Block for Residual Physics

In MATLAB, we typically integrate residual physics with Simulink. The analytical model is implemented as a Simulink subsystem, and the residual network is implemented as a custom block (e.g., via MATLAB Function block) that takes state and torque and outputs a correction.

% Dimensions and time step
nq = 7;
nx = 2 * nq;
nu = nq;
dt = 0.01;

% Example: one-step hybrid prediction in MATLAB
function x_next = hybrid_step(x, tau, dt, residual_params)
    % x: nx-by-1, tau: nu-by-1
    x_phys = f_phys(x, tau, dt);       % call to analytic discrete-time model
    r = residual_net(x, tau, residual_params);
    x_next = x_phys + r;
end

function x_next = f_phys(x, tau, dt)
    % Very simple placeholder; replace with your rigid-body dynamics.
    nq = numel(x) / 2;
    q = x(1:nq);
    qdot = x(nq+1:end);
    q_next = q + dt * qdot;
    qdot_next = zeros(size(qdot));
    x_next = [q_next; qdot_next];
end

function r = residual_net(x, tau, params)
    % params.W1: [H, nx+nu], params.b1: [H, 1]
    % params.W2: [nx, H], params.b2: [nx, 1]
    in = [x; tau];
    h = tanh(params.W1 * in + params.b1);
    r = params.W2 * h + params.b2;
end

% In Simulink:
%  - Implement f_phys as a subsystem.
%  - Implement residual_net as a MATLAB Function block using the above code.
%  - Sum the outputs to produce x_next.
      

9.2 Wolfram Mathematica: Residual Correction of Dynamics

Wolfram Mathematica can directly integrate residual-corrected dynamics using NDSolve. We treat the residual as a function rTheta[x, u] obtained, for example, by fitting a neural network externally and importing it as an interpolating function.

(* State: x = {q1, ..., qn, q1dot, ..., qndot} *)
nq = 2;
nx = 2 * nq;

(* Placeholder analytic continuous-time dynamics fPhys[x_, u_] *)
fPhys[x_, u_] := Module[{q, qdot, dqdt, dqdotdt},
  q = x[[1 ;; nq]];
  qdot = x[[nq + 1 ;; nx]];
  dqdt = qdot;
  dqdotdt = ConstantArray[0.0, nq];
  Join[dqdt, dqdotdt]
];

(* Residual dynamics, assumed given as an interpolating function *)
rTheta[x_, u_] := ConstantArray[0.0, nx]; (* replace with learned residual *)

(* Hybrid continuous-time dynamics *)
fHybrid[x_, u_] := fPhys[x, u] + rTheta[x, u];

(* Example control input as a function of time *)
uFunc[t_] := {Sin[t], Cos[t]};

(* ODE system for NDSolve *)
odeSystem = Table[
   x[i]'[t] == fHybrid[Table[x[j][t], {j, 1, nx}], uFunc[t]][[i]],
   {i, 1, nx}
];

(* Initial conditions *)
x0 = ConstantArray[0.0, nx];
ic = Table[x[i][0] == x0[[i]], {i, 1, nx}];

sol = NDSolve[
   Join[odeSystem, ic],
   Table[x[i], {i, 1, nx}],
   {t, 0, 5}
];

(* Evaluate trajectory at t = 5 *)
xFinal = Table[x[i][5] /. sol[[1]], {i, 1, nx}];
      

In practice, rTheta could be obtained by fitting a surrogate model to data, exporting it as a Mathematica function (e.g., via Import from ONNX or manually as an interpolating grid).

Problem 1 (Residual vs Pure: Optimal Functions).
Let \( \mathbf{f}_{\text{real}} \) and \( \mathbf{f}_{\text{phys}} \) be as defined above. Show that if both residual physics and pure learning use squared loss with infinite capacity and data, then the optimal hybrid model \( \mathbf{f}_{\theta^\star}^{\text{res}} \) and optimal pure model \( \mathbf{f}_{\theta^\star}^{\text{pure}} \) satisfy

\[ \mathbf{f}_{\theta^\star}^{\text{res}}(\mathbf{x},\boldsymbol{\tau}) = \mathbf{f}_{\text{phys}}(\mathbf{x},\boldsymbol{\tau}) + \mathbb{E}[\mathbf{e}(\mathbf{x},\boldsymbol{\tau}) \mid \mathbf{x},\boldsymbol{\tau}], \]

\[ \mathbf{f}_{\theta^\star}^{\text{pure}}(\mathbf{x},\boldsymbol{\tau}) = \mathbb{E}[\mathbf{f}_{\text{real}}(\mathbf{x},\boldsymbol{\tau}) \mid \mathbf{x},\boldsymbol{\tau}], \]

and conclude that the two differ by the deterministic component of the physical model.

Solution.
From Sections 3 and 4, we already know that the optimal residual and pure models correspond to conditional means of \( \mathbf{e} \) and \( \mathbf{x}_{t+1}^{\text{real}} \), respectively. Using \( \mathbf{f}_{\text{real}} = \mathbf{f}_{\text{phys}} + \mathbf{e} \), we have

\[ \mathbf{f}_{\theta^\star}^{\text{res}} = \mathbf{f}_{\text{phys}} + \mathbb{E}[\mathbf{e} \mid \mathbf{x},\boldsymbol{\tau}], \]

\[ \mathbf{f}_{\theta^\star}^{\text{pure}} = \mathbb{E}[\mathbf{f}_{\text{phys}} + \mathbf{e} \mid \mathbf{x},\boldsymbol{\tau}] = \mathbf{f}_{\text{phys}} + \mathbb{E}[\mathbf{e} \mid \mathbf{x},\boldsymbol{\tau}], \]

whenever \( \mathbf{f}_{\text{phys}} \) is deterministic given \( (\mathbf{x},\boldsymbol{\tau}) \). Therefore in the infinite-data, infinite-capacity limit, both residual and pure learning converge to the same conditional mean model; the difference lies in finite-sample behavior and implicit regularization.

Problem 2 (Sample Complexity Intuition).
Suppose we have a parametric residual model class \( \{\mathbf{r}_\theta\} \) with \( d_{\text{res}} \) free parameters and a pure model class \( \{\mathbf{f}_\theta^{\text{pure}}\} \) with \( d_{\text{pure}} \gg d_{\text{res}} \). Under a classical learning-theoretic bound, the generalization error scales as \( \mathcal{O}\big(\sqrt{d/N}\big) \). Explain why, for the same number of real samples \( N \), the residual model is likely to generalize better than the pure learning model when the physical model is reasonably accurate.

Solution.
The generalization error bound of order \( \sqrt{d/N} \) suggests that the estimation error grows with \( \sqrt{d} \). Residual modeling keeps \( d_{\text{res}} \) small because much of the structure is hard-coded in \( \mathbf{f}_{\text{phys}} \). Pure learning requires a larger function class, with \( d_{\text{pure}} \gg d_{\text{res}} \), to express the same behaviors. Hence, for the same \( N \), the residual model typically has a smaller estimation error, provided that the physical prior is not badly misspecified. In addition, the data do not have to explain features already captured by \( \mathbf{f}_{\text{phys}} \), further reducing effective complexity.

Problem 3 (Long-Horizon Error Bound).
Assume the conditions of Section 5: \( \mathbf{f}_{\text{real}} \) and \( \mathbf{f}_\theta \) are Lipschitz in \( \mathbf{x} \) with constant \( L \) and the one-step prediction error is bounded by \( \delta \). Prove by induction on \( k \) that

\[ \big\| \mathbf{x}_{t+k}^{\text{real}} - \mathbf{x}_{t+k}^{\theta} \big\| \leq \delta \, \frac{L^k - 1}{L - 1}. \]

Solution.
For \( k = 1 \), the bound reduces to \( \|\mathbf{x}_{t+1}^{\text{real}} - \mathbf{x}_{t+1}^\theta\| \leq \delta \), which holds by assumption. Suppose it holds for \( k \). For \( k+1 \) we have

\[ \begin{aligned} \big\| \mathbf{x}_{t+k+1}^{\text{real}} - \mathbf{x}_{t+k+1}^{\theta} \big\| &= \big\| \mathbf{f}_{\text{real}}(\mathbf{x}_{t+k}^{\text{real}},\boldsymbol{\tau}_{t+k}) - \mathbf{f}_{\theta}(\mathbf{x}_{t+k}^{\theta},\boldsymbol{\tau}_{t+k}) \big\| \\ &\leq \big\| \mathbf{f}_{\text{real}}(\mathbf{x}_{t+k}^{\text{real}},\boldsymbol{\tau}_{t+k}) - \mathbf{f}_{\text{real}}(\mathbf{x}_{t+k}^{\theta},\boldsymbol{\tau}_{t+k}) \big\| \\ &\quad + \big\| \mathbf{f}_{\text{real}}(\mathbf{x}_{t+k}^{\theta},\boldsymbol{\tau}_{t+k}) - \mathbf{f}_{\theta}(\mathbf{x}_{t+k}^{\theta},\boldsymbol{\tau}_{t+k}) \big\| \\ &\leq L \, \big\| \mathbf{x}_{t+k}^{\text{real}} - \mathbf{x}_{t+k}^{\theta} \big\| + \delta. \end{aligned} \]

Applying the induction hypothesis gives

\[ \big\| \mathbf{x}_{t+k+1}^{\text{real}} - \mathbf{x}_{t+k+1}^{\theta} \big\| \leq L \, \delta \frac{L^k - 1}{L - 1} + \delta = \delta \frac{L^{k+1} - 1}{L - 1}, \]

completing the induction.

Problem 4 (Residual Magnitude Regularization).
Consider the regularized residual objective

\[ L_{\text{total}}(\theta) = L_{\text{res}}(\theta) + \lambda \, \mathbb{E}\big[ \|\mathbf{r}_\theta(\mathbf{x},\boldsymbol{\tau})\|^2 \big]. \]

Derive the first-order optimality condition for \( \mathbf{r}_\theta \) in functional form (ignoring parameterization), and interpret the effect of \( \lambda > 0 \).

Solution.
Treat \( \mathbf{r}_\theta \) as a free function \( \mathbf{r} \). The population objective is

\[ \mathcal{L}_{\text{total}}[\mathbf{r}] = \mathbb{E}\big[ \|\mathbf{x}_{t+1}^{\text{real}} - \mathbf{f}_{\text{phys}}(\mathbf{x},\boldsymbol{\tau}) - \mathbf{r}(\mathbf{x},\boldsymbol{\tau})\|^2 \big] + \lambda \mathbb{E}\big[ \|\mathbf{r}(\mathbf{x},\boldsymbol{\tau})\|^2 \big]. \]

For fixed \( (\mathbf{x},\boldsymbol{\tau}) \), differentiating with respect to \( \mathbf{r} \) and setting the derivative to zero yields

\[ -2 \, \mathbb{E}\big[ \mathbf{x}_{t+1}^{\text{real}} - \mathbf{f}_{\text{phys}} - \mathbf{r}^\star \mid \mathbf{x},\boldsymbol{\tau} \big] + 2 \lambda \, \mathbf{r}^\star(\mathbf{x},\boldsymbol{\tau}) = \mathbf{0}. \]

Solving for \( \mathbf{r}^\star \) gives

\[ \mathbf{r}^\star(\mathbf{x},\boldsymbol{\tau}) = \frac{1}{1 + \lambda} \, \mathbb{E}\big[ \mathbf{e}(\mathbf{x},\boldsymbol{\tau}) \mid \mathbf{x},\boldsymbol{\tau} \big], \]

so regularization shrinks the residual toward zero by a factor \( 1/(1+\lambda) \). This formalizes the intuition that nonzero \( \lambda \) encourages trust in the physical model unless the data strongly indicate otherwise.

In this lesson we formalized the distinction between residual physics and pure learning for robot dynamics in the sim-to-real context. Starting from standard rigid-body dynamics, we derived the residual formulation and showed that, under squared loss, residual models estimate the conditional mean of the modeling error on top of an analytical prior. We compared this to pure learning, which must learn the entire dynamics map, and discussed implications for sample complexity and long-horizon error. We then implemented residual dynamics in Python, C++, Java, MATLAB/Simulink, and Wolfram Mathematica, emphasizing how these abstractions integrate with existing physics engines. Finally, we developed problem sets that connect training objectives, regularization, and multi-step error propagation to the theoretical foundation of residual physics vs pure learning.