Alternating optimization of design and stress for stress-constrained topology optimization

Design of structures with local stresses upper-bounded by a critical stress value is of paramount importance in engineering. To this end, the incorporation of stress constraints has been an important field of study in topology optimization of continuum structures (Duysinx and Sigmund 1998). Over the past two decades, three computational challenges have been recognized (Le et al. 2010; Holmberg et al. 2013), and solutions for some of them have been proposed:

• The “singularity” problem — the feasible design space contains degenerate sub-spaces of a lower dimension (Kirsch 1990; Rozvany 2001). The globally optimal solution, which often locates in the degenerate sub-spaces, is not accessible to nonlinear programming algorithms. It has been shown that the problem of degenerate sub-spaces can be alleviated by relaxing the stress constraints, i.e., the 𝜖-relaxation originally developed for trusses (Cheng and Guo 1997) and its variants for continuum structures (Duysinx and Bendsøe 1998; Bruggi 2008; Le et al. 2010).

• The local nature of stress constraints — the stress limit applies to every material point in the domain. This results in a large number of constraints. An often used solution strategy is to approximate these local constraints by a global one which can then be more efficiently addressed, e.g., the p-norm (Duysinx and Sigmund 1998) and Kreisselmeier-Steinhauser (KS) function (Yang and Chen 1996).

• The highly nonlinear dependence of stress on design variables. Especially at stress concentration regions, stresses are sensitive to density changes in neighborhoods. This leads to convergence problems: a large number of iterations, fluctuations in the objective and constraints, and a suboptimal objective value which potentially could be further reduced.

This last challenge is coupled with solutions of the first two. For instance, the aggregation for reducing the number of constraints may further exacerbate the nonlinearity. It is found in the literature that research has been mostly focusing on the first two challenges by reformulating the optimization problem, and effective alternative approaches have been proposed, e.g., Verbart et al. (2016, 2017), Wang and Qian (2018).

Stress-constrained topology optimization problems (and more general structural optimization) are typically solved by using sequential convex programming. Notable algorithms include the Convex Linearization method (CONLIN) (Fleury 1989), the Method of Moving Asymptotes (MMA) (Svanberg 1987) and its globally convergent version GCMMA (Svanberg 2002). These algorithms make use of (first order) approximations of objective and constraint functions. An assessment of these optimization algorithms, as well as general methods such as the primal-dual interior point methods (Forsgren and Gill 1998) and Sequential Quadratic Programming (SQP) (Boggs and Tolle 1995), is presented by Rojas-Labanda and Stolpe (2015). The benchmark problems in the comparative study include compliance/volume minimization and mechanism design. Unfortunately stress constraints are not included.

An alternative solution strategy to stress-constrained topology optimization is to incorporate local stress constraints in the objective function using an augmented Lagrangian formulation (Pereira et al. 2004). It has been used for topology optimization based on density (Fancello 2006; da Silva and Cardoso 2017; da Silva et al. 2019) and level sets (James et al. 2012; Emmendoerfer and Fancello 2014). Very recently (da Silva et al. 2020) demonstrated that this formulation allows handling very large problems in 3D manufacturing-tolerant topology optimization, with hundreds of millions of stress constraints. Also aiming for 3D large scale optimization, Senhora et al. (2020) proposed to modify both the penalty and objective function terms of the augmented Lagrangian function, leading to consistent solutions under mesh refinement and driving the mass minimization towards black and white solutions. Giraldo-Londoño and Paulino (2020) applied this to handle multiple classical failure criteria with a unified yielding function.

In this paper we introduce an interpretation of local stresses as optimization variables, using an augmented Lagrangian formulation. We consider auxiliary stresses as optimization variables, in addition to the design variables (i.e., densities) representing the material distribution. The stress limit is then imposed upon the auxiliary stresses as an upper bound. Given a material distribution and boundary conditions, the actual stress is computed by a finite element analysis. The auxiliary stress is related to the actual stress by an equality constraint. This equality constraint is then incorporated into the objective function as linear and quadratic terms using an augmented Lagrangian form.

This reformulation offers some conceptual benefits. Firstly, thanks to the extra set of optimization variables, it opens a larger space for the optimization to search. It even allows the optimization to explore the infeasible region in the conventional formulation — the actual stress could be larger than the prescribed stress limit. This can be beneficial for rapid convergence at the beginning iterations. Secondly, the stress limit is regarded as side or bound constraints for the auxiliary stress variables, and bound constraints are readily handled by optimization algorithms.

To efficiently solve the optimization problem with two sets of variables, we show that the formulation is separable regarding the auxiliary stresses and structural variables (i.e., densities). This gives rise to a variant of augmented Lagrangian solvers known as Alternating Direction Method of Multipliers (ADMM). Within each iteration, the structural variables, stress variables, and Lagrangian multipliers are alternatingly updated. This solver converges rapidly at the first few iterations. Our tests on continuum structures demonstrate reduced objective values in comparison with using the MMA family for solving the conventional formulation.

ADMM algorithms have provable convergence for solving convex optimization problems (Boyd et al. 2011). Their applicability to nonconvex problems has also been investigated, e.g., Diamond et al. (2018). For nonconvex problems, while ADMM is not guaranteed to converge, in practice it can often find a reasonable objective value with small computational cost (Kanno and Kitayama 2018). It has been recently applied in structural optimization for problems with a few or hundreds of design variables. Kanno and Kitayama (2018) used ADMM as an effective heuristic for mixed-integer nonlinear structural optimization. Palanduz and Groenwold (2020) analyzed the applicability of a subset of ADMM-type algorithms for optimal structural design, in combination with a novel scaling method and quadratic approximations of the primal problem. Eckstein and Bertsekas (1992) pointed out that ADMM is an application of the Douglas-Rachford splitting algorithm (Douglas and Rachford 1956), which is in turn an instance of the proximal point algorithm. Recently, Evgrafov and Sigmund (2020) proposed to split stress variables into positively and negatively semi-definite parts, and applied this splitting for compliance minimization in the vanishing volume ratio limit. In this paper we adopt ADMM to solve stress-constrained problems with a large number of variables and constraints, and compare its performance with conventional formulations and solvers.

Our main contributions in this work include the following: (i) an interpretation of local stresses as optimization variables, and (ii) an extensive comparison of different solution strategies for stress-constrained topology optimization. This intuitive interpretation leads to equations that turn out to be mathematically equivalent to those from the augmented Lagrangian formulation (e.g., Giraldo-Londoño and Paulino 2020). Using simple truss examples, we validate that the augmented Lagrangian formulation allows reaching the global optimum that is located in the degenerate sub-space, similarly to a unified aggregation and relaxation approach (Verbart et al. 2017). Tested on continuum structures, our comparison reveals that the augmented Lagrangian formulation can achieve superior solutions in terms of the objective value. Furthermore, we propose a hybrid approach which based on our comparison is most favorable regarding the objective value and computational efficiency.

We start with a conventional density-based formulation for stress-constrained topology optimization (Section 2.1). We then introduce an interpretation of stresses as optimization variables (Section 2.2). An optimization strategy to this reformulation is explained in Section 2.3.

2.1 Stress-constrained compliance minimization

Using a finite element discretization of the design domain, the density distribution is represented by a vector, ρ. Let σ_e denote the stress tensor for generic element e when the structure is under a constant external load, F. We introduce a variable, α_e, to indicate a scalar measure of the stress tensor, i.e.,

The stress-constrained compliance minimization is then written as

This conventional formulation is denoted as \(\mathcal {Q}_{0}\) to differentiate from a reformulation to be introduced soon. In the objective, U(ρ) represents the displacement vector. In (3), vector v is composed of the volume per element, and γ is a prescribed volume ratio. The stress constraint is encoded by (4), where \(\sigma _{\lim }\) is a prescribed limit on the scalar stress measure. At last, (5) restricts the pseudo density of each element to between 0 (empty) and 1 (solid).

The displacement vector, U(ρ), is computed by solving a static equilibrium equation,

where K(ρ) is the assembled stiffness matrix. The element stiffness matrix is calculated using the Solid Isotropic Material with Penalization (SIMP) model,

where K₀ is the stiffness matrix of a solid element with Young’s modulus E₀. E_e is the interpolated Young’s modulus. E_min is a small value to prevent the assembled stiffness matrix (K) from becoming singular. k is a penalization parameter (k = 3).

2.2 Interpreting stresses as optimization variables

An innovative aspect of our approach is to interpret stress α_e as an auxiliary optimization variable, and consequently (1) as an equality constraint rather than a definition. To efficiently handle the equality constraint, the objective function is updated by incorporating an augmented Lagrangian form of the equality constraint,

where λ is a Lagrange multiplier and μ is a penalty parameter. The conventional optimization problem is reformulated as

Intuitively, interpreting stresses as optimization variables enlarges the optimization space, giving extra flexibility for searching an optimal solution. It relaxes the stress constraint in the sense that it permits, at the beginning of the optimization process, the actual stress (i.e., \(\tilde {\sigma }_{e}\), computed from ρ) being larger than the upper bound. This is beneficial for a rapid decrease of the objective at the onset of the optimization process.

The augmented Lagrangian formulation can be solved by simultaneously optimizing the two sets of variables, by using for instance, MMA. To efficiently solve this formulation, in the following we introduce a method to decouple the two sets of optimization variables.

2.3 Alternating optimization

While two sets of optimization variables are employed, the objective, f_c(ρ), and inequality constraints, i.e., (3), (4) and (5), involve either ρ or α, but not both. In other words, each of these functions can be separated into two parts: one part that is dependent on ρ, and the other part that is dependent on α. In each of these functions, one of these two parts is null.

This so-called separability is essential for a variant of augmented Lagrangian methods, known as the alternating direction method of multipliers or ADMM. ADMM uses partial updates for the variables, similar to the Gauss-Seidel method for solving linear equations. This partial update scheme allows an efficient search of the optimal solution in the enlarged optimization space.

By making use of the separability, the density and stress variables are updated alternatingly. ADMM consists of the iterations

This is followed by updating the Lagrange multiplier,

Each of the two sub-problems itself is an optimization problem. In initialization, we set the density equal to the prescribed volume ratio, i.e., \(\rho ^{[0]}_{e}=\gamma , \forall e\), and compute \( {\alpha }^{[0]} = \tilde { {\sigma }}( {\rho }^{[0]})\). The initial value of Lagrange multipliers is \({\lambda }^{[0]}_{e} = 1, \forall e\). In the first sub-problem, the density variables are updated, while the stress variables and Lagrange multipliers are fixed. We solve this sub-problem by using MMA. Alternatingly, in the second sub-problem, the stress variables are optimized with updated densities (ρ^[i]) and Lagrange multipliers from the last iteration (λ^[i− 1]). In fact, (11) is a quadratic programming problem with respect to α, and can be reformulated as

Considering the bound constraint, the solution of the second sub-problem is

Before we proceed to the implementation details of continuum structure optimization, let us examine the validity of the proposed approach on a simple two-bar truss optimization problem. The example is illustrated in Fig. 1 (Stolpe 2003; Verbart et al. 2017). The problem is to minimize its mass subjected to a stress limit, \(\sigma _{\lim }\). The design variables are the cross-sectional areas A = [A₁,A₂]^T. The cross-sectional area is upper bounded by \(A_{\max \limits }=2\). Both members have a Young’s modulus E. ρ_j and L_j denote the density and length of member j, respectively. The stress in the members is calculated by

A two-bar truss example (Stolpe 2003) is used to validate the proposed approach. The optimization problem is to minimize the total mass by varying the cross-sectional areas, A₁ and A₂, under stress limit \(\sigma _{\lim }\)

Full size image

The initial problem formulation, denoted by \(\mathcal {P}_{0}\), is written as,

Figure 2a illustrates the solution space of the problem. It is composed of the polygon region BCDEF and the line segment EG, with the global optimum located at G. The global optimum is located in a 1-dimensional sub-space within the 2-dimensional solution space. Without special treatments, the optimization reaches a local minimum in the polygon region BCDEF.

Solution space of the two-bar truss problem in Fig. 1, corresponding the initial formulation \(\mathcal {P}_{0}\) (a) and the p-mean aggregated formulation \(\mathcal {P}_{0}^{M}\) (b)

Full size image

The unified approach proposed by Verbart et al. (2017) is included here as a reference for addressing the singularity problem. The optimization problem is modified as

Here, G^L denotes the aggregated constraint, based on the p-mean function

with

Illustrated in Fig. 2b, the p-mean aggregated formulation, \(\mathcal {P}_{0}^{M}\), enlarges the solution space to additionally include the shaded region \(E\widetilde {FG}\). Curve \(\widetilde {FG}\) corresponds to the boundary obtained when p = 8 in the p-mean function (Verbart et al. 2017).

3.1 Augmented Lagrangian formulation and alternating optimization

Following Section 2.2, we introduce an auxiliary variable α_j, and relate it to the stress by α_j = σ~_j, with \(\Tilde {\sigma }_{j} = \frac {A_{j}}{A_{\max \limits }}\frac {\sigma _{j}}{\sigma _{\lim }}\)^{Footnote 1}. The equality constraint α_j = σ~_j is then incorporated in the objective using an augmented Lagrangian form,

The optimization problem is reformulated as

ADMM consists of alternatingly solving two sub-problems, i.e.,

s.t. (26).

s.t. (27).

And the Lagrange multiplier is updated by

In our test, we initialize λ^[0] = (1,1)^T and μ = 0.5. μ is multiplied by 1.05 after every 5 iterations.

3.2 Results and discussion

We compare six solution strategies for solving three formulations, \(\mathcal {P}_{0}\), \(\mathcal {P}_{0}^{M}\), and \(\mathcal {P}_{1}\). The results are shown in Fig. 3. In all cases, the initial values of A₁ and A₂ equal to \(A_{\max \limits }=2\), corresponding to C in the solution space. The global optimum is located at G(1,0). Using MMA to solve \(\mathcal {P}_{0}\) leads to a local minimum point F(0,1) in the polygon region (Fig. 3a). With the p-mean aggregated formulation \(\mathcal {P}_{0}^{M}\), both MMA and GCMMA are able to converge to the global optimum (Fig. 3b and c). Figure 3d and e show results of solving the augmented Lagrangian formulation \(\mathcal {P}_{1}\) using MMA and ADMM, respectively. Lastly, Fig. 3f shows the results of a hybrid approach: we use ADMM to solve \(\mathcal {P}_{1}\), and, after five iterations, switch to MMA for solving \(\mathcal {P}_{0}^{M}\).

Comparison of convergence histories. The optimization starts at point C(2,2). Except in (a) where MMA is used to solve \(\mathcal {P}_{0}\), in other cases the global optimum, G(1,0), is reached

Full size image

From Fig. 3, it can be observed that the global optimum (G) is correctly identified using the aggregated constraint formulation \(\mathcal {P}_{0}^{M}\) and the augmented Lagrangian formulation \(\mathcal {P}_{1}\). This is also confirmed in Table 1. The minimal mass is 0.60. In the plots of the augmented Lagrangian formulation (Fig. 3d and e), the path leading to the global optimum travels beyond the feasible region. This is not unexpected since the equality constraint \( {\alpha }= {\tilde {\sigma }}\) is transformed to penalty terms in the objective, and thus the stress constraint is not strictly enforced during the process. Yet, as the process converges, the equality condition and consequently the stress constraint are satisfied.

Discussion

The augmented Lagrangian formulation reduces constraints into a single scalar value. On this regard, it plays the role of constraints aggregation. The augmented Lagrangian formulation is also able to achieve what constraint aggregation schemes can offer, i.e., allowing an expansion of the search space and thus making the global optimum accessible (Verbart et al. 2017). We refer to Verbart et al. (2017) for a clever use of constraint aggregation schemes for stress-constrained optimization.

The convergence of the total mass is plotted in Fig. 4. ADMM is capable of rapidly reducing the objective at the first few iterations. As it approaches the final solution, it progresses slowly and leads to more iterations. This behavior is known and has been discussed, for instance, in Boyd et al. (2011). This motivates the hybrid approach which uses ADMM to quickly find a good start solution. The computation time is summarized in Table 1. It suggests that the hybrid approach, ADMM (\(\mathcal {P}_{1}\)) + MMA (\(\mathcal {P}_{0}\)), is promising in terms of the number of iterations and wall-clock time.

The convergence of the total mass using different solution strategies

Full size image

In the context of continuum structure optimization which is highly nonlinear and involves a large number of design variables, as will be shown in later sections, the augmented Lagrangian formulation demonstrates not only benefit in the computation time but more importantly in the optimality of the optimized solution.

4.1 Density operations

We adopt density operations that are commonly applied in density-based approaches, i.e., a smoothing operator to avoid checkerboard patterns (i.e., regions of alternating solid and empty elements), and a smoothed Heaviside operator for improving convergence towards black-and-white designs (Guest et al. 2004; Wang et al. 2011).

4.1.1 Density smoothing

Instead of taking ρ as the design variable, a new variable, ϕ ∈ [0,1], is introduced as the design variable. A smoothed version of ϕ is obtained by computing the weighted average,

where v_i is the volume of an element, and the weighting function is defined as

where r is a predefined filter radius, x_e and x_i are respectively the centroid of element e and an element (i) that is in close vicinity, i ∈ S_e = {i|∥x_i − x_e∥≤ r}.

4.1.2 Heaviside projection

The smoothed Heaviside projection is written as

where β controls the sharpness of the step function and η is a projection threshold (η = 0.5). To avoid numerical instability, we start with β = 1 and double its value every 50 iterations.

4.2 Stress relaxation

We follow the qp-relaxation scheme proposed by Le et al. (2010) to resolve the stress singularity problem. The scalar measure of per element stress is computed by

where q is a relaxation parameter (q = 0.5). \(\bar {\sigma }_{e}\) refers to the von Mises stress of element e, defined by

where V is a symmetric matrix

and σ_e is the Voigt notation of a stress tensor, calculated by

Here, D₀ is the elasticity tensor in the Voigt notation for a solid element, B_c is the strain-displacement matrix of the element centroid and u_e is the displacement vector of element e.

4.3 Sensitivity analysis

We present the derivatives of the augmented Lagrangian form, L_μ(ρ,α,λ), regarding density variables (ϕ_e). Since α and λ are not functions of the density variables, \(\frac {\partial L_{\mu }( {\rho }, {\alpha }, {\lambda })}{\partial \phi _{e}}\) is computed by

Detailed sensitivity analysis is presented in Appendix 1.

4.4 Algorithm

Algorithm 1 details the process of using ADMM to solve the augmented Lagrangian form. The algorithm takes the prescribed volume fraction (γ) and stress limit (\(\sigma _{\lim }\)) as input parameters. The output is the material distribution (ρ). We prescribe the maximum number of ADMM iterations I_ADMM = 1000 or I_ADMM = 20. The later is used in the hybrid solution which uses ADMM for \(\mathcal {Q}_{1}\), followed by using MMA for \(\mathcal {Q}_{0}\) (see Appendix 2 for solving \(\mathcal {Q}_{0}\)).

The proposed method is implemented using MATLAB 2018b, and tested on a desktop PC with an Inter(R) Xeon(R) W-2123 CPU @3.60GHz and 64GB RAM. The material model has a Young’s modulus of E = 1.0 and Poisson’s ratio of ν = 0.3. Bi-linear square elements are used for stress analysis. The solution strategies for the conventional and the augmented Lagrangian formulation for compliance minimization are quantitatively compared using ten examples. The same termination criteria are applied: the maximum change in design variables in successive iterations is smaller than a threshold, 𝜖 = 10^− 2, and that both the volume and stress constraints are satisfied, or the maximum number of iterations (1000) is reached. For each solution strategy, all examples use the same parameters, which are selected empirically based on multiple tests using the L-shaped beam. Important parameters in MMA include move, asyincr, and asydecr. Based on tests, move = 0.1, asyincr = 1.2 and asydecr = 0.7 are selected.

5.1 L-shaped beam

The first example is an L-shaped beam as depicted in Fig. 5 (left). The design domain is specified by DH = DW = 150, Dh = Dw = 90. The top edge of the domain is fixed while a distributed force (F = 1) is applied downwards on the right, over a short span of 4 elements horizontally. Figure 5 (right) visualizes the von Mises stress distribution of a compliance minimized structure which is optimized without stress constraints, under volume ratio γ = 0.3. The visualization indicates a stress concentration at the corner. The maximum stress is 1.35. In the test of stress-constrained optimization, we set a stress limit of \(\sigma _{\lim }=0.5\).

Design domain of the L-shaped beam (left) and visualization of the von Mises stress distribution of a structure optimized without stress constraints (right)

Full size image

The stress-constrained optimization problem is solved using two formulations, in total six algorithms. Fig. 6 illustrates the optimized structural layouts and the corresponding von Mises stress distributions. The first row is obtained using the conventional formulation with densities as design variables (\(\mathcal {Q}_{0}\)), while the second row corresponds to the proposed formulation with two sets of optimization variables (\(\mathcal {Q}_{1}\)). Note that in the last approach, ADMM (\(\mathcal {Q}_{1}\)) + MMA (\(\mathcal {Q}_{0}\)), we use MMA to solve the conventional formulation rather than the new formulation which has a doubled number of variables.

Optimized structural layouts and stress distributions on the L-shaped beam using different optimization methods

Full size image

Table 2 summarizes the statistics of using different optimization methods. The methods are grouped into two categories according to the formulation of the optimization problem, \(\mathcal {Q}_{0}\) and \(\mathcal {Q}_{1}\). The methods are comparable in terms of constraints on the total volume, V, and maximum stress, \(\max \limits (\tilde { {\sigma }})\). Of interest here are the objective value and computation time of these six solution strategies, which are visualized using bar graphs in Fig. 7. In the former category, GCMMA leads to a smaller compliance than MMA (296.61 vs 301.53), at the price of a longer computing time (1464.6 vs 759.6). GCMMA requires a smaller number of iterations than MMA. However, on average a GCMMA iteration is more costly than an MMA iteration. GCMMA involves an inner loop which requires solving the state equation. The last row (# Solves) indicates how many times the state (and adjoint) equation is solved. The hybrid approach of GCMMA and MMA achieves a comparable compliance (295.91).

Bar graph of the objective value and computation time for minimizing compliance of the L-shaped beam. The numbers 1 to 6 represent the different solution strategies in the order as they appear in Table 2, i.e., MMA (\(\mathcal {Q}_{0}\)), GCMMA (\(\mathcal {Q}_{0}\)), GCMMA (\(\mathcal {Q}_{0}\))+MMA (\(\mathcal {Q}_{0}\)), MMA (\(\mathcal {Q}_{1}\)), ADMM (\(\mathcal {Q}_{1}\)), ADMM (\(\mathcal {Q}_{1}\))+MMA (\(\mathcal {Q}_{0}\))

Full size image

The objectives from the category of \(\mathcal {Q}_{1}\) outperform those from the conventional formulation \(\mathcal {Q}_{0}\). ADMM (\(\mathcal {Q}_{1}\)) achieves a compliance of 281.99, which is 6.48% smaller than that by MMA (\(\mathcal {Q}_{0}\)). On the last column, the hybrid approach arrives at a slightly larger compliance, with a reduced total computation time.

Figure 8 plots the convergence of the volume and stress constraint (\(\max \limits (\tilde { {\sigma }})-\sigma _{\lim }\)) over iteration. It is found that ADMM is able to quickly meet the volume constraint, in less than 50 iterations. The stress constraints are gradually satisfied, since the penalty term increases during the optimization progress. MMA in general exhibits less severe fluctuations compared ADMM. This motivates the hybrid approach, ADMM (\(\mathcal {Q}_{1}\)) + MMA (\(\mathcal {Q}_{0}\)). The convergence in compliance for all six solutions is plotted together in Fig. 9. A few fluctuations in these plots are due to the continuation of the projection parameter (β) in (33).

Plots of convergence of the volume and stress constraint (\(\max \limits (\tilde { {\sigma }})-\sigma _{\lim }\)) for optimizing the L-shaped beam

Full size image

Plots of compliance convergence for optimizing the L-shaped beam

Full size image

5.2 Step-shaped structure

The second set of experiments is performed on a step-shaped structure shown in Fig. 10 (left). The design domain is descried by DH = DW = 150 and Dh = Dw = 50. A horizontally distributed force is applied on the top of 4 elements in the top right of the domain. Figure 10 (right) visualizes the von Mises stress distribution of a compliance-minimized structure optimized in the absence of stress constraints, with a target volume ratio of γ = 0.25. Three stress concentration regions can be observed. The maximum von Mises stress is 1.16. For this example, we test three stress limits: \(\sigma _{\lim }\in \{0.6, 0.7, 0.8\}\).

Design domain of a step-shaped structure (left) and visualization of the von Mises stress distribution of a structure optimized in the absence of stress constraints (right)

Full size image

Table 3 summarizes the statistics of tests on this example. In all cases, the constraints on the total volume and maximum stress are satisfied. From the top block to the bottom one, as the allowed maximum stress increases, the objective value decreases. This trend in compliance is better visualized from the bar graphs in Fig. 11. There is however no clear pattern in the computation time for different stress limits.

Bar graph of the objective value and computation time for minimizing compliance of the step-shaped structure. The numbers 1 to 6 represent the different solution strategies in the order as they appear in Table 3, i.e., MMA (\(\mathcal {Q}_{0}\)), GCMMA (\(\mathcal {Q}_{0}\)), GCMMA (\(\mathcal {Q}_{0}\))+MMA (\(\mathcal {Q}_{0}\)), MMA (\(\mathcal {Q}_{1}\)), ADMM (\(\mathcal {Q}_{1}\)), ADMM (\(\mathcal {Q}_{1}\))+MMA (\(\mathcal {Q}_{0}\)). The triple indicates the increasing stress limits, from left to right, 0.6, 0.7 and 0.8

Full size image

Comparing the final objectives from the six algorithms, it is observed that the augmented Lagrangian formulation in general results in smaller compliance values. ADMM is able to reduce the compliance value from MMA (\(\mathcal {Q}_{0}\)) by 8.78% (\(\sigma _{\lim }=0.6\)), 7.73% (\(\sigma _{\lim }=0.7\)) and 3.16% (\(\sigma _{\lim }=0.8\)), while for the hybrid strategy, ADMM (\(\mathcal {Q}_{1}\)) + MMA (\(\mathcal {Q}_{0}\)), the reduction is 7.06%, 6.44% and 4.57%, respectively. This reduction in objective is in agreement with results from the previous example. In these tests, the hybrid approach have comparable to slightly reduced computation time than that of MMA (\(\mathcal {Q}_{0}\)).

As a side note, we noticed that GCMMA was sensitive to the stress limit. As the stress limit became smaller, the constraints were not satisfied after 1000 iterations under a default setting. A parameter in GCMMA was thus modified for achieving satisfactory results. The optimization framework of GCMMA consists of “outer” and “inner” iterations. Within each outer iteration, depending on whether the approximation function is conservative or not, a number of inner iterations are carried out. The more inner iterations, the more conservative the approximation function. We found that an increase of the inner iterations in GCMMA from 5 to 20 was able to cope with the tight stress limit in our example.

Figure 12 shows the optimized structural layouts and stress distributions using different solvers, under the stress limit \(\sigma _{\lim }=0.6\). The overall layouts are similar. In all cases, smooth geometry is developed around the sharp corners in the design domain, alleviating stress concentration. Interestingly, comparing (d)–(f) with (a), we find that the vertical substructure near the stress concentration region on the bottom corner (indicated by an arrow in (a)) is closer to the domain boundary in (d)–(f) than in (a). Note that in the structure optimized without stress constraints, shown in Fig. 10 (right), the vertical substructure conforms to the domain boundary with no gap. Similar trend can be observed in the L-shaped designs (Fig. 6).

Optimized structural layouts and stress distributions in the step-shaped design domain using different optimization methods

Full size image

5.3 Parameters and discussion

5.3.1 More examples

In addition to the four tests above, six other examples are tested. The specifications and detailed statistics can be found in Appendix 3. To compare results of different examples, we take MMA (\(\mathcal {Q}_{0}\)) as a reference, and normalize the compliance and computation time from other solution strategies by that of the reference. The average and standard deviation of the ten examples are summarized in Fig. 13. From the left, it is observed that ADMM (\(\mathcal {Q}_{1}\)) on average attains the smallest compliance, followed by the hybrid approach, ADMM (\(\mathcal {Q}_{1}\)) + MMA (\(\mathcal {Q}_{0}\)). On the right, the comparison shows that the computation time of the hybrid approach is marginally smaller than that of MMA (\(\mathcal {Q}_{0}\)).

Comparison of the different solution strategies using the average and standard deviation from 10 examples. The numbers 1 to 6 represent the solution strategies in the order as they appear in other comparison figures and tables, i.e., MMA (\(\mathcal {Q}_{0}\)), GCMMA (\(\mathcal {Q}_{0}\)), GCMMA (\(\mathcal {Q}_{0}\))+MMA (\(\mathcal {Q}_{0}\)), MMA (\(\mathcal {Q}_{1}\)), ADMM (\(\mathcal {Q}_{1}\)), ADMM (\(\mathcal {Q}_{1}\))+MMA (\(\mathcal {Q}_{0}\))

Full size image

5.3.2 Lagrange multiplier λ and penalty parameter μ

Augmented Lagrangian formulation involves two parameters: the Lagrange multiplier λ and the penalty parameter μ in (9). These two parameters effectively control the discrepancy between α and \(\tilde { {\sigma }}\). To achieve a smooth convergence process, in general it is beneficial to start with small λ and μ, and thus in the beginning iterations the combined objective is dominated by the compliance. We choose λ_e = 1 and μ = 0.5, such that the compliance (f_c(ρ)) accounts for more than 99.9% of the objective (L_μ(ρ,α,λ)) after the very first iteration. To penalize the discrepancy, μ is multiplied by 1.05 after every 5 iterations.

The evolution of the Lagrange multiplier is depicted in Fig. 14 (right) for three iterations. The actual stress and stress design variable are visualized in the second and third column, respectively. The magnitude of the Lagrange multiplier is large in regions where the stress design variable significantly deviates from the actual stress. As the optimization converges, the Lagrange multiplier takes values close to zero (see (f)).

Visualization of the material distribution (ρ), actual stress (\(\tilde {\sigma }\)), stress design variable (α) and Lagrange multiplier λ at three different iterations

Full size image

5.3.3 Initial iteration number

In the ADMM (\(\mathcal {Q}_{1}\)) + MMA (\(\mathcal {Q}_{0}\)) approach, a fixed number of ADMM iterations is performed. This value is chosen based on analyzing the stability in the compliance history. We use the relative standard deviation (RSD) of the compliance value over a period of five iterations as an index for evaluating stability,

with \(\bar {f_{j}} = \frac {1}{5}{\sum }_{i=j-2}^{i=j+2}{f_{i}}\). From a number of tests we found that it took 12 to 25 iterations for the relative standard deviation to become less than 0.2%. Thus in our examples we chose to run 20 ADMM iterations.

In both examples the maximum stress is of the same magnitude to the range of density variables. We have performed additional tests by increasing the external forces by two orders of magnitude. Consequently the density variables and the stress variables are no longer of the same magnitude. In these tests, solving the new formulation by alternating optimization of design and stress (ADMM) is able to converge, while simultaneously solving for both sets of variables (MMA (\(\mathcal {Q}_{1}\))) does not work well. Some scaling of the variables, possibly problem-dependent, might be needed for the simultaneous optimization.

5.3.4 Mass minimization

Augmented Lagrangian formulation for stress-constrained mass minimization has been studied in multiple recent works targeting 3D large scale optimization (da Silva et al. 2020; Senhora et al. 2020). We also performed a number of tests using different solution strategies for mass minimization. The optimized layout and stress distribution of the L-shape beam are shown in Fig. 15, while the statistics are reported in Table 4. In this example, ADMM achieves a slightly smaller volume fraction. However we note that this involves empirically selecting parameters, and the same set of parameters does not necessarily lead to consistent advantages in other examples. The objective of augmented Lagrangian formulation is written as

It consists of three terms, a scaled mass, the Lagrangian term, and the penalty term. The scaling factor κ is introduced to balance the three terms, following Senhora et al. (2020). The best set of parameters that fits different examples in the context of mass minimization remains to be found. We leave this as important future work.

Final layouts of L-shaped structure on stress-based mass minimization topology optimization

Full size image

In this paper, we have presented an interpretation of local stresses as optimization variables. The introduction of auxiliary stress variables enlarges the optimization space. It leads to an augmented Lagrangian formulation which can be solved by the alternating direction method of multipliers. Using simple truss examples, we have validated that the augmented Lagrangian formulation allows reaching the global optimum that is located in the degenerate sub-space. Numerical results on continuum structures suggest that the augmented Lagrangian formulation, solved using different solution strategies, leads to compliance values that are up to 9.1% smaller than that from the conventional formulation. In particular, using ADMM for solving the augmented Lagrangian formulation continued by MMA for the conventional formulation is most favorable in terms of both the objective and computation time.

The departure point of our approach is to introduce a second set of optimization variables corresponding to local (stress) constraints. We envision this is applicable to other topology optimization problems with local constraints, such as the minimum feature scale (Zhou et al. 2015) and porosity (Wu et al. 2018). The computational benefit for reformulating such problems is yet to be investigated.