Let the robot state be \( x(t) \in \mathbb{R}^n \) and control input be \( u(t) \in \mathbb{R}^m \). The real system evolves as \( \dot{x} = f_{\text{real}}(x,u,t) \), while the simulator uses \( \dot{x} = f_{\text{sim}}(x,u,t;\theta) \) with parameters \( \theta \) (masses, frictions, delays, contact constants, sensor models).

The sim-to-real gap is the trajectory discrepancy \( e(t) = x_{\text{real}}(t) - x_{\text{sim}}(t) \). We want to (i) quantify the gap, (ii) bound it analytically, and (iii) validate whether a simulator or digital twin is “good enough” for a task.

flowchart TD
  A["Specify task + metrics"] --> B["Build simulator model f_sim(x,u;theta)"]
  B --> C["Collect real data under same u(t)"]
  C --> D["Compute mismatch e(t)=x_real-x_sim"]
  D --> E["Estimate/adjust parameters theta_hat"]
  E --> F["Re-simulate with theta_hat"]
  F --> G["Validation tests pass?"]
  G -->|no| E
  G -->|yes| H["Use sim/twin for design + deployment"]
        

The key idea: treat simulation as a mathematical model and apply the same rigor we use in linear control for model uncertainty and validation.

Write the difference between real and simulated vector fields as \( \Delta f(x,u,t) = f_{\text{real}}(x,u,t) - f_{\text{sim}}(x,u,t;\theta) \). We decompose

\[ \Delta f = \Delta f_{\text{param}} + \Delta f_{\text{struct}} + \Delta f_{\text{num}} + \Delta f_{\text{sens}} . \]

• \( \Delta f_{\text{param}} \): wrong physical parameters (mass, inertia, friction).
• \( \Delta f_{\text{struct}} \): missing physics (e.g., unmodeled compliance, stiction).
• \( \Delta f_{\text{num}} \): numerical artifacts (time step, solver tolerances).
• \( \Delta f_{\text{sens}} \): sensor/actuator model mismatch (noise, delay, saturation).

In an introductory setting, we emphasize parameter and numerical errors, because these can be corrected using data and careful simulation settings without invoking advanced robotics.

flowchart LR
  R["Real robot"] -->|unknown| F1["f_real(x,u)"]
  S["Simulator"] -->|uses| F2["f_sim(x,u;theta)"]
  F1 --> DF["Delta f = f_real - f_sim"]
  F2 --> DF
  DF --> P["Param error"]
  DF --> ST["Structural error"]
  DF --> N["Numerical error"]
  DF --> SE["Sensor/actuator error"]
        

Assume both systems start from the same state \( x(0)=x_0 \), and inputs are identical. Let \( f_{\text{sim}} \) be Lipschitz in \( x \) with constant \( L \), i.e.,

\[ \| f_{\text{sim}}(x_1,u,t;\theta) - f_{\text{sim}}(x_2,u,t;\theta) \| \le L \|x_1-x_2\|. \]

Then the gap dynamics satisfy:

\[ \dot{e}(t)= f_{\text{real}}(x_{\text{real}},u,t)-f_{\text{sim}}(x_{\text{sim}},u,t;\theta) = \Delta f(x_{\text{real}},u,t) + \\ \big(f_{\text{sim}}(x_{\text{real}},u,t;\theta)-f_{\text{sim}}(x_{\text{sim}},u,t;\theta)\big) \]

Proposition (Gap bound). If \( \|\Delta f(x,u,t)\| \le \delta \) for all relevant \( (x,u,t) \), then

\[ \|e(t)\| \le \frac{\delta}{L}\left(e^{Lt}-1\right). \]

Proof. Take norms and apply the Lipschitz inequality:

\[ \|\dot{e}(t)\| \le \|\Delta f(x_{\text{real}},u,t)\| + L\|e(t)\| \le \delta + L\|e(t)\|. \]

Define \( g(t)=\|e(t)\| \). Then \( \dot{g}(t) \le \delta + L g(t) \) with \( g(0)=0 \). Multiply by \( e^{-Lt} \):

\[ \frac{d}{dt}\big(e^{-Lt}g(t)\big) \le \delta e^{-Lt}. \]

Integrate from \( 0 \) to \( t \):

\[ e^{-Lt}g(t) - g(0) \le \delta \int_0^t e^{-Ls}\,ds = \delta \frac{1-e^{-Lt}}{L}. \]

Since \( g(0)=0 \),

\[ g(t) \le \frac{\delta}{L}(e^{Lt}-1). \]

This shows that even small per-step modeling error \( \delta \) can grow exponentially in time when dynamics are sensitive (large \( L \)). Thus validation must focus on the task horizon and sensitivity.

For many robots near an operating point, a simulator and real plant can be approximated as LTI systems (students already know state-space).

Real: \( \dot{x}=A x + B u \), simulated: \( \dot{x}=(A+\Delta A)x + (B+\Delta B)u \).

The gap state satisfies

\[ \dot{e} = A e + \Delta A\,x_{\text{sim}} + \Delta B\,u. \]

If \( A \) is stable (all eigenvalues with negative real part), then for bounded signals \( x_{\text{sim}},u \) the error remains bounded. Using the induced norm and the resolvent of \( A \),

\[ \|e(t)\| \le \int_0^t \|e^{A(t-s)}\| \big(\|\Delta A\|\,\|x_{\text{sim}}(s)\| + \|\Delta B\|\,\|u(s)\|\big)\,ds. \]

Corollary. If \( \|e^{At}\| \le M e^{-\alpha t} \) for some \( M\ge 1,\alpha>0 \), and \( \|x_{\text{sim}}(t)\|\le X_{\max} \), \( \|u(t)\|\le U_{\max} \), then

\[ \sup_{t\ge 0}\|e(t)\| \le \frac{M}{\alpha} \left(\|\Delta A\|X_{\max} + \|\Delta B\|U_{\max}\right). \]

This aligns with robustness ideas from linear control: stability margins and gain bounds translate to sim-to-real reliability.

Suppose the simulator depends on a parameter vector \( \theta \in \mathbb{R}^p \) (e.g., friction coefficient, actuator constant). We collect real data \( \{x_k,u_k\}_{k=0}^{N} \) sampled at \( T_s \).

A simple discrete-time model (first-order Euler) is \( x_{k+1} \approx x_k + T_s f_{\text{sim}}(x_k,u_k;\theta) \). Define residuals \( r_k(\theta)= x_{k+1}-x_k - T_s f_{\text{sim}}(x_k,u_k;\theta) \).

We estimate parameters by least squares:

\[ \hat{\theta} = \arg\min_{\theta}\; J(\theta),\quad J(\theta)=\sum_{k=0}^{N-1}\|r_k(\theta)\|^2. \]

If the model is linear in parameters, \( f_{\text{sim}}(x_k,u_k;\theta)=\Phi_k \theta \) with regressor \( \Phi_k \in \mathbb{R}^{n\times p} \), then stacking data yields \( y = \Phi \theta + \varepsilon \) where \( y = \frac{x_{k+1}-x_k}{T_s} \). The LS solution is closed-form:

\[ \hat{\theta} = (\Phi^\top \Phi)^{-1}\Phi^\top y. \]

Proposition (Unbiasedness and covariance). If \( \varepsilon \) is zero-mean with covariance \( \sigma^2 I \), and \( \Phi^\top \Phi \) is invertible, then

\[ \mathbb{E}[\hat{\theta}] = \theta,\qquad \operatorname{Cov}(\hat{\theta}) = \sigma^2(\Phi^\top \Phi)^{-1}. \]

Proof. Substitute \( y=\Phi\theta+\varepsilon \):

\[ \hat{\theta}=(\Phi^\top\Phi)^{-1}\Phi^\top(\Phi\theta+\varepsilon) =\theta + (\Phi^\top\Phi)^{-1}\Phi^\top\varepsilon. \]

Taking expectation gives \( \mathbb{E}[\hat{\theta}]=\theta \). For covariance:

\[ \operatorname{Cov}(\hat{\theta}) =(\Phi^\top\Phi)^{-1}\Phi^\top \operatorname{Cov}(\varepsilon)\Phi(\Phi^\top\Phi)^{-1} =\sigma^2(\Phi^\top\Phi)^{-1}. \]

After fitting \( \hat{\theta} \), we validate using fresh data. Define prediction error on a validation set: \( e_k = x^{\text{val}}_{k} - x^{\text{sim}}_{k}(\hat{\theta}) \).

Task-level error norms.

\[ \text{RMSE} = \sqrt{\frac{1}{N}\sum_{k=1}^{N}\|e_k\|^2},\qquad \text{MaxErr} = \max_{1\le k\le N}\|e_k\|. \]

A common validation question: “Is the simulator consistent with data within tolerance \( \epsilon \)?” We can formulate

\[ H_0:\; \mathbb{E}[\|e_k\|^2] \le \epsilon^2 \quad \text{vs.}\quad H_1:\; \mathbb{E}[\|e_k\|^2] > \epsilon^2. \]

If we assume \( \|e_k\|^2 \) are i.i.d. with sample mean \( \bar{z} \) and sample variance \( s^2 \), then by the one-sided t-test,

\[ t = \frac{\bar{z}-\epsilon^2}{s/\sqrt{N}} \]

Reject \( H_0 \) at level \( \alpha \) if \( t > t_{1-\alpha,\,N-1} \). This yields a statistically grounded “pass/fail” validation.

Note: validation is always relative to a task and horizon. The bound in Section 3 often motivates choosing short-horizon metrics for unstable or contact-rich tasks.

Consider a 1D actuator with viscous friction: \( \dot{v} = \frac{1}{m}(u - b v) \). We identify \( b \) from real data, then validate.

import numpy as np
from numpy.linalg import lstsq

# "Real" data generation (unknown to identification)
np.random.seed(0)
m_true = 2.0
b_true = 0.8
Ts = 0.01
N = 2000

u = 2.0*np.sin(0.5*np.arange(N)*Ts)  # known input
v = np.zeros(N)
for k in range(N-1):
    vdot = (u[k] - b_true*v[k]) / m_true
    v[k+1] = v[k] + Ts*vdot
v_meas = v + 0.02*np.random.randn(N)  # noisy measurement

# Identification model: v_{k+1} - v_k = Ts*(u_k/m - (b/m) v_k)
# Fix m to known nominal (from datasheet), estimate b by LS.
m_nom = 2.0
y = (v_meas[1:] - v_meas[:-1]) / Ts - u[:-1]/m_nom
Phi = (-v_meas[:-1] / m_nom).reshape(-1,1)
b_hat, _, _, _ = lstsq(Phi, y, rcond=None)
b_hat = float(b_hat)

print("Estimated b:", b_hat, "True b:", b_true)

# Validation on fresh input
Nv = 800
u_val = 1.5*np.cos(0.7*np.arange(Nv)*Ts)
v_real = np.zeros(Nv)
v_sim  = np.zeros(Nv)

for k in range(Nv-1):
    # real (unknown to sim)
    v_real[k+1] = v_real[k] + Ts*(u_val[k] - b_true*v_real[k])/m_true
    # sim with estimated friction
    v_sim[k+1]  = v_sim[k]  + Ts*(u_val[k] - b_hat*v_sim[k])/m_nom

e = v_real - v_sim
rmse = np.sqrt(np.mean(e**2))
maxerr = np.max(np.abs(e))
print("Validation RMSE:", rmse, "Max error:", maxerr)

# t-test against tolerance epsilon
epsilon = 0.05
z = e**2
zbar = np.mean(z)
s = np.std(z, ddof=1)
t_stat = (zbar - epsilon**2) / (s/np.sqrt(Nv))
print("t statistic:", t_stat)
      

This minimal example mirrors a simulator tuning loop: choose a parametric model, fit \( \theta \), then validate on unseen trajectories.

#include <Eigen/Dense>
#include <iostream>
#include <vector>

int main(){
    double Ts = 0.01;
    double m_nom = 2.0;
    int N = 2000;

    std::vector<double> u(N), v_meas(N);
    // Assume u and v_meas are filled from logs (CSV/ROS bag).
    // For demo, create dummy signals:
    for(int k=0;k<N;k++){
        u[k] = 2.0*std::sin(0.5*k*Ts);
        v_meas[k] = 0.5*std::sin(0.5*k*Ts) + 0.02*((rand()%1000)/1000.0 - 0.5);
    }

    Eigen::VectorXd y(N-1);
    Eigen::MatrixXd Phi(N-1,1);

    for(int k=0;k<N-1;k++){
        double dv = (v_meas[k+1]-v_meas[k])/Ts;
        y(k) = dv - u[k]/m_nom;
        Phi(k,0) = -v_meas[k]/m_nom;
    }

    // LS: minimize ||Phi*b - y||^2
    double b_hat = (Phi.transpose()*Phi).ldlt().solve(Phi.transpose()*y)(0);
    std::cout << "Estimated b = " << b_hat << std::endl;
    return 0;
}
      

Eigen provides stable solvers for the normal equations. In practical robot stacks, the same logic is embedded in calibration nodes.

import org.ejml.simple.SimpleMatrix;
import java.util.Random;

public class FrictionLS {
    public static void main(String[] args){
        double Ts = 0.01;
        double mNom = 2.0;
        int N = 2000;
        double[] u = new double[N];
        double[] v = new double[N];
        Random rnd = new Random(0);

        for(int k=0;k<N;k++){
            u[k] = 2.0*Math.sin(0.5*k*Ts);
            v[k] = 0.5*Math.sin(0.5*k*Ts) + 0.02*(rnd.nextDouble()-0.5);
        }

        SimpleMatrix Phi = new SimpleMatrix(N-1,1);
        SimpleMatrix y   = new SimpleMatrix(N-1,1);

        for(int k=0;k<N-1;k++){
            double dv = (v[k+1]-v[k])/Ts;
            y.set(k,0, dv - u[k]/mNom);
            Phi.set(k,0, -v[k]/mNom);
        }

        // b_hat = (Phi^T Phi)^-1 Phi^T y
        double bHat = Phi.transpose().mult(Phi).invert().mult(Phi.transpose()).mult(y).get(0,0);
        System.out.println("Estimated b = " + bHat);
    }
}
      

EJML is a lightweight linear algebra option for robot software on JVM-based stacks.

Matlab code for LS fitting and validation:

Ts = 0.01; N = 2000; m_nom = 2.0;
t = (0:N-1)*Ts;
u = 2*sin(0.5*t);

% Load measured velocity v_meas from experiment/log
% Here we mock it:
v_meas = 0.5*sin(0.5*t) + 0.02*randn(size(t));

dv = diff(v_meas)/Ts;
y = dv - u(1:end-1)/m_nom;
Phi = -v_meas(1:end-1)/m_nom;

b_hat = (Phi*Phi')\(Phi*y');  % scalar LS
disp(["Estimated b = ", num2str(b_hat)]);

% Validation on a step input
Nv = 800; tv = (0:Nv-1)*Ts;
u_val = ones(1,Nv);  % step

v_real = zeros(1,Nv); v_sim = zeros(1,Nv);
b_true = 0.8; m_true = 2.0;

for k=1:Nv-1
    v_real(k+1) = v_real(k) + Ts*(u_val(k)-b_true*v_real(k))/m_true;
    v_sim(k+1)  = v_sim(k)  + Ts*(u_val(k)-b_hat *v_sim(k))/m_nom;
end

e = v_real - v_sim;
rmse = sqrt(mean(e.^2));
maxerr = max(abs(e));
disp(["RMSE = ", num2str(rmse), ", MaxErr = ", num2str(maxerr)]);
      

In Simulink, build two blocks: (1) “Real Plant” and (2) “Sim Model”. Feed identical \( u(t) \), log both outputs, and compute RMSE using a “MATLAB Function” block. This mirrors digital twin validation by co-simulation.

Problem 1 (Gap Bound Derivation): Assume \( \|\Delta f(x,u,t)\|\le \delta \) and \( f_{\text{sim}} \) is Lipschitz with constant \( L \). Starting from \( e(0)=0 \), derive an upper bound on \( \|e(t)\| \).

Solution: From Section 3:

\[ \|\dot{e}\| \le \delta + L\|e\|. \]

Let \( g(t)=\|e(t)\| \), then \( \dot{g} \le \delta + L g \), \( g(0)=0 \). Multiplying by \( e^{-Lt} \) and integrating yields

\[ g(t)\le \frac{\delta}{L}(e^{Lt}-1), \]

hence \( \|e(t)\| \) satisfies the same bound.

Problem 2 (Stable Linear Gap): Let \( \dot{x}=Ax+Bu \) and simulated \( \dot{x}=(A+\Delta A)x + (B+\Delta B)u \). Assuming \( A \) is stable and \( \|e^{At}\|\le M e^{-\alpha t} \), show that \( \sup_{t\ge 0}\|e(t)\| \) is bounded by a constant proportional to \( \|\Delta A\| \) and \( \|\Delta B\| \).

Solution: The error dynamics are \( \dot{e}=Ae+\Delta A x_{\text{sim}}+\Delta B u \). Using variation of constants:

\[ e(t)=\int_0^t e^{A(t-s)}(\Delta A x_{\text{sim}}(s)+\Delta B u(s))\,ds. \]

Taking norms and the exponential bound:

\[ \|e(t)\| \le \int_0^t M e^{-\alpha(t-s)}(\|\Delta A\|X_{\max}+\|\Delta B\|U_{\max})\,ds \le \frac{M}{\alpha}(\|\Delta A\|X_{\max}+\|\Delta B\|U_{\max}), \]

which is finite for stable \( A \).

Problem 3 (Least Squares Properties): For \( y=\Phi\theta+\varepsilon \) with \( \mathbb{E}[\varepsilon]=0 \) and \( \operatorname{Cov}(\varepsilon)=\sigma^2 I \), prove that the LS estimate \( \hat{\theta}=(\Phi^\top\Phi)^{-1}\Phi^\top y \) is unbiased and compute its covariance.

Solution: Substitute:

\[ \hat{\theta}=\theta + (\Phi^\top\Phi)^{-1}\Phi^\top\varepsilon. \]

Taking expectation gives unbiasedness. Covariance follows as

\[ \operatorname{Cov}(\hat{\theta})=\sigma^2(\Phi^\top\Phi)^{-1}. \]

Problem 4 (Validation Threshold Test): Suppose on a validation set you obtain squared errors \( z_k=\|e_k\|^2 \). Design a one-sided statistical test for \( H_0:\mathbb{E}[z_k]\le \epsilon^2 \).

Solution: Use the t-statistic

\[ t=\frac{\bar{z}-\epsilon^2}{s/\sqrt{N}}, \]

where \( \bar{z} \) is the sample mean and \( s \) is the sample standard deviation. Reject at level \( \alpha \) if \( t > t_{1-\alpha,N-1} \).

Problem 5 (Design Choice): A simulator matches real data well for slow trajectories but fails for fast motions. Using the bound in Section 3, give two plausible explanations involving \( L \) and \( \delta \).

Solution: (i) Fast motions increase sensitivity, effectively increasing the Lipschitz constant \( L \), causing exponential amplification of small errors. (ii) At high speeds, unmodeled effects (aerodynamic drag, actuator saturation, friction regime changes) raise \( \delta \), enlarging the bound.

We modeled the sim-to-real gap as a structured difference between real and simulated dynamics, derived trajectory error bounds via Lipschitz arguments, and specialized them to linearized LTI models aligned with your linear control background. We then showed how data-driven parameter identification reduces the gap and how statistical validation decides if a simulator or digital twin is reliable for a given task.