V Yaghoubi

S Marelli

B Sudret

T Abrahamsson

Sparse polynomial chaos expansions of frequency response functions using stochastic frequency transformation

Keywords:

Frequency Response Functions (FRFs) are important for assessing the behavior of stochastic linear dynamic systems. For large systems, their evaluations are time-consuming even for a single simulation. In such cases, uncertainty quantification by crude Monte-Carlo (MC) simulation is not feasible. In this paper, we propose the use of sparse adaptive Polynomial Chaos Expansions (PCE) as a surrogate of the full model. To overcome known limitations of PCE when applied to FRF simulation, we propose a frequency transformation strategy that maximizes the similarity between FRFs prior to the calculation of the PCE surrogate. This strategy results in lower-order PCEs for each frequency. Principal component analysis is then employed to reduce the number of random outputs. The proposed approach is applied to two case studies: a simple 2-Degree Of Freedom (DOF) system and a 6-DOF system with 16 random inputs. The accuracy assessment of the results indicates that the proposed approach can predict single FRFs accurately. Besides, it is shown that the first two moments of the FRFs obtained by the PCE converge to the reference results faster than with the MC methods.

Introduction

Interest towards working with large engineering systems is increasing recently, but long simulation time is one of the main limiting factors. Although the development of the computational power of modern computers has been very fast in recent years, increasing model complexity, more precise description of model properties and more detailed representation of the system geometry still result in considerable execution time and memory usage. Model reduction [START_REF] Vakilzadeh | An improved modal approach for model reduction based on input-output relation[END_REF][START_REF] Rahrovani | Modal dominancy analysis based on modal contribution to frequency response function 2-norm[END_REF], efficient simulation [START_REF] Yaghoubi | An efficient simulation method for structures with local nonlinearity[END_REF][START_REF] Avitabile | Efficient techniques for forced response involving linear modal components interconnected by discrete nonlinear connection elements[END_REF][START_REF] Liu | An efficient backward Euler time-integration method for nonlinear dynamic analysis of structures[END_REF] and parallel simulation methods [START_REF] Yaghoubi | A parallel solution method for structural dynamic response analysis[END_REF][START_REF] Tak | High scalable non-overlapping domain decomposition method using a direct method for finite element analysis[END_REF] are different strategies to address this issue.

Consequently, uncertainty propagation in these systems cannot be carried out by classical approaches such as crude Monte-Carlo (MC) simulation. More advanced methods such as stochastic model reduction [START_REF] Amsallem | An online method for interpolating linear parametric reducedorder models[END_REF] or surrogate modeling (Frangos, Marzouk, Willcox, and van Bloemen Waanders, Frangos et al.) are required to replace the computationally expensive model with an approximation that can reproduce the essential features faster. Of interest here are surrogate models. They can be created intrusively or non-intrusively. In intrusive approaches, the equations of the system are modified such that one explicit function relates the stochastic properties of the system responses to the random inputs. The perturbation method [START_REF] Schuëller | Uncertain linear systems in dynamics: Retrospective and recent developments by stochastic approaches[END_REF]) is a classical tool used for this purpose but it is only accurate when the random inputs have small Coefficients Of Variation (COV), say COV<5%. An alternative method is intrusive Polynomial Chaos Expansion (PCE) [START_REF] Ghanem | Stochastic finite elements: a spectral approach[END_REF]. It was first introduced for Gaussian input random variables [START_REF] Wiener | The homogeneous chaos[END_REF] and then extended to the other types of random variables leading to generalized polynomial chaos [START_REF] Xiu | The Wiener-Askey polynomial chaos for stochastic differential equations[END_REF][START_REF] Soize | Physical systems with random uncertainties: chaos representations with arbitrary probability measure[END_REF].

In non-intrusive approaches, already existing deterministic codes are evaluated at several sample points selected over the parameter space. This selection depends on the methods employed to build the surrogate model, namely regression [START_REF] Blatman | An adaptive algorithm to build up sparse polynomial chaos expansions for stochastic finite element analysis[END_REF][START_REF] Berveiller | Stochastic finite element: a non intrusive approach by regression[END_REF] or projection methods [START_REF] Gilli | Uncertainty quantification for criticality problems using non-intrusive and adaptive polynomial chaos techniques[END_REF][START_REF] Knio | A stochastic projection method for fluid flow: I. basic formulation[END_REF]. Kriging [START_REF] Fricker | Probabilistic uncertainty analysis of an FRF of a structure using a Gaussian process emulator[END_REF][START_REF] Jones | Efficient global optimization of expensive black-box functions[END_REF] and non-intrusive PCE (Blatman and Sudret, 2011a) or combination thereof [START_REF] Kersaudy | A new surrogate modeling technique combining Kriging and polynomial chaos expansions-application to uncertainty analysis in computational dosimetry[END_REF][START_REF] Schöbi | Polynomial-chaos-based Kriging[END_REF] are examples of the non-intrusive approaches. The major drawback of PCE methods, both intrusive and non-intrusive, is the large number of unknown coefficients in problems with large parameter spaces, which is referred to as the curse of dimensionality [START_REF] Sudret | Uncertainty propagation and sensitivity analysis in mechanical models -contributions to structural reliability and stochastic spectral methods[END_REF]. Sparse [START_REF] Blatman | Sparse polynomial chaos expansions and adaptive stochastic finite elements using a regression approach[END_REF] and adaptive sparse (Blatman and Sudret, 2011b) PCE have been developed to dramatically reduce the computational cost in this scenario.

To propagate and quantify the uncertainty in a Quantity of Interest (QoI) of a system, its response should be monitored all over the parameter space. This response could be calculated in time, frequency or modal domain. For dynamic systems, the frequency response is important because it provides information over a frequency range with a clear physical interpretation. This is the main reason of the recent focus on Frequency Response Functions (FRF) for uncertainty quantification of dynamic systems and their surrogates [START_REF] Fricker | Probabilistic uncertainty analysis of an FRF of a structure using a Gaussian process emulator[END_REF][START_REF] Goller | An interpolation scheme for the approximation of dynamical systems[END_REF][START_REF] Kundu | A hybrid spectral and metamodeling approach for the stochastic finite element analysis of structural dynamic systems[END_REF][START_REF] Adhikari | Doubly spectral stochastic finite-element method for linear structural dynamics[END_REF][START_REF] Chatterjee | A bi-level approximation tool for the computation of FRFs in stochastic dynamic systems[END_REF].

Several attempts have been made to find a surrogate model for the FRF by using modal properties or random eigenvalue problems. [START_REF] Pichler | A mode-based meta-model for the frequency response functions of uncertain structural systems[END_REF] proposed a mode-based meta-model for the frequency response functions of stochastic structural systems. [START_REF] Yu | Hermite polynomial chaos expansion method for stochastic frequency response estimation considering modal intermixing[END_REF] used Hermite polynomials to solve the random eigenvalue problem and then employed modal assurance criteria (MAC) to detect the phenomenon of modal intermixing. [START_REF] Manan | Prediction of uncertain frequency response function bounds using polynomial chaos expansion[END_REF] used non-intrusive PCE to find the modal properties of a system and predict the bounds for stochastic FRFs. They implemented the method on models with one or two parameters and COV ≤ 2%.

Very few and recent papers addressed the direct implementation of PCE on the frequency responses of systems. [START_REF] Pagnacco | Polynomial chaos for modeling multimodal dynamical systems-investigations on a single degree of freedom system[END_REF] investigated the use of PCE for modeling multimodal dynamic systems using the intrusive approach by studying a single degree of freedom (DOF) system. They showed that the direct use of the polynomial chaos results in some spurious peaks and proposed to use multi-element PCE to model the stochastic frequency response but, to the knowledge of the authors, they did not publish anything on more complex systems yet. Jacquelin et al. (2015b) studied a 2-DOF system to investigate the possibility of direct implementation of PCE for the moments of the FRFs and they also reported the problem of spurious peaks. They showed that the PCE converges slowly on the resonance parts. They accelerate the convergence of the first two statistical moments by using Aitken's method and its generalizations (Jacquelin et al., 2015a).

In general, there are two main difficulties to make the PCE surrogate model directly for the FRFs: their non-smooth behavior over the frequency axis and the frequency shift of the eigenfrequencies due to uncertainties in the parameters. This results in very high-order PCEs even for the FRFs obtained from cases with 1 or 2 DOFs. The main contribution of this work is to propose a method that can solve both problems.

The proposed approach consists of two steps. First, the FRFs are transformed via a stochastic frequency transformation such that their associated eigenfrequencies are aligned in the transformed frequency axis, called scaled frequency. Then, PCE is performed on the scaled frequency axis.

The advantage of this procedure is the fact that after the transformation, the behavior of the FRFs at each scaled frequency is smooth enough to be surrogated with low-order PCEs. However, since PCE is made for each scaled frequency, this approach results in a very large number of random outputs. To solve this issue, an efficient version of principal component analysis is employed.

Moreover, the problem of the curse of dimensionality is resolved here by means of adaptive sparse PCEs.

The outline of the paper is as follows. In Section 2, the required equations for deriving the FRFs of a system are presented. In Section 3, all appropriate mathematics for approximating a model by PCE are presented. The main challenges for building PCEs for FRFs are elaborated and the proposed solutions are presented. In Section 4, the method is applied to two case studies, a simple case and a case with a relatively large number of input parameters.

Frequency response function (FRF)

Consider the spatially-discretized governing second-order equation of motion of a structure as

M q + V q + Kq = f (t) (1)
where for an n-DOF system with n u system inputs and n y system outputs, q(t) ∈ R n is the displacement vector, f (t) is the external load vector which is governed by a Boolean transformation of stimuli vector f (t) = P u u(t); with u(t) ∈ R nu . Real positive-definite symmetric matrices M , V , K ∈ R n×n are mass, damping and stiffness matrices, respectively. The state-space realization of the equation of motion in Eq. ( 1) can be written as

ẋ(t) = Ax(t) + Bu(t), y(t) = Cx(t) + Du(t) (2)
where

A ∈ R 2n×2n , B ∈ R 2n×nu , C ∈ R ny×2n , and D ∈ R ny×nu . x T (t) = [q(t) T , qT (t)] ∈ R 2n is
the state vector, and y(t) ∈ R ny is the system output. A and B are related to mass, damping and stiffness as follows

A =   0 I -M -1 V -M -1 K   , B =   0 M -1 P u   . (3)
The output matrix C, which has application dependent elements, linearly maps the states to the output y and D is the associated direct throughput matrix. The frequency response of the model ( 2) can be written as

H(jω) = C(jωI -A) -1 B + D, (4)
where

H = [H 1 , H 2 , • • • , H nu×ny ] T ∈ C (ny×nu)×1
, ∀ω and j = √ -1. (•) T stands for the transpose of the matrix. It should be mentioned that the eigenvalues of A are the poles of the system. They are complex values and their imaginary parts can be approximated as the frequencies, in rad/s, at which the maximum amplitude occurs.

Methodology

Polynomial chaos expansions

Let M be a computational model with M -dimensional random inputs X={X 1 , X 2 , ..., X M } T and a scalar output Y . Further, let us denote the joint probability distribution function (PDF) of the random inputs by f X (x) defined in the probability space (Ω,F , P).

Assume that the system response Y = M(X) is a second-order random variable, i.e. E Y 2 < +∞ and therefore it belongs to the Hilbert space H = L 2 f X (R M , R) of f X -square integrable functions of X with respect to the inner product:

E [ψ(X)φ(X)] = D X ψ(x)φ(x)f X (x)dx (5)
where D X is the support of X. Further assume that the input variables are independent, i.e.

f X (x) = M i=1 f Xi (x i ).
Then the generalized polynomial chaos representation of Y reads [START_REF] Xiu | The Wiener-Askey polynomial chaos for stochastic differential equations[END_REF]:

Y = α∈N M ũα ψ α (X) (6)
in which ũα is a set of unknown deterministic coefficients, α = (α 1 , α 2 , ..., α M ) is a multi-index set which indicates the polynomial degree of ψ α (X) in each of the M input variables. ψ α s are multivariate orthonormal polynomials with respect to the joint PDF f X (x), i.e. :

E [ψ α (X)ψ β (X)] = D X ψ α (x)ψ β (x)f X (x)dx = δ αβ (7)
where δ αβ is the Kronecker delta. Since the input variables are assumed to be independent, these multivariate polynomials can be constructed by a tensorization of univariate orthonormal polynomials with respect to the marginal PDFs, i.e.

ψ α (X) = M i=1 ψ (i)
αi (X i ). For instance, if the inputs are standard normal or uniform variables, the corresponding univariate polynomials are Hermite or Legendre polynomials, respectively.

In practice, the infinite series in Eq. ( 6) has to be truncated. Given a maximum polynomial degree p, the standard truncation scheme includes all polynomials corresponding to the set A M,p = {α ∈ N M : |α| ≤ p}, where |α| = M i=1 α i is the total degree of polynomial ψ α . The cardinality of the set A M,p = M +p p = P increases rapidly by increasing the number of parameters M and the order of polynomials p. However, it can be controlled with suitable truncation strategies such as q-norm hyperbolic truncation [START_REF] Blatman | An adaptive algorithm to build up sparse polynomial chaos expansions for stochastic finite element analysis[END_REF], that drastically reduce the number of unknowns when M is large.

The estimation of the vector of coefficient ũα can be done non-intrusively by projection [START_REF] Ghanem | Stochastic seismic soil-structure interaction using the homogeneous chaos expansion[END_REF][START_REF] Ghiocel | Stochastic finite element analysis of seismic soil-structure interaction[END_REF] or least square regression methods [START_REF] Blatman | An adaptive algorithm to build up sparse polynomial chaos expansions for stochastic finite element analysis[END_REF][START_REF] Berveiller | Stochastic finite element: a non intrusive approach by regression[END_REF]. The latter is based on minimizing the truncation error via least square as follows:

Y = M(X) = α∈A M,p ũα ψ α (X) + ≡ Ũ T Ψ(X) + (8)
This can be formulated as

Û = arg min E Ũ T Ψ(X) -M(X) 2 . ( 9
)
Let X = {x (1) , x (2) , ..., x (N ED ) } and Y = {y (1) = M(x (1) ), y (2) = M(x (2) ), ..., y (N ED ) = M(x (N ED ) )} be an experimental design with N ED space-filling samples of X and the corresponding system responses, respectively. Then, the minimization problem (9) admits a closed form solution

Û = (Ψ T Ψ) -1 Ψ T Y, (10)
in which Ψ is the matrix containing the evaluations of the Hilbertian bases, that is Ψ ij = ψ αj (x (i) ), i = 1, 2, ..., N ED , j = 1, 2, ..., P .

The accuracy of PCE will be improved by reducing the effect of over-fitting in least square regression. This can be done by using sparse adaptive regression algorithms proposed in [START_REF] Hastie | Forward stagewise regression and the monotone lasso[END_REF][START_REF] Efron | Least angle regression[END_REF]. In particular, the Least Angle Regression (LAR) algorithm has been demonstrated to be effective in the context of PCE by Blatman and Sudret (2011a).

Vector-valued response

In the case of vector-valued response, i.e. Y ∈ R N , N > 1, the presented approach may be applied componentwise. This can make the algorithm computationally cumbersome for models with large number of random outputs. To decrease the computational cost, one can extract the main statistical features of the vector random response by Principal Component Analysis (PCA). The concept has been adapted to the context of PCE by [START_REF] Blatman | Sparse polynomial chaos expansions of vector-valued response quantities[END_REF].

To perform sample-based PCA, let us rewrite the Y as the combination of its mean Ȳ and covariance matrix as follows:

Y = Ȳ + N i=1 u i v T i (11)
where the v i 's are the eigenvectors of the covariance matrix:

COV (Y) = E (Y -Ȳ) T (Y -Ȳ) = [v 1 , ..., v N ]      l 1 . . . 0 . . . 0 . . . l N           v T 1 . . . v T N      T ( 12
)
and the u i 's are vectors such that

u i = (Y -Ȳ)v i . ( 13
)
One can approximate Y by the N -term truncation:

Y = Ȳ + N i=1 u i v T i , N N. ( 14
)
Since Ȳ and v T i are the mean and the eigenvectors of the system responses, respectively, they are independent of the realization. Therefore, PCE can be applied directly on the N N auxiliary variables u i . Besides, acknowledging the fact that the PCA is an invertible transform, the original output can be retrieved directly from Eq. ( 14) for every new prediction of u.

Vector-valued data with extremely large output size

Assume Y ∈ R N ED ×N , has an extremely large N and N N ED . Then COV (Y) ∈ R N ×N is exceptionally large and solving the eigenvalue problem numerically may be unfeasible. To address this issue, the following well-known theorem and the associated corollary is presented.

Theorem 1. (Singular value decomposition) Let A ∈ R n×N , n < N and rank(A) = n then there exist two orthogonal matrices, U and V and two diagonal matrices S and Σ such that

A = U SV T = U   Σ 0 0 0   V T in which U ∈ R n×n , S ∈ R n×N , Σ ∈ R n×n and V ∈ R N ×N . Furthermore, this
decomposition can be written as eigenvalue decomposition as AA T U = U Σ and A T AV = V S.

Corollary 1. The nonzero eigenvalues of A T A and AA T are equal. Furthermore, U and V are related to each other by

U = AV S -1 (15)
The proof of the theorem can be found in any matrix analysis book, e.g. [START_REF] Laub | Matrix Analysis For Scientists And Engineers[END_REF] and the corollary directly follows from the theorem.

Therefore, instead of the eigenvalue calculation of (Y -Ȳ) T (Y -Ȳ) ∈ R N ×N , which may be an

extremely large matrix, one can consider (Y -Ȳ)(Y -Ȳ) T ∈ R N ED ×N ED , which is much smaller.
The associated eigenvectors can be transformed to the ones in Eq. ( 14) through Eq. ( 15).

Stochastic frequency transformation

In this section, the method of stochastic frequency transformation is developed to address the challenge of frequency shift at eigenfrequencies due to uncertainty in the parameters. The idea is basically to apply a transformation to the system responses to maximize their similarity before building PCE, as first proposed by [START_REF] Mai | Polynomial chaos expansions for damped oscillators[END_REF]. Here, the technique is extended and adopted into the frequency domain to obtain PCEs of the FRFs.

To this end, the following algorithm is proposed. First, an experimental design X and the corresponding model responses Y are evaluated. Each system response will be called a trajectory in the remainder of this paper. Let the frequency range of interest be discretized to n ω equidistant

frequencies

Ω d = [ω 1 , ω 2 , ..., ω nω ].
Then, the required system responses are matrices H(Ω d ) ∈

C nuny×nω and F ∈ R nuny×n sf . The matrix H(Ω d ) is obtained by evaluating Eq. ( 4) at frequencies Ω d . The matrix F consists of all the peak and valley frequencies of the system's input-output relations for one system realization, as follows:

F =         ω 1 ω p1 ω 1 m1 ω p2 . . . ω pn p-1 ω 1 mn p -1 ω pn p ω nω ω 1 ω p1 ω 2 m1 ω p2 . . . ω pn p-1 ω 2 mn p -1 ω pn p ω nω . . . . . . . . . . . . . . . . . . . . . . . . . . . ω 1 ω p1 ω nu×ny m1 ω p2 . . . ω pn p -1 ω nu×ny mn p -1 ω pn p ω nω         =         F 1 F 2 . . . F nu×ny         (16)
in which n p is the number of eigenvalues of the system. Furthermore, {ω pi i = 1, 2, ..., n p } are the resonant frequencies and {ω l mi i = 1, 2, ..., n p -1, l = 1, 2, ..., n u × n y } are frequencies between each two consecutive resonant frequencies at which the minimum amplitude occurs. Throughout the paper, these important frequencies, shown by red asterisks in Figure 1 for a typical frequency response, will be referred to as selected frequencies. Their number n sf is assumed to be constant across different realizations of the system inputs. {F i , i = 1, 2, • • • , n u ×n y } includes all the selected frequencies for the i th input-output relation.

For the next step of the algorithm, let x (ref ) be selected randomly among the sample points in the ED to have its associated trajectory as the reference, i.e. :

H ref = H(x (ref ) , Ω d ), F ref = F (x (ref ) ; ω).
Then, the other trajectories are transformed in the frequency axis so as to have the peaks and valleys as close to the corresponding locations in the reference trajectory as possible i.e. :

T k i = T (k) i (ω, ν (k) i ) = {ν (k) i = f (ω)|F i (x (k) ; ν (k) ) = F ref i } ( 17
) where i = 1, 2, • • • , n u × n y , k = 1, 2, • • • , N ED
and ν is the transformed frequency axis called scaled frequency. The transform T k i consists of a continuous piecewise linear transform of the intervals between the identified selected frequencies that align them to the corresponding ones of the reference trajectory as follows,

T k i : ν (k) i,l = a (k) ω l + b (k) F (k) i (j) ≤ ω l ≤ F (k) i (j + 1) (18)
where

a (k) = F ref i (j) -F ref i (j + 1) F (k) i (j) -F (k) i (j + 1) , b (k) = F ref i (j)F (k) i (j + 1) -F ref i (j + 1)F (k) i (j + 1) F (k) i (j) -F (k) i (j + 1) , j = 1, 2, • • • , n sf -1 and l = 1, 2, • • • , n ω .
This transformation results in the FRFs which are similar to the reference one in the scaled frequency domain:

Hi (x (k) , N (k) i ) = H i (x (k) , Ω d ) • T k i (19)
where the set

N (k) i = {ν (k) i,1 , ν (k) i,2 , • • • , ν (k)
i,nω } consists of the discretized scaled frequencies which are non-equidistantly spread over the frequency range of interest.

Figures 2a and2b illustrate the FRFs of a 2-DOF system versus frequency and scaled frequency,

respectively. An example of such a transform used for transforming the FRFs of a 2-DOF system is presented in Figure 3.

One should notice that since N (k) i contains the non-equidistant scaled frequencies, a final interpolation is required to obtain a common discretized scaled frequency

N ref = {ν ref 1 , ν ref 2 , • • • , ν ref nω }
between the reference and all other trajectories. To reduce interpolation error in the system response, small frequency steps should be selected. The proposed approach for preprocessing the FRFs is summarized in Algorithm 1.

Algorithm 1 Data preprocessing: continuous piecewise-linear transformation k) ; ω)] T using Eq. ( 16) k) ; Ω d )] T using Eq. ( 4)

1: Input: X = {x (1) , x (2) , ..., x (N ED ) } 2: H ref =H(x (r) , Ω d ), F ref = F (x (r) ; ω), for a random r ∈ [1, ..., N ED ] 3: for k = 1 to N ED do 4: F (x (k) ; ω)=[F 1 (x (k) ; ω)F 2 (x (k) ; ω), • • • , F nu×ny (x (
5: H(x (k) ; Ω d )=[H 1 (x (k) ; Ω d ), H 2 (x (k) ; Ω d ), • • • , H nu×ny (x (
6: for i = 1 to n u × n y do 7:
Evaluate T k i using Eq. ( 18)

8: Hi (x (k) , N (k) i )=H i (x (k) , Ω d ) • T k i 9: Hi (x (k) , N ref )=interpolate( Hi (x (k) , N (k) i ), N (k) i , N ref ) 10:
end for 11:

H(x (k) ; N ref )=[ H1 (x (k) ; N ref ), H2 (x (k) ; N ref ), • • • , Hnu×ny (x (k) ; N ref )] T 12: end for 13: F={vect(F (x (1) ; ω)), vect(F (x (2) ; ω)), ..., vect(F (x (N ED ) ; ω))} T 14: H(N ref )= {vec( H(x (1) , N ref )), vec( H(x (2) , N ref )), ..., vec( H(x (N ED ) , N ref ))} T , 15: Output: F, G R =real( H(N ref )), G I =imag( H(N ref ))

Polynomial chaos representation

The non-smooth behavior of the FRFs makes their surrogation by polynomials a problematic task.

To solve this issue, one PCE could be calculated for each scaled frequency. This means that to compute PCEs for the FRFs, two sets of PCEs are required. The first set is to predict the selected frequencies, collected in the matrix F (16), which are required for performing the stochastic transformation as explained in Section 3.3. This matrix includes eigenfrequencies of the system, therefore by obtaining this set of PCEs, the problem of random eigenvalue calculation is solved by the use of PCE as a byproduct. This problem has been addressed in some recent works, e.g. [START_REF] Pichler | A meta-modeling technique for the natural frequencies based on the approximation of the characteristic polynomial[END_REF]. Since the number of random outputs for this set is not very large, PCE can be applied to each of the selected frequencies separately, i.e.

for i = 1, 2, • • • , n u × n y and j = 1, 2, • • • × n sf Fi (j) = α∈A M,p f i α (j)ψ α (x). ( 20
)
The second set of PCE is for the system response at each individual scaled frequency. To this nω×nu×ny) , defined in Algorithm 1, be a matrix of trajectories at scaled frequencies. Since the FRFs have complex-valued responses, whereas the PCEs are defined for real-valued functions only1 , separate PCEs need to be performed for real and imaginary parts of the FRFs. Therefore, the matrix:

end, let H(N ref ) ∈ C N ED ×(
G = {G R , G I } = {real( H(N ref )), imag( H(N ref ))} ∈ R N ED ×(2×nω×nu×ny)
is the response matrix for which the PCE should be built.

The number of random outputs for this set, N = 2 × n ω × n u × n y , can be extremely large. As discussed in Section 3.2, the PCEs are therefore applied directly to the principal components of G,

yielding: ĜR = ḠR + N j=1 α∈A M,p (u R α ψ α (x)) j v R T j , (21)
ĜI = ḠI + N j=1 α∈A M,p (u I α ψ α (x)) j v I T j , (22)
where u R α and u I α are respectively the vectors of coefficients of the PCEs made for the real and imaginary parts of the FRFs.

Surrogate response prediction

To predict the surrogate model response at a new sample point x (0) , several steps need to be taken to transform the PCE predictors in Eqs. ( 21) and ( 22) from the scaled frequency axis ν to the original frequency axis ω. The matrices ĜR and ĜI are obtained by evaluating the second set of PCEs in Eqs. ( 21) and ( 22), respectively. Then, the FRFs at the scaled frequencies can be obtained at the new sample point by the inverse vectorization of Ĥ(N ref ) = ĜR + j ĜI where j = √ -1. To obtain the FRF at the original frequency ω the following transformation is used,

Ĥi (x (0) , Ω (0) i ) = Ĥi (x (0) , N ref ) • (T 0 i ) -1 , i = 1, 2, • • • , n u × n y (23)
where T 0 i is obtained by evaluating Eq. ( 18) at F (0) = F (x (0) ; ω) which is the matrix of selected frequencies at the new sample point x (0) evaluated by Eq. ( 20). Besides, Ω (0) i is a set of discretized frequencies which are non-equidistantly spread over the frequency rage of interest. In order to provide the frequency response at the desired frequencies Ω d , interpolation is inevitable. The algorithm for predicting the system response at a new sample point is briefly presented in Algorithm 2.

Examples

Introduction

In this section, the proposed method will be applied to two case studies. The first one is a simple 2-DOF system to illustrate how the method works. The second one is a system with a relatively large (16-dimensional), parameter space.

To assess the accuracy of the surrogate models quantitatively, two measures are defined. The first one aims at measuring the relative difference between the selected frequencies:

∆ω = |ω ex -ω approx | ω ex × 100. ( 24
)
Algorithm 2 Predicting system responses ĜR = ĜR (x (0) , N ref ) using Eq. ( 21).

1: Input: x (0) = x (l) , l = 1,
3: ĜI = ĜI (x (0) , N ref ) using Eq. ( 22). Evaluate F0 i = Fi (x (0) ; ω) using Eq. ( 20)

8:

Evaluate T 0 i using Eq. ( 18)

9:
Ĥi (x (0) , Ω

i )= Ĥi (x (0) , N ref ) • (T 0 i ) -1 10: Ĥi (x (0) , Ω d )=interpolate( Ĥi (x (0) , Ω (0) i ), Ω (0)
i , Ω d ) 11: end for

12: Output: Ĥ(Ω d )={ Ĥ1 (x (0) , Ω d ), Ĥ2 (x (0) , Ω d ), • • • , Ĥnu×ny (x (0) , Ω d )}
The second one aims at measuring the relative difference between vectorial data, such as one FRF or the mean and standard deviation of several FRFs. This error is based on the root mean square (rms) error of the logarithm of the vectors, defined as:

Error(•) = rms(log(|(•) ex |) -log(|(•) approx )|) rms(log(|(•) ex |)) × 100, (25)
in which (•) is the vector of interest. (•) ex and (•) approx represent results obtained by running the true and surrogate models, respectively. For the mean and standard deviation of the data, the reference results are obtained by evaluating the true model at 10,000 MC samples and the approximations are calculated by the PCE surrogate at the same 10,000 points.

4.2 Simple 2-DOF system (Jacquelin et al., 2015b) As the first example, the simple 2-DOF system shown in Figure 4 is selected to highlight the steps of the proposed method. In this system, stiffness is assumed to be uncertain

k = k(1 + δ k ξ) ( 26
)
where ξ is a standard normal random variable. Other properties of the system are listed in Table 1. The system has one input force f at mass 1, two physical outputs q 1 and q 2 and thus, two FRFs.

The FRFs of the system are obtained in the range of 10 to 35 Hz with a frequency step of 0.01 Hz, as shown in Figure 1. The selected frequencies are also shown in the figure with red asterisks.

m m k c k c q 2 q 1 f Figure 4: 2-DOF system Table 1: 2-DOF system's charactristic Characteristics m (kg) k (Nm -1 ) c (Nm -1 s -1 ) δ k Value 1 15000 1 5%
40 points are sampled in the parameter space using Latin Hypercube Sampling (LHS) to form an Experimental Design (ED) X and the model is evaluated at these points to find the system responses of interest, namely the FRFs, H, and the selected frequencies F .

Figure 2a shows the FRFs of the system evaluated at X . To find the transformed FRFs, i.e.

FRFs at scaled frequencies ν, one trajectory was selected randomly as the reference and the others were scaled such that their peaks and valleys were at the same scaled frequencies as that of the reference trajectory. The transformed FRFs and the corresponding continuous piecewise-linear transformations are respectively shown in Figure 2b and Figure 3.

The next step is to find a suitable basis and the associated coefficients for the PCE. In this case, since the random variable is Gaussian, the basis of the polynomial chaos consists of Hermite polynomials. The LARS algorithm (Blatman and Sudret, 2011b) is employed here to calculated a sparse PCE with adaptive degree.

The first set of expansion consists of 10 PCEs to surrogate the selected frequencies F , shown in Figure 1 by red asterisks. As the second set of expansions, PCE is made for the dominant components of G as explained in Section 3.2.1. To do so, the N largest principal components are selected such that the sum of their associated eigenvalues amounts to 99% of the sum of all the eigenvalues, i.e.

N i=1 λ i = 0.99 N i=1 λ i in which λ i 's are the eigenvalues of the covariance matrix of either G R or G I . By this truncation, the number of random outputs is reduced from 2501 × 2 × 2 to 6 components, namely 3 components for the real part and 3 components for the imaginary part.

The spectra of the eigenvalues of the covariance matrix of the G R and G I are displayed in Figure 5. All the PCEs used to make this surrogate model, including those for the selected frequencies and those for the dominant components have orders between 3 and 8. The efficiency of the proposed approach is assessed by comparing the prediction accuracy of the surrogate model on a large reference validation set (10,000 samples calculated with the full model).

# Eigenvalues

PCE estimate of the mean and standard deviation of the surrogate model are compared to their MC estimators on ED of increasing size. The resulting convergence curves are given in Figure 6.

They indicate that the PCE converges faster to the reference results for the mean and standard deviation. Their estimates are approximately two and one order of magnitude more accurate than those from the MC estimators. In addition, one can conclude that 40 points are enough for the ED in this example, since for larger sizes the accuracy does not improve significantly.

It is worth mentioning that the COV of the parameters and the level of damping are among the criteria that can affect the size of the ED. Therefore, larger COV and lower levels of damping are not obstacles for the proposed method provided a sufficiently large ED is used. The 10,000 model evaluations used to produce the convergence curves in Figure 6 are also used to provide a detailed validation of the performance of the PCE surrogates on various quantities of interests. The results are presented in the following figures. In Figure 7, the selected frequencies obtained by the PCE are shown versus the true ones. The results show that the PCE model accurately predicts the selected frequencies. To assess this accuracy quantitatively, the error defined in Eq. ( 24) is used for each selected frequency and their histograms are presented. The excellent agreement shown in Figure 7 is confirmed by the error distribution in Figure 8. To show the feasibility of the proposed method to estimate the statistics of the FRFs, the results obtained here are compared to their counterparts in two of the most recent works available in the literature. The first study (Jacquelin et al., 2015a) directly uses high-order PCE for estimating the first two moments of the FRF, whereas the second method (Jacquelin et al., 2015b) proposes to use Aitken's transformation in conjunction with PCEs. Both methods use PCEs of order 50 and tend to produce spurious peaks around the resonance region. The use of Aitken's transformation slightly improves convergence. Their results for the mean and standard deviation are shown in Figure 13 and Figure 14 respectively. For comparison, the results from our approach in Figure 11a and Figure 12a are reproduced in Figure 13c and Figure 14c respectively, with a scaling similar to the other panels. They indicate that the stochastic frequency transformation approach proposed here significantly improves the estimation accuracy of the PCE surrogate, as no spurious peaks are visible in this case.

# ED points

6-DOF system: large parameter space

The second example is chosen to illustrate the application of the proposed method to a problem with a relatively large parameter space. The system, shown in Figure 17, consists of 10 springs and 6 masses which are modeled by random variables with lognormal distributions. Their mean values are listed in Table 2. The uncertainty on the springs (resp. masses) has a COV = 10% (resp.

COV = 5%). The damping matrix is V = 0.01 M , where M is the matrix of the mean value of the system masses. The system has one input force at mass 6 and 6 system outputs, one for each mass. The FRF of the system is evaluated at a frequency range from 1 to 25 rad/s with the step of 0.01π rad/s.

In this example, the ED consists of 1,000 points sampled from the parameter space using LHS.

The marginal distributions of the input vectors X consists of lognormal distributions. Therefore, the chosen PCE basis consists of Hermite polynomials on the reduced variable Z = ln(X). Eq. ( 6) thus, can be written as

Y = M(X) = α∈A M,p
ũα ψ α (ln(X)).

The LARS algorithm has been employed to build sparse PCEs with adaptive degree for both the selected frequencies and the principal components of the scaled FRF. For the second set of PCEs, PCA has been performed and the dominant components are selected such that N i=1 λ i = 0.999 N i=1 λ i . This truncation reduced the number of random outputs from 761 × 6 × 2 to 82 components. Since the dimension of the input parameter space is large, to reduce the unknown coefficients of the PCEs and avoid the curse of dimensionality, a hyperbolic truncation with q-norm of 0.7 was used before the LARS algorithm. Besides, only polynomials up to rank 2 were selected here (i.e. polynomials that depend at most on 2 of the 16 parameters). It should be mentioned that all the PCEs used for the surrogate model have eventually maximum degrees less than 10.

To start building PCEs, the basis should be selected. Since X ∼ LN (λ, ξ) to find suitable basis, one isoprobabilistic transform is employed here as follows, Z = ln(X) ∼ N (λ, ξ), therefore, the Hermite polynomials can be used as basis for Z, and Eq. ( 6) can be written as

Y = M(X) = α∈A M,p ūα ψ α (ln(X)).
In order to select a proper ED, the convergence of the mean and standard deviation of the FRFs w.r.t the reference ones has been studied by enlarging the ED. They are shown in Figures 18 and19, respectively, for the mean and standard deviation at two outputs. The results for the other outputs are presented in A. They indicate that the mean converges faster than when using Monte-Carlo simulation but the standard deviations converges slower due to the error introduced by interpolation and PCA. Besides, it can be inferred that 1,000 points are enough for the ED.Therefore, 1,000 sample points are selected from the parameter space using LHS. In order to assess the accuracy of the surrogate model in estimating various quantities of interests, the same 10,000 points used as the reference validation set to study the convergence are used here. Figure 20 illustrates some of the predicted selected frequencies (only eigenfrequencies) versus the true ones so that the accuracy of the surrogate model in this step can be inferred. While the overall accuracy is very good for all frequencies, it tends to degrade somewhat at higher frequencies.

m 1 m 2 m 3 m 6 m 5 m 4 k 1 c 1 k 2 c 2 k 3 c 3 k 4 c 4 k 5 c 5 k 6 c 6 k 7 c 7 k 8 c 8 k 9 c 9
Besides, at all the validation points the FRFs are calculated by both the true model and the surrogate model. Plots of the individual FRFs are reported in B. In order to assess the error one of the true model in the root-mean-square sense. This error is evaluated using Eq. ( 25 In general, it can be concluded that the proposed method predicts the eigenfrequencies and FRFs with excellent precision, although the detailed peaks are generally captured less accurately. This is the main reason of slow convergence of the standard deviation process shown in Figure 19.

One of the sources for this error is the interpolation. One can reduce it by reducing the frequency step. In addition, close to the resonance frequency, the amplitude is driven by damping ratio, therefore another reason for this error is the low damping ratio. In this example, damping about 0.01-0.1% of critical damping was assumed for each mode.

Conclusions

A novel method to build a surrogate model directly for the FRFs of stochastic linear dynamic systems based on sparse PCE has been proposed. To this end, there were two major challenges which have been addressed in this paper: the frequency shifts in the selected frequencies of the FRFs, i.e. peaks and valleys, due to the uncertainty in the parameters of the system and the non-smooth behavior of the FRFs. These can lead to very high-order PCEs even for the FRFs obtained from cases with 1 or 2 DOFs. We thus propose a stochastic frequency-transformation as a preprocessing step before building PCEs. This transformation scales the FRFs in the frequency horizon so that their selected frequencies become aligned. Although this preprocessing step results in one extra set of PCEs, they do not require any additional full model evaluations. After the transformation, FRFs are very similar and low-order PCEs can be built for each frequency. This leads to an enormously large number of random outputs. An efficient implementation of principal component analysis has been used to alleviate this issue. Moreover, the problem of curse of dimensionality of PCEs in cases with large parameter space was resolved by employing the LARS algorithm to build sparse PCEs together with adaptive degree.

Successfully applied to two case studies, the proposed method shows its capability of accurately 1) predicting individual FRFs, 2) estimating the mean and standard deviation of the FRFs, and 3) estimating the eigenfrequencies of the system and their statistics. Some errors have been observed

at the peaks, likely caused by a combination of interpolation error and very small damping ratio.

Further investigation needs to be made to consider the effect of damping on the accuracy of the FRF prediction.

A Convergence analysis

In this appendix, the results of the convergence analysis of the statistics of the FRFs of the 6-DOF system is presented. The reference results are obtained by running the true model at 10,000 points sampled randomly from the parameter space.

A.2 Standard deviation of the FRFs C Moments of the FRFs of the 6-DOF system

In this appendix, the statistics of the FRFs are compared. They are obtained by evaluating the surrogate model and the true model at 10,000 points sampled from the parameter space.

C.1 Mean of the FRFs

D Individual FRFs comparison of the 6-DOF system

In this appendix, individual FRFs are compared for two particular cases. One FRF that has average error comparing to the true model and one FRF with the maximum error among 10,000 realizations.

D.1 Typical FRFs

Frequency (Hz)

Frequency

  Figure 1: FRFs of the 2-DOF system presented in Figure 4. The selected frequencies F and the associated notations are illustrated with asterisks ( * ).

Frequency

  FRFs after frequency transformation.

Figure 2 :FrequencyFigure 3 :

 23 Figure 2: Several realizations of the FRFs of the 2-DOF system at first system output.

4:

  Ĥ(x (0) , N ref ) = ĜR +j ĜI 5: Construct Ĥ(x (0) , N ref ) from Ĥ(x(0) , N ref ) by inverse vectorization operation 6: for i = 1 to n u × n y do 7:

Figure 5 :

 5 Figure 5: Spectrum of the eigenvalues of the covariance matrix for the 2-DOF system; evaluated for both real part G R and imaginary part G I

Figure 6 :

 6 Figure 6: Convergence plot of the statistics of the FRFs obtained by the PCE ( * ) and the true model (×) with increasing ED size. The reference results were obtained by 10,000 MC simulations of the true model.

Figure 7 :Figure 8 :

 78 Figure 7: Selected frequencies predicted by the surrogate model versus by the true model. upper row: eigenfrequencies, lower row: frequencies where minimum amplitude occurs, see matrix (16) for notations.

Figure 9 :Figure 10 :

 910 Figure 9: All the FRFs obtained by evaluating the true and the surrogate model at 10,000 MC samples.

FrequencyFigure 11 :FrequencyFigure 12 :

 1112 Figure 11: Mean of the FRFs evaluated over 10,000 sample points, by the true model (red) and by the surrogate model (black).

Frequency

  Direct use of PCE with order 50. Shown in Jacquelin et al. (2015a,b) PCE with order 50 and Aitken's transformation. Proposed in Jacquelin et al. (2015b) The proposed method in this work. Sparse PCEs with order less than 8 and stochastic frequency transformation

Figure 13 :Frequency

 13 Figure 13: Average of the FRFs at the first mass, evaluated by PCE: comparison between the proposed method and state of the art methods. The red line is for the reference result and the black line is for the surrogate models

Figure 14 :FrequencyFigure 15 :FRFFigure 16 :

 141516 Figure 14: Standard deviation of the FRFs at the first mass, evaluated by PCE: comparison between the proposed method and state of the art methods. The red line is for the reference result and the black line is for the surrogate models

Figure 18 :Figure 19 :

 1819 Figure 17: The 6-DOF system

Figure 20 :Figure 21 :

 2021 Figure 20: The eigenfrequencies predicted by the surrogate model versus obtained by the true model, evaluated at 10,000 MC samples.

Frequency

Figure 22 :

 22 Figure 22: Typical FRFs predicted by the surrogate model at 2 outputs, evaluated by the true model (red line) and the surrogate model (black line).

FrequencyFigure 23 :

 23 Figure 23: The worst FRFs (out of 10,000) predicted by the surrogate model at 2 outputs, evaluated by the true model (red line) and the surrogate model (black line).

FrequencyFigure 24 :FrequencyFigure 25 :

 2425 Figure 24: Mean of the FRF of the 6-DOF system at 2 outputs, evaluated at 10,000 MC sample points by the true model (red line) and the surrogate model (black line).

Figure 27 :Frequency

 27 Figure 27: Convergence plot of the standard deviation of the FRFs at 4 outputs obtained by the PCE (black * ) and the true model (red ×) by enlarging the experimental design. The reference results were obtained by 10,000 Monte-Carlo simulation of the true model

Figure 28 :Frequency

 28 Figure 28: Envelope of the FRF of the 6-DOF system, evaluated at 10,000 MC sample points by the true model and the surrogate model. The black lines are the averaged FRFs.

Figure 29 :

 29 Figure 29: Envelope of the FRF of the 6-DOF system, evaluated at 10,000 MC sample points by the true model and the surrogate model. The black lines are the averaged FRFs.

FrequencyFigure 30 :Figure 31 :

 3031 Figure 30: Mean of the FRF of the 6-DOF system at 4 outputs, evaluated at 10,000 MC sample points by the true model (red line) and the surrogate model (black line).

Figure 32 :Figure 33 :

 3233 Figure 32: Typical FRFs predicted by the surrogate model at 4 outputs, evaluated by the true model (red line) and the surrogate model (black line).

  2, ..., N ED and H ref , F ref , Ḡ and N ref

	2:

Limited literature is available on the use of PCE for complex-valued functions, see e.g.[START_REF] Soize | Physical systems with random uncertainties: chaos representations with arbitrary probability measure[END_REF].

Acknowledgment

The first author would like to express his gratitude to the ETH Zurich for its kind host during this work.


Baldwin Dumortier

Emmanuel Vincent

Madalina Deaconu

RECURSIVE BAYESIAN ESTIMATION OF THE ACOUSTIC NOISE EMITTED BY WIND FARMS

Keywords: Audio source separation, nonlinear Kalman filtering, uncertainties, wind turbine noise

Wind turbine noise is often annoying for humans living in close proximity to a wind farm. Reliably estimating the intensity of wind turbine noise is a necessary step towards quantifying and reducing annoyance, but it is challenging because of the overlap with background noise sources. Current approaches involve measurements with on/off turbine cycles and acoustic simulations, which are expensive and unreliable. This raises the problem of separating the noise of wind turbines from that of background noise sources and coping with the uncertainties associated with the source separation output.

In this paper we propose to assist a black-box source separation system with a model of wind turbine noise emission and propagation in a recursive Bayesian estimation framework. We validate our approach on real data with simulated uncertainties using different nonlinear Kalman filters.

INTRODUCTION

Wind energy is one of the most used sustainable energies in the world. The stochastic nature of wind raises numerous research problems including control [START_REF] Guo | An approximate wind turbine control system model for wind farm power control[END_REF], short-term wind prediction [START_REF] Holland | Location robust estimation of predictive Weibull parameters in short-term wind speed forecasting[END_REF][START_REF] Kuh | Applications of complex augmented kernels to wind profile prediction[END_REF], resource assessment [START_REF] Bahlke | Stochastic load scheduling for risk-limiting economic dispatch in smart microgrids[END_REF] or power curve estimation [START_REF] Lydia | Advanced algorithms for wind turbine power curve modeling[END_REF][START_REF] Liu | An improved interpolation method for wind power curves[END_REF]. Among the limiting factors of the electric production, acoustic annoyance has been somewhat left behind in the prediction and control field. Yet, complaints from inhabitants can lead to severe penalties for the owners that must ensure the compliance with the applicable law of the country where the farm is set. In all countries acoustic annoyance is quantified by either the total noise level resulting from all acoustic sources which is called ambient noise, or by the acoustic emergence [START_REF] Schild | State of the art and new perspectives for the development of noise regulation of wind farms[END_REF]. The acoustic emergence is the difference between the ambient noise level and the background noise level due to other noise sources besides the turbines. It can be seen as the acoustic gain brought by the use of wind turbines.

The acoustic emergence is difficult to estimate because both the propagation of the sound of the wind turbines [START_REF] Guarnaccia | Wind turbine noise: theoretical and experimental study[END_REF] and the background noise rapidly fluctuate in an outdoor context. Nowadays, numerical simulation and measurement campaigns are used for that purpose but the calculated emergence is often inaccurate because of the sensitiveness of the acoustic variables to the fine topology and meteorological conditions [START_REF] Tickell | Wind turbine generator noise prediction-comparison of computer models[END_REF][START_REF] Van Den | Effects of the wind profile at night on wind turbine sound[END_REF]. To address this problem, the iEar project 1 introduced an online control system [START_REF] Dumortier | Acoustic control of wind farms[END_REF] based on permanent sound level meters that estimate the intensity of wind turbine noise and background noise in real time by means of a proprietary audio source separation system using audio input only. The output of this source separation system is not perfect, though, since wind turbine noise and background noise often exhibit similar spectral and spatial characteristics that make them difficult to separate. Uncertainties associated with source separation have received some interest for speech and speaker recognition [START_REF] Kolossa | Robust Speech Recognition of Uncertain or Missing Data[END_REF][START_REF] Delcroix | Cluster-based dynamic variance adaptation for interconnecting speech enhancement pre-processor and speech recognizer[END_REF][START_REF] Kallasjoki | Estimating uncertainty to improve exemplar-based feature enhancement for noise robust speech recognition[END_REF][START_REF] Tran | Nonparametric uncertainty estimation and propagation for noise robust ASR[END_REF][START_REF] Abdelaziz | Uncertainty propagation through deep neural networks[END_REF], but they have not been studied for environmental sounds yet, to the best of our knowledge.

In this paper, we propose a recursive Bayesian estimation framework for estimating the intensity of wind turbine noise and background noise that jointly takes into account the output of a source separation system (considered here as a black box) and a model of wind turbine noise emission and propagation. Our approach has the advantage of considering simultaneously the sound level meter measurements, the noise emission values provided by turbine manufacturers, and the results of audio source separation and acoustic simulation. We evaluate the suitability of different nonlinear Kalman filters to solve the estimation problem. Since the ground truth is unknown in real conditions, we validate our approach on real data obtained from an acoustic compliance study of a selected wind farm to which we added simulated uncertainties.

The rest of this paper is organised as follows. In Section 2, we present the background and the notations. Then we present the proposed recursive Bayesian model in Section 3. In Section 4, we present the conducted experiments with nonlinear Kalman filters tested on semi-simulated data. Finally perspectives and discussions are given in Section 5.

BACKGROUND AND NOTATIONS

Acoustic variables

We assume that acoustical engineers have set sound level meters in different locations indexed by j ∈ {1, . . . , J}. Each sound level meter measures the ambient noise b t,j , i.e. the total noise level at location j and time t in A-weighted decibels (dBA). Each time t corresponds to a period of 10 min.

The ambient noise is the sum of the wind turbine noise l t,j and the background noise r t,j at location j and time t. These two signals are assumed to be uncorrelated and to sum up on a linear scale (W/m2 ) [START_REF] Hansen | Fundamentals of acoustics[END_REF] which gives in dBA: b t,j = 10 log 10 10 l t,j 10 + 10 r t,j 10

.

(1)

Furthermore, the acoustic emergence e t,j is defined as

e t,j = b t,j -r t,j . (2)
This quantity represents the increase of the noise level due to the wind turbines and it has the meaning of an a posteriori signal-to-noise ratio.

Acoustic model

Besides acoustic measurements, acoustic simulations are conducted to predict the wind turbine noise l t,j with the help of the acoustic noise emission curves of the turbine manufacturer. These simulations are conducted by means of specific software such as Datakustik cadnaA 2 , that uses ray-tracing to simulate the propagation. To reduce the number of simulations, the variables c t and o t that respectively represent the turbine command variables (pitch angle, etc) and the meteorological variables (wind speed and direction), are reduced to a small number of number of discrete values denoted by C = {c 1 , ..., c N } and O = {o 1 , ..., o N } and the following acoustic formula is generally considered [START_REF] Rogers | Wind turbine acoustic noise[END_REF]:

l t,j 10 log 10 I i 10 a o t ,i,j 10
10 x c t ,o t ,i 10
. ( 3
)
where i ∈ {1, . . . , I} indexes the wind turbines, x ct,ot,i is the noise emitted by turbine i in dBA that both depends on the command c t in use and on the meteorological conditions o t , and a ot,i,j is an attenuation factor between wind turbine i and sound level meter j which depends on the meteorological conditions o t only.

Measurements and uncertainties

Finally, we assume that a source separation system separates the wind turbine noise and background noise signals and produces an estimate of the background noise level at each location and time. Neither the acoustic model nor the sound level measurements and the source separation outputs are exact: they suffer from uncertainties. To distinguish them from the ground truth values, we adopt the following notations. b * t,j and r • * t,j respectively denote the sound level measurement of b t,j and the audio source separation estimate of r t,j . They are stacked across locations into J × 1 vectors b * t and r * • t . x ct,ot,i , a ot,i,j , and r t,j respectively denote the real values of turbine noise emission, propagation, and background noise. They are stacked across turbines and/or locations into an I × 1 vector x t , an IJ × 1 vector a t , and a J × 1 vector r t . We denote by

z t = [b * T t r * • T t ]
T the overall 2J × 1 observation vector and by

y t = [x T t a T t r T t ]
T the overall (IJ +I+J)×1 state vector. Finally, w v denotes the uncertainty related to a given variable v, that is the deviation from the ground truth.

BAYESIAN ESTIMATION

The aim of the proposed Bayesian approach is to refine the prior distribution resulting from the acoustic model thanks to the observations that comprise of the measured ambient noise and the outputs of the source separation system. The posterior distribution is then inferred thanks to the Bayes theorem [START_REF] Bishop | Pattern Recognition and Machine Learning[END_REF]. In order to define the likelihood and the prior distribution, we model the uncertainties of the system as Gaussian random variables. This assumption is suitable since wind turbine manufacturers, sound level meter manufacturers, and simulation software systematically quantify uncertainties using standard deviations on the dBA scale.

Prior

For the prior, we use a first order auto-regressive moving average with external output vectors (ARMAX) model [START_REF] Bertsekas | Dynamic Programming: Deterministic and Stochastic Models[END_REF] for the emission, the propagation and the background noise:

y t =   x t a t r t   =   F x x t-1 + G x m xt + w xt F a a t-1 + G a m at + w at F r r t-1 + G r m rt + w rt   . (4)
This model takes into account both the temporal correlations (via the matrices F x , F a , and F r ) and the mean values m xt , m at , and m rt derived from the acoustic studies (via the matrices G x , G a , and G r ). The matrices F x , F a , F r , G x , G a , and G r can be determined by learning in principle, even though we haven't tried this yet and we chose a simple random walk model for r t (F x = 0) as well as simple model for the propagation (G x = 0 and G a = 0).

Likelihood

For the likelihood, we consider the following model: b * t,j = 10 log 10 I i 10 a t,i,j 10 10

x t,i 10 + 10

r t,j 10
+ w b * t,j (5)
r * • t,j = r t,j + w r • t,j + w b * t,j . (6)
We model the ambient noise b * t,j measured by the sound level meter at location j as a Gaussian deviation from the real ambient noise b t,j expressed as a nonlinear function of the state variables with (1)- [START_REF] Kuh | Applications of complex augmented kernels to wind profile prediction[END_REF]. Similarly, we model the background noise r * • t,j estimated by the source separation algorithm as a Gaussian deviation from the true r t,j . The source separation uncertainty w r * • t,j , and the measurement uncertainty w b * t,j are strongly correlated with each other since source separation is performed from the same sound level meter input. Therefore the overall observation covariance is given by:

Σ z = σ * 2 b I J σ * 2 b I J σ * 2 b I J (σ * 2 b + σ •2 r )I J (7)
with I J the identity matrix of size J.

Recursive Bayesian model

The proposed formulation now fits with the recursive Bayesian filtering framework. We can summarize ( 4)-( 6) with the general recursive Bayesian estimation model

(prior) y t = Fy t-1 + Gm y t + w y t (8) (likelihood) z t = h(y t ) + w zt (9)
with h the nonlinear function in ( 5)-( 6

with z 1:t = {z 1 , . . . , z t } the set of observations up to time t.

EXPERIMENTS AND RESULTS

Data

We validated our approach on data collected by VENATHEC SAS on a wind farm located in the north-east of France. For privacy reasons, the name and exact location of the wind farm are not mentioned. The wind farm is composed of 6 wind turbines (Vestas V112 -3.3 MW) and 5 measurement locations were set up. In this paper, we consider 19 h of measurements that correspond to 115 time frames of 10 min. Since the ground truth is unknown in real conditions, we generated the test data as follows. For background noise, we considered the real series r t measured when the turbines are stopped. For noise emission and propagation, we considered the mean values m xt and m at provided by the turbine manufacturer and cadnaA, respectively, and generated x t and a t by adding random Gaussian uncertainties w xt and w at . We then generated observations z t using ( 9) with random Gaussian measurement and separation uncertainties w b * t and w r • t . Such a simulation procedure is common in the wind turbine acoustics community. The uncertainties were assumed to be temporally and spatially uncorrelated and independent of i and j. Hence, we omit indices i and j in the following.

We generated two datasets. The first dataset contains 125 test cases with all 6 turbines and 5 measurement locations. The standard deviations σ x , σ a , and σ r • of the uncertainties w xt , w at , and w r • t are chosen in {0.5, 1.5, 2.5, 3.5, 4.5} and all combinations of values are considered. These are plausible values encountered in practice. In the second dataset, σ x , σ a , and σ r • are fixed to 3.5, but subsets of 1 to 6 turbines and 1 to 5 measurement locations are selected, resulting in 30 test cases. For both datasets, the standard deviation σ * b of the measurement uncertainty w b * is fixed to 1.5 dB, that is the standardized value for Type 1 sound level meters [START_REF]1: 2002 Electroacoustics-Sound level meters-Part 1[END_REF].

Tested nonlinear Kalman filters

We processed these data using various nonlinear Kalman filters with the following values for the model parameters. We set F x = I, F a = I, F r = 0, G x = 0, G a = 0, and G r = I. The covariance of w y t was defined as

Σ y =   σ 2 x I J 0 0 0 σ 2 a I IJ 0 0 0 Σ r   . (12)
with Σ r the covariance matrix of r t -r t-1 measured when the turbines are stopped. The values of m xt , m at , σ x , σ a , and σ r • were set as above. The assumptions behind these settings are also common in the wind turbine acoustics community. The Kalman filter [START_REF] Kailath | Linear estimation[END_REF] implements the general recursive Bayesian filter in (10)- [START_REF] Dumortier | Acoustic control of wind farms[END_REF]. It provides update rules for the posterior mean ŷt and covariance Σy t of the state vector y t given the sequence of observations z 1:t up to time t. In the case of a linear function h, these update rules are exact.

In our specific problem, the function h is nonlinear, hence we evaluated three different approximations : the extended Kalman filter [START_REF] Haykin | Kalman filtering and neural networks[END_REF], the unscented Kalman filter [START_REF] Julier | New extension of the Kalman filter to nonlinear systems[END_REF] and the central difference Kalman filter [START_REF] Van Der Merwe | Sigma-point Kalman filters for probabilistic inference in dynamic state-space models[END_REF]. The first one is the most common approximation: it relies on first-order Taylor development of h around a single point, namely the estimated state mean. When the covariance is large or h is strongly nonlinear it poorly approximates the posterior distribution and the estimation can be weak. The two other filters are ingenious improvements of the extended Kalman filter based on statistical linearisation. They approximate the prior distribution with a number of well-chosen points that are propagated through the nonlinearity to capture the mean and covariance. The two filters slightly differ in the choice of the points and the way they reconstruct the mean and covariance.

Evaluation metrics

We evaluated the results by means of the root mean square error (RMSE) between the estimated background noise value rt and the ground truth r t and the log-likelihood of the ground truth r t given the posterior distribution of the background noise, assumed to be Gaussian with mean rt (subvector of ŷt ) and covariance Σrt (submatrix of Σy t ). We also computed the RMSE and the log-likelihood of the emergence, where the ground truth e t was obtained from y t and the posterior mean êt and covariance Σet were obtained by propagating ŷt and Σy t through the nonlinearity via the nonlinear Kalman filter under test.

Results

Tables 1 and2 report the average RMSE on the background noise estimate and the emergence estimate. The former is similar to the intrinsic standard deviation of the sound level meter (σ * b = 1.5 dB), while the latter is much smaller due to the fact that the measurement uncertainty w b * t on b t and r t cancels when subtracting those two quantities from each other. It is also much smaller than the average uncertainty of the source separation system (2.5 dB). The central difference Kalman filter appears to perform best in terms of RMSE, however the posterior covariance is underestimated so that the unscented Kalman filter yields a better log-likelihood. Figure 1 shows the dependency of the RMSE on the uncertainties in the data and on the number of wind turbines and measurement locations. As expected, thanks to the fusion of the model uncertainties and the separation uncertainties operated by the nonlinear Kalman filter, the RMSE on the emergence estimate decreases with the uncertainties in the data and it is always smaller than the standard deviation of the separation system, except for (σ x , σ a , σ r • ) = (3.5, 3.5, 0.5) and (σ x , σ a , σ r • ) = (4.5, 4.5, 0.5) where the central difference approximation is poor. Also, the RMSE increases with the number of wind turbines and decreases with the number of measurement locations.

CONCLUSION

In this paper, we proposed a recursive Bayesian estimation framework adapted to the context of wind farms in order to improve a black-box source separation system. The proposed framework leverages the noise emission values provided by turbine manufacturers and the results of acoustic simulation. Our experiments showed that the central difference Kalman filter provides the lowest average RMSE on the emergence estimate and that this average RMSE is almost always smaller than the uncertainty of the source separation system alone. Future work includes addressing the poor performance of the central difference Kalman filter when the separation uncertainties are small and the model uncertainties are large, for instance by using iterated Kalman [START_REF] Sibley | The iterated sigma point Kalman filter with applications to long range stereo[END_REF] or particle filters [START_REF] Ristic | Beyond the Kalman filter: Particle filters for tracking applications[END_REF].

ACKNOWLEDGEMENT

This work was done with the collaboration and the support of VENATHEC SAS. In particular, the authors thank Jeremy Schild and Patrice Cornu from Venathec SAS for their implication. Experiments presented in this paper were carried out using the Grid'5000 testbed, supported by a scientific interest group hosted by Inria and including CNRS, RENATER and several Universities as well as other organizations (see https://www.grid5000.fr).

Fig. 1 .

 1 RMSE of the emergence estimate as a function of (a) the standard deviations σ x and σ a of the model uncertainties and the standard deviation σ r • of the source separation system on the first dataset and (b) the number of measurement locations and wind turbines on the second dataset. The central difference Kalman filter is used.

Table 1 .

 1 Average RMSE and log-likelihood of the background noise estimate on the first dataset.

	Filter	RMSE log-likelihood
	source separation estimate (without Kalman filter)	3.24	2.50 × 10 5
	extended filter	1.63	2.42 × 10 5
	unscented filter	1.63	2.41 × 10 5
	central difference filter	1.61	2.40 × 10 5
	Filter	RMSE log-likelihood
	source separation estimate (without Kalman filter)	2.87	4.2675 × 10 4
	extended filter	1.01	-0.69 × 10 5
	unscented filter	0.79	1.54 × 10 5
	central difference filter	0.76	-0.1 × 10 5

Table 2 .

 2 Average RMSE and log-likelihood of the emergence estimate on the first dataset.

http://www.venathec.com/fr/iear

http://www.datakustik.com/en/products/cadnaa