Contents

I have a concrete cantilever beam. I know the material, the load, and the geometry, but none of them are exact: concrete modulus varies batch to batch, the load has tolerance, construction gives me slightly different dimensions each time. So the deflection, stress, and rotation are not single numbers. They are distributions.

The questions I want to answer: How spread out are my outputs? Which input uncertainty matters most? Where should I spend effort to reduce risk?

This post walks through the full process: define the uncertain inputs, build a cheap polynomial surrogate (PCE), extract sensitivity rankings (Sobol indices), and interpret what the numbers actually mean for design decisions.

Interactive: PCE dashboard

Parameter Mean CoV $E$ (GPa) $q$ (kN/m) $L$ (m) $b_f$ (mm) $h_w$ (mm)
First-order Sobol indices $S_i$ (%)
Top PCE terms by variance contribution

Step 1: the structural model

A cantilever beam, fixed at one end, carries a uniformly distributed load $q$ over its span $L$. I want three outputs:

\[\delta = \frac{qL^4}{8EI} \tag{1}\] \[\sigma = \frac{qL^2\,c}{2I} \tag{2}\] \[\theta = \frac{qL^3}{6EI} \tag{3}\]

where $\delta$ is the tip deflection, $\sigma$ is the maximum bending stress at the wall, and $\theta$ is the tip rotation. The cross-section is a concrete T-beam (flange width $b_f$, flange thickness $t_f = 100$ mm, web height $h_w$, web thickness $t_w = 200$ mm). For simplification we neglect reinforcement and use gross section properties. The second moment of area $I$ follows from the parallel-axis theorem:

\[\bar{y} = \frac{A_w\,\tfrac{h_w}{2} + A_f(h_w + \tfrac{t_f}{2})}{A_w + A_f}, \qquad I = \sum_i \left(\frac{b_i h_i^3}{12} + A_i\,d_i^2\right) \tag{4}\]

where $d_i$ is each rectangle’s centroid distance from $\bar{y}$.

One thing to notice: $\sigma$ does not depend on $E$. This will show up clearly in the Sobol matrix.

Step 2: what is uncertain and by how much?

Parameter Symbol Distribution Mean CoV Why this distribution?
Elastic modulus $E$ Lognormal 30 GPa 15% Must be positive, concrete is highly variable
Distributed load $q$ Normal 25 kN/m 10% Load tolerance, symmetric around nominal
Span length $L$ Normal 4 m 2% Construction tolerance, small
Flange width $b_f$ Normal 600 mm 3% Formwork tolerance
Web height $h_w$ Normal 500 mm 3% Formwork tolerance

Setting CoV = 0 in the dashboard makes a parameter deterministic. Try it to isolate the effect of individual inputs.

Step 3: build the polynomial surrogate (PCE)

Why not just run Monte Carlo directly? For this simple beam formula, I could. But if the model were a finite element analysis taking 10 minutes per run, 10,000 MC samples would take 70 days. PCE trains on a handful of runs and gives me a cheap polynomial that replaces the model.

Each uncertain input $X_i$ is mapped to a standard normal variable $\xi_i$:

  • Normal inputs: $\xi_i = (X_i - \mu_i) / \sigma_i$
  • Lognormal $E$: $\xi_1 = (\ln E - \lambda)/\zeta$ where $\zeta^2 = \ln(1 + \text{CoV}_E^2)$, $\lambda = \ln\mu_E - \zeta^2/2$

The output is then approximated as a polynomial in these $\xi_i$:

\[Y(\boldsymbol{\xi}) \approx \sum_{|\boldsymbol{\alpha}| \leq p} c_{\boldsymbol{\alpha}}\, \Psi_{\boldsymbol{\alpha}}(\boldsymbol{\xi}) \tag{5}\]

where $\Psi_{\boldsymbol{\alpha}}$ is a product of univariate Hermite polynomials, one per input. The multi-index $\boldsymbol{\alpha} = (\alpha_1, \ldots, \alpha_5)$ controls which polynomial degree applies to each variable.

How many terms? For $d = 5$ inputs and polynomial order $p$:

\[P = \binom{d+p}{p} \tag{6}\]
Order $p$ 1 2 3 4
Terms $P$ 6 21 56 126

Step 4: train it

I generate $N = 3P$ sample points using Latin Hypercube Sampling (LHS) in standard normal space, evaluate the beam model at each, and solve a least-squares system:

\[\hat{\mathbf{c}} = (\boldsymbol{\Psi}^T\boldsymbol{\Psi})^{-1}\boldsymbol{\Psi}^T\mathbf{y} \tag{7}\]

The design matrix $\boldsymbol{\Psi}$ has one row per sample, one column per basis term. The same matrix works for all three outputs ($\delta$, $\sigma$, $\theta$): only the right-hand side $\mathbf{y}$ changes.

For order $p = 3$: 56 terms, 168 training samples. That is 168 model evaluations to build a surrogate that replaces all future model calls.

Step 5: what does the polynomial look like?

For $p = 2$ the expansion has 21 terms, grouped by role:

\[Y \approx \underbrace{c_{\mathbf{0}}}_{mean} + \underbrace{\sum_{i} c_i\,\xi_i}_{linear} + \underbrace{\sum_{i} c_{ii}(\xi_i^2 - 1)}_{curvature} + \underbrace{\sum_{i<j} c_{ij}\,\xi_i\xi_j}_{interactions} \tag{8}\]
  • $c_0$: the mean output (close to the deterministic value, but not identical due to nonlinearity)
  • Linear terms $c_i \xi_i$: capture how each input shifts the output up or down
  • Curvature terms $c_{ii}(\xi_i^2 - 1)$: capture nonlinear effects like $\delta \propto 1/E$
  • Interaction terms $c_{ij}\xi_i\xi_j$: capture coupled effects between inputs

Open the “Top PCE terms” panel in the dashboard to see the actual coefficient values.

Step 6: extract Sobol sensitivity indices

Once I have the coefficients, I get sensitivity indices for free. The first-order Sobol index $S_i$ tells me what fraction of the output variance is caused by input $i$ alone:

\[S_i = \frac{\sum_{\boldsymbol{\alpha}:\,\alpha_i>0,\;\alpha_{j\neq i}=0} \boldsymbol{\alpha}!\,c_{\boldsymbol{\alpha}}^2}{\sum_{\boldsymbol{\alpha}\neq\mathbf{0}} \boldsymbol{\alpha}!\,c_{\boldsymbol{\alpha}}^2} \tag{9}\]

where $\boldsymbol{\alpha}! = \prod_i \alpha_i!$. The total index $S_i^T$ adds all terms involving $\xi_i$ (including interactions). When $S_i \approx S_i^T$, interactions are weak.

How to read the results

“My deterministic deflection is 4.84 mm but the mean is 4.98 mm. Why are they different?”

Because $\delta \propto 1/E$ is a convex function. When $E$ is uncertain, the average of $1/E$ is larger than $1/(\text{average } E)$ (Jensen’s inequality). The 3% bias is real and a deterministic analysis would miss it.

“The CoV of my output (21%) is larger than any single input CoV. How?”

Multiple uncertainties compound. Even though $E$ has 15% CoV and $q$ has 10%, the deflection sees $q \cdot L^4 / E$, so the variances add up (roughly: $\sqrt{0.15^2 + 0.10^2 + (4 \times 0.02)^2 + \ldots} \approx 0.20$).

“I want to reduce deflection uncertainty. Where should I invest?”

Look at the Sobol matrix, $\delta$ row. $E$ accounts for ~50% of the variance. If I can tighten the concrete modulus specification (better QC, more testing), I cut the largest contributor. Tightening $L$ (construction tolerance) barely helps because it only contributes ~14% despite the $L^4$ amplification.

“Why does $E$ not appear in the stress row?”

Because $\sigma = qL^2 c / (2I)$ has no $E$ in it. Stress depends on geometry and load, not stiffness. This is a well-known structural mechanics result, and the Sobol matrix confirms it numerically: $S_E = 0.0\%$ for $\sigma$.

“Why is flange width $b_f$ negligible for everything?”

Two reasons: (1) its CoV is only 3%, and (2) $b_f$ enters $I$ roughly linearly (through the flange area term). Compare with $h_w$, which enters as $h_w^3$ in the web contribution to $I$, making it ~12% despite the same 3% CoV.

“How do I know the PCE surrogate is accurate?”

The comparison table shows analytical (Monte Carlo) vs PCE mean and standard deviation. At order $p = 3$ with default parameters, both mean and std errors are typically below 0.5%. The density curves in the plot should overlap almost perfectly. If they diverge, increase the PCE order.

“What happens if I set all CoVs to zero except one?”

The problem reduces to a single uncertain input. The Sobol matrix will show 100% for that input and 0% for all others. The density plot becomes the 1D distribution of the output driven by that one parameter alone. This is a useful sanity check.

Does this make physical sense?

A quick engineering check on the default results:

  • $\delta_0 = 4.84$ mm for a 4 m concrete cantilever under 25 kN/m: $L/\delta \approx 826$, reasonable for a stiff T-section.
  • $\sigma_0 = 13.2$ MPa: well below concrete compressive strength (~30 MPa), consistent with a serviceability-governed design.
  • $E$ dominates deflection uncertainty: concrete modulus is notoriously variable (15% CoV is standard for C30 concrete). This matches field experience where deflection predictions have the widest scatter.
  • $q$ dominates stress uncertainty: makes sense because stress is linear in $q$, and load uncertainty (10%) is the largest remaining contributor once $E$ drops out.
  • $L^4$ amplification: the span enters the deflection formula to the fourth power, so even 2% construction tolerance on $L$ generates ~8% sensitivity in $\delta$. This is why span tolerance matters more than you might expect.

Python

import numpy as np
from scipy.stats import norm
from scipy.stats.qmc import LatinHypercube
from math import factorial

TF, TW = 0.1, 0.2  # m

def hermite(n, x):
    if n == 0: return np.ones_like(x, dtype=float)
    if n == 1: return x.astype(float).copy()
    h0, h1 = np.ones_like(x, dtype=float), x.astype(float).copy()
    for k in range(2, n + 1):
        h0, h1 = h1, x * h1 - (k - 1) * h0
    return h1

def multi_index(d, p):
    out = []
    def rec(depth, rem, cur):
        if depth == d:
            out.append(tuple(cur)); return
        for k in range(rem + 1):
            cur[depth] = k; rec(depth + 1, rem - k, cur)
    rec(0, p, [0] * d)
    return out

def compute_I(bf, hw):
    Af, Aw = bf * TF, TW * hw
    yb = (Aw * hw / 2 + Af * (hw + TF / 2)) / (Af + Aw)
    I = (TW * hw**3 / 12 + Aw * (hw / 2 - yb)**2
         + bf * TF**3 / 12 + Af * (hw + TF / 2 - yb)**2)
    c = max(yb, hw + TF - yb)
    return I, c

def f_delta(E, q, L, bf, hw):
    I, _ = compute_I(bf, hw); return q * L**4 / (8 * E * I)

def f_sigma(E, q, L, bf, hw):
    I, c = compute_I(bf, hw); return q * L**2 * c / (2 * I)

def f_theta(E, q, L, bf, hw):
    I, _ = compute_I(bf, hw); return q * L**3 / (6 * E * I)

# Parameters: (mean, cov, lognormal?)
params = [
    (30e9, 0.15, True), (25e3, 0.10, False),
    (4.0, 0.02, False), (0.6, 0.03, False), (0.5, 0.03, False),
]
names = ['E', 'q', 'L', 'b_f', 'h_w']
d = 5

def to_physical(xi):
    X = np.empty_like(xi)
    for i, (mu, cov, ln) in enumerate(params):
        if ln:
            z2 = np.log(1 + cov**2)
            X[:, i] = np.exp(np.log(mu) - z2 / 2 + np.sqrt(z2) * xi[:, i])
        else:
            X[:, i] = mu + mu * cov * xi[:, i]
    return X

p_order = 3
alphas = multi_index(d, p_order)
P = len(alphas)
N = 3 * P

xi_train = norm.ppf(LatinHypercube(d, seed=42).random(N))
X_train = to_physical(xi_train)

models = [
    ('delta (mm)', f_delta, 1e3),
    ('sigma (MPa)', f_sigma, 1e-6),
    ('theta (mrad)', f_theta, 1e3),
]

for label, fn, scale in models:
    y = np.array([fn(*X_train[k]) for k in range(N)]) * scale

    Psi = np.ones((N, P))
    for j, a in enumerate(alphas):
        for i in range(d):
            if a[i] > 0: Psi[:, j] *= hermite(a[i], xi_train[:, i])

    c, *_ = np.linalg.lstsq(Psi, y, rcond=None)

    norms = [np.prod([factorial(a) for a in alpha]) for alpha in alphas]
    var_pce = sum(norms[j] * c[j]**2 for j in range(1, P))

    S, ST = np.zeros(d), np.zeros(d)
    for j in range(1, P):
        w = norms[j] * c[j]**2
        active = [i for i in range(d) if alphas[j][i] > 0]
        for i in active: ST[i] += w
        if len(active) == 1: S[active[0]] += w
    S /= var_pce; ST /= var_pce

    print(f"\n{label}:  mean={c[0]:.3f},  std={np.sqrt(var_pce):.3f}")
    for i in range(d):
        print(f"  {names[i]:>3s}:  S={S[i]:.1%},  ST={ST[i]:.1%}")

References

[1] Ghanem, R. G. & Spanos, P. D. (1991). Stochastic Finite Elements: A Spectral Approach. Springer.

[2] Xiu, D. & Karniadakis, G. E. (2002). The Wiener-Askey polynomial chaos for stochastic differential equations. SIAM Journal on Scientific Computing, 24(2), 619-644.

[3] Sudret, B. (2008). Global sensitivity analysis using polynomial chaos expansions. Reliability Engineering & System Safety, 93(7), 964-979.

[4] Sobol, I. M. (1993). Sensitivity estimates for nonlinear mathematical models. Mathematical Modelling and Computational Experiments, 1(4), 407-414.