Generating minimal Pareto sets in multi-objective topology optimisation: an application to the wing box structural layout

Multi-objective topology optimisation problems are often tackled by compromising the cost functions according to the designer’s knowledge. Such an approach however has clear limitations and usually requires information which especially at the preliminary design stage could be unavailable. This paper proposes an alternative multi-objective approach for the generation of minimal Pareto sets in combination with density-based topology optimisation. Optimised solutions are generated integrating a recently revised method for a posteriori articulation of preferences with the Method of Moving Asymptotes. The methodology is first tested on an academic two-dimensional structure and eventually employed to optimise a full-scale aerospace structure with the support of the commercial software Altair OptiStructⓇ. For the academic benchmark, the optimised layouts with respect to static and dynamic objectives are visualised on the Pareto frontier and reported with the corresponding density distribution. Results show a progressive and consistent transition between the two extreme single-objective layouts and confirm that the minimum number of evaluations was required to fill the smart Pareto front. The multi-objective strategy is then coupled with Altair OptiStruct and used to optimise a full-scale wing box, with the clear purpose to fill a gap in multi-objective topology optimisation applied to the wing primary structure. The proposed methodology proved that it can generate efficiently non-dominated optimised configurations, at a computational cost that is mainly driven by the model complexity. This strategy is particularly indicated for the preliminary design phase, as it releases the designer from the burden to assign preferences. Furthermore, the ease of integration into a commercial design tool makes it available for industrial applications.


Introduction
Since the conception of the SIMP (Solid Isotropic Material with Penalisation) method (Bendsøe 1989;Zhou and Rozvany 1991;Mlejnek 1992), density-based topology optimisation has been successfully applied to a wide variety  (Zhu et al. 2016) and multi-physics applications (Deaton and Grandhi 2014), becoming a generally accepted approach both in the academia (Rozvany 2008) and for the development of industrial software (Choi et al. 2016). As a result, over the last decades, a considerable effort has been devoted to extend the method's capabilities to deal with the increasing complexity introduced by the new challenges of topology optimisation. Among them, the definition of novel strategies for multi-objective problems.
Unlike the case of single-objective optimisation, in the presence of many, usually conflicting, criteria, it is more appropriate to generate a set of optimal solutions which constitute the so-called Pareto set. Each of these solutions represents a trade-off between the objectives. According to Messac et al. (2003), an optimisation method can generate an effective Pareto set if it is able to (1) produce points which are evenly distributed in the design space; (2) identify all the optimal solutions (necessary condition of optimality); (3) retain only non-dominated points (sufficient condition of optimality); and finally, (4) be easily applicable. With these criteria in mind, numerous algorithms have been assessed. Metaheuristics, which are out of the scope of this work, for multi-objective structural optimisation problems, were reviewed by Zavala et al. (2014), while Marler and Arora (2004) presented a survey of continuous non-linear multi-objective methods.
In the latter, the authors introduced a classification of these methods depending on how the decision maker is required to express the relative importance or weight of each cost function, usually referred to as designer's preference. Some methods require introducing these attributes before performing the optimisation. For this reason, they are indicated as methods with a priori articulation of preferences. Both Compromise Programming (Chen et al. 1999) and Physical Programming (Messac 1996) belong to this class. A detailed analysis of each method is available at the referenced articles. However, two aspects must be retained in this context. Firstly, the ability to explore efficiently the design space cannot be achieved by consistently changing the weights. Secondly, the accuracy of the global utility function, the relationship embedding all preferences, depends on information which at the early design stage could be limited or unavailable.
In spite of these weaknesses, most of the literature concerning multi-objective problems for SIMP topology optimisation implement methods with a priori articulation of preferences. Luo et al. (2006) combined compromise programming with sequential convex programming method to optimise a missile frame under multiple load cases, considering compliance and eigenfrequency as objectives. A single optimised structure was obtained by applying a multilevel sequential approach consisting in the consecutive optimisation of two single-objective problems. On the same line, Hongwei et al. (2008) and Peng et al. (2018) used the weighted sum method to optimise respectively a vehicle frame under multiple load cases and a connection structure with respect to compliance and the first three eigenfrequencies. Zhang et al. (2018) obtained an optimised wind turbine brake pad by applying compromise programming in case of thermal and structural cost functions. Due to the high geometrical complexity associated, the authors of these works took advantage of the commercial software Altair OptiStruct. Stanford and Ifju (2009) also used compromise programming to perform SIMP topology optimisation of the flexible wing of a micro-aerial vehicle for aerodynamic objectives. Unlike the previously cited applications, the authors were able to produce a Pareto set by modifying the weights assigned to each objective. Chen et al. (2010) and Marck et al. (2012) also presented a set of optimised solutions implementing the weighted sum method to optimise the layout of a heat exchanger for structural and temperature-related criteria. Similarly, Aulig et al. (2014) built a set of compromise solutions for an automotive component under static and crash loads. A slightly different approach is proposed by Lin et al. (2010). In this case, the authors combined physical programming and convex programming methods to improve the design of compliant mechanisms with respect to displacement and strain energy. Modifying the designer's priorities, they were able to produce a set of optimised layouts.
To overcome the limitations deriving from the need to assign a relative importance to each criterion, other methods have been developed. Some of them do not require the designer to express preferences before running the optimisation. Rather, a palette of optimal solutions is generated and proposed to the decision maker which only at this point can choose among a set of best alternatives. Marler and Arora (2004) refer to these methods as a posteriori articulation of preferences or generatefirst-choose-later approaches. The normalised Normal Constraint (NC) method (Messac et al. 2003) belongs to this second group, and according to its authors it is able to produce evenly spaced Pareto frontiers , which are independent of the scale of the individual objectives. Thus far, there is no application of these methods to SIMP topology optimisation. Recently, the NC method has been used in combination with the bi-directional evolutionary structural optimisation (BESO) method. Munk et al. (2017a) proposed an updated version of the normal constraint method, modified to take advantage of the BESO ability to handle multiple constraints. Compared with the original method, this variant proved to be more efficient producing well-distributed optimal solutions and avoiding additional optimisation runs resulting in redundant points. Following a similar approach, the updated version is coupled with the Method of Moving Asymptotes (MMA) (Svanberg 1987) to perform multi-objective SIMP topology optimisation. The integration is first tested on a bidimensional benchmark and later used to optimise the internal structure of a full-scale wing. The rest of this paper is organised as follows. Section 2 explains the difference between the original NC method and the updated version, outlining the integration with the MMA algorithm. Section 3 is dedicated to the formulation of the optimisation problem and its numerical implementation with the SIMP method. In Section 4, the approach is applied to a bidimensional academic benchmark and the corresponding results are discussed. Section 5 first reviews the other works concerning the application of single-and multi-objective topology optimisations to the wing structural layout. Then, it illustrates the implementation of the multi-objective strategy, presented in this paper, with the commercial software Altair OptiStruct. The optimised wing box layouts and the corresponding Pareto frontier are reported and discussed. Concluding remarks are provided in Section 6.

The updated Smart Normal Constraint method
Recently, Munk et al. (2017a) proposed to generate the Pareto set for multi-objective BESO by implementing what they called updated Smart Normal Constraint method, abbreviated as updated-SNC or uSNC in the rest of this paper.
The normalised Normal Constraint (NCC) method introduced by Messac, Yahaya, and Mattson (2003) is a variant of the original version proposed earlier by the same authors (Ismail-Yahaya and Messac 2002). In the following years, the method underwent a series of changes, first with the introduction of the smart Pareto filter by ; later, Boyce and Mattson (2008) increased its efficiency by discarding redundant points associated with disjointed frontiers. Haddock et al. (2008) reformulated the same filter in the form of additional linear constraints through the introduction of the Practically Insignificant Trade-off (PIT) region concept. Eventually, the evolution into the smart Normal Constraint (SNC) method was due to Hancock and Mattson (2013) which refined this strategy embedding the smart Pareto filter directly into the optimisation process and avoiding simulations resulting in discarded points. The updated version suggested by Munk et al. (2017a) differs from the SNC method for the introduction of an additional normal constraint which reduces the feasible design space for each iteration, as illustrated in Fig. 1 for the case of a bi-objective problem. The generic multi-objective optimisation problem (MOP) can be formulated as: Being μ(x) the vector of objective functions, x the density design vector, m the number of objectives, and g(x) and h(x) the inequality and equality constraints, respectively. The SNC method turns the original problem (1) into a single-objective problem (SOP) subject to m − 1 additional normal constraints (2).
If n ap is the number of approximation points U defined on the utopia line, the Pareto frontier is built by solving the SOP associated with each of the approximation points. The where μ j * is the j -th anchor point, obtained as solution of the SOP with respect to the j -th objective. With reference to Fig. 1, being m = 2, the only additional normal constraint is defined on the approximation point U l currently under examination and produces the red-shaded feasible region. By contrast, the updated-SNC version modifies the search area, by introducing two normal constraints which are defined on the points adjacent to the active approximation point U l . The MOP in this case is converted into the equivalent SOP in (4), with the search area filled with grey stripes (Fig. 1).
The authors of the updated-SNC showed that the restricted feasible region avoids the issue of generating redundant points, a problem that affects the original version of the SNC, increasing the overall efficiency of the method (Munk et al. 2017a).

Integration with the MMA
The solution of the SOP (4) presented by Munk et al. (2017a) employs an iterative scheme where Lagrange multipliers are updated upon satisfaction of the additional normal constraints. This approach takes advantage of the BESO ability to satisfy simultaneously multiple constraints. For the SIMP method, a similar strategy would be possible implementing the optimality criteria (Bendsøe 1995). This heuristic scheme proved to be easily applicable and extremely efficient when a single constraint on the material resources is introduced. However, from (4), it is clear that in this case the required approach has to deal with multiple constraints. Although such an extension of the optimality criteria has been suggested by Yin and Yang (2001), other solutions, like the sequential linear programming (SLP) and the Method of Moving Asymptotes (MMA) (Svanberg 1987), are largely employed in multi-constrained topology optimisation problems due to their superior convergence qualities. In this paper, the MMA method is used in conjunction with SIMP topology optimisation. The two normal constraints already defined in (4) can be passed directly into the MMA algorithm in the form of scalar products. Furthermore, as this method requires gradient information both for the objectives and constraints, also the derivatives in (5) must be provided.
Where the subscript e indicates the element in the discretisation. In order to guarantee that all quantities are of the same order of magnitude, each objective is normalised with respect to the extreme values obtained from the anchor points previously calculated.

Problem formulation and sensitivity analysis
A bi-objective problem for minimisation of structural compliance (C) and maximisation of the fundamental frequency (f 1 ) is defined to test the updated-SNC in combination with SIMP method (in the rest of the paper indicated as uSNC-SIMP). Following the general formulation in (1), the problem is reported in (6).
With V * f the target volume fraction, v e volume of the e-th finite element, NE the number of elements in the discretisation, and x min = 0.001 the minimum density.
Following the passages illustrated in Section (2), the initial bi-objective problem is re-written as a SOP for compliance minimisation, as showed in (7).
For the bi-objective problem in (6), the derivatives in (5) can be explicitly expressed as follows: To calculate the derivatives of the two objective functions, an interpolation scheme must be defined. The basic version employed by SIMP method for stiffness and material density in the case of isotropic material is reported in (9)-(10), as also suggested by Huang et al. (2010).
Being q = 1 for the material density and p ≥ 3 for the elastic modulus, the values that guarantee an appropriate penalisation for intermediate densities. A value of p = 3 will be used throughout the paper. The superscript "0" refers to values of solid material.

Compliance sensitivity
Compliance sensitivity number can be easily derived making use of the interpolation scheme enunciated in (10) and under the hypothesis of design independent loads. In (11), it is expressed for the single element.
where u e is the element displacement vector and [K 0 e ] the element stiffness matrix for solid material. This expression is used in (7) to calculate the objective's derivative and in (8) to obtain the first component of the normal constraint.

Eigenfrequency sensitivity
Mathematical derivation of the eigenfrequency sensitivity number has been thoroughly obtained by Adhiakri (1999), in a paper which also includes the contribution of damping. The dynamic behavior of an undamped continuum structure is described by the eigenvalue problem in (12).
With [K] and [M] the global stiffness and mass matrix, ω n the angular frequency associated with the n-th eigenmode, and φ n the corresponding eigenvector. These are related to each other through the Rayleigh quotient below.
Differentiating (13) with respect to the element density x e , the general expression for the eigenfrequency sensitivity in (14) is obtained.
Assuming that the element eigenvector φ e n is normalised with respect to the mass matrix [M]. Similar to (11), derivatives of the stiffness and mass matrices depend upon the interpolation scheme used for the material density and the elastic modulus. In contrast to the compliance minimisation, optimisation involving eigenfrequencies presents two main issues. On the one hand, localised eigenmodes occur in low-density elements when their stiffness/mass ratio becomes too small. Pedersen (2000) proposed a modified interpolation scheme with a reduced penalisation factor for low-density elements. Similarly, Huang et al. (2010) suggested a novel interpolation scheme which converges to the basic SIMP when x min tends to 0. Frequency coalescence, on the other hand, occurs when adjacent eigenvalues which are originally separate converge to the same value during the optimisation process, resulting in a duplicate eigenfrequency. Bendsøe (1995) suggested to prevent it by implementing a bounded formulation of the optimisation problem. Both the first and second issues are addressed, for the interested reader, in a paper by Du and Olhoff (2007). With respect to the first issue, in this paper, the interpolation scheme proposed by Huang and reported in (15) and (16) has been adopted, with a value p = 3 as suggested by its authors.
It follows that the eigenfrequency sensitivity is defined by (14) with the derivatives calculated according to the modified interpolation scheme in (15) and (16). It is employed in (8) to calculate the second component of the normal constraint.

Sensitivity filtering
The existence of the solution to the problem in (2) is guaranteed by implementing the original version of the mesh-independent sensitivity filter (Sigmund 2007). It modifies the sensitivity of each element with respect to the neighbor elements located within a circular area defined by the radius r min . In this way, the discontinuous change occurring between the element boundaries is smoothed over the surrounding area. The filtered sensitivity for the generic objective μ and the element "e" is defined in (17).
where x e represents the pseudo-density associated to the central element and k is the index used to identify any other element within the radius considered for the filtering procedure {k| dist (e, k) ≤ r min }. It follows that the weight factor H k is calculated as in (18).
While the sensitivity filter has been proved to be an effective strategy to achieve regularisation, it must be mentioned that its purely heuristic character introduces an artificial distortion of the optimised results.

Numerical implementation
The whole software framework consists of a set of Matlab scripts which integrate SIMP topology optimisation, the MMA algorithm, and the updated-SNC version.
The approach is first applied to the two-dimensional benchmark taken from Proos et al. (2001) and also reused for comparison by Munk et al. (2017a). It consists in the simply supported structure of Fig. 2, loaded symmetrically on the top edge. To remain coherent with the cited papers, the bi-dimensional domain is discretised into 80×50 quadrilateral elements with two translational Fig. 2 Simply supported two-dimensional benchmark structure used in two previous articles (Proos et al. 2001;Munk et al. 2017a). It has a thickness t = 10 mm, and is discretised with 80 membrane elements in x and 50 elements along y degrees of freedom per node, in x and y directions. Despite its symmetry, the entire structure is solved because some modes could be non-symmetric. Isotropic material properties are assigned according to the following values, ρ = 7 · 10 −6 kg/mm 3 , E = 200 GPa, and ν = 0.3. This structure was selected as a benchmark because the optimised layout with respect to the fundamental frequency preserves material where the load is applied, ensuring a finite value of compliance. Furthermore, an analysis of the baseline (full solid) structure shows that the fundamental frequency f 1 is unimodal and equal to 600.4 Hz (a value of 591.1 Hz is obtained from OptiStruct using shell elements under the hypothesis of plane-stress).

Single-objective results
Single-objective results are presented as they identify the anchor points of the smart Pareto set generated by the implementation of the updated-SNC. It follows that they are used as extreme configurations for the multi-objective optimisation. They are obtained by solving problem (6) separately for each objective. A target volume fraction V * f = 0.7 has been imposed. In Fig. 3, the optimised layout for compliance minimisation is reported. The target volume of 70% of the original material is satisfied. A value of C = 6.3 N mm is obtained for compliance, while the fundamental frequency raised, with respect to the baseline case, to f 1 = 681.5 Hz. Numerical results are in agreement with those provided by OptiStruct, C = 6.4 N mm and f 1 = 681.2 Hz. The optimised layout for maximisation of the first eigenfrequency is given in Fig. 4. The layout converged to a compliance value C = 7.5 N mm with the fundamental frequency increased to f 1 = 715.3 Hz. Results, in this case, showed some discrepancy with respect to those provided by OptiStruct (C = 7.9 N mm and f 1 = 697.3 Hz). This difference is ascribed to the interpolation scheme employed for the calculation of

Multi-objective results
This section presents the multi-objective results, relative to the bi-dimensional benchmark in Fig. 2, obtained with the uSNC-SIMP approach proposed in this paper. The authors of the updated-SNC demonstrated how the method is superior to the SNC version in building minimal Pareto set with evenly distributed and non-redundant points. Here, the same features will be assessed against the results presented in the two works used as references (Proos et al. 2001;Munk et al. 2017a). Non-dimensional Pareto curves are compared in Fig. 5. Overall, the Pareto frontier is correclty intercepted. The uSNC method could find, in both cases, more and better distributed solutions compared with the global criterion method employed by Proos et al. (2001). Restricting the comparison to the applications of the uSNC, the results obtained in this paper show that some of the optimal points slightly differ in the frequency or compliance value. This difference can be ascribed, on the one hand, to the ability of SIMP to define a wider design space with respect to BESO, on the other, to the interpolation scheme used to optimise the dynamic objective. The smart Pareto set obtained by solving problem (7) is reported in Fig. 6. For each point, the corresponding PIT region is marked in red. As described in detail in the original paper (Haddock et al. 2008), the PIT region is defined through two parameters. The objective variation a m denotes the minimum significant change in the objective and was here set at 5%. Its curvature κ, which controls the PIT "sharpness", was set to 0.4. Both values correspond to those used in Munk et al. (2017a). With these settings, the algorithm could identify ten non-dominated points. The fact that they belong to a minimal Pareto set can be deduced The non-dimensional Pareto frontier obtained for the case of simply supported structure with the uSNC-SIMP approach (blue dots).
Results proposed for the same test case by two other authors are compared. A set of points was obtained using the global criterion method (in black), the other employing the uSNC-BESO. The utopia point in this case is located at the top-left corner by observing that no point falls inside the PIT region of an adjacent solution. In its original implementation, with the bi-directional evolutionary topology optimisation, for the same test case, the same number of optimal points was found. In addition, for the results discussed here, each optimisation run was successful in identifying a point on the Pareto frontier, confirming the ability of the updated-SNC version to fill the Pareto front with the minimum number of calculations. Figure 6 proves also that the entire design space has been explored systematically as all PIT regions are partially superimposed or in contact with each other. This is further corroborated by the smooth and consistent transition visible in the sequence of optimised layout corresponding to each optimal point. The results presented in Fig. 7 are obtained for a mesh grid of 80 × 50 elements, which corresponds to the mesh size of the works used as reference. In order to verify if the generation of the Pareto frontier can be somehow affected by the mesh size, the optimisation was repeated for a coarser mesh (40 × 25 elements) and a finer mesh (160 × 100 elements). In the case of the coarser mesh, the update SNC identified eight optimal points; the corresponding optimised layouts are reported in Fig. 8. In total, four single-objective optimisation runs resulted in dominated solutions, which suggests a difficulty for the method to compromise between the different objective values when smaller changes in the material layout are restricted by the mesh size. By contrast, the finer mesh produced eleven optimised layouts (Fig. 9); this time with no redundant calculations. The optimised layouts overall confirm the trade-off solutions already found with intermediate mesh size. A minimal amount of grey material can be observed in some of the layouts, based on the discreteness index defined in (19). Although grey elements are not desired, the SIMP algorithm produces final structures with some grey  As it shows that minor difference between the topologies produces a noticeable variation in the objectives, indicating that small design changes can result in an appreciable effect on performance. These results confirm the correct implementation of the method and lay the foundations for the following application.

Application to the wing box layout
Multi-objective optimisation attracts great interest among the aerospace community because in this field the presence of multiple conflicting objectives and stringent constraints accentuates the need of trade-off solutions. As a result, the ability to choose among a limited number of optimised designs appears to be crucial, especially at the preliminary design stage.

Topology optimisation of the wing structure
Aerospace applications of topology optimisation have been reviewed by Zhu et al. (2016), while more recently Munk et al. (2012) emphasised on feasibility of the optimised layout and its compliance with aviation regulations. Although the number of publications in this area is growing constantly, the attention here is restricted to those works addressing the optimisation of the wing primary structure, in a clear attempt to rethink the traditional spar-rib design or, at least, to achieve some improvement in terms of weight reduction and in-flight performance. The first example of topology optimisation applied to the wing structural layout (Balabanov and Haftka 1995) employed the ground structure approach to solve the minimum compliance problem for the wing of a high-speed civil transport aircraft.
Results, limited to a component-level, outlined the layout of the primary structure. On the opposite end of the resolution scale, Aage et al. (2017) obtained probably one of the most impressive results. The authors of this study performed SIMP topology optimisation of the Nasa Common Research Model (CRM) wing (Vassberg et al. 2008), achieving a level of mesh refinement without precedents for density-based topology optimisation. The interpretation suggested the presence of curved ribs in the inboard wing, a feature also observed in another similar study (Crescenti et al. 2018). In between, a variety of approaches and benchmarks have been presented. Maute and Reich (2006) integrated SIMP topology optimisation into a multidisciplinary framework to perform single-objective optimisation of quasi three-dimensional adaptive wing sections. Although some simplifications were introduced, the authors pointed out the importance of a multidisciplinary approach to account for the change in the structural response. A similar design was employed by Walker et al. (2015) to find the minimum compliance layout using Altair Optistruct. Eves et al. (2009) andOktay et al. (2014) went further, building a CFD topological framework. The former improved the design of a blended wing-body aircraft; the latter solved the minimum compliance problem for a three-dimensional rectangular wing. Munk et al. (2017b) presented a study where the dynamic stability of a plate wing first, and the CRM wing later, was optimised by maximising the largest frequency gap among the first ten modes. Optimisation with respect to dynamic stability was also investigated by Stanford and Beran (2011) for a plate wing having different volume fractions and sweep angles. Other examples of topology optimisation applied to the wing structural design adopt the Level Set method (Wang et al. 2003;Allaire et al. 2004). Gomes and Suleman (2008) maximised the aileron reversal speed of a wing torsion box using Level-Set topology optimisation and re-designing the upper skin. Brampton et al. (2012) and Dunning et al. (2014) optimised the primary structure of a three-dimensional full-scale wing with the Level-Set method, considering structural compliance as response and modifying load cases and volume fraction.
Unlike the works cited so far, only two papers relative to multi-objective topology optimisation of the wing design have been reviewed. Stanford and Ifju (2009) carried out multi-objective topology optimisation of a membrane micro air vehicle wing for aerodynamic objectives. Trade-off solutions were constructed compromising two, out of three, responses. Sleesongsom and Bureerat (2013) applied the ground structure method in conjunction with a populationbased algorithm to perform topology and size optimisation of an unswept, untapered wing box. The resulting Pareto frontier was obtained for a combination of three aeroelastic and structural objectives. Although this specific case does not rely on a compromise method, it becomes unpractical when applied to density-based topology optimisation with a significantly higher number of design variables.
The reviewed articles highlight a deficiency in multiobjective applications. While a significant effort has been made to integrate topology optimisation into a multidisciplinary framework for aerodynamic coupling, little effort has been made to integrate multi-objective algorithms. In addition, because of the difficulty to handle complex geometries, about half of the reviewed works use simplified designs, consisting in two-dimensional or tridimensional rectangular wings. In an attempt to overcome these limitations and to demonstrate the superiority of a multi-objective approach, the present work will be applied to the bi-objective optimisation of the CRM wing box.

Software framework
The Altair package (2017) is employed in combination with Matlab. HyperMesh is used for the pre-processing and OptiStruct for the structural optimisation. Different from what was done in Section 4, Matlab is only used to implement the uSNC method and to update the HyperMesh input file containing the normal constraints. Once the optimisation in OptiStruct converges, the objective vector μ(x) is updated. Depending on whether the new point is optimal or dominated, the Pareto frontier is also updated and a new approximation point selected. The two neighbouring points U l−1 and U l+1 are substituted in (20)-(21) to start the new search. The jig-shape approach (Maute and Allen 2004) is adopted for the optimisation; as a result, aerodynamics and topology optimisation are not coupled. This choice is justified as more attention here is paid to the implementation of the method. The optimisation is restricted to the bi-objective problem defined in (6) because compliance (static) and frequency (dynamic) are conflicting and vital considerations in aircraft wing design. The MOP is transformed into the single-objective problem (7) with the two additional conditions.
OptiStruct also offers the possibility of selecting among three optimisation algorithms. The Method of Feasible Design (MFD) is the default option; however, the Sequential Linear Programming (SLP) method has been chosen as it allows for the introduction of the equality constraint associated with the volume fraction, which in this case is set to V * f = 0.5.

Wing computational model
The geometry of the computational model is represented by the wing box of the CRM model. Data are available from Vassberg et al. (2008), with the exception of the spars location that is assumed at 25% and 75% of the chord for the front and rear spars, respectively. The model consists of two components. The outer surface is treated as non-design skin, coloured in blue in Fig. 10; this prevents the material from being removed, preserving the original aerodynamic shape and a sealed volume. The second component, greycoloured, is the internal region, subject to the material redistribution. The non-design skin is discretised using quadrilateral shell elements, to which was assigned a thickness property of t = 1 mm. As a result, the skin provides little contribution to the overall wing box stiffness. The entire model was discretised using 1.6·10 5 elements and 1.2 · 10 5 degrees of freedom. A length control technique is implemented in OptiStruct by defining a minimum member size which according to the Altair User Guide is between two and three times the average element size. The same isotropic material properties are assigned both to the skin and the internal volume, Young's modulus E = 73 GPa, density ρ = 2.7 · 10 −6 kg/mm 3 and Poisson coefficient ν = 0.3. The root section of the wing box is fully clamped (6dof constrained per node). The aerodynamic load is applied in the form of a pressure distribution (Fig. 11) on the upper surface only and was obtained for cruise conditions Mach = 0.85, Re = 4 · 10 7 per reference chord and a C L = 0.5. No concentrated loads are considered for Where the bar symbol indicates that both compliance and frequency are normalised with respect to their extreme values, while χ is the ratio between the first and second components of the normal vector defined in (3). The two expressions are thus associated with the corresponding responses and subcases and finally assigned to two design constraints. Regular convergence in OptiStruct is achieved when the convergence criteria are satisfied for two consecutive iterations. The convergence criteria used here are indicated as OBJTOL and CONTOL as defined by the Altair User's Guide (2017). OBJTOL checks the relative change of the objective function between two consecutive iterations, and it has been set at 10 −3 . CONTOL expresses the percentage of constraint violation and is set equal to 0.1%. A third parameter, DELTOP, sets the move limits. It is used by OptiStruct to define the upper and lower bounds of the design variable changes during each iteration. Smaller values usually lead to smoother convergence but a larger number of iterations. DELTOP has been set at 0.1 for this model.

Wing box: multi-objective results
As for the academic two-dimensional structure, the two anchor point solutions of the SOPs for compliance minimisation and fundamental frequency maximisation provide essential information to normalise the objectives and to identify the boundaries of the entire design space. The minimum compliance problem produced the optimised layout reported in Fig. 16a, with C min = 2.964 · 10 7 N mm and f 1 = 2.404 Hz. The maximum fundamental frequency resulted instead in the structure of Fig. 12, with C = 1.932 · 10 8 N mm and a fundamental frequency f max = 3.372 Hz. This configuration was coherent with those found in two previous studies (Stanford and Beran 2011;Munk et al. 2017b), Starting from these two anchor points, the other optimised solutions found on the Pareto frontier resulted in a compliance value of C ∈ (C min , 10 8 N mm), but inconsistent frequency values f > f max . As a result, the anchor point for frequency maximisation became a dominated solution. In order to explain this incongruence, the frequency anchor point was further investigated. Its compliance convergence history (blue curve in Fig. 13) is compared with that of two other simulations. First, the frequency maximisation problem was solved under an additional constraint on the compliance value (C < 10 8 N mm) producing the red curve in Fig. 13. Such a value was chosen as the threshold because all the optimised solutions previously found stayed below this limit. Although the starting point is the same as the blue curve (material density is initialised by default at 0.6 in both cases), the optimisation produced a new anchor point with a frequency f 1 = 4.975 Hz, which is higher than the previous one. The corresponding layout is reported in Fig. 16h. For the third comparison (green curve), instead of introducing an empirical threshold for compliance, the material density was initialised at 1.0 . As a result, the optimisation started from a lower compliance value but without constraints it was Fig. 13 Compliance convergence history for the frequency maximisation. Baseline solution (blue), constrained optimisation C < 10 8 N mm (red), constrained optimisation with move limit equal to 0.01 (cyan), and baseline problem with material density initialised at 1.0 (green) potentially free to explore the design space. As expected, convergence was slower but attained the same value as the red curve. This comparison pointed out as by starting from a less stiff structure (MatInit= 0.6) led the optimisation to a local minimum of compliance, preventing the algorithm from further increasing the fundamental frequency. Forcing the solution towards lower compliance values or simply moving the starting point far from this area allowed the optimiser to find a non-dominated solution. Although all curves attained regular convergence, they did not show the expected flatness. To justify this fact, a fourth case (cyan) has been reported in Fig. 13. It reproduces the same case as the red curve (compliance-constrained) with a reduced DELTOP parameter set at 0.01. The reduced move limit clearly smoothes the convergence achieving the same compliance value, but at the expenses of the iteration number, 54 against the 14 of the red curve. Despite its smoothness, the cyan curve presents a jump in the final part, which is also recognisable in the other cases, due to the continuation strategy automatically implemented by OptiStruct to achieve better discrete results. Following these considerations, the new anchor point for frequency maximisation was employed for the generation of the minimal optimal set, along with the minimum compliance solution whose internal features are visible in Fig. 14.
The resulting Pareto frontier is reported in Fig. 15, where each optimal solution is associated to the corresponding PIT region (a m = 5% and κ = 0.4) and the "datum" point, full wing box with volume fraction set at 0.5, is visualised for reference. Under the guide of the uSNC, OptiStruct was able to identify eight optimal solutions and only in two cases the optimisation process resulted in a dominated layout. This produced the two gaps between points (c)-(d) and points (e)-(f) on the Pareto frontier. In both cases, the dominated points were slightly in contact with the PIT region of the neighbour point, suggesting that an optimal solution could still be achieved by refining the mesh or, alternatively, by setting a tighter convergence threshold for the normal constraints. This explanation is also supported by looking at the material distributions in Fig. 16, where only density values x ≥ 0.5 are reported to highlight the elements that contribute the most to structural stiffness. For the minimum compliance layout (a), material is distributed along the span to counteract the bending effect induced by the aerodynamic load. On the contrary, the maximum eigenfrequency design (h) tends to remove material from the tip to avoid the "lumped mass" effect which otherwise would decrease the fundamental frequency. Interestingly, the progression between these two configurations follows a smooth transition which also guarantees a connected structure, in contrast to that shown in Fig. 12. The absence of any rapid change in the structural layout in correspondence of the gaps on the Pareto frontier discards the possibility of a discontinuity in the distribution of optimal solutions. The multi-objective optimisation is performed on a desktop computer having 4 cores and the CPU time for the generation of a new optimised layout in OptiStruct takes between 6 and 8 min. For a bi-objective problem, as the one solved in this study, the additional constraints do not seem to affect the CPU time, which is mainly driven by the model complexity. One could expect, however, that in the case of many objectives the benefits of the update may be outweighed by the increase in computation required by the higher number of constraints.

Conclusions
This work presented a strategy to deal with multiple objectives in SIMP topology optimisation, offering an alternative to the most common compromise programming, which is employed in almost the totality of the reviewed literature. Starting from a recently revised version of the Smart Normal Constraint method, the authors adapted this new version to the Method of Moving Asymptotes, which is universally recognised as an efficient and reliable optimisation algorithm for density-based topology optimisation. The optimisation strategy was first tested on a two-dimensional academic test case for a bi-objective problem, with a twofold purpose: to define a reference for the industrial application and, above all, to obtain a direct comparison with two previous works. The nondimensional Pareto frontier showed a strong agreement with the results already available in the literature and full coherence in terms of the method's efficiency. It must be noticed however that the high value of volume fraction used for the 2D benchmark, to match the previous results, covers the possible occurrence of multiple eigenfrequencies, so further tests are necessary to provide exhaustive conclusions on the method. Once the approach was proved for the 2D structure, the same method was used on a real aerospace application, consisting of the wing box of the CRM model. In order to cope with the increased geometric complexity, the SIMP module in Matlab was replaced by the commercial software OptiStruct, showing how the methodology can be promptly integrated into a professional engineering tool for industrial design. Based on the information provided by the uSNC-SIMP multi-objective approach, the solution corresponding to a local minimum was discarded and the search focused on an area where more points were found. While the set of optimised results showed a consistent evolution of the layouts and provided a new insight into the design of the wing box internal structure, the absence of detailed features in the internal optimised layout suggests the need for the use of a finer mesh, here limited by computational resources.