Eduardo Haro Sandoval

Floriane Anstett-Collin

Michel Basset

Floriane Anstett

Sensitivity study of dynamic systems using polynomial chaos

Keywords: Global sensitivity analysis, dynamic system, polynomial chaos expansion, sensitivity functions

Global sensitivity has mainly been analyzed in static models, though most physical systems can be described by differential equations. Very few approaches have been proposed for the sensitivity of dynamic models and the only ones are local. Nevertheless, it would be of great interest to consider the entire uncertainty range of parameters since they can vary within large intervals depending on their meaning. Other advantage of global analysis is that the sensitivity indices of a given parameter are evaluated while all the other parameters can be varied. In this way, the relative variability of each parameter is taken into account, revealing any possible interactions. This paper presents the global sensitivity analysis for dynamic models with an original approach based on the polynomial chaos (PC) expansion of the output. The evaluation of the PC expansion of the output is less expensive compared to direct simulations. Moreover, at each time instant, the coefficients of the PC decomposition convey the parameter sensitivity and then a sensitivity function can be obtained. The PC coefficients are determined using non-intrusive methods. The proposed approach is illustrated with some well-known dynamic systems.

Introduction

Very often, the equations of a model involve unknown parameters which must be estimated from experimental data. A number of parameters are estimated with more or less precision, which may lead to unacceptable uncertainty on the model output. Among all the parameters, however, only few have a small or insignificant influence on the model response and therefore do not need to be determined precisely. On the other hand, some parameters are decisive for the model response and thus influence its uncertainty significantly. These parameters may require additional measurement data in order to be estimated with relatively high accuracy. To prepare and plan the experiments, it is necessary to distinguish the parameters with an insignificant influence on the response uncertainty, so as to set them at their nominal value in their interval of variation, thanks to the sensitivity analysis. Numerous studies have focused on the sensitivity analysis for static nonlinear models, for example, [START_REF] Blatman | Sparse polynomial chaos expansion and adaptive stochastic finite element using a regresion approach[END_REF][START_REF] Borgonovo | A new uncertainty importance measure[END_REF][START_REF] Cukier | Nonlinear sensitivity analysis of multiparameter model systems[END_REF][START_REF] Frey | Identification and review of sensitivity analysis methods[END_REF][START_REF] Saltelli | Sensitivity analysis for importance assessment[END_REF][START_REF] Saltelli | Sensitivity Analysis in Practice. A Guide to Assessing Scientific Models[END_REF][START_REF] Sobol | Sensitivity estimates for nonlinear mathematical models[END_REF]. The approaches may be local or global. Local approaches help to determine the impact of a small parameter variation around a nominal value [START_REF] Turyani | Local methods in sensitivity analysis[END_REF]. Global approaches also allow the determination of the same impact but by varying the parameter in its entire range of variation. Global methods are often based on the analysis of the output variance and are known as ANOVA (ANalysis Of VAriance) techniques [START_REF] Homma | Importance measures in global sensitivity analysis of model output[END_REF][START_REF] Saltelli | A quantitative model independent method for global sensitivity analysis of model output[END_REF][START_REF] Sobol | Sensitivity estimates for nonlinear mathematical models[END_REF]. More recently, sensitivity moment-independent methods have been used, where emulation model is used to compute densitybased sensitivity measure [START_REF] Borgonovo | Model emulation and moment-independent sensitivity analysis: An application to environmental modelling[END_REF]. The emulator is the one of [START_REF] Ratto | Using recursive algorithms for the efficient identification of smoothing spline anova models[END_REF].

The model function is split into a sum of functions of increasing dimension [START_REF] Sobol | Sensitivity estimates for nonlinear mathematical models[END_REF]. This decomposition, known as High Dimensional Model Representation (HDMR), performs the separation of the effects of different input parameters, which are transmitted in the decomposition of the variance. The present study exclusively focuses on global approaches. In order to quantify the contribution of a parameter to the output variance, a sensitivity index is calculated, often analytically when the model function is known and relatively simple. However, some models may be complex with a high number of parameters so that analytical calculations of the sensitivity indices become time consuming or even impossible. It is therefore necessary to estimate them [START_REF] Cukier | Nonlinear sensitivity analysis of multiparameter model systems[END_REF][START_REF] Jacques | Sensitivity analysis in presence of model uncertainty and correlated inputs[END_REF][START_REF] Mara | Application of global sensitivity analysis of model output to building thermal simulations[END_REF][START_REF] Mckay | Evaluating prediction uncertainty in simulation models[END_REF][START_REF] Sobol | Sensitivity estimates for nonlinear mathematical models[END_REF]. Very often, they are computed using Monte Carlo simulations, but for computationally demanding models, this can become intractable. To overcome this drawback, the model of interest is replaced by an analytical approximation, called metamodel, for example, by polynomials which are less expensive. The sensitivity indices are then obtained straightforwardly from the algebraic expression of the coefficients of the polynomial expansion. The polynomial chaos (PC) decomposition is an example of such metamodelling. The PC expansion appeared in the 1930's as an effective means to represent stochastic processes in mechanics [START_REF] Wiener | The homogeneous chaos[END_REF]. It is based on a probabilistic framework and represents amounts from stochastic spectral expansions of orthogonal polynomials [START_REF] Ghanem | Propagation of probabilistic uncertainty in complex physical systems using a stochastic finite element approach[END_REF][START_REF] Ghanem | Stochastic finite elements -A spectral approach[END_REF]. It has recently been used in an original manner for sensitivity analysis purposes in [START_REF] Blatman | Sparse polynomial chaos expansion and adaptive stochastic finite element using a regresion approach[END_REF][START_REF] Crestaux | Polynomial chaos expansion for sensitivity analysis[END_REF][START_REF] Sudret | Global sensitivity analysis using polynomial chaos expansion[END_REF][START_REF] Blatman | Efficient computation of global sensitivity indices using sparse polynomial chaos expansions[END_REF][START_REF] Buzzard | Global sensitivity analysis using sparse grid interpolation and polynomial chaos[END_REF][START_REF] Rahman | Global sensitivity analysis by polynomial dimensional decomposition[END_REF]. The great advantage of PC-based sensitivity approaches is that the full randomness of the response is contained in the set of the expansion coefficients. On the other hand, the analysis of influential parameters is also important for dynamic models since most physical systems (biological, mechanical, electrical and so on) can be described by differential equations. Very few approaches have been proposed in the literature for the sensitivity analysis of dynamic models and the proposed ones are based on local derivatives or on one-at-a-time approaches [START_REF] Pauw | Practical aspects of sensitivity function approximation for dynamic models[END_REF][START_REF] Perry | Principal components analysis in sensitivity studies of dynamic systems[END_REF]. However, for some applications, mechanical or biological ones for instance, it can be of great importance to consider the entire uncertainty range of parameters since they can vary within large intervals depending on their meaning. Another advantage of global sensitivity analysis is that the sensitivity estimates of individual parameters are evaluated while all the other parameters are varied. In this way, the relative variability of each parameter is taken into account, thus revealing any existing interactions. The global sensitivity analysis for dynamic models is addressed in this paper. In [START_REF] Witteveen | Modeling arbitrary uncertainties using Gram-Schmidt polynomial chaos[END_REF][START_REF] Xiu | The Wiener-Askey polynomials chaos for stochastic differential equations[END_REF], the PC expansion for stochastic differential equations has been studied to represent the model output and to get its statistic properties, but the parameter sensitivity has not been dealt with. Based on these studies, an original approach using the PC decomposition of the output is investigated here, to calculate the parameter sensitivity for dynamic models. At each time instant, the PC coefficients of the decomposition convey the parameter sensitivity and then a sensitivity function of each parameter can be obtained from the algebraic expression of the coefficients. The PC coefficients are determined either by regression or projection techniques which have the advantage of being nonintrusive methods. The proposed approach is illustrated with the well-known mass-spring-damper and DC motor systems. The outline of this paper is as follows. The sensitivity functions for dynamic systems are presented in Section 2 and the PC expansion for the output of a differential equation in Section 3. Section 4 is focused on the determination of the PC coefficients. Section 5 proposes a PC-based approach to the estimation of the sensitivity functions, which approach is summed up in Section 6. Finally, Section 7 presents an analytical test case to show the convergence of the numerical results. Moreover, the provided approach is applied on some representative dynamic physical systems.

Global sensitivity analysis

Consider the following stochastic differential equation:

L (t, ω, p(ω); y) = f (t, ω, p(ω)) (1)
where L is a linear or non-linear differential operator and p(ω) = (p 1 (ω), . . . , p n (ω)) with p i (ω), i = 1, . . . , n, the n unknown parameters, considered as uniformly random and independent variables, defined on the unit cube K. The stochastic variable ω is used to indicate the randomness of the input variable p. For the sake of simplicity, ω will be omitted in the following and p stands for p(ω).

The solution y = y(t, p), corresponding to the output of the model, can be decomposed into summands of increasing dimension [START_REF] Sobol | Sensitivity estimates for nonlinear mathematical models[END_REF], at each time instant:

y(t, p) = f 0 (t) + n i=1 f i (t, p i ) + n-1 i=1 n i<j f ij (t, p i , p j ) + . . . + f 1...n (t, p 1 , . . . , p n ) (2)
where y(t, p) ∈ R the model output is assumed continuous, derivative and square-integrable. The term f 0 (t) is the mean value of the output at each time instant:

f 0 (t) = K n y(t, p)dp (3)
The summands of equation ( 2) are given by:

f i (t, p i ) = E[y(t, p)|p i ] -f 0 (t) f ij (t, p i , p j ) = E[y(t, p)|p i , p j ] -f i (t) -f j (t) -f 0 (t) . . . (4)
where

1 0 f i1...is (t, p i1 , . . . , p is )dp i k = 0, (5)
with k ∈ {i 1 , . . . , i s } and 1 ≤ i 1 ≤ . . . ≤ i s ≤ n. Due to equation ( 5), the summands are orthogonal to each other:

K n f i1...is (t, p i1 , . . . , p is )f j1...jr (t, p j1 , . . . , p jr )dp = 0 (6)
for {i 1 , . . . , i s } = {j 1 . . . , j r }.

There are infinite possible decompositions but only one is satisfying Sobol' orthonormality condition [START_REF] Saltelli | Sensitivity analysis for importance assessment[END_REF]. Moreover, the variance of the output, denoted V (t), is given by:

V (t) = K n y 2 (t, p) -f 2 0 (t) dp (7)
The decomposition (2) leads to the following decomposition of the variance V (t):

V (t) = n i=1 V i (t) + n-1 i=1 n i<j V ij (t) + . . . + V 1...n (t) (8)
with:

V i (t) = V [E[y(t, p)|p i ]] V ij (t) = V [E[y(t, p)|p i , p j ]] -V i (t) -V j (t) . . . V 1...n (t) = V (t) - n i=1 V i (t) - 1≤i<j≤n V ij (t) -• • • - 1≤i1<...<in-1≤n V i1...in-1 (t) (9) where V [E[y(t, p)|p i ]] (resp. V [E[y(t, p)|p i , p j ]]
) is the variance of the conditional expectation of y(t, p) when p i is set (resp. p i and p j are set).

Since the output y(t, p) varies with time, there is one value for the classic Sobol sensitivity indices at each time instant, thus leading to sensitivity functions. In the same manner as for static models, the sensitivity functions of parameter p i are obtained by renormalising [START_REF] Homma | Importance measures in global sensitivity analysis of model output[END_REF] with the total variance V (t). Thus, the first order sensitivity function, denoted S i (t), is defined as follows:

S i (t) = V i (t) V (t) (10)
The first order sensitivity function S i (t) represents the main effect of the parameter p i which corresponds to its contribution alone. The value of S i (t), at each time instant, lies between 0 and 1. The closer to 1 its value is, the more parameter p i contributes to the total variance of the output. The sensitivity functions of higher orders, denoted S i1...is (t), are defined as:

S i1...is (t) = V i1...is (t) V (t) , 1 ≤ i 1 ≤ . . . ≤ i s ≤ n (11)
The sensitivity functions S i1...is (t) represent the collective contribution of the parameters p i1 , . . . , p is , which corresponds to the influence of the interactions of these parameters. By normalising equation ( 8) with the total variance V (t), the following expression is obtained:

1 = n i=1 S i (t) + n-1 i=1 n i<j S ij (t) + . . . + S 1...n (t) (12)
The total sensitivity function of parameter p i , denoted S T i (t), represents the contribution of p i alone in addition to the contribution due to its interaction with the other parameters [START_REF] Homma | Importance measures in global sensitivity analysis of model output[END_REF]. This total sensitivity function is defined by:

S T i (t) = S i (t) + n j=1 j =i S ij (t) + n j=1 j =i n k=j+1 k =i j =i S ijk (t) + . . . + S ijk...n (t) (13)
The previous sensitivity functions can be calculated using Monte Carlo simulations. However, the number of model evaluations can become very large and the calculations intractable, especially as they have to be performed at each time instant. To overcome this problem, the model response is substituted by its PC expansion, which is an analytical representation in an orthonormal polynomial basis. Therefore, evaluating the PC expansion of the output is less expensive compared to direct simulations. Furthermore, the sensitivity functions can be evaluated as analytical expressions of the PC coefficients. In the next section, the PC expansion of the model output is detailed.

PC expansion of the model output

The beginnings of the polynomial chaos can be traced back to Wiener [START_REF] Wiener | The homogeneous chaos[END_REF]. He suggested that the spectral expansion of Hermite polynomials in terms of Gaussian random variables could be used to represent certain stochastic processes. It was then shown that the homogeneous chaos expansion could be used to approximate any function in the Hilbert space of square-integrable functions. As a result, any second-order random variable with finite variance can be represented by a spectral expansion standing on a basis of Hermite polynomials. In [START_REF] Witteveen | Modeling arbitrary uncertainties using Gram-Schmidt polynomial chaos[END_REF][START_REF] Xiu | The Wiener-Askey polynomials chaos for stochastic differential equations[END_REF], this is extended to differential equations. In this case, any second-order random variable X(t, ω) can be expanded as:

X (t, ω) = ∞ j=0 x j (t) Φ j (ξ (ω)) (14)
where x j (t) are unknown deterministic coefficients, ξ(ω) = (ξ 1 (ω), . . . , ξ n (ω)) the n standard normal independent random variables (zero mean and unit standard deviation) and Φ j the multivariate Hermite polynomials. The basis is orthogonal with respect to the Gaussian measurement. As previously, for the sake of simplicity, ω will be omitted and ξ stands for ξ(ω).

In practice, the PC expansion ( 14) is truncated to a finite number M + 1 of terms, i.e. the coefficients to be calculated. The degree d of the polynomials is related to M and to the number n of random variables ξ by:

M + 1 = (n + d)! n!d! (15)
It is worth noting that the required degree for the polynomials is not known a priori and must be tested to obtain the best accuracy. The multivariate polynomial Φ j is given by the tensor product of the corresponding one-dimensional Hermite polynomials φ a j k :

Φ j (ξ) = n k=1 φ a j k (ξ k ) (16)
with a j k the degree of the one-dimensional Hermite polynomials, such that

|a j | = n k=1 a j k ≤ d, j = 0, . . . , M .
The Hermite expansion above has been used effectively for equations with Gaussian inputs. But, according to the theorem of Cameron and Martin [START_REF] Xiu | Modeling uncertainty in flow simulations via generalized polynomial chaos[END_REF], this function converges for arbitrary random processes. For example, normalised Legendre polynomials can be associated with a uniform distribution. Table 1 sums up the type of distribution and the associated polynomials.

Assuming that the model output y(t, p) is square-integrable and the parameters independent, y(t, p) can also be decomposed in PC: As stated in Section 2, the model parameters p are uniformly independent random variables. Therefore, the one-dimensional polynomials φ a j k in ( 17) are Legendre polynomials, defined by the following recurrence equation:

y(t, p) ≈ M j=0 α j (t)Φ j (p), Φ j (p) = n k=1 φ a j k (p k ) (17)
φ 0 (p k ) = 1 φ i+1 (p k ) = 2i + 1 i + 1 p k φ i (p k ) - i i + 1 φ i-1 (p k ), for i > 1 (18)
Moreover, Table 1 emphasizes that the Legendre polynomials are orthogonal with respect to a uniform probability measurement over [-1; 1], while the parameters p k , k = 1, . . . , n, lies in [0; 1]. Consequently, the parameters are first mapped onto the interval [-1; 1], before completing the PC expansion.

The deterministic coefficients α j (t), j = 1, . . . , M , are unknown. The next section is focused on their determination.

Determination of the PC coefficients

Once the structure of the PC expansion is obtained for the model output, as explained in the previous section, the deterministic coefficients α j (t) must be computed. There exist two types of methods for this, the intrusive and the non-intrusive ones. Historically, the intrusive Stochastic Galerkin (SG) method was used [START_REF] Field | Numerical methods to estimate the coefficients of the polynomial chaos expansion[END_REF][START_REF] Xiu | The Wiener-Askey polynomials chaos for stochastic differential equations[END_REF][START_REF] Xiu | Modeling uncertainty in flow simulations via generalized polynomial chaos[END_REF]. A Galerkin projection helps to minimize the error of the truncated expansion, and all the resulting coupled equations are solved for the coefficients of the expansion equations. SG method is well suited for the treatment of ordinary and partial differential equations as it has the ability to cope with a large nonlinear dependence of the solution on data from random patterns. However, SG method is intrusive in the sense that it must solve a system of coupled equations which require robust and efficient algorithms and the modification of an existing deterministic code. Often, the form of equations or code used to solve the deterministic equations is complex and it makes the implementation of the PC difficult, if not impossible. Therefore, non-intrusive methods have been proposed based on regression technique [START_REF] Berveiller | Stochastic finite element: a non intrusive approach by regression[END_REF][START_REF] Sudret | Application of a stochastic finite element procedure to reliability analysis[END_REF] or on projection method [START_REF] Crestaux | Polynomial chaos expansion for sensitivity analysis[END_REF][START_REF] Sudret | Global sensitivity analysis using polynomial chaos expansion[END_REF]. Instead of building prior expansions of uncertainty sources and of finding the coefficients by solving the corresponding equations, these approaches use interpolation methods and project a set of deterministic simulations, which are defined using samples of parameters chosen carefully, with a polynomial basis. These approaches are very useful when trying to quantify uncertainty in complex deterministic models which cannot easily be modified. Similarly to SG methods, non-intrusive methods reach a fast convergence when the solutions are sufficiently smooth in the random space.

Regression method

The PC coefficients are obtained through the regression of the mean-square optimization problem. In the following, consider a sample of size N for the output y and the parameters p, at each time instant. The k-th realisation of p and of y, at time t, is denoted, respectively, by p (k) and y (k) (t, p (k) ). The mean-square minimization helps to solve the following problem:

α(t) = arg min 1 N N k=1   y (k) (t, p (k) ) - M j=0 α j (t)Φ j p (k)   2 (19)
where α(t) = (α 0 (t), . . . , αM (t)) is the vector of PC coefficients. The notation α(t) is employed to distinguish the estimated value of α(t) from its true one. Denoting by Φ the matrix whose elements are given by Φ kj = Φ j (p (k) ), k = 1, . . . , N and j = 0, . . . , M , by Φ T its transpose and by y(t) = (y (1) (t, p (1) ), . . . , y (N ) (t, p (N ) ) the vector containing the N realisations of y, at time t, the solution of ( 19) is:

α(t) = Φ T Φ -1 Φ T y(t) (20)
It is generally advisable to use an over-sampling to determine N , resulting in a least squares solution for the over-determined system. It can be calculated by N = 2 (M + 1) [START_REF] Eldred | Recent advances in non-intrusive polynomial chaos and stochastic collocation methods for uncertainty analysis and design[END_REF], or N = (n -1) (M + 1) [START_REF] Sudret | Global sensitivity analysis using polynomial chaos expansion[END_REF]. As stated in [START_REF] Sudret | Global sensitivity analysis using polynomial chaos expansion[END_REF][START_REF]Non linear non intrusive stochastic finite element method -application to a fracture mechanics problem[END_REF], taking more points does not improve the accuracy of the results.

Projection method

The projection method takes advantage of the orthogonality properties of the PC representation computing the expansion coefficients from some evaluations of the model. Multiplying both sides of equation ( 17) by Φ l (p) and taking the inner product, it reads:

y(t, p), Φ l (p) ≈ M j=0 α j (t) Φ j (p) , Φ l (p) (21)
where •, • is the inner product in the Hilbert space. Orthogonality property is now used. Indeed, polynomials Φ j (p) form an orthogonal basis so that:

Φ j (p), Φ l (p) = δ jl ( 22
)
where δ jl is the Delta Kronecker symbol, δ jl = 1 if j = l and 0 otherwise. Thus, from [START_REF] Blatman | Efficient computation of global sensitivity indices using sparse polynomial chaos expansions[END_REF], the coefficients αj (t) can be obtained:

αj (t) = y(t, p), Φ j (p) Φ j (p) , Φ j (p) , ∀j (23)
In the latter expression the denominator is already known and the numerator can be computed with the multidimensional integral:

y(t, p), Φ j (p) = K n y(t, p)Φ j (p) h(p) dp (24)
where h(p) is the joint probability density function of p. In order to obtain the j-th coefficients, it can be shown that M + 1 multidimensional integrals must be solved. There exist different methods to numerically compute [START_REF] Pauw | Practical aspects of sensitivity function approximation for dynamic models[END_REF], random sampling, quadrature product and Smolyak sparse grids.

Random sampling

The Random sampling is the Monte Carlo approach. Here the integral evaluation is equivalent to compute the mean of the function product (left member of equation ( 24)) for each term in the expansion. To have appropriate estimates a total-order chaos expansion must be used and the samples must have the density of the weighting function. In practice, this method can be intractable because the number of samples required for acceptable convergence can be large.

Quadrature

Consider the use of a numerical integration technique to estimate the coefficients of the expansion. The integral of (24) can be computed as a finite weighted summation evaluated at special points (the integration points). The theory of quadrature approaches is based on the use of orthogonal polynomials of the PC basis. This numerical integration method can be expressed in the form:

D n y(t, p)Φ j (p) h(p) dp ≈ N k=1 y t, p (k) Φ j p (k) W (k) (25)
with p (k) being the integration points, W (k) the weights and N the number of integration points. These methods require the smoothness of the function.

The Gauss-Hermite integration of this equation requires N n computations, where n is the number of random variables. As for the number of the integration points, a minimal Gaussian quadrature order of N = d + 1 will be required to obtain a good accuracy in the coefficients, where d is the degree of the polynomials [START_REF] Crestaux | Polynomial chaos expansion for sensitivity analysis[END_REF][START_REF] Eldred | Recent advances in non-intrusive polynomial chaos and stochastic collocation methods for uncertainty analysis and design[END_REF][START_REF] Field | Numerical methods to estimate the coefficients of the polynomial chaos expansion[END_REF][START_REF] Sudret | Global sensitivity analysis using polynomial chaos expansion[END_REF].

Smolyak sparse grids

When the number of parameters is very large and the number of evaluations of the integral (25) becomes impractical, a solution is to use Smolyak cubature [START_REF] Crestaux | Polynomial chaos expansion for sensitivity analysis[END_REF][START_REF] Eldred | Recent advances in non-intrusive polynomial chaos and stochastic collocation methods for uncertainty analysis and design[END_REF]. Smolyak formulas are linear combinations of the products expressed in equation [START_REF] Perry | Principal components analysis in sensitivity studies of dynamic systems[END_REF], which only takes into account products with relatively few integration points. These methods are based on the differences between quadratures of successive levels. The one-dimensional quadrature formula of level l, with points p and weights W , respectively, writes:

D n y(t, p)Φ j (p) h(p) dp ≈ N k=1 y t, p (k) Φ j p (k,l) W (k,l) ≈ U 1 l (26)
Setting U 1 0 = 0 and with i ≥ 1 an increasing sequence of integration formulas, the difference quadrature formula is defined as:

∆ 1 i≥1 = U 1 i -U 1 i-1 (27)
With a set of multi-indices i = (i 1 , . . . , i n ), the n-dimensional difference formula is the tensor product of one-dimensional difference quadrature. Then, the ndimensional Smolyak's cubature formula of level l is constructed by the sum of tensor products of difference quadratures, over a set of multi-indices i:

D n y(t, p)Φ j (p) h(p) dp ≈ |i|≤l+n-1 ∆ i1 ⊗ • • • ⊗ ∆ in (28)

PC-based sensitivity functions

Consider the PC expansion of the model output [START_REF] Ghanem | Propagation of probabilistic uncertainty in complex physical systems using a stochastic finite element approach[END_REF] with the coefficients computed with [START_REF] Crestaux | Polynomial chaos expansion for sensitivity analysis[END_REF] or [START_REF] Rahman | Global sensitivity analysis by polynomial dimensional decomposition[END_REF]. Due to the orthogonality of the basis, it can be shown that the mean ŷ(t) and the variance V (t) of the output are, respectively, given by: ŷ(t) = α0 (t) (29)

V (t) = M j=1 α2 j (t) Φ j (p), Φ j (p) (30)
As explained in [START_REF] Sudret | Uncertainty propagation and sensitivity analysis in mechanical models -Contributions to structural reliability and stochastic spectral methods, habilitation à diriger des recherches[END_REF], to compute the sensitivity function S i (t) (10) of parameter p i , it is necessary to reorder the expansion [START_REF] Ghanem | Propagation of probabilistic uncertainty in complex physical systems using a stochastic finite element approach[END_REF] in order to separate the different contributions -single and collective -of each parameter, as (2). To do so, I k1,...,ks must define the set of multi-indices j such that :

I k1,...,ks = {0 ≤ a j k ≤ d, a j k = 0 ∀k ∈ {1, . . . , n}\{k 1 , . . . , k s }} (31)
where a j k is the one-dimensional polynomial degree. The set I i corresponds to the polynomials depending on parameter p i only, with the exception of all the others. Using this notation, the expansion ( 17) can be written as:

y (t, p) ≈ α 0 (t)Φ 0 (p) + n i=1 j∈Ii α j (t)Φ j (p i ) + 1≤i1≤i2≤n j∈Ii 1 ,i 2 α j (t)Φ j (p i1 , p i2 ) + • • • + 1≤i1≤...≤i2≤n j∈Ii 1 ,...,is α j (t)Φ j (p i1 , . . . , p is ) + • • • (32)
The first-order sensitivity function can be obtained straightforwardly from [START_REF] Eldred | Recent advances in non-intrusive polynomial chaos and stochastic collocation methods for uncertainty analysis and design[END_REF] with the estimated coefficients αi (t) (equation ( 19) or ( 23)). Thus, the estimated first-order sensitivity function Ŝi (t) of parameter p i is given by:

Ŝi (t) = j∈Ii α2 j (t) Φ j (p i ), Φ j (p i ) M j=1 α2 j (t) Φ j (p), Φ j (p) (33)
The estimated sensitivity functions of a higher order, Ŝi1,...,is (t), reflecting the interactions effect, can be obtained in the same manner:

Ŝi1,...,is (t) = j∈Ii 1 ,...,is α2 j (t) Φ j (p i1 , . . . , p is ), Φ j (p i1 , . . . , p is ) M j=1 α2 j (t) Φ j (p), Φ j (p) (34)
Finally, the estimated total sensitivity function, ŜT i (t) is given by:

ŜT i (t) = j∈I i + α2 j (t) Φ j (p i ), Φ j (p i ) M j=1 α2 j (t) Φ j (p), Φ j (p) (35)
with I i + the set of multi-indices j defined by:

I i + = {0 ≤ a j ≤ p, a j i = 0, ∀k = i} (36)
The set I i + corresponds to polynomials Φ i depending on parameter p i and possibly on the other parameters.

The great advantage of using the PC-based approach to the computation of the sensitivity functions is that the calculation is straightforward once the PC expansion is obtained. Then, this approach will be efficient if the PC expansion is accurate and inexpensive.

The next section gives a summary of the proposed approach.

Summary of the proposed approach

The computation of the sensivity functions of the parameters for the dynamic model ( 2), requires the following steps: a) PC decomposition 1. Polynomial type Choose the polynomial type associated to the parameter distribution, see Table 1.

If necessary, map the parameters onto the required interval given in Table 1, in order to get an orthogonal basis for the chosen polynomials.

PC degree

Set the required degree d for the PC.

PC truncation

Compute the number M + 1 of terms in the expansion [START_REF] Ghanem | Propagation of probabilistic uncertainty in complex physical systems using a stochastic finite element approach[END_REF], according to (15).

PC coefficients

Compute the PC coefficients according to [START_REF] Crestaux | Polynomial chaos expansion for sensitivity analysis[END_REF] or [START_REF] Rahman | Global sensitivity analysis by polynomial dimensional decomposition[END_REF].

b) Sensitivity functions

Once the PC decomposition of the output is completed, compute the sensitivity functions according to [START_REF]Non linear non intrusive stochastic finite element method -application to a fracture mechanics problem[END_REF], [START_REF] Sudret | Uncertainty propagation and sensitivity analysis in mechanical models -Contributions to structural reliability and stochastic spectral methods, habilitation à diriger des recherches[END_REF] and [START_REF] Baudin | Polynômes de chaos sous Scilab via librarie NISP[END_REF], for each parameter.

Application examples

A simple test case is provided in order to compute analytically the sensitivity indices and to show the convergence of the numerical results obtained using the PC decomposition. Then, two examples of dynamic system are given to illustrate the contribution of this work.

Analytical test case

Consider the solution y(t) of a simple first-order differential equation of the form ẏ(t) = p 1 :

y(t) = p 1 t + p 2 (37)
The parameters p 1 and p 2 are assumed to be independent uniformly distributed over [c 1 ; d 1 ] and [c 2 ; d 2 ], respectively.

Analytical sensitivity indices

The mean value ȳ(t) of the output is given by (3) and can be calculated analytically:

ȳ(t) = c 2 + d 2 2 + (c 1 + d 1 )t 2 (38)
The output variance V (t) is given by ( 7) and its analytical expression is:

V (t) = (c 2 -d 2 ) 2 12 + (c 1 -d 1 ) 2 t 2 12 ( 39
)
The first-order sensitivity functions S p1 (t) and S p2 (t) of the parameters p 1 and p 2 are given by [START_REF] Saltelli | A quantitative model independent method for global sensitivity analysis of model output[END_REF] and their analytical expressions are, respectively:

S p1 (t) = (c 1 -d 1 ) 2 t 2 (c 2 -d 2 ) 2 + (c 1 -d 1 ) 2 t 2 (40)
S p2 (t) = (c 2 -d 2 ) 2 (c 2 -d 2 ) 2 + (c 1 -d 1 ) 2 t 2 (41)
It can be worth noting that model (37) is additive. Thus, the second-order index S p1p2 (t) is zero and the total indices, S T p1 (t) and S T p2 (t), are equal to the first-order indices, S p1 (t) and S p2 (t), respectively.

PC-based sensitivity indices

Here, the number of parameters is n = 2 and the degree used for the polynomials is d = 1 as the output is linear in the parameters. Thus, the number of coefficients to compute in the decomposition is M +1 = 3 [START_REF] Mckay | Evaluating prediction uncertainty in simulation models[END_REF]. The multivariate polynomials Φ j are given by ( 16):

Φ j (p 1 , p 2 ) = 2 k=1 φ a j k (p k ) (42)
where

|a j | = 2 k=1
a j k ≤ 1, j = 0, . . . , 2, and φ a j k are one-dimensional Legendre polynomials. The mean value ŷ(t), the variance V (t) and the sensitivity functions Ŝp1 (t) and Ŝp2 (t) are given, respectively, by ( 29), ( 30), ( 33), ( 34) and ( 35).

In the following, the numerical values are c 1 = 0.76, d 1 = 0.84, c 2 = 0.9 and d 2 = 1.1. The figure 1 represents the evolution of the output variance computed analytically vs using PC decomposition (top figure). The relative error, shown in figure 1 (down), is given, in percentage, by

100|V (t) -V (t)| V (t)
. This error does not exceed 6%, showing good accuracy of the results. The figures 2 and 3 represent, respectively, the sensitivity functions of p 1 and p 2 , computed analytically vs using PC decomposition (top) and the relative error of the sensitivity functions (down). These figures show that the relative errors, computed as previously, do not exceed 3% for p 1 and 7.5% for p 2 , showing satisfactory accuracy of the PC decomposition.

Mass-spring-damper system

Also to illustrate the contribution of this paper, an ideal mass-spring-damper is considered, with a mass m, a spring constant b and a viscous damping coefficient k, subject to a constant force F , as represented in figure 4. This system is described by the following differential equation:

mÿ(t) + b ẏ(t) + ky(t) = F (43)
where y(t), the output of the system, represents the mass position. This equation can be represented by a state space system. Defining the states as x 1 (t) = y(t) and x 2 (t) = ẏ(t), it reads:

ẋ1 (t) ẋ2 (t) = 0 1 - k m - b m x 1 (t) x 2 (t) + 0 1 m F ( 44
)
Figure 5 is the simulation when a step signal is applied on F , for the fixed values m = 1, b = 1 and k = 4.5.

The sensitivity of the parameters m, b and k on the mass position y(t) is considered. The proposed approach is applied in order to obtain the sensitivity functions.

p i -E[p i ] (b i -a i )/2 , ∀i ∈ {1, 2, 3} (45)
where b i and a i are, respectively, the maximum and minimum values of the parameter interval and E[p i ] the expectation of p i .

PC degree

Here, the number of parameters is n = 3 and the degree used for the polynomials is d = 4. This degree of expansion was sufficient due to the smoothness of the function, as shown in figure 6.

PC truncation

Consequently, the number of terms in the decomposition (the number of coefficients to be calculated) is M + 1 = 35, according to (15).

PC coefficients

After substituting (45) for p i , the multivariate polynomials Φ j are given by ( 16):

Φ j (p 1 , p 2 , p 3 ) = 3 k=1 φ a j k (p k ) (46)
where |a j | = 3 k=1 a j k ≤ 4, j = 0, . . . , 34, and φ a j k are one-dimensional Legendre polynomials. The 35 coefficients α j (t) are determined with the projection method because the given dynamic system is relatively smooth. Therefore, the coefficients are given by [START_REF] Rahman | Global sensitivity analysis by polynomial dimensional decomposition[END_REF]. The software Scilab toolbox NISP [START_REF] Baudin | Polynômes de chaos sous Scilab via librarie NISP[END_REF] helps to compute them. The error y(t) -

Sensitivity functions

The next step is to compute the sensitivity functions for the three parameters according to [START_REF]Non linear non intrusive stochastic finite element method -application to a fracture mechanics problem[END_REF], [START_REF] Sudret | Uncertainty propagation and sensitivity analysis in mechanical models -Contributions to structural reliability and stochastic spectral methods, habilitation à diriger des recherches[END_REF] and [START_REF] Baudin | Polynômes de chaos sous Scilab via librarie NISP[END_REF]. Figures 9, 10 and 11 represent, respectively, the first-order, the higher order and the total sensitivity functions. Several instants of time were chosen to cover the transient and the steady-state of the system: the functions were calculated every 0.2 seconds. For the first seconds of the transient, the most influential parameters on the mass position are obviously the mass and the damper, because they depend on the acceleration and the speed (see figure 9). When the speed reaches its maximum value and the acceleration is zero, approximately at 0.6 seconds, the mass and the spring become equally influent. From second 2 up to the steady state, the spring is clearly the most relevant parameter because the velocity and the acceleration tend to zero. The second order indices S 13 (mass-spring) and S 23 (damperspring) appear more clearly in the middle of the simulation (figure 10), when the position is already significant and the velocity and acceleration are still not zero.

DC motor system

Another example of a known dynamic system is the permanent magnet DC motor (figure 12). The parameters of this system are the internal resistance R, the inductance L, the motor constant k m , the counter electromotive constant k a , the mechanical inertia J and the mechanical damping b. The input is the voltage V applied to the motor and the output is the shaft speed ω(t). The torque is given by Γ = k m i(t) and the counter electromotive voltage by V b = k a ω(t). We assume that the motor is unloaded, so the equations of this system are:

V = Ri(t) + L di(t) dt + k a ω(t) J ω(t) = k m i(t) -bω(t)
Defining the states as x 1 (t) = i(t) and x 2 (t) = ω(t), the state space representation of this system is:

ẋ1 (t) ẋ2 (t) = -R L -ka L k m J - b J x 1 (t) x 2 (t) + 1 L 0 V (47)

PC degree

Here, the number of parameters is n = 6 and the degree used for the polynomials is d = 5. This degree is sufficient due to the smoothness of the function, as shown in figure 13. The error between the original system and the PC approximation is obtained to verify the accuracy of the PC expansion, as shown in figure 14.

PC truncation

The number of terms in the decomposition is M + 1 = 462.

PC coefficients

The multivariate polynomials Φ j are given by ( 16), with n = 6 and where In practice, these parameters can be neglected or settled to their nominal value. The inertia J, the resistance R and the motor constant k m are important during the transient. Then, their influence disappears when reaching the steady state. The parameter k a is the most influential parameter in the steady state, since it is a constant which is multiplied by the speed of the shaft.

|a j | =

Conclusion

This paper has presented the problem of the global sensitivity of dynamic systems which proves important when parametric uncertainty falls within a wide range for several parameters and when the sensitivity changes significantly from one operation point to another. The polynomial chaos expansion is used to calculate the sensitivity indices of the system. The expansion coefficients are calculated either through the technique of regression or projection, which are non-intrusive methods. As the sensitivity of a dynamic system changes over time, the sensitivity indices are calculated at different time instants. In this way, the global sensitivity functions of the dynamic system are obtained. A major advantage of this approach is that, if the function to be represented with the PC is smooth, the number of coefficients that must be calculated is small. So, computation time is small too, compared with the Monte Carlo simulation. The approach presents good results for relative simple dynamic models, as presented in this work. Further work will be carried out on more complex models (partial differential equations, for instance) to attest its efficiency and test its accuracy of PC decomposition in long term.

tFigure 1 :

 1 Figure 1: Evolution of the output variance (top) -Relative error for the output variance V (t) (down)

Figure 4 :Figure 5 :

 45 Figure 4: Mass-spring damper system

34 j=0α

 34 j (t)Φ j (p), represented in figure6, is low in the transient and stabilizes in the steady state, showing a good approximation of the output with the PC decomposition. Moreover, figures 7 and 8 show the PDF of the output, at different time instants.

Figure 6 :Figure 7 :Figure 8 :

 678 Figure 6: Error between the output and its PC-based approximation

Figure 9 :Figure 10 :

 910 Figure 9: 1st-order sensitivity functions

Figure 11 :Figure 12 :

 1112 Figure 11: Total sensitivity functions

Figure 13 represents

 13 Figure 13 represents the evolution of ω(t), when the step signal V = 9 is applied for the fixed values L = 0.3, R = 15, J = 0.2, b = 0.2, k m = 7.5 and k a = 3.

Figure 14 :

 14 Figure 14: Error between the output and its PC-based approximation

6 k=1 a j k ≤ 5 ,

 65 j = 0, . . . , 461, and φ a j k are one-dimensional Legendre polynomials. The error y(t) -461 j=0 α j (t)Φ j (p), represented in figure14, is less than 0.15%, showing a good approximation of the output with the PC decomposition.

Figure 15 :Figure 16 :

 1516 Figure 15: 1st-order sensitivity functions

  E[y(t, p)|p i ] (resp. E[y(t, p)|p i , p j ]) is the conditional expectation of y(t, p) when p i is set (resp. p i and p j are set). The integral of each summand f i1...is (t, p i1 , . . . , p is ) is zero:

Table 1 :

 1 Wiener-Askey polynomials

		Distribution	Polynomials	Support
	Continuous	Gaussian	Hermite	R
		Uniform	Legendre	[-1; 1]
		Gamma	Laguerre	[0; +∞[
		Beta	Jacobi	[-1; 1]
	Discrete	Poission	Charlier	{0, 1, . . . , }
		Negative binomial	Miexner	{0, 1, . . . , }
		Binomial	Krawtchouk {0, 1, . . . , N }
		Hypogeometric	Hahn	{0, 1, . . . , N }

  Figure 2: Sensitivity function Sp 1 (t) (top) -Relative error for Sp 1 (t) (down)

		1.5
										Analytical Sp1(t)
										PC Sp1(t)
	1 ( t)	1
	Sp	0.5
		0	2	4	6	8	10	12	14	16	18	20
	1 ( t)						t
	Re lat ive e r r or for Sp	0 2 4 6 8 10	2	4	6	8	10	12	14	16	18	20
							t

Acknowledgements

The authors wish to thank Michael Baudin (Consortium Scilab) for his valuable guidance in the use of software NISP.


Meiling Luo
email: meiling.luo@ranplan.co.uk

Nikolai Lebedev

Guillaume Villemaud

Guillaume De La Roche
email: guillaume.delaroche@mindspeed.com

Jie Zhang
email: jie.zhang@sheffield.ac.uk

Jean-Marie Gorce
email: jean-marie.gorce@insa-lyon.fr

On Predicting Large Scale Fading Characteristics with the MR-FDPF Method

The MR-FDPF (Multi-Resolution Frequency Domain Partial Flow) method is proven to be a fast and efficient method to simulate radio wave propagation. It is a deterministic model which can provide an accurate radio coverage prediction. In reality, radio channels have the nature of randomness due to e.g. moving people or air flow. Thus they can not be rigorously simulated by a pure deterministic model. However, it is believed that some statistics can be extracted from deterministic models and these statistics can be very useful to describe radio channels in reality. In this paper, large scale fading statistical characteristics are extracted based on the MR-FDPF method. They are validated by comparison to both the theoretical result and measurement. The match also demonstrates that MR-FDPF is capable of simulating large scale fading.

I. INTRODUCTION

Since the radio propagation channel is the foundation of wireless communications, a good understanding of the radio channel is necessary for reliable system design [START_REF] Rappaport | Wireless communications: principles and practice[END_REF]. For instance, an accurate knowledge of channel characteristics is required for cellular operators to optimize the coverage and maintain the interference at the lowest possible level.

Typically, the characteristics of the radio channel can be obtained by the following ways:

• Channel measurements: Conducting channel measurements is the most reliable way to know the radio channel, but it is not always possible to conduct extensive channel measurements for every radio system since measurement campaigns are usually time consuming and costly.

• Empirical channel models: Empirical models [START_REF] Rappaport | Wireless communications: principles and practice[END_REF] are constructed either based on simplifying assumptions concerning the physical geometry of the propagation environments or based on a best fit to measurements data conducted in a typical environment. However, empirical models suffer from a lack of accuracy as they do not take into account specific propagation environments.

• Deterministic channel models: Deterministic channel models typically include ray optical methods [START_REF] Balanis | Advanced Engineering Electromagnetics[END_REF] and finite difference methods [START_REF] Yee | Numerical solution of initial boundary value problems involving maxwell's equations in isotropic media[END_REF]. They are based on the theory of electromagnetic propagation. Deterministic models are known for their relatively high accuracy since they intrinsically take into account the propagation environment. For that reason, they are also called site-specific models [START_REF] Molisch | Wireless communications[END_REF]. However, deterministic models have the disadvantage of high computational load.

Among the three, deterministic channel models are a good choice considering the trade off between accuracy and cost. Therefore, in this paper, the MR-FDPF method [START_REF] Gorce | Deterministic approach for fast simulations of indoor radio wave propagation[END_REF][6] [START_REF] Luo | Estimating channel fading statistics based on radio wave propagation predicted with deterministic mrfdpf method[END_REF][8], a finite difference method, is chosen as the propagation model. The MR-FDPF method is proven to be fast and efficient in providing the radio coverage prediction. The provided coverage has a high level of accuracy because it takes into account the specific propagation environment, and also because of its high spatial resolution. Since the MR-FDPF method is a deterministic channel model, the signal strength provided by it is deterministic, but not stochastic. However, in reality radio channels have the nature of randomness due to e.g. moving people or air flow. Therefore, even while possessing the property of high accuracy, the MR-FDPF method can not describe the random aspect of radio channels. However, it is believed that some useful statistics, like large scale fading statistics, can be extracted from the deterministic model. In this paper, we investigate this point and validate it by both the theoretical result and measurement.

The rest of the paper is organized as follows. Section II presents the concept of large scale fading. Then, in section III large scale fading characteristics are extracted based on the MR-FDPF method, followed by experimental evaluation in section IV. Finally, conclusion is given in section V.

II. LARGE SCALE FADING CHARACTERISTICS

When expressed in dB, the instantaneous path loss can be considered as the sum of the mean path loss, shadow fading and small scale fading as follows [START_REF] Molisch | Wireless communications[END_REF]

P L(d) = L(d) + X σ + F (1)
where P L(d), L(d), X σ and F denote the instantaneous path loss, the mean path loss, shadow fading and small scale fading associated with the Tx-Rx (transmitter and receiver) separation distance d, respectively. Unlike small scale fading which characterizes the rapid fluctuations of the received signal strength over very short travel distances (a few wavelengths), large scale fading, referring to the total effect combining the mean path loss and shadow fading shown in equation ( 1), characterizes signal strength variation over large distances [START_REF] Rappaport | Wireless communications: principles and practice[END_REF].

Typically the mean path loss L(d) is deterministic and is log dependent on the Tx-Rx separation distance d as follows [START_REF] Rappaport | Wireless communications: principles and practice[END_REF] L

(d) = L 0 + 10n • log 10 (d) (2)
where L 0 is a constant which accounts for system losses and n is the path loss exponent depending on the specific propagation environment. For instance in free space propagation n = 2.

The shadow fading X σ is a zero-mean Gaussian distributed random variable (in dB) with standard deviation σ. Small scale fading is typically described by Rayleigh distribution for NLOS (non-line-of-sight) propagation and Rice distribution for LOS (line-of-sight) propagation [START_REF] Rappaport | Wireless communications: principles and practice[END_REF].

Large scale fading characteristics are very useful for determining the coverage area or evaluating system performance. For instance, with large scale fading characteristics, system designers can determine how large the coverage area of a transmitter is and what fade margin is required to achieve a certain level of edge reliability, e.g., 90% edge reliability.

III. EXTRACTION OF LARGE SCALE FADING CHARACTERISTICS BASED ON THE MR-FDPF METHOD

In this section, a brief introduction of the MR-FDPF method is provided. We then investigate whether the MR-FDPF method is capable of simulating shadow fading phenomenon and how to extract fading characteristics based on the MR-FDPF method.

A. Multi-Resolution Frequency Domain Partial Flow method

MR-FDPF method is the frequency domain version of (time domain) Partial Flow method [START_REF] Gorce | Deterministic approach for fast simulations of indoor radio wave propagation[END_REF]. The Partial Flow method is a numerical solution of Maxwell's equations based on the cellular automata formalism [START_REF] Chopard | Lattice boltzmann method for wave propagation in urban microcells[END_REF] [START_REF] Chopard | Cellular automata modeling of physical systems[END_REF]. Detailed information about the Partial Flow method can be found in [START_REF] Luthi | Lattice wave automata: from radio wave to fracture propagation[END_REF] [START_REF] Chopard | Lattice boltzmann computations and applications to physics[END_REF]. For the reason of simplification, the Partial Flow method is transfered into frequency domain by taking the Fourier transform. Moreover, the computational load of the frequency domain Partial Flow method can be further reduced by introducing the multi-resolution structure [START_REF] Gorce | Deterministic approach for fast simulations of indoor radio wave propagation[END_REF].

Considering the high computational load, most of the time MR-FDPF method is restricted to 2 dimensions (2D). This introduces an additional approximation error when compared to a real 3D propagation. Therefore, a calibration process seems to be imperative [START_REF] Runser | Assessment of a new indoor propagation prediction method based on a multi-resolution algorithm[END_REF].

The calibration can be done in two steps. The first step is to estimate the constant offset as follows

∆Ψ = 1 m m k=0 (Ψ mes (k) -Ψ sim (k)) (3)
where Ψ mes (k) and Ψ sim (k) are the mean powers from measurements and simulations, respectively, and m is the number of samples. A constant offset always exists because of the numerical sources used in the MR-FDPF method, compared to the real transmitters in reality. An attenuation coefficient of the air α air is introduced by MR-FDPF in order to compensate the inaccuracy of the 2D model to simulate a real 3D propagation [START_REF] De La Roche | On predicting in-building wifi coverage with a fast discrete approach[END_REF]. Besides, the set of parameters (n mat , α mat ) of materials, where n mat and α mat are the reflection factor and the attenuation coefficient of materials respectively, should also be calibrated since we never know exactly the material and the width of walls. Therefore, the second step of the calibration process is to estimate the α air and (n mat , α mat ) by minimizing the cost-function Q defined by the root mean square error (RMSE) between measurements and predictions

Q = RM SE = 1 m m k=0 Ψ mes (k) -Ψ pred (k) 2 (4)
where

Ψ pred (k) = Ψ sim (k) + ∆Ψ (5)
are the mean powers from predictions. The minimization process is solved by the direct search algorithm "DIRECT" by Jones et Al. in [START_REF] Jones | Lipschitzian optimization without the lipschitz constant[END_REF]. In Fig. 1, radio coverage provided by the MR-FDPF method is presented, obtained from a narrowband signal at 2.45 GHz with 0 dBm transmitter power.

B. Capability of Simulating Shadow Fading Phenomenon

From a theoretical point of view, small scale fading and shadow fading can be averaged out over an area, e.g. a square area, of a proper size. In other words, neither small scale fading nor shadow fading can be observed if the averaging area is large enough. The obtained path loss curve after averaging over a large enough area will be very smooth because it only contains the mean path loss.

If MR-FDPF is capable of simulating the fading characteristics, it should show the averaging effect described above. Therefore, we investigate this by doing averaging over square areas of different sizes, from 5λ × 5λ to 50λ × 50λ, with a step of 5λ, where λ denotes the wavelength. The averaging effect from MR-FDPF simulation is shown in Fig. 2. In this figure, a shadowing effect is evidenced over 5λ × 5λ area, and 10λ × 10λ area too, but there is a trend: the larger the size of the averaging area is, the less the shadowing effect can be observed. When the size of the area is about 40λ × 40λ, the shadowing effect almost can not be observed, which means that the shadowing has been averaged out.

Fig. 2 demonstrates that MR-FDPF is capable of providing the shadow fading phenomenon.

C. Extraction of Large Scale Fading Characteristics

In this subsection, we describe the extraction of the large scale characteristics, namely the mean path loss and the shadow fading. In order to get the large scale characteristics, first we should average out the small scale fading. The obtained path loss after averaging out the small scale fading is usually called the local mean path loss. The most critical point for averaging out the small scale fading is the determination of the size of the averaging area.

According to [START_REF] Saunders | Antennas and propagation for wireless communication systems[END_REF][16], if the samples are expressed in dB, the number of samples associated with 90% confidence inter- val is N = 85. Provided that for MR-FDPF a square area can be chosen, the number of samples associated with each side of the square is thus √ N . In addition, for a Rayleigh distributed signal envelope, the uncorrelated distance for two adjacent samples is 0.38λ [START_REF] Parsons | The mobile radio propagation channel[END_REF]. Finally, we obtain the minimum length of the side of the square is 0.38λ × √ N ≈ 3.8λ. On the other side, if it is a LOS propagation, the received signal envelope may be Rice distributed rather than Rayleigh, then in this case a smaller sample size may be sufficient [START_REF] Parsons | The mobile radio propagation channel[END_REF].

Since the performance of MR-FDPF will be verified by comparison to the measurement conducted at Stanford by Dr. N. Czink [START_REF] Czink | radio measurement campaign: Measurement documentation[END_REF] which will be described in section IV, we run MR-FDPF simulation with the same scenario as in the measurement. The simulation frequency is 2.45GHz, and there are 8 transmitters and 8 receivers (i.e. 64 links in total). For more detailed scenario information, please refer to section According to the description above, the local mean path loss averaged over an area of 3.8λ×3.8λ is obtained shown in Fig. 3 by the red stars.

The local mean path loss includes both the mean path loss and the shadow fading. Since the mean path loss depends on the log distance shown in equation ( 2), we choose to get the L 0 and n by using the curve fitting tool of Matlab such that the difference between the local mean path loss and the estimated mean path loss is minimized in a mean square error sense. Substituting the estimated L 0 = 50.26 and n = 1.592 into equation ( 2), we obtain the mean path loss

L(d) = 50.26 + 10 × 1.592 • log 10 (d) (6)
Here n = 1.592 indicates that there exists a waveguide effect in the propagation. The obtained mean path loss is shown by the black continuous curve in Fig. 3. Finally, the shadow fading can be obtained by just subtracting the mean path loss from the local mean path loss, which is denoted by the blue circles in Fig. 3. In order to validate the shadow fading we have extracted, we compare its CDF (cumulative distribution function) with that of normal distribution. The comparison result is shown in Moreover, we also compute the standard deviation of the extracted shadow fading, which is 5.87 dB.

IV. EXPERIMENTAL EVALUATION

The results from the measurement campaign corresponding to "I2I moving receivers" scenario conducted at Stanford [START_REF] Czink | radio measurement campaign: Measurement documentation[END_REF] are exploited to verify the performance of the proposed approach based on the MR-FDPF.

A. The office scenario

The scenario was a typical 16 x 34 m office space made of 30 cubicles and 7 small separated rooms. 8 transmitters and 8 receivers were distributed in the office as illustrated in Fig. 5. All the transmitters were fixed in their locations, whereas all the receivers were randomly moved inside their cubicles when conducting the measurements. All the transmitters and receivers were equipped with omnidirectional antennas. Four materials were mainly used in the office, i.e., concrete for the main walls, plaster for the internal walls, glass for the external glass wall and wood for the cubicles located in the central part of the office.

B. Measurement setup 8 x 8 MIMO channels at a center frequency of 2.45 GHz were measured simultaneously with a RUSK MEDAV channel sounder [START_REF]Rusk medav channel sounders[END_REF]. In the measurement, 120 time blocks covering a total time of 32 seconds and 220 frequency bins covering a bandwidth of 70 MHz were recorded. For large scale fading study, narrowband measurements are sufficient, hence the measurement data is restricted to the center frequency. Note that the averaging over the 120 time blocks implies the averaging over space since the receivers were moving when recording.

C. Measurement results

In this subsection, we extract large scale fading characteristics from the measurement. First of all, we obtain the local mean power by averaging over the 120 time blocks, which is shown by the red stars in Fig. 6. In the same way as in the simulation, the deterministic mean path loss is obtained by curve fitting

L(d) = 47.25 + 10 × 1.442 • log 10 (d) (7)
The mean path loss is shown in Fig. 6 by the black continuous curve. Once we get the mean path loss, the shadow fading is obtained by subtracting the mean path loss from the local mean path loss, which is denoted by the blue circles in Fig. 6. The CDF comparison between the extracted shadow fading and normal distribution is shown in Fig. 7. This good match demonstrates that the proposed approach for extracting the shadow fading is effective. The standard deviation of the extracted shadow fading from the measurement is 7.66 dB.

D. Match between the simulation and measurement

Compare the parameters extracted from the MR-FDPF method and measurement in Tab. I, we can see that they fit each other very well. The path loss exponents n both from the MR-FDPF method and measurement are smaller than 2, which indicates that there exists a waveguide effect in the propagation. The standard deviation of shadowing from simulation is slightly lower than that from measurement. We suggest that this is due to the unmodeled furniture which would contribute more to fading. The match between the simulation and measurement demonstrates that the MR-FDPF method is capable of simulating large scale fading characteristics.

V. CONCLUSION AND PERSPECTIVES

Although deterministic channel models possess the advantage of high accuracy due to their ability to take into account the specific propagation environment, they can not rigorously simulate realistic radio channels presenting random phenomenon. However, some useful statistics can be extracted from deterministic channel models. Basically, deterministic channel models combine all propagation phenomenons, i.e., the deterministic mean path loss, shadow fading and small scale fading. Thus, in order to extract fading statistics, we should first try to separate them. In this paper, large scale fading characteristics are extracted based on the MR-FDPF method. The extracted fading characteristics are validated by comparison to both the theoretical results and measurement.

The match demonstrates that MR-FDPF method is capable of simulating large scale fading characteristics.

In the future, more measurements would be conducted in order to validate the extraction approach based on MR-FDPF in other scenarios and/or frequency bands. Furthermore, although this paper is mainly based on the MR-FDPF method, it is worth investigating whether our approach can be applied to other deterministic channel models, e.g. ray tracing models to extract the fading statistics.

Fig. 1 .Fig. 2 .

 12 Fig. 1. radio coverage prediction provided by the MR-FDPF method at 2.45GHz.

Fig. 3 .

 3 Fig. 3. The obtained local mean path loss, mean path loss and shadowing from the MR-FDPF method.

Fig. 4 .

 4 Fig. 4. CDF comparison between the extracted shadow fading from the MR-FDPF method and normal distribution.

Fig. 4 .

 4 Fig. 4. From the figure, we can see that the extracted shadow fading and the theoretical result are very similar.Moreover, we also compute the standard deviation of the extracted shadow fading, which is 5.87 dB.

Fig. 5 .Fig. 6 .

 56 Fig. 5. The measurement scenario.

Fig. 7 .

 7 Fig. 7. CDF comparison of the extracted shadowing and normal distribution from measurement.

TABLE I

 I

	COMPARISON BETWEEN THE SIMULATION AND MEASUREMENT
		n	The standard deviation σ
	MR-FDPF	1.59	5.87 dB
	Measurement	1.44	7.66 dB

ACKNOWLEDGMENT

The authors would like to thank Dr. Nicolai Czink from FTW of Vienna for performing the measurements and providing the data. This work is funded by the FP7 IPLAN Project.


Qianqian Li

Christian A Rottmair

Robert F Singer

CNT Reinforced Light Metal Composites Produced by Melt Stirring and by High Pressure Die Casting

Keywords: A. Carbon nanotubes (CNTs), A. Metal-matrix composites (MMCs), A. Nano composites, B. Mechanical properties, E. Casting

Light metal matrix composites are of great interest due to their potential for reducing CO 2 emission through lightweight design e.g. in the automotive sector. Carbon nanotubes can be considered as ideal reinforcements, due to their high strength, high aspect ratio and thermomechanic properties. In this research, CNT reinforced light metal composites were produced by melt stirring and by high pressure die casting, which can be both easily scaled up. The light metal composites showed significantly improved mechanical properties already at small CNT contents. The influence of CNT concentration on the composites was also studied.

Introduction

The issue of CO 2 emission has become more and more critical during the last decades. One of the main sources of CO 2 emission is from transportation. One liter of petrol consumption induces 2.34 kg CO 2 and the petrol consumption is directly connected to the average weight of a car. European Commission has proposed to reduce the average CO 2 emission from new cars to 130 g/km by 2012, which means a requirement of 5.6 l petrol / 100 km for a car. To fulfil this regulation, one of the solutions is to reduce the weight of a car. An average weight for a car of less than 800 kg is required. Therefore the use of light metals is becoming important.

Currently lightweight metal alloys are used in relatively small quantities in automotive applications due to their relatively low strength which limits their potential applications. To overcome the limitations, much research has been done on producing light metal matrix composites to improve the mechanical properties by adding different reinforcements.

The reinforcing mechanisms considered to be most important in nano particle reinforced metal matrix composites are load transfer, Orowan mechanism and thermal mismatch. In order to obtain a high strengthening effect,

• Load transfer requires i) a high aspect ratio of the reinforcement ii) a much stronger reinforcement so the load can be efficiently transferred;

• Orowan mechanism requires a small interparticle spacing, i.e. for a certain volume fraction of the reinforcements, smaller particles yield smaller interparticle spacing as the number of particles increases; Additionally rod shaped reinforcements cause more strengthening compared to spherical reinforcements [1];

• Thermal mismatch needs a high difference of the coefficient of thermal expansion (CTE) between the reinforcements and the matrix to create higher dislocation density around the reinforcements.

To sum up, an ideal reinforcement to obtain a good strengthening effect in light metal matrix should have high strength, high aspect ratio, small diameter, a rod shape and a different CTE from light metal matrix, additionally with other good properties such as low density to keep the composites as light as their alloys. Accordingly carbon nanotubes (CNTs) can be considered as ideal reinforcements, due to their: favourable geometrical properties (high aspect ratio, high specific surface area), good physical properties and excellent mechanical properties.

Theoretical simulation [2][3][4][5] and experimental results [6][7][8][9] have indicated their extraordinary strength (about 150 GPa) and Young's modulus (about 1 TPa).

The main reason for introducing CNTs into light metal matrix composites is enhancement of the mechanical properties but also the tailoring of density and hardness. Recent studies also mentioned the modified electrical conductivity, tribological and friction behaviour, and coefficient of thermal expansion. Because of the difficulties in introducing CNTs into metal melts, most work used powder metallurgical (PM) approaches [10][11][12][13][14] including hot pressing [15][16], plasma spray forming [17][18] and hot extrusion [19][20][21][22][23][24]. Generally speaking, the property improvement achieved in CNT/metal composites fabricated by PM techniques has not fulfilled the expectations. This is commonly ascribed to the fact that CNTs are inhomogeneously dispersed and/or agglomerated at the grain boundaries of the metal matrix. A technology capable of dispersing CNTs homogenously in the metal melts is expected to overcome the problems listed above.

If few CNT are present, their reinforcing effect increases with increasing volume fraction. At elevated CNT volume fractions, homogeneous dispersion becomes difficult as CNTs begin to agglomerate. Hence, most research groups [START_REF] Esawi | Carbon nanotube-reinforced aluminium strips[END_REF][START_REF] Deng | Thermal expansion behaviors of aluminum composite reinforced with carbon nanotubes[END_REF][START_REF] Zhou | Fabrication and tribological properties of carbon nanotubes reinforced Al composites prepared by pressureless infiltration technique[END_REF] pointed out that there is an optimum CNT volume fraction which produces maximum enhancement of mechanical properties.

In this work, carbon nanotube reinforced light metal composites were produced by two different methods. One is to produce CNT reinforced Mg AZ91 composites by a two-step process, which has been described in our previous papers [28][START_REF] Li | Magnesium: 8th International Conference on Magnesium Alloys and their Applications[END_REF]. The two-step process includes: step one is a pre-dispersing procedure of CNTs on the Mg, in order to break down big agglomerates; step two is the fabrication of CNT/Mg alloy composite by a melt stirring technique. At step one a block copolymer is used as a dispersing agent to pre-disperse CNTs on Mg alloy chips. The block copolymer was chosen because it has already been proven in previous re-search that it can improve the dispersion of CNT in ethanol [START_REF] Li | [END_REF][31]. In our previous paper, we have already shown that by adding a small amount of CNTs (0.1 wt%), the mechanical properties of the composites have been improved. In this paper, we continue our research and use different amount of CNTs to produce the Mg composites.

Another method is to produce CNT reinforced Al alloy composites by high pressure die casting process. High pressure die casting is an efficient and economic process manufacturing complex shaped components in high quantities. Furthermore, the high-speed casting process provides a high sheer speed to disperse carbon nanotubes homogenously into the Al alloy melt. In this paper, we introduce high pressure die casting to produce CNT/Al alloy composites and investigate the resulting enhancement of mechanical properties.

Experimental

Melt stirring

The block copolymer Disperbyk-2150 (BYK Chemie GmbH) was first dissolved in ethanol in a small beaker. Then multiwall CNTs (MWNTs, mass ratio to the block copolymer 1:1, Bay-tubes® C 150P) were added to the as-prepared solution. This mixture was put at roomtemperature into an ultrasonic bath for 15 minutes. Then it was stirred for 30 minutes at 250 rpm. After adding Mg alloy chips (AZ91 D, ECKA), the suspension was further stirred at 250 rpm inside a fume cupboard to evaporate ethanol and homogenize the mixture. After the mixture was dried, the MWNT coated chips were placed in a cylindrical sample crucible as show in Figure 1 [28]. This crucible was placed into an oven used in our previous paper and heated up to 670 °C under an inert gas atmosphere to avoid oxidation. When the Mg alloy chips were molten, the liquid was mechanically stirred at 350 rpm for 30 minutes to further disperse MWNTs. After stirring, the molten MWNT/Mg composite was poured into a mould. The cooled sample was machined to cylindrical shaped specimens (diameter 5 mm × height 7 mm) for subsequent compression tests. According to the described procedure, AZ91 composites with different contents of MWNTs (0.1 wt%, 0.5 wt% and 1 wt%) were produced. Reference samples were made using exactly the same procedure but from pure AZ91.

After the pre-dispersion step scanning electron microscopy (SEM) was used to observe the microstructure on the surface of CNT coated Mg chips. The resulting samples were tested by compression testing to determine the compression at failure, compressive yield strength and ultimate compressive strength. Tests were conducted at ambient temperature using standard tensile/compression testing equipment (100 kN load cell). Testing was performed at a constant strain rate of 0.01 second -1 . The molten Al alloy (Al239D, AlSi10Mg) was pushed into the die by the piston with a high velocity. Because of the small cross section at the entrance, the melt is accelerated to a high speed (as fast as 250 m/s), which leads to a turbulent flow in the die and could help the dis-persion of MWNTs [32]. After cooling, the samples were cut to dog-bone shapes for tensile testing. Tests were conducted at ambient temperature using standard testing equipment (Instron 4505, 100kN). Testing was performed at a constant cross-head speed of 5 mm/second.

X-ray diffraction (XRD, Philips X' Pert)) was conducted on the pure Al alloy and the MWNT reinforced Al alloy to check for carbide formation in the composites.

Results and Discussion

Melt stirring

For Step 1, the pre-dispersion of CNTs on the Mg alloy chips, SEM analysis was used to study the microstructure of the raw MWNTs and the MWNT coated Mg alloy chips. Figure 3 (a) and (b) are SEM images of the MWNTs as received. From the images, it can be observed clearly that the raw MWNTs are agglomerated in big bundles. Figure 3 (c)-(f) exhibit the SEM images of dispersed MWNTs on Mg alloy chips. Individual MWNTs can be found on the Mg surface after the pre-dispersion process. In the case of 0.1 wt% and 0.5 wt% of MWNTs, no agglomerates can be observed. When increasing the amount of MWNTs, agglomerates of entangled MWNTs with a size of up to 2 µm can be observed (Figure 3f).

Therefore, we can conclude that for a small amount of MWNTs, during step one, MWNTs can be homogenously dispersed on the Mg chips. For the samples with increased amount of MWNTs, agglomerates can still be found. The compression at failure and the ultimate compressive strength increases significantly for the MWNT/AZ91 composites compared to the AZ91 alloy by adding a small amount of CNTs, as discussed in our previous paper [START_REF] Li | Magnesium: 8th International Conference on Magnesium Alloys and their Applications[END_REF]. An increase of the ultimate compressive strength of 20 % has been reached by adding 0.1 wt% MWNTs. An even higher increase was observed for the compression at failure of the MWNT/AZ91 composites, which increased about 40 % compared to the AZ91. The improvement of the mechanical properties of the composites is contributed to excellent properties of carbon nanotubes. Moreover, combining the SEM observation with the results of our mechanical testing, we attribute the improvement of the mechanical properties to the good dispersion of MWNTs in the Mg metal matrix.

Compared to AZ91, adding MWNTs did not change the 0.2 % yield strength as dramatically as other properties. The value of the 0.2 % yield strength for AZ91 is 130 MPa. For the 0.1 wt% MWNT/AZ91 composite, the 0.2 % yield strength is 125 MPa and does not significantly change by increasing the amount of MWNTs, regarding the scattering of the data.

However, because of the limitation of the measured sample geometry, i.e. the top surface of the cylinder sample was not perfectly parallel to the bottom surface, and the short gauge length (7 mm in the case) in the compression test, the 0.2 % yield strength may not be very accurate. 1 % or 2 % yield strength will be measured for comparison in our further study.

When increasing the amount of the MWNTs, there was no obvious effect on the mechanical properties. This can be explained by the fabrication process. It has been observed that with an increasing amount of MWNTs, not all MWNTs could be separated into individual MWNTs during the pre-dispersion process (see Figure 3). Larger agglomerates may be incorporated into the composite. Hence, the mechanical properties remain the same when increasing the amount of MWNTs, as the large agglomerates do not enhance the properties of the resulting composites.

Furthermore, during the experiments it could be observed that by increasing the amount of MWNTs, a larger amount separated to the top or bottom of the melt and did not reside in the melt. Therefore it could be possible that the real amount of MWNTs embedded in the com-

ACCEPTED MANUSCRIPT

posites is similar for the nominal 0.1, 0.5 and 1.0 wt.% MWNT reinforced composites. As a result, the mechanical properties are similar for all the composites. Especially the elongation at fracture has also increased together with the tensile strength, which is different from by using traditional reinforcements such as carbon fibres [33]. In order to check for carbide formation in our composites XRD measurements were conducted on pure Al alloy and the MWNT/Al composites [34]. As can be seen in Figure 6 The lack of carbide formation indicates that the enhanced mechanical properties of the composite can be attributed to the added carbon nanotubes. In the future further characterisation by scanning electron microscopy and transmission electron microscopy would be desirable.

High pressure die casting

Combining the mechanical property improvements and the XRD results, we attribute the enhancement of all mechanical properties to the advanced mechanical properties of the carbon nanotubes.

Conclusions

In this study, CNT reinforced light metal composites were produced successfully by two different methods: melt stirring and high pressure die casting. It demonstrated that CNT reinforced light metal composites can be produced by straight-forward procedures that can be easily scaled up for industrial applications. These composites hold great potential for the use in lightweight construction and through this for the reduction of CO 2 emission in areas like transportation.

By adding only a small amount of CNTs, the mechanical properties of the light metal matrix composites have all been improved. For the melt stirring method, research was carried out following our previous study [START_REF] Zhou | Fabrication and tribological properties of carbon nanotubes reinforced Al composites prepared by pressureless infiltration technique[END_REF]. Compared to the AZ91 Mg alloy, the compression at failure of the MWNT/Mg composites could be improved up to 40 % and the ultimate compressive strength could be improved up to 20 %. For high pressure die casting, tensile testing showed that by adding only 0.05 wt% of MWNTs, the elongation at fracture has increased by 27% compared to pure Al alloy and the tensile strength of CNT/Al composites increased 8% compared to pure Al alloy. XRD results indicate that there is no carbide formation in our CNT/Al composites. We attribute the improvement of the mechanical properties to the excellent properties of CNTs and relatively good dispersion of CNTs in the Mg matrix.

However, no obvious change was found for different amounts of MWNTs. This can be explained by agglomerates still existing in composites of high MWNT content and the difference between embedded amount of MWNTs in the composites and the nominal amount due

Figure 1 .Figure 2 .

 12 Figure 1. Melt stirring equipment [28] Figure 2. High pressure die casting

Figure 3 .

 3 Figure 3. . SEM images of (a), (b) raw MWNTs as received; (c-e) MWNTs dispersed on the surface of the AZ91 chips, where different nominal amounts of MWNTs were added. (f) With increasing weight fraction of MWNTs larger agglomerates can be observed, which are non-existent in the case of the 0.1 wt% and 0.5 wt% MWNT/AZ91 chips.

Figure 4 .

 4 Figure 4. (a) Typical compression stress-strain curves of different MWNT AZ91 alloy composites (0.1 wt.%, 0.5 wt.% and 1 wt.%) and AZ91. Comparison of (b) compression at failure; (c) 0.2% yield strength and (d) ultimate compressive strength averaged over 24 samples.

Figure 5 (

 5 Figure 5(a) shows typical stress-strain curves of the MWNT/Al alloy composites and the pure

Figure 5 .

 5 Figure 5. (a) Typical tensile stress-strain curves of MWNT Al composites (two samples) and pure Al alloy (two samples). Comparison of (b) elongation at fracture and (c) tensile strength as average of obtained 14 samples. (d) dog-bone shaped sample.

  , the XRD results of Al alloy and MWNT/Al composites show complete overlapping. The intensity of the Al and Si peaks shows also no changes. The results indicate that there is no carbide formation in our MWNT/Al composites, which may be due to the low content of CNTs (only 0.05 wt%) in the composites. It can be also explained that in our die casting process, the Al melt has only limited contact time with the CNTs at high temperature (several seconds) before the composite solidifies in the relatively cold die (200 °C).

Figure 6 .

 6 Figure 6. XRD scans of pure Al alloy and MWNT-Al alloy.

Acknowledgements

A financial grant of the Bavarian State Ministry of Science is gratefully acknowledged. We would like to thank Christian Rauber, Natalie Kömpel, Jens Schaufler, Matthias Rübner, Jörg Komma and Sergiu P. Albu for their kind technical assistant. We are indebted to Bayer AG, BYK Chemie GmbH and ECKA Granulate GmbH & Co. KG. for their generous offers of the materials used in this study. We greatly appreciate Dr. Jan Schwerdtfeger for fruitful discussions and suggestions.