The study presents deflation constraints that enable a systematic exploration of the design space during the design of composite structures. By incorporating the deflation constraints, gradient-based optimizersbecomeabletofindmultiplelocaloptimaoverthedesignspace.Thestudypresentstheidea behinddeflationusingasimplesinefunction,whereallrootswithinanintervalcanbesystematically found. Next, the novel deflation constraints are presented: hypersphere, hypercube and hypercuboid; consisting of a combination of Gaussian and sigmoid functions. As a test case, the developed constraints are applied to the optimization of a double-cosine function, where all the 13 minima points could be found with 24 deflation constraints. It is shown that a new optimum is encountered after each deflation constraint is added, with the optimization subsequently re-started from the same initial point, or resumed from the last found minimum, being the latter the recommended approach. Thenewdeflationconstraintsarethenusedinheuristic-baseddirectsearchmethods,whereagenetic algorithm optimizer is able to find new optimum individuals for straight-fiber composites. Lastly, variable-stiffness composites were designed with the deflation constraints applied to the multimodal optimization problem of recovering fiber orientations from a set of optimum lamination parameters. 1. Introduction Structuraloptimizationisafairlyvastsubjectinitselfthatis tightly associated with the design of lightweight structures. Composite materials have been a major enabler of these designs, owing to their superior specific material properties ascomparedtometals,andtotheirsignificantlylargerdesign space enabled by their anisotropy. Both aeronautical and space industries have been continuously developing new conceptsforcompositestopursuethesebenefits,andFigure 1showstheevolutionoftheuseofcompositesinaeronautics overthelast60decades.NotethatBoeing’s787Dreamliner was the first large commercial aircraft have 50% of the structural weight in composites (1). 1.1. Types of design space and optimization composites Ghiasietal.reviewedindetaildifferentcompositeoptimiza- tion methods used for constant (3) and variable (4) stiffness laminates, classifying the techniques into four categories, being in order of relevance: gradient-based, direct search and heuristic, specialized and hybrid methods. Gradient- based methods can find a local minimum generally faster than all the other methods, limited to contiguous problems with first or second derivatives (3). For problems that are largeindimension,withusuallymorethan30variables(5), gradient-basedmethodsmaybetheonlyviableoptimization method (6). However, gradient-based optimizers can find the global optimum with high accuracy only in unimodal designspaces(7;8),whereasinmultimodalspacesthefinal solution largely depends on the initial point (3). A convex design space is a subset of a unimodal design space, as illustrated in Figure 2, and convexity requires that all line ∗Corresponding author S.G.P.Castro@tudelft.nl (S.G.P. Castro) ORCID(s):0000-0001-9711-0991 (S.G.P. Castro) 1The two authors had an equal contribution to this manuscript.segments connecting any two points in the function are above the function and do not intersect it. Hence, not all unimodalfunctionsareconvex,butallmultimodalfunctions are non-convex (6). As the design space becomes increas- ingly multimodal or non-convex, finding a global minimum becomesmoredifficultandwouldrequirespecialtechniques that can evaluate multiple regions of the design space to findthepointthatbestminimizesthefunction.Theseglobal search schemes involve different sampling or starting point approaches (9; 10) to increase the likelihood to converge to a global optimum, despite it is never guaranteed that a global optimum will be found, not even that the global search will converge onto a solution that is different from those previously found (6). In commercial software such as Altair’s Optistruct, such global searches can be activated with the "dglobal" user input (11). Directsearchmethodsdonotrequirederivatives,andcan be more appropriately applied for composite lay-up design, handling a mixture of continuous and discrete variables. They can find the global optimum of multimodal objective functions, although with a significantly lower rate of con- vergence when compared to gradient-based methods (3). Genetic algorithms (GA) has been the most popular class of direct search method, according to Ghiasi et al. (3), with simulated annealing ranking second. Haftka and his group have pioneered the application of GA in composite design (12; 13; 14). António (15) proposed a hierarchical GA for multimodaloptimizationofhybridcompositeswithmultiple solutions.Thesunfloweralgorithmisanotherrelevantglobal optimizeroriginallyappliedtoinverseproblemofstructural damage detection in laminated composite plates, proposed by Gomes (16). Variable stiffeness (VS) laminated composites can fur- ther increase the design space and hence the potential for better structural performance of composites (17), manufac- turable by automated fiber placement (18) and continuous SS Bangera, SGP Castro Preprint Page 1 of 16 Deflation constraints for global optimization of composite structures 1950 1970 2000 2010 2020GFRP and balsawood - Sidewalls - Bag racks - GalleysAircra� interior Boeing 737/777: - Composite ﬂoor - Horizontal/ver�cal stabilizer - Empennage1stall-composite planes CFRP: - Fuselage - Wing - 50% of compositesAirbus A350: 53% of composites Phönix FS-24Long EZBusiness jets - Variable angle tow high-aspect ra�o Wing: NASA - lighter, more structurally eﬃcient and more ﬂexibleBoeing 707 Boeing 787 Figure 1: Utilization of fiber-reinforced composites in aircraft over the last decades, from Wang et al. (2). twoshearing(19).However,thelargerdesignspacecreated byVSlaminatescreatesnewchallengesfortheoptimization, due to highly non-convex design spaces associated with the variable fiber angles (20). Figure 2: Unimodal and multimodal functions, from Martins & Ning (6). Even in the design of constant stiffness laminated com- posites, the use of ply-angles as design variables is often associated with non-convenxity (21). Miki (22) proposed the use of lamination parameters, introduced by Tsai & Pagano in 1968 (23), as design variables to render a convex optimization. Scardaoni & Montemurro (24) proved non- convexity of the feasible domain for both anisotropic and orthotropic-membrane laminates, even when those are pa- rameterizedusinglaminationparameters.However,approx- imations that lead to convex feasible domains have been proposed, such as the one for the in-plane and flexural lamination parameters developed by Fukunaga & Sekine (25; 26). Even when the optimization of the laminationparameters is performed in a convex design space, the next step of retrieving the fiber angles involves a multimodal optimizationthatisubiquitouslydoneusingGAalgorithms. Inthepresentstudy,itisshownthatthisretrievalcanbedone with gradient-based optimizers when deflation constraints are used. 1.2. The deflation method 1.2.1. Deflation of scalar-valued functions The first instance of the deflation methods were applied for systematic root finding in nonlinear functions (27), provid- ing an easy visualization of the concept. Assume 𝑝(𝑥)to be a scalar-valued nonlinear function having 𝑛multiple roots 𝑥1,𝑥2,...,𝑥𝑛,witheachrootfoundusinganiterativemethod such as the Newton-Raphson (28). After having evaluated the initial root of 𝑝(𝑥), more roots can be systematically found by considering the following deflated function (29): 𝑞(𝑥)=𝑝(𝑥) 𝑛∏ 𝑖=1(𝑥−𝑥𝑖), (1) for which the already obtained roots can be effectively removed by the multiplicative term in the denominator. Consider the sine function 𝑝(𝑥) =𝑠𝑖𝑛(𝜋𝑥)with 7 roots within−𝜋≤𝑥≤+𝜋. Figure 3 depicts 𝑝(𝑥)and the two newfunctions 𝑞1(𝑥)and𝑞2(𝑥),createdaftertwoconsecutive deflationsteps,respectivelyafterfindingtherootsat 𝑥=−2 and𝑥=+1, given as: 𝑞1(𝑥) =sin(𝜋𝑥) 𝑥−(−2) 𝑞2(𝑥) =sin(𝜋𝑥) (𝑥−(−2))(𝑥−(+1))(2) SS Bangera, SGP Castro Preprint Page 2 of 16 Deflation constraints for global optimization of composite structures Figure 3 shows that the deflated functions have the spe- cificrootscompletelyremoved,i.e.thedeflatedfunctionsno longercrosszeroatthosepoints,whilekeepingalltheother roots unchanged. This property of the deflation technique can be exploited to remove selected minima points from an arbitrary objective function, as demonstrated later in the present study. Figure 3: Deflated sine function obtained using Equation 1 1.2.2. Deflation applied to nonlinear equations The idea behind polynomial deflation was generalized to partial differential equations (PDE) by Farrell et al. (29), extending the concept of deflation matrix 𝑀and deflation operator𝑚thathadbeenpreviouslyintroducedbyBrown& Gearhart (30), also with the intent of solving PDEs. Farrell provedthatconvergencetodifferentsolutionsfromthesame starting point was achieved after deflation. Consider the residuals for a system of 𝑛nonlinear equations 𝐹(𝑥)such that: 𝐹(𝑥)=0 (3) After solving the system of equations once, solution 𝑥1 is obtained, such that a deflation matrix can be calculated using a deflation operator: 𝑀(𝑥;𝑥1)=𝑚(𝑥;𝑥1) 𝑚(𝑥;𝑥1)=||𝑥−𝑥1||−𝑝+𝜎(4) where istheidentitymatrix( 𝑛×𝑛);𝑝isthepolestrengthor powerthatdictatestherateatwhichthefunctionapproaches infinity;and 𝜎isanoffsetparameterthatisreachedwhenthe norm distance ||𝑥−𝑥1||→∞(31). The deflated nonlinear system of equations 𝐺(𝑥)becomes (31): 𝐺(𝑥)≡𝑀(𝑥;𝑥1)𝐹(𝑥)=0 (5) The system in Equation 5 can be solved to obtain a new solution𝑥2. Here𝐺satisfies the following properties (31):1. The two system of equations, 𝐹(𝑥)=0and𝐺(𝑥)=0 both have the same solutions for all 𝑥≠𝑥1. 2. Withtheknownsolution 𝑥1,𝐺willnotconvergeagain to𝑥1under the assumption lim𝑥→𝑥1‖𝐺(𝑥)‖>0 After having found ̃𝐾multiple solutions, the method in (29)proposesmultiplyingthedeflationmatricesandsolving the following set of equations: ̃𝐾∏ 𝑘=1𝑀(𝑥;𝑥𝑘)𝐹(𝑥)=0. (6) Farrelletal.inarecentstudy(32)havefurtherextended this method for semi-smooth equations. Deflation does not guarantee that all the solutions to a problem are found; however, it provides a systematic method to explore a pro- gressively larger number of nonlinear solutions. 1.2.3. Deflation in global search Inthecontextofsizingoptimization,theideatoapplydefla- tion as a new constraint has been independently developed by the authors, with the only literature reference using a similar approach for topology optimization being the non- peerreviewedpre-printbyTarek&Huang(33).Considerthe nonlinear programming (NLP) problem for which multiple solutions need to be evaluated: minimize 𝑥∈𝑛𝑓(𝑥) subject to𝑐(𝑥)=0 𝑙≤𝑥≤𝑢(7) where𝑙and𝑢are respectively the lower and upper bounds of the design variable 𝑥. This constrained minimization problem can be converted into an unconstrained problem by expressing the objective function as the Lagrangian, expressed as (33): (𝑥,𝜆,𝑧+,𝑧−)=𝑓(𝑥)+𝑐(𝑥)𝑇𝜆+(𝑥−𝑢)𝑇𝑧+−(𝑥−𝑙)𝑇𝑧− (8) where𝜆∈𝑚is the Lagrange multiplier vector associ- ated with the equality constraints, 𝑧−,𝑧+∈𝑛 +is the Lagrange multiplier vectors associated with the lower and upper bound. The Karush-Kuhn-Tucker (KKT) conditions are satisfied when (33): ∇𝑥(𝑥,𝜆,𝑧+,𝑧−)=0 𝑐(𝑥)=0 𝑙≤𝑥≤𝑢 𝑧+≥0 𝑧−≥0 (𝑥−𝑢)𝑇𝑧+=0 (𝑥−𝑙)𝑇𝑧−=0(9) A point𝑥that satisfies these conditions is therefore known as a KKT point, and finding multiple solutions to SS Bangera, SGP Castro Preprint Page 3 of 16 Deflation constraints for global optimization of composite structures the NLP consists of finding multiple KKT points. In the work of Tarek & Huang (33) two versions of the deflation constraint were proposed. First, an additional variable 𝑦is introduced into the optimization problem along with the deflation constraint 𝑚. With𝑦≥0, the deflation constraint is given as (33): 𝑚(𝑥;𝑥1)=||𝑥−𝑥1||−𝑝+𝜎≤𝑦 (10) where𝑥1isthefirstoptimumpointthatisfoundpreviously, 𝑝 isapower(usuallyvariesbetween2-4),andanoffsetterm 𝜎. After this constraint becomes active, it is considered that if thenewfoundoptimumpoint 𝑥2isassociatedwithvariable 𝑦such that𝑦is finite and has a value in exact arithmetic, then𝑥1≠𝑥2(sincelim𝒙→𝒙1𝑚(𝑥;𝑥1) = ∞). It can further beprovedthateveryKKTpointofthenewNLPwithafinite value of𝑦is a KKT point of the original formulation (33). A proof of how ( 𝑥∗,𝑦∗) would be a regular KKT point for the modified problem after the introduction of the deflation constraint is provided in Appendix A. If it is assumed that 𝑚 >0and𝑦 > 𝜎, the deflation constraint would be equivalent to the distance constraint as follows (33): ||𝑥−𝑥1||𝑝≥𝑧 (11) where𝑧= 1∕(𝑦−𝜎). Thus, it becomes clear that the de- flation constraint essentially puts a constraint on the known solutions.Ifatanypoint,theoptimizerapproaches 𝑧=0or 𝑦=∞,thenthedeflationoperatorisnotdistancingthenew points from the known optimum points sufficiently. Hence, in this case, the main hyperparameters that can be varied to ensure that this distancing occurs are the offset 𝜎and the power𝑝, as shown in Equation 10. The first deflation constraint defined in Equation 10 can be modified by omitting the additional variable 𝑦, and replacingitwithalargefiniteconstant 𝑀,suchthatasecond deflation constraint is obtained as (33): 𝑚(𝑥;𝑥1)=||𝑥−𝑥1||−𝑝+𝜎≤𝑀 (12) Because there are no additional variables to this second deflation constraint, a proof similar to what is given in Appendix A can be used to conclude that the newly found point𝑥∗is a regular KKT point to Equation 7. So far, Equation 10 and Equation 12 have been evaluated with a single known solution, whereas for ̃𝐾known solutions the deflationconstraintcanbeexpressedasasummation,given by (33): ̃𝐾∑ 𝑘=1𝑚(𝑥;𝑥𝑘)≤𝑦 (13)1.3. Contributions of the present work Giventhechallengesrelatedtoperformingsystematicglobal search in multimodal design spaces, the present study pro- poses a novel deflation constraint that is included in the optimizationproblembyminimizingtheLagrangianinstead of the main objective function, thus being compatible with any algorithm that allows constrained optimization. The methodology of Tarek & Huang (33), similarly to the methodology initially attempted by Bangera (34), does notperformadequatelyinalargerdesignspaceafterthefirst few deflation constraints are added. This lack of robustness ismainlyduetotheirdeflationconstraintnotbeingdifferen- tiableatthedeflatedpoints,andduetotherelativelyarbitrary offsetparameter 𝜎.Thedeflationconstrainthereinproposed addressestheseshortcomings,enablingarobustexploration of larger design spaces by guaranteeing the discovery of new distinct minima when the optimization is re-started or resumed, after a new deflation constraint is added. This new constraint is applied without altering the optimization algorithm, making it suitable to both gradient-based and gradient-free optimization methods. 2. Methodology 2.1. Interior-point optimization algorithm Interior-point algorithms stems from interior penalty meth- ods that associate a penalization term with the constraints, being a key difference the fact that the constraints are not directly penalised. In interior-point, the penalization acts on slack variables that are added to the constraints, with the penalty term increasing as the optimizer moves towards theboundaryoftheconstraineddomain.Theformulationis written as (6): minimize 𝑥,𝑠𝑓(𝑥)−𝜇𝑏∑𝑛𝑔 𝑗=1ln𝑠𝑗 subject to ℎ(𝑥)=0 𝑔(𝑥)+𝑠=0(14) where𝜇𝑏is a barrier parameter. The resulting formulation turns the inequality constraint into an equality constraint with the addition of the slack variables 𝑠. Newton’s method can be applied to solve the KKT sys- tem of equations in Equation 14, where the logarithm term associated with the slack variables is only defined for pos- itive𝑠values, acting as a barrier for negative 𝑠values (see Figure 4). Owing to the positive 𝑠values,𝑔(𝑥∗)<0at the solution, hence satisfying the inequality constraints (6). The interior-point formulation of Equation 14 is equiv- alent to the original constrained problem when 𝜇𝑏→0. Hence, a sequence of solutions needs to be obtained such that𝜇𝑏→0. A constrained problem can be reformulated as an unconstrained problem by utilizing the Lagrangian function along with associated Lagrange variables, or La- grangemultipliers(6).TheLagrangianfortheinterior-point optimization problem can be written as (6): (𝑥,𝜆,𝜎,𝑠)=𝑓(𝑥)+𝜇𝑏𝑒⊤ln𝑠+ℎ(𝑥)⊤𝜆+(𝑔(𝑥)+𝑠)⊤𝜎, (15) SS Bangera, SGP Castro Preprint Page 4 of 16 Deflation constraints for global optimization of composite structures Figure4: PlotofinversebarrierandLogarithmicbarrierpenalty functions (6) where,ln𝑠is a vector having 𝑛𝑔components being the logarithms of each value of 𝑠;𝑒= [1,...,1]is a vector of length𝑛𝑔containing ones, introduced to express the sum in vector form; 𝜆is the Lagrange variable vector associated withtheequalityconstraints;and 𝜎istheLagrangevariable associated with the inequality constraints. The KKT condi- tionscanbederivedbytakingthederivativeswithrespectto 𝑥,𝜆,𝜎and𝑠, leading to: ∇𝑓(𝑥)+𝐽ℎ(𝑥)⊤𝜆+𝐽𝑔(𝑥)⊤𝜎=0 ℎ=0 𝑔+𝑠=0 −𝜇𝑏𝑒+𝑆𝜎=0(16) where𝑆is a diagonal matrix with its diagonal values given by the slack variable vector. With this set of residual equa- tions,Newton’smethodcanbeapplied.TakingtheJacobian oftheequationsinEquation16,thefollowinglinearsystem is obtained: ⎡ ⎢ ⎢ ⎢⎣𝐻𝐿(𝑥)𝐽ℎ(𝑥)⊤𝐽𝑔(𝑥)⊤0 𝐽ℎ(𝑥) 0 0 0 𝐽𝑔(𝑥) 0 0  0 0 𝑆Σ⎤ ⎥ ⎥ ⎥⎦⎡ ⎢ ⎢ ⎢⎣𝑝𝑥 𝑝𝜆 𝑝𝜎 𝑝𝑠⎤ ⎥ ⎥ ⎥⎦=−⎡ ⎢ ⎢ ⎢⎣∇𝑥(𝑥,𝜆,𝜎) ℎ(𝑥) 𝑔(𝑥)+𝑠 𝑆𝜎−𝜇𝑏𝑒⎤ ⎥ ⎥ ⎥⎦(17) whereΣis a diagonal matrix whose entries are given by the values of vector 𝜎; and is the identity matrix. For numerical efficiency, the system can be made symmetric after multiplying the last rows in Equation 17 by 𝑆−1, which can be calculated as 𝑆−1 𝑘𝑘= 1∕𝑠𝑘. This results in the following symmetric linear system that can be directly solved: ⎡ ⎢ ⎢ ⎢⎣𝐻𝐿(𝑥)𝐽ℎ(𝑥)⊤𝐽𝑔(𝑥)⊤0 𝐽ℎ(𝑥) 0 0 0 𝐽𝑔(𝑥) 0 0  0 0 𝑆−1Σ⎤ ⎥ ⎥ ⎥⎦⎡ ⎢ ⎢ ⎢⎣𝑝𝑥 𝑝𝜆 𝑝𝜎 𝑝𝑠⎤ ⎥ ⎥ ⎥⎦=−⎡ ⎢ ⎢ ⎢⎣∇𝑥(𝑥,𝜆,𝜎) ℎ(𝑥) 𝑔(𝑥)+𝑠 𝜎−𝜇𝑏𝑆−1𝑒⎤ ⎥ ⎥ ⎥⎦(18) The optimization algorithm developed in the present study fulfils the requirement of handling versatile inputs, beingapplicabletoconstrainedandunconstrainedoptimiza- tion problems, while handling both equality and inequalityconstraints. Therefore, generally grouped in four types of optimization problems: 1. Fullyunconstrained,whichcanbesolvedbyapplying Newton’s method to the simple unconstrained condi- tions of optimality. 2. Only with equality constraints that can be solved by employing the method of sequential quadratic pro- gramming,whichinvolvesasimilarsetofequationsas theinterior-pointmethod(Equation18),butomitsthe equations associated with the inequality constraints. 3. Onlywithinequalityconstraints,whichcanbesolved using the interior-point method (Equation 18), but omitting the equations associated with the equality constraints. 4. With both equality and inequality constraints, which can be solved using the interior-point method. 2.2. Deflation constraint Consider the deflation constraint of Equation 4 evaluated at the design point 𝑥=𝑥∗: 𝑚(𝑥;𝑥∗)=||𝑥−𝑥∗||−𝑝+𝜎 (19) Figure 5 plots Equation 19 using 𝑥∗= 0and𝑝= 2. Here, the offset parameter 𝜎influences the ability to overcome saddle points created by the deflation constraint when us- ing gradient-based optimizers, but the use of 𝜎keeps the influence of the deflation constraint even far from 𝑥=𝑥∗. When||𝑥−𝑥∗||→0, the constraint behaves as 𝑚(𝑥;𝑥∗)→ ∞, which is a discontinuity of the deflation constraint that prevents its robust use in gradient-based optimizers. 4 3 2 1 0 1 2 3 4 x01234567m(x) x* Figure 5: Plot of deflation function given by Equation 19, also proposed by Tarek & Huang (33) Aiming to reach a differentiable deflation constraint function that has no influence outside the proximity of 𝑥=𝑥∗, the present study proposes the use of the following SS Bangera, SGP Castro Preprint Page 5 of 16 Deflation constraints for global optimization of composite structures 2 1 0 1 2 x012345m(x) x*a Magnified (a) Gaussian deflation function (Equation 20) Magnified (b)Gaussiandeflationfunctionwiththesigmoidterms(Equation23) Figure 6: Comparison of the proposed Gaussian deflation with and without the sigmoid terms. Gaussianfunctionasoneofthemaincomponentsofthenew deflation function: 𝑚(𝑥;𝑥∗)=𝑎𝑒−(𝑥−𝑥∗)2 2𝑐2 (20) where𝑎is the maximum amplitude; and 𝑐is the standard deviation.Itisproposedtoutilisethetotalspanofthedesired deflation divided by 6 as the standard deviation, such that 99.7% of the Gaussian distribution is covered. Figure 6a depicts the Gaussian distribution using 𝑎= 5and𝑐= 1∕3. For both the original deflation function, depicted in Figure 5, and the proposed Gaussian deflation function depicted in Figure 6a, it can be observed that the deflation constraint does not terminate to zero beyond the deflationspan −1≤(𝑥−𝑥∗)≤+1,noritbecomesnegative, making it inappropriate to be used as a new equality or inequalityconstraint.Thisisbecausethedeflationconstraint function still has a residual value 𝑚(𝑥;𝑥∗)>0that ex- tends way beyond the desired deflation span. Therefore, it is suggested to use a differentiable decay function, being a sigmoid, that multiplies to the Gaussian function with the aim to terminate the aforementioned residual value entirely beyond the desired deflation span. Two one-sided sigmoid functions are used, given as (35): 𝑆𝑖𝑔𝑚𝑜𝑖𝑑𝐿(𝑥)=1 1+𝑒−𝐾(𝑥−𝑥∗+𝑏)(21) 𝑆𝑖𝑔𝑚𝑜𝑖𝑑𝑅(𝑥)=1 1+𝑒𝐾(𝑥−𝑥∗−𝑏)(22)where𝑏ishalfoftherequireddeflationspan; 𝐾isaparam- eter that controls the steepness of how the sigmoid function goes from1to0at the required point of termination, i.e. at 𝑥−𝑥∗=𝑏.Withthetwoone-sidedsigmoidfunctions,Equa- tion 21 and Equation 22, it is possible to make 𝑚(𝑥;𝑥∗)→ 0respectively for negatives and positives values of ||𝑥− 𝑥∗||; with the−𝑏term becoming +𝑏and𝐾becoming−𝐾, fromEquation22toEquation21.Hence,withinthedesired deflation span, the product of the two sigmoid functions results in1, whereas outside the span this product results in 0. Figure 6b shows how the Gaussian function is changed bythemultiplyingsigmoidterms,resultinginthefollowing deflation function: 𝑚(𝑥;𝑥∗)=𝑎𝑒−(𝑥−𝑥∗)2 2𝑐2⋅𝑆𝑖𝑔𝑚𝑜𝑖𝑑𝐿(𝑥)⋅𝑆𝑖𝑔𝑚𝑜𝑖𝑑𝑅(𝑥)(23) where𝑐=𝑏∕3. In Figure 6b, 𝑏= 1and𝐾= 100and the resulting deflation function terminates to zero beyond ||𝑥−𝑥∗||= 1. Noticeably, the deflation function value does not become zero immediately after crossing 1, but this can be adjusted using a higher value of 𝐾for a steeper termination, or a parameter 𝑏that is slightly smaller than half of the required deflation span. Note that, the associated Gaussianandsigmoidfunctionsincreasewiththenumberof dimensions being deflated in an optimization problem. Three deflation constraint methods are herein proposed: hypersphere,hypercubeandhypercuboid;beingthelasttwo basedonEquation23,withtheiradvantagesanddrawbacks discussed in the following. SS Bangera, SGP Castro Preprint Page 6 of 16 Deflation constraints for global optimization of composite structures 2.2.1. Hypersphere approach The hypersphere or n-sphere approach is a distance-based approach in which the deflation span region is defined by a hypersphere. The main advantage is that only one scalar parameter controls the size of the deflated region, i.e. the radius of action 𝑟that is the same over all dimensions undergoing deflation within the optimization problem. The distance constraint is given as: ||𝑥−𝑥∗||≤𝑟 (24) which whether the current design point 𝑥falls within the requireddeflationspan.Followedbythat,thedeflationcon- straint of Equation 20 can be directly applied, without re- quiringthesigmoiddecayfunction.Deflationisnotapplied when||𝑥−𝑥∗||>𝑟. 2.2.2. Hypercube approach Here, a hypercube is used to represent the required deflated span, where each side of the hypercube has the same length represented by the scalar 𝑟. The distance constraint is given as: 𝑥𝑖−𝑥∗ 𝑖≤𝑟 (25) where𝑖=1,2,…,𝑛;with𝑛beingthenumberofdimensions oftheproblem.Foran 𝑛−dimensionalproblem,thedeflation function becomes: 𝑚(𝑥;𝑥∗)=𝑛∏ 𝑖=1( 𝑎𝑒−(𝑥𝑖−𝑥∗ 𝑖)2 2𝑐2⋅𝑆𝑖𝑔𝑚𝑜𝑖𝑑𝐿(𝑥𝑖) ⋅𝑆𝑖𝑔𝑚𝑜𝑖𝑑𝑅(𝑥𝑖))(26) Figure7aplotsthehypercubedeflationfunctionofEqua- tion 26 for a 2-dimensional problem, showing the even spread across all dimensions, here using a value of 𝑟=4.5. 2.2.3. Hypercuboid approach When the deflated region is represented by a hypercuboid, a different length is attributed to each side of the cuboid, being the number of sides equal to the number of dimen- sions being deflated within the optimization problem. The 𝑛−dimensional vector of deflation spans is 𝑟= {𝑟1,𝑟2,…,𝑟𝑖,…,𝑟𝑛}; such that the distance constraint becomes: 𝑥𝑖−𝑥∗ 𝑖≤𝑟𝑖 (27) where𝑖= 1,2,…,𝑛. The expression for the deflation function becomes: 𝑚(𝑥;𝑥∗)=𝑛∏ 𝑖=1( 𝑎𝑒−(𝑥𝑖−𝑥∗ 𝑖)2 2𝑐2 𝑖⋅1 1+𝑒𝐾(𝑥𝑖−𝑥∗ 𝑖−𝑟𝑖) ⋅1 1+𝑒−𝐾(𝑥𝑖−𝑥∗ 𝑖+𝑟𝑖))(28)with𝑐𝑖=𝑟𝑖∕3for𝑖= 1,2,…,𝑛, being compatible with thedeflationspanvector 𝑟.Figure7billustratesEquation28 for a2−dimensional problem, where deflation constraint is unevenacrossdifferentdimensions,hereusing 𝑟={2,4.5}. 2.2.4. Advantages and drawbacks The idea behind having a deflation constraint is to have a simple approach that can be applied to any optimization method intuitively seamlessly to find different solutions in a multimodal design space. The three proposed deflation constraints can be applied to any optimization problem; however, there are certain cases in which some perform better than the others. Thehypersphereandhypercubeapproachhaveonlyone distanceparameterinthedeflationproblem,i.e.thescalar 𝑟, beingthereforeadequatewhenalldesignvariablesshouldbe equally deflated from existing minima points. Even though it is simpler to have a single parameter to control the size of the deflated region, it increases the chances of masking local minima throughout the design space. This scenario is shown in Figure 8, which could happen in any of the threeproposeddeflationconstraints.Thefiguredisplaystwo local minima in green that were masked due to 𝑟being too large. Moreover, Figure 8 illustrates that at the point of terminationofthedeflationconstraint,theobjectivefunction intersectsthedeflatedregioncreatingfalse-minimadepicted in red. These false-minima become possible minima to be discoveredbytheoptimizerinsubsequentsearchesafterthe deflation constraint is applied, and can be overcome when further deflation constraints are added. Masking of minima pointsshouldbealwaysavoided,whichispossiblebyusing smaller values of 𝑟, with the side effect of increasing the number of deflations that are necessary to move away from a local minimum region to the next. However, when it is desired that different design vari- ablesaredeflateddifferently,thehypercuboidapproachpro- videsflexibilityinselectingtherequired 𝑟𝑖valuesalongeach dimension. 3. Results and discussions Theproposeddeflationconstraintsareappliedacrossvarious problems to evaluate the capabilities of the method. First, a test case study based on a double-cosine function is pre- sented. The next two case studies represent different sizing problems encountered in the design of composites. 3.1. Test case: double-cosine function Thedouble-cosinefunctionofEquation29isusedtoverify whether the proposed deflation constraint can be used to obtain all the 13 minima points within a bounded region of the objective function. 𝑓(𝑥1,𝑥2)=−cos(𝑥1𝜋)cos(𝑥2𝜋) (29) The results below were generated for the three defla- tion constraint approaches herein proposed: hypersphere, SS Bangera, SGP Castro Preprint Page 7 of 16 Deflation constraints for global optimization of composite structures (a) (b) Figure 7: (a) Plot of Gaussian deflation function distribution using hypecube approach (Equation 25) (b)Plot of Gaussian deflation function distribution using hypercuboid approach (Equation 28) x*r Figure 8: Masking of local minima (green) due to a large value of𝑟, and creation of false minima (red) at the boundaries of the deflated region. hypercube and hypercuboid. Here, the same starting point for the optimizer is used for all cases, even after deflation is performed. However, as shown in Appendix B, it was observed with the double-cosine function that consistently fewer deflation iterations are necessary if the starting point afterdeflationisbesettothelatestfoundsolution,suchthat the search for a new solution proceeds from there. 3.1.1. Hypersphere Thehyperspheredeflationconstraintisappliedtothedouble- cosinefunctionwithintheboundedregion −2.5𝜋≤𝑥1,𝑥2≤ +2.5𝜋and using as starting point 𝑥1,𝑥2= −5.5,−5.5. The radius of the deflation hypersphere is 𝑟=𝜋. Figure 9 showsaplotofthedouble-cosinefunction,wherethedashed regions represent negative values; the red dot depicts the startingpointoftheoptimizer;whereastheblackdotsarethe obtainedsolutions.Afterdeflatingthefunction29times,all the13minimawithintheboundspacehavebeendiscovered. Note the false minima points are created at the end of the regions intersecting with the deflation radius, represented as red dots in Figure 8. The presence of false minima is not a problem for the overall optimization result given that 6 4 2 0 2 4 6 x16 4 2 0246x2 x0(-5.50,-5.50)Deflated 29 times using hypersphere approachFigure 9: Hypersphere approach applied to locate minima of the double-cosine function (Equation 29) enough additional deflations are carried on, because when new deflation constraints are added the false minima points areeventuallyovercome.Thevalueofdeflationradius 𝑟=𝜋 has shown to be small enough, preventing local minima of being masked, e.g. the green dots in Figure 8. Hence, to guaranteetherobustnessofthemethod,itisrecommendedto usesmallvaluesforthedeflationradius 𝑟,despiteapossible increase in the number of deflation iterations needed to discover the required solutions. SS Bangera, SGP Castro Preprint Page 8 of 16 Deflation constraints for global optimization of composite structures 6 4 2 0 2 4 6 x16 4 2 0246x2 x0(-5.50,-5.50)Deflated 51 times using hypercube approach Figure 10: Hypercube approach applied to locate minima of the double-cosine function (Equation 29) 3.1.2. Hypercube Forthehypercube,theboundsandstartingpointarethesame as for the hypersphere, with the side of the hypercube set to𝑟=𝑝𝑖. As seen in Figure 10, 51 deflation iterations are required to find all 13 minima points, as compared to 29 of the hypersphere. However, this difference is not verified when the starting point for the new optimization, after each deflation step, is set to be the latest minimum point (see Appendix B), with both the hypersphere and the hypercube approachesrequiring24deflationiterationstofindallthe13 minima points. 3.1.3. Hypercuboid 6 4 2 0 2 4 6 x14 2 024x2 x0(-5.85,-3.24)Deflated 64 times using hypercuboid approach Figure 11: Hypercuboid approach applied to locate minima of the double-cosine function (Equation 29)For the hypercuboid, the double-cosine function is skewed in one dimension, as follows: 𝑓(𝑥1,𝑥2)=−cos(𝑥1𝜋)cos(2𝑥2𝜋) (30) whereitcanbeseenthateachsub-spaceinthewholedesign space is now shaped like a rectangle, making a proper test for the hypercuboid approach. The design space is bounded such that−2.5𝜋≤𝑥1≤+2.5𝜋, and−1.5𝜋≤𝑥2≤+1.5𝜋. The dimensions used for the 2-sided deflation hypercuboid (i.e. a rectangle in two-dimensions) are, 𝑟1=𝜋, along 𝑥1; and𝑟2=𝜋∕2, along𝑥2. As seen in Figure 11, the functionhastobedeflated64timestoyieldallthe13minima within the bounded region. When the starting point after eachdeflationissettobethelatestminimumpoint,asshown in Appendix B, the hypercuboid approach required only 44 deflations to find all the 13 minima. 3.2. Application to composite design Case study 1 presents a genetic algorithm-based optimiza- tionseekingminimummassandconstrainedbybucklingand strength,wherethelaminateisparameterizedusingdiscrete ply orientations. Case study 2 presents a gradient-based optimization, where the laminate is parameterized using laminationparameters(LP)andthetotalthickness,typically applied in the design of variable-stiffness (VS) laminates. 3.2.1. Case study 1: deflation in gradient-free discrete-based optimization 4590Y X 0 ab NxNy Nx Ny Figure 12: The geometry of laminated plate and applied loads (12). Riche et al. (12) presents optimum benchmark cases for composite plates under bi-axial compression. As illustrated in Figure 12, 𝑁𝑥and𝑁𝑦are the loads acting along the 𝑥 and𝑦axes,and𝜆representsaloadamplitudeparameterthat affects all loads. The rectangular plates with length 𝑎and width𝑏are simply supported. Discrete fiber angles of each layer are optimized, aiming minimum weight for a given design load, and the laminate is considered to have variable numberofplieswhileremainingsymmetricandbalanced.It is assumed that only angles 0◦,90◦, and±45◦can be used, and that all plies have the same thickness ℎ𝑝𝑙𝑦. SS Bangera, SGP Castro Preprint Page 9 of 16 Deflation constraints for global optimization of composite structures Toaccountforavariablenumberofpliesinthisdiscrete optimization, the thickness value associated with each ply ℎ𝑝𝑙𝑦, becomes an additional design variable. Schläpfer (36) introducedtheconceptofghostlayers,whichessentiallyisa layer that carries information about material properties and alsotheplyfiberorientationbutcanbeassociatedwithazero thicknessvalue.Here,everylayerconstitutingthelaminates are allowed to become ghost layers. Two constraints are used in this optimization problem, represented by their respective load factors, being: critical buckling (𝜆𝑐𝑏); and material failure ( 𝜆𝑐𝑠). These are calcu- lated as detailed in Appendix C. In the optimization, these constraintsarenormalizedusingthedesignloadfactor( 𝜆𝑡), such that the optimization problem becomes: minimize 𝜃,𝑡∈𝑛(𝑗∑ 𝑖=1ℎ𝑖)𝑎𝑏 subject to−𝑀𝑆𝑐𝑏≤0 −𝑀𝑆𝑐𝑠≤0(31) where(𝑗∑ 𝑖=1ℎ𝑖)𝑎𝑏isthetotalvolumeofthelaminate;defining the margins of safety for buckling 𝑀𝑆𝑐𝑏and failure𝑀𝑆𝑐𝑠 as: 𝑀𝑆𝑐𝑏=(𝜆𝑐𝑏∕𝜆𝑡)−1 𝑀𝑆𝑐𝑠=𝜆𝑁𝑦=(𝜆𝑐𝑠∕𝜆𝑡)−1(32) Ageneticalgorithm(GA)optimizerisselectedgiventhe discretenatureoftheoptimizationofEquation31.Adopting the ghost layer approach, the thickness of each ply can either be zero or ℎ𝑝𝑙𝑦, and each ply angle can assume one of the three values; 0◦,90◦or±45◦. For implementation, the GA present within the Python module pymoo (37) is utilized. The deflation constraint herein proposed is applied to the optimizer to support the claim that the developed methodology is applicable in any optimization scheme that supports inequality constraints, even heuristic optimization schemes based on discrete variables. For this optimizer, the conversion from continuous to discrete variables can be done within the optimization, directly affecting the value of the objective function. In gradient-based methods this con- version is done after the optimization process is complete, in a post-processing step, usually affecting the objective and constraints. Regarding the main GA parameters, the population size herein utilized is 20, and the number of generations is 30, being the termination criterion. Table 1 shows 4 stacking sequences with 48 plies with their corresponding load factors, extracted from Riche et al. (12). Their optimization focused on maximizing the bucklingandfailureloadfactorsforafixedamountofplies. Note the lowest load factor of 13,518.66corresponding to the failure constraint. Here, with the intent to allow the optimizer to remove some layers, the design load factor is set to𝜆𝑡= 10,000. A total of 24 independent and discreteLoad factor Stacking sequence Buckling Failure (902,±454,04,±45,04,±45,02)𝑠 14,168.12 13,518.66 (±453,02,±452,02,902,04±45,02)𝑠14,134.76 13,518.66 (902,±453,02,±45,02,±45,04,±45,02)𝑠14,013.71 13,518.66 (±452,02,±452,902,04,±45,02,±45,02)𝑠13,662.61 13,518.66 Table 1 Results for the buckling load maximization problem obtained by Riche et al. (12). (48 plies, 𝑎=20in•𝑏=5in,𝑁𝑥=1lb, 𝑁𝑦=0.125lb). Margin of Safety Stacking sequence PliesBuckling Failure (904,06,902,±456)𝑠 487.44% 15.85% (±453,04,±452,02,±452,02)𝑠 446.49% 12.94% (±454,04,±45,02,±45,02,±45)𝑠4410.3% 12.94% (902,06,904,02,904,±453)𝑠 483.84% 9.21% (904,04±452,902,02,904,02,902)𝑠4813.73 % 5.155 % (902,±45,02,±45,02,±454,04)𝑠442.55% 14.2% Table 2 Results for the buckling load maximization problem obtained using the GA optimizer with deflation. (6 solutions generated using hypersphere approach with 𝑝=3,𝜎=1and𝑟=1) design variables are utilized, being 12 for the ply angles and 12 for the thickness of each ply. Note that due to the assumptionofsymmetricandbalancedlaminated,thissetup will allow a maximum of 48 plies in the laminate. The solutions obtained are displayed in Table 2. The first solution is without deflation, resulting in a laminate with 48 plies, not changing from the baseline value. After deflation is applied, new constraints are added to the GA to enablefindingnovelsolutions,wherefivenewsolutionsare herein reported. The third solution is the most promising, combining the smallest amount of plies herein encountered with relatively high margins of safety. Hence,applyingdeflationtothisGAoptimizationproves that the deflation constraint works even in heuristic algo- rithmsforoptimizationanditcreatesmeanstoforcefurther exploration of the design space. 3.2.2. Case study 2: deflation enabling gradient-based retrieval of fiber angles Setoodeh et al. (39) discuss a method for generating curvi- linear fiber paths for the manufacturing of VS laminates using AFP, where part of the work consists of finding the fiberangleswhenthelaminationparameters(LP)arealready known. Here, the deflation constraints are applied to enable the use of gradient-based methods to solve the multimodal optimization problem of retrieving the fiber angles. A balanced and symmetric layup is assumed, and the objective function is given by: SS Bangera, SGP Castro Preprint Page 10 of 16 Deflation constraints for global optimization of composite structures Figure 13: Cantilever plate with uniform load (38). minimize 𝜃|𝑉−𝑉∗| subject to0◦⩽𝜃⩽90◦(33) which minimizes the least square distance between the known LPs 𝑉∗and those calculated 𝑉based on the fiber anglesbeingoptimized.ThecalculatedLPsusetherelations inEquationD.7.TheoptimizationofEquation33isalready non-convex even for a single design variable, as shown in Figure 14. 0 20 40 60 80 0.5 1.0 1.5 2.0 2.5|V V*| Non-convex nature of design space Figure 14: Variation of the objective function plotted against 𝜃for a node along the laminate The optimal set of lamination parameters 𝑉∗is found through a compliance minimization problem, where the feasibility constraints of LPs allow for a convex design space according to Fukunaga and Sekune (25; 26). The minimization problem is formulated as:Normalized compliance Method used ̄𝐶∗ ̄𝐶 %Difference Baseline from reference (39) 0.0374 0.0389 4.00 Gradient-based optimizer 0.0389 0.0414 6.42 Genetic Algorithm 0.00392 0.0419 6.88 Table 3 Normalized compliance obtained with optimum LPs and re- trieved stacking sequence. minimize 𝑉𝑖1 2𝑁𝑇 𝑖⋅𝐴−1(𝑉𝑖)⋅𝑁𝑖 subject to2𝑉2 1(1−𝑉3)+2𝑉2 2(1+𝑉3)+ 𝑉2 3+𝑉2 4−4𝑉1𝑉2𝑉4⩽1 𝑉2 1+𝑉2 2⩽1 −1⩽𝑉3⩽1(34) where𝑁𝑖isthevectoroftheresultantforcesforthe 𝑖𝑡ℎnode; 𝐴is the in-plane stiffness matrix which is a function of 𝑉𝑖 calculated with Equation D.8; and 𝑉𝑖= {𝑉𝑖1,𝑉𝑖2,𝑉𝑖3,𝑉𝑖4}, the vector of in-plane LPs of the 𝑖𝑡ℎnode. For balanced VS laminates, LPs 𝑉2=𝑉4= 0. The material properties used forthepresentnumericalevaluationsare: 𝐸11=181.0𝐺𝑃𝑎, 𝐸22=10.3𝐺𝑃𝑎,𝐺12=7.17𝐺𝑃𝑎,and𝜈12=0.28. The results for the retrieved fiber angles can be verified byplottingthedistributionofLPsonceagainandcomparing it with the optimal LPs. In this case a non-dimensional compliancetermisalsoutilisedtonumericallycomparethe results of the retrieved angles and optimal LPs. The non- dimensional compliance term is given as (38): ̄𝐶=𝐸22ℎ𝑏3𝐶 𝑞2 0𝑎5(35) where𝐸2is the transverse modulus of elasticity; ℎ,𝑎and 𝑏are respectively the thickness, length and width of the laminate;𝑞0is the uniformly distributed load; and 𝐶is the compliance obtained after optimization, given as: 𝐶=1 2𝐹𝑇⋅𝑈 (36) where𝐹and𝑈arerespectivelythevectorsofexternalforces and displacements. Forthesakeofbenchmarking,astandardGAalgorithm (37) is also used to retrieve the fiber angles. Table 3 shows normalized compliance values obtained with the optimum LPs (̄𝐶∗); and with the retrieved fiber angles ( ̄𝐶). For the baselineresult,Setoodehetal.(39)achievedfiberanglesthat led to just 4% difference in compliance, and this is because the authors utilized a maximum curvature constraint for the curved fibers of (𝜅= 3.333𝑚−1)that further constrained theminimumcomplianceoptimizationofEquation34.Fig- ure15illustratesthefiberangledistributionscorresponding to Table 3. SS Bangera, SGP Castro Preprint Page 11 of 16 Deflation constraints for global optimization of composite structures (a)+𝜃by Setoodeh et al. (39) (b)−𝜃by Setoodeh et al.(39) (c)+𝜃using gradient-based optimization with deflation (d)−𝜃using gradient-based optimization with deflation (e)+𝜃using GA (f)−𝜃using GA Figure 15: Distribution of fiber angles of balanced laminate. ( 𝑎∕𝑏=3,and31×11nodes) The results for the gradient-based method using the deflation constraint are obtained after deflating the design space 7 times, such that a total of 8 solutions are ob- tained after finding the retrieved angles that best minimizes the objective function of Equation 33. In Table 3, it can be seen that the percentage difference between the non- dimensional compliances obtained with 𝑉∗and𝑉is ap- proximately6.68%,achievingaslightlybetterresultthanthe benchmark GA that has a difference of 6.88%. Figure16showsthedistributionsoflaminationsparam- eters obtained with the gradient-based method constrained by deflation, whereas Figure 17 shows these distributions obtained with the GA. Again, the good agreement shows that the deflation constraints enabled global search by the gradient-based optimizer. Figure 18 shows the convergence of the normalized compliance obtained obtained with each optimization, after progressive deflation steps, indicating convergence towards a global optimum. 4. Conclusion The present study proposed a novel set of deflation con- straints applicable to gradient-based and heuristics-based optimizationalgorithms.Virtually,theproposedconstraints canbeusedinanyoptimizationalgorithmthatiscompatible with constraints. Deflation constraints force the optimizer to look for other optimal points, even if the optimization is restarted from the same initial guess. Such behavior makes it even possible for gradient-based optimizers to exploremultimodal design spaces and ultimately find a global min- ima, after an enough number of deflation constraints are progressively added. The three deflation constraint schemes herein proposed: hypersphere, hypercube, hypercuboid; were demonstrated in detail for a double-cosine function. Thereafter, two case studies related to composite design and optimization were investigated, where the first case study demonstrated the minimizationoftheweightofaplateforagiventargetdesign load,inwhichdeflationisusedtocreatenewconstraintsfor a GA algorithm. The second case study covered the opti- mizationofLPsfollowedbytheretrivaloffiberangles,with deflation applied only to the gradient based optimizer. The obtained results were promising and showed the versatility ofthedevelopedmethodology,uncoveringnewpathsforits application composite design. Future studies could focus in design cases involving multi-objective optimizations, investigating how the use of deflation constraints could improve the Pareto fronts ob- tained with gradient-based optimizers. A. Deflation constraint - KKT proof This appendix provides a proof on how the modified prob- lem, after the introduction of the deflation constraint shown in Equation 10, provides a solution ( 𝑥∗,𝑦∗) that also satis- fies the stationary conditions of the original problem. The modified problem is given as (33): SS Bangera, SGP Castro Preprint Page 12 of 16 Deflation constraints for global optimization of composite structures Lamination Parameter = V1 0.857 0.571 0.286 0.000 0.286 0.571 0.857 (a)𝑉1∗ Lamination Parameter = V3 0.857 0.571 0.286 0.000 0.286 0.571 0.857(b)𝑉3∗ V1from retrieved balanced stack (deflation) 0.857 0.571 0.286 0.000 0.286 0.571 0.857 (c)∼ 𝑉1 V3from retrieved balanced stack (deflation) 0.857 0.571 0.286 0.000 0.286 0.571 0.857 (d)∼ 𝑉3 Figure 16: Distribution of Lamination parameters obtained by applying gradient-based optimization method with upper images representing optimum LP values and lower distributions obtained with retrieved angles shown in Figure 15 (c) and (d). ( 𝑎∕𝑏=3, and31×11nodes) minimize 𝑥∈𝑛𝑓(𝑥) subject to𝑐(𝑥)=0 𝑚(𝑥;𝑥1)≤𝑦 𝑙≤𝑥≤𝑢(A.1) It is to be proved that, 𝑥∗is a regular KKT point to the problem given in Equation 7 and that 𝑥∗≠𝑥1. Let a new Lagrange multiplier 𝜂be introduced in addition to the ones in Equation 8 associated with the deflation constraint. Therefore, the stationary conditions now will be: ∇𝑥𝑓(𝑥)+∇𝑥𝑐(𝑥)𝑇𝜆+𝑧++𝑧−+𝜂∇𝑥𝑚(𝑥;𝑥1)=0 𝜂=0 (A.2) Itcanbesaidthat 𝜂wouldbe0atanyKKTpoint,hence the stationary conditions for Equation 7 would be satisfied. SincetheconstraintsgiveninEquation7areasubsettothose in Equation A.1, 𝑥∗must be a feasible point in the original problem (Equation 7) as well, and the complementary con- ditions of those constraints must be satisfied. Since the obtained 𝑦∗is finite, (𝑥∗,𝑦∗) is feasible to the deflation constraint problem in Equation A.1 and 𝑚is bounded from below, then 𝑚(𝑥∗;𝑥1)must be finite, which implies that(𝑥∗≠𝑥1. This proof can also be generalized to inequality-constrained NLPs, hence the proposed deflation constraint approach is a generic and non-invasive way ofutilizing the deflation method for an optimization problem (33). B. Using last optimum as starting point Afteranewdeflationconstraintisadded,theoptimizationis restartedusingasstartpointthelastfoundminimum,which resulted in a decrease of the number of deflations required to locate all minima points of the double-cosine function. Figures 19, 20 and 21 show the evaluated points. C. Buckling and failure constraints The critical buckling load factor can be calculated with the following analytical expression (12): 𝜆𝑏(𝑚,𝑛) 𝜋2=𝐷11(𝑚 𝑎)4+2(𝐷12+2𝐷66)(𝑚𝑛 𝑎𝑏)2+𝐷22(𝑛 𝑏)4 (𝑚 𝑎)2𝑁𝑥+(𝑛 𝑏)2𝑁𝑦 (C.3) where𝑚,𝑛are respectively the half waves along the length and width of the plate; and 𝐷𝑖𝑗are the components of the bendingstiffnessmatrix 𝐷𝐷𝐷ofthelaminate,fromtheclassical laminatedplatetheory(1).Thevaluesof 𝑚and𝑛needtobe foundsuchthatitminimizes 𝜆𝑏toyieldthecriticalbuckling load𝜆𝑐𝑏;andtheyvarywiththenumberofplies,constitutive properties, plate dimensions and the load case. The material failure constraints are given by the prin- cipal allowable strains for each ply. The calculation of the SS Bangera, SGP Castro Preprint Page 13 of 16 Deflation constraints for global optimization of composite structures Lamination Parameter = V1 0.857 0.571 0.286 0.000 0.286 0.571 0.857 (a)𝑉1∗ Lamination Parameter = V3 0.857 0.571 0.286 0.000 0.286 0.571 0.857(b)𝑉3∗ V1from retrieved balanced stack (GA) 0.857 0.571 0.286 0.000 0.286 0.571 0.857 (c)∼ 𝑉1 V3from retrieved balanced stack (GA) 0.857 0.571 0.286 0.000 0.286 0.571 0.857 (d)∼ 𝑉3 Figure 17: Distribution of Lamination parameters obtained by applying GA optimization method with upper images representing optimum LP values and lower distributions obtained with retrieved angles shown in Figure 15 (e) and (f). ( 𝑎∕𝑏=3,and31×11 nodes) 0 2 4 6 8 10 Number of solutions0.04160.04170.04180.04190.0420Normalized compliance Convergence of compliance Figure 18: Convergence of Normalized Compliance with in- crease in number of solutions for deflation principal strains are done using the following relations: 𝜆𝑁𝑥=𝐴11𝜖𝑥+𝐴12𝜖𝑦 𝜆𝑁𝑦=𝐴12𝜖𝑥+𝐴22𝜖𝑦(C.4) 𝜖𝑖 1=cos2𝜃𝑖𝜖𝑥+sin2𝜃𝑖𝜖𝑦 𝜖𝑖 2=sin2𝜃𝑖𝜖𝑥+cos2𝜃𝑖𝜖𝑦 𝛾𝑖 12=sin2𝜃𝑖(𝜖𝑦−𝜖𝑥)(C.5) where𝜃𝑖referstothefiberangleassociatedwiththe 𝑖𝑡ℎlayer, with𝜖𝑥and𝜖𝑦beingtheglobalstrainsoftheplate.Theterm 6 4 2 0 2 4 6 x16 4 2 0246x2 x0(-5.78,5.78)Deflated 24 times using hypersphere approachFigure 19: Hypersphere approach applied to locate minima of the double-cosine function (Equation 29) with updated start points 𝜆in Equation C.4 is the load amplitude and it becomes 𝜆𝑐𝑠 for the lowest value of 𝜆such that only one of the principal strainsinoneofthelayersexceedstheallowablestrain.Inthe casehereinconsidered, 𝛾𝑥𝑦isnotpresentbecausethelayers SS Bangera, SGP Castro Preprint Page 14 of 16 Deflation constraints for global optimization of composite structures 6 4 2 0 2 4 6 x16 4 2 0246x2 x0(3.64,-2.64)Deflated 24 times using hypercube approach Figure 20: Hypercube approach applied to locate minima of the double-cosine function (Equation 29) with updated start points 6 4 2 0 2 4 6 x14 2 024x2 x0(-2.64,1.45)Deflated 44 times using hypercuboid approach Figure 21: Hypercuboid approach applied to locate minima of the double-cosine function (Equation 29) with updated start points of the laminate are always balanced and symmetric. The allowable strains for the adopted material, graphite epoxy composite, are given as: 𝜖𝑖 1≤0.008 𝜖𝑖 2≤0.029 𝛾𝑖 12≤0.015(C.6) D. Lamination parameters In this case study a few additional aspects need to be covered for enabling the reader’s understanding. It is im- portant to know the relation for lamination parameters tounderstand how it relates to the design variables and also the objective function. Lamination parameters were first introduced by Tsai et al. (23; 40) to represent the laminate layup configuration in a compact form. The transformation properties are derived using trigonometric relations (sines and cosines) in terms of multiple angles. The Lamination Parameters are non-dimensional and can in turn be used to obtain the in-plane and bending properties (ABD matrices) as they are related. A great advantage of using lamination parametersisthatthenumberofdesignvariablesisreduced and is independent of the number of layers (20; 39). The 12 variables are given as: (𝑉1𝐴,𝑉2𝐴,𝑉3𝐴,𝑉4𝐴)=∫1 2 1 2(cos2𝜃,sin2𝜃,cos4𝜃,sin4𝜃)𝑑𝑧, (𝑉1𝐵,𝑉2𝐵,𝑉3𝐵,𝑉4𝐵)= 4∫1 2 1 2̄ 𝑧(cos2𝜃,sin2𝜃,cos4𝜃,sin4𝜃)𝑑̄ 𝑧, (𝑉1𝐷,𝑉2𝐷,𝑉3𝐷,𝑉4𝐷)= 12 ∫1 2 −1 2̄ 𝑧2(cos2𝜃,sin2𝜃,cos4𝜃,sin4𝜃)𝑑𝑧 (D.7) where𝑉𝑖𝐴,𝑉𝑖𝐵and𝑉𝑖𝐷are in-plane, coupling and flexural laminationparameters. ̄ 𝑧isthenormalizedzcoordinate (̄ 𝑧= 𝑧∕ℎthroughthethicknessand 𝜃isthefiberorientationangle of that layer. The [A], [B], and [D] Matrices can be written as a linear function of the lamination parameters and the parameters of Tsai and Pagano as (23), 𝐴=ℎ(Γ0+Γ1𝑉1𝐴+Γ2𝑉2𝐴+Γ3𝑉3𝐴+Γ4𝑉4𝐴) 𝐵=ℎ2∕4(Γ1𝑉1𝐵+Γ2𝑉2𝐵+Γ3𝑉3𝐵+Γ4𝑉4𝐵) 𝐷=ℎ3∕12(Γ0+Γ1𝑉1𝐷+Γ2𝑉2𝐷+Γ3𝑉3𝐷+Γ4𝑉4𝐷) (D.8) where the material invariant matrices Γ𝑖’s are given by the parametersofTsaiandPaganowhicharematerialinvariants. CRediT authorship contribution statement Sankalp S. Bangera: Methodology, investigation, val- idation, formal analysis, software, writing - original draft. SaulloG.P.Castro: Conceptualization,supervision,method- ology, software, writing - original draft, data curation. Data Availability Statement The data that support the findings of this study are availablefromthecorrespondingauthor,S.G.P.Castro,upon reasonable request.