using TopOpt, Zygote, Flux
using NNlib: leakyreluNeural 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
Problem Definition
E = 1.0 # Young's modulus
v = 0.3 # Poisson's ratio
f = 1.0 # downward force
els = (160, 40)
problem = PointLoadCantilever(Val{:Linear}, 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_paramsOptimizer 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.01Density Functions
function todensities(w; filter=true)
if filter
PseudoDensities(proj.(cheqfilter(nn(NNParams(w))).x))
else
PseudoDensities(proj.(nn(NNParams(w)).x))
end
endOptimization 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.7635, constr=0.3622, eps=1.0 Epoch 2: p=1.2000000000000002, α=21.0, obj=0.8306, constr=0.3151, eps=1.0 Epoch 3: p=1.3000000000000003, α=31.0, obj=0.9326, constr=0.2488, eps=1.0 Epoch 4: p=1.4000000000000004, α=41.0, obj=1.1173, constr=0.1515, eps=1.0 Epoch 5: p=1.5000000000000004, α=51.0, obj=1.4888, constr=0.0123, eps=1.0 Epoch 6: p=1.6000000000000005, α=61.0, obj=2.0739, constr=-0.1209, eps=1.0 Epoch 7: p=1.7000000000000006, α=71.0, obj=2.2265, constr=-0.1223, eps=1.0 Epoch 8: p=1.8000000000000007, α=81.0, obj=2.1216, constr=-0.0711, eps=1.0 Epoch 9: p=1.9000000000000008, α=91.0, obj=1.9379, constr=-0.0004, eps=1.0 Epoch 10: p=2.000000000000001, α=100.0, obj=1.7964, constr=0.0637, eps=1.0 Epoch 11: p=2.100000000000001, α=100.0, obj=1.7672, constr=0.0969, eps=1.0 Epoch 12: p=2.200000000000001, α=100.0, obj=1.8219, constr=0.1052, eps=1.0 Epoch 13: p=2.300000000000001, α=100.0, obj=1.9459, constr=0.0954, eps=1.0 Epoch 14: p=2.4000000000000012, α=100.0, obj=2.1449, constr=0.0716, eps=1.0 Epoch 15: p=2.5000000000000013, α=100.0, obj=2.4291, constr=0.0374, eps=1.0 Epoch 16: p=2.6000000000000014, α=100.0, obj=2.7764, constr=0.0016, eps=1.0 Epoch 17: p=2.7000000000000015, α=100.0, obj=3.0751, constr=-0.0195, eps=1.0 Epoch 18: p=2.8000000000000016, α=100.0, obj=3.2458, constr=-0.0209, eps=1.0 Epoch 19: p=2.9000000000000017, α=100.0, obj=3.2774, constr=-0.0052, eps=1.0 Epoch 20: p=3.0, α=100.0, obj=3.2178, constr=0.0211, eps=1.0 Epoch 21: p=3.0, α=100.0, obj=2.956, constr=0.0493, eps=1.0 Epoch 22: p=3.0, α=100.0, obj=2.7854, constr=0.0686, eps=1.0 Epoch 23: p=3.0, α=100.0, obj=2.7037, constr=0.0759, eps=1.0 Epoch 24: p=3.0, α=100.0, obj=2.6933, constr=0.0721, eps=1.0 Epoch 25: p=3.0, α=100.0, obj=2.7445, constr=0.0585, eps=1.0 Epoch 26: p=3.0, α=100.0, obj=2.8502, constr=0.0376, eps=1.0 Epoch 27: p=3.0, α=100.0, obj=2.9729, constr=0.0161, eps=1.0 Epoch 28: p=3.0, α=100.0, obj=3.049, constr=0.0028, eps=1.0 Epoch 29: p=3.0, α=100.0, obj=3.0373, constr=0.0016, eps=1.0 Epoch 30: p=3.0, α=100.0, obj=2.9395, constr=0.0119, eps=1.0 Epoch 31: p=3.0, α=100.0, obj=2.7982, constr=0.0293, eps=1.0 Epoch 32: p=3.0, α=100.0, obj=2.6681, constr=0.0467, eps=1.0 Epoch 33: p=3.0, α=100.0, obj=2.5914, constr=0.0566, eps=1.0 Epoch 34: p=3.0, α=100.0, obj=2.573, constr=0.0571, eps=1.0 Epoch 35: p=3.0, α=100.0, obj=2.6092, constr=0.0488, eps=1.0 Epoch 36: p=3.0, α=100.0, obj=2.6901, constr=0.0346, eps=1.0 Epoch 37: p=3.0, α=100.0, obj=2.7779, constr=0.021, eps=1.0 Epoch 38: p=3.0, α=100.0, obj=2.8135, constr=0.0153, eps=1.0 Epoch 39: p=3.0, α=100.0, obj=2.767, constr=0.0201, eps=1.0 Epoch 40: p=3.0, α=100.0, obj=2.6638, constr=0.0324, eps=1.0 Epoch 41: p=3.0, α=100.0, obj=2.5625, constr=0.0455, eps=1.0 Epoch 42: p=3.0, α=100.0, obj=2.5085, constr=0.0522, eps=1.0 Epoch 43: p=3.0, α=100.0, obj=2.5142, constr=0.0502, eps=1.0 Epoch 44: p=3.0, α=100.0, obj=2.5718, constr=0.041, eps=1.0 Epoch 45: p=3.0, α=100.0, obj=2.6508, constr=0.0299, eps=1.0 Epoch 46: p=3.0, α=100.0, obj=2.694, constr=0.024, eps=1.0 Epoch 47: p=3.0, α=100.0, obj=2.6646, constr=0.0269, eps=1.0 Epoch 48: p=3.0, α=100.0, obj=2.5857, constr=0.0361, eps=1.0 Epoch 49: p=3.0, α=100.0, obj=2.5128, constr=0.0452, eps=1.0 Epoch 50: p=3.0, α=100.0, obj=2.4838, constr=0.0485, 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.4837692825480158 constr_final(w) = 0.048451364046329104 eps = 1.0
1.0Visualization
using CairoMakie
topology = todensities(w; filter=true)
fig = visualize(problem; topology=topology)Precompiling packages... 7501.2 ms ✓ QuartoNotebookWorkerMakieExt (serial) 1 dependency successfully precompiled in 8 seconds Precompiling packages... 5227.8 ms ✓ QuartoNotebookWorkerCairoMakieExt (serial) 1 dependency successfully precompiled in 6 seconds ┌ 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
