# Benchmark Problems And Setup

This page describes the empirical benchmark behind the claims in {doc}`benchmark_results`: 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}
:name: alg-benchmark-protocol

\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

```bash
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/`.