# Wind Turbine Quantification and Reduction of Uncertainties Based on a Data-Driven Data Assimilation Approach

Adrien Hirvoas, Clémentine Prieur, Élise Arnaud, Fabien Caleyron, Miguel Munoz Zuniga

► To cite this version:

Adrien Hirvoas, Clémentine Prieur, Élise Arnaud, Fabien Caleyron, Miguel Munoz Zuniga. Wind Turbine Quantification and Reduction of Uncertainties Based on a Data-Driven Data Assimilation Approach. Journal of Renewable and Sustainable Energy, 2022, 14 (5), pp.053303. 10.1063/5.0086255. hal-03855143v2

HAL Id: hal-03855143

https://ifp.hal.science/hal-03855143v2

Submitted on 16 Nov 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Wind turbine quantification and reduction of uncertainties based on a data-driven data assimilation approach

Adrien Hirvoas,¹,² Clémentine Prieur,² Élise Arnaud,² Fabien Caleyron,¹ and Miguel Munoz Zuniga³

¹) IFP Energies nouvelles, Rond-point de l'échangeur de Solaize, BP 3, 69360 Solaize, France

²) Univ. Grenoble Alpes, CNRS, Inria, Grenoble INP* LJK, 3800 GrenobleFrance

* Institute of Engineering Univ. Grenoble Alpes

³) IFP Energies nouvelles, 1 et 4 avenue de Bois-Préau, 92852 Rueil-MalmaisonFrance

(*Electronic mail: adrien.hirvoas@gmail.com)

(Dated: June 15, 2022)
In this paper, we propose a procedure for quantifying and reducing uncertainties impacting
numerical simulations involved in the estimation of the fatigue of a wind turbine struc-
ture. The present study generalizes a previous work carried out by the authors proposing
to quantify and to reduce uncertainties affecting the properties of a wind turbine model by
combining a global sensitivity analysis and a recursive Bayesian filtering a pproach. We
extend the procedure to include the uncertainties involved in the modeling of a synthetic
wind field. Unlike the model properties having a static or slow time-variant behavior, the
parameters related to the external solicitation have a non-explicit dynamic behavior which
must be taken into account during the recursive inference. A non-parametric data-driven
approach to approximate the non-explicit dynamic of the inflow related parameters is used.
More precisely, we focus on data assimilation methods combining a nearest neighbor or
analog sampler with a stochastic filtering method such as the ensemble Kalman filter. The
so-called data-driven data assimilation approach is used to recursively reduce the uncer-
tainties that affect the parameters related to both model properties and wind field. For the
approximation of the non-explicit dynamic of the wind inflow related parameters, in-situ
observations, obtained from a Light Detection And Ranging system and a cup-anemometer
device, are used. For the data-assimilation procedure, synthetic data simulated from the
aero-servo-elastic numerical model are considered. The next investigations will be to ver-
ify the procedure with real in-situ data.
1

# I. INTRODUCTION

2

3

4

5

6

7

8

9

10

11

12

A major challenge in wind energy industry is to propose robust designs withstanding unknown environmental conditions. Design standards (IEC, 2019) are mainly based on dynamic load simulations describing the structural behavior of the wind turbine under different wind and operational conditions weighted by their probability of occurrence. Most of the time the number of wind scenarios considered during the conception phase is moderate and far from exploring the set of environmental conditions. Moreover, the dynamic response of the structure and its lifetime can be affected by some uncertainties or evolution in the wind turbine properties. Consequently, the prediction of the operating wind turbine lifetime by taking into account all the inherent uncertainty is crucial. In that context, the quantification and reduction of uncertainties involved in the aero-servo-elastic numerical models play an important role to determine the effective fatigue loads of the turbine.

13

14

15

16

17

18

19

20

21

22

23

24

25

The approach introduced in this paper generalizes the one in Hirvoas et al. (2021) by taking into account the uncertainties affecting the parameters related to the wind inflow. It relies on a complete framework including a global sensitivity analysis, an identifiability analysis, and a recursive Bayesian inference approach. First, a surrogate based global sensitivity analysis through the estimation of Sobol' indices allows one to determine the most relevant input parameters in the variability of the fatigue loads of a wind turbine. After assessing the identifiability properties of these influential parameters, a second objective is to reduce their uncertainty by using an ensemble Kalman filter. Data assimilation allows one to gather all the information obtained from real time measurements of the physical system and from the numerical model. The procedure is closely related to the industrial concept of digital twin which consists in combining measurements from the wind turbine with a numerical model to build a digital equivalent of the real-world structure. However, unlike the model properties having a static or slow time-variant behavior, the parameters related to the external conditions have a dynamic that has to be learnt from data.

26

27

28

29

30

31

For design certification, offshore wind turbines in pilot farms are more and more monitored thanks to a large number of sensors. In that context, the measured data can be efficiently used in order to learn the non-explicit dynamic behavior of the wind parameters needed for numerical simulation. In the present work, we focus on non-parametric learning strategies. In the literature, several non-parametric methods have been developed such as regression machine learning (Brunton, Proctor, and Kutz, 2016), echo state networks (Pathak et al., 2018) or more recently
32 residual neural networks (Bocquet, Farchi, and Malartic, 2020). Our study investigates an analog
33 forecasting method relying on the principle of the nearest neighbors (Lorenz, 1969). Analogue
34 methods have the advantage of being computationally inexpensive, so that ensemble forecasts,
35 used in data assimilation, are easy to make. The aforementioned non-parametric procedure has
36 been firstly coupled with data assimilation filtering schemes in (Tandeo et al., 2015) and further
37 detailed by Lguensat et al. (2017). In the present work, we propose an algorithm, developed in
38 (Hirvoas et al., 2021), interfacing Python library AnDA<fnref n="1" /> combining analog forecasting with en-
39 ensemble data assimilation. The algorithm we propose takes profit of the parallelization capabilities
40 of high performance computing architectures which allows for example to evaluate the real-time
41 damage of an operating wind turbine using a digital twin.

42 The outline of this paper is as follows. Firstly, Section II describes the different uncertainties
43 involved in the framework of this study. In Section III, the theoretical framework of data-driven
44 data assimilation with a specific focus on the ensemble Kalman filtering scheme coupled with the
45 analog forecasting strategy is detailed. Finally, results of an application of this complete procedure
46 of uncertainty quantification and reduction to a reference wind turbine are presented in Section IV.

## II. CONTEXT

### A. Uncertainty in wind turbine modeling

49 Before their exploitation, wind turbine rotors are designed thanks to a site classification strat-
50 egy. It relies on design standard classes characterized by the reference turbulence intensity $I_{ref}$,
51 defined as the mean turbulence intensity expected at 15 m/s mean wind speed and the reference
52 wind $\bar{u}_{ref}$, defined as the extreme 10-minute average wind speed with a recurrence period of 50
53 years. In the IEC-61400-1 standard (IEC, 2019), two safety classes are considered. The first one,
54 named as normal safety class, allows one to cover most applications by giving specific values for
55 $I_{ref}$ and $\bar{u}_{ref}$.
56 In Table I, the corresponding values for the nine categories of the normal safety are
given. The proposed parameter values are supposed to represent many different sites and conse-
57 quently do not give a precise representation of a specific site. The second category is mentioned
58 as a special safety class S which allows to consider site-specific values for the wind speed and
59 turbulence terms.

<footnote n="1">¹ see https://github.com/ptandeo/AnDA</footnote>
<table><caption>Table I. Safety class design classification of the wind turbines: the normal safety class containing nine categories from I-A to III-C and the special safety class S (IEC, 2019) giving the reference turbulence intensity $I_{ref}$ and the reference wind speed $\bar{u}_{ref}$</caption><thead><tr><th colspan="2">Wind turbine<br>Generators class</th><th>I</th><th>II</th><th>III</th><th>S</th></tr></thead><tbody><tr><td colspan="2">$\bar{u}_{ref}$ [m/s]</td><td>50</td><td>42.5</td><td>37.5</td><td></td></tr><tr><td rowspan="3">Turbulence<br>class</td><td>A</td><td rowspan="3">$I_{ref}$ [-]</td><td>0.16</td><td></td><td rowspan="3">Site-specific<br>values</td></tr><tr><td>B</td><td>0.14</td><td></td></tr><tr><td>C</td><td>0.12</td><td></td></tr></tbody></table>

60

61

62

63

64

65

66

67

68

69

For each class, the design relies on numerical aero-servo-elastic simulations under different environmental and operational conditions. They allow one to estimate the ultimate and fatigue loads in order to certify the structural integrity. Nevertheless, operating wind turbines experience real wind and operational conditions that can differ from the ones mentioned in the design standard classes. Consequently, there is a need for an estimation of the fatigue life of the structure based on the real wind solicitation seen by the structure. Moreover, the wind turbine itself can present some uncertainties or evolution in its mechanical properties (defaults appearance, degradation with time) that will affect the dynamic response of the structure and its lifetime.

70

71

72

73

74

75

76

77

As a consequence, these aero-servo-elastic numerical models involve many uncertain and potentially variable over time parameters. The ubiquitous uncertainty may be found in the parameters of the wind turbine numerical model as well as in the external conditions. To ensure the tracking of fatigue and defaults of an operating wind turbine structure, it is important to quantify the impact of these uncertainties on predictions and then to reduce them based on the combination of measurements and model predictions. For that purpose, the field of uncertainty quantification is well-adapted. Hereafter, we propose to determine the sources of uncertainties affecting the wind field parameters and the wind turbine numerical model properties.

78

# 1. Uncertainties in wind field modeling

79

80

First, the uncertainty of wind field parameters has to be determined. In our context, these parameters are used to characterize a synthetic three-dimensional turbulent wind field based on
81 the Kaimal spectrum (Kaimal et al., 1972) having a one-sided power spectral density defined as:

$$S_k(f) = \frac{4\sigma_k^2 \frac{L_k}{\bar{u}}}{(1 + 6f \frac{L_k}{\bar{u}})^{\frac{5}{3}}}, \quad 82$$

83 where f is the frequency of occurrence, the subscript $k \in \{u, v, w\}$ represents the turbulent longitudinal, crosswise or vertical components, $L_k$ is the Kaimal length scale, u is the longitudinal mean wind speed at hub height, and $\sigma_k$ is the standard deviation of the wind speed.

86 The wind inflow over the swept area is generated based on a grid of points thanks to an exponential spatial coherence method (Jonkman, 2009). The related coherence function for the longitudinal wind component of two distinct points i and j separated by a distance $\Delta r$ on a plan perpendicular to the wind direction is defined as:

$$\coth_{i,j}(f) = \exp \left( -a \left( \frac{\Delta r}{z_m} \right)^\gamma \sqrt{ \left( \left( \frac{f \Delta r}{\bar{u}_m} \right)^2 + \left( \frac{b' \Delta r}{L_u} \right)^2 \right) } \right), \quad (1)$$

91 where $z_m$ and $\bar{u}_m$ are respectively the mean height of the two points and the mean of the wind speeds of the two points, a and $b'$ are respectively the input coherence decrement and offset parameter, and $\gamma$ is the coherence exponent.

94 Eight input parameters related to the wind field have been identified to be tainted by uncertainties, see Table II. We have considered the mean and the standard deviation of the wind speed at hub height, the vertical wind shear exponent, the mean wind inflow direction relative to the wind turbine in terms of vertical or horizontal inflow angles, and the longitudinal turbulence length scale parameter. Moreover, we have supposed as unknown the input coherence decrement and offset parameter.

100 In an operational context, some information on the mean and standard deviation of the wind speed at hub height can be obtained from 10-minute data measured from a nacelle mounted anemometer. Nevertheless, these measurements are known to be very perturbed and never fully describe the parameters of interest due mainly to the wake effect of the rotor and the non-exact transfer function used to retrieve them. In this work, we assume that the 10-minute mean and standard deviation free wind speed can be obtained from the 10-minute data obtained from the anemometer modulo an additive error term. So that the mean free wind speed at hub height can be
obtained from the anemometer as:

$$ \bar{u} = \bar{u}_{scada} + \Delta\bar{u}, \qquad (108) $$

where $\bar{u}_{scada}$ is the 10-minute mean wind speed obtained from the anemometer mounted on the wind turbine nacelle and $\Delta\bar{u}$ is an additive error assumed to follow the distribution defined in Table II.

In a similar manner, the free wind speed standard deviation can be obtained from the measurement obtained by the anemometer mounted on the nacelle of the wind turbine as:

$$ \sigma_u = \sigma_{scada} + \Delta\sigma_u, \qquad (114) $$

where $\sigma_{scada}$ is the 10-minute standard deviation wind speed obtained from the nacelle anemometer of the wind turbine nacelle and $\Delta\sigma_u$ is an additive error assumed to follow the distribution defined in Table II.

Unless having high frequency supervisory control and data acquisition system (SCADA) data, no information can be obtained on the other parameters. An investigation of the distribution of the uncertainty affecting these remaining wind inflow parameters has to be properly made. The vertical wind shear is modeled with the following power law that uses a shear coefficient $\alpha$:

$$ \bar{u}_z = \bar{u} \times \left( \frac{z}{z_{hub}} \right)^\alpha, \qquad (122) $$

where $\bar{u}$ is the prescribed hub-height mean wind velocity, $z$ is the vertical distance from the ground surface, $z_{hub}$ is the hub height, and $\alpha$ is the vertical wind shear coefficient. We adapt the **the mean and the standard deviation of the Gaussian distribution** proposed by Dimitrov, Natarajan, and Kelly (2015) for the 10-minute vertical wind shear exponent, such that:

$$ \mu_{\alpha} = 0.088(\ln(\bar{u}_{scada}) - 1), \qquad (2) $$

$$ \sigma_{\alpha} = 1/\bar{u}_{scada} $$

Table II summarizes the wind-inflow parameters that we consider unknown and their respective uncertainty modeling. In particular, we defined the probability distributions of the parameters used in the exponential coherence model defined in Equation (1).
<table><caption>Table II. Wind field parameters - uncertainties affecting the inputs of the wind turbine model. $\mathcal{U}$: uniform distribution and $\mathcal{G}$: Gaussian distribution.</caption><thead><tr><th>Input</th><th>Variable</th><th>Unit</th><th>Distribution</th><th colspan="2">Parameters</th><th>REF</th></tr></thead><tbody><tr><td>Error of hub mean wind speed SCADA vs undisturbed inflow</td><td>Δū</td><td>[m/s]</td><td>U</td><td>Min: −0.15·ū<sub>scada</sub> Max: 0.15·ū<sub>scada</sub></td><td></td><td>Expert knowledge</td></tr><tr><td>Error of hub standard deviation SCADA vs undisturbed inflow</td><td>Δσ<sub>u</sub></td><td>[m/s]</td><td>U</td><td>Min: −0.2·σ<sub>scada</sub> Max: 0.2·σ<sub>scada</sub></td><td></td><td>Expert knowledge</td></tr><tr><td>Vertical wind inflow angle</td><td>φ<sub>v</sub></td><td>[°]</td><td>U</td><td>Min: 0 Max: 10</td><td></td><td>Expert knowledge</td></tr><tr><td>Horizontal wind inflow angle</td><td>φ<sub>h</sub></td><td>[°]</td><td>U</td><td>Min: −15 Max: 15</td><td></td><td>Expert knowledge</td></tr><tr><td>Longitudinal turbulence length scale</td><td>Λ<sub>u</sub></td><td>[m]</td><td>U</td><td>Min: 20 Max: 170</td><td></td><td>(Dimitrov, Natarajan, and Mann, 2017) (Solari and Piccardo, 2001)</td></tr><tr><td>Decrement parameter of coherence model</td><td>a</td><td>[-]</td><td>U</td><td>Min: 1.5 Max: 26</td><td></td><td>(Robertson et al., 2019a)</td></tr><tr><td>Offset parameter of coherence model</td><td>b'</td><td>[-]</td><td>U</td><td>Min: 0 Max: 0.17</td><td></td><td>(Robertson et al., 2019a) (Saranyasoontorn, Manuel, and Veers, 2004)</td></tr><tr><td>Vertical wind shear exponent</td><td>α</td><td>[-]</td><td>G</td><td>μ = μ<sub>α</sub> σ = σ<sub>α</sub>, see Equation (2)</td><td></td><td>(Dimitrov, Natarajan, and Kelly, 2015)</td></tr></tbody></table>

131

## 2. Uncertainties in aero-servo-elastic numerical model

132

133

134

135

136

137

138

139

140

141

142

143

144

145

146

147

148

149

150

151

Moreover, as suggested in Hirvoas et al. (2021), a total of twelve parameters can be considered as uncertain in the aero-servo-elastic wind turbine numerical model properties. All these input parameters are assumed to be independent of one another with Gaussian or truncated Gaussian distributions obtained from expert knowledge or literature. Considering the support structural properties of the turbine model, we have selected six parameters: as nacelle mass and center of mass, tower Rayleigh damping, inertial nacelle and drive-train torsion stiffness. Lastly, the geometry of the tower, resulting from fabrication tolerances, has been also included in these uncertainties by uniformly scaling the distributed tower thickness. The probability distribution of this last mentioned parameter is determined by changing the first fore-aft tower frequency mode by ±10% of its nominal value. The uncertainties in blade structural properties have been represented using five parameters. The blade structural responses have led to the definition of the uncertainty range. Indeed, the frequency of the edge-wise (EW) and flap-wise (FW) modes are changed about 10% each from their reference value. These modifications of the frequency modes are done by uniformly scaling the associated stiffness and the distributed blade mass of all blades. For anomaly diagnosis, blade mass imbalance effects have been also included by applying a different mass factor value to each blade. One blade's mass property is modified to be a value that is higher than the nominal value, and another one modified to a lower value. The third blade remains unchanged at the nominal value. Finally, for the individual blade pitch error, a constant offset angle is applied to two of the blades, respectively above and below the nominal value. These different parameters are considered independent from each other. Table III gathers information about the probability
152 distribution of each of these paremeters.

<table><caption>Table III. Model parameters - uncertainties affecting the inputs of the wind turbine model. *G*: Gaussian distribution and *FG*: Truncated Gaussian distribution where *a* and *b* are the cut-off parameters.</caption><thead><tr><th>Input</th><th>Variable</th><th>Unit</th><th>Distribution</th><th>Parameters</th><th>REF</th></tr></thead><tbody><tr><td>Nacelle mass</td><td>N<sub>mass</sub></td><td>[kg]</td><td>G</td><td>μ = 6.90e+04 σ = 2.30e+03</td><td>(Witcher, 2017)</td></tr><tr><td>Nacelle center of mass</td><td>N<sub>CMx</sub></td><td>[m]</td><td>G</td><td>μ = 1.00 σ = 3.35e-02</td><td>(Robertson et al., 2019b)</td></tr><tr><td>Tower thickness adding coefficient</td><td>e</td><td>[-]</td><td>G</td><td>μ = 0.00 7.00e-01</td><td>Expert knowledge ±10% 1 FA</td></tr><tr><td>Tower Rayleigh damping</td><td>β<sub>TR</sub></td><td>[-]</td><td>FG</td><td>μ = 2.55 σ = 0.82<br>a = μ - 3σ b = μ + 3σ</td><td>(Koukoura, 2014)</td></tr><tr><td>Inertial nacelle</td><td>I<sub>zz</sub></td><td>[kg·m<sup>2</sup>]</td><td>G</td><td>μ = 7.00e+05 σ = 2.33e+04</td><td>Expert knowledge ±10% μ</td></tr><tr><td>Drive-train torsional stiffness</td><td>K<sub>D</sub></td><td>[N·m<sup>2</sup>/rad]</td><td>G</td><td>μ = 9.08e+09 σ = 3.03e+07</td><td>(Holierhoek et al., 2010)</td></tr><tr><td>Blade flap wise stiffness</td><td>α<sub>BF</sub></td><td>[N·m<sup>2</sup>]</td><td>G</td><td>μ = 1.00 σ = 3.33e-02</td><td>Expert knowledge ~ ±10% 1 FW</td></tr><tr><td>Blade edge wise stiffness</td><td>α<sub>BE</sub></td><td>[N·m<sup>2</sup>]</td><td>G</td><td>μ = 1.00 σ = 3.33e-02</td><td>Expert knowledge ~ ±10% 1 EW</td></tr><tr><td>Blade mass coefficient</td><td>α<sub>mass</sub></td><td>[-]</td><td>G</td><td>μ = 1.00 σ = 1.67e-02</td><td>(Witcher, 2017)</td></tr><tr><td>Blade Rayleigh damping</td><td>β<sub>BR</sub></td><td>[-]</td><td>FG</td><td>μ = 1.55 σ = 4.83e-01<br>a = μ - 3σ b = μ + 3σ</td><td>(Robertson et al., 2019b)</td></tr><tr><td>Blade mass imbalance</td><td>η<sub>B</sub></td><td>[-]</td><td>G</td><td>μ = 2.50e-02 σ = 8.33e-03</td><td>(Robertson et al., 2019b)</td></tr><tr><td>Individual pitch error</td><td>Ω</td><td>[°]</td><td>G</td><td>μ = 0.10 σ = 3.33e-02</td><td>(Simms et al., 2001)</td></tr></tbody></table>

153

154

## B. Methodology for uncertainty quantification

In the monitoring context of an operating wind turbine, one of the major challenges is to predict the remaining lifetime of the structure. Hence, the current study focuses on a complete framework first quantifying and then reducing in a recursive fashion the uncertainties affecting the damage loads obtained from an aero-servo-elastic simulation. Hereafter, we will focus on the estimation of the effective damage equivalent load (DEL) describing the fatigue behavior of the wind turbine at some specific locations. The DEL is obtained by considering the internal loads and is defined as a virtual load amplitude that would create, in reference regular cycles, the same damage as the considered irregular load history.

The aim of the work in this article is to generalize the complete methodology proposed in (Hirvoas et al., 2021) for quantifying and reducing the uncertainties affecting a wind turbine numerical model by handling wind turbine model properties in addition to wind inflow uncertainties,
167 respectively denoted by $\mathbf{x}^1$ and $\mathbf{x}^2$ in Figure 1.

<figure><img src="image_2.png" /><figcaption>Figure 1. Wind turbine modeling framework, where $\mathbf{x}$ is the extended vector gathering both uncertainties from wind inflow parameters and model properties.</figcaption></figure>

The procedure relies on a global sensitivity analysis (GSA) based on Sobol' index estimation and a recursive Bayesian inference procedure to reduce the uncertainties, similarly as in (Hirvoas et al., 2021). In order to alleviate the computational cost of index estimation during the sensitivity analysis of the fatigue loads, the aero-servo-elastic time-consuming numerical model is approximated by a surrogate. A major challenge in building such surrogate model relies on the fact that the turbulent wind inflow realization causes variations in the quantities of interest obtained from the model. Thus, to take into account the inherent variability on the turbine response induced by different turbulent wind field realizations, the approach focuses on the use of heteroscedastic Gaussian process regression models. Then, a recursive reduction of the influent parameter uncertainties based on an ensemble Kalman filter is proposed. This data assimilation filtering method is computationally efficient with high-performance computing tools which is a major advantage for online calibration of time-consuming codes, such as aero-servo-elastic wind turbine models. Nevertheless, a challenge in this kind of inverse problem is to determine whether the measurements are sufficient to unambiguously determine the parameters that generated the observations, i.e., identifiability properties. In that context, GSA is also proposed to detect non identifiable parameters considering the current measurements.

The main contribution of the presented work is the inference of parameters involved in both the model properties of the wind turbine having a static or slow evolution and the short-term wind inflow varying at each inference iteration of 10-minute. To take into account the non-explicit dynamics of the parameters related to the wind inflow in the recursive inference procedure, the study relies on a data-driven approach combining a K-nearest neighbors with an ensemble Kalman filtering scheme. In the next section, we propose to describe this data-driven procedure used in our model calibration strategy.
## III. DATA-DRIVEN DATA ASSIMILATION

### A. Data assimilation

State-space model (SSM) is a useful framework to perform recursive inference strategy such as sequential data assimilation techniques (Bertino, Evensen, and Wackernagel, 2003; Durbin and Koopman, 2012; Hirvoas et al., 2021). In order to take into account the information obtained from the SCADA system of the wind turbine, we consider the SSM formulation involving forcing variables defined $\forall k \in \mathbb{N}^*$ as:

$$ \mathbf{x}_{(k)} = \mathcal{M}_{(k-1,k)}(\mathbf{x}_{(k-1)}) + \boldsymbol{\varepsilon}_{(k)}^{m}, \quad (3) $$

$$ \mathbf{y}_{(k)} = \mathcal{H}_{(k)}(\mathbf{x}_{(k)}, \mathbf{v}_{(k)}) + \boldsymbol{\varepsilon}_{(k)}^{o}, \quad (4) $$

where $\mathbf{y}_{(k)}$ corresponds to the observation at step $k$ and $\mathbf{x}_{(k)}$ is a $p$-dimensional vector representing the hidden-state variables which depict the input parameters in the aero-servo-elastic numerical model such that:

$$ x = \begin{pmatrix} x^1 \\ x^2 \end{pmatrix} $$

with $\mathbf{x}_{(k-1)}^1$ and $\mathbf{x}_{(k-1)}^2$ respectively the uncertain parameters for the wind inflow conditions, described in Table II, and the model properties, described in Table III. The model denoted by $\mathcal{M}$ (potentially nonlinear) allows us to describe the dynamic behavior of the hidden process. The model error $\boldsymbol{\varepsilon}_{(k)}^m$ is supposed to be a Gaussian white noise of zero mean and of covariance $\mathbf{Q}_{(k)}$, modeling the uncertainties related to the dynamics model structure. The propagator $\mathcal{H}$ relates the hidden-state vector to the measured observations and contains some forcing variables $\mathbf{v}_{(k)}$, e.g., mean wind speed obtained from the anemometer of the wind turbine. The sources of errors in the observation model defined in Equation (4) are reflected by the Gaussian white noise of zero mean and of covariance $\mathbf{R}_{(k)}$, denoted by $\boldsymbol{\varepsilon}_{(k)}^o$, and assumed to be independent of the model error $\boldsymbol{\varepsilon}_{(k)}^m$. This SSM formulation can be represented thanks to the directed graph given below.
$$ (\mathcal{M}_{(k-1,k)}, \mathbf{Q}_{(k)}) $$

$$ \begin{array}{ccccccc} \text{Hidden-state} & \cdots & \rightarrow & \mathbf{x}_{(k-1)} & \rightarrow & \mathbf{x}_{(k)} & \rightarrow & \mathbf{x}_{(k+1)} \rightarrow \cdots \\ & & & \downarrow & & \downarrow & (\mathcal{H}_{(k)}, \mathbf{v}_{(k)}, \mathbf{R}_{(k)}) & \downarrow \\ \text{Observations} & & & \mathbf{y}_{(k-1)} & & \mathbf{y}_{(k)} & & \mathbf{y}_{(k+1)} \end{array} $$

## B. Data-driven data assimilation

In many situations, the dynamical model $\mathcal{M}$ is numerically intractable or unknown. In the literature different studies have been conducted to emulate this propagator, used in Equation (3), from historical data. Several surrogate techniques have been employed for the reconstruction of nonlinear dynamics model of chaotic system. Authors in (Tandeo et al., 2015) propose a K-nearest neighbors based method, also known as the analog strategy in meteorology or geoscience community. Nevertheless, it has been argued that methods relying on a K-nearest neighbors technique are plagued by the curse-of-dimensionality, i.e., fails in very high dimensional applications (Friedman, 1997; Chen, 2009). Consequently, other non-parametric surrogate modeling approaches have been investigated to learn the underlying dynamics by using for example regression machine learning (Brunton, Proctor, and Kutz, 2016), echo state networks (Pathak et al., 2018) or more recently residual neural networks (Bocquet, Farchi, and Malartic, 2020).

Due to the limited dimension of our inference problem, we have decided to investigate and to use the analog forecasting strategy coupled with data assimilation proposed in (Tandeo et al., 2015; Hamilton, Berry, and Sauer, 2016; Lguensat et al., 2017). Analog forecasting is related to the notion of atmospheric predictability introduced by Lorenz (1969). Later, this approach has been widely used in several atmospheric, oceanic, and climate studies (Toth, 1989; Alexander et al., 2017; Ayet and Tandeo, 2018). Hereafter, we detail the principle of the analog forecasting technique.

The main idea of the methodology is to substitute the dynamical model in Equation (3) by a data-driven model relying on an analog forecasting operator, denoted by $\mathcal{A}$, such as :

$$ \forall k \in \mathbb{N}^*, \begin{cases} \mathbf{x}_{(k)} = \mathcal{A}_{(k-1,k)}(\mathbf{x}_{(k-1)}) + \boldsymbol{\varepsilon}_{(k)}^m. \\ \mathbf{y}_{(k)} = \mathcal{H}_{(k)}(\mathbf{x}_{(k)}, \mathbf{v}_{(k)}) + \boldsymbol{\varepsilon}_{(k)}^o \end{cases}. $$
The analog forecasting principle consists in searching for one or several similar situations of the current hidden-state vector that occurred in historical trajectories of the system of interest, then retrieving the corresponding successors of these situations, and finally assume that the forecast of the hidden-state can be retrieved from these successors. Consequently, this strategy requires the existence of a representative catalog of historical data, denoted by $\mathcal{C}$. The reference catalog is formed by pairs of consecutive hidden-state vectors, separated by the same lag (Fablet et al., 2017). The first component of each pair is named as the analog (denoted by $\mathbf{a}$) while the corresponding state is referred to as the successor (noted as $\mathbf{s}$). The corresponding representative dataset of hidden-state sequences can be written as:

$$ \mathcal{C} = \{(\mathbf{a}_i, \mathbf{s}_i), i = [1 \cdots P]\}, \text{ with } P \in \mathbb{N}^*. $$

This historical catalog can be constructed using observational data recorded using in-situ sensors as well as from numerical simulations. Based on this database, the analog forecasting operator $\mathcal{A}$ is a non-parametric data-driven sampling of the state from iteration $k-1$ to iteration $k$. Three analog forecasting operators have been originally proposed by the authors in Lguensat et al. (2017). They are all based on nearest neighbors of the hidden-state in the reference catalog $\mathcal{C}$ weighted thanks to a kernel function. Among the different kernels, Chau, Ailliot, and Monbet (2021) propose to use a tricube kernel which has a compact support and is smooth at its boundary. Throughout this article, as selected by Lguensat et al. (2017), a radial basis function (also known as Gaussian kernel, squared exponential kernel, or exponentiated quadratic) is considered and defined as:

$$ g(\mathbf{u}, \mathbf{v}) = \exp \left( -\lambda \|\mathbf{u} - \mathbf{v}\|^2 \right), \qquad (5) $$

where $(\mathbf{u}, \mathbf{v})$ are two distinct variables in the hidden-state space, $\lambda$ is a scale parameter, and $\|\cdot\|$ is the Euclidean distance or any other relevant distance function for our application. The kernel choice is case dependent. The Gaussian kernel used hereafter is isotropic and parameterized by $\lambda$ allowing to easily control the bandwidth.

Let us denote by $\{\mathbf{a}_n\}_{n \in \mathcal{J}}$ the K-nearest neighbors (also known as analog situations) of a given hidden-state at iteration $k-1$, where $\mathcal{J} = \{i_1, \cdots, i_K\}$ contains the K indices of these situations. From the reference catalog $\mathcal{C}$, one can retrieve the corresponding successors $\{\mathbf{s}_n\}_{n \in \mathcal{J}}$. Then for every pair of analog and successor $(\mathbf{a}_n, \mathbf{s}_n)_{n \in \mathcal{J}}$, a normalized kernel weight $(\omega_n)_{n \in \mathcal{J}}$ can be
assigned such that:

$$ \omega_n = \frac{g(\mathbf{x}_{(k-1)}, \mathbf{a}_n)}{\sum_{j=1}^{K} g(\mathbf{x}_{(k-1)}, \mathbf{a}_{i_j})} $$

This term provides more importance to pairs that are best suited according to the kernel function for the estimation of the hidden-state $\mathbf{x}_{(k)}$ in the K-nearest neighbors obtained from the catalog. Nevertheless, the parametrization of this weight is highly dependent on the kernel function. Moreover in the context of the Gaussian kernel as defined in Equation (5), the normalized kernel weight involves the choice of the number of nearest neighbors K and the scale parameter $\lambda$. Two common strategies in the statistic field are used for the estimation of K: either a distance threshold in order to consider the nearest neighbors which respect it, or an expert based number of analogs (Peterson, 2009). In our work, we consider the last strategy for simplicity. As proposed by Lguensat et al. (2017), the scale parameter can be fixed following the adaptive rule defined as:

$$ \lambda = \frac{1}{\text{md}(\mathbf{x}_{(k-1)})}, $$

where md($\mathbf{x}_{(k-1)}$) is the median distance between the hidden-state at iteration $k-1$ and its K nearest neighbors. Nevertheless, a more sophisticated procedure not used hereafter, based on a cross-validation procedure, can be employed to optimize the choice of these hyper-parameters.

Three analog forecasting operators $\mathcal{A}$ have been defined in (Lguensat et al., 2017). The complexity of model studied and the available computational resources are the two main constraints that will drive the choice of one forecasting operator over the others. For example in situations facing some extreme values of the hidden-state based on the available catalog, the locally-constant gives poor results due to the fact that the forecasting estimate is held in the range of K-nearest neighbors. In that context, the locally-incremental and the locally-linear forecasting operators are much more efficient (Lguensat et al., 2017). A graphical representation of the three different analog forecasting operators for a 2-dimensional hidden-state is given in Figure 2. In this example, the underlying dynamics model has a simple polynomial form and the analogs are obtained by using a normal distribution sampling centered on the real value of the hidden-state at iteration $k-1$. In this research work, we have decided to focus on the locally-linear operator.

The locally-linear forecasting operator consists of performing a weighted least square linear regression between the K-nearest neighbors and their corresponding successors in the catalog $\mathcal{C}$. The multivariate linear regression provides a slope matrix of size $p \times p$ denoted by $\alpha$, a vector
intercept of size $p \times 1$ designated hereafter by $\beta$, and residuals defined as the following vectors
$\forall j \in [1, \dots, K], \mathbf{s}_{ij} = (\alpha \mathbf{a}_{ij} + \beta)$. The Gaussian sampling resorts to:

$$ \mathbf{x}_{(k)}^f \sim \mathcal{N}(\mu_{(k)}^{\text{LL}}, \Sigma_{(k)}^{\text{LL}}), $$

where the mean forecast is $\mu_{(k)}^{\text{LL}} = \alpha \mathbf{x}_{(k-1)} + \beta$, and $\Sigma_{(k)}^{\text{LL}}$ is the weighted empirical covariance of the residuals thanks to the normalized kernel weight $(\omega_n)_{n \in \mathcal{J}}$.

<figure><img src="image_3.png" /><figcaption>Figure 2. Analog forecasting operator strategies. The hidden-state real values $\mathbf{x}_{(k-1)}$ and its forecast $\mathbf{x}_{(k)}$ are represented by full circles. Analogs are displayed in colored down-pointing triangles and successors in up-pointing triangles with their equivalent colors. Each triangle size is proportional to the normalized kernel weight. The ellipsoids in black and red represent respectively the 95 % confidence intervals of the hidden state distribution before and after the analog forecasting strategy.</figcaption></figure>

## C. Non-parametric EnKF method

Hereafter, we propose to describe the data assimilation framework coupled with the analog forecasting method firstly proposed by Tandeo et al. (2014) and further detailed in (Lguensat et al., 2017). Data assimilation methods allow us to combine all the sources of information obtained from a physical model and observations. In particular, sequential data assimilation techniques, also known as filtering approaches consist of estimating the filtering posterior distribution of the current hidden-state knowing past and present observations $p_{\mathbf{X}_{(k)}|\mathbf{Y}_{(1:k)}}(\mathbf{x}_{(k)}|\mathbf{y}_{(1:k)})$ where $\mathbf{Y}_{(1:k)} = [\mathbf{Y}_{(1)}, \dots, \mathbf{Y}_{(k)}]$. Different methods are available in order to compute the filtering distribution of interest. In the context of linear Gaussian state-space models, Kalman filter methods can
be considered to provide the exact filtering methods (Kalman, 1960; Brown, 1986; Harvey, 1990; Haykin, 2004; Wells, 2013). Nevertheless in real applications, the linear assumption is often un- realistic and more sophisticated Kalman-based approaches have to be used (Julier and Uhlmann, 1997; Evensen, 2009). In particular, the ensemble Kalman filter (EnKF) which is a Monte Carlo variant relying on an ensemble of members to represent the statistics. This sequential Monte Carlo filter, introduced by Evensen (1994), is widely used in data assimilation applications to take into account the nonlinearities in the state-space formulation and to handle the high dimensional prob- lems (Houtekamer and Mitchell, 2001; Snyder and Zhang, 2003; Aanonsen et al., 2009). The principle of the EnKF is to sequentially update all members of the ensemble by means of a cor- rection term relying on the Kalman gain which allows one to blend the model responses and the observations at a given iteration, see Evensen (2003). Due to the fact that this approach is based on an ensemble, it is inherently well-adapted to parallelization which is a crucial advantage with the current high-performance computing architectures for the inference of time-consuming numerical models (Houtekamer, He, and Mitchell, 2014).

Thus, we present the formulation of a non-parametric EnKF method, also known as analog EnKF (AnEnKF), see (Tandeo et al., 2014; Lguensat et al., 2017). The procedure is similar to the stochastic ensemble Kalman recursion (Evensen, 2009). Nevertheless, the main difference of the AnEnKF occurs for the forecast step where the non-parametric data-driven sampling, i.e., the analog forecasting operator, is used instead of the dynamic model $\mathcal{M}$ in Equation (3). The analog ensemble Kalman filter applies one of the three analog forecast sampling strategies, presented in (Lguensat et al., 2017), to each analysis member of the ensemble to generate a forecast term at each iteration. Then, the equations used in the procedure are equivalent to the EnKF strategy. At each iteration during the analysis step, each forecast member of the ensemble is corrected by computing $\mathbf{x}_{(k)}^{a(i)} = \mathbf{x}_{(k)}^{f(i)} + \mathbf{K}_{(k)} (\mathbf{y}_{(k)}^{(i)} - \mathcal{H}_{(k)}(\mathbf{x}_{(k)}^{f(i)}, \mathbf{v}_{(k)}))$ where $\mathbf{K}_{(k)} = \mathbf{P}_{(k)}^{f} \mathbf{H}_{(k)}^{T} (\mathbf{R}_{(k)} + \mathbf{H}_{(k)} \mathbf{P}_{(k)}^{f} \mathbf{H}_{(k)}^{T})^{-1}$ is named as the Kalman Gain with $\mathbf{P}_{(k)}^{f}$ the forecast covariance matrix and $\mathbf{H}_{(k)}$ the observation operator. Due to the nonlinearity of the model $\mathcal{H}_{(k)}$, the terms $\mathbf{P}_{(k)}^{f} \mathbf{H}_{(k)}^{T}$ and $\mathbf{H}_{(k)} \mathbf{P}_{(k)}^{f} \mathbf{H}_{(k)}^{T}$ are respectively empirically estimated based on the ensemble members. The ensemble Kalman filter coupled with the analog forecasting strategy is detailed in Algorithm 1.
**Algorithm 1** Ensemble Kalman Filter with analog forecast methodology, so-called AnEnKF.

1: Input:

2: number of members in the ensemble $N_{ens}$;

3: number of inference iterations $T$;

4: catalog $\mathcal{C}$ and number of nearest neighbors $K$;

5: prior guess of the parameter vector $\mathbf{x}_b$ and prior parameter covariance matrix $\mathbf{P}_b$.

6: Initialisation step:

7: **for** $i = 1$ to $N_{ens}$ **do**

8: $\mathbf{x}_{(0)}^{(i)} = \mathbf{x}_b + \varepsilon^b$ with, $\varepsilon^b \sim \mathcal{N}(0, \mathbf{P}_b)$

9: **for** $k = 1$ to $T$ **do**

10: **Forecast step:**

11: **for** $i = 1$ to $N_{ens}$ **do**

12: $$ \mathbf{x}_{(k)}^{(i)} = \mathcal{A}_{(k-1,k)}(\mathbf{x}_{(k-1)}^{(i)}) + \varepsilon_{(k)}^{(i)}, \text{ where,} $$ 

**Locally-linear analog operator:** $\mathcal{A}_{(k-1,k)}(\mathbf{x}_{(k-1)}^{(i)}) := \mu_{(k)}^{\text{LL}} = \alpha \mathbf{x}_{(k-1)} + \beta$ 

and $\varepsilon_{(k)}^{(i)} \sim \Sigma_{(k)}^{\text{LL}}$

where $(\mathbf{a}_n, \mathbf{s}_n)_{n \in \mathcal{I}}$ (with $\mathcal{I} = \{i_1, \dots, i_K\}$) are the K-pairs of analog and successor for the i-th analysis member of the ensemble at iteration $k-1$ and $\text{cov}_\omega$ is the weighted covariance.

13: **Update step:**

$$ \mathbf{P}_{(k)}^{f} \mathbf{H}_{(k)}^{T} = \frac{1}{N_{ens}-1} \sum_{i=1}^{N_{ens}} (\mathbf{x}_{(k)}^{(i)} - \bar{\mathbf{x}}_{(k)}^{f}) (\mathcal{H}_{(k)}(\mathbf{x}_{(k)}^{(i)}, \mathbf{v}_{(k)}) - \mathcal{H}_{(k)}(\bar{\mathbf{x}}_{(k)}^{f}, \mathbf{v}_{(k)}))^{T} \\ \mathbf{H}_{(k)} \mathbf{P}_{(k)}^{f} \mathbf{H}_{(k)}^{T} = \frac{1}{N_{ens}-1} \sum_{i=1}^{N_{ens}} (\mathcal{H}_{(k)}(\mathbf{x}_{(k)}^{(i)}, \mathbf{v}_{(k)}) - \mathcal{H}_{(k)}(\bar{\mathbf{x}}_{(k)}^{f}, \mathbf{v}_{(k)})) (\mathcal{H}_{(k)}(\mathbf{x}_{(k)}^{(i)}, \mathbf{v}_{(k)}) - \mathcal{H}_{(k)}(\bar{\mathbf{x}}_{(k)}^{f}, \mathbf{v}_{(k)}))^{T} \\ \mathbf{K}_{(k)} = \mathbf{P}_{(k)}^{f} \mathbf{H}_{(k)}^{T} (\mathbf{R}_{(k)} + \mathbf{H}_{(k)} \mathbf{P}_{(k)}^{f} \mathbf{H}_{(k)}^{T})^{-1} $$

with $\bar{\mathbf{x}}_{(k)}^f$ the mean of the forecast members of the ensemble

14: **for** $i = 1$ to $N_{ens}$ **do**

15: $$ \mathbf{y}_{(k)}^{(i)} = \mathbf{y}_{(k)} + \mathbf{e}_{(k)}^{(i)} \text{ with, } \mathbf{e}_{(k)}^{(i)} \sim \mathcal{N}(0, \mathbf{R}_{(k)}) $$ 
 $$ \mathbf{x}_{(k)}^{(i)} = \mathbf{x}_{(k)}^{(i)} + \mathbf{K}_{(k)} (\mathbf{y}_{(k)}^{(i)} - \mathcal{H}_{(k)}(\mathbf{x}_{(k)}^{(i)}, \mathbf{v}_{(k)})) $$

## IV. NUMERICAL RESULTS

In this section, the numerical results of the proposed methodology to quantify and reduce the uncertainties based on global sensitivity analysis and a data-driven data assimilation approach are presented in the context of an industrial operating wind turbine. The two categories of param-
eters investigated in this application are the wind turbine model properties and the wind-inflow conditions. In the sensitivity analysis of the fatigue loads of the wind turbine, we assume that the 10-minute mean and standard deviation obtained from the SCADA are respectively equal to 10 m/s and 1.4 m/s.

## A. Case description

For the purpose of this work, the considered model is a numerical representation of a reference 2MW onshore horizontal-axis wind turbine based on the open-source aero-servo-elastic software FAST developed by the National Renewable Energy Laboratory (NREL) (Jonkman, Buhl Jr. et al., 2005). This numerical code employs a combined modal and multibody dynamics formulation which allows one to consider a limited degree of freedom number for the structure. Moreover, the aerodynamic model relies on the blade-element momentum theory coupled with some corrections, e.g., dynamic stall. The generation of the synthetic turbulent wind field solicitation uses a Kaimal turbulence model with an exponential spatial coherence method thanks to the TurbSim software (Jonkman, 2009). Some specifications of the turbine are presented in Table IV.

<table><caption>Table IV. Reference wind turbine specifications.</caption><thead><tr><td>Quantity</td><td>Value</td></tr></thead><tbody><tr><td>Number of blades</td><td>3</td></tr><tr><td>Rated power</td><td>2.0 MW</td></tr><tr><td>Rotor speed range</td><td>8.5 – 17.1 rpm (±16 %)</td></tr><tr><td>Rated wind speed</td><td>13 m/s</td></tr><tr><td>Cut-in wind speed</td><td>3.0 m/s</td></tr><tr><td>Cut-out wind speed</td><td>25 m/s</td></tr><tr><td>Rotor radius</td><td>41 m</td></tr><tr><td>Hub height</td><td>80 m</td></tr></tbody></table>

The in situ data used to assess the performances of our procedure are based on a specific measurement campaign of eight months from the French national project SMARTEOLE. For that purpose, the wind turbine has a SCADA gathering 10-minute statistics about the external conditions at the level of the nacelle hub, e.g., wind speed or direction, and also information on the turbine operation, e.g., generator speed, generated power. Alongside, a Light Detection And Ranging (LI-DAR) system is placed on top of the wind turbine nacelle in order to measure the upstream wind flow conditions. A graphical representation of the monitoring system configuration is proposed
in Figure 3. In the study, we suppose that the wind speed at hub height reconstructed from the LIDAR system is the free wind to be applied on the servo-aero-elastic model through the synthetic turbulence wind field. Lastly, bi-axial measuring devices are located at mid and top tower height position. From these sensors, we can record four functional acceleration time series. Then, the power spectral density (PSD) of each measured acceleration time series is computed using Welch's method (Welch, 1967).

<figure><img src="image_4.png" /><figcaption>Figure 3. Monitoring system configuration for the reference wind turbine.</figcaption></figure>

365

366

367

## B. Global sensitivity analysis on fatigue loads

368

369

370

371

372

373

374

375

376

377

378

379

To quantify the importance of each input parameter on the variability of the fatigue loads obtained from the aero-servo-elastic numerical model, a global sensitivity analysis (GSA) based on Sobol' index estimation has been investigated. We focus our interest on total Sobol' sensitivity indices (Sobol', 1990). The total Sobol' index associated with each input parameter represents the amount of variance due to the quantity of interest alone or in interaction with any other subset of parameters. It allows one to quantify the part of variation in the damage equivalent load that could be reduced if the parameter was to be fixed to a single value. To alleviate the computational cost in the sensitivity index estimation, heteroscedastic Gaussian process (GP) models (Ginsbourger et al., 2008) are built independently for each Damage Equivalent Load (DEL). The notion of damage equivalent load (DEL) (Veldkamp, 2006) is often used and is defined as a virtual stress amplitude that would create the same fatigue damage as a particular load history considering a specific number of regular cycles. These short-term fatigue estimations are computed from the load time
series, obtained from the aero-servo-elastic model, by applying a rain-flow counting algorithm and the Palmgren-Miner rule for linear fatigue damage accumulation (Hansen, 2015). Fitting such surrogate model to the load behavior of a wind turbine requires a design of experiments covering the range of variation in all parameters. In that context, we rely on a Latin Hypercube Sampling (LHS) with a geometrical criterion maximizing the minimum distance between the design points (i.e., the sample points obtained from the LHS) (Damblin, Couplet, and Iooss, 2013). To testify the accuracy of the fitted surrogate model for each output of interest described in Table V, an augmented LHS of size 200 has been generated. Then, ten different turbulent inflow realizations are generated using the Kaimal spectrum with an exponential spatial coherence model for each point of the DOE, from which the empirical mean and standard deviation of the fatigue loads are estimated. The heteroscedastic property of the GP, as described in (Hirvoas et al., 2021), allows one to capture the global fatigue behavior of the turbine but also to estimate the inherent variability due to different turbulent wind field realizations. This study leads to a total number of 11,960 aero-servo-elastic numerical model evaluations.

Eight different model quantities of interest are considered for describing the fatigue behavior of the wind turbine, see Table V. For each output, the total effect Sobol indices are estimated using the corresponding heteroscedastic Gaussian process metamodel based on the estimator proposed by Jansen (1999) and implemented in the function sobolGP of the R package sensitivity (Iooss et al., 2019). The prediction performance of each Gaussian process is quantified thanks to a validation set, not seen during the training phase, and the predictivity coefficient $Q^2$ (Marrel et al., 2008). In the presented application, the surrogate model have $Q^2$ factors over at least 0.8. The GSA estimation approach relies on the complete conditional predictive distribution of the metamodel which allows one to evaluate the uncertainty in the estimation due to the Monte Carlo procedure and the surrogate approximation, see (Hirvoas et al., 2021).
<table><caption>Table V. Wind turbine model fatigue load outputs with their Wöhler exponent *m* (Meyers and Chawla, 2008).</caption><thead><tr><th>Quantity of interest</th><th>m</th></tr></thead><tbody><tr><td>DEL blade root in-plane bending moment</td><td>10</td></tr><tr><td>DEL blade root out-of-plane bending moment</td><td>10</td></tr><tr><td>DEL tower bottom fore-aft bending moment</td><td>3</td></tr><tr><td>DEL tower bottom side-to-side bending moment</td><td>3</td></tr><tr><td>DEL tower top side-to-side bending moment</td><td>3</td></tr><tr><td>DEL tower top fore-aft force</td><td>3</td></tr><tr><td>DEL shaft torsional moment</td><td>3</td></tr></tbody></table>

For the estimation procedure, two distinct LHSs with a maximin criterion have been generated. The uncertainty related to the kriging approximation is quantified by using 100 samples from the conditional distribution of the predictor based on the learning sample. Moreover, the uncertainty due to Monte Carlo integration was estimated with a bootstrap procedure with a sample size of 100, see Efron (1981) for further details in bootstrapping strategy. The estimated total Sobol' indices for the considered quantities of interest with their corresponding 95% confidence intervals are presented in Figure 4. Most of the outputs have a large total Sobol' index for the errors relative to the wind speed $\Delta\bar{u}$ and $\Delta\sigma_u$. These input parameters have an important impact on the variability of fatigue loads obtained from our aero-servo-elastic numerical model. The vertical wind shear coefficient $\alpha$ has also a clear impact in particular for the torsional moment of the shaft and the out-of-plane bending moment of the blade. The noticeable effect of the wind shear for rotating components can be explained by the fact that they will face cyclic changes in wind velocity if wind shear is considered. Eight other parameters describing the wind inflow conditions or the wind turbine model properties have total Sobol' indices higher to the arbitrary threshold (set to 5.0e-02) and can be considered as influential. The arbitrary threshold is used to discriminate efficiently sensitive and insensitive input parameters. For clarity, these parameters are underlined in Figure 4. In particular, we can notice that model property parameters related to tower thickness, lineic mass and mass imbalance related to the blades (e, $\alpha_{mass}$, and $\eta_B$) have a non-negligible influence on fatigue load variance of the considered wind turbine components. The remaining parameters can be fixed to any specific value in their range of variability without affecting the considered fatigue loads. In our research work, these parameters having a negligible influence on fatigue loads are fixed to the mean of their statistical distribution, see Table II and Table III.

After assessing the sensitivity analysis of the fatigue load of some critical components of the
wind turbine structure, one major challenge is to reduce the uncertainties affecting the most influ-
ential input parameters.

<figure><img src="image_5.png" /></figure>

Figure 4. Estimation of total Sobol' indices (y-axis) with their 95% confidence interval corresponding to each of the 20 parameters (x-axis) for the different fatigue loads. The dashed line corresponds to a threshold arbitrarily chosen to $5.0e-2$. <u>Underlined total Sobol' indices represent indices higher to the threshold value.</u> Confidence intervals (CI) are obtained by taking into account the uncertainties due to both the metamodel and the Monte Carlo estimation. The number of samples for the conditional Gaussian process, in order to quantify the uncertainty of the kriging approximation, was set to 100. The uncertainty due to Monte Carlo integration was computed with a bootstrap procedure with a sample size of 100.

C. Identifiability study

A major issue for parameter estimation problem is the identifiability. In this context, Dobre
*et al.* (2012) highlights that nullity of total sensitivity index for a specific input parameter implies
its non-identifiability from the measured output used during the recursive inference procedure.
Consequently, we perform a GSA on the measured outputs in order to determine which parameters
cannot be inferred with the current sensors on the wind turbine. In our industrial application, six
measured outputs are considered, see Table VI.

For the acceleration outputs, we are mainly interested in their response in the frequency-domain
by using the power spectral density (PSD). When performing GSA, discretized PSD series involve a substantial dimensionality and a high degree of redundancy. To overcome this issue, the different discretized PSD outputs have been reduced using a Principal Component Analysis (PCA) (Wold, Esbensen, and Geladi, 1987). This dimensionality reduction approach allows the functional output expansion in a new reduced space spanned by the most significant directions in terms of variance. Then, a method based on PCA and GSA with a GP model is used to compute an aggregated Sobol' index for each input parameter of the model. Indeed, the multivariate sensitivity analysis decomposes the discretized output into the non-correlated principal components and computes sensitivity indices on each principal component. Finally, the overall effect of each parameter can be summarized by the aggregated Sobol' index equal to the weighted sum of the indices over the principal component, with weights proportional to the inertia associated with the components (Lamboni, Monod, and Makowski, 2011).

The proposed index synthesizes the influence of the parameter on the whole discretized functional output. Table VII summarizes the estimated Sobol' indices for the scalar 10-minutes average observations and aggregated Sobol' indices for the discretized PSD series. In this sensitivity analysis, the input parameters having total Sobol' index values under a threshold set at $1e-02$ are considered as non-identifiable from the measured output.

<table><caption>Table VI. Observations from our reference wind turbine used in the data assimilation procedure.</caption><thead><tr><td>Observation</td><td>Unit</td></tr></thead><tbody><tr><td>10-minute mean power production</td><td>[kW]</td></tr><tr><td>10-minute mean rotor speed</td><td>[rpm]</td></tr><tr><td>Tower middle fore-aft acceleration's PSD</td><td>[dB]</td></tr><tr><td>Tower middle side-to-side acceleration's PSD</td><td>[dB]</td></tr><tr><td>Tower top fore-aft acceleration's PSD</td><td>[dB]</td></tr><tr><td>Tower top side to side acceleration's PSD</td><td>[dB]</td></tr></tbody></table>

<footnote>457</footnote>

<footnote>458</footnote>
<table><caption>Table VII. Total Sobol' and aggregated Sobol' indices for each output used during the recursive inference procedure. Due to the GSA on fatigue loads, the eight wind parameters have been whittled down to six, and the twelve model parameters have been whittled down to three. Estimated Sobol' indices higher than the arbitrary threshold, set at 1e-02, are underlined.</caption><thead><tr><th colspan="2">Measured output</th><th>Δμ̄<br>[m/s]</th><th>Δσ<sub>u</sub><br>[m/s]</th><th>φ<sub>h</sub><br>[°]</th><th>Λ<sub>u</sub><br>[m]</th><th>a<br>[-]</th><th>α<br>[-]</th><th>e<br>[%]</th><th>α<sub>mass</sub><br>[%]</th><th>η<sub>B</sub><br>[%]</th></tr></thead><tbody><tr><td rowspan="2">Sobol' index</td><td>10-minute mean power production</td><td>9.81e-01</td><td>4.29e-04</td><td>1.71e-02</td><td>1.30e-04</td><td>3.70e-04</td><td>1.50e-02</td><td>3.84e-05</td><td>3.83e-04</td><td>5.23e-05</td></tr><tr><td>10-minute mean rotor speed</td><td>9.75e-01</td><td>3.30e-03</td><td>1.87e-02</td><td>9.43e-04</td><td>1.61e-03</td><td>1.62e-02</td><td>1.03e-04</td><td>7.56e-04</td><td>7.34e-05</td></tr><tr><td rowspan="4">Aggregated Sobol' index</td><td>Tower middle fore-aft acceleration's PSD</td><td>1.44e-01</td><td>2.49e-01</td><td>1.00e-02</td><td>1.77e-01</td><td>3.70e-01</td><td>1.33e-02</td><td>4.58e-02</td><td>5.82e-03</td><td>3.48e-03</td></tr><tr><td>Tower middle side-to-side acceleration's PSD</td><td>2.04e-01</td><td>2.51e-01</td><td>1.09e-02</td><td>1.92e-01</td><td>3.00e-01</td><td>1.33e-02</td><td>4.49e-02</td><td>4.86e-03</td><td>3.42e-03</td></tr><tr><td>Tower top fore-aft acceleration's PSD</td><td>3.12e-01</td><td>2.16e-01</td><td>1.87e-02</td><td>1.75e-01</td><td>2.69e-01</td><td>9.59e-03</td><td>3.36e-02</td><td>8.49e-03</td><td>7.01e-03</td></tr><tr><td>Tower top side to side acceleration's PSD</td><td>2.84e-01</td><td>1.87e-01</td><td>1.18e-02</td><td>1.76e-01</td><td>2.50e-01</td><td>1.21e-02</td><td>8.33e-02</td><td>5.52e-03</td><td>2.38e-02</td></tr></tbody></table>

According to the GSA, the coefficient related to the blade mass coefficient $\alpha_{mass}$ is not identifiable with the current observations. Consequently, the model parameter properties remaining for the inference procedure are the tower thickness coefficient $e$, and the mass imbalance factor $\eta_B$. Moreover, all the influent parameters related to the wind field remain candidates for the recursive inference strategy.

### D. Recursive inference strategy based on AnEnKF approach

The current in situ wind data availability or quality from the LIDAR system does not allow a proper extraction of the mean flow angle $\phi_h$, the longitudinal turbulence length scale $\Lambda_u$, and the decrement parameter of the coherence model $a$. Consequently, only the five remaining parameters having an influential effect on the fatigue behavior of the structure and potentially identifiable are considered during the recursive inference procedure. These input parameters and their corresponding prior Gaussian distributions are detailed in Table VIII. Their corresponding reference variable in the augmented state vector is also specified.

For assessing the performance of the AnEnKF for our recursive inference procedure, we rely on pseudo-experimental numerical tests. They consist in performing forward aero-servo-elastic simulations considering known values of the input parameters, and then adding a Gaussian noise of known variance to the simulated measurements. In our study, the simulated data are perturbed by considering a covariance matrix such as the obtained standard deviation is equivalent to a 10% signal-to-noise ratio. The pseudo-simulated responses of the wind turbine structure are generated
<table><caption>Table VIII. A-priori Gaussian distribution $\mathcal{G}$ for each of the considered input parameters.</caption><thead><tr><th>Input parameter</th><th>Variable</th><th>Distribution</th><th>Initial prior</th><th>State</th></tr></thead><tbody><tr><td>Tower thickness</td><td>e</td><td>G</td><td>μ = 0.00 σ = 7.00e-01</td><td rowspan="2">x<sup>1</sup></td></tr><tr><td>Blade mass imbalance</td><td>η<sub>B</sub></td><td>G</td><td>μ = 2.50e-02 σ = 8.33e-03</td></tr><tr><td>Error mean of the wind speed at hub height</td><td>Δu&#773;</td><td>G</td><td>μ = 0.00 σ = 9.11e-01</td><td rowspan="2">x<sup>2</sup></td></tr><tr><td>Error standard deviation of the wind speed at hub height</td><td>Δσ<sub>u</sub></td><td>G</td><td>μ = 0.00 σ = 9.70e-02</td></tr><tr><td>Vertical wind shear exponent</td><td>α</td><td>G</td><td>μ = 1.30e-01 σ = 2.90e-01</td><td></td></tr></tbody></table>

using the wind inflow conditions obtained from the nacelle mounted LIDAR for a specific day and the mean values of the model properties described in Table III. The noisy pseudo-experimental outputs used to recursively update the wind turbine model are 10-minute mean power production and rotor speed, and the PSD of the acceleration time series obtained for side to side and fore-aft at the two different tower positions. Our recursive inference problem using a filtering-based estimation procedure can be considered as a state estimation problem for the following augmented system:

$$ \forall k \in \mathbb{N}^*, \left\{ \begin{array}{l} \mathbf{x}_{(k)} = \begin{bmatrix} \mathbf{x}_{(k)}^1 \\ \mathbf{x}_{(k)}^2 \end{bmatrix} = \begin{pmatrix} \mathbf{x}_{(k-1)}^1 \\ \mathcal{A}_{(k-1,k)}(\mathbf{x}_{(k-1)}^2) \end{pmatrix} + \boldsymbol{\varepsilon}_{(k)}^m . \\ \mathbf{y}_{(k)} = \mathcal{H}_{(k)}(\mathbf{x}_{(k)}, \mathbf{v}_{(k)}) + \boldsymbol{\varepsilon}_{(k)}^o \end{array} \right. $$

where $\mathbf{x}_{(k-1)}^1$ and $\mathbf{x}_{(k-1)}^2$ are respectively the whittled-down set of five uncertain parameters for the model properties and the wind inflow conditions at iteration $k-1$ as described in Table VIII, $\mathcal{A}_{(k-1,k)}$ is the analog forecasting operator as detailed in Section III, $\mathbf{v}_{(k)}$ is the forcing vector corresponding to the 10-minute mean and standard deviation wind speed obtained from the SCADA system, and $\mathcal{H}_{(k)}$ is the combination of the aero-servo-elastic model FAST and the turbulent wind field generation software TurbSim.

For the initialization of the EnKF approach, independent Gaussian distributions are assumed to be the initial prior for each of the input parameters, see Table VIII. The initial error covariance matrix of the input parameters, denoted by $P_b$, is thus assumed to be diagonal. To create the catalog, we rely on the measurements obtained from both the SCADA system and the LIDAR installed on the onshore wind turbine. A data pretreatment has been performed in order to find any corrupted observations. The obtained database consists in both 4,735 analog situations to be
compared to the current parameters related to the wind inflow and their corresponding successors at a 10-minute interval.

Figure 5 shows the results of the identification of the considered input parameters by applying the AnEnKF approach with the locally-linear forecasting operator using $N = 500$ members and $K = 50$ nearest-neighbors. It can be noticed that the augmented state vector is well reconstructed by using this non-parametric data assimilation procedure which allows to emulate the dynamical model from a dataset. Indeed, the mean of the empirical distribution obtained from the members of the ensemble is close to the true hidden-state for every parameter. A major advantage of the procedure is the confidence intervals obtained at each inference iteration allowing us to give information about the difficulty to retrieve the value of the input parameters from the measured outputs.

We can notice that the uncertainties, represented by the ensemble, have drastically reduced trough the iterations. The confidence intervals take into account the errors related to the model and the observations.

<figure><img src="image_6.png" /></figure>

Figure 5. Iteration evolution of the posteriori estimates of the input parameters. Results obtained by running the AnEnKF procedure with $N = 500$ members of the ensemble used for the estimation and considering pseudo-experimental numerical observations.
512 V. CONCLUSION

In the present work, we extend a procedure to quantify and reduce the uncertainties affecting the fatigue load estimation of a wind turbine numerical model. The fatigue loads encountered by a wind turbine structure are function of the parameters describing the turbulent wind field, the structural properties, and the control system. The study aims at taking into account these parameters used as input to aero-servo-elastic fatigue load simulations of an operating wind turbine. The procedure relies on a global sensitivity analysis and a recursive Bayesian inference method. A major challenge during the recursive inference procedure is the dynamic behavior of the inflow-related parameters. Unfortunately, the underlying dynamic behavior of these parameters is not explicitly known. To overcome this issue, a combination of the implicit analog forecasting of the dynamics with the ensemble Kalman filtering scheme is investigated.

Finally, we demonstrate the applicability and performance of the procedure using a numerical representation of a reference wind turbine. The study leads to the following main conclusions. The global sensitivity analysis based on heteroscedastic Gaussian processes for the estimation of Sobol' indices shows that parameters related both to the wind and the structure have an influence on the fatigue loads of a wind turbine structure. The presented metamodeling approach is an efficient way to capture the inherent stochasticity of aero-servo-elastic simulations due to the turbulent inflow realization leading to variations in the quantities of interest. After determining the most influential parameters in terms of fatigue loads variability, an identifiability study based on a global sensitivity analysis is performed to assess if these parameters can be inferred from the current sensors. The sensitivity analysis is based on the estimation of the so-called aggregated Sobol' indices involving a principal component analysis in order to take into account the functional behavior of the measured outputs. Finally, the ensemble Kalman filtering method coupled with the analog forecasting strategy used in this study is very suitable for carrying the recursive inference of parameters related to the wind field solicitation and the wind turbine numerical description.

Further research should focus on the quality of the catalog used for the analog forecasting strategy. Additionally, other types of kernels in the forecasting operator have to be studied. Lastly, the hyperparameters used in the K-nearest neighbors method and the chosen kernel function could be optimized for each member of the ensemble Kalman filtering procedure by using a cross-validation approach. From an industrial perspective, the proposed AnEnKF methodology has to be performed using measured acceleration time-series obtained from the sensor devices of the onshore wind turbine.
bine.

544 ACKNOWLEDGMENTS

545 The author acknowledge SMARTEOLE project partners for the use of experimental data from
546 national project SMARTEOLE (ANR-14-CE05-0034) measurement campaigns.
REFERENCES

Aanonsen, S. I., Nævdal, G., Oliver, D. S., Reynolds, A. C., Vallès, B., *et al.*, "The ensemble kalman filter in reservoir engineering-a review," *Spe Journal* **14*, 393-412 (2009).

Alexander, R., Zhao, Z., Székely, E., and Giannakis, D., "Kernel analog forecasting of tropical intraseasonal oscillations," *Journal of the Atmospheric Sciences* **74**, 1321-1342 (2017).

Ayet, A. and Tandeo, P., "Nowcasting solar irradiance using an analog method and geostationary satellite images," *Solar Energy* **164**, 301-315 (2018).

Bertino, L., Evensen, G., and Wackernagel, H., "Sequential data assimilation techniques in oceanography," *International Statistical Review* **71**, 223-241 (2003).

Bocquet, M., Farchi, A., and Malartic, Q., "Online learning of both state and dynamics using ensemble kalman filters," *arXiv preprint arXiv:2006.03859* (2020).

Brown, S. D., "The kalman filter in analytical chemistry," *Analytica chimica acta* **181**, 1-26 (1986).

Brunton, S. L., Proctor, J. L., and Kutz, J. N., "Discovering governing equations from data by sparse identification of nonlinear dynamical systems," *Proceedings of the national academy of sciences* **113**, 3932-3937 (2016).

Chau, T. T. T., Ailliot, P., and Monbet, V., "An algorithm for non-parametric estimation in state-space models," *Computational Statistics & Data Analysis* **153**, 107062 (2021).

Chen, L., "Curse of dimensionality," in *Encyclopedia of Database Systems*, edited by L. LIU and M. T. ÖZSU (Springer US, Boston, MA, 2009) pp. 545-546.

Damblin, G., Couplet, M., and Iooss, B., "Numerical studies of space-filling designs: optimization of latin hypercube samples and subprojection properties," *Journal of Simulation* **7**, 276-289 (2013).

Dimitrov, N., Natarajan, A., and Kelly, M., "Model of wind shear conditional on turbulence and its impact on wind turbine loads," *Wind Energy* **18**, 1917-1931 (2015).

Dimitrov, N., Natarajan, A., and Mann, J., "Effects of normal and extreme turbulence spectral parameters on wind turbine loads," *Renewable Energy* **101**, 1180-1193 (2017).

Dobre, S., Bastogne, T., Profeta, C., Barberi-Heyob, M., and Richard, A., "Limits of variance-based sensitivity analysis for non-identifiability testing in high dimensional dynamic models," *Automatica* **48**, 2740-2749 (2012).

Durbin, J. and Koopman, S. J., *Time series analysis by state space methods* (Oxford university press, 2012).
Efron, B., "Nonparametric estimates of standard error: the jackknife, the bootstrap and other methods," Biometrika **68**, 589-599 (1981).

Evensen, G., "Sequential data assimilation with a nonlinear quasi-geostrophic model using monte carlo methods to forecast error statistics," Journal of Geophysical Research: Oceans **99**, 10143-10162 (1994).

Evensen, G., "The ensemble kalman filter: Theoretical formulation and practical implementation," Ocean dynamics **53**, 343-367 (2003).

Evensen, G., *Data assimilation: the ensemble Kalman filter* (Springer Science and Business Media, 2009).

Fablet, R., Viet, P., Lguensat, R., and Chapron, B., "Data-driven assimilation of irregularly-sampled image time series," in *2017 IEEE International Conference on Image Processing (ICIP)* (IEEE, 2017) pp. 4302-4306.

Friedman, J. H., "On bias, variance, 0/1-loss, and the curse-of-dimensionality," Data mining and knowledge discovery **1**, 55-77 (1997).

Ginsbourger, D., Picheny, V., Roustant, O., and Richet, Y., "Kriging with heterogeneous nugget effect for the approximation of noisy simulators with tunable fidelity (krigeage avec effet de pépite hétérogène pour l'approximation de simulateurs bruités à fidélité réglable)," (2008).

Hamilton, F., Berry, T., and Sauer, T., "Ensemble kalman filtering without a model," Physical Review X **6**, 011021 (2016).

Hansen, M., *Aerodynamics of wind turbines* (Routledge, 2015).

Harvey, A. C., *Forecasting, structural time series models and the Kalman filter* (Cambridge university press, 1990).

Haykin, S., *Kalman filtering and neural networks*, Vol. 47 (John Wiley & Sons, 2004).

Hirvoas, A., Prieur, C., Arnaud, E., Caleyron, F., and Munoz Zuniga, M., "Quantification and reduction of uncertainties in a wind turbine numerical model based on a global sensitivity analysis and a recursive bayesian inference approach," International Journal for Numerical Methods in Engineering (2021), https://doi.org/10.1002/nme.6630, https://onlinelibrary.wiley.com/doi/pdf/10.1002/nme.6630.

Holierhoek, J., Korterink, H., van de Pieterman, R., Rademakers, L., and Lekou, D., "Recommended Practices for Measuring in Situ the 'Loads' on Drive Train, Pitch System and Yaw System." Energy Research Center of the Netherlands (ECN) (2010).

Houtekamer, P. L., He, B., and Mitchell, H. L., "Parallel implementation of an ensemble kalman
filter,” Monthly Weather Review **142**, 1163–1182 (2014).

Houtekamer, P. L. and Mitchell, H. L., “A sequential ensemble kalman filter for atmospheric data assimilation,” Monthly Weather Review **129**, 123–137 (2001).

IEC, I. E. C., “Iec 61400-1: 2019: Wind energy generation systems-part 1: Design requirements,” (2019).

Iooss, B., Janon, A., Pujol, G., with contributions from Baptiste Broto,, Boumhaout, K., Veiga, S. D., Delage, T., Fruth, J., Gilquin, L., Guillaume, J., Le Gratiet, L., Lemaitre, P., Nelson, B. L., Monari, F., Oomen, R., Rakovec, O., Ramos, B., Roustant, O., Song, E., Staum, J., Sueur, R., Touati, T., and Weber, F., sensitivity: Global Sensitivity Analysis of Model Outputs (2019), r package version 1.16.0.

Jansen, M. J., “Analysis of variance designs for model output,” Computer Physics Communications **117**, 35–43 (1999).

Jonkman, B. J., “Turbsim user’s guide: Version 1.50,” Tech. Rep. (National Renewable Energy Lab (NREL), Golden, CO, USA, 2009).

Jonkman, J. M. , Buhl Jr., M. L., et al., “Fast user’s guide,” Golden, CO: National Renewable Energy Laboratory **365**, 366 (2005).

Julier, S. J. and Uhlmann, J. K., “New extension of the kalman filter to nonlinear systems,” in *Signal processing, sensor fusion, and target recognition VI*, Vol. 3068 (International Society for Optics and Photonics, 1997) pp. 182–193.

Kaimal, J. C., Wyngaard, J. C., Izumi, Y., and Coté, O., “Spectral characteristics of surface-layer turbulence,” Quarterly Journal of the Royal Meteorological Society **98**, 563–589 (1972).

Kalman, R. E., “A new approach to linear filtering and prediction problems,” (1960).

Koukoura, C., *Validated Loads Prediction Models for Offshore Wind Turbines for Enhanced Component Reliability*, Ph.D. thesis, Technical University of Denmark (2014).

Lamboni, M., Monod, H., and Makowski, D., “Multivariate sensitivity analysis to measure global contribution of input factors in dynamic models,” Reliability Engineering and System Safety **96**, 450–459 (2011).

Lguensat, R., Tandeo, P., Ailliot, P., Pulido, M., and Fablet, R., “The analog data assimilation,” Monthly Weather Review **145**, 4093–4107 (2017).

Lorenz, E. N., “Atmospheric predictability as revealed by naturally occurring analogues,” Journal of the Atmospheric Sciences **26**, 636–646 (1969).

Marrel, A., Iooss, B., Van Dorpe, F., and Volkova, E., “An efficient methodology for modeling
complex computer codes with gaussian processes,” Computational Statistics and Data Analysis **52**, 4731–4744 (2008).

Meyers, M. A. and Chawla, K. K., *Mechanical behavior of materials* (Cambridge university press, 2008).

Pathak, J., Hunt, B., Girvan, M., Lu, Z., and Ott, E., “Model-free prediction of large spatiotemporally chaotic systems from data: A reservoir computing approach,” *Physical Review Letters* **120** (2018), 10.1103/PhysRevLett.120.024102.

Peterson, L. E., “K-nearest neighbor,” *Scholarpedia* **4**, 1883 (2009).

Robertson, A. N., Shaler, K., Sethuraman, L., and Jonkman, J. M., “Sensitivity analysis of the effect of wind characteristics and turbine properties on wind turbine loads,” *Wind Energy Science (Online)* **4** (2019a).

Robertson, A. N., Shaler, K., Sethuraman, L., and Jonkman, J. M., “Sensitivity of uncertainty in wind characteristics and wind turbine properties on wind turbine extreme and fatigue loads,” *Wind Energy Science Discussions*, 1–41 (2019b).

Saranyasoontorn, K., Manuel, L., and Veers, P., “On estimation of coherence in inflow turbulence based on field measurements,” in *42nd AIAA Aerospace Sciences Meeting and Exhibit* (2004) p. 1002.

Simms, D., Schreck, S., Hand, M., and Fingersh, L. J., “Nrel unsteady aerodynamics experiment in the nasa-ames wind tunnel: a comparison of predictions to measurements,” *Tech. Rep.* (National Renewable Energy Lab., Golden, CO (US), 2001).

Snyder, C. and Zhang, F., “Assimilation of simulated doppler radar observations with an ensemble kalman filter.” *Monthly Weather Review* **131** (2003).

Sobol', I. M., “On sensitivity estimation for nonlinear mathematical models,” *Matematicheskoe modelirovanie* **2**, 112–118 (1990).

Solari, G. and Piccardo, G., “Probabilistic 3-d turbulence modeling for gust buffeting of structures,” *Probabilistic Engineering Mechanics* **16**, 73–86 (2001).

Tandeo, P., Ailliot, P., Fablet, R., Ruiz, J., Rousseau, F., and Chapron, B., “The analog ensemble kalman filter and smoother,” (2014).

Tandeo, P., Ailliot, P., Ruiz, J., Hannart, A., Chapron, B., Cuzol, A., Monbet, V., Easton, R., and Fablet, R., “Combining analog method and ensemble data assimilation: application to the lorenz-63 chaotic system,” in *Machine learning and data mining approaches to climate science* (Springer, 2015) pp. 3–12.
Toth, Z., “Long-range weather forecasting using an analog approach,” Journal of climate **2**, 594–607 (1989).

Veldkamp, H. F., “Chances in wind energy: a probabilistic approach to wind turbine fatigue design,” (2006).

Welch, P., “The use of fast fourier transform for the estimation of power spectra: a method based on time averaging over short, modified periodograms,” IEEE Transactions on audio and electroacoustics **15**, 70–73 (1967).

Wells, C., *The Kalman filter in finance*, Vol. 32 (Springer Science and Business Media, 2013).

Witcher, D., “Uncertainty Quantification Techniques in Wind Turbine,” (2017).

Wold, S., Esbensen, K., and Geladi, P., “Principal component analysis,” Chemometrics and intelligent laboratory systems **2**, 37–52 (1987).