Neural Network Topology Optimization with Adam Optimizer

Description

This tutorial demonstrates neural network-parametrized topology optimization using Flux.jl’s Adam optimizer. Unlike the IPOPT-based approach, this method uses first-order gradient descent with adaptive learning rates. A continuation strategy progressively increases the SIMP penalty (1→3) and constraint aggregation weight (1→100) to converge to binary designs.

The network uses 6 hidden layers with sigmoid output for density parametrization.

Setup

using TopOpt, Zygote, Flux
using NNlib: leakyrelu
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.

Problem Definition

E = 1.0 # Young's modulus
v = 0.3 # Poisson's ratio
f = 1.0 # downward force
els = (160, 40)

problem = PointLoadCantilever(els, (1.0, 1.0), E, v, f)

V = 0.5       # volume fraction
xmin = 1e-6   # minimum density
rmin = 3.0    # filter radius

# Continuation parameters: penalty ramps 1 → 3 over the epochs, and the
# constraint weight ramps from 1 to 100 so the volume constraint is
# progressively enforced.
p = 1.0           # initial SIMP penalty
delta_p = 0.1     # penalty increment per epoch
p_max = 3.0       # maximum penalty

alpha = 1.0       # initial constraint weight
delta_alpha = 10.0
alpha_max = 100   # maximum constraint weight

# Initial solver and compliance for C0 (normalization constant)
penalty = PowerPenaltyFun(p)
solver = FEASolver(DirectSolver, problem; xmin, penalty)
cheqfilter = DensityFilterFun(solver; rmin)
comp = ComplianceFun(solver)
volfrac = VolumeFun(solver)
C0 = comp(cheqfilter(PseudoDensities(fill(V, getncells(problem)))))

Neural Network

# 6-layer MLP with sigmoid output (scalar density in [0, 1])
m = 20
act = leakyrelu
nn = NeuralNetworkFun(
    Chain(
        Dense(2, m, act; init=Flux.glorot_normal),
        Dense(m, m, act; init=Flux.glorot_normal),
        Dense(m, m, act; init=Flux.glorot_normal),
        Dense(m, m, act; init=Flux.glorot_normal),
        Dense(m, m, act; init=Flux.glorot_normal),
        Dense(m, 1, sigmoid; init=Flux.glorot_normal),
    ),
    problem;
    scale=true,
)
w0 = nn.init_params

Optimizer Setup

# Adam with a moderate learning rate; gradient clipping avoids spikes.
alg = Flux.Optimise.Adam(0.01)
clip_alg = Flux.Optimise.ClipValue(1.0)
w = copy(w0)

# Heaviside projection with nonzero steepness pushes densities toward 0/1.
proj = HeavisideProjectionFun(2.0)

# Termination criteria
eps = Inf
eps_star = 0.05    # target intermediate-density fraction
maxiter = 50
constr_tol = 0.01

Density Functions

function todensities(w; filter=true)
    if filter
        PseudoDensities(proj.(cheqfilter(nn(NNParams(w))).x))
    else
        PseudoDensities(proj.(nn(NNParams(w)).x))
    end
end

Optimization Loop

epoch = 1
while true
    global epoch, p, alpha, eps, w
    epoch > maxiter && break
    eps < eps_star && break

    # Ramp penalty and constraint weight
    global p = min(p + delta_p, p_max)
    global alpha = min(alpha + delta_alpha, alpha_max)

    # Update penalty on the existing solver (mutates in-place, so the
    # compliance and volume functions pick up the new penalty automatically)
    setpenalty!(solver, p)

    obj = w -> comp(todensities(w; filter=true)) / C0
    constr = w -> volfrac(todensities(w; filter=false)) / V - 1
    combined_obj = w -> obj(w) + alpha * constr(w)^2

    # Adam step with gradient clipping
    Δ = Zygote.gradient(combined_obj, w)[1]
    Δ = Flux.Optimise.apply!(clip_alg, w, Δ)
    Flux.Optimise.update!(alg, w, Δ)

    # Check convergence
    x = nn(NNParams(w)).x
    global eps = sum(0.05 .< x .< 0.95) / length(x)
    violation = abs(constr(w))

    println("Epoch $epoch: p=$p, α=$alpha, obj=$(round(obj(w); digits=4)), constr=$(round(constr(w); digits=4)), eps=$(round(eps; digits=3))")
    global epoch += 1
end

┌ Warning: Layer with Float32 parameters got Float64 input.
│   The input will be converted, but any earlier layers may be very slow.
│   layer = Dense(2 => 20, leakyrelu)  # 60 parameters
│   summary(x) = "2-element Vector{Float64}"
└ @ Flux ~/.julia/packages/Flux/hrg9M/src/layers/stateless.jl:60
Epoch 1: p=1.1, α=11.0, obj=0.8621, constr=0.2566, eps=1.0
Epoch 2: p=1.2000000000000002, α=21.0, obj=1.1108, constr=0.1019, eps=1.0
Epoch 3: p=1.3000000000000003, α=31.0, obj=1.6165, constr=-0.0923, eps=1.0
Epoch 4: p=1.4000000000000004, α=41.0, obj=1.704, constr=-0.0875, eps=1.0
Epoch 5: p=1.5000000000000004, α=51.0, obj=1.5855, constr=-0.016, eps=1.0
Epoch 6: p=1.6000000000000005, α=61.0, obj=1.4474, constr=0.0647, eps=1.0
Epoch 7: p=1.7000000000000006, α=71.0, obj=1.4271, constr=0.0976, eps=1.0
Epoch 8: p=1.8000000000000007, α=81.0, obj=1.5263, constr=0.086, eps=1.0
Epoch 9: p=1.9000000000000008, α=91.0, obj=1.7263, constr=0.0459, eps=1.0
Epoch 10: p=2.000000000000001, α=100.0, obj=2.0227, constr=-0.0083, eps=1.0
Epoch 11: p=2.100000000000001, α=100.0, obj=2.252, constr=-0.0345, eps=1.0
Epoch 12: p=2.200000000000001, α=100.0, obj=2.3162, constr=-0.0246, eps=1.0
Epoch 13: p=2.300000000000001, α=100.0, obj=2.2522, constr=0.009, eps=1.0
Epoch 14: p=2.4000000000000012, α=100.0, obj=2.1573, constr=0.048, eps=1.0
Epoch 15: p=2.5000000000000013, α=100.0, obj=2.1397, constr=0.0701, eps=1.0
Epoch 16: p=2.6000000000000014, α=100.0, obj=2.2178, constr=0.0721, eps=1.0
Epoch 17: p=2.7000000000000015, α=100.0, obj=2.3966, constr=0.0572, eps=1.0
Epoch 18: p=2.8000000000000016, α=100.0, obj=2.6723, constr=0.0322, eps=1.0
Epoch 19: p=2.9000000000000017, α=100.0, obj=2.9621, constr=0.0132, eps=1.0
Epoch 20: p=3.0, α=100.0, obj=3.1319, constr=0.0126, eps=1.0
Epoch 21: p=3.0, α=100.0, obj=2.9522, constr=0.0283, eps=1.0
Epoch 22: p=3.0, α=100.0, obj=2.7388, constr=0.0504, eps=1.0
Epoch 23: p=3.0, α=100.0, obj=2.6064, constr=0.0648, eps=1.0
Epoch 24: p=3.0, α=100.0, obj=2.5826, constr=0.0658, eps=1.0
Epoch 25: p=3.0, α=100.0, obj=2.6563, constr=0.0547, eps=1.0
Epoch 26: p=3.0, α=100.0, obj=2.7932, constr=0.0376, eps=1.0
Epoch 27: p=3.0, α=100.0, obj=2.8967, constr=0.0257, eps=1.0
Epoch 28: p=3.0, α=100.0, obj=2.8694, constr=0.027, eps=1.0
Epoch 29: p=3.0, α=100.0, obj=2.7402, constr=0.039, eps=1.0
Epoch 30: p=3.0, α=100.0, obj=2.6112, constr=0.0523, eps=1.0
Epoch 31: p=3.0, α=100.0, obj=2.5514, constr=0.0579, eps=1.0
Epoch 32: p=3.0, α=100.0, obj=2.572, constr=0.053, eps=1.0
Epoch 33: p=3.0, α=100.0, obj=2.6565, constr=0.0402, eps=1.0
Epoch 34: p=3.0, α=100.0, obj=2.7414, constr=0.0278, eps=1.0
Epoch 35: p=3.0, α=100.0, obj=2.7382, constr=0.0259, eps=1.0
Epoch 36: p=3.0, α=100.0, obj=2.65, constr=0.0345, eps=1.0
Epoch 37: p=3.0, α=100.0, obj=2.5507, constr=0.0455, eps=1.0
Epoch 38: p=3.0, α=100.0, obj=2.506, constr=0.0498, eps=1.0
Epoch 39: p=3.0, α=100.0, obj=2.5297, constr=0.0442, eps=1.0
Epoch 40: p=3.0, α=100.0, obj=2.5921, constr=0.0331, eps=1.0
Epoch 41: p=3.0, α=100.0, obj=2.6182, constr=0.0268, eps=1.0
Epoch 42: p=3.0, α=100.0, obj=2.551, constr=0.0323, eps=1.0
Epoch 43: p=3.0, α=100.0, obj=2.4376, constr=0.0445, eps=1.0
Epoch 44: p=3.0, α=100.0, obj=2.4026, constr=0.0456, eps=1.0
Epoch 45: p=3.0, α=100.0, obj=2.4615, constr=0.0329, eps=1.0
Epoch 46: p=3.0, α=100.0, obj=2.4724, constr=0.0265, eps=1.0
Epoch 47: p=3.0, α=100.0, obj=2.3411, constr=0.0387, eps=1.0
Epoch 48: p=3.0, α=100.0, obj=2.2489, constr=0.0457, eps=1.0
Epoch 49: p=3.0, α=100.0, obj=2.327, constr=0.0274, eps=1.0
Epoch 50: p=3.0, α=100.0, obj=2.2687, constr=0.0281, eps=1.0

Results

obj_final = w -> comp(todensities(w; filter=true)) / C0
constr_final = w -> volfrac(todensities(w; filter=false)) / V - 1
@show obj_final(w)
@show constr_final(w)
@show eps

┌ Warning: Layer with Float32 parameters got Float64 input.
│   The input will be converted, but any earlier layers may be very slow.
│   layer = Dense(2 => 20, leakyrelu)  # 60 parameters
│   summary(x) = "2-element Vector{Float64}"
└ @ Flux ~/.julia/packages/Flux/hrg9M/src/layers/stateless.jl:60
obj_final(w) = 2.2687309682818877
constr_final(w) = 0.028102551417824362
eps = 1.0

1.0

Visualization

topology = todensities(w; filter=true)
fig = visualize(problem; static=true, topology=topology)
display_app(fig)

┌ Warning: Layer with Float32 parameters got Float64 input.
│   The input will be converted, but any earlier layers may be very slow.
│   layer = Dense(2 => 20, leakyrelu)  # 60 parameters
│   summary(x) = "2-element Vector{Float64}"
└ @ Flux ~/.julia/packages/Flux/hrg9M/src/layers/stateless.jl:60

undeformed mesh

load arrows

support arrows

Figure 1: Adam-optimized neural network design (6-layer MLP)