2. Python ecosystem for crashworthiness surrogate modeling

Crashworthiness optimization sits at the intersection of nonlinear finite element analysis (FEA), design of experiments (DOE), and statistical learning. A single LS-DYNA frontal impact run on a full vehicle body-in-white may consume 6–24 wall-clock hours on 32–64 cores, so a surrogate model that predicts intrusion, peak acceleration, or mass within milliseconds is not a convenience — it is the only way to make multi-objective optimization tractable. This chapter maps the Python toolchain that supports that workflow, from raw solver output to a trained, validated, and queried metamodel.

2.1 The Layered Stack

It helps to think of the ecosystem as four cooperating layers rather than a flat list of libraries. The numerical core (NumPy, SciPy) provides array algebra, linear solvers, and interpolation. The data layer (pandas, plus solver-specific parsers) handles tabular design matrices and time-history curves. The learning layer (scikit-learn, PyTorch, GPy, SMT) fits and validates the surrogate. The orchestration layer (SALib, job schedulers, solver Python APIs) generates DOE samples, dispatches simulations, and collects responses.

2.2 Numerical Core: NumPy and SciPy

Every surrogate is ultimately a function mapping a design vector x ∈ R^n to a response vector y ∈ R^m. NumPy supplies the ndarray that stores the design matrix X (shape n_samples × n_features) and the response matrix Y. For Kriging, the training step requires solving a linear system K α = y where K is the n × n covariance matrix; numpy.linalg.solve or a Cholesky factorization via scipy.linalg.cho_factor handles this. SciPy adds the optimization machinery — scipy.optimize.minimize with L-BFGS-B is the standard choice for maximum-likelihood estimation of Kriging hyperparameters, and scipy.optimize.differential_evolution is useful when the likelihood surface is multimodal.

A practical detail: crash responses often span several orders of magnitude (e.g., mass in kg, intrusion in mm, acceleration in g). Always normalize X and Y to zero mean and unit variance before fitting, otherwise the hyperparameter optimizer is dominated by the largest-magnitude feature.

2.3 Data Layer: pandas and Solver Output Parsing

pandas is the natural container for the design matrix and responses. A typical crashworthiness dataset looks like this:

import pandas as pd

df = pd.read_csv("doe_results.csv")
# columns: thickness_B_pillar, thickness_rail, material_yield,
#          mass, intrusion_mm, peak_accel_g, energy_absorbed_kJ
X = df[["thickness_B_pillar", "thickness_rail", "material_yield"]].values
Y = df[["mass", "intrusion_mm", "peak_accel_g"]].values
print(df.describe())

Getting data out of the solver is the harder half. LS-DYNA writes binary d3plot and binout files; the lasso-python package (formerly lasso) reads both and exposes nodal displacement, element stress, and rigid-wall forces as NumPy arrays. Radioss output is typically read through pyradioss or by parsing T01/T02 ASCII files. Abaqus users can call abaqus python scripts that use the odbAccess module to extract field output from an .odb, then dump to CSV for downstream pandas ingestion. A robust pattern is to keep a thin parser module per solver that returns a tidy DataFrame with one row per design point and one column per scalar response.

2.4 Learning Layer: Choosing the Right Surrogate

The choice of learner depends on sample budget and response smoothness. With 20–80 FEA runs, Kriging (Gaussian process regression) is usually the best default because it provides both a prediction and an uncertainty estimate — essential for infill sampling. SMT (Surrogate Modeling Toolbox) is the most crashworthiness-oriented option: it offers ordinary Kriging, universal Kriging, gradient-enhanced Kriging (GEK), and co-Kriging for multi-fidelity data, all with a scikit-learn-compatible API. GPy gives finer control over kernel composition and is preferred when you need custom covariance structures or sparse GP approximations for larger datasets. scikit-learn covers the broader baseline — GaussianProcessRegressor, RandomForestRegressor, GradientBoostingRegressor, and MLPRegressor — and is the right tool for cross-validation utilities like KFold and metrics like r2_score. PyTorch enters when the response is a field (e.g., full deformation map) or when you want a physics-informed neural network that embeds the governing momentum equation as a soft constraint.

2.5 Orchestration: SALib and Solver Integration

SALib generates the DOE and performs sensitivity analysis. A Latin hypercube or Sobol sequence is created with salib.sample.latin.sample or salib.sample.sobol.sample, then each row is written into a solver input deck. After the runs complete, salib.analyze.sobol.analyze returns first-order and total-order indices, telling you which thickness or material parameters actually drive intrusion — this prunes the design space before surrogate training.

from SALib.sample import latin
from SALib.analyze import sobol

problem = {
    "num_vars": 3,
    "names": ["t_bpillar", "t_rail", "yield"],
    "bounds": [[1.2, 2.4], [1.5, 3.0], [280, 420]],
}
X = latin.sample(problem, 64, seed=42)
# ... run 64 LS-DYNA jobs, collect intrusion_mm ...
Si = sobol.analyze(problem, Y_intrusion, print_to_console=False)
print(Si["S1"], Si["ST"])

The integration loop itself is usually a Python driver that (1) reads a template .k or .rad file, (2) substitutes parameter values, (3) launches the solver via subprocess or a scheduler API, (4) waits for termination, and (5) parses the result. For LS-DYNA, lasso-python can also write modified keyword files programmatically, avoiding brittle string replacement.

2.6 A Minimal End-to-End Skeleton

Putting the layers together, a reproducible surrogate pipeline looks like this:

import numpy as np, pandas as pd
from SALib.sample import latin
from smt.surrogate_models import KRG
from sklearn.model_selection import KFold
from sklearn.metrics import r2_score

# 1. DOE
X = latin.sample(problem, 64, seed=42)
# 2. Simulate (external) -> df with X and Y
df = pd.read_csv("runs.csv")
# 3. Fit Kriging per response
kf = KFold(n_splits=5, shuffle=True, random_state=0)
for col in ["intrusion_mm", "peak_accel_g"]:
    y = df[col].values
    scores = []
    for tr, te in kf.split(X):
        sm = KRG(theta0=[1e-2]*3, print_global=False)
        sm.set_training_values(X[tr], y[tr])
        sm.train()
        scores.append(r2_score(y[te], sm.predict_values(X[te])))
    print(col, "CV R2 =", np.mean(scores))

An R² above 0.9 on held-out folds is a reasonable acceptance threshold for intrusion and mass; peak acceleration is noisier and often lands between 0.75 and 0.88, which is acceptable if the surrogate is used for ranking rather than absolute prediction. If R² is poor, the usual causes are an under-sampled design space, a discontinuous response (e.g., a buckling mode switch), or unnormalized inputs. Adding 10–20 infill points near the region of highest Kriging variance, then refitting, typically recovers most of the lost accuracy.

The next step in any real project is to wrap this pipeline so that the surrogate can be queried inside an optimizer — but that belongs to the optimization chapters. Here the essential point is that NumPy and SciPy do the algebra, pandas and solver parsers do the data plumbing, SMT/GPy/scikit-learn/PyTorch do the learning, and SALib plus a subprocess driver close the loop back to LS-DYNA, Radioss, or Abaqus.