Robust structural optimization in presence of manufacturing uncertainties through a boundary-perturbation method

Most manufacturing processes are inevitably characterized by process tolerances that ultimately affect the way a component behaves and complies with the design requirements. These uncertainties determine the real performance of a structure, with their impact growing with increasing deviations from the nominal values. This work introduces a simple approach, applicable to both static and dynamic cases, to conduct robust structural topology optimization in presence of manufacturing uncertainties. This approach, based on the level set method, makes use of a computationally efficient boundary-perturbation technique to describe over- and under-etching errors. Compared to the existing methods, it does not require a frequent re-initialization of the level set function, nor does it require a mapping between the etched structures and the nominal one. Moreover, compared to the standard case with uniform uncertainty, the technique presented in this work allows dealing with arbitrary spatially varying errors without increasing the computational cost.


Introduction
While most topology optimization techniques deal with deterministic geometries, real structures are usually characterized by a number of uncertainties. These can be caused by, e.g., material properties, loading and boundary conditions, manufacturing errors, and others. In cases where uncertainties are likely to arise, a structure optimized using a deterministic topology optimization approach (Bendsoe and Sigmund 2003) will generally lead to sub-optimal performance or, at worst, failure. The aim of Robust Topology Optimization (RTO, Beyer and Sendhoff (2007)) is to optimize the required performance measure while minimizing its variance due to external uncertainties that might affect the structure. Several approaches exist for robust topology optimization. These methods are usually classified into non-probabilistic and probabilistic, based on how they represent and manage uncertainties (Beyer and Sendhoff 2007).
Replacing the stochastic quantities with their average values (Asadpoure et al. 2011) is arguably the simplest and most straightforward probabilistic approach to RTO. On the other hand, the most complete probabilistic approach would be to take the higher-order statistics into account. For instance, Asadpoure et al. (2011) developed a probabilistic method, applied to truss structures, to deal with Gaussian uncertainties that affect the stiffness matrix, e.g., material properties. The aim is to exploit the higher-order statistics to obtain a layout that is less sensitive to the considered uncertainty. As an alternative, the Reliability-Based Topology Optimization (RBTO) (Kim and Kwak 1996;Kharmanda et al. 2004;Allen and Maute 2005) aims at optimizing a given performance, while including in the constraint the effect of the uncertainty to satisfy a reliability criterion with a prescribed confidence level. A much different approach is the Risk-Averse Topology Optimization (RATO) (Conti et al. 2011;Martínez-Frutos et al. 2018;Martínez-Frutos and Ortigosa 2021), which aims at minimizing a risk measure that somehow quantifies the loss in performances due to the uncertainty. However, in presence of non-Gaussian uncertainties, or when dealing with manufacturing uncertainties, a probabilistic approach will result in too complex of an analytical formulation to be of practical use for most non-trivial problems.
A possible alternative is to resort to a Monte Carlo method, where a sufficiently high number of realizations is simulated at each step to approximate the statistics (Rubinstein and Kroese 2016), thus trading off a higher simulation cost with lower analytical complexity. In the field of topology optimization, this has been used by Schevenels et al. (2011) to optimize structures in the presence of Gaussian spatially varying manufacturing errors. However, due to the computational complexity involved, it quickly becomes prohibitive for problems characterized by a high number of design variables.
Another approach is to arbitrarily choose a finite number of cases that are reasonably believed to be representative of the real behavior of the structure. This allows approximating the influence of several uncertainties by assigning a probability to each possible case. An example of this is the multi-load approach, suitable in presence of uncertainties in the load's position and magnitude. Guest and Igusa (2008), for example, applied this concept to truss structures subjected to uncertain loads and extended the method to consider uncertain nodal locations as well.
In the field of Density-based Topology Optimization (DTO) (Bendsøe and Kikuchi 1988;Bendsoe and Sigmund 2003), Sigmund (2009);Wang et al. (2011b) introduced erosion, intermediate, and dilation filtering stages to model the over-and under-etched structures. These projections are used for the design of compliant mechanisms and heat conduction under a minimum length scale control (see Lazarov et al. (2016) for a comprehensive review). Similarly, Andreasen et al. (2020) implemented a non-probabilistic worst-case formulation for the compliance minimization problem, whereas Wang et al. (2011a); Elesin et al. (2012) applied these projections to dynamic problems. A modified version of this technique can be used to manage Gaussian spatially varying manufacturing errors (Schevenels et al. 2011).
Methods to deal with load uncertainties are well-established in the field of Level Set Topology Optimization (LSTO) (Allaire et al. 2002;Wang et al. 2003;Allaire et al. 2004) as well. For instance, Dunning et al. (2011a) introduced a method, equivalent to the multi-load approach, to minimize the expected compliance in presence of uncertainties both in the load's magnitude and direction. On the other hand, Dunning and Kim (2013) proposed the minimization of the expected value and the variance of the compliance in presence of an uncertain load magnitude.
If we consider the optimization in the case of other sources of uncertainties, however, far fewer publications exist that make use of a Robust Level Set Topology Optimization (RLSTO) approach. Compared to density-based techniques, level-set approaches are intrinsically able to model the boundary of the structure without the need for any projection filter, thus providing a straightforward method to describe manufacturing errors (Chen and Chen 2011). Among these are the works of Jang et al. (2012); Andreasen et al. (2020), who dealt with the etching problem in a static case, and Li et al. (2019), who showed how to tackle the maximization of the first natural frequency in presence of possible over-etching. These methods make use of the signed-distance property to model uniform under-etched and over-etched structures using contours of the level-set function at offset levels. In this sense, a re-initialization scheme is required to maintain the signed-distance property, even though it may cause problems in the optimization process (Jang et al. 2012). Moreover, mapping approaches (Chen and Chen 2011;Li et al. 2019) are required to map the perturbed domain to the original one.
This paper presents a simple Robust Level Set Topology Optimization (RLSTO) approach, based on a worst-case formulation, to deal with both uniform and spatially varying manufacturing uncertainties, with examples spanning both the static (compliance minimization and compliant mechanism) and dynamic (frequency assignment) domains. The etching error is treated as an unknown (but bounded) uncertainty, which is characterized by an unspecified probability distribution bounded between predefined lower and upper values. In order to consider the etching error in the optimization, this method makes use of two additional structures: the reduced (over-etched) and the expanded (under-etched) realizations. They are efficiently obtained through the Level Set Method (LSM) (Osher and Sethian 1988;Sethian 1999;Osher and Fedkiw 2006) by perturbing the nominal interface under a velocity field that represents the manufacturing uncertainty, either uniform or spatially varying. Therefore, the computational cost associated with both cases is the same. In addition, this level-set approach directly relates the manufacturing error with the velocity field used for the perturbation, thus allowing the description of a variety of different etching fields without using any projection filters (Sigmund 2009;Wang et al. 2011b;Schevenels et al. 2011). By properly defining the perturbing velocity field, there is the possibility to apply the etching error only to specific regions of the design domain, thus leaving the other ones untouched. Moreover, since this method does not rely on the signed-distance property, a frequent re-initialization of the level-set function is not required, thus further improving the overall computational efficiency. In addition, using the Least-Squares Interpolation (LSI, Dunning et al. (2011b)) to compute the shape sensitivities, the proposed approach does not require any mapping between the original and perturbed domains (Chen and Chen 2011;Li et al. 2019).
The extension to other physical domains is also straightforward, since the method acts on the level-set function itself, making it highly independent of the underlying physics description. It is important to point out that in other physical domains, such as those involving fluid-structure interaction, the three realizations may not be sufficient to describe the behavior of the structure. Nevertheless, the boundary-perturbation technique presented herein still applies. In fact, if required, more than three realizations can be obtained at a low computational cost, and successively combined to formulate an optimization problem more suitable to the particular case.
It is worth mentioning that there are other techniques (such as polynomial chaos expansion), to perform a robust optimization in the presence of manufacturing uncertainties (Tootkaboni et al. 2012;Lazarov et al. 2012;Chen and Chen 2011;Zhang and Kang 2017). These methods are used to analytically describe the manufacturing uncertainty field, thus allowing the computation of the response statistics and their sensitivities. These techniques, however, increase the analytical complexity of the problem.
The rest of the paper is organized as follows. Section 2 briefly presents the necessary LSM background and the deterministic LSTO process, followed by the analytical description of our boundary-perturbation RLSTO in the case of manufacturing uncertainties. Numerical examples regarding both the static and dynamic domains, including a case in which a spatially varying manufacturing error is considered, are presented in Sect. 3. The results are then validated in Sect. 4 using COMSOL Multiphysics©. Section 5 closes with the conclusions.

Level-set topology optimization
In this Section, we briefly illustrate the level-set method. After the initial overview of deterministic optimization, the boundary-perturbation RLSTO formulation is presented.

Level-set function
The level-set function (Osher and Sethian 1988;Sethian 1999;Osher and Fedkiw 2006) is used to implicitly define the boundary Γ of a N-dimensional structure as the zero contour of a (N +1)-dimensional implicit function Φ ( Fig. 1): where Ω is the region belonging to the structure.
In most problem settings, the level set function Φ is a signed-distance function, characterized by a unitary gradient, so that:

Deterministic topology optimization
Among the advantages of the LSM are its capability of implicitly describing the interface of a structure, the computational convenience of propagating the boundary only (instead of the whole domain), and the absence of filtering.
A generic topology optimization problem can be stated as, where g represents the vector of inequality constraints with a bound g u .
Problem (3) is solved iteratively by propagating the boundary at each step, according to the Hamilton-Jacobi equation, where V n (x) indicates the optimal velocity field normal to Γ.
At each iteration, the Marching Squares algorithm (Lorensen and Cline 1987;Maple 2003) is used to discretize the interface Γ , thus allowing to compute the equivalent area (or volume) fractions in each element. The Finite Element Method (FEM) is then applied to compute the solution to the structural problem. Here, the structural matrices, i.e. stiffness K and mass M, are assembled using the ersatz material approximation with a linear interpolation scheme for the material properties (Allaire et al. 2004;Dambrine and Kateb 2010). Subsequently, the Sequential Linear Programming (SLP) formulation (Dunning and Kim 2015) is used to solve the sub-optimization problem that allows obtaining the optimal boundary velocities V n (x) . This is done using IPOPT (Interior Point OPTimizer, Wächter and Biegler (2006)), an open-source library for large-scale nonlinear constrained optimization. Equation (4) is then propagated in its discrete form, where i is a grid node, and k is the iteration number. Throughout this work, a fixed grid made of unitary square elements is used. The Hamilton-Jacobi Weighted Essentially Non-Oscillatory method (HJ-WENO, Osher and Fedkiw (2006)) is implemented for spatial discretization, whereas a Forward Euler approach is used for time discretization. Finally, the Fast Marching Method (FMM, Osher and Fedkiw (2006)) is implemented to extend the boundary velocity to a neighboring band of Γ . In this work, the level-set function is re-initialized every twenty iterations, to satisfy Eq. (2).
In Eq. (5), Δt = 1 is used to satisfy the Courant-Friedrichs-Lewy (CFL) condition for numerical stability (Osher and Fedkiw 2006): For the sub-optimization problem, the shape derivatives of f (Γ) and g(Γ) are required. Considering the objective function, for example, where s f is the shape sensitivity of f with respect to a boundary movement. After a discretization of Γ , this becomes: where l i is the length of the portion of the boundary associated with the i-th boundary point p i , and s f ,i is the shape sensitivity of f with respect to a displacement of p i . The shape derivatives of g(Γ) are obtained in a similar fashion.
In this work, the shape sensitivity of each boundary point is computed using the Least-Squares Interpolation (LSI, Dunning et al. (2011b)), which weights the analytical sensitivities, i.e. the derivatives of the objective function and the constraints with respect to the area fractions k . Other methods, such as boundary perturbation (Kambampati et al. 2021), are also commonly used.
The algorithm described in the previous paragraphs is repeated iteratively until convergence is reached, i.e., all constraints are satisfied and the objective function difference between the last five iterations is sufficiently small.

Robust topology optimization through boundary perturbation
This Section describes how to implement the boundary-perturbation RLSTO. Assuming a constant uncertainty along the whole boundary, in a direction normal to Γ , it is possible to define three structures (Wang et al. 2011b): • Nominal: blueprint structure, without uncertainties; • Expanded or under-etched: a structure whose features are thicker than in the nominal case, of a quantity ; • Reduced or over-etched: a structure whose features are thinner than in the nominal case, of a quantity .
In the rest of the paper, these three variations will be denoted by the subscripts n, e, and r, respectively. The etching error is treated as an unknown-but-bounded uncertainty, with upper and lower limits ± . A possible value for can be the maximum expected deviation, for which 3 is commonly used. Note that is not required to be constant, and a more generic uncertainty field (x) can easily be implemented as well. This is investigated in Sect. 3.5. The boundary propagated during optimization is that of the nominal structure. The expansion of the nominal structure may lead to features merging with each other if they are closer than 2 . Conversely, the reduced structure may see some very thin ( < 2 ) features disappear completely.
To used to describe uniform etching errors, and it requires reinitialization schemes to maintain the signed-distance property of the level-set function (Jang et al. 2012). However, these schemes may introduce inconsistencies that ultimately affect the convergence of the algorithm.
The boundary-perturbation technique presented herein allows us to efficiently overcome these limitations. The method is illustrated in Fig. 3: the reduced and expanded structures are obtained by propagating the nominal boundary, according to Eq. (5), under a normal boundary velocity field ± (x) that represents the manufacturing uncertainty. This approach does not require a frequent re-initialization of Φ , thus reducing the cost associated with this operation. Furthermore, as mentioned above, both uniform and nonuniform velocity fields can be handled at the same computational cost. In this way, different etching conditions can be straightforwardly managed by properly choosing the velocity field (x) . For instance, the etching error can be applied only to particular regions, leaving the others untouched.
After obtaining the three structure realizations, the Finite Element Analysis (FEA) can be carried out as usual. Then, the properties of the three structures can be combined as appropriate to formulate the robust optimal problem. As far as this work is concerned, the robust optimization problem takes the form of: where the objective function is a linear combination of the costs obtained by analyzing the nominal, the reduced, and the expanded structures. The parameters w n , w r , and w e are weights for the linear combination: their values lie between 0 and 1, and their sum is equal to 1. In this way, by adjusting the weights, the proposed objective function can describe different optimization problems. In particular, if w n = 1 , w r = 0 , and w e = 0 , Problem (9) falls back into Problem (3). On the other hand, if w n = 0 , w r = 1 , and w e = 0 are used, the usual formulation for worst-case compliance minimization problems is obtained (Wang et al. 2011b;Andreasen et al. 2020).
In level-set approaches, to compute the shape sensitivities of the nominal boundary points, mapping techniques are usually implemented to map the perturbed structures (reduced and expanded) onto the nominal one. For instance, Chen and Chen (2011) successfully applied the continuum mechanics large deformation theory to find point correspondence. However, this technique inevitably increases the numerical complexity of the problem.
In this work, we show how to compute shape sensitivities without using any mapping between the initial and perturbed domains. Firstly, the analytical sensitivities are obtained through a linear combination of the derivatives of each of the three cases, with respect to k : After that, the LSI method is used to compute the shape sensitivities of the boundary points of the nominal interface, so that each boundary point is associated with a sensitivity that considers the contributions of all three structures. Therefore, there is no need to map the boundary points of the reduced and the expanded structures onto the nominal ones.

Numerical examples
In this Section, we present the results obtained when applying the boundary-perturbation RLSTO to four different problems.
First, a compliance minimization problem is shown. Second, the design of a gripper, a compliant mechanism, is carried out. Third, a frequency assignment problem is presented. This is relevant, for example, for MEMS resonators (Wu et al. 2020) whose resonant frequency should be (9) min Γ w n f n (Γ) + w r f r (Γ) + w e f e (Γ) . . g l ≤ g(Γ) ≤ g u  accurate and precise despite possible manufacturing uncertainties. Fourth, the frequency assignment problem is solved under a spatially varying manufacturing error. The simulation parameters used in the following examples, as well as the material properties, are presented in Table 1.
All the results are obtained using OpenLSTO (Kambampati et al. 2018).

Compliance minimization
We start by considering the result of deterministic optimization as a baseline for comparison. For the case of minimum compliance, the optimal problem can be stated as: where K is the stiffness matrix, u is the vector of nodal displacements, f is the vector of external forces, C is the compliance, and A is the total area fraction of the structure. In this example, the upper limit of the area fraction is A max = 40% , while the magnitude of the external force ( Fig. 4a) is equal to 20 kN. Figure 4 shows the problem setting and initial design, while Fig. 5 shows the final (optimal) layout of the structure.
We now introduce the manufacturing uncertainties. In a compliance problem, the reduced structure is representative of the worst-case scenario. For instance, Fig. 6 shows the effect of a uniform manufacturing error = 0.5 mm on the nominal structure in Fig. 5. Note that the thinnest parts of the deterministic structure (Fig. 5) disappear completely in the reduced version (Fig. 6). On the other hand, the expanded structure is representative of the best-case scenario, thus it is removed from the robust optimization process. This not only reduces the computational load but also considers a more conservative case.
The robust optimal problem is stated as:  where 0 ≤ w ≤ 1 . If w = 0 , the objective function is equivalent to the one studied in Andreasen et al. (2020); if w = 1 , the deterministic case is obtained. The area constraint is imposed on the nominal structure. Figure 7 shows, for example, the results obtained with w = 0.5 . Table 2 compares the results obtained in the three cases: w = 0 , w = 0.5 , and w = 1 . In particular, note that w = 1 leads to the same solution (Fig. 5) obtained by solving the deterministic problem (Problem (11)).
Notice that, as opposed to the deterministic case (Fig. 6), no features of the robust structure have disappeared during the reduction (Fig. 7). In addition, even though the nominal compliance C n is slightly increased, the boundary-perturbation RLSTO with w < 1 leads to a lower reduced compliance C r and a smaller absolute variation ΔC . By using a weight w strictly between 0 and 1, both C n and C r are included in the objective function. On the other hand, using w = 0 , only C r is minimized, meaning that C n (and consequently ΔC ) are not directly considered in the optimization. This could (12) lead the algorithm to converge to a point of local optimality characterized by higher C r and ΔC . For instance, the values of C r and ΔC obtained with w = 0 are higher than those associated with w = 0.5 (Table 2). In the case of C r being more important than C n , a value of w close to zero should be used. In such problems, to avoid the issue of local optimality, it is either possible to try different initial conditions or to slightly increase the value of w.

Compliant mechanism
Similarly to the problem in Sect. 3.1, we start by considering the results of a deterministic case as a baseline for comparison. In particular, Fig. 8 shows the problem setting and initial design of a gripper. Due to the symmetry of the problem, only the bottom half of the mechanism is considered. The input force (green arrow) is represented by the vector f , whereas the output dummy force (red arrow) is (a) Nominal structure.
(b) Reduced structure. Fig. 7 Result of the boundary-perturbation RLSTO under a uniform manufacturing error = 0.5 mm and w = 0.5 . As opposed to Fig. 6, no feature disappears in this case (a) Design space and boundary conditions. The blue regions are assumed to be always full. The green and red arrows represent respectively the input and output forces.
(b) Initial design. where ME is the mechanical efficiency of the mechanism (Lau et al. 2001). The workpiece is modeled as a spring (Sigmund 1997) placed in the application point of the output dummy force l . Regarding the area fraction constraint, an upper limit of 20% is chosen. The sensitivity of ME is (Lau et al. 2001): where K e,k is the stiffness matrix of the k-th element, and is the solution of the adjoint equation Figure 9 shows the optimal layout along with its reduced version.
In presence of a uniform over-etching error = 0.5 mm, some features and hinges of the nominal layout disappear. This may lead to a drastic variation in the behavior of the mechanism. The RLSTO approach presented in this work aims at preventing both these problems. The robust optimization problem is stated as follows: where the objective function J is the linear combination of the mechanical efficiency of the nominal ( ME n ) and reduced ( ME r ) structures. If w = 1 , this formulation falls back into the deterministic one. On the other hand, if w = 0 , only the reduced realization is considered, yielding the optimization problem used by Wang et al. (2011b).
The results obtained with A max = 20% and w = 0 are reported in Fig. 10. As compared to the deterministic layout ( Fig. 9), no features of the robust structure disappear in the reduced version. (16) (b) Reduced structure. Fig. 9 Optimal solution of the deterministic mechanism synthesis, with A max = 20% . Due to a uniform manufacturing error = 0.5 mm, some features of the nominal structure disappears (a) Nominal structure.
(b) Reduced structure. Fig. 10 Result of the boundary-perturbation RLSTO under a uniform manufacturing error = 0.5 mm and w = 0 . As opposed to Fig. 9, no feature disappears in this case Page 9 of 13 120

Deterministic frequency assignment
The goal is to assign a prescribed value to the first angular frequency 0 of a structure fixed at both ends. The optimal problem is stated as: where K is the stiffness matrix, M is the mass matrix, the target frequency is 0,target = 16000 rad/s ( f 0,target = 8 kHz), and X is the corresponding mode shape. The structure is subjected to a maximum area fraction A max = 50%.
Provided that the mode shapes are mass-normalized, the first constraint-the characteristic equation-can be removed using the adjoint method with the Lagrange multiplier : Under the assumption of mass-normalized mode shapes, the derivative of the first natural frequency with respect to the area fraction k is (Bendsoe and Sigmund 2003): where K e,k and M e,k are respectively the element stiffness and mass matrices of the k-th element. Figure 11 shows the problem setting and the initial design, whereas the optimal design is shown in Fig. 12.

Robust frequency assignment
Under a deterministic frequency assignment setting, there is no control over the variability of the natural frequencies with respect to small topology perturbations. Thus, even a small shape change can lead to large eigenfrequencies variations. Table 3 shows the impact of manufacturing uncertainties ( = 0.5 mm) on the first natural frequency of the geometry obtained through a deterministic optimization (Fig. 12).
A RLSTO approach allows to minimize this effect. The problem formulation becomes: where 0,n , 0,r , and 0,e are respectively the first natural frequency of the nominal, reduced, and expanded structure.    Table 3 Effect of a uniform manufacturing error = 0.5 mm on the natural frequency of the structure shown in Fig. 12 f [kHz] Δf ∕f [-] Nominal structure 8.00 -Reduced structure 7.87 − 1.53% Expanded structure 8.11 + 1.44% The target frequency is once again 0,target , and the area constraint is enforced on the nominal structure. Figure 13 shows the optimal design obtained through the RLSTO method, while the numerical results are shown in Table 4. The average deviation from the nominal case is reduced to 0.39% (with a worst-case deviation of 0.47% ), as opposed to 1.49% in the deterministic case (with a worst-case 1.53% variation). The optimizer converges to a solution with no holes inside the structure (pockets). Intuitively, having a lower value of the perimeter reduces the variation of the matrices M and K following a perturbation, thus reducing the overall eigenvalue variation.

Robust frequency assignment in presence of a non-uniform manufacturing error
We show how the method can be adapted to cases where the over-and under-etching uncertainties are not constant. For our example problem, we define an etching error that is radially distributed in space (shown in Fig. 14): which is scaled by the maximum absolute error max .
For simplicity and ease of comparison, we consider the same problem setting above (Sects. 3.3 and 3.4). The robust optimal solution is shown in Fig. 15. The results of the optimizations presented above are summarized in Table 5.
Compared to the case with a uniform etching uncertainty (Table 3), the radial etching error produces a higher variation in the natural frequency of the deterministic structure Fig. 13 Result of the boundary-perturbation RLSTO applied on a frequency assignment problem in the presence of a uniform manufacturing tolerance = 0.5 mm with w = 0.8 . To be compared with the deterministic design in Fig. 12 Table 4 Effect of a uniform manufacturing error = 0.5 mm on the natural frequency of the structure shown in Fig. 13 To be compared with the much larger frequency variations obtained in the deterministic case,    (Table 5a). Nevertheless, the Boundary-Perturbation RLSTO can still reduce the frequency deviation in both the over-and under-etched structures. In fact, the target natural frequency is still matched, and the average frequency variation has been more than halved (from 2.12% to 0.99% ). Comparing the layouts shown in Figs. 12 and 15, the latter has thicker features close to the boundaries, where the error is higher, but is left mostly unchanged close to the center of the design domain, where the error goes to zero.

Validation of the numerical examples
In this Section, the previously obtained layouts (Sect. 3) are validated using COMSOL Multiphysics©. The geometry is discretized using an unstructured body-fitted mesh of triangular elements. The plane stress approximation is used along with a geometrically linear formulation. After enforcing the boundary conditions, a static step is used to compute the deformed configurations under the external load, whereas the natural frequencies are computed by solving an eigenfrequency step. The deformed configurations of the structures obtained in Sects. 3.1 and 3.2 are reported respectively in Figs. 16 and 17. The mode shapes of the structures shown in Sects. 3.3, 3.4, and 3.5 are shown in Fig. 18. In Table 6, the results obtained with OpenLSTO presented in Sect. 3 are compared with those obtained from COMSOL Multiphysics© simulations. The maximum discrepancy between COMSOL Multiphysics© and OpenLSTO for the case of compliance minimization, gripper optimization, and frequency assignment are 2.70%, 1.47%, and 0.22%, respectively.

Conclusion
In this paper, a new RLSTO formulation that can deal with manufacturing uncertainties was introduced. The method makes use of two additional structural realizationsexpanded and reduced-which can be efficiently obtained through the Level Set equation. Linear combinations of the three structures can then be used to define both the objective function and the problem sensitivities, without the need for any mapping technique.
The effectiveness of the approach was shown by applying it to two different classes of problems: static and dynamics. In particular, the compliance minimization problem and the synthesis of a gripper mechanism make it easy to compare the method to the existing ones. For the dynamic case, two variations (uniform and non-uniform etching conditions) of the problem were presented, to demonstrate the generality of this approach.
Since this method relies on the LSM and the LSI, no additional effort is required to compute the sensitivities expressions, meaning that this approach does not affect the analytical complexity of the optimal problem. Therefore, in the future, this method could be easily adapted to consider other physical phenomena, even multi-physics ones. If required by the specific problem, more than three realizations could be used to obtain alternative formulations, more suitable for cases other than the compliance minimization (stress-related problems or fluid domain optimizations). Through the boundary-perturbation technique, the additional structures can be obtained with a low computational cost, by changing the velocity field used to perturb the nominal level-set function.   Fig. 10 76.01 % 74.91 % 1.47 % (c) Natural frequencies of the results of the frequency assignment (Sects. 3.3, 3.4, and 3.5). Fig. 12 8.00 kHz 7.99 kHz 0.13 % Fig. 13 8.00 kHz 7.99 kHz 0.13 % Fig. 15 8.00 kHz 7.99 kHz 0.13 %