Simulation Optimization

Overview

Simulation optimization finds the best input parameters for a system that can only be evaluated through simulation. Unlike classical optimization where the objective function has a closed form, here each function evaluation requires running a stochastic simulation — expensive, noisy, and often black-box.

┌────────────┐    parameters    ┌────────────┐    noisy output
│ Optimization│───────────────►│ Simulation  │──────────────►  statistics
│  Algorithm  │◄───────────────│   Model     │
└────────────┘   performance    └────────────┘

Ranking and Selection

Given k alternative system configurations, identify the best (highest or lowest mean performance) with a statistical guarantee.

Indifference-Zone Procedures

Guarantee: select the best system with probability >= 1 - alpha, whenever the best is at least delta better than the second best.

Procedure (Rinott's two-stage):

1. Take n_0 replications of each system
2. Compute sample means and variances
3. Calculate total replications needed per system: n_i = max(n_0, ceil((h·S_i / delta)²))
4. Take additional replications as needed
5. Select system with best sample mean

// Select the system with the smallest mean from simulation results
PROCEDURE RANKING_AND_SELECTION(means, std_errs, delta) → index
    // Simple: return index of minimum mean
    // In practice: apply Rinott or KN procedure for guarantees
    best_index ← 0
    FOR i ← 1 TO LENGTH(means) - 1 DO
        IF means[i] < means[best_index] THEN
            best_index ← i
    RETURN best_index

Procedures for Many Systems

• KN procedure (Kim-Nelson): Sequential elimination. After each batch of replications, eliminate systems that are statistically dominated.
• OCBA (Optimal Computing Budget Allocation): Allocate replications proportionally to uncertainty about which system is best.

Response Surface Methodology (RSM)

Approximate the simulation response with a surrogate model, then optimize the surrogate.

Workflow

1. Screening — Identify influential factors (fractional factorial, Morris method)
2. First-order model — Fit y = beta_0 + sum(beta_i·x_i) + epsilon using a 2^k factorial design
3. Steepest ascent/descent — Follow the gradient of the fitted model
4. Second-order model — Near the optimum, fit y = beta_0 + sum(beta_i·x_i) + sum(beta_ii·x_i²) + sum(beta_ij·x_i·x_j) using central composite or Box-Behnken design
5. Optimize — Find the stationary point of the quadratic surface

Metamodel-Based Optimization

Modern alternatives to classical RSM:

• Kriging / Gaussian Process (GP) — Interpolating surrogate with uncertainty estimates. Enables Bayesian optimization.
• Radial Basis Functions — Flexible interpolation for high-dimensional spaces.
• Neural network surrogates — For very complex response surfaces with large training budgets.

Bayesian Optimization

Use a GP surrogate with an acquisition function to balance exploration and exploitation:

• Expected Improvement (EI): E[max(f(x) - f_best, 0)]
• Upper Confidence Bound (UCB): mu(x) + kappa·sigma(x)
• Knowledge Gradient: value of information from one additional sample

Particularly effective when simulation is very expensive (minutes to hours per run) and the parameter space is moderate-dimensional (< ~20).

Simulation-Based Optimization

Stochastic Approximation

Robbins-Monro: iteratively update parameters along noisy gradient estimates:

theta_{n+1} = theta_n - a_n · g_hat(theta_n)

where g_hat is a noisy gradient estimate and {a_n} is a diminishing step-size sequence (a_n -> 0, sum(a_n) = inf, sum(a_n²) < inf).

Stochastic Gradient Methods

• Finite differences: Estimate gradient by perturbing each parameter. Cost: 2d simulations per iteration for d parameters.
• Simultaneous Perturbation (SPSA): Perturb all parameters at once with random direction. Cost: 2 simulations per iteration regardless of d.
• Likelihood ratio / score function: Differentiate the expectation analytically through the probability distribution.

Metaheuristics

For combinatorial or non-smooth problems:

• Genetic algorithms — Population of solutions evolved through selection, crossover, mutation
• Simulated annealing — Random perturbations accepted with probability depending on temperature schedule
• Particle swarm — Solutions move through parameter space influenced by personal and global best
• Cross-entropy method — Iteratively update sampling distribution toward elite solutions

Reinforcement Learning for Simulation

RL treats the simulation as an environment. An agent learns a policy mapping states to actions to maximize cumulative reward.

Simulation as RL Environment

┌─────────┐  action   ┌────────────┐  next state, reward
│ RL Agent│──────────►│ Simulation │──────────────────────►
│ (Policy)│◄──────────│Environment │
└─────────┘ observation└────────────┘

Advantages over classical simulation optimization:

• Handles sequential decision-making naturally
• Learns adaptive policies, not just static parameter settings
• Generalizes across states (neural network function approximation)

Challenges:

• Sample efficiency — millions of simulation episodes may be needed
• Sim-to-real gap — policy trained in simulation may fail in reality
• Reward shaping — defining the right reward signal is critical

Applications

• Inventory management — RL agent learns reorder policies
• Traffic signal control — Optimize signal timing based on real-time state
• Manufacturing scheduling — Sequence jobs across machines
• Robot control — Learn locomotion or manipulation in physics simulator

Digital Twins

A digital twin is a live virtual replica of a physical asset, process, or system, continuously synchronized with real-world data.

Architecture

┌─────────────────────────────────────────────┐
│              Physical System                 │
│  sensors ──► data stream ──► actuator cmds   │
└──────┬──────────────────────────────┬────────┘
       │ IoT / SCADA                  ▲
       ▼                              │
┌──────────────────────────────────────────────┐
│              Digital Twin Platform           │
│  ┌──────────┐  ┌──────────┐  ┌───────────┐  │
│  │  Data    │  │Simulation│  │ Analytics │  │
│  │Ingestion │─►│  Model   │─►│& Decision │  │
│  └──────────┘  └──────────┘  └───────────┘  │
│  ┌──────────┐  ┌──────────┐  ┌───────────┐  │
│  │  State   │  │ What-If  │  │Visualization│ │
│  │Estimation│  │ Scenarios│  │    & UI   │  │
│  └──────────┘  └──────────┘  └───────────┘  │
└──────────────────────────────────────────────┘

Synchronization

The twin must reflect current physical state:

• State estimation — Kalman filters, particle filters, or data assimilation techniques fuse sensor data with model predictions
• Parameter updating — Online calibration adjusts model parameters as the physical system drifts
• Update frequency — Ranges from real-time (milliseconds) for control to periodic (hourly/daily) for strategic planning

Predictive Maintenance

A primary digital twin application. The twin tracks degradation state and predicts remaining useful life (RUL):

1. Ingest sensor data (vibration, temperature, current)
2. Update degradation model (physics-based, data-driven, or hybrid)
3. Simulate forward to predict failure time distribution
4. Optimize maintenance schedule: minimize cost = f(downtime, repair, spare parts)

// Simple degradation model with digital twin update
STRUCTURE DegradationTwin
    degradation: real       // current estimated degradation level
    drift_rate: real        // degradation rate parameter
    threshold: real         // failure threshold
    process_noise: real
    measurement_noise: real
    estimation_var: real    // Kalman filter variance

// Kalman filter update with new sensor reading
PROCEDURE UPDATE(twin, measurement, dt)
    // Predict
    predicted ← twin.degradation + twin.drift_rate * dt
    predicted_var ← twin.estimation_var + twin.process_noise * dt

    // Update
    kalman_gain ← predicted_var / (predicted_var + twin.measurement_noise)
    twin.degradation ← predicted + kalman_gain * (measurement - predicted)
    twin.estimation_var ← (1.0 - kalman_gain) * predicted_var

// Estimate remaining useful life
PROCEDURE PREDICT_RUL(twin) → real
    remaining ← twin.threshold - twin.degradation
    IF remaining ≤ 0.0 THEN RETURN 0.0
    RETURN remaining / twin.drift_rate

Parallel and Distributed Simulation

Parallel Discrete-Event Simulation (PDES)

Partition the simulation model into Logical Processes (LPs) executing on separate processors.

Conservative approach (Chandy-Misra-Bryant):

• LPs only process events when safe (no earlier event can arrive from other LPs)
• Uses null messages or lookahead to avoid deadlock
• No rollback but may have idle time waiting for safety guarantees

Optimistic approach (Time Warp):

• LPs process events speculatively
• If a straggler event arrives with earlier timestamp, rollback via state restoration and anti-messages
• Global Virtual Time (GVT) provides a floor below which no rollback can occur — enables fossil collection

Scaling Considerations

Cloud and HPC Simulation

• Embarrassingly parallel: Independent replications across nodes (trivial parallelism)
• Parameter sweeps: Each node runs a different parameter combination
• Large-scale PDES: Federated simulations using HLA (High Level Architecture) standard
• GPU acceleration: Thousands of agents or Monte Carlo paths computed in parallel on GPU

Distributed Simulation Standards

• HLA (IEEE 1516): Federates communicate through a Runtime Infrastructure (RTI). Used in defense and large-scale training simulations.
• DIS (IEEE 1278): Protocol for real-time military simulations. Simpler than HLA, uses PDUs (Protocol Data Units) over UDP.

Choosing an Approach