Put it to work.

A space-filling design (sampling) spreads runs evenly before anything is known. Once a first surrogate exists, it can say where the next run would teach it most. Choosing runs one batch at a time this way is sequential design, or active learning.

The loop

1. Run an initial space-filling design, a few runs per input dimension.
2. Fit a surrogate that knows its uncertainty: a Gaussian process, or an ensemble.
3. Score a large set of candidate designs with an acquisition function, and pick the best.
4. Run the solver there, add the result, refit. Stop when the budget is spent or the score stops improving.

What to ask for

• Uncertainty sampling: pick the candidate with the largest σ. It improves the model everywhere, which is what you want for a general-purpose surrogate.
• Expected improvement: pick the candidate where the chance of beating the best design so far, weighted by how much, is highest. It balances exploiting good regions against exploring uncertain ones, and it is the usual choice for Bayesian optimisation.

import numpy as np
from scipy.stats import qmc
from sklearn.gaussian_process import GaussianProcessRegressor
from sklearn.gaussian_process.kernels import RBF, ConstantKernel

def solver(X):
    # Stands in for one expensive CFD run per row.
    a, s = X[:, 0], X[:, 1]
    return np.sin(3 * a) * np.cos(2 * s) + 0.5 * a

rng = np.random.default_rng(0)
lhs = qmc.LatinHypercube(d=2, rng=rng)
X = lhs.random(8)               # 8 initial runs
y = solver(X)
X_cand = lhs.random(2000)       # designs we could run next
X_test = rng.random((500, 2))
kernel = ConstantKernel() * RBF([0.3, 0.3])

for step in range(12):          # 12 more runs, one by one
    gp = GaussianProcessRegressor(kernel, normalize_y=True)
    gp.fit(X, y)
    mean, std = gp.predict(X_cand, return_std=True)
    # The candidate where the model is least sure.
    x_next = X_cand[[std.argmax()]]
    X = np.vstack([X, x_next])
    y = np.append(y, solver(x_next))
    err = gp.predict(X_test) - solver(X_test)
    rmse = np.sqrt(np.mean(err**2))
    print(step, len(y), round(rmse, 3), round(std.max(), 3))

In this run the test RMSE falls from 0.046 to 0.001 as 12 runs are added, and the largest σ falls with it: σ is a usable stopping signal. For expected improvement, batches and constraints, use BoTorch ↗ rather than writing your own.

In practice

• Batches. On a cluster you run 10 or 50 simulations at once. Picking the top 10 σ puts them side by side; batch methods add one point, pretend its result is known, update σ, and pick the next.
• Neural surrogates. Use the spread of a 5-member ensemble as σ. Retraining five networks per round is the main cost; fine-tune from the previous round instead of starting over.
• When it does not help. If the response is smooth and the budget is large, a good space-filling design is nearly as good and much simpler. Active learning pays when runs are expensive and the interesting behaviour is local.

Further on this site: Active operator learning chooses simulations for neural operators with cheap uncertainty estimates; PhysicsNeMo added active-learning workflows in its 26.05 release.

A neural surrogate is differentiable. Backpropagation (backpropagation and autograd) gives the gradient of the loss with respect to the weights; the same machinery gives the gradient of the prediction with respect to the inputs. If the inputs are design parameters, that gradient says how to change the design to improve it, at the cost of one backward pass. Solvers can give such gradients too, through adjoint methods, but each costs about as much as another simulation.

Optimise the inputs, not the weights

The recipe turns the training loop inside out. Freeze the network. Make the design a tensor with requires_grad=True and hand it to an optimiser. The loss is the objective to minimise, for example −CL/CD to maximise lift-to-drag. After each step, clamp the design back into the training box.

import torch
from torch import nn

torch.manual_seed(0)

def solver(d):
    # Stands in for one CFD run: returns (CL, CD).
    a, s = d[..., 0], d[..., 1]
    cl = 0.09 * (a + 1.5) * torch.cos(torch.deg2rad(s))
    cd = 0.008 + 0.002 * (s / 30) ** 2 + 0.045 * cl**2
    return torch.stack([cl, cd], -1)

# 1. A small surrogate trained on 45 runs inside the box.
lo, hi = torch.tensor([-2.0, 20.0]), torch.tensor([8.0, 36.0])
X = lo + (hi - lo) * torch.rand(45, 2)
scale = torch.tensor([1.0, 20.0])  # CL, CD similar sizes
net = nn.Sequential(nn.Linear(2, 64), nn.SiLU(),
                    nn.Linear(64, 64), nn.SiLU(),
                    nn.Linear(64, 2))
def surrogate(d):
    return net((d - lo) / (hi - lo)) / scale
opt = torch.optim.Adam(net.parameters(), lr=3e-3)
for _ in range(3000):
    opt.zero_grad()
    loss = nn.functional.mse_loss(surrogate(X) * scale,
                                  solver(X) * scale)
    loss.backward()
    opt.step()
net.requires_grad_(False)          # freeze the surrogate

# 2. Optimise the design: alpha and sweep are now
# the parameters being trained.
design = torch.tensor([2.0, 30.0], requires_grad=True)
opt = torch.optim.Adam([design], lr=0.05)
for _ in range(500):
    cl, cd = surrogate(design)
    loss = -cl / cd                # maximise lift-to-drag
    opt.zero_grad()
    loss.backward()
    opt.step()
    with torch.no_grad():
        design.clamp_(lo, hi)      # stay inside the box

# 3. Re-simulate the winner before believing it.
with torch.no_grad():
    cl, cd = surrogate(design)
    cl_s, cd_s = solver(design)
print(design.detach().numpy().round(2))
print(round(float(cl / cd), 2), round(float(cl_s / cd_s), 2))

A typical run ends near α ≈ 3.7°, Λ = 20°, with a predicted L/D of about 25 and the solver agreeing to within about 0.1. Note where the optimum landed: on the edge of the box, at the smallest sweep allowed. That is common, and it is exactly where the surrogate has the fewest neighbouring runs.

Optimisers find the surrogate's mistakes

An optimiser searches for the input with the best predicted value. If the surrogate is wrong somewhere in a favourable direction, the optimiser will find that place, because it looks better than anything real. The more freedom the design has, the more such places exist. Guards, from simplest to strongest:

• Bounds and a trust region. Keep the design inside the training box, and limit each optimisation round to a small region around known runs.
• Penalise uncertainty. Optimise the predicted value plus k times the ensemble spread (for a minimisation), so uncertain regions look worse rather than better.
• Robust objectives. Optimise the worst case inside a conformal interval (uncertainty), so the chosen design is good even if the surrogate is off by its calibrated error.
• Close the loop. Simulate the optimum, add it to the data, retrain, optimise again: active learning aimed at the optimum.

Gradient-free optimisers (differential evolution, genetic algorithms, Bayesian optimisation) work with any surrogate, including non-differentiable ones, and escape local optima more easily; with a surrogate their many evaluations are affordable. Gradients win when there are tens or hundreds of design variables.

Further on this site: Robust predict-then-optimize builds conformal bounds into the optimisation; CANTO optimises CAD parameters through a surrogate and verifies every design with CFD, which is the subject of learning straight from CAD.