---
kind: concept
slug: surrogate-modeling
title: Surrogate modeling
summary: >-
  A learned stand-in for an expensive simulator or physical system, so design sweeps and gradient estimates run cheaply. Variants differ in what they predict: one objective, a field, per-node stress, a material law with built-in thermodynamics, a gradient-only digital model, a time series of modal amplitudes, or a score for finished runs that saves no simulator calls.
tags: [surrogate, simulation, optimization, cfd, physics-informed, reduced-order-model]
papers: ["2609.16066", "2507.05291", "2104.13386", "2412.13802", "2005.12183", "2609.17101"]
related: [neural-operator, fourier-neural-operator, graph-neural-network, physics-informed-neural-network, finite-element-method, hardware-aware-training, metaheuristic-optimization, self-attention, surrogate-gradient, fuzz-testing, automatic-differentiation, dynamic-mode-decomposition, digital-twin, dimensionality-reduction]
status: draft
updated_at: 2026-11-01
---

## Overview

A surrogate model learns a map from inputs to a system's output, so that each
new query costs one model evaluation and not a simulation run or a device
measurement [arxiv:2609.16066 §1] [arxiv:2507.05291 §1]. Its accuracy is
bounded by the reference data it was trained on. The first two variants below
train on solver output and keep the solver in the loop for checking final
designs [arxiv:2609.16066 §2.3] [arxiv:2609.16066 §4.3]
[arxiv:2507.05291 §4.1]. The fourth is different: a learned scorer for runs
that have already been simulated, which saves no simulator calls
[arxiv:2412.13802 §5.1.2].

The main design choice is what the surrogate predicts. A classical regressor
predicts one objective directly [arxiv:2609.16066 §4.1.3]. A field surrogate
predicts a whole field and computes the objectives afterwards
[arxiv:2609.16066 §1] [arxiv:2609.16066 §3.2]. A graph surrogate predicts a
stress value at every node of a mesh [arxiv:2507.05291 §3]. A material
surrogate predicts a stored-energy function, and stress follows from its
derivatives [arxiv:2005.12183 §3.3]. A further use replaces the derivatives of
a physical device rather than its output [arxiv:2104.13386 app. S1]. That use
is not the same idea as a [surrogate gradient](/concepts/surrogate-gradient),
which replaces a non-differentiable step inside a network during training of
spiking networks. Another use scores a finished run, so that a search can rank
it [arxiv:2412.13802 §4.2.1]. A last use is over time rather than over
design: a few spatial modes from a transient solver are kept, and a small
network predicts how their amplitudes evolve [arxiv:2609.17101 §4]
[arxiv:2609.17101 §5].

## Variants

### Field surrogate with a branch-trunk neural operator

The surrogate maps the 28 blade-geometry variables and a query point (r, z) on
a meridional plane to five basic flow fields: static pressure, static
temperature, absolute and relative squared velocity, and mass flux
[arxiv:2609.16066 §2.1] [arxiv:2609.16066 §3.3]. Performance metrics such as
efficiency, pressure ratio and mass flow are then computed from the predicted
fields by averaging, not learned as separate outputs [arxiv:2609.16066 §3.2]
[arxiv:2609.16066 tab. 6]. The model is a branch-trunk neural operator
(TNO) whose trunk uses Galerkin-type attention blocks instead of softmax
attention [arxiv:2609.16066 §3.1]. Card: [Panoramic aerodynamic performance
prediction for turbomachinery cascades](/papers/2609.16066).

Reported results, all author-reported on one compressor (Rotor 37) with a
validation split only: mean field error 0.083% against 0.572% for an FNO
baseline [arxiv:2609.16066 tab. 4]; one forward pass in about 0.9 ms, against
roughly 800 s per CFD sample on the authors' CPU [arxiv:2609.16066 tab. 3]
[arxiv:2609.16066 §2.2]. Training the model took 40.5 minutes on the authors'
GPU [arxiv:2609.16066 tab. 3].

### Graph surrogate for local stress on a periodic micro-structure

The surrogate maps a periodic micro-structure mesh and its mean stress to the
stress at each mesh node, so local stress is available without a full
micro-scale solve for each macro point [arxiv:2507.05291 §1]
[arxiv:2507.05291 §3]. The mesh becomes a graph whose nodes carry the mean
stress, the position and a boundary label, and whose edges carry the length
between neighbours. Extra edges join opposite sides of the cell, which is how
periodicity enters [arxiv:2507.05291 §3.1]. The network is a message-passing
graph neural network with 10 steps and a latent size of 128
[arxiv:2507.05291 §3.2]. Training adds a penalty on the discrete stress
divergence, weighted by λ, so that predictions approximately satisfy local
equilibrium [arxiv:2507.05291 §3.3]. The linear dataset has 10,000 periodic
meshes and the hyperelastic one 9,795 [arxiv:2507.05291 §4.1]
[arxiv:2507.05291 §5.1]. Card: [Physics-Informed Graph Neural Networks to
Reconstruct Local Fields Considering Finite Strain Hyperelasticity](/papers/2507.05291).

Reported results, author-reported on single values with no spread shown: on
linear elastic meshes, NMSE is 1.08e-2 with the divergence penalty, against
1.13e-2 for a periodic graph network without it and 2.33e-2 without periodic
edges [arxiv:2507.05291 tab. 1]. On hyperelastic meshes, NMSE is 2.18e-2
against 5.36e-2 for the periodic graph network [arxiv:2507.05291 tab. 2]. Its
inference is said to be about 500 times faster than finite-element solves on
a CPU, but that compares GPU inference with CPU solves and excludes training
[arxiv:2507.05291 §5.3.1].

### Differentiable digital model for physical systems

Here the surrogate is a neural network fitted to input-output pairs measured
on a physical device. The device runs the forward pass, and the surrogate
supplies the derivatives in the backward pass. The loss is computed from the
measured outputs, so the surrogate's errors affect the gradient direction and
not the loss value [arxiv:2104.13386 app. S1]. This avoids estimating the
device's own derivatives by finite differences, which would cost one physical
call per input and parameter [arxiv:2104.13386 app. S1]. The digital model in
typical runs is trained on about 10^4 uniformly sampled inputs
[arxiv:2104.13386 app. S2]. Card: [Deep physical neural networks enabled by a
backpropagation algorithm for arbitrary physical systems](/papers/2104.13386).

Reported results, author-reported on three physical systems. A five-stage
optical chain reaches 93% on vowel classification, against about 40% at best
for training entirely in simulation [arxiv:2104.13386 fig. 2]
[arxiv:2104.13386 fig. 3]. On MNIST, a mechanical plate reaches 87%, an
electronic circuit 93% and a hybrid optical network 97%
[arxiv:2104.13386 fig. 4]. The paper's in-simulation baseline includes a
fitted noise model [arxiv:2104.13386 app. S2].

### Learned violation scorer for scenario search

Here the model does not stand in for the simulator. A Transformer encoder reads
the per-frame states of the vehicles in a finished simulation, seven features
per vehicle, and outputs one scalar that the paper calls a violation
probability [arxiv:2412.13802 §4.1] [arxiv:2412.13802 §4.2.1]
[arxiv:2412.13802 tab. 2]. That scalar is one of four objectives in a
multi-objective genetic search over driving scenarios. The others are a
road-feature classifier, the minimum distance to another vehicle, and the
number of unique violations [arxiv:2412.13802 §4.2.1]. The scorer is trained
on 1,000 short scenarios of two minutes each, with binary cross-entropy and
Adam [arxiv:2412.13802 §5.1.2]. Card: [SimADFuzz: Simulation-Feedback Fuzz
Testing for Autonomous Driving Systems](/papers/2412.13802).

Reported results, author-reported on one driving system (InterFuser) in CARLA
on one map: 35 unique violations in six hours, against 26, 18 and 3 for the
three baselines [arxiv:2412.13802 §5.2.2]. The paper lists the single system
and map as a threat to external validity [arxiv:2412.13802 §5.3]. The scorer
does not reduce cost. Every candidate is still simulated, and the search is
bounded by time budgets rather than by a count of simulator calls
[arxiv:2412.13802 §5.1.2] [arxiv:2412.13802 §5.2].

### Thermodynamics-constrained network for a material law

Here the surrogate replaces a hand-made constitutive model, the rule that gives
stress from strain at a material point. Such models need calibration, and
micromechanical simulations are too costly to run inside large finite-element
analyses [arxiv:2005.12183 §1]. The network predicts a Helmholtz free-energy
function and the increments of internal variables. Stress and dissipation are
not network outputs. They are derivatives of the energy network, so the first
law holds by construction [arxiv:2005.12183 §3.3] [eq. 13]. The training data
are generated from a reference thermodynamic model, not from measurements
[arxiv:2005.12183 §4.2]. The energy network uses a squared-argument ELU, whose
second derivative is nonzero everywhere. Activations whose second derivative
vanishes over part of their range fit the stress output poorly
[arxiv:2005.12183 §3.2] [tab. 2]. The scope is strain-rate-independent,
small-strain processes, with isothermal experiments [arxiv:2005.12183 §2.3]
[arxiv:2005.12183 §2.4]. Card: [Thermodynamics-based Artificial Neural Networks
for constitutive modeling](/papers/2005.12183).

Reported results, author-reported on synthetic data from 1D and 3D
elasto-plastic reference models. In 1D, a standard network with a matched
number of hyper-parameters gives the wrong energy and dissipation, and
generalizes worse on cyclic and random paths [arxiv:2005.12183 §5.1.4]. Beyond
the training range of strain increments, accuracy drops, but the predicted
dissipation stays positive [arxiv:2005.12183 §5.1.3]. In 3D, the network agrees
with the reference for uniaxial cyclic loading, with less accuracy at the
largest increments. The authors report larger errors for it in the mean and
deviatoric stress [arxiv:2005.12183 §5.2.2] [arxiv:2005.12183 §5.2.3]. The
paper states that code is available on request and gives no repository
[arxiv:2005.12183 §1].

### Reduced-order DMD model with a neural amplitude model

Here the surrogate is a time-stepping model of a transient solver. The
training snapshots come from a finite-difference solver for the one-dimensional
viscous Burgers equation with a discontinuous initial state, at Reynolds numbers
10², 10³ and 10⁴ [arxiv:2609.17101 §2]. Each case gets its own model, and
no model is fitted across Reynolds numbers [arxiv:2609.17101 tab. 1].
Randomized dynamic mode decomposition gives a set of spatial modes and their
amplitudes. The number of modes is chosen by a search over all ranks, which
keeps the rank with the smallest time-averaged error [arxiv:2609.17101 §4]
[arxiv:2609.17101 alg. 1] [arxiv:2609.17101 eq. 29]. The randomized SVD inside
it draws a random test matrix [arxiv:2609.17101 alg. 2]. A cascade-forward
network with 10 hidden units then models the time series of each amplitude,
using past estimator outputs and delayed amplitude inputs [arxiv:2609.17101 §5]
[arxiv:2609.17101 fig. 3] [arxiv:2609.17101 eq. 32]. The modes and the
estimated coefficients together form the data model. The paper does not
require a Galerkin projection of the governing equations [arxiv:2609.17101
abstract]. Card: [Randomized-DMD digital twin data models](/papers/2609.17101).

Reported results, author-reported and in-sample only. The DMD error is 3.2e-7
at Re 10² with 15 modes, 4.4e-7 at Re 10³ with 20 modes and 5.4e-8 at Re 10⁴
with 20 modes [arxiv:2609.17101 tab. 1]. The data-model error, grouped by
coefficient, is between 1.3e-4 and 8.1e-4 across the three cases
[arxiv:2609.17101 tab. 2] [arxiv:2609.17101 tab. 3] [arxiv:2609.17101 tab. 4].
The paper names no external baseline and reports no run time for the
full-order solver. Offline work takes at most about 2 s and the online stage
8 to 17 s, on hardware that is not stated [arxiv:2609.17101 §6]
[arxiv:2609.17101 fig. 13]. A speed-up over the solver is therefore not
shown.

## When to use which

The papers compare against their own baselines and, for the first variant,
against classical regressors on one task. The rows below say what the reported
evidence supports, and nothing more.

| Situation | Use | Why |
|---|---|---|
| You need one scalar objective and have few simulator runs (about 500) | Try a classical regressor first, then compare | At 500 samples SVR has the lowest pressure-ratio error, while TNO has the lowest efficiency error (0.823‰ against 1.040‰ for SVR) [arxiv:2609.16066 tab. 8] |
| You need many metrics from one run, such as efficiency, pressure ratio and mass flow | A field surrogate, then derive metrics | One predicted field feeds several formulas, so sensitivity and optimisation reuse the same model [arxiv:2609.16066 §3.2] [arxiv:2609.16066 §4.2] |
| You have about 1500 or more samples on a fixed design parameterisation | A field surrogate is a candidate | At 1500 samples and above, TNO has the lowest error of the baselines on efficiency, pressure ratio and mass flow [arxiv:2609.16066 tab. 8] |
| An optimisation loop needs thousands of evaluations | A surrogate with a GA or NSGA-II search | The paper runs GA and NSGA-II on the surrogate and reports 121 s for the optimisation [arxiv:2609.16066 §4.3] |
| You will act on the optimum | Verify the optimum with the simulator | Optimal designs were re-simulated with CFD; the largest pressure-loss error in the multi-objective table was 14.7‰ [arxiv:2609.16066 tab. 10] |
| You need local stress at each node of a periodic micro-structure, because the mean stress alone is not enough (for example for fracture or fatigue) | A graph surrogate, if you can afford about 10,000 solves for training and the cell is periodic | Reduced models keep only the average field, which loses local stress [arxiv:2507.05291 §1]. Its linear and hyperelastic datasets have about 10,000 meshes each [arxiv:2507.05291 §4.1] [arxiv:2507.05291 §5.1] |
| You need a fast stand-in for a micro-scale solve, and the speed-up figure matters | Use the graph surrogate's timing with care | The 500x figure compares GPU inference with CPU solves and excludes training [arxiv:2507.05291 §5.3.1] |
| The device can run the forward pass, but you cannot differentiate it | Physics-aware training with a digital model | Gradients come from the model and the loss from measured outputs [arxiv:2104.13386 app. S1]. The model needs about 10^4 measured samples [arxiv:2104.13386 app. S2], and at 200 Hz one million device runs take about 5,000 s [arxiv:2104.13386 app. S2] |
| You need a score for a finished run, to rank candidates in a search | A learned scorer on the run's own trajectory | The scorer reads the vehicle states of each run and enters the search as one objective [arxiv:2412.13802 §4.2.1] |
| You need a rate-independent, small-strain material law, and a reference model can generate its training data | A thermodynamics-constrained network that predicts an energy function and derives stress | The first law holds by construction, and the authors report better accuracy and generalization than a standard network in 1D [arxiv:2005.12183 §5.1.4]. In 3D the reported stress errors are larger in places [arxiv:2005.12183 §5.2.3], and the paper's data are synthetic [arxiv:2005.12183 §4.2] |
| You need fewer simulator runs | Not this variant; use a simulator-bound surrogate above | Every candidate is still simulated before the scorer runs [arxiv:2412.13802 §5.1.2] |
| You need a cheap time-stepping stand-in for a transient solver, at one fixed setting, and in-sample error is enough for now | A reduced-order DMD model with a neural amplitude model, built per setting | The paper reports low in-sample errors on three Burgers cases, but tests no held-out time window, initial condition or Reynolds number [arxiv:2609.17101 §6] [arxiv:2609.17101 tab. 1] |
| You need a model that covers several settings, such as Reynolds numbers | Not shown by this variant; the paper trains one model per case | Each case has its own mode count and estimator orders, and no model is tested across cases [arxiv:2609.17101 tab. 1] [arxiv:2609.17101 tab. 2] |

## Implementation sketch

The loop below covers the two simulator-bound variants: sample the design
space, run the simulator, fit the model, search on the model, and verify the
result with the simulator [arxiv:2609.16066 §2.3] [arxiv:2609.16066 §4.3]. The
graph variant fits the same loop with a graph input and a per-node target
[arxiv:2507.05291 §3]. The details of the architecture are not stated in full
for the field variant, so this is a skeleton, not a reimplementation. The
physics-aware variant does not follow this loop. It trains on device
measurements and uses the model only for gradients [arxiv:2104.13386 app. S1].
The learned-scorer variant adds one scoring call after each simulation and
removes no simulation [arxiv:2412.13802 §5.1.2]. The material-law variant
follows neither loop. It trains on data from a reference material model, and
its stress comes from differentiating the energy network [arxiv:2005.12183
§3.3] [arxiv:2005.12183 §4.2]. The reduced-order variant follows neither
loop. It decomposes solver snapshots in time and fits a network to the
amplitude series, with no design search and no per-candidate simulation
[arxiv:2609.17101 §4] [arxiv:2609.17101 §5].

```pseudocode
# Surrogate loop for a simulator-bound design task  [§2.3] [§3.3] [§4.3]
# Input:  design-variable box; CFD simulator
# Output: surrogate model and a CFD-verified design

def surrogate_loop(design_box, simulator):
    # Sample designs spread across the box (Latin hypercube in the paper)
    designs = latin_hypercube(design_box, n=2900)     # [§2.3]
    # Run the simulator once per design (about 800 s each in the paper)
    fields = [simulator(d) for d in designs]          # [§2.2]
    # Reduce each simulation to the training target
    # (circumferential mass-flux average in the paper)  [§3.2] [§3.3]
    targets = [average_fields(f) for f in fields]
    # Fit the surrogate on a train split, check on a validation split
    model = fit(designs[:2500], targets[:2500])       # [§3.3]
    # Search on the surrogate, not the simulator (GA or NSGA-II)  [§4.3]
    candidate = search(lambda d: metrics(model(d)))
    # Verify the chosen design with the simulator before using it  [§4.3]
    return candidate, simulator(candidate)
```

## Contradictions

No paper disagrees with another in the reported numbers, but the metrics are
not comparable. The field error in percent [arxiv:2609.16066 tab. 4], the NMSE
of the graph variant [arxiv:2507.05291 tab. 1], the classification accuracy
of the physical variant [arxiv:2104.13386 fig. 2] and the unique-violation
count of the scorer variant [arxiv:2412.13802 §5.2.2] and the DMD and
data-model errors of the reduced-order variant [arxiv:2609.17101 eq. 27]
[arxiv:2609.17101 tab. 1] measure different things on different tasks. Read
them within each paper only.

Inside the cited papers, several points disagree or do not match:

- Inside the field-surrogate paper, the text and the printed equation describe
  the branch-trunk combination differently, so check that step before
  implementing [arxiv:2609.16066 §3.1] [arxiv:2609.16066 eq. 1].
- Inside the graph paper, the linear divergence value is 4.56e-4 in the text
  and 3.23e-4 in the table, with no explanation given
  [arxiv:2507.05291 §4.3] [arxiv:2507.05291 tab. 1].
- Inside the physical-variant paper, the in-simulation baseline is described
  as "barely distinguishable from random guessing" in the figure caption, and
  as reaching about 40% at its best in the text [arxiv:2104.13386 fig. 3]
  [arxiv:2104.13386 fig. 2].
- Inside the scorer paper, the six-hour margin over DriveFuzz is 35 − 18 = 17
  from the table, but the text and conclusion say 27 more
  [arxiv:2412.13802 §5.2.2] [arxiv:2412.13802 §7].
- Inside the reduced-order paper, the snapshot count is 301 in the text
  [arxiv:2609.17101 §2] but 300 in the caption of its Table 1
  [arxiv:2609.17101 tab. 1]. The text claims the method has the best
  correlation, yet the table lists 1.0000 for every case, so the metric does
  not separate them [arxiv:2609.17101 §6] [arxiv:2609.17101 tab. 1].
- Inside the scorer paper, the output is described as a violation probability,
  but the text does not say whether it is squashed before it is used as a
  fitness term [arxiv:2412.13802 §4.2.1].

## Papers

- [Panoramic aerodynamic performance prediction for turbomachinery cascades](/papers/2609.16066): a transformer-style neural operator predicts compressor flow fields from blade geometry, and metrics are computed from the fields.
- [Physics-Informed Graph Neural Networks to Reconstruct Local Fields Considering Finite Strain Hyperelasticity](/papers/2507.05291): a graph network predicts the stress at every node of a periodic micro-structure, with a divergence penalty for equilibrium.
- [Deep physical neural networks enabled by a backpropagation algorithm for arbitrary physical systems](/papers/2104.13386): a digital model fitted to measured device data supplies the gradients for training physical networks.
- [SimADFuzz: Simulation-Feedback Fuzz Testing for Autonomous Driving Systems](/papers/2412.13802): a genetic fuzzer whose Transformer scorer of simulated runs is one of four selection objectives.
- [Thermodynamics-based Artificial Neural Networks for constitutive modeling](/papers/2005.12183): a network predicts a free-energy function, and stress and dissipation follow by differentiation, so the first law holds by construction.
- [Randomized-DMD digital twin data models](/papers/2609.17101): a randomized DMD gives spatial modes, and a cascade-forward network predicts their amplitudes over time, for one Burgers setting at a time.