Shape optimization of trusses

Shape optimization adjusts the positions of nodes (vertices) in a fixed truss topology to improve structural performance.

In this notebook we optimize a bridge truss to minimize compliance (maximize stiffness) under downward loads, while keeping the total volume of material constant. The design variables are the vertical positions of the top chord nodes.

from math import sqrt
from scipy.optimize import bisect
import torch
import matplotlib.pyplot as plt

from torchfem import Truss
from torchfem.materials import IsotropicElasticity1D

torch.set_default_dtype(torch.double)

# Create material
material = IsotropicElasticity1D(500.0)

Truss definition

We define a bridge truss with $A = 17$ nodes per chord (top and bottom), giving a $17 \times 2$ node grid. The connectivity consists of:

• Horizontal members connecting adjacent nodes along each chord
• Vertical members connecting the top and bottom chords at each column
• Diagonal members that alternate direction depending on position, forming the characteristic V-shaped cross-bracing

Boundary conditions:

• Left end: pinned support (both directions fixed)
• Right end: roller support (only vertical direction fixed)
• Downward point loads $f = -0.1$ applied at all intermediate bottom chord nodes

In the initial configuration, the top chord is straight and flat. The optimization will reshape the top chord into an efficient arch.

# Dimensions
A = 17
B = 2

# Nodes
n1 = torch.linspace(0.0, 5.0, A)
n2 = torch.linspace(0.0, 0.5, B)
n1, n2 = torch.stack(torch.meshgrid(n1, n2, indexing="xy"))
nodes = torch.stack([n1.ravel(), n2.ravel()], dim=1)

# Elements
connections = []
for i in range(A - 1):
    for j in range(B):
        connections.append([i + j * A, i + 1 + j * A])
for i in range(A):
    for j in range(B - 1):
        connections.append([i + j * A, i + A + j * A])
for i in range(A - 1):
    for j in range(B - 1):
        if i >= (A - 1) / 2:
            connections.append([i + j * A, i + 1 + A + j * A])
        else:
            connections.append([i + 1 + j * A, i + A + j * A])
elements = torch.tensor(connections)

# Truss
bridge = Truss(nodes.clone(), elements, material)

# Forces at bottom edge
bridge.forces[1 : A - 1, 1] = -0.1

# Constraints by the supports
bridge.constraints[0, :] = True
bridge.constraints[A - 1, 1] = True

# Plot
bridge.plot()

Method of Moving Asymptotes (MMA)

The Method of Moving Asymptotes (Svanberg, 1987) is one of the most widely used approximation methods in structural optimization. The optimization belongs to the family of sequential convex approximation methods: at each iteration it builds a convex approximation of each function around the current point, solves that approximation exactly, and moves to the solution.

Approximation

At iteration $k$, MMA replaces the original function $f$ by a separable, convex approximation:

$$\tilde{f}^{(k)}(\mathbf{x}) = f(\mathbf{x}^{(k)}) - \sum_{j} \left(\frac{p_j^{(k)}}{U_j^{(k)} - x_j^{(k)}} + \frac{q_j^{(k)}}{x_j^{(k)} - L_j^{(k)}}\right) + \sum_{j} \left(\frac{p_j^{(k)}}{U_j^{(k)} - x_j} + \frac{q_j^{(k)}}{x_j - L_j^{(k)}}\right)$$

where the moving asymptotes $L_j^{(k)} < x_j < U_j^{(k)}$ act as poles that push $x_j$ away from the asymptote boundaries. The coefficients $p_j^{(k)}, q_j^{(k)} \geq 0$ are derived from the gradient $\nabla f(\mathbf{x}^{(k)})$:

$$p_j^{(k)} = \begin{cases} (U_j^{(k)} - x_j^{(k)})^2 \cdot \frac{\partial f}{\partial x_j} & \text{if } \frac{\partial f}{\partial x_j} \geq 0 \\ 0 & \text{otherwise} \end{cases}, \qquad q_j^{(k)} = \begin{cases} -(x_j^{(k)} - L_j^{(k)})^2 \cdot \frac{\partial f}{\partial x_j} & \text{if } \frac{\partial f}{\partial x_j} < 0 \\ 0 & \text{otherwise} \end{cases}$$

This guarantees $\tilde{f}^{(k)}$ is convex and agrees with $f$ to first order at $\mathbf{x}^{(k)}$.

Asymptote update heuristic

The asymptotes are updated each iteration based on whether each design variable is oscillating (sign change in consecutive steps) or not:

• Oscillating $\Rightarrow$ contract the asymptotes (factor $s < 1$) to slow convergence and improve stability
• Non-oscillating $\Rightarrow$ expand the asymptotes (factor $s^{-1/2} > 1$) to accelerate convergence

Implementation

The MMA function computes the MMA approximation for a given function in three steps:

1. Gradient computation - torch.autograd.grad(func(x_lin), x_lin) performs one forward pass through the function and one backward pass to obtain $\nabla f(\mathbf{x}^{(k)})$. This works for any differentiable func, including the full FEM solve.
2. Coefficient assembly - $p_j^{(k)}$ and $q_j^{(k)}$ are assembled from the gradient and the current asymptotes.
3. Approximation closure - the returned approximation callable evaluates $\tilde{f}(\mathbf{x})$ at any point $L_j^{(k)} < x_j < U_j^{(k)}$.

def MMA(func, x_k, L_k, U_k):
    x_lin = x_k.clone().requires_grad_()
    grads = torch.autograd.grad(func(x_lin), x_lin)[0]
    pg = grads >= 0.0
    ng = grads < 0.0
    f_k = func(x_k)

    def approximation(x):
        p = torch.zeros_like(grads)
        p[pg] = (U_k[pg] - x_k[pg]) ** 2 * grads[pg]
        q = torch.zeros_like(grads)
        q[ng] = -((x_k[ng] - L_k[ng]) ** 2) * grads[ng]
        return (
            f_k
            - torch.sum(p / (U_k - x_k) + q / (x_k - L_k))
            + torch.sum(p / (U_k - x) + q / (x - L_k))
        )

    return approximation, grads

Helper function to compute bar lengths

The volume constraint requires the total volume $V = \mathbf{a} \cdot \mathbf{l}(\mathbf{x})$, where $\mathbf{a}$ is the vector of cross-sectional areas and $\mathbf{l}(\mathbf{x})$ is the vector of element lengths. The lengths are an explicit function of the node positions $\mathbf{x}$:

$$l_e(\mathbf{x}) = \|\mathbf{x}_{n_2} - \mathbf{x}_{n_1}\|_2$$

Since this expression is built with PyTorch operations, its gradient $\frac{\partial l_e}{\partial \mathbf{x}}$ is available for free via autograd - no manual differentiation needed.

def compute_lengths(truss):
    start_nodes = truss.nodes[truss.elements[:, 0]]
    end_nodes = truss.nodes[truss.elements[:, 1]]
    dx = end_nodes - start_nodes
    return torch.linalg.norm(dx, dim=-1)

Design variables and bounds

The design variables $\mathbf{x} \in \mathbb{R}^{A}$ are the vertical coordinates of the top chord nodes. We extract them using a boolean mask that selects exactly these degrees of freedom from the full node tensor.

The bounds $[\mathbf{x}^-, \mathbf{x}^+]$ allow nodes to move 0.4 units downward or 0.5 units upward from their initial positions. The asymmetric range gives more room to form an arch while preventing the top chord from crossing the bottom chord.

$V_0$ records the initial total volume (sum of area × length over all elements) so the volume constraint $\mathbf{a} \cdot \mathbf{l}(\mathbf{x}) \leq V_0$ can be enforced throughout the optimization.

# Mask for design variables.
mask = torch.zeros_like(nodes, dtype=torch.bool)
mask[A : 2 * A, 1] = True

# Limits on design variables
x_0 = nodes[mask].ravel()
x_min = x_0 - 0.4
x_max = x_0 + 0.5

# Compute current volume
l = compute_lengths(bridge)
V0 = torch.inner(bridge.areas, l)

Optimization problem

We minimize the structural compliance $C(\mathbf{x}) = \mathbf{f} \cdot \mathbf{u}(\mathbf{x})$, which is the work done by the external forces. Minimizing compliance is equivalent to maximizing stiffness - a stiffer structure deflects less under the same load, so a smaller $\mathbf{u}$ leads to smaller $C$.

$$\min_{\mathbf{x}} \quad C(\mathbf{x}) = \mathbf{f} \cdot \mathbf{u}(\mathbf{x}) \\\text{s.t.} \quad \mathbf{a} \cdot \mathbf{l}(\mathbf{x}) - V_0 \le 0\\\quad \quad \mathbf{x} \in [\mathbf{x}^{-}, \mathbf{x}^{+}]$$

Passing differentiable_parameters=x to bridge.solve() tells torch-fem to keep the full computation graph through the entire FEM pipeline (stiffness matrix assembly, boundary condition application, and the linear solve) as a differentiable function of x. PyTorch autograd can then backpropagate through this graph to compute $\frac{\partial C}{\partial \mathbf{x}}$ without any manual sensitivity derivation. This is the central advantage of using torch-fem for gradient-based shape optimization.

def f(x):
    # Update nodes
    bridge.nodes[mask] = x
    # Solve truss with updated nodes
    u_k, f_k, _, _, _ = bridge.solve(differentiable_parameters=x)
    # Return compliance
    return torch.inner(f_k.ravel(), u_k.ravel())

def g(x):
    # Update nodes
    bridge.nodes[mask] = x
    # Return constraint function
    return torch.inner(bridge.areas, compute_lengths(bridge)) - V0

Optimization loop

Each iteration of the loop follows the MMA workflow:

1. Update asymptotes — for the first two iterations the asymptotes are initialized from the variable ranges. From iteration 3 onward the heuristic contraction/expansion rule is applied based on whether each variable is oscillating.
2. Build approximations — MMA(f, ...) and MMA(g, ...) each call torch.autograd.grad once to compute $\nabla C$ and $\nabla g$ at the current point, then assemble the convex surrogates.
3. Solve the subproblem analytically — because $\tilde{f}$ and $\tilde{g}$ are separable, the Lagrangian $\tilde{f} + \mu\tilde{g}$ can be minimized analytically in each variable $x_j$ for a given multiplier $\mu$. The four cases correspond to the four sign combinations of the gradients of f and g.
4. Bisect for $\mu^*$ — the optimal Lagrange multiplier satisfies the dual feasibility $\tilde{g}(\mathbf{x}^*(\mu)) = 0$. Since the dual function is monotone in $\mu$, a simple bisection over $[10^{-10}, 1]$ is sufficient.
5. Update $\mathbf{x}^{(k+1)} = \mathbf{x}^*(\mu^*)$.

s = 0.8

# Set up lists for L, U, x
L = []
U = []
x = [x_0]

# Define the lower bound and upper bound of "mu"
mu_lower = torch.tensor([1e-10])
mu_upper = torch.tensor([1.0])

for k in range(100):
    # Update asymptotes with heuristic procedure (see Exercise 04)
    if k <= 1:
        L.append(x[k] - s * (x_max - x_min))
        U.append(x[k] + s * (x_max - x_min))
    else:
        L_k = torch.zeros_like(L[k - 1])
        U_k = torch.zeros_like(U[k - 1])
        # Shrink all oscillating asymptotes
        osci = (x[k] - x[k - 1]) * (x[k - 1] - x[k - 2]) < 0.0
        L_k[osci] = x[k][osci] - s * (x[k - 1][osci] - L[k - 1][osci])
        U_k[osci] = x[k][osci] + s * (U[k - 1][osci] - x[k - 1][osci])
        # Expand all non-oscillating asymptotes
        L_k[~osci] = x[k][~osci] - 1 / sqrt(s) * (x[k - 1][~osci] - L[k - 1][~osci])
        U_k[~osci] = x[k][~osci] + 1 / sqrt(s) * (U[k - 1][~osci] - x[k - 1][~osci])
        L.append(L_k)
        U.append(U_k)

    # Compute lower and upper move limit in this step
    x_min_k = torch.max(x_min, 0.9 * L[k] + 0.1 * x[k])
    x_max_k = torch.min(x_max, 0.9 * U[k] + 0.1 * x[k])

    # Compute the current approximation function and save gradients
    f_tilde, f_grad = MMA(f, x[k], L[k], U[k])
    g_tilde, g_grad = MMA(g, x[k], L[k], U[k])

    # Define x_star by minimizing the Lagrangian w. r. t. x analytically
    def x_star(mu):
        x_hat = torch.zeros_like(f_grad)
        # Case 1
        cond1 = (f_grad > 0.0) & (g_grad > 0.0)
        x_hat[cond1] = x_min[cond1]
        # Case 2
        cond2 = (f_grad < 0.0) & (g_grad < 0.0)
        x_hat[cond2] = x_max[cond2]
        # Case 3
        cond3 = (f_grad > 0.0) & (g_grad < 0.0)
        roots = torch.sqrt(
            (-mu * g_grad[cond3] * (x[k][cond3] - L[k][cond3]) ** 2)
            / (U[k][cond3] - x[k][cond3]) ** 2
            / f_grad[cond3]
        )
        x_hat[cond3] = (U[k][cond3] * roots + L[k][cond3]) / (1 + roots)
        # Case 4
        cond4 = (f_grad < 0.0) & (g_grad > 0.0)
        roots = torch.sqrt(
            (-mu * g_grad[cond4] * (U[k][cond4] - x[k][cond4]) ** 2)
            / (x[k][cond4] - L[k][cond4]) ** 2
            / f_grad[cond4]
        )
        x_hat[cond4] = (U[k][cond4] + L[k][cond4] * roots) / (1 + roots)
        return torch.clamp(x_hat, x_min_k, x_max_k)

    # Gradient of dual function
    def dual_grad(mu):
        return g_tilde(x_star(mu)).detach()

    # Solve for optimal "mu" with bisection method
    mu_star = bisect(dual_grad, mu_lower, mu_upper)

    # Compute current optimal point with dual solution
    x.append(x_star(mu_star).detach())

Results

Each curve in the following plot shows the history of one design variable (vertical position of one top-chord node) over the 100 iterations. It shows, that the design variables have converged to a stationary solution.

# Plot the development of design variables
plt.plot(torch.stack(x).detach())
plt.xlabel("Iteration")
plt.ylabel("Values $x_i$")
plt.show()

The optimized truss forms a characteristic arch shape: the top chord rises toward the center. This result - discovered purely by minimizing compliance with no geometric prior - matches classical results from bridge engineering and confirms that gradient-based shape optimization with torch-fem can recover physically meaningful optimal forms.

# Plot the optimized bridge
bridge.plot(node_labels=False)