# A unified mixed-variables framework for manufacturing-constrained optimization of composite laminates applied to conventional and double-double layups

Alejandro Garcia-Pis, Christian Gogu, Laurent Michel, Stéphane Grihon, Simone Coniglio

► To cite this version:

Alejandro Garcia-Pis, Christian Gogu, Laurent Michel, Stéphane Grihon, Simone Coniglio. A unified mixed-variables framework for manufacturing-constrained optimization of composite laminates applied to conventional and double-double layups. Structural and Multidisciplinary Optimization, 2025, 68 (9), pp.176. 10.1007/s00158-025-04101-z. hal-05466289

HAL Id: hal-05466289

https://ut3-toulouseinp.hal.science/hal-05466289v1

Submitted on 19 Jan 2026

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# A unified mixed-variables framework for manufacturing-constrained optimization of composite laminates applied to conventional and double-double layups

Alejandro Garcia-Pis<fnref n="1" />,<fnref n="2" /> · Christian Gogu<fnref n="1" /> · Laurent Michel<fnref n="1" /> · Stéphane Grihon<fnref n="2" />,<fnref n="3" /> · Simone Coniglio<fnref n="2" />,<fnref n="3" />

Received: 19 February 2025 / Revised: 19 June 2025 / Accepted: 27 July 2025 / Published online: 10 September 2025

© The Author(s) 2025

## Abstract

Structural optimization may imply in many instances the presence of categorical, discrete, and continuous variables. This is particularly the case for composite structures due to the stacking sequence definition which involves not only structural considerations, but also a set of complex manufacturing- and cost-related constraints. While approaches ranging from relaxation of all variables to evolutionary methods have been proposed in the literature, it is still challenging to deal consistently with all manufacturing and structural constraints, especially in high-dimensional problems (hundred or more design variables). This work proposes a mixed-variables optimization formulation that can integrate a large variety of manufacturing and structural constraints in a monolithic approach based on a modified Bi-level Outer Approximation algorithm. The proposed framework allowed the treatment of composite optimization problems accounting for catalog selection of stacking sequences, manufacturing compatibility of adjacent elements, as well as variability control schemes. Various optimization formulations have been proposed tailored both for conventional laminates and for double-double laminates. The proposed approaches were tested on both an academic use-case, a multi-cell panel and, on an industrial use-case, a vertical tail plane front spar.

**Keywords** Structural optimization · Mixed-variable optimization · Stacking sequence optimization

## Acronyms

CaF Categorical formulation

CC Catalog compliant

BLOA Bi-level outer approximation

DD Double-double laminates

FEA Finite-element analysis

GA Genetic algorithm

MVF Mixed-variable formulation

MMA Method of moving asymptotes

NGSA Next-generation structural analysis

LPF Lamination parameter formulation

QUAD Conventional laminates

VTP Vertical tail plane

Responsible editor: Marco Montemurro.

<figure><img src="image_2.png" /></figure>

Christian Gogu
christian.gogu@gmail.com

## 1 Introduction

Composite structures are becoming more and more predominant in aerospace applications. The lightweight properties of composite materials make them more appealing than conventional metals for many applications. One of their advantages lies in the possibility to tailor their mechanical properties to the specific loading, notably through the ply layup. This opens up vast opportunities for optimized structures, specific to each location on the aircraft. However, the production of this type of structure is linked to a series of specific manufacturing constraints that have to be addressed, such that the optimized structures can indeed be materialized in practice (Bruyneel 2008). These constraints are traditionally a set of rules intended to maximize the manufacturing quality and structural performance of these materials, of which an extensive review is provided in Irisarri et al. (2014), and which can be divided into two main groups:

* *Laminate level constraints.* These are the constraints that affect the individual laminates of a structure. This group includes aspects, such as: layup balance, contiguity of identical plies, surface ply restrictions (e.g., having 45°

<footnote n="1">¹ Institut Clement Ader (ICA), CNRS, UPS, ISAE-SUPAERO, INSA, Mines Albi, Université de Toulouse, Toulouse, France</footnote>

<footnote n="2">² Airbus Operations SAS, Toulouse, France</footnote>

<footnote n="3">³ Institute of Technology IRT Saint Exupéry, Toulouse, France</footnote>
plies as surface plies), rules regarding ply orientation variations, etc.

* *Compatibility constraints.* These are those affecting the structure at a global level. This group includes blending and drop-off constraints. Blending states that when two laminates are placed next to each other in the structure, all the plies in the thinnest must continue in the thickest. The drop-off constraint sets a limit on how many plies in thickness two contiguous laminates can differ.

The complexity of modeling the effect of plies orientation and positioning on the mechanical properties limits the optimization strategies to three main groups: the parametrization-based, the categorical, and, more recently, the penalization-based formulations (Venkataraman and Haftka 1999). The parametrization strategies relax the problem to the point where only continuous variables are considered. The lamination parameters are a common choice [e.g., Miki (1985)] as they provide a continuous design space. Even if locally and under certain projections, this design space can be convex, it is not in general, as proven in Picchi Scardaoni and Montemurro (2022). However, the solution of these relaxation strategies cannot always be transformed into a unique manufacturable stacking sequence matching exactly the continuous values. The manufacturing rules and the discrete nature of the stacking sequence choice in most applications lead to non-unique or approximate matches to the relaxed parameters. This led to the development of the bi-stage approaches, where in a first step, the optimization with respect to the continuous parameters is performed, and then, a gradient-free method is used to find a stacking sequence with equivalent parameters close to those. This led to the development of the bi-stage approaches, where in a first step, the optimization with respect to the continuous parameters is performed, and then, a gradient-free method is used to find a stacking sequence with equivalent parameters close to those. In Jing et al. (2019), this is done through a shuffling of the plies, while the compatibility is ensured through a shared ply identification scheme. In Picchi Scardaoni and Montemurro (2022), the search propagation direction is used to generate a series of compatible stackings from the continuous results. Additionally, some frameworks provide a bi-level approach, as that in Montemurro et al. (2023), where the two optimization loops are nested and solved through metaheuristic methods. Besides lamination parameters, another parameterization commonly used in optimization is the polar parametrization. They have a direct physical meaning and are related to tensor invariants. Several examples of its application can be found in Montemurro and Catapano (2017); Catapano and Montemurro (2014); Picchi Scardaoni et al. (2021). The categorical formulations and optimization strategies are those that directly use the stacking sequence as variables. These formulations are commonly

addressed with heuristic solvers. Several examples can be found in Nikbakt et al. (2018). One of these methodologies is the Stacking Sequence Table approach from Irisarri et al. (2014). This approach uses a genetic algorithm to find the series of stacking sequences that optimizes the mass of the structure, accounting for the manufacturing and structural constraints. Finally, there are also multiple works on the use of relaxation strategies to ease the optimization in the categorical variables domain. The Discrete Material Optimization (DMO) (Lund and Stegmann 2005) approach is similar to the Solid Isotropic Material with Penalization (SIMP) method of Bendsoe (1989) in that for every ply, a weighted combination of orientations is set and constraints or penalization are imposed in the problem to make this combination converge to a categorical choice. This method was extended in Sorensen and Lund (2013); Sorensen et al. (2014) to include a relaxed variable for the presence of individual plies, allowing the number of plies in the structure to vary and to consider manufacturing constraints, developing the Discrete Material and Thickness Optimization (DMTO).

Double-double composites (DDs) (Tsai 2021) are a layup methodology that relies on the repetition of sub-elements called super-plies to generate the stacking sequence rather than a set of manufacturing or stacking rules as in the case of usual composite laminates, referred to here as QUAD. A DD layup is composed of the repetition of a super-ply, which is a set of 4 plies with orientations [Φ / − Ψ / − Φ / Ψ], where Φ and Ψ are two given angles. A DD laminate can present several advantages over QUAD laminates. In terms of manufacturing, one advantage is the ease of generating a thickness variation while maintaining the same super-ply in the structure through a “card-sliding” approach. Also, it is claimed in Tsai et al. (2017) that different super-plies can be placed together through a staggered transition region, effectively removing the blending constraint common in the conventional manufacturing. The DD approach thus has the potential to generate thinner, and thus, lighter laminates than QUAD for the same loading, thanks to its simpler manufacturing rules and its adaptability (Vermes et al. 2021). Among the drawbacks of the DD approach, we can mention that it is unable to cover the entire lamination parameters space (although QUAD also has significant limitations). Note that DD also faces manufacturability challenges, though they are somewhat different from those faced for QUAD. A main challenge for DD laminates is due to the fact that they are still a relatively recent design concept and thus do not have the same industrial maturity as QUAD laminates.

Given the potential of DD laminates, there has been a rising interest in recent years regarding the optimization of composite structures based on this manufacturing technology. The work presented in Wang et al. (2023) showed several applications of the SIMP method for DD composites optimization. To address the identification of

<figure><img src="image_3.png" /></figure>
the ideal DD for a given loading, an optimization strategy based on minimizing the mass and material-effort threshold was developed in Kappel (2022). The use of several machine learning algorithms to predict this optimal DD was explored in Zhang et al. (2022). In Wang et al. (2024), deep-learning techniques are used to propose optimal DD grid/skin considering only buckling constraints. In Garofano et al. (2023), the fuselage section frames were optimized through an iterative evaluation process, where the reserve factors for the Tsai-Wu first ply failure criterion were used to define the optimal thickness and super-ply for each cell of the section.

In industry applications, many times the design is also constrained by the design chain. Suppliers and manufacturers usually agree on the use of a catalog of materials, limiting the number of laminates that can be chosen for the designs. This is done to provide a cost-effective, high-quality set of materials, which can go through an extensive testing and validation campaign, typically required for certification of aeronautical structures.

Based on the previous literature review, it appears that while optimization of DD laminates has been previously investigated, there is still a gap in terms of optimization frameworks that can simultaneously account for a relatively large number of variables and a large variety of constraints, including local failure, global stiffness and local buckling, as well as industrial constraints that would involve categorical choices from a catalog. In terms of optimization of QUAD composites, there also appears to be a gap in terms of frameworks that can address such problems while simultaneously accounting for manufacturing constraints, notably compatibility constraints.

Accordingly, the objectives of this work are to propose and investigate a new algorithmic framework able to deal with DD and QUAD laminates' optimization involving the industrial context constraints reviewed above. To this purpose, the new framework seeks to:

* select the optimal choices from a stacking sequence catalog

* account for structural and manufacturing constraints

* control the variability of choices in the solution.

We will consider optimization with both QUAD and DD laminates for several reasons. First, both face manufacturing challenges, though somewhat different in nature; the main challenge for QUAD is the consideration of compatibility constraints, while the main challenge for DD is its low industrial maturity. Second, there seems to be a literature gap in terms of comparing on the same structure the optimization of a QUAD laminate with manufacturing constraints and the optimization of a DD laminate with manufacturing constraints. We thus seek to contribute to filling these gaps.

The rest of the paper is structured as follows. In Sect. 2, the optimization problem is introduced through an academic use-case, including the formulations of the different failure criteria. In Sect. 3, we propose different formulations for solving the considered optimization problem in the different cases, as well as a brief description of the optimization formulations and solvers used for this analysis. The results and trade-off study are presented in Sect. 4. In addition, an industrial case replicating the most interesting results on the sizing process of a composite vertical tail plane spar. Finally, Sect. 5 provides concluding remarks and outlooks. In addition, the Appendix presents more in depth information about the DD properties, Appendix A, the analytical formulae used to compute the stress distribution, Appendix B, the catalog used in the QUAD use cases, Appendix C, and some additional representations of the results, Appendix D.

## 2 Optimization problem

This section presents the academic use-case, which will serve as a benchmark for the various formulations investigated to solve the layup optimization problem under industrial and manufacturing constraints. Note that an actual industrial application will be presented in the last section of this article.

The application test case considered concerns the minimization of the mass of a fully one-sided clamped multi-cell composite stiffened panel supporting the equivalent loading of a punctual force, F, applied at the center of its free edge (see Fig. 1). This use-case can be seen as representative of a fuselage panel optimization, as in Allaire and Delgado (2016). Each cell is defined through its laminate layup. For QUAD composites, this would be the layup from a catalog such as the one in Appendix C. For DD composites, this would be the super-ply angles and number of repetitions.

The panel is rectangular with a size of $L \times H$ and counts $n_p$ cells. As illustrated in Fig. 2, it is divided into $n_{px}$ on the x-axis and into $n_{py}$ on the y-axis, so the total number of cells $n_p = n_{px}n_{py}$. The total mass of the panel, m, can be computed as

$$m(t) = \rho \sum_{i=1}^{n_p} t_i \frac{H}{n_{py}} \frac{L}{n_{px}}, \quad (1)$$

where $\rho$ is the density of the material and $t_i$ is the thickness of the i-th cell, which are grouped in the vector t. At the separation between the cells in the y-axis, as well as at the top and bottom of the panel, a series of longitudinal stringers are present, but are not considered here among the sizing variables.
<figure><img src="image_4.png" /><figcaption>Fig. 1 Multi-cell clamped panel scheme ($n_{px} = 2, n_{py} = 4$)</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Fig. 2 Panel lateral view ($n_{px} = 3, n_{py} = 6$)</figcaption></figure>

We consider in this mass optimization problem 3 sets of constraints: those dealing with the first-ply failure criteria, $g_{fpf}$, and those related to the buckling criteria. The latter is itself composed of two different criteria, the pure shear buckling and the compression-shear buckling criteria, $g_{shearbuck}$ and $g_{mixbuck}$. Note that accounting for buckling constraints in composites optimization significantly complicates the optimization problems. This is why buckling is either not considered or considered in very simplified ways in most studies on composite optimization available in the literature. However, this failure criterion significantly affects the stacking sequence selection, while not the case for the first-ply failure. This is why, it is essential to consider them in industrial applications.

In this work, the blending constraint among different DD super-plies was not considered, as it is said in Tsai et al. (2017) that patching the contiguous super-plies in an interlocking manner is enough to ensure the correct stress transmission, but a general criterion for blending was not provided. In the case of QUAD composites, the manufacturing constraints of balanced, pseudo-symmetry, having at least 8% plies in each direction, skins with fixed orientations, and no more than 3 contiguous plies with the same orientation are ensured by the catalog shown in the

Appendix C. The blending and drop-off constraints are also considered, and their implementation is explained in Sect. 3.

To extract the values of the stresses, analytical expressions based on the linear tension distribution assumption are used. In Appendix B, the formulation used to compute the geometrical parameters is described. The values of the resulting stresses depend on the cell thickness and the chosen material, as they control the first and second moments of inertia of the sections. As the focus of this study was to compare different optimization strategies and not a detailed structural analysis of the multi-cell panel, these expressions were further simplified by setting the longitudinal stringers as static elements regarding the optimization problem, being set to $E_s$ and $S_s$.

In this particular scenario, with the applied force, each of the cells is sustaining a stress with normal and shear components, the first one being parallel to the x-axis and the second one in the plane transversal to the panel. This is consistent with the membrane behavior assumption used in the Classical Laminate Theory (CLT), which was used to compute laminate properties.

The normal stress can be approximated in our model following the expression:

$$ \sigma_x(x, y) = \frac{F(L - x)(y - y_c(x))}{I(x)}, \qquad (2) $$

where $y_c(x)$ is the position of the cell center on the y-axis. In a similar way, the shear stress at any given point can be computed as

$$ \tau(x, y) = \frac{Q(x, y)F}{I(x)t(x, y)}, \qquad (3) $$

where $Q(x, y)$ and $I(x)$ are the first and second moments of inertia. In both of these stresses, there are both effects of the thickness of the individual cells and their LSS, through the stiffness of each cell. Both of these variables affect the
stress redistribution as shown in expressions (31) and (32) in Appendix B (i.e., $\sigma_x = \sigma_x(E, t)$, $\tau = \tau(E, t)$).

The next subsection further defines the optimization constraints considered based on these stresses.

## 2.1 Failure criteria

We consider here a classical first-ply failure, which can be expressed for an individual cell through the 2D Tsai–Wu criterion

$$g_{fpf} = F_{xx}\hat{\sigma}_x^2 + 2F_{xy}\hat{\sigma}_x\hat{\sigma}_y + F_{yy}\hat{\sigma}_y^2 + F_s\hat{\tau}^2 + F_x\hat{\sigma}_x + F_y\hat{\sigma}_y - 1, \quad (4)$$

where

$$\left\{ \begin{aligned} F_{xx} &= \frac{1}{X_c X_t} \\ F_{yy} &= \frac{1}{Y_c Y_t} \\ F_x &= \frac{1}{X_t} - \frac{1}{X_c} \\ F_y &= \frac{1}{Y_t} - \frac{1}{Y_c} \\ F_s &= \frac{1}{S^2} \\ F_{xy} &= -\frac{1}{2\sqrt{X_c X_t Y_c Y_t}}, \end{aligned} \right. \quad (5)$$

where $X_t$ and $X_c$ are the stress limits for traction and compression in the direction of the fibers, $Y_t$ and $Y_c$ are the limits in the perpendicular direction, and $S$ is the shear stress limit. The values of these quantities used in this analysis can be found in Table 7.

The compression and shear loads may also induce buckling of the cells of the panel. These constraints were implemented here using the analytical expressions provided by Herencia et al. (2007) for a fully supported panel. Assuming that all cells are simply supported, the critical stress fluxes in both shear and compression, $N_x^{crit}$ and $N_{xy}^{crit}$, can be computed for each one of them based on their mechanical properties through the Nemeth coefficients. This assumption is based on the real case of a fuselage panel, where stiffeners are located between these cells, providing support and commonly modeled using these formulas. There is an interaction between these two types of buckling, and a destabilizing effect of the shear loading on compression buckling may appear. This is referred to as mixed-buckling. These phenomena are taken into account in the following formulation of the buckling constraints:

$$g_{shearbuck} = \frac{1}{RF_{shearbuck}} - 1 = \frac{N_{xy}}{N_{xy}^{crit}} - 1 \quad (6)$$

$$g_{mixbuck} = \frac{1}{RF_{mixbuck}} - 1 = \frac{-N_x}{N_x^{crit}} + \left( \frac{N_{xy}}{N_{xy}^{crit}} \right)^2 - 1, \quad (7)$$

where the loads for each cell can be computed as the product of the stresses and the thicknesses, $[N_x, 0, N_{xy}] = [\sigma_x t, 0, \tau t]$.

## 3 Optimization strategies and formulations

In this section, different formulations for solving the previously stated mass minimization problem are proposed. These formulations vary in terms of variables and objective function definition, in terms of optimization approach (e.g., continuous, discrete, and mixed), as well as in terms of associated optimization algorithms. The idea behind testing these multiple formulations is to investigate the advantages and disadvantages of these different strategies in the context of DD structure optimization and give a comparative insight about the accuracy of our solver. Additionally, a variability control scheme is proposed, based on the mixed-variable formulation. Furthermore, the necessary modifications to the formulation to implement the mixed-variable approach to a QUAD structure accounting for the manufacturability constraints are presented.

In the previously defined problem, for the case of DD, the definition of each cell is inherently a set of two discrete variables, a choice for the super-ply, and a number of repetitions of it. As explained in Sect. A, the normalized properties of a homogenized DD depend only on the super-ply and, therefore, its actual properties on the number of repetitions. To ensure that it is a homogenized laminate, a minimal number of repetitions must be ensured. In an industrial-like context, these designs were therefore ruled by a categorical choice of super-plies and a discrete number of repetitions for each cell. The optimal designs are then supposed to be non-continuous. However, as stated before, a fully non-continuous strategy might, in many cases, prove intractable for the optimization as the design space would be too restricted. Also, in this particular case, even a full relaxation might be acceptable thanks to the properties of DDs, as these are not restricted to a given number of orientations, a continuous design space can be defined for the super-ply, and the discrete number of repetitions can be easily translated into a continuous value. An example of this would be the *ceil* function, which retrieves the closest integer bigger than a given real value

$$t_{discrete} = \text{ceil}(t/t_{sply})t_{sply}, \quad (8)$$

where $t_{discrete}$ is the discrete solution obtained by lamination formulation and $t_{sply}$ is the thickness of a given super-ply. Then, in this case, full relaxation is simple, allowing for the comparison of different formulations ranging all the way from this fully continuous one to a fully categorical one. Note, however, that the compatibility among adjacent super-plies is not considered in this work, due to the low maturity
of the DD concept leading to a lack of consensual definition of the blending requirements for DD laminates.

In the case of QUAD composites, the problem formulation is not as straightforward. As mentioned in Sect. 1, there are some strategies dealing with this problem using fully relaxed variables (e.g., lamination parameters and laminate thickness as design variables), but these strategies are also reliant on a second optimization step where a gradient-free method is used to look for the actual stacking sequence, which is not trivial and can present scalability limitations. To surpass this limitation, while still considering manufacturing constraints, a mixed-variable bi-level approach will also be presented here.

This section is structured in the following way. To begin with, three formulations and associated optimization strategies are presented for DD composites. First, a full-continuous strategy is presented. Then, a mixed-variable strategy, Bi-Level Outer Approximation, is introduced. Finally, a full-categorical approach is outlined. The continuous lamination parameters' formulation and the mixed-variable one are then combined with a variability control scheme to ensure adherence to a given number of materials in the solution. Finally, we show that small modifications of the mixed-variable formulation allow us to address QUAD composites optimization with manufacturing constraints.

## 3.1 Lamination parameters' formulation for DD

With DD laminates, the problem presented in Sect. 2 can be relaxed to a full-continuous one. Each cell definition can be expressed in terms of continuous variables. Following the typical approach for composite structural optimization, the super-ply is defined with two continuous lamination parameters, $\xi_1$ and $\xi_2$, as expressed in Zhao et al. (2023); Tsai et al. (2017). The laminate thickness is also approximated as a continuous variable, denoted $t$. In reality, this thickness is a discrete variable, a number of repetitions multiplied times the thickness of a super-ply. This relaxation was chosen for convenience, but others, such as the polar method, could also be applied

$$ \left\{ \begin{array}{l} \min_{\xi_1, \xi_2, t} m(t) \\ \text{s.t.} \\ g_{f1}(\xi_1, \xi_2) \le 0 \\ g_{f2}(\xi_1, \xi_2) \le 0 \\ g_{f3}(\xi_1, \xi_2) \le 0 \\ g_{fpf}(\xi_1, \xi_2, t) \le 0 \\ g_{shearbuck}(\xi_1, \xi_2, t) \le 0 \\ g_{mixbuck}(\xi_1, \xi_2, t) \le 0 \\ -1 \le \xi_1 \le 1 \\ -1 \le \xi_2 \le 1 \\ t \le t \le t. \end{array} \right. \qquad (9) $$

Note that when reasoning in lamination parameters space, not all values of these parameters lead to feasible DD

laminates. The additional constraints $g_{f1}$, $g_{f2}$ and $g_{f3}$ are thus introduced as defining an approximation of the feasible space for the combinations of $\xi_1$ and $\xi_2$ that can be translated into a DD super-ply (cf. Zhao et al. 2023)

$$ g_{f1}(\xi_1, \xi_2) = \xi_2 - (1 - 2\xi_1)^2 \qquad (10) $$

$$ g_{f2}(\xi_1, \xi_2) = \xi_2 - (1 + 2\xi_1)^2 \qquad (11) $$

$$ g_{f3}(\xi_1, \xi_2) = (2\xi_1^2 - 1) - \xi_2. \qquad (12) $$

The problem, expressed in this formulation, could be solved by any continuous solver. For this analysis, the Method of Moving Asymptotes (MMA) (Svanberg 1987) was chosen. This method, widely used in the field of structural optimization, uses an inverse first-order approximation and a set of asymptotes to facilitate the resolution of the problem. The continuous solution thickness solution $t$ is post-processed using the ceiling function (8). This is a conservative approach with respect to the constraint function, which is monotonically decreasing with respect to the thickness variables. Since the $\Phi_i$ and $\Psi_i$ angles can indeed vary continuously to obtain the exact values of the optimized lamination parameters, following the identity:

$$ \left\{ \begin{array}{l} \xi_{1i} = \frac{1}{2}(\cos(2\Phi_i) + \cos(2\Psi_i)) \\ \xi_{2i} = \frac{1}{2}(\cos(4\Phi_i) + \cos(4\Psi_i)); \end{array} \right. \qquad (13) $$

for any given cell $i$, the continuous optimum $(\xi_1^*, \xi_2^*)$ will be accepted as the solution of the lamination formulation.

Note, however, that as explained earlier, this flexibility might not always be acceptable from a manufacturing standpoint, since some suppliers of DD fabrics may choose to only offer a catalog of orientations out of which the customer should choose. In this case, the optimized super-plies would need to be compared with the catalog and the closest catalog choice selected, which may or may not be conservative.

## 3.2 Mixed-variable formulation for DD

Another possible formulation for the considered optimization problem implies the use of both categorical and continuous optimization variables.

In this work, a new Bi-level Outer Approximation strategy was adapted, based on the findings of Barjhoux (2020) and the modifications of Coniglio (2023). This family of methods solves the problem by dividing it into two levels according to the variables' relaxation. The inner level, where the continuous ones are optimized for a fixed value of the non-continuous ones, is used to generate a set of hyper-planes. These hyper-planes form an approximation of the categorical problem that is solved iteratively using
Mixed Integer Linear Programming (MILP) until conver-
gence is reached. A more detailed description can be found
in Barjhoux (2020). One of these modifications introduced
in the algorithm was the addition of convexification meth-
ods to enforce that the objective function seen by the algo-
rithm always behaves in a convex way in the **α** design
space. Two methods were introduced: a direct convexifica-
tion method and an adaptive convexification method. In
the first one, all hyper-planes are deformed using a given
factor. The values of the objective function, *f*(**α**), for the
values of **α** that satisfy the variable bounds: Σ<sub>*i*=1</sub><sup>*n*<sub>*c*</sub></sup> α<sub>*i*,*j*</sub> = 1
∀*i* = 1, 2,..., *n*<sub>*e*</sub> and **α** ∈ {0, 1}<sup>*n*<sub>*e*</sub>×<sup>*n*<sub>*c*</sub></sup></sup> are not altered, but the
behavior in the intermediate points is more convex. The
expression of the modified objective function, *f*<sup>*</sup>(**α**), fol-
lowing this strategy is

$$
f^*(\alpha) = f(\alpha) + \kappa r(\alpha). \tag{14}
$$

The second method, adaptive convexification, is the result
of a modification of the outer approximation concept, where
the hyper-planes are now computed in a different way to
ensure the approximation is convex. These new hyper-planes
are denoted $\tilde{f}_k^*(\alpha)$ where, as before, $k$ indicates the choice
matrix value on which it is computed and $m$, the points
stored in the memory for the approximation. The expres-
sion of this new hyper-plane is

$$
\tilde{f}_{k}^{*}(\boldsymbol{\alpha}) = f^{*}(\boldsymbol{\alpha}_{k}) + \chi_{k} \cdot (\boldsymbol{\alpha} - \boldsymbol{\alpha}_{k}), \quad (15)
$$

where the slope is computed as
$$
\chi_k = \frac{df^*}{d\alpha}\bigg|_{\alpha=\alpha_k} + \sum_{s=1}^{m} (\alpha_s - \alpha_k)\gamma_{s,k}.
$$
These coefficients $\gamma_{s,k}$ are
computed at each iteration. Their value is so that the approx-
imation satisfies for all already explored points to be an
underestimation. This is done by imposing that

$$
\tilde{f}_k^*(\alpha_l) = \min(f^*(\alpha_l) - \delta, \overline{f}_k^*(\alpha_l)) \\
\qquad \forall l \neq k = 1, 2, \dots, m \tag{16}
$$

with $\delta \ge 0$ being the convexity margin, a constant indicating
the intensity of the convexification effect. This condition
leads to the definition of the following system of equations
to define the $\gamma_{s,k}$ coefficients:

$$
\begin{align}
& \sum_{s=1}^{m} (\alpha_s - \alpha_k) \cdot (\alpha_l - \alpha_k) \gamma_{s,k} \nonumber \\
&= \min(f_l^*(\alpha_l) - f_k^*(\alpha_l)) \nonumber \\
& \quad - \delta - \frac{df^*}{d\alpha} \bigg|_{\alpha=\alpha_k} : (\alpha_l - \alpha_k), 0 \tag{17} \\
& \forall l \neq k = 1, 2, \dots, m. \nonumber
\end{align}
$$

While the α-Under Estimator is fully capable on its own
to render a highly non-convex problem convex enough to
be solved by Bi-Level OA, the same cannot be said for the
Adaptive Convexification. The latest mechanism needs
information from previous points to provide a correct

convexification. This means that to use it efficiently, at least,
the initial point must provide a hyper-plane that underesti-
mates the value of the function. Both of them trigger design
space exploration in the algorithm by providing a more con-
vex objective. This is in both cases related to an increase
in the number of evaluations needed for convergence when
increasing κ or δ.

The optimization problem formulation is introduced in (18).
It involves a catalog, C<sub>S</sub>, that contains n<sub>C<sub>S</sub></sub> choices for the super-
plies. This choice variable is the categorical variable in this
problem. The choice matrix α<sub>S</sub> represents, in this case, which
super-ply is used in each cell. As in the continuous formula-
tion, the thickness of the laminate, t, is our other variable and
is considered continuous. Treating the thickness as a continu-
ous variable greatly facilitates the resolution of the optimiza-
tion problem through relaxation, since the thickness is typi-
cally a sizing variable which allows a given super-ply choice
to become feasible (with respect to the first-ply failure and
buckling constraints). Again, the ceiling function is employed
to post-process the solution to obtain a discrete solution

$$
\left\{
\begin{array}{l}
\min_{t, \alpha_S} m(t) \\
\text{s.t.} \\
\mathbf{g}_{\text{fpf}}(t, \boldsymbol{\alpha}_S) \le \mathbf{0} \\
\mathbf{g}_{\text{shearbuck}}(t, \boldsymbol{\alpha}_S) \le \mathbf{0} \\
\mathbf{g}_{\text{mixbuck}}(t, \boldsymbol{\alpha}_S) \le \mathbf{0} \\
\underline{t} \le t \le \bar{\underline{t}} \\
\sum_{j=1}^{n_{C_S}} \alpha_{S,j} = 1 \quad \forall i = 1, 2, \dots, n_p \\
\boldsymbol{\alpha}_S \in \{0, 1\}^{n_p \times n_{C_S}}
\end{array}
\right.
\tag{18}
$$

3.3 Full-categorical formulation for DD

In this section, we introduce a catalog-based full-categorical
formulation of the problem. An evolutionary method, using
a genetic algorithm (GA) (Golberg 1989), was chosen as the
solver for this formulation.

The formulation of this full-categorical problem is as
follows:

$$
\left\{
\begin{array}{l}
\min_{\alpha_S, \alpha_R} m(\alpha_R) \\
\qquad \text{s.t.} \\
\quad \mathbf{g}_{\text{fpf}}(\alpha_S, \alpha_R) \le \mathbf{0} \\
\quad \mathbf{g}_{\text{shearbuck}}(\alpha_S, \alpha_R) \le \mathbf{0} \\
\quad \mathbf{g}_{\text{mixbuck}}(\alpha_S, \alpha_R) \le \mathbf{0} \\
\sum_{j=1}^{n_{C_S}} \alpha_{S,j} = 1 \quad \forall i = 1, 2, \dots, n_p \\
\qquad \alpha_S \in \{0, 1\}^{n_p \times n_{C_S}} \\
\sum_{j=1}^{n_{C_R}} \alpha_{R,j} = 1 \quad \forall i = 1, 2, \dots, n_p \\
\qquad \alpha_R \in \{0, 1\}^{n_p \times n_{C_R}},
\end{array}
\right.
\tag{19}
$$

where $\alpha_S$ is the choice matrix for the super-plies of each
cell and $C_S$ is the catalog of super-ply choices, of size $n_{C_S}$
as used in the mixed-variable formulation (Problem 18). In
this case, the thickness is led by the number of repetitions of categorical choice, expressed by the variable $\alpha_R$, the choice matrix for the repetitions of each cell. The catalog for this variable, $C_R$, has a size $n_{C_R}$. In this formulation, the total thickness of each cell, $t_i(\alpha_{Ri}) = r_i(\alpha_{Ri})t_{sply}$, where $r_i(\alpha_{Ri})$ is a function that takes the value of the number of repetitions for the i-th cell from the choice matrix and $t_{sply}$ is the thickness of a super-ply.

GA is a well-known evolutionary method that can address categorical problems (Golberg 1989). The implementation of GA used in this analysis was the one available within the Pymoo Python Library (Blank and Deb 2020).

## 3.4 Optimal cell clustering approach for DD

In the previous problem formulations, the optimal solution may include as many different super-plies as the number of cells in the panel. This, though possible, may prove extremely costly to be implemented in actual practice. Too many different orientations in the same panel may result in increased manufacturing costs, whether or not the catalog compliance is mandatory. To mitigate this issue, a clustering optimization strategy is also proposed and implemented in this section. Instead of finding the optimal laminate for each of the cells, the new problem first groups the cells into regions sharing the same super-ply, to be optimized. Then, for each cell, its thickness is also optimized individually. This approach can potentially be implemented to any of the previously detailed formulations. In this work, it was implemented for both the lamination parameter formulation (9), and the mixed-variable formulation (18).

This modification forces the new formulations to have additional categorical variables. A choice matrix $\Gamma$ of size $n_p \times n_{C_Γ}$ is introduced to indicate which cells are associated with the different groups; in the same manner, the $\alpha_S$ variable in Formulation (18) indicates their super-ply choice. These groups of cells or clusters constitute the catalog of clusters, $C_Γ$, of size $n_{C_Γ}$. As in Problem (18), $\Gamma$ is optimized by the solver to find the optimal grouping of cells using the same super-plies. To guide the BLOA to not get stuck in solutions using just fewer regions than provided, a constraint is added to the formulation, $g_{distribution}$. Otherwise, the algorithm may arrive at intermediate choices where one of the regions is empty, and, as there would be no retrieved sensitivities with respect to its related materials, it would be ignored in the future. The expression for $g_{distribution}$ is

$$g_{distribution}(\Gamma) = 1 - \sum_{i} \Gamma_{i,j} \forall j, \qquad (20)$$

where i refers to the element that is being clustered and j to each of the $n_{C_Γ}$ clusters.

When modifying the lamination parameter formulation, the previously presented variables representing the lamination parameters, $\xi_1$ and $\xi_2$, need to have their dimensions modified. In this new formulation, they are no longer assigned to each cell, but to each region. Therefore, each of them has a size equivalent to the imposed number of different regions in the panel, which is $n_{C_Γ}$. The thickness $t$ for the variable remains the same as in the full-continuous Formulation (9). The corresponding problem formulation is provided in (21)

$$ \left\{ \begin{array}{l} \min_{\xi_1, \xi_2, t, \Gamma} m(t) \\ \text{s.t.} \\ g_{f1}(\xi_1, \xi_2) \le 0 \\ g_{f2}(\xi_1, \xi_2) \le 0 \\ g_{f3}(\xi_1, \xi_2) \le 0 \\ g_{fpf}(\xi_1, \xi_2, t, \Gamma) \le 0 \\ g_{shearbuck}(\xi_1, \xi_2, t, \Gamma) \le 0 \\ g_{mixbuck}(\xi_1, \xi_2, t, \Gamma) \le 0 \\ g_{distribution}(\Gamma) \le 0 \\ -1 \le \xi_1 \le 1 \\ -1 \le \xi_2 \le 1 \\ t \le t \le \bar{t} \\ \sum_{j=1}^{n_{C_Γ}} \Gamma_{i,j} = 1 \forall i = 1, 2, \dots, n_p \\ \Gamma \in \{0, 1\}^{n_p \times n_{C_Γ}}. \end{array} \right. \qquad (21) $$

Note that these modifications have made Problem (21) a mixed-categorical one that can be addressed, as in previous examples, using the BLOA solver.

When adding this clustering matrix into the mixed-variable formulation, the definition of the choice matrix $\alpha_S$ needs to be modified. In this implementation, it no longer relates each individual cell with a super-ply from the catalog, but each of these clusters of cells with a catalog choice. As in the previous mixed-variable case, the thickness variable remains a continuous one. The new formulation is expressed in (22)

$$ \left\{ \begin{array}{l} \min_{t, \alpha_S, \Gamma} m(t) \\ \text{s.t.} \\ g_{fpf}(t, \alpha_S, \Gamma) \le 0 \\ g_{shearbuck}(t, \alpha_S, \Gamma) \le 0 \\ g_{mixbuck}(t, \alpha_S, \Gamma) \le 0 \\ g_{distribution}(\Gamma) \le 0 \\ t \le t \le \bar{t} \\ \sum_{j=1}^{n_{C_S}} \alpha_{S_{i,j}} = 1 \forall i = 1, 2, \dots, n_{C_Γ} \\ \alpha_S \in \{0, 1\}^{n_{C_Γ} \times n_{C_S}} \\ \sum_{j=1}^{n_{C_Γ}} \Gamma_{i,j} = 1 \forall i = 1, 2, \dots, n_p \\ \Gamma \in \{0, 1\}^{n_p \times n_{C_Γ}}. \end{array} \right. \qquad (22) $$
## 3.5 Mixed-variable formulation for QUAD

In this section, we seek to extend the previous mixed-variable formulation also to the optimization of QUAD laminates based on catalog choices for the possible ply layups, while simultaneously accounting for manufacturing constraints, notably layup compatibility constraints between adjacent panel cells, which are rarely considered in the existing literature. This optimization problem may seem a priori more prone to be solved via a categorical scheme, as in this case, the catalog choice, representing a stacking sequence, imposes a thickness on the cell. However, the use of the Bi-level OA scheme in this case provides certain advantages: the Bi-level OA is known to be more efficient than heuristic methods, and it can integrate categorical constraints directly in the upper level, which cannot be done with relaxation schemes. This, as later shown, allows the introduction of a mathematical formulation for the compatibility constraint.

Even if using the same solver, the Formulation (18) needs to be modified. Both categorical and continuous variables need to be slightly redefined with respect to the DD case. The continuous variable is the ply thickness of each cell, as now the thickness has to be related to the stacking sequence, and the categorical variable is the stacking sequence for each cell. A catalog of stacking sequences is defined beforehand. This affects the formulation in two ways: there must be a relation between the chosen stacking sequences and the thicknesses of the cells, and the solution must satisfy the manufacturing constraints, including compatibility constraints.

Instead of dealing with the thickness of each cell, $t_i$, the ply thickness, $t_{ply}$, will be used as an optimization variable. If the number of plies of a given cell is a function of the choice of stacking sequence for that cell, now indicated with the choice matrix $\alpha_L$, such as $n_{plyi}(\alpha_{Li})$, the relation between ply and laminate thickness for the i-th cell can be expressed as

$$t_i = t_{ply} n_{plyi}(\alpha_{Li}). \qquad (23)$$

To ensure that the optimizer finds a solution where the discrepancy between the computed ply thickness and the objective one, $t_{plytrue}$, which is a parameter for the optimization, a penalty is added to the objective. This term is computed as proportional to the discrepancy, $\epsilon$, expressed as

$$\epsilon = \left( \frac{t_{ply} - t_{plytrue}}{t_{plytrue}} \right)^2. \qquad (24)$$

Regarding the manufacturing constraints, as previously described, these can be divided into two groups: the laminate-related and compatibility-related constraints. The first group is already taken into account in the mixed-variable formulation through the catalog definition. All possible

choices in the catalog must be feasible from a manufacturing point of view. The second group is dealt with through a compatibility constraint. This group contains the blending constraint, ensuring that all plies in the thinner of two contiguous cells continue in the other one; and drop-off, limiting the difference in number of plies among adjacent cells. To build this constraint, a pre-optimization analysis of the catalog is performed. This analysis compares the different choices and identifies which ones are compatible. From this analysis, a square and symmetric matrix referred to as the compatibility matrix, $C$, is generated, such that

$$C_{i,j} = \begin{cases} 1 & \text{if } i \text{ and } j \text{ choices are compatible} \\ 0 & \text{otherwise.} \end{cases} \qquad (25)$$

This allows the formulation of a quadratic constraint for the upper level of BLOA. For each individual cell, a number $n_{adj}$ of contiguous cells can be identified. For a given material choice, it is possible to identify both the material choice corresponding to the cell in the vector $\alpha_{Li}$ and the material choices for the adjacent cells in the matrix $\alpha_{Ladj}$. With this, it is possible to define the constraint for the compatibility as

$$g_{comp} = \sum_{i=1}^{n_p} (n_{adj_i} - \sum_{j}^{n_{adj_i}} \alpha_{Li} C \alpha_{Lj}^T) \leq 0, \qquad (26)$$

where, if cells i and j are compatible, the term in the second addition is equal to 1, and 0 otherwise. This way, all the adjacent cells for every cell in the structure must be compatible to satisfy this constraint. This constraint is linearized in the Bi-Level OA outer level.

With these modifications in place, the formulation of the MVF problem for a QUAD structure is

$$\left\{ \begin{array}{l} \min_{t_{ply}, \alpha_L} m(t_{ply}, \alpha_L) + \lambda \epsilon(t_{ply}) \\ \text{s.t.} \\ \mathbf{g}_{fpf}(t_{ply}, \alpha_L) \leq \mathbf{0} \\ \mathbf{g}_{shearbuck}(t_{ply}, \alpha_L) \leq \mathbf{0} \\ \mathbf{g}_{mixbuck}(t_{ply}, \alpha_L) \leq \mathbf{0} \\ \mathbf{g}_{comp}(\alpha_L) \leq \mathbf{0} \\ \underline{t_{ply}} \leq t_{ply} \leq \overline{t_{ply}} \\ \sum_{j=1}^{n_{c_L}} \alpha_{Li,j} = 1 \forall i = 1, 2, \dots, n_p \\ \alpha_L \in \{0, 1\}^{n_p \times n_{c_L}}, \end{array} \right. \qquad (27)$$

where $\lambda$ is the penalty factor for the thickness discrepancy and $C_L$ is the catalog of stacking sequences of size $n_{c_L}$. It is important to remark here that, as in classical penalization methods (e.g., Bendsoe 1989), a progressively higher penalty factor leads to a better value of the discrepancy in the result. This is the reason why, in this case, the formulation was solved iteratively with different values of $\lambda$.
## 4 Results

In this section, the optimization results obtained for the different formulations shown in Sect. 3 for the composite panel use-case are provided. Two different configurations of the problem presented in Sect. 2 were solved to investigate the influence of the number of variables on the performance of the different solvers, one using DD laminates and the other using QUAD laminates. As explained in the previous section, the DD case was used to compare the mixed-variable formulation (MVF) (18), using the BLOA as the solver against the lamination parameter formulation (LPF) (9), and the full categorical formulation (CaF) (19), and show the benefits of variability control schemes in formulations (21) and (22). The second case should serve as an example of how to apply MVF to QUAD structures, with formulations (27), and allow a comparison between these laminate concepts and their associated manufacturing techniques. The clustering strategies were also tested in these two configurations. In addition, the implementation of the Bi-level OA and MVF in an industrial use-case is also presented via a vertical stabilizer spar optimization.

In both configurations, the overall size of the panel, $H = 6$ m and $L = 10$ m, is kept constant. The loading is the main difference between the two configurations. Both are divided into 50 cells, i.e., $n_{px} = 5$ and $n_{py} = 10$. In the DD case, the applied force is 500 kN and the values of the stiffener properties were $S_s = 400$ mm² and $E_s = 43.7$ GPa. These values were chosen to provide stiffeners with approximately 10% of the area of the average panel and approximately the same stiffness. As the applicability of the QUAD composites is defined by the material and stacking sequence catalog used, the use-case has to be adapted to our QUAD catalog. In the QUAD case, the applied force is 20 kN and the values of the stiffener properties were $S_s = 100$ mm² and $E_s = 10$ GPa. Table 1 summarizes the naming of each of the strategies studied in this section. An extra formulation, named LPF(1DD), is added to the comparison to serve as a reference for assessing

the advantage of using different super-plies in each cell. This optimization will be done using the LPF previously explained, but forcing all cells to share the same super-ply and thickness.

The uni-directional ply properties used to compute the mechanical properties of the panel are provided in Table 7 in Appendix A.1, representing those of IM7/977-3 composite material as defined in Tsai et al. (2017). The thickness of each super-ply, containing 4 plies of the same thickness, $t_{sply}$, was 0.125 mm and the density of the material $\rho = 1.4 \times 10^{-6}$ kg/mm³.

For the DD, a catalog for the super-plies, $C_s$, needed to be defined. The catalog used in this analysis contained all possible super-ply combinations when mixing orientations from 0° to 90° with a step of 15°. This makes a total of $n_{C_s} = 28$ choices in this catalog, taking into account that property-wise there is no difference between super-plies [$\Phi, \Psi$] and [$\Psi, \Phi$]. These super-plies are shown in Table 2. This catalog of super-plies has been selected by testing different values of this sweep and analyzing how close the solutions of the MVF become to those of the LPF. This approach aimed to find a good compromise between the number of choices in the catalog and the proximity between the solutions of both methods. In an industrial application, it is important to reduce the size of the catalog to control the quality and cost of the choices. This catalog contains the matrices $A^*$ and $D^*$ of every super-ply for a homogenized configuration consisting of 15 repetitions of the super-ply, leading to a total thickness of 1.875 mm. These quantities are constant with respect to the thickness or number of repetitions of the laminate as long as the number of repetitions is higher than 4 or 6, i.e., the laminate can be considered homogenized. This condition was imposed as the lower bound for the thickness in the formulations. Additionally, the CaF formulation needs a catalog for the number of repetitions, $C_R$. The chosen catalog is a list of all discrete values between 4 and 70, meaning that $n_{C_R} = 66$. In the case of QUAD, the catalog contains a set of possible stacking sequences extracted from previous industrial applications, which is detailed in the Appendix, C, and the common 0.127 mm plies with the same mechanical

<table><caption>**Table 1** Summary of formulation naming (acronyms) and description (equations, solvers, and major remarks)</caption><thead><tr><th>Name</th><th>Acronym</th><th>Formulation</th><th>Solver</th><th>Remarks</th></tr></thead><tbody><tr><td>Lamination parameters' formulation</td><td>LPF</td><td>Eqs. (9)</td><td>MMA</td><td>–</td></tr><tr><td>Lamination parameters' formulation with a unique DD for all cells</td><td>LPF(1 DD)</td><td>Eqs. (9)</td><td>MMA</td><td>A single DD laminate is used for all cells</td></tr><tr><td>Mixed-variable formulation</td><td>MVF</td><td>Eqs. (18)</td><td>Bi-level OA + MMA</td><td>–</td></tr><tr><td>Categorical formulation</td><td>CaF</td><td>Eqs. (19)</td><td>GA</td><td>–</td></tr><tr><td>Lamination parameters' formulation with optimal clustering</td><td>LPF + CL#</td><td>Eqs. (21)</td><td>Bi-level OA + MMA</td><td>Implements a clustering of # regions</td></tr><tr><td>Mixed-variable formulation with optimal clustering</td><td>MVF + CL#</td><td>Eqs. (22)</td><td>Bi-level OA + MMA</td><td>Implements a clustering of # regions</td></tr></tbody></table>
<table><caption>**Table 2** DD super-plies catalog represented by the combinations of Φ and Ψ expressed in degrees</caption><tr><td>0,0</td><td>15,15</td><td>30,30</td><td>45,45</td><td>60,60</td><td>75,75</td><td>90,90</td></tr><tr><td>0,15</td><td>15,30</td><td>30,45</td><td>45,60</td><td>60,75</td><td>75,90</td><td></td></tr><tr><td>0,30</td><td>15,45</td><td>30,60</td><td>45,75</td><td>60,90</td><td></td><td></td></tr><tr><td>0,45</td><td>15,60</td><td>30,75</td><td>45,90</td><td></td><td></td><td></td></tr><tr><td>0,60</td><td>15,75</td><td>30,90</td><td></td><td></td><td></td><td></td></tr><tr><td>0,75</td><td>15,90</td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>0,90</td><td></td><td></td><td></td><td></td><td></td><td></td></tr></table>

properties as those used in DD. The compatibility matrix is computed for this catalog, accounting for blending and a maximum of 8 plies total difference between adjacent laminates. It is important to note that this particular catalog contains non-symmetric choices. This is due to the fact that these pseudo-symmetric choices behave close to symmetric ones, and previous experiences proved their utility.

All formulations were coded using the GEMSEO library for Python (Gallard et al. 2018), which contains the MMA (Svanberg 1987), and BLOA solvers and provides an environment where multidisciplinary problems can easily be implemented and solved through different strategies.

An important parameter to track in this analysis was the number of stress analyses (cf. Sect. 2) that were performed by the solver. In an industrial application, this typically involves a Finite-Element Analysis (FEA). The efficiency of a solver is related, in these applications, to the number of calls to this model that it requires to converge.

The BLOA algorithm hyperparameters were tuned in each case to obtain the lightest solution. In particular, for some tests of the DD configuration, the adaptive convexification (see Sect. 3) was used to foster design space exploration in the DD case and the conventional convexification in the QUAD panel case. The values for the convexity margin, δ, and convexity constant, κ, were chosen after a manual tuning. In terms of the solver for the subproblem, MMA was used. The results are identified in the next paragraphs as MVF. For the relaxed formulation, the MMA solver was used as well, and these solutions are identified as LPF. In the case of the genetic algorithm, CaF, the Pymoo (Blank and Deb 2020) library GA solver was used, and the solver was run several times over the problem with the same hyperparameters. The number of evaluations for each run is specified through the population size and the number of generations, and then, this run is repeated to reduce the statistical variability.

## 4.1 DD panel

The optimum masses obtained by each of the solvers and formulations are shown in Table 3, with the number of model evaluations in each of the cases. In this particular application, the time needed for an evaluation was around

0.2 s, including the retrieval of the derivatives (on a PC with AMD Ryzen 5 PRO 5650U processor and 16GB of RAM). As different formulations propose different grades of relaxation, some transformations are needed to obtain a manufacturable solution.

In these tests, all solvers started from the super-ply [45,45] and the maximum laminate thickness. The lower and upper bounds for the thickness were 0.5 and 100 mm. LPF and MVF, in this case with $κ = 20$, provided directly the results shown in Figs. 3 and 4. In both cases, the algorithms arrive at the convergence criteria, variation of the objective function for the LPF, and collapse of the upper and lower bound for MVF. This remains true for the rest of the results in this paper. These are identified as relaxed results, as they include a relaxed value for the thickness, as directly obtained by the algorithm. These values can be rounded up to the closest discrete thickness value compliant with a given number of super-plies, increasing the mass conservatively, following the expression (8). For the MVF results, this led directly to a Catalog-Compliant (CC) design with a mass of 407.34 kg. On the other hand, the solution provided by LPF presented super-plies that were not included in $C_S$. Discretizing the thicknesses respecting these super-plies led to a No-catalog Compliant (No-CC) design with a mass of 401.52 kg. The small difference in mass proves that even if the continuous formulation can find better super-plies than those present in the catalog $C_S$, the gains in terms of mass might not be as substantial to compensate for the

<table><caption>**Table 3** Optimization results in DD panel case. The relaxed optimum mass refers to the mass of the solution with a continuous value of the thickness. The discretized optima are either CC or No-CC. CC indicates that the solution is built only from catalog choices, while No-CC states the opposite</caption><thead><tr><th>Formulation</th><th>Relaxed<br>optimum<br>mass</th><th>No-CC<br>optimum<br>mass</th><th>CC optimum<br>mass</th><th>Evaluations</th></tr></thead><tbody><tr><td>LPF(1 DD)</td><td>826.38 kg</td><td>829.75 kg</td><td>829.75 kg</td><td>28</td></tr><tr><td>LPF</td><td>396.5 kg</td><td>401.52 kg</td><td>—</td><td>85</td></tr><tr><td>MVF</td><td>403.4 kg</td><td>—</td><td>407.34 kg</td><td>31,999</td></tr><tr><td>CaF (best<br>run)</td><td>—</td><td>—</td><td>431.14 kg</td><td>160,000*</td></tr></tbody></table>
<figure><img src="image_6.png" /><figcaption>**Fig. 4** Super-ply and cell thickness distribution for the DD panel using the MVF formulation. The solution mass is 403.4 kg</figcaption></figure>

<figure><img src="image_7.png" /></figure>

cost increase related to using non-standardized materials. In terms of evaluations, the continuous solver is significantly more efficient than the mixed-categorical one. A discussion about the scalability of BLOA in comparison with other methods can be found in Barjhoux (2020). Finally, the CaF was unable to retrieve the same or a similar design to the others, using the same super-plies as shown in Fig. 5, even at the cost of significantly more evaluations, which highlights how challenging the optimization problem is. These evaluations accounted for a population size of 400 and 400 generations. The results of GA shown in Table 3 represent the best of the 50 repetitions in terms of mass, 431.14 kg. The mean of the runs was 437.44 kg and the standard deviation of the distribution was 3.08.

To quantify the gains obtained by defining different laminates in each cell of the panel, a reference panel was defined. This reference panel (LPF(1 DD)) was computed by introducing the continuous formulation constraints to impose the same super-ply and laminate thickness in all cells. It is possible to see that the mass of this reference is more than twice that of the results of LPF and MVF.

Note that all our formulations can appropriately handle load redistribution effects and buckling constraints. This represents a step forward compared to existing DD optimization approaches in the literature (e.g., Garofano et al. 2023; Wang et al. 2023), which only consider strength constraints that do not need to account for load redistribution.

The results of implementing the optimal cell clustering scheme will be presented next and compared to the non-clustered MVF and LPF results. To reduce the effect on the manufacturing cost, the maximum number of different super-plies in the optimal design is fixed before the optimization.

For the previously shown configuration, different numbers of clusters have been tested. This number was fixed to 1, 2, and 3, indicated as CL1, CL2, and CL3, respectively. The results in terms of mass and evaluations are shown in Tables 4 and 5. As in the previous MVF solution, adaptive convexification was used to aid in the convergence and was also fine-tuned. For the clustering based on LPF, $\delta = 20$ was selected, and for those based on MVF, $\delta = 30$ was chosen. The reader may realize that the restrictions imposed by the clustering do not just reduce the variability by removing
<figure><img src="image_8.png" /><figcaption>**Table 4** Clustering results in the DD panel case with LPF</figcaption></figure>

<table><caption>**Table 5** Clustering results in the DD panel case with MVF</caption><thead><tr><th>Formulation</th><th>Relaxed optimum mass</th><th>No-CC optimum mass</th><th>Evaluations</th></tr></thead><tbody><tr><td>LPF</td><td>396.5 kg</td><td>401.52 kg</td><td>85</td></tr><tr><td>LPF+CL3</td><td>407.72 kg</td><td>412.66 kg</td><td>5300</td></tr><tr><td>LPF+CL2</td><td>422.52 kg</td><td>427.78 kg</td><td>625</td></tr><tr><td>LPF+CL1</td><td>448.28 kg</td><td>453.82</td><td>17</td></tr></tbody></table>

<table><thead><tr><th>Formulation</th><th>Relaxed optimum mass</th><th>CC optimum mass</th><th>Evaluations</th></tr></thead><tbody><tr><td>MVF</td><td>403.4 kg</td><td>407.34 kg</td><td>31,999</td></tr><tr><td>MVF+CL3</td><td>403.4 kg</td><td>408.88 kg</td><td>14,366</td></tr><tr><td>MVF+CL2</td><td>415.76 kg</td><td>421.48 kg</td><td>3227</td></tr><tr><td>MVF+CL1</td><td>448.76 kg</td><td>453.6 kg</td><td>70</td></tr></tbody></table>

choices from the final result while slightly increasing the mass, but also, in combination with the BLOA algorithm, reduce the number of evaluations. As this algorithm can become more costly in more complex cases, even if less than other mixed-categorical choices for the same accuracy, this modification in the formulation may be beneficial in cases where a Catalog-Compliant solution is required. However, tests in more complex problems may be required to validate such affirmation. The solution for a 3-region clustering applied to the MVF formulation is shown in Fig. 6.

As expected, the clustering with just one choice for all cells resulted in heavier panels than almost all other solutions found for this problem. However, they were still almost two times lighter than the reference case [i.e., LPF(1DD)], proving that a substantial part of the gains in this optimization came from the choice of an optimal thickness for each cell. From there, when increasing the number of regions

in the panel, the mass decreased. In a real application, this would also generate an increase in terms of manufacturing cost, leading to a trade-off between mass and cost. This clustering approach provided a simple way to control the variability in this application to help make this trade-off.

When comparing the solutions with CL1 and those obtained by the LPF(1DD), a remark regarding the optimization of DD structures can be drawn. It is possible to optimize a structure using only one super-ply definition, potentially leading to a lower manufacturing cost. This solution is simpler to produce and design than that obtained with QUAD laminates, as no blending constraint needs to be addressed to ensure compatibility among cells. Optimizing the number of repetitions per cell and the overall super-ply resulted in a 45% reduction in mass with respect to the LPF(1DD) baseline. If the optimization also considers the optimization of super-plies at each cell, the gains in mass can arrive at 50%.

## 4.2 QUAD panel

As the previous results showed, even if the MVF was not the most efficient in the DD case, it retrieved catalog-compliant solutions that were very close to the relaxed ones. After also introducing manufacturing constraints and in particular compatibility constraints as proposed in formulation (27), we can also solve realistic QUAD laminates optimization problems. The results for this case using the MVF with and without the compatibility constraint are shown in Table 6. A reminder here that the loadings, in this case, are not the same as in the previous DD cases due to the differences in available catalogs. For the optimizations associated with these results, a convexity constant $κ = 50$ was found to be the most efficient in terms of computational cost and exploration. Also, regarding the values of the penalty factor $λ$, a series of two values was used to ease convergence as explained in Sect. 3, 0.5 and 1. The initial guess for the choice was the first option
<figure><img src="image_9.png" /><figcaption>**Table 6** MVF results for the QUAD panel with and without compatibility constraint</figcaption></figure>

<table><thead><tr><th>Compatibility con-<br>straint</th><th>Relaxed<br>optimum<br>mass</th><th>CC optimum mass</th><th>Evaluations</th></tr></thead><tbody><tr><td>g<sub>comp</sub> active</td><td>165.24 kg</td><td>162.14 kg</td><td>3942</td></tr><tr><td>g<sub>comp</sub> not active</td><td>164.36 kg</td><td>158.96 kg</td><td>8737</td></tr></tbody></table>

of the catalog and the maximum possible ply thickness in all cells. The bounds for the ply thickness were 0.1 and 10 mm.

As expected, a restriction in the catalog choices, through the blending and drop-off constraint, leads to a higher optimal mass, but the relative difference is quite small. The significant differences between the chosen QUAD laminates when accounting for or not for compatibility constraints show that the problem of accounting for them is not trivial and cannot be simply solved by a post-processing step starting from the solution without compatibility constraints. Note also that the computational costs between the two solutions are quite different and are related to how the compatibility constraint is implemented. As it affects the solution of the MILP outer loop problem, it reduces or limits the potential solutions seen by the solver, reducing the exploration. This is seen in Table 6 as a reduction of the number of evaluations. The optimization without the blending constraint led to some inconsistencies, as several incompatible stacking sequences were placed next to each other. In the optimization with the blending constraint, this situation does not occur. All cells are blending compatible, as shown in Fig. 7. In Fig. 8a, b, transversal cuts at the positions $x = 1000$ mm and $y = 250$ mm are provided. These cuts show the continuity of all plies in the thinner cells in the contiguous ones. The same can be said of any other cut in the solution. Another important remark from these figures is the simultaneous drop-offs from the symmetry axis and from intermediate positions of the different plies. These drop-offs cannot be

managed with other gradient-based methods from the literature, proving another advantage of the catalog-based approach. It is relevant to highlight that the shown connections are just one of the several possibilities to represent a blended configuration among the given laminates. Several other configurations equally compatible with the blending constraint can be identified in both cases. The aim of the current algorithm does not revolve around these connections as a result of the optimization, but on the fact that there are feasible connections.

Finally, we carried out a direct comparison between DD and QUAD laminates for the latter use-case. A new design with the same geometry, stiffeners, and loading was optimized with the DD catalog using a super-ply as thick as the individual plies of the QUAD. The results are shown in Fig. 9a. The DD solution has a mass of only 139.54 kg, once discretized, which shows that DD leads to a significantly lighter panel than QUAD for this particular case. For this particular result, it is important to mention that compatibility constraints were considered for the QUAD but not for the DD case, considering claims by Tsai et al. (2017) that compatibility constraints are not necessary for DD. This claim would have to be validated experimentally in this case before the DD solution can indeed be considered the lightest, all else being equal. To provide an easier comparison, the panel was also optimized using MVF+CL1, shown in Fig. 9b. This resulted in a smaller mass reduction, to 159.18 kg, but in a configuration which should be simpler to manufacture, as no blending among different super-plies is concerned. Finally, note that some additional insights into the results obtained in this section are provided in Appendix D.

## 4.3 Industrial use-case: short-range aircraft vertical tail plane front spar

In addition to the previously shown work over the multi-cell stiffened panel, which can be considered representative of a
<figure><img src="image_10.png" /><figcaption>**Fig. 7** Stacking and cell thickness distribution for the QUAD panel with blending constraint. The solution mass is 165.24 kg. For stacking sequence details, see Appendix C</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>(a) One of the possible schemes of the connections among plies in a transversal vertical cut at $x = 1000\,\text{mm}$ of the blending compatible solution.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>(b) One of the possible schemes of the connections among plies in a transversal horizontal cut at $y = 250\,\text{mm}$ of the blending compatible solution.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>(a) Super-ply and cell thickness distribution for the QUAD use case using DD materials and the MVF formulation. The solution mass is 139.54kg.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>(b) Super-ply and cell thickness distribution for the QUAD use case using DD materials and the MVF+CL1 formulation. The solution mass is 159.18kg.</figcaption></figure>

<caption>**Fig. 9** Solutions to the QUAD use-case using DD materials</caption>

fuselage panel, a more comprehensive industrial application case, using actual industrial processes, was analyzed to show how the MVF scheme can be easily integrated into these.

The addressed problem is the design of a multi-cell spar that resembles that of a short-range aircraft Vertical Tail Plane (VTP). The structure is analyzed using the Next
Generation Structural Analysis (NGSA) tool from Airbus, which is able to combine FEA, analytical, and empirical formulae to provide the different reserve factors of the structure. The spar is divided into 9 panels and is supposed to be clamped to the left side (cf. Fig. 10). The tool allows the modification of this structure by introducing manholes and stiffeners, which are taken into account in the computation of the reserve factors of each panel. The addition of stiffeners divides each panel into a given number of subpanels. The optimization is performed at the subpanel level, equivalent to the cells in the previous case. The material choice and thickness can be modified for each of them. As shown in Fig. 10, a configuration with 11 subpanels was obtained by placing 2 stiffeners in the first panel from the left of the figure (see the red lines in the figure). In the first 3 panels, different manholes were also placed. Note that the rectangular hole at 45° represents the modeling of actual manholes of more complex geometries. Similarly to the previously presented use-case, the objective of this optimization is to define the design that minimizes the mass of the spar while satisfying relevant structural constraints. As in the academic case, the stiffeners' definition is assumed fixed. In this case, the structural constraints are computed by NGSA, which is also able to retrieve the sensitivities of these reserve factors with respect to the ply orientations. In addition to the global and buckling constraints, which experience showed they were dominant over other strength criteria for this application, the solver is able to introduce a damage tolerance criterion in each panel and take into account the manholes located in the first three panels from the left in this damage assessment as an empirical modification over the classical strength constraints. This damage tolerance computation to account for free edge impacts or impacts on the interior of the manholes is empirically derived from experimental campaigns performed by Airbus.

The configuration shown in Fig. 10 was optimized considering QUAD and DD materials. The first investigation evaluated the ability of the MVF formulation to find a QUAD-based solution that accounts for all the failure criteria of the tool, using the same catalog as in previous examples (C), but removing those stacking sequences where one orientation is present in less than 8% of the plies, to adhere to this more strict industrial constraint. To account for the discrepancy between the relaxed and discrete formulation, a knock-down factor $k_D \le 1$ was added to the formulation

<figure><img src="image_15.png" /><figcaption>Fig. 10 Spar scheme. The black lines separate the panels, and the red ones indicate the position of the stiffeners and separate the subpanels</figcaption></figure>

for QUAD, otherwise identical to the formulation (27). The second investigation sought a DD solution to compare it with QUAD to evaluate the gains.

The MVF formulation, a convexification factor $kappa = 30$, a set of two penalties, 0.5 and 1, and an interval of 5 iterations without improvement to reduce the trust region was applied to this problem. The solver was able to retrieve a blending-compatible design of 30.06 kg. Similarly to the previous panel, the cells closer to the root present higher thickness, while the presence of holes in the second, fourth, and fifth cells also made them less stiff by design, requiring more thickness than those around them. This is also influenced by the damage tolerance criteria. Sudden changes in the optimal thickness are found in certain areas to cope with different loadings, which are not as linear as in the academic use-case. This is not a failure of the optimizer or a bad solution, as it helps to minimize the mass and is blending compatible, but this sudden change may increase assembly costs. To solve this particular problem, an extra constraint was introduced, profiting from the capability of BLOA to manage categorical constraints. This constraint imposes a continuous evolution of the thickness in the spar, forcing the algorithm to select categorical choices, such that the evolution of the number of plies along the horizontal axis does not present sudden changes. It is expressed in the following formula:

$$ 
\begin{aligned}
g_{\text{slope}} &= -(n_{\text{plyi}}(\alpha_L) - n_{\text{plyi+n}_y}(\alpha_L)) \\
(n_{\text{plyi+n}_y}(\alpha_L) - n_{\text{plyi+2n}_y}(\alpha_L)) \le 0 \\
\forall i = 1, \dots, n_p - 2 * n_y,
\end{aligned}
\quad (28) $$

where $i$ is a given subpanel of the total $n_p$ and $n_y$ is the number of subpanels in a vertical cut. This constraint, in addition to $g_{\text{comp}}$, ensures that the solution presents a smoother thickness transition. The result, with this constraint taken into account, is shown in Fig. 11. As the new design space is more restrictive, the optimal discretized mass was increased to 31.74 kg. The blending constraint is still respected as shown in Fig. 12.

The MVF for the DD case was launched from an initial [90,90] super-ply with maximum thickness in all the cells, an adaptive convexification margin of $\delta = 30$, and the same trust region reduction scheme as in the QUAD case, but with a minimum of just 1 change in the choice matrix. As expected for such loading and considering only buckling constraints, global and local, the super-ply [45,45] was selected for all the cells. Aside from the mass reduction from 30.06 to 26.52 kg, the DD configuration also presents a simpler to build VTP Spar where all cells can be produced from the same super-ply using the card-slicing technique mentioned in the bibliography. In this case, this was a result of the conventional formulation, but
it can also be enforced using the optimal cell clustering approach if needed.

Finally, note that, as for the DD case, some additional insights into the results obtained in this section are provided in Appendix D.

# 5 Conclusions and perspectives

In this article, a catalog-based mixed-variable bi-level optimization framework for composite materials was presented, the BLOA. This method is based on the Bi-Level Outer Approximation solver and was shown capable of identifying the optimal choice out of a stacking sequence catalog in several use cases, both considering QUAD and DD laminate concepts.

<figure><img src="image_16.png" /><figcaption>**Fig. 11** Stacking and cell thickness distribution for the 11 subpanels Spar in QUAD with slope constraint $g_{slope}$. The solution mass is 31.74 kg. For stacking sequence details, see Appendix C</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>**Fig. 12** One of the possible schemes of the connections among plies for the 11 subpanel Spar case with slope constraint $g_{slope}$ active</figcaption></figure>
The proposed method included two convexification strategies that extend its applicability. These strategies trigger the exploration of the design space even when it is non-convex. However, this overexploration also leads to an increase in the computational time. Furthermore, these strategies define from which initial guesses the algorithm is able to identify the global optimum.

To further improve the manufacturability and costs of the optimal solutions with respect to stacking sequence layup variability from cell to cell, an optimal clustering strategy was explored. Both the parametrization and the mixed-variable methods were tested in combination with it, achieving the desired reduction of variability in the results and additionally reducing the computational cost of the mixed-variable formulation. Nevertheless, there are trade-offs to be considered when using this method between the number of materials present in the solution and the optimal mass.

The performance of the proposed approaches was tested against standard optimization techniques in the frame of the optimization of a DD multi-cell stiffened panel. Its accuracy was comparable to that of classical parametrization approaches, while the combinatorial exploration was better than that of classical metaheuristic approaches. At the same time, it was possible to see that while the proposed algorithm provided a directly manufacturable solution, the classical parametrization method allowed for a quicker design from which a post-process can be launched to identify the manufacturable solution. Also, the classical methods are more suited for exploration outside of catalog choices, when the freedom to do so is present.

This mixed-variable strategy was also extended to QUAD laminates, for which a compatibility constraint based on a catalog pre-processing was introduced. The use of the BLOA and this constraint assured that the provided solutions are compatible in terms of blending and drop-offs at the same time, as they satisfy the structural constraints. In addition, this strategy allowed a comparison between DD and QUAD, where it was shown that the former can achieve lighter structures both when considering a simple to manufacture single super-ply structure and even more when considering more complex multi-super-ply structures in the cases where the nature of the loading varies in the different cells. These results were consolidated in an industrial application to a vertical tail plane spar. At the same time, it is important to note that in this work, no blending rules were considered for the DD laminates due to low maturity of DD and lack of consensus, which may render some DD optima overop timistic. Nevertheless, even when the same DD super-ply was considered over the entire panel and only thickness was varied, the optimal DD laminate still showed lower mass than the optimal, manufacturable QUAD solution.

While the proposed method is more efficient and shows better scalability when solving the combinatorial problems

considered here than classical metaheuristics, it can still be challenging to address large industrial problems, and more refining is needed. Challenges arise in particular in the handling of the large number of hyper-planes associated with larger number of design variables. On the other hand, the method is adaptable to different types of problems due to its mixed-variables structure, as shown in this work by addressing different manufacturing techniques. In the future, its extension to other composite material concepts (e.g., curve-fibers, tri-axial or quasi-trivial laminates) will be worth exploring. Furthermore, the modifications made to the BLOA algorithm allow it to consider a wide range of constraints, not needing to build a convex categorical problem. In this regard, further constraints, such as fatigue, damage, crashworthiness, dynamic responses, or thermal behavior, could be integrated in the future. In addition, a comparison against other state-of-the-art formulations for the QUAD formulation should be considered to robustify the results of this paper.

## Appendix A: Double-double composites

This section presents the methods used to compute the properties of double-double (DD) laminates and the characteristics of the mechanical properties of homogenized DDs.

### Homogenization

One of the claimed advantages of the DDs in comparison to QUADs is that the thickness of homogenized laminates can be smaller than in the case of QUAD composites. Homogeneity is considered achieved when the bending and stiffness properties of a laminate are decoupled (Tsai et al. 2017).

The conditions for homogeneity of DD in terms of components of the stiffness matrix are as defined in Tsai (2021)

$$ 
\left\{ 
\begin{array}{l} 
B_{i,j}^* / \text{trace}(\mathbf{A}^*) \le 0.02 \quad \forall i, j \\ 
A_{i,j}^* - D_{ij}^* / \text{trace}(\mathbf{A}^*) \le 0.02 \quad \forall i, j \\ 
D_{1,6}^* / \text{trace}(\mathbf{A}^*) \le 0.02 \\ 
D_{2,6}^* / \text{trace}(\mathbf{A}^*) \le 0.02. 
\end{array} 
\right. 
\qquad (29) 
$$

The * in the equations (29) denotes that these matrices were normalized with the thickness of the laminate, t, as

$$ A^* = \frac{1}{t} A, B^* = \frac{2}{t^2} B, D^* = \frac{12}{t^3} D. \qquad (30) $$

Using the CLT to compute the laminate properties of different super-plies for varying numbers of repetitions, the number of plies to achieve homogenization for different super-plies can be computed. For this analysis, as well as for the proposed optimization, the ply properties are those
<table><caption>**Table 7** IM7/977-3 mechanical properties as extracted from Tsai et al. (2017)</caption><thead><tr><th>E<sub>x</sub>[GPa]</th><th>E<sub>y</sub>[GPa]</th><th>G[GPa]</th><th>ν[−]</th><th>X<sub>c</sub>[MPa]</th><th>X<sub>t</sub>[MPa]</th><th>Y<sub>t</sub>[MPa]</th><th>Y<sub>c</sub>[MPa]</th><th>S[MPa]</th></tr></thead><tbody><tr><td>191</td><td>9.9</td><td>7.8</td><td>0.35</td><td>3250</td><td>1600</td><td>62</td><td>98</td><td>75</td></tr></tbody></table>

<figure><img src="image_18.png" /><figcaption>**Fig. 13** Repetitions for homogenized material according to expression (29)</figcaption></figure>

proposed in Tsai et al. (2017) and summarized in Table 7 for the IM7/977-3 composite material.

When analyzing the number of repetitions needed to homogenize a given super-ply, according to (29), it presents a symmetric behavior regarding the orientations Ψ and Φ as can be seen in Fig. 13. Note also that there is a trend to increase the number of repetitions as the difference between the positive and negative orientations in the super-ply grows. In the super-ply [0,15], the maximum difference in orientations is 30° (between the negative and positive 15° orientations), and it is homogenized with 2 repetitions, while for the super-ply [0,75], this difference is up to 150° and it is homogenized with 6 repetitions. The average number of repetitions to have a homogenized material is 4, meaning that the homogenized laminate has 16 plies. It is possible to notice that this homogeneity is, for most super-plies, an approximation based on the threshold defined in Tsai (2021). Mathematical homogeneity is achieved only in symmetrical DDs, as all of them are already balanced. For many industry applications, the approximate homogeneity provided by the threshold appears however to be satisfactory.

## DD homogenized properties study

Once the material is homogenized, the mechanical properties of a laminate with a given super-ply are defined by its number of repetitions. The thickness-normalized properties only depend on the super-ply choice. This means that for a given

loading, the super-ply choice defines these properties and, therefore, the laminate thickness needed to sustain that load. A simple parametric study showed that the loading is enough to define the optimal super-ply and number of repetitions. This study was performed by evaluating the material properties, regarding pure compression and pure shear buckling behaviors for different numbers of repetitions of each super-ply. CLT was used in this analysis to compute the material properties. Compression and shear buckling reserve factors were computed using the expressions from Herencia et al. (2007). The reserve factors for a supported 0.6 × 10 m² panel of each of the super-plies used in the previous analysis were computed. The loading for pure shear was a flux of N<sub>xy</sub> = 25 N/mm and that for pure compression was N<sub>x</sub> = 167 N/mm. Figure 14a, b show the number of repetitions needed to obtain a greater than one reserve factor in the cases of pure shear and pure compression respectively.

## Appendix B: Analytical expressions for the stress redistribution

The analytical formulation of the stress redistribution was based on the computation of the moments of inertia, accounting for the thickness and material variations. This formulation is based on beam stresses computation and it uses the thickness and stiffness of each cell to compute the stress distribution on a given transversal section, being a transversal section the set of cells that share the same x-coordinates. With this information from each cell forming the section, its center of inertia, and the first and second moments of inertia are computed. These are then used to compute the stresses along the section.

The transversal area of the j-th element of the i-th section is computed as a multiplication of its height, H/n<sub>py</sub>, times its thickness t<sub>i,j</sub>

$$S_{i,j} = \frac{H}{n_{py}} t_{i,j} \qquad (31)$$

The stiffness, E<sub>i,j</sub>, can be extracted from the A* matrix associated to the material chosen for this cell using the following expression:

$$E_{i,j} = A_{1,1i,j}^* \left( 1 - \frac{A_{1,2i,j}^* A_{2,2i,j}^*}{A_{1,1i,j}^* A_{2,2i,j}^*} \right) \qquad (32)$$

The center of inertia is shared for all cells in the same i-th transversal section. For each of these sets of elements, the vertical position of the center of inertia is computed as
<figure><img src="image_19.png" /><figcaption>(a) Repetitions to sustain pure shear loading</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>(b) Repetitions to sustain pure compression loading</figcaption></figure>

<caption>Fig. 14 Double-double homogenized properties study</caption>

$$y_{c_i} = \frac{\sum_{j=1}^{n_{py}} E_{i,j} y_{i,j} S_{i,j}}{\sum_{j=1}^{n_{py}} E_{i,j} S_{i,j}}, \qquad (33)$$

where $y_{i,j}$ is the y-coordinate of the center of the cell in the position $i, j$.

The stiffeners do not play any role in the position of the center of inertia, as they are assumed to be all equivalent in terms of mass and stiffness and evenly distributed on the panel. However, their effect is visible, in addition to that of the cells, in the computation of the average stiffness of the panel

$$E_{0_i} = \frac{\sum_{j=1}^{n_{py}} E_{i,j} S_{i,j} + n_s E_s S_s}{\sum_{j=1}^{n_{py}} S_{i,j} + n_s S_s}, \qquad (34)$$

where $n_s$ is the number of stiffeners, which for a given panel description is always $n_s = n_{py} + 1$, $E_s$ is their stiffness, and $S_s$ their transversal section, both are considered fixed parameters during the optimization.

Once the center of inertia is computed for all the sets of equal x-axis coordinates, the first and second moments of inertia can be computed for all cells. The cumulative first moment of inertia for any rectangular cell of the i-th section, taking into account the stiffeners, takes the following form:

$$Q_{i,j} = Q_{i,j-1} + S_{i,j}(y_{i,j} - y_{c_i}) \frac{E_{i,j}}{E_{0_i}} + S_s(y_{s,j} - y_{c_i}) \frac{E_s}{E_{0_i}}, \qquad (35)$$

where each time you add a panel contribution, the precedent stiffener contribution is also included.

The second moment of inertia is also expressed as the addition of the contributions of the different cells and stiffeners

$$I_i = \sum_{j=1}^{n_{py}} I_{i,j}^c + \sum_{k=1}^{n_s} S_s(y_k - y_{c_i})^2 \frac{E_s}{E_{0_i}}, \qquad (36)$$

where $I_{i,j}^c$ is the contribution of each cell of the i-th section to the moment of inertia. It is computed as

$$I_{i,j}^c = \left( \frac{\left( \frac{H}{n_{py}} \right)^3 t_{i,j}}{12} + S_{i,j}(y_{i,j} - y_{c_i})^2 \right) \frac{E_{i,j}}{E_{0_i}}. \qquad (37)$$

## Appendix C: QUAD stacking sequence catalog and compatibility

The following table represents the stacking sequences of the catalog used in these analyses. This catalog is a subset of a list of industrially explored configurations, selected through engineering knowledge. It is important to notice that some of the choices represent some compromises of the laminate MCs common to industrial applications. An example is the presence of pseudo-symmetric laminates. These are laminates that, while not being geometrically symmetric, present properties that are numerically close to those of symmetric laminates (Table 8).
<table><caption>**Table 8** Stacking sequence catalog</caption><thead><tr><th>Material ID</th><th>Stacking</th></tr></thead><tbody><tr><td>Material1</td><td>45/-45/-45/45/90/-45/0/45/90/45/-45/-45/45</td></tr><tr><td>Material2</td><td>45/-45/-45/45/90/-45/0/45/90/45/-45/-45/45</td></tr><tr><td>Material3</td><td>45/-45/-45/45/90/-45/0/90/0/45/90/45/-45/-45/45</td></tr><tr><td>Material4</td><td>45/-45/-45/45/90/-45/45/0/90/0/45/-45/90/45/-45/-45/45</td></tr><tr><td>Material5</td><td>45/-45/-45/45/90/-45/45/0/-45/90/45/0/45/-45/90/45/-45/-45/45</td></tr><tr><td>Material6</td><td>45/-45/-45/45/90/-45/45/0/-45/90/45/0/45/-45/90/45/-45/-45/45</td></tr><tr><td>Material7</td><td>45/-45/-45/45/90/-45/45/0/45/-45/90/45/-45/0/45/-45/90/45/-45/-45/45</td></tr><tr><td>Material8</td><td>45/-45/-45/45/90/-45/45/0/45/-45/-45/90/45/-45/45/0/45/-45/90/45/-45/-45/45</td></tr><tr><td>Material9</td><td>45/-45/-45/45/90/-45/45/0/45/-45/90/-45/90/45/90/-45/45/0/45/-45/90/45/-45/-45/45</td></tr><tr><td>Material10</td><td>45/-45/-45/45/90/-45/45/0/45/-45/90/-45/90/0/45/90/-45/45/0/45/-45/90/45/-45/-45/45</td></tr><tr><td>Material11</td><td>45/-45/-45/45/90/-45/45/0/45/-45/90/45/-45/90/-45/45/0/45/-45/90/45/-45/-45/45</td></tr><tr><td>Material12</td><td>45/-45/-45/45/90/-45/45/0/45/-45/90/45/-45/90/-45/45/0/45/-45/90/45/-45/-45/45</td></tr><tr><td>Material13</td><td>45/-45/-45/45/90/-45/45/0/45/-45/90/45/-45/90/-45/0/45/90/-45/45/0/45/-45/90/45/-45/-45/45</td></tr><tr><td>Material14</td><td>45/-45/-45/45/90/-45/45/0/45/-45/90/45/-45/90/-45/45/0/45/-45/90/-45/45/90/-45/45/0/45/-45/90/45/-45/-45/45</td></tr><tr><td>Material15</td><td>45/-45/-45/45/90/-45/45/0/45/-45/90/45/-45/90/-45/45/0/45/-45/90/-45/45/90/45/-45/45/0/45/-45/90/45/-45/-45/45</td></tr><tr><td>Material16</td><td>45/-45/-45/45/90/-45/45/0/45/-45/90/45/-45/90/-45/45/0/-45/45/0/45/-45/90/-45/45/90/45/-45/45/0/45/-45/90/45/-45/-45/45</td></tr><tr><td>Material17</td><td>45/-45/-45/45/90/-45/45/0/45/-45/90/45/-45/90/-45/45/0/-45/0/45/0/45/-45/90/-45/45/90/45/-45/45/0/-45/90/45/-45/45</td></tr><tr><td>Material18</td><td>45/-45/-45/45/90/-45/45/0/-45/90/0/45/90/45/-45/90/45/-45/-45/45</td></tr><tr><td>Material19</td><td>45/-45/-45/45/90/-45/45/0/45/-45/90/0/45/45/0/45/-45/90/45/-45/-45/45</td></tr><tr><td>Material20</td><td>45/-45/-45/45/90/-45/45/0/45/-45/90/45/0/-45/90/-45/45/0/45/-45/90/45/-45/-45/45</td></tr><tr><td>Material21</td><td>45/-45/-45/45/90/-45/45/0/45/-45/90/45/0/0/-45/90/-45/45/0/45/-45/90/45/-45/-45/45</td></tr><tr><td>Material22</td><td>45/-45/-45/45/90/-45/45/0/45/-45/90/45/-45/0/0/-45/45/90/-45/45/0/45/-45/90/45/-45/-45/45</td></tr><tr><td>Material23</td><td>45/-45/-45/45/90/-45/45/0/45/-45/90/45/-45/0/0/-45/45/90/45/-45/45/0/45/-45/90/45/-45/-45/45</td></tr><tr><td>Material24</td><td>45/-45/-45/45/90/-45/45/0/45/-45/90/45/-45/0/0/0/-45/45/90/45/-45/45/0/45/-45/90/45/-45/-45/45</td></tr><tr><td>Material25</td><td>45/-45/-45/45/90/-45/45/0/45/-45/90/45/-45/45/0/0/0/-45/-45/45/90/45/-45/45/0/45/-45/90/45/-45/-45/45</td></tr><tr><td>Material26</td><td>45/-45/-45/45/90/-45/45/0/45/-45/90/45/-45/-45/45/0/0/0/45/-45/45/90/45/-45/45/0/45/-45/90/45/-45/-45/45</td></tr><tr><td>Material27</td><td>45/-45/-45/45/90/-45/45/0/45/-45/90/45/-45/90/-45/45/0/45/-45/90/45/-45/45/0/45/-45/90/45/-45/-45/45</td></tr><tr><td>Material28</td><td>45/-45/-45/45/90/-45/0/0/0/45/90/45/-45/-45/45</td></tr><tr><td>Material29</td><td>45/-45/-45/45/90/-45/45/0/0/0/45/-45/90/45/-45/-45/45</td></tr><tr><td>Material30</td><td>45/-45/-45/45/90/-45/45/0/90/0/45/-45/90/45/-45/-45/45</td></tr><tr><td>Material31</td><td>45/-45/-45/45/90/-45/45/0/90/0/90/0/45/-45/90/45/-45/-45/45</td></tr><tr><td>Material32</td><td>45/-45/-45/45/90/-45/45/0/45/-45/90/0/90/0/-45/45/0/45/-45/90/45/-45/-45/45</td></tr><tr><td>Material33</td><td>45/-45/-45/45/90/-45/45/0/45/-45/90/45/0/-45/90/45/-45/45/0/45/-45/90/45/-45/-45/45</td></tr><tr><td>Material34</td><td>45/-45/-45/45/90/-45/45/0/45/-45/90/45/0/0/-45/90/45/-45/45/0/45/-45/90/45/-45/-45/45</td></tr><tr><td>Material35</td><td>45/-45/-45/45/90/-45/45/0/45/-45/90/45/0/0/0/-45/90/45/-45/45/0/45/-45/90/45/-45/-45/45</td></tr><tr><td>Material36</td><td>45/-45/-45/45/90/-45/45/0/45/-45/90/45/0/-45/90/45/0/-45/45/90/45/0/-45/90/45/-45/-45/45</td></tr><tr><td>Material37</td><td>45/-45/-45/45/90/-45/45/0/45/-45/90/45/0/-45/90/45/0/-45/45/90/45/0/-45/90/45/-45/-45/45</td></tr><tr><td>Material38</td><td>45/-45/-45/45/90/-45/45/0/45/-45/90/45/0/-45/90/45/0/-45/45/90/45/0/-45/90/45/-45/-45/45</td></tr><tr><td>Material39</td><td>45/-45/-45/45/90/-45/45/0/45/-45/90/45/0/-45/90/45/0/-45/45/90/45/0/-45/90/45/-45/-45/45</td></tr><tr><td>Material40</td><td>45/-45/-45/45/90/-45/45/0/45/-45/90/45/0/-45/90/45/0/-45/90/45/0/-45/45/90/45/0/-45/90/45/-45/-45/45</td></tr></tbody></table>

**Appendix D: Additional figures**

This section contains a series of figures providing additional insights into the results presented in Sect. 4.

Figure 15 shows our reference result for the comparison among the different formulations. It was obtained using LPF and forcing al cells to have the same suer-ply and thickness (LPF(1DD)).
<figure><img src="image_21.png" /><figcaption>**Fig. 15** Super-ply and cell thickness distribution for the DD panel using the LPF(1DD) formulation. The solution mass is 826.38 kg. This reference was obtained imposing the same super-ply and thickness in all cells</figcaption></figure>

Figures 16, 17, and 18 represent the solutions obtained when applying the optimal clustering strategy over the

LPF. The figures represent the solutions with three, two, and one cluster, respectively.

<figure><img src="image_22.png" /><figcaption>**Fig. 16** Super-ply and cell thickness distribution for the DD panel using the LPF+CL3 formulation. The solution mass is 407.72 kg</figcaption></figure>

<figure><img src="image_23.png" /></figure>
<figure><img src="image_24.png" /><figcaption>Fig. 17 Super-ply and cell thickness distribution for the DD panel using the LPF+CL2 formulation. The solution mass is 422.52 kg</figcaption></figure>

<figure><img src="image_25.png" /><figcaption>Fig. 18 Super-ply and cell thickness distribution for the DD panel using the LPF+CL1 formulation. The solution mass is 448.28 kg</figcaption></figure>
Similarly, Figs. 19, and 20 represent the solutions obtained when applying the optimal clustering strategy

over the MVF with two and one cluster, respectively. They complement Fig. 6 shown in Sect. 4.

<figure><img src="image_26.png" /><figcaption>Fig. 19 Super-ply and cell thickness distribution for the DD panel using the MVF+CL2 formulation. The solution mass is 415.76 kg</figcaption></figure>

<figure><img src="image_27.png" /><figcaption>Fig. 20 Super-ply and cell thickness distribution for the DD panel using the MVF+CL1 formulation. The solution mass is 448.76 kg</figcaption></figure>

<figure><img src="image_28.png" /></figure>
Figure 21 represents the solution obtained by MVF when the blending constraint (26) is not considered. It can be seen that some incompatible laminates were placed together, such as “material29” and “material2” or “material3”.

Figure 22 represents the solution obtained by the MVF over the VTP spar without considering the constraint in Eq. (28). It is possible to see by observing the thicknesses

evolution and the blending scheme in Fig. 23 that it presents some irregularities in its design, even if compatible, due to the non-monotonous behavior of the loading. This is the reason that led to the implementation of constraint (28).

Figure 24 shows the results for the VTP Spar considering the DD technique. As in the previous case, the effect of the different loading can be seen in the thickness

<figure><img src="image_29.png" /><figcaption>Fig. 21 Stacking and cell thickness distribution for the QUAD panel without blending constraint. The solution mass is 164.36 kg. For stacking sequence details, see Appendix C</figcaption></figure>

<figure><img src="image_30.png" /><figcaption>Fig. 22 Stacking and cell thickness distribution for the 11 subpanels Spar in QUAD. The solution mass is 30.06 kg. For stacking sequence details, see Appendix C</figcaption></figure>
<figure><img src="image_31.png" /><figcaption>**Fig. 24** Super-plies and cell thickness distribution for the 11 subpanels Spar with DD material. The solution mass is 26.52 kg</figcaption></figure>

<figure><img src="image_32.png" /></figure>

distribution, but due to the nature of DD, it is not as critical for manufacturing as in QUAD.

**Author Contributions** Alejandro Garcia Pis: conceptualization, formal analysis, methodology, investigation, visualization, and writing—original draft. Christian Gogu: conceptualization, investigation, validation, supervision, funding acquisition, and writing—review & editing. Stephane Grihon: conceptualization, investigation, validation, supervision, funding acquisition, and writing—review & editing. Simone Coniglio: conceptualization, investigation, validation, supervision, and writing—review & editing.

**Funding** Open access funding provided by ISAE-SUPAERO.

**Data Availability** To reproduce the data and results of the problem presented in Sect. 2 with the redistribution scheme shown in Appendix B and the material properties computed as expressed in Appendix A, the open source GEMSEO (Gallard et al. 2018) library was used in conjunction with its Pymoo and MMA plugins https://gitlab.com/gemseo/dev/gemseo-mma. The BLOA algorithm with the modifications described in Sect. 3 is available in the GEMSEO library at

https://gitlab.com/gemseo/dev/gemseo-bilevel-outer-approximation. The hyperparameters used in the BLOA formulation are detailed for each case in Sect. 4. The NGSA tool used in the wing spar use-case is an industrial closed source code which cannot be shared.

## Declarations

**Conflict of interest** The authors declare that they have no conflict of interest.

**Open Access** This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

<figure><img src="image_33.png" /></figure>
## References

Allaire G, Delgado G (2016) Stacking sequence and shape optimization of laminated composite plates via a level-set method. J Mech Phys Solids 97:168-196

Barjhoux PJ (2020) Towards efficient solutions for large scale structural optimization problems with categorical and continuous mixed design variables. PhD thesis. http://www.theses.fr/2020E SAE0004.

Bendsøe M (1989) Optimal shape design as a material distribution problem. Struct Optim 1:193-202. https://doi.org/10.1007/BF016 50949

Blank J, Deb K (2020) Pymoo: multi-objective optimization in python. IEEE Access 8:89497-89509. https://doi.org/10.1109/ACCESS. 2020.2990567

Bruyneel M (2008) Optimization of laminated composite structures: problems, solution procedures and applications. Composite materials research progress. NOVA Science Publishers, Hauppauge, NY, pp 51-107

Catapano A, Montemurro M (2014) Homogenisation of core properties. A multi-scale approach for the optimum design of sandwich plates with honeycomb core: part I. Compos Struct 118:664-676

Coniglio S (2023) Non-convex, mixed categorical-continuous variable structural optimization applied to metallic stiffened panel sizing, 15th World Congress on Structural and Multidisciplinary Optimization, 5-9 June 2023, Cork, Ireland

Enrique Herencia J, Weaver PM, Friswell MI (2007) Optimization of long anisotropic laminated fiber composite panels with t-shaped stiffeners. AIAA J 45(10):2497-2509. https://doi.org/10.2514/1. 26321

Gallard F, Vanaret C, Guénot D, Gachelin V, Lafage R, Pauwels B, Barjhoux PJ, Gazaix A (2018) Gems: a python library for automation of multidisciplinary design optimization process generation. In: AIAA/ASCE/AHS/ASC Structures. Structural Dynamics and Materials Conference, p 0657

Garofano A, Sellitto A, Di Caprio F, Riccio A (2023) On the use of double-double design philosophy in the redesign of composite fuselage barrel frame components. Polym Compos 45(5):4250-4265. https://doi.org/10.1002/pc.28056

Golberg DE (1989) Genetic algorithms in search, optimization, and machine learning. Addison Wesley, Reading

Irisarri FX, Lasseigne A, Leroy FH, Le Riche R (2014) Optimal design of laminated composite structures with ply drops using stacking sequence tables. Compos Struct 107:559-569. https://doi.org/10. 1016/j.compstruct.2013.08.030

Jing Z, Chen J, Sun Q (2019) Constrained-manufacturable stacking sequence design optimization using an improved global shared-layer blending method and its 98-line Matlab code. Struct Multidiscip Optim 59(2):539-575. https://doi.org/10.1007/s00158-018- 2083-3. (ISSN 1615-1488.)

Kappel E (2022) Double-double laminates for aerospace applications-finding best laminates for given load sets. Compos C Open Access 8:100244. https://doi.org/10.1016/j.jcomc.2022.100244

Lund E, Stegmann J (2005) On structural optimization of composite shell structures using a discrete constitutive parametrization. Wind Energy 8(1):109-124. https://doi.org/10.1002/we.132

Miki M (1985) Design of laminated fibrous composite plates with required flexural stiffness. ASTM International, West Conshohocken, PA

Montemurro M, Catapano A (2017) On the effective integration of manufacturability constraints within the multi-scale methodology

for designing variable angle-tow laminates. Compos Struct 161:145-159

Montemurro M, Fiordilino GA, Carrera E (2023) Multi-level optimization of composite structures through a global-local modelling approach based on high-order theories. Comput Struct 275:106932. https://doi.org/10.1016/j.compstruc.2022.106932

Nikbakt S, Kamarian S, Shakeri M (2018) A review on optimization of composite structures part I: laminated composites. Compos Struct 195:158-185. https://doi.org/10.1016/j.compstruct.2018.03.063

Picchi Scardaoni M, Montemurro M (2022) The design of blended laminates regardless of the stack: the search propagation direction. Compos Commun 36:101385. https://doi.org/10.1016/j. coco.2022.101385

Picchi Scardaoni M, Montemurro M, Panettieri E, Catapano A (2021) New blending constraints and a stack-recovery strategy for the multi-scale design of composite laminates. Struct Multidiscip Optim 63(2):741-766. https://doi.org/10.1007/ s00158-020-02725-x

Sørensen S, Lund E (2013) Topology and thickness optimization of laminated composites including manufacturing constraints. Struct Multidiscip Optim 48:08. https://doi.org/10.1007/ s00158-013-0904-y

Sørensen S, Sørensen R, Lund E (2014) DMTO: a method for discrete material and thickness optimization of laminated composite structures. Struct Multidiscip Optim 50:25-47

Svanberg K (1987) The method of moving asymptotes-a new method for structural optimization. Int J Numer Methods Eng 24(2):359-373. https://doi.org/10.1002/nme.1620240207

Tsai SW (2021) Double-double: new family of composite laminates. AIAA J 59(11):4293-4305. https://doi.org/10.2514/1.J060659

Tsai SW, Melo JDD, Sihn S, Arteiro A, Rainsberger R (2017) Composite laminates: theory and practice of analysis, design and automated layup. Stanford Aeronautics & Astronautics, Stanford, CA

Venkataraman S, Haftka R (1999) Optimization of composite panels—a review. In: 14th Technical Conference on American Society for Composites, Fairborn, OH, pp 479-488

Vermes B, Tsai SW, Massard T, Springer GS, Czigany T (2021) Design of laminates by a novel double-double layup. Thin Walled Struct 165:107954

Wang Y, Wang D, Zhong Y, Rosen DW, Li S, Tsai SW (2023) Topology optimization of double-double (DD) composite laminates considering stress control. Comput Methods Appl Mech Eng 414:116191. https://doi.org/10.1016/j.cma.2023.116191

Wang Z, Zheng C, Wang L, Gu J, Jing L, Lai X (2024) Optimal design of composite grid/skin structures based on deep learning and double-double layup strategy. Aerosp Sci Technol 147:109030

Zhang Z, Zhang Z, Di Caprio F, Gu GX (2022) Machine learning for accelerating the design process of double-double composite structures. Compos Struct 285:115233. https://doi.org/10.1016/j. compstruct.2022.115233

Zhao K, Kennedy D, Miravete A, Tsai SW, Featherston CA, Liu X (2023) Defining the design space for double-double laminates by considering homogenization criterion. AIAA J 61(7):3190-3203. https://doi.org/10.2514/1.J062639

**Publisher's Note** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.