Simulation-Based Testing: From Monte Carlo to Autonomous Vehicles

Simulation-Based Testing: From Monte Carlo to Autonomous Vehicles

Testing in production is expensive. Testing against reality is sometimes impossible. A self-driving car cannot accumulate millions of edge-case miles before launch. A satellite cannot be rebooted mid-orbit. A financial exchange cannot be stress-tested with real market crashes on demand. Simulation fills the gap — and doing it right is an engineering discipline of its own.

This guide covers the core concepts and practical tools engineers use to build simulation-based test suites: discrete event simulation, Monte Carlo methods, robotics simulators like Gazebo, autonomous vehicle environments like CARLA, property-based testing layered over simulations, and deliberate fault injection.

Simulation vs Emulation vs Virtualization

These terms are often used interchangeably and shouldn't be.

Virtualization creates an isolated instance of the same hardware and OS. A VM running Ubuntu on a Linux host is virtualization. The guest and host share the same instruction set. Tests that run in VMs are testing real software in a sandboxed environment — there is no abstraction of the system under test.

Emulation replaces one system's behavior with another at the hardware or firmware level. QEMU emulates ARM hardware on x86. A Game Boy emulator mimics the exact CPU, memory bus, and display controller. Emulation is faithful to the original system but not to physical reality — timing, power consumption, and analog behavior are typically not modeled.

Simulation abstracts a system into a model. The model captures the behaviors that matter for your test and ignores the rest. A traffic simulation does not emulate CPU instructions — it models vehicle physics, road topology, and driver behavior at a level useful for testing a routing algorithm. Simulation trades fidelity for speed and controllability.

For software testing purposes, the relevant distinction is: can you inject arbitrary conditions into your test environment? With simulation you can. With a real system you often cannot.

Discrete Event Simulation Testing

Discrete event simulation (DES) models a system as a sequence of events occurring at points in time. Between events, nothing changes. This makes it extremely efficient for testing systems like queues, message brokers, schedulers, and network protocols.

Python's simpy library is the standard tool for DES in Python.

import simpy
import random
from collections import defaultdict

RANDOM_SEED = 42
random.seed(RANDOM_SEED)

class MessageBroker:
    def __init__(self, env, processing_time_mean=0.5, queue_capacity=100):
        self.env = env
        self.queue = simpy.Resource(env, capacity=1)
        self.processing_time_mean = processing_time_mean
        self.queue_capacity = queue_capacity
        self.metrics = defaultdict(list)
        self.dropped = 0
        self.processed = 0

    def process_message(self, message_id, arrival_time):
        if len(self.queue.queue) >= self.queue_capacity:
            self.dropped += 1
            return

        with self.queue.request() as req:
            yield req
            wait_time = self.env.now - arrival_time
            self.metrics["wait_times"].append(wait_time)

            proc_time = random.expovariate(1.0 / self.processing_time_mean)
            yield self.env.timeout(proc_time)
            self.processed += 1
            self.metrics["processing_times"].append(proc_time)


def message_generator(env, broker, arrival_rate):
    message_id = 0
    while True:
        yield env.timeout(random.expovariate(arrival_rate))
        arrival_time = env.now
        env.process(broker.process_message(message_id, arrival_time))
        message_id += 1


def test_broker_under_load():
    env = simpy.Environment()
    broker = MessageBroker(env, processing_time_mean=0.5, queue_capacity=50)

    # Arrival rate slightly exceeds processing capacity — stress condition
    env.process(message_generator(env, broker, arrival_rate=2.2))
    env.run(until=1000)

    drop_rate = broker.dropped / (broker.processed + broker.dropped)
    avg_wait = sum(broker.metrics["wait_times"]) / len(broker.metrics["wait_times"])

    assert drop_rate < 0.05, f"Drop rate {drop_rate:.2%} exceeds 5% SLA"
    assert avg_wait < 2.0, f"Average wait {avg_wait:.2f}s exceeds 2s SLA"
    print(f"Processed: {broker.processed}, Dropped: {broker.dropped}")
    print(f"Drop rate: {drop_rate:.2%}, Avg wait: {avg_wait:.3f}s")


test_broker_under_load()

The key insight here is that you can test SLAs that would require weeks of real traffic in seconds of simulation time. You can also run the same simulation under different load profiles — 10% overload, 50% overload, sustained burst — by changing a single parameter.

Monte Carlo Simulation in Testing

Monte Carlo methods replace deterministic inputs with random sampling to explore the probability distribution of outcomes. Where DES tests a specific sequence of events, Monte Carlo tests what happens across thousands of random scenarios.

This is especially powerful for testing numerical algorithms, financial models, retry logic, and any system where you care about tail behavior.

import random
import math
from dataclasses import dataclass
from typing import List

@dataclass
class RetryConfig:
    max_attempts: int
    base_delay: float
    max_delay: float
    jitter: bool

def simulate_retry_sequence(config: RetryConfig, failure_probability: float) -> dict:
    """Simulate one request with exponential backoff retry."""
    total_time = 0.0
    for attempt in range(config.max_attempts):
        if random.random() > failure_probability:
            return {"success": True, "attempts": attempt + 1, "total_time": total_time}

        delay = min(config.base_delay * (2 ** attempt), config.max_delay)
        if config.jitter:
            delay *= random.uniform(0.5, 1.5)
        total_time += delay

    return {"success": False, "attempts": config.max_attempts, "total_time": total_time}


def monte_carlo_retry_test(config: RetryConfig, failure_prob: float, trials: int = 10_000):
    results = [simulate_retry_sequence(config, failure_prob) for _ in range(trials)]

    successes = [r for r in results if r["success"]]
    failures = [r for r in results if not r["success"]]
    success_rate = len(successes) / trials
    avg_time_on_success = sum(r["total_time"] for r in successes) / max(len(successes), 1)
    p99_time = sorted(r["total_time"] for r in successes)[int(len(successes) * 0.99)] if successes else 0

    return {
        "success_rate": success_rate,
        "avg_time_on_success": avg_time_on_success,
        "p99_time": p99_time,
        "total_failures": len(failures),
    }


def test_retry_policy_at_50_percent_failure():
    config = RetryConfig(
        max_attempts=5,
        base_delay=0.1,
        max_delay=10.0,
        jitter=True,
    )
    stats = monte_carlo_retry_test(config, failure_prob=0.5, trials=50_000)

    # With 5 attempts and 50% failure rate, success probability = 1 - 0.5^5 = 96.875%
    assert stats["success_rate"] > 0.96, f"Success rate {stats['success_rate']:.3f} too low"
    assert stats["p99_time"] < 30.0, f"P99 latency {stats['p99_time']:.1f}s exceeds budget"
    print(f"Success rate: {stats['success_rate']:.3%}")
    print(f"P99 time on success: {stats['p99_time']:.2f}s")


test_retry_policy_at_50_percent_failure()

Monte Carlo testing catches policies that look fine on average but have catastrophic tail behavior. A retry config with no jitter produces thundering herd effects that only appear under concurrent load — Monte Carlo with correlated failures will surface this before production does.

Gazebo for Robotics Testing

Gazebo (now Gz Sim) is the standard simulation environment for ROS-based robotics. It provides physics simulation, sensor modeling, and plugin APIs that let you write tests against virtual robots running real ROS 2 node code.

A typical test scenario looks like this:

# test_navigation_recovery.py — ROS 2 + Gazebo integration test
import pytest
import rclpy
from rclpy.node import Node
from geometry_msgs.msg import PoseWithCovarianceStamped, Twist
from nav_msgs.msg import OccupancyGrid
import subprocess
import time

@pytest.fixture(scope="module")
def ros_context():
    rclpy.init()
    yield
    rclpy.shutdown()

@pytest.fixture(scope="module")
def gazebo_world():
    """Launch Gazebo with a world containing known obstacles."""
    proc = subprocess.Popen([
        "ros2", "launch", "my_robot_bringup",
        "sim.launch.py", "world:=obstacle_course.world"
    ])
    time.sleep(5)  # Wait for Gazebo and ROS nodes to initialize
    yield proc
    proc.terminate()
    proc.wait()


class NavigationTestNode(Node):
    def __init__(self):
        super().__init__("navigation_test_node")
        self.goal_reached = False
        self.collision_detected = False
        self.pose_sub = self.create_subscription(
            PoseWithCovarianceStamped,
            "/amcl_pose",
            self._pose_callback,
            10,
        )
        self.cmd_pub = self.create_publisher(Twist, "/cmd_vel", 10)

    def _pose_callback(self, msg):
        x = msg.pose.pose.position.x
        y = msg.pose.pose.position.y
        # Goal zone: 5m corridor end
        if abs(x - 5.0) < 0.3 and abs(y - 0.0) < 0.3:
            self.goal_reached = True


def test_robot_navigates_to_goal(ros_context, gazebo_world):
    node = NavigationTestNode()

    # Send navigation goal via action client (simplified here)
    import subprocess
    subprocess.run([
        "ros2", "action", "send_goal",
        "/navigate_to_pose",
        "nav2_msgs/action/NavigateToPose",
        "{pose: {header: {frame_id: map}, pose: {position: {x: 5.0, y: 0.0, z: 0.0}}}}"
    ])

    deadline = time.time() + 60  # 60 second timeout
    while time.time() < deadline:
        rclpy.spin_once(node, timeout_sec=0.1)
        if node.goal_reached:
            break

    assert node.goal_reached, "Robot failed to reach goal within 60 seconds"
    assert not node.collision_detected, "Robot collided with obstacle during navigation"
    node.destroy_node()

The power of Gazebo testing is repeatability. A physical robot navigating a corridor will encounter sensor noise, floor imperfections, and lighting variation that you cannot control. In Gazebo, every run starts from the same initial state. You can introduce specific failure modes — add a dynamic obstacle, introduce GPS spoofing, simulate a stuck wheel — and verify the recovery behavior every time.

CARLA for Autonomous Vehicle Testing

CARLA is an open-source simulator built on Unreal Engine specifically for autonomous driving research and testing. It provides a Python API for scenario scripting, sensor configuration, and traffic control.

import carla
import time
import math

def test_emergency_braking_scenario():
    client = carla.Client("localhost", 2000)
    client.set_timeout(10.0)
    world = client.get_world()

    # Deterministic simulation: fixed delta time
    settings = world.get_settings()
    settings.synchronous_mode = True
    settings.fixed_delta_seconds = 0.05
    world.apply_settings(settings)

    blueprint_library = world.get_blueprint_library()

    # Spawn ego vehicle
    ego_bp = blueprint_library.find("vehicle.tesla.model3")
    spawn_point = carla.Transform(carla.Location(x=0, y=0, z=0.5))
    ego_vehicle = world.spawn_actor(ego_bp, spawn_point)

    # Attach collision sensor
    collision_bp = blueprint_library.find("sensor.other.collision")
    collision_sensor = world.spawn_actor(
        collision_bp,
        carla.Transform(),
        attach_to=ego_vehicle,
    )
    collisions = []
    collision_sensor.listen(lambda event: collisions.append(event))

    # Spawn obstacle 30m ahead, stationary
    obstacle_bp = blueprint_library.find("vehicle.audi.a2")
    obstacle_spawn = carla.Transform(carla.Location(x=30, y=0, z=0.5))
    obstacle = world.spawn_actor(obstacle_bp, obstacle_spawn)

    # Apply throttle to ego vehicle — AEB should trigger
    ego_vehicle.apply_control(carla.VehicleControl(throttle=0.8, brake=0.0))

    # Run simulation for 5 seconds
    for _ in range(100):  # 100 ticks × 0.05s = 5s
        world.tick()

    final_velocity = ego_vehicle.get_velocity()
    speed = math.sqrt(final_velocity.x**2 + final_velocity.y**2)
    ego_location = ego_vehicle.get_location()

    # Cleanup
    collision_sensor.destroy()
    ego_vehicle.destroy()
    obstacle.destroy()

    assert len(collisions) == 0, f"Collision detected — AEB failed ({len(collisions)} events)"
    assert speed < 1.0, f"Vehicle still moving at {speed:.1f} m/s — AEB did not engage"
    print(f"AEB test passed. Final speed: {speed:.2f} m/s, position: x={ego_location.x:.1f}")


test_emergency_braking_scenario()

CARLA's synchronous mode is critical for testing: it ensures the simulation advances exactly one tick per world.tick() call, giving you full determinism. Without this, sensor callbacks and physics updates arrive at non-deterministic times, making assertions unreliable.

Property-Based Testing with Simulations

Property-based testing (Hypothesis in Python) generates random inputs and checks that invariants hold. When combined with simulation, it becomes a powerful way to explore the input space of a system model.

from hypothesis import given, settings, assume
from hypothesis import strategies as st
import math

def pid_controller_step(setpoint, current, prev_error, integral, kp, ki, kd, dt):
    error = setpoint - current
    integral += error * dt
    derivative = (error - prev_error) / dt
    output = kp * error + ki * integral + kd * derivative
    return output, error, integral


def simulate_pid(kp, ki, kd, setpoint, initial_value, steps=200, dt=0.05):
    current = initial_value
    prev_error = setpoint - current
    integral = 0.0
    history = [current]

    for _ in range(steps):
        output, prev_error, integral = pid_controller_step(
            setpoint, current, prev_error, integral, kp, ki, kd, dt
        )
        # First-order system response
        current += (output - current) * dt
        history.append(current)

    return history


@given(
    kp=st.floats(min_value=0.1, max_value=5.0),
    ki=st.floats(min_value=0.0, max_value=1.0),
    kd=st.floats(min_value=0.0, max_value=2.0),
    setpoint=st.floats(min_value=-10.0, max_value=10.0),
)
@settings(max_examples=500)
def test_pid_stability_property(kp, ki, kd, setpoint):
    assume(not math.isnan(kp) and not math.isnan(ki) and not math.isnan(kd))

    history = simulate_pid(kp, ki, kd, setpoint, initial_value=0.0)

    # Property 1: No NaN or Inf values (numerical stability)
    assert all(math.isfinite(v) for v in history), "PID produced non-finite output"

    # Property 2: Output bounded to [-100, 100] (actuator saturation check)
    assert all(abs(v) <= 100.0 for v in history), f"Output exceeded bounds: {max(abs(v) for v in history):.1f}"

    # Property 3: Converges within 10% of setpoint by end of simulation
    final_values = history[-20:]
    tolerance = max(abs(setpoint) * 0.1, 0.5)
    final_error = abs(sum(final_values) / len(final_values) - setpoint)
    # Only assert convergence for stable-gain combinations
    if kp < 2.0 and ki < 0.5:
        assert final_error < tolerance, f"PID did not converge: final error {final_error:.3f}, tolerance {tolerance:.3f}"

This approach finds gain combinations that cause instability, overflow, or divergence — issues that would only appear in field deployment if you only test with hand-picked parameters.

Fault Injection

Fault injection is the deliberate introduction of failures into a running system or simulation to verify that error-handling, recovery, and degraded-mode behaviors work correctly. This is distinct from testing the happy path — you are specifically testing what happens when things go wrong.

import random
import time
from unittest.mock import patch, MagicMock
from contextlib import contextmanager

class FaultInjector:
    def __init__(self):
        self.active_faults = {}

    @contextmanager
    def network_partition(self, failure_rate=1.0, latency_ms=None):
        """Inject network failures into any requests made within this context."""
        original_send = None  # Would patch actual HTTP library in real usage

        def faulty_send(*args, **kwargs):
            if random.random() < failure_rate:
                raise ConnectionError("Simulated network partition")
            if latency_ms:
                time.sleep(latency_ms / 1000.0)
            return original_send(*args, **kwargs)

        yield faulty_send

    @contextmanager
    def disk_full(self, target_module, write_function="write"):
        """Simulate disk-full errors on file writes."""
        with patch(f"{target_module}.{write_function}") as mock_write:
            mock_write.side_effect = OSError(28, "No space left on device")
            yield mock_write

    @contextmanager
    def memory_pressure(self, allocation_failure_rate=0.1):
        """Randomly fail memory allocations."""
        original_list = list

        def faulty_list(*args, **kwargs):
            if args and isinstance(args[0], range) and len(args[0]) > 1000:
                if random.random() < allocation_failure_rate:
                    raise MemoryError("Simulated OOM")
            return original_list(*args, **kwargs)

        with patch("builtins.list", faulty_list):
            yield


class DataPipeline:
    def __init__(self, storage_backend):
        self.storage = storage_backend
        self.checkpoint = None

    def process_batch(self, records):
        results = []
        for i, record in enumerate(records):
            try:
                processed = self._transform(record)
                self.storage.write(processed)
                results.append(processed)
                if i % 10 == 0:
                    self.checkpoint = i  # Save progress
            except OSError:
                # Should resume from checkpoint, not lose all progress
                return {"processed": results, "checkpoint": self.checkpoint, "error": "disk_full"}
        return {"processed": results, "checkpoint": len(records), "error": None}

    def _transform(self, record):
        return {**record, "processed": True}


def test_pipeline_resumes_from_checkpoint_after_disk_full():
    injector = FaultInjector()
    mock_storage = MagicMock()
    pipeline = DataPipeline(mock_storage)

    records = [{"id": i, "value": i * 2} for i in range(25)]

    # Storage fails after 15 writes
    write_count = [0]
    def counting_write(data):
        write_count[0] += 1
        if write_count[0] > 15:
            raise OSError(28, "No space left on device")

    mock_storage.write.side_effect = counting_write
    result = pipeline.process_batch(records)

    assert result["error"] == "disk_full"
    assert result["checkpoint"] is not None, "Pipeline must save checkpoint before failure"
    assert result["checkpoint"] >= 10, f"Checkpoint {result['checkpoint']} is too early — lost too much progress"
    assert len(result["processed"]) > 0, "Pipeline returned no results before failure"

    print(f"Checkpoint at record {result['checkpoint']}, {len(result['processed'])} records processed before failure")


test_pipeline_resumes_from_checkpoint_after_disk_full()

Fault injection testing should cover the full taxonomy of failures: network (partition, latency, packet loss), storage (disk full, permission denied, corruption), process (OOM, timeout, crash), and dependency (downstream service unavailable, returning malformed data). For each failure mode, the question is not "does the system fail?" but "does it fail in a controlled, recoverable way?"

Putting It Together

Simulation-based testing occupies a specific niche in your test pyramid. Unit tests verify logic. Integration tests verify contracts. Simulation tests verify behavior under conditions you cannot reliably reproduce in any other environment.

The practical rule: reach for simulation when:

  1. The real system is too expensive or dangerous to test against — vehicles, satellites, medical devices, financial markets.
  2. You need statistical confidence about tail behavior — Monte Carlo for retry policies, SLA compliance under load.
  3. You need determinism for debugging — replaying the exact sensor sequence that caused a robotics failure.
  4. You need to test failure modes deliberately — fault injection requires control that production systems don't provide.

The tools covered here — simpy for DES, Hypothesis for property-based exploration, Gazebo and CARLA for physical simulations, and manual fault injection — each address a different class of problem. The common thread is that you are building a model, and the quality of your tests depends directly on the fidelity of that model to the system it represents.

Model too simple: your tests pass but the system fails in production on unmodeled behaviors. Model too complex: the simulation itself has bugs and becomes another system to maintain. The engineering judgment is knowing which simplifications are safe.

Start with the behaviors that keep you up at night, and simulate those first.

Start now free