A robot simulation is a computational surrogate of a physical robot and its environment. In control terms, the physical robot is a dynamical system (the plant) and the simulator is an approximate model of that plant, driven by the same inputs.

Let the physical robot state be \( x(t) \in \mathbb{R}^n \), control input be \( u(t) \in \mathbb{R}^m \), and measured output be \( y(t) \in \mathbb{R}^p \). A general continuous-time plant can be written as

\[ \dot{x}(t) = f\big(x(t),u(t),w(t)\big), \quad y(t) = h\big(x(t)\big) + v(t), \]

where \( w(t) \) represents unmodeled disturbances (contacts, friction variability, aerodynamic drag, etc.) and \( v(t) \) is sensor noise. A simulator replaces \( f \) and \( h \) with approximations \( \hat{f} \), \( \hat{h} \):

\[ \dot{\hat{x}}(t) = \hat{f}\big(\hat{x}(t),u(t)\big), \quad \hat{y}(t) = \hat{h}\big(\hat{x}(t)\big). \]

The simulator may be deterministic or stochastic, and may include rigid-body dynamics, actuation models, sensor models, and a physics engine for contact and constraints.

flowchart TD
  A["Design goal or control idea"] --> B["Build model of robot + environment"]
  B --> C["Run simulation experiments"]
  C --> D["Analyze failures and performance"]
  D --> E["Refine controller / design"]
  E --> C
  C --> F["Deploy on physical robot"]
        

Define the simulation error as \( e(t) := x(t) - \hat{x}(t) \). Subtracting the two systems gives

\[ \dot{e}(t) = f\big(x(t),u(t),w(t)\big) - \hat{f}\big(\hat{x}(t),u(t)\big). \]

Add and subtract \( f(\hat{x}(t),u(t),0) \):

\[ \dot{e}(t) = \underbrace{f\big(x(t),u(t),w(t)\big) - f\big(\hat{x}(t),u(t),0\big)}_{\text{state + disturbance mismatch}} + \underbrace{f\big(\hat{x}(t),u(t),0\big) - \hat{f}\big(\hat{x}(t),u(t)\big)}_{\text{model mismatch}}. \]

Assume (a standard control assumption) that \( f \) is Lipschitz in state with constant \( L > 0 \) and that disturbances are bounded: \( \|w(t)\| \le w_{\max} \). Also assume the model mismatch is bounded: \( \|f(\hat{x},u,0)-\hat{f}(\hat{x},u)\| \le \delta \). Then

\[ \|\dot{e}(t)\| \le L\|e(t)\| + L_w w_{\max} + \delta \]

for some disturbance gain \( L_w \). By Grönwall’s inequality,

\[ \|e(t)\| \le e^{Lt}\|e(0)\| + \int_0^t e^{L(t-s)}(L_w w_{\max}+\delta)\,ds = e^{Lt}\|e(0)\| + \frac{e^{Lt}-1}{L}\,(L_w w_{\max}+\delta). \]

Interpretation: even with perfect initial alignment \( e(0)=0 \), mismatch accumulates with time and is proportional to two unavoidable terms: (i) real-world disturbances \( w_{\max} \) and (ii) imperfect modeling \( \delta \).

Simulation errors come from (a) modeling error and (b) numerical error. Modeling error is structural: wrong friction law, inaccurate mass distribution, simplified contact model, etc. Numerical error is algorithmic: finite time steps, iterative solvers, and discretization artifacts.

flowchart RL
  P["Physical robot"] -->|unknown disturbances| W["w(t)"]
  P -->|true dynamics f| X["x(t)"]
  S["Simulator"] -->|approx dynamics f_hat| XH["x_hat(t)"]
  X --> E["error e(t)=x-x_hat"]
  XH --> E
  W --> E
  M["model mismatch delta"] --> E
  N["numerical error"] --> E
        

3.1 Numerical integration error (discretization)

Simulators usually integrate continuous dynamics in discrete time with step \( h \). For simplicity, consider the general ODE \( \dot{x} = f(x,u) \). Explicit Euler gives

\[ \hat{x}_{k+1} = \hat{x}_k + h\,f(\hat{x}_k,u_k). \]

If \( f \) is continuously differentiable, the local truncation error is \( O(h^2) \), so the global error after \( T/h \) steps is \( O(h) \).

Proof sketch (global error order for Euler):

By Taylor expansion of the true solution about \( t_k \),

\[ x(t_{k+1}) = x(t_k) + h\,f(x(t_k),u_k) + \frac{h^2}{2}\,\ddot{x}(\xi_k), \]

where \( \xi_k \in (t_k,t_{k+1}) \). Subtract Euler’s update:

\[ x(t_{k+1})-\hat{x}_{k+1} = \big(x(t_k)-\hat{x}_k\big) + h\big(f(x(t_k),u_k)-f(\hat{x}_k,u_k)\big) + O(h^2). \]

Using Lipschitz continuity: \( \|f(x)-f(\hat{x})\| \le L\|x-\hat{x}\| \),

\[ \|e_{k+1}\| \le (1+hL)\|e_k\| + Ch^2. \]

Unrolling this recursion for \( k \approx T/h \) yields \( \|e_k\| \le O(h) \). Hence Euler is first-order accurate.

3.2 Stability limits of discrete simulation

Students already know linear systems. Consider a stable continuous-time system \( \dot{x}=Ax \) with eigenvalues \( \lambda_i(A) \) satisfying \( \Re(\lambda_i)<0 \). Euler discretization gives

\[ \hat{x}_{k+1} = (I+hA)\hat{x}_k. \]

Discrete-time stability requires spectral radius \( \rho(I+hA) < 1 \). For each eigenvalue \( \lambda \), the mapped eigenvalue is \( 1+h\lambda \). Thus stability requires

\[ |1+h\lambda| < 1 \quad \text{for all eigenvalues } \lambda. \]

For real negative \( \lambda \), this becomes \( -2 < h\lambda < 0 \Rightarrow 0 < h < 2/|\lambda| \). So a simulator can become unstable purely due to a large timestep, even when the physical robot is stable.

3.3 Contact and friction complexity

Many robots interact with the environment through impacts and Coulomb-like friction. A common ideal model is

\[ f_{\text{fric}} = -\mu N\,\operatorname{sign}(v), \]

but real contacts exhibit stiction, velocity-dependent friction, deformation, and microslip. Simulators replace these with smooth or piecewise models to remain numerically solvable. This introduces structural mismatch \( \delta \) in Section 2.

4.1 Safety and cost

Real robots are expensive and can be dangerous. In simulation, we can explore wide regions of state space without physical risk. If a controller is parameterized by \( \theta \) and performance is measured by a cost functional \( J(\theta) \), we can evaluate many candidates:

\[ \theta^\star_{\text{sim}} = \arg\min_{\theta} J_{\text{sim}}(\theta), \]

before testing only a small, safe subset on hardware.

4.2 Scalability and reproducibility

Simulation enables parallel experimentation. If \( J_{\text{sim}}(\theta;\omega) \) depends on random environment parameters \( \omega \) (terrain, friction, sensor noise), then we can estimate expected cost with Monte Carlo:

\[ \mathbb{E}[J_{\text{sim}}(\theta)] \approx \frac{1}{N}\sum_{i=1}^N J_{\text{sim}}(\theta;\omega_i). \]

4.3 Exploring uncertainty (domain randomization)

Since real-world parameters are uncertain, a robust design solves

\[ \theta^\star_{\text{rob}} = \arg\min_{\theta}\; \mathbb{E}_{\omega}[J_{\text{sim}}(\theta;\omega)]. \]

This reduces sensitivity to mismatch \( \delta \) when transferring to hardware.

4.4 “Ground-truth” access

In simulation we can access latent variables (exact state, contact forces) that are hard to measure on real robots. This is essential for debugging and learning.

To keep within previously introduced topics, we simulate a 1D mass–spring–damper, which is a standard actuator/load model in robotics. The continuous dynamics are

\[ m\ddot{q} + b\dot{q} + k q = u(t), \]

written as a first-order system with state \( x = \begin{bmatrix} q \\ \dot{q} \end{bmatrix} \):

\[ \dot{x} = \begin{bmatrix} 0 & 1 \\ -k/m & -b/m \end{bmatrix}x + \begin{bmatrix} 0 \\ 1/m \end{bmatrix}u. \]

We compare Euler simulation for different timesteps to show numerical limits.

5.1 Python (NumPy / SciPy)

import numpy as np

m, b, k = 1.0, 0.6, 4.0
A = np.array([[0.0, 1.0], [-k/m, -b/m]])
B = np.array([[0.0], [1.0/m]])

def euler_sim(h, T=5.0):
    N = int(T/h)
    x = np.zeros((2, N+1))
    u = lambda t: 1.0  # unit step
    for i in range(N):
        x[:, i+1] = x[:, i] + h*(A@x[:, i] + (B*u(i*h)).ravel())
    t = np.linspace(0, T, N+1)
    return t, x

for h in [0.05, 0.2, 0.5]:
    t, x = euler_sim(h)
    print("h =", h, "final q =", x[0, -1])
      

Related robotics simulation ecosystems in Python include PyBullet, MuJoCo Python API, Isaac Sim Python, and ROS2 Gazebo interfaces, which we will overview later.

5.2 C++ (from-scratch Euler; libraries: Bullet, ODE, DART, Gazebo)

#include <iostream>
#include <vector>

int main() {
    double m=1.0, b=0.6, k=4.0;
    double h=0.2, T=5.0;
    int N = (int)(T/h);

    std::vector<double> q(N+1,0.0), dq(N+1,0.0);
    auto u = [](double t){ return 1.0; };

    for(int i=0;i<N;i++){
        double ddq = (u(i*h) - b*dq[i] - k*q[i])/m;
        dq[i+1] = dq[i] + h*ddq;
        q[i+1]  = q[i]  + h*dq[i];
    }
    std::cout << "final q = " << q[N] << std::endl;
    return 0;
}
      

5.3 Java (simple ODE loop; robotics libs: ROSJava, jMonkeyEngine Physics)

public class MassSpringDamper {
    public static void main(String[] args){
        double m=1.0, b=0.6, k=4.0;
        double h=0.2, T=5.0;
        int N = (int)(T/h);

        double q=0.0, dq=0.0;
        for(int i=0;i<N;i++){
            double u = 1.0;
            double ddq = (u - b*dq - k*q)/m;
            dq += h*ddq;
            q  += h*dq;
        }
        System.out.println("final q = " + q);
    }
}
      

5.4 Matlab / Simulink (script + toolbox notes)

m = 1.0; b = 0.6; k = 4.0;
A = [0 1; -k/m -b/m];
B = [0; 1/m];

h = 0.2; T = 5.0;
N = floor(T/h);

x = zeros(2, N+1);
u = @(t) 1.0;

for i = 1:N
    x(:,i+1) = x(:,i) + h*(A*x(:,i) + B*u((i-1)*h));
end
disp(['final q = ', num2str(x(1,end))]);
      

In Simulink, the same system is a second-order block diagram. For robotics-scale models, Simscape Multibody and the Robotics System Toolbox provide standard components and URDF import; those are used in later lessons.

Problem 1 (Bounding Simulation Error): Suppose the physical robot follows \( \dot{x} = f(x,u,w) \) and the simulator uses \( \dot{\hat{x}} = \hat{f}(\hat{x},u) \). Assume \( f \) is Lipschitz in \( x \) with constant \( L \), disturbances satisfy \( \|w\|\le w_{\max} \), and model mismatch satisfies \( \|f(\hat{x},u,0)-\hat{f}(\hat{x},u)\|\le \delta \). Show that \( \|e(t)\| \le \frac{e^{Lt}-1}{L}(L_w w_{\max}+\delta) \) when \( e(0)=0 \).

Solution:

From Section 2, \( \|\dot{e}(t)\| \le L\|e(t)\| + L_w w_{\max} + \delta \). With \( e(0)=0 \), apply Grönwall:

\[ \|e(t)\| \le \int_0^t e^{L(t-s)}(L_w w_{\max}+\delta)\,ds = \frac{e^{Lt}-1}{L}(L_w w_{\max}+\delta). \]

Problem 2 (Euler Stability Condition): Let a stable scalar system be \( \dot{x} = \lambda x \) with \( \lambda < 0 \). Explicit Euler gives \( x_{k+1}=(1+h\lambda)x_k \). Derive the step-size condition for discrete stability.

Solution:

Stability requires \( |1+h\lambda|<1 \). Because \( \lambda < 0 \), this is equivalent to

\[ -1 < 1+h\lambda < 1 \;\Rightarrow\; -2 < h\lambda < 0 \;\Rightarrow\; 0 < h < \frac{2}{|\lambda|}. \]

Problem 3 (Order of Global Error): For Euler simulation of \( \dot{x}=f(x) \), show that a local truncation error \( O(h^2) \) implies global error \( O(h) \).

Solution:

Let one-step error be \( \tau_k = x(t_{k+1}) - x(t_k) - h f(x(t_k)) \), with \( \|\tau_k\|\le Ch^2 \). Then

\[ e_{k+1} = e_k + h\big(f(x(t_k))-f(\hat{x}_k)\big) + \tau_k. \]

Using Lipschitz continuity, \( \|e_{k+1}\|\le (1+hL)\|e_k\|+Ch^2 \). Iterating for \( k=T/h \) steps yields a sum of \( O(h^2) \) terms over \( O(1/h) \) steps, giving total \( O(h) \).

Problem 4 (Why Randomize Parameters?): Consider optimizing a controller parameter \( \theta \) in simulation. Explain using formulas why minimizing only \( J_{\text{sim}}(\theta;\omega_0) \) for a single nominal environment \( \omega_0 \) can fail in real robots, and why minimizing \( \mathbb{E}_{\omega}[J_{\text{sim}}(\theta;\omega)] \) can help.

Solution:

Real deployment corresponds to a distribution of conditions \( \omega \sim p(\omega) \). A controller tuned to a single \( \omega_0 \) solves

\[ \theta^\star_{\text{nom}}=\arg\min_{\theta} J_{\text{sim}}(\theta;\omega_0), \]

which need not minimize performance under \( p(\omega) \). Minimizing the expectation

\[ \theta^\star_{\text{rob}}=\arg\min_{\theta}\mathbb{E}_{\omega\sim p(\omega)}[J_{\text{sim}}(\theta;\omega)] \]

reduces sensitivity to mismatch by spreading optimization over likely real conditions.

We modeled simulation as an approximate dynamical surrogate of a physical robot, derived bounds on sim-to-real error, and highlighted two major limitations: structural modeling mismatch and numerical discretization error. We proved Euler’s global error order and its timestep stability restriction, then discussed benefits such as safety, scalability, uncertainty exploration, and access to ground-truth variables. The mini-lab connected these ideas to executable simulations in Python, C++, Java, and Matlab.