Batoul

Donghai

Daniel

Montassar, Menouar, Pierre, Shiqi Jim Hannah

Yanfeng

Alaa Abbas

Bilal

Fatima

Keywords: Uncertainty quantification, Backward uncertainty propagation, sensitivity analysis, Sobol method

Nearly three years have passed since I came to Toulouse, three years with some dark days, but many light days.

Résumé

Contents

Acknowledgment iii

Résumé vii

Table of Contents x Figures xiii

List of . 85 3.15 The difference between the EFDCs obtained experimentally (blue) and the EFDCs obtained using the first order matrix approximation models of A, B, C, and D (red). Different values of (d, v, w) are considered: (a) (10 nm, 6 v, 6 µm) , (b) (100 nm, 8 v, 40 µm), (c) (50 nm, 6 v, 6 µm), and (d) (10 nm, 15 v, 40 µm). . . . . . . . . . . . . . . . . . . . . . . . . . . . 88 3.16 The difference between the EFDCs obtained experimentally (blue) and the EFDCs obtained using the second order matrix approximation models of A, B, C, and D (red). Different values of (d, v, w) are considered: (a) (10 nm, 6 v, 6 µm)

General Introduction

Nowadays, no one can doubt the basic role of modeling in any scientific process. In short, a model is a systematic description of the relationship between input and output. It aims to imitate, translate, or predict the behavior of real systems. However, undesirable disturbances may prevent a model from achieving its aim perfectly. Indeed, scientists can not be completely accurate during the modeling process, while constructing the model and collecting the input data. This leads to the presence of uncertainty, indicating a state of being unsure about the correctness of the performance of the model. Consequently, it is becoming no longer acceptable to submit any scientific project without providing a comprehensive assessment of the reliability and validity of the results under the effect of uncertainty. For that, studying uncertainty is becoming of a great interest in almost all disciplines, considering it an essential procedure for robust modeling. This manuscript presents our work in developing new methods to deal with uncertainty, and applying such methods in different domains.

The input of a model usually consists of variables and parameters. The variables are said to be uncertain if it is not sure that the values given to run the model are the actual true ones. Mainly, variable uncertainty appears due to imperfect measurements, inherent variability, and incomplete data collection. In addition, the model is said to have parameter uncertainty if its parameters are not surely characterizing the real system. Such uncertainty arises due to poor calibrations, imprecise estimations, or bad curve fittings.

Moreover, the model is said to have structural uncertainty if we cannot be confident that the form of the model is accurately imitating the real studied system. Usually, structural uncertainty appears due to ambiguity in the definitions of the given concepts, limitations of the acquired knowledge in the studied domain, or difficulties in some systems to be represented as equations or codes.

Whenever there is uncertainty in the input or the model structure, the output of course will be uncertain, hence the obtained results can not be trusted. Moreover, this output uncertainty may end up with severe consequences, especially in some sensitive domains like risk assessment and decision making. A simple example in this manner is the uncertainty in the construction of an airplane. This may happen due to uncertainty in the global process from the pre-design step up to the final manufacturing and assembling.

Such uncertainty, combined to unpredicted meteorological configurations, can be a reason for airplane crashes and crises. As a consequence, scientists insist that uncertainty cannot be tolerated or ignored, and it should be studied carefully.

The methods developed in this manner are of different concerns, depending on the source of the uncertainty (input, parameter, or model structure) and on the aim of the modeler. Some methods are concerned just with structural uncertainty. Other methods are dedicated for parameter and input variable uncertainty. Generally, these methods can be classified into three groups:

1. Structural Uncertainty Assessment: Methods in this group are applied only in case of structural uncertainty. They seek an optimal model representation of the real system with a reduced structural uncertainty as much as possible.

2. Uncertainty Quantification: Methods in this group are mainly applied in the case of input uncertainty. Their goal is to find a quantitative characterization of the uncertainty.

Two types of uncertainty quantification methods exist: forward uncertainty propagation and backward uncertainty propagation. In the forward uncertainty propagation the uncertainty of the output is to be quantified by propagating the uncertainty of the input. In the backward uncertainty propagation, the input uncertainty is to be determined starting from the given output uncertainty.

3. Sensitivity Analysis: Methods in this group are applied in the case of input uncertainty. Their aim is to find which input elements have the genuine output impact. This helps in indicating whether an uncertain input element will cause a significant uncertainty in the output or not.

Note that the methods of sensitivity analysis are even applied in the case of absence of uncertainty. Indeed, their main goal is to detect the effect of each input on the output by studying the output variation with respect to the variation of the inputs. However, in an uncertainty study a sensitivity analysis is first done to detect the influence of each uncertain input according to its variation on its range of uncertainty. Then, for a model having several inputs, only the most influencing inputs are taken into account in the uncertainty study while the others are fixed at some specific values.

Relative to the methods of the groups considered above, the backward propagation of uncertainty is the one with least consideration in literature. In addition, rare studies have conducted input variable uncertainty, knowing that this uncertainty is frequently encountered especially in problems where the input variable is an experimental data. For that, in this thesis we focus on deriving a new backward uncertainty propagation method applied mainly for input variable uncertainty. In parallel to this, we followed a new manner while applying Sobol method which is an already existing sensitivity analysis method.

The applications are done on two real models in order to keep the context of our work and the obtained conclusions realistic. This work is presented in this dissertation which is organized as follows:

In Chapter 1 we give a general review of the main methods in structure uncertainty assessment, uncertainty quantification, and sensitivity analysis. This will provide the necessary background concepts needed to study and understand uncertainty. There will be a clear focus on the uncertainty quantification and sensitivity analysis methods since these two topics are the main concern of our work. The methods are thoroughly described with their applicability and limitations.

In Chapter 2 we present a new derived backward uncertainty propagation method. The aim of this method is to determine the input uncertainty starting from the given data of the uncertain output. The main idea is to partition the output uncertainty between the inputs using the probabilistic representation of uncertainty. This partition helps in generating a nonlinear system of equations whose unknowns are the uncertainties of the inputs. The system can be solved solved numerically as a non linear least square problem, and the input uncertainty is obtained. The method is mainly applied in case of having input variable uncertainty, especially for problems with inputs coming from experimental data. Different examples are also presented in this chapter in order to see how the method is applied in reality. In general, the method is simple, however the partition of the output uncertainty becomes more complicated in the case of having a big number of inputs with a complex form of the model. In this case, sensitivity analysis can be used to detect the most impacting inputs, so that the backward propagation can be restricted to these important inputs. In this work we consider one of the sensitivity analysis methods, called Sobol method, and develop the way of applying it and analyzing its results. These ideas are presented in detail in the next chapters.

In Chapter 3 we present our first application of the Sobol sensitivity method. The aim of this application is to examine the performance of the Sobol method in case of having a model whose explicit form is unknown, plus having a limited number of data points to apply the sensitivity method numerically. The model considered in this application is from the domain of force spectroscopy, in which we study the sensitivity of an experimental curve called Electrostatic Force Distance Curve (EFDC). This curve is obtained by a microscopic scanning technique, which uses a very thin tip to scan surfaces. The EFDC plots the electrostatic force between this tip and a scanned dielectric. The EFDCs for different experimental settings are given as experimental data, then using this we study the sensitivity of this curve with respect to the variation of the settings which are considered as inputs. To derive the conclusions concerning the performance of the Sobol method we use Design Of Experiment (DOE). DOE is a methodology for designing experiments that allows, by some special plots, the analysis of the effect of the experimental factors on the response. For that we validate the sensitivity results of the Sobol method by the plots of DOE. All these notions are presented in detail in this chapter as well as the obtained conclusions and consequences.

In Chapter 4 we present a different framework for the application of the Sobol sensitivity method from the already existing ones. Usually with the Sobol method the sensitivity is studied by computing for each input a sensitivity index that reflects the effect of this input on the output. Our idea in this chapter is to extend this by studying the evolution of these indexes with respect to an outside factor or an experimental condition like time, distance, and temperature. The aim of this extension is to detect the most convenient conditions at which conclusions about the impact of each input can be derived. In addition, studying the variation of Sobol indexes helps in giving more information about the inputs, which helps also in the backward propagation of uncertainty. In this manner we consider two different models.

The first model is from the domain of computer vision, which is a programming representation of a 3D reconstruction method called Shape-From-Template (SFT). This method uses a single 2D image and a 3D template to recover a deformed 3D surface. We study the sensitivity of this model with respect to the depth of the surface in front of the camera, its orientation, and the focal length of the camera by which the image is taken. The sensitivity indexes are then computed and analyzed as a variation of the depth. This helps in revealing how, at each depth, the position of the surface affects the quality of the reconstruction. All these specific points, the description of the SFT method, its sensitivity study and the conclusions derived are presented in this chapter.

The second considered model is a model for charge transport in dielectrics. This model is usually modeled using a set of partial differential equations however in this work we consider it as a black box model. We study its sensitivity with respect to four of its main inputs. The sensitivity indexes are studied under the variation of three experimental conditions: the temperature, the time and the intensity of the applied electric field.

The aim of this study is to find the experimental conditions at which each input has the significant impact on the output. This highlights the experimental conditions that should be followed in order to acquire a data suitable for estimating each input. The details of these ideas and the results obtained are all presented in this chapter.

Finally we close up with the Conclusion chapter that gives a full summary of the work with the conclusions drawn and the future perspectives.

Introduction

Methods for studying uncertainty are of great necessity in all disciplines. These methods are categorized into three main groups: structural uncertainty assessment, uncertainty quantification, and sensitivity analysis. In this chapter, we provide a literature review of these three groups. First, we start by giving the general notation of a model that will be used throughout this manuscript. Next, we briefly present the concept of the structural uncertainty assessment. Then we focus on the main methods of uncertainty quantification and sensitivity analysis. To ease the explanation of these methods we introduce the notion of the input uncertainty representation. Then we start by reviewing some forward uncertainty propagation methods, which are divided into two groups, probabilistic methods and non probabilistic methods. After that, we continue by reviewing some sensitivity analysis methods which are also divided into two main groups: local and global. For each presented method, we give its basic idea and then discuss its applicability and limitations.

General notation of a model

Lastly, we present the state of the art of the backward uncertainty propagation methods and then we finish off with a conclusion.

General notation of a model

Any model can be expressed formally as

F (x, α) = y (1.1)
The symbols x, y and α refer to the model input variable, output, and parameter. The model input variable x is a vector of n components in D, where D is the domain of F and It is important to note that, throughout this manuscript the model structure F is assumed to be deterministic and not stochastic i.e. it produces the same output when it is run exactly with the same input. In addition, the notation of the model that will be used is F (x) = y, the parameter symbol α is removed for simplicity as x plays the same role.

D ⊆ R n (n ≥ 1). x is
In the next section we give a brief overview of the concepts of structural uncertainty and the methods used in this manner.

Structural uncertainty assessment

Structural uncertainty usually represents the lack of confidence that the structure of the constructed model reflects adequately the studied real system. In such cases, a modeler could not be sure that the output would be correct, even if the true values of all the inputs and parameters are known. To clarify this type of uncertainty, consider for example the case study of having the corresponding data points of a studied system and the aim is to find its structure. As shown in Fig. 1.1, the data points can be fitted by a logarithmic function and by a 4th root function. In this case it is not sure which formula is the true representation of the model, and hence there is a structural uncertainty.

Figure 1.1: The data fitted to a logarithmic function and to a 4th root function.

During modeling, several factors may lead to structural uncertainty. One of these factors is the simplifications and scientific judgments made when constructing and interpreting the model. Other factors are the incomplete understanding of the system under study and the inappropriate equations used to express this system. Even though these causes were taken in consideration, the definite elimination of structural uncertainty was impossible, however the emphasis was to reduce this uncertainty as much as possible. Methods of structural uncertainty assessment try to find the optimal model structure that best represents the real studied system.

Multi-model analysis [Lu et al., 2013] is mainly applied in this domain, in which plausible models are built representing all the possible emulations of the true structure. Then a strategy is used to derive the aspect of the final model from the set of possible models.

Model Selection [Leeb et Pötscher, 2005] and Model Averaging [START_REF] Strong | Structural uncertainty in health economic decision models[END_REF] are the two most popular and broad approaches used in this context. In Model Selection, an optimal model from the set of models is selected according to some criteria provided by experts. Some criteria that are proposed in [START_REF] Bojke | Characterizing structural uncertainty in decision analytic models: A review and application of methods[END_REF] include Residual Mean Squared Error, Finite-Prediction-Error, Minimum Variance Criteria, and subjective probabilities. On the other hand, in Model Averaging, instead of choosing one single model, the weighted average of the proposed models is taken. In this case, a suitable weight is assigned to each plausible model according to how much it matches reality.

Another method to cope structural uncertainty more flexibly is known as "Parameterization of structural uncertainty" [START_REF] Strong | Managing structural uncertainty in health economic decision models: a discrepancy approach[END_REF]. The idea is to describe the uncertainty of the model by introducing new uncertain parameters, for instance correction factors, boolean elements, or exogenous variables. Thus a single general model is constructed in such a way that every other plausible model is considered as a particular case of this general model. In such a situation, the general model is considered to be studied under the effect of only input uncertainty with no structural uncertainty.

Although these strategies seem straightforward and really helpful in reducing and studying structural uncertainty, they still have some limitations. Indeed, using Model Selection may be disadvantageous in several cases. This is because selecting one model will probably discard important eventualities from other alternative models. On the other hand, the Model Averaging strategy allows the collection of all possible models, however, when using large models with with highly computational cost, it becomes difficult to find the average. Parameterization of structural uncertainty is practical when dealing with structural uncertainty, however not all structural uncertainties are that easy to be represented by a parameter. As a conclusion, a modeler should be cautious while choosing the appropriate approach to reduce the structural uncertainty.

In this section, a general idea about methods that deal with structural uncertainty was presented. Although different methods exist in this manner, but none of them can guarantee that the true model can be attained. Thus the goal was to achieve the best level of confidence while choosing the structure of the model.

In the following sections, we continue with the methods of uncertainty quantification and sensitivity analysis, which are the main focus of our work. However in the sequel we will assume that we are dealing with an exact true model that has no structural uncertainty.

Input uncertainty representation

For the sake of clarity, it is important to note that methods that study input uncertainty rely first on finding a representation of the uncertainty. Several ways have been proposed.

However, the most practical and used one is the probabilistic way. In the following paragraph, in each method of uncertainty quantification and sensitivity analysis the associated representation way of uncertainty will be introduced. In exception, the probabilistic way is introduced here since its notion is used in most methods and in Chapter 2.

The uncertainty at an input point a = (a 1 , • • • , a n ) is represented probabilistically by a random vector (vector of random variables) which will be denoted by

X = (X 1 , • • • , X n ).
Each X i represents the uncertainty at a i . The uncertainties at a i 's are assumed to be independent, and hence the random variables X 1 , • • • , X n are mutually independent. The realizations of each random variable X i are the uncertain values x i supplied to F that are supposed to be equal to a i . These realizations are collected in a set denoted by Ω i . Then the collection of the realizations of the random vector x is the set Ω

= Ω 1 × • • • × Ω n .
To clarify this notation of uncertainty consider the following example. Let

F (x 1 , x 2 ) = πx 1 x 2 2 (1.2)
be the function that gives the volume of some liquid in a cylinder. The first input x 1 represents the height attained by the liquid in the cylinder in mm, and the second input

x 2 represents the radius of the cylinder also in mm. Suppose that for a certain liquid the values of x 1 and x 2 are 9.762 and 3.91 respectively, but these true exact values are not known. Thus to get the volume of the liquid, measurements for x 1 and x 2 should be done.

The measurement process is pictured in Fig. 1.2.

Note that with the measurements of Fig. 1.2 we can not be sure about the exact values of x 1 and x 2 for the given liquid. This implies that there is input uncertainty for both Figure 1.2: The measurements of the height x 1 and the radius of the cylinder x 2 .

variables x 1 and x 2 . Let X 1 be the random variable representing the uncertainty at x 1 , and let X 2 be the random variable representing the uncertainty at x 2 . Then the input uncertainty of F at the point (9.762, 3.91), which we suppose is not known, is represented by the random vector (X 1 , X 2 ). Note that the realizations of the random vector X 1 are all the possible true values for x 1 and they are collected in a set Ω 1 . Similarly for X 2 , its realizations are all the possible true values of x 2 and they are collected in a set Ω 2 .

From Fig. 1.2, one can define the sets Ω 1 and Ω 2 by Ω 1 = [9. 5,10] and Ω 2 = [3.5, 4]. We associate here intervals to the sets Ω 1 and Ω 2 since x 1 and x 2 take real values. So any real value between 9.5 and 10 is a possible true value for x 1 , and this gives an interval, and similarly for x 2 its possible true values are collected in an interval.

Actually what we present here is a very simple example of uncertainty that one may face.

However in big problems, uncertainty could be much more complicated and it cannot be eliminated even with highly accurate measurements. Now we continue with the probabilistic representation of the input uncertainty. Each random variable X i , as a representation of the input uncertainty, is associated with a probability distribution which specifies the probability of each possible value x i to be the true value a i . Here we distinguish two cases:

1. If Ω i is not finite, X i is a continuous random variable and its probability distribution is characterized by a probability density function. A simple example of this continuous case, is the uncertainty presented in the inputs of the example of Fig. 1.2. Realizations in this example are real numbers taken between two limits and Ω i is an interval, thus it is not finite. In this case each X i is characterized by a probability density function. . 1.3(a) shows the probability density function of the height presented in the example above.

Fig

2.

If Ω i is finite, X i is a discrete random variable and its probability distribution is characterized by a probability mass function. An example of this discrete case, is an uncertainty in an input which represents the number of items sold by a store per year. Realizations in this example are natural numbers and could not be real numbers, and so the Ω i is for sure finite. In this case the random variable is characterized by a mass function (see Fig. In this work, we assume that the random variables X 1 , • • • , X n are continuous, however the techniques presented are equally applicable for discrete distributions. The probability density functions of X 1 , • • • , X n are denoted by p 1 (x 1 ), • • • , p n (x n ). Accordingly, the probability distribution of the random vector X is p X (x) = n i=1 p i (x i ). Usually, the probability distributions of the random variables are selected based on either prior data or subjective judgments of experts. In [START_REF] Hammonds | An Introductory Guide to Uncertainty Analysis in Environmental and Health Risk Assessment[END_REF], the authors proposed some simple guidelines to derive the appropriate continuous distribution. These guidelines state that when the data are limited and the uncertainty range is relatively small, a uniform distribution can be used and the associated support interval is characterized by the uncertainty range. If there is more knowledge about a most likely value or midpoint, in addition to the range of the uncertainty, a triangular distribution may be assigned. When the range of the uncertainty is very large, a log-uniform or log-triangular distribution may be more appropriate than the uniform or the triangular distribution.

The assumption of normal, log-normal, or empirical distributions usually depends on the availability of the relevant data, where a fitting process is usually used to guess such probability density functions. In addition to this, the authors in [START_REF] Hammonds | An Introductory Guide to Uncertainty Analysis in Environmental and Health Risk Assessment[END_REF] indicated that other continuous distributions can be also used such as Gamma, Beta, and Poisson. Note that, analogous guidelines can be derived for the discrete case.

With this probabilistic representation of the input uncertainty at a, the corresponding output uncertainty at F (a) is represented by a random variable denoted by Y . This random variable is defined as Y = F (X 1 , • • • , X n ), and its probability density function is denoted by p Y . Note that the random variable Y and the random variables X 1 , • • • , X n are dependent, as y is a function of the other random variables.

The best estimate of the true value of a is given by the mean of X which will be denoted by µ

X = (µ 1 , • • • , µ n ).
Similarly, the best estimate of the true value of F (a) is given by the mean of Y which will be denoted by µ Y .

Concerning the quantity of uncertainty at a, it is usually represented by the variance of

X, denoted by V = (V 1 , • • • , V n ).
However, for some specific distributions, other statistical parameters can be used to represent the quantity of uncertainty. For instance, for a uniformly distributed random variable, the radius of the support vector could be used to represent the quantity of uncertainty since it gives the dispersion of the values around the expected value. Similarly, the quantity of the output uncertainty is usually represented by the variance, which will be denoted by V ar(Y ). However, the quantity of the output uncertainty can also be represented by other statistical parameters depending on the associated distribution. In the following , whenever the probabilistic point of view is used to represent uncertainty, the representative of the quantity of the uncertainty will be stated explicitly.

Forward uncertainty propagation

Forward Uncertainty Propagation is performed to investigate the uncertainty in the model's output that is generated from the uncertainty in the model's input [Marino et al., 2008]. The idea is to associate a representation for the input uncertainty, then accordingly try to find the output uncertainty in the same representation type. Methods in this group are classified as either probabilistic and non probabilistic, according to the way the input uncertainty is represented. Probabilistic methods use the probabilistic point of view to represent the input uncertainty (see section 1.4). Then through propagation, the probabilistic presentation of the model output is to be determined. Mainly, we seek the mean µ y and the variance V ar(Y ) of the the random variable Y which represents the output uncertainty. The non probabilistic methods, however, use a non probabilistic forms to represent the input uncertainty. Then, according to the used form, the associated output representation is to be determined. In the following two subsections, classical methods from both categories are presented.

Probabilistic methods

In the probabilistic uncertainty propagation methods, each uncertain input is represented by a random variable and it is specified by a probability density function. Then, model output, as a function of the model input, is also a random variable whose statistical moments are to be determined.

values {x k = (x (k) 1 , • • • , x (k) n )} k=1•••M
are drawn randomly from the distributions of X i , according to their probability density functions. These sample points are then run by the model F to obtain their corresponding output values. The obtained values of the output are then collected to find its statistical characteristics. For instance, the computation of the expectation and variance of the output Y is done using the approximation formulas:

µ Y = 1 M M k=1 F (X (k) ) (1.3) V ar(Y ) = 1 M -1 M k=1 (F (X (k) ) -µ Y ) 2 (1.4)
On the other hand, the distribution of Y can also be derived simply by using the obtained output data in a fitting process.

As it can be seen, the method is simple from a theoretical point of view, universally applicable and does not require any assumptions on the model form as linearity or continuity. Moreover, the number of sample points M needed in the simulation is generally independent of n, the size of the input vector x [START_REF] Helton | [END_REF]. However, it is important to note that the method converges to the exact stochastic solution as the number of samples goes to infinity, so thousands or millions of samples may be required to obtain accurate estimations [Iaccarino, 2009]. This could be problematic in the case of computationally expensive models and/or in the case of important size of input vector.

Several In Fig. 1.5, the Latin Hyper cube sampling shows some gaps and clusters like the random sampling because the sample size is too small. Note that as the number of samples in-creases, the number of rows and columns increases, and hence there will be more coverage of the range of inputs.

To conclude, the Monte Carlo method is a practical method for propagating uncertainty, however its performance depends on the number of samples available. With the development of high performance software, computation of big samples becomes more easy.

However it remains a problem for the models with high computation cost.

Spectral Methods

Spectral methods represent an alternative strategy for uncertainty propagation. Their basic idea is to write the model as an infinite sum of some basis functions. Indeed, this expansion eases the derivation of the moments (mean and variance) of the output random variable Y . Different basis functions have been used for such expansion, depending on the distribution of the input random variables X 1 , • • • , X n . However the most used basis are polynomials. Several approaches with different polynomial basis have been used [Gilli, 2013]. The most known one in this manner is the Polynomial Chaos Expansion (PCE), defined using multidimensional orthogonal polynomials as representative basis. Here, we give a general overview about the Polynomial Chaos Expansion, considering it as a typical illustration of the expansions of the spectral methods.

The first proposed Polynomial Chaos expansion employed the Hermite polynomials in terms of Gaussian random variables to generate the expanded series. According to [Lee et Chen, 2009], its expression is:

u = a 0 H 0 + ∞ i 1 =1 a i 1 H 1 (φ i 1 ) + ∞ i 1 =1 i 1 i 2 =1 a i 1 i 2 H 2 (φ i 1 , φ i 2 ) + ∞ i 1 =1 i 1 i 2 =1 i 2 i 3 =1 a i 1 i 2 i 3 H 3 (φ i 1 , φ i 2 , φ i 3 ) + ... (1.5)
for an arbitrary random variable u, where {φ i 1 } ∞ i=1 is a set of standard normal variables, H i is a generic element in the set of multidimensional Hermite polynomials of order i, and a i are coefficients to be determined. For convenience, the expression was rewritten in a more compact way:

u = ∞ i=0 b i Ψ i (φ) (1.6)
where b i and Ψ i (φ) correspond to a i 1 i 2 ...ip and H p (φ i 1 , φ i 2 , ..., φ ip ) respectively. Note that the orthogonality property with respect to the standard normal probability density function of the Hermite polynomials implies that

E[Ψ i Ψ j ] = E[Ψ 2 i ]δ ij and E[Ψ i ] = 0 f or i = 0 (1.7)
where E represents the expectation. Hence, the set {Ψ i } forms an orthogonal basis of the space of functions having normally distributed variables.

In practice, when there are n uncertain inputs standard normally distributed, the output response can be approximated by n-dimensional PCE, truncated at some order p [Lee et

Chen , 2009]. In this case, the number of terms in PCE becomes P + 1 where P is given as

P = p s=1 (n + s -1)! s!(n -1)! (1.8)
Thus, the model output approximation is given by

Ỹ = P i=0 b i Ψ i (X) (1.9)
with X = (X 1 , X 2 , ..., X n ). The derivation of the coefficients can be carried out analytically [Schick, 2011] [START_REF] Sepahvand | Uncertainty quantification in stochastic systems using polynomial chaos expansion[END_REF].

Unfortunately, this generalization maintained the property that only identical random variables can be involved in the same expansion. Thus, later, the approach was extended by assuming that the input vector components are independent distinct random variables, where a multi-dimensional basis is constructed as a simple product of the corresponding constructed one dimensional orthogonal polynomials.

Once the expansion of the model function is obtained, the moments (mean and variance)

of the output can be derived. While the accuracy of the generalized PC approach can be improved by increasing the polynomial order of truncation, it should also be noted that as the number of inputs and the expansion order increase, the number of unknown coefficients to be determined increases exponentially, thereby increasing the computational costs [Kewlani et al., 2012]. Thus this method is suitable for models with a small number of uncertain inputs.

Perturbation Method

The Perturbation method is an alternative way for uncertainty propagation, based on local expansion of the model function [Sudret, 2007]. The idea is to consider the truncated Taylor expansion of the model F in the neighborhood of µ X the mean of the input random vector X:

F (X) = F (µ X ) + n i=1 ∂F ∂X i X=µ (X i -µ i ) + 1 2 n i=1 n j=1 ∂ 2 F ∂X i ∂X j X=µ (X i -µ i )(X j -µ j ) + o( X -µ X 2 ) (1.10)
Then the expectation of Y is:

µ Y = E[Y ] = E[F (X)] ≈ F (µ) + n i=1 ∂F ∂X i X=µ E[(X i -µ i )]+ n i=1 N j=1 ∂ 2 F ∂X i ∂X j X=µ E[(X i -µ i )(X j -µ j )] (1.11) Note that E[(X i -µ i )] = 0 for every i, and E[(X i -µ i )(X j -µ j )] = Cov(X i , X j ).
However, the assumption of independence between the inputs uncertainties implies that Cov(X i , X j ) = 0 for all i = j. Hence, the expectation of Y is simplified into:

µ Y = E[F (X)] ≈ F (µ) + n i=1 ∂ 2 F ∂ 2 X i X=µ V i (1.12)
As it can be seen, the expectation of the uncertain output is approximated by two terms.

The first term is F (µ X ), the value of the model F at the mean of the uncertain input, which is called the first order approximation of µ Y . The second term is a second order correcting term which depends on the variances of the inputs and the partial derivatives of the model form F . In a similar manner, the variance of Y can be approximated. So starting from the formula:

V ar(Y ) = E[(Y -E(Y )) 2 ] ≈ E[(Y -F (µ)) 2 ] (1.13)
Using the Taylor expansion, then V ar(Y ) is approximated by:

V ar(Y ) ≈ E     n i=1 ∂F ∂X i X=µ (X i -µ i )   2   ≈ n i=1 n j=1 ∂F ∂X i X=µ ∂F ∂X j X=µ E[(X i -µ i )(X j -µ j )]
(1.14)

As E[(X i -µ i )(X j -µ j )] = 0 for all i = j, V ar(Y ) ends up with:

V ar(Y ) ≈ n i=1   ∂F ∂X i X=µ   2 V i (1.15)
An interpretation of this approximation implies that the variance of the response is the sum of contributions of each input, where each contribution is a mix of the variance of this input and the gradient of the response with respect to this input.

This method appears quite general and it is applied at a low computational cost, especially if the gradient of the model response is available. However, it can be applied only for models with small uncertainties, due to the local nature of the Taylor expansion approximation.

Note that the choosing between spectral method and perturbation method for uncer-tainty propagation, is the same as choosing between high accuracy and less cost. Spectral methods rely on complete expansion of the model while the perturbation method relies on the local expansion of the model. The first costs more, but give more accurate results than the second. Thus in applications one should compromise between cost and accuracy.

In this subsection, three different uncertainty propagation methods were revised. The common point between these methods is that they rely on the probabilistic representation of uncertainty. Their main goal is to find the probabilistic moments of the output uncertainty. This is done either numerically by simulation (Monte Carlo method) or by using an expansion of the output formula (spectral and perturbation methods). In general, these methods are considered simple from a theoretical point of view. In the next section, non probabilistic methods for uncertainty propagation are presented. The concept of some of these methods may be considered as generalization of the probabilistic approach, however some other methods have completely different notions.

Non-Probabilistic methods

As discussed in the previous section, the probabilistic uncertainty propagation methods are usually applied when the information about the uncertain inputs are sufficient to construct a probabilistic distribution. However, if the information of the uncertain element is insufficient, the non-probabilistic approaches can be used [START_REF] Gao | Hybrid probabilistic interval analysis of bar structures with uncertainty using a mixed perturbation monte-carlo method[END_REF]. Below, alternative approaches to the probabilistic approach of uncertainty presentation are discussed, including Interval theory, Fuzzy theory, Possibility theory, and Evidence theory.

The way the uncertainty is propagated using these uncertainty presentation approaches is also discussed.

Interval Analysis

Interval analysis [Moore et al., 2009] is one of the simplest ways to propagate uncertainty in data-poor situations. In interval analysis, it is assumed that nothing is known about the uncertain input except that it lies within certain bounds. Each uncertain input is represented by an interval. An interval here refers to a close set in R, which includes the possible values of a number. An interval is usually expressed as

[a, b] = {c ∈ R, a ≤ c ≤ b},
where the real numbers a and b are the lower and upper limits of the interval. According to the author of [P Swiler et al., 2009], the problem of uncertainty propagation can be turned into an interval analysis problem: given the inputs' uncertainties represented by intervals, what is the corresponding interval of the output uncertainty?

The interval of the output uncertainty is determined by finding the infimum and supremum attained by the model when the uncertain inputs vary entirely over their associated intervals. When the model is a simple expression with simple arithmetic operations (+, -, ×, ÷), then the output interval is determined by extending these elementary arithmetic operations to intervals. Such extension is as follows

I op J = {c op d such that c ∈ I, d ∈ J} (1.16)
where I and J are two intervals, and op refers to one of the arithmetic operations. However, the determination of the output interval becomes much more complicated if the model has a complex form. Non linear optimization problems might be used to determine the upper and the lower limits of the output interval. This probably requires a large number of model evaluations. Furthermore, most optimization solvers are local, and thus the global optima is not guaranteed. Which means that it is not easy to find the infimum and supremum. Hence, to solve interval analysis problems properly, global methods must be used, and usually these approaches can be very expensive.

In brief, interval analysis is suitable for propagating uncertainty in problems where the model has an elementary form, with a small number of inputs.

Fuzzy Set Theory

Another approach to represent the uncertainty of the input is the Fuzzy Set Theory. In the Fuzzy Set Theory, an uncertain input is treated as fuzzy number, and its corresponding uncertainty is characterized by membership functions. Such membership functions associate a weight between 0 and 1 to every possible value of the uncertain input. Then, to propagate the input uncertainty and find the output uncertainty, the output membership function is to be determined. To this end, some procedures are carried. Here, we give a quick review of the notions of the Fuzzy Set theory as it is used in uncertainty propagation.

The starting point will be the definition of the widely used term "Fuzzy Set". Let U be a universe set of τ values (elements), set A is called a fuzzy set if it is composed of ordered pairs in the following form:

A = {(τ, Φ A (τ )) | τ ∈ U, Φ A (τ ) ∈ [0, 1] } (1.17)
where Φ A (τ ) is the degree of membership of τ in A. The function Φ A is called the membership function of A.

A special case of fuzzy sets are the so-called fuzzy numbers. According to [START_REF] Schulz | [END_REF], a Fuzzy number A is a fuzzy set satisfying the following :

1. The membership function Φ A is piecewise continuous, and the elements τ are real

numbers, i.e. A = {(τ, Φ A (τ )) | τ ∈ U ⊆ R, Φ A (τ ) ∈ [0, 1] }.
2. A is normal, i.e. there exists at least one (τ, Φ A (τ )) ∈ A such that Φ A (τ ) = 1.

3. A is a convex set, i.e. for any (τ 1 , Φ A (τ 1 )), (τ 2 , Φ A (τ 2 )), and (τ 3 , Φ A (τ 3 )) ∈ A, the following implication holds:

τ 1 < τ 3 < τ 2 ⇒ Φ A (τ 3 ) ≥ min{Φ A (τ 1 ), Φ A (τ 2 )} (1.18)
This simple concept of fuzzy numbers permits its application in the domain of uncertainty propagation. Indeed, an uncertainty at one real valued input point a i of point a = (a 1 , • • • , a n ) is represented by a fuzzy number A i with a specific membership function Φ i . The elements τ of A i are the uncertain values supplied to F that are supposed to be equal to a i . The membership function Φ i is usually constructed based on the available information about the uncertain input. In this manner, a general followed guidance is:

the closer Φ i (τ ) is to 1, the more the element τ is accepted as a true value for a i , and the closer it is to 0, the less it is accepted. Several forms of the membership function of a fuzzy numbers has been developed and used in the domain of uncertainty [Wierman, 2010], including triangular and trapezoidal shaped membership functions. Illustrations After associating to each uncertain input a fuzzy number, the propagation of uncertainty is performed using the extended principle of fuzzy set theory [Maskey et al., 2004]. According to this principle, the membership function of the output y is given by:

Φ y (y) =        sup{min(Φ 1 (x 1 ), • • • , Φ n (x n )) y = f (x 1 , ..., x n ) 0 if no (x 1 , • • • , x n ) exist such that f (x 1 , • • • , x N ) = y (1.19)
Thus, the uncertain output is represented by a fuzzy number {(y, Φ y (y))}. Its membership function is determined at each possible value y of F (a) by applying the formula (1.19). This indicates the spread of the output uncertainty as well as the most probable true value F (a).

This method of uncertainty propagation has been applied in various domains, examples are in [START_REF] Schulz | [END_REF]Maskey et al., 2004], where usually the number of uncertain parameters is small, or the model form is monotone with respect to the inputs. In cases where more complex models have to operate on fuzzy numbers, the above procedure results in nonlinear numerical optimization problems at each possible value of the output, and hence this would be computationally expensive.

The Dempster-Shafer Theory

The Dempster-Shafer Theory(DST), also known as evidence theory, is an uncertainty propagation method used when the available information is mostly provided by experts.

Unlike probability theory, in DST there are two measures of likelihood and not a single probability distribution function. These two measures are called belief and plausibility.

The following paragraph is a description of the derivation of such measures as well as their role in uncertainty propagation. Readers interested in more details about DST are referred to [Baraldi et Zio, 2010].

As a first step for quantifying the input uncertainty in DST, a finite set Ω i is assigned to each uncertain input x i at the uncertain point a i consisting of all the uncertain values supplied to F that are supposed to be equal to a i . Then a mass function is associated with each set Ω i , called the Basic Belief Assignment (BBA) and denoted by B i . The BBA is a mapping B i : P(Ω i ) -→ [0, 1] satisfying:

B i (φ) = 0 and A⊂Ω i B i (A) = 1 (1.20)
where P(Ω i ) is the power-set of Ω i i.e. the set of all subsets of Ω i . For every A ⊂ Ω i , the value B i (A) indicates how likely the true value of the uncertain input falls within the subset A. Each A ⊂ Ω i with B i (A) > 0 is called a focal element of B i . Note that, as

A⊂Ω i B i (A) = 1
, the BBA function has a finite number of focal elements. Moreover, a BBA function is completely defined by these focal elements and their associated masses according to [Limbourg, 2008].

In DST, the function B i is not the fundamental measure of likelihood. Rather, this BBA is used to derive the two measures of likelihood: the belief and the plausibility. According to [START_REF] Helton | An exploration of alternative approaches to the representation of uncertainty in model predictions[END_REF], the belief, Bel(A), and the plausibility, P l(A), for a subset A ⊂ Ω i are defined by:

Bel i (A) = C⊆A B i (C) (1.21)
and

P l i (A) = C∩A =∅ B i (C) (1.22)
Observe that these two summations can be simply computed since B i has a finite number of focal elements. Further more, the belief Bel i (A) is commonly considered as a lower bound of the probability that the true value of x i is within A. On the other hand, the plausibility P l i (A) is considered as an upper bound for this probability. Thus, together, belief and plausibility define an interval-valued probability distribution [Bel i (A), P l i (A)],

and not a single probability distribution [START_REF] Giunta | [END_REF]. Now, as the uncertainty of each input x i is characterized by Bel i and P l i , it is now the question about the uncertainty of the output. In fact, to propagate the uncertainty using DST theory, it is required to find the belief and plausibility functions of y, denoted Bel y and P l y respectively. In [START_REF] Helton | An exploration of alternative approaches to the representation of uncertainty in model predictions[END_REF], the author indicated a practical way to find both Bel y and P l y . First a generalization of the BBA denoted by B x is defined on

Ω = Ω 1 × • • • × Ω n by: B x (A) = n i=1 B i (A i ) (1.23)
where

A = A 1 × • • • × A n ⊆ Ω and each A i ⊆ Ω i .
Then, the image of Ω by the model function F is by definition

F (Ω) = {F (x), for every x = (x 1 , • • • , x N ) ∈ Ω 1 × Ω 2 ... × Ω N } (1.24)
This F (Ω) represents the set of all possible values of the uncertain output y. Accordingly, for every A ⊆ F (Ω), the two measure functions of its belief and plausibility are defined as follows:

Bel y (A) = Bel y (F -1 (A)) = C⊆F -1 (A) B x (C) (1.25)
and

P l y (A) = P l y (F -1 (A)) = C∩F -1 (A) =∅ B x (C) (1.26)
where F -1 (A) = {x ∈ Ω; F (x) ∈ A}. These functions Bel y and P l y characterize the uncertainty in the output y.

Despite the simplicity of the concept of this method, its application in practice depends on the number of the uncertain inputs involved in the analysis. In addition, it depends on the number of possible values of each uncertain element. Moreover, the evaluation of Ω, its image by F , and then the two measures become more complex as the number of uncertain inputs increase.

Possibility Theory

Possibility theory provides another alternative to probability theory for uncertainty rep-resentation. As Dempster-Shafer Theory, possibility theory involves two measures for likelihood, a necessity and a possibility. The following is a summary of how necessity and possibility are derived, and their contribution in uncertainty propagation.

In any modeling process, where the available data is incomplete to assign a random variable for an uncertain input x i at an input point a i in a probabilistic mode, possibility theory can be used instead of probability theory. A possibility function P os i is defined on Ω i , the set of all possible values x i . This function assigns to each element ξ in Ω i a degree in the interval [0, 1]. Such a function specifies which values are possible to be the true ones and which are less possible. It also differentiates between the values that are surprising to be true and the values that are expected to be true. More precisely, if for some ξ ∈ Ω i , P os i (ξ) = 0, this means that ξ is impossible and totally surprising to be true. On the other hand, if P os i (ξ) = 1, then ξ is considered as totally possible and unsurprising to be true, however it could not be true at all. Indeed, P os i (ξ) = 1 is a much weaker property than saying probability is 1 [Ripamonti et al., 2013]. Analogously to the Dempster-Shafer Theory, possibility theory introduces two likelihood measures, possibility and necessity.

For any A ⊆ Ω i , the possibility and necessity for A are defined as follows [Ripamonti et al., 2013]:

P i (A) = sup ξ∈A P os i (ξ) (1.27)
and

N i (A) = 1 -P i (A c ) = 1 -sup ξ / ∈A P os i (ξ) (1.28)
where A c = Ω i \ A. After finding the P i and N i of each uncertain input, then, in a similar way to Dempster-Shafer method, the possibility and necessity of y are derived. These two measures P y and N y characterize the possibility of each of the possible values of the outcome. Moreover, as the number of the inputs increase, it become more expensive to compute such measures.

In this section several non probabilistic approaches for forward uncertainty propagation were presented. The idea is first to find a representation of the input uncertainty and then find output uncertainty according to this representation. In this manner, some methods adopt simple representations like interval theory, and other methods use more compli-cated concepts like fuzzy set theory. In addition, some non probabilistic methods use a generalization of the probabilistic representation to propagate uncertainty. Nevertheless the probabilistic approaches are still considered to be more practical than the non probabilistic methods, for that, probabilistic representation is more preferred in applications.

In the next section we give the state of art of the backward propagation of uncertainty and see how it differs from the forward propagation.

Backward uncertainty propagation

The basic idea of the Backward uncertainty propagation is to start from some output data in order to quantify the input uncertainty [Chantrasmi et Iaccarino, 2012] Even more, there are no clear guidelines in literature for the applications done in this manner. Maybe this is because the backward propagation problem can be an ill posed problem in some cases, since different partitions of the input uncertainties may give the same output uncertainty. This makes the problem more challenging.

The first attempt for a backward propagation method was to find an inverse of the model and then apply a forward propagation procedure. In this way the input uncertainty can be derived by propagating the output uncertainty as in [Chen et al., 2015]. One of the big obstacles for this approach is that most models cannot be easily inverted, especially complicated black box models. In some applications, modelers try to find a meta model of the initial model, in which this meta model can be inverted. For instance, in [START_REF] Baumgärtel | Inverse uncertainty propagation for demand driven data acquisition[END_REF] the model is approximated by a Gaussian process. Then using this Gaussian process the mean and the variance of the input uncertainty are computed. Another example in [Chantrasmi et Iaccarino, 2012], where a non continuous model is approximated by a fraction of polynomials. Then using this approximation, the backward propagation is done. Unfortunately, not all models can be easily approximated into simpler ones and

Backward uncertainty propagation

get their inverse and then apply a forward propagation of uncertainty. Thus a backward propagation method adapted to deal with the model itself, not its inverse, is really needed.

In the few other applications done in backward propagation domain, the main concept that is used is the Bayesian Inference. In statistics, inference refers to the process of deducing the properties of a probability distribution by the analysis of data [Keeping, 1962]. Bayesian inference, which is one of the approaches of statistical inference, is used to quantify the input uncertainties in a probabilistic framework based on the available output data. However, it is also required to have at least a basic idea about the input uncertainty. In the following paragraphs, we give the general formalism of the Bayesian inference used in the backward uncertainty propagation.

Bayesian Inference

The basic idea of the Bayesian inference is to use the output data and prior information about the input uncertainty in order to give a probabilistic characterization of the input uncertainty. In this method, quantification of the input uncertainty is done as an update of the already acquired information.

First, the prior information of the input uncertainty is used to drive a probabilistic distribution of the random variable X called prior distribution. This prior distribution represents the probability distribution of the input uncertainty before any output data is collected. Usually, this distribution is extracted from physical constraints, expert knowledge, and previous experimental data [Nagel, 2017]. For instance, uniform distributions are often chosen for inputs that can be bounded from above and below, as having some physical constraints. Gaussian or lognormal distributions are often used for parameters that are unbounded or strictly positive. Alternatively, the Principle of Maximum Entropy [Jaynes, 1957] provides an objective method to determine suitable prior distributions that yield optimal representation of the uncertainty given the available information.

The backward propagation of uncertainty is done by updating the prior distribution using the output data. The new updated distribution is called the Posterior distribution. It is defined as the probability distribution of X conditioned to the given data represented by Y . Its formula is given by the well known Bay's Law:

p X (x|y) = P XY (x, y) p Y (y) = p Y (y|x)p X (x) p Y (y) (1.29)
In the above equation, p Y (y|x)is known as likelihood function. This function is considered as the connection between the output and the input uncertainty. p Y (y) is the distribution of Y , and it can be computed using the following formula: [START_REF] Jaynes | Information theory and statistical mechanics[END_REF] where dp refers to integration over the probability space of X. An illustration of this backward propagation method is given in Fig. 1.7 : In Fig. 1.7, both the prior and the posterior distributions are normal distributions.

p Y (y) = p Y (y|x)p X (x)dx (1.
However in reality the posterior is a complex probability distribution even if the prior distribution is a simple one [Nagel, 2017]. Moreover, the posterior and the prior distributions need not to be of the same family. If the prior and the posterior distributions are in the same family, then they are are called conjugate distributions.

Note that the posterior probability density function holds all the information of the input uncertainty. So for instance, the mean value which represents the best estimate of the uncertain input is computed by :

E[X|y] = xp X (x|y)dx (1.31)
Even the covariance matrix of X = (X 1 , • • • , X n ) can be computed using the formula:

Cov(X|y) = (x -E[X|y])(x -E[X|y]) T p X (x|y)dx (1.32)
where (x-E[X|y]) T refers to the transpose of the vector (x-E[X|y]). By the assumption that the X i 's are mutually independent, the covariance matrix is a diagonal matrix, where the diagonal elements are the associated variances representing the quantity of the input uncertainty. Even though formulae (1.31) and (1.32) are used in the derivation of input uncertainty, the computation is usually done numerically without finding the explicit formula of the posterior distribution.

As it can be seen, Bayesian inference for backward propagation depends mainly on the prior distribution. This initial guess of the input uncertainty cannot be easily established.

In chapter 2, we introduce a new backward propagation method that does not rely on the prior distribution of the input uncertainty. Instead the input uncertainty is computed by decomposing the variance of Y in terms of the statistical moments of the X i 's. Note that the problem of backward propagation becomes more difficult as the number of inputs increases. In such a case, sensitivity analysis is used in order to detect the most important inputs and hence involve them in the uncertainty study. In the next section we give the general idea of sensitivity analysis and then we detail the most common sensitivity methods. The focus will be on the Sobol method, which constitutes the cornerstone of the work of chapter 3 and 4.

Sensitivity analysis

The modeling of complex systems usually requires a large number of inputs. Carrying an uncertainty analysis involving all these inputs would be a real burden. However, in most real world problems, only a limited number of inputs happens to influence the response significantly [Sudret, 2007]. Thus, it is essential to determine which inputs contribute most to the output variability and which of them are insignificant so that they can be ignored during investigation. In this manner, sensitivity analysis has gained considerable attention, as it assesses how variations in the model output can be apportioned to different input sources [Marino et al., 2008]. So, sensitivity analysis indicates how much each uncertain input contributes to the output uncertainty. This helps in detecting the essential inputs.

Sensitivity analysis was firstly performed qualitatively [Gan et al., 2015], in which a statement of confidence as "low", "moderate", "high" was given to describe the impact of an uncertain input on the output. Unfortunately, such statements are difficult to be interpreted practically. Later, quantitative methods appeared substantially, giving a new form for the sensitivity analysis. In quantitative sensitivity analysis methods, an index is assigned to each input, reporting how much the output responds to the changes in the values of this input. In fact, these methods enable analysts to rank the inputs according to their impact on the output, and thus improving the state of knowledge in order to reduce output uncertainty more effectively.

Quantitative sensitivity analysis methods are divided into two main categories: local and global [START_REF] Saltelli | Sensitivity analysis wiley series in probability and statistics[END_REF]. Local sensitivity analysis methods derive the sensitivity index of each input by computing or approximating the partial derivative of the model function in a specific neighborhood of the input. On the other hand, global sensitivity analysis methods derive the sensitivity indexes by allowing inputs to vary over the whole range of their possible values [Tong, 2007]. In the next two subsections, methods from both categories are represented. Recall that the model is assumed to be deterministic, i.e. it gives identical results when it is run with exactly the same set of input values.

Local sensitivity analysis

Local sensitivity analysis concentrates on the local impact of the inputs on the model's output [Cacuci, 2003]. It is based on the computation of the partial derivatives of the model with respect to each input at some specific value of the input. Here we will assign the value at which the partial derivative is computed by a, so this is to say the sensitivity of the model F is studied locally at a. The idea behind using the partial derivative to study sensitivity comes from the fact that: varying a single input x i in the vicinity of a i , while keeping other inputs fixed, will provide a comprehensive assessment of the effect of

x i on the model's output. The partial derivative then is called the sensitivity index of F at a i . Numerous techniques have been developed to compute the gradient of a model efficiently, here the main ideas of the direct method, indirect method, and the automatic differentiation methods will be represented consecutively.

When the considered model is expressed as an explicit algebraic equation, a symbolic computation of the partial derivatives of this equation can produce the sensitivity index.

Such a method is usually referred to as the direct method. Moreover, to ensure a more comparable concept, normalized sensitivity coefficients can be computed instead. The idea is just to multiply the partial derivative by the ratio of the value of an input over the value of the output. So the normalized sensitivity index of F at a i is:

S i = ∂F ∂x i (a) × a i F (a) (1.33)
The advantage of the normalized sensitivity coefficients is that they are dimensionless, this facilitates identifying the most sensitive inputs. Note that these methods cannot be applied to models with non specified mathematical equations, like computer codes. In addition, it may become troublesome when having models with complicated expressions.

In such situations, symbolic computation of the partial derivatives becomes inefficient even with the use of the modern computer algebra systems such as Mathematica and Maple. In these cases, numerical methods are much preferred.

One of the most applicable numerical methods for deriving partial derivatives is the finite difference method, also known as the indirect method or the Brute Force method. The model is solved first at the specified value of the inputs, then a perturbation is added just to one input element and the model is resolved. This generates sensitivity indexes in the following form

S i = F (a 1 , • • • , a i + ∆x i , • • • , a n ) -F (a 1 , • • • , a i , • • • , a n ) ∆x i (1.34)
Thus for computing sensitivity indexes for n inputs, this method requires at least n + 1 model evaluations. Hence, one should be cautious in the case when n is large.

Despite the large applicability of the indirect method, it is usually associated with the challenge of selecting the input step size ∆x i . As the indirect finite difference method intended a local investigation of sensitivity, then by choosing a large step size this aim will not be satisfied. On the other hand, if the step size is too small, the difference between the original and perturbed solutions can be so small, leading to serious numerical errors in the solution.

To overcome such limitations of the numerical methods, more sophisticated techniques have been instituted, named under Automatic Differentiation (AD) [Sandu, 1997]. In general, AD is a way to give numerical values of the derivatives while using algebraic manipulation steps for the computation, thus giving more accurate values than the numerical methods. Moreover, unlike the direct method which uses symbolic computation just for models formulated under explicit mathematical expressions, AD techniques differentiate computer codes. Indeed, in AD a compiler analyses the code of the model, then adds some instructions to the code, which are needed to compute derivatives, in a manner that reduces complexity and saves computational time. Then, the new expanded code can automatically evaluate the partial derivatives of the output with respect to the inputs with minimum human effort. The basic idea of such a process lies within the fact that any computer program, no matter how complicated, performs a sequence of binary ( +, -, ×, ÷) or/and unitary (sin, power, square root, log) operations. Thus, by applying the chain rule successively to these operations, derivatives of arbitrary order can be computed automatically, exactly up to machine error, and using more arithmetic operations than the original program.

Two main modes are usually used in the AD techniques, the forward mode and the backward mode. In the forward mode, chain rule is applied from inside to outside while in the backward mode, chain rule is applied from outside to inside. So for example, if some function Υ(t) = g(h(t)), then by chain rule ∂Υ ∂t = ∂g ∂h ∂h ∂t , thus in the forward mode ∂h ∂t is computed first, while ∂g ∂h is firstly computed in the backward mode. Although, this method is much more effective than other differentiation methods, it is given less attention. This is mostly because it is poorly understood, plus it is frequently confused with the better known symbolic and numerical differentiation methods.

In addition to the above three differentiation techniques, several alternative methods exist to find local sensitivity indexes, such as the Green function methods and polynomial approximation. However all such methods only investigate the behavior of the model in a small neighborhood of a i , which may not take into account all the possible values of the uncertain input a i . In addition, these methods do not taking into consideration the effect of inputs' interactions on the output. In order to generalize the concept of sensitivity analysis and overcome the limitation of local sensitivity analysis, the notion of global sensitivity analysis methods were developed. The next section details thid notion for better investigation of sensitivity analysis.

Global sensitivity analysis

In global sensitivity analysis methods, inputs are varied simultaneously over their entire range of possible values [Bokov, 2012]. The effects on the output of both individual inputs and interactions between them are assessed in terms of sensitivity indexes. Several approaches have been dedicated to derive such indexes, each based on different concepts.

In general, all these approaches are based on the probabilistic framework in representing uncertainty. So the uncertain output is represented by Y random variable, and each uncertain input is represented by random variable X i .

The elementary techniques for global sensitivity analysis are specified for the linear case.

For instance, if the studied model depends linearly on its input, then Correlation coefficients (CC), Partial correlation coefficients (PCC), or Standardized Regression coefficients (SRC) are used as sensitivity indexes. Note that such linearity can be easily detected by examining scatter plots, in which two dimensional graphics are plotted representing the variation of Y with respect to X i . Actually, these plots give a full understanding of the relationship between y and each x i including monotonicity and dependency.

The Correlation Coefficient (CC), also called Pearson Correlation Coefficient, provides a measure of the strength of the linear relationship between X i and Y [START_REF] Helton | Survey of sampling-based methods for uncertainty and sensitivity analysis[END_REF].

According to [Marino et al., 2008], it is defined as follows:

C(X i , Y ) = Cov(X i , Y ) V ar(X i ) V ar(Y ) (1.35)
The CC has a value between -1 and +1. A positive value of CC indicates that X i and Y either increase together or decrease together. However, a negative value of CC indicates that X i and Y tend to move in opposite directions [START_REF] Helton | Survey of sampling-based methods for uncertainty and sensitivity analysis[END_REF]. On the other hand, values of CC close to 0 indicate a weak (linear) relationship between X i and Y whereas values close to -1/1 indicate the relationship is strong. Thus, taking the absolute value of the CC would be an appropriate indication of the effect of X i on Y i , and hence the sensitivity of Y with respect to X i .

Although CC considers the effect of X i on Y , it does not take into account the possible effects on Y due to other uncertain inputs. A more thorough representation is the Partial correlation coefficients (PCC) [Marino et al., 2008]. The Partial correlation characterizes the linear relationship between the input X i and the output Y after the linear effects of the remaining inputs on Y are discounted. This is done using a sequence of regression models. First, the following two regression models are introduced:

Xi = a 0 + n p=1 p =i a p X p and Ŷ = c 0 + n p=1 p =i c p X p (1.36)
Then the CC between the two residuals X i -Xi and Y -Ŷ is the PCC of X i . Indeed these two residuals remove the linear effect of the other inputs on Y . So PCC represents the sensitivity of Y due to only X i .

Another convenient indication for sensitivity assuming linearity is the Standardized Regression Coefficient (SRC) [START_REF] Saltelli | Sensitivity analysis of model output: an investigation of new techniques[END_REF]. Practically, a least square procedure is usually used to construct a regression model describing the relation between the output and the inputs:

Y = b 0 + n i=1 b i X i (1.37)
Then the regression coefficients b 0 , • • • , b n can characterize the influence of each input on Y . However, a drawback of these coefficients is that they depend on the units in which the X i 's and Y are expressed. To cope with this problem, a normalized form of the regression is adopted, having the following expression :

Y -µ Y V ar(Y ) = n i=1 (b i V i V ar(Y ) ) X i -µ i V i (1.38)
where µ i and µ Y are the means (averages) of the X i and Y values respectively, V ar(Y ) and V i are the variances of the Y and X i respectively. The coefficient b i

V i V ar(Y ) is called the Standardized Regression Coefficient (SRC) of X i .
According to the author of [START_REF] Helton | [END_REF], the absolute values of the SRCs can be used to provide a measure of the inputs' importance.

While the above three coefficients CC, PCC, and SRC are well behaved sensitivity indexes in the linearity case, they perform poorly in the nonlinear case [START_REF] Helton | [END_REF].

However, if the relationship between the input and the output is still monotone, a rank transformation can be used to reduce the effect of non linearity [START_REF] Helton | [END_REF].

In rank transformation, the sampled data are replaced by their corresponding ranks, and then the usual regression and correlation procedures are performed on these ranks. Specifically, the smallest value of each variable is ranked by 1, the next smallest value is ranked by 2, and so on up to the largest value, which is ranked by the value of the sample size. To clarify this, consider the following example. A model with two inputs X 1 and

X 2 , having samples {(α 1 , β 1 ), (α 2 , β 2 ), (α 3 , β 3 ), (α 4 , β 4 )}, such that α 4 < α 1 < α 3 < α 2 and β 3 < β 2 < β 4 < β 1 . Then, their corresponding ranking is {(2, 4), (4, 2), (3, 1), (1, 3)}.
Considering the corresponding outcomes :

X =           α 1 β 1 α 2 β 2 α 3 β 3 α 4 β 4           -→           F (α 1 , β 1 ) = Y 1 F (α 2 , β 2 ) = Y 2 F (α 3 , β 3 ) = Y 3 F (α 4 , β 4 ) = Y 4           (1.39) such that Y 3 < Y 1 < Y 4 < Y 2 , the new ranked data is           2 4 4 2 3 1 1 3           -→           2 4 1 3           (1.40)
The analysis is then performed with these ranks used as the values for the input and output variables. The sensitivity indexes of such analysis are the Ranked Correlation Coefficients (RCCs), the Partial Ranked Correlation Coefficients (PRCCs), and the Standardized Ranked Regression Coefficients (SRRCs) instead of CCs, PCCs, and SRCs, respectively [Pereira et Broed, 2006]. As for the linear case, the absolute values of such coefficients may give a convenient degree of the importance of each input. Note that, the use of rank transformed data results in an analysis based on the strength of monotonic relationships rather than on the strength of linear relationships [START_REF] Helton | [END_REF].

Since linearity and monotonicity are specific for a limited number of models, more sophisticated methods have been investigated to derive sensitivity indexes for general cases.

The most popular approaches in this manner are the ANOVA methods, an abbreviation of Analysis of Variance methods. The basic idea of such methods is to decompose the total variance of the output into a sum of partial variances. Each partial variance corresponds to an input or group of inputs. Such decomposition indicates the contributions of each input variable or group of inputs in the output variance. The sensitivity indexes are derived as the ratios of the partial variances over the total variance of the output.

In the following subsection, one of the main ANOVA methods, called Sobol method, is presented in detail, as it is a cornerstone in the applications presented in chapters 3 and 4.

Sobol Method

A very powerful sensitivity analysis technique that is gaining popularity in many fields is the Sobol method [Sobol, 1993]. It is based on the decomposition of the variance of the output into a sum of partial variances, each contributing either to an individual input or a group of inputs. Then sensitivity indexes, called Sobol indexes, are introduced as the fraction of each of these partial variances over the total variance. The following paragraphs form a detailed description of the concept behind the formulation of Sobol indexes, as well as the strategies applied to compute them.

Consider the random variable

Y = F (X 1 , • • • , X n ), where X = (X 1 , X 2 , ..., X n ) is a
random vector of independent inputs. Each random variable X i is characterized by probability density function p i (x i ). In the sequel, and for simplicity, the notation dp i will be adopted instead of p i (x i )dx i when integrating over the probability space associated to X i . Moreover, dp will represent the product of all the density measures n i=1 dp i , and the notation dp ∼i will refer to the product of all density measures accept dp i .

A starting point for the Sobol method was to consider the ANOVA decomposition of F [Sobol, 1993]. Indeed, whenever F is integrable, it can be written as summands of increasing dimensions as:

F (X 1 , X 2 , ..., X n ) = F 0 + n i=1 F i (X i ) + n 1≤i<j≤n F ij (X i , X j ) + • • • + F 1•••n (X) (1.41)
where F 0 is a constant, and

F i 1 •••is (X i 1 , X i 2 , ..., X is )dp i k = 0 for 1 ≤ k ≤ s (1.42)
The classical properties for this decomposition, as stated in [Sudret, 2007], are:

1. The constant F 0 is equal to the mean of Y :

F 0 = F (X)dp (1.43)
2. The summands are orthogonal to each other in the following sense:

for {i 1 , i 2 , ..., i s } = {j 1 , j 2 , ..., j t },

F i 1 i 2 ...is (X i 1 , X i 2 , ..., X is )F j 1 j 2 ...jt (X j 1 , X j 2 , ..., X jt )dp = 0 (1.44)
Accordingly, the terms of the above decomposition can be determined as follows:

F 0 = F (X)dp = E[Y ] (1.45) F i (X i ) = F (X)dp ∼i -F 0 = E[Y /X i ] -E[Y ] (1.46) F ij (X i , X j ) = F (X)dp ∼i,j -F i (X i ) -F j (X j ) -F 0 = E[Y /X i , X j ] -E[Y /X i ] -E[Y /X j ] + E[Y ]
(1.47)

Analogously one can proceed for the higher order terms. Now, squaring both sides of the decomposition formula (1.41), and then integrating over all the input spaces and by the orthogonality property, one ends up by:

F 2 (X)dp = F 2 0 + n i=1 F 2 i (X i )dp i + n 1≤i<j≤n F 2 ij (X i , X j )dp i dp j + • • • + F 2 12•••n (X)dp 1 dp 2 • • • dp n (1.48) Note that F 2 (X)dp -F 2 0 (1.49) is the variance Y , denoted V ar(Y ). Moreover F 2 i (X i )dp i is the variance of F i (X i ) and it is denoted by V ar(F i ). Also, F 2 ij (X i , X j )dp i dp j is the variance of V ar(F ij (X i , X j )
) and it is denoted by V ar(F ij ), and so follows for the the higher order terms.

Thus the variance of Y , can be written as a sum of partial variances:

V ar(Y ) = n i=1 V ar(F i ) + n 1≤i<j≤n V ar(F ij ) + • • • + V ar(F 1•••n ) (1.50)
According to this variance decomposition, the author in [Sobol, 1993] defined Sobol indexes as:

S i 1 •••is = V ar(F i 1 •••is ) V ar(Y ) (1.51)
The First-order Sobol index, S i , estimates the main effect of X i on Y . On the other hand, the Higher-order Sobol indexes, S i 1 •••is , estimate the corresponding effect of interaction between the inputs X i 1 , • • • , X is on Y . Note that, by dividing both sides in (1.50) by

V ar(Y ), we find that the sum of the all first and higher order Sobol indexes is equal to 1. Fig. 1.8 is an illustration of the Sobol indexes for model with two variables.

To illustrate Sobol's point of view of using variance decomposition to get the impact of each input on the output, consider the following clarification. If one predicted the influence of a certain input X i on Y , he could measure the variation of Y while fixing X i at a specific possible value a and keeping the other inputs to vary randomly. If the variation of Y stays the same as when all inputs vary, this means that X i has no influence on the output. However if the variation is smaller than that when all inputs vary, this means that X i is influencing the output. Note that, the scientific description of the variation of Y while fixing X i at a specific possible value a is nothing but V ar(Y /X i = a). Due to uncertainty, no special value can be given for X i . Thus it is preferable to take the mean of this variance over the whole space of possible values of X i to get better understanding of the influence of X i . In other words,

E[V ar(Y /X i )] represent a reference of the sensitivity of Y to X i . Indeed, comparing with the values of V ar(Y ), small values of E[V ar(Y /X i )]
indicates that X i highly affects Y , while large values indicates that X i has no strong effect on Y .

Knowing that:

V ar(Y ) = E[V ar(Y /X i )] + V ar(E[Y /X i ]) (1.52) then V ar(E[Y /X i ]) can
F i ) is exactly V ar(E[Y /X i ]), as F i (X i ) is E[Y /X i ] -E[Y ]
. Similar deduction can be done for the higher order terms.

In addition to the above effective indication for sensitivity, other proposed interesting sensitivity indexes are the Total Sobol indexes S T [Homma et Saltelli, 1996]. Such indexes describe the total contribution of an input X i , including all its interactions, on the output Y . So for the input X i , the total Sobol index S T i is defined as the sum of all Sobol sensitivity indexes involving X i , and it can be written as:

S T i = V ar(F i ) + j =i V ar(F ij ) + • • • + V ar(F 1•••n ) V = {i}⊆J S J (1.53)
where this last sum is overall J ⊆ {1, • • • , n} such that i ∈ J. The advantage of this Total

Sobol index is that it can be numerically computed with no need to compute all second and higher order Sobol indexes.

In general, in applications where Sobol indexes are applied to study sensitivity, the understudy models may have complex forms. Thus, computing Sobol indexes, which requires, as we see, integrating the model function, may be a hard mission. In fact, practically, only numerical methods are used. More precisely, according to the distribution each input, a sampling is done to end up with a sample space

{X (k) = (X (k) 1 , • • • , X (k) n )} k=1,••• ,M .
Then using this sample, different formulae can be used to approximate the the Sobol indexes [START_REF] Saltelli | Variance based sensitivity analysis of model output. design and estimator for the total sensitivity index[END_REF]. For example, if the mean of Y is estimated by:

F0 ≈ E[Y ] = 1 M M k=1 F (X (k) ) (1.54)
then, the total variance V can be estimated by:

V ar(Y ) ≈ 1 M M k=1 F 2 (X (k) ) -F0 2 (1.55)
On the other hand, the partial variances are usually approximated by another set of formulas. Noting that, in most applications the analyst often computes the first-order Sobol indexes, and sometimes the second order ones corresponding to the interaction of every two inputs. The following is one of the approximating formulas mentioned in [START_REF] Saltelli | Variance based sensitivity analysis of model output. design and estimator for the total sensitivity index[END_REF], it depends on approximating the partial variances using two different samples:

{X (k) = (X (k) 1 , X (k) 2 , • • • , X (k) n )} k=1,••• ,M , { X(k) = ( X(k) 1 , X(k) 2 , • • • , X(k) n )} k=1,••• ,M (1.56) Considering V ar(F i ) = V ar(E[Y /X i ]), one may write V ar(F i ) = E 2 [Y /X i ]dp i - E[Y /X i ]dp i 2 (1.57) But E[Y /X i ]dp i = F (X)dp ∼i dp i = F (X)dp = F 0 (1.58)
and so this part can be approximated by F0 in (1.54). On the other hand, E[Y /X i ] can also be written as

F (X)dp ∼i (1.59)
where X i is considered in this integral as a constant and integrating is over the dummy variable X ∼i . Thus, it is possible to write:

E 2 [Y /X i ]dp i =   F (X i , X ∼i )dp ∼i F (X i , X∼i )dp ∼i   dp i (1.60)
Rearranging this formula gives:

E 2 [Y /X i ]dp i = F (X i , X ∼i )F (X i , X∼i )dpdp ∼i (1.61)
Hence, the proposed estimator for V ar(F i )is written as:

V ar(F i ) ≈ 1 M M k=1 F (X (k) i , X (k) ∼i )F (X (k) i , X(k) ∼i ) -F0 2 (1.62)
where X i , X ∼i , X∼i are taken from the two different samples defined above. Similar approximation formulas are also defined to any higher order Sobol indexes as well as to the total order Sobol indexes. Thus, with these numerical approximations, Sobol indexes can be efficiently computed.

However, one should pay attention to the computational cost of such approximations, especially for models with a large number of inputs.

In this section the concepts of sensitivity analysis and some common sensitivity methods were revised. Two main groups of sensitivity methods were considered: the local methods and the global methods. The local sensitivity methods use mainly the partial derivative as an indication of the sensitivity of the output with respect to each input. Thus this kind of methods provides information only at the base point where the partial derivative is computed and do not take into account the rest of the variation ranges of the model inputs. On the other hand, global sensitivity methods give a more thorough comprehension of the sensitivity of the output with respect to the inputs, as they take into account the interaction and the global ranges of variation of the inputs. These global methods rely on the probabilistic point of view to describe the variation of each input. In this manner, the ANOVA method, Sobol method is considered as one of the strongest methods, and we rely mainly on it in our applications where we develop its usage. Note that all sensitivity analysis methods aim to detect the influence of the inputs of a model on its output. This in fact plays an important role in an uncertainty analysis, since it restricts the study to the significant inputs and this simplifies the uncertainty study. In particular this helps a lot in the backward propagation of uncertainty, since as the number of the uncertain inputs increases the complexity of the problem increases.

Conclusion

In this chapter, several methods that deal with uncertainty during modeling are revised.

First the concepts of the structural uncertainty assessment methods are briefly reviewed.

Then different methods for forward uncertainty propagation are described in detail. Then, the state of art of the backward propagation methods is presented. A little consideration of these methods has been taken in account in literature. Actually, this was our main motivation for deriving a new backward propagation method which is presented in chapter 2. In addition, various sensitivity analysis techniques have been reported and discussed in this chapter. The main focus was on the Sobol method, which is the main method used in our work presented in chapters 3 and 4.

CHAPTER 2

Introduction

One of the remaining challenges in modeling under uncertainty is the inability to quantify the input uncertainty given the output uncertainty. In most problems which rely on simulations and experiments, different outcomes are obtained even when the experiment is carried out at a supposed same value of the input. This is usually due to a hidden disturbance in the values of the input, which is indeed an unknown input uncertainty. So given the data of the obtained uncertain output, how the input uncertainty, that causes this output uncertainty, can be quantified and apportioned between different elements of the input. This type of problems has been rarely tackled in the domain of studying uncertainty, one cause for that would be the difficulty of such problems. Another reason could be the inability to guarantee solutions for such problems that can be ill posed in some cases since several solutions for the input uncertainty may correspond to the same output uncertainty. In this chapter we present a new backward uncertainty propagation method, concerned with finding the input uncertainty starting from the data of the uncertain output.

The method uses the probabilistic point of view to represent uncertainty. In addition, it consists of two main steps. Its basic idea is to solve a nonlinear least square problem whose residues are defined using the formula that expresses the variance of the output in terms of the statistical moments of the inputs. The solution of this problem gives directly the quantification of the input uncertainty. Unlike most previously proposed backward propagation methods, this method does not rely on any prior data or information about the input uncertainty. In addition, it is applied for both cases, whether the input uncertainty is represented by a normal distribution or a uniform distribution. Even more it is applied in the case when the uncertain inputs are not represented by the same distribution.

In following paragraphs, we first introduce again the problem of the backward propagation using the probabilistic notation. Then the two main steps of the method are presented. After that some numerical examples are described to prove the validity of the method. The chapter ends with the derived conclusions as well as the obtained notes. Concerning the quantity of the input uncertainty, the representation is more specific.

Problem definition and notation

Indeed, in uncertainty analysis, the uncertain inputs are mostly considered as either normally distributed or uniformly distributed, depending on the way measurements and experiments are done. For that we restrict our method here to these two types of distributions. Accordingly, if X i has a normal distribution, the quantity of the uncertainty is best represented by the variance of X i denoted by V i . However, if X i has a uniform distribution, then X i is associated by a support interval centered at its mean µ i . So let [µ i -r i , µ i + r i ] be the support interval of X i if it is uniformly distributed. In this case, the quantity of uncertainty is best represented by the radius of the support interval r i , which indicates the dispersion of the values around µ i .

For that we consider the two sets I 1 and I 2 defined by:

I 1 = {i ∈ N | 1 ≤ i ≤ n and X i has normal distribution} (2.1) I 2 = {i ∈ N | 1 ≤ i ≤ n and X i has uniform distribution} (2.2)
Let V = {V i } i∈I 1 be the set of the variances of all normally distributed X i 's and let R = {r i } i∈I 2 be the set of the radii of the support intervals of all uniformly distributed X i 's. Thus the elements of V and R represent the quantity of the input uncertainty. In our study case, the input uncertainty is unknown, and so the elements of the two sets V and R are unknown.

So to do a backward propagation of uncertainty, we should start with the given data of the uncertain output in order to quantify the input uncertainty. Using the notation defined in this section, the backward propagation problem can be defined as follow: starting from the values of µ X , µ Y and V ar(Y ) the elements of V and R are to be determined.

Backward propagation: Variance Partitioning

The idea of the presented backward propagation method is to consider the elements of V and R as unknowns. Then, we solve a nonlinear least square problem to find these unknowns. The residues of the least square problem are generated using the formula of V ar(Y ). The execution of this method is done in two main steps. The first step is to write the output variance V ar(Y ) in terms of the elements of V and R. The second step is to generate the least square problem using the derived expression of V ar(Y ).

Step 1: Output variance in terms of V and R

To derive an explicit expression of the output variance V ar(Y ) in terms of the elements of V and R the formula of the model F is used. In this manner, two different cases are distinguished: F is a multivariate polynomial in x 1 , • • • , x n and F is a smooth non polynomial function.

CASE 1:

If F is a multivariate polynomial in x 1 , • • • , x n , then it can be written in the following form:

F (x) = K j=1 q j Q j (x 1 , • • • , x n ) (2.3) where each of Q 1 , • • • , Q K is a monomial in x 1 , • • • ,
x n and K is the number of terms of F . In addition, the q j 's are the coefficients associated with the monomials Q j 's. To clarify this notation, consider an example of F a polynomial of three inputs defined by

F (x 1 , x 2 , x 3 ) = x 2 1 +4x 1 x 2 +x 2 x 2
3 . Then K = 3 and q 1 = 1, q 2 = 4 and q 3 = 1. Furthermore, the monomials are

Q 1 = x 2 1 , Q 2 = x 1 x 2 , Q 3 = x 2 x 2 3 .
According to the formulation of F in (2.3), the expression of the random variable Y in terms of the random vector X is given by:

Y = F (X) = K j=1 q j Q j (X 1 , • • • , X n ) (2.4) So Y is the sum of K different random variables each of the form Q j (X 1 , • • • , X n ),
where

Q j is a monomial of X 1 , • • • , X n .
Note that these K random variables are obtained from various multiplications of the random variables

X 1 , • • • , X n . For simplicity, each random variable Q j (X 1 , • • • , X n ) is used in the sequel as Q j .
Recall that, the variance of the sum of different random variables is the sum of the covariances of each couple of these random variables. Applying this to the Y in (2.4), then V ar(Y ) is written as:

V ar(Y ) = V ar K j=1 q j Q j = K j,k=1 q j q k Cov Q j , Q k (2.5)
However, according to the simplified formula, each covariance is decomposed into:

Cov Q j , Q k = E[Q j * Q k ] -E[Q j ]E[Q k ] (2.6)
where E [.] stands for the expected value of a random variable. Note that, since X 1 , • • • , X n are mutually independent, then Q j and Q k are independent if they have no X i in common with any degree. But if Q j and Q k are independent, the product of their expectation is equal to the expectation of their product, and so their covariance is equal to zero.

However, if Q j and Q k have at least one of the variables X 1 , • • • , X n in common with any degree, they are dependent. In addition, their covariance is expressed in terms of the non central moments of the variables X 1 , • • • , X n . To explain how this is done, consider a generic case of two dependent random variables

Q k and Q j such that Q k = X u l X v i and Q j = X s h X t i , where u, v, s, t ∈ N * . Q k and Q j are dependent since X i is a common random variable in their expressions. Using the fact that X 1 , • • • , X n are mutually independent, the expectations of Q k * Q j , Q k , and Q j are decomposed into: E[Q j * Q k ] = E[X u l X v i * X s h X t i ] = E[X v+t i ] * E[X u l ] * E[X s h ] E[Q j ] = E[X u l X v i ] = E[X v i ] * E[X u l ] E[Q k ] = E[X s h X t i ] = E[X t i ] * E[X s h ]
(2.7)

Recall that any expectation of the form E[Z z ], where Z is a random variable, is called the non central moment of Z of order z. We say non central because the random variable is not centered by its mean, so the central moment of Z of order z has the form

E[(Z -E[Z]) z ].
So the expectations which appeared on the right hand side in (2.7) are the non central moments of the random variables X l , X h and X i with different orders.

Substituting the expectation decomposition of (2.7) in the formula (2.6) implies that the covariance of the two dependent random variables Q k and Q j is written in terms of the non central moments of X 1 , • • • , X n . Hence, any non zero covariance in (2.5), which is in fact the covariance of two dependent random variables, is written in terms of the non

central moments of X 1 , • • • , X n . Consequently, V ar(Y ) is expressed in terms of the non central moments of X 1 , • • • , X n .
Note that, for any X i , the first order non central moment is the expected value µ i .

Furthermore, the higher order moments are defined as follows:

• If X i has a normal distribution: the non central moments are expressed as polynomials in terms of the first and second moments, the mean µ i and the variance V i , using the Moment Generating Function. See Appendix A for details.

• If X i has a uniform distribution: the t-th non central moment of X i is defined by the following formula:

E[X t i ] = 1 2r i µ i +r i µ i -r i x t i dp i = 1 t + 1 (µ i + r i ) t+1 -(µ i -r i ) t+1 2r i (2.8)
This expression of E[X t i ] can be simplified into:

E[X t i ] = 1 t + 1 t k=0 (µ i + r i ) k (µ i -r i ) t-k
(2.9)

So the moments of a uniformly distributed X i are expressed in terms of µ i and r i .

Accordingly, the moments of

X 1 , • • • , X n are expressed in terms of the elements of µ X = (µ 1 , • • • , µ n ),
V, and R.

So to sum up, the non zero covariances in (2.5) are written in terms of the non central moments of X 1 , • • • , X n . In addition, the moments of X 1 , • • • , X n are expressed in terms of the elements of µ X , V, and R. Thus V ar(Y ), which is known, is expressed in terms of µ X and the elements of V and R that we are searching for.

CASE 2:

If F is not a multivariate polynomial, then by Weierstrass Theorem F can be approximated by a multivariate polynomial [Reimer, 2012]. This well known theorem states that any continuous function on a closed and bounded domain in R n is approximated by a multivariate polynomial. So whenever the model F is continuous, we may select a bounded and closed subset of D and approximate F to a polynomial on this set. In our method, a bounded and closed subset of D is sufficient to do a backward propagation, since we are using a finite number of input points of F .

So theoretically, F can be approximated by a polynomial, however it is important to see how this is done practically. One of the simplest ways is to use the multivariate version of Taylor's expansion. However such polynomial expansion is done locally, in a small neighborhood of an input point. In addition, it requires the computation of the partial derivatives of F , which are used as the coefficients of the polynomial. More sophisticated and general approximations are the multivariate polynomial interpolation methods. These methods use a finite set of points of the model to construct a polynomial approximation that matches the model at the given points. They mainly rely on some polynomial basis to construct the approximations. For instance the Lagrangian interpolation, which is initially defined for the case of one variable, is generalized to the multivariate case. In [Sauer et Xu, 1995], the algorithms of two methods to find the coefficients of the polynomial approximation using multivariate Lagrangian interpolation are presented. In [Duchoň, 2011] a multivariate polynomial approximation using the Bernstein basis is proposed as a generalization of the univariate case.

Note that approximating a function by a polynomial in the univariate case is much simpler than the multivariate case. For that most methods are first developed for the univariate case, and then generalized to the multivariate case.

Thus, whenever the model F is continuous, it is approximated by a polynomial, and so we are again in case 1.

Therefore as a conclusion of the above two cases, whenever F is a continuous function, the variance of the output can be expressed in terms of µ X , V and R as

V ar(Y ) = G µ X , V, R (2.10)
for some algebraic function G. µ X which is the best estimate of the input point is actually the input point a. However the elements of V and R are unknown. The expression (2.10)

is the key point in generating the residues of the least square problem that is solved in step 2. Figure 2.1 is a flowchart that summarizes all the procedures of Step 1. In the flowchart, the term poly(.) refers to a polynomial of the inside indicated variables, and M i refers to the moments of the random variable X i of any order.

Step 2: Solving nonlinear least square problem

In order to do a backward propagation of uncertainty we should find the elements of V and R. So here we have n unknowns. However, we just have one equation involving these unknowns which is (2.10). To cope wih this problem we consider the uncertainty at different input points in D.

Note that, input uncertainty usually comes from some inaccurate measurement instruments or improper experimental procedures. So the uncertainty at different input points is due to carrying the same experiments and procedures with different configurations.

Thus, it can be assumed that the input uncertainty is a function of the input point. To Y = F (X) As it can be seen in this example the input uncertainty is a linear function of the input value W m . So with this example it is clear how input uncertainty can be expressed as a function of the input point. Accordingly, we define the input uncertainty at any point

F polynomial F non polynomial Y = q j Q j (X 1 , • • • , X n ) V ar(Y ) = q j q k Cov Q j , Q k Cov Q j , Q k = 0 Cov Q j , Q k = poly(M 1 , • • • , M n ) V ar(Y ) = q j q k poly(M 1 , • • • , M n ) M i = poly(µ i , V i ) M i = poly(µ i , r i ) V ar(Y ) = G µ X , V, R Case 1 Case 2 approximation Q j , Q k independent Q j , Q k dependent X i normal X i uniform
x = (x 1 , • • • , x n ) ∈ D as a functions of x:        V(x) = {V i (x i ) | i ∈ I 1 } R(x) = {r i (x i ) | i ∈ I 2 } (2.11)
So to find the input uncertainty at all the input points, we should find the function representations V i 's and r i 's.

In this manner, a first approximation of these functions can be derived, by assuming these functions as constant. This means that the quantities of the uncertainty at different input points are equal and independent of the chosen input point i.e. V(x) = V and R(x) = R.

It is to say that the same quantity of uncertainty arises whenever the measurements or the experiments are carried out. So to find the input uncertainty we should find the elements of V and R, which are n constant quantities corresponding to the input uncertainty for any x ∈ D. For that we select n different points in D denoted a (1) , • • • , a (n) . Then by step 1, the corresponding output variances of the selected points are decomposed into :

                 V ar(Y a (1) ) = G a (1) , V, R . . . V ar(Y a (n) ) = G a (n) , V, R
(2.12) This is a system of n unknowns, and its solution gives the values of the elements of V and R. To solve (2.12) we consider the least square problem: min

V,R n i=1 ∇ i || 2 (2.13)
whose residues that should be minimized are

∇ i = V ar(Y a (i) ) -G a (i) , V, R (2.14)
By solving this least square problem we obtain the values of the elements of V and R,

and hence the values of the input uncertainty. Note that the problem (2.13) is non linear, thus it is best to be solved numerically.

To be satisfied by this first approximation of the input uncertainty, which is represented by the obtained values of the elements of V and R, we use the following test. We select new points from D the domain of F and a tolerance, ε = 10 -6 for instance. Then we compute the residues of the new points using the obtained values of the elements of V and R and the formula: If for some x ∈ D, we get ∇ > ε, then the obtained uncertainty is not promising. For that we resolve a least square problem, however this time we relax the assumption that the uncertainty at different input points is constant. We assume that the uncertainty is a linear function of the input point. So for an arbitrary input point

∇ x = V ar(Y x ) -G x, V,
x = (x 1 , • • • , x n ) ∈ D

Applications

In each of the following two examples we perform a reverse computation of the input uncertainty. So we take the model F , which has in both cases two variables, and we specify the input uncertainty by the associated random variable of each input. Then we run a Monte Carlo simulation in order to get the data of the uncertain output. After that, our method is applied on the obtained data to find again the input uncertainty. The results obtained are then compared to the initially chosen values of the inputs uncertainty.

In addition, we present the convergence of the results as a function of the number of the generated data points used in the Monte Carlo simulation.

First example

Consider the model F to be a polynomial of two variables defined as:

F (x 1 , x 2 ) = x 1 * x 2 (2.18)
The input uncertainty of x 1 is assumed to have a normal distribution with variance V = 0.7, and the input uncertainty of x 2 is assumed to have a uniform distribution such that the radius of its support interval is r = 0.25. Then, using these assigned values the output data is generated at two input points a = (0, -2) and b = (3, 6). Then, based on the obtained data V ar(Y a ) and V ar(Y b ) are computed. To perform step 1 of the method, we use the expression of F to derive an expression of the output variance. Indeed, the output variance is written as:

V ar(Y ) = V ar(X 1 * X 2 ) = E[X 2 1 * X 2 2 ]-E[X 1 ] 2 * E[X 2 ] 2 = E[X 2 1 ] * E[X 2 2 ]-µ 2 1 * µ 2 2 (2.19)
Using the moments of the normal and the uniform distributions, and with some computations, V ar(Y ) becomes:

V ar(Y ) = 1 3 (V + x 1 )(3x 2 2 + r 2 ) -x 2 1 * x 2 2 (2.20)
In step 2 of the method, we use the expression of the output variance to derive the residues of the least square problem. Thus, here we first assume that the uncertainty is constant over the input range. Then we derive the residues ∇ a and ∇ b using equation (2.20):

∇ a = V ar(Y a ) - 1 3 (V )(3(-2) 2 + r 2 ) -0 ∇ b = V ar(Y b ) - 1 3 (V + 3)(3(6) 2 + r 2 ) -(3) 2 * (6) 2
(2.21)

Then we solve the least square problem using these two residues. The results obtained for V and r as a function of the number of sample points are presented in Fig. 2.3. As it can be seen from Fig. 2.3, the method detects the input uncertainty as we get very near values. However as the number of the data points increases the accuracy of the obtained values increases.

In the following we consider another example which is already considered in a previous publication to compare the results of our method to already existing methods.

Second example

In this example we consider a model presented in [Chen et al., 2015], that is defined by:

       F 1 (x 1 , x 2 ) = x 2 1 + e x 2 F 2 (x 1 , x 2 ) = sin(x 1 ) + x 2 (2.22)
So this model takes two inputs x 1 and x 2 and gives two outputs. According to [Chen et al., 2015], the uncertainties of both inputs are considered to be uniformly distributed and it is given that at the point a = (1, 0.7) the radii of the associated intervals of distributions are r 1 = 0.867 and r 2 = 0.6922. Starting from this, we generate the data of the uncertain output, and then this data is used by our method to find again r 1 and r 2 . In [Chen et al., 2015] the method used to detect the input uncertainty is the reversal workflow solver, which is restricted to models assembled as a computational workflow.

Concerning the first step of our method, we find a polynomial approximation of order 4

for each of e x 2 and sin(x 1 ). Then using these approximations we decompose the output variance in terms of x 1 , x 2 , r 1 and r 2 .

Concerning the second step of our method, we consider only the input point a given above.

This is because we have two outputs, giving two equations, with two unknowns r 1 and r 2 .

Thus substituting one point in the decomposition is enough to find the solution.

The following plots represent the values obtained of r 1 and r 2 as a function of the number of samples of the output data used: So our method detects the input uncertainty by giving very near values for r 1 and r 2 . The results obtained validate the efficiency of our method, however the accuracy depends on the number of sample points used.

Conclusion

In this chapter, we establish a new method for the backward propagation of uncertainty.

The method aims at quantifying the input uncertainty starting from the data of the uncertain output. We use the probabilistic representation in order to derive the input uncertainty. To clarify the concepts of the method we split it into two steps. In a first step we derive an expression of the variance of the output in terms of the statistical moments of the inputs. Then using this expression, in a second step, we construct a system of equations whose unknowns are the uncertainties of the inputs. The obtained system is nonlinear and it is solved numerically as a least square problem giving a quantification of the input uncertainty. Two applications were considered to validate the method, and the results showed that the method gives the values of the input uncertainty. However, the accuracy of the results depends on the number of the sample points used in the data of the uncertain output. In general, the method is considered as simple as it relies on some concepts from probability theory, however the decomposition of the output variance becomes more difficult as the number of the variables increase. In such a case, sensitivity analysis can be used in order to detect the most important inputs and hence involve them in the uncertainty study. In the next two chapters we consider one of the most important sensitivity analysis methods, called Sobol method, and we develop the way it is applied.

CHAPTER 3

Introduction

Sensitivity analysis methods are mainly used to detect the impact of the inputs on the output in a model. This assists in detecting the most influencing inputs. Accordingly it is possible to determine whether an uncertain input will cause a significant uncertainty in the output or not. Thus in an uncertainty analysis one may restrict the study to only the most influencing uncertain inputs. Several sensitivity analysis methods were revised in chapter 1, in which the Sobol method is considered as the strongest among them.

The Sobol method is based on the derivation of sensitivity indexes, called Sobol indexes, that indicate the impact of the inputs on the output. The computation of these indexes is done either symbolically according to the formula (1.51) or numerically using a simulation technique. To compute the indexes symbolically the expression of the model F should be given, and to compute them numerically many data samples of the model are needed.

The goal of this chapter is to examine the performance of the Sobol method in case of having no evidence about the expression of the model F plus having a limited number of data points coming from experiments.

To reach our goal we apply the Sobol method to study the sensitivity of an experimental curve, called Electrostatic Force Distance Curve (EFDC) [Villeneuve-Faure et al., 2014],

with respect to some experimental factors. The following sections present in detail this sensitivity study. First, we introduce the concept of the EFDC and how it is obtained.

After that, in section 3.3 we indicate the experimental factors that are considered as the input of the model and that their impact on the EFDC is studied by the Sobol method.

Then in section 3.4, we describe the experimental procedures and materials that are used to acquire the data of the EFDCs. This is important in order to understand the sensitivity results from a physics point of view. In order to apply the Sobol method, the output of the considered model should be a scalar. For that, the EFDC is fitted by a 4 parameter logistic law in section 3.5. This enables us to investigate the sensitivity of the EFDC by studying the sensitivity of the associated four logistic parameters and the sensitivity results obtained are presented in section 3.6. Then, to validate the results we present an effect detecting technique used in experiments called Design of Experiment (DOE) which is introduced in section 3.7. Lastly in section 3.8, we construct an approximation model, called the matrix model, for the four logistic parameters giving an approximating formula for the electrostatic force.

What is an EFDC

To best describe EFDCs, it is important to know first how these curves are traced. Indeed, the main technique used in this manner is the Atomic Force Microscopy(AFM). AFM is a kind of Scanning Probe Microscopes (SPM), which are microscopic techniques designed to detect the local properties of materials, at a microscopic scale, by scanning with a probe.

AFM has been initially developed as an imaging technique to study the topography of surfaces and materials. However, its ability to detect and measure the interaction forces between the probe and the sample surface makes it a powerful tool to disclose much more than just the chemical and mechanical properties of the sample. In this manner, AFM is used to find the electrostatic force caused by charge accumulation in dielectrics. In this section we explain how AFM is used to compute the electrostatic force and how the EFDCs are extracted, readers interested in more details about the general procedures of AFM are referred to [Cappella et Dietler, 1999].

The AFM process

The AFM operator consists of three main parts: a probe, a detector, and a scanner (see Fig. 3.1):

• The probe is formed of a very sharp thin tip (typically less than 5 µm tall and often less than 10nm in diameter at the apex [Cappella et Dietler, 1999]) attached to the free-swinging end of a small spring-like cantilever that is usually 100 -500µm long.

• The detector records the deflection and the motion of the cantilever as the tip scans the sample. Usually, a laser beam with a photo-detecto are used as a detector (as in Fig. 3.1), however other deflection detecting methods can be used as Piezoelectric detection [Giessibl, 1998], Optical Interferometry [START_REF] Rugar | Improved fiberoptic interferometer for atomic force microscopy[END_REF], and Scanning

Tunneling Microscope [START_REF] Binnig | Atomic force microscope[END_REF]].

• The scanner controls the probe-sample displacement both vertically and laterally, in order to allow the probe to scan the sample in all directions. As illustrated in Fig. 3.1, the drive is attached to the sample, however it can be also attached to the cantilever while keeping the sample fixed.

There are two principle modes in which AFM operates: the contact mode and the dynamic mode [START_REF] Binnig | Atomic force microscope[END_REF]. This depends on how the probe scans the sample. For acquiring the EFDCs, the contact mode is used, in which the tip makes a soft contact with the sample at a particular location of the sample surface. In this mode, the process starts by the tip and the sample at a rest state in which a large distance separates them so that no interaction forces are presented (Fig. 3.2(a)). Then the scanner starts to move vertically upward, bringing the sample very close to the tip.

As the sample approaches the tip, the cantilever remains at an equilibrium state, until the sample comes close enough to the tip so that the tip experiences the attractive Van der Waals force. So the tip snaps into surface forming the jump to contact point (Fig.

3.2(b))

, and the cantilever bends slightly towards the surface. As the scanner continues to move upward, the cantilever deflects away from the surface (Fig. 3.2(c)). Whenever the Van der Waals force force is detected, implying that the tip is in contact with the surface, the scanner begins to retract. However, the interaction forces between the tip and the sample hamper the retraction, making a gradual withdrawal of the tip from the sample (Fig. 3.2(d)). Lastly, the tip withdraws and loses its contact with the sample (Fig. 3.2(e)). These different phases of the AFM are illustrated in the Fig. 3.2.

During all these phases, the deflection of the cantilever is recorded by the detector. Note that this deflection is due to the net vertical forces acting on the tip, which are in fact the interaction forces between the tip and the sample. So the quantity of these forces at each moment can be computed using the Hooke's law:

F = k δ c
where k is the spring constant of the cantilever defined by the manufacturer, and δ c is the deflection of the cantilever from its equilibrium position calculated using the records of the detector. Plotting F as a variation of the sample-tip distance gives the well known Force Distance Curve (FDC). The FDC provides direct measurement of interaction forces between the AFM tip and the sample surface. Note that the tip sample interaction forces may include other forces than just the Van der Waals force. This depends on the materials used as well as the supplement settings added to the AFM setup. For instance, if the tip is grounded, and a dielectric sample is scanned by the tip, in which charges are injected on the surface of the sample or an electrode is buried in the sample, which induces an electrostatic force.

Indeed, in the case of a buried electrode, when a voltage is supplied to the electrode, the electrical potential difference between the electrode and the tip induces an electrostatic force. However in the case of injected charges, the charge density is trapped in the dielectric layers, and this charge density induces electrostatic force on the tip. This induction of the electrostatic force in both cases affects the tip-sample interaction during the different phases of the AFM operation. This idea plays a basic role in deriving the EFDC, which is our concern in this work.

Acquiring the EFDC

In fact, EFDCs are obtained by two steps. The first step is to record the FDC over the considered sample. The second step is to record the FDC again after either injecting some charges in the sample or supplying a buried electrode in the sample by a potential.

Then the EFDC is obtained as the difference between the two obtained FDCs in the approach and contact stage [Villeneuve-Faure et al., 2014]. This implies that the EFDC represents only the action of the electrostatic force on the tip-sample position, without the consideration of the other forces. Fig. 3.4 shows the two obtained FDCs while scanning an oxy-nitride dielectric. The black curve corresponds to a blank FDC, however the gray one corresponds to the FDC after injecting charges on the surface of the oxy-nitride sample.

Charges are injected on the surface of this dielectric for obtaining the second FDC. The insert in the graph is the difference between the two FDCs curves which represents the EFDC. So as we see the EFDC is based on the measurements of the electrostatic force between the AFM tip and the sample. In the next section we indicate the motivation for studying the sensitivity of the EFDC from the physics point of view. This enables us to indicate the experimental factors that are of scientific interest in this manner. These factors will be the inputs of our considered model, in which their impact is to be studied using the Sobol method.

The model inputs

All solid dielectrics have the property of being able to accumulate electrical charges under electrical stress beyond an electric field threshold [Villeneuve-Faure et al., 2014]. The accumulation of electrical charges induces a local increase of the electric field, and this may lead to a failure of the system containing the dielectric and/or a premature of dielectric breakdown [Normand et al., 2003]. Consequently, it is important to quantify the density of accumulated charges and their localization in dielectric layers to improve the reliability of devices and systems.

Several methods for charge detection in insulators have been proposed; see [Rezende et al., 2009] for a review of these methods. In the same manner, in [Villeneuve-Faure et al., 2014] a recent method has been proposed that totally relies on the Electrostatic Force Distance Curve (EFDC), from which the method takes its name as it is also called the EFDC method. This method aims to detect the 3D localization of charges across dielectrics [Boularas et al., 2016]. It is specialized for thin dielectric films with a thickness of less than 200 nm.

It has been demonstrated experimentally that EFDCs are sensitive to charge localization in three dimensions [START_REF] Alhossen | Sensitivity analysis of the electrostatic force distance curve using sobol's method and design of experiments[END_REF]. So studying their sensitivity can be a first step in detecting the charge spatial localization in dielectrics. For that, we study the sensitivity of the EFDC with respect to the charge position and density in dielectrics using the global sensitivity analysis method of the Sobol method.

In our case study, the considered EFDCs are obtained by AFM scanning a dielectric with a buried electrode. Indeed, from an experimental point of view the potential is easier to be controlled than charge distribution. Accordingly, the electrical potential localization and density can be represented by three elements: the width of the electrode w, the depth of electrode d, which in turn represents the electrode's position, and v the potential applied on the electrode. So the inputs of our model are w, v and d of which we study the impact on the EFDC using Sobol method. In the next section, we give more details about the experimental procedures and the materials used to derive the EFDCs that we use in our sensitivity study. To manufacture these samples, a 270nm-thick SiN x dielectric layer was deposited using High Frequency Plasma Enhanced Chemical Vapor Deposition [Zaghloul et al., 2010] over a highly doped silicon substrate. For the embedded aluminum electrode fabrication, a lift-off process was used. First, a 2.5µm-thick N-LOF photoresist was deposited and patterned by photolithography. After this step, SiN x was chemically etched to a depth of 70nm and filled with aluminum. The combination of etching and electrode deposition with the same photoresist layer ensures an intimate metal/dielectric contact and a small surface roughness. Finally, the electrodes were embedded at different depths d depositing a SiN x cover layer of different thickness.

Experimental procedures

The AFM measurements were done using a Bruker Multimode 8 apparatus. To avoid capillarity parasitic effect all measurements were done under N 2 atmosphere after drying the samples at 100 • C for 15 min to remove the water layer. In addition, these measurements were done using a Pt-coated silicon tip provided by Bruker (SCM-PIT) for which the spring constant of the cantilever was calibrated combining photodiode sensitivity measurement and thermal tune mode. Moreover, all measurements were performed at the middle of the electrode width. Moreover, during the experiments the cantilever remains parallel to electrode to minimize parasitic effect [Negoescu et Axinte, 2007a].

The different values taken by w, v and d to obtain the samples are:

• w = 6 µm, 20 µm, 40 µm.

• v = 4 V, 6 V, 8 V, 10 V, 15 V.

• d = 10 nm, 50 nm, 100 nm.

Accordingly, 45 different samples were scanned. Implying that 45 different EFDCs were obtained. These EFDC are used in our sensitivity study.

Fig. 3.6 shows some EFDCs obtained for different configurations according to the above explained experimental procedures. In general, these curves start from their negative initial values and then undergo a logarithmic growth. There is rapid progress at the beginning, but continue to progress slowly until they approach their horizontal asymptotes. As it can be seen from these two figures, the variation of d and w, respectively, leads to an apparent modification in the plot of the EFDC. Even though this influence is clearly visible from these plots, its quantification is not straightforward. So quantifying this influence could be a key point in using the EFDC for detecting the charge position in 3D [Palleau et al., 2010]. With this in mind, we use the Sobol sensitivity method to thoroughly understand and quantify this influence. However, to apply the Sobol sensitivity method, the output should be a scalar. For that, in the next section a representative model of the EFDC is introduced, in order to prepare the requirements for applying the Sobol sensitivity method.

EFDC as a logistic law

Recall that from chapter 1, the Sobol method is a global sensitivity analysis applied to study the sensitivity of a model with respect to its input variable. In this manner, the output of the studied model should be a scalar. So for studying the sensitivity of the EFDC with respect to d, v, and w, it is needed to find a representative model that takes d, v, and w and link them to the EFDC, while giving a scalar output.

Since we are dealing with graphs, at first sight one would think of consulting analytical geometry to find the link between the graphical features of the curve and the values of d, v, and w, for instance take the curvature as an output. However, the data of the curves is derived from experimental procedures, and hence it is surely infected by random noise.

This leads to EFDCs with almost perturbed plots, and so it would be problematic to carry a geometric procedure directly. For that, a more systematic way to deal with such noisy curves is used.

The clue here is that these curves have a sigmoid shape compatible with the logistic law curves. So, according to [Villeneuve-Faure et al., 2014], the best experimental data fitting for an EFDC matches a 4 parameter logistic law (or 4PL for short), which is defined by the following expression: The main advantage of fitting by the 4PL law compared to more general laws (like polynomials) is that the 4PL fits accurately the EFDC over a large distance range. The plots in 3.8 compare the regression residues of fitting the EFDC by 4PL and by a polynomial for different degrees. In the ordinate, we have the residuals after fitting (using least-square) the experimental EFDCs with polynomial expression and 4PL. In the abscissa, we have the degree of the polynomial used in the fitting.

f (z) = D + A -D 1 + z C B ( 3
Residuals coming from polynomial regression are always higher than those coming from logistic Law regression. This demonstrates that the logistic Law is better than polynomial expression to fit the EFDCs. Regression using polynomial order higher than 13 has not been performed since the fitting algorithm was not able to converge. Moreover, using polynomial expression of order more than 13 to fit the EFDCs curves does not make sense physically. For that, we choose the 4PL as a representative formula for the EFDC.

After fitting the experimental EFDCs by 4PL curves, for each value of the triplet variables Figure 3.9: The model Fitted-4PL

First order Sobol indexes

The general formula of the first order Sobol index of an arbitrary input x i is given by:

S i = V ar(F i ) V ar(Y ) (3.2)
where V ar(F i ) is defined by:

F i (X i ) = E[Y /X i ] -E[Y ] (3.3)
For each output A, B, C, and D and for each of the inputs d, v, and w, the first order Sobol indexes are computed by substituting the above formulas according to a discrete random variables. So for instance, for the output A and the input d we have:

S d = V ar(A d ) V ar(A) (3.4)
where According to these first results, the applied bias v seems to be the main influence for the parameters on of EFDC curve [START_REF] Alhossen | Sensitivity analysis of the electrostatic force distance curve using sobol's method and design of experiments[END_REF]. This highlights the Sobol method interest in this sensitivity study

V ar(A d ) = V ar(E[A/d] -E[A]) = V ar(E[A/d]) = V ar(E[A/d = 10], E[A/d = 50], E[A/d = 100]) = 1 3 E[A/d = 10] -E[A] 2 + E[A/d = 50] -E[A] 2 + E[A/d = 100] -E[A]

Second order Sobol indexes

The general formula of the first order Sobol index of an arbitrary input x i is given by: to have no effect on A as all second order Sobol indexes are approximately zero.

S ij = V ar(F ij ) V ar(Y ) (3.6)
To confirm the results of the Sobol indexes as computed by the formulas of the discrete random variables, we use Design of Experiment (DOE) that is an influence detecting method specialized for experimental model. In the next section we introduce briefly the notion of DOE, then we present the results obtained by using DOE to detect the effects of d, v, and w on A, B, C, and D as stated in [START_REF] Alhossen | Sensitivity analysis of the electrostatic force distance curve using sobol's method and design of experiments[END_REF].

Design Of Experiment (DOE)

To validate the results obtained above using Sobol method we use a graphical way from Design Of Experiments (DOE). In general, DOE is an effective strategy to examine the behavior of a simulation model when changing its input values in order to detect their effect [Wakeland et al., 2004]. This helps in understanding the input-output relationship and in identifying the relative importance of the inputs. DOE has a property of expressing the results graphically, which allows scientists to directly analyze of the effect of each input. Hence, applying a DOE technique to the present case study, allows the identification of the impact of each of the inputs d, v, and w on the logistic parameters (outputs) A, B, C, and D.

The basics of DOE

To start, first we introduce some fundamental terms used when dealing with DOE techniques as mentioned in [Telford, 2007] Several techniques of DOE have been developed [Cavazzuti, 2012], each is applied to some particular problems depending on the required purpose of the experiment. Among those, Factorial Design is considered as one of the most efficient experimental designs for assessing the effects of the experimental factors on the response. In particular, a Full Factorial Design of experiment allows not only the detection the individual effects of the factors on the response ( first order Sobol indexes), but also the interaction between the factors and its effect on the response ( second order Sobol indexes). By this DOE technique, all the possible combinations of the selected levels of the factors are taken in account.

Although this collection could be costly for complex experiments, however it enables the broad investigation of the effects of the factors as well as their interactions. Thus, results

obtained by Full Factorial DOE are expected to be more robust than the results obtained by any fraction Factorial DOE, in which a fractional number of combinations of the factor levels is taken.

The data executed from the full Factorial DOE is usually interpreted graphically using the Factorial Plots. These plots mainly depend on showing the contrasts of the averages of the response at different configurations of the factors to derive conclusions. The two frequently used factorial plots are the Main Effect Plots and the Interaction Plots.

The main effect plots are the factorial plots responsible for visualizing the individual effects of the factors on the response. The interaction plots are the factorial plots responsible for visualizing the effects of the interaction of the factors on the response. In the next subsections we detail how these plots are done and the associated factorial plots of our study case.

The main effect plots

The main effect plots are formed by plotting the different averages of the response while fixing one of the factors at a certain level. To see this practically, consider the Fig. 3.12 which displays the main effect plots of the logistic parameters A, B, C, and D with respect to the factors d, v, and w. Analyzing the main effect plots of Fig. 3.12 we see:

• Fig. 3.12(a) shows that parameter A (i.e. the minimum value of the electrostatic force) (i) increases as the depth d increases, (ii) decreases quadratically as the bias v increases and (iii) decreases slowly as the width w increases.

• Fig. 3.12(b) shows that parameter B (curve bending): (i) increases slowly as the depth d increases, (ii) decreases with as the bias v increases and (iii) is quite constant with respect to the width w.

• Fig. 3.12(c) shows that parameter C increases slightly as the depth d increases and the applied voltage v increases. However C seems to vary as the the width w increases. But it is important to notice that the variation of C remains small until the voltage reaches around 12v.

• Fig. 3.12(d) shows that parameter D (the maximum value of the electrostatic force) increases as the depth d increases, and decreases as the bias v and the electrode width w increase.

In general, these results validate the results obtained by Sobol sensitivity method. As it appears, A is mostly affected by d and v, and this is compatible with the bar diagram 3.10. Also, compatible results obtained for D as it approximately equally affected by d, v, and w. Moreover, the main plots shows that B is mostly affected by v and this is in agreement with the obtained Sobol index of v which is around 0.8. Similarly for C, which appears to be mostly affected by v. Now we proceed to check the compatibility between the interaction plots and the second order Sobol indexes.

The interaction effect plots

The interaction effect plots indicate whether the effect of a factor on the response changes depending on the setting of another factor. These plots display by lines all the main effect plots of one factor at every level of another factor. The connecting lines in every plot lead to the understanding of how the interactions between the factors affect the response. If the lines are parallel, this means that the effect of the indicated factor on the response is identical across all the levels of the other factor. Hence, the effect does not depend on the level of the other factor, and so no interaction occurs. However, when the lines are not parallel, the effect of one factor depends on the setting of the other factor, and so in this case there is an interaction effect.

To see this practically, consider the matrix plot in Fig 3 .13 which represents the interaction effects matrix plot of the parameter A. This is a multi-plot per figure, displaying the original main effect plots of A on the diagonals, and all of the two-factor interaction effect plots of A on the off-diagonal positions. For instance, the first plot in the first row presents the main effect plot of d on A. On the other hand, the second plot in the first row presents six traces for the effect of d on A, each trace corresponds to a level of v.

Analogously the other plots can be seen.

In general, despite the slight non parallelism in some interaction plots of A, most traces possess approximately parallel lines. This shows that no serious interaction between d, v,

and

Matrix Model

To continue investigating the relationship between the logistic parameters A, B, C, and , 2007b] that allows us to find a multivariate polynomial approximation using little experimental data. This is favorable in our case since we have a limited number of data points for A, B, C, and D (only 45).

D
In addition, the coefficients in this model are computed by formulas similar to those that express the mean average in the discrete case, which we have acquired from the computation of the Sobol indexes and the effects in DOE.

The matrix model is expressed with different orders, depending on how many contributions of the inputs are needed in the approximation formula. So for instance, the first order matrix model considers only the individual contributions of the inputs and it is formed of the sum of univariate polynomials of the inputs. Second order Matricidal model includes the second order interaction of the inputs and it is formed of the first order matrix formula plus bivariate polynomials of the inputs. Analogously higher order matrix models are defined.

Since three factors are involved in the experiments of our study case (d, v, and w) then a general response symbolized by y is considered, written as y = M(d, v, w) for some unknown formula M. M here represents the unknown expression of the model to be approximated which is in our case the logistic parameters A, B, C, and D.

The following subsections represent the results obtained when applying a first and a second order matrix approximations for A, B, C, and D in terms d, v and w [START_REF] Alhossen | Sensitivity analysis of the electrostatic force distance curve using sobol's method and design of experiments[END_REF]. This includes a confirmation of the results obtained by DOE and Sobol method.

First order matrix model

The First order matrix Model M 1 which approximates M is defined by:

M 1 (d, v, w) = ȳ + I d 1 (d) I d 2 (d) • • • I d l (d)           E d 1 E d 2 . . . E d l           + I v 1 (v) I v 2 (v) • • • I vm (v)           E v 1 E v 2 . . . E vm           + I w 1 (w) I w 2 (w) • • • I wn (w)           E w 1 E w 2 . . . E wn           (3.10)
where ȳ is the total average of the response, and l, m, and n are the number of levels of d, v, and w respectively. I d , I v , and I w are polynomials associated to the levels of the factors and called indicators. For instance, the indicator of the factor d at level i, denoted by I d i (d), is a polynomial in d and it is given by:

I d i (d) = k =i d -d k d i -d k (3.11)
where d k is the kth level of d in the given data. Moreover, in formula (3.10) E d , E v , and E w are constants associated to the levels of the factors and called the first order effects.

Practically, the first order effect of the factor d at the ith level, denoted by E d i , can be easily computed using the following formula:

E d i = Mean of y when d is fixed at d i -Total mean (3.12)
According to [Pillet, 1997] As it can be seen in Fig. 3.15, a good agreement between the two plots of the EFDCs is found but not for all cases. Indeed, disagreement is observed for the maximum attained asymptotic (Fig. 3.15(c)) or at the bending of the curve (Figure 7.d). These two discrepancies are related directly to D and the couple B-C respectively. They show that the approximation model is missing some information. By referring to the results of the interaction plots of DOE (Fig. 3.14) and the second order Sobol indexes (Fig. 3.11), this defect can be justified. Indeed, the interaction plots demonstrate the importance of the interaction between d, v, and w and its influence on the logistic parameters especially B and D. Thus, to avoid this drawback in the approximation, the interaction factors should be involved in the formula. This is done by considering a second order matrix model.

Second order matrix model

The second order matrix model add second order interactions of the inputs to the first order matrix model. This is done by adding to formula (3.10) new terms of the following form:

I d 1 (d) I d 2 (d) • • • I d l (d)           E d 1 v 1 E d 1 v 2 • • • E d 1 vm E d 2 v 1 • • • . . . . . . . . . . . . E d l v 1 E d l v 2 • • • E d l vm                     I v 1 (v) I v 2 (v) . . . I vm (v)           (3.13)
where E d i v j is the interaction effect between d and v at levels i and j respectively. Its formula is given by:

E d i v j = Mean of y when d and v are fixed at d i and v j -Total mean -E d i -E v j (3.14)
Analogously to equation (3.14), the interaction effects E d i w j and E v i w j are defined. As it can be seen from Fig. 3.16, the approximated EFDCs computed when using the second order model exhibit a very good agreement with the experimental EFDCs. These results are promising as they are compatible with what obtained suing DOE and Sobol method.

Usually in approximating expressions, attention should be paid to the cost of each interaction term added to formula (3.10). So compensation usually is taken between the efficiency of the constructed model and the computational cost. For the EFDC measurements, the price of getting the 45 trials, compared to the results of the second order matrix Model, is acceptable according to the experts. The second order Matricidal model permits a precise computation of EFDCs for whatever values of (d, v, w), which will be powerful for charges localization investigation in thin dielectric film [Boularas et al., 2016].

Conclusion

In this chapter we check the use of Sobol indexes if they are computed using the formulas of discrete random variables. Usually when the Sobol method is applied, Sobol indexes are computed numerically and the sampling is taken according to continuous random variables (normal, uniform, • • • ). In this chapter we show that the Sobol method may give robust results even for sampling of discrete random variables.

To this end, we study the sensitivity of the experimental curve EFDC from the domain Moreover, we propose an approximation model, called matrix model, of the logistic parameters of A, B, C, and D in terms of d, v, and w. Using this approximation we show the importance of taking into account the interaction between the inputs, while comparing a first order approximation by a second order one relative to the experimental results.

The final words concerning the application done in this chapter are to take advantage of the obtained results in the domain of force spectroscopy. According to experts, the constructed matrix model can be used to characterize the charge localization as it in-cludes the three main variables: the depth d and the width w, hence the position and volume of the charge, and the voltage v, meaning the density of the charge.

In the next chapter we will continue with the Sobol method, but we extend its basic application into studying the variation of the Sobol indexes with respect to some active factors or experimental conditions. Accordingly, new conclusions and scientific consequences will be derived. in a particular manner and we examine its performance in the discrete case. However here we consider it in a general manner and we study the variation of the Sobol indexes with respect to some external factors or experimental conditions. The aim of studying this variation is to reveal the optimal settings at which conclusions and information about the inputs can be derived. This helps in the case of uncertainty study as it indicates the best configurations for analyzing and quantifying of the input uncertainty. In this manner, we consider two different models and we study the evolution of their corresponding Sobol indexes with respect to some factors, in order to see how conclusions about the inputs can be derived.

The first considered model is from the domain of computer vision. It represents mainly the functioning of a 3D reconstruction method called Shape-From-Template (SfT) [START_REF] Bartoli | Shape-from-template[END_REF]. This method uses a single 2D image and a 3D template to reconstruct a 3D deformed surface. An important factor in the 3D reconstruction process is the depth of the 3D deformed surface in front of the camera, so for that we study the evolution of the Sobol indexes concerning the SfT method as a function of the depth. To keep the concepts clear, we start by recalling in section 4.2 the mathematical model of image formation by a camera and the parameters involved in this process. Then in section 4.3, we present the details of the 3D reconstruction method SfT. In section 4.4, we introduce the inputs of our considered model according to the concern of the SfT method. After that in section 4.5, we explain how we construct our model while including the numerical solution of the SfT method and getting a scalar output to be adapted to apply Sobol method. Keeping in mind that the solution of the SfT is theoretical and that in reality 3D reconstruction is always accompanied by noise, for that we modify the constructed model by adding noise to imitate reality. The procedures done are detailed in section 4.6. To conclude, we describe in section 4.7 how Sobol indexes are computed as a variation of the depth and the conclusions derived from these variations.

The second model considered in this chapter is a charge transport model for dielectrics. This model describes how a dielectric may conduct charges under the effect of high electric fields. We deal with it here as a black box, and we study the variation of its corresponding Sobol indexes as a function of time, temperature and the strength of the applied electric field. In section 4.9, we define the inputs and the outputs of the considered charge transport model. In this manner, we will see how studying the variation of the Sobol indexes helps in deriving the best configuration for applying an optimization process to find the inputs.

Image formation: from 3D into 2D

In this section, we review the mathematical formulation of a camera in order to form the 2D image of a 3D object. First of all, it is important to note that the cameras used in our daily life are accompanied with lenses to focus light. However, all cameras follow the same concept as the pinhole camera in image formation. Thus the pinhole camera model is presented with the associated notations that are used throughout this manuscript.

A pinhole camera is a black box that is punctured from one side forming a small hole.

The rays of light coming from the outside world pass through the hole and fall on the opposite side of the box, forming a 2D image of the 3D outside scene. Observe that the real pinhole image is an upside down image of the scene. However the virtual image, which is usually seen on a photograph or on a computer screen, corresponds to the projection of the scene onto a hypothetical plane. This hypothetical plane is situated in front of the camera at the same distance from the hole to the opposite wall on which the image is actually formed (see Fig. 4.2).

The distance from the hole to the opposite wall on which the image is formed is called the focal length of the camera, and it is denoted throughout this manuscript by f . In addition, the 'hole' of the camera is called the camera center and it is denoted by C, whereas the To see geometrically how the 3D object is projected into the image plane, we follow the procedure explained in [Moons et al., 2008]. Consider an orthogonal reference frame

(D) = {(x = X M t, y = Y M t, z = Z M t) | t ∈ R} (4.1)
In addition, (D) intersects the image plane, and let m be the point of intersection. As m ∈ (D) then the coordinates of m in I are (X M t m , Y M t m , Z M t m ) for some t m ∈ R. But m also belongs to the image plane, thus its third coordinate should be equal to the focal length f , and so Z M • t m = f . This implies that t m = f Z M , and so the coordinates of

m in I are (f X M Z M , f Y M Z M ), and the coordinates of m in J are ( f Z M X M , f Z M Y M ).
Thus the image of the 3D point M is the 2D point m whose coordinates in the image plane are

( f Z M X M , f Z M Y M ).
Applying this to all the points (X, Y, Z) in the 3D scene, the 2D image is formed on the image plane and the corresponding coordinates of the points are (f X Z , f Y Z ). Accordingly, we define the projection function of a camera Π ∈ C ∞ (R 3 , R 2 ) defined for any M = (X, Y, Z) by: The change of units between J and P is done according to the following transformation:

Π(M ) = f Z (X, Y). ( 4
1 - → U = α u -→ U p 1 - → V = α v -→ V p (4.3)
for some positive values α u and α v . Let the coordinates of O in the pixel coordinate system P be (U po , V po ), then a point m having coordinates (f X Z , f Y Z ) in J, its pixel coordinates are derived as follows:

--→ O p m = --→ O p O + --→ Om = U po -→ U p + V po -→ V p + f Z X - → U + f Z Y - → V = U po -→ U p + V po -→ V p + f Z Xα u -→ U p + f Z Yα v -→ V p = U po + α u f X Z -→ U p + V po + α v f Y Z -→ V p = m 1 -→ U p + m 2 -→ V p (4.4)
So the coordinates of m in the frame P are U po + α u f X Z and V po + α v f Y Z . These coordinates can be derived in a compact way using a matrix expression:

Z        m 1 m 2 1        =        α u f 0 U po 0 α v f V po 0 0 1               X Y Z        (4.5)
The above matrix is called the characteristic matrix of a camera and it is usually denoted by K. In addition, α u f, α v f, U po , and V po are called the internal parameters of the camera.

In this section we see how a 2D image of a 3D object is formed using the camera parameters. Each point in the 3D object undergoes a simple geometric transformation to get its pixel coordinates in the 2D image:

M = (X, Y, Z) 3D camera frame m = (f X Z , f Y Z ) 2D image plane m = (U po + α u f X Z , V po + α v f Y Z )
2D pixel frame

In the next section, the reverse operation is considered i.e. we go from a 2D image into recovering the 3D object. This is a much more complicated process than the simple image formation.

Shape-from-Template: 3D reconstruction

3D reconstruction is the process of reproducing the shape and the appearance of a 3D object starting from acquired information. Mainly, the starting point is a given 2D image of the scene. So for a 2D point in the image plane having coordinates (a, b), 3D

reconstruction is to recover its corresponding 3D point (X, Y, Z) in the scene. Note that if the focal length of the camera f and the depth of the scene Z are known, then the 3D reconstruction problem can be directly solved. This is because the remaining unknowns are X and Y, and they can be derived using f and Z by:

X = Z f a and Y = Z f b (4.6)
So the 3D reconstruction problem can be defined as a problem of finding the focal length f and the depth Z. Note that, this is not an easy mission and it cannot be solved by a single 2D image. For that, this topic has occupied a broad research area in the domain of computer vision. Various methods exist in this manner [Moons et al., 2008]. Simple methods start by using multiple 2D images (at least two) to drive the depth geometrically, while other methods use videos instead. In some problems, the used camera is calibrated, i.e. its intrinsics are known, and so its focal length is known. In this case, the 3D reconstruction problem is reduced to the problem of finding the depth Z only. In some other problems, the scene is a non rigid object, i.e. it deforms and changes its shape with time.

This makes the problem more complicated in reconstructing the deformed object.

In this work our concern is the method Shape-from-Template (SfT for short) [START_REF] Bartoli | Shape-from-template[END_REF][START_REF] Salzmann | [END_REF]Moreno-Noguer et Porta, 2016 ;Agudo et Moreno-Noguer, 2017]. This method recovers a 3D deformed surface using one 2D image and a 3D template of the surface. This template is a 3D reference shape that under a specific deformation gives the 3D object. An example in this manner is the template • L ⊂ R 2 is the known flattening of the 3D template T into R 2 .

In addition, the labeled arrows represent functions that are characterized as follows:

• Π ∈ C ∞ (R 3 , R 2 ) is the known camera projection function.

• Ψ ∈ C 1 (T , R 3 ) is the unknown isometric deformation function.

• ∆ ∈ C 1 (L, R 3 ) is the known template embedding function, it is an invertible function that maps the 3D template to its flattening.

• ∆ -1 ∈ C 1 (T , R 3 ) is the known flattening function of the 3D template.

• η ∈ C 1 (L, R 2 ) is the known warp function mapping the flattened template into the 2D image.

• ϕ ∈ C 1 (L, R 3 ) is the unknown embedding into the deformed surface.

• γ ∈ C 1 (L, R + * ) is the depth function, mapping each point in the flattened template to the depth of its point in the 3D deformed surface.

Remember that to do a reconstruction, the depth Z is to be determined. According to the notation provided in Fig. 4.6, the function γ is to be determined. The SfT method reconstructs the 3D deformed surface based on a single image and a 3D template. This gives two constraints: the reprojection constraint and the deformation constraint.

The reprojection constraint represents the compatibility of the projection of the deformed surface with the given image data. This is expressed by:

Π • Ψ • ∆ = η (4.7)
The deformation constraint represents the restrictions of the deformation done on the 3D template. Here the focus is on isometric deformations, which means that the distance between any two points on the surface is preserved. This implies that:

J Ψ J Ψ = I 3 (4.8)
where J Ψ is the Jacobian matrix of the unknown function Ψ and I 3 is the identity matrix of size 3. The symbol raised on the Jacobian matrix refers to its transpose.

Note that the above two constraints form a first order system of partial differential equations whose main unknown is the function Ψ:

       Π • Ψ • ∆ = η (reprojection constraint) J Ψ J Ψ = I 3 (deformation constraint) (4.9)
It is a first order system since it includes J Ψ the first derivative Ψ. In addition, system (4.9) is a system of nonlinear partial differential equations of dimension 5: two equations from the reconstruction constraint and three equations from the deformation constraint.

This system is the starting point of the SfT method.

Indeed, the SfT method proceeds in finding the unknown depth function γ by reformulating system (4.9) two times, and then by doing a change of variable. The idea of these reformulations is to introduce γ as the unknown of the system instead of Ψ, and to get a simpler system with a lower dimension.

First reformulation

The first reformulation is done by using the equality:

ϕ = Ψ • ∆ (4.10)
The idea is to introduce ϕ that is defined on the 2D parameterized space L which is also the domain of the depth function γ. By (4.10) the reconstruction constraint becomes:

Π • ϕ = η (4.11)
Now, differentiating both sides of (4.10) gives:

J ϕ = (J Ψ • ∆)J ∆ (4.12)
Then multiplying each side of equation ( 4.12) by its transpose and using the deformation constraint implies:

J ϕ J ϕ = J ∆ (J Ψ • ∆) (J Ψ • ∆)J ∆ = J ∆ I 3 J ∆ = J ∆ J ∆ (4.13)
Thus system (4.9) becomes:

       Π • ϕ = η (reprojection constraint) J ϕ J ϕ = J ∆ J ∆ (deformation constraint) (4.14)
This new system (4.14) is also a non linear system of partial differential equations however of dimension 4: two equations from the first equality and two equations from the second equality. In addition, its main unknown here is ϕ instead of Ψ.

Second reformulation

In the second reformulation the depth function γ is introduced. Note that, whenever a point in the 2D image is multiplied by Z f , the first two coordinates of its corresponding 3D point, X and Y, are obtained. So if a point in the image of the function η is multiplied by Z f , the coordinates X and Y of the corresponding point in the deformed surface S are obtained. To express this symbolically, consider a new function η ∈ C 1 (L, J × {f }) defined by:

η = [η, f ] (4.15)
where f is the known focal length of the camera. Then the depth function γ is introduced using ϕ by the equality:

ϕ = 1 f γ η (4.16)
To substitute this new formulation of ϕ in the system (4.14), equation (4.16) is differentiated on both sides:

J ϕ = 1 f ηJ γ + γJ η (4.17)
So multiplying each side of (4.17) by its transpose gives:

1 f 2 η 2 2 J γ J γ + γ 2 J η J η + γ(J γ η J η + J η ηJ γ ) = J ϕ J ϕ (4.18)
Using the deformation constraint of (4.14), a new system is obtained, defined by:

η 2 2 J γ J γ + γ 2 J η J η + γ(J γ η J η + J η ηJ γ ) = f 2 J ∆ J ∆ (4.19)
This system is of dimension three. In addition, its main unknown is the depth function γ.

Solving this system gives the solution of the reconstruction problem. However, the non linearity presented in the system, appeared as γ multiplied by J γ , harden its solving. For that, a change of variable is done to eliminate the mixed terms of γ and J γ .

Change of variable

First, to reduce the form of the expression in (4.19), the authors in [START_REF] Bartoli | Shape-from-template[END_REF] introduce ρ = ||η|| 2 . Note that ρ ∈ C 1 (L, R + * ) since the third component of η is the focal length f = 0. In addition the derivative of ρ is: 4.20) This implies that:

J ρ = 1 ρ η J η (
ρJ ρ = η J η (4.21)
Substituting ρ in (4.19) gives:

ρ 2 J γ J γ + γ 2 J η J η + ρ γ J γ J ρ + J ρ J γ = f 2 J ∆ J ∆ (4.22)
Proceeding into the change of the variable, the unknown γ is substituted by λ ∈ C 1 (L, R)

defined by:

λ = ργ f (4.23)
This implies that:

γ = f λ ρ and J γ = f ρ 2 (ρJ λ -λJ ρ ) (4.24)
Substituting γ and J γ in terms of λ, ρ, J λ and J ρ in (4.22) gives:

f 2 ρ 2 ρ 2 J λ J λ + λ 2 J ρ J ρ -λρ(J λ J ρ + J ρ J λ ) + f 2 λ 2 ρ 2 J η J η+ f 2 λ ρ 2 ρJ λ J ρ -λJ ρ J ρ + ρJ ρ J λ -λJ ρ J ρ = f 2 J ∆ J ∆ (4.25)
With simple symbolic computation, equation (4.25) is simplified into:

J λ J λ + 1 ρ 2 J η J η -J ρ J ρ λ 2 = J ∆ J ∆ (4.26)
Which is equivalent to the equation:

J λ J λ + ξλ 2 = J ∆ J ∆ where ξ = 1 ρ 2 J η J η -J ρ J ρ (4.27)
The form of the equation in (4.27) is much simpler than the equation of (4.19). In addition, it has λ as the main unknown. But λ = ργ f , so finding λ is equivalent to find the depth function γ as ρ is known. In the next paragraph, we give the solution of the equation in (4.27).

Analytical solution

The equation in (4.27) has a unique solution, and it can be obtained analytically. To this end first, J λ is substituted by an independent vector function β . So the equation (4.27) becomes:

β β + ξλ 2 = J ∆ J ∆ (4.28)
Rearrangement of the equation (4.28) gives:

β βξ -1 = J ∆ J ∆ ξ -1 -λ 2 I 2 (4.29)
According to [START_REF] Bartoli | Shape-from-template[END_REF], equation ( 4.29) has a unique solution for λ given by:

λ = λ 2 J ∆ J ∆ ξ -1 (4.30)
where λ 2 (.) refers to the second eigenvalue of the associated matrix. Substituting λ in equation (4.28) gives two possible solutions for β:

β = ± λ 1 (Λ)V 1 (Λ) (4.31)
where λ 1 (.) and v 1 (.) refer to the first eigenvalue and eigenvector of the associated matrix and

Λ = J ∆ J ∆ -λ 2 J ∆ J ∆ ξ -1 ξ (4.32)

Numerical solution

In [START_REF] Bartoli | Shape-from-template[END_REF], a numerical algorithm is proposed to find directly the 3D reconstructed surface using the analytical solution obtained by the SfT method. Practically, a set of points in the 2D flattened template L is selected. Then the depth of the corresponding 3D points in the deformed surface S are computed. This is done point by point independently using the formula:

γ(p) = λ 2   J ∆ (p)J ∆ (p) J η (p)J η(p) - 1 ||η(p)|| 2 2 J η (p)η(p)η (p)J η(p) -1   (4.33)
where p is a selected point in L. So as it appears, the numerical solution depends mainly on the warp function η presented in η, its Jacobian matrix, and the Jacobian matrix of the embedding function ∆. Accordingly, a Matlab algorithm called SfT_BGCC12I has been written to find the 3D reconstructed point M ∈ S corresponding to the 2D point p ∈ L. This algorithm executes in the following manner:

The inputs:

• The coordinates of η(p), which are the coordinates of some point m in the image J .

• The 2 × 2 Jacobian matrix J η (p).

• The 2 × 2 matrix J ∆ (p)J ∆ (p). This matrix is taken by default as I 2 , the identity matrix of dimension 2.

• The 3 × 3 matrix K -1 ), which is the inverse of the characteristic matrix of the camera K. This input is needed if the considered image is taken in pixels and not in mm.

The output:

• The 3D reconstructed point M ∈ S, which is in fact the image of p by ϕ.

• The two solutions of the 3 × 2 Jacobian matrix J ϕ (p) denoted by J 1 and J 2 . Two solutions exist for J ϕ since β in (4.29) has two possible solutions as in (4.31). J 1 and J 2 represent the two possible tangents to the 3D deformed surface S at M .

• The two solutions of the normal vector to the deformed surface S at the reconstructed point M . These two vectors, denoted by N 1 and N 2 , are obtained by the cross-product of the two columns of J 1 and J 2 respectively, followed by normalization.

The Matlab algorithm SfT_BGCC12I, presented schematically in Fig. 4.7, is the basic reference that we rely on to study the sensitivity of the SfT method. In the next sections we indicate the inputs that are considered in our sensitivity study and then we describe the model that we use to study its sensitivity. This model should take the indicated inputs and it should involve SfT_BGCC12I which represents the concept of the SfT method.

The model inputs

It was detected experimentally that the quality of the 3D reconstruction by the SfT method depends on the depth Z of the deformed surface in front of the camera. In addition, according to the [Chhatkuli et al., 2017], two other parameters may also affect the results of the SfT reconstruction method: the focal length of the camera f and the orientation of the 3D deformed surface in front of the camera. The effect of these two parameters appears either as an independent defect or as an interaction with Z. For that, we select Z, f , and the orientation of the 3D deformed surface to be the inputs of our considered model.

In the next section we present how we construct this model, which is used to study the sensitivity of the SfT method with respect to these indicated inputs.

The model ZfT_SfT

The main inputs of the SfT_BGCC12I algorithm are the 2D point η(p) and J η (p), so at least the depth does not appear as an input, in fact as output. Since our interest is the sensitivity of the SfT method with respect to the depth, f , and the orientation, another model that takes these inputs should be constructed, while involving the concept of the SfT method. For that, we build a new model ZfT_SfT whose inputs are: the depth Z, the focal length f , and the orientation of the 3D surface.

The model ZfT_SfT is described in the following steps:

1. We fix as an initialization two coordinates (X 0 , Y 0 ) of a 3D point M in the deformed surface. Then the given input depth denoted Z G is the third coordinate of M , So M = (X 0 , Y 0 , Z G ). This point is considered as a reference and it is the point to be recovered by the SfT reconstruction method.

2. The given focal length f is used to construct the intrinsic matrix K of the camera.

3. The orientation of the 3D deformed surface is given as an angle θ which is used to indicate the direction of the normal vector to the deformed surface at M . Indeed, the vector M is rotated around the Z direction by the angle θ. The new obtained vector is normalized and then taken as the normal the normal vector N at a point M ∈ S. This is illustrated in Fig. 4.8. J Π (M ) can be simply computed since the projection function of the camera Π is known. So we still have to find J ϕ (p), which can be determined using the normal vector N . Indeed, the Jacobian matrix J ϕ (p) is a 3 × 2 matrix, whose two columns, denoted by U 1 and U 2 , are the generating vectors of the tangent plane at ϕ(p) = M in S. But N is orthogonal to this tangent plane, so N is orthogonal to U 1 and U 2 .

Thus, by rotating N by 90 • we obtain U 1 . Then U 2 is obtained as the cross product of U 1 and N , i.e. U 2 = U 1 ∧ N . This implies that J ϕ (p) is determined using N , and consequently J η (p) is found.

6. η(p), J η (p) and K are now ready, then the algorithm SfT_BGCC12I is run to reconstruct the point M . Note that in this model we take J ∆ (p)J ∆ (p) = I 2 , the identity matrix of dimension 2, as the default case. Let M denote the new reconstructed 3D point, and let N 1 and N 2 denote the two obtained normal vectors from the SfT_BGCC12I algorithm.

7. In this last step the error which forms the output of the model ZfT_SfT is computed. Indeed, two errors are taken into consideration:

• The 3D reconstruction error defined by:

E r = ||M -M || 2 (4.36)
• The normal error defined by:

E N = min{arccos( N • N 1 ), arccos( N • N 2 )} (4.37)
The model ZfT_SfT defined by the above steps is summarized in the flow chart of Fig. 4.9.

The model ZfT_SfT presented matches our needs to study the sensitivity of the SfT Figure 4.9: A scheme of the model ZfT_SfT.

method as it has a scalar output and it takes the indicated inputs. However, while doing simulations, the obtained errors are very small, and almost negligible. This is because all the steps and the procedures are a direct implementation of the theoretical concepts, and this does not reflect reality. Indeed, images taken in reality are exposed to noise whatever the focal length and the resolution are. For that, we propose a method to add a reasonable amount of noise on η(p) and J η (p) simultaneously while running the ZfT_SfT model.

The method of adding noise is explained in the next section, and the modified model of ZfT_SfT that is associated with noise will be denoted by ZfT_SfT'. This model takes the same inputs as ZfT_SfT and gives the same outputs, so it is the one used in our sensitivity study.

Imitating reality by adding noise

To imitate reality while studying the sensitivity of the SfT method by the ZfT_SfT model, we add Gaussian noise on the computed η(p) and J η (p) in the body of the model before running the SfT_BGCC12I algorithm. This Gaussian noise should be added to η(p) and J η (p) simultaneously and in a compatible manner as it happens in reality. For that, we develop a method that derives a correspondence between the standard deviation of the Gaussian noise σ and the error in η(p) and J η (p). This correspondence appears as a plot of the variation of the error in η(p) and J η (p) as a function of σ. This helps in indicating the appropriate value of the σ according to a required amount of error in η(p) and J η (p). The following is a detailed description of the method used to derive the correspondence:

1. We create a 2D regular grid representing the pixel coordinates (U 5. To detect how much error is obtained on the warp function and its derivative due to the Gaussian noise associated with the selected σ, we derive an affine approximation H N of the warp function between I N ref and I N def . This derivation is done using the method DIRT explained in [Bartoli, 2008], which needs an initial guess of

H N defined locally in I ref .
For that, we select a point q ∈ I ref , and we define a rectangular patch D ⊂ I ref around q. Then we derive a local affine transformation G around q that matches H on D. To find this G we use the Taylor Expansion of H around q which is defined by:

H(q + h) = H(q) + J H (q)h (4.38)
So we define G for any q ∈ D by substituting h in (4.38) by q -q, this implies that:

G(q ) = J H (q)q -J H (q)q + H(q) (4.39)

Note that G is well defined according to (4.39), since J H (q) is a 2 × 2 matrix and q and q are in D ⊂ I ref ⊂ R 2 . Once G is found, H N is computed using the DIRT method and the error of the warp and its derivative can be derived. Indeed, H N is an affine map between I N ref and I N def , approximating the warp function. Thus its formula consists of the value of the warp and its first derivative.

6. To derive the error of the warp function and its derivative we consider the point q N ∈ I N ref corresponding to the point q ∈ I ref .

Then we select two different points in the neighborhood of q N and denote them by q N and q N corresponding to the points q and q in I ref .

Then the error is taken as

RM S = 1 3 H N (q N ) -H(q) 2 2 + H N (q N ) -H(q ) 2 2 + H N (q N ) -H(q ) 2 2 (4.40)
Accordingly, the error associated to the selected σ is derived. Note that in the above procedure only one point q is selected to derive the local approximation, but in practice we choose a set of points of interest, and then for each point of interest we derive the associated RM S as above. Then totally we take the average of all the RM S obtained to get a reasonable value of the error of the warp function and its derivative.

With the above steps we derive a correspondence between the σ of the Gaussian noise added to the gray levels of the 2D image and the error obtained on η and J η . In our case study, the selected value of σ is 1.1155, since according to experts this value gives a reasonable error on η and J η . Accordingly, we update the model ZfT_SfT, to obtain a new model ZfT_SfT' whose sensitivity is to be studied. Fig. 4.10 represents a scheme of the model ZfT_SfT'.

In the following section we present the framework we follow when applying the Sobol sensitivity method.

Sensitivity analysis of SfT

Recall that our aim is to see how information about the inputs can be extracted by studying the variation of Sobol indexes of these inputs, with respect to some factors. For that, the sensitivity of the model ZfT_SfT' is not studied by just applying basic Sobol method. Indeed, we extend the application to study the evolution of the Sobol indexes of Z, f , and θ as a function of the depth Z. The reason behind choosing the variation factor to be the depth Z is technical. Actually, a basic role of the depth has been detected experimentally when applying the SfT method, in which three different results obtained when varying Z. So studying the variation of the Sobol indexes of Z, f , and θ as a function of the depth Z would give information about the three depth positions that are affecting the 3D reconstruction by the SfT method. In this way, we keep our study realistic and the conclusions are derived with scientific significance.

To this end, we first compute the Sobol indexes of Z, f , and θ corresponding to the output E r . Then, in a next step, we compute the Sobol indexes of Z, f , and θ corresponding to the output E N . Each of the inputs Z, f , and θ is considered to vary as a uniformly distributed random variable with a specified interval:

• For θ we consider three different cases with three different intervals: [15 -The depth Z as uniformly distributed over the interval [Z 0 -5, Z 0 + 5] with Z 0 varying in the interval [44mm, 2040mm].

-The focal length f as uniformly distributed over the interval [f 0 -5, f 0 + 5]

where f 0 = Z 0 -20mm.

The computation of the Sobol indexes is done numerically using the formulas (1.54),

(1.55),and (1.62) of chapter 1. We use the Monte Carlo simulation method with 4000 samples to compute each Sobol index.

We start by the computation of the first order Sobol indexes. Algorithm 1 represents the steps we follow in this computation.

Algorithm 1: First order Sobol indexes of Z, f , and θ as a function of Z with [15 • , 30 • ] as a sampling space for θ For Z 0 = 40mm upto 2040mm

f 0 = Z 0 -20mm sample spaces: Z: [Z 0 -5, Z 0 + 5] f : [f 0 -5, f 0 + 5] θ: [15 • , 30 • ] sample1 = {(Z k , f k , θ k )} k=1,••• ,4000 sample2 = {( Zk , fk , θk )} k=1,••• ,4000 Ēr = 1 4000 4 k=1 000E r (Z k , f k , θ k ) V ar(E r ) = 1 4000 4 k=1 000E 2 r (Z k , f k , θ k ) -Ē2 r V ar(E rZ ) = 1 4000 4000 k=1 E r (Z k , f k , θ k )E r (Z k , fk , θk ) -Ē2 r V ar(E rf ) = 1 4000 4000 k=1 E r (Z k , f k , θ k )E r ( Zk , f k , θk ) -Ē2 r V ar(E rθ ) = 1 4000 4000 k=1 E r (Z k , f k , θ k )E r ( Zk , fk , θ k ) -Ē2 r S Z = V ar(E rZ ) V ar(Er) (first order Sobol index of Z) S f = V ar(E rf ) V ar(Er) (first order Sobol index of f ) S N = V ar(E rθ ) V ar(Er) (first order Sobol index of θ)
Similarly, the first order Sobol indexes are computed for the cases of θ sampled on the intervals [30 For θ ∈ [15 • , 30 • ], the reconstruction error E r (Fig. 4.11(a)) is highly sensitive to Z for depth greater than 500 mm, but has very little sensitivity to θ. In addition, for all the values of the depth, the influence of the focal length f is constant with a small value of around 0.2. To conclude, it is possible to extract two separate phases:

• Depth less than 500 mm: in this range E r is mostly affected by the variation of Z and θ, but the influence of Z increases while the influence of N decreases.

• Depth greater than 500 mm: in this range E r is highly affected by the variation of Z while the effects of θ and f are almost negligible.

For θ ∈ [30 • , 45 • ], the reconstruction error E r (Fig. 4.11(b)) is also highly sensitive to Z especially for depth greater than 500 mm. In addition, the sensitivity of E r with respect to the focal length f is very small (about 0.2%), however the sensitivity of E r with respect to θ varies over the range of the depth. Indeed, S N increases slightly for depth less than 200 mm, after that it decreases gradually to reach S f at 700 mm, and then almost vanishes at depth 2000 mm. This enables us to extract three different phases:

• Depth less than 200 mm: in this range E r is mostly affected by the variation of Z, but the influence of θ increases. The influence of f is negligible.

• Depth in [200 mm, 700 mm]: in this range E r is affected by the variation of Z and N , but the influence of N decreases. The influence of f is negligible.

• Depth greater than 700 mm: in this range E r is highly affected by the variation of Z. The effects of N and f are negligible.

For θ ∈ [45 • , 60 • ], the variation of S Z ,S f and S N (Fig. 4.11(c)) is analogous to the previous case. Accordingly, three different phases can be extracted:

• Depth less than 250 mm: in this range E r is mostly affected by the variation of Z, but the influence of θ increases. The influence of f is negligible.

• Depth in [250 mm, 1000 mm]: in this range E r is affected by the variation of Z and N , but the influence of N decreases. The influence of f is negligible.

• Depth greater than 1000 mm: in this range E r is highly affected by the variation of Z. The effects of N and f are negligible.

Concerning E N , for the different cases of θ (Fig. 4.11(d), (e), and (f)), the evolution of the Sobol indexes is approximately the same. All over the range of the depth, θ has the most influence on the E N . The effects of S Z and S f are very small (around 0.15 mm), except for depth less than 250 mm where S Z has a slightly more effect than S f .

These results obtained from the first order Sobol indexes are compromising, however to check if there is an interaction between Z, f , and θ affecting E r and E N , we compute the Total Sobol indexes. We apply the same numerical procedures as for the First order indexes, with 4000 samples in a Monte Carlo simulation. We consider three different sampling intervals for θ: [15 given in [START_REF] Saltelli | Variance based sensitivity analysis of model output. design and estimator for the total sensitivity index[END_REF] for the computation of the total Sobol indexes. Algorithm 2 gives the details of the computation in case of sampling interval [15 • , 30 • ] for θ.

Similarly, the total Sobol indexes are computed for the other sampling intervals of θ and for the output E N . Fig. 4.12 presents the results obtained.

As it can be seen, the total Sobol indexes are approximately equal to the first order Sobol indexes for both outputs E R and E r , this implies that the second order Sobol indexes are almost zero. Thus we conclude that there is no interaction between Z, f , and θ affecting

Z 0 = 40mm upto 2040mm f 0 = Z 0 -20mm sample spaces: Z: [Z 0 -5, Z 0 + 5] f : [f 0 -5, f 0 + 5] θ: [15 • , 30 • ] sample1 = {(Z k , f k , θ k )} k=1,••• ,4000 sample2 = {( Zk , fk , θk )} k=1,••• ,4000 Ēr = 1 4000 4 k=1 000E r (Z k , f k , θ k ) V ar(E r ) = 1 4000 4 k=1 000E 2 r (Z k , f k , θ k ) -Ē2 r V ar(E r∼Z ) = 1 4000 4000 k=1 E r (Z k , f k , θ k )E r ( Zk , f k , θ k ) -Ē2 r V ar(E rf ) = 1 4000 4000 k=1 E r (Z k , f k , θ k )E r (Z k , fk , θ k ) -Ē2 r V ar(E rθ ) = 1 4000 4000 k=1 E r (Z k , f k , θ k )E r (Z k , f k , θk ) -Ē2 r T S Z = V ar(Er)-V ar(E rZ ) V ar(Er) (total order Sobol index of Z) T S f = V ar(Er)-V ar(E rf ) V ar(Er) (total order Sobol index of f ) T S N = V ar(Er)-V ar(E rθ )
V ar(Er)

(total order Sobol index of θ)

Charge transport model as a black box

The charge transport model considered in our work is described in [Le Roy et al., 2003].

Here the model is represented as a black box, having four inputs and two outputs (see The scientific idea behind studying such variation is indicated as follows. Indeed, most of the parameters (inputs) of a charge transport model cannot be determined by independent experiments and it is a heavy task to estimate their values that best fit the experimental data. In this manner, optimization algorithms play an important role in systematizing this part of the modeling activity. However, to facilitate the convergence of the optimization algorithms, it is important to quantify the effect of each input on the output in order to limit the optimization to the most influential inputs. Plus, it is important to know the best experimental configuration at which the data can be collected to be used in an optimization process. In the next, section we present in detail the sensitivity study carried out and the obtained results, and then discuss their significance.

Sensitivity analysis of charge transport model

To carry out a sensitivity study we consider each of the inputs as a uniformly distributed random variable on a given range. Then, accordingly, the first order Sobol indexes are computed according to the numerical method defined in chapter 1. The ranges of the inputs are defined by their lower and upper bound which are indicated in According to experts, these ranges of variation are chosen for several reasons. Firstly, to be certain to keep the physical sense to the conditions. Secondly, to have a large range of inputs in order to assume a broad and consistent representation of the output data.

Lastly to have tractable computation. The computation of the Sobol indexes is done under the variation of the three experimental factors: the temperature T , strength of the applied electric field E, and the time of application t. The way this variation is applied is described in the next paragraph.

The variation protocols

Three different protocols are applied to the dielectric material used in order to drive the data used in the computation of the Sobol indexes:

_ First protocol (red curve): Sobol indexes are estimated using experimental data obtained by using a DC electric field of 30kV.mm -1 applied for charging and discharging times of 20 min. The sensitivity analysis is carried out considering that experimental data can be acquired over a temperature range of [0, 90

• C].
_ Second protocol (blue curve): Sobol indexes are estimated using the same material under a temperature of 40 • C and for charging and discharging times of 20 min. The sensitivity analysis is performed considering an applied electric field varying over the range [10, 80kV.mm -1 ].

_ Third protocol (green curve): Sobol indexes are estimated under a temperature of 40 • C and an applied electric field of 30kV.mm -1 . The sensitivity analysis is performed considering charging and discharging times varying over the range [1, 60min].

According to these protocols, the variation of the first order Sobol indexes is studied.

In the analysis of the results, we consider the indexes relatively in % to facilitate the comparison between the different protocols. Also, we consider the influence of a given input on charge or current density as negligible if its associated index does not exceed 20% (hatched area on the figures). Indeed, a Sobol index below 20% shows that the chosen experiment protocol does not give sufficient information to estimate the selected input with an optimization algorithm. does not influence the current density much. For this output, Sobol indexes are below 10% whatever the protocols used. On the other hand, Sobol indexes exceed 50%, meaning a great influence, on the charge density at low temperature (below 30 • C) or in charging at short time (less than 10 min). For the both cases, it means that the impact on barrier for injection on the deposited charge is important. The fact that this parameter is influential at the beginning of polarization is in phase with the experimental observation. Indeed, when an electric field is applied, charges are injected in the vicinity of the electrodes.

The presence of these charges close to the electrode induces a decrease of the electric field at the interface over time and so a decrease of the injection flux. So, the influence over longer times is less important. Roughly, it corresponds to space charge limited process, which also explains why the barrier to injection is not strongly influential on the external current, which corresponds to the space-averaged trapped current [START_REF] Baudoin | Bipolar charge transport model with trapping and recombination: an analysis of the current versus applied electric field characteristic in steady state conditions[END_REF].

Finally, Fig. 4.14 also shows that the barrier height to injection does not impact the charge density at T = 40 • C, irrespective of the field. Sobol indexes are always below 10%. According to the results, it seems difficult to find a suitable experimental protocol for the optimization purpose (most of the results are in the hatched area). The temperature seems to be the most impacting protocol factor for the mobility for both outputs: charge and current density. A temperature higher than 70 • C allows us to achieve a Sobol index higher than 20%. Conduction takes place via a constant effective mobility µ, leading to the transport of carriers through shallow levels that are related to the structural disorder of the dielectric.

Sensitivity analysis of the mobility

Deep trapping, mainly due to chemical disorder in the material, is described using a unique level of deep traps for each kind of carrier. Charges have a certain probability to escape from deep traps by overcoming a de-trapping barrier w tr . Based on this physical description, results show that this model gives more importance to the charges in shallow traps than in deep traps at high temperature. Indeed, for a temperature higher than 70 • C, the effective mobility has a greater impact on the outputs (Sobol indexes exceed 20%) which helps to estimate it well by the optimization process. Perhaps this increase is due to the fact that at high temperature the fraction of charges in shallow traps is higher.

Sensitivity analysis of the de-trapping barrier height

The results related to the deep trap depth, or de-trapping barrier height w tr are summarized in Fig. given applied electric field of 30kV.mm -1 . However, this input does not affect the current density so much. Sobol indexes are below 10% whatever the protocols used except at high temperature. In general, long charging times are preferred for improving sensibility to the de-trapping coefficient.

Sensitivity analysis of the trapping coefficient

Fig. 4.17 concerns the influence of trapping coefficient B on current and charge density.

Very clearly here, this input has little effect on the charge or current density. Whatever the protocols used, Sobol indexes are always below 10%. It is not very easy to explain this feature because trapping and de-trapping phenomena are obviously linked by nature.

According to experts, an explanation could be provided for the obtained difference between the effect of the trapping an de-trapping coefficient. One reason could be the chosen ranges used for each input. Even though these ranges are chosen in a way to keep the physical sense, this may affect the response of each of the input differently. Another option for this difference could be the fixed trap density used in this model. Indeed, this may represent a very low density of defects (3.2 × 10 14 /cm 2 ), and this may limit the role of the trapped charges in the net charge distribution and in the current density.

Parameter optimization

By analyzing the results obtained from the variation of the Sobol indexes a strategy of study can be designed for parameter optimization for the charge transport model.

Indeed, optimization algorithms are used to find a set of parameters able to minimize the deviations between experimental data and simulation data, as shown in Fig. 4.18.

Experimental data are the net density of charge as measured by the pulsed electro-acoustic method -PEA -and external charging and discharging current measurements [Liu et al., 1993]. Simulated data are those produced by the designed charge transport model. Based on the parameter sensitivity analysis it is possible to find suitable experimental conditions to obtain optimized estimation of the model parameters used in our charge transport model. The three experimental protocols used with the chosen dielectric material, in film form of thickness D = 200µm, give the following guidelines to provide a good approach to the model parameters: Then, the obtained experimental results could be inserted into an optimization algorithm in order to find the new set of parameters. Unfortunately, no straightforward optimal conditions appear for identifying the trapping coefficient. According to experts, analysis is in progress to understand why Sobol indexes are so low in this case. Recombination processes, and electroluminescence as its pending experimental information, are not incorporated in the model. This could be a route to resolve the issue.

Conclusion

In this chapter we study how the variation of the Sobol indexes helps in deriving conclusions concerning the inputs. In this manner we consider two applications, the first is for a model from the domain of computer vision and the second is for a charge transport model.

In the first application, we were able to get information concerning the inputs while studying the variation of Sobol indexes. In detail we consider a 3D reconstruction method which seems to be sensitive to the depth of the 3D surface in front of the camera. We adopt a model involving this method, which gives a scalar output and imitates reality by adding noise to images. Then we examine the sensitivity of this model with respect to three inputs: the focal length of the camera, the depth of the 3D deformed surface, and the orientation of the surface. Sobol indexes are computed numerically and they vary as a function of the depth. Plotting the obtained results allows us to visualize clearly how the effect of the depth can be partitioned into different domains.

Considering the charge transport model, it consists of four inputs and two outputs. We study the Sobol indexes of the four inputs under three different protocols. Each protocol accounts for the variation of one of the experimental factor: the temperature, strength of the applied electric field, and the time of application. Analyzing the results obtained allows us to conclude the most suitable experimental configuration for each input to collect the data associated with its estimation process.

Final words, extending the application of the Sobol method by considering the examination of the variation of Sobol indexes has proved its efficiency. Indeed, more than just information and conclusions concerning the inputs can be derived. Appropriate setting for doing an uncertainty study can be also deduced from studying this variation and this would help a lot in the domain of the backward propagation of uncertainty. In addition, the most favorable ranges for approximating the parameters of a model in an optimization process can be also detected by this variation study.

study. We focus in this manner on the Sobol sensitivity method which is considered as one of the strongest sensitivity methods. It performs the sensitivity study by computing sensitivity indexes, called Sobol indexes, in a probabilistic manner.

In chapter 2 we establish a new method for backward uncertainty propagation. Its aim is to quantify the input uncertainty starting from the data of an uncertain output. The method consists of two main steps. The first step is to write the output variance in terms of the moments of the inputs. The second step is to solve a nonlinear least square problem generated by the expression of step 1. By solving this least square problem numerically the input uncertainty is obtained. Applications show that the method gives very near approximations for the values of the input uncertainty. However the accuracy increases as the number of samples used increases. This issue could be one of our perspectives for future work.

In chapter 3 we show that the Sobol method can give reliable results even when applied in the discrete case. In general, the Sobol sensitivity method is applied to determine the effect of each of the inputs on the output represented by the Sobol indexes. Usually inputs in this method are considered to vary as continuous random variables in order to compute the corresponding indexes. In this chapter, we expand this idea and we applied Sobol methods for inputs with a discrete representation. To this end, we considered a model from the domain of force spectroscopy and we studied its sensitivity by deriving the associated Sobol indexes. The model is a representation of an experimental curve called the Electrostatic Force distance curve (EFDC). The data was limited in this case study (only 45 samples) due to the cost of each experimental run. For that, we performed the computation of the Sobol indexes using the formulas of discrete random variables.

The obtained results of the Sobol method were confirmed using Design of Experiment (DOE), where a total agreement was noted. In addition to this, we invested the obtained sensitivity results in constructing an approximating formula that describes the EFDC.

The approximated EFDCs are plotted and compared to the experimental EFDCs. The results obtained demonstrate a precise compatibility between the two curves.

In addition to their contribution in the area of sensitivity analysis, these results allow the EFDC to be a key starting point to detect charge localization in dielectrics. A next step for this work would check if a similar derivation of the Sobol method can be done for discrete random variables. In addition, future work in the force spectroscopy domain would progress in deriving a systematic strategy for detecting the charge localization using the EFDC.

In chapter 4 we continued in developing the framework for the applications of the Sobol method and we showed how different information and conclusions can be derived for the inputs by studying the variation of their corresponding Sobol indexes with respect to some active factors of the model or some experimental settings. In addition, we demonstrated that studying this variation also allows the determination of the most suitable configuration for estimating the inputs. This can help in the quantification of the input uncertainty in a backward propagation process. These deductions are the consequence of two different applications considered in this chapter.

In the first application we were able to get information concerning the inputs while studying the variation of the Sobol indexes. In detail we considered a 3D reconstruction method which seems to be sensitive to the depth of the 3D surface in front of the camera. We adopted a new model involving this method, that at the same time gives a scalar output and imitates reality by adding noise to images. The strategy we derived to add noise in the model depends mainly on the gray level of the images. After that we examined the sensitivity of the model with respect to three inputs: the focal length of the camera, the depth of the 3D deformed surface, and the orientation of the surface. Sobol indexes were computed numerically and they varied as a function of the depth, since the depth is a significant factor in the 3D reconstruction process as it affects the quality of the obtained reconstruction according to experts. Plotting the variation of the Sobol indexes with respect to the depth allowed us to visualize clearly how the effect of the depth can be partitioned into different domains.

In the second application, we considered a charge transport model, it consists of four inputs and two outputs. We studied the Sobol indexes of the four inputs under three different protocols. Each protocol accounts for the variation of one of the experimental factors: the temperature, the strength of the applied electric field, and the time period at which the electric field is applied. From this study, we were able to extract for each input the most suitable experimental configuration in order to collect the data associated with the inputs estimation process. However, only one of the inputs seems not to be sensitive in any of the protocols. The future work will focus on how experimental data can be best derived to be used in an optimization process to approximate this input.

For the variation of Sobol indexes, we concentrated in this work on only first order Sobol indexes. Our next step would be to study the variation of the second order Sobol indexes in case of having interaction between the inputs (non zero second order Sobol indexes).

The aim of studying this variation is to try to derive a conclusion or configuration so that the interaction of the inputs can be minimized.

From here several goals can be set as future work for the short and the long term. As short term goals, we could seek improvements in the performance of the derived backward propagation method, mainly concerning the number of samples needed. We could also consider the case study of backward propagation of uncertainty for models that are not defined explicitly by functions or if the defined function is not continuous. On the other hand, we can study the variation of the second order Sobol indexes searching for an indication about minimizing the interaction between inputs. This helps in decreasing the over all uncertainty of the model by minimizing the interaction between uncertain inputs.

As long term goals, one would consider the case in which the structural uncertainty is taken into account in the backward propagation of uncertainty. This means that the output uncertainty is not only due to input uncertainty, the structural uncertainty plays a role in producing output uncertainty. So how a backward propagation of uncertainty can partition the output uncertainty between the inputs and the structure of the model.

In general, these goals are ideas that can contribute in solving some of the remaining issues in the domain of uncertainty. .4)

m 1 = µ m 2 = µ + V m 3 = µ + 3µV m 4 = µ + 6µ 2 V + 3V 2 m 5 = µ + 10µ 3 V + 15µV 2 m 6 = µ + 15µ 4 V + 45µ 2 V 2 + 15V 3 m 7 = µ + 21µ 5 V + 105µ 3 V 2 + 105µV 3 m 8 = µ + 28µ 6 V + 210µ 4 V 2 + 420µ 2 V 3 + 105V 4 m 9 = µ + 36µ 7 V + 378µ 5 V 2 + 1260µ 3 V 3 + 945µV 4 m 10 = µ + 45µ 8 V + 630µ 6 V 2 + 3150µ 4 V 3 + 4725µ 2 V 4 + 945V 5 (A

  Dans de nombreuses disciplines, les approches permettant d'étudier et de quantifier l'influence de données incertaines sont devenues une nécessité. Bien que la propagation directe d'incertitudes ait été largement étudiée, la propagation inverse d'incertitudes demeure un vaste sujet d'étude, sans méthode standardisée. Dans cette thèse, une nouvelle méthode de propagation inverse d'incertitude est présentée. Le but de cette méthode est de déterminer l'incertitude d'entrée à partir de données de sortie considérées comme incertaines. Parallèlement, les méthodes d'analyse de sensibilité sont également très utilisées pour déterminer l'influence des entrées sur la sortie lors d'un processus de modélisation. Ces approches permettent d'isoler les entrées les plus significatives, c'est a dire les plus influentes, qu'il est nécessaires de tester lors d'une analyse d'incertitudes. Dans ce travail, nous approfondierons tout d'abord la méthode d'analyse de sensibilité de Sobol, qui est l'une des méthodes d'analyse de sensibilité globale les plus efficaces. Cette méthode repose sur le calcul d'indices de sensibilité, appelés indices de Sobol, qui représentent l'effet des données d'entrées (vues comme des variables aléatoires continues) sur la sortie. Nous démontrerons ensuite que la méthode de Sobol donne des résultats fiables même lorsqu'elle est appliquée dans le cas discret. Puis, nous étendrons le cadre d'application de la méthode de Sobol afin de répondre à la problèmatique de propagation inverse d'incertitudes. Enfin, nous proposerons une nouvelle approche de la méthode de Sobol qui permet d'étudier la variation des indices de sensibilité par rapport à certains facteurs du modèle ou à certaines conditions expérimentales. Nous montrerons que les résultats obtenus lors de ces études permettent d'illustrer les différentes caractéristiques des données d'entrée. Pour conclure, nous exposerons comment ces résultats permettent d'indiquer les meilleures conditions expérimentales pour lesquelles l'estimation des paramètres peut être efficacement réalisée. Mots clés: Quantification d'incertitude, Propagation inverse d'incertitude, Analyse de sensibilité, Méthode de Sobol.

Figures 1 . 1

 11 The data fitted to a logarithmic function and to a 4th root function. . . . . 1.2 The measurements of the height x 1 and the radius of the cylinder x 2 . . . . 1.3 Probabilistic presentation of input uncertainty: (a) continuous distribution, (b) discrete distribution . . . . . . . . . . . . . . . . . . . . . . . . . 1.4 Probabilistic propagation of input uncertainty. . . . . . . . . . . . . . . . . 1.5 Scatter plot of sampling techniques: (a) Random sampling, (b) Latin Hyper cube sampling. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.6 Fuzzy membership functions: (a) Triangular, (b) Trapezoidal. . . . . . . 1.7 Bayesian inference for backward uncertainty propagation: prior distribution and output data are used to derive a posterior distribution. . . . . . 1.8 Output variance decomposition for model of two variables representing the concept of the Sobol's method for the derivation of the sensitivity indexes. 2.1 Flowchart that summarizes Step 1. . . . . . . . . . . . . . . . . . . . . . . 2.2 Uncertain real weight (W R ) in terms of the measured weight (W m ). . . . . 2.3 The obtained values of V and r as a function of the number of sample points used. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.4 The obtained values of r 1 and r 2 as a function of the number of sample points used. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.1 Schematic of the AFM parts: 1.Piezoscanner 2.Sample 3.Tip 4.Cantilever 5. Laser emitter 6.Photo-detector. . . . . . . . . . . . . . . . . . . . . . 3.2 The tip-sample positioning during the different phases of AFM. . . . . . . 3.3 A typical AFM Force Distance Curve showing the approach stage, the contact stage, and the retract stage. The labels (a), (b), (c), (d), and (e) refer to the phases of the tip-sample positioning of Fig. 3.2. . . . . . . . . 3.4 The FDCs while scanning an oxy-nitride dielectric (before and after the effect of the electrostatic force). Their difference produces the EFDC. This figure is extracted from [Villeneuve-Faure et al., 2014]. . . . . . . . . . . . 3.5 Sample scheme (electrode is represented in black). . . . . . . . . . . . . . . 3.6 Evolution of EFDC as function of electrodes parameters: (a) Width w (bias v and depth d are fixed at 15v and 100nm respectively) and (b) Depth d (bias v and width w are fixed at 15v and 20µm respectively). . . . 3.7 The illustration of the logistic curve f (z). . . . . . . . . . . . . . . . . . . . 3.8 Residuals of the polynomial and logistic regressions. . . . . . . . . . . . . . 3.9 The model Fitted-4PL . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.10 The first order Sobol indexes of d, v, and w. The indexes S d , S v , and S w represent the influence of d, v, and w, respectively, on A, B, C, and D. . . 3.11 The second order Sobol indexes of d, v, and w. The indexes S dw , S dv , and S vw represent the influence of the interactions between d and w, d and v, and v and w, respectively, on A, B, C, and D. . . . . . . . . . . . . . . . . 79 3.12 The main effect plots of depth d, width w and applied voltage v for each of the 4PL parameters: (a) A, (b) B, (c) C and (d) D. . . . . . . . . . . 82 3.13 The Interaction Effects Matrix plot of the parameter A. . . . . . . . . . . . 84 3.14 one-way interaction plots for the parameters: (a) B, (b) C, and (d) D, in the order : d-w plots, d-v plots and v-w plots. . . . . . . . . . . . . . . .

  , (b) (100 nm, 8 v, 40 µm), (c) (50 nm, 6 v, 6 µm), and (d) (10 nm, 15 v, 40 µm). . . . . . . . . . . . . . . . . . . . . . . . . . . . 90 4.1 The pinhole camera. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 95 4.2 Real vs virtual image by a pinhole camera. . . . . . . . . . . . . . . . . . 96 4.3 Geometric model of the 2D image formation by a pinhole camera. . . . . . 96 4.4 The pixel coordinate system P = {O p ; U p , V p } vs the image plane frame J = {O; U, V}. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 97 4.5 A paper template deformed into a cylinder. . . . . . . . . . . . . . . . . . 100 4.6 The 3D reconstruction model in SfT. . . . . . . . . . . . . . . . . . . . . . 101 4.7 A scheme of the algorithm SfT_BGCC12I. . . . . . . . . . . . . . . . . . 108 4.8 The normal at the point M determined by the angle θ. . . . . . . . . . . . 109 4.9 A scheme of the model ZfT_SfT. . . . . . . . . . . . . . . . . . . . . . . 111 4.10 A schematic of the model ZfT_SfT'. . . . . . . . . . . . . . . . . . . . . . 114 4.11 The variation of the first order Sobol indexes S Z and S f and S N as a function of the depth : (a) For E r with θ ∈ [15 • , 30 • ] , (b) For E r with θ ∈ [30 • , 45 • ], (c) For E r with θ ∈ [45 • , 60 • ], (d) For E N with θ ∈ [15 • , 30 • ], (e) For E N with θ ∈ [30 • , 45 • ], (f) For E N with θ ∈ [45 • , 60 • ] . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 116 4.12 The variation of the total order Sobol indexes T S Z and T S f and T S N as a function of the depth : (a) For E r with θ ∈ [15 • , 30 • ] , (b) For E r with θ ∈ [30 • , 45 • ], (c) For E r with θ ∈ [45 • , 60 • ], (d) For E N with θ ∈ [15 • , 30 • ], (e) For E N with θ ∈ [30 • , 45 • ], (f) For E N with θ ∈ [45 • , 60 • ] . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 120 4.13 The charge transport model. . . . . . . . . . . . . . . . . . . . . . . . . . . 121 4.14 Evolution of the first order Sobol indexes of the injection height barrier w. 124 4.15 Evolution of the first order Sobol indexes of the mobility µ. . . . . . . . . . 125 4.16 Evolution of the first order Sobol indexes of the de-trapping barrier height w tr . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 127 4.17 Evolution of the first order Sobol indexes of the trapping coefficient B. . . 128 4.18 Principle of the optimization technique. . . . . . . . . . . . . . . . . . . . . 129

  the part of the model that varies in D at each model run to generate a new output. The parameter α is a vector in R m (m ≥ 0). α is the part of the model that defines its characteristics. It does not change with each model run, however it is given an initialization value once at the beginning of the model use. Both x and α are called the input of the model. The F in the above notation represents the model's structure, it could be mathematical equation(s), computer code(s), or visual representation(s). The output y is the response of the input by F . It is considered here as a scalar, y ∈ R, since same results hold for a multi-scalar output, by considering each component alone. In a modeling process, the parts of the model that can be uncertain are the input x, α and/or the model structure F .

Figure 1

 1 Figure 1.3: Probabilistic presentation of input uncertainty: (a) continuous distribution, (b) discrete distribution .

  Figure 1.4: Probabilistic propagation of input uncertainty.

  of the rows and one of the columns. This ensures more coverage of the range of the inputs than the case of just random sampling. Fig.1.5 shows the difference between the random sampling (basic Monte Carlo) and the Latin Hyper cube sampling for eight samples for two inputs.

Figure 1 . 5 :

 15 Figure 1.5: Scatter plot of sampling techniques: (a) Random sampling, (b) Latin Hyper cube sampling.

Figure 1

 1 Figure 1.6: Fuzzy membership functions: (a) Triangular, (b) Trapezoidal.

  . Due to confusion, it should be noted that backward propagation differs from the well known topics Parameter Calibration and Inverse Problem. Indeed, Parameter Calibration and Inverse Problem are concerned in estimating the parameters of a model, and not the uncertainty, starting from output data. Methods for backward propagation are the least developed among the methods that are derived to study uncertainty [Chen et al., 2015].

Figure 1 . 7 :

 17 Figure 1.7: Bayesian inference for backward uncertainty propagation: prior distribution and output data are used to derive a posterior distribution.

Figure 1

 1 Figure 1.8: Output variance decomposition for model of two variables representing the concept of the Sobol's method for the derivation of the sensitivity indexes.

  Let a = (a 1 , • • • , a n ) be an input point in D the domain of F . Given the data of several samples of the output, coming from some experiments or measurements, in which unequal outcome values correspond to the input point a. This variety in the outcomes indicates that the output is uncertain, i.e. we are not sure what the true value of F (a) is. Since the model is assumed to be deterministic, this implies that the only source of this output uncertainty is the input. It is to say that the values given to F to produce F (a) are not surely equal to a. Therefore, there is an input uncertainty at a causing an output uncertainty at F (a). The aim of our new method is to find this input uncertainty using the given output data.In this work, the probabilistic point of view is used to represent the uncertainty, as described in section 1.4 of chapter 1. The input uncertainty at a is represented by the random vector X = (X 1 , • • • , X n ) and the output uncertainty at F (a) is represented by the random variable Y . The expectation of X, denoted by µ X = (µ 1 , • • • , µ n ), is equal to the best estimate of the true input value. Since the runs of the model are supposed to be performed at the same input point (a) then we take µ X = a. On the other hand, the expectation of Y , denoted by µ Y , which represents the best estimate of F (a) is computed from the given data. The quantity of the output uncertainty at F (a) is represented here by the variance of Y which is denoted by V ar(Y ) and it is computed from the given output data.

  clarify the idea of taking the input uncertainty as a function of the input point consider the following simple example: Consider a balance that measures the human's weight, labeled by ±0.01 W m , where W m here refers to the measured weight. This label means that this balance has an uncertainty defined by the interval [-0.01, +0.01] × W m . Clearly, one can see that this uncertainty is a function of the measured weight W m . So for a man weighing 70 kg on this balance, his real weight is a value in the interval [70 -0.7, 70 + 0.7] and not exactly 70. Thus the uncertainty here is 0.7kg. However, for a man weighing 100 kg, the real weight is in the interval [100-1, 100+1], and so the uncertainty here is 1kg. This evolution of uncertainty is illustrated graphically in Fig.2.2.

Figure 2 . 1 :

 21 Figure 2.1: Flowchart that summarizes Step 1.

Figure 2 . 2 :

 22 Figure 2.2: Uncertain real weight (W R ) in terms of the measured weight (W m ).

  R . If for any x ∈ D, ∇ x ≤ ε , then we are done and the quantity of the input uncertainty is the obtained values of the elements of V and R.

Figure 2

 2 Figure 2.3: The obtained values of V and r as a function of the number of sample points used.

Figure 2 . 4 :

 24 Figure 2.4: The obtained values of r 1 and r 2 as a function of the number of sample points used.

Figure 3

 3 Figure 3.1: Schematic of the AFM parts: 1.Piezoscanner 2.Sample 3.Tip 4.Cantilever 5. Laser emitter 6.Photo-detector.

Figure 3 . 2 :

 32 Figure 3.2: The tip-sample positioning during the different phases of AFM.

  Fig. 3.3 shows a typical FDC illustrated with the associated phases of the tip-sample positioning.

Figure 3

 3 Figure 3.3: A typical AFM Force Distance Curve showing the approach stage, the contact stage, and the retract stage. The labels (a), (b), (c), (d), and (e) refer to the phases of the tip-sample positioning of Fig. 3.2.

Figure 3 . 4 :

 34 Figure 3.4: The FDCs while scanning an oxy-nitride dielectric (before and after the effect of the electrostatic force). Their difference produces the EFDC. This figure is extracted from [Villeneuve-Faure et al., 2014].

  As stated in the previous section, the considered EFDCs are obtained by AFM scanning a dielectric with buried electrodes. In detail, the sample structure of the considered dielectric consists of aluminum electrodes embedded in a SiN x layer as shown in Fig.3.5. Since the considered inputs are w, v and d, different samples are designed by varying w, v and d, then their EFDCs are recorded.

Figure 3 . 5 :

 35 Figure 3.5: Sample scheme (electrode is represented in black).

Fig. 3

 3 Fig. 3.6(a) displays EFDCs for different values of the electrode's width w. However Fig. 3.6(b) displays EFDCs for different values of the electrode's depth d.

Figure 3

 3 Figure 3.6: Evolution of EFDC as function of electrodes parameters: (a) Width w (bias v and depth d are fixed at 15v and 100nm respectively) and (b) Depth d (bias v and width w are fixed at 15v and 20µm respectively).

  .1) where the four parameters A, B, C, and D characterize the curve's sigmoid shape. As shown in Fig. 3.7, parameters A and D are the curve's minimum and maximum value respectively, either taken at some values of z or attained asymptotically by the curve.The parameter B is called the Hill slope, which is responsible for the steepness of the bending of the curve. C is called the mean response point, it is the z value at which the curve is at its midway between the min and max. Note that if B ≤ 1 the first bending of the curve almost disappears. In fact, this is the case of the EFDC as shown in Fig.3.6.

Figure 3 . 7 :

 37 Figure 3.7: The illustration of the logistic curve f (z).

Figure 3

 3 Figure 3.8: Residuals of the polynomial and logistic regressions.

  (d, v, w), there will be corresponding values for the logistic parameters A, B, C, and D.Hence, studying the effects of the d, v, w on the EFDCs can be exchanged by studying the effect of (d, v, w) on each of logistic parameters A, B, C, and D, which are scalars.Starting from this point, the mission is to study the sensitivity of the logistic parameters (A, B, C, D) with respect to(d, v, w). For that we consider the model whose inputs are d, v, and w, and outputs are A, B, C, and D. This model is called a Fitted-4PL and it is schematized in Fig.3.9.

  first order Sobol indexes are computed. The results obtained are represented as a bar diagram in Fig. 3.10. For each output A, B, C, and D there are three bars representing the values of the first order Sobol indexes of the three inputs d, v, and w. As depicted in Fig. 3.10, A (i.e. maximum electrostatic force) is mainly influenced by applied bias v and electrode depth d. Contrary to A, D (i.e. electrostatic long range force)is influenced mainly by the electrode width w and equally in the same amount by the applied bias v and the electrode depth d. Moreover, B and C (i.e. curve bending and shape) are more influenced by the applied bias v, however as less influence is by the electrode depth d on B and electrode width w on C.

Figure 3 . 10 :

 310 Figure 3.10: The first order Sobol indexes of d, v, and w. The indexes S d , S v , and S w represent the influence of d, v, and w, respectively, on A, B, C, and D.

Figure 3 .

 3 Figure 3.11: The second order Sobol indexes of d, v, and w. The indexes S dw , S dv , and S vw represent the influence of the interactions between d and w, d and v, and v and w, respectively, on A, B, C, and D.

  . The Factors are the inputs of the experimental model, they are the elements that are varied when the experiment is conducted. A factor may assume at least two distinct values during an experiment, such values are called Levels. So, Levels are the different possible taken values of the factors. The Response is the output of the experimental model; it describes how the system responds under a certain configuration of the input factors. Thus, concerning our study case, the factors are the three inputs d, v, and w. According to the conducted EFDC measurements, the levels of w are 6 µm, 20 µm and 40 µm, the levels of v are 4 V, 6 V, 8 V, V10 V and 15 V, and lastly the levels of d are 10 nm, 50 nm and 100 nm. In addition, the response in this case is the four logistic parameters A, B, C, and D. To proceed in studying the factor-response model with DOE, one should select a DOE technique that holds behind the procedures in which the experimental data is collected, as well as the tests required to analyze this data.

Fig. 3 .

 3 Fig. 3.12 represents the mean of A, B, C, and D as function of the levels of d, v, and w. At each level, the corresponding factor is fixed, and then the mean of the response is computed and plotted as a dot. After performing this at all the levels of the factors, the plotted points of each factor (d, v, and w) are joined by a line. Then accordingly, the dependence of A, B, C, and D (linear, quadratic, • • • ) to each of d, v, and w is derived.

Figure 3

 3 Figure 3.12: The main effect plots of depth d, width w and applied voltage v for each of the 4PL parameters: (a) A, (b) B, (c) C and (d) D.

  Fig. 3.14(a) shows the interaction effects plots of the parameter B. Contrary to A, more crossed lines appear especially in the d-v plot and in the v-w plot. Nevertheless these lines possess a slight non parallelism. For the d-w plot, all the lines are quite similar

Fig. 3 .Figure 3

 33 Fig. 3.14(b) shows the interaction effects plots of the parameter C. Crossed lines appear especially in the d-w plot and in the v-w plot. For the d-v plot, all the lines are quite similar indicating a poor influence of their interaction on B. This is also compatible with obtained second order Sobol indexes associated to C, where S dv ≈ 0.1 and S dw ≈ S vw ≈ 0.2.

  and the experimental variables d, v and w, we proceed in finding an approximation formula of each of A, B, C, and D in terms of d, v and w. This helps also in finding an approximating expression of the electrostatic force in terms of d, v and w by substituting A, B, C, and D in the formula of the logistic law (3.1). The proposed model that we use is known as the matrix Model [Negoescu et Axinte

Figure 3

 3 Figure 3.15: The difference between the EFDCs obtained experimentally (blue) and the EFDCs obtained using the first order matrix approximation models of A, B, C, and D (red). Different values of (d, v, w) are considered: (a) (10 nm, 6 v, 6 µm) , (b) (100 nm, 8 v, 40 µm), (c) (50 nm, 6 v, 6 µm), and (d) (10 nm, 15 v, 40 µm).

Fig. 3 .

 3 Fig. 3.16 represents the EFDCs computed using the second order matrix model for approximating the parameters A, B, C, and D of the logistic law curve, compared to the experimental EFDCs for the same values of (d, v, w) used in Fig. 3.15.

Figure 3

 3 Figure 3.16: The difference between the EFDCs obtained experimentally (blue) and the EFDCs obtained using the second order matrix approximation models of A, B, C, and D (red). Different values of (d, v, w) are considered: (a) (10 nm, 6 v, 6 µm), (b) (100 nm, 8 v, 40 µm), (c) (50 nm, 6 v, 6 µm), and (d) (10 nm, 15 v, 40 µm).

  of force spectroscopy. The inputs of the model are the experimental factors d, v, and w. The output of the model considered should be scalar to allow the application of the Sobol method. For that, the EFDC is fitted by 4PL logistic law, which allows the EFDC to be represented exclusively by the four logistic parameters A, B, C, and D. Then the sensitivity of each of A, B, C, and D with respect to d, v, and w is studied by deriving the corresponding Sobol index. By definition, the inputs d, v, and w can take continuous values, however due to the cost of each experiment run, limited data is given corresponding to the limited discrete values taken by d, v, and w . This gives only 45 different combinations. So we deal with the inputs d, v, and w as discrete variables and we compute the corresponding Sobol indexes using the formulas of discrete random variables. Two kinds of Sobol indexes are considered as the first order Sobol indexes and the second order Sobol indexes. The results are then validated using the factorial plots of Design Of Experiments (DOE). The main effect plots of DOE confirm the results obtained by the first order Sobol indexes. In addition, the interaction effect plots also confirm the results obtained by the second order Sobol indexes.

Fig. 4 .

 4 1 is an illustration of the pinhole camera.

Figure 4 .

 4 Figure 4.1: The pinhole camera.

Figure 4 .

 4 Figure 4.2: Real vs virtual image by a pinhole camera.

I

  = {C; X, Y, Z} centered at the camera center C, and whose positive Z axis is directed towards the scene. This frame induces another orthogonal frame J = {O; U, V} on the image plane as shown in Fig. 4.3. So any point in the image plane with coordinates (a, b) in J, has (a, b, f ) coordinates in the camera frame I.

Figure 4 . 3 :

 43 Figure 4.3: Geometric model of the 2D image formation by a pinhole camera.

  .2) Usually, images we use in daily life have a common measurement unit called 'pixel'. However the unit used in the orthogonal frame J = {O; U, V} is millimeter (mm). For that, a new frame P = {O p ; U p , V p } is attached to the image plane with a new coordinate system called the pixel coordinate system.

Figure 4 . 4 :

 44 Figure 4.4: The pixel coordinate system P = {O p ; U p , V p } vs the image plane frame J = {O; U, V}.

Figure 4 .

 4 Figure 4.6: The 3D reconstruction model in SfT.

Figure 4 .

 4 Figure 4.7: A scheme of the algorithm SfT_BGCC12I.

Figure 4 . 8 :

 48 Figure 4.8: The normal at the point M determined by the angle θ.

Figure 4 .

 4 Figure 4.10: A schematic of the model ZfT_SfT'.

Figure 4 .

 4 Figure 4.11: The variation of the first order Sobol indexes S Z and S f and S N as a function of the depth : (a) For E r with θ ∈ [15 • , 30 • ] , (b) For E r with θ ∈ [30 • , 45 • ], (c) For E r with θ ∈ [45 • , 60 • ], (d) For E N with θ ∈ [15 • , 30 • ], (e) For E N with θ ∈ [30 • , 45 • ], (f) For E N with θ ∈ [45 • , 60 • ] .

E

  r and E N . So the effect of Z, f , and θ on the reconstruction error E r and normal error E N are only individual effects. Accordingly, the three different ranges of the depth can be derived from the sensitivity reconstruction error E r to the individual effects of Z, f and θ. So conclusions and information concerning the input Z are derived by considering the variation of the Sobol indexes of Z, f and θ concerning the SfT method. In the next section, we present another way in which conclusions about the inputs are detected by studying the variation of the Sobol indexes. We consider a new model corresponding to charge transport of dielectric and we study the variation of Sobol indexes of its input as a function of three different experimental factors. Algorithm 2: Total order Sobol indexes of Z, f , and θ as a function of Z with [15 • , 30 • ] as a sampling space for θ For

Figure 4 .

 4 Figure 4.12: The variation of the total order Sobol indexes T S Z and T S f and T S N as a function of the depth : (a) For E r with θ ∈ [15 • , 30 • ] , (b) For E r with θ ∈ [30 • , 45 • ], (c) For E r with θ ∈ [45 • , 60 • ], (d) For E N with θ ∈ [15 • , 30 • ], (e) For E N with θ ∈ [30 • , 45 • ], (f) For E N with θ ∈ [45 • , 60 • ] .

Fig. 4 .

 4 Fig. 4.13). The inputs are the barrier height for injection w, mobility µ, the trapping

Figures 4 .

 4 Fig. 4.14 concerns the influence of barrier height for injection w on the current and charge density. It appears that, for the model and protocol considered, the barrier for injection

Figure 4 .

 4 Figure 4.14: Evolution of the first order Sobol indexes of the injection height barrier w.

Fig. 4 .

 4 Fig. 4.15 concerns the influence of the mobility µ on the current and the charge density.

Figure 4 .

 4 Figure 4.15: Evolution of the first order Sobol indexes of the mobility µ.

  4.16. For temperatures higher than 50 • C, the influence of the charge trapping coefficient decreases considering the output Y 1 charge density. The same happens for the current density Y 2 for temperature above 80 • C. From room temperature to 70 • C, the charge density is impacted by the release of charges from deep traps, while for a temperature higher than 70 • C the charge density is linked to the mobility of the charges in the shallow traps. For a temperature below 20 • C, the charge density is only related to the injection phenomena, charges tend to be close to the electrodes and to remain there. The influence of de-trapping barrier height on the charge density increases over time to reach 70% at one hour of charging time for a given temperature of 40 • C and a

Figure 4 .

 4 Figure 4.16: Evolution of the first order Sobol indexes of the de-trapping barrier height w tr .

Figure 4 .

 4 Figure 4.17: Evolution of the first order Sobol indexes of the trapping coefficient B.

Figure 4 .

 4 Figure 4.18: Principle of the optimization technique.

Backward Propagation method: Variance Partitioning Contents 2

.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47 2.2 Problem definition and notation . . . . . . . . . . . . . . . . . 48 2.3 Backward propagation: Variance Partitioning . . . . . . . . . 50 2.3.1 Step 1: Output variance in terms of V and R . . . . . . . . . . 50 2.3.2 Step 2: Solving nonlinear least square problem . . . . . . . . . 54 2.4 Applications . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59 2.4.1 First example . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59 2.4.2 Second example . . . . . . . . . . . . . . . . . . . . . . . . . . . 60 2.5 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61

Sobol Indexes in the Discrete Case Contents 3

  The AFM process . . . . . . . . . . . . . . . . . . . . . . . . .

	3.2.1

.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2 What is an EFDC . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2.2 Acquiring the EFDC . . . . . . . . . . . . . . . . . . . . . . . . 3.3 The model inputs . . . . . . . . . . . . . . . . . . . . . . . . . . 3.4 Experimental procedures . . . . . . . . . . . . . . . . . . . . . . 3.5 EFDC as a logistic law . . . . . . . . . . . . . . . . . . . . . . . 3.6 Sobol indexes of the logistic parameters . . . . . . . . . . . . . 3.6.1 First order Sobol indexes . . . . . . . . . . . . . . . . . . . . . 3.6.2 Second order Sobol indexes . . . . . . . . . . . . . . . . . . . . 3.7 Design Of Experiment (DOE) . . . . . . . . . . . . . . . . . . . 3.7.1 The basics of DOE . . . . . . . . . . . . . . . . . . . . . . . . . 3.7.2 The main effect plots . . . . . . . . . . . . . . . . . . . . . . . . 3.7.3 The interaction effect plots . . . . . . . . . . . . . . . . . . . . 3.8 Matrix Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.8.1 First order matrix model . . . . . . . . . . . . . . . . . . . . . 3.8.2 Second order matrix model . . . . . . . . . . . . . . . . . . . . 3.9 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

  ref , V ref ) of points in an initial 2D reference image I ref . Then, using the internal parameters of the camera, the grid is converted into mm. This obtained grid is considered as a 3D template T ref whose points have coordinates (X ref , Y ref , 0). 2. A 3D transformation is done to T ref consisting of a rotation and a translation. Then the grid obtained formed a deformed surface S def , whose points have coordinates:(X def , Y def , Z def ). After that, using the camera internal parameters, the deformed surface is projected in to a pixel image I def defined by the points (U def , V def ) 3. According to[Faugeras et al., 2004], a homography H exists which maps the points of I ref into I def . Thus we compute the H following the normalized direct linear transformation algorithm proposed in[START_REF] Hartley | [END_REF].

4. For a selected value of σ, a Gaussian noise is inserted to the gray level of the two images I ref and I def . The obtained noisy images are called I N ref and I N def respectively.

  • , 30 • ] , [30 • , 45 • ] , and [45 • , 60 • ]. In each of these cases we consider:

Table 4

 4 .1.

	Inputs	Notation	Units	Lower bound	Upper bound
	Barrier height for injection, w	X 1	eV	1.10	1.20
	Mobility, µ	X 2	m 2 .V -1 .s -1	10 -14	10 -12
	Trapping coefficient, B	X 3	s -1	5×10 -4	10
	De-trapping barrier height, w tr	X 4	eV	0.73	1

Table 4 .

 4

1: The range of variation of the four inputs.

1.5 Forward uncertainty propagation

1.8 Conclusion

citations pages20, 32, et

List of Algorithms

First order Sobol indexes of Z, f , and θ as a function of Z with [15 Note that we try to keep this manuscript self-contained as it may include some concepts from statistics (chapters 1 and 2), force spectroscopy (chapter 3), and computer vision (chapter 4). We try to make the information presented completely sufficient to understand how uncertainty is studied. However for further details, readers are invited to consult the references that will be mentioned in each chapter.

the corresponding input uncertainty is defined by:

(2.15)

where I 1 and I 2 are as defined in 2.1, and A 1 , • • • , A n , B 1 , • • • , B n are constants to be determined. Thus the expression of the output variance at the input point x becomes a

, and it can be written as:

So, to determine the input uncertainty it is enough to find

are 2n unknowns. For that we select different 2n points a (1) , • • • , a (2n) from D, and we define new residues by:

Then by solving the associated least square problem with the new residues, we get the

This gives the input uncertainty as a linear function of the input point using the expressions in (2.15). In fact this can be considered as a first order approximation of the uncertainty of the inputs. However, for most experimental instruments the uncertainty presented can be assumed to take constant values or at most linear since the manufacturing error of the instruments is usually constant or linear.

In the following section, three different examples for the application of our method are presented. In the first example the model is a polynomial function with mixed distributions of the uncertain inputs. In the second example, we consider a model presented ??,

where the input uncertainty is computed by inverting the model. The results obtained by our method are compatible with the results presented in the above reference.

Note that the model Fitted-4PL is not defined explicitly, however by the fitting of the EFDC by 4PL, data points of Fitted-4PL can be acquired.

Sobol indexes of the logistic parameters

Recall that sensitivity analysis studies how the variation in the output can be apportioned to the variation of the inputs. This leads to the determination of how the output is dependent on each of the inputs. Thus, sensitivity analysis allows the identification of the input (or set of inputs) that has the greatest influence on the output. In this manner, our aim is to find qualitative information about the influence of each of where V ar(F ij ) is the defined by:

For each output A, B, C, and D and for each of the inputs d, v, and w, the second order Sobol indexes are computed by substituting the above formulas according to a discrete random variables. So for instance for the output A and the second order Sobol index of the inputs d and w is:

where Apparent effect is for the interaction between v and w on D, its associated second order

Sobol index is about 0.4. The output C is slightly affected by the interaction between the inputs where its associated second order Sobol indexes are all about 0.2. Similarly the output B is slightly affected by the interaction between the inputs where its associated second order Sobol indexes are all about 0.07. The interaction between the inputs appears projection formula in order to form a system of partial differential equations. Some reformulations and a change of variables are done to the system, allowing it to be solved locally.

Then the obtained solution of the system is found analytically and it gives the depth of the deformed surface that is the main unknown in the 3D reconstruction problem. Both isometric and conformal deformations are treated by the SfT method, however we focus here on the method presented in [START_REF] Bartoli | Shape-from-template[END_REF] which is an analytical, isometric, and local method. The following paragraphs give a detailed presentation of the concept of the SfT method in the case of an isometric deformation as described in [START_REF] Bartoli | Shape-from-template[END_REF].

Problem definition

The reconstruction problem to be solved by SfT is summarized in Fig. 4.6 with the notation In Fig. 4.6 the sketched forms are described in the following manner:

• S ⊂ R 3 is the unknown 3D deformed surface.

• J ⊂ R 2 is the given 2D image of S .

• T ⊂ R 3 is the given 3D template of S .

Conclusion and Perspectives

In this dissertation our aim was to make a contribution to the domain of uncertainty analysis especially for the backward uncertainty propagation. To this end, we studied different methods in this domain. We concentrated the most on the Sobol sensitivity analysis method where we developed its application framework. In addition, we derived a new backward propagation method that, unlike previous methods, does not rely on any preliminary approximations of the input uncertainty. The fruitful results obtained allowed us to prepare and publish different papers in different scientific journals. In the following we give a compact summary of what we did and how the conclusions were derived, and at the end we set our short and long term perspectives for the future work.

In

Distributed Random Variables

Let X be a normally distributed random variable with mean µ and variance denoted V .

Then the moment generating function of X is a function M X : [0, +∞) → R defined by:

With easy symbolic computation, the integral in A.1 can be solved and the moment generating function ends by:

Recall that whenever the moment generating function exists in some neighborhood of 0, the moments of the random variable are expressed in terms of the derivatives of the moment function at t = 0 [Heathcote, 2000]. So if m n denotes the n-th non central moment of X, then m n is written as:

Accordingly, the non central moments of X for any order can be directly found using A.3.

To ease the computation proces of the derivatives of the Cited page 70