Design-Based Simulations: Methods & Applications

• Design-based simulations are methodologies that treat observed outcomes as fixed while resampling design elements to generate controlled simulation data.
• They integrate design of experiments, adaptive sampling, and surrogate modeling to efficiently bridge expensive computational models with finite-sample inference.
• These techniques enable rapid optimization and robust analysis across econometrics, engineering design, and complex simulation-driven workflows.

Searching arXiv for recent and foundational papers on "design-based simulations" and closely related simulation-driven design frameworks. Design-based simulations are procedures in which the realized outcomes are treated as fixed, and randomness in the simulated data is generated solely by resampling elements of the research design—typically treatment assignments or exogenous shocks (Ferman, 11 Mar 2026). In a broader simulation-and-design literature, the design problem is the choice of locations in parameter space at which simulations are to be run, the choice of fidelity levels, or the choice of geometry and operating conditions from which an emulator, surrogate, or optimizer is constructed [(Schneider et al., 2010); (Savage et al., 2023)]. Across these uses, design-based simulations connect expensive computational models to finite-sample inference, adaptive sampling, inverse design, and interactive engineering workflows. This suggests a family of methodologies organized around deliberate control of what is held fixed, what is varied, and how simulation budget is allocated.

1. Scope and conceptual foundations

In the narrow econometric sense, a design-based simulation starts from an observed dataset obeying a model such as $Y_i=\alpha+\beta T_i+\varepsilon_i$, fixes the realized vector $Y=(Y_1,\dots,Y_N)$, imposes a null hypothesis such as $H_0:\beta=0$, and repeatedly draws new treatment vectors $T^*$ from a known or assumed assignment distribution. The empirical rejection frequency of the resulting test statistics is then interpreted as an estimate of finite-sample size or power in the original application (Ferman, 11 Mar 2026).

In the broader emulator literature, the central object is not the treatment assignment alone but the simulation campaign itself. “Simulation design” is the choice of $n$ design points $$\Theta=\{\theta_1,\dots,\theta_n\}\subset\mathbb{R}^p$$ at which to run the full simulator so that an emulator can predict the simulator output $y(\theta)$ at arbitrary $\theta\in\Theta$. Under that formulation, the problem is to minimize emulator error given a fixed budget $n$ (Schneider et al., 2010).

A third formulation appears in Bayesian optimal design. There, one chooses a design $\xi\in\Xi$ to maximize an expected utility

$Y=(Y_1,\dots,Y_N)$0

and the optimal design is $Y=(Y_1,\dots,Y_N)$1. When $Y=(Y_1,\dots,Y_N)$2 is unavailable in closed form but simulation from the model is possible, approximate Bayesian computation can be incorporated into the design algorithm (Hainy et al., 2013).

Taken together, these formulations show that design-based simulations are not limited to one discipline. They include resampling-based inferential assessments, emulator-oriented placement of expensive runs, and utility-driven design search under intractable likelihoods.

2. Planning simulation studies and placing simulation runs

A recurring theme is that simulation studies themselves should be designed rather than merely tabulated. In Design and Analysis of Experiments, the planning stage uses fishbone diagrams to enumerate candidate factors, then chooses factors and levels, full or fractional factorial designs, and, when robustness is the target, Taguchi Robust Parameter Design. Analysis is subsequently carried out via ANOVA, main-effect and interaction plots, half-normal plots, and signal-to-noise ratios such as $Y=(Y_1,\dots,Y_N)$3, $Y=(Y_1,\dots,Y_N)$4, and $Y=(Y_1,\dots,Y_N)$5 (Chipman et al., 2021).

For emulator construction, Orthogonal-Array-Based Latin Hypercube Sampling provides a standard design in which runs are spread across one- and two-dimensional projections. Schneider and Knox proposed a Fisher-whitened modification, OALHSFS, that restricts points to a “Fisher sphere.” In the six-parameter CMB TT power-spectrum case study, OALHSFS halved the number $Y=(Y_1,\dots,Y_N)$6 of runs needed to reach a given emulator RMSE, and for $Y=(Y_1,\dots,Y_N)$7 gave $Y=(Y_1,\dots,Y_N)$8 smaller errors than standard OALHS. In that same study, quadratic polynomial interpolation performed nearly as well as Gaussian-process interpolation for $Y=(Y_1,\dots,Y_N)$9, with interpolation errors $H_0:\beta=0$0 over $H_0:\beta=0$1 up to 2000; the “Emu CMB” implementation achieved $H_0:\beta=0$2 interpolation errors at all $H_0:\beta=0$3 up to 5000 after a log-space principal-component decomposition and quadratic interpolation of the first 20 modes (Schneider et al., 2010).

The same design logic extends to likelihood-free settings. Müller’s MCMC formulation treats the joint integrand $H_0:\beta=0$4 as an unnormalized density whose marginal in $H_0:\beta=0$5 is proportional to $H_0:\beta=0$6. Simulated-annealing-style augmentation with $H_0:\beta=0$7 independent replicates concentrates the marginal law of $H_0:\beta=0$8 around the optimizer, while two-stage ABC plus design-MCMC replaces prior draws of $H_0:\beta=0$9 with draws from an ABC approximation to $T^*$0 (Hainy et al., 2013).

These strands share a common methodological commitment: simulation runs are scarce, and their placement determines what can later be inferred, interpolated, or optimized.

3. Surrogate models and learned simulators

A major contemporary form of design-based simulation replaces expensive solvers with learned surrogates. In “Graph Neural Network Based Surrogate Model of Physics Simulations for Geometry Design” (Wong et al., 2023), any 2D or 3D unstructured simulation mesh is cast as a graph $T^*$1, with nodes corresponding to mesh points and edges corresponding to bidirectional mesh connectivity. Node features include relative coordinates, norms of those offsets, one-hot problem flags, and topological information such as node degree; edge features include relative displacement vectors and their $T^*$2 norms. The model uses an encoder–processor–decoder architecture with message passing, graph-level latent aggregation $T^*$3, and node-level outputs $T^*$4. On additive-manufacturing feature design, the surrogate obtained Train/Test1 median relative $T^*$5-norm error $T^*$6 with a range $T^*$7–$T^*$8, while high-fidelity LPBF simulation required $T^*$9–$n$0 hours per geometry and the GNN surrogate required $n$1 second per geometry, yielding a $n$2–$n$3 speedup. On the airfoil task, pressure prediction had Test median $n$4, drag $n$5 had $n$6, lift $n$7 had $n$8, and $n$9 for $$\Theta=\{\theta_1,\dots,\theta_n\}\subset\mathbb{R}^p$$0 and $$\Theta=\{\theta_1,\dots,\theta_n\}\subset\mathbb{R}^p$$1.

SimuLearn adopts a similar replacement strategy for morphing materials, but couples finite-element data generation to a graph-convolutional neural simulator. Abaqus generates 4,377 training samples, each taking $$\Theta=\{\theta_1,\dots,\theta_n\}\subset\mathbb{R}^p$$2 min on an 8 core, 5 GHz CPU; the learned simulator then predicts in $$\Theta=\{\theta_1,\dots,\theta_n\}\subset\mathbb{R}^p$$3 s per simulation on CPU, with a $$\Theta=\{\theta_1,\dots,\theta_n\}\subset\mathbb{R}^p$$4 speedup versus Abaqus FEA at $$\Theta=\{\theta_1,\dots,\theta_n\}\subset\mathbb{R}^p$$5 s, mean vertex error $$\Theta=\{\theta_1,\dots,\theta_n\}\subset\mathbb{R}^p$$6 mm, and physical-prototype deviation $$\Theta=\{\theta_1,\dots,\theta_n\}\subset\mathbb{R}^p$$7, corresponding to overall $$\Theta=\{\theta_1,\dots,\theta_n\}\subset\mathbb{R}^p$$8 accuracy. The framework supports forward, inverse, and hybrid CAD workflows in Rhino + Grasshopper, including 100 parallel design variants in $$\Theta=\{\theta_1,\dots,\theta_n\}\subset\mathbb{R}^p$$9 s in inverse mode (Yang et al., 2020).

In structural engineering, “Learning to simulate and design for structural engineering” (Chang et al., 2020) uses a pre-trained differentiable structural simulator, NeuralSim, together with a learned proposal network, NeuralSizer. NeuralSim predicts story drifts from graph-encoded buildings in 6.8 ms per structure, compared with RSA simulation at $y(\theta)$0 s, for a $y(\theta)$1 speedup and $y(\theta)$2 relative accuracy on test. NeuralSizer then proposes cross-sections under mass, drift-ratio, variety, and entropy constraints, and its designs are comparable to those optimized by genetic algorithm while satisfying all constraints.

The significance of these systems is not merely acceleration. They change the granularity of feasible design interaction: from hours or days per geometry to sub-second or millisecond evaluation, which enables iterative design loops that would otherwise be computationally impractical.

4. Differentiable optimization, active learning, and multi-fidelity querying

A stronger integration of design and simulation occurs when gradients are propagated through the simulator itself. In “Physical Design using Differentiable Learned Simulators” (Allen et al., 2022), the inverse-design objective is posed as

$y(\theta)$3

where $y(\theta)$4 is a ground-truth simulator, $y(\theta)$5 is a learned GNN approximation, and $y(\theta)$6 maps design parameters and uncontrollable conditions into an initial state that is then rolled forward for $y(\theta)$7 steps. Because $y(\theta)$8, $y(\theta)$9, and $\theta\in\Theta$0 are fully differentiable, gradients are obtained by back-propagation through time. On 2D fluid tools, GD-M was approximately $\theta\in\Theta$1–$\theta\in\Theta$2 better than CEM-M and $\theta\in\Theta$3–$\theta\in\Theta$4 better than CEM-S when evaluated with the ground-truth simulator. On the airfoil task, GD-M with a single model reached $\theta\in\Theta$5, an ensemble of 5 reached $\theta\in\Theta$6, and the DAFoam specialized adjoint solver achieved $\theta\in\Theta$7; runtimes were 21 s for a single model, 62 s for the ensemble, and 1021 s on 8-core CPU for DAFoam.

When the expensive object is not a trajectory simulator but a portfolio of candidate runs, active learning can select which simulations to perform. In offshore riser design, pool-based active learning represents 526 candidate loading cases as an unlabeled pool $\theta\in\Theta$8, initializes with $\theta\in\Theta$9 simulations, and repeatedly queries the point with largest GP predictive uncertainty. With a total budget $n$0, the method used approximately $n$1 fewer runs than the full $n$2–526 set, achieved $n$3 for empty-fill DNVUF-201 and $n$4 N for $n$5, and delivered a speed-up factor $n$6. Random sampling at 100 runs produced $n$7 larger errors (Elsas et al., 2020).

A related strategy is continuous multi-fidelity optimization. DARTS models reactor performance $n$8 jointly over geometry, operating conditions, and continuous mesh fidelities $n$9, and then optimizes a cost-adjusted UCB acquisition that trades off expected gain at high fidelity against evaluation cost and posterior correlation to the highest-fidelity point. After approximately 64 h of wall-clock budget, DARTS returned an optimal helical-tube design with $\xi\in\Xi$0, $\xi\in\Xi$1, $\xi\in\Xi$2, and pulsed-flow $\xi\in\Xi$3, giving $\xi\in\Xi$4. Experimental RTDs from a 3D-printed reactor yielded $\xi\in\Xi$5, in excellent agreement with the predicted RTD (Savage et al., 2023).

High-fidelity fluid–structure interaction design uses a similar GP logic but places stronger emphasis on verification and bridging. The workflow in “A Principled Approach to Design Using High Fidelity Fluid-Structure Interaction Simulations” proceeds from solver verification on Turek–Hron FSI3, through bridge simulations that vary elasticity under simplified geometries, to constrained Bayesian optimization for a UUAV sail plane. The final optimization used 8 initial and 8 BO iterations, for 16 FSI calls total, reducing drag from $\xi\in\Xi$6 N at the mid-range baseline to $\xi\in\Xi$7 N at $\xi\in\Xi$8, a $\xi\in\Xi$9 reduction while satisfying $Y=(Y_1,\dots,Y_N)$00 N, $Y=(Y_1,\dots,Y_N)$01 mm, and $Y=(Y_1,\dots,Y_N)$02 kPa (Wu et al., 2020).

At the solver level, open-source inverse design for 3D nanostructures has also been pushed into this design-based regime. FDTDX implements FDTD in JAX with reverse-mode automatic differentiation derived from time reversibility in Maxwell’s equations. By storing only the PML boundary slices and reconstructing interior fields during the backward pass, memory scales as $Y=(Y_1,\dots,Y_N)$03 rather than the naive $Y=(Y_1,\dots,Y_N)$04. In the silicon waveguide bend example, the best device reached $Y=(Y_1,\dots,Y_N)$05 dB attenuation, or $Y=(Y_1,\dots,Y_N)$06 efficiency; in the 3D polymer stitching device, coupling was maintained between $Y=(Y_1,\dots,Y_N)$07 dB and $Y=(Y_1,\dots,Y_N)$08 dB over random $Y=(Y_1,\dots,Y_N)$09 (Mahlau et al., 2024).

5. Domain-specific embodiments

Classical CFD-based engineering design remains an important embodiment of design-based simulation. For hydraulic Cross-Flow turbines, the methodology in (Mehr et al., 2017) separates a one-step system-level design from a three-step detail-design using ANSYS CFX: nozzle geometry optimization, runner parameter optimization, and performance evaluation under variable loads. The resulting turbine achieves peak hydraulic efficiency $Y=(Y_1,\dots,Y_N)$10 and peak overall efficiency $Y=(Y_1,\dots,Y_N)$11, with $Y=(Y_1,\dots,Y_N)$12 over $Y=(Y_1,\dots,Y_N)$13 head variation and down to $Y=(Y_1,\dots,Y_N)$14 of nominal flow. Because the optimized geometry is parameterized by dimensionless similarity, efficiency maps $Y=(Y_1,\dots,Y_N)$15 can be reused across a turbine family.

Nanophotonic design problems expose the same issues in a different physical regime. The benchmark framework of (Kim et al., 2023) formulates parametric structure design over film stacks, nanocones, nanowires, nanospheres, film-plus-cone structures, and combinatorial material-block systems, with FDTD simulations at low, medium, and high fidelity. The framework explicitly supports a discretized-search-space mode, a surrogate-model mode, and a direct-simulation mode. A recommended workflow is coarse search at low fidelity, screening at medium fidelity, and final ranking at high fidelity. Example high-fidelity runtimes ranged from $Y=(Y_1,\dots,Y_N)$16 s for a three-layer film to $Y=(Y_1,\dots,Y_N)$17 s for close-packed nanospheres.

Detector design offers another long-standing use. In plastic scintillator studies, GEANT4 simulations vary the volume-to-area ratio $Y=(Y_1,\dots,Y_N)$18, coating thickness $Y=(Y_1,\dots,Y_N)$19, reflectivity $Y=(Y_1,\dots,Y_N)$20, surface roughness, and scintillation decay profile. Standard coating with $Y=(Y_1,\dots,Y_N)$21 boosts collection sufficiently that the overall detected fraction $Y=(Y_1,\dots,Y_N)$22 is approximately $Y=(Y_1,\dots,Y_N)$23, $Y=(Y_1,\dots,Y_N)$24, and $Y=(Y_1,\dots,Y_N)$25 for $Y=(Y_1,\dots,Y_N)$26, whereas without coating these values fall by factors $Y=(Y_1,\dots,Y_N)$27–$Y=(Y_1,\dots,Y_N)$28. Small volumes yield fast and narrow pulses; large volumes maximize deposited energy at high energies but suffer strong optical losses (Ros et al., 2018). In atmospheric muon tomography, an end-to-end CORSIKA+GEANT4 pipeline simulates ground-level muon flux, detector response, and object interaction. Using a $Y=(Y_1,\dots,Y_N)$29 threshold, lead can be differentiated from aluminum in 4.9 days by absorption and 9.9 days by scattering, and the absorption method gives the best results (Rengifo et al., 2024).

A newer application is LLM-based social simulation. “The Silicon Society Cookbook” varies seven design parameters: number of agents, base LLM, network topology, homophily, survey-in-context, biased news agent, and persona proportions. The reported design space has a “hybrid geometry”: some factors behave additively, while others show strong interactions. The base LLM is the dominant factor for both stylistic realism and opinion-shift direction, with single-factor $Y=(Y_1,\dots,Y_N)$30 for BERT accuracy and $Y=(Y_1,\dots,Y_N)$31 for Net Consensus Change. Survey context raises detectability additively, but its effect on consensus drift depends on model choice, and the $Y=(Y_1,\dots,Y_N)$32 interaction on Majority Follow Rate is negligible for $Y=(Y_1,\dots,Y_N)$33 and strong for $Y=(Y_1,\dots,Y_N)$34 (Bück-Kaeffer et al., 30 Apr 2026).

These examples indicate that design-based simulations are not tied to one solver family or one ontology of “design.” The design object may be a nozzle curvature, a mesh fidelity, a detector coating, a treatment assignment rule, or the base model in a simulated social network.

6. Validity, misalignment, and recurring limitations

A common misconception is that any simulation that “looks like” the application is automatically informative about inference or design validity. The econometric literature rejects that view explicitly. In shift-share designs, standard simulations that fix outcomes and resample shocks can rely on a working data-generating process that is not aligned with the true one. They can confound true treatment effects with error dependence and thereby overstate inference distortions due to spatial correlation. Alternative procedures therefore either fix estimated residual structure while resampling shocks or jointly simulate shocks and spatially correlated errors (Ferman, 11 Mar 2026).

The engineering literature reports analogous limits in generalization. The GNN surrogate for geometry design performs well in distribution, but out-of-distribution brackets and stems produce errors between $Y=(Y_1,\dots,Y_N)$35 and $Y=(Y_1,\dots,Y_N)$36, and the paper explicitly lists reduced accuracy on heavily out-of-distribution shapes or extreme operating points as a limitation. SimuLearn likewise notes that collisions are not modeled, that accuracy is limited by FEA fidelity and training-set coverage, and that larger out-of-distribution geometries may degrade performance (Wong et al., 2023, Yang et al., 2020).

Multi-fidelity and active-learning methods also impose their own conditions. In nanophotonic FDTD benchmarks, lower fidelity produces noisier, biased spectra, whereas high fidelity smooths out interference fringes and converges to reference. In offshore riser design, GP training scales as $Y=(Y_1,\dots,Y_N)$37, even though $Y=(Y_1,\dots,Y_N)$38 remained cheap in the reported study. The quality of uncertainty-driven selection also depends on kernel adequacy and domain-specific preprocessing (Kim et al., 2023, Elsas et al., 2020).

A broader implication is that design-based simulations are most reliable when three alignments hold simultaneously: the simulated data-generating process matches the inferential target, the design space covers the operating regime of interest, and the fidelity hierarchy preserves the ranking or qualitative structure needed for optimization. Where these alignments fail, the literature typically recommends either redesigning the simulation protocol, enriching the training distribution, or validating the final design with a higher-fidelity or physical experiment.