First-Order Bilevel Topology Optimization for Fast Mechanical Design

Several emerging techniques over the past decade are contributing to the new trend of lightweight, low-cost mechanical designs. Commodity-level additive manufacturing has been made possible using 3D printers and user-friendly modeling software, allowing a novice to prototype devices, consumer products, and low-end robots. While the appearance of a mechanical part is intuitive to perceive and modify, other essential properties, such as structural robustness, fragility, vibration frequencies, and manufacturing cost, are relatively hard to either visualize or manage. These properties are highly related to the ultimate performance criteria, including flexibility, robotic maneuverability, and safety. For example, the mass and inertia of a collaborative robot should be carefully controlled to avoid harm to the human when they collide. For manipulators on spaceships, the length of each part should be optimized to maximize the end-effector reachability under strict weight limits. Therefore, many mechanical designers rely on TO methods to search for the optimal structure of mechanical parts by tuning their center-of-mass [1], inertia [2], infill pattern [3], or material distribution [4].

Unfortunately, the current TO algorithms [6] pose a major computational bottleneck. Indeed, TO algorithms rely on FEA to predict the relationship between the extrinsic (boundary loads) and intrinsic status (stress, strain, and energy densities) of a mechanical part. Such predictions are made by solving systems of (non-)linear equations during each iteration of TO. The design space and resulting system size are often gigantic when designing fine-grained material shapes. For example, an $128\times 256$ design space illustrated in Figure 1 incurs $10^{4}$ decision variables. Prior works make full use of high-end GPU [7] or many-core CPU [8, 9] to solve such problems. Even after such aggressive optimization, the total time spent on FEA is typically at the level of hours. As a result, mechanical engineers cannot modify or preview their designs without waiting for tens of minutes or even hours, which significantly slow down the loops of design refinement.

Main Results: We propose a low-cost algorithm that solves a large subclass of TO problems known as Solid Isotropic Material with Penalization (SIMP) [10]. The SIMP model is widely used to search for fine-grained geometric shapes by optimizing the infill levels. Our main observation is that, in conventional TO solvers, the computation resources are predominated by solving the system of equations to the machine precision. We argue that pursuing such precision for intermediary TO iterations is unnecessary as one only expects the system to be exactly solved on final convergence. Therefore, we propose to maintain an approximate solution that is refined over iterations. We show that such a TO solver corresponds to a bilevel reformulation of the SIMP problems (Section III). We further propose a first-order algorithm to solve these bilevel problems with guaranteed convergence to a first-order critical point (Section IV) of moderate accuracy. Compared with conventional TO algorithms [6], our method incurs a much lower iterative cost of $\mathcal{O}(E\log E)$ ($E$ is the number of decision variables), allowing the mechanical designers to interactively preview the TO results and accelerate the loops of design refinement. We then discuss several useful extensions including different preconditioning schemes for refining the approximate solutions (Section IV-A) and a GPU-based implementation that achieves an asymptotic iterative cost of $\mathcal{O}(\log E)$ (Section IV-C). We highlight the effectiveness of our method through a row of 2D benchmark problems from [11, 12, 13] (Section V), where our iterative cost can be as low as $11-40$ms when handling a resolution of up to $110$k grid cells.

We review representative TO formulations, numerical solvers for TO problems, and bilevel optimizations.

TO Formulations & Applications: Generally speaking, TO involves all mathematical programming problems that use geometric shapes as decision variables. A TO problem has three essential component: geometric representations, governing system of equations, and design objectives/constraints. The SIMP model [10] uses density-based representation where the geometric shape is discretized using FEA and decision variables are the infill levels of discrete elements. There are other popular representations. For example, truss optimization [14] uses the decision variables to model the discrete beam existence. Infill pattern optimization [3] searches for the local arrangements of shape interiors. Compared with the SIMP model, these formulations typically induce much fewer decision variables due to simplicity of shapes or locality.

In terms of governing system of equations, the SIMP model assumes linearized elastic models, which means that the extrinsic and intrinsic properties of a mechanical part are linearly related. Linear models make fairly accurate predictions under the small-deformation assumption, which is the case with most mechanical design problems with rigid parts. For large deformations, there are extensions to the SIMP model that use nonlinear elastic models [15, 16] or even elastoplastic models [17]. However, the use of these nonlinear, non-convex models induces optimization problems of a more complex class than the SIMP model [18], for which our analysis does not apply. In a similar fashion as SIMP, truss optimization [19] has small- and large-deformation variants. We speculate that similar ideas and analysis apply for truss optimization problems under small-deformations.

TO formulations deviate from each other mostly due to their design objectives and constraints. The SIMP model minimizes the total potential energy under a single external load and a total-volume constraint. Standard SIMP model does not have constraints to limit deformations, but there are extensions to include strain/stress constraints [20]. Unfortunately, stress-constrained formulations also induce a more complex problem class, where our analysis does not apply. Yet another extension to the SIMP model [21, 22, 23] is to consider multiple or uncertain external loads. Such extension is favorable because a mechanical system can undergo a variety of unknown extrinsic situations. We speculate it is possible to adapt our method to these cases by sampling the uncertain situations and use our method as underlying solvers for each sampled scenario.

Numerical TO Solvers: Three categories of numerical approaches have been extensively used: truncated gradient method, general-purpose method, and hierarchical method. The reference implementation [10] of the SIMP solver uses the heuristic, truncated gradient method [24]. This method first computes the exact gradient and then updates the design along the negative gradient direction using heuristic step sizes. Although the convergence of these heuristic schemes is hard to analyze, they perform reasonably well in practice and are thus widely used, e.g. in [12]. General-purpose methods are off-the-shelf optimization algorithms (e.g., sequential quadratic programming (SQP)/interior point method [25], and augmented Lagrangian method (ALM) [26]). The generality of these solvers makes it easier to experiment with different TO formulations. Without problem-specific optimization, however, they typically scale poorly to the high-dimensional decision space of practical problems. In particular, the (Globally Convergent) Method of Moving Asymptotes (GCMMA) [27] is a special form of SQP that can handle a large set of box constraints, making it particularly suitable for SIMP-type problems. Finally, solving problem in a hierarchical approach is an attempting way to tackle high-dimensional problems. [28] [28] used a quadtree to focus computational resources on the material boundaries. The multigrid method [12, 29] constructs a hierarchy of grids to accelerate FEA system solves.

Bilevel Optimization [18]: When a constraint of an optimization problem (high-level problem) involves another optimization (low-level problem), the problem is considered bilevel. In our formulation, the high-level problem optimizes the material infill levels, while the low-level problem performs the FEA analysis. Bilevel problems frequently arise in robotics for time optimality [30], trajectory search [31], and task-and-motion planning [32]. As a special case considered in our method, if the low-level problem is strictly convex, then its solution is unique and the low-level problem can be replaced by a well-defined solution mapping function. This treatment essentially reduces the bilevel to a single-level optimization. Otherwise, low-level problem can be non-convex and the solution map is multi-valued, where additional mechanism is needed to compare and select low-level solutions (e.g., a maximum- or minimum-value function can be used and the associated problem is NP-hard [33]). It is noteworthy that oftentimes an arbitrary solution of the low-level problem suffice, in which case penalty method [34] can be used to approximate one of the solutions on the low-level. Bilevel formulation has been used in TO community for years, but prior works resort to off-the-shelf solvers [35] or machine learning [36] to (approximately) solve the underlying optimization problem, while we utilize the special structure of the SIMP model to design our efficient, first-order method.

In this section we introduce the SIMP model using notations in [12]. (see Table II in our appendix for a list of symbols) This model uses FEA to discretize a mechanical part governed by the linear elastic constitutive law. If the material is undergoing an internal displacement $u(x):\Omega\to\mathbb{R}^{2}$ where $\Omega$ is the material domain, then the internal potential energy is accumulated according to the constitutive law as follows:

where $P[\bullet,\bullet]$ is the potential energy functional, $k_{e}$ is the spatially varying stiffness tensor parameterized by infinite-dimensional decision variable $v(x):\Omega\to\mathbb{R}$ that determines the infill levels. FEA discretizes the material domain $\Omega$ into a set of $E$ elements each spanning the sub-domain $\Omega_{e}$. The elements are connected by a set of $N$ nodes. By restricting $P[\bullet]$ to a certain finite-dimensional functional space (the shape space), the potential energy can be written as the following sum over elements:

Here $K_{e}$ is the stiffness matrix of the $e$th element, $K$ is the stiffness matrix assembled from $K_{e}$. With a slight abuse of notation, we denote $P(\bullet,\bullet)$ as a finite-dimensional potential energy function, $u\in\mathbb{R}^{2N}$ without bracket as a vector of 2D nodal displacements ($u\in\mathbb{R}^{3N}$ in 3D), and $v\in\mathbb{R}^{E}$ without bracket as a vector of elementwise infill levels. The SIMP model is compatible with all kinds of FEA discretization scheme and the case with regular grid is illustrated in Figure 2.

Suppose an external force $f$ is exerted on the mechanical part, then the total (internal+external) potential energy is:

SIMP model works with arbitrary force setup, but external forces are only applied to the boundary in practice. For example, if only the $i$th boundary node is under a force $f_{i}$, then we have $f=e_{i}f_{i}$ with $e_{i}$ being the unit vector. The equilibrium state is derived by minimizing $P_{f}$ with respect to $u$, giving $u_{f}=K^{-1}(v)f$ (Figure 2a). As the name suggests the two variables $u_{f},f$ are linearly related under the linear elastic constitutive law. The goal of TO is to minimize the induced internal potential energy due to $f$:

where $X$ is a compact, convex set bounding $v$ away from potentially singular configurations. Several widely used stiffness matrix parameterization of $K(v)$ includes the linear law $K(v)=\sum_{e=1}^{E}v_{e}K_{e}$ or the power law $K(v)=\sum_{e=1}^{E}v_{e}^{\eta}K_{e}$ where $\eta>0$ is some constant. In both cases, we need to choose $X$ such that the spectrum of $K$ is bounded from both sides:

for any $v\in X$, where $\rho(\bullet)$ is the eigenvalues of a matrix. A typical choice of $X$ to this end is:

where $v_{e}$ is the infill level of the $e$th element, $\underline{v}$ is the minimal elementwise infill level, $\bar{v}$ is the total infill level, and $\mathbb{1},\mathbb{0}$ are the all-one and all-zero vectors, respectively. In this work, we use the following more general assumption on the stiffness matrix parameterization and the shape of $X$:

All the conventional methods such as the (GC)MMA [38, 27], SQP [5, 25], and [12] for solving Equation 1 involve computing the exact gradient via sensitivity analysis (Figure 2b):

where we have used the derivatives of the matrix inverse. The major bottleneck lies in the computation of $u_{f}$ that involves exact matrix inversion. Here we assume that $\frac{\partial^{2}{K}}{\partial{v}^{2}}$ is a third-order tensor of size $|u|\times|u|\times|v|$ and its left- or right-multiplication means contraction along the first and second dimension. Note that the convergence of GCMMA and SQP solvers rely on the line-search scheme that involves the computation of objective function values, which in turn requires matrix inversion. Practical methods [12] and [10] avoid the line-search schemes using heuristic step sizes. Although working well in practice, the theoretical convergence of these heuristic rules is difficult to establish.

In this section, we first propose our core framework to solve SIMP problems that allows inexact matrix inversions. We then discuss extensions to accelerate convergence via preconditioning (Section IV-A) and fast projection operators (Section IV-B). Finally, we briefly discuss parallel algorithm implementations (Section IV-C). Inspired by the recent advent of first-order bilevel optimization algorithms [39, 40, 41], we propose to reformulation Equation 1 as the following bilevel optimization:

The low-level part of Equation 3 is a least-square problem in $u$ and, since $K(v)$ is always positive definite when $v\in X$, the low-level solution is unique. Plugging the low-level solution into the high-level objective function and Equation 1 is recovered. The first-order bilevel optimization solves Equation 3 by time-integrating the discrete dynamics system as described in Algorithm 1 and we denote this algorithm as First-Order Bilevel Topology Optimization (FBTO). The first line of Algorithm 1 is a single damped Jacobi iteration for refining the low-level solution using an adaptive step size of $\beta_{k}$. But instead of performing multiple Jacobi iterations until convergence, we go ahead to use the inexact result after one iteration for sensitivity analysis. Finally, we use a projected gradient descend to update the infill levels with an adaptive step size of $\alpha_{k}$.

Here $\text{Proj}_{X}(\bullet)$ is the projection onto the convex set $X$ under Euclidean distance. This scheme has linear iterative cost of $\mathcal{O}(E)$ as compared with conventional method that has superlinear iterative cost due to sparse matrix inversions (some operations of Algorithm 1 have a complexity of $\mathcal{O}(N)$ but we know that $E=\Theta(N)$).

As our main result, we show that Algorithm 1 is convergent under the following choices of parameters.

We prove Theorem IV.6 in Section VIII where we further show that the high-level optimality error scales as $\mathcal{O}(k^{m-1})$, so that $m$ should be as small as possible for the best convergence rate. Direct calculation leads to the constraint of $m>3/4$ from the first two conditions of Assumption IV.5, so the convergence rate can be arbitrarily close to $\mathcal{O}(k^{-1/4})$ as measured by the following relative projection error [42]:

Although our analysis of convergence uses a similar technique as that for the two-timescale method [40], our result is only single-timescale. Indeed, the low-level step size $\beta_{k}$ can be constant and users only need to tune a decaying step size for $\alpha_{k}$. We will further show that certain versions of our algorithm allow a large choice of $\beta_{k}=1$ (see Section IX-E). In addition, our formula for choosing $\alpha_{k}$ does not rely on the maximal number of iterations of the algorithm.

In this section, we compare our most efficient Algorithm 3 with other methods that can provide convergence guarantee. We implement our method using native C++ on a laptop machine with a 2.5G 6-core Intel i7 CPU and a 1.48G Nvidia GTX 1060 mobile GPU having 1280 cores. Our implementation makes full use of the cores on CPU via OpenMP and GPU via Cuda and we implement the regular grid FEA discretization. In all the examples, we apply a $7\times 7$ separable Gaussian kernel $\mathbb{C}$ as our material filter (see Remark III.2), which can be implemented efficiently on GPU as a product of multiple linear kernels. We found that a $7\times 7$ kernel achieves the best balance between material thickness and boundary sharpness. Unless otherwise stated, we choose Algorithm 3 with $\underline{v}=0.1,\eta=3,\alpha_{k}=\alpha_{0}k^{-3/4},\beta_{k}=1,D=20$ for our Krylov-preconditioner in all the experiments where $\alpha_{0}$ varies across different experiments. We justify these parameter choices in Section IX-E. For the three algorithms in Section IV, we terminate when the relative change to infill levels over two iterations ($\|v_{k+1}-v_{k}\|_{\infty}$) is smaller than $10^{-4}$ and the absolute error of the linear system solve ($\|K(v_{k})u_{k}-f\|_{\infty}$) is smaller than $10^{-2}$, which can be efficiently checked on GPU by using a reduction operator that costs $\mathcal{O}(\log E)$. The two conditions are measuring the optimality of the high- and low-level problems, respectively. The later is otherwise known as the tracking error, which measures how closely the low-level solution $u_{k}$ follows its optima as the high-level solution $v_{k}$ changes. We choose 7 benchmark problems from [11, 12, 13] to profile our method as illustrated in Figure 3 and summarized in Table I. Our iterative cost is at most $40$ms when handling a grid resolution of $192\times 576$, which allows mechanical designers to quickly preview the results.

We present a novel approach to solve a large subclass of TO problems known as the SIMP model. Unlike classical optimization solvers that perform the exact sensitivity analysis during each iteration, we propose to use inexact sensitivity analysis that is refined over iterations. We show that this method corresponds to the first-order algorithm for a bilevel reformulation of TO problems. Inspired by the recent analytical result [39], we show that this approach has guaranteed convergence to a first-order critical point of the original problem. Our new approach leads to low-iterative cost and allows interactive result preview for mechanical designers. We finally discuss extensions to use preconditioners and fine-grained GPU parallelism. Experiments on several 2D benchmarks from [11, 12, 13] highlight the computational advantage of our method.

An inherent limitation of our method is the requirement to estimate step size parameters $\alpha_{k},\beta_{k}$. We have shown that $\beta_{k}$ is fixed and our results are not sensitive to $\alpha_{0}$, but users still need to give a rough estimation of $\alpha_{0}$ that can affect practical convergence speed. Further, our first-order solver is based on plain gradient descendent steps, for which several acceleration schemes are available, e.g. Nesterov [54] and Anderson [55] accelerations. The analysis of these techniques are left as future works. Finally, when additional constraints are incorporated into the polytopic set $X$ or non-polytopic constraints are considered, our analysis needs to be adjusted and the efficacy of Algorithm 4 might be lost.

We compare the performance of FBTO and PGD for all the benchmarks in Figure 3 and the results are summarized in Figure 11. Note that PGD uses an exact low-level solver so it requires less iterations then FBTO, while FBTO uses more iterations with a much lower iterative cost. For fairness we compare them based on computational time. PGD outperforms or performs similarly to FBTO on two of the problems (Figure 11df), where PGD converges rapidly in $5-10$ iterations while FBTO requires thousands of iterations. In all other cases, FBTO converges faster to an approximate solution.

We provide proof of Theorem IV.6. We start from some immediate observations on the low- and high-level objective functions. Then, we prove the rule of error propagation for the low-level optimality, i.e., the low-level different between $u_{k}$ and its optimal solution $u_{f}(v_{k})$. Next, we choose parameters to let the difference diminish as the number of iterations increase. Finally, we focus on the high-level objective function and show that the difference between $v_{k}$ and a local optima of $l(v)$ would also diminish. Our analysis resembles the recent analysis on two-timescale bilevel optimization, but we use problem-specific treatment to handle our novel gradient estimation for the high-level problem without matrix inversion.