1. OOD Detection Baseline for LS-DYNA Simulation Outputs

Why a Mahalanobis Baseline for LS-DYNA Outputs

An LS-DYNA run produces a high-dimensional, correlated feature vector: peak nodal displacement, internal energy, hourglass energy, contact force, kinetic energy at a reference time, and so on. These features are not independent — a heavier impact raises both kinetic energy and contact force together. A naive per-feature z-score threshold will therefore either miss joint anomalies or fire constantly on normal correlated variation. The Mahalanobis distance solves this by measuring how far a run sits from the nominal population in the metric of the nominal covariance, so correlated features are handled jointly and the score is scale-free.

For a run feature vector x ∈ ℝd, with nominal mean μ and covariance Σ estimated from a reference set of nominal runs, the squared Mahalanobis distance is:

D²(x) = (x − μ)ᵀ Σ⁻¹ (x − μ)

Under a Gaussian assumption, D² follows a chi-squared distribution with d degrees of freedom, which gives a principled threshold: flag a run when D² > χ²d, 1−α. This is the baseline we implement. It is deliberately simple and interpretable — the goal of this first case is a reference detector that later fallback logic can be calibrated against, not a state-of-the-art model.

Feature Contract and Reference Set

Before any math, fix the feature contract. Every run must emit the same ordered vector, extracted at the same simulation time, with the same units. A mismatch here is the most common cause of a detector that "works" on the reference set and fails in production. For this baseline we use six features per run:

The reference set should contain at least 30–50 nominal runs (varying mesh seed, small material scatter, solver version) so that Σ is well-conditioned. With d = 6 and n = 40, the sample covariance is estimable; if n is close to d, regularize with Σ + λI.

Implementation

The detector is a small, self-contained module. It fits on the reference set, persists μ, Σ⁻¹, and the threshold, then scores new runs. The code below is complete and runnable with NumPy/SciPy.

import numpy as np
from scipy.stats import chi2

FEATURES = ["peak_disp", "int_energy", "hg_energy",
            "contact_force", "ke_20ms", "added_mass"]

class MahalanobisOOD:
    def __init__(self, alpha=0.99, reg=1e-6):
        self.alpha = alpha
        self.reg = reg
        self.mu = None
        self.prec = None      # Σ⁻¹
        self.threshold = None

    def fit(self, X_ref):
        X_ref = np.asarray(X_ref, dtype=float)
        n, d = X_ref.shape
        self.mu = X_ref.mean(axis=0)
        cov = np.cov(X_ref, rowvar=False)
        cov += self.reg * np.eye(d)          # regularize
        self.prec = np.linalg.inv(cov)
        self.threshold = chi2.ppf(self.alpha, df=d)
        return self

    def score(self, x):
        x = np.asarray(x, dtype=float)
        diff = x - self.mu
        return float(diff @ self.prec @ diff)

    def predict(self, x):
        d2 = self.score(x)
        return d2 > self.threshold, d2

# --- example reference set (illustrative values, not measured) ---
rng = np.random.default_rng(0)
n_ref = 40
# correlated nominal population: contact force tracks internal energy
base = rng.normal(0, 1, size=(n_ref, 6))
base[:, 3] = 0.8 * base[:, 1] + 0.6 * base[:, 3]   # induce correlation
X_ref = base * np.array([2.0, 50.0, 1.5, 300.0, 20.0, 0.02]) \
        + np.array([12.0, 800.0, 20.0, 4500.0, 150.0, 0.10])

det = MahalanobisOOD(alpha=0.99).fit(X_ref)
print("threshold D²:", round(det.threshold, 3))

# nominal run and an off-nominal run (contact force far above nominal)
x_nom = X_ref[0]
x_ood = X_ref[0].copy(); x_ood[3] *= 2.5
for name, x in [("nominal", x_nom), ("off-nominal", x_ood)]:
    flag, d2 = det.predict(x)
    print(f"{name:12s} D²={d2:8.2f}  flag={flag}")

Expected output shape (values depend on the random seed): the nominal run scores below the threshold and is not flagged; the off-nominal run, whose contact force is 2.5× nominal, produces a D² well above the threshold and is flagged. The exact numbers are illustrative — the point is the ordering and the flag decision, which is what the fallback logic consumes.

Reading the Score: Which Feature Drove the Flag

A single scalar is useful for triggering but useless for diagnosis. Decompose the score into per-feature contributions to see what pushed a run out of the nominal cloud. The contribution of feature i is diff_i · (Σ⁻¹ diff)_i; summing over i recovers D². This tells an engineer whether the anomaly is a contact-force spike (likely a modeling or contact-definition issue) or an energy imbalance (likely hourglass or mass-scaling).

def contributions(det, x):
    diff = np.asarray(x) - det.mu
    contrib = diff * (det.prec @ diff)
    return dict(zip(FEATURES, np.round(contrib, 2)))

print(contributions(det, x_ood))

Calibration and Validation

Two checks decide whether the baseline is trustworthy before it feeds any fallback logic.

• False-positive rate on held-out nominal runs. Split the reference set (e.g., 30 fit / 10 hold-out). With α = 0.99, expect roughly 1% of nominal runs to exceed τ. If the hold-out rate is far higher, the Gaussian assumption is strained — inspect feature histograms for heavy tails and consider a log transform on energy features.
• Separation on known off-nominal runs. Assemble a small set of deliberately perturbed runs (e.g., doubled contact force, hourglass energy at 8% of internal energy, mass scaling above the run's tolerance). A useful baseline flags most of these while keeping the nominal false-positive rate near 1 − α. If a perturbation is not flagged, that feature is under-weighted by Σ⁻¹ — check whether it is nearly collinear with another feature.

Common Failure Symptoms

The output of this case is a persisted (μ, Σ⁻¹, τ) triple plus a score(x) function. That triple is the reference the fallback trigger will be calibrated against: the trigger's job is to decide what to do when D² > τ, while this baseline's job is to produce a defensible, reproducible D² in the first place.