2. Design of experiments for CAE

In CAE-driven sensitivity analysis and response surface optimization, the design of experiments (DoE) is the single decision that most strongly constrains what you can later learn. A solver run may take minutes to hours; the sampling plan determines whether a modest budget of runs yields a well-conditioned surrogate or a matrix of redundant, clustered points. This chapter treats DoE as an engineering allocation problem: choose sample locations so that the information gained per simulation is maximized, given the physics, the parameter ranges, and the intended downstream model.

Why CAE Sampling Differs from Physical Experimentation

Classical DoE was built for physical trials where replication estimates pure error and randomization protects against nuisance factors. In simulation, the "experiment" is deterministic: rerunning identical inputs returns identical outputs, so replication buys nothing except solver-noise diagnostics. There is no block structure, no randomization constraint, and no cost asymmetry between factor levels. What remains is the core problem — where to place points in a k-dimensional box — plus two CAE-specific realities: runs are expensive and often batched, and the response may be strongly nonlinear or discontinuous (contact status changes, buckling, mesh-dependent stress concentrations).

Consequently, CAE DoE practice favors designs that (1) cover the interior of the domain rather than only corners, (2) remain well-conditioned when some runs fail and must be dropped, and (3) can be extended sequentially without discarding prior results.

Screening Designs: Finding the Few Factors That Matter

When a model carries 10–50 input parameters but only a handful dominate the response, screening precedes any response surface work. The objective is ranking, not accurate prediction, so the design must separate main effects with the fewest runs.

Plackett–Burman designs are the workhorse: for N runs (a multiple of 4), up to N−1 factors can be estimated, with main effects aliased against two-factor interactions. A 12-run PB design screens 11 factors at two levels. The practical caveat is severe: if a strong interaction exists, its effect is folded into a main effect and the ranking can mislead. Use PB only when interaction effects are expected to be small relative to main effects, and always confirm the shortlist with a follow-up design.

Fractional factorial (2^k−p) designs offer a more transparent alias structure. A resolution IV design keeps main effects clear of two-factor interactions (though interactions alias with each other); resolution V keeps both main effects and two-factor interactions clear. For a 7-factor study, a 2⁷⁻³_IV design needs 16 runs and cleanly ranks main effects — a defensible screening choice when you can afford slightly more than PB.

Two CAE-specific cautions apply to both. First, two-level screening assumes monotonic behavior across the range; if a parameter has an interior optimum (e.g., rib thickness where mass and stiffness trade off), a two-level design will report a near-zero main effect and wrongly discard it. Second, screening ranges must be physically credible — a range that drives the solver into contact chatter or mesh distortion produces failed runs that corrupt the effect estimates.

Space-Filling Latin Hypercube Sampling

Once the factor shortlist is fixed and the goal shifts to building a surrogate (response surface, Kriging, or polynomial chaos), the design should fill the domain. Latin hypercube sampling (LHS) stratifies each of the k factors into N equal-probability intervals and places exactly one sample per interval per factor, then randomly pairs the strata across factors. This guarantees one-dimensional projection uniformity: every factor is sampled across its full range, with no clustering.

Plain LHS, however, says nothing about multi-dimensional spacing — two points can be arbitrarily close in the joint space. The standard remedy is maximin LHS: generate many candidate LHS matrices (e.g., 1,000), compute the minimum pairwise Euclidean distance in the normalized [0,1]^k space, and retain the candidate with the largest minimum distance. This is cheap to implement and dramatically improves conditioning. An alternative criterion, centered L₂ discrepancy, measures uniformity of the full point set; minimizing it produces designs that are slightly more regular than maximin but can be less robust to failed runs.

# Maximin LHS for k=5 factors, N=60 runs (Python, scipy + numpy)
import numpy as np
from scipy.stats import qmc

k, N, candidates = 5, 60, 1000
best_d, best_X = -1.0, None
for seed in range(candidates):
    X = qmc.LatinHypercube(d=k, seed=seed).random(N)   # rows in [0,1]^k
    d = qmc.distance.pdist(X).min()                     # min pairwise distance
    if d > best_d:
        best_d, best_X = d, X
# Map normalized samples to physical ranges
lo = np.array([2.0, 1.5, 3.0, 0.8, 20.0])   # e.g. thickness, fillet, ... 
hi = np.array([6.0, 4.0, 8.0, 2.5, 120.0])
X_phys = lo + best_X * (hi - lo)
print("maximin min-distance:", best_d)

For a five-factor, 60-run study this typically yields a minimum normalized distance around 0.15–0.20, versus roughly 0.05–0.08 for unoptimized LHS — a threefold improvement in spacing that directly reduces surrogate prediction variance between samples. Note that LHS is not sequentially extensible: adding runs requires regenerating the whole matrix, so plan the final budget before launching the batch.

Optimal Designs: When the Criterion Must Match the Model

When the surrogate form is known in advance — a quadratic response surface, a Kriging covariance, or a polynomial chaos basis — the design can be optimized against that model's information matrix. For a linear model y = Xβ + ε, the Fisher information is X^TX, and three classical criteria follow:

For CAE response surfaces intended for prediction, I-optimal is usually the better choice than D-optimal, because D-optimality concentrates runs at the design extremes where leverage is highest, leaving the interior — where the optimizer will search — thinly sampled. A practical workflow: define a candidate set (e.g., a 3^k or 5^k grid, or a large maximin LHS), specify the quadratic model, and run a coordinate-exchange or Fedorov exchange algorithm to select N rows minimizing the average prediction variance.

Choosing Among the Three Families

The decision is driven by the question being asked, not by software defaults. Screening designs answer "which factors matter?" and tolerate aliasing in exchange for few runs. Space-filling LHS answers "what does the response surface look like across the whole domain?" and is model-agnostic. Optimal designs answer "where should points go to best fit this specific model?" and require committing to the model form up front.

The scatter comparison above makes the trade-off concrete. Maximin LHS spreads points across the interior; Plackett–Burman collapses onto two levels per factor, which is efficient for ranking but useless for curvature; D-optimal for a quadratic model pushes points to corners and edge midpoints, maximizing leverage but leaving the center sparse. If the optimizer will later search the interior, that sparsity is a liability — which is why I-optimal or a hybrid (LHS augmented with optimal points) is often the pragmatic CAE choice.

Practical Execution and Failure Symptoms

A defensible CAE DoE workflow proceeds as follows. Normalize all factors to [0,1] using physically justified bounds, and record the mapping so results can be de-normalized. Generate the design, then inspect it before launching runs: check the minimum pairwise distance, the per-factor histogram (should be flat for LHS), and the correlation matrix between factors (should be near zero). Launch the batch, and log solver status per run. If more than roughly 10% of runs fail, do not simply drop them — the surviving design is no longer space-filling, and the surrogate will be biased toward the feasible region. Instead, tighten the offending factor bounds and regenerate.

Common failure symptoms and their causes: a surrogate that predicts poorly in the interior despite low leave-one-out error usually indicates a screening design was reused for surface fitting (two-level points cannot resolve curvature). A design matrix with high pairwise factor correlation indicates a botched LHS implementation where strata were not independently permuted. Runs that cluster at one corner after de-normalization point to a unit-conversion error in the bounds mapping. Each of these is detectable from the design matrix alone, before spending a single solver hour — which is precisely why design inspection belongs in the workflow, not after the results arrive.