Written and maintained byCASRAI Editorial Board
Surrogate models—also referred to as metamodels, response surface approximations, or computer experiment emulators—are data-driven mathematical models designed to approximate computationally demanding numerical simulations. In computational engineering and scientific research, executing high-fidelity solvers such as computational fluid dynamics (CFD), finite element analysis (FEA), and multi-physics partial differential equation (PDE) integrations often requires hours or days per design iteration. A surrogate model learns the underlying mapping between input parameters and simulation outputs from a structured set of training runs, allowing researchers to evaluate new design configurations in milliseconds.
What it is used for in research
Surrogate modeling is widely adopted across scientific computing, aerospace engineering, materials science, and physical systems modeling to replace repetitive calls to expensive simulators. Key research applications include:
- Design space exploration and global optimization: Optimization algorithms (such as genetic algorithms, particle swarm optimization, or Bayesian optimization) require thousands of objective function evaluations. Substituting a surrogate model enables rapid convergence without overburdening high-performance computing (HPC) clusters.
- Uncertainty quantification (UQ) and sensitivity analysis: Calculating Sobol sensitivity indices or propagating parameter uncertainties via Monte Carlo sampling demands tens of thousands of evaluations, which is computationally prohibitive with raw PDE solvers but trivial with a trained emulator.
- Digital twins and real-time control: Physical systems require fast feedback loops for hardware-in-the-loop testing and real-time monitoring, where response times must remain below tens of milliseconds.
- Multi-fidelity data integration: When limited high-fidelity simulation budgets exist, multi-fidelity surrogates bridge sparse high-accuracy data with abundant low-fidelity approximations (such as coarse-grid simulations).
Unlike simple polynomial curve fitting implemented through baseline routines like NumPy polyfit, modern scientific surrogate modeling relies on spatial statistics, Gaussian process regression (Kriging), and machine learning methods capable of capturing complex non-linear response surfaces and quantifying prediction uncertainty.
Install or access (official
The Surrogate Modeling Toolbox (SMT) is an open-h Aerospace Lab), ISAE-SUPAERO, NASA Glenn Research Center, and the University of Michigan Multidisciplinary Design Optimization (MDO) Lab. SMT provides specialized surrogate modeling methods, sampling techniques, and gradient-enhanced algorithms tailored for engineering design
SMT requires Python 3.8 or newer along with core scientific libraries: NumPy, SciPy, scikit-learn, pyDOE, and Cython. Certain internal operations (including radial basis function distance evaluations) rely on compiled C++ extensions.
Install the latest release from the official Python Package Index (PyPI):
pip install smt
To verify the installation and check the installed SMT version from an interactive Python shell or terminal session, execute:
import smt
print(smt.__version__)
For research projects tracking code reproducibility, ensure all dependencies and version pins are logged in a requirements file alongside the appropriate open
Step-by-step task (numbered, with menu names or commands exactly as in the current docs)
This workflow demonstrates constructing, training, and evaluating a Kriging Gaussian process surrogate model on a two-dimensional non-linear benchmark using SMT.
- Configure design variable bounds and generate samples via Latin Hypercube Sampling:
Use SMT’s
LHSclass to produce space-filling sample configurations across the parameter domain. Unlike unconstrained uniform sampling discussed in random sampling in Python, Latin Hypercube Sampling ensures stratified coverage across each dimension.import numpy as np from smt.sampling_methods import LHS from smt.surrogate_models import KRG # Define bounds: 2 design variables, x0 in [0, 4], x1 in [0, 3] xlimits = np.array([[0.0, 4.0], [0.0, 3.0]]) # Initialize LHS with Enhanced Stochastic Evolutionary (ESE) optimization sampling = LHS(xlimits=xlimits, criterion='ese', random_state=42) num_samples = 30 xt = sampling(num_samples) print("Training input shape:", xt.shape) # (30, 2) - Evaluate the simulation objective function:
Compute response values
ytfor each sampled design point. In real applications, this step invokes a numerical simulation executable or finite element batch script.# Analytical non-linear response surface for demonstration def simulation_solver(x): # Non-linear function with coupling between variables return ((x[:, 0] - 2.0)**2 + (x[:, 1] - 1.5)**2 + np.sin(2.0 * x[:, 0]) * np.cos(2.0 * x[:, 1])).reshape(-1, 1) yt = simulation_solver(xt) print("Training output shape:", yt.shape) # (30, 1) - Instantiate the Kriging surrogate model:
Create a
KRGsurrogate model instance. Specify the polynomial trend component (poly='constant') and the correlation kernel (corr='matern52'for second-order differentiability orcorr='squar_exp'for infinite smoothness).# Initialize Kriging model sm = KRG(poly='constant', corr='matern52', print_global=False) # Supply training inputs and target responses sm.set_training_values(xt, yt) - Train the surrogate model:
Invoke train() to optimize hyperparameters (such as correlation lengths θ)
# Optimize hyperparameters and fit the Gaussian process sm.train() print("Kriging model training completed successfully.") - Evaluate predictions and compute uncertainty on new test points:
Generate an independent test grid to measure predictive accuracy and compute spatial prediction variances.
# Create evaluation grid eval_sampling = LHS(xlimits=xlimits, random_state=99) x_test = eval_sampling(100) # Predict mean response values and prediction variances y_pred = sm.predict_values(x_test) var_pred = sm.predict_variances(x_test) y_true = simulation_solver(x_test) # Calculate Root Mean Squared Error (RMSE) rmse = np.sqrt(np.mean((y_true - y_pred)**2)) r2 = 1.0 - (np.sum((y_true - y_pred)**2) / np.sum((y_true - np.mean(y_true))**2)) print(f"Validation RMSE: {rmse:.4f}") print(f"Validation R2 score: {r2:.4f}")
Output and how to read it
The Kriging model outputs two primary arrays when evaluated on unobserved query coordinates:
- Predicted response (
predict_values): A 2D NumPy array of shape(n_eval, n_outputs)representing the surrogate’s best estimate ŷ(x). When plotted against the true simulator outputs, points should cluster tightly along the 45° line of identity. An R² value above 0.95 indicates that the surrogate captures the dominant physical trends of the simulation. - Prediction variance (
predict_variances): A 2D NumPy array of shape(n_eval, n_outputs)representing the spatial epistemic uncertainty σ²(x). Because Kriging is an exact interpolator for deterministic computer experiments, σ²(x) = 0 at all training sample points x_t. The variance increases as query points move further into unexplored regions of the design space. - Confidence bounds: The local standard error σ(x) = √(σ²(x)) allows researchers to construct 95% confidence intervals: ŷ(x) ± 1.96 × σ(x). These intervals directly feed active learning, Bayesian optimization acquisition functions (e.g., Expected Improvement), and adaptive sampling workflows.
When running multi-dimensional sweeps or archiving large validation matrices, store sample points and surrogate prediction arrays using standard scientific containers like h5py for structured array storage.
Troubleshooting
| Issue | Cause | Resolution |
|---|---|---|
| Ill-conditioned covariance matrix / Cholesky failure | Two sample points in xt are positioned too close together in the design space, or correlation parameter values θ become excessively large, producing a near-singular covariance matrix. |
Ensure design samples generated with LHS enforce a minimum Euclidean distance threshold. If training on legacy simulation datasets, add a small regularization nugget (e.g., theta0=[1e-2]) or filter duplicate configurations prior to calling set_training_values. |
| High-dimensional training slowdown (d ≥ 15) | Standard Kriging requires optimizing d hyperparameter dimensions and inverting an n × n covariance matrix during likelihood evaluation, scaling poorly when design spaces exceed 15–20 parameters. | Switch from KRG to Kriging with Partial Least Squares (from smt.surrogate_models import KPLS). KPLS projects the high-dimensional variable space onto a lower-dimensional principal latent subspace while maintaining variance information. |
| Spurious oscillations near non-smooth responses | Using the squared exponential kernel (corr='squar_exp') imposes infinite differentiability assumptions, causing Gibbs-like overshoot near step changes or boundary layer transitions. |
Change the correlation kernel to Matérn 5/2 (corr='matern52') or Matérn 3/2 (corr='matern32'). Matérn kernels relax smoothness constraints and handle non-analytic physical gradients robustly without numerical divergence. |
| Shape mismatch errors during set_training_values | Passing 1D arrays of shape (n,) causes SMT internal matrix operations to fail. |
Ensure inputs are strictly 2D arrays: xt must have shape (n_samples, n_features) and yt must have shape (n_samples, n_outputs). Use x.reshape(-1, 1) for single-variable problems. |
Free alternatives (table)
Researchers select surrogate modeling packages based on their mathematical algorithms, programming language ecosystem, and high-performance computing (HPC) requirements.
| Software | License | Primary Algorithms | Key Strength | Research Use Case |
|---|---|---|---|---|
| SMT (Surrogate Modeling Toolbox) | BSD 3-Clause | Kriging, KPLS, GEKPLS, RBF, Polynomial Chaos (QPCE) | Gradient-enhanced models, partial least squares, and mixed-variable support | Multidisciplinary design optimization (MDO), aerospace aerodynamics, and engineering simulations |
| Dakota | GNU LGPL | Gaussian Process Models, Polynomial Chaos Expansions (PCE), MARS | Enterprise HPC clustering, iterative coupling, and industrial-scale UQ | National laboratory engineering physics, reactor safety, and complex multi-physics analysis |
| OpenTURNS | GNU LGPL | Kriging, Polynomial Chaos, Copula probabilistic modeling | Rigorous mathematical reliability analysis and sensitivity index calculation | Civil engineering structures, nuclear safety assessments, and environmental risk modeling |
| BoTorch / GPyTorch | BSD 3-Clause | Gaussian Processes, Sparse GPs, Deep GPs | PyTorch GPU acceleration, modular Bayesian optimization acquisition functions | High-throughput experimental chemistry, automated drug discovery, and machine learning hyperparameter tuning |
| Surrogates.jl | MIT | Kriging, Radial Basis Functions, Neural Network Surrogates | Native Julia compilation speed and seamless integration with SciML differential equations | Scientific machine learning, dynamic physical system emulation, and real-time control |
| scikit-learn (GaussianProcessRegressor) | BSD 3-Clause | Gaussian Process Regression with customizable composite kernels | Universal adoption, simple API, and native interoperability with Python data science stacks | Low-dimensional baseline metamodeling, exploratory data analysis, and teaching laboratories |
What is the difference between a surrogate model and a reduced-order model (ROM)?
A surrogate model (black-box metamodel) treats the simulation purely as an input-output mapping, fitting statistical or machine learning functions to sampled coordinates without direct reference to the governing equations during evaluation. In contrast, a reduced-order model (such as Proper Orthogonal Decomposition or Galerkin projection) projects the actual governing physical equations and discretization meshes onto a reduced-basis subspace, preserving the underlying conservation laws.
How many training samples are needed to train a reliable surrogate model?
A standard empirical heuristic for Gaussian process and Kriging surrogates is 10 to 20 training samples per input dimension (10d to 20d) for moderately non-linear behaviors. Highly non-linear physical phenomena or high-dimensional spaces require adaptive sampling, where initial Latin Hypercube samples are augmented iteratively at locations exhibiting high prediction variance.
When should researchers choose Kriging over deep neural network surrogates?
Kriging is preferred when the computational budget allows only small-to-moderate sample sizes (N ≤ 2,000), because it provides exact interpolation and rigorous analytical variance bounds. Deep neural networks (or Fourier Neural Operators) require significantly larger training datasets (N ≥ 10,000) to generalize effectively, but excel when modeling full 2D/3D spatial field distributions rather than scalar engineering metrics.
Why does prediction variance drop to zero at training sample points?
Standard Kriging models deterministic computer experiments where evaluating the same input parameter set will always yield the exact same simulation result. Because there is no measurement noise, the model interpolates training points exactly, and epistemic uncertainty at those coordinates is zero.
How does Latin Hypercube Sampling compare to random sampling?
Latin Hypercube Sampling divides each parameter dimension into N equally probable intervals and places exactly one sample in each interval, ensuring uniform marginal coverage across all variables. Pseudo-random sampling frequently creates spatial clustering and leaves large unobserved voids in multi-dimensional parameter spaces.
Sources (official docs)
- SMT Documentation (ReadTheDocs): https://smt.readthedocs.io/en/latest/ (checked 2026-11-08)
- Bouhlel, M. A., Hwang, J. T., Bartoli, N., Lafage, R., Morlier, J., and Martins, J. R. R. A. (2019). “A Python surrogate modeling framework with derivatives”. Advances in Engineering Software, 135, 102662. DOI: 10.1016/j.advengsoft.2019.03.005 (checked 2026-11-08)
- Saves, P., Lafage, R., Bartoli, N., Diouane, Y., Bussemaker, J., Lefebvre, T., Hwang, J. T., Morlier, J., and Martins, J. R. R. A. (2024). “SMT 2.0: A Python package for surrogate modeling with mixed variables and hierarchical spaces”. Advances in Engineering Software, 188, 103571. DOI: 10.1016/j.advengsoft.2023.103571 (checked 2026-11-08)
- Sandia National Laboratories Dakota Project: https://dakota.sandia.gov/ (checked 2026-11-08)
- OpenTURNS Uncertainty Quantification Platform: https://openturns.github.io/ (checked 2026-11-08)
Research-administration guidance, standards updates and independent tool reviews.