Shape Optimization and Topological Optimization

From Formula to Form: A Practical Task

An engineer is designing a load-bearing beam of an aircraft wing. Requirements: withstand a given load and weigh as little as possible. Where to remove material? The shape of the wing is optimized precisely as the solution of a variational problem—to minimize the volume (mass) under strength constraints. This is shape optimization. Every day, Boeing, Airbus, and Toyota use it in their calculations. Without variational calculus, modern aircraft engineering would be impossible.

Formulation of the Shape Optimization Problem

Given a domain $\Omega \subset \mathbb{R}^n$—the "body"—with boundary $\Gamma = \partial\Omega$. A partial differential equation (deformation, flow, heat) is solved on $\Omega$. The problem: find a shape $\Omega$ minimizing a functional $J(\Omega)$.

Examples of functionals:

• $J(\Omega) = \int_\Omega u ,dx$ (average displacement under load)
• $J(\Omega) = \int_{\partial\Omega} q^2 ,dS$ (heat flux through the surface)
• $J(\Omega) = \max_{x \in \Omega} \sigma(x)$ (maximum stress)
• $J(\Omega) = |\Omega|$ (volume of material) under strength constraints

A challenge: the region $\Omega$ is an infinite-dimensional object. How to take a "derivative" with respect to the shape?

Method of Boundary Variations (Hadamard's Formula)

A family of domains: $\Omega_t$, deformed along a velocity vector $V(x)$ on the boundary:

$ x_t = x + t \cdot V(x), \quad \Omega_t = { x_t : x \in \Omega } $

Shape derivative (Hadamard's formula): $\left.\frac{dJ}{dt}\right|{t=0} = \int{\partial\Omega} j(x) \cdot (V \cdot n) , dS$

where $n$ is the outward normal, $j(x)$ is the "sensitivity" to shape deformation. The derivative depends only on the normal component of velocity!

Optimality condition: for the optimal shape, $j(x) = \text{const}$ on $\partial\Omega$ (or equal to the Lagrange multiplier under an isoperimetric constraint).

Example: For a diffusion problem ($-\Delta u = f$ in $\Omega$, $u = 0$ on $\partial\Omega$), $J = \int_\Omega u , dx$: the shape derivative $\frac{dJ}{dt} = -\int_{\partial\Omega} \left( \frac{\partial u}{\partial n} \right)^2 V \cdot n , dS$. Optimality condition: $\left(\frac{\partial u}{\partial n} \right)^2 = \text{const}$ on $\partial\Omega$—the flux through the optimal boundary is constant.

Topological Optimization: The SIMP Method

Classical shape optimization alters only the boundary (shape). Topological optimization also changes topology—new holes may appear or old ones be eliminated.

Density approach: $\rho(x) \in [0, 1]$—"material density" at each point. $\rho = 1$: material present. $\rho = 0$: void.

SIMP method (Solid Isotropic Material with Penalization):

Modulus of elasticity: $E(\rho) = \rho^p \cdot E_0$ ($p \approx 3$—penalty on intermediate densities)

Task: $\min_{\rho \in [0,1]} \int_\Omega \rho , d\Omega$ (material volume) subject to compliance constraints (structure stiffness).

Equivalent: $\min_{(u, \rho)} \int_\Omega E(\rho) , \varepsilon(u):\varepsilon(u) , d\Omega$ with $\int_\Omega \rho , d\Omega \le V_0$, plus the equilibrium equation.

Algorithm:

1. Set initial distribution $\rho = V_0/|\Omega|$ (homogeneous)
2. FEM calculation of displacement field $u$
3. Compute "sensitivity": $\frac{\partial C}{\partial\rho_i} = -p \rho_i^{p-1} \varepsilon_i : \varepsilon_i$
4. Update $\rho$ (optimality criterion or gradient method)
5. Apply filter for regularization
6. Repeat until convergence to a "black-and-white" distribution ($\rho \approx 0$ or $1$)

Result—a "lacy" structure with minimal mass for the given stiffness. Similar to the bones of birds!

Aerodynamic Optimization: Adjoint Method

Task: Find the shape of a wing (section profile) minimizing aerodynamic drag under a given lift force.

Physics: the velocity field $u(x)$ satisfies the Navier-Stokes equations (or Euler for inviscid fluid). The goal: $J(\Omega) = \int_\Omega f \cdot u , dV$ (power of drag).

Problem: $J$ depends on $u$, which depends on $\Omega$. The gradient $\frac{dJ}{d\Omega}$ is needed for optimization, but brute-force calculation requires $N$ computations for $N$ shape parameters.

Adjoint method: introduce an adjoint field $p$ (solution of the "reverse" problem), then:

$ \frac{dJ}{d\Omega} = \text{function of } (u, p) \text{ on } \partial\Omega $

Instead of $N$ computations—just 2 (the direct and adjoint problems). This is a revolutionary approach.

Application: Airbus uses the adjoint method to optimize A380 wing profiles. NASA optimizes turbine blades. Fuel savings—2–5% due to optimal shaping.

Full Analysis: Topological Optimization of a Beam

Problem: Beam $2 \times 1$ m, fixed at both ends, load applied in the center from below. Material volume $50%$ of the total. Minimize maximum deflection.

Initial state: $\rho = 0.5$ everywhere. Stiffness is uniform.

After 5 iterations: "struts" appear from the supports to the loaded point, material flows out of unloaded regions.

After 50 iterations: black-and-white structure. Triangular struts stretching from the load center to the anchor points. Looks like a bridge truss!

Result: stiffness increased by 3.5 times at the same mass compared to uniform distribution. Such a design is impossible without topological optimization.

Theory of Optimal Control

The Pontryagin maximum principle (1956) generalizes the Euler-Lagrange equation to problems with control constraints. For the problem: $\dot{x} = f(x, u, t)$, $u \in U$, minimize $J = \int_0^T L(x, u, t) , dt + g(x(T))$. The Hamiltonian: $H(x,p,u,t) = p^\top f(x,u,t) - L(x,u,t)$. The maximum principle: optimal $u^(t) \in \arg\max_{u \in U} H(x^, p^*, u, t)$, where $p$ is the adjoint variable (Lagrange multiplier) satisfying $\dot{p} = -\partial H/\partial x$.

Bang-bang Control

With linear dependence of $H$ on $u$ and compact $U$, the optimal control takes boundary values of $U$: $u^$ lies "on the boundary" of the admissible set. This is bang-bang control—switching between extreme modes. For example, in the minimum-time problem for a linear system: $u^ = \operatorname{sign}(B^\top P(t) x(t))$—switching between "full throttle" and "full brake" at moments when the switching function changes sign.

Numerical Methods of Optimal Control

• Indirect methods: derive necessary conditions (maximum principle), solve the boundary value problem for the ODE system
• Direct methods: discretization and conversion to NLP—IPOPT, SNOPT
• Pseudospectral methods: representing the trajectory by a high-order polynomial (GPOPS-II, DIDO)
• Differential Dynamic Programming (DDP): iterative second-order method, basis for modern RL and trajectory planning algorithms

Applications

• Aerospace industry: optimal rocket launch trajectories (fuel minimization to reach a given orbit), interplanetary missions (Voyager, New Horizons, Cassini)
• Robotics: motion planning for manipulators and drones
• Finance: optimal consumption and investment (Merton model)
• Energy: optimal real-time control of power plants
• Medicine: optimization of drug dosage, chemotherapy scheduling
• Epidemiology: optimal vaccination and social distancing strategies

Modern control theory is a direct descendant of variational calculus, uniting its tools with the theory of differential equations and numerical methods.