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

Handling stress constraints is an important topic in topology optimization. In this paper, we introduce an interpretation of stresses as optimization variables, leading to an augmented Lagrangian formulation. This formulation takes two sets of optimization variables, i.e., an auxiliary stress variable per element, in addition to a density variable as in conventional density-based approaches. The auxiliary stress is related to the actual stress (i.e., computed by its definition) by an equality constraint. When the equality constraint is strictly satisfied, an upper bound imposed on the auxiliary stress design variable equivalently applies to the actual stress. The equality constraint is incorporated into the objective function as linear and quadratic terms using an augmented Lagrangian form. We further show that this formulation is separable regarding its two sets of variables. This gives rise to an efficient augmented Lagrangian solver known as the alternating direction method of multipliers (ADMM). In each iteration, the density variables, auxiliary stress variables, and Lagrange multipliers are alternatingly updated. The introduction of auxiliary stress variables enlarges the search space. We demonstrate the effectiveness and efficiency of the proposed formulation and solution strategy using simple truss examples and a dozen of continuum structure optimization settings.


Introduction
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 subspaces, 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. (2016Verbart et al. ( , 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 stressconstrained 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.

Alternating optimization of design and stress
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.

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., α e = σe .
( 5 ) This conventional formulation is denoted as 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 σ 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, K e (ρ e ) = E e (ρ e )K 0 , (7) where K 0 is the stiffness matrix of a solid element with Young's modulus E 0 .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).

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 (3), ( 4) and (5).
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., σ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.

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 3) and ( 5).
(11) 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., ρ [0]   e = γ, ∀e, and compute α [0] = σ (ρ [0] ).The initial value of Lagrange multipliers is λ [0]  e = 1, ∀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 min Considering the bound constraint, the solution of the second sub-problem is

Test on a two-bar truss example
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, σ lim .
The design variables are the cross-sectional areas A = [A 1 , A 2 ] T .The cross-sectional area is upper bounded by A max = 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 The initial problem formulation, denoted by 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 .
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 Illustrated in Fig. 2b, the p-mean aggregated formulation, P M 0 , enlarges the solution space to additionally include the shaded region E F G. Curve F G corresponds to the boundary obtained when p = 8 in the p-mean function (Verbart et al. 2017).
Fig. 1 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 1 and A 2 , under stress limit σ lim

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 σj = 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 subproblems, i.e., (26).

Results and discussion
We compare six solution strategies for solving three formulations, P 0 , P M 0 , and P 1 .The results are shown in Fig. 3.In all cases, the initial values of A 1 and A 2 equal to A max = 2, corresponding to C in the solution space.The global optimum is located at G(1, 0).Using MMA to solve 1 The scaling by Amax is included for convergence to the global optimum.
P 0 leads to a local minimum point F (0, 1) in the polygon region (Fig. 3a).With the p-mean aggregated formulation P M 0 , 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 P 1 using MMA and ADMM, respectively.Lastly, Fig. 3f shows the results of a hybrid approach: we use ADMM to solve P 1 , and, after five iterations, switch to MMA for solving P M 0 .From Fig. 3, it can be observed that the global optimum (G) is correctly identified using the aggregated constraint formulation P M 0 and the augmented Lagrangian formulation 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 α = σ 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 (P 1 ) + MMA (P 0 ), is promising in terms of the number of iterations and wallclock time.
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 Implementation of stress-constrained continuum structure optimization

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-andwhite designs (Guest et al. 2004;Wang et al. 2011).

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}.

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.

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 Fig. 4 The convergence of the total mass using different solution strategies where q is a relaxation parameter (q = 0.5).σ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 0 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.

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, is computed by Detailed sensitivity analysis is presented in Appendix 1.

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 (σ 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 Q 1 , followed by using MMA for Q 0 (see Appendix 2 for solving 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.

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 σ lim = 0.5.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 (Q 0 ), while the second row corresponds to the proposed formulation with two sets of optimization variables (Q 1 ).Note that in the last approach, ADMM (Q 1 ) + MMA (Q 0 ), we use MMA to solve the conventional formulation rather than the new formulation which has a doubled number of variables.
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, Q 0 and Q 1 .The methods are comparable in terms of constraints on the total volume, V , and maximum stress, max( σ ).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.6vs 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  2 Statistics of minimizing the compliance of an L-shaped beam.In the last row, #Solves refers to the number of solving the linear systems, including both the state equation and adjoint analysis  GCMMA and MMA achieves a comparable compliance (295.91).
The objectives from the category of Q 1 outperform those from the conventional formulation Q 0 .ADMM (Q 1 ) achieves a compliance of 281.99, which is 6.48% smaller than that by MMA (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( σ ) − σ 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 (Q 1 ) + MMA (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).

Step-shaped structure
The second set of experiments is performed on a stepshaped 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 complianceminimized 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: σ lim ∈ {0.6, 0.7, 0.8}.
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.
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 (Q 0 ) by 8.78% (σ lim = 0.6), 7.73% (σ lim = 0.7) and 3.16% (σ lim = 0.8), while for the hybrid strategy, ADMM (Q 1 ) + MMA (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 (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  3 Data statistics for optimizing a step-shaped structure under three different stress limits Fig. 11 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 . The triple indicates the increasing stress limits, from left to right, 0.6, 0.7 and 0.8 Fig. 12 Optimized structural layouts and stress distributions in the step-shaped design domain using different optimization methods 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 σ 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).

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 (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 (Q 1 ) on average attains the smallest compliance, followed by the hybrid approach, ADMM (Q 1 ) + MMA (Q 0 ).On the right, the comparison shows that the computation time of the hybrid approach is marginally smaller than that of MMA (Q 0 ).

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 σ .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

Initial iteration number
In the ADMM (Q 1 ) + MMA (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, 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 (Q 1 )) does not work well.Some scaling of the variables, possibly problem-dependent, might be needed for the simultaneous optimization.

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.

Conclusion
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 subspace.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.

Appendix 1. Sensitivity analysis
The first term on the right-hand side in (38), ∂fc(ρ)  ∂φe , is the same as in conventional density approaches (e.g., Wang et al. 2011) With the definition of σe given by ( 34), the derivative ∂ σi ∂ρ j has two cases,

Derivative of the von Mises stress w.r.t. the stress components
The von Mises stress σi is computed from the Voigt notation of stress tensor σ i = [σ i,x , σ i,y , τ i,xy ] T by ( 35) and ( 36).Thus, the derivative of the von Mises stress with respect to its stress components is where

Derivative of the stress components w.r.t. the density
From the stress form, (37), the derivative of stress is given by It requires the derivative of the displacement field regarding the density.This is obtained by differentiating (6), Since the load vector F is assumed to be constant, the above equation leads to where ∂K ∂ρ j is derived from the function of K in (7).By substituting (A1.8) into (A1.6),we get where the subscript () i indicates extracting the entries that correspond to the degrees of freedom of element i.

Appendix 3. More examples
The design domain and boundary conditions of six more tests are illustrated in Fig. 16 on the left of each example.The three C-shaped structures share the same design domain but with different boundary conditions.The concentrated forces are applied over a short span of 4 elements horizontally.Fixations are applied either to the entire boundary edge or ten elements.In (d) and (f), a distributed force with a sum of F = 5 is applied on the top edge.Figure16, on the right of each example, shows the stress distribution in the structural layout optimized without stress constraints.The stress limits used in stress-constrained compliance minimization are 0.4, 0.4, 0.4, 0.4, 2.5, 2, respectively.The results of different solution strategies are summarized in Table 5, while the optimized structures are shown in Figs. 17,18,19,20,21 and 22.

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

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

Fig. 5 Fig. 6
Fig. 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)

Fig. 8 Fig. 9
Fig. 8 Plots of convergence of the volume and stress constraint (max( σ ) − σ lim ) for optimizing the L-shaped beam

Fig. 10
Fig. 10 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)

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

Fig.
Fig.In each example, from left to illustration of the design domain and boundary conditions, the optimized material distribution without stress constraints, and the corresponding stress distribution

Fig. 18
Fig. 18 Optimized structural layouts and stress distributions in the C-shaped structure 2 using different solution strategies

Table 1
Statistics in solving the two-bar truss optimization problem

Table 4
Statistics of stress-constrained mass minimization on the L-shaped beam.γ is the volume ratio, i.e., γ = ρ T v/1 T v , and thus is left out here.The second and third terms involves ∂ σ ∂φe .Consider element i.Using the chain rule, ∂ σi ∂φe is calculated byAlternating optimization of design and stress for stress-constrained...
The derivative of the maximum stress regarding to φ e is stated below:

Table 5
Statistics for optimizing six different structures/boundary conditions