Simultaneous Analysis and Design (SAND)

Description

Nested topology optimization solves the equilibrium equation K(x) u = f for the displacements u at every design iterate, then differentiates through that solve. Simultaneous analysis and design (SAND) instead treats the displacements u as optimization variables and enforces the equilibrium equation as a constraint:

\[ \min_{x, u}\; f^\mathsf{T} u \quad\text{s.t.}\quad K(x)\,u = f,\quad 0 \le x_e \le 1 \]

At a feasible point K u = f, the objective fᵀu equals the compliance uᵀKu, so this is the same physical problem — but the linear solve is folded into the optimization, and the equilibrium is enforced by a nonlinear solver rather than a nested factorization.

TopOpt.jl exposes the differentiable pieces needed for SAND: ElementKFun maps densities to element stiffness matrices, AssembleKFun assembles them into the global matrix, and apply_boundary_with_meandiag! applies the Dirichlet boundary conditions — all differentiable through Zygote.

Setup

using TopOpt, LinearAlgebra
using WGLMakie
WGLMakie.activate!(; resize_to=:parent)
using Bonito
if haskey(ENV, "QUARTO_PROJECT_DIR")
    Bonito.Page(exportable=true, offline=true)
else
    Bonito.browser_display()
end
display_app(app) = display(app)

WGLMakie.activate!(; resize_to=:parent) selects the browser renderer and fills the Quarto output column. Bonito.Page(exportable=true, offline=true) embeds the assets needed by visualize(...; static=true) in the Quarto output, so the visualization does not require a running Julia process.

The differentiable assembly

E, ν, f = 1.0, 0.3, 1.0
problem = PointLoadCantilever((20, 10), (1.0, 1.0), E, ν, f)
ch = problem.ch
nel = getncells(problem)

solver = FEASolver(DirectSolver, problem; xmin=0.001, penalty=PowerPenaltyFun(3.0))
solver.vars .= 1.0
solver()

fvec = copy(solver.globalinfo.f)          # assembled load (BC applied)
ndofs = length(fvec)
element_k = ElementKFun(solver)            # x -> element stiffness matrices
assemble_k = AssembleKFun(problem)         # element matrices -> global K

K(x) = apply_boundary_with_meandiag!(assemble_k(element_k(PseudoDensities(x))), ch)

We verify these blocks reproduce the solver’s own equilibrium: assembling K(x) and solving K u = f matches the nested displacement field, and the compliance fᵀu matches ComplianceFun.

x = ones(nel)
Kx = K(x)
u_nested = copy(solver.u)
u_direct = Kx \ fvec
comp = ComplianceFun(solver)
println("assembly reproduces equilibrium: |K u - f| = $(norm(Kx * u_direct - fvec))")
println("nested vs assembled u:          |Δu| = $(norm(u_nested - u_direct))")
println("f'u = $(dot(fvec, u_direct)),  ComplianceFun = $(comp(PseudoDensities(x)))")

assembly reproduces equilibrium: |K u - f| = 1.7491395730072561e-13
nested vs assembled u:          |Δu| = 0.0
f'u = 35.29121684231742,  ComplianceFun = 35.29121684232222

The augmented Lagrangian method

The equilibrium is a vector of ndofs equality constraints r(x, u) = K(x) u - f = 0. The augmented Lagrangian folds them into the objective:

\[ L(x,u;\lambda,\rho) = f^\mathsf{T} u + \lambda^\mathsf{T} r + \frac{\rho}{2}\,r^\mathsf{T} r, \qquad \lambda \leftarrow \lambda + \rho\, r,\quad \rho \leftarrow \gamma\rho. \]

For a fixed design (analysis only) the subproblem is a convex quadratic in u, so it can be solved exactly. Below we run the multiplier/penalty loop and show the equilibrium residual collapse to zero, recovering the nested displacement.

function al_subproblem(λ, ρ)
    # stationarity of L w.r.t. u:  K(λ + ρ (K u - f)) = -f   =>   u = K⁻¹f - K⁻¹(K⁻¹f + λ)/ρ
    w = Kx \ fvec
    return w .- (Kx \ (w .+ λ)) ./ ρ
end

λ = zeros(ndofs)
ρ = 1.0
u = zeros(ndofs)
for k in 1:6
    global λ, ρ, u
    u = al_subproblem(λ, ρ)
    r = Kx * u - fvec
    println("AL iter $k: |r| = $(norm(r)),  |u - u_nested| = $(norm(u - u_nested))")
    λ .+= ρ .* r
    ρ *= 10.0
end

AL iter 1: |r| = 277.208878401867,  |u - u_nested| = 647437.9746524752
AL iter 2: |r| = 4.415424300535528e-11,  |u - u_nested| = 5.086145683582884e-9
AL iter 3: |r| = 1.650512649755094e-13,  |u - u_nested| = 1.4029480611646794e-12
AL iter 4: |r| = 1.6058522997914846e-13,  |u - u_nested| = 1.657937226649138e-12
AL iter 5: |r| = 1.6338611980131567e-13,  |u - u_nested| = 3.5075010148938525e-12
AL iter 6: |r| = 1.6388752745859102e-13,  |u - u_nested| = 2.4198905860534385e-12

After the first multiplier update the equilibrium is satisfied to machine precision — in exact arithmetic the augmented Lagrangian converges in a single multiplier update for a quadratic equality, and in practice a few iterations suffice.

Feasibility of the recovered solution

The displacement u is an optimization variable, so we verify the equilibrium constraint directly: the residual K u - f must vanish, and u must match the nested solver’s displacement.

r_final = Kx * u - fvec
println("equilibrium residual |K u - f| = $(norm(r_final))")
println("recovered vs nested |Δu|     = $(norm(u - u_nested))")
println("feasible: $(norm(r_final) < 1e-8)")

equilibrium residual |K u - f| = 1.6388752745859102e-13
recovered vs nested |Δu|     = 2.4198905860534385e-12
feasible: true

fig = visualize(problem; static=true, topology=ones(nel), u=u, default_exagg_scale=0.5)
display_app(fig)

undeformed mesh

load arrows

support arrows

Figure 1: Deformed cantilever with the displacement field recovered by the augmented Lagrangian (the SAND solution)

The residual is at machine precision, so the SAND analysis solution satisfies K u = f exactly and reproduces the nested displacement field.

From analysis to design

The same augmented Lagrangian applies to the full SAND problem: add the volume constraint and optimize over (x, u) jointly, with K(x) u - f = 0 enforced as the augmented-Lagrangian term. TopOpt.jl’s first-order augmented Lagrangian solver, PercivalAlg from NonconvexPercival, does exactly this.

model = Model(z -> dot(fvec, z[nel+1:end]))
addvar!(model, vcat(zeros(nel), fill(-1e3, ndofs)), vcat(ones(nel), fill(1e3, ndofs)))
add_ineq_constraint!(model, z -> sum(z[1:nel]) / nel - 0.5)
add_eq_constraint!(model, z -> K(z[1:nel]) * z[nel+1:end] - fvec; dim=ndofs)
z_sand = optimize(model, PercivalAlg(), vcat(fill(0.5, nel), solver.u);
    options=PercivalOptions(; max_iter=25, subsolver_max_eval=250)).minimizer
x_sand = z_sand[1:nel]
u_sand = z_sand[nel+1:end]

[ Info:   iter        fx    normgp    normcx         μ     normy    sumc     inner_status        iter_type  
[ Info:      0   3.5e+01   2.3e+00   8.7e-01   1.0e+01   2.2e+01       5
[ Info:      1   3.5e+01   2.3e+00   8.7e-01   1.0e+02   2.2e+01       9      first_order         update_μ
[ Info:      2   3.5e+01   1.5e+00   1.5e-01   1.0e+02   1.4e+01     740         max_eval         update_y
[ Info:      3   3.5e+01   1.5e+00   1.5e-01   1.0e+03   1.4e+01     744         max_eval         update_μ
[ Info:      4   3.5e+01   3.6e+01   1.5e-01   1.0e+03   1.5e+02     748         max_eval         update_y
[ Info:      5   3.5e+01   3.6e+01   1.5e-01   1.0e+04   1.5e+02     752         max_eval         update_μ
[ Info:      6   3.5e+01   3.9e+02   1.5e-01   1.0e+04   1.6e+03     756         max_eval         update_y
[ Info:      7   3.5e+01   3.9e+02   1.5e-01   1.0e+05   1.6e+03     760         max_eval         update_μ
[ Info:      8   3.5e+01   3.9e+03   1.5e-01   1.0e+05   1.7e+04     764         max_eval         update_y
[ Info:      9   3.5e+01   3.9e+03   1.5e-01   1.0e+06   1.7e+04     768         max_eval         update_μ
[ Info:     10   3.5e+01   1.1e+04   1.5e-01   1.0e+06   1.7e+05     772         max_eval         update_y
[ Info:     11   3.5e+01   1.1e+04   1.5e-01   1.0e+07   1.7e+05     776         max_eval         update_μ
[ Info:     12   3.5e+01   1.8e+04   1.5e-01   1.0e+07   1.7e+06     780         max_eval         update_y
[ Info:     13   3.5e+01   1.8e+04   1.5e-01   1.0e+08   1.7e+06     784         max_eval         update_μ
[ Info:     14   3.5e+01   2.0e+04   1.5e-01   1.0e+08   1.7e+07     788         max_eval         update_y
[ Info:     15   3.5e+01   2.0e+04   1.5e-01   1.0e+09   1.7e+07     792         max_eval         update_μ
[ Info:     16   3.5e+01   2.0e+04   1.5e-01   1.0e+10   1.7e+07     796         max_eval         update_μ
[ Info:     17   3.5e+01   2.0e+04   1.5e-01   1.0e+11   1.7e+07     800         max_eval         update_μ
[ Info:     18   3.5e+01   2.0e+04   1.5e-01   1.0e+12   1.7e+07     804         max_eval         update_μ
[ Info:     19   3.5e+01   2.0e+04   1.5e-01   1.0e+13   1.7e+07     808         max_eval         update_μ
[ Info:     20   3.5e+01   2.0e+04   1.5e-01   1.0e+14   1.7e+07     812         max_eval         update_μ
[ Info:     21   3.5e+01   2.0e+04   1.5e-01   1.0e+15   1.7e+07     816         max_eval         update_μ
[ Info:     22   3.5e+01   2.0e+04   1.5e-01   1.0e+16   1.7e+07     820         max_eval         update_μ
[ Info:     23   3.5e+01   2.0e+04   1.5e-01   1.0e+17   1.7e+07     824         max_eval         update_μ

462-element Vector{Float64}:
   0.01636111637855903
   0.01636111637855903
  -1.2034417636527244
  -0.8734758156064086
  -0.7610502677316888
  -0.5034453400585096
   0.01636111637855903
   0.01636111637855903
  -2.1570522686847937
  -1.5825553304712823
   ⋮
 -24.659715534830095
  10.237010150581538
 -26.84677520156709
  10.504408455323539
 -29.070766052588272
  10.401467890042532
 -31.072628013871874
  10.831331360594
 -33.529567751262064

Feasibility of the SAND design

The optimization now jointly returns a design x_sand and a displacement u_sand. Because u is a decision variable, its feasibility is not guaranteed — we must check the equilibrium residual K(x) u - f and compare the objective fᵀu against the nested compliance (which satisfies equilibrium by construction).

r_sand = K(x_sand) * u_sand - fvec
solver.vars .= x_sand
solver()
u_nested_sand = copy(solver.u)
println("volume fraction         = $(round(sum(x_sand) / nel, digits=4))")
println("equilibrium residual    = $(round(norm(r_sand), digits=4))  (relative $(round(norm(r_sand) / norm(fvec), digits=4)))")
println("SAND objective  f'u     = $(round(dot(fvec, u_sand), digits=3))")
println("nested compliance at x  = $(round(dot(fvec, u_nested_sand), digits=3))")
println("feasible (|K u - f| ≈ 0) = $(norm(r_sand) < 1e-6)")

volume fraction         = 0.5266
equilibrium residual    = 0.1464  (relative 0.1464)
SAND objective  f'u     = 34.854
nested compliance at x  = 860.004
feasible (|K u - f| ≈ 0) = false

The recovered u is only approximate: the problem is nonconvex and the solver may not find a solution that satisfies all constraints perfectly, so a significant residual K(x)u - f means the design is not actually feasible. Consequently fᵀu no longer equals the nested compliance (which satisfies equilibrium by construction).

Discussion

• SAND exposes u as a design variable and enforces K u = f as a constraint, avoiding a nested solve per iterate at the cost of a larger, nonconvex problem.
• The building blocks (ElementKFun, AssembleKFun, apply_boundary_with_meandiag!) are differentiable, so the equilibrium residual and its gradient are available through Zygote.
• The augmented Lagrangian is the natural method for the equilibrium constraint; for the convex analysis-only case it converges immediately, while the full SAND problem relies on first-order solvers such as NonconvexPercival.PercivalAlg.