Efficient buckling constrained topology optimization using reduced order modeling

We present an efficient computational approach to continuum topology optimization with linearized buckling constraints, using Reduced Order Models (ROM). A reanalysis technique is employed to generate basis vectors that subsequently are used to significantly reduce the size of the generalized eigenvalue problems. We demonstrate the efficacy of this approach by optimizing for stiffness with buckling constraints and show results for several test cases. Based on our findings, we conclude that the ROM can potentially save significant computational effort without compromising the quality of the results.


Introduction
Topology Optimization (TO) is a powerful computational tool for designing structural components. The archetype TO design problem is compliance minimization under a volume constraint, which favors slender designs in states of pure compression and tension. Unfortunately, these structures are inherently prone to buckle, which renders them impractical under realistic loading conditions. A common way of accommodating for stability is to constrain the critical buckling load factor (BLF) such that safety against buckling exists (Neves et al. 1995, Gao and Ma 2015, Ferrari and Sigmund 2019, Dalklint et al. 2021. The TO problem is solved iteratively using a nested approach, wherein the implicitly enforced equilibrium equations are solved exactly in each design iteration. In practice, several hundred design iterations are required to find a minimum, and the linear system solves constitute the majority of the computational cost. One approach to reduce this cost is using Reduced Order Models (ROM), wherein the equilibrium equations are solved in a reduced space much smaller than the original high-fidelity space. The efficiency of a ROM is ultimately determined by the cost of generating the reduced space versus the loss of accuracy. One ROM approach that has been integrated into TO is to base the reduced space on previous solutions to the equilibrium equations, and update the space as new solution vectors become available (Gogu 2015, Choi et al. 2019).
An attractive ROM technique for structural optimization problems is reanalysis, in which the system matrix's factorization is stored and reused to build a reduced space for subsequent design iterations. Kirsch (1991) introduced a particular type of reanalysis known as Combined Approximation (CA), initially intended for approximating linear static displacements. This method has been showed to be equivalent to a preconditioned conjugate gradient method for linear systems of equations (Kirsch et al. 2002), and is known for its robustness and efficiency. The scope of CA was extended to include nonlinear, dynamic and eigenvalue problems (Kirsch 2008). In the context of integrating CA into TO, Amir et al. (2009) introduced reanalysis for topology optimization, and solved compliance minimization problems. Senne et al. (2022) extended reanalysis to TO problems with geometric and material non-linearities. Bogomolny (2010) used reanalysis for vibration problems.
When buckling is included in the optimization, the computational burden increases drastically, because the BLFs are obtained from a generalized eigenvalue problem. Eigenvalue problems are notoriously expensive to solve and scale poorly with the problem size. This cost becomes prohibitive for medium to large problems, and motivates the search for efficient numerical methods.
Recently, efforts have been made to extend TO with buckling constraints to large scale problems. Dunning et al. (2016) used a Block-Jacobi CG solver with a shift-andinvert strategy based on information from previous design iterations. Ferrari and Sigmund (2020) propose a multigrid approach wherein the multigrid preconditioner from the initial equilibrium problem is reused to estimate the BLFs and to filter out nonphysical, local modes. This reduced the cost for solving the eigenvalue problems significantly, while retaining sufficient accuracy in the eigenvalues. Also worth noting is the work by Bian and Fang (2017), where an assembly free method was used to solve the large-scale eigenvalue problem, and the work by Kang et al. (2020), where the vibration modes were estimated concurrently with the progress of optimization, improving the accuracy sequentially.
To the authors knowledge, reanalysis has not been applied to buckling constrained topology optimization. In this paper, we show that ROM using CA can be effectively integrated into TO with buckling constraints.

Preliminaries
We consider the finite element discretized linear system in combination with appropriate boundary conditions to describe our mechanical system. In (1), L ∈ ℝ n×n is the symmetric and positive definite linear stiffness matrix, o ∈ ℝ n is the reference displacement vector and o ∈ ℝ n is the load vector, where n is the number of degrees-of-freedom. To solve (1), we utilize the Cholesky decomposition of L in combination with forward/backward substitutions, i.e.
where ∈ ℝ n×n is upper triangular. After solving (1) for o using (2), we compute the symmetric and indefinite stress stiffness matrix G ( o ) ∈ ℝ n×n , which enables the generalized eigenvalue problem to be solved for the eigenpairs ( j , j ) ∈ (ℝ × ℝ n ) 1 . In (3), we assume T i L j = ij and a descending sorting of the eigenvalues. Solving the generalized eigenvalue problem in (3) for the BLFs j = 1 j , is referred to as linearized buckling analysis. The smallest BLF, i.e. largest j , defines the critical buckling load, j = 1 j o , at which a small perturbation results in that the structure buckles in the shape of j .

Reduced order modeling
We seek an approximation to the displacement field in (1) using a ROM, where = 1 , 2 , … s ∈ ℝ n×s contains the basis vectors and = [y 1 , y 2 , … y s ] ∈ ℝ s the reduced model coefficients. The coefficients are found by replacing o in (1) with ̃ o defined in (4), and projecting onto the basis i by premultiplying with T . The success of the ROM relies on the choice of the basis vectors, that is . In this work the basis vectors are generated using CA, wherein the factorization of the stiffness matrix is systematically stored and reused for subsequent design iterations. Quantities at the current design iteration and past design iterations will be distinguished using an overbar, for example L o = o and L o = o . Overbar quantities are understood to belong to the same, past, design iteration and are constants with respect to the current design. Approximate quantities are distinguished from their exact counterpart by • , for example ̃ o ≈ o denotes an approximate solution to (1).

Reanalysis of the linear equilibrium problem
Let us start by considering the approximate ̃ o solution to (1) and let L be the stiffness matrix from a past design iteration. In the CA approach, we utilize the Cholesky factorization of L , to compute o i.e.
where Δ L ≡ L − L is the stiffness change due to a design update. Applying −1 L to both sides of (5) yields 1 We pose the eigenvalue problem like (3) since G is indefinite, but emphasize that the original eigenvalue problem is where j are the BLFs. Now, we expand the inverse of ( + −1 L Δ L ) in a power series 2 Hence, the solution to (6) is which we truncate at k = s ≪ n , such that ̃ o = , cf. (4). The basis vectors in are obtained recursively from (8) as After the basis vectors are generated, we solve the reduced s × s problem for , and insert the result into (4). In practice, the magnitudes of i rapidly decrease due to repeated application of −1 L Δ , resulting in a poorly conditioned system (10). To mitigate this issue we introduce normalization and L -orthogonalization of the basis vectors wherefore the final basis vectors = 1 , 2 , … s are obtained recursively as Due to the orthogonalization, T L = , the reduced system in (10) boils down to = = T o , such that the approximate solution becomes An alternative, numerically equivalent, approach to (12) is to use the factorization of L as a preconditioner in a Preconditioned Conjugate Gradients (PCG) procedure, cf. e.g. Kirsch et al. (2002) and Amir et al. (2012). However, in this work the eigenproblem (3) is solved with ROM so for consistency we use the ROM approach for the static problem (1) as well.

Reanalysis of the buckling problem
Next, we seek an approximation of each eigenpair to (3). We do this by approximating the eigenmodes, after which the eigenvalues are given by the Rayleigh quotient. The eigenmodes are approximated as where j and j are analogous to the quantities defined in (4) and j ∈ ℕ n . To generate the basis vectors we utilize the additive decomposition of L in (3), i.e.
where G = G (̃ o ) . Premultiplication of (14) by −1 L yields and, using the result (7), we find from (15) The basis vectors j for mode j are obtained by truncating the infinite sum at term k = s ≪ n , replacing j on the righthand side with j and ignoring the arbitrary scaling factor of −( k ) −1 , i.e.
We scale the basis vectors to ensure that their relative magnitudes do not differ significantly, analogously to the static problem. Additionally we employ orthogonalization of the basis vectors with respect to lower order modes ̃ l by introducing ij = ij + ∑ j−1 l=1 l̃ l and requiring ̃ T l L ij = 0 for all i ≤ s, l < j . Thereby ensuring that the orthogonality condition of the eigenmodes is retained, that is ̃ T l L̃ j = lj . Due to the orthogonalization, we also avoid that the ̃j approximates the lower order l , l < j [compare Kirsch (2008) and Bogomolny (2010)]. The old eigenmodes, l , may be chosen instead of the current approximation, allowing for independent solution of the eigenproblems. However, in this work the latter approach is taken. The basis generation can be defined recursively as 2 This inverse exists as the (7) power series if

Renardy and Rogers (2006).
where j is the jth eigenvector corresponding to the matrix L . Now, inserting the (13) approximation into (3) and premultiplying by j = 1j , 2j , … sj gives a reduced generalized eigenvalue problem The system (19) has s solutions, and we seek the one which most accurately approximates the jth eigenpair ( j , j ) . Since the eigenpairs are solved in order of descending magnitude and j are L -orthogonal to ̃ l l < j , we choose eigenpair with largest magnitude eigenvalue. To summarize, for each sought eigenpair, we first generate a set of basis vectors using (18), and then solve (19) for the largest magnitude (̃j,̃ j ).

Design representation, filtration and thresholding
The goal of our topology optimization is to optimally distribute a linear elastic material in the design domain which is quantified by the piece-wise constant non-dimensional volume fraction ∈ ℝ n with 0 ≤ i ≤ 1 , where n is the number of elements. A well-posed optimization problem is obtained using restriction, via the Helmholtz PDE-filter, cf. Lazarov and Sigmund (2011). In this way, the filtered densities are obtained from where and appear in e.g. Lazarov and Sigmund (2011);Wallin et al. (2020) and T takes the nodal densities −1 ∈ ℝ n to the averaged, element-based densities ∈ ℝ n with 0 ≤ ν i ≤ 1 . To limit regions wherein ν i ∈ (0, 1) , we use thresholding [cf. Guest et al. (2004); Wang et al. (2014)] and penalization (Bendsøe (1989)), such that In (21), and are scalars defined such that lim →∞ H , ( ) = u s ( − ) , where u s is the unit step function.
Increasing q ∈ ℝ + enforces increasing levels of penalization, and o ∈ ℝ + is the ersatz material stiffness scaling. Also, e L and e G denote element stiffness matrices of L and G , respectively. Note that we use different density interpolations for G and L to limit the occurrence of spurious eigenmodes in low density regions (cf. Bendsøe MP and Sigmund (2003); Dalklint et al. (2020)).

Aggregation
In this work a lower bound constraint is enforced on the lowest magnitude BLF, that is the largest magnitude (or critical) i , which we denote c . To address the non-smoothness due to coinciding eigenvalues, we approximate c using aggregation in the form of a p-norm, cf. Torii and De Faria (2017). Therefore, the constraint on the critical BLF is stated as where ⋆ ∈ ℝ is the upper bound, p ∈ ℝ and we consider n ∈ ℕ n BLFs.

Optimization formulation
The topology optimization problem we solve is which is a standard compliance minimization problem subject to a lower bound constraint on the critical buckling load factor with maximum allowed volume fraction ξ . The (23) problem is solved in a nested approach wherein equilibrium equations (1) and (3) are enforced implicitly.

Sensitivity analysis
When computing the sensitivity of a non self-adjoint displacement dependent function, the implicit sensitivity of the displacements with respect to the design variables appears. In this work, we annihilate this implicit sensitivity through the adjoint method. Without going into details, the sensitivities of the BLFs are (Ferrari and Sigmund 2019) where it is assumed that (1) and (3) are fulfilled. An important consequence of (24) is that for each BLF we must compute an adjoint j which requires solving a linear system with L . The j can be estimated the same way o is estimated from the ROM by generating a new set of basis vectors using (11) and solving (12), where in both cases o is replaced by j . In this work, since we approximate o and j , we must instead use when estimating the adjoints.

Consistent sensitivity analysis
Since the sensitivity (24) assumes that (1) and (3) are solved exactly, it is not consistent when they are solved approximately. Amir et al. (2009) and Bogomolny (2010) account for the inconsistencies by replacing j with ̃j and introducing additional adjoint vectors corresponding to the steps in (18). This consistent sensitivity approach is computationally expensive and has shown to not offer any benefit over estimating (24) with ̃ o and (̃j,̃ j ) . Indeed, if the approximations are accurate the sensitivities are sufficiently accurate nonetheless. In this work the latter approach will be taken, and j in (24) will be estimated from the ROM (12). (24)

Numerical examples
To investigate the performance of our ROM, we study the optimization problem (23) for the three geometries depicted in Figs. 1a, 3a and 5a. The design domains are discretized using 180 × 420 quadrilateral, plane stress, bi-linear finite elements.
We first solve (23) without enforcing buckling constraints, that is we solve a stiffness optimization problem, which yields the critical BLF o . In our subsequent problems we let ⋆ = o ∕s where s is the safety factor, set to 3 unless stated otherwise. The stiffness optimized designs are used as initial guesses when solving the full (23) optimization problem, see Figs. 1b, 3b and 5b.
The design updates are computed using MMA with default parameters, see Svanberg (1987). Young's modulus and Poisson's ratio are 2 × 10 5 and 0.3, and the ersatz material stiffness scaling is 0 = 10 −6 . The material penalization is initially q = 2 , cf. (21), but is after iteration 30 incremented by 0.25 every 20 iterations until q = 6 . Initially we set p = 8 in (22), but after iteration 30 p is incremented by 0.50 every 20 iterations until p = 16 , and n = 6 unless stated otherwise. The filter length-scale parameter is l o = 3 2 √ 3 . We set the surface filter length-scale to l s = l 0 , everywhere except at nodes subjected to boundary conditions, where l s = 0 . Finally, = 6 and = 0.5 . All parameters are dimensionless. We stop the iteration when ‖ ‖ i − i+1 ‖ ‖2 ≤ 0.01 , or if the number of iterations reaches 400.
We distinguish between two different approaches when solving for equilibrium and adjoints. In the first (reference) approach we start by computing the Choleksky factorization of L , after which we solve (1) for o , (3) for the ( j , j ) and (24) for the adjoint j . The reference approach is summarized in Algorithm 1.

Algorithm 1 Reference approach
Solve for initial displacements 1: Solve for adjoints 4: for i ← 1, nµ do 5: In the second (approximate) approach we first determine if the factorization should be updated. If so, we follow the reference approach, otherwise we solve (12) for ̃ o , (19) for the (̃j,̃ j ) , and (12) for the ̃ j . A new factorization of L is performed if 1) p or q are updated or 2) the difference between the current physical densities i and the densities in the last factorization update is too large. We choose the cosine of the angle between the designs as a measure of their similarity, that is we update the factorization if where we use the threshold value of c min = 0.9996 . The approximate approach is summarized in Algorithm 2.

15:
end for 16: end if

Spire
The first geometry we study is shown in Fig. 1a. In Fig. 1b the stiffness maximized design is depicted. For this problem we use 8, 4 and 8 basis vectors for the linear equilibrium, eigenproblem and the adjoints, respectively. The designs obtained when constraining the critical BLF using the reference and approximate approach, shown in Fig. 1c, d, are very similar. This is also seen in Table 1 showing the value of the objective and constraints computed in a postprocessing step using the reference approach (see Algorithm 1). The constraints are satisfied for both designs, and the differences in the objective is minute. Table 1 also shows the number of matrix factorizations and triangular solves without and with ROM. We see that the number of factorizations is significantly reduced using ROM, while the difference in triangular solves is very small. Figure 2 shows the evolution of i ∕ ⋆ without and with ROM. We note that the evolution is similar for both approaches.
To verify that our ROM performs well for cases where a good initial guess is not available, the spire problem is solved starting from a homogeneous initial guess. The same number of basis vectors as in the previous example is used. The resulting designs are very similar to each other, but are slightly different from those in Fig. 1. The post-process analysis shows that they have the same structural qualities, and the reduction in computational effort is similar to the previous example: the number of factorizations was reduced from 352 to 62, and the number of triangular solves from 31344 to 28623.

V-structure
The second geometry we study is shown in Fig. 3a. For this problem we use 8, 4 and 8 basis vectors. Again, the designs obtained without and with ROM, seen in Fig. 3c, d, are similar. The objective and constraints at the end designs are provided in Table 2. The constraints are satisfied for both approaches, and the difference in objective is, again, very small. Table 2 also shows the number of matrix factorizations and triangular solves. We note that again the number of factorizations is significantly smaller and that the number of triangular solves is slightly smaller for the approximate approach. Figure 4 shows the evolution of i ∕ ⋆ , which is similar for the two approaches.

Column
The last geometry we study is shown in Fig. 5a. For this problem, the number of basis vectors for the static problem, eigenproblem, and adjoint was increased from 8, 4 and 8 to 10, 6 and 10, respectively. When constraining the smallest BLF, the designs shown in Fig. 5c, d are obtained without  and with ROM, respectively. For Fig. 5c design, the buckling constraint is satisfied and the p-norm approximating the smallest BLF renders undershooting. However, for Fig. 5d design, the constraint on the p-norm of the BLFs is slightly violated. Even so, due to the undershooting, the critical BLF still satisfies the upper bound constraint. The evolution of i ∕ ⋆ for Fig. 5c, d designs are shown in Fig. 6, and we note that, like in previous examples, the graphs are very similar.
Additionally, Table 3 shows the number of matrix factorizations and triangular solves without and with ROM. When using ROM the number of factorizations is reduced significantly, by over 75% , whereas the total number of triangular solves increased by about 37% . While this increase is noticeably larger compared to the previous examples, it still represents a small portion of the total computational effort which is dominated by the matrix factorizations.

Mode switching
As a final example, we demonstrate that our ROM can accurately approximate the BLFs when they switch order and coincide. For this reason, we solve the spire problem with a  higher safety factor of 6, starting from the compliance initial guess. We also increase the number of buckling factors to 10 to properly take into account the grouping of BLFs near the critical BLF. The number of basis vectors is increased to 10, 6 and 10, respectively.
The resulting designs, which can be found in Fig. 7, are again very similar. Table 4 shows the value of the objective and constraints, which are, in terms of stiffness and buckling resistance, almost identical. This example shows again that the ROM approach provides a significant speedup, in this case by 70% , and the number of triangular solves is increased only slightly. Figure 8 shows the evolution of the i ∕ ⋆ . We can clearly see that the ROM is able to resolve both coinciding BLFs and mode switching, see for example modes 6 and 7, and modes 8, 9, 10. Figure 9 shows the relative residuals of the eigenpairs, which we define by the quotient Of course, when using the reference approach the residuals are determined by the tolerance of the eigenvalue solver, so residuals are only shown at design iterations where the approximate approach is used.

Residuals
The relative residuals for the column problem are, towards the end of the optimization, in the range 10 −4 − 10 −2 and the residuals directly after a factorization update are higher compared to the other problems. Thus, to reduce the residuals and improve the quality of the approximation, the number of basis vectors could be increased further. For the spire problem the residuals are lower in general, but also have a larger spread which can be attributed to the spread in the BLFs themselves. For the V-structure problem, the residuals are comparatively lower, 10 −6 − 10 −3 , and are less scattered. Figure 9 shows that the geometry has an influence on the ROM's efficacy. While the method works well for all considered problems, some cases require more basis vectors  and/or factorizations to reach the same level of accuracy as others.

Discussion
The examples show that designs produced without and with ROM are very similar. In our implementation, we did not compute the consistent sensitivities that account for the inaccuracies of the ROM, nor did we compute the adjoints with full accuracy. Nevertheless, we still found designs very   Fig. 9 Relative residuals of the eigenvalue problem for a the spire problem, b the V-structure problem and c the column problem. Residuals are only shown at design iterations where the approximate approach is used similar to those generated by the reference approach. We demonstrate that our ROM can, without modifications, be applied to problems where a good initial guess is not available. We also show that the ROM can accurately resolve both mode switching and coinciding eigenmodes. Compared to the reference approach, the ROM approach reduced the number of factorizations significantly, and with exception to the column example, the number of triangular solves did not change significantly. Since the triangular solves only constitute a small portion of the total effort compared to the matrix factorizations, the computational effort is reduced significantly by the proposed ROM approach. When the required number of BLFs grows, the savings might decrease due to the increasing cost of approximating the adjoints.
The number of basis vectors and when to update the factorization was chosen heuristically. Even so, we were able to reduce the computational effort without fine-tuning the parameters. In some cases, the ROM was inaccurate, but adding more basis vectors improved the accuracy without significantly increasing the costs. Based on our numerical experiments, the optimal number of basis vectors depends on the problem and on the design. Ideally, the number of basis vectors and the frequency of factorization should be adaptive, such that the approximation is both efficient and sufficiently accurate.
The ROM-solver for the eigenvalue problem has several advantages over traditional eigenvalue solvers that iterate over a subspace: the number of basis vectors for each eigenpair can be chosen independently. Thus, the effort can be balanced based on how each eigenpair affects the problem-higher accuracy can be enforced for more important eigenpairs. Furthermore, the approximate solution of eigenpairs is embarrassingly parallelizable as the eigenpairs can be solved independently of each other, if the old eigenmodes are chosen in the orthogonalization step. This too should be considered in future work.
The proposed ROM requires a direct solve using the system matrix, limiting its use to problems which can be efficiently solved using direct solvers. Solving problems beyond this size requires other approaches, such as highly tailored multigrid preconditioned subspace iteration methods. For this we refer to the work by Peetz and Elbanna (2021), or the multigrid approach by Ferrari and Sigmund (2020).

Concluding remarks
In this paper, we have shown that a reanalysis-based ROM can be applied to topology optimization considering linearized buckling, and that it has clear potential to reduce the computational cost. Our examples indicate that the ROM works for a range of problems, and that more effort needs to be spent on developing adaptive schemes-we presume that further computational savings are possible if the number of basis vectors and the frequency of factorizations are adapted throughout the optimization process.