Lightweight Design of Variable-Stiffness Cylinders with Reduced Imperfection Sensitivity Enabled by Continuous Tow Shearing and Machine Learning

The present study investigates how to apply continuous tow shearing (CTS) in a manufacturable design parameterization to obtain reduced imperfection sensitivity in lightweight, cylindrical shell designs. The asymptotic nonlinear method developed by Koiter is applied to predict the post-buckled stiffness, whose index is constrained to be positive in the optimal design, together with a minimum design load. The performance of three machine learning methods, namely, Support Vector Machine, Kriging, and Random Forest, are compared as drivers to the optimization towards lightweight designs. The new methodology consists of contributions in the areas of problem modeling, the selection of machine learning strategies, and an optimization formulation that results in optimal designs around the compromise frontier between mass and stiffness. The proposed ML-based framework proved to be able to solve the inverse problem for which a target design load is given as input, returning as output lightweight designs with reduced imperfection sensitivity. The results obtained are compatible with the existing literature where hoop-oriented reinforcements were added to obtain reduced imperfection sensitivity in composite cylinders.


Introduction
The buckling performance of cylindrical shells is recognized to be highly sensitive to geometric, load, and material imperfections, a behavior that was identified already since the first attempts to correlate experimental with theoretical predictions. The use of straight-fiber laminated composite materials has enabled the maximization of the mass-normalized loadcarrying capacity of these shells, while customarily applying conservative knock-down factors to reduce the theoretical load-carrying capacity. Recent studies present measured imperfections for filament-wound composite cylinders [1], showing a thickness pattern that highly depends on the design and that significantly affects the buckling strength. Furthermore, other defects, such as those created during service due to low-velocity impacts [2][3][4], can also lead to a significant knock-down on the theoretical performance of cylindrical shells.
The application of novel composite manufacturing technologies such as continuous tow shearing (CTS) can enable a larger design space by means of fiber steering and thickness variation. This larger design space encompasses composite cylindrical shell designs that showed an imperfection-insensitive behavior and which have the potential to undergo buckling without collapse.
Imperfection-sensitive structures such as cylindrical shells are known for their large number of nonlinear paths and nearby bifurcation points. During the experimental evaluation of such structures, only one nonlinear path is preferably chosen according to the localization event produced by geometric, load, and material imperfections. The high sensitivity of the chosen nonlinear path with respect to this localization event was called, by Thomson and Virgin, "spatial chaos" [5], justifying the large scatter of experimental data that is usually obtained from the experimental assessment of cylindrical shells, and other imperfection-sensitive shells.
Recently, Groh et al. [6] proposed a representative and simple model based on rigid links supported by transverse springs to numerically demonstrate the concept of spatial chaos, where small changes in the geometric imperfections result in a large scatter of the maximum load-carrying capacity. The NASA-8007 guideline [7], which is derived from the collection of experimental results from Seide, Weingarten, and Morgan [8,9], reports a large number of experimental results for cylindrical shells having different radius-per-thickness ratios, where the aforementioned scatter can be clearly observed. The guideline defines a lower-bound curve with respect to the scattered data that is frequently used to calculate knock-down factors (KDF), and these KDF can be directly used in combination with the theoretical linear buckling load of cylindrical shells, as discussed by Hilburger [10]. Such a resource is useful while designing shells for a given required compression load and while taking the imperfection sensitivity into account.
White et al. [34] proposed the use of Koiter's asymptotic method to achieve reduced imperfection sensitivity in variable-stiffness cylindrical panels, solving the numerical problem by means of generalized differential and integral quadrature. A genetic algorithm optimizer was used to find optimal straight-and curved-fiber composite layouts for the curved panels under investigation. The best post-buckling performance achieved by the curved-fiber designs showed virtually no drop in axial stiffness while entering the postbuckling regime. The work of White et al. [34] inspired the present study on the lightweight design of cylindrical shells, aiming at reduced imperfection sensitivity.
The main scientific contributions of the present study are: • a new parameterization scheme for variable-stiffness cylindrical shells manufactured using continuous tow shearing (CTS), which result in a vast design space that is explored in order to attain designs with reduced imperfection sensitivity; • the comparison of the performance of three classes of machine learning strategies with multiple kernel models, in a total of eighteen distinct instances; • an algorithm to calculate the optimal design that meets the target levels of mass, buckling load, and post-buckling stiffness. The inverse problem is formulated within the scope of global optimization for mixed-integer variables.
The study is concluded by showing the potential of the algorithm to solve the inverse problem related to the proposed design scheme, where for each design load, the optimal design that is lightweight and has a reduced imperfection sensitivity can be obtained efficiently.
The paper is organized as follows. Section 2 presents the characteristics of the computational model of a cylinder and explains the key features considered in the analysis. The linear buckling is discussed in Section 3, and the asymptotic expansion is discussed in Section 4. Section 5 outlines the characteristics of the Support Vector Machine, Kriging, and Random Forest algorithms, and presents the results obtained from the adjustment of the kernel parameters. Section 6 details the proposed algorithm that computes the inverse problem and discusses the results obtained from numerical experiments. The conclusions are summarized in Section 7.

Design Parameterization with Circumferentially Oriented Thickness Patterns
The proposed parameterization consists of circumferentially oriented thickness patterns created from the nonlinear steer-thickness coupling that is enabled by the continuous tow shearing (CTS) manufacturing process. Figure 1 schematically shows the design parameterization in an example with two regions c 2 (n = 2) and three regions c 1 . Between the regions c 1 and c 2 , there is a transition region t. The tow angle θ(x) at the regions c 1 and c 2 is, respectively, θ 1 and θ 2 . In the continuous tow shearing (CTS) process, the robot head moves in such a way that the deposited tow width remains constant, such that any tow angle different than zero will generate a tow of higher thickness, where the local thickness can be calculated with as shown in [35][36][37][38]. This is true if we assume that the shift direction of the deposition head, where it moves to start the next tow, is the circumferential direction of the cylinder. The transition regions t that are located between the two regions c 1 and c 2 are made by steering the fibers with a defined steering radius r CTS ≥ r CTSmin , being r CTSmin = 50 mm [35], as inspired by the work of Ummels and Castro [39] on overlap-stiffened structures suggesting pre-determined steering angle values. In the present study, it is assumed that the fiber angles in the transition region change linearly with x. Hence, the length t can be calculated as If n regions with length c 2 are used between n + 1 regions with length c 1 , the following simple equation for the total cylinder length L can be derived: A nondimensional quantity c 2ratio can be conveniently used to express the ratio between regions c 1 and c 2 : such that 0 ≤ c 2ratio ≤ 1. The value c 2max is the maximum length of region c 2 , and is determined assuming c 1 = 0 in Equation (3), leading to Assuming that the whole cylinder length is covered with transition regions t, implying c 1 = c 2 = 0, the maximum number of c 2 regions, or the maximum value for n, can be calculated for the proposed design parameterization as The special case c 2ratio = 0 and n > 0 leads to a feasible design consisting of a combination of regions c 1 and transition regions t. The case with c 2ratio = 1 and n > 0 also leads to a feasible design combining transition regions t with regions c 2 . When c 2ratio = 0 and n = 0, a constant-stiffness cylinder with stacking sequence ±θ 1 is generated, whereas c 2ratio = 1 and n = 0 leads to a constant-stiffness cylinder with stacking sequence ±θ 2 .
In summary, the proposed parameterization takes, as input: • n: the number of regions c 2 ; • r CTS : the CTS steering radius in the transition zone; • c 2ratio : the ratio between the design length c 2 with respect to c 2max , which is calculated as per Equation (5); • θ 1 and θ 2 : the respective angles at regions c 1 and c 2 .
With these inputs, the following design parameters are calculated: • The length of the transition regions t, according to Equation (2); • The length c 2 , based on c 2ratio and c 2max ; • The length c 1 , according to Equation (3).
Furthermore, the adopted strategy to control the finite element mesh size takes, as input: • n x t : the number of nodes along x, along the length of each transition zone, given by t; • n y : the number of nodes along the circumferential direction, which is constant along the length. If n y is not given, then it is calculated assuming a desired aspect ratio of n y n x = 1; • the maximum aspect ratio of the elements in the circumferential direction, with respect to the axial direction: n y n x max .
If this value is exceeded, more nodes in the x direction are added to the regions c 1 and c 2 to keep n y n x ≤ n y n x max .
Note that Equation (7) does not affect the transition regions with length t, which have their number of nodes along the length determined by n x t .
With these mesh inputs, the number of nodes along the lengths c 1 and c 2 are then calculated. The workflow is shown in Algorithm 1.
The proposed design parameterization creates a tailored circumferentially oriented thickness pattern, and enables the proposition of an optimization formulation. Setting the main objective as a lightweight design aiming at minimum mass m, with a given design load P design , a positive post-buckling stiffness b I > 0, and the linear buckling load of the cylinder P critical , the resulting optimization problem can be formally stated as min (r CTS ,n,c 2ratio ,θ 1 ,θ 2 ) m (9) subject to Algorithm 1: The modeling phase /* Initial setup */ n ← number of regions c 2 ; r CTS ← steering radius; c 2ratio ← ratio in agreement with Equation (5); θ 1 ← tow angle at c 1 ; θ 2 ← tow angle at c 2 ; /* Intermediary computation */ t ← length of the transition regions (Equation (2)); c 2 ← length (Equations (4) and (5)); c 1 ← length (Equation (3) Although the formulation expressed by Equations (9) and (10) may lead to optimal designs, there is no guarantee that a feasible design exists for an arbitrary P design and b I > 0. Therefore, a new multi-objective formulation is proposed: min (r CTS ,n,c 2ratio ,θ 1 ,θ 2 ) where m re f , P re f , and b re f are reference values for m, P critical , and b I , respectively. The weighting values α 1 , α 2 , and α 3 represent the preference in optimizing each criteria. An advantage of this formulation (Equations (11) and (12)) is the ability to contemplate some constraints as objectives to be minimized. It avoids the issue that the set of inequalities (Equation (10)) may, in some circumstances, lead to an empty set of solutions. Furthermore, the optimization process will provide feasible designs even when an infeasible set of reference values m re f , P re f and b re f is given.
Although finding suitable reference values for the new objective function (Equation (11)) is not trivial, the current proposition consists in the usage of machine learning strategies as a heuristic to provide useful estimates. Those parameters will then be used as part of a global optimization procedure, in the domain of mixed-integer variables, which will refine the search of the global optimum.
Further details on this topic can be found in Section 6.

Linear Buckling Constraint
The linear buckling load P critical is compared to the design load P design in the optimization scheme described by Equation (10). In order to calculate the linear buckling load of the proposed variable-stiffness design, the single-curvature Bogner-Fox-Schmit-Castro (SC-BFSC) finite element is selected. First presented by Wang et al. [40], the SC-BFSC element is based on the plate version presented by Castro and Jansen [24], which is a C1-contiguous and confirming element obtained by taking tensor products of cubic Hermite splines [41]. With 4 nodes per element and 10 degrees of freedom per node, the BFSC approximates the in-plane and out-of-plane displacements using third-order polynomials, leading to a fast convergence for cases with variable stiffness. Figure 2 illustrates a SC-BFSC element and the global coordinate system xyz used, where coordinate y is curvilinear following the circumferential perimeter, such that at y = 2πr, the path closes on itself. The displacements along x, y, z are, respectively, u, v, w, and they are approximated within each element as: where S S S u,v,w i and u e u e u ei are the shape functions and the 10 degrees of freedom of the ith node of the SC-BFSC element, being in the following order: u, ∂u/∂x, ∂u/∂y, v, ∂v/∂x, ∂v/∂y, w, ∂w/∂x, ∂w/∂y, ∂ 2 w/∂x∂y . For the SC-BFSC element, the same shape functions of the plate BFSC element [24] can be used: with the cubic Hermite functions H i , H x i , H y i , H xy i calculated using natural coordinates [42][43][44]: where x , y are, respectively, the finite element dimensions along x, y, as represented in Figure 2. Using the proposed nodal connectivity for the SC-BFSC element, the nodal degrees of freedom u e u e u ei and the respective shape functions S S S u,v,w i are concatenated as: u e u e u e = u e u e u e1 u e u e u e2 u e u e u e3 u e u e u e4 with S S S u , S S S v , and S S S w being matrices of shape 1 × 40.
The total potential energy of the cylindrical shell under axial compression can be represented by Equation (17), where Ω represents the shell domain and δΩ the shell boundaries, i.e., the loaded edges. The first term represents the strain energy based on equivalent single-layer theories [45,46], while the second term represents the work conducted by external forces at the boundaries N N N, multiplied by a scalar λ, with N N N expressed in force per length units.
The strain energy functional of the entire shell continuum is built from the individual contributions of all finite elements φ e : where the membrane forces are N N N = N xx , N yy , N xy and the distributed moments are The integration limits x 1 ≤ x ≤ x 2 and y 1 ≤ y ≤ y 4 define the domain of one finite element Ω e . The membrane ε ε ε and rotational κ κ κ strains are assumed to follow Sanders-type kinematics [47], also discussed in references [48,49] for cylindrical shells, and in references [50][51][52][53] for conical shells: with (·), x = ∂(·)/∂x used as a compact notation for partial derivatives. At the bifurcation point, the following state of equilibrium exists, considering all n e elements: Expressing the displacements within one element in terms of nodal coordinates u e u e u e , as in Equation (13), δε ε ε and δκ κ κ can be calculated from Equation (19) as: where B m B m B m and B b B b B b are defined as: noting that the partial derivatives of S S S u,v,w are directly calculated from the shape functions of Equation (14) in terms of the natural coordinates ξ, η, using the Jacobian relations ∂/∂x = x /2∂/∂ξ and ∂/∂y = y /2∂/∂η. Applying the neutral equilibrium criterion δ 2 φ e = 0 [53] leads to: The first integral of Equation (22) becomes the constitutive stiffness matrix of the element, calculated using the constitutive relations from classical laminated plate theory [46], Note that the geometric non-linearity appears in the constitutive stiffness matrix due to w ,x , w ,y , and v in Equation (21). Therefore, the linear constitutive stiffness matrix of a finite element K e K e K e is calculated assuming The second integral of Equation (22) becomes the geometric stiffness matrix of the finite element K G0e K G0e K G0e , capturing the geometrically nonlinear effects of a linear pre-buckling membrane stresses N 0 N 0 N 0 = N 0xx , N 0yy , N 0xy on the membrane stiffness. Noting that δ 2 κ κ κ = 0 0 0 [24,53], the equation for K G0e K G0e K G0e becomes: The contributions of all n e finite elements are added to build the global constitutive stiffness matrix K K K and the geometric stiffness matrix K G0 K G0 K G0 of the system: The linear pre-buckling stress field of one finite element N 0xx , N 0yy , N 0xy is calculated from the corresponding nodal displacements u 0e u 0e u 0e as: where u 0e u 0e u 0e is directly extracted from the full pre-buckling displacement vector u 0 u 0 u 0 that can be obtained from a linear static analysis, derived from the equilibrium of Equation (20): with f 0 f 0 f 0 representing any general pre-buckling force, here calculated based on a uniformly distributed axially compressive unit force P unit = 1. Assuming that there is a stress state N N N = λN 0 N 0 N 0 that leads to the condition δ 2 φ = 0, the problem consists of finding the value of λ, such that: which holds true for any variation δu u u, such that the required condition for the equality of Equation (28) is: Equation (29) represents a linear buckling eigenvalue problem, solved numerically in the present study by means of the locally optimal block preconditioned conjugate gradient method [54], implemented in SciPy [55]. The lowest eigenvalue is the critical linear buckling load multiplier λ c , from which the critical compressive load can be calculated with P critical = λ c .
The integration of K e K e K e and K G0 e K G0 e K G0 e over the finite element domains are performed numerically using standard Gauss-quadrature with 4 × 4 integration points per element. The authors verified that this amount of integration points leads to a converged behavior even for variable-stiffness filament-wound cylinders [40]. For each integration point, the local shell constitutive properties are calculated using the smeared approach proposed by Castro et al. [36] based on a constant-volume constraint, verified against a discreet thickness modeling by Vertonghen and Castro [37].

Post-Buckling Stiffness
The post-buckling stiffness b I , also known as b f actor , is herein calculated using the proposed SC-BFSC element, based on Koiter's asymptotic approach. Given a total potential energy functional φ[u u u, λ] that depends on displacements u u u and a scalar load multiplier λ, a pre-buckling static equilibrium with solution u 0 u 0 u 0 that corresponds to a load level λ 0 can be described as per Equation (30), where the notation φ δu u u is used instead of δφ to conveniently express the functional variation as a tensor product between the Fréchet derivative φ and the variation of the vector containing all degrees of freedom δu u u [24,25,56]. Assume that there exists at least one point of equilibrium that intersects [u u u(λ), λ] at a bifurcation point [u c u c u c , λ c ], such that u c u c u c = λ c u 0 u 0 u 0 , with λ c representing the critical buckling load, or critical bifurcation load, such that: Koiter [11] proposed expressing u u u − u c u c u c and λ − λ c using asymptotic expansions to represent the difference between the current displacements and the displacements at the bifurcation point with the corresponding load increment λ − λ c : where: 1. ξ is a scalar parameter; 2.
u I u I u I is a first-order field, taken directly from one or a linear combination of multiple linear buckling modes. Vector u I u I u I is customarily re-scaled by dividing with the maximum normal displacement amplitude and multiplying by the plate or shell thickness; 3.
u I I u I I u I I is a second-order field that provides a correction to the first-order field; 4. the third-order field u I I I u I I I u I I I , and higher fields, are assumed to have a negligible contribution; 5.
a I and b I are, respectively, the first-and second-order load parameters to be determined.
Equation (31) is a reduced-order model (ROM) relating the load λ and displacement u u u around the equilibrium point [u c u c u c , λ c ]. Note that Equation (31) consists of a reduced-order model to calculate displacements u u u based on a pre-buckled state u c u c u c with known first-and second-order fields u I u I u I and u I I u I I u I I . The coefficient ξ can be determined for each known load level λ, after the calculation of the coefficients a I and b I .
The expression given by Equation (31) can be applied to the expanded total potential energy functional of Equation (32) [24,25], leading to: The resulting formula is shown in Equation (33), where the terms multiplying ξ 2 and ξ 3 are collected. Note that this expansion concerns only a single mode, and the reader is referred to Castro and Jansen [24,25] for details about the multi-modal expansion adopting the same notation. It follows that where φ c is the tangent stiffness matrix; the calculation of and φ iv c is discussed in detail by Castro and Jansen [24,25]. For the expanded equilibrium to be stationary, each term in Equation (33) must vanish separately. Assuming δu u u = u I u I u I in Equation (33), the expressions for a I and b I of Equation (31) can be obtained, as given, respectively, in Equations (34) and (35): The second-order fields u I I u I I u I I needed for the calculation of b I can be obtained by first calculating a non-orthogonal second-order fieldū I Ī u I Ī u I I . The calculation ofū I Ī u I Ī u I I can be conducted directly from Equation (33), knowing that the term multiplying ξ 2 must vanish separately: The orthogonal second-order field vector u I I u I I u I I in the single-modal asymptotic expansion can be obtained after a Gram-Schmidt orthogonalization [57] operation, used to remove the components ofū I Ī u I Ī u I I that are parallel to the linear buckling mode used in the single-modal expansion named u I u I u I , as formulated in Equation (37).
u I I u I I u I I =ū I Ī u I Ī u I I − u I u I u I ū I Ī u I Ī u I I , u I u I u I u I u I u I , u I u I u I This topic concludes the discussion of the modeling strategy, with the exposition that follows focusing on the application and performance evaluation of three machine learning (ML) strategies applied on the optimization problem of Equation (11).

Machine Learning Strategies for Meta-Modeling
An interesting characteristic of machine learning algorithms is their ability to generalize the training experience and provide an unexpected output that best fulfills the objective function, predicting future events or scenarios that are not explicitly mapped in the training process, and that can, therefore, reach results that are unexpected and non-intuitive. Among a myriad of alternatives for implementing a learning strategy, a comprehensive overview of popular strategies and resources is presented by Russell [58].
In the present work, the viability of the Support Vector Machine (Section 5.2), the Kriging surrogate (Section 5.3), and the Random Forest (Section 5.4) algorithms are investigated as learning strategies to meta-model and optimize cylinder design. Numerical calculations were performed by several packages of the R software platform [59]. The discussion starts with the design of the experiments (Section 5.1).

Design of Experiment
The steering radius r CTS , the number n of c 2 regions (Equation (6)), the ratio c 2ratio (Equation (4)), the tow angle θ 1 at region 1, and the tow angle θ 2 at region 2 constitute the five design variables v = [r CTS , n, c 2ratio , θ 1 , where The real-valued variables v 4 and v 5 are limited to two decimal places, in order to be closer to a viable solution for manufacturing. The response in terms of buckling load P critical (v) and post-buckling stiffness b I (v) are evaluated throughout the design space by means of a Design of Experiment (DOE) where a set of 2000 feasible points is chosen from the design envelope by means of a Latin hypercube sampling methodology, as available in the LHS toolbox [60]. The DOE is publicly available in a dataset published by the authors [61]. The finite element model is evaluated at each point of the DOE, and the results are presented below. Figure 3 shows the mass versus b I obtained from every evaluation within the DOE, where no clear correlation between the variables can be seen. Figure 4 plots mass versus P critical , showing a Pareto frontier at the upper limit of P critical for increasing values of mass.  The effectiveness of meta-modeling the b I , mass, and P critical indices by means of SVM, Kriging, and Random Forest methodologies are discussed in Sections 5.2-5.4, respectively. The numerical results, obtained from a set of n sample = 2000 feasible points, are divided into two groups: the first group of n training = 1600 points is used in the training phase, to adjust the models, and the second group of n test = 400 points is used in the test phase, to infer the quality of the models.
The relative error E(v) of the meta-model is used to quantify the deviance between the actual f (v) and estimatedf (v) values of a design v, and is evaluated as while providing estimates for the post-buckling stiffness b I , the mass, and the buckling load P critical . There are several boxplot diagrams presented in Figures 5-16, whose data are calculated following Equation (40). The results summarized in Tables 1-9 are also obtained by using Equation (40). As the sample points are randomly chosen in the design space, from the project perspective, it is important to identify the points along the frontier that offer a positive post-buckling stiffness, herein represented by the criteria b I > 0. This discussion will be expanded in Section 6.

Support Vector Machine
The standard SVM is formulated as a classifier whose decision function is represented by a hyper-plane that maximizes the distance of separated samples from different classes. The methodology ultimately finds a linear combination of features that characterizes or separates two or more classes of objects or events (as described by y i ). Since the convex optimization problem is solved by a deterministic algorithm, for a specific input dataset, the same optimal hyper-plane parameter is obtained, as a solution of Equations (41)- (43).
Formally, the minimization problem is given by subject to where and m ≤ n training is the number of active constraints. A detailed discussion about the SVM method and variants can be found in Awad and Khanna [62] and Salcedo-Sanz and Rojo-Álvarez [63].
In the present study, the SVM methodology is evaluated considering the characteristics of seven different kernels (as expressed by φ in Equation (42)): Laplace, RBF, Polynomial, Vanilla, Hyperbolic Tangent, Bessel, and ANOVA. The executions are performed thought the Kernlab toolbox [64]. A set of 2000 feasible points is chosen from the design envelope using a Latin hypercube sampling methodology that is available in the LHS toolbox [60]. Then, the set is divided into two groups: the first group of 1600 points is used in the training phase, to adjust the models, and the second group of 400 points is used in the test phase, to infer the quality of the models. As a result, the performance of each kernel will be evaluated against a set of 400 designs not used in the training phase.
The minimum error, the median error, and the maximum error in the estimation of the b I are shown in Table 1. The error dispersions (Equation (40)) of all SVM methods are compared in Figure 5.
The median of the error dispersion is used as a metric to identify the best adjustment. According to this criteria, the first model, using the Laplace kernel, is found to be the best b I predictor.  For the mass estimation, the minimum error, the median error, and the maximum error are shown in Table 2. The error dispersion of the mass estimates using all SVM methods are compared in Figure 6. Based on the median error, the Laplace kernel is the best b I predictor.
For the estimation of P critical , the minimum error, the median error, and the maximum error are shown in Table 3. The error dispersion of the P critical estimate is shown in Figure 7. The Laplace kernel is the best P critical predictor, based on the median error.
According to the metric of the median error, the SVM combined with the Laplace kernel is the best predictor for the design space created under the proposed parameterization of variable-stiffness cylinders.

Kriging Surrogate
Haeri and Fadaee [65] present a reliability analysis of laminated composites using a Kriging surrogate model, where the performance is compared with results obtained from a radial basis neural network. The comparison results confirm the computational cost efficiency and accuracy of Kriging model for reliability analyses of laminated composites. The Kriging method [66,67] treats the function of interest as a realization of a stochastic process y(v) by means of the equation: where f (v) is a linear regression of the data, with k regressors modeling the drift of the process mean over the domain. The second part, Z(v), is a model of a Gaussian and stationary random process with zero mean and covariance between two observations, v 1 and v 2 . The authors of [68] present an overview of common spatial correlation functions R(v 1 , v 2 ) for approximating a deterministic computer model and the influence of parameter selection on the properties of these functions.
In the present study, the Kriging methodology is evaluated considering the characteristics of five different spatial correlation kernels: Gaussian, Matern 5 2 , Matern 3 2 , Exponential, and Power-Exponential. The performance of each model is evaluated against a set of 400 testing designs not used in the training phase. The executions are performed through the DiceKriging toolbox [69].
A set of 2000 feasible points is chosen from the design envelope using a Latin hypercube sampling methodology, available in the LHS toolbox [60]. Then, the set is divided into two groups: the first group of 1600 points is used in the training phase, to adjust the models, and the second group of 400 points is used in the test phase, to infer the quality of the models.
For the b I estimation using Kriging with different kernels, the minimum error, the median error, and the maximum error as evaluated by Equation (40) are shown in Table 4. The error dispersion of all Kriging methods are compared in Figure 8. The median of the error dispersion is used as a metric to identify the best adjustment. According to this criteria, the first model, using the Gaussian kernel, is found to be the best b I predictor.  The minimum error, the median error, and the maximum error in the estimation of the mass m are shown in Table 5. The error dispersion of the mass estimates using all Kriging methods are compared in Figure 9, where the first model, based on the Gauss kernel, is the best mass predictor according to the metric of the median error.  For the estimation of P critical , the errors are shown in Table 6. The error dispersion of the P critical estimate is shown in Figure 10. The Matern 3 2 kernel is the best P critical predictor, based on the median error; the other kernels lead to similar results, and therefore, they can be chosen without compromising the analysis.
It noteworthy that there is no single option with superior performance on all three datasets (mass, b I , and P critical ). In this case, the performance is dependent on the data under consideration.  Figure 10. Error dispersion of the P critical estimate using distinct Kriging kernel options.

Random Forest
The Random Forest [70] algorithm is a supervised learning procedure which operates according to the divide-and-conquer principle: sample fractions of the data, grow a randomized tree predictor on each small piece, and then aggregate these predictors together. It can be applied to a wide range of prediction problems and have few parameters to tune.
A Random Forest consists in a predictor composed by a collection of M randomized regression trees. For the j-th tree in the family, the predicted value at the query point x is denoted by m n (x; Θ j , D n ), where Θ 1 , . . . , Θ M are independent random variables.
A key result of such methodology is that as the number of trees increases, all sequences of generalization error converge almost surely. This result explains why Random Forests do not overfit as more trees are added, but produce a limiting value of the generalization error.
In the present study, the Random Forest is evaluated considering the number of trees to grow (tree ∈ {50, 100, 200}) and the number of variables randomly sampled as candidates at each split (split ∈ {3, 5}). The performance of each model is evaluated against a set of n test = 400 designs not used in the training phase. The executions are performed thought the randomForest toolbox in R [71].
The minimum error, the median error, and the maximum error in the estimation of the b I , as evaluated by Equation (40), are shown in Table 7. The error dispersion of all Random Forest methods are compared in Figure 11. The median of the error dispersion is used as a metric to identify the best adjustment. According to this criteria, the fifth model, using the 100 trees and 5 splits, is found to be the best b I predictor. Random Forest Error Figure 11. Error dispersion of the b I estimate using distinct Random Forest parameters.
The errors in the estimation of the mass are shown in Table 8, while the dispersion using all Random Forest methods are compared in Figure 12. The parameters tree = 200 and split = 5 provided the best mass predictor, based on the median error.
For the estimation of P critical , the error data is given in Table 9, with the error dispersion shown in Figure 13 where, based on the median error, the parameters tree = 200 and split = 3 provided the best P critical predictor.   Different values for trees and splits resulted in small performance deviations. A similar profile is observed in all the datasets considered. As a result, this method is seen as a robust choice for meta-modeling which does not require specific knowledge about parameter tuning in order to obtain good results.

The Inverse Problem
Raissi et al. [72] proposed a physics-informed neural network for solving inverse problems after learning the unknown model parameters, by means of gradient-based learning, of the Navier-Stokes and Korteweg-de Vries equations. Haghighat et al. [73] applied a similar approach to solid mechanic problems. Yan et al. [74] applied extreme machine learning to the analysis of composite structures, demonstrating how the framework can be applied to inverse problems with the aim of finding optimal designs.
A typical challenge that occurs in the design phase is the ability to satisfy a trade-off between mass and stiffness. In order to meet the changing need for various requirements, the possibility of exploring the design space and the optimal compromise between variables are important aspects that are discussed now.
In the present study, the direct problem consists in calculating the mass m, the critical buckling load P critical , and the post-buckling stiffness b I from the input design variables v 1 , . . . , v 5 (Equation (38)). The inverse problem to be solved with the proposed machine learning framework consists of finding the values of the design variables v = [v 1 , . . . , v 5 ] for a given m, P critical , and b I . The solution of the inverse problem is computed as outlined in Algorithm 3 below, where the relative error between the index estimate and the actual index value is computed through Equation (40). The model with the smallest median of error values is chosen as representative of each category: SVM, Kriging, and Random Forest.
While evaluating b I , the error dispersion of the selected methods are compared, as shown in Figure 14. Note that the Random Forest presented the lowest error median; thus, it is used in the computations of the inverse problem discussed later. For the mass evaluation, the error dispersion of the selected methods are compared, as is shown in Figure 15. The Kriging algorithm is used in the computations discussed next because it presented the lowest median error.
For the evaluation of P critical , the error dispersion of the selected methods are compared, as is shown in Figure 16. The Kriging algorithm was selected as the best meta-modeling technique to compute P critical , although the best results of each of the three groups presented similar performances. In this case, the Random Forest or the SVM are used without compromising the result, and both methods show comparable computational costs. The steps performed in adjusting the models for a given dataset-the training phase-are shown in Algorithm 2.
Given a set of target valuesb,m, andp, representing the desired b I , mass, and P critical , the inverse problem is computed by means of the minimization of subject to where α 1 , α 2 , α 3 , and α 4 are weighting parameters, indicating preference in meeting each goal. Although the design values v j that solve the inverse problem are not known at the beginning of the optimization process, an useful estimate is provided by the meta-model G j .

Algorithm 2:
The training phase /* The index of the design variables */ j ← 1, . . . , 5; /* The number of DOE sample points */ n sample ← 2000; n training ← 0.8 × n sample ; Generate a Design of Experiment v i,j ∈ R n training ×5 ; Evaluate the FEM model at v i,j designs; /* The estimate of the index performance b I by means of machine learning */ Adjust a meta-model F b such that b ← F b (v i,j ); /* The estimate of the mass m index by means of machine learning */ Adjust a meta-model F m such that m ← F m (v i,j ); /* The estimate of the P critical index */ Adjust a meta-model F p such that p ← F p (v i,j ); /* A new design estimate by means of machine learning */ for j ∈ {1, . . . , 5} do Adjust a meta-model G j such that v j = G j (b, m, p, v i,k ), i ∈ 1, . . . , n training , k ∈ 1, . . . , 5, k = j; It is worth mentioning that Equation (46) retains the concepts provided by Equation (11), and extends the analysis by considering further information provided by G j . Furthermore, the inequalities explicitly described by Equation (10) are implicit in the learning process of the F b , F m , F p , and G j meta-models. The remaining box constraints (Equation (12)) will also be checked by Algorithm 3. As a result, the proposition of Equations (11) and (12) also constitutes a contribution of the present study.
After the training phase is accomplished (Algorithm 2), the computation of a new design v * j , j = 1, . . . , 5 that solves the inverse problem is given by Algorithm 3. The effectiveness of this methodology in providing feasible designs is discussed in the following.
A Design of Experiments of n sample = 2000 points is generated by means of a Latin hypercube sampling algorithm [60]. The optimization phase is driven by a differential evolution strategy implemented in the DEoptim toolbox [75,76]. The mass m, the postbuckling stiffness b I , and the buckling load P critical of each point are shown in Figure 17, where a correlation between the upper value of P critical and the mass can be noticed. However, there is no apparent correlation between mass and b I , nor is there one between P critical and b I . As a result, it would be of practical interest for the design optimization of the present study to identify designs on this Pareto-like frontier between m and P critical that contains a positive b I . Using no further information about feasible designs, an arbitrary set of experiments are evaluated by this inverse problem algorithm. The entire process of modeling (Algorithm 1), training (Algorithm 2), and optimization (Algorithm 3) is summarized in Figure 18.
The designs obtained as a response to the inverse problem are evaluated by the finite element model, resulting in the values of m and P critical shown in Figure 19. Those results resemble the trend found in the DOE when mass and P critical are compared. Therefore, the meta-models herein applied have proven effective in evaluating new designs in a region of interest.
The comparison between the mass and the post-buckling stiffness b I is presented in Figure 20. The feasible designs contain a positive b I . The methodology was effective in identifying feasible designs in two ranges of mass (m < 2 kg and m > 2.5 kg). The comparison of b I and P critical is shown in Figure 21. Feasible designs (b I > 0) are found in the ranges P critical < 10,000 and P critical > 25,000.  A similar trend is found in the comparison of Figures 20 and 21. Since the mass m and buckling load P critical are chosen in a correlated fashion, this similarity was expected. The design variables and the corresponding b I , m, and P critical are presented in Table 10. Two decimal digits are considered when calculating the values of the design variables v 4 and v 5 throughout the optimization process. This constraint addresses a manufacturing requirement. Moreover, the v 2 design variable has integral values. The corresponding results are shown in columns 2, 4, and 5 of Table 10, respectively.
The correlation between m and P critical in those designs can be verified as shown in Figure 19.
It is worth mentioning the design presented in the second line of Table 10, whose mass is 1.707 kg and b I > 0, confirms that the inverse problem is effectively solved in the discovery of feasible designs with reduced mass requirements. Designs with reduced mass values of 1.448 kg and 1.690 kg, and with positive post-buckling stiffness b I > 0, were also found in lines 12 and 13 of the table, respectively.
The three best optima of Table 10 have the following optimum mass values: m = 1.448 kg, m = 1.510 kg, and m = 1.690 kg; these are illsutrated, respectively, in Figures 22-24. The thickness pattern of these optima have a close resemblance with the imperfection-insensitive hybrid shells proposed by Wagner et al. [77], where hoop-oriented belts were installed on top of a monolithic cylindrical shell to achieve the imperfection-insensitive behavior described by the authors. Furthermore, Lincoln et al. [78] also reported imperfectioninsensitive designs with hoop-oriented thickness patterns created by CTS, and the encountered optima reported therein also resemble what the inverse ML-driven solver found in the present study. The nonlinear behavior of the optimum configuration with m = 1.448 kg was evaluated using the ABAQUS finite element solver using a mesh of 240 general-purpose elements S4R around the circumference, which have 4 nodes and reduced integration, while keeping the element aspect ratio close to 1. The results are compared against two constantstiffness cases, where the cylinders were produced by setting θ 1 = θ 2 and equal to 0 • and 45 • , respectively. For the latter, the thickness of the plies is higher by a factor of 1/ cos 45 • due to the nonlinear steering-thickness coupling produced by the CTS manufacturing process, as previously explained.

Imperfection-Sensitivity Analysis
The imperfection sensitivity of the optimal design obtained after solving the aforementioned inverse problem is evaluated using single-perturbation load imperfections (SPLI) [79]. The performance is compared with two designs: one design named "Case 0 • ", having two plies with all the fibers oriented parallel to the axial direction; and another named "Case ±45 • ", with two plies oriented at ±45 • . Note that, in the latter design, the thickness is increased by a factor of 1/ cos 45, assuming manufacturing by means of CTS. All analyses are performed using ABAQUS and 240 S4R elements around the circumference, as detailed in Castro et al. [80]. The name boundary conditions used for the optimization are herein realized by fixing the bottom-and top-edge nodes in the y and z directions, and by fixing the translation in the x direction for all nodes laying in the middle cross-section of the cylindrical shell. Note that this renders a load-controlled compression analysis, where the nodes at the edges are allowed to move separately from each other in the axial direction. The results of the imperfection sensitivity analysis are presented in Table 11, where the optimum cylinder herein obtained (design ID = 12 of Table 10) shows a reduced imperfection sensitivity when compared to Case 0 • and Case ±45 • . In summary, a new methodology was proposed that encompasses modeling (Algorithm 1), training (Algorithm 2), and optimizing (Algorithm 3). Numerical computation confirmed the feasibility of the proposal, whose main advantages are: • effective modeling strategy that explores the design space; • uses popular machine learning methods without needing to adjust special parameters; • optimization proposition that leads from the Design of Experiments ( Figure 17) to optimal designs on the boundary of the compromise between mass and stiffness (Figures 19-21).
Eighteen machine learning instances were evaluated on three datasets. No single strategy demonstrated superior performance in all cases (Figures 14-16). On the other hand, sub-optimal machine learning strategies were also effective in supporting the optimization process. This confirms the need to select a reliable workflow ( Figure 18) that will not degrade due to poor machine learning performance.
As a result, the authors believe that the present methodology constitutes a contribution to the design phase of a new project.

Conclusions
The present study proposed a design parameterization compatible with the continuous tow shearing (CTS) manufacturing process that is able to produce lightweight variable-stiffness cylindrical structures with tailored buckling load and reduced imperfection sensitivity.
The new methodology consists of contributions in the areas of problem modeling, the selection of machine learning strategies, and the proposition of an optimization formulation that results in optimal designs around the compromise frontier between mass and stiffness.
The parameterization developed herein took into account the nonlinear steeringthickness variation, which resulted in hoop-oriented regions of variable thickness for steering angles larger than zero. The maximum number of regions of larger thicknesses depends on the minimum fiber-steering radius, here assumed to be 50 mm, reportedly achievable by the CTS process.
An optimization goal to minimize the mass of cylindrical shells while keeping their post-buckling stiffness positive was proposed. The investigated problem consisted of a cylindrical shell under load-controlled axial compression, where both edges are simply supported and allowed to deform unevenly in the axial direction under the load application. The shells should be able to withstand a given design load, which is then compared with its linear buckling load. The post-buckling stiffness and linear buckling load are assessed by means of an enriched finite element model, using a formulation based on 10 degrees of freedom per node and a third-order interpolation of all the displacement field variables. For the post-buckling stiffness, a displacement-based formulation of Koiter's method was adopted.
The machine learning models enabled a new design perspective based on the solution of a highly nonlinear inverse problem, and this proposition was effective in determining the input variables from a pre-determined design load. It is worth noting that all three of the methods performed well when using a suitable kernel choice, and the election of the best methodology would change if another criteria is considered, such as computational time or maximum error.
In addition to the quality of the designs provided by the methodology, as can be seen in Table 10, the formulation does not require additional parameter adjustments. Therefore, the investigation of the design space became straightforward.
Further studies will focus on applying Koiter's method from a nonlinear pre-buckling state while looking for fully imperfection-insensitive shell designs that are capable of reaching post-buckled stiffness in load-and displacement-controlled applications. Furthermore, a detailed reliability-based evaluation of the obtained shells will be performed to map how robust the obtained positive post-buckling stiffness is in terms of angle and thickness imperfections, as well as variations in material properties. The evaluation of buckling behavior after an imperfection created post-impact, which was investigated, for instance, by Maziz et al. for filament-wound hybrid composite pipes [4], would be an interesting next step.