Topology optimization with linearized buckling criteria in 250 lines of Matlab

Abstract

We present a 250-line Matlab code for topology optimization for linearized buckling criteria. The code is conceived to handle stiffness, volume and buckling load factors (BLFs) either as the objective function or as constraints. We use the Kreisselmeier-Steinhauser aggregation function in order to reduce multiple objectives (viz. constraints) to a single, differentiable one. Then, the problem is sequentially approximated by using MMA-like expansions and an OC-like scheme is tailored to update the variables. The inspection of the stress stiffness matrix leads to a vectorized implementation for its efficient construction and for the sensitivity analysis of the BLFs. This, coupled with the efficiency improvements already presented by Ferrari and Sigmund (Struct Multidiscip Optim 62:2211–2228, 2020a), cuts all the computational bottlenecks associated with setting up the buckling analysis and allows buckling topology optimization problems of an interesting size to be solved on a laptop. The efficiency and flexibility of the code are demonstrated over a few structural design examples and some ideas are given for possible extensions.

Appendices

Appendix: 1. Notes

Let $X^{e} = \left [{x^{e}_{i}},{y^{e}_{i}}\right ]^{k}_{i=1}$ be the array collecting the coordinates of a k-noded element and $N(\xi ,\zeta ) = \left [ N_{1}(\xi ,\zeta ), N_{2}(\xi ,\zeta ), \ldots , N_{k}(\xi ,\zeta )\right ]$ be the array collecting the shape functions, where (ξ,ζ) are the logical coordinates (Bathe 1982). The (x,y)-gradient of the shape functions is computed as ∇_(x,y)N = J^− 1∇_(ξ,ζ)N, where the (ξ,ζ)-gradient and Jacobian are

$$ \begin{aligned} \nabla_{(\xi, \zeta)}N & = \left[ \begin{array}{ccc} \partial_{\xi} N_{1}, \ldots, \partial_{\xi} N_{k} \\ \partial_{\zeta} N_{1}, \ldots, \partial_{\zeta} N_{k} \end{array} \right] \\ J(\xi, \eta) & = \left[ \begin{array}{cc} {\sum}_{k}\partial_{\xi} N_{k} x_{k} & {\sum}_{k}\partial_{\xi} N_{k} y_{k} \\ {\sum}_{k}\partial_{\zeta} N_{k} x_{k} & {\sum}_{k}\partial_{\zeta} N_{k} y_{k} \end{array} \right] = \nabla_{(\xi, \zeta)}N X^{e} \end{aligned} $$

(29)

The discretized deformation gradient (B₁) and strain-displacement operator (B₀) are then recovered as

$$ \begin{aligned} B_{1} & = \left[ \begin{array}{c} \nabla_{(x, y)}N \otimes [ 1, 0 ]^{T} \\ \nabla_{(x, y)}N \otimes [ 0, 1 ]^{T} \end{array} \right] \\ B_{0} & = L (\nabla_{(x, y)}N \otimes I_{2} ) \end{aligned} $$

(30)

where L is the placement matrix (here shown for a $\mathcal {Q}_{4}$ element)

$$ L = \left[ \begin{array}{cccc} 1 & 0 & 0 & 0 \\ 0 & 0 & 0 & 1 \\ 0 & 1 & 1 & 0 \end{array} \right] $$

and I₂ is the identity matrix of order 2.

The operations above have a remarkably compact implementation in the code of Appendix C, for the particular case of a $\mathcal {Q}_{4}$ element (see lines extracted below).

The instructions above could be further simplified, since we are considering a uniform discretization (i.e., J is constant). However, the current definition is general enough to be easily extended to higher order elements by only modifying the definition of dN.

Appendix: 2. Details on the re-design rule and on the KS aggregation function

Let us have the two sets of functions {g_0(i)(x)}_i= 1…,q and {g_j(x)}_j= 1…,s and the following optimization problem

$$ \begin{cases} & \underset{\mathbf{x}\in\left[ 0, 1 \right]^{m}}{\min} \underset{i=1,\ldots, q}{\max} g_{0(i)}(\mathbf{x}) \\ \text{s.t.} & g_{j}\left( \mathbf{x} \right) \leq 0 , \qquad\qquad j = 1, \ldots, s \end{cases} $$

(31)

We replace $J^{KS}_{0}\left (\mathbf {x} \right ) := J^{KS}[g_{0(i)}\left (\mathbf {x} \right )]$ to the “$\max \limits$” operator in the objective, and $J^{KS}_{1}\left (\mathbf {x} \right ) := J^{KS}\left [g_{j}\left (\mathbf {x} \right )\right ]$ to the set of constraints, obtaining the single-objective, single-constraint problem

$$ \begin{cases} & \underset{\mathbf{x}\in\left[ 0, 1 \right]^{m}}{\min} J^{KS}_{0}\left( \mathbf{x} \right) \\ \text{s.t.} & J^{KS}_{1}\left( \mathbf{x} \right) \leq 0 \end{cases} $$

(32)

which is analogous to (13). At the current design point (x_k = ξ) the objective and constraint functions are expanded in terms of intervening variables y_e(x_e), e = 1,…m, each selected such to build a convex approximation. Choosing an MMA-like form for the intervening variables (Svanberg 1987), the objective (i = 0) and constraint (i = 1) are expanded as

$$ \begin{array}{@{}rcl@{}} J^{KS}_{i}(\mathbf{x}) &\approx& J^{KS}_{i}(\boldsymbol{\xi}) + {\sum}^{m}_{e=1} \left[ (U_{e}-\xi_{e})^{2}{p^{i}_{e}}(\boldsymbol{\xi}) \left( \frac{1}{U_{e}-x_{e}} - \frac{1}{U_{e}-\xi_{e}} \right)\right. \\ &&\left.- (\xi_{e}-L_{e})^{2}{q^{i}_{e}}(\boldsymbol{\xi}) \left( \frac{1}{x_{e}-L_{e}} - \frac{1}{\xi_{e}-L_{e}} \right) \right] \\ && =\underbrace{J^{KS}_{i}(\boldsymbol{\xi}) - {\sum}^{m}_{e=1} \left[ (U_{e}-\xi_{e}) {p^{i}_{e}}(\boldsymbol{\xi}) - (\xi_{e}-L_{e}) {q^{i}_{e}}(\boldsymbol{\xi})\right]}_{\widehat{J^{KS}_{i}}(\boldsymbol{\xi})} \\&&+ {\sum}^{m}_{e=1} \left[ \frac{(U_{e}-\xi_{e})^{2}{p^{i}_{e}}(\boldsymbol{\xi})} {U_{e}-x_{e}} - \frac{(\xi_{e}-L_{e})^{2}{q^{i}_{e}}(\boldsymbol{\xi})}{x_{e}-L_{e}} \right] \\ & = &\widehat{J^{KS}_{i}}(\boldsymbol{\xi}) + {\sum}^{m}_{e=1} \left[\frac{(U_{e}-\xi_{e})^{2} {p^{i}_{e}}(\boldsymbol{\xi})}{U_{e}-x_{e}} - \frac{(\xi_{e}-L_{e})^{2} {q^{i}_{e}}(\boldsymbol{\xi})}{x_{e}-L_{e}} \right] \end{array} $$

(33)

where ${p^{i}_{e}}(\boldsymbol {\xi }) = \max \limits \{\partial _{e}J^{KS}_{i}(\boldsymbol {\xi }),0\} \geq 0$ and ${q^{i}_{e}}(\boldsymbol {\xi }) = \min \limits \{\partial _{e}J^{KS}_{i}(\boldsymbol {\xi }),0\} \leq 0$ and (L_e,U_e), e = 1,…,m are the moving asymptotes, such that L_e < (x_e,ξ_e) < U_e.

Once (33) is substituted into (32), the Lagrangian of the problem reads

$$ \begin{array}{@{}rcl@{}} \ell(\mathbf{x}, \lambda ) &=& \widehat{J^{KS}_{0}}(\boldsymbol{\xi}) + \kappa \widehat{J^{KS}_{1}}(\boldsymbol{\xi}) + {\sum}^{m}_{e=1} \left[ \frac{(U_{e}-\xi_{e})^{2} ({p^{0}_{e}}(\boldsymbol{\xi})+\kappa {p^{1}_{e}}(\boldsymbol{\xi})) }{U_{e}-x_{e}}\right. \\ &&\left.- \frac{(\xi_{e}-L_{e})^{2} ({q^{0}_{e}}(\boldsymbol{\xi})+\kappa {q^{1}_{e}}(\boldsymbol{\xi}))}{x_{e}-L_{e}} \right] \end{array} $$

(34)

where κ ≥ 0 is the Lagrange multiplier. The stationarity condition of the Lagrangian with respect to x_e

$$ \partial_{e} \ell(\mathbf{x}, \lambda ) = \frac{(U_{e}-\xi_{e})^{2} ({p^{0}_{e}}(\boldsymbol{\xi})+\kappa {p^{1}_{e}}(\boldsymbol{\xi})) }{(U_{e}-x_{e})^{2}} + \frac{(\xi_{e}-L_{e})^{2} ({q^{0}_{e}}(\boldsymbol{\xi})+\kappa {q^{1}_{e}}(\boldsymbol{\xi}))}{(x_{e}-L_{e})^{2}} = 0 $$

(35)

can be solved explicitly, giving the primal update map

$$ x_{e}(\kappa) = \max\left\{\delta_{-},\min\left\{ \delta_{+}, \frac{L_{e}(U_{e}-\xi_{e})\sqrt{{p^{0}_{e}}(\boldsymbol{\xi}) + \kappa {p^{1}_{e}}(\boldsymbol{\xi})} + U_{e}(\xi_{e}-L_{e})\sqrt{-{q^{0}_{e}}(\boldsymbol{\xi}) - \kappa {q^{1}_{e}}(\boldsymbol{\xi})}} {(U_{e}-\xi_{e})\sqrt{{p^{0}_{e}}(\boldsymbol{\xi}) + \kappa {p^{1}_{e}}(\boldsymbol{\xi})} + (\xi_{e}-L_{e})\sqrt{-{q^{0}_{e}}(\boldsymbol{\xi}) - \kappa {q^{1}_{e}}(\boldsymbol{\xi})}} \right\}\right\} $$

(36)

where 0 < δ₋ < δ₊ < 1 are the adaptive move limits, depending on the asymptotes and we have highlighted the dependence on the dual variable κ. This latter can be computed from the following equation, obtained by formally substituting (36) into (34) and imposing stationarity with respect to κ

$$ \psi(\kappa ) = \widehat{J^{KS}_{1}}(\boldsymbol{\xi}) + {\sum}^{m}_{e=1} \left[ \frac{(U_{e} - \xi_{e})^{2} {p^{1}_{e}}(\boldsymbol{\xi}) }{U_{e}-x_{e}(\kappa)} - \frac{(\xi_{e} - L_{e})^{2} {q^{1}_{e}}(\boldsymbol{\xi})}{x_{e}(\kappa)-L_{e}} \right] \!= 0 $$

(37)

The re-design procedure is implemented in the routine ocUpdate, listed at the bottom of this section. The input parameters are the current re-design step (loop), the current design point, objective sensitivity, constraint value and sensitivity (xT, dg0, g1, dg1), the design variables at the two previous iterations (xOld, xOld1), and the current asymptotes (as). The input ocPar=[move,asReduce,asRelax] collects the steplenght (move) and the two numbers used for tightening and relaxing the asymptotes according to the smoothness of the optimization history (asReduce,asRelax).

The evolution of the asymptotes and their link with the adaptive bounds δ₋ and δ₊ (called xL and xU) are conceptually identical to that suggested by Svanberg (1987) (see lines 4–13). Following Guest et al. (2011) we allow the tightening of the initial asymptotes to accomodate large β values, by passing this parameter to the ocUpdate routine. The user can also reset the asymptotes whenever continuation is applied on the penalization and projection parameters, by passing the additional variable restartAsy.

The primal update (36) and stationarity condition (37) are defined on lines 20–23, based on the positive variables p0, q0, p1 and q1, representing $(U_{e}-\xi _{e})^{2}{p^{i}_{e}}$ and $-(\xi _{e}-L_{e})^{2}{q^{i}_{e}}$, i = 0, 1, respectively (lines 17–18). The dual variable κ is computed between lines 26 and 33 accounting for three possible situations: given the search window $[0, \bar {\kappa }]$ (e.g., $\bar {\kappa } = 10^{6}$)

• 1. If $\psi (0)\psi (\bar {\kappa }) < 0 \rightarrow$ the Lagrange multiplier is within $(0, \bar {\kappa })$, and we compute it by the built-in Matlab function fzero, applying a version of the Brent’s algorithm (Brent 1973) (lines 27–28);

• 2. $\psi (0)<0 \rightarrow$ the constraint is not active within the current local expansion. Thus, we set κ = 0 and x_e = x_e(0) (line 30);

• 3. $\psi (\bar {\kappa })>0 \rightarrow$ the constraint cannot be fulfilled within the current local expansion (i.e., we cannot find a feasible solution). In this case, we set $\kappa = \bar {\kappa }$ and $x_{e} = x_{e}(\bar {\kappa })$ (line 32);

We stress that the present update rule can still be interpreted as an “OC-like” scheme since, for the simple case of a single constraint we are concerned with, the primal map (36) can be written out explicitly.

In our testing we have generally observed good behavior of this update rule. However, we acknowledge that the rather heuristic update for the special cases 2 and 3 may lead to oscillations in the convergence history, or even breakdowns. Since a robust optimizer is beyond the reach of a 35-line optimization routine and beyond the purpose of this paper, the users who might face bad behavior of the ocUpdate are recommended to replace it with the MMA. As an example, considering the classic Matlab version of the MMA “mmasub” (Svanberg 2007), the buckling maximization problem, with compliance and volume constraint, can be solved by replacing lines 224–226 with the following

if the user still wants to consider a single, KS aggregated constraint. Alternatively, the following call

can be used for treating the two constraints separately. The parameters a0MMA,aMMA,ccMMA and ddMMA can be selected as recommended in Svanberg (2007).

1.1 A2.1 Some properties of the KS function

Given the set of functions {g_i(x)}_i= 1…,q, where each $g_{i}\colon \mathbb {R}^{m}\rightarrow \mathbb {R}$ is not necessarily smooth, we can build the smooth Kreisselmeier-Steinhauser (KS) aggregation function (Kreisselmeier and Steinhauser 1979)

$$ J^{KS}[g_{i}](\mathbf{x}) = g_{\ast}(\mathbf{x}) + \frac{1}{\rho} \ln\left( \sum\limits^{q}_{i=1} e^{\rho\left( g_{i}(\mathbf{x}) - g_{\ast}(\mathbf{x}) \right)} \right) $$

(38)

depending on the parameter $\rho \in [1, \infty )$ and where $g_{\ast }(\mathbf {x}) = \max \limits _{i=1,\ldots ,q}\{g_{i}(\mathbf {x})\}$. Equation (38) is a convex function if and only if all its arguments g_i, i = 1,…,q are convex and fulfill $g_{\ast }(\mathbf {x}) \leq J^{KS}[g_{i}](\mathbf {x}) \leq g_{\ast }(\mathbf {x}) + \ln (\rho ^{-1} q_{a})$, where the right inequality is tight at points where q_a functions simultaneously attain the maximum value (Wrenn 1989; Raspanti et al. 2000).

With the aggregation of several constraint functions in mind, the following considerations are in order. At a given point x, if all g_i(x) > 0, then g_∗(x) > 0 and (g_i − g_∗)(x) ≤ 0 for all i. If all g_i(x) < 0, then g_∗(x) < 0 and (g_i − g_∗)(x) ≤ 0 for all i. In both cases J^KS[g_i](x) ≥ g_∗(x) and, if all the terms are negative this means that J_KS[g_i](x) is closer to zero than all g_i(x). Thus, it gives an upper bound to the constraint that is closer to become active. The interesting case is when some g_i(x) > 0 and some other g_j(x) < 0. Obviously, g_∗(x) > 0 and J^KS[g_i](x) gives an upper bound to the positive maximum; in other words, to the most violated constraint.

Appendix: 3. Matlab code