A generalized approach for robust topology optimization using the first-order second-moment method for arbitrary response functions

The paper presents a rigorous formulation of adjoint systems to be solved for a robust design optimization using the first-order second-moment method. This formulation allows to apply the method for any objective function, which is demonstrated by considering deformation at certain point and maximum stress as objectives subjected to random material stiffness and geometry, respectively. The presented approach requires the solution of at most three additional adjoint systems per uncertain system response, when compared to the deterministic case. Hence, the number of adjoint systems to be solved is independent of the number of random variables. This comes at the expense of accuracy, since the objective functions are assumed to be linear with respect to random parameters. However, the application to two standard cases and the validation with Monte Carlo simulations show that the approach is still able to find robust designs.


Introduction
Design optimization under uncertainty has become a broad field in structural optimization. The fact that the variability of input parameters affects the optimal design has motivated a large variety of approaches embedding uncertainty quantification into design optimization. An overview of different approaches is, for instance, given by Park et al. (2006), Schuëller and Jensen (2008), and more recently by Kanno (2020).
The current paper focuses on a probabilistic consideration of uncertainties, which, depending on the problem formulation, is referred to as robust design optimization (RDO) or reliability-based design optimization (RBDO). While RBDO approaches constrain the probability of failure, RDO approaches typically aim at minimizing the mean and the standard deviation of an objective function. Both types of problems can be solved with sampling-based approaches, such as Monte Carlo (Schevenels et al. 2011) or stochastic collocation methods (Lazarov et al. 2012a). However, even when boosting such an approach with surrogate models, the computational cost is significantly higher than in a deterministic optimization, especially for high-dimensional optimization problems.
Probabilistic methods that are based on a Taylor series expansion offer more computational efficiency. For RBDO problems, the First-Order Reliability Methods, also referred to a Hashofer-Lind method, is often employed. It is based on an approximation of the limit state function, which divides the random space into a failure region and a safe region. Finding the relevant limit state(s) requires the solution of an optimization problem. Due to the computational cost caused by embedding this optimization into the design optimization, a sequential approach is often used (Schuëller and Valdebenito 2010).
Here, a Taylor expansion at the mean vector of the random parameters is employed to estimated the stochastic moments and their gradients. When considering only the linear therms, this approach is referred to as first-order second-moment (FOSM) method (Elishakoff et al. 1987), mean centered FOSM (Der Kiureghian 2022) or mean value FOSM approach (Haldar and Mahadevan 1999). 1 A similar approach is the perturbation technique, where the state equation (i.e., static equilibrium) is expanded at the mean vector instead of the objective function. In any case, the approaches based on a Taylor expansion at the mean vector require only one evaluation of the objective function. However, determining the gradients of the objective and the gradients of the stochastic moments can be computationally expensive.
The first approach for RDO using the adjoint method in conjunction with the perturbation technique was presented by Doltsinis and Kang (2004); Doltsinis et al. (2005). Here, a Taylor series expansion of the equilibrium condition is carried out, which yields the adjoint systems to be solved to get all information for the gradients of mean and variance of the objective function. The number of system of equations to be solved and hence, the computational cost, increases quadratically with the number of random parameters.
In topology optimization, one design and often one random variable is usually related to each finite element. Different approaches for topology optimization emerged over the last decades (Sigmund and Maute 2013). In this work, the authors focus on the well-known solid isotropic material penalization (SIMP) approach (Bendsøe 1989). Lazarov et al. (2012b) applied the perturbation technique to topology optimization, considering random geometry. The geometric randomness is modeled via spatially scattering projection threshold, which is considered as a random field. Hence, the number of random parameters obtained from discretizing the random fields equals the number of finite elements. In order to reduce computational cost, the number of random parameters are reduced using discrete Karhunen-Loève transformation.
Kriegesmann and Lüdeker (2019) used the first-order second-moment (FOSM) method for robust design optimization, which is also based on a Taylor expansion. In difference to the perturbation technique, here the objective function itself is expanded. By applying the adjoint method directly to the variance obtained by the FOSM method, only one additional adjoint system needs to be solved per iteration. Hence, the computational cost is less than two times the deterministic analysis. However, the approach of Kriegesmann and Lüdeker (2019) is only applicable for the compliance as objective obtained from a linear analysis.
The current paper provides a generalized framework for RDO using FOSM. Thereby, the approach can easily be applied to different objective functions. This is demonstrated by considering a displacement (other than load introduction point) and maximum stress as objective functions subjected to random material stiffness and geometry, respectively. Stresses are of special interest, due to their local nature which causes sensitivities with respect to small, local variations. da Silva and Cardoso (2017) considered stress constraints in robust topology optimization (RTO), but subjected to spatially scattering material stiffness and applied the augmented Lagrangian method to handle the high number of stress constraints. Random geometry is also considered in conjunction with stresses (see da Silva et al. 2019), but here the geometric uncertainty is modeled by an uniform variation of the projection threshold (compare to Wang et al. 2011). In the current paper, the focus is on a maximum stress approximation subjected to spatially scattering projection threshold.
The work is organized as follows. In Sect. 2, the generalized first-order approach is described. Sect. 3 briefly recapitulates the state-of-the-art topology optimization framework, stress calculation, and optimization algorithm. Followed by numerical results of the two mentions use cases in Sect. 4. A conclusion is drawn in Sect. 5. Details on derivatives and stress aggregation are given in the Appendix.

Generalized first-order approach for RDO
Structural optimization is almost always strongly connected to numerical methods (e.g., finite elements), which approximate the solution of partial differential equations (PDE) representing a physical state variable u (deformation, velocity, temperature, etc.). These methods result in a non-linear algebraic system of equations R = 0 , called state equilibrium. A general structural optimization setup seeks to minimize a certain objective function f regarding the inequality constraints c i , lower and upper bounds for the design variables v , and the state equation.
This problem can be solved efficiently by gradient-based optimization algorithms, if analytical sensitivity information is provided (Sigmund 2011).
Many programs or software libraries exist for solution preparation of the extremely problem depended equilibrium R . Therefore this constraint is usually not handed to the optimization algorithm, but is considered within the gradient information of any other design response Φ . The later may be any objective f or constraint c i , which depend on the physical state variable u . In topology optimization one usually prefers an adjoint sensitivity analysis (see Michaleris et al. 1994;van Keulen et al. 2005) The Lagrange multiplier Φ can be computed by the solution of a single system of equations: Thereby the optimization algorithm is restricted to solutions u(v) that satisfy R = 0.

RDO using FOSM
RDO is considering a probabilistic system g( , , ) depending on design, state, and random variables . In structural mechanics g may represent typical system responses like compliance, stress, or deformation.
FOSM is an efficient possibility to perform a rough probabilistic analysis, that only approximates mean g and variance Here, is the mean of the random vector and its symmetric covariance matrix. RDO utilizing FOSM considers at least one probabilistic design response Φ g in problem (1), which depends on the mean g and/or the variance 2 g of the probabilistic function. The total sensitivity of Φ g with respect to the design variables has to be determined for the gradient-based optimization algorithm and can be computed as A typical choice for the probabilistic design response Φ g considered in RDO is given by a weighted sum of the mean g and the standard deviation g : Then, Eq. (6) can then be specified to Dv .

Generalized sensitivty analysis of FOSM approximations
The focus is now on the terms D g ∕Dv and D 2 g ∕Dv . The sensitivity of the approximated mean equals the sensitivity of the probabilistic function and can be computed equivalently to any deterministic response function using Eqs. (2) and (3). Now, the discrete equilibrium and therefore also the state variable may also depend on the random variables, i.e., The sensitivities with respect to the random variables within Eq. (5) are determined in the same way as the sensitivity with respect to the design variables: Hence, the approximated variance (5) can be rewritten: The direct differentiation of the approximated variance involves the derivative of the Lagrange multiplier g with respect to the design variables. In order to circumvent its computation, we introduce a new residuum, defined as The sensitivity of the approximated variance with respect to the design variables is obtained by introducing the following Lagrangian: Here, we introduced two additional Lagrange multipliers u and . Differentiation of L yields In order to avoid the direct differentiation of the unknown terms Du∕Dv and D g ∕Dv in Eq. (14), the Lagrange multipliers u and are derived from the following adjoint system. 2 (10) The second row is solved first to eliminate the second column of the first row. Hence, DR∕Du is the only matrix that has to be factorized for the entire probabilistic sensitivity analysis. Furthermore, this Then, the sensitivity of the approximated variance with respect to the design variables reads: An overview of the required partial derivatives (details are given in Appendix) for the total sensitivity analysis of the FOSM approximated variance is summarized in Table 1.
Depending on the actual choice of the probabilistic function and the random variables, some of the second-order derivatives may be non-sparse large matrices. None of those has to be factorized or inverted. Nonetheless, it has to be emphasized, that even storing these matrices and/or performing the according multiplications to set up the adjoint system (15) may be computational extremely demanding and can even end up dominating the computing time for the entire procedure.
The equations provided so far include the case that the distribution of the random variables may be a function of the design variables. For the FOSM approximation it means that and may be functions of v . For a detailed discussion on the design dependency of the random parameter distribution refer to Kriegesmann (2020).

Topology optimization and structural response functions
In the remainder of the paper, the first-order approach for robust design optimization is applied to topology optimization. In this section, a brief recapitulation of the considered and established methods is given. Additionally, two different combinations of structural response function and random parameter are selected for demonstration in Sect. 4 and described in the following.

Topology optimization
The design variables considered are the pseudo elemental densities v i ∶= e , each related to one element of the design domain. Those densities are then filtered using the commonly known variable filter (described by Bourdin (2001) and Bruns and Tortorelli (2001). These, filtered densities ̃e are projected to the variable ̄e , using the simplified formulation by Wang et al. (2011) originating from (Xu et al. 2010).
Projected densities ̄e are also referred to as physical densities, since they are used to scale the elemental stiffness matrix K e = E e K 0 e by employing the modified SIMP scheme (Sigmund 2007).
with E min = 1e − 9 and p = 3 . K 0 e is the elemental stiffness matrix for unit Young's modulus. E 0 is the actual material stiffness parameter to be varied or predefined.

Optimization algorithm
For all results shown, we apply the MATLAB implementation of the method of moving asymptotes (MMA) for the design variable update (Svanberg 1987). To ensure convergence for stress-based optimization, we use external move limits. Based on numerical experience, we find a design variable change in the range of ±0.05 is a good compromise between speed and stable convergence. Furthermore, the internal move limits are set as = 0.01 , = 1.2 , and = 0.7.

Nodal displacement as objective function
The first structural response considered is the nodal displacement u k , which can be expressed as Footnote 2 (continued) factorization does not require any further computation, since it has been computed to obtain the physical state variables already. with e T k being the kth unit vector, where k is the degree of freedom corresponding to the desired output displacement to be minimized.
Since u k is only explicitly depending on u the only nonzero partial derivative from Table 1 is given in Eq. 20. This is true for all considerable random parameters.

Compliance as objective function
The often-used structural response compliance c can also be employed as objective function and is included in this section due to its popularity.
with f being the load vector.
The compliance c is only explicitly depending on u and f . Thus, the only non-zero partial derivative from Table 1 are The later is only to be considered if random load is studied. The fact that compliance is self-adjoint makes it so popular and easy to implement. Making use of this fact, the proposed generalized framework could be further simplified, increasing the efficiencies for a compliance objective. As derived in Kriegesmann and Lüdeker (2019) this means, the state equation and one adjoint system need to be solved, whereas for the general approach two additional adjoint systems are required. Nonetheless, the main purpose of this work is generalization. Simplification is left to the interested reader. Still, compliance can easily be applied to the proposed approach.

Elemental von Mises stress as objective function
Elemental von Mises stresses q e are calculated at the element's centroid with M being the constant von Mises matrix. Assuming plane strain/stress state B e is the strain-displacement matrix and C e = E e C 0 is the constitutive of element e, where C 0 is similar to K 0 e the constitutive matrix for unit Young's modulus. E e contains the scale factor for actual material stiffness E 0 and the RAMPinterpolation scheme 3 (Stolpe and Svanberg 2001).
with an interpolation parameter of p = −0.5 . Corresponding partial derivatives of q e are given in Appendix . Note, for stresses all required derivatives stated in Table 1 are nonzero and need to be determined.

Stress aggregation
Typically, the maximum stress is the quantity of interest in stress-based optimization, but since the max function is nondifferentiable, the maximum stress is approximated by the upper bound KS function (Kreisselmeier and Steinhauser 1979).
where q max = max(q) . Due to a rapidly changing design and maximum stress value during optimization, the aggregation parameter is updated every iteration, to ensure a constant approximation quality. As proposed by Jansen et al. (2013) the numerator in = 15∕q max is predefined and kept constant during optimization. The aggregation parameter is updated accordingly. Note, during the first iterations, the predefined value is set to 40 and the aggregation parameter is kept constant until ∕q max < 15 . Values are chosen based on numerical experiments and are a good trade-off between a stable, converging optimization and a good maximum approximation. Initially the value is increased, since it does not compromise the convergence and helps to avoid local minima, where, for e.g., stress singularities at sharp corners are not rounded of. All needed derivatives of the KS function can be found in Appendix . Also here, all required derivatives stated in Table 1 are non-zero and needs to be determined.

Numerical results
In this section, the well-known inverter (see Fig. 1) and the common benchmark example for stress-related problems, the L-Beam (see Fig. 2), are chosen for demonstration purpose. The proposed generalized FOSM approach is applied to a displacement objective considering random Young's modulus and maximum stress with variations in the projections threshold (17), respectively. In both cases, robust designs are obtained. Validation is performed by means of the Monte Carlo method. Shared parameters for the depicted demonstrators are given in Table 2. For all results shown, the volume is constrained to a fraction of the whole design domain. The inverter's volume fraction is constrained to 35% and the L-Beam's is 40%. In the following section, a linear finite element model ( R = Ku − f ) is considered and a regular mesh with quadrilateral, bi-linear, isoparametric continuum elements, and isotropic material is chosen. Derivatives of R are given in Appendix. Plain stress condition is assumed.

Inverter and spatially scattering Young's modulus
The inverter is a popular example also in RTO (see (Lazarov et al. 2012a) and da Silva et al. (2019) for instance). The optimization objective is the minimization of the horizontal displacement u h (see Fig. 1). The design should be symmetric, but still take into account an asymmetric spatially scattering material stiffness. For that, the whole model is calculated, but symmetry is enforced for the design variables during optimization. Due to the asymmetric random field the vertical displacement u v ≠ 0 , which is the reason we included its standard deviation v in the robust objective: with h and h being the mean and standard deviation of u h , respectively. The weighting factor = 10 for all inverter examples. Note that can be interpreted as the reliability index related to the chosen distribution type and the probability of failure, as explained in Kriegesmann and Lüdeker (2019) in Sect. 3.3. Here, the value is increased to provoke a visual effect due to the low coefficient of variation, as will be explained in the following. Furthermore, the mean value and standard deviation are only approximated values, which might differ from the real counterparts. Thus leaving room to adjust for finding a good trade-off between final robustness and mean value. Usually, the displacement in horizontal direction has a positive value at first, approaches zero, and then starts to become negative until converging to some value of u h < 0 . Numerical experiments showed if the initial weighting factor is set too high, the MMA is not able to reduce u h below zero. Thus, the initial value of int = ∕1000 is kept constant for the first 10 iteration and then increased in each step until it reaches its predefined value 40 iterations later. For the deterministic case only u h is to be minimized.  Load input is set to F in = 2 . k in = 1 and k out = 0.01 are the input and output spring stiffness, respectively. The optimization objective and constraint values are scaled by a factor of 100 to gain well-posed function values. As in Silva and Cardoso (2016), da Silva and Cardoso (2017), and Lazarov et al. (2012a), the spatially scattering elemental Young's modulus is modeled as a random field using the parameters given in Table 3.
From an engineering perspective, scattering material stiffness might especially be relevant for additive manufacturing. Thillaithevan et al. (2022) showed that tolerances could lead to spatially degenerated material properties for microscale geometry. In the case of fiber-reinforced composites, local varying fiber orientations lead to a significant impact on the material stiffness as investigated by Rauter and Lammering (2020) or in the field of civil engineering variations in soil stiffness can even be measured in areas considered as homogeneous (see Grabe 1994).
The resulting physical density fields are depicted in Fig. 3. Figure 3 shows the deterministically optimized inverter with bulk material agglomerated at the input node, which is connected to the rest of the structure via two hinges. The robust designs avoid this by directly connecting the input node.
Removing v from the robust objective yields designs similar to the deterministic one. Thus, obviously avoiding those two hinges and directly connecting the input node helps to reduce the variation in vertical displacements. This is further proven by the reduced standard deviation of horizontal displacement v (see Table 4). As can be observed, v is reduced by 48.4% comparing deterministic and FOSM postprocessing values. Monte Carlo benchmark shows an improvement of 67.2. Note, Monte Carlo optimization is performed with 1000 realizations. The increased robustness comes at the expense of around 4% less horizontal displacement, which is to be expected in robust design optimization.
Evaluating the proposed approach, the approximation quality of FOSM is very good, if optimization and postprocessing values for the horizontal displacement values h and h are compared. The approximation of v on the other hand is close to a factor of three worse, but still a more robust design is obtained. Note, the coefficient of variation of the deterministic design CoV v = v ∕ h ≈ 0.5% , which is a quite low scatter, compared to the chosen Young's modulus variation CoV E = 10%.

L-Beam and spatially scattering projection threshold
The L-Beam is a standard benchmark for stress-based topology optimization (to name a few, see (da Silva and Cardoso 2017; Giraldo-Londoño and Paulino 2021). For the L-Beam example, the maximum aggregated stress q KS is considered as objective (as described in Sec. 3.5 and 3.6). The beam is loaded by a force F = 10 , which is distributed along the horizontal and vertical edge (see Fig. 2) to help reduce stress peaks. The robust objective is given by with = 3 . Note that the weighting factor had to be relaxed due to numerical reasons. To properly scale the optimization   Table 5. Figure 4 shows the obtained designs and corresponding nominal von Mises stress plots. FOSM shows the lowest nominal stress, otherwise only minor differences are spotted. Note, nominal stresses plotted represent only postprocessing values which are not considered during optimization. In general, the nominal objective value should be lowest for the deterministic case. From the optimizer's perspective this is also the case. Considering column "Optimization" in Table 6 it can be observed that the objective value for the deterministic case is 2.08 and for FOSM it is slightly increased to 2.08 + 3 ⋅ 0.001 . The higher value in 4a can thus   be explained by a local minima. This may happen due to the aggregation and is thus not completely unexpected nor avoidable. It is also important to notice that stresses are very sensitive to little changes in geometry. For instance, the top strut running from the loaded tip to left is slightly thicker for the robust designs. The sensitivity of stresses with respect to small geometry changes can also be seen in Table 6. A uniform distribution of in the interval [0.45, 0.55] corresponds to a coefficient of variation of CoV = 6% , resulting in a CoV̂q ≈ 22% , where q is the non-aggregated true maximum von Mises stress q max . It is observed, that FOSM decreases the standard deviation of the true maximum stress q by more than 12%. Also the mean value is slightly reduced, which still yields a CoV̂q ≈ 20% . For Monte Carlo a significant reduction of q by more than 84% is achieved, resulting in a CoV̂q ≈ 4% . Comparing the standard deviation estimated by the FOSM approach with the Monte Carlo postprocessing of the same design (see Table 6) reveals that the linear approximation is very inaccurate in this case. Consequently, the RDO with FOSM method by far not exploit the same optimization potential as the RDO using Monte Carlo sampling. Still, the FOSM approach provides a more robust design than the deterministic optimization.
Note that the computational time is very high if aggregated stresses are considered for the proposed generalized FOSM approach. It even exceeds Monte Carlo run time for very fine discretized models. For Monte Carlo, 1000 models are solved in parallel on 24 cores, thus making a direct comparison unfair, but still only three adjont systems need to be solved for the proposed generalized FOSM approach per iteration. For coarse discretizations FOSM is much faster in any case. The reason for this increasingly high computational time can be traced down to gradient calculation of aggregated stresses (compare Appendix). Many required partial derivatives are not only non-zero and some nonsparse ( 2 q e ∕ u i u j in (52) for instance), but also simple matrix operations on those non-sparse and big matricies take up much computing time. Generally, calculation of secondorder partial derivatives for aggregated stresses takes up most of the computational time.

Concluding remarks
A generalized approach for robust topology optimization is derived based on a linear Taylor series expansion, namely generalized FOSM, and applied to two commonly used objective functions, nodal displacement and maximum aggregated von Mises stress. In both cases, a robust design could be found with the proposed method. As expected, due to the inaccurate linear approximation, performance is below Monte Carlo benchmark.
For a displacement objective subjected to spatially scattering material stiffness, a low sensitivity with respect to the random quantity is observed. Here, a good robust design with a significant increase in robustness is achieved.
The stress objective showed a high sensitivity with respect to random geometry variations. Even though the standard deviation is reduced, the Monte Carlo benchmark showed a much better robustness. Additionally, the von Mises stress objective requires all second-order partial derivatives to be set up, which results in a very high computational effort due to big non-sparse matrices and operations on those. Especially the aggregation and corresponding gradient calculation takes up a significant amount of computational time.
Nonetheless, for single valued or low number aggregated objective functions, the proposed generalized FOSM approach is computationally very efficient, since only two additional adjoint system needs to be solved compared to the deterministic case. For self-adjoint objectives, e.g., compliance, the run time is improved. Furthermore, the proposed approach is independent of the number of random variables, which is a big benefit compared to other robust topology optimization approaches.
In future work improvements for stress-based robust optimization are the main focus. First, the run time for secondorder derivative calculation needs be reduced or possibly bypassed completely. Second, since stresses are highly nonlinear, a reciprocal extension seems promising to increase approximation accuracy for FOSM (compare to Kriegesmann and Lüdeker (2021).

Derivatives for robust topology optimization using linear analyses
This appendix summarizes the derivatives required for applying the FOSM-based RDO approach to topology optimization. Here, we consider linear analyses and standard element densities as design variables, resulting in the well-known residuum in (36). It is assumed that the distribution of random parameters is independent of the design variables and that the load f is design independent.
Furthermore, partial derivatives with respect to random variables ( and E ) are required. Considering random , firstorder derivatives are determined below. Note, due to the way the proposed generalized approach is formulated (see Sec. 2), only explicitly depending terms are differentiated.
with ̄ ∕ being the derivative of Eq. 17 with respect to .
Second-order derivatives are obtained as follows: Kriegesmann and Lüdeker (2019) already derived the second-order partial derivative of K , but Eq. (28) on page 274 is incorrect. The second-order partial derivative of K in (40) cannot directly be determined by differentiating K∕ with respect to ̄ , since ̄ is just a simple substitution in terms of differentiation. Thus, K is differentiated with respect to the filtered variable ̃j and rewritten as follows: In line two of the above equation, the partial derivative ̄ ∕̃ is placed outside the brackets causing the inverse ̃ ∕̄ for the second therm to occur. Derivatives for random are obtained accordingly, without the need to differentiate K∕ E with respect to ̃ , since it is directly differentiable with respect to E and ̄ ..
So far derivatives with respect to the physical densities ̄ are derived for the objective function f. Applying the chain rule yields the derivatives with respect to the design variables .

Stress derivatives
Here, all required derivatives of the von Mises elemental equivalent stresses are given with respect to the projected design variables ̄ and the random parameter (compare to (17). Dependencies are as follows: First-order derivatives are given below. Note, due to the way the proposed generalized approach is formulated (see Sec. 2), only explicitly depending terms are differentiated. The second-order partial derivative of C e in (57) is determined the same way as the second-order derivative of K (for details see (43).

General KS aggregation
First-and second-order derivatives of the upper bound KS function (compare to Eq. 27) are as follows: where f 0 is considered constant during differentiation as derived in Kranz et al. (2021). The first-order derivative of f KS with respect to an arbitrary variable x 1 , assuming the function value f i = f (x 1 ) and substituting Σ e = ∑ n e i=1 exp (fi−f0) yields The second derivative with respect to another arbitrary variable x 2 , which can be the same variable as x 1 or a different one, again assuming the function value f i = f (x 1 , x 2 ) yields All partial and final derivatives are checked by means of finite differences.