RasmahRasmah

Robust (Polynomial Chaos) Topology Optimization

println("Julia ", VERSION)
Julia 1.12.6

Overview

Real loads are never known exactly. Robust topology optimization designs a structure that performs well across a range of loads, not just one. Here the applied load is a smooth function f(ξ) of a standard-normal random variable ξ ~ N(0,1), so the compliance

\[C(\xi) = f(\xi)^{\mathsf T} K(x)^{-1} f(\xi)\]

is itself a random response. Its Hermite polynomial chaos expansion (PCE) coefficients (Ghanem & Spanos 1991; Xiu & Karniadakis 2002)

\[a_k = \frac{\mathbb{E}[C\,H_k]}{k!} = \frac{1}{k!}\sum_q w_q\, C(\xi_q)\, H_k(\xi_q)\]

are evaluated by Gauss–Hermite quadrature, giving the mean a₀ and the variance Σ_{k≥1} aₖ² k!. The optimizer minimizes the robust objective mean(C) + κ·std(C), trading off average performance against its variability.

using Rasmah
using LinearAlgebra

E, ν = 1e3, 0.3
m = tetrahedralize_box(2.0, 0.5, 0.5, 8, 2, 2)
U = FESpace(m, ReferenceFE(3, 1); valuetype = VectorValue{3})
n = ndofs(U); nn = size(m.nodes, 2)

# clamp the left face, pull the right face in x
fixed = Int[]; f0 = zeros(n)
for i in 1:nn
    m.nodes[1, i] < 1e-12 && append!(fixed, 3i - 2, 3i - 1, 3i)
    m.nodes[1, i] > 2 - 1e-12 && (f0[3i - 2] = 1.0)
end

# random load magnitude: f(ξ) = f0 (1 + 0.2 ξ + 0.1 ξ²) — a quartic compliance
fload = ξ -> f0 .* (1.0 + 0.2ξ + 0.1ξ^2)
#2 (generic function with 1 method)

pce_robust_topopt projects C(ξ) onto the Hermite chaos at nq Gauss–Hermite nodes and minimizes mean + κ·std. With κ = 0 it reduces to the deterministic compliance; a positive κ (risk aversion) penalizes load variability.

x, hist = pce_robust_topopt(m, E, ν, fixed, fload, 0.5, 1.0; p = 4, nq = 8, max_iter = 40)

println("robust objective: ", round(hist[1]; digits = 3), " → ", round(hist[end]; digits = 3))
println("final volume frac = ", round(sum(x) / length(x); digits = 3))

vtk(m; field = x, fieldname = "material density", edges = true)
robust objective: 10.189 → 2.68
final volume frac = 0.501

Verification

For a linear load f(ξ) = f₀ + f₁ ξ the compliance C(ξ) = C₀₀ + 2C₀₁ξ + C₁₁ξ² is quadratic in ξ, whose Hermite coefficients are a₀ = C₀₀ + C₁₁, a₁ = 2C₀₁, a₂ = C₁₁. The mean and variance then reduce to the closed forms C₀₀ + C₁₁ and 4C₀₁² + 2C₁₁² — the same statistics the driver recovers by quadrature:

C00, C01, C11 = 1.0, 0.3, 0.5
a  = [C00 + C11, 2C01, C11]            # Hermite coefficients
μ  = pce_compliance_mean(a)
σ² = pce_compliance_variance(a)
println("PCE mean = ", round(μ; digits = 6), "   (analytic C₀₀ + C₁₁ = ", C00 + C11, ")")
println("PCE var  = ", round(σ²; digits = 6), "   (analytic 4C₀₁² + 2C₁₁² = ", 4C01^2 + 2C11^2, ")")
PCE mean = 1.5   (analytic C₀₀ + C₁₁ = 1.5)
PCE var  = 0.86   (analytic 4C₀₁² + 2C₁₁² = 0.86)

Next steps