Benchmark Problems And Setup

This page describes the empirical benchmark behind the claims in Benchmark Results And Discussion: the seven engineering problems and their datasets, the surrogate configurations, the sampling and scoring protocol, and the acceptance criteria. Everything here can be reproduced with the scripts in BM/.

Note

The benchmark tests the package’s central premise: a few hundred samples of an expensive, deterministic engineering model are enough to train a surrogate that reproduces it accurately on unseen inputs, and a single interface makes it cheap to find which surrogate family suits a problem.

Datasets at a glance

All datasets are generated by BM/datasets/generate_datasets.py from deterministic physics models. They are therefore noise-free: any test error is approximation error of the surrogate, not measurement noise. Inputs are sampled with a Latin hypercube design (seed 2026). Inputs marked with * in the tables below are sampled log-uniformly. Each row is assigned once to a fixed train or test split (Bernoulli, \(P(\text{test}) = 0.2\)), so all models are always scored on the same held-out points.

Dataset

Discipline

Physics model

\(d\)

Outputs

Train / test rows

aero_naca4_airfoil_panel

Aerodynamics

Linear-vortex panel method, Prandtl–Glauert, viscous drag estimate

6

5

1229 / 271

aero_finite_wing_liftingline

Aerodynamics

Prandtl lifting line, 24-term Glauert series

6

5

1630 / 370

struct_10bar_truss_fem

Structures

Linear truss FEM, direct stiffness

14

15

2427 / 573

struct_thick_cylinder_lame

Structures

Lamé thick-walled cylinder, closed ends

7

8

2427 / 573

heat_pin_fin

Heat transfer

1-D fin with convective tip

6

5

2427 / 573

heat_2d_plate_conduction_fdm

Heat transfer

2-D steady conduction, 41 × 41 finite differences

7

7

665 / 135

heat_counterflow_hx_entu

Heat transfer

Counter-flow heat exchanger, ε-NTU

7

6

2427 / 573

Together the seven problems have 51 outputs. They cover smooth closed-form responses, discretized PDE solutions, a high-dimensional finite-element model and some deliberately non-smooth targets.

Problem descriptions

Airfoil aerodynamics (NACA 4-digit panel method)

Inputs: maximum camber \(m \in [0, 0.06]\), camber position \(p \in [0.2, 0.6]\), thickness \(t \in [0.06, 0.20]\), angle of attack \(\alpha \in [-4^\circ, 8^\circ]\), Mach number \(M \in [0.05, 0.5]\), Reynolds number \(Re^* \in [2 \times 10^5, 10^7]\). Outputs: \(C_L\), \(C_D\), \(C_{M,c/4}\), \(L/D\), \(C_{p,\min}\).

The geometry is discretized into 160 linear-strength vortex panels with a Kutta condition. The incompressible coefficients are corrected with the Prandtl–Glauert rule,

\[ C_L = \frac{C_{L,0}}{\sqrt{1 - M^2}}, \qquad C_M = \frac{C_{M,0}}{\sqrt{1 - M^2}}, \]

and drag is an engineering estimate built from a turbulent flat-plate skin friction and a thickness form factor:

\[ C_D = 2\,c_f\,\mathrm{FF}\,(1 + 0.06\,C_L^2), \qquad c_f = \frac{0.455}{(\log_{10} Re)^{2.58}}, \qquad \mathrm{FF} = 1 + 2t + 60t^4. \]

Only subcritical samples (\(C_{p,\min} > C_{p,\mathrm{crit}}(M)\)) are kept. The responses are smooth but strongly coupled, and \(L/D\) is a ratio of two outputs.

Finite wing (Prandtl lifting line)

Inputs: aspect ratio \(AR \in [4, 14]\), taper ratio \(\lambda \in [0.2, 1]\), linear washout \(\in [-6^\circ, 0^\circ]\), root incidence \(\in [-2^\circ, 10^\circ]\), section lift slope \(a_0 \in [5.6, 6.8]\), zero-lift angle \(\alpha_{L0} \in [-4^\circ, 0^\circ]\). Outputs: \(C_L\), \(C_{D_i}\), \(c_{l,\max}\) (section maximum), \(\eta_{c_{l,\max}}\) (spanwise location of that maximum), root bending coefficient.

The circulation is expanded as \(\Gamma(\theta) = 2bV_\infty \sum_{n=1}^{24} A_n \sin(n\theta)\), and the monoplane equation is solved at collocation points by least squares:

\[ \sum_{n=1}^{24} A_n \sin(n\theta)\left(\frac{n\,\mu(\theta)}{\sin\theta} + 1\right) = \mu(\theta)\,\alpha(\theta), \qquad \mu = \frac{c\,a_0}{4b}, \]
\[ C_L = \pi\,AR\,A_1, \qquad C_{D_i} = \frac{C_L^2}{\pi\,AR}(1 + \delta), \qquad \delta = \sum_{n \ge 2} n\left(\frac{A_n}{A_1}\right)^2. \]

\(\eta_{c_{l,\max}}\) is the location of a maximum on a discrete spanwise grid, so it is piecewise constant in the inputs: a deliberately hard, non-smooth target.

10-bar truss (linear finite elements)

Inputs: member areas \(A_1, \dots, A_{10} \in [0.645, 225.8]\ \mathrm{cm}^2\), Young’s modulus \(E \in [60, 210]\) GPa, nodal loads \(P_2, P_4 \in [200, 600]\) kN, density \(\rho \in [2700, 7850]\ \mathrm{kg/m^3}\). Outputs: mass, maximum displacement, vertical displacements of nodes 2 and 4, maximum absolute stress, member stresses \(\sigma_1, \dots, \sigma_{10}\).

The classic cantilever truss (bays of 9.144 m, two pinned nodes) is solved with the direct stiffness method:

\[ K(A, E)\,u = F, \qquad \sigma_e = \frac{E}{L_e}\begin{bmatrix} -c_e & -s_e & c_e & s_e \end{bmatrix} u_e, \qquad m = \rho \sum_{e=1}^{10} A_e L_e. \]

With 14 inputs this is the highest-dimensional problem. Stresses and displacements scale roughly like \(1/A\) in the members that carry the load, so small areas produce heavy, signed tails. The mass, in contrast, is bilinear in \((\rho, A_e)\).

Thick-walled cylinder (Lamé)

Inputs: inner radius \(r_i \in [20, 200]\) mm, wall ratio \(r_o/r_i \in [1.1, 3]\), internal and external pressures \(p_i \in [5, 150]\), \(p_o \in [0, 30]\) MPa, \(E \in [70, 210]\) GPa, Poisson ratio \(\nu \in [0.25, 0.35]\), yield stress \(\sigma_y \in [250, 900]\) MPa. Outputs: hoop stress (inner and outer), axial stress, von Mises and Tresca stress at the bore, radial displacements (inner and outer), von Mises safety factor.

This is the exact closed-form solution for a closed-end cylinder:

\[ A = \frac{p_i r_i^2 - p_o r_o^2}{r_o^2 - r_i^2}, \qquad B = \frac{(p_i - p_o)\,r_i^2 r_o^2}{r_o^2 - r_i^2}, \qquad \sigma_\theta(r) = A + \frac{B}{r^2}, \quad \sigma_r(r) = A - \frac{B}{r^2}, \quad \sigma_a = A, \]
\[ \sigma_{\mathrm{vM}} = \sqrt{\tfrac{1}{2}\bigl[(\sigma_\theta - \sigma_r)^2 + (\sigma_r - \sigma_a)^2 + (\sigma_a - \sigma_\theta)^2\bigr]}, \qquad SF = \frac{\sigma_y}{\sigma_{\mathrm{vM}}}. \]

The responses are rational functions of the inputs. At the bore, \(\sigma_\theta - \sigma_a = \sigma_a - \sigma_r = B/r_i^2\), so

\[ \sigma_{\mathrm{vM}} = \sqrt{3}\,\frac{B}{r_i^2} \;\propto\; (p_i - p_o), \qquad \sigma_{\mathrm{Tresca}} = \frac{2}{\sqrt{3}}\,\sigma_{\mathrm{vM}}. \]

The equivalent stresses therefore vanish, and the safety factor grows without bound, when the two pressure ranges overlap (\(p_i \approx p_o\)). That makes safety_factor_vm heavy-tailed: from 0.26 to 2193 in the test set.

Pin fin (1-D fin theory)

Inputs: diameter \(D \in [2, 20]\) mm, length \(L \in [10, 150]\) mm, conductivity \(k^* \in [15, 400]\) W/(m·K), convection coefficient \(h^* \in [5, 500]\) W/(m²·K), base and ambient temperatures. Outputs: heat rate, tip temperature, fin efficiency, fin effectiveness, Biot number.

With \(m = \sqrt{hP/(kA_c)}\) and \(\theta_b = T_b - T_\infty\), the exact solution for a fin with a convective tip is

\[ q = \sqrt{hPkA_c}\,\theta_b\, \frac{\sinh mL + \frac{h}{mk}\cosh mL}{\cosh mL + \frac{h}{mk}\sinh mL}, \qquad T_{\mathrm{tip}} = T_\infty + \frac{\theta_b}{\cosh mL + \frac{h}{mk}\sinh mL}. \]

The responses are smooth, and several are close to power laws in the log-sampled inputs.

Plate conduction (2-D finite differences)

Inputs: width \(W \in [20, 200]\) mm, \(k^* \in [5, 200]\) W/(m·K), volumetric generation \(q'''^* \in [10^{-3}, 1]\) MW/m³, top convection \(h^* \in [10, 2000]\) W/(m²·K), ambient and side temperatures. Outputs: maximum, mean and centre temperatures, location of the hottest point \((x, y)_{T_{\max}}\), convective heat to the fluid, net heat through the fixed-temperature sides.

The steady heat equation \(\nabla^2 T + q'''/k = 0\) is solved on a 41 × 41 grid with fixed temperatures at the left and right edges, a convective top edge (\(-k\,\partial T/\partial n = h(T - T_\infty)\)) and an insulated bottom edge. As for the wing, the location of the hottest point jumps between grid nodes, so \(x_{T_{\max}}\) and \(y_{T_{\max}}\) are non-smooth by construction.

Counter-flow heat exchanger (ε-NTU)

Inputs: mass flow rates \(\dot m_h^*, \dot m_c^* \in [0.05, 5]\) kg/s, specific heats \(c_{p,h}, c_{p,c} \in [1000, 4200]\) J/(kg·K), conductance \(UA^* \in [100, 20000]\) W/K, inlet temperatures. Outputs: duty \(Q\), outlet temperatures, effectiveness, NTU, LMTD.

With \(C = \dot m c_p\), \(C_r = C_{\min}/C_{\max}\) and \(\mathrm{NTU} = UA/C_{\min}\):

\[ \varepsilon = \frac{1 - e^{-\mathrm{NTU}(1 - C_r)}}{1 - C_r\,e^{-\mathrm{NTU}(1 - C_r)}}, \qquad Q = \varepsilon\,C_{\min}\,(T_{h,\mathrm{in}} - T_{c,\mathrm{in}}). \]

The responses are smooth, but the effectiveness saturates towards 1 at high NTU and the min/max switch between the two streams makes them non-smooth along \(C_h = C_c\).

Surrogate configurations

Each output is fitted by its own single-output surrogate (51 surrogates per model and training sample). The configurations are defined in build_model() in BM/benchmark.py:

Model

Configuration

Hyperparameter tuning

LS

Full quadratic basis (degree 2, with interactions), ridge \(\lambda_2 = 10^{-6}\)

none

RBF

Cubic kernel \(\phi(r) = r^3\) with a linear polynomial tail

none (the cubic kernel has no shape parameter)

Kriging

Gaussian correlation (\(p = 2\)), nugget \(10^{-10}\), anisotropic \(\theta \in [10^{-3}, 10^{2}]^d\)

maximum likelihood with L-BFGS-B, 1 + 2 starts

Ensemble

WeightedEnsemble of the three models above

inverse-MSE weights from 3-fold cross-validation

All four use the package’s default input standardization. The configurations are deliberately off the shelf: nothing is tuned per dataset, so the results show what a user gets from a single, generic setting.

Protocol

Pre-processing

Following the dataset README:

  • Inputs sampled log-uniformly (marked * above) are transformed with \(\log_{10}\).

  • Outputs that are strictly positive and heavy-tailed (training-set skewness \(> 2\)) are fitted as \(\log y\). Of the 51 outputs, 13 qualify: NTU, Q_kW, T_center_C, T_max_C, biot_number, fin_effectiveness, heat_rate_W, max_abs_stress_MPa, max_disp_mm, safety_factor_vm, stress1_MPa, tresca_inner_MPa and von_mises_inner_MPa.

  • Targets are standardized to zero mean and unit variance before fitting.

Metrics are computed in the original output units, except for the 13 log-fitted outputs, which are scored in log space. Otherwise a handful of extreme test values would dominate \(R^2\) and the range normalization. For example, two test points of safety_factor_vm exceed 2000 while its median is 3.9. --score-space original scores everything in raw units instead.

Sampling and training

Training sets are nested random subsamples of the fixed training split. For each seed \(s\), one random permutation of the training rows is drawn, and the first \(n\) rows form the training set of size \(n\). All models and all outputs see the same rows for a given \((s, n)\), so model comparisons are paired, and a larger \(n\) only ever adds points.

Algorithm 9 Benchmark protocol for one dataset

Require: \((X_{\mathrm{tr}}, Y_{\mathrm{tr}})\), \((X_{\mathrm{te}}, Y_{\mathrm{te}})\) with \(Y \in \mathbb{R}^{N \times q}\), models \(\mathcal{M}\), seeds \(\mathcal{S} = \{0, \dots, 4\}\), sizes \(\mathcal{N} = \{50, 100, 200, 400\}\)
Ensure: Metrics \(\bigl(R^2, \mathrm{NRMSE}, \mathrm{NMAX}\bigr)_{m, s, n, j}\) for every model, seed, size and output
1:for \(s \in \mathcal{S}\) do
2:\(\pi^{(s)} \gets\) random permutation of \(\{1, \dots, N_{\mathrm{tr}}\}\) seeded with \(s\)
3:for \(n \in \mathcal{N}, \; j = 1, \dots, q, \; m \in \mathcal{M}\) do▷ independent jobs, run in parallel
4:\(I \gets \{\pi^{(s)}_1, \dots, \pi^{(s)}_n\}\)
5:\(t \gets g_j\bigl(Y_{\mathrm{tr}}[I, j]\bigr), \quad z \gets (t - \bar{t}) / \operatorname{sd}(t)\) ▷ \(g_j = \log\) or identity
6:\(\widehat{f} \gets \mathcal{A}_m\bigl(X_{\mathrm{tr}}[I],\, z\bigr)\) ▷ seed \(s\) also drives Kriging restarts and CV folds
7:\(\widehat{t}_{\mathrm{te}} \gets \bar{t} + \operatorname{sd}(t)\, \widehat{f}(X_{\mathrm{te}})\)
8:score \(\widehat{t}_{\mathrm{te}}\) against \(g_j\bigl(Y_{\mathrm{te}}[:, j]\bigr)\) (log-fitted outputs) or \(g_j^{-1}(\widehat{t}_{\mathrm{te}})\) against \(Y_{\mathrm{te}}[:, j]\)
9:end for
10:end for
LaTeX source
\begin{algorithm}
\caption{Benchmark protocol for one dataset}
\begin{algorithmic}[1]
\Require $(X_{\mathrm{tr}}, Y_{\mathrm{tr}})$, $(X_{\mathrm{te}}, Y_{\mathrm{te}})$ with $Y \in \mathbb{R}^{N \times q}$, models $\mathcal{M}$, seeds $\mathcal{S} = \{0, \dots, 4\}$, sizes $\mathcal{N} = \{50, 100, 200, 400\}$
\Ensure Metrics $\bigl(R^2, \mathrm{NRMSE}, \mathrm{NMAX}\bigr)_{m, s, n, j}$ for every model, seed, size and output
\For{$s \in \mathcal{S}$}
    \State $\pi^{(s)} \gets$ random permutation of $\{1, \dots, N_{\mathrm{tr}}\}$ seeded with $s$
    \For{$n \in \mathcal{N}, \; j = 1, \dots, q, \; m \in \mathcal{M}$} \Comment{independent jobs, run in parallel}
        \State $I \gets \{\pi^{(s)}_1, \dots, \pi^{(s)}_n\}$
        \State $t \gets g_j\bigl(Y_{\mathrm{tr}}[I, j]\bigr), \quad z \gets (t - \bar{t}) / \operatorname{sd}(t)$ \Comment{$g_j = \log$ or identity}
        \State $\widehat{f} \gets \mathcal{A}_m\bigl(X_{\mathrm{tr}}[I],\, z\bigr)$ \Comment{seed $s$ also drives Kriging restarts and CV folds}
        \State $\widehat{t}_{\mathrm{te}} \gets \bar{t} + \operatorname{sd}(t)\, \widehat{f}(X_{\mathrm{te}})$
        \State score $\widehat{t}_{\mathrm{te}}$ against $g_j\bigl(Y_{\mathrm{te}}[:, j]\bigr)$ (log-fitted outputs) or $g_j^{-1}(\widehat{t}_{\mathrm{te}})$ against $Y_{\mathrm{te}}[:, j]$
    \EndFor
\EndFor
\end{algorithmic}
\end{algorithm}

Metrics

For an output with test values \(y_i\) and predictions \(\widehat{y}_i\), \(i = 1, \dots, n_{\mathrm{te}}\), and range \(\Delta y = \max_i y_i - \min_i y_i\):

\[ R^2 = 1 - \frac{\sum_i (y_i - \widehat{y}_i)^2}{\sum_i (y_i - \bar{y})^2}, \qquad \mathrm{NRMSE} = \frac{1}{\Delta y}\sqrt{\frac{1}{n_{\mathrm{te}}}\sum_i (y_i - \widehat{y}_i)^2}, \qquad \mathrm{NMAX} = \frac{\max_i |y_i - \widehat{y}_i|}{\Delta y}. \]

\(R^2\) measures the explained variance, NRMSE the typical error relative to the output’s range, and NMAX the worst-case error, which matters in optimization because an optimizer is drawn towards regions where the surrogate is wrong.

Acceptance criteria

An output passes when the median over the five seeds, at the largest training size, satisfies all three thresholds:

\[ \operatorname{median}_s R^2 \ge 0.90, \qquad \operatorname{median}_s \mathrm{NRMSE} \le 0.10, \qquad \operatorname{median}_s \mathrm{NMAX} \le 0.30. \]

A model passes a dataset only if every output of that dataset passes. Multi-output surrogates often look fine on average but miss one output, and this rule doesn’t let an average hide that. The thresholds are command-line options (--r2-min, --nrmse-max, --nmax-max).

Execution

The reported run has 7 datasets × 51 outputs × 4 models × 5 seeds × 4 sizes = 4080 fits. Every fit is an independent job in a ProcessPoolExecutor with 24 workers. Each worker is limited to one BLAS thread so that the processes don’t oversubscribe the CPU, and the most expensive jobs are scheduled first. None of the 4080 fits failed.

The run took 117 minutes of wall time (46.7 CPU-hours). About 97% of that was spent fitting Kriging and the ensemble at \(n = 400\). The default sizes are therefore 50/100/200, and \(n = 400\) is opt-in.

Reproducing the results

pip install -e ".[bench]"

# the reported run (about 2 h on 24 cores); omit --train-sizes for a run of minutes
python BM/benchmark.py --train-sizes 50 100 200 400

# rebuild every figure on the results page from BM/results/raw_results.csv
python BM/make_doc_figures.py

BM/benchmark.py writes raw_results.csv (one row per fit), summary.csv (mean, s.d., median and quartiles over seeds), criteria.csv (per-output pass/fail) and report.md to BM/results/.