Ludovic Coelho
email: ludovic.coelho@onera.fr

Christian Fagiano
email: christian.fagiano@onera.fr

Julien Cédric
email: cedric.julien@onera.fr

Dmas

Nicolo Fabbiane
email: nicolo.fabbiane@onera.fr

Didier Lucor
email: didier.lucor@lisn.upsaclay.fr

MULTI-SCALE APPROACH FOR RELIABILITY-BASED DESIGN OPTIMIZATION WITH METAMODEL UPSCALING

Keywords: Reliability-based design optimization, multiscale modeling, uncertainty, Inverse problem, Surrogate models, composite material

For multi-scale materials, the interplay of material and design uncertainties and reliability-based design optimization (RBDO) is complex and very dependent on the chosen modeling scale. Uncertainty quantification and management are often introduced at lower scales of the material, while a more macroscopic scale is the preferred design space at which optimization is performed. How the coupling between the different scales is handled strongly affects the efficiency of the overall model and optimization. This work proposes a new iterative methodology that combines a low-dimensional macroscopic design space with gradient information to perform accurate optimization and a highdimensional lower-scale space where design variables uncertainties are modeled and upscaled. An inverse problem is solved at each iteration of the optimization process to identify the lower-scale configuration that meets the macroscopic properties in terms of some statistical description. This is only achievable thanks to efficient metamodel upscaling. The proposed approach is tested on the optimization of a composite plate subjected to buckling with uncertain ply angles. A particular orthonormal basis is constructed with Fourier chaos expansion for the metamodel upscaling, which provides a very efficient closed-form expression of the lamination parameters statistics. The results demonstrate a drastic improvement in the reliability compared to the deterministic optimized design and a significant computational gain compared to the approach of directly optimizing ply angles via a genetic algorithm.

Introduction

The structural design process often results in an optimization problem for which the structural performance is maximized within a set of constraints imposed by the materials and general behavior of the structure. The result of the optimization is strongly dependent on the type of modeling that is considered (e.g., atomistic vs. continuum modeling or deterministic vs. stochastic), in part because the nature of materials is broad. Design variables may exist at different scales of the material and also because it needs to be better understood how to account for the many materials, structural or modeling uncertainties of the system. Complex nonlinear materials such as newly architected materials (e.g. polymeric nano-composite, fiber-reinforced polymers) or porous materials often cover a wide range of different length scales, i.e. the scales of the heterogeneities from nano-, micro-, to mesoscale, whose problematic optimization due to large design spaces strongly influences the macroscale physical and mechanical properties. Once focus is tied to a particular scale, the level of details and the objective function's regularity may impeach the optimal design's performance. In topology optimization, multi-scale approaches are proposed to account for scale-related effects of materials and structures [START_REF] Zhang | Scale-related topology optimization of cellular materials and structures[END_REF] and simultaneously achieve the layout updating and the material change by macro-micro concurrent design [START_REF] Long | Concurrent topology optimization of composite macrostructure and microstructure constructed by constituent phases of distinct Poisson's ratios for maximum frequency[END_REF]. [START_REF] Gao | Topology optimization for multiscale design of porous composites with multi-domain microstructures[END_REF] proposes a multi-scale method for porous composites where the overall distribution in the macrostructure and the material microstructures are optimized. For simple isotropic cases, methods of inverse homogenization were also introduced [START_REF] Sigmund | Materials with prescribed constitutive parameters: An inverse homogenization problem[END_REF] to find the internal topology of a base (truss) cell with given homogenized coefficients. In the field of metamaterials, another inverse problem must be solved with the same idea of matching macro physical properties with microstructure optimization. More general than periodic truss-, plate-, or shell-based architectures, [START_REF] Kumar | Inverse-designed spinodoid metamaterials[END_REF] takes advantage of recent machine learning capabilities to achieve the robust inverse design of the non-periodic class of spinodoid topologies, used to represent biomimetic artificial bone architectures. Laminate composite optimization involves many design variables at a mesoscale and requires evolutionary algorithms [START_REF] Venkataraman | Optimization of composite panels -a review[END_REF], still struggling because of the multimodal nature of the functions of interest. Homogenization procedures can be used to introduce a macroscopic space in which the functions of interest are mostly regular. Bi-level approaches are commonly used, relying on a first optimization of the laminate in the macroscopic space (e.g., via gradient-based methods) while the stacking sequence is then retrieved in a second step (mostly via genetic algorithms) [START_REF] Macquart | Derivation and application of blending constraints in lamination parameter space for composite optimisation[END_REF][START_REF] Scardaoni | New blending constraints and a stack-recovery strategy for the multi-scale design of composite laminates[END_REF]. [START_REF] Ferreira | Hierarchical optimization of laminated fiber reinforced composites[END_REF] use a hierarchical optimization to simultaneously design macro-and microstructural levels of plies orientations and fiber volume fractions and also in terms of fiber cross-sectional shapes, respectively. A clear separation of scales is not always possible, and it is sometimes unclear (both from a physical and a computational point of view) at which scale(s) the optimization should take place. Moreover, while the forward problem, i.e., mapping the macroscale properties from finer (e.g., micro-or mesostructure) properties, is well defined, the inverse problemidentifying a small-scale structural topology or configuration that meets the mechanical property requirements -is ill-posed (multiple solutions can have the exact same properties). Moreover, the idea of a homogenized description of the system [START_REF] Charalambakis | Homogenization Techniques and Micromechanics. A Survey and Perspectives[END_REF][START_REF] Gineau | Macroscopic model of fluid structure interaction in cylinder arrangement using theory of mixture[END_REF], which seems computationally useful, has to be handled with care in order to benefit to the optimization process.

In addition to the various scales considerations, for some materials, as manufacturing processes are complex, some uncertainties (e.g., parameters, geometry, material and mechanical properties, molecular interactions, operating conditions, lack of repeatability) may creep at different scales into the system (material uncertainties, structural uncertainties, modeling uncertainties) and have a significant influence on the overall structural performance. Design variables and other system parameters may also contain some inherent random quantities. Machine learning (ML) methods, defined here in a broad sense, have gained particular interest in the computational mechanics community towards developing surrogate models with uncertainty quantification as ML is effective in solving complex inverse and ill posed problems while also suitable for optimization query applications, [START_REF] Cheng | Machine learning with data assimilation and uncertainty quantification for dynamical systems: a review[END_REF]. For deterministic optimization, such uncertainties are often dealt with by introducing simplifying hypotheses, for instance, using safety factors, considering only average or extreme values. However, such approaches are inefficient and can lead to conservative and inefficient designs or optimistic designs with poor reliability [START_REF] Beck | A comparison of deterministic, reliability-based and risk-based structural optimization under uncertainty[END_REF]. The research topic of optimization under uncertainty (OUU) is vast and presents important computational challenges in terms of implementing efficient numerical procedures. A crucial initial step is carefully crafting the OUU problem formulation, conditioned by the evolution of the system from a deterministic design problem to a design under uncertain conditions. Stochastic formulations of the objective and constraint functions must be carefully described.

Robustness is habitually associated with considering uncertainties in the objective function, and reliability-based optimization with introducing uncertainties in the constraint functions. Reliability turns out to be an increasingly critical issue in structural design [START_REF] Kuschel | Two basic problems in reliability-based structural optimization[END_REF][START_REF] Reddy | Reliability-based structural optimization: a simplified safety index approach[END_REF][START_REF] Lopez | Reliability-based design optimization strategies based on FORM: a review[END_REF][START_REF] Moustapha | Quantile-based optimization under uncertainties using adaptive kriging surrogate models[END_REF]. In this work, the focus is made on design variable uncertainties for a Reliability-Based Design Optimization (RBDO). The challenge is to find an optimal reliable solution such that the probability of the design not satisfying constraints does not exceed a chosen limit failure probability. It is an active multidisciplinary research topic, and some reviews can be found in [START_REF] Aoues | Benchmark study of numerical methods for reliability-based design optimization[END_REF][START_REF] Yao | Review of uncertainty-based multidisciplinary design optimization methods for aerospace vehicles[END_REF][START_REF] Lelièvre | On the consideration of uncertainty in design: optimization -reliability -robustness[END_REF][START_REF] Acar | Modeling, analysis, and optimization under uncertainties: a review[END_REF]. Taking account of the uncertainty rely on repeated evaluations of the mechanical model, mainly in the reliability analysis step. The range of applications may be limited due to this optimization's time-consuming computation. Surrogate models have been used in order to alleviate the computation time. For the RBDO task, surrogate models relying, for instance, on Polynomial Chaos Expansion were considered in [START_REF] Suryawanshi | Reliability based optimization in aeroelastic stability problems using polynomial chaos based metamodels[END_REF], [START_REF] López | Reliability-based design optimization of composite stiffened panels in post-buckling regime[END_REF], or Kriging in [START_REF] Dubourg | Reliability-based design optimization using kriging surrogates and subset simulation[END_REF], [START_REF] Li | A local kriging approximation method using MPP for reliability-based design optimization[END_REF]. Another interesting approach is the one of [START_REF] Moustapha | Quantile-based optimization under uncertainties using adaptive kriging surrogate models[END_REF], where a new quantile-based formulation is proposed motivated by the relatively high target failure probabilities that can be accepted in the car body design field.

In the multidisciplinary optimization field, some works have been done in optimization under uncertainty using gradient-based algorithms and surrogate model strategies to reduce the time of convergence [START_REF] Lucor | Stochastic design optimization: Application to reacting flows[END_REF][START_REF] Andrieu | Gradient-based simulation optimization under probability constraints[END_REF][START_REF] Dubourg | Reliability-based design optimization using kriging surrogates and subset simulation[END_REF][START_REF] Díaz | Efficient methodologies for reliability-based design optimization of composite panels[END_REF][START_REF] Fang | A gradient-based uncertainty optimization framework utilizing dimensional adaptive polynomial chaos expansion[END_REF]. In multi-scale framework, recent works [START_REF] Duan | Reliability-based multi-scale design optimization of composite frames considering structural compliance and manufacturing constraints[END_REF], Jung et al., 2021[START_REF] Liu | A general integrated procedure for uncertainty-based design optimization of multilevel systems by hierarchical decomposition framework[END_REF] consider "concurrent" design variables and uncertainties at different scales. No matter their proposed formulation, all of these works point to the severe computational expense due to the dimensionality of the uncertainties to be propagated/quantified and the function evaluations required to assess the reliability. Many of them rely on the use of surrogate models in order to alleviate this cost. Different sequential strategies and algorithmic loops are proposed to decouple and alleviate the optimization, uncertainty quantification, and system response evaluations. [START_REF] Liu | Learning macroscopic internal variables and history dependence from microscopic models[END_REF] use a deep neural network to learn macroscopic internal variables and history dependence of a polycrystalline from microscopic models. It gives access to multi-scale accuracy for a computational cost comparable to a conventional empirical constitutive relation. Other surrogate techniques can be used. For instance, [START_REF] Ghasemi | Uncertainties propagation in metamodel-based probabilistic optimization of CNT/polymer composite structure using stochastic multi-scale modeling[END_REF] created a kriging metamodel for the multi-scale uncertainty propagation model of CNT/polymer structure. Based on a similar concept, [START_REF] Omairey | Multiscale surrogate-based framework for reliability analysis of unidirectional FRP composites[END_REF] used a surrogate model to establish the relationship between microscale uncertainties and macro-scale material property uncertainties, thereby reducing computational costs and improving computational efficiency. This is applied to RBDO of unidirectional fiber-reinforced polymer (FRP) composite laminates in [START_REF] Omairey | Multi-scale reliability-based design optimisation framework for fibre-reinforced composite laminates[END_REF] where uncertainties are taken into account at different scales of the composite. In their case, the smallest scale design space is optimized by evolutionary algorithms (i.e., particle swarm optimization). A similar strategy was used in [START_REF] Liu | Reliability-based design optimization of composite battery box based on modified particle swarm optimization algorithm[END_REF] with the optimization method, and a different surrogate modeling strategy.

Our research proposes a multi-scale RBDO methodology where the meaningful design variables are (discrete) material properties at a mesoscopic scale. The multi-scale approach is exploited to use an efficient gradient-based optimization method in the homogenized space. The originality of this work is that the multi-scale strategy is used at each iteration of the optimization to propagate the uncertainty inherited from the mesoscopic design variables. Indeed, the parametric uncertainties affect both design variables at meso-and macro-scale. Moreover, our approach is sequential: gradient-based optimization is pursued in a "deterministic" fashion and provides a homogenized macroscale optimum design that we treat as a "mean" design. A critical step is the one of solving an inverse problem of identifying a corresponding uncertain mesoscale design that matches our mean target design. A surrogate model makes quick statistical homogenization possible and efficiently delivers the needed macroscopic parameter statistics. Finally, the proposed methodology can handle complex geometries and several constraints and converges in a decent computational time.

The core of the present manuscript is Section 2, where the multi-scale and surrogate-based strategy to accelerate the solution of the inverse problem is presented. Hence, the proposed approach is applied to the optimization of a laminated composite material subject to uncertainties at the mesoscale: the test-case is introduced in Section 3, and the main results are presented in Section 4. Section 5 discusses the sensitivity of the optimization process to the initial design and the to the computation of the gradient of the failure probability. Section 6 will conclude the paper and provide some perspectives.

Multi-scale reliability-based design optimization approach

In general RBDO problems, both design variables and other system parameters can contain deterministic and/or random quantities, here represented as random variables. For simplicity, θ is defined as the vector of design variables, and p describes the environmental parameters. Moreover, subscripts r and d will refer to random and deterministic quantities, respectively. These quantities can be more compactly written as vectors as follows1 :

θ = Θ r θ d , p = P r p d , R = Θ r P r , d = θ d p d . (1)
Reliability analysis quantifies structural safety considering the random nature of all phenomena affecting the structural system. A (or more) performance function g(R, d) characterizes the system response. We emphasize that this function is also random due to its dependency on R. The design region is divided into two classes of domains:

Failure domain:

F = {R, d | g(R, d) ⩾ 0} , (2)
and the Safety domain:

S = {R, d | g(R, d) < 0} . (3)
The boundary between failure and safety domains is the limit state surface, which generally is a hypersurface in the n-dimensional space of random variables R. Other (deterministic) constraints may exist as well. According to this, the failure probability P is formulated as:

P (g(R, d) ⩾ 0) = • • • F π R dR, (4)
where π R is the joint probability density function of the random variables. Nevertheless, except in some particular cases, the integral expression cannot be computed analytically. Indeed, this is often out of reach because of the nonlinearity of π R , the number of random variables which can be large, and the lack of knowledge about the exact position of the failure domain.

Many different formulations and computational algorithms exist for RBDO problems. For instance, the double loop approaches consider the reliability constraints within the optimization loop [START_REF] Nikolaidis | Reliability based optimization: A safety index approach[END_REF]]. Mono-level approaches exist where the probabilistic constraint is approximated with deterministic values, converting the double loop into a single loop [START_REF] Kuschel | Two basic problems in reliability-based structural optimization[END_REF]. Finally, uncoupled approaches solve sequentially deterministic optimization procedures with a reliability analysis at the end of each optimization [START_REF] Du | Sequential optimization and reliability assessment method for efficient probabilistic design[END_REF].

In this work, a double loop approach is considered to control the uncertainties on the failure probability during the whole optimization process. Moreover, because some design variables are considered random, we will optimize some of their moments. This will be made easier by making some simplifying assumptions about the joint distribution of the random variables designating this part of the design. For instance, we can assume that Θ r follows a generic class of distributions for which one (or more) parameters are not fixed (for instance, the mean value for a Gaussian distribution).

In this case, the formulation of the problem is taken from [START_REF] Díaz | Efficient methodologies for reliability-based design optimization of composite panels[END_REF] and can be formalized as:

min µΘ r ,θ d f (µ R , d) subject to: h i (µ R , d) ≤ 0 i = 1, ..., n d P (g j (R, d) ⩾ 0) ≤ P max j j = 1, ..., n p (5)
where f is the deterministic cost function to be minimized, µ R is the vector gathering the random variables mean values, d the deterministic parameters, µ Θr the mean of random design variables, θ d the deterministic design variables, h i are the n d deterministic constraints, and P are the n p failure probabilities of the limit state functions g j , which have to be below the maximum failure probabilities P max j .

In this work, the design variables only bear a random component Θ r due to complex manufacturing processes, and the environmental parameters do not contain any random variable. For the sake of simplicity, the functions f , g or h will now be written as function of Θ or µ Θ .

Multi-scale formulation

In this work, we propose a multi-scale RBDO approach. Without lack of generality, we will consider two adjacent scales only, i.e., meso and macroscale, for which we know a nonlinear homogenization process mapping the mesoscale to the macroscale description of the material. The design variables that must be optimized for manufacturing and to provide reliable performance are the mesoscale variables, called Θ. We wish to account for the intrinsic uncertainties of these variables. We assume that we have a good model for these uncertainties. Therefore, the uncertainty of the mesoscopic variables is mostly known and modeled with a probability density function π Θ given some parameters.

On the other hand, the uncertainties associated with the continuous macroscale variables, called v, are unknown and must be quantified. Moreover, it can be complex to model due to various correlations. As we will see, this will directly impact the calculation of the failure probability for reliability optimization in this space. Despite this complication, it remains more efficient to take advantage of the homogenized space where gradient-based algorithms can accelerate the convergence of the optimization problem. Our idea is to take an iterative approach, repeatedly moving from one scale to the other during the optimization, in order to take advantage of facilitating properties at all scales and deal with the uncertainty quantification upscaling. Nevertheless, computational burdens will remain at each iteration of the optimization process, such as evaluating the probability of failure and its gradient.

The macro (homogenized) space is utilized for the global optimization process to rely on a gradient algorithm. This strategy, exploiting both spaces for calculating the failure probability, is illustrated in Fig. 1. The objective function and the constraints are defined, and the optimizer improves the design in this macroscopic space -cf. the red line in Fig. 1. The original mesospace is going to be used to help in evaluating the failure probability because the distribution of the discrete design variables is known. Nevertheless, the inverse mapping from the continuous design variables to the discrete one is far from being trivial. It is a highly multimodal optimization problem where it is necessary to use metaheuristic optimization methods. In this study, the solution of the inverse problem -illustrated by the dotted pink line in Fig. 1 -is solved thanks to a genetic algorithm. Roughly speaking, once the genetic algorithm identifies a potent set of discrete mesoscopic variables, it is conceptually straightforward to propagate the uncertainty to the homogenized space, then to the model in order to calculate the failure probability (.e.g, by the Monte Carlo method).

Then the gradient optimizer can propose another design point and repeat the process until convergence is reached.

More specifically, the formulation of our multi-scale reliability-based design optimization taking into account design variables uncertainties can be described by Eq.( 6).

min µ Θ f µ v(Θ|µ Θ ) subject to: h i µ v(Θ|µ Θ ) ≤ 0, i = 1, . . . , n d P j (g j (v = H(Θ)) ⩾ 0) ≤ P max j , j = 1, . . . , n p (6)
where Θ are the random design variables with a joint probability density function π Θ|σΘ given their variances, µ Θ are the design parameters mean values to optimize, h i are deterministic nonlinear constraints (emerging from the homogenization process), g j are limit state functions, H is the nonlinear mapping2 from the mesoscopic space to the macroscopic space, used to obtain v for a given µ Θ , and P j are the failure probabilities evaluated in the macroscopic space.

Because of the design variables uncertainties in both spaces, it is important to efficiently propagate forward uncertainties from the original mesospace to the target macrospace, but also to be able to solve an inverse mapping that satisfies certain statistical constraints: e.g., given a mean macroscale design, be able to identify the original uncertain mesoscale design that generates a population of macroscale designs with the corresponding mean solution. A trial-and-error approach to resolve this inverse problem would be highly costly. Therefore, we propose to build a surrogate mapping Ĥ that is constructed in place of H in order to rapidly access needed statistics of v for a given population of variables Θ. The type of deployed surrogate modeling is obviously problem-dependant but does not negatively impact the formulation as long as it is sufficiently accurate. Moreover, a strategy is proposed to approximate the failure probability and its gradient efficiently. Technical details are presented in the following subsections.

Surrogate-based assistance strategy for the inverse problem resolution

In the following, we consider the uncertainty of the input parameter Θ as being defined as:

Θ = µ Θ + σ Θ X, (7)
where µ Θ is the unknown mean vector of Θ, σ Θ is the known vector standard deviation, and X is chosen as an independent and identically distributed (i.i.d.) random vector, described by an independent joint Gaussian distribution π X with zero mean and unit variance. For a given state in the homogeneous space provided by µ v , the inverse problem optimization to be solved in the stochastic context is written as:

µ v target := µ v min µΘ ∥µ v target -µ H(Θ|µΘ) ∥, (8)
where µ H(Θ|µΘ) is the mean of the random quantity representing the image through the mapping H from the mesoscopic to the macroscopic space. For the inverse problem resolution, the distance is here minimized with respect to the mean of the mesoscopic and macroscopic parameters for simplicity3 . The distance is here minimized -instead of being made null -to account for a very general scenario in which there does not exist an exact corresponding mesoscale design matching the target. This situation is common in cases where the mean vector of the mesoscale design variables can only take some discrete values from a set of finite configurations. It remains that the several computations of the µ H(Θ|µΘ) vector (needed for the inverse problem) and higher-order correlation statistics (needed for the failure probabilities and their gradients) call for the use of a mapping approximation Ĥ. In this work, we will rely on generalized polynomial surrogate models, referred here under the generic naming of polynomial chaos expansions (PCE) [START_REF] Xiu | The wiener-askey polynomial chaos for stochastic differential equations[END_REF] and used under various forms for mechanical systems, e.g. [START_REF] Bijl | Uncertainty Quantification in Computational Fluid Dynamics[END_REF][START_REF] Garroussi | Tackling random fields non-linearities with unsupervised clustering of polynomial chaos expansion in latent space: application to global sensitivity analysis of river flooding[END_REF]. It is an effective tool that guarantees exponential convergence with

g(v) > 0 v 1 v 2 Θ 1 Θ 2 Θ 3 v k v k+1 Θ k Θ k+1
macro space meso space

g ( v ) = 0 v k ≈ Ĥ ( Θ k ) μ Θ k μ Θ k+1 μ v k μ v k+1 μ Θ 1 μ Θ 2 μ Θ N v k+1 ≈ Ĥ ( Θ k+1 )
Inverse problem metamodel Failure domain

Figure 1: Cartoon illustrating the proposed sequential multi-scale RBDO approach. Design variables bear some uncertainties in the meso space. For a given design, surrogate homogenization allows very efficient and rich uncertainty quantification (bottom large blue arrow) of the mapped image in the macro space of the lower dimension. Gradient RBDO optimization is then carried out on the mean value (solid red arrow). An inverse problem is solved to identify the best corresponding mesoscale design (dashed pink arrow). This involves the combination of metaheuristic optimizations and forward metamodeling (top large blue arrow).

The whole process is re-iterated until convergence.

increasing expansion order for sufficiently smooth functionals. Moreover, once the approximation is obtained, the access to the statistics is very straightforward from the expansion coefficients. Nevertheless, the plain version of this method is inaccurate for models which are highly nonlinear or are of very high dimensionality. Other options are possible, including conditioned Gaussian processes like Kriging [START_REF] Rasmussen | Gaussian processes for machine learning[END_REF], which is employed for larger experimental areas such as Bayesian optimization, approximating deterministic function, and machine learning. It combines a regression model with a stationary Gaussian process error model. The advantage is that it provides a stochastic error bound.

Assuming the quantity of interest H is a second-order random variable, then it can be expanded into an infinite series of polynomials:

Ĥ(Θ(X)) = ∞ i ξ i ϕ i (X), (9)
where ϕ i must be orthonormal polynomials with respect of the measure of X, ξ i are the expansion coefficients. In practice, the expansion must be truncated, and a choice must be made on the number of retained terms for sufficient convergence. Thanks to the orthonormality of the approximation space, the deterministic coefficients are obtained as projections of the function of interest onto each member of the approximation basis.

With the Gaussian distribution hypothesis on X, Hermite polynomials can be used [START_REF] Wiener | The homogeneous chaos[END_REF]. Later, another basis (made of polynomials of trigonometric functions) is introduced, which is much more efficient for the application. This surrogate model is very convenient to approximate the statistics of H(Θ(X)|µ Θ ) needed for the inverse problem and the computation of the failure probability and its gradient.

Here, we have focused on metamodeling the upscaling process alone. Obviously, such an approach could be deployed at several other levels. For instance, the evaluation of the problem constraints could be costly and rely on some type of high-order finite elements model and may also benefit from metamodeling. This will not be pursued in this paper.

Computation of failure probability and its gradient

Once the mesoscopic set variables are identified at a given step of the RBDO process, one can proceed with the uncertainty quantification given the distribution π Θ . The failure probability in Eq.( 6) can be approximated by the Monte Carlo method (MC) as follows:

P g(v = Ĥ(Θ)) ⩾ 0 = 1 n M C n M C i=1 I F v = Ĥ(Θ (i) ) , (10)
where Θ (i) , i = 1, . . . , n M C is a large sample of n M C independent copies of the random vector Θ and I F is the failure indicator:

I F v = Ĥ(Θ (i) ) = 1 if g v = Ĥ(Θ (i) ) ⩾ 0 0 otherwise. ( 11
)
In order to assess the accuracy of the estimator and deduces a reasonable sample size to use, it is interesting to resort to the coefficient of variation of the estimator defined as follows:

CoV = σ P = 1 -P n M C P . ( 12
)
The sampling-based RBDO approach relying on the Monte Carlo simulation is general but requires intensive computation. This is particularly crucial for reliability-based optimization where the threshold is low, and many samples are needed for these rare events, making the point for surrogate modeling. The computational effort is even more strenuous to calculate the sensitivities of the probabilistic response.

Since the macroscopic variables are used for the smoothness of the functions of interest, a gradient algorithm is used.

The challenging part resides in the computation of the failure probability gradient. Several approaches exist in the literature, such as the finite differences, the score function approach, or the Smooth Perturbation Analysis methods [START_REF] Fu | Smoothed perturbation analysis derivative estimation for markov chains[END_REF]Hu, 1994, Rubinstein, 1986], well summarized for Monte-Carlo approaches in machine learning in [START_REF] Mohamed | Monte carlo gradient estimation in machine learning[END_REF]. Finite difference methods are straightforward for computing failure probability gradients but may produce inaccurate sensitivities for too small sample size due to the statistical noise, or because an appropriate step size (or perturbation size) for each design variable needs to be determined, which is often arbitrary.

When the uncertainty of macroscopic parameters can be modeled with a parametric probability density function, the score function approach [START_REF] Rubinstein | The score function approach for sensitivity analysis of computer simulation models[END_REF] can be used for the failure probability sensitivity computation. This approach was brought to the structural reliability community by [START_REF] Wu | Computational methods for efficient structural reliability and reliability sensitivity analysis[END_REF]. Assuming that the joint PDF π v of the continuous parameters is continuously differentiable with respect to µ vi and that the integration range V, i.e., the design space domain, does not depend on µ vi , the partial derivative of the failure probability recasts as follows:

∂P ∂µ vi = V I F (v) ∂π v (v, µ v ) ∂µ vi dv. (13)
Then, in order to compute this integral more efficiently as an expectation, [START_REF] Rubinstein | The score function approach for sensitivity analysis of computer simulation models[END_REF] proposed to use an importance sampling trick. Given a sample V = v (i) = Ĥ Θ (i) , i = 1, . . . , n M C of n M C copies of the random vector V coming from the uncertainty of the mesoscopic parameters, the following estimator

∂P ∂µ vi ≡ 1 n M C n M C i=1 I F (v (i) ) π v (v (i) , µ vi ) ∂π v (v (i) , µ v ) ∂µ vi , (14)
is unbiased and asymptotically convergent according to the central limit theorem. The advantage is that the failure probability gradient is estimated using the same sample used for the failure probability estimation. However, the estimator variance of the score function approach is sensitive to the function probability function [START_REF] Mohamed | Monte carlo gradient estimation in machine learning[END_REF]. [START_REF] Zhu | A new sampling-based RBDO method via score function with reweighting scheme and application to vehicle designs[END_REF] proposed a reweighting scheme to improve the accuracy and decrease the variance of the score function.

µ Θ [0] Initial guess µ [0] v = µ H(Θ)
Objective and its gradient f [0] , ∇f [0] Constraints and its gradient h [0] , ∇h [0] , P [0] , ∇P [0] Mesoscopic scale (discrete)

Macroscopic scale (continuous)

Gradient optimizer f [j-1] , ∇f [j-1] , h [j-1] , ∇h [j-1] P [j-1] , ∇P [j-1] Next design point µ

[j] v target ← µ [j] v
Inverse problem Eq.( 8) Surrogate model assistance

Parameters retrieval µ Θ [j] Sample generation µ Θ [j] + σ Θ X (i) , i = 1, . . . , n MC } Homogenization transport Eq.(9) Surrogate model assistance Direct statistics µ v * , σ v * and sample V [j] = v * (i) , i = 1, . . . , n MC Evaluation of determin- istic objective/constraint
and their gradients f [j] (µ v * ), ∇f [j] (µ v * ) h [j] (µ v * ), ∇h [j] (µ v * )

Evaluation of failure probability P [j] with Monte Carlo Eq.( 10) and its gradient ∇P [j] with Score Function Eq.( 14) Depending on the application, when π v is not known a priori, statistical tests can be made on a sample to know if it can be fitted with a parametric law. In the case of a non-parametric distribution, the gradient of the failure probability may be approximated by centered finite differences. One of the drawbacks is the number of simulations that is of the order of O(2n v n M C ) for a random vector of n v dimension. One simple way to reduce the variance of the estimator is to use Common Random Numbers (CRN) [Royset andPolak, 2004b,a, Taflanidis, 2007].

Converge ? Done j = 1 yes Θ (j) = µ Θ (j) + σ Θ X X ∼ N (0, 1) j = j + 1
The whole multi-scale RBDO approach proposed in Section 2 is summarized in the flowchart of Fig. 2. The optimization convergence is reached when all the constraints are respected, and the norm of the macroscopic design variables between two iterations is below ϵ = 0.01. 3 Application to a composite material: test-case of the buckling of a laminated plate with uncertain layers orientations Laminated composite materials are widely used in aerospace, automotive, shipping, and civil engineering thanks to their great material properties for high specific strength and stiffness. In an RBDO framework, Conceição António and Hoffbauer [2017] proposed a new methodology using a genetic algorithm. Here, the loads imposed on the structure are sources of uncertainties. Subsequently, a sensitivity analysis was performed on the optimal structure to study the influence of parameters and design variables on the structural response. The most influential parameters were found to be the modulus of longitudinal elasticity and the ply orientations. However, the uncertainty on ply orientations was not considered during the RBDO. [START_REF] Scarth | Reliability-based aeroelastic design of composite plate wings using a stability margin[END_REF] investigated the minimization of the probability of the flutter instability, modeled as Gaussian processes, in a simple, composite-plate wing, with random design variables, i.e., the ply orientations. Even if the resulting design is, in some sense, an optimized solution with respect to the onset of the flutter instability, their proposed procedure is far from being representative of the real aircraft design process since it is a mono-objective optimization where the aeroelastic stability is taken into account as an objective instead of a constraint. Moreover, a genetic algorithm was used to solve the optimization, which can be limited to complex high-dimensional problems.

In the following, we will rely on gradient-based algorithms for our laminate optimization in the macroscopic space. The multi-scale approach is applied to reliability-based design optimization of a composite plate subjected to a compressive load, where the ply orientations design variables are uncertain due to some given dispersion around design values, as represented in Eq.( 7). The objective is to increase the stiffness of the wing while remaining reliable with respect to the buckling phenomenon. The buckling response has a modal behavior even in the macroscopic space. Since Scarth et al.

[2014] and [START_REF] Scarth | Reliability-based aeroelastic design of composite plate wings using a stability margin[END_REF] show the interest in taking into account the uncertainty of the orientations in an aeroelastic design process in order to avoid the existing discontinuity zones, the buckling is an interesting constraint to apply a RBDO process to show this new methodology.

Composite model

We consider a simplified composite wing, represented as a flat rectangular cantilever plate. The plate dimensions and the applied load direction are shown in Fig. 3. The composite laminate is t = 2-mm thick, which accounts for a total of 16 plies, each stacked at a specific θ i orientation with respect to the global coordinate system. Available orientations to be chosen from are, in general, uniformly distributed over [-75 o : ∆ inc θ : 90 o ], with ∆ inc θ the angular increment (e.g., 15 o ). Table 1 shows the dimension and material properties.

In classical lamination theory [START_REF] Tsai | Introduction to composite materials[END_REF], the macroscale constitutive equation relating applied bending moments to the curvature of a symmetrically laminated plate may be written as

M = Dκ, ( 15
)
where M is a vector of resultant bending moments, and κ is the vector of plate curvatures. The out-of-plane stiffness matrix D can be written as a linear combination of the Tsai-Pagano parameters U, the lamination parameters v D and the thickness t [Tsai andPagano, 1968, Miki and[START_REF] Miki | Optimum design of laminated composite plates using lamination parameters[END_REF]

       D 11 D 22 D 12 D 66 D 16 D 26        = t 3 12        1 v D 1 v D 3 0 0 1 -v D 1 v D 3 0 0 0 0 -v D 3 1 0 0 0 -v D 3 0 1 0 v D 2 /2 v D 4 0 0 0 v D 2 /2 -v D 4 0 0             U 1 U 2 U 3 U 4 U 5      ( 16
)
where the Tsai-Pagano parameters U only depend on material properties and are defined as:

     U 1 U 2 U 3 U 4 U 5      = 1 8      3 3 2 4 4 -4 0 0 1 1 -2 -4 1 1 -6 -4 1 1 -2 4         Q 11 Q 22 Q 12 Q 66    (17)
where Q ij are the reduced stiffness components for unidirectional lamina:

Q 11 = E 1 1 -ν 12 ν 21 , Q 12 = ν 12 E 1 1 -ν 12 ν 21 , Q 22 = E 2 1 -ν 12 ν 21 , Q 66 = G 12 . (18)
In Eq.( 18), E 11 , E 22 , G 12 and ν 12 are the lamina longitudinal, transverse and shear moduli, and Poisson's ratio respectively and their values are available in Table 1. The out-of-plane laminations parameters are defined as:

v D = H(θ) = 12 t 3 N k z 3 k -z 3 k-1 3 [cos(2θ k ), sin(2θ k ), cos(4θ k ), sin(4θ k )] dz (19)
where θ k, k∈[1,N ] denotes the k th ply orientations and z k is the coordinate of the k th ply. The laminate has a discrete set of plies with orientations [θ 1 , θ 2 , ..., θ n ]. The lamination parameters are defined as functions of the ply orientations and are very convenient for representing a laminate with a small number of variables. They are defined in a convex space and, therefore, can be used as macrospace design variables, allowing efficient gradient-based optimization.

Buckling RBDO

The objective is to maximize the bending stiffness D 11 ,

D 11 (µ v ) = U 1 + µ v D 1 U 2 + µ v D 3 U 3 , (20)
while remaining reliable with respect to the buckling phenomenon g, with:

g(v D ) = λ crit -min (λ) (21)
λ = π 2 D 11 m 4 a 4 + (D 12 + 2D 66 ) m 2 a 2 n 2 b 2 + D 22 n 4 b 4 m 2 a 2 N x ,
where λ crit is the buckling limit criterion, N x is the compressive load, a and b are the plate dimensions, D is the bending stiffness matrix, m and n are the number of half-wavelengths in the x and y directions.

The orientations variance σ 2

Θ is set to 2°. Moreover, the lamination parameters v D 2 and v D 4 are set to zero during the optimization process to follow the orthotropic nature of the laminate. The formulation of the gradient RBDO problem is written in the following form:

min µΘ -D 11 µ v D (Θ|µΘ)
subject to:

h LP (µ v D (Θ|µΘ) ) ≤ 0 P g = P g(v D = H(Θ)) ⩾ 0 ≤ P max (22)
where µ Θ are the mean ply orientation design variables, v D are the macroscopic lamination parameters with their mean µ v D , h LP is the compatibility constraint defined by [START_REF] Miki | Optimum design of laminated composite plates using lamination parameters[END_REF] for an orthotropic laminate, H is the mapping between ply orientations and lamination parameters found in Eq.( 19) and P max = 1% is the maximum failure probability. The failure probability P is approximated via the Monte Carlo method with a sample size n M C = 200000.

The normalized objective D 11 , the constraints h LP and g are illustrated in Fig. 4, with the obtained deterministic optimization solution in the lamination parameters space, represented by a small black star.

Fourier Chaos Expansion for lamination parameters

As explained in section 2.2, an inverse problem must be solved at each iteration, which is accelerated thanks to a surrogate model, to retrieve a stacking sequence. With the trigonometric nature of the lamination parameters, the orthonormal basis used for the surrogate model is different from the classical one using Hermite polynomials.

In this work, we construct a surrogate model of lamination parameters based on trigonometric polynomials that we name the Fourier Chaos Expansion (FCE). This type of decomposition was sketched out in the literature where it was shown that it is more efficient to quantify uncertainties of such form than classical expansions, such as the ones relying on Hermite or Legendre polynomials (depending on the type of random variables) [START_REF] Ko | Mixing layer growth response to inflow forcing with random phase shift[END_REF]. In this work, we deepen the FCE formalism and computation, such that it then becomes ideally suited to LPs since analytical formulation of the expansion coefficients e i in Eq.( 23) is possible and detailed in a further appendix. The lamination parameters H are expanded into a series of polynomials: where ψ i are Fourier orthogonal polynomials which depend on the standard deviation σ Θ and e i are the expansion coefficients.

Ĥ(X) = ∞ i e i ψ i (σ Θ X) (23) Basis Fourier chaos polynomials number ψ n i ( X = σ Θ X) 0 1 1 Z 11 sin( X) 2 Z 21 cos( X) -Z 22 3 Z 31 sin(2 X) -Z 32 sin( X) 4 Z 41 cos(2 X) -Z 42 cos( X) -Z 43 5 Z 51 sin(3 X) -Z 52 sin(2 X) -Z 53 sin( X) 6 Z 61 cos(3 X) -Z 62 cos(2 X) -Z 63 cos( X) -Z 64
The orthonormal properties of the FCE are analogous to the standard PCE basis. Nevertheless, the bases to be orthogonalized are not classical polynomials basis, e.g., u i = 1, x, x 2 , ...; instead, the development begins with the Fourier basis, that is, the set u i = {1, sin(nx), cos(nx)}, where n = 1, ..., ∞. An FCE representation was rapidly sketched in [START_REF] Millman | Airfoil pitch-and-plunge bifurcation behavior with fourier chaos expansions[END_REF] to obtain the probability distribution of an airfoil pitch angle where oscillatory motion was involved. In our case, we have formalized the approach to extend it to an arbitrary order of accuracy and extended it to random variables with different wavenumbers. The trigonometric polynomials are functions of scaled random variables ψ i ( X = σ Θ X), and are numerically constructed thanks to the Gram-Schmidt orthogonalization method [Y. K. [START_REF] Wong | An application of orthogonalization process to the theory of least squares[END_REF]. This method has, for instance, been used in [START_REF] Navarro | Polynomial chaos expansion for general multivariate distributions with correlated variables[END_REF] to construct arbitrary polynomial chaos and applied in aeroelasticity for the critical flutter velocity uncertainty quantification [START_REF] Nitschke | Influence of stochastic perturbations of composite laminate layups on the aeroelastic flutter of a cantilevered plate wing[END_REF].

In order to orthogonalize the polynomials with respect to the distribution π Θ , the Gram-Schmidt algorithm calculates the coefficients of the polynomials using the inner product to ensure that each polynomial is orthogonal to all of its predecessors:

ψ 0 ( X) = 1 ψ i ( X) = u i ( X) - i-1 k=0 C ik ψ k ( X), ( 24
)
where u i are the set of Fourier polynomials in (cos(n X), sin(n X)) and deterministic quantities C ik must be computed as:

C ik = E u i ( X)ψ k ( X) E ψ k ( X)ψ k ( X) (25)
The generic form of the obtained orthonormal polynomials are shown in Table 2. The details about the computational approach to efficiently evaluate the C ik coefficients can be found in Appendix ??.

With this surrogate, the statistics of the lamination parameters are readily available and very accurate. In the following, we choose that their probability density function π v are very well modeled by correlated Gaussian laws, which is useful for the failure probability gradient. This is a reasonable assumption; indeed [START_REF] Kriegesmann | Closed-form probabilistic analysis of lamination parameters for composite structures[END_REF] discusses that the distribution of lamination parameters tends asymptotically to Gaussian, increasing the number of plies in the stacking sequence.

Results

In our approach, a gradient-based optimizer is deployed in the macroscopic space in order to take advantage of the convergence properties. Indeed, it is more efficient to take advantage of the homogenized space because analytic gradients are available for the objective D 11 and the deterministic constraint h LP . Concerning the computation of the failure probability gradient, it was explained in Section 2.3. For the inverse problem step, a metaheuristic algorithm is The gradient-based MMA method [START_REF] Svanberg | The method of moving asymptotes: a new method for structural optimization[END_REF] is used, with the non-monotonic approximation of the GCMMA [START_REF] Svanberg | A class of globally convergent optimization methods based on conservative convex separable approximations[END_REF], in the macroscopic space, without making an exhaustive comparison with other potential optimization algorithms candidates. Table 3 presents the MMA parameters used in this work. Regarding the inverse problem, the genetic algorithm optimizer developed by [START_REF] Vicente | Stacking sequence retrieval of large composite structures in bi-step optimization strategies using mechanical constraints[END_REF], which followed the formulation devised by [START_REF] Irisarri | Optimal design of laminated composite structures with ply drops using stacking sequence tables[END_REF] known as SST, is used to retrieve the stacking sequence at each iteration. The parameters used for the genetic algorithm are shown in Table 4. The optimization is run through Algorithm 1, which follows the strategy detailed in Fig. (2). i) , i = 1, . . . , nMC 6: P [0] := P g V [0] > 0 7: Optimize := True 8: while Optimize do 9: j := j+1 10:

Algorithm 1 Multi-scale RBDO 1: maxit = 50, nMC = 200000, j := 0 2: Θ [0] := µΘ [0] + σ Θ X, X ∼ N (0, 1) 3: µ [0] v = µ H(Θ [0] ) 4: f [0] := D11(µ [0] v ), ∇f [0] := ∇D11(µ [0] v ) 5: V [0] := v (i) = Ĥ Θ [0],(
µ [j] v := GradientBasedOptimizer f [j-1] ∇f [j-1] , P [j-1] g , ∇P [j-1] g 11: µΘ [j] := InverseProblem µ [j] v 12: Θ [j] := µΘ [j] + σΘX 13: V [j] := v (i) = Ĥ Θ [j],(i) i = 1, . . . , nMC 14
: P [j] := P g V [j] > 0 15: µ v * , σv * := FCE Θ [j]
16:

f [j] := D11(µ v * ), ∇f [j] := ∇D11(µ v * ) 17: πV ∼ N (µ v * , σv * ) 18: ∇P [j] := ScoreFunction πV , V [j]
19: Optimize := (j < maxit) and ∥µ [j] v -µ

[j-1] v ∥ > ϵ or P [j]
g > P max 20: end while The first results were obtained with the hypothesis that the lamination parameters followed Gaussian distributions modeled at each iteration, thanks to the accurate statistics provided by the FCE. With this setup, the score function approach is performed for the gradient failure probability computation at each iteration. In Fig. 5, the optimization path (5a) and the close-up (5b) around the final design (cyan star) are shown with the shaded area corresponding to the failure domain. The green points show the designs proposed by the MMA algorithm during the iterative process, and the blue ones are the design points retrieved from the inverse problem solutions. The probability density function around the final design point is shown in the close-up. The final RBDO solution is more reliable and at a different location than the deterministic optimal.

The optimization convergence is shown in Fig. 5c with the objective D 11 and the probability plotted through the 19 iterations. Fig. 5d shows the design evolution in the lamination parameter and orientation spaces. The convergence is relatively slow toward the end, but that is a feature of the MMA algorithm. We notice that while the maximum probability constraint is not violated by the final design, the obtained value is lower than the threshold. This is, in fact, due to the discrete nature of the ply orientations. Indeed, although the MMA optimizer is designed to propose a solution closer to the probability threshold, in practice, the inverse problem resolution nudges the solution a bit, inducing slight deviations (cf. pink segments in Fig. 5b). Therefore, the final result at convergence could remain overly conservative due to the nature of the application.

Additionally, the buckling PDFs corresponding to each RBDO and deterministic optimized designs are shown in Fig. 6 for comparison. The PDF for the deterministic case is obtained by uncertainty quantification around the stacking sequence design. The deterministic design leads to poor reliability.

Discussion

In the following, we first test the robustness of the proposed numerical approach. In particular, we are interested in the sensitivity of the method to -the choice of the initial guess (to determine whether or not there is a correlation between the initial design point and the convergence of the optimization) and -the choice of different methods to evaluate the failure probability gradients. We then compare its performance with more classical optimization methods.

Study of the impact of the initial design

Optimizations were made with different initial points in the design space to check the convergence. In Fig. 7a, the initial points are spread over the macroscopic design space. Some of these initializations do not converge to a reliable solution. These initializations are mainly close to the limit-state function or in the failure domain. In Fig. 7b, the optimization paths are shown for the not converged case, and the initial point is the one surrounded by a square in Fig. 7a.

The Gaussian hypothesis for the lamination parameters can lead to an error in the probability gradients. Moreover, when the optimizer goes too deep into the failure zone, he could have difficulties returning to the safety domain since the sensitivity is close to zero. A first solution to prevent this problem is to do a multi-start optimization, which means running 3 or 4 optimizations with different initialization, mainly in the safety domain. In order to verify that the non-convergence of the optimization for some initial conditions might be due to the sensitivity of the probability in the failure domain, another approach is tested for which a deterministic constraint relative to the buckling is added. This is done without accounting for the uncertainty associated with this new component in the optimization formulation as defined in Eq.( 26) in order to control the optimization. subject to:

   h LP (µ v D (Θ|µΘ) ) ≤ 0 g(µ v D ) ≤ 0 P g = P g(v D = H(Θ)) ⩾ 0 ≤ P max (26)
The same convergence study has been done with this formulation. In Fig. 8a, the RBDOs converge to a reliable design except for one initialization. Therefore, the potential to add the deterministic constraint in a gradient-based method is shown for this case.

Gradient probability comparison

In the reference optimization in Section 3, the Gaussian hypothesis for the lamination parameters is made, which can lead to an error in the failure probability gradient evaluation. For this reason, a study is made to compare the computation of failure probability gradients with the score function approach and the centered finite differences for different levels of Gaussianity. For each Gaussianity level, the same lamination parameters and orientations are chosen for both sensitivity approaches. For the score function approach, the procedure explained in Section 2.3 is followed.

Concerning the centered finite differences, at each perturbation made around the lamination parameters, a stacking sequence is retrieved around the original stacking sequence set. All sensitivities computations were repeated twenty times with a sample size of 200000 to compute the mean and standard deviation for each sensitivity approach. In Fig. 10, the score function approach is validated for a set of lamination parameters well modeled by a Gaussian law. We emphasize that the minimal CFD step size is imposed in relation to the Monte Carlo estimator error for a given sample. In Fig. 11, the chosen lamination parameter sample does not follow a Gaussian trend. In this case, the Gaussian hypothesis of the score function is misleading, and there is a huge discrepancy for the v D 3 direction. More generally, we also notice that the score function approach is quite sensitive along the v D 1 direction, for both cases, with a higher standard deviation than the centered finite differences method. This is most likely because, in this example, the limit-state function of the constraint is parallel to the v D 1 direction, and the gradient value is close to null. Consequently, the optimization of Eq.( 22) is run with the same initialization. However, the failure probability gradients are now computed with centered finite differences without assuming Gaussianity. The optimization path is different (see Fig. 9a); however, the convergence is quite similar at the end with almost the same number of iterations and the same lamination parameters, as shown in Table 5. Nevertheless, the stacking sequences may be different because of the non-unicity of the inverse problem. In any case, their mechanical response is very similar since they have the same lamination parameter values. This study can lead to a better approach of the multi-scale RBDO. Instead of choosing one of the two approaches for the failure probability sensitivity, a hybrid approach is proposed. At each iteration, a Henze-Zirkler Multivariate Normality Test [START_REF] Henze | A class of invariant consistent tests for multivariate normality[END_REF]] is computed to determine whether the macroscopic sample can be approximated with a Gaussian law or not to know if the score function approach or the finite CFD SF

Figure 10: Comparison of the failure probability gradients obtained from score function (SF) or centered finite differences (CFD) approach. The lamination parameters sample (not represented) seems to follow a gaussian distribution in this case. On the left figure: the yellow dashed contour represents the fitted Gaussian law, and the blue contour is the non-parametric law obtained from a kernel-density estimate using Gaussian kernels.

differences are used for the sensitivity. With this approach, the optimization converges with fewer iterations and almost the same design of lamination parameters. In this case, the score function approach is used for half of the iteration.

Comparison with standard evolutionary optimization algorithm

The proposed method is compared to a more standard one often used in this context. Indeed, various multi-scale approaches rely on direct methods, such as evolutionary algorithms, at the mesoscale design space. In the following, a genetic algorithm is used as the reference approach to solve the optimization problem in the mesoscopic space (i.e., the orientation space in this case). Therefore, the macroscopic space is not exploited.

Results for the new multi-scale RBDO approach combined with different methods for the failure probability sensitivity are now compared with the result obtained by the genetic algorithm. Details of the results are found in Table 5. First, the final designs in the lamination parameters space are similar for all the methods. However, in terms of computational time, it is quite different. The reference optimization with the genetic algorithm took three times longer than the proposed method with the hybrid approach for failure probability sensitivity. Overall, the multi-scale approach, exploiting the lamination parameters, is much more efficient regarding the computational time. Nevertheless, if the macroscopic parameters (i.e., the lamination parameters in this case) can not be modeled with a parametric distribution (i.e., the Gaussian one), using only finite differences may take as much computational time as the genetic algorithm if the design variable dimension increases. Step size . We notice that the shape places much more weight under the unfeasible region.

Method

In this case, failure probability gradients differ depending on the approaches.

Conclusion

This paper proposes a new RBDO approach for multi-scale optimization exploiting macroscopic design space and a lower scale space with relevant uncertain design variables. It is an iterative approach, repeatedly focusing from one scale to the other during the optimization process in order to take advantage of each space attributes. The smoother macroscopic design space with regular gradient information allows fast gradient-based optimization, while at the lower scale, the uncertain design variables are modeled and upscaled. An inverse problem is solved at each iteration to identify a corresponding uncertain mesoscale design that matches the mean macroscopic design and gives the capability to perform efficient reliability analysis. The whole approach shows some efficiency in terms of computation time.

The strenght of this approach is demonstrated for composite laminate optimization for which the ply orientation design variables are considered uncertain. Moreover, a Fourier chaos expansion has been proposed to efficiently and accurately propagate the effect of ply orientations uncertainties from the mesoscopic (orientations) to the macroscopic (lamination parameters) space, which is used for the inverse problem and to model the joint probability density function of lamination parameters. The multi-scale RBDO approach has been applied to a buckling optimization of a composite plate. Since the optimizer needs the failure probability sensitivity, two existing methods have been tested: the score function approach and the centered finite differences. The score function approach is computationally more efficient but requires a parametric distribution of the random variables. Therefore, a hybrid approach has also been proposed where a statistical test of lamination parameters is performed to know if the score function or centered finite differences are used for the failure probability gradient. The proposed RBDO method combined with the different approaches for the failure probability sensitivity is compared with the reference composite RBDO using a genetic algorithm. Regarding the computational time, the proposed RBDO method looks better for this application. Nevertheless, a critical point concerns the initial guess. If the latter is in the failure domain or the optimizer leads the design to the failure domain during the optimization, it could have some issues converging well.

Future works will concern the application of the developed methodology to aeroelastic optimization with flutter instability. To do so, a surrogate model strategy must be implemented to compute the failure probability within a reasonable time. Future developments will include the parametrization of the composite stack thickness (i.e., number of plies) and spatially non-uniform distributions of stiffness in the structure.

Figure 2 :

 2 Figure 2: Flowchart of the proposed methodology for multi-scale RBDO.

Figure 3 :

 3 Figure 3: Geometry of the wing.

  Response of the normalized objective D11 with the deterministic design solution. The unfeasible region has been shaded to reveal the boundary between failure and safety domains.

Figure 4 :

 4 Figure 4: Objective function and constraints in the 2D lamination parameters macroscopic design space.

  Zoom around the final design. The joint PDF of the uncertainties associated with the final design is plotted with a blue color-map contour. Convergence of the objective f , failure probability P.

  convergence of the design variables represented in both design spaces.

Figure 5 :

 5 Figure 5: Optimization path with the probability sensitivity computed with the score function approach.

Figure 6 :

 6 Figure 6: Buckling constraint PDFs for deterministic and reliability-based optimized designs.

Figure 7 :

 7 Figure 7: Convergence study with different initializations.

  Sensitivity study to initial design guess: regular lattice in the macroscopic space + the reference optimization. Detailed results are presented for the square symbol case. Optimization path of a initial design from Fig.8a.

Figure 8 :

 8 Figure 8: Convergence study with different initializations with the addition of the deterministic constraint.

  2D lamination parameters macroscopic space.

  Convergence of the objective f , failure probability P and the design variables.

Figure 9 :

 9 Figure 9: Optimization path with the probability sensitivity computed with the centered finite differences.

Figure 11 :

 11 Figure11: Same caption as the previous figure but this time, the lamination parameters sample clearly does not follow a Gaussian distribution (cf. left figure). We notice that the shape places much more weight under the unfeasible region. In this case, failure probability gradients differ depending on the approaches.

Table 2 :

 2 Example of generic orthonormal Fourier basis first terms.

Table 4 :

 4 Parameters of the genetic algorithm.

	Parameters Value
	albefa	0.3
	move	0.5
	asyinit	0.25
	asydecr	1.1
	asyincr	0.5
	Table 3: MMA parameters.
	Parameters		Value
	Initial population size		160
	Population size per generation 40
	Probability of crossover	0.75
	Probability of mutation		0.5
	used.

Table 5 :

 5 Comparison of the deterministic and the various RBDO approaches: the proposed multi-scale RBDOs with various ways of computing the failure probability gradients, a direct RBDO with genetic algorithm and the deterministic optimization followed by uncertainty quantification at the final design.

								1.50
			Unfeasible region										CFD
	0.32							1.25						SF
								1.00
	0.31							0.75
	v3					Pf	v1	0.50
	0.30
								0.25
	0.29							0.00
								0.25
	0.28	0.48	0.49	0.50	0.51	0.52		0.000 0.50	0.005	0.010	0.015	0.020	0.025	0.030	0.000	0.005	0.010	0.015	0.020	0.025	0.030
				v1							Step size

Capital symbols emphasize the random nature of some of the components.

In this section, we will keep this operator simple for the sake of clarity, but as we will see in coming sections, it can take quite complex nonlinear forms.

We note that other formulations involving higher-order moments or distributions could be developed in this framework, depending on various physical hypotheses.