Léon, Lucet, Emilia, Mickaelle Sully

Pr. LEMAIRE Arnaud

Pr Maurice

Thierry Tison

Jérémy Lebon

Guénhaël Le Quilliec

Piotr Breitkopf

Rajan Filomeno Coelho

Pierre Villon

Filomeno Coelho

A two-pronged approach for springback variability assessment using sparse polynomial chaos expansion and multi-level simulations

Keywords: Springback variability assessment, Sparse polynomial chaos expansion, Semi-analytical bending-under-tension model Uncertainty Quantification, Model Sensitivity, Springback Variability Assessment, sensitivity-constrained Design of Experiment, Sparse Polynomial Chaos Expansion

Now comes the time to write the last page of this thesis work, bringing by the way my student journey to an end. Here are a few words full of gratitude to those who contributed to make this journey a success, so rich in lessons, in both scientific and personal

1.1. Introduction

Introduction

Numerical models are extensively used in science disciplines such as physics, chemistry, computational engineering, biological and social sciences, etc. These are aimed at providing a better understanding of real world phenomena. In mechanical engineering, they are devoted to play a central role in the analysis and design of structures and processes. They provide an invaluable tool to the engineer to the search for the best performing structures and processes. In most of real-life applications, the performances are assessed simultaneously by multiple criteria (cost, environmental impact, safety, robustness). It has been now commonly admitted that uncertainties are inherent of this process and may have a non-negligible influence on the design performances.

In this thesis work we explore two issues the engineer has to face: the quantification of the impact of uncertainties on the design (Uncertainty quantification), and its incorporation into an optimization process with multiple criteria (Multiobjective Optimization under Uncertainty).

Uncertainty Quantification (UQ) Along the elaboration process of the numerical models, approximations are made. They result from the implementation into a computer code of an imperfect abstraction (the mathematical model) of the reality (the real experiment). Figure 1.1 proposes to identify three sources of approximations involved in the construction of a numerical model. One of them is issued from the elaboration of the mathematical model of real experiment. It describes the underlying physics of the observed phenomenon. The induced errors are due to simplifications or lack of knowledge that the mathematical abstraction of the physical system may suffer from (simplified physics error) but also to the choice of the underlying mathematical formulation (Galerkin formulation instead of strong FEM formulation, etc). The implementation of the mathematical model into a computer code also generates its own type of errors. It encompasses the necessary approximation/discretization errors as well as the numerical and round-off errors, and programming errors. Finally, a third source of errors may be emphasized by repeatedly performing the real experiment under the same configuration:

1.1. Introduction different results may though be observed as a result of errors in measurements, of the evolution in time and space of material properties, the environmental conditions, etc. All these errors laid end to end, may lead simulation results to represent the underlying physics with a low level of fidelity. The question of the accuracy of the simulation results must be inevitably confronted. This topic is addressed in the literature by the Verification1 and Validation domain2 (V&V) [OTH04] as well as Uncertainties Quantification (UQ). This thesis work is concerned with the latter domain where the errors are preferentially referred to "uncertainties". Different theories to represent uncertainties are provided in the literature (fuzzy logic, confidence theory, probabilistic,... ). This thesis focuses on the probabilistic viewpoint: uncertainty sources are described by random variables characterized by a probability density function 3 .

UQ provides the engineer with valuable tools to assess the impact of the designed process design and sensitivity analysis, reliability based design and variability analysis:

• Robust design and Sensitivity analysis are aimed at quantifying the contribution of each random input variable to the variability of the response [SAA + 10]. One may distinguish between two types of approaches addressing the sensitivity analysis from two different points of view. The local sensitivity analysis focuses on small variations of the input parameters around a nominal value and studies how the model response is affected; the global sensitivity analysis focuses on quantifying the output uncertainty due to changes of input parameters (taken either separately or in combination) in their entire range of variation.

• The reliability of a system is defined as the assessment of its failure probability, that is the nonsatisfaction of its expected performance [LCM09]. As a system is usually composed of subsystems, the failure of the system may be caused by different scenarii caused by the failure of one or more components to be identified.

• Finally, the variability study aims at completely characterizing the output probability density function.

Achieving these tasks firstly requires the definition of the probabilistic model. It consists in identifying the uncertain input data and to model them according to their respective (joint) probability density function (Identification). Then, to assess their impact onto the model output, one has to transform the modeled uncertainties on the input variables to output statistical (e.g. statistical moments) or probabilistic quantities of interest (e.g. probability density function). This process forms the propagation of uncertainties. Two approaches may then be distinguished: intrusive methods consist in adapting the governing equations of the deterministic model to uncertainty propagation; and non-intrusive methods, which make use of a series of calls to the deterministic model to propagate the uncertainties. In this thesis, we consider the non-intrusive uncertainty propagation scheme enabling the use of complex in-house or commercial simulation codes. In this context, the UQ process is illustrated in Fig1.2.

Introduction

Figure 1.2: Steps for non-intrusive uncertainty quantification

Once the stochastic model has been identified, a Design Of Experiment (DoE) is built on the stochastic variables with respect to their probability density function. For each point in the DoE, deterministic calls to the numerical model are performed. A post-treatment phase allows to retrieve the statistical (histogram, mean, variance, etc.) and/or the probabilistic (quantiles, probability of failure, etc.) data of interest. The uncertainty propagation phase is computationally expensive as a single call to the numerical model may already be costly and numerous calls may be necessary to assess the stochastic and/or probabilistic quantities of interest.

Optimization Under Uncertainty (OUU) The second issue addressed in this thesis concerns the multiplicity of pontentially competing criteria required to assess the performances of a structure.

The incorporation of uncertainties [START_REF] Eldred | Formulations for surrogatebased optimization under uncertainty[END_REF] transposes the concepts of robustness, and reliability into the optimization context.

• Robust design optimization is aimed at identifying the most performing responses with a particular insight on their sensitivity to random perturbation of design variables or parameters.

• Reliability-based design optimization incorporates the failure probability as a constraint of the optimization problem. The assessment of robustness and reliability measure usually require the computation of statistical data (mean/variance) and or probabilistic data such quantiles and probability of failure. The robustness and the reliability enters the optimization problems either as objective of constraints and may then be assessed each time a new solution is provided by the optimizer. The optimization and uncertainty quantification loops may then be nested which leads to high computational costs. Strategies to overcome this limitation have ben developed. One way to circumvent it consists in multifidelity approaches whose numerical challenges are addressed in more details in the next section.

Numerical challenges

As in practice both OUU and UQ are performed by repeated calls to the (deterministic) numerical model, the computational expense may rapidly become unaffordable. Multifidelity approaches address this issue by combining in a hierarchical manner two levels of numerical models: a "high-fidelity" model which is characterized by a high accuracy but also a high computational cost is combined with "low-fidelity" model less accurate but also less computationally costly.

On the development of metamodel in multifidelity approaches

Categorization of surrogate models Non-intrusive schemes based on surrogate-based approaches have been proposed in optimization and uncertainty propagation domains:

• on the UQ side, the Chapter 2 provides an overview of computational methods and shows that the Polynomial Chaos Expansion (PCE) has been widely extended and adapted to deal with computationally costly and high-dimensional black-box functions;

• on the optimization side, optimization strategies for computationally costly and high-dimensional black-box simulations [SW10] as well as metamodeling techniques for optimization [WS07] are proposed.

In the following the original computational model will be referred as the "high-fidelity" model, and the surrogate model as the "low-fidelity" model as it is less accurate than the original numerical model.

Numerical challenges

One may classify these approaches as "multi-fidelity" or "variable-fidelity" when physics-based and non-physics surrogates are considered:

• physics-based surrogates built from partial convergence of the "high -fidelity" model [BBM + 00], those based on simplified physics [START_REF] De | Multivariate modelling of variability in sheet metal forming[END_REF], those based on a coarse discretization [ALG + 99, VHS02, SLZ + 11];

• interpolating or regression-based response surface from a Design of Experiment (DoE) evaluated with "high-fidelity" simulations. When the derivatives are available, one may build a first order consistent approximation [WG95]. When only the function evaluation is available, general approaches, such as general least square, polynomial response surface [START_REF] Myers | Response surface methodology: process and product optimization using designed experiments[END_REF], kriging [Ste99, Cre90], moving least squares [START_REF] Lancaster | Surfaces generated by moving least squares methods[END_REF]NTV92], radial basis function [B + 00] may be used to build a global approximation from the available set of points.

The non-physics based approaches do not perform well in high dimensions ("curse of dimensionality").

Their accuracy depends on the underlying "high-fidelity" model and may be sensitive to the type and to the size of the DoE used. Moreover, the computational expense related to the training phase depends on the cost of the underlying "high-fidelity model". Finally tuning such metamodel is not straightforward (e.g. excessive smoothing, overfiting). Physics-based approaches do not suffer from these issues but have to be specifically built for one application, thus may not be always available.

In the UQ domain, during the last two decades, the Polynomial Chaos Expansion (PCE) has been widely used to propagate uncertainties (Chapter 2). Broadly speaking, the PCE [Gha91] is a metamodel that is intended to give an approximation of the stochastic behavior of a function y (scalar random process). When a single random variable ξ is considered, the P th -order PCE of y is defined as:

y = P -1 i=0 γ i Ψ i (ξ) (1.1)
where ξ S are standardized random variables, and γ = {γ} P -1 i=0 the set of PCE coefficients which have to be computed during the training phase.

On the validation of surrogate based approaches Twofold challenge has to be tackled:

• ensuring the convergence of the multifidelity surrogates,

• ensuring the accuracy/stability of the underlying "high-fidelity" model.

Convergence of hierarchical approaches

In the optimization field, the key point is to guarantee that the surrogate-based optimization process converges to an optimum of the original model. In the UQ field, the key point is to guarantee that the surrogate-based UQ approach converges in terms of statistical measures of the system responses (e.g. the output variances, sensitivities, robustness measures, probability measures)

In the optimization field, rigorous approaches combining in a hierarchical way a "low-" and "high-'' fidelity model are referred to surrogate-based optimization and/or model management framework [ALG + 01].

The premise of this methodology is to claim that a surrogate based on a physics-based low-fidelity model 1.2. Numerical challenges and an interpolant of the discrepancy may provide a more cost-effective approximation of the high fidelity model. A proof of convergence is ensured provided that one employs a sufficiently rigorous verification (e.g. trust region methodology) and that the surrogate model satisfies first-order [AL01] and sometimes even second order [EGC + 04] consistency conditions with the underlying high-fidelity model. In the trust region methodology, the "low-fidelity" model values are corrected to fit the "high-fidelity" model values (up to the second order if necessary) using either additive, multiplicative, or combined multiplicativeadditive correction functions on a set of collocation points. In addition, to guarantee that the progress made with the "lower-fidelity" model also leads to an improvement with the "higher-fidelity" model, the "lower-fidelity" model is regularly updated using systematical calls to the "high-fidelity" model.

In the UQ field, the use of non-physics based metamodel is usually not discussed or is highlighted a posteriori by assessing error measures on the statistical or probability data of interest. Moreover, in some cases, adaptive sampling is used in order to update the metamodels and control its accuracy in a particular space region [DSB11b, WCS + 13]. It is only recently that [START_REF] Ng | Multifidelity uncertainty quantification using nonintrusive polynomial chaos and stochastic collocation[END_REF] investigates the extension of the classic multi-fidelity optimization concepts to the uncertainty quantification fields. The "low-fidelity" model values are corrected using combined additive-multiplicative correction function to match the highfidelity values on some collocation points, and the non-intrusive PCE is trained on the low-fidelity model to propagate the uncertainty.

On the accuracy/stability of the "high-fidelity model". In the optimization field, ensuring the accuracy and/or the stability of the "high-fidelity" model and of the gradients in the whole domain of interest belongs to the "good practice rules" (e.g. the minimum step size of the finite difference scheme is assessed a priori). In UQ, these instabilities may be treated as epistemic uncertainties through diverse approaches such as interval analysis, possibility theory, evidence theory or probability theory [START_REF] Helton | Representation of analysis results involving aleatory and epistemic uncertainty[END_REF][START_REF] Helton | Latin hypercube sampling and the propagation of uncertainty in analyses of complex systems[END_REF].

When considering the high-fidelity model as a black-box, no direct mapping is provided from the input variables ξ = [ξ 1 , ξ 2 , . . . , ξ M ] to the derivatives of a selected output function y, except in seldom cases where the black-box itself provides it. Thus a numerical scheme has to be used in order to compute the sensitivities (i.e the derivatives) of the model. In this thesis, we are interested in the model output stability with respect to the variation of the variable ξ i , i ∈ {1, . . . , M} assessed by the following non dimensional finite difference scheme :

µ i = ∆y(ξ i ) ∆ξ i × ξ i nom y(ξ nom i ) (1.2)
where ξ nom i , i ∈ {1, . . . , M} is a nominal value.

∆y(ξ i ) = y ξ nom i + ∆ξ i 2 -y ξ nom i - ∆ξ i 2 . (1.3)
When decreasing the order of magnitude of the perturbation (-log(∆ξ) increasing), the non-dimensional sensitivity µ i computed for the "high-fidelity" model exhibits subsequent behaviors illustrated in a qualitative manner 4 on Fig. 1.4:

1.2. Numerical challenges 1. firstly, for "large" variation of ∆ξ i , the variation µ i reveals the non-linear behavior of the model.

No brutal variation of µ i is observed and the model is considered as trustworthy.

2. secondly, µ i stabilizes around a constant value μi where the model may be considered as linear.

3. thirdly, on reaching the threshold , µ i becomes unstable.

4. finally, the threshold shows the model sensitivity limit: for this range of variation, the model is not sensitive anymore. In both optimization and UQ domains, comparing the model error with the variation magnitude of the output function is of paramount importance to ensure the validity of the approach. To do so, one may arbitrarily define a resolution threshold δ * i on ξ i determined from an acceptable value of µ i arbitrarily defined. δ * i is then defined as.

δ * i = argmin ∆ξi µ i (∆ξ i ) > µ * i , i ∈ {1, . . . , M}. (1.4)
When evaluating two sample points ξ

(1) i and ξ

(2) i one considers model response trustworthy when ∆ξ =

|ξ (1) i -ξ (2)
i | ≥ δ * i and erroneous otherwise. In this work, we focus on the observation that the resolution threshold and the variation range of the training data may be close. The number of achievable simulations by the "high-fidelity" model is not only limited by the cost of an individual simulation but also by its intrinsic resolution. Thus one has to face two issues: firstly, an insufficient number of "high-fidelity" simulations may harm the accuracy of the response surface, secondly, a too high number of simulations may introduce numerical noise which also directly leads to an inaccurate surrogate model.

MultiObjective Optimization under Uncertainty (MOOU)

In real-life applications the solution of the optimization problem involving possibly competing criteria is not unique: a set of best compromise solutions (called Pareto set) is sought rather than a single solution.

Contributions of the thesis

The Pareto set P * is then defined as the set of all non-dominated solutions in the design variable space Ω:

P * = {x * ∈ Ω | x ∈ Ω; y(x) ≽ y(x *
)} and the Pareto front defines its representation in the objective function space. The symbol "≽" defines the dominance concept between the vectors y(x) and y(x * ) such as ∀i ∈ {1, ..., m}, y i (x) ≤ y i (x * ) and ∃i ∈ {1, ..., m} | y i (x) < y i (x * ). The Moreover, in most of the studies (e.g. [YCL + 11]), the influence of uncertainties is rarely considered using a metamodel approach. The formulation of the optimization problem has to be adapted and the notion of robustness and reliability has to be defined in conjunction with the dominance notion. Thus, the development of a non-intrusive stochastic metamodel strategy, specially when a complete probabilistic description might be useful to compute the probability of failure (e.g. in Reliability Based Optimization) is a challenged we address in this thesis work.

Contributions of the thesis

The contribution of this thesis is twofold.

• The issue of Uncertainty Quantification with high fidelity model characterized by inadequate resolution for uncertainty propagation is illustrated in Chapter 3 and in Chapter 4. It highlights that a careful attention has to be given to the validation of the high-fidelity model to perform uncertainty quantification. In fact the variation range of the stochastic variables may be smaller than the resolution to the model.

• The incorporation of uncertainty into a multi-objective optimization in a non intrusive framework is addressed in Chapter 5, based on the hierarchical combination of two metamodels: one for 1.3. Contributions of the thesis the deterministic variation of design variables and the other to assess the variability of the model response for each set of design variables.

In the following sections we discuss both issues.

A multifidelity approach for springback variability assessment

In this section we illustrate the first contribution of this thesis using an example of deep drawing of metal parts. The deep drawing process is a manufacturing process which aims at permanently changing the shape of a metal sheet through the action of a moving punch forming the metal against a motionless die (Fig. 1.6(a)). When the loading is removed, additional deformations appear and the so called springback phenomenon occurs (Fig. 1.6(b)). To perform a variability study of the post springback shape, we propose a two-pronged approach .

A first ingredient consists in defining a physics-based "low-fidelity" model with a finer resolution to replace the "high-fidelity" model for small variation of parameters. A custom sparse stochastic surrogate (Sparse PCE, where only the most relevant term of the PCE are retained) is then used to perform the variability study onto the newly "low-fidelity" "high-resolution" model.

The simulation of the springback process provides a noisy numerical behavior against small thickness variations. The FEM model may be considered as trustworthy for a resolution which represents a tolerance on the thickness around 1 20 of its nominal value. The resolution of the numerical model is thus clearly insufficient to perform a variability study on the tolerance range. Three aspects of the FE implementation of the deep drawing process may be responsible for these numerical instabilities: the coarsity of the mesh, the incorrect stress integration through the thickness and the contact algorithm. Improving the FE model needs (among other) to refine the discretization of the model in every direction preserving the same aspect ratio which leads to excessive computational cost. We choose to combine both models: the FE "highfidelity" model by a "lower-fidelity" physics-based metamodel with both features: a lower computational cost and a higher resolution. Such a metamodel has been recently proposed by [START_REF] Le Quilliec | Semi-analytical approach for plane strain sheet metal forming using a bending-under-tension numerical model[END_REF]. It describes the deep drawing process 1.3. Contributions of the thesis of a 2D metal sheet as a 2D plain strain Bending-Under-Tension process (B-U-T) in a semi-analytical framework combining an analytical approach with a FE model (assuming negligible shear stress).

The contact occurring during the deep drawing process is modeled here using an analytical approach inducing consequently no numerical instabilities. Up to 200 integrations points are considered. In this configuration, the impact the numerical noise on the springback shape parameters is drastically reduced.

We have trained the PCE model on the B-U-T physics-based surrogate in order to perform a variability study of the springback shape parameters. The coefficients of the PCE are computed using a regression based approach from a Design of Experiment (DoE) defined by a standard Latin-Hypercube Sampling.

We propose a stochastic model representing up to 8 independent random variables for the springback shape parameter study.

A well known limitation of the PCE lies in the computational cost to compute the full set of coefficients when the number of variables or the degree increases. We here choose to apply the Least Angle Regression Stagewise (LARS) method [BS08c] on an 7 th order PCE and the empirical error estimate as the stopping criterion.

The savings in computational cost is significant since in each case only a fourth of the coefficients are needed to reach the convergence on each of the function representative of the springback shape with a preserved accuracy. This approach demonstrates that the use of simplified physics-based model allows for accurate UQ and acceptable computational costs.This approach is not limited to 2D and opens the way to the use of other types of physics-based metamodels such as one-step or PGD/POD.

Adapting the sampling to the model resolution: Fat-Latin Hypercube sampling

To train the stochastic metamodel, a compromise has to be found in the distribution and the size of the set of training points. In fact, a too small number of simulations leads to an excessive smoothing of the response surface while a too high number of simulations may introduce numerical noise leading by the way to the same lack of accuracy of the stochastic metamodel.

In Chapter 3, the compromise is addressed. A modified LHS scheme is proposed in order to take into account the resolution of the model. It provides an upper bound on the sampling density. A coherent metamodel scheme based on the PCE is built taking into account this upper bound on sampling density, and the lower bound given by the regression approach to compute the PCE coefficients.

The main idea is to build around each sampling point a restricted area free from other samples while preserving the LHS property. The characteristic size of each area is parameterized by the a priori identified model resolution δ * i , {i = 1, . . . , M} in each of the M dimensions. Depending on the chosen norm, different shapes of the restricted area are considered. Here, an illustration is given for the L ∞ and the L 2 norms. The "Fat-LHS" scheme provides a dense sampling along with a maximum number S ub of samples respecting the model resolution. With this limited number of samplings, we train a sparse version of PCE where only the most significant terms are kept.

Contributions of the thesis

The efficiency of the Fat-LHS is illustrated on the variability assessment of springback of a 2D deep drawn metal sheet. It exhibits a higher convergence rate to the two first statistical moments than a similar size LHS (without taking into account the model resolution). On a higher dimension test case (8D), the convergence results for the different truncation strategies are compared. The sampling is then coupled with an adapted metamodel strategy in order to efficiently propagate the uncertainties at low computational costs.

Towards Multi-objective Optimization Under Uncertainty

In this section, we address the incorporation of uncertainty treatment into a multi-objective optimization process using surrogate-models in a non-intrusive framework. After a brief recall on the basics of the multi-objective optimization under uncertainty, we present the following original contributions:

• the development of a non-intrusive hierarchical stochastic metamodels based on Moving-Least Squares (MLS) and PCE;

• its application to multi-objective reliability based optimization of space truss structures.

Evolutionary population-based algorithms are designed to efficiently address multi-objective optimization problems these optimization problems due to their ability to possibly provide many of the Pareto-optimal solutions in a single algorithm iteration [START_REF] Fogel | The advantages of evolutionary computation[END_REF]And02]. We use Non Sorted Genetic Algorithm II [Deb02] for 3 reasons: low computational complexity, elitist approach and ability to preserve diversity in the population.

ξ = [ξ 1 , . . . , ξ M ] are the M independent stochastic variables described by their probability density function. x represent the set deterministic design variables. y(x, ξ) = [y 1 (x, ξ), . . . , y m (x, ξ)] is the vector of the m stochastic objective functions The challenge consists in the computation of the non-dominance probability for each set of solution produced by the algorithm. A straightforward approach would consist in building a polynomial chaos expansion of each provided solution. The statistical quantities and non-dominance probability needed are then computed exclusively using this approximation.

We propose a hierarchical metamodeling approach. It consists in a surrogate model providing an original mapping from the mixed deterministic-stochastic variable space to the response function space.

The first step consists in building a metamodel of the PCE coefficients:

x → γi (x), i = 0, . . . , P -1 (1.5)
and then to reconstruct the PCE for each design variable provided by the algorithm:

(x, ξ) → y(x, ξ) = P -1 i=0 γi (x)Ψ i (ξ) (1.6)
resulting in an increased computational efficiency.

We validate the proposed approach by comparing the obtained stochastic Pareto set with a Monte Carlo sampling under different probability levels. An explicit two-variable test case allows to analytically compute the response surface and to compare it with those obtained by an MLS approximation of the PCE coefficients. A perfect agreement is observed. A second test case deals with the multi-objective sizing optimization of three 3D linear elastic truss structures. In most of the situations the combination of second order MLS interpolation of the PCE coefficients provides the best results. The influence of the 1.4. Thesis outline PCE order is investigated. In most of the cases the second order approximation is better, highlighting a possible overfiting phenomenon when the approximation order increases or a sampling size is too small.

Finally, the sizing optimization of trusses is verified using Monte Carlo Simulations.

Thesis outline

The original contributions of this thesis are organized in the following manner:

• Chapter 2 presents a literature review with a special insight firstly on general UQ methodologies and then a focus on the spectral methods for uncertainty propagation. For the latter approach the numerous challenges are put into evidence as well as the concepts proposed in the literature.

• when modifying the experiment parameters or when changing the software used to predict the springback values. However the stochastic analysis of metal forming processes requires both a high precision and low cost numerical models in order to take into account very small perturbations on inputs (physical as well as process parameters) and to allow for numerous repeated analysis in a reasonable time. Chapter 5 proposes a metamodel based Multi-Objective Optimization under uncertainty based on a hierarchical metamodel approach with respect to the deterministic and random input parameters.

Introduction and general concepts

In the past decades a copious amount of work has been devoted to the characterization of the impact of uncertainties on the behavior of a mechanical system described by numerical models. A first step in their quantification may consists in the categorization of the uncertainties. This topic is subjected to many discussions by statisticians and engineers, but two concepts are clearly distinguished and commonly admitted : epistemic uncertainties are reducible uncertainties by means of gathering more data or by refining the model; on the opposite, random or aleatory uncertainties are irreducible uncertainties, intrinsic to the components or to the environment of the real experiment. However, the classification of uncertainties into one of these two categories may depend on the context of the study but is of paramount importance. [START_REF] Kiureghian | Aleatory or epistemic? does it matter?[END_REF] highlights that a correct a priori identification of the types of uncertainties allows to identify (hence to elaborate) reasonable strategies in order to a priori reduce the uncertainties at the end of the line: a wrong categorization may influence the results from several orders of magnitude.

To quantitatively assess their influence, one may classify the methodologies into two groups: intrusive and non-intrusive methodologies. The formers require mathematical developments of the governing equations to produce, most of the time, semi-analytical solutions for the stochastic analysis. This class of methods present a major drawback: their implementation tends to become complex and analytically cumbersome in case of highly non linear physical systems. One may then prefer the use non-intrusive techniques which consider the numerical simulation process as a black-box. This section only deals with the latter class of methodologies described on a probabilistic framework.

First and foremost, a brief recall on the probabilistic description of uncertainties is provided together with a classification and a short review of the non-intrusive approaches for UQ. It essentially aims at pointing out the notations and to provide the reader with an overview of non intrusive UQ.

The second section forms the heart of this chapter. It is dedicated to the class of spectral methodologies taking their roots in functional evaluation of the stochastic data of interest. A special focus is provided on the use of the Polynomial Chaos Expansion (PCE) in a non-intrusive context. Since the groundbreaking work of [Gha91], this methodology has encountered a growing interest and has been used in many different fields. This section highlights the evolution of the Polynomial Chaos Expansion from its initial formulation as Hermite-Homogeneous Polynomial Chaos to the most recent advances. The goal of this section is not to systematically provide a deep insight into each theoretical advances but more to guide the interesting reader by providing him in an organized manner some fundamental references to the field.

Mathematical framework

Probabilistic space

The observation of a random phenomenon provides uncertain outcomes which nevertheless may follow a regular distribution for a large number of trials. The set of all possible outcomes, denoted Ω is called the "elementary sample space". The result of one trial, is a subset ω ∈ Ω called "elementary event". The set of all admissible events associated to aprobability measure p F forms a σ-algebra associated with Ω denoted by F. Finally, the triplet (Ω, F, p F ) forms the "probability space".

ξ : F → B(R) ω → ξ(ω).
(2.1)

In this thesis we assume that the uncertainties can be modeled using a finite number random variables gathered into a M -sized random vector ξ = {ξ 1 , . . . , ξ M }. We thus define the random vector ξ as a measurable function defined on the σ-algebra F with values in a Borel set of R M denoted B(R M ).

ξ : F → B(R M ) ω → ξ(ω).
(2.

2)

The dimensions of the random vector ξ are denoted using subscript, while their realizations are denoted by superscripts: ξ (j) i denotes the j th realization of the i th component of the random vector

ξ = {ξ 1 , . . . , ξ M }.
Let ξ be any of the components of ξ and τ one realization of ξ. In this thesis work, we consider that ξ is completely characterized by a continuous probability density function (pdf) p ξ (ξ) whose integral over the sample space R is called the cumulative probability density function (cdf). The mean value of the random variable ξ is given by:

µ = E[ξ] = R τ dp ξ (τ ) (2.3)
where {ξ = τ } = {ω ∈ Ω|ξ(ω) = {τ }} and dp(ξ = τ ) = p ξ (τ ) = p F ({ξ = τ }).

The variance V of the random variable is given by:

V = R (τ -µ) 2 dp ξ (τ ) (2.4)
The n-th order moment is given by

E[ξ n ] = R τ n dp ξ (τ ) (2.5)
In this thesis work, we only consider real random variables with finite second order moment (E[ξ 2 ] < ∞).

We denote the corresponding vectorial space by L 2 (B(R), T , P ξ ).

Considering now two random variables ξ = {ξ 1 , ξ 2 }, the joint probability density function is naturally denoted p ξ1,ξ2 or more concisely p ξ . The covariance between these two random variables is denoted

C(ξ) = E [(ξ 1 -µ ξ1 ) × (ξ 2 -µ ξ2 )] . (2.6)
The correlation factor ρ ξ1,ξ2 is given by:

ρ ξ1,ξ2 = C(ξ) σ ξ1 σ ξ2 .
(2.7)

Introduction and general concepts

A random field A random field y(x, ω) [Adl81] may be defined as a collection of random variables indexed by a continuous parameter x defined on a bounded set D ∈ R d with value in R (Eq. 2.8). A one-dimensional random field (D ∈ R) is usually called a stochastic process.

y :R d × Ω → R (x, ω) → y(x, ω).
(2.8)

For a chosen x 0 ∈ D, y(x 0 , .) is a random variable. One denotes by y(x 0 ) = {y(x 0 , ω 1 ), . . . , y(x 0 , ω S )} the set of S realizations of the random field. Its mean is then defined as:

µ(x 0 ) = Ω y(x 0 , ω)dp y (x 0 , ω) (2.9)
The continuity of a random field is often characterized by its mean-square continuity:

E[∥y(x, .)∥ 2 ] < ∞, ∀x ∈ D lim x→x0 E[∥y(x, .) -y(x 0 , .)∥ 2 ] = 0 (2.10)
The matrix-valued covariance function of the random field taken as the covariance function between two distinct x 1 , x 2 :

C(x 1 , x 2 ) = E [(y(x 1 , .) -µ(x 1 )) × (y(x 2 , .) -µ(x 2 ))] (2.11)
The correlation factor ρ ξ1,ξ2 is given by: ρ y(x1,.),y(x2,.) = C(x 1 , x 2 ) σ y(x1,.) σ y(x1,.)

(2.12)

Common concepts for probabilistic approaches

Building a stochastic model usually consists in two steps. The first one is concerned with the proper identification of the random inputs by a finite number M of random variables (such as parameters of the system, material properties, etc.) or processes (random time dependent loading,etc.) [START_REF] Franklin | Numerical simulation of stationary and non-stationary gaussian random processes[END_REF]. For Gaussian random variables and processes, the identification is straightforward as they are completely determined by the two first statistical moments (mean and covariance). This is however not the case in general for non Gaussian processes and their identification still remains an open research area [YS88, WC94, SG02a, PPS02, LCS07]

The second one consists in identifying the (joint) probability distribution of the random variables, fields of interest. Two approaches may be distinguished:

• A direct approach consists in directly constructing the probability function p ξ based on information theory (i.e using only the information at hand on ξ).

• An indirect approach consists in introducing a measurable mapping h between the random vector y(ξ) of unknown probability density function (p.d.f) p y to the random vector ξ whose p.d.f p ξ is known: y(ξ) = h(ξ). Then, p y (y) is the transformation of p ξ (ξ) by h. This method is advantageous each sample associated with the probability p ξ (ξ (i) ), i = 1, . . . , S such as

S i=1 p ξ (ξ (i) ) = 1.
The entropy measure is given by:

T (p ξ (ξ (1) , . . . , ξ (S) )) = -K S i p ξ (ξ (i) )ln(p ξ (ξ (i) )) (2.13)
where K is a positive constant. To solve the inference problem stated in the direct approach, the only possible unbiased assignment consists in searching the probability distribution which maximizes the entropy with regards to the information at hand. The classical optimization problem to solve may then be stated as: max

p ξ
T (p ξ (ξ (1) , . . . , ξ (S) ))

s.t E[y(ξ)] = S i=1 y(ξ (i) ) S i=1 p ξ (ξ (i) ) = 1.
(

The use of Lagrange multipliers λ 0 and λ 1 leads to the following optimal solution:

p ξ (ξ) = e -λ0-λ1ξ (2.15)
In the multi-variate case where ξ = {ξ 1 , . . . , ξ M } the maximum entropy probability distribution generalizes to [START_REF] Jaynes | Information theory and statistical mechanics[END_REF]:

p ξ (ξ) = n i=1 exp(-[λ 0 + λ 1 ξ (i) 1 + . . . + λ M ξ (i) M ) (2.16)
The set of Lagrangian coefficients, may be retrieved by solving the set of equation obtained by replacing the expression of p ξ (ξ) in the expressed constraints.

Maximum likelihood Let us consider the following problem (also referred as "forward problem") [START_REF] Millar | Maximum likelihood estimation and inference: with examples in R, SAS and ADMB[END_REF] y ≈ f (ξ) + ϵ (2.17)

where ξ = {ξ 1 , . . . , ξ M } represents here a vector of model parameters, y = {y 1 , . . . , y r } is a vector of output of interests seen as realizations of random variables, f is a (non-linear) function R M → R r and ϵ ∈ R r a random vector encompassing all possible kind of uncertainties.

Let us assume that a set of observed output is available. The goal of Maximum Likelihood estimation is to retrieve the input parameters which has most likely generated the observed sample. To do so, one has to maximize the following probability p ξ (y|ξ), which is called the Likelihood.

Introduction and general concepts

A trivial (and common) assumption consists in assuming that the random vector ϵ ∈ R r representing the components of the error produced are independent and identically distributed characterized by a normal probability density function ( E[ϵ i ] = 0, i = 1, . . . , r, the covariance matrix of the error is C = σ 2 I r×r where I r×r is the r × r identity matrix). In this case p ξ (y|ξ) ∈ N (f (ξ)) and the likelihood function becomes

L(ξ) = r i=1 p ϵi (ϵ i )) L(ξ) = r i=1 p ϵi (y i -f i (ξ)) = r i=1 exp - (y i -f i (ξ)) 2 2σ 2 = exp - r i=1 (y i -f i (ξ)) 2 2σ 2 (2.18)
Searching for the maximum of this function is equivalent to finding the minimum of r i=1 (y i -f i (ξ)) 2 which may be done by any appropriate optimization algorithm or by considering a classical minimum distance leasts square problem. In the general case ϵ is not supposed to be Gaussian, the Likelihood function may be computed using Monte Carlo Simulations which may be too expensive. Another more efficient consists in using a kernel approximation method.

Bayesian Inference This approach allows to determine the prior probability distribution of input parameters ξ = {ξ 1 , . . . , ξ M } when the only information at hand concerns the output of the model y = {y 1 (ξ), . . . , y J (ξ)}. This process is called inference. In the context of the Bayes' framework, the knowledge about the true values of the parameters before and after having observed the data are described using probability. In the general case one may use inference only in order to refine the probability distribution of the input knowing the outputs, but we present here this approach in the general case where the whole probability density function of the input parameters are searched.

The Bayes' theorem states that:

p(ξ|y) = p(y|ξ)p ξ (ξ) p(y|ξ)p ξ (ξ)dξ (2.19)
where p(ξ|y) is the posterior distribution of the parameters values (obtained after having observed the data), p(y|ξ) = L(ξ) is the likelihood function and p ξ (ξ) is the prior probability density function (obtained before having observed the data) of the parameters values and p(y|ξ)p ξ (ξ)dξ is a constant (as p(ξ|y) is a distribution on ξ, y is a constant) Thus a more concise form of the Bayes theorem may be obtained by:

p(ξ|y) ∝ (y|ξ)p ξ (ξ) (2.20)
In the general case, the described inference process consists in refining an a priori given set of probabilities (prior models) in the light of the observed outputs.

Introduction and general concepts

In the case where the input probabilities are unknown, the Bayes' inference approach methodology may often be used in conjunction with the Maximum Entropy Principle. In general, the Bayes theorem offers a robust mean for inference, however vague prior information on the parameter value the posterior remains improper, which may cause serious problems for the verification of complex models [START_REF] Geyer | Practical markov chain monte carlo[END_REF].

However, if a sufficient amount of information on the prior is available, the Bayes' inference has been successfully use in the literature to assess the confidence in model prediction by comparing the model output with experimental data.

A first direct approach in order to compute the posterior distribution consists in using Monte Carlo Sampling on the prior distribution. The sample produced is thus non optimal with regards to the posterior distribution (many sample points may located in the tail of the posterior probability) and thus the chosen set of parameters with higher posterior probability may not be relevant. The use of Markov Chain Monte Carlo (MCMC) [START_REF] Se Ahmed | Markov chain monte carlo: Stochastic simulation for bayesian inference[END_REF] based on the Metropolis-Hasting's algorithm [START_REF] Chib | Understanding the metropolis-hastings algorithm[END_REF] may allow to increase the computational efficiency. MCMC are used to directly sampling the posterior distribution, to build an optimal sample with regards to the posterior distribution. Moreover, it eliminates the need of computing the constant p(y|ξ)p ξ (ξ)dξ.

A classification of indirect approaches

Simulation based methodology Monte Carlo Samplings (MCS) is the most direct approach for uncertainty quantification. In the literature, the Monte Carlo simulation is widely used in order to assess the two first statistical moment of a multi-dimensional random variable. These quantities are extracted from S independent realizations of random inputs. These are generated from their a priori prescribed probability density function. For each realization, the data are fixed and the problem becomes deterministic: the is obtained from a set of deterministic calls to the considered model. From this set of deterministic evaluations, statistical information are extracted (mean, variance, histograms). MCS is then straightforward to apply but a large number of realizations are needed e.g. the mean value typically converges in 1/ √ S. Slow convergence rates may incur excessive computational costs, specially when one deterministic simulation is already computationally costly. To accelerate this computational cost, different techniques have been developed in the literature. Among them, Quasi Monte Carlo Simulation, Stratified Sampling and the Latin Hypercube Sampling.

Moment equations

In this approach, the objective functional is replaced by its Taylor series expansion around the mean value (of the input random variable). The mean and the variance of the output function of interest are computed using the moments of the input variables. This leads to a more efficient integration of the first and second order moment. The proof is done below only for the 2 first statistical moments using a first order (y FO ) and a second order (y SO ) Taylor series expansion of the output function.

Let y the output function of interest where ξ = {ξ 1 , . . . , ξ M } with mean ξ = { ξ1 , . . . , ξM } and let us define its first and second order Taylor series approximation by y FO and by y SO [START_REF] Putko | Approach for uncertainty propagation and robust design in cfd using sensitivity derivatives[END_REF][START_REF] Martinelli | On the use of second-order derivatives and metamodel-based monte-carlo for uncertainty estimation in aerodynamics[END_REF]. It is respectively assumed that y ∈ C 1 and y ∈ C 2 .

y FO (ξ) = y( ξ) + M i=1 ∂y ∂ξ i (ξ i -ξi ) (2.21)

Introduction and general concepts

y SO (ξ) = y( ξ) + M i=1 ∂y ∂ξ i (ξ i -ξi ) + 1 2! M j=1 M i=1 ∂ 2 y ∂ξ i ∂ξ j (ξ i -ξi )(ξ j -ξj ). (2.22)
After some algebra, the expected value from the mean (unbiased) and from the variance of the output y are approximated by:

ȳFO = y FO ( ξ) (2.23) σ 2 (y FO ) = M i=1 ∂y FO ∂ξ i σ(ξ i ) 2 (2.24)
and by

ȳSO = y SO ( ξ) + 1 2! M j=1 M i=1 ∂ 2 y ∂ξ i ∂ξ j σ(ξ i )σ(ξ j ) (2.25) σ 2 (y SO ) = M i=1 ∂y SO ∂ξ i σ(ξ i ) 2 + 1 2! M j=1 M i=1 ∂ 2 y ∂ξ i ∂ξ j σ(ξ i )σ(ξ j ) 2 (2.26)
In some special cases, where the first and second order derivatives are available, the integration of the first and second order moment may be realized analytically as shown in equations Eq.2.23,2.25 and Eq.2.24, 2.26. However, in most of the case, finite difference schemes are required which leads to unaffordable computational costs when M increases. Moreover, the issue of the truncation the Taylor series expansion is still an open issue as in general no fundamental result proves the monotonic convergence of the Taylor series expansion.

Methods dedicated to reliability studies Performing reliability analysis of a mechanical system or process results in determining the probability of a particular event ξ failure leading to the failure of the structure (the "failure event"). This failure event may usually be expressed as the negative value (by convention) of a criterion G(ξ failure ). The probability of failure P failure of the system is then the probability that G(ξ failure ) < 0. The set D f = {ξ failure |G(ξ failure ) < 0} defines the failure domain. The function Ξ → ∂D f = {ξ failure |G(ξ failure ) = 0} defines the border of the failure domain also called limit state function.

First and Second Order Reliability Method These methods consist in approximating the failure domain G by a simpler domain whose probability may be computed analytically. A preliminary step to the construction of the simpler domain consists in transforming the observed random variable ξ into a vector of independent centered normalized Gaussian random variables. In this particular space, the most probable failure point ξ * also called "conception point" has to be identified. In the transformed space the conception point is the nearest point to the origin (this is due to the "bell symmetric" shape of the Gaussian probability density function.) Computing this point may be done by solving the following optimization problem min(∥ξ∥ 2 )

s.t G(ξ failure ) ≤ 0 . (2.27)

Introduction and general concepts

Many algorithms have been proposed in the literature to solve this problem and find the ξ * which may be referred in the literature as the conception point, the beta point of most probable point (MPP).

Once the MPP has been identified, the FORM consists in building a linear approximation of the limit state function points by defining a hyperplane passing through the design point and orthogonal to the vector ξ * . The probability of failure is then simply approximated by:

P failure ≈ Φ(-β), (2.28)
where Φ is the standard Gaussian cumulative distribution function and β is the signed distance between the origin and the conception point. The SORM improves the FORM methodology by proposing a second order approximation of the limit state function defined by a paraboloid tangent to the limit state at the conception point. The curvatures of the paraboloid at the design point are computed as the eigenvalues of the Hessian matrix of the limit state function. In this context, the failure probability may then be approximated by :

P failure ≈ Φ(-β) M i=1 (1 -α i β) -1/2 (2.29)
where the α i are the principal curvatures of the paraboloid. For relatively small stochastic dimensions and smooth continuous limit state these methodologies provide relatively accurate results with reasonable computational time. However, they are not adapted to discontinuous or multi-modal limit state functions and may lead in these cases to an erroneous estimation of the probability of failure. Moreover, they do not intrinsically provide an error estimator on the probability of failure. Finally, these methodologies require the need of derivatives which are not always available when the function of interest is given in an implicit way (as it is the case in non intrusive approaches). Moreover, in high dimension the computation of the Hessian matrix for the SORM may be prohibitively costly [START_REF] Zhao | A general procedure for first/second-order reliabilitymethod (form/sorm)[END_REF]. [START_REF] A D Kiureghian | Efficient algorithm for second-order reliability analysis[END_REF] proposes a more efficient iterative algorithm to compute the curvatures at the design point without the need of the Hessian matrix.

Sampling based approach To keep the control on the error committed on the estimated probability of failure, direct integration techniques to assess the reliability of the system may also be used.

Introduction and general concepts

The probability of failure is then written in its integral form as:

P failure = G ξdP ξ (ξ) = E[1 G (ξ)], (2.30)
where 1 G is the indicator function on G.

The computation of the above integral using brute Monte Carlo simulation becomes rapidly prohibitive in the transformed variable space as most of the sampling will concentrate around the origin of the domain.

It may be shown that for low P failure , the coefficient of variation C V of the MCS estimator converges in subset simulation [START_REF] Au | Estimation of small failure probabilities in high dimensions by subset simulation[END_REF]. These advanced methodologies are extended by [START_REF] Gi Schuëller | A critical appraisal of reliability estimation procedures for high dimensions[END_REF] to be efficient in high dimensions.

Metamodel based approaches A metamodel is a response surface coupled with a sampling strategy which is aimed at describing the behavior of a complex computational model while being less expensive to evaluate. In the context of structural reliability a number of them has been developed in the last decade in order to circumvent the costs and accuracy issues of the proposed FORM and SORM methodologies.

In the early 90's, [START_REF] Bucher | A fast and efficient response surface approach for structural reliability problems[END_REF][START_REF] Faravelli | Response-surface approach for reliability analysis[END_REF][START_REF] Wong | Slope reliability and response surface method[END_REF] proposes an adaptive second order polynomial response surface to interpolate the limit state function. To increase the efficiency and the accuracy of the metamodel, statistical information on the basic variables are used in order to update the obtained metamodel. Monte

Carlo simulations are then used in conjunction of the metamodel in order to assess the desired reliability

estimates. An obvious limitation of this methodology lies in the fact the limit state function has to be smooth enough and continuous. In some practical problems this is clearly not the case: the performance function may highlight multiple design points, and multiple regions that make significant contributions to the failure probability. This question is addressed by [START_REF] Gupta | An improved response surface method for the determination of failure probability and importance measures[END_REF] who proposes an algorithm to identify the multiple points and regions and to interpolate these points using a second-order polynomial response surface whose coefficients are determined using a least square analysis. A measure of the sensitivity of the reliability index is also provided. However, for non-smooth limit state function, the use of second order polynomial response surface may show insufficient accuracy in approximating the limit state function, and large errors may be observed in the computation of the sensitivity of the reliability index, particularly.

[YC04a], [KKC10] used Moving Least Squares giving higher weight to the experimental points closer to the design point and allowing the response surface function (RSF) to be closer to the limit state function. The authors in [DSB11b] take advantage of the kriging intrinsic features in order to propagate the approximation error on the limit-state surfaces to the failure probabilities estimates providing thus an empirical error measure. In [START_REF] Li | Evaluation of failure probability via surrogate models[END_REF] the authors argue that the straightforward sampling on surrogates models may lead to biased results and thus proposes a hybrid approach based on a large sampling of the metamodel in the probability space and refining this sampling using the numerical model in some region of interest. In [START_REF] Sudret | Meta-models for structural reliability and uncertainty quantification[END_REF][START_REF] Dubourg | Metamodel-based importance sampling for structural reliability analysis[END_REF] the author proposed an efficient scheme to solve the problem of the potential biasedness in the estimation of a probability of failure due to a direct substitution. It combines the use of an importance sampling strategy guided by a kriging metamodel replacing the indicator function

Spectral representation of stochastic field through Polynomial Chaos Expansion

of the failure domain by a probabilistic function.

Spectral representation of stochastic field through Polynomial Chaos Expansion

Introduction

This section reviews in more details the construction of the Polynomial Chaos Expansion (PCE). This methodology has been widely used as an alternative to the Monte Carlo Simulation in order to propagate the uncertainty through a deterministic numerical model. From its initial form, many obstacles have to be circumvented for an efficient use. In this section we briefly review these obstacles and the circumventing solutions.

Functional evaluation of a random variable by orthogonal polynomials

Functional evaluation of random variable

In practice, random variables of interest are seen as the output of a deterministic numerical model y(ξ).

They often belong to the set of square integrable functions denoted by L 2 (R M , dP ξ ) and defined by:

L 2 (R M , dP ξ ) = {y : ξ → y(ξ) ∈ R; E[y 2 ] := R M y(ξ) 2 dP ξ (ξ) < ∞} (2.31)
where Ξ ∈ R M since in practice one often deals with a finite number of second order (with finite second order moment) random variables ξ = {ξ 1 , . . . , ξ M }.

When this space is equipped with the following inner product, it becomes an Hilbert space:

< u, v > L 2 (R M ,dP ξ ) = E[uv] = R M u(ξ)v(ξ)dP ξ (ξ) (2.32)
An Hilbertian basis {B i } i∈I of L 2 (R M , dP ξ ) is a complete set of orthonormal functions verifying the following properties:

< B i , B j > L 2 (R M ,dP ξ ) = δ ij (2.33)
where δ i , j is the Kronecker symbol and

∀u ∈ L 2 (R M , dP ξ ), < B i , u >= 0, ∀i ∈ I ⇒ u = 0 (2.34)
Finally, each function y ∈ L 2 (R M , dP ξ ) admits an unique decomposition on this Hilbertian basis:

y(ξ) = i∈I γ i B i (ξ) (2.35)
where γ i , i ∈ I are the coordinate of y projected on the orthonormal basis {B i } i∈I :

2.2. Spectral representation of stochastic field through Polynomial Chaos Expansion

γ i =< γ, B i > L 2 (R M ,dP ξ ) = E[y, B i ] = R M y(ξ)B i (ξ)dP ξ (ξ) (2.36)
Several choices have been proposed in the literature in order to construct the Hilbertian basis. Among the most popular, one may find the Karhunen Loève (KL) expansion and the polynomial expansion methods.

Link between random variable and orthogonal polynomials

In probability theory, since the preliminary work of Wiener on Brownian motion, the relationship between orthogonal polynomials and random variables has been widely studied [START_REF] Schoutens | Stochastic Processes and Orthogonal Polynomials[END_REF].

Let P = {H N (ξ), N ∈ N} be a set of monovariate polynomials whose maximal degree is N . Let Φ be a real positive measure. P forms a set of orthogonal polynomials with respect to the measure Φ if:

Ω H n (ξ)H m (ξ)dΦ(ξ) = Kδ nm , n,m∈ N (2.37)
where δ nm is the Kronecker product and K is a non zero constant (if K=1, P is orthonormal). A well known property of real orthogonal polynomial lies in the fact that they all may be built using a three terms recurrence relationship on the form:

-ξH n (ξ) = b n H n+1 + γ n H n (ξ) + c n H n-1 (ξ), n∈ N (2.38)
where b n , c n ̸ = 0 and c n /b n-1 > 0

In the space L 2 (R, dP ξ ) of continuous second order random variable, the following classical polynomials form an orthogonal basis: they are referred in the literature to the classical orthogonal polynomials of the Askey-scheme:

Pdf type Density Function Orthogonal polynomial Support

Uniform

1 ]-1,1[ (ξ)/2 Legendre: P k (ξ) [-1; 1] Gaussian 1 √ (2π) e -ξ 2 /2
Hermite:

H k (ξ) [-∞; ∞] Gamma ξ a e -ξ 1 R + (ξ) Laguerre: L k a (ξ) [0; +∞] Exponential e -x Gen. Laguerre L (α) n (x) [0; +∞] Beta 1 ]-1;1[ (ξ) (1-ξ) a (1+ξ) b B(a)B(b) Jacobi: J k (ξ) [-1; 1]
Table 2.1: Classical orthogonal polynomial of the Askey scheme with respect to their probability density function

The resulting expansion of the random variable of interest is readily:

y(ω) = i∈J γ i H i (ξ) (2.39)
where J ⊂ N, {H i } i∈J is an infinite series of one of the Askey-scheme polynomial, ξ S the associated

Spectral representation of stochastic field through Polynomial Chaos Expansion

standard variable, {γ i } i∈J the coordinate of the function on the orthogonal polynomial basis.

Construction of Wiener Polynomial Chaos Expansions (PCEs)

N. Wiener [START_REF] Wiener | The homogeneous chaos[END_REF] proposes for standard Gaussian variable the use of Hermitian polynomial as an Hilbertian basis to construct an approximation of the random variable y(ξ) ∈ L 2 (Ξ, dP ξ ). In the following we describe the construction of the original PCE. It is usually named Hermite homogeneous PCE in the literature.

Univariate case

Let ξ be any component of ξ supposed to be a Gaussian standardized random variable (which may be obtained using an iso-probabilist transformation [LCM09]) the expansion is readily given by Eq.2.39

where {H i } i∈N are the mono-variate Hermite polynomials. The first 5 Hermite's polynomials are given below:

H 0 (ξ) = 1 H 1 (ξ) = ξ 2 -1 H 2 (ξ) = ξ 3 -3ξ H 3 (ξ) = ξ 4 -6ξ 2 + 3 H 4 (ξ) = ξ 5 -10ξ 3 + 15ξ
A fundamental result for the Hermite homogeneous PCE is established by the Cameron Martin theorem [START_REF] Cameron | The orthogonal development of non-linear functionals in series of fourier-hermite functionals[END_REF]. It states that any second-order functional of the Brownian motion (Gaussian random field) can be expanded as a mean-square convergent series in terms of infinite-dimensional Hermite polynomial in Gaussian variables. Thus, the exact (in a mean square sense) PCE of the random variable y(ξ) is readily:

y(ξ) = i∈N γ i H i (ξ) (2.40)
where H i is the i th degree Hermite polynomial. Exploiting the polynomial orthogonality, the coefficients of the expansion are defined as:

∀i ∈ N, γ i = < y(ξ), H i (ξ) > ∥H i ∥ 2 (2.41)

Multi-variate case

In the multi-variate case ξ = {ξ 1 , ξ 2 , . . . , ξ M } with M coordinates, according to the independence assumption, the joint probability density function f ξ (ξ) may be decomposed on a product of the marginal probability density functions f ξ (i) (Eq.2.42):

f (ξ) = M i=1 f ξi (ξ i ) (2.42)
Given the natural inner product for arbitrary function φ with respect to each of the marginal proba-31 2.2. Spectral representation of stochastic field through Polynomial Chaos Expansion

bility function f ξ (ξ) defined on D < φ 1 , φ 2 >= D φ 1 (ξ)φ 2 (ξ)f ξ (ξ)dξ (2.43)
Using the tensor product on these mono-variate polynomials one may obtain an infinite set of multi-variate polynomials (with a preserved orthogonality property)

H = {H α , α ∈ N M } where α = {α 1 , . . . , α M } ∈ N M is a multi-index set. H α (ξ) = ⊗ M i=1 H αi (ξ i ) (2.44)
According to the theorem of Cameron Martin [START_REF] Cameron | Transformations of wiener integrals under translations[END_REF], the exact polynomial expansion of the functional y is readily

y(ξ) = α∈N M γ α H α (ξ), (2.45)
where {γ α }, α ∈ N M are the coefficients of the PCE to be identified

For practical use, one may truncate the full set of tensor product polynomials in order to only retain a finite set of polynomial terms.

Classical truncation scheme Among all {H α , α ∈ N M } the classical truncation scheme [Gha91, BSL06a, SDK00a] retains only the multi-variate polynomial terms whose degree does not exceed an arbitrarily fixed N leading to the following multi-index set:

A M q = {α ∈ N M , ||α|| q ≤ p}, (2.46)
where ||α|| q = M i=1 α q i 1/q and q = 1. Then, the truncated model may be written as:

y A M q (ξ) ≈ α∈A M q γ α H α (ξ).
(2.47) Fig. 2.2 gives an illustration of this scheme for a 2-variate 7 th order PCE. The number P of coefficients in the PCE is given by

Spectral representation of stochastic field through Polynomial Chaos Expansion

P = N k=0 C k M +k+1 = (N + M )! N !M ! (2.48)
and increases exponentially both with N and M .

A list of the multi-variate Hermitian monomials for M = 2 and N = 3 is provided below:

i α = {α 1 , α 2 } H i 0 {0, 0} H 0 (ξ 1 )H 0 (ξ 2 ) = 1 1 {1, 0} H 1 (ξ 1 )H 0 (ξ 2 ) = ξ 1 2 {0, 1} H 0 (ξ 1 )H 1 (ξ 2 ) = ξ 2 3 {2, 0} H 2 (ξ 1 )H 0 (ξ 2 ) = ξ 2 1 -1 4 {1, 1} H 1 (ξ 1 )H 1 (ξ 2 ) = ξ 1 ξ 2 5 {0, 2} H 0 (ξ 1 )H 2 (ξ 2 ) = ξ 2 2 -1 6 {3, 0} H 3 (ξ 1 )H 0 (ξ 2 ) = ξ 3 1 -1ξ 1 7 {2, 1} H 2 (ξ 1 )H 1 (ξ 2 ) = (ξ 2 1 -1)ξ 2 8 {1, 2} H 1 (ξ 1 )H 2 (ξ 2 ) = ξ 1 (ξ 2 2 -1) 9 {0, 3} H 0 (ξ 1 )H 3 (ξ 2 ) = ξ 3 2 -3ξ 2
(i) = y(ω i ) = {y(x 1 , ω i ), . . . , y(x N , ω i ))} ⊤ , i = 1, .
. . , S be a real column vector in R N . Let {ȳ 1 , . . . , ȳN } ⊤ the vector of expectations:

∀i ∈ {1, . . . , N}, ȳi = 1 S S j=1 y(x i , w j ). (2.51)
Karhunen-Loève decomposition The full Karhunen-Loève expansion [START_REF] Loeve | Probability theory[END_REF] is aimed at efficiently reducing the statistical complexity of the random field by achieving its decomposition as a denumerable sum of product functions of the form:

y(ω) = ȳ + N i=1 λ i κ i (ω)ϕ i (2.52)
decoupling the stochastic Ω and deterministic D spaces.

The suited decomposition may be obtained by orthogonal projection of the stochastic field of interest onto separable Hilbert spaces. The convenience of the Karhunen Loève (KL) expansion stays precisely in the fact that KL expansion is toptimal with regards to the total mean square error. To construct, the Hilbertian basis one performs an eigenvalue decomposition of the covariance operator which may be defined as

C(x i , x j ) = D y(x i , ω)y(x j , ω)dp ω (ω), (i, j) ∈ {1, . . . , N} 2 . (2.53)
Replacing y by its collocation vector y, an unbiased estimator of the covariance operator may be obtained by:

C(x i , x j ) = 1 S -1 S i=1 (y (i) -ȳi ) ⊤ (y (j) -ȳj ), (i, j) ∈ {1, . . . , N} 2 .
(2.54)

By definition, the covariance kernel function (since D is bounded), symmetric and positive definite.

Spectral representation of stochastic field through Polynomial Chaos Expansion

It admits a decomposition of the form:

C = N i=1 λ i ϕ i ϕ i ⊤ = N i=1 φ i φ i ⊤ (2.55)
where {λ i } N i=1 are the eigenvalues, and {ϕ i } N i=1 the eigenfunctions. They constitute the elements of the Hilbertian basis

{φ i } N i=1 = { √ λ i ϕ i } N i=1
. From 2.55, one obtained the following decomposition y KL of y(ω):

y KL (ω) = ȳ + φφ ⊤ (y(ω) -ȳ) = ȳ + N i=1 φ i κ i (ω) (2.56) where we have introduced {κ i } N i=1 = {φ ⊤ i (y(ω) -ȳ)} N i=1 .
Out of this definition one may deduce the following properties:

• κ i are second order random variables • E[κ i ]= 0 • E[κ i κ j ]= 0 i ̸ = j ∈ N
Consequently, the covariance matrix of κ is the identity matrix insuring the independence between the vector components. Moreover, since the underlying random field y(ω) is generally not gaussian, the κ variables are in general non Gaussian variables. Due to the symmetry and to the positiveness of the covariance kernel, the eigenfunctions ϕ i form a complete set and exhibit an orthogonality property

φφ ⊤ = I N ×N , (2.57)
where I N ×N is the N × N identity matrix.

Regarding the eigenvalues {λ i } N i=1 they form a decreasing sequence of positive values conventionally ordered as follows: λ 1 > λ 2 > λ 3 > . . . → 0. Moreover, with the additional assumption of mean-square continuity of the random field, this sequence is convergent. From these assertion one may prove that the Karhunen Loève expansion is unique.

In the following we show that the Karhunen Loève expansion is optimal in the mean square sense, i.e increasing the truncation the order of the expansion reduces the total mean square error.

Let us first introduce the following additional notations:

-Let n < N and let us consider {φ} n i=1 the n first eigenvectors, we denote them φ t , -Let Y = {y(ω 1 ), . . . , y(ω S )} and

-C = S i=1 (y (i) -ȳ)(y (i) -ȳ) ⊤ = Y Y ⊤ .

Spectral representation of stochastic field through Polynomial Chaos Expansion

A truncated expansion y KLt (ω) of y(ω) is readily

y KLt (ω) = ȳ + φ t φ ⊤ t (y(ω) -ȳ) . (2.58)
The truncation error ϵ is given by:

ϵ = S i=1 ∥y(ω i ) -y KLt (ω i )∥ 2 = y (i) -ȳ -φ t φ ⊤ t (y (i) -ȳ) ⊤ y (i) -ȳ -φ t φ ⊤ t (y (i) -ȳ) = S i=1 (y (i) -ȳ)(I -φ t φ ⊤ t ) ⊤ (I -φ t φ ⊤ t )(y (i) -ȳ) = S i=1 (y (i) -ȳ)(I -φ t φ ⊤ t )(y (i) -ȳ) = tr Y ⊤ I -φ t φ ⊤ t Y = tr C I -φ t φ ⊤ t = tr (C) -tr φ ⊤ t Cφ t = N i=1 λ i - n i=1 λ i = N i=n+1 λ i
Thus when increasing the order of the expansion one automatically reduces the total mean square error. In practice one truncates the expansion to a finite number of summands leading to :

y(ω) = ȳ(ω) + n i=1 φ i κ i (ω) (2.59)
The truncated order may be chosen such that for some n ≤ N ,

n i=1 λ i N i=1 λ i is closed enough to one.
Reduced PCE expansion Let the random n variables κ j be collected in to a random vector κ with value in R n and let us approximate them using a PCE of order P and dimension n:

κ j (γ) = P α,|α|=1 γ α H α (ξ) (2.60) where |α| = M j=1 α j .
Injecting this reduced expansion into the KL decomposition of the field y, one obtains:

2.2. Spectral representation of stochastic field through Polynomial Chaos Expansion

y(ω) = ȳ + n i=1 ⎛ ⎝ P αi,|αi|=1 γ αi H αi (ξ) ⎞ ⎠ φ i . (2.61)
We note that ξ is valued in R n where n is chosen sufficiently large in order to limit the accuracy loss due to the dimension reduction process. Considering the difference between the tails of the true function and the PCE to the upper 1% fractile, in some cases non monotonic convergence are observed: for some non-gaussian random variables increasing the number of terms may not improve the distribution tails making the use of the PCE approaches risky for reliability analysis. Moreover, the authors in [START_REF] Field | Convergence properties of polynomial chaos approximations for l2 random variables[END_REF] provide mathematical proofs showing that higherorder statistical moments (i.e., greater than two) computed using homogeneous hermitian PCE may not always converge with the number of terms in the Hermite-PCE series. For example, it is shown that the third absolute moment of the PC approximation for a lognormal random variable does converge, while moments of order four and higher for uniform random variables do not converge. It has been previously demonstrated through numerical study in [START_REF] Field | On the accuracy of the polynomial chaos approximation[END_REF]) who proved that a lack of convergence in the higherorder moments can have a profound effect on the rate of convergence of the tails of the distribution of the PC approximation.

On the balance between computational costs and accuracy

The highest challenge consists then in the computation of the PCE coefficients. One may then have to face many different problems:

• in high dimension, the number of terms to compute increases exponentially with the number of variables and with the order of the PCE which may in some cases improve some metrics. This problematic is called the "curse of dimensionality". To tackle this issue, a special care has to be simulatneously given to the reduction of the number of calls of the exact model in order to compute 2.2. Spectral representation of stochastic field through Polynomial Chaos Expansion a fixed number of coefficient, and to adapt the polynomial basis in order to only capture the basis elements the most significant with regards to the stochastic function to be expanded.

• limited set of data or erroneous set of data introduces [START_REF] Jr | A probabilistic approach to uncertainty quantification with limited information[END_REF] its own type of uncertainties. It is then of paramount importance to build a PCE with regards to the data at hand and to be able to evaluate the influence of missing and/or erroneous data on the PCE prediction capability. XH05, Xiu07a, Xiu09a, Xiu09b], showed (without providing mathematical proofs) that in many cases a higher convergence rate may be obtained using the orthogonal polynomials of the Askey-scheme in lieu of the Hermite polynomial.

However, for these generalized polynomial basis, despite the demonstration of their efficiency in many cases no proof of convergence has been established until very recently. The authors in [START_REF] O G Ernst | On the convergence of generalized polynomial chaos expansions[END_REF] notably

show that an arbitrary random variable with finite variance may be expanded in generalized polynomial chaos expansions if the underlying probability measure is uniquely determined by its moment and continuous giving by the way the missing theoretical insight to the gPC.

arbitrary PCE (aPCE)

Finally, arbitrary Polynomial Chaos expansion (aPC) [START_REF] Oladyshkin | Data-driven uncertainty quantification using the arbitrary polynomial chaos expansion[END_REF][START_REF] Riesz | Sur le problème des moments et le théorème de parseval correspondant[END_REF] proposed the construction of a PCE basis using only the statistical moments of the random variable computed from the data at hand. It

shows that the construction of a polynomial basis of degree d needs the knowledge of a finite number (2d) of the statistical moments of the continuous random variable to be expanded without requiring any knowledge of a probability density function. The degree of the expansion is also required to be not greater than the number of available sampling points. The stochastic analysis is handled using only the data at hand, hence avoiding the subjective association of a parametric probability function to the limited set of data. The complexity of the stochastic analysis is thus aligned with the reliability and

Spectral representation of stochastic field through Polynomial Chaos Expansion

the detail level of the statistical information at hand. However, when limited or erroneous data sets are available epistemic uncertainties are introduced in the model. The proposed methodology does not offer any intrinsic property in order to assess the impact of these kind of uncertainties. External methods to assess their impact are then proposed by the authors such as jacknife and bootstrap. Finally, the authors highlight an improved convergence rate of this approach with regards to the gPC approach on the same limited set of data, but no convergence proof of such an approach has been provided yet. At the current state of the art, no work has been proposed yet in order to extend the aPC technique for correlated random variable.

Multi-Element Generalized PCE

The previously proposed generalizations, the Hilbertian basis used for the projection of the random field over the whole stochastic space is made of continuous polynomial functions: they do not allow to treat non continuous probability density function. To address this limitation, the authors in [WK05, WK06a, This methodology shows its limitation in high dimension when no other solution may exists in order to build a partition of the stochastic space.

Computing the PCE coefficients in adequacy with computational ressources

Another well addressed limitation of the PCEs (Hermite, generalized, or arbitrary) appears when the stochastic dimension increases. As shown by Eq. 2.48, when the dimension increases the number of terms in the PCE increases exponentially and may rapidly lead to a situation where the construction of the PCE metamodel itself may become too computationally expensive. This situation is often referred to as an aspect of the "curse of dimensionality". To circumvent this curse, a handfull of methodologies are proposed in the literature proposing to retain only the number of polynomial in the expansion which are significant to describe the response variability.

Low rank index set truncation

The most intuitive and direct approach relies on the "sparsity of effects principle". It states that a system is usually dominated by main effects and low-order interactions. Q-norm generalizes the classical truncation scheme by varying 0 ≤ q ≤ 1 [BS08a]. Fig. 2.3 illustrates a typical truncated index set different for q = 0.6 and a 7 th order 2-variate PCE.

(a) q = 0.6 The set of active polynomials in the PCE decomposition is decreased when q decreases. Fig. 2.4(a)

illustrates the evolution of the number of 2-variate polynomial terms in linear scale against q values, and Having performed this a priori truncation scheme, the PCE coefficients may then be computed using either a projection scheme exploiting the orthogonality of the PCE basis, or using a collocation scheme.

The cost reduction is done using so called "sparse approaches" were the sparsity either characterizes the Design of Experiment (DoE) used to compute the projection coordinates of the random data on the orthogonal polynomial basis (PCE coefficients), or the polynomial basis itself.

γ α = < y, H α (•) > < H α (•), H α (•) > = R M y(ξ)H α (ξ)Φ(ξ)dξ R M H α (ξ) 2 H α (ξ)dξ = E[y(.)H α (•)] E[H 2 α ]
(2.62)

Each tensor product involves the computation of a multi-dimensional integral over the support range of the multivariate probability measure. Let us assume that the random variables ξ = {ξ 1 , . . . , ξ M } are all independent standard Gaussian random variables, hence with M i=1 φ(ξ i ) = Φ(ξ) is a M -multivariate standardized Gaussian probability density function obtained from the product of the M one-dimensional standardized Gaussian probability density function.

The integral equation on the denominator is known and easy to compute in the Hermitian case:

E[H 2 α ] = < H α , H α > = R H 2 α (ξ)H α (ξ)dξ = Ξ1 . . . Ξ M Ψ 2 j (ξ 1 , . . . , ξ M )φ(ξ 1 ) . . . φ(ξ M )dξ 1 . . . dξ M = Ξ1 . . . ΞM M j=1 H 2 αj (ξ j ))φ(ξ 1 ) . . . φ(ξ M )dξ 1 . . . dξ M = M j=1 Ξj H 2 αj (ξ j )φ(ξ j )dξ j = M j=1 α j ! (2.63)
Only the integral equation on the numerator is still to be computed

R M y(ω)H α (ξ)dξ (2.64)
The key point is to choose the collocation points. In one dimension the choice is straight forward. A one dimensional quadrature rule is aimed at approximating an integrand of the form R y(ξ)φ(ξ)dξ by a finite summand number of S terms composed by the product of "weights" {w 1 , . . . , w S } and of evaluations of the output function y(ξ) at some quadrature abscissas {ξ (1) , . . . , ξ (S) } ∈ R.

U(y) = S i=1 y(ξ (i) )w i (2.65)
For example the Gaussian Quadrature rule abscissas are composed of the zeros of the polynomials which are orthogonal to the probability density weighting function and referred in the Wiener-Askey scheme. Numerous studies show that the Gaussian quadrature rule is usually the best choice.

The tensor product of one dimensional quadrature rule generalizes this approach in higher dimension.

Let ξ = {ξ 1 , . . . , ξ M }, be a M dimensional vector and i ∈ {1, . . . , M} an index referring to the variables dimensions. For each dimension i, let {ξ

i , . . . , ξ

} be a sequence of abscissas for quadrature on R. For f ∈ C 0 (R) and n = 1 let us introduce a sequence of one dimensional quadrature operators. where {w k 1 . . . w k M } and {ξ k 1 , . . . , ξ k M } are respectively the integration weights and integration points.

U i (y) = mi j=1 y(ξ i j )w i j (2.
The integration order K in each direction is assessed by considering the following assertion: if the response y(ξ) was polynomial of order p the terms to be integrated (y i ξ)Ψ j (ξ)) would be of maximal order 2p and a K = p + 1 order integration scheme would give the exact solution. The integration weights and points are computed according to the probabilistic measure considered. In this context for each direction, the integration weights are computed in the following way:

w k = < H k , H k > H ′ k (ξ k )H k-1 (ξ k )
with H k the k th order mono-variate Hermitian polynomial.

(2.68)

Considering now the integration points, it is shown that they have to be taken as the roots of the maximal possible order (2p) monomial Hermitian polynomial appearing in the terms to be integrated.

This approach has been used in a number of recent papers notably [GS93, MH03, MK05, XK03a] and some investigations towards an error estimation has been analyzed in [START_REF] Babuška | A stochastic collocation method for elliptic partial differential equations with random input data[END_REF].

The main limitation of this approach is the exponential growth of the computational cost with the number of variables. If s collocation points are chosen in each of M -dimension, the number of total calls to the exact model is S = s M . Realizing this number of simulations becomes rapidly unfeasible particularly when high computational efforts are demanded for one evaluation of y(ξ) and thus limit the use of the tensor product of integration rules to lower dimensions.

Sparse approaches based on Smolyak's algorithm In projection based approaches, the challenge consists in computing at low computational costs the integral Eq.2.64 when the number of random variables become moderately large. Smolyak [?] proposes to consider sparse tensor product of quadrature rules.

The Smolyak sparse grid quadrature rule is aimed at drastically reducing the number of collocation points while preserving a high level of accuracy. It consists in tensor product of quadrature rules retaining only those where a small number of points are used by the following rules:

For U 0 = 0 and for i ≥ 1, and |i| = i 1 + . . . , i n , the Smolyak quadrature formula is defined by:

A(w, n) = |i|≤w+n (∆ i1 ⊗ . . . ⊗ ∆ in ) (2.69)
where

∆ i = U i -U i-1 (2.70)

Spectral representation of stochastic field through Polynomial Chaos Expansion

which may be equivalently written as:

A(w, n) = w+1≤|i|≤w+n (-1) w+n-|i| n -1 w + n -|i| (U i1 ⊗ . . . ⊗ U in ) (2.71)
For each level index i, the number of points considered in the one dimensional quadrature is ruled by a growth function which may be linear or non linear:

Clenshaw-Curtis rules m = ⎧ ⎨ ⎩ 1i = 1 2 i-1 + 1i > 1 (2.72) Gauss-Patterson rules m = 2 i -1, (2.73)
Gausisan rules m = 2i -1.

(2.74)

As an extension of this idea, [START_REF] P G Constantine | Sparse pseudospectral approximation method[END_REF] re-examine the Smolyak's algorithm in order to reconstruct the coefficients that naturally corresponds to the sparse grid integration rule. The number of terms in the obtained expansion is consistent with the number of points in the sparse grid integration rule.

The key point of the proposed approach is to separately compute the coefficients of the tensor product polynomial expansion for each tensor grid in the sparse grid. Then, the linear combination of the tensor weight is used to linearly combine the coefficients of each tensor expansion. This method produces a pointwise equivalent polynomial surrogate to the one constructed from a linear combination of tensor product Lagrange polynomial.

Cubature rules Cubature rules are specifically adapted to the computation of integrand in high dimension. As in Eq.2.66 they aim at transforming a multi-dimensional integral in a weighted sum, without though using the tensor structure of the multivariate stochastic space. An extensive review of cubature rules may be found in [START_REF] Cools | An encyclopaedia of cubature formulas[END_REF][START_REF] Haber | Numerical evaluation of multiple integrals[END_REF]. Cubature rules are characterized by a degree "d" for which the equation 2.66 is exact if the integrand is any multi-variate polynomials of degree at most d but not d + 1. As an example, a set of d = 2 cubature rule discussed in [START_REF] Xiu | Numerical integration formulas of degree two[END_REF] are generalized for arbitrary integration domain with arbitrary probability measures. Due to the excellent consistency between the number of polynomials in the expansion and the number of required terms, these methods may be highly efficient for smooth (the integration degree remains quite low) high dimensional problems. However these cubature rules may only be built using homogeneous random variables. Moreover, the interpolation accuracy can not systematically be reduced. These limiting assumptions make them superseded in the literature and sparse grid integration rules are preferred.

Regression approach to compute PCE parameters Classical Collocation Approach Another scheme referred as collocation (or regression scheme) consists in performing least squares regression on the polynomial chaos coefficients. For each design point

x, a set of S scattered data in the standardized random variable space Ξ = {ξ (1) , . . . , ξ (S) } are sampled.

Spectral representation of stochastic field through Polynomial Chaos Expansion

Then on each of these points the mechanical model (high fidelity simulations) PCE coefficients γ are obtained by minimizing the least square residual ϵ LS

ϵ LS = ⎧ ⎪ ⎪ ⎨ ⎪ ⎪ ⎩ y(ω 1 ) . . . y(ω S ) ⎫ ⎪ ⎪ ⎬ ⎪ ⎪ ⎭ - ⎡ ⎢ ⎢ ⎣ H 1 (ξ (1) ) • • • H j (ξ (1) ) . . . . . . H 1 (ξ (S) ) • • • H j (ξ (S) ) ⎤ ⎥ ⎥ ⎦ ⎧ ⎪ ⎪ ⎨ ⎪ ⎪ ⎩ γ 1 . . . γ j ⎫ ⎪ ⎪ ⎬ ⎪ ⎪ ⎭ 2 .
(2.75)

It has to be noticed that the results issued from the collocation method is quite sensible to the samples configuration of DOE points: for different sampling, different accuracies on the output can be obtained.

Moreover, the number of sampling points taken to solve the system, also influences the accuracy on the outputs. An empirical rule shows that the number of numerical experiments Q giving the best accuracy has to be chosen such as

Q ≈ (M -1) × P ([BSL06b]).
Improved Collocation method Collocation methods can be improved by performing a smarter sampling of the collocation points [START_REF] Berveiller | Stochastic finite element : a non intrusive approach by regression[END_REF]. Another alternative for the choice of the collocation points consists in constructing them from the roots of the P + 1 th order mono-variate Hermitian polynomial, where P is the maximal order of the polynomial chaos. Having constructed the M -uplets formed by all the possible combinations of the P + 1 roots for each stochastic variables ξ, only P of them are chosen to solve the system. As ξ follows a multi-normal standardized gaussian distribution, the P M-uplets EHJT04b, HTF09, HCMF08]) to efficiently select the most relevant polynomial basis elements of the full polynomial chaos expansion. The combination of these methods are particularly interesting since they may be use to build iteratively a sparse polynomial chaos metamodel using a small number approximation.

retained
By "small", the authors the authors mean not greater than the number of terms in the expansion, or even possibly significantly smaller. The selection of the most "relevant" polynomials coefficients are performed using statistical correlation measures between the current models and the exact models evaluated on a reference population of small size. The authors also build an adaptive algorithm to stop adding terms in the polynomial expansion when overfitting is detected. This algorithm is applied to a rigid frame with 21 stochastic variables. The LARS adaptive algorithm is used to build a sparse polynomial representation of the displacement of the structure. [HCMF08]

Adressing the limited number and erroneous training data

An alternative approach to compute the PCE coefficients, is to characterize them as a random variables.

Some approaches have been developed in very recent years in the literature tackling this problem most of the time in an intrusive manner. A very interesting feature of this proposed approach is that they set the computation of the PCE coefficients as a stochastic inverse problem. The coefficients of the PCE are 2.2. Spectral representation of stochastic field through Polynomial Chaos Expansion seen as random variables. Following the same idea as for aPC, classical direct methods for (maximum likelihood, maximum entropy, Bayes' theorem) are applied in order to assess the most probable values of these coefficients with regards to the data at hand. This approach has been investigated in the framework of stochastic inverse problem where the characteristics of random fields have to be identified using some of their realizations. Most of the time the complexity of the random data of interest is reduced using the KL decomposition:

y(ω) = n i=1 κ (i) (ω)φ i .
(2.76)

The random variables κ = {κ 1 , . . . , κ n } are then expanded onto a PCE leading to the following approximation:

y(ω) = ȳ + n i=1 ⎛ ⎝ P αi,|αi|=1 γ αi Ψ αi (ξ) ⎞ ⎠ φ i .
(2.77)

The coefficients to be identified and seen as random variables are the Γ = {γ} n i=1 . In a set of papers, [DGS + 05b, DGS06, DSG + 06, DSG07] investigate the use of the maximum likelihood in order to identify the coefficients Γ of the reduced PCE expansion . The likelihood function is then written as a function of the coefficients:

L(Γ) = N i=1 p κ (κ (i) |Γ) (2.78)
where Γ = {γ α , |α| ≤ q}. Then identifying the γ set of coefficients lead to the maximization of L(γ):

Γ * = argmax Γ L(Γ).
(2.79)

The authors in [DGS + 05b] propose the use of the characteristic function to describe each probability density function. Such a methodology may be costly in high dimension and is not always applicable. In a general way, the computation of the probabilities involve Monte Carlo simulations, thus high computational costs. Moreover, another issue consists in performing efficiently the optimization problem on the likelihood function. The authors in [PSBP07, PFS10, BSS07] specially address these issues.

To circumvent this issue, an alternative method has been proposed quite recently and encounter a growing success. It consists in considering the Bayesian approach in order to estimate the posterior probability of the coefficients knowing some evaluation of the κ variables. The Bayes' rule is used: If no information is available on the set of γ coefficients, the prior distribution p γ (γ) a non-informative probability distribution [AGS10] may be arbitrary defined. As usual, the likelihood

p(γ|y) = p(y|γ)p γ (γ) p(y|γ)p γ (γ)dγ (
L(γ) = n j=1 p κ (κ (j) |γ) (2.82)
may be performed using Monte Carlo simulation or more efficiently using the kernel density method [START_REF] Scott | Multivariate density estimation: theory, practice, and visualization[END_REF]. This framework has firstly been proposed in [START_REF] Ghanem | On the construction and analysis of stochastic models: characterization and propagation of the errors associated with limited data[END_REF] and further developed in [AGS10, MDMV12, BSS10].

Adressing simultaneously the two problematics

Finally, [Soi10b, SD10] proposes a general framework based on the stochastic identification of the PCE coefficients of a reduced random field expansion aiming at extending the originally proposed approach to high-dimension (several millions of coefficients) polynomial chaos expansions with random coefficients for non-Gaussian tensor-valued random fields using partial and limited experimental data. The proposed methodology consists in the following steps:

• introducing a family of prior probability models,

• identifying an optimal prior model in the constructed family using the experimental data,

• constructing a statistical reduced order optimal prior model,

• constructing the polynomial chaos expansion with deterministic vector-valued coefficients of the reduced order optimal prior model and finally,

• constructing the probability distribution of random coefficients of the polynomial chaos expansion and identifying the parameters using experimental data.

2.2.7

Why is PCE efficient to perform UQ analysis?

PCE for robust design

Having computed the set of PC coefficients, a number of convenient analytical features allows to make the use of these techniques particularly attractive for global and local sensitivity analysis, and design under uncertainty in a general way. In fact, the analytical moments of the responses may be obtained in closed form:

E[y(ω)] = α,|α|≤P E[γ α H α (ξ)] = γ 0 σ 2 [y(ω)] = α,|α|≤P E[γ α H α (ξ)] 2 -E[y(ω)] 2 = α,|α|≤P γ 2 α E[H 2 ]
(2.83)

Local sensitivity analysis Introducing non probabilistic variables such as design variable for example gives the following expression:

2.2. Spectral representation of stochastic field through Polynomial Chaos Expansion

y(ξ) ≈ α,|α|≤P γ α (x)H α (ξ) (2.84)
The sensitivity of the mean and variance responses may be obtained:

dµ dx = d dx γ 0 (x) (2.85) dσ 2 dx = α,|α|≤P d dx γ α (x) 2 (2.86)
where d dx γ α is defined as (the theorem of uniform convergence being verified):

d dx γ α = < dy dx , H α > < H 2 α > (2.87)
Due to the independence property, the coefficients computed in the last equation may be interpreted either as the nonprobabilistic sensitivities of the chaos coefficients for the response expansion, or as the chaos coefficients of an expansion for the nonprobabilistic sensitivities of the response. The resulting expansion are valid only for a particular set of nonprobabilistic variables and must be in a first approach computed again each time the nonprobabilistic variable are modified.

Another alternative to compute local sensitivities consists in building the stochastic expansion over both stochastic ξ and deterministic x variables. Considering a bounded domain x upper bound ≤ x ≤

x lower bound with no probabilistic content the PCE expansion then becomes:

y(x, ω) ≈ α,|α|≤P γ α H α (x, ξ) (2.88)
where Ψ may be built over the tensor product of optimal orthogonal polynomials with regards to the probability measure for each component of ξ and by Legendre polynomial (orthogonal with respect to the uniform distribution) for each component of x. In this case, the statistical moment of interest are computed using the scalar product only on the probabilistic variables:

E[y] ξ (x) = α,|α|≤P γ α E ξ [H α (ξ, x)] σ 2 [y] ξ (x) = α,|α|≤P E ξ [Ψ i (ξ, x)γ i ] 2 -E ξ [y(ξ, x)] 2 (2.89)
Considering this approach one advantage is that the obtained expansion is valid for the full nonprobabilistic variable range and the sensitivity of the statistical moment does not need the sensitivity of the response to be computed. However, by adding another type of variable, the dimension of the PCE is increased which makes its construction more computationally expensive.

PCE for sensitivity study

Global sensitivity analysis Recent work from [START_REF] Sudret | Global sensitivity analysis using polynomial chaos expansions[END_REF] and [START_REF] Crestaux | Polynomial chaos expansion for sensitivity analysis[END_REF] introduces the use of PCE to perform global sensitivity analysis. Two types of methodologies are exposed in the literature.

Spectral representation of stochastic field through Polynomial Chaos Expansion

with the following property:

u⊆{1,2,...,M } u̸ =∅ V u = V. (2.

96)

S u denotes the s-order sensitivity index where s=card(u) and measures the sensitivity of the variance of y due only to the interactions of the variables ξ u . One may count 2 M -1 sensitivity indices. As this number may increase rapidly with the number of variable, one defined the total sensitivity indices S Ti by:

S Ti = u∋i S u .
(2.97)

These express the total sensitivity of the variance of y due to the variable ξ i alone, and all its interaction with the other variables.

To compute the Sobol indices sampling methods may be used but rapidly become unaffordable.

Moreover one may exhibit a straightforward relationship between the polynomial chaos coefficients and the Sobol indices [START_REF] Sudret | Global sensitivity analysis using polynomial chaos expansions[END_REF]. In fact, due to the orthogonality property of the {Ψ} P -1 i=0 the total variance D may be approximated by:

V = R M y 2 (ξ)p(ξ)dξ -ȳ2 = α,|α|≤P γ 2 α E[H 2 α ] (2.98)
The Sobol decomposition of y is approximated by the Sobol decomposition of y PCE . The elements y PCE u of the Sobol decomposition of y PCE may be simply obtained by:

y PCE u (ξ u ) = β∈βu γ β H β (ξ u ) (2.99)
where the selected indices are

β u = {β ∈ {1, . . . , P }, H β (ξ) = s i=1 H α β i (ξ ui ), α β i > 0} (2.100)
The approximated conditional variances may then be expressed as:

Ṽu = β∈βu γ 2 β E[H 2 β ] (2.101)
and thus the sensitivity indices as:

S u ≈ Su = β∈βu γ 2 β E[H 2 β ] P β=1 γ 2 β E[H 2 β ]
(2.102)

PCE for reliability study

As mentioned earlier, the PCE is not of suitable used to approximate distribution tails and thus may be adapted in order to compute low probability of failure. Different strategies have been adopted in order to adapt the PCE to the computation of small probabilites. A first attempt has been proposed by the authors in [Gha91,Gha99b] to sample the orthogonal residuals in order to improve the prediction

Conclusion

accuracy of the distribution tails. As a first step [START_REF] Paffrath | Adapted polynomial chaos expansion for failure detection[END_REF] proposes the use of an adapted PCE to the computation of the failure probability. It defines two methodologies: a "shifted" PCE and a "windowed" PCE inspired from the space discretization in ME-gPC.

An open issue: Extension to correlated non Gaussian random variables

The polynomial chaos expansion extends the random process of interest on a set of polynomials depending on standards random variables. However the observed random variables (denoted ξ obs ) are usually not standards and may be correlated. One may then transform the random variable into standard non correlated random variables and apply the stochastic expansion into the transformed probability space.

Two alternatives are briefly formulated.

One alternative consists in numerically generating optimal orthogonal polynomials and to compute their Gauss points and weights using one of the approaches described in the literature. These approaches produce optimal orthogonal polynomials for an arbitrary probability density function, but a gap in theoretical advances stays wide open concerning the assessment of joint probability distributions of correlated random variables with arbitrary probability density functions.

Another alternative consists in using an non-linear transformation in order to transform non-normal correlated distributions on normal independent distributions. These transformations are usually non linear. Among them one may count, the Nataf transformation, the Rosenblatt transformation, and the Box-Cox transformation.

The input of the Nataf transformation are the marginal density functions and the correlation matrix, it does not need the knowledge of the joint probability density function. The traditional Nataf distribution may be applied if and only if the copula relative to the random variables ξ 1 , . . . , ξ M is normal. However, a generalized Nataf transform has been defined and relax the condition on the copula to the elliptical cases.

The Rosenblatt transformation does not require any condition on the copula type of the random variables but requires the knowledge of the cumulative distribution of conditional random variables, which is usually not provided. The (generalized-)Nataf transformation is more suitable for the common case.

Finally a last alternative has been proposed by [START_REF] Soize | Physical systems with random uncertainties: chaos representations with arbitrary probability measure[END_REF] to generalize the polynomial chaos expansion.

It proposes a theoretical insight to construct multivariate orthogonal basis when the ξ variables are dependent. In this case, L 2 (R M , dp R M ) has no more a tensor product structure.

Conclusion

This chapter offered an overview of the different categories of approaches developed in the literature in order to quantify the impact of uncertainties. A special focus has been provided on the efficient use of the Polynomial Chaos Expansion in Uncertainty Quantification in a probabilistic and non-intrusive framework. Since the preliminary work of Wiener on the homogenous Hermite PCE many developments have been made in order to alleviate some of its initial limitations. (Multi-Element) generalized Polynomial Chaos and arbitrary Polynomial Chaos Expansion mainly addresses the issue of the approximation of non smooth stochastic behaviors exploiting the flexibility in the choice of the polynomial basis. Sparse

Conclusion

approaches based on adaptive strategies investigates the computational cost of the construction of the PCE in high dimensional space aiming at circumventing the curse of dimensionality. Based on the sparsity of effect principle, these strategies aimed at identifying in the expansion terms that are the most significant with regards to the response variability. Finally, an emerging set of methodologies consists in considering the coefficients of the PCE themselves as random variable and to identify them by solving an inverse stochastic problem. These approaches notably allows to quantify the impact of missing data on the accuracy of the produced expansion. Significant contributions within this approach also allows to tackle high dimensional problems. This chapter also points out the limitations of these methodologies opening the way to future work. An open issue has still to be addressed concerning the spectral representation of correlated random variables. A way to deal with this limitation would be to consider the tensor decomposition of high-dimensional stochastic fields [START_REF] Giraldi | Tensor-based methods for numerical homogenization from high-resolution images[END_REF].

Chapter 3

A physics-based metamodel approach for springback variability assessment

In this chapter, we address the stochastic analysis of the U-shaped deep drawing process of a metal sheet using a regression based Polynomial Chaos approach. In this context, very small perturbations of the input variables provide the regression points used to build the PCE. We show that a careful attention has to be paid to the resolution5 of the "high-fidelity" model used to produce those points. In the following, the resolution of the model is of the same order of magnitude than the range of variation of the stochastic variable. We show in this case, that one possible solution to perform the stochastic analysis, consists in defining an intermediate surrogate model. We here consider a physics based metamodel which highlights two advantages. The first one, is to improve the model resolution, and the second one, is to be less computationally costly. This model may then be used in order to build a Polynomial Chaos of the response.

Introduction

In order to accurately assess the springback variability of a formed metal sheet one has to take into account uncertainties linked with the deep drawing process (blank-holder forces, tool radius tolerances, punch speed, etc.) as well as those linked with the physical model parameters (material properties, thickness of the metal sheet, etc.). The assessment of the final springback shape variability relies on the combination of two numerical tools: a numerical model of the deep drawing process and a stochastic process to propagate the uncertainties into this model.

In the context of variability study, the numerical model has to be computed for very small variations of the input parameters. Comparing the model error with the variation magnitude of the output function is of paramount importance to ensure the validity of the variability study. The Finite Element Method (FEM) is a standard tool to construct the numerical model which is used consequently to perform sensitivity analysis on springback of sheet metal forming [1] or to perform reliability analysis (coupled to classical reliability methodology) of the deep drawing process [2]. However, the FEM shows important limitations regarding the variability issues. First of all, to model the highly non-linear phenomena involved (large strains, plasticity, frictional contact, etc.) with a high sufficient accuracy for small variations of the input parameters, the model has to be refined in every direction (decreasing the mesh size while preserving a good enough aspect ratio and increasing the number of integration points). This rapidly leads to unaffordable computational costs, specially when a high number of calls is needed. To circumvent this cost issue, alternative approaches have been proposed: in a sensibility analysis context using an analytical model of the springback [3]; in an optimization context using a "one step" inverse approach for the analysis and optimum process design of deep drawn industrial metal part [4], this model is also used in [5] to train a Moving Least Square surrogate to efficiently perform the optimization. A second limitation is that FEM generates intrinsic errors such as a discretization Int J Mater Form error (due to mesh size) and computational errors (roundoff, quadrature, ...) [6] (page 93). Contact description and through-thickness integration may play a leading role [7].

Common computational methods propagate the uncertainties defined on the random parameters to the shape of the formed metal sheet based on a sampling of the random input parameters (Monte Carlo simulation, Importance Sampling, etc.) and thus require a relatively high number of calls to the underlying computational model. The authors in [3] use Monte Carlo simulations to perform a sensitivity analysis of the springback with regards to process parameters, [2] and [8] use FORM methodology coupled with an enhanced adaptive Monte Carlo methodology to assess the reliability of the deep drawing process. The use of stochastic metamodels combined with advanced sampling methods offers an alternative to the crude sampling methods. The authors in [9] proposes the use of linear and quadratic response surfaces combined with Monte Carlo Simulation to perform the reliability assessment of sheet metal forming process. In [10], a second order Polynomial Chaos Expansion (PCE) is used as a stochastic response surface to assess the reliability of the formed metal sheet with regards to tolerance criteria. Howewer, in these studies, the number of random variables taken into account by the stochastic response surface is limited.

In this paper, we propose an original two-pronged approach to accurately propagate the uncertainties into a "high resolution" model at low computational cost. Our approach is based on a combination of a stochastic surrogate (PCE) and a physics-based Reduced Order Model (ROM). We demonstrate, that in the case of variability study, as the order of magnitude of the variation range decreases the full scale model prediction becomes unstable. We investigate then regularization features of a "high resolution" physics based ROM, allowing us to reach numerical stability for smaller variations of the input variables. Then, a stochastic surrogate model is trained on the full parameter variation range and it is then used to perform the variability study at low computational cost.

In the second section, we introduce the Polynomial Chaos Expansion (PCE) as a stochastic response surface for the uncertainty propagation. We highlight the need of a numerical model characterized by a high precision for small perturbations of the input variables as well as low computational cost. The third section quantifies the limitations of the FEM modeling with regards to the variability problematic by considering "very" small variations of the input parameters around a nominal value. The two following sections introduce the ingredients of the proposed approach: the semi-analytical deep drawing model based on a plain strain Bending-Under-Tension model (denoted semi-analytical B-U-T in the sequel) and the sparse custom version of the PCE (inspired from [11]) for an enhanced accuracy of the variability study at low computational cost. The fourth section introduces the semi-analytical B-U-T and shows how it allows to alleviate the FEM limitations by side stepping the main numerical noise sources (contact description and through-thickness integration). The fifth section discusses a custom sparse version of the PCE, putting into evidence the need to create a sparse polynomial chaos expansion when the number of random variables increases. Finally, the sixth section of the paper illustrates the methodology on the deep drawing process of U-shaped metal sheet with 8 random variables to assess the variability of the springback shape parameters considering a relatively large number of random variables chosen among the most influential ones [1,10].

Assessing the springback variability using polynomial chaos expansion

The purpose of this paragraph is to provide the necessary evidence for the following rather than giving a comprehensive introduction to the PCE which may be found in [12,13].

The PCE is a stochastic metamodel, that is intended to give an intrinsic representation of the stochastic behavior of a function y (scalar random variable) that is defined as a function of an input random vector ξ = {ξ 1 , ξ 2 , . . . , ξ M } with M coordinates and with prescribed probability density function f ξ (ξ ) [START_REF] Wang | Improved two-point function approximations for design optimization[END_REF]. Its accuracy to predict the variability of the scalar output function y highly depends on the underlying "high fidelity model".

Building the multi-variate Hermitian PCE

Assuming that ξ has independent Gaussian components and that y is a second order random variable (E[y 2 ] < ∞), then according to the Cameron Martin theorem, generalized in M dimensions [15] an exact PCE of y is given by

y(ξ ) = ∞ j =0 γ j j (ξ ), (1)
where { j (ξ ), j = 0, 1, . . . , ∞} are the multivariate polynomials orthogonal with respect to the natural inner product

< k , l >= k (ξ ) l (ξ ) M i=1 f ξ i (ξ i )dξ i (2)
with f ξ i (ξ i ) the marginal distribution of the i th component of the random vector ξ . For M independent Gaussian random variables, multivariate Hermitian polynomials exhibit the orthogonality property. Given a degree N , the multi-variate set of M N +1 Int J Mater Form Hermitian polynomials are built using the tensor product of mono-variate Hermitian polynomials. Practically, only a finite P < M N +1 number of terms is retained yielding to the reduced expansion

y(ξ ) ≈ P -1 j =0 γ j j (ξ ), (3)
where

P = (N + M)! N !M! . ( 4
)
In order to properly interpret the sparse PCE algorithm, an appropriate scheme for numbering the multi-variate polynomials of degree d ≤ N needs to be introduced. The scheme used is illustrated in Table 1 on an example of a 3 rd order Hermitian multi-variate PCE with 2 random variables ξ 1 and ξ 2 .

Computing the coefficients of the multi-variate PCE

The PCE accuracy highly depends on the precise computation of its P coefficients. An intrusive Galerkin type approach has been first proposed by [START_REF] Wang | Improved two-point function approximations for design optimization[END_REF]. Projections based non-intrusive methods take benefit from the orthogonal properties of the multivariate polynomials of the expansion. Equivalent stochastic collocation is based on a Lagrangian interpolation in the stochastic space [16]. These methods are the most accurate but also become rapidly unaffordable in relatievly high dimensions.

The regression based approach (emphasized in this paper) consists in finding the best set of PCE coefficients γ = [γ 0 , . . . , γ P -1 ] by minimizing the residual error in the least square sense:

γ = argmin(∥y(ξ ) -(ξ )γ ⊤ ∥ 2 ) (5) yielding to γ = ( (ξ ) (ξ ) ⊤ ) -1 (ξ )y(ξ ) (6)
Table 1 Illustration of the numbering scheme used for a 3rd order PCE with 2 random variables Bi-variate monomials (ξ 1 , ξ 2 ) where y = y(ξ (1) ), . . . , y(ξ

1 ξ 2 ξ 2 2 -1 ξ 3 2 -3ξ 2 ξ 1 ξ 1 ξ 2 ξ 1 ξ 2 2 -1 ξ 2 1 -1 ξ 2 1 -1 ξ 2 ξ 3 1 -
(Q) ) ⊤ , (7)
= ⎡ ⎢ ⎣ 0 (ξ (1) ) . . . P -1 (ξ (1) ) . . . . . .

0 (ξ (Q) ) . . . P -1 (ξ (Q) ) ⎤ ⎥ ⎦ (8)
with ξ (i) , i ∈ {1, ..., Q} representing Q > P samples of the M dimensional random vector ξ . These samples are generated using a Q sized M dimensional standard normal LHS sampling.

The optimal Q number of samples needed to assess the coefficients with a good accuracy is still an open research area, but an empirical rule proposes Q ≥ (M -1) × P [17]. Once the set of P coefficients {γ 0 , ..., γ P -1} has been determined, one may compute the statistical moments of y analytically avoiding Monte Carlo simulations. The first two moments are given by:

E(y) =γ 0 (9) σ 2 (y) = P -1 j =1 E 2 j γ 2 j ( 10
)
For an accurate assessment of the statistical moments, a special attention has to be given to the computations of the P coefficients (9,10). Thus, according to Eq. 6, the underlying model has to be as accurate as possible. Moreover, computing the whole set of coefficients may require a high number of calls Q = (M -1)×P to the underlying model as P (4) increases exponentially with the number of variables M and the degree N .

In the following section we introduce a test case and we quantify the precision requirements for the simulation model. We show that using typical FEM simulations may lead to inaccurate results for a reasonable computational cost when small perturbations are involved.

Inpact of through-thickness integration on the springback prediction

Test case

To illustrate the issue of model variability with regards to small perturbations, we choose to model a 2D deep drawn U-shaped metal sheet. The example used here corresponds to the B3 benchmark test proposed in the conference Numisheet '93 [18] (blank made of mild steel, blank holder force = 2.45 kN). We propose here to use the legacy software ABAQUS v6.10. The model we use is directly based on the implemented example [START_REF] Kang | An efficient response surface method using moving least squares approximation for structural reliability analysis[END_REF].

The Fig. 1 gives the overall geometrical configuration of the deep drawing process. The problem is symmetric (Fig. 1), thus only a half structure is modeled using appropriate symmetry boundary conditions. The blank is modeled using a single row of 175 first-order shell elements (S4R) with Simpson integration rules (7 integration points across the thickness). As the problem is essentially in plane strain state (the width of the blank is 35 mm and its thickness nominal value is 0.8 mm), corresponding boundary conditions are applied on each node. The blank is made of mild steel modeled as an elasticplastic material. Isotropic elasticity and the Hill 48 [START_REF] Ghanem | Adaptative polynomial chaos expansions applied to statistics of extremes in non linear random vibration[END_REF] anisotropic yield criterion for the plasticity are considered. The following values are used: -Young's modulus: 206 GPa, -Poisson's ratio: 0.3, -Yield strength: 167 MPa, -Anisotropic yield criterion: r 00 = 1.79, r 45 = 1.51, r 90 = 2.27.

Int J Mater Form

The tools (punch, blank holder and die) are modeled as rigid body surfaces. The contact occurring during forming phase is modeled using contact pairs. The punch velocity is taken here as 15 m/s and its displacement is s = 70 mm. The blank holder force is defined as F b = 2.45 kN and a mass of 5 kg is attached.

The whole deep drawing process is simulated in two steps. The forming phase is modeled using the common dynamic explicit approach to solve the problem in a reasonable computational time. During this period, the blank holder force is applied with a smooth ramp to minimize the inertia effect and the punch velocity using a triangle step definition starting and ending with 0 velocity and reaching the 15 m/s with the half run.

To investigate the sensitivity limit of this model, we choose to uniformly vary the thickness of the blank from 0.7 mm to 0.9 mm with an 0.002 mm increment size. For each thickness value the springback shape parameters, the curvature ρ, the angles β 1 and β 2 are measured as shown in Fig. 2. The results of the simulations are depicted in Fig. 3 for 7 integration points across the thickness.

The observed responses (Fig. 3(left)) highlight numerical instabilities for small variations (0.04 mm) of the thickness parameter. A deeper insight into the sensibilities (Fig. 3(right)) shows that the model may be used for thickness variations up to an order of x ≈ 10 -5 m which is insufficient to perform a variability study. Moreover, to perform a variability study, an unaffordable number of calls to the fine FEM model may be necessary.

Investigation on the numerical instability on FEM modeling

We model a section of the metal sheet submitted to a typical 2D deep drawing process undergoing bending-unbending loading path and focus on the in-plane stress σ tt distribution across the thickness.

When increasing the number of integration points across the thickness, the in-plane stress σ tt profile through the section reaches numerical convergence shown in figure Fig. 4(left) when the number of integration points through the section is increased from 2 to around 200 (Fig. 4(right)). The convergence is assessed using the following mean square error

Err = I i=1 σ i -σ ref i 2 I i=1 σ ref i 2 , ( 11
)
where I is the number of integration points and σ ref i corresponds to the stress profile obtained with 200 integration points through the section.

This conclusion is in conformity with [7] where adaptive integration was investigated reducing the number of integration points (from 50 to 11). Moreover, among other numerical instability sources one may identify the coarsity Int J Mater Form Fig. 3 Numerical instability for FEM simulations of the deep drawing process of 2-D U-shaped metal sheet under small variations of the thickness (m). On the left side, the obtained responses are depicted, as on the right side the numerical sensitivities of the model with regards to different order of magnitude of thickness variation is plotted of the mesh. To reach a higher accuracy, the mesh of the FEM model has to be refined in every direction (preserving the same aspect ratio). Another limitation to FEM simulations of deep drawing process is contact modeling, also known to induce a non negligible numerical noise and low convergence rates. This leads to the conclusion that using a full scale refined FEM model to perform the variability study becomes rapidly unaffordable.

In the following two sections we successively describe the two ingredients of the proposed two-pronged approach: the semi-analytical B-U-T model (introduced in [21]) to alleviate the main FEM limitations for small variations of input parameters, and the sparse PCE to propagate the uncertainties at low computational cost.

A semi-analytical bending-under-tension model for deep drawing applications

In this section, we describe the first ingredient of the twopronged approach, namely, the physics based reduced order model. We choose a semi-analytical approach introduced in [21] based on a B-U-T model, and highlight how it decreases numerical instabilities while preserving a reasonable precision and low computational costs.

Semi-analytical bending-under-tension model (B-U-T)

The semi-analytical B-U-T model considers the deep drawing process of a 2D U-shaped metal sheet as a 2D plain Int J Mater Form Fig. 5 Finite element description of one slice of the metal sheet strain Bending-Under-Tension (B-U-T) forming process with negligible shear stress (the width of the sheet being sufficiently large). It is based on a semi-analytical numerical approach which combined an analytical approach with finite element modeling taking benefit from material laws independence and avoiding time consuming and low convergence issues such as contact modeling and high number of degrees of freedom. This model may be constructed in three steps [21]:

1. The first step consists in identifying a finite (usually small) number of regions of the metal sheet with an homogeneous loading state in the length direction. For the U-shaped sheet (as shown in the Fig. 5(left)) 5 regions are identified. Knowing the behavior of a single typical section of one of these regions is sufficient to deduce the behavior of the whole region as the loading state is supposed homogeneous in each region. This section (or slice) may be modeled by a handful of 2D/3D solid element or even a single shell element as shown in Fig. 5. 2. The second step consists in defining the loading path that the region is subjected to and to divide it into a sequence of loading states. For each sequence, particular boundary conditions are automatically applied (the reader is invited to refer to [21] for more details on the boundary conditions for each loading state.).

3. Finally, in the last step the whole model is integrated and the springback shape is reconstructed.

Assessment of the numerical noise in the springback shape parameters

Considering the semi-analytical B-U-T model, only very few finite elements are needed to model the deep drawing and springback phases. Moreover, the contact modeling of metal sheet with the punch and/or with the die does not induce numerical noise: the contact is modeled using an analytical uniform pressure applied on the lower or upper side of each slice. In this semi-analytical model no friction is directly applied, but its main effect, the induced tensile force, is taken into account [21]. To highlight the influence of the through-thickness integration, we use the B-U-T model with the configuration described in Table 2 and compare the springback shape parameters obtained for small variations of the thickness. The same range of thickness values as in Fig. 3 is considered. Four nodes plate elements model each slice of the B-U-T model. P 10 , P 50 , P 400 plates use respectively 10, 50, 400 integration points across the sections. As shown in Fig. 6, the through-thickness integration noise is non-negligible for the B-U-T model. For a large range of variation, numerical instabilities may be observed using 10 integration points, and a smooth response may be observed Fig. 7 Evolution of ρ/ (x) with regards to different order of magnitude of thickness variation x, around the nominal value of 0.8 mm using 50 and 400 integration points. When focusing to a smaller range of variation a smooth response is only observed with 400 integration points.

As only one slice is considered at a time, considering as many integrations points through the sections becomes computationally affordable on the contrary to a full FEM. Moreover a variability study is possible for a variation range [10 -4 ; 10 -7 ] m uncovered by the refined FEM model (Fig. 7). We thus retain the B-U-T model to construct the custom PCE surrogate introduced in the next section.

Sparse polynomial chaos expansion approach using least angle regression stagewise algorithm

In this section, the second ingredient of the two-pronged approach, namely, the sparse PCE is described. When using the regression scheme ( 6) in a high dimensional space, the computational cost to compute the full PCE rapidly becomes unaffordable. Strategies to truncate the PCE are needed in order to reduce the number of terms by retaining the most significant ones (Sparsity of effect principle).

Some methods inspired from model selection schemes have been proposed to create a sparse approximation of the PCE [11]. Among them l 1 penalization methods are the most popular. A review of these methods may be found in [22] and also in [START_REF] Donglai | Optimization and tolerance prediction of sheet metal forming process using response surface model[END_REF]. These studies highlight that among the great number of methods, Least Angle Regression Stagewise (LARS) is the most efficient method [24].

We build here a custom algorithm based on the Least Angle Regression Stagewise.

Let us first define the J th -sparse approximation of the output function y by:

ỹJ A (ξ ) = j ∈A γ j j (ξ ) (12)
where A is a sparse index set with card(A) = J and J ≤ P . Then A only contains the set of J indices taken among { j , j ∈ {0, . . . , P -1}}. The J th residual vector is then defined as:

r J = y -ỹJ A (13)
and the corresponding J th correlation vector as

cJ (ỹ) = ⊤ y -ỹJ A (14)
The proposed algorithm proceeds as follows:

Offline phase:

1 Build a stochastic design of experiments. Q samples on the hypercube [0, 1] M are performed using a "space filling" Latin Hypercube Sampling. The corresponding Int J Mater Form normal standard realizations are then retrieved using an iso-probabilist transformation (using the invert multivariate Gaussian PDF):

ξ = {ξ (1) , . . . , ξ (M) } = ⎛ ⎜ ⎜ ⎝ ξ (1) 1 ξ (1) 2 ... ξ (1) M ξ (2) 1 ξ (2) 2 ... ξ (2)
M . . . . . . . . . . . . ξ (Q) 1 ξ (Q) 2 ... ξ (Q) M ⎞ ⎟ ⎟ ⎠
At this stage, ξ is a M-dimensional standard gaussian random variable: ξ ∈ N M (0, 1) represented by Q samples. Finally, the small random variations of the parameters are modeled by scaling each random variable to the desired means µ = [µ 1 , . . . , µ M ] and standard deviations σ = [σ 1 , . . . , σ M ] according to the probabilistic model (see Table 3 as an example). We denote ξ µ,σ , the set of random variables matching the probabilistic models requirements. 2 For each Qsized sample {ξ (1) µ,σ , . . . , ξ (M) µ,σ } the corresponding evaluations of the "high fidelity" model are stored in y (7).

3 Choose an arbitrary N (high) degree PCE composed of P orthogonal multivariate monomials { 0 , P -1 }. In the case where independent random Gaussian variables are considered as, multivariate Hermitian polynomials satisfy the orthogonality property (2) with regards to the measures

f ξ i (ξ i ) = 1 √ 2π e -1 2 ξ 2
i , i = 1, . . . , M. Build the matrix (8) by evaluating each orthogonal multivariate monomial { 0 , P -1 } on the Qsized sampling {ξ (1) , . . . , ξ (M) }.

Online phase:

1 Initialize the coefficients γ 0 , ..., γ P -1 = 0 which sets the current residual r equals to y obtained in the offline phase. 2 Compute the correlation vector c 0 between each element of and y. Retain the predictor ψ j * where j * = arg max| c 0 |. The model then becomes y = γ j ψ j = γ ⊤ (1) ψ. 3 Update γ j to γ * j = γ j + ϵ * (1) , where ϵ * (1) is the LARS step: at this step, another predictor ψ l * has as much correlation with the current residual as does ψ j * (see [24] for numerical computations of ϵ * ). Add ψ * l to the current basis:

y = γ * j ψ * + γ * l ψ * l = γ ⊤ (2) ψ. 4 Update jointly γ (2) = {γ j * , γ l * } following the direction u (2) (∥u (2)
∥ = 1) defined by the joint least-square coefficient on the current residual:

γ * (2) = γ (2) + ϵ * (2) × u (2)
. At this step another predictor ψ k is found to have much correlation with the current residual and is added to the model. 5 Repeat step 4 until m = min(P ; Q) predictors have been entered or until the empirical error (15) has reached the desired threshold.

Each step of the algorithm allows us to add one term to the basis chosen according to its correlation with the current residual. The best obtained model y A * is chosen to give a low enough empirical error

Err emp = E ( fA -f ) 2 σ (f ) 2 . ( 15
)
Once the best model is obtained, the sparse coefficients values are computed according to Eq. 6. If the number Q of available simulations is too low to reach the condition Q ≥ (M -1) × P then the design of experiment may be enriched and the previous algorithm run again. As all computations are analytical, running the algorithm as many times as necessary is considered to have a negligible cost.

The combination of LARS and PCE allows us to propagate the uncertainties through the model at low computational costs. Nevertheless the accuracy of the results highly depends on the accuracy of the "high fidelity" model which has to be accurate for small random variations on the input parameters.

Illustration of the two-pronged B-U-T/sparse PCE approach

In this section, we combine the ingredients of the twopronged approach on the springback variability assessment of a 2D deep drawn U-shaped metal sheet introduced in section "Test case". We demonstrate the validity of the proposed approach firstly when a single random variable (blank initial thickness) is considered (mono-variate case) and then the whole probabilistic model is considered (multi-variate case) in which the standard variation has been arbitrarily fixed at 3 % of the mean value for each variable so that the Int J Mater Form variation of each variable implies a variation in the input parameters with a comparable order of magnitude.

Mono-variate case using full PCE

As an introductory example, we consider the thickness as the unique random normal variable (according to the first line of Table 3) and use Hermitian PCE to assess the stochastic behavior of the springback shape parameters ρ, β 1 and β 2 . As shown in Fig. 8, for each of these responses, the convergence is reached for the mean and standard deviation respectively for a PCE of order 5, 5 and 7. A 7 th order polynomial chaos expansion may be chosen for all the responses. This expansion contains 8 terms and thus Q = 56 calls to the model are needed to compute the whole set of coefficients (6). Additionally, the Fig. 8 also compares the values for the mean and the standard deviation obtained using a post-treatment of the coefficients (circles) (using the equations 9 and 10) and 10 4 Monte Carlo simulations (line) on the PCE surrogate. A good agreement is observed.

The Fig. 9(left) compares the 7 th order PCE surrogate to simulations using the "high-fidelity" B-U-T model around the nominal value of the thickness x = 8.10 -4 m. It shows that the responses most affected by the thickness variation are ρ and β 1 . For these responses, a good fit is observed. Considering the response β 2 , a slight error may be observed specially when comparing the responses values corresponding to low probable values of the thickness. These localized errors on non probable values have a limited influence on the evaluation of the mean and standard deviation as noticed in Fig. 8.

The right side of the Fig. 9 illustrates good agreement of the probability density function obtained with 10 5 Monte Carlo simulations on the PCE metamodel with the one obtained by sampling the B-U-T model directly. In order to illustrate the performance of the approach for an increasing number of variables, we consider now the full set of normal random variables described in Table 3.

In this study, we choose an a priori N = 5 th order PCE with 8 variables yielding to P = 1287 polynomial terms. Q = (M -1) × P ≈ 9000 simulations would be necessary to compute the whole set of coefficients using the collocation scheme (6). We apply the LARS strategy proposed Int J Mater Form Fig. 10 Evolution of the empirical error when increasing the number of PCE terms using LARS algorithm in the previous section, with a stopping criterion on the empirical error. We stop the algorithm when the criterion error reaches 10 -4 . The Fig. 10 exhibits the convergence rate of the empirical error as the terms are added step by step during the LARS procedure. The error threshold 10 -4 is reached for the number of PCE terms given in the second column of the Table 4 for ρ, β 1 , β 2 .

Considering this empirical error, the Fig. 11 illustrates for each response the retained polynomials (numbered according to the procedure described in the first section.)

A more general insight is given in Table 5 which exhibits the number of polynomial terms with regards to their degree in the LARS PCE expansion. We show a posteriori that considering this level on the empirical error, a 4 th order polynomial approximation is enough to assess the springback variability of responses β 1 and β 2 while for ρ, only a sparse 3 rd PCE approximation is needed.

Considering the obtained histograms (using 1, 000 Monte Carlo simulations on the B-U-T model and the PCE model) in Fig. 12 a good agreement for each response is observed. Comparing with the obtained histograms in 11 Representation of the sparsity of the polynomial basis after the LARS algorithm (see section "Sparse polynomial chaos expansion approach using least angle regression stagewise algorithm") has been applied for each response ρ, β 1 , β 2 . A point corresponds to the presence in the basis of one polynomial according to the numbering scheme illustrated in Table 1 Fig. 9(left) (when only the thickness was considered as a random variable), we note that according to our stochastic model, the same variables ρ and β 1 are the most affected by randomness on the input parameters. Moreover, the mean value has moved for the responses 1/ρ, β 1 , β 2 respectively from 4.36 mm -1 , 100.83, 91.74 to 9.6 mm -1 , 112.28, 91.11. Except for β 2 , whose variation range is small compared to ρ and β 1 , considering the full stochastic model highly affects the variability responses: the mean computed in the mono-variate case corresponds to low probable values in the multi-variate case, and the type of the p.d.f for these responses is highly modified. An exponential type distribution is observed in the mono-variate case as with the full stochastic probabilistic model, a nearly symmetric distribution characterizes the responses variability.

Conclusion and prospects

In this paper we combine a semi-analytical Bending-Under-Tension model and a custom Polynomial Chaos Expansion to accurately assess the springback parameters for small variations on the input parameters. The B-U-T model allows us to circumvent typical cost issues and numerical instability from full FEM simulations (contact modeling, through thickness integration). The use of such a model allows us to reach a sufficient numerical stability for small variations of the random parameters. Using these high resolution outputs, we are able to accurately train a custom stochastic surrogate to efficiently propagate the uncertainties through the model. Then, this approach allows us to accurately assess the springback variability when multiple random variables are taken into account with a limited budget. In a more general way, our approach demonstrates that the use of simplified physics based model for large strain forming process allows to reduce the numerical instability and makes possible an accurate and low cost variability study. The approach is of course not limited to 2D plain strain and sparse PCE could be combined with other types of physics-based metamodels such as one-step or POD/PGD approaches could possibly presenting similar smoothing properties in 3D. In the current state of development we have demonstrated the validity of our the approach from numerical point of view using standard benchmarks. The comparison with experiment requires an implementation of the method within a specialized metal forming framework, which is beyond the scope of the current study dedicated to stochastic modeling. Moreover, further work is needed for sensitivity analysis performed using the polynomial chaos expansion to a priori identify the most influent input variables. Adaptive features may also be included in the sparse construction algorithm to keep on reducing the number of terms used, and increasing its prediction accuracy.

Chapter 4

Adapting the sampling to the model resolution: Fat Latin Hypercube sampling.

In this section we propose an alternative methodology to solve the same problematic: the stochastic analysis of the U-shaped deep drawing process of a metal sheet using a regression based Polynomial Chaos approach.

In many cases, the model resolution is not of the same order of magnitude of the range of variation of the random variable. In this case, the use of an intermediate non-physics based metamodel is not mandatory. One may solve the problem by focusing on the design of experiment. We defined a custom Latin Hypercube Sampling (LHS) taking the model resolution into account. A restricted area whose shape is parameterized by the resolution of the model is defined around each sampling point. The resulting LHS consists in sampling points, each characterized by a restricted area free of any other sampling points. This limits the total number of available points. We then propose to build a consistent regression-based PCE with the remaining number of samples points. This way this methodology aligns the stochastic model complexity with the limited number of trustworthy data at hand. This approach has been described into a paper submitted in july 2013 to the journal CMAME and is currently under review.

The paper is given in its submitted version.

Fat Latin Hypercube Sampling and efficient Sparse Polynomial Chaos Expansion for uncertainty propagation on finite precision models -Application to 2D deep drawing process -

Introduction

A profuse literature review reveals rigorous approaches to reduce the computational expense, combining in a hierarchical way a 'high-fidelity" (costly) model with a "lower-fidelity" model (less accurate but also less expensive) to perform parametric studies dedicated either to the search for optimal design [1,2,3] or to the characterization of system variability or finally to sensitivity studies [4,5,6]. In an optimization context, rigorous convergence proofs may be established as long as the "lower-fidelity" model is consistent (generally to the first [7] or second order [8]) with the "higher-fidelity" model. With regards to the nature of the metamodel involved one may classify these approaches as "multi-fidelity" or "variable-fidelity" when physical based surrogate model is involved [9,10,11,12,13] or in surrogate-based approaches when non-physics interpolating or regression-based metamodels built from a Design of Experiment (DoE) are considered. [START_REF] Wang | Improved two-point function approximations for design optimization[END_REF]15,16].

These latter approaches are widely spread in stochastic analysis. For example, [17] uses a combination of Kriging model and subset simulations in order to efficiently assess structural failure probability at low computational costs. [18] uses Monte Carlo simulation and linear response surface to detect the most significant variable and to give an approximation of the probabilistic response. [START_REF] Kang | An efficient response surface method using moving least squares approximation for structural reliability analysis[END_REF] uses the moving least square instead of the classical quadratic order response surface to better fit the limit state function and perform reliability analysis of the sheet metal forming process. [START_REF] Ghanem | Adaptative polynomial chaos expansions applied to statistics of extremes in non linear random vibration[END_REF]21,22,6] demonstrates the advantages of adaptive scheme for Polynomial Chaos Expansion (PCE) to perform robust, reliability, and sensitivity analysis. A second order PCE is used in [START_REF] Donglai | Optimization and tolerance prediction of sheet metal forming process using response surface model[END_REF] to assess the variability of the tolerance prediction of the formed metal sheet submitted to random parameters. In the field of metal forming applications, a classical approach consists of using Monte Carlo simulation on a quadratic polynomial response surface method to quantify probabilistic characteristics (mean and standard deviation), e.g. [24] uses this classical approach to estimate the variability of the shape and dimensional errors in net-shape metal forming. In this approach, the exactness of the statistical quantity mainly depends on the capability of the metamodel of picking up the relationship between the dependent explanatory random variables. The training data may thus play a leading role.

In this work we focus our attention on the quality of the "high-fidelity" training data. We claim that there may exist a threshold on the magnitude of variation of the input variables below which the "high-fidelity" simulations may not be trustworthy. Thus, the number of achievable simulations using "high-fidelity" model is not only limited by the cost of a single simulation but more importantly by its intrinsic sensitivity to small perturbations. In this case, we face two competing issues. On one hand, a too small number of "high-fidelity" simulations harms the accuracy of the response surface. On the other hand, a too high number of simulations introduces numerical noise which also directly leads to a noisy response surface.

In the present paper, we address both issues simultaneously in order to propose a coherent PCE scheme taking into account the upper bound on sampling density given by sensitivity consideration, and the lower bound given by the regression approach to compute PCE coefficients.

The rest of the paper is organized as follows. In section 2, we investigate a sampling strategy taking into account the model precision leading to a non-noisy reduced set of sampling. Then, in section 3 we compare three sparse methodologies applied to PCE to efficiently and accurately propagate the uncertainty with the few number of remaining simulations. Finally, in section 4 we demonstrate the efficiency of the proposed methodology considering the deep drawing process of a 2D-U-shaped metal sheet as a test case.

Sampling schemes taking into account of model sensitivities

In the scope of this work, the FE model is referred to as the "high-fidelity" model and is used to train the analytical "lower-fidelity" model, namely the PCE. Firstly, we show that "very" small random perturbations on the input parameters ξ = [ξ (1) , ξ (2) , . . . , ξ (M) ] around a nominal value ξ nom = [ξ (1) nom , ξ(2) nom , . . . , ξ (M) nom ] yield noisy training data set. We characterize the model output stability by the non-dimensional sensitivity for each considered variable ξ (i) , i ∈ {1, . . . , M} using the finite difference scheme

µ i = ∆y(ξ (i) ) ∆ξ (i) × ξ (i) nom y(ξ (i) nom ) (1)
where

∆y(ξ (i) ) = y ξ (i) nom + ∆ξ(i) 2 -y ξ (i) nom - ∆ξ (i) 2 . ( 2
)
When decreasing the order of magnitude of the perturbation (-log(∆ξ) increasing), the nondimensional sensitivity µ i computed for the "high-fidelity" model exhibits different behaviors:

1. Firstly, for "large" variation of ∆ξ i , the variation µ i reveals the non-linear behavior of the model. No brutal variation of µ i is observed and the model is considered as trustworthy.

2. Secondly, µ i stabilizes around a constant value μi where the model may be considered as linear.

3. Thirdly, on reaching the threshold , µ i becomes unstable. Within the corresponding range of variation noisy data estranged from physical reality is generated. These data points have to be discarded to train a metamodel.

4. Finally, the threshold shows the model sensitivity limit: for this range of variation, the model is not sensitive anymore.

The model output is considered to be unstable when the following criterion is satisfied:

|µ i -μi | > µ * i 2 . (3)
The Fig. 1 illustrates the tendencies obtained for µ i against the order of magnitude of the variation range in ∆ξ i . The results are issued from actual computations shown for a real case (section 4). We thus focus our attention on avoiding noisy training data which may lead to an inaccurate response surface for small variation of the input parameters. We propose to modify the traditional Latin Hypercube Sampling scheme in order to take into account model sensitivities.

2.1. A Fat-Latin Hypercube Sampling taking into account model sensitivities Standard LHS. LHS ( [START_REF] Mckay | A comparison of three methods for selecting values of input variables in the analysis of output from a computer code[END_REF]) is an efficient technique to generate joint probability distributions by distributing samples in equiprobable bins. Let Ξ = ξ (1) , ξ (2) , . . . ,

ξ (M) = ⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ ξ 1 ξ 2 . . . ξ S ⎫ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎬ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎭ = ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ ξ (1) 1 ξ (2) 1 . . . ξ (M) 1 ξ (1) 2 ξ (2) 2 . . . ξ (M) 2 . . . . . . . . . . . . ξ (1) S ξ (2)
S . . . ξ (M) s ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ (4)
an S -sized and M-dimensional DoE and F j , j = 1, . . . , M the joint cumulative distribution (assumed to be known) of each of the random variables. Let Ξ U , the special case where F is an uniform joint distribution for all variables.

A possible 2D sampling for 2 independent uniform distributions laws is given in Fig. 2(a) and for 2 independent normal distributions in Fig. 2(b). By construction there is exactly one observation per row in each of the M directions. The LHS advantages [START_REF] Stein | Large sample properties of simulations using latin hypercube sampling[END_REF] are:

• as long as the number of samples S is large compared to the number of variables M, LHS eventually provides estimators with lower variances for any function with finite variance,

• in any case S -sized LHS does not perform worse than (S -1)-sized crude Monte Carlo.

However, LHS shows also some limitations:

• The error estimates may not be improved by iteratively increasing the number of samples: the resulting sampling is not a LHS anymore (see [START_REF] Qian | Construction of nested space-filling designs[END_REF][START_REF] Qian | Nested latin hypercube designs[END_REF][START_REF] Qian | Nested space-filling designs for computer experiments with two levels of accuracy[END_REF] for Nested LHS).

• When used for training response surfaces, a space-filling optimal design is interesting in order to sample the design space with a minimum number of response evaluations. When using LHS there is a risk that some of the random samples form a cluster to the detriment of some unexplored part of the design space. To circumvent these issues, some strategies may be found [START_REF] Wyss | A user's guide to lhs: Sandia's latin hypercube sampling software[END_REF][START_REF] Dalbey | Fast generation of spacefilling latin hypercube sample designs[END_REF][START_REF] Dalbey | Generating a maximally spaced set of bins to fill for high-dimensional spacefilling latin hypercube sampling[END_REF].

In the following we propose to alleviate another fundamental limitation: when performing stochastic studies small variations of the random input parameters result in noisy responses.

Fat-LHS. In the present paper, we propose to build a restricted area (free of other sampling points) around each sampling point is defined. The shape of this restricted area is parameterized by δ * i (Eq.3) and may be defined as follows:

δ * i = argmin ∆ξ (i) µ i (∆ξ (i) ) > µ * i , i ∈ {1, . . . , M}. (5)
Depending on the chosen norm, different shapes are obtained for the restricted area (Fig. 3). Fig. 3 This approach coupled with LHS requirements permits us to filter the spacial sensitivity but limits the maximum number of samples available for PCE determination. In the remainder of the paper we denote this upper bound on the number of samples at hand by S ub .

Implementation of the Fat-LHS

We propose here an algorithm to identify the maximum number of points to be sampled for a given µ i , i ∈ {1, . . . , M}. We assume that for small variations of the input random parameters the shape and size of the restricted area stays identical for each sampling point. The general idea of the procedure is to start with a given LHS with a prescribed density probability on each random variable and then to:

1. identify the restricted area for each sampled point (Fig. 4 Find and suppress the illegal neighbors

-Let i = 1 while card(s set (i)) > card(s set (i + 1)) do for s = 1 : card(s set (i)) do -Compute: D s = sign ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ ∥ξ (1) s -ξ (1)
1 ∥µ (1) . . . ∥ξ (M) s ξ (M) 1 ∥µ (M) ∥ξ (1) sξ (2) 1 ∥µ (1) . . .

∥ξ (M) s -ξ (M) 2 ∥ -µ (M) . . . . . . . . .
∥ξ (1) sξ (1) S ∥µ (1) . . .

∥ξ (M) s -ξ (M) S ∥ -µ (M) ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ if ∃ j| D (s) ( j, :) = [1, 1, . . . , 1] m times then -s set (i + 1) = s set (i) -{ j} -Ξ {-i} = ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ ξ (1) 1 ξ (2) 1 . . . ξ (M) 1 ξ (1) 2 ξ (2) 2 . . . ξ (M) 2 . . . . . . . . . . . . ξ (1) i-1 ξ (2) i-1 . . . ξ (M) i-1 ξ (1) i+1 ξ (2) i+1 . . . ξ (M) i+1 . . . . . . . . . . . . ξ (1) S ξ (2) S . . . ξ (M) S ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ .
end if end for i=i+1 end while Re-build an LHS sampling -Considering the number of remaining points, adapt the "bin" size to recover the equiprobability property.

-Identify the bins with more than one sampling and the empty bins -Randomly distribute the points from the over occupied bins to the empty bins

Hybrid Q-norm+LARS Polynomial Chaos Expansion (PCE) scheme

In this section we introduce the theoretical aspects for the construction of a sparse PCE decreasing its computational cost and thus allowing us to consider only a limited number of sampling points at hand.

Multivariate PCE

The PCE [START_REF] Ghanem | Stochastic finite Elements: A spectral approach[END_REF] is a stochastic metamodel, that is intended to give an approximation of the stochastic behavior of a functional y (scalar random process) that is defined as a function of an input random vector ξ = {ξ 1 , ξ 2 , . . . , ξ M } with M coordinates. We assume that y is a second order random variable (E(y 2 ) < ∞) and that the probability density function f ξ (ξ) may be decomposed on a product of the marginal probability density functions f ξ (i) (Eq.6):

f (ξ) = M i=1 f ξ i (ξ i ) (6)
Given the natural inner product for arbitrary function φ with respect to each of the marginal probability function f ξ (ξ) defined on

D < φ 1 , φ 2 >= D φ 1 (ξ)φ 2 (ξ) f ξ (ξ)dξ (7)
one may define an infinite set of mono-variate orthogonal polynomials ϕ = {ϕ k , k ∈ N} verifying < ϕ j , ϕ k >= δ jk . Hermitian polynomials respect this condition for Gaussian random variables. For other types of random variables, different orthogonal polynomials may be retained (Table 1) leading to the generalized PCE or Wiener-Askey scheme [START_REF] Xiu | The wiener-askey polynomial chaos for stochastic differential equations[END_REF].

Probabilistic measure Orthogonal Polynomial

Uniform: 1 ]-1,1[ (ξ)/2 Legendre: P k (ξ) Gaussian: 1 √ (2π) e -ξ 2 /2
Hermite: H e k (ξ) Gamma: Using the tensor product on these mono-variate polynomials one may obtain an infinite set of multi-variate polynomials (with a preserved orthogonality property) ψ = {ψ α , α ∈ N M } where α ∈ N M is a multi-index set.

ξ a e -ξ 1 R + (ξ) Laguerre: L k a (ξ) Beta: 1 ]-1;1[ (ξ) (1-ξ) a (1+ξ) b B(a)B(b) Jacobi: J k (ξ)
According to the theorem of Cameron Martin [START_REF] Cameron | Transformations of wiener integrals under translations[END_REF], the exact polynomial expansion of the functional y is

y(ξ) = α∈N M γ α ψ α (ξ). ( 8
)
where {γ α }, α ∈ N M are the coefficients of the PCE to be identified

Truncating multi-variate polynomials expansion

For practical use, one may truncate the full set of tensor product polynomials in order to only retain a finite set of polynomial terms. For Fat-LHS sampling scheme, we need an economical PCE scheme requiring less than S lb ≤ S ub samples. In the following we revisit three classical truncation schemes needed for the construction of our hybrid approach.

Classical truncation scheme. Among all {ψ α , α ∈ N M } the classical truncation scheme [START_REF] Ghanem | Stochastic finite Elements: A spectral approach[END_REF]21,[START_REF] Sudret | Stochastic finite elements methods and reliability-a state of the art[END_REF] retains only the multi-variate polynomial terms whose degree does not exceed an arbitrarily fixed N leading to the following multi-index set:

A M q = {α ∈ N M , ||α|| q ≤ p}, (9)
where

||α|| q = ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ M i=1 α q i ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ 1/q
and q = 1. Then, the truncated model may be written as:

y A M q (ξ) ≈ α∈A M q γ α ψ α (ξ). (10)
Fig. 5 gives an illustration of this scheme for a 2-variate 7 th order PCE. The number P of coefficients in the PCE is given by

P = N k=0 C k M+k+1 = (N + M)! N!M! (11)
and increases exponentially both with N and M. So does the number of "high-fidelity" function evaluations needed to compute the number of PCE coefficients: whatever the method used, at least S = P + 1 samples are necessary.

Q-norm based truncation. This approach relies on the "sparsity of effects principle": a system is usually dominated by main effects and low-order interactions. Q-norm generalizes the classical truncation scheme by varying 0 ≤ q ≤ 1 [START_REF] Blatman | Sparse polynomial chaos expansions and adaptive stochastic finite elements using a regression approach[END_REF]. Fig. 6 illustrates a typical truncated index set different for q = 0.6 and a 7 th order 2-variate PCE. The set of active polynomials in the PCE decomposition is decreased when q decreases. Fig. 7(a) illustrates the evolution of the number of 2-variate polynomial terms in linear scale against q values, and Fig. 7 7: Number of polynomials terms in q-truncated PCE with regards to the q truncation parameter for 2 and 8 variables LARS truncation scheme. The Least Angle Regression Stagewise algorithm [START_REF] Hastie | The Elements of Statistical Learning, Data Mining, Inference and Prediction[END_REF][START_REF] Blatman | Sparse polynomial chaos expansions and adaptive stochastic finite elements using a regression approach[END_REF] is issued from the variable selection community [START_REF] Hesterberg | Least angle and l 1 penalized regression: A review[END_REF]. Roughly speaking, it iteratively adds to the current model the polynomial terms which are the most correlated with the residual response (Algorithm 2).

y A M LARS (k) = α∈A M LARS (k) γ α ψ A M LARS (k) α (12)
where A M LARS (k) is the corresponding multi-index set whose cardinal is k. The current residual response r (k) and correlation vector c(k) 0 are respectively defined by

r (k) = y -y A M LARS (k) (13)
and

c(k) 0 = (y -E[y])(r (k) -E[r (k) ]) (y -E[y]) 2 × (r (k) -E[r (k) ]) 2 . ( 14
)

Algorithm 2 Least Angle Regression Algorithm

Offline phase

• Build a stochastic S -sized and M-dimensional design of experiments according to the stochastic model denoted Ξ.

For each Ssized sample {ξ 1 , . . . , ξ M } store the corresponding evaluations of the "highfidelity" model in y.

• Build a (truncated) polynomial chaos basis Ψ A M q and evaluate it at each sample.

Online phase

• Initialize the coefficients γ α which sets the current residual r equals to y obtained in the offline phase. LARS (k) following the direction u (k) (∥u (k) ∥ = 1) defined by the joint least-square coefficient on the current residual. At this step another predictor ψ θ * is found to have much correlation with the current residual and is added to the model.

• Repeat the previous step until m = min(P; Q) predictors have been entered or until a previously chosen error estimate has reached a minimal value.

Combining Q-norm and LARS

Considering the limited Fat-LHS sampling, we need to find the optimal sparse index set A * such that the error produced on the resulting approximation model y A * is as low as possible. To perform this optimization task, we combine in an iterative manner the Q-norm and LARS truncations. We index by A M q the set of polynomials Ψ q-truncated obtained by a Q-norm. From this set, one may apply the LARS which selects the most correlated polynomial. We thus index by A M q+LARS the sparser set of polynomials Ψ q+k obtained after a Q-norm and k steps using LARS. Fig. 9 illustrates this method by showing a sparse set of active polynomials obtained for q = 0.6 and k = 8.We note that card(A M q+k ) ≤ card(A M q ).

(a) LARS+Q-norm (q = 0.6, k = 8)

Figure 9: Illustration of combined Q-norm + LARS truncation for a 7 th PCE order

We thus face a combinatorial optimization problem:

⎧ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎩ argmin N,q,k Error(A M q+k ) s.t. S lb ≤ S ≤ S ub (15)
where Error is an estimator of the PCE quality that we describe in the section 3.4, and S is the available number of samples. Eq. 15 may be solved using any appropriate algorithm (genetic approach, simulated annealing, etc.). The strategy we adopt here is inspired by [22] and is shown in Algorithm 3.

Algorithm 3 Optimization of Q-norm + LARS parameters

• Arbitrarily choose a set N = {N 1 , . . . , N n } of PCE orders, (N n possibly high).

• Store the response obtained using an Fat-LHS sampling in a vector y = [y 1 , . . . , y N ] for idx n = 1 : n do • Choose a set of significant q values q ∈ {q (i) |i ∈ 1, . . . , Q} for idx q = 1 : Q do -Compute the N th idx n order full polynomial basis Ψ full . -Truncate the full polynomial basis using the q (idx q ) norm giving Ψ q-truncated -Let P remain = card(A M q ) be the number of remaining polynomials after truncation -Perform a V-fold cross validation as follows, with K=2. for idx P = 1 : min(P remain , N) do -Divide the sampling in 2 populations of equal size ξ test and ξ verif for v = 1 : V do -Compute the LARS Algorithm on the P 1 population -Verify the results on the P 2 population by computing the chosen error estimate. end for -Retain the best LARS step k * according to the selected error criterion end for end for end for • Retain the best model with the best N * , q * , k * according to the selected error criterion

Error evaluation of the polynomial expansion

According to Cameron Martin's theorem [START_REF] Cameron | Transformations of wiener integrals under translations[END_REF], when truncating the multi-index set, one may not reach the convergence to the exact solution in the L 2 sense. We assess the results for an another classical error estimate called the corrected Leave-One-Out (LOO) error estimate taking into account the overfitting phenomenon [START_REF] Blatman | Sparse polynomial chaos expansions and adaptive stochastic finite elements using a regression approach[END_REF].

Global error estimate. An estimation of the exact L 2 error is given by the empirical error:

Err emp = 1 N N i=1 (y(ξ i ) -y A M q+k (ξ i )) 2 (16)
However, this estimator is known to under-predict the exact L 2 error: when increasing the complexity of the PCE, the empirical error is systematically reduced, as the exact L 2 error may increase (overfitting phenomenon). By construction, the Leave-One-Out error (LOO error) [22] may be less sensitive to the overfitting. It relies on the computation of the predicted residual

∆ (i) = y(ξ i ) -y A M q+k -i (ξ i ), (17)
for each evaluation ξ i , i ∈ {1, . . . , N}, where y A M q+k -i (ξ i ) denotes the approximated model evaluated in ξ i trained in Ξ/{ξ i }. The LOO error is then computed as

Err LOO = 1 N N i=1 ∆ (i) 2 . (18)
In the general case, the computation of the predicted residuals is a greedy process. In our case, these may be analytically computed

∆ (i) = y(ξ i ) -y A M q+k (ξ i ) 1 -h i (19)
where h i is the i th diagonal term of the

Ψ A M q+k Ψ A M q+k ⊤ Ψ A M q+k -1 Ψ A M q+k ⊤ matrix.
Finally, we compute the absolute LOO as

Err LOO = 1 N N i=1 ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ y A M q+k (ξ i ) -y A M q+k (ξ i ) 1 -h i ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ 2 . ( 20
)
and its relative counterpart

ϵ LOO = Err LOO σ(y) 2 . (21)

Computing the coefficients of the truncated multi-variate PCE

To compute the coefficients of the PCE, intrusive Galerkin type approach has been proposed by [START_REF] Ghanem | Stochastic finite Elements: A spectral approach[END_REF]. Among non intrusive, projection based methods take advantage of the orthogonal properties of the multivariate polynomials of the expansion. Stochastic collocation is based on a Lagrangian interpolation in the stochastic space. It may be proved that this method is equivalent to the former [START_REF] Sudret | Uncertainty propagation and sensitivity analysis in mechanical models -Contributions to structural reliability and stochastic spectral methods, habilitation à diriger des recherches[END_REF].

The regression based approach (on which we focus in this paper) consists in solving the overdetermined system of equations, where each Ξ (i) , i ∈ {1, ..., S } represent S > P samples of the random vector Ξ = {ξ 1 , ..., ξ M }. The optimal number of realizations needed to assess the coefficients with a good accuracy is still an open research area, but [21] proposes an empirical rule (Eq.22). In the following, we consider this empirical requirement as a lower bound in the simulation requirements (S ≥ S lb ) to build a PCE.

S lb = (M -1) × P.

The set of coefficients may be computed as

γ = argmin(∥y(ξ) -Ψ(ξ)γ ⊤ ∥ 2 ) ( 23
) yielding γ = (Ψ(ξ)Ψ(ξ) ⊤ ) -1 Ψ(ξ)y(ξ) (24)
with ξ (i) , i ∈ {1, ..., S } representing S > P samples of the M dimensional random vector ξ. These samples are generated using a S -sized M-dimensional Fat-LHS.

The overdetermined system of equations (Eq.24) may be solved using a singular value decomposition of the ΨΨ ⊤ matrix.

Once the set of P coefficients {γ α } α∈A has been determined, one may compute the statistical moments of y analytically avoiding Monte Carlo simulations. The first two moments are given by: E(y) =γ 0 (25)

σ 2 (y) = α∈A-{0} E(Ψ 2 α )γ 2 α ( 26
)
4. Results and discussions

Experiment description

We model the variability of the springback parameter of a 2D deep drawn U-shaped metal sheet from the B3 Numisheet '93 [START_REF]Numisheet'93 2nd International Conference, Numerical Simulatoin of 3-D Sheet Metal Forming Process -Verification of Simulation with Experiment[END_REF] benchmark.

The overall geometrical configuration of the deep drawing process is illustrated in Fig. 10. The process is modeled using a legacy software [START_REF] Lstc | LS-Dyna Keyword User's Manual[END_REF] (Fig. 10) using appropriate symmetry boundary conditions. A single row of 175 first-order shell elements is used to model the blank with Simpson integration rule with 10 integration points across the thickness. As the problem is essentially in plane strain state (the width of the blank is 35 mm and its thickness nominal value is 0.8 mm), corresponding boundary conditions are applied on each node. The blank is made of mild steel modeled as an elastic-plastic material. Isotropic elasticity and the Swift isotropic hardening law are considered σ = K 0 (ϵ 0 + ϵ p ) n 0 .

Die

Figure 11: Representation of the Swift hardening law for the parameters described in Table 2 The value of the geometrical, material, loading and contact parameters are summarized in The tools (punch, blank holder and die) are modeled as rigid body surfaces. The punch velocity is taken here as 15 m/s and its displacement is s = 70 mm. The blank holder force is defined as F b = 2.45 kN. The whole deep drawing process is simulated in two steps.

1. The forming phase is modeled using the explicit dynamic approach to solve the problem in a reasonable computational time. During this period, the blank holder force is applied with a smooth ramp to minimize the inertia effect and the punch velocity using a triangle step definition starting and ending with 0 velocity and reaching the 15 m/s with the half run. The contact occurring during forming phase is modeled using contact pairs.

2. The springback phase is modeled using an implicit approach. At the end at this phase, the springback shape parameters (output functions of interest), the curvature ρ, the angles β 1 and β 2 are measured as shown in Fig. 12.

Sensitivity analysis

Results of the sensitivity according to the method presented in section 2 are given in Table 3.

Stochastic model

In Table 4, we identify the set of independent random variables considered in the model. If the variation range of the parameters may be considered as realistic, the Gaussian hypothesis is only illustrative.

The mean values correspond to the nominal values, and the standard deviations are adjusted so that ξ min = ξ mean -3σ and ξ max = ξ mean + 3σ.

2D validation test case

In this paragraph, we consider only 2 random variables: the thickness of the blank and the Young modulus. If we consider the identified sensitivity in section 4.2, the sampling methodology developed in section 2.1 allows us to consider S ub = 343 samples respecting the sensitivity criterion Eq.3. The Fig. 15 is obtained using the classical LHS and Fat-LHS. For the same number of samples a non-negligible bias in the mean value is observed when the number of sampling becomes high, the classical LHS mean and standard deviation converges to the values produced by the Fat-LHS for a smaller number of sampled points.

. Comparison of sparse PCE strategies

We here compare the convergence results of the different truncation strategies for different 2variate PCE of increasing order. We consider a limited number of 457 simulations issued from the Fat-LHS previously described.

We choose a polynomial order N ∈ {1, . . . , n} for the PCE approximation. For each PCE order N we apply 3 truncation strategies:

1. We select the Q-norm parameters such that S lb ≤ S ≤ S ub . For all the possible q parameters we compute the approximate model y A M q , and retain the one with the lowest LOO error. 2. We use the LARS based algorithm on the classical truncation scheme and retain the best model.

3. We combine the Q-norm and LARS approach and retain the best approximate model.

Fig. 16 shows the evolution of the LOO corrected error with regards to the number of terms contained in the best PCE approximation. Each point refers to the best model obtained during the truncation for different PCE order. LARS alone provides the worst results, most of the time less accurate and more costly than the two other methods. Comparing the Q-norm and the combined LARS+Q-norm we observe similar results in terms of accuracy. However, a slight advantage may be given to the Combined LARS+Q-norm as it gives similar accuracy for a sparse PCE expansion. In addition, we note that due to smooth training data, we finally obtain the best results for truncated low order PCE. The histograms (Fig. 17) illustrate the variability obtained for the best retained model for each response obtained using the Algorithm 3. A good agreement is observed with the exact "highfidelity" simulations as the relative error in mean is close to 0 and in standard deviation lower than 1% as depicted in table 5.

a) ρ, N = 3 (b) ρ, N = 4 (c) ρ, N = 5 (d) ρ, N = 6 (e) β 1 , N = 3 (f) β 1 , N = 4 (g) β 1 , N = 5 (h) β 1 , N = 6 (i) β 2 , N = 3 (j) β 2 , N = 4 (k) β 2 , N = 5 (l) β 2 , N = 6

Conclusions and prospects

In this paper, we highlighted a fundamental limitation of the surrogate-based approach for uncertainty propagation. We showed that when using non linear FE scheme as "high-fidelity" simulation, small variations of the random input parameters may lead to noisy input training which alters the accuracy of the training data set and may distort the statistical quantities of interest. We illustrate this claim using the non-linear FEM simulation (involving contact/friction and material non linearities) of the springback of a 2D deep drawing process of U shaped metal sheet. We propose a sampling methodology called Fat-LHS allowing us to filter noisy data preserving their LHS property. This heuristic strategy provides the maximum number of simulations available considering the finite model sensitivity. We then use this limited number of non-noisy samples to build a PCE in order to propagate the uncertainty. But, the low number of samples leads us to consider sparse strategies to make affordable possible identification of the PCE terms. We compare three different methodologies to build a sparse PCE (LARS, Q-norm and LARS+Q-norm) and retain the best possible PCE for each of them. The comparison of the results shows that generally the Q-norm+LARS hybrid is more efficient. We obtain the best results for truncated low order sparse PCE leading to unbiased estimation. Further work is required to

• economically compute the model sensitivity threshold,

• generate more space-filling LHS design (in this paper, only the min-max strategy has been tested)

Finally, to assess the PCE accuracy, we use a LOO corrected PCE error in order to

• assess the goodness of fit of the PCE

• and to limit the overfilling phenomenon simultaneously.

Thus it is difficult to separate which part of the inaccuracy comes from the model misspecification and which part comes from the overfitting phenomenon. A formulation of an overfitting measure for PCE approximation inspired from [START_REF] Bilger | Measuring overfitting and mispecification in nonlinear models[END_REF] may open a new way to efficiently select the most significant polynomials terms in a sparse PCE expansion.

Chapter 5

Towards multi-objective optimization under uncertainty

As a natural step forward, incorporating the treatment of uncertainties into an optimization process may provide the engineer with additional valuable information concerning for example, the robustness and the reliability performances of the system to design. Most of the time such studies are made in a single objective context. However, in many real-life cases, the performance of a system is defined by multiple criteria, some of them competing.

In this section we thus propose an efficient metamodel-based strategy in order to tackle the issue of multi-objective optimization under uncertainty. Most of the time, evolutionary algorithms offer a reasonable way to solve this multiobjective optimization problem. However, these populations based algorithms are known to require a large number of evaluation of the numerical models. Inserting the treatment of uncertainties as an inner loop of the multiobjective optimization process may make the process become unaffordable. Moreover, the concepts defined in single objective optimization are not directly appliable in multiobjective optimization as they usually do not take into account the intrinsic multiobjective nature of the problem.

Both issues are addressed in this section. An original hierarchical metamodel approach is provided in order to perform a multiobjective optimization task at affordable computational costs. It is then applied in an original formulation to take uncertainties into account in a multiobjective context as well as an efficient approach to solve it based on a hierarchical approach.

The following paper has been published in the journal Structural Multidisciplinary Optimization.

As a second author in this paper, I particularly focus my attention on the definition of the hierarchical metamodel strategy for uncertainty propagation. (x) J x ( )

J x ( ) = i w i (x (i) , x) T (x (i) ) -γ (i) , (x)
x (i) w i x (i) x

w i = w x (i) -x) r , w w (s) = -s + s ≤ s ≤ , s ≥ ,
f i = f i -f i f i -f i f i f i f i 4.
M N N N N R k Q > Q P Q = k Q P k Q = N A = × × = . N B = × × × = N C = α = α = b ξ i ∂b ∂ξ i = j ∂b ∂ψ j ∂ψ j ∂ξ i , ψ j ξ i ψ j (ξ ) = M l= H κ l ( j) (ξ l ),

Conclusions

In this thesis we investigated the following question in the scope of variability studies when small variations are concerned: how far may we trust the "high-fidelity" models? We have shown on a classical metal forming test problem than when using the non-linear FE scheme as "high-fidelity" simulations, small variations of the random input parameters lead to noisy training data and may distort the statistical data of interest.

When using non linear FE scheme as "high-fidelity" simulation, small variations of the random input parameters may lead to noisy input training data which alters the accuracy of the training data set and may distort the statistical quantities of interest. We have illustrated this claim using the non-linear FEM simulation (involving contact/friction and material non linearities) of the springback of a 2D deep drawing process of U shaped metal sheet.

We have introduced two methodologies in order to bypass this limitation. The first methodology consists in combining a physical reduced order metamodel (semi-analytical Bending-Under-Tension model) and a custom Polynomial Chaos Expansion to accurately assess the springback parameters for small variations on the input parameters. The B-U-T model allowed us to circumvent typical cost issues and numerical instability from full FEM simulations (contact modeling, through thickness integration). The use of such a model has allowed us to reach a sufficient numerical stability for small variations of the random parameters. Using these high resolution outputs, we have been able to accurately train a custom stochastic surrogate to efficiently propagate the uncertainties through the model. Then, this approach allowed us to accurately assess the springback variability when multiple random variables are taken into account with a limited budget. However, physics-based metamodels with higher resolution are not always available. To circumvent this issue, one opportunity is to take into account the model sensitivity in the sampling scheme. We proposed a modified Latin Hypercube Sampling methodology called Fat-LHS allowing to filter noisy data and to preserve their Latin Hypercube Sampling property. This heuristic strategy provides the maximum number of simulations available considering the finite model sensitivity. We then used this limited number of non-noisy samples to build a PCE in order to propagate the uncertainty. But, the low number of samples has lead us to consider sparse strategies to make affordable possible identification of the PCE terms. We have compared three different methodologies to build a sparse PCE (LARS, Q-norm and LARS+Q-norm) and have retained the best possible PCE for each of them. The comparison of the results has shown that generally the Q-norm+LARS hybrid is more efficient. We have obtained the best results for truncated low order sparse PCE.

Moreover, the combination of optimization and uncertainty quantification has been addressed. In most of the studies, these both fields have been combined when only one objective function is considered. However, most of real life problem involved the simultaneous optimization of possibly competing objectives. Then two issues have to be tackled: the definition of a formulation of MOOU, and its resolution. In this thesis, we particularly focused on the resolution phase. Most of the time this task is achieved using a nested combination of metamodels constructed separately either in the deterministic design variables space or in the stochastic variables space. Here, we have proposed a computationally efficient non-intrusive procedure to construct a metamodel in both deterministic and stochastic spaces. It is based on a moving least squares interpolation of polynomial chaos coefficients. Those metamodels have 6.2. Prospects been thoroughly validated and incorporated into an original multi-objective reliability-based formulation, and successfully applied to the cost-effective structural design optimization of space trusses.

Prospects

To deal with real-life multi-objective optimization problem under uncertainty, one has to take into account the eventually high dimensionality of the problem and to develop robust methodologies to this purpose.

The examples treated in this work are characterized by small numbers of variables and rather smooth responses. For more complex problems, higher degrees of PCE and MLS bases might be required, which would considerably increase the size of the training database, hence the overall CPU cost. In that case, considering sparse PCE to discard from the PCE coefficient matrix the terms that bring no significant contribution to the random response might be an option, in addition, to a finer updating process during the optimization; Moreover, tackling wider design spaces, improving the adaptive strategy of the stochastic metamodels should also be considered, for example by including a prediction error procedure to validate the surrogate model accuracy during the optimization. Another approach could imply screening or space reduction techniques in a MOOU framework.

Finally, focusing on the UQ part of the problem, our approach demonstrates that the use of simplified physics based model for large strain forming process allows to reduce the numerical instability and makes possible an accurate and low cost variability study. The approach is of course not limited to 2D plain strain and sparse PCE opens the way to a combination with other types of physics-based metamodels such as one-step or POD/PGD approaches presenting similar smoothing properties in 3D.

The consistency between the different fidelities levels in the approach have not been addressed as [START_REF] Ng | Multifidelity uncertainty quantification using nonintrusive polynomial chaos and stochastic collocation[END_REF] recently investigates this issue. In the current state of development we have demonstrated the validity of our the approach from numerical point of view using standard benchmarks. The comparison with experiment requires an implementation of the method within a specialized metal forming framework, which is beyond the scope of the current study dedicated to stochastic modeling, but may be done in further studies.

Figure 1 . 1 :

 11 Figure 1.1: A possible illustration for uncertainty sources

Figure 1 .Figure 1 . 3 :

 113 Figure 1.3 plots a straightforward manner to perform optimization under uncertainty in a non intrusive context.

Figure 1 . 4 :

 14 Figure 1.4: Typical sensitivity results issued from actual computation

Fig. 1 .

 1 5 illustrates Pareto fronts for two different values of α. The best performing solutions are obtained for a minimization of the functions y 1 and y 2 .

Figure 1 . 5 :

 15 Figure 1.5: Illustration of dominated and non dominated solutions. The red vector of objectives is dominated by 3 green vectors. The green vectors are not dominated and form the Pareto front.

  (a) Geometrical configuration of the modeled Numisheet'93 benchmark (b) Illustration of springback shape [Car09]

Figure 1 . 6 :

 16 Figure 1.6: Springback phenomenon for a 2D deep drawn U-shaped metal sheet

2 Figure 1 . 7 :

 217 Figure 1.7: Sensitivity restricted area shape around two sampling points for L ∞ and L 2 in 2D

1. 3 .

 3 Contributions of the thesis Instead of comparing each objective function separately, one may take into account of the intrinsic multi-objective nature of the problem by checking the probabilistic non-dominance ([Tei01]) P non-dominance ≡ P [y(x, ξ) ≻ ζ] ≥ α assessing that the vector of the objective functions y(x, ξ) = [y 1 (x, ξ), y 2 (x, ξ)] should dominate the set of quantiles ζ with a minimum probability level (user defined) α while satisfying constraints (such as reliability for example). The Fig.1.8 illustrates Pareto fronts for two different values of α. The best performing solutions are obtained for a minimization of the functions y 1 and y 2 . When the values of α increases, the stochastic Pareto front represented by the non dominated quantiles ζ moves farther from the deterministic one being more and more conservative.

Figure 1 . 8 :

 18 Figure 1.8: Stochastic Pareto fronts

  Chapters 3, 4 and 5 refer either to our published articles (Chapter 4, Chapter 5) or to papers under review (Chapter 3). Chapter 3 and Chapter 4 rise the question of the trust to give to the nominal model in the context of uncertainty propagation for variability study. Answers are proposed in the field of metal forming where highly non linear phenomena such as material non linearities, contact friction, springback, . . . occur, and which makes in practice the stability of the numerical model difficult to control. In fact, when considering the springback prediciton of the 2D U-shaped deep drawn metal sheet, a large variety of results is observed [MNOW93]

  Chapter 3 addresses this issue by presenting an original two-pronged approach based on the combination of a semi-analytical model dedicated to plain strain deep drawing based on a Bending-Under-Tension numerical model (B-U-T model) to accurately predict the influence of small random perturbations around a nominal solution estimated with a full scale Finite Element Model (FEM). A custom sparse variant of the Polynomial Chaos Expansion (PCE) is used to model the propagation of uncertainties through this model at low computational cost. In Chapter 4, a particular attention is given to the definition of an adapted Design of Experiment (DoE) taking the model sensitivity into account which limits the number of sampling points. The construction of an adaptive sparse PCE based on the limited set of data is investigated.

Figure 2 . 1 :

 21 Figure 2.1: Illustration of FORM/SORM methodologies

1 √S

 1 . For example, for P failure = 10 -4 , and a targeted C V = 10 -1 , 10 6 samples may be needed. The Importance Sampling methodologies[AB99][ER93] allow to improve the efficiency of the crude MCS by centering the samples around the conception point following a translated probability density function p translated = p ξ (ξξ * ). Directional sampling [Bje88], [Mel94, NE00], line sampling [KPS04],

Figure 2 . 2 :

 22 Figure 2.2: Illustration of the classical truncation scheme for a 7 th order PCE

  WK06c, WK06b] introduced Multi-Element generalized Polynomial Chaos Expansion (MEgPC) as an extension of the gPC dedicated to the assessment of non continuous probability density function. It is based on a decomposition of the random space into local elements, and subsequently implements gPC locally within the individual elements. An error control has been developed for the ME-gPC in [WK09].

Figure 2 . 3 :

 23 Figure 2.3: Illustration of Q-norm truncation with different values of the truncation parameter q for a 7 th order PCE

Fig. 2 .Figure 2 . 4 :

 224 Fig.2.4(b) shows it for an 8-variate polynomial in log scale.

2. 2 . 6 . 2 40 2. 2 .

 262402 Adressing the curse of dimensionality using sparse approaches Non-intrusive spectral projection The projection method uses the orthogonality property of the polynomial basis elements on the Hilbert space of finite variance random variables, {H α (ξ)} |α|≤P to compute the coefficients of the expansion {γ α } |α|≤P : Spectral representation of stochastic field through Polynomial Chaos Expansion

  2.80) p(γ|y) ∝ p(y|γ)p γ (γ) (2.81) where p(γ|y) is the posterior distribution of the parameters values (obtained after having observed the data), p(y|γ) = L(γ) is the likelihood function and p γ (γ) is the prior probability density function (obtained before having observed the data) of the parameters values. The full characterization of the posterior probability density function allows to investigate the impact of data limitations associated with the calibration of parameters on the overall predictive accuracy of the 2.2. Spectral representation of stochastic field through Polynomial Chaos Expansion PCE. Moreover, the Gaussian variables of the PC basis may be seen as the part of the approximation modeling the aleatory uncertainty, as the random PC coefficients represent the epistemic uncertainty.

Fig. 1

 1 Fig. 1 Geometrical configuration of the modeled Numisheet'93 benchmark (values in mm)

Fig. 2

 2 Fig. 2 Definition of springback parameters, ρ, β 1 and β 2

Fig. 4

 4 Fig.4 Evolution of the stress across the section for different number of integration points (thickness 0.8 mm)

Fig. 6

 6 Fig.6 Evolution of ρ, β 1 and β 2 with regards to small thickness variation for different numbers of through-thickness integration points

Fig. 8 Fig. 9

 89 Fig. 8 Convergence in Mean and Standard deviation with regards to the polynomial degree for ρ, β 1 , β 2

Figure 1 :

 1 Figure 1: Typical sensitivity results issued from actual computation 3

Figure 2 :

 2 Figure 2: LHS sampling illustration for uniform distribution

  (a) describes L ∞ related restricted area, the lengths of the border of the hypercube being δ * i , i ∈ {1, . . . , M}. Fig.3(b) defines an elliptic sensitivity region restricted area, where the δ * i , i ∈ {1, . . . , M} define the lengths of the axes.

2 Figure 3 :

 23 Figure 3: Sensitivity restricted area shape around two sampling points for L ∞ and L 2 in 2D

  (a)), 2. iteratively remove the illegal neighbors (Fig.4(b)) 3. re-adapt the bin's size to recover the equiprobability property (this makes the sampling loosing in empty bins and over-occupied bins (Fig.4(c))), 4. reconstruct an LHS sampling by randomly allocating the samples from over-occupied bins to the empty bins (Fig.4(d)).

Figure 4 :Algorithm 1

 41 Figure 4: Illustration of the Fat Latin Hypercube Sampling procedure

Figure 5 :

 5 Figure 5: Illustration of the classical truncation scheme for a 7 th order PCE

Figure 6 :

 6 Figure 6: Illustration of Q-norm truncation with different values of the truncation parameter q for a 7 th order PCE

Fig. 8 24 Figure 8 :

 8248 Fig.8illustrates truncated index obtained after k = 24 iterations of the LARS algorithm applied to a 2-variate 7 th order polynomial chaos expansion.

Figure 10 :

 10 Figure 10: Geometrical configuration of the modeled Numisheet'93 benchmark (values in mm).

Figure

  Figure 12: Definition of springback parameters, ρ, β 1 and β 2 .

Figure 13 :

 13 Figure 13: Numerical instability for FEM simulations of the deep drawing process of 2-D U-shaped metal sheet under small variations of the Thickness and Young modulus (MPa) for the ρ response. On the left side, the obtained responses are depicted, as on the right side the numerical sensitivities of the model with regards to different order of magnitude of thickness variation is plotted.

Figure 14 :

 14 Figure 14: Illustration of LHS samplings obtained with taking (a) and not taking (b) into account the sensitivity constraint for 343 samplings. In red appears the illegal neighbors.

Figure 15 :

 15 Figure 15: Comparison of the mean and the standard deviation evolution for different sampling size of classical LHS with the Fat-LHS

(

Figure 16 :

 16 Figure 16: Comparison of the evolution of the LOO error corrected for different polynomial degree with regards to springback parameter ρ (a), β 1 (b), β 2 (c). The parameters considered as random variable are depicted in table 4

H 3 Table 5 : 5 Figure 17 :

 35517 Figure 17: Histograms obtained for the modified cross validation strategy for the 8-variate case using the ϵ LOO corrected

  γ s ) s {γ s , . . . } ξ H j (ξ ) j H (ξ ) = ; H (ξ ) = ξ ; H (ξ ) = ξ -; H (ξ ) = ξξ ; H (ξ ) = ξξ + . , . . . , ξ M ) = P-j= γ s j ψ j (ξ ), ξ , . . . , ξ M ) = M k= H κ k (ξ k ), {κ , . . . , κ M } κ k ≥ ∀k ( ) γ s , {ξ ( ) , . . . , ξ (Q) } Q ( ) = [s(ξ ( ) ) s(ξ (Q) )] T

  f . . . γ g p ] P × (m + p) P ×Q = . . . p [Q×(m+ p)] ξ s( , ξ ) = γ s ( ) + P-j= γ s j ( )ψ j (ξ ). {ξ ( ) , . . . , ξ (Q) } ( ) = ( ) ( , ). γ s j ------------→ ( ) ----------→ ( , ξ ). ) ≈ γ (x) = T (x) (x), T (x) = [ x x . . . ].

∈ζ≡

  f , σ f , . . . , µ f m , σ f m T . Robust Multi-Objective Genetic Algorithm Concurrent Subspace Optimization * multiobjective robust solution of type I δ B( , δ) |B( , δ)| ⎧ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎩ f i ( ) ≡ |B( , δ)| ′ ∈B( ,δ) f i ( ′ )d ′ , i = , . . . , m f (µ ), . . . , µ f m (µ ) , P g j ( , , ) ≥ ≥ R j , j = , . . . , J, h k ( ) ≥ , k = , . . . , K , (L) ≤ µ ≤ (U ) , (L) ≤ ≤ (U ) , = ζ , . . . , ζ m T P ≡ P ( , ξ ) ≻ ζ ≥ α , P ≡ P ( , ξ ) ≤ ≥ α , ⎪ ⎪ ⎪ ⎩ ,η ζ = ζ , . . . , ζ m T P ≡ P ( , ξ ) ≻ ζ ≥ α , P ≡ P ( , ξ ) ≤ ≥ α , ζ i = µ f i ( , ξ ) + η i .σ f i ( , ξ ) , µ f i ( , ξ ) σ f i ( , ξ ) ∀i ∈ I = { , . . . , m} : f A i P ( , ξ ) ≻ ζ = P { f ( , ξ ) ≤ ζ } ∧ • • • ∧ { f m ( , ξ ) ≤ ζ m } = P { f ( , ξ )ζ ≤ } ∧ . . . = P i= ,...,m { f i ( , ξ )ζ i } ≤ . P ≡ P ( , ξ ) ≤ = P {g ( , ξ ) ≤ } ∧ • • • ∧ g p ( , ξ ) ≤ = P j= ,..., p g j ( , ξ ) ≤ . {. . . } C ξ g j ≤ , j = , . . . , p → g i ( , ξ )ζ i ) g ( , ξ ) = ρ

F

  ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ ,η ζ i ≡ γ f i ( ) + η i ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ ,η γ f i ( ) + η i ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ µ i ( , ξ ) = , ξ ) = ⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ b ( , ξ ) ξ )ỹ( , ξ ) direct approximation of the statistical outputs E σ Ẽy ( , ξ ) σy ( , ξ ) = F ( + ξ ), F = F ( + ξ ),

6 General

 6 =i H κ l ( j) (ξ l ) ⎞ ⎠ ∂ H κ i ( j) ∂ξ i . κ i ( j) = ∂ H ∂ξ i =i H κ l ( j) (ξ l ) ⎞ ⎠ κ i ( j)H κ i ( j)-(ξ i ) . =k H κ l ( j) (ξ l ) ⎤ ⎦ κ k ( j)H κ k ( j)-(ξ k ) =k H κ l ( j) (ξ l ) ⎤ ⎦ κ k ( j)H κ k ( j)-(ξ k ) ( j)H κ k ( j)-) κ k ( j) = Chapter

  = {ξ 1 , . . . , ξ M }. Let us assume now that ξ has been evaluated in S samples denoted ξ(1) , . . . , ξ(S)

	2.1. Introduction and general concepts
	when experiments are available. It also allows to develop non-intrusive approaches by defining the
	mapping function h independently from the physical model.
	2.1.2.1 An overview of direct approaches
	Maximum Entropy Principle The Maximum Entropy Principle [Jay82, Gif08] is based on the defini-
	tion of a unique measure which quantifies the amount of uncertainties represented by a discrete probability
	density function, namely the entropy measure. It follows the natural intuition that a "broad" probability
	distribution represents more uncertainty than a "sharp" one. Let ξ be any component of the random vec-
	tor ξ

Table 2 .

 2 may be built over the tensor product of optimal orthogonal polynomials with regards to the probability measure for each ξ i , i = 1, . . . , M ξ and by Legendre polynomial (orthogonal with respect to the uniform distribution) for the x i , i = 1, . . . , M x variables[START_REF] Ms Eldred | Recent advances in non-intrusive polynomial chaos and stochastic collocation methods for uncertainty analysis and design[END_REF]. The cost of the latter expansion may become rapidly prohibitive as it increases exponentially with the number M ξ + M x of variables (see

	γ α (x)Ψ α (ξ)	(2.49)
	α∈N M

2: Bi-variate Hermitian Polynomials Basis elements up to the order p = 3 Note: One may imagine different ways to classify the obtained multi-variate monomials. When a M dimensional expansion of order N is searched, one may classify the multi-variate monomial by making the value of the M -digits of α evolving according to the "rain drop" algorithm. The previous tabular provides an example with M = 2 and N = 3.

2.2.3.3 PCE for stochastic field representation

PCE for random field The extension to the approximation of random process as may introduce an additional non probabilistic variables such as design variable for example gives the following expression: 2.2. Spectral representation of stochastic field through Polynomial Chaos Expansion y(x, ξ) ≈ Let y(x, ω) : D × Ω → R be a mean-square continuous random field (i.e. with continuous covariance function) where D denotes an open, bounded domain with values in R d , and Ω is the set of elementary events.

Let us give a set of points x 1 , . . . , x N ∈ R d and ω 1 , . . . , ω S a set of S elementary events. Let y

  According to the Cameron Martin theorem, the PCE converges in a mean square sense for L 2 random

	2.2.4 On the limitation of Wiener Hermite PCE for non Gaussian probability
	measure
	2.2.4.1 Convergence for non Gaussian random variable

variables and stationary stochastic processes. However, as only a few terms of the PCE are retained in many applications this statement is sometimes irrelevant.

[START_REF] Kennedy | Solution to the practical problem of moments using nonclassical orthogonal polynomials, with applications for probabilistic analysis[END_REF]

particularly shows on analytical random non Gaussian variables that when increasing the number of terms of the PCE, some metrics are systematically improved: among them, the mean square accuracy, the relative error accuracy, the probability of the PCE approximation of taking non physical values, the probability that the PCE is equal in probability to the true function, the difference between the tails of the true function and the PCE to the upper 1% fractile, the relative errors in variance and kurtosis. When increasing the number of terms all the considered metrics are improved in a monotonic way with different convergence rate, some of them being highlighting particularly slow convergence rates.

  66)2.2. Spectral representation of stochastic field through Polynomial Chaos Expansionwith m i ∈ N. The Gauss quadrature method allows to compute multivariate integrals by transforming it into a weighted sum. In the M -dimension case such a weighted sum is written:

			K		K
	R M	y(ξ)Ψ j (ξ)Φ(ξ)dξ =	k1	. . .	kM	w k1 . . . w kM y(ξ k1 , . . . , ξ kM )	(2.67)

Table 2

 2

	Geometrical, material, loading and contact parameter	Geometry		Material		Loading		Contact
	of the U-shaped B-U-T model
		L s	300 mm	E	70.5 GPa	F b	300 kN	µ	0.15
		h 0	0.8 mm	ν	0.342	s	60 mm
		W s	1 mm	ρ	2700 kg/m 3
		r p	10 mm	H	1.5 GPa
		W d	62 mm
		r d	10 mm

Table 3

 3 Stochastic input normal random variables for the springback shape parameter study

	Parameters	Mean	Std dev
	Blank thickness	0.8 mm	3 %
	Young's modulus	70.5 GPa	3 %
	Yield Strength	0.180 GPa	3 %
	Poisson's ratio	0.342	3 %
	Friction coefficient	1.50e-1	3 %
	Radius of the punch	10 mm	3 %
	Radius of the matrix	10 mm	3 %
	Clamp force	600 kN	3 %

Table 4

 4 Empirical error and number of retained coefficients using LARS procedure

	Springback parameters	Err emp	Nb of terms
			LARS	Full PCE
	ρ	1,042e-4	112	1287
	β 1	1,009e-4	381	1287
	β 2	1,05e-4	383	1287

Table 5

 5 Proportion of terms sorted according to their degree in LARS based PCE for each response ρ, β 1 , β 2 Histograms of the springback shape responses obtained using the sparse PCE model

	Polynomial degree	Number of terms in PCE
		ρ	β 1	β 2
	0	1	1	1
	1	8	6	8
	2	3 3	28	31
	3	6 9	92	97
	4	0	254	245
	5	0	0	0

Table 1 :

 1 Some orthogonal polynomial types with respect to different continuous probability density function types

  (b) shows it for an 8-variate polynomial in log scale.

	Number of PCE terms	15 20 25 30 35 40 45	2 3 rd order PCE nd order PCE 4 th order PCE 5 th order PCE 6 th order PCE 7 th order PCE 8 th order PCE				Number of PCE terms	10 10 10 10	5 2 3 4	2 nd order PCE 3 rd order PCE 4 th order PCE 5 th order PCE 6 th order PCE 7 th order PCE 8 th order PCE
		10
		0 5	0.2	0.4 q norm	0.6 parameter	0.8	1	10	0 1	0.2	0.4 q norm	0.6 parameter	0.8	1
				(a) 2-variate					(b) 8-variate
	Figure

  • Compute the correlation vector c(k) 0 . Retain the predictor ψ α * where α * = argmax|c(k) 0 |. The model becomesy A M LARS (k) . • Update γ α to γ * α = γ α + ϵ *(1) , where ϵ * (1) is the LARS step where another predictor ψ β * has as much correlation with the current residual as does ψ α * (see Tibshirani and al. 2009 for numerical computations of ϵ * ). Add β * to the current index set of retained polynomial and update y A M LARS (k) . • Update jointly {γ α } α∈A M

Table 2 .

 2

	Geometry	Material	Loading Contact
	L s : 300 mm	E : 71 GPa	F b : 300 N µ : 0.15
	h 0 : 0.81 mm W s : 1 mm	ν : 0.342 ρ : 2700 kg/m 3	s : 60 mm
	r p : 10 mm K 0 : 576.79 MPa
	W d : 62 mm	ϵ 0 : 0.3593
	r d : 10 mm	n 0 : 0.01658

Table 2 :

 2 Geometrical, material, loading and contact parameter of the U-shaped B-U-T model .

  12: Definition of springback parameters, ρ, β 1 and β 2 .

	Variables	ρ	Responses β 1	β 2
	Thickness	1e-5 1e-5 1e-5
	Young's Modulus	1e6 1e6 1e6
	K 0	1e6 1e6 1e6
	ϵ 0	1e-2 1e-2 1e-2
	ν 0	1e-2 1e-2 1e-2
	Poisson's coefficient 1e-2 1e-2 1e-2
	Friction coefficient 1e-2 1e-2 1e-2
	Clamp force	1e1 1e1 1e1

Table 3 :

 3 Sensitivity threshold obtained for µ *

i = 2

Table 4 :

 4 Full stochastic model under study for the deep drawing process application4.4.1. Validation of Fat-LHSIn order to illustrate the efficiency of the Fat-LHS, we propose to compare the value of the two first statistical moments of each response and illustrates for two different sampling methodologies: 20 a classical LHS sampling and the proposed sensitivity constrained LHS sampling. We illustrate the results only on the 1/ρ response, the interpretation remaining unchanged for the other outputs.

	ξ	min	max	E[ξ]
	Thickness (h 0 )	0.805	0.815	0.8 mm
	Young's Modulus (E b )	70.5	71.5	71 GPa
	K 0	575.79 577.79 576.79 MPa
	ϵ 0	0.3493 0.3693	0.3593 -
	ν 0	0.015	0.017	0.01658
	Poisson's Coefficients (ν)	0.325	0.335	0.33
	Friction coefficients (µ)	0.155	0.170	0.162
	Clamp force (F)	34.5e3	35.5e3	35e3 kN

Verification aims at assessing the accuracy of the solution in the context of the study by comparing in some useful sense the consistency of the numerically obtained results with some trustworthy experimental results at hand (Do we solve right the equations?)

Validation consists in providing a quantitative description of the distance between the real experiment, the numerical model or the real-life system it is aimed to represent (Do we solve the right equation?).

The notions of random variables and probability density function are described in more details in Chapter 2

For a quantitative illustration see Fig.7in Chapter 3 or Fig.13in Chapter 4

The resolution of a model has been defined on p.8.

J. Lebon ( ) • G. Le Quilliec • P. Breitkopf • P. Villon Laboratoire Roberval, UMR 7337, Université de Technologie de Compiègne, BP 20529, 60205 Compiègne cedex, France e-mail: jeremlebon974@gmail.com J. Lebon • R. Filomeno Coelho BATir Department CP 194/2, Université Libre de Bruxelles (ULB), 50, avenue Franklin Roosevelt, 1050 Brussels Belgium

Acknowledgements

Acknowledgments This research was conducted as part of the OASIS project, supported by OSEO within the contract FUI no. F1012003Z. The first author is grateful to Innoviris for its support to BB2B project entitled multicriteria optimization with uncertainty quantification applied to the building industry. The authors also acknowledge the support of Labex MS2T.

Acknowledgements

This research was conducted as part of the OASIS project, supported by OSEO within the contract FUI no. F1012003Z. The first author is grateful to Innoviris for its support to BB2B project entitled multicriteria optimization with uncertainty quantification applied to the building industry. The authors also acknowledge the support of Labex MS2T.

, Frédéric, and my beloved partner Sabine, nothing would have been possible. Finally, I would like to gratefully acknowledge the support provided for this work under the BB2B project from the Innoviris institution (Bruxelles Capitale).

Chapter 2

Non-intrusive Polynomial Chaos Expansion for uncertainty quantification 2.2. Spectral representation of stochastic field through Polynomial Chaos Expansion

• The regression based methods based on a linear regression of the output random vector on the output vector. Among these methodologies, one may count the Pearson correlation coefficient and the Partial Correlation Coefficient. These are useful to quantify the effect of random input onto the random output if the stochastic model is nearly linear, but fails when in most of the case when the model is non linear.

• The variance based methodologies are more general. Broadly speaking they are based on a decomposition of the variance of the output of interest into summands First of all let us give an overview of the Sobol decomposition. We consider ξ = {ξ 1 , . . . , ξ M } ∈ R M , M independent identically distributed random variables and y(ξ) an output function. The Sobol decomposition of the y function is given by:

where u is a set of integers, ξ u = (ξ u1 , . . . , ξ us ) with s=card(u). This summand counts 2 M elements. Each of the terms of this summand (except the mean value y 0 = ȳ) verifies for any ξ ui the following property: Ω y u (ξ u )p(ξ i )dξ i = 0 ∀u ∋ i which implies the orthogonality of the functions y u :

and thus the uniqueness of the Sobol decomposition. Each Sobol function y u may be computed according to:

where ξ ∼u denotes the ξ vector without the u components.

We denote y 0 the mean of the output y, V the total variance of the output y,

and V u the conditional variances of the Sobol functions y u .

Thanks to the orthogonality of the decomposition, we have:

The Sobol sensitivity indices are then defined as: