4. Design of experiments for expensive simulations

When a single simulation run costs minutes to hours of wall-clock time, the sampling plan becomes the most consequential decision in the entire surrogate-modeling workflow. A poorly placed set of design points wastes budget on redundant information; a well-designed set extracts maximum coverage from a small number of runs. This section covers the three dominant strategies for expensive black-box functions: space-filling designs, Latin hypercube sampling (LHS), and adaptive sampling.

Why Space-Filling Matters

Classical designs such as full factorials or central composite designs were developed for physical experiments where replication is cheap and the response surface is assumed smooth and low-order. For a crash simulation or a CFD sweep, neither assumption holds. The response can be highly nonlinear, the input space may have 5–20 dimensions, and replication is prohibitively expensive. The goal shifts from estimating variance to covering the input domain uniformly so that no region is left unobserved.

Space-filling designs aim to distribute n points across a d-dimensional hypercube such that the minimum pairwise distance is maximized (maximin criterion) or the discrepancy from a uniform distribution is minimized. The maximin criterion is the most common practical choice:

maximize  min_{i≠j}  ||x_i − x_j||₂
subject to  x_i ∈ [0,1]^d,  i = 1..n

This is a combinatorial optimization problem. For small n and d, exchange algorithms or simulated annealing can solve it directly. For larger budgets, LHS provides a structured starting point that already guarantees one-dimensional stratification.

Latin Hypercube Sampling

LHS partitions each of the d dimensions into n equal-probability intervals and places exactly one point in each interval per dimension. The pairing of coordinates across dimensions is then randomized. This guarantees that the projection of the design onto any single axis is perfectly stratified, which is a significant improvement over simple random sampling.

However, plain LHS does not guarantee good coverage in higher-dimensional projections. Two points can be very close in the full space while still satisfying the one-dimensional stratification. The standard remedy is to optimize the LHS by maximizing the minimum distance or minimizing a correlation-based criterion. A widely used implementation is the lhs package in R or scipy.stats.qmc.LatinHypercube in Python with an optimization step:

from scipy.stats.qmc import LatinHypercube
import numpy as np

sampler = LatinHypercube(d=6, seed=42)
sample = sampler.random(n=40)

# Maximin optimization via random swaps
def maximin(sample):
    from scipy.spatial.distance import pdist
    return pdist(sample).min()

best = sample.copy()
best_score = maximin(best)
for _ in range(5000):
    i, j = np.random.choice(len(sample), 2, replace=False)
    candidate = best.copy()
    candidate[i], candidate[j] = candidate[j], candidate[i]
    score = maximin(candidate)
    if score > best_score:
        best, best_score = candidate, score

In practice, 40–80 points for a 6-dimensional problem is a reasonable starting budget. The maximin-optimized LHS typically improves the minimum pairwise distance by 15–30% over plain LHS, which translates directly into better surrogate accuracy near the boundaries.

Rule of thumb: Use at least 10·d points for an initial LHS when the response is expected to be smooth. For highly nonlinear responses (e.g., crash pulses, combustion), start with 20·d and plan for adaptive refinement.

Adaptive Sampling Strategies

Space-filling designs are one-shot: all points are chosen before any simulation is run. Adaptive sampling, by contrast, uses the surrogate model's own uncertainty estimate to decide where to sample next. The workflow is iterative:

  1. Build an initial space-filling design (e.g., 20·d LHS points).
  2. Fit a surrogate model — typically a Gaussian process (Kriging) because it provides both a prediction and a variance estimate.
  3. Select the next point by maximizing an acquisition function.
  4. Run the expensive simulation at that point, augment the dataset, and refit.
  5. Repeat until the budget is exhausted or the maximum predicted variance falls below a threshold.

The acquisition function balances exploration (sampling where uncertainty is high) and exploitation (sampling where the prediction is promising). Three common choices are:

Acquisition FunctionFormulaBehavior
Maximum varianceσ(x)Pure exploration; good for global surrogate accuracy
Expected improvement (EI)E[max(0, y* − y(x))]Balances improvement and uncertainty; good for optimization
Lower confidence bound (LCB)μ(x) − κσ(x)Tunable κ controls exploration/exploitation trade-off

For surrogate accuracy rather than optimization, maximum variance is often the most direct choice. For design optimization, EI is preferred because it naturally concentrates samples near the optimum while still exploring uncertain regions.

Practical Workflow and Diagnostics

A reproducible adaptive sampling loop for a 6-dimensional crashworthiness simulation might proceed as follows. The inputs are normalized to [0,1]^6: front rail thickness, material yield stress, bumper beam gauge, crash speed, barrier offset, and spot-weld density. The output is peak deceleration (g).

# Pseudocode for adaptive Kriging loop
X = maximin_lhs(n=120, d=6)
y = run_simulation(X)          # expensive: ~45 min per run
for iteration in range(30):
    gp = fit_gaussian_process(X, y)
    candidates = sobol_sequence(10000, d=6)
    sigma = gp.predict_variance(candidates)
    x_next = candidates[argmax(sigma)]
    y_next = run_simulation(x_next)
    X = vstack([X, x_next])
    y = append(y, y_next)
    if sigma.max() < 0.05 * y.std():
        break

Observable diagnostics during the loop include the maximum predicted variance (should decrease monotonically), the leave-one-out cross-validation R² (should approach 0.95+ for a well-behaved response), and the distribution of new sample locations (should spread across the domain rather than cluster). If the maximum variance plateaus, the surrogate has likely reached the noise floor of the simulation or the kernel hyperparameters need re-tuning.

Common failure mode: If the initial LHS is too small, the Gaussian process hyperparameters are poorly estimated, and the variance field is misleading. The adaptive loop then samples in the wrong regions. Always verify that the initial design has at least 10·d points and that the length-scale estimates are stable before trusting the acquisition function.

Choosing Between Strategies

The decision depends on the goal. If the objective is a globally accurate surrogate for visualization or sensitivity analysis, a maximin-optimized LHS with 20·d points is often sufficient and requires no iteration. If the objective is optimization or reliability assessment, adaptive sampling with EI or LCB will typically find better designs with 30–50% fewer simulation runs. In practice, a hybrid approach works well: start with a maximin LHS for robustness, then switch to adaptive sampling once the surrogate has stabilized.

Sampling Strategy Decision Flow for Expensive Simulations Initial Design Maximin LHS n = 10·d to 20·d Fit Surrogate Gaussian Process (Kriging) Acquisition Function max σ(x) or EI(x) or LCB(x) Run Simulation at x_next, augment data iterate Convergence Check σ_max < 0.05·σ_y ? Stop Surrogate ready for use no Initial design Surrogate fit Acquisition Simulation run Convergence

The diagram above summarizes the adaptive sampling loop. The initial maximin LHS provides a robust starting point; the Gaussian process supplies both prediction and uncertainty; the acquisition function selects the next simulation point; and the convergence check terminates the loop when the maximum predicted variance drops below a fraction of the response standard deviation. This workflow is directly applicable to any expensive simulation where a surrogate model is needed for optimization, sensitivity analysis, or reliability assessment.

← Previous 📖 Contents Next →