Sensitivity and Uncertainty of the FLORIS Model Applied on the Lillgrund Wind Farm

Abstract

Wind farms experience significant efficiency losses due to the aerodynamic interaction between turbines. A possible control technique to minimize these losses is yaw-based wake steering. This paper investigates the potential for improved performance of the Lillgrund wind farm through a detailed calibration of a low-fidelity engineering model aimed specifically at yaw-based wake steering. The importance of each model parameter is assessed through a sensitivity analysis. This work shows that the model is overparameterized as at least one model parameter can be excluded from the calibration. The performance of the calibrated model is tested through an uncertainty analysis, which showed that the model has a significant bias but low uncertainty when comparing the predicted wake losses with measured wake losses. The model is used to optimize the annual energy production of the Lillgrund wind farm by determining yaw angles for specific inflow conditions. A significant energy gain is found when the optimal yaw angles are calculated deterministically. However, the energy gain decreases drastically when uncertainty in input conditions is included. More robust yaw angles can be obtained when the input uncertainty is taken into account during the optimization, which yields an energy gain of approximately 3.4%.

1. Introduction

Global wind power capacity has increased significantly from about 17 GW in 2000 to 623 GW in 2019 [1], and this trend is expected to continue as part of the sustainable energy transition. Currently, 95% of the 623 GW is from the onshore sector, which has increased by a factor of four in the past decade. In comparison, the installed capacity of offshore wind has increased by a factor of 13 [1], where numerous wind turbines are generally clustered in wind farms. However, clustering of turbines can result in significant reductions of the wind farm efficiency, e.g., Barthelmie et al. [2] show wake losses of 10–20%, due to the wake interaction of the wind turbines. The wakes behind wind turbines are characterized by a velocity deficit and increased turbulence compared to the free stream flow that the upstream turbine experiences [3]. The velocity deficit results in lower power production for the downstream wind turbines, while the increased turbulence generally leads to higher fatigue loads.

The adverse wake effects can potentially be mitigated through layout optimization [4] or various wind farm control methods. The oldest and also widely studied control technique is axial induction control [5]. The axial induction and thereby thrust of a given turbine can be reduced, which decreases the velocity deficit for a downstream turbine. This yields a lower electrical power production for the upstream turbine but potentially increase the electrical power production of the downstream turbine [6]. Induction control will generally lead to gains when applying low-fidelity models but have also been shown to result in power losses when higher-fidelity models are applied [7], while recent full-scale control test using a calibrated low fidelity model have shown potential benefit [8]. A second control technique is re-positioning of floating offshore wind turbines, which is still at an early stage of development, yet looks promising in terms of power gain [9,10,11].

Thirdly, wake effects can potentially be reduced by intentionally misaligning turbines with respect to the incoming wind direction [12,13]. This control technique is called yaw-based wake steering, where the yaw misalignment introduces a lateral force component, which results in a deflected wake. The misaligned wind turbine operates sub-optimal and produces less power than otherwise, but the downstream turbine may produce more power to compensate and even surpass the total power production as a result of the deflected wake. Numerous studies using low-fidelity models [14,15,16,17,18] and high-fidelity ones [9,19,20] have shown increases in the power production of wind farms. Furthermore, wind tunnel experiments reported positive gains [21,22], and even recent field tests have shown promising results [23,24].

A summary of the current state-of-the-art research on wind farm control is given in [10], with particular focus on the re-positioning, induction-based, and yaw based control strategies and the associated results derived from field tests, wind tunnel tests, as well as a wide range of low-, mid-, and high-fidelity model results. Kheirabadi and Nagamune [10] concluded that particular yaw based wind farm control has the potential to achieve increased power production. Recent investigations show the benefit and potential necessity of using dynamic control to improve wind farms operation [25,26,27,28,29]. However, despite the promising potential it is still paramount to diligently validate both the various control strategies and the models used before they can be used in full-scale applications [30].

Hence, the present work focuses on yaw-based wake steering and including the uncertainty and sensitivities in the process. This paper investigates the potential for improved performance of the Lillgrund wind farm, which has significant wake losses due to a very compact layout [31]. The low-fidelity engineering FLOw Redirection and Induction in Steady-state (FLORIS) model [32,33] aimed specifically at yaw-based wake steering is applied and calibrated. Low-fidelity models like FLORIS require calibration, which can be done using flow data from, e.g., lidar measurements [34] or high-fidelity simulations [14,35,36,37,38,39,40]. An alternative to calibrating low-fidelity engineering wake models is to build surrogate models based on high-fidelity simulations but requires numerous high-fidelity simulations [41].

In this research, a calibration method is applied using power measurements from Supervisory Control Data Acquisition (SCADA) data from the Lillgrund wind farm, Sweden. Furthermore, a sensitivity analysis is performed to assess the importance of each model parameter. Finally, the FLORIS model is used to estimate the potential Annual Energy Production (AEP) gain using wake steering and it is investigated whether the robustness of the optimized yaw angles can be increased by including uncertainty in the estimation.

The paper is organized as follows. Section 2 introduces the data set used for the calibration, the engineering model used to simulate the wind farm, and the calibration and optimization methods. The next section, Section 3, first introduces and discusses the sensitivity analysis and its results. Afterward, the accuracy of the model is assessed using an uncertainty analysis. Finally, the results of the yaw optimization are discussed. At last, the conclusions of this study are described in Section 4.

2. Materials and Methods

2.1. Data Set

The analysis is performed on the Lillgrund wind farm, which contains 48 Siemens SWT-2.3-93 turbines [42]. The turbines have a rotor diameter of D = 92.6 m, a hub height of 65 m, and a rated power of 2.3 MW. The corresponding rated wind speed is around 12 m/s. The wind farm layout is shown in Figure 1a and it can be seen in Figure 1 that the farm is densely packed with the smallest turbine spacing being 3.3D. All aligned wind directions ${\theta }_{\infty }$ and the corresponding turbine spacings are summarized in

Table 1. The alignment wind directions per turbine spacing.

Spacing (D)	${\mathit{\theta }}_{\mathit{\infty }}$ (deg)
3.3	120 & 300
4.3	42 & 222
4.8	0 & 180
5.9	75 & 255
7.1	156 & 336
8.6	90 & 270

.

The analysis is performed using SCADA data with 10-min statistical values from the 48 turbines, which record electrical power, rotational speed, local wind speed, pitch angle, and nacelle position measurements. The SCADA data measurement campaign was conducted from 1 January 2008 until 31 December 2012. Furthermore, meteorological measurements are used from a nearby meteorological mast.

The SCADA data is carefully prepared to ensure the validity of the used measurements by removing unrealistic values from the data set. Examples of unrealistic values are wind speeds above 50 m/s and pitch angles larger than 90°. Additional outliers are identified and removed using the Local Outlier Factor (LOF) method [43]. The LOF algorithm computes a score that indicates the abnormality of observations based on the neighboring observations. The score indicates the local density deviation of a given data point with respect to its neighbors. Removing outliers and filtering out unrealistic values yields a data set that mainly contains valid and useful data points, which capture the general behavior including natural variability. After the LOF method is applied and the data is filtered, a total of 4,592,626 10-min average power measurements are left for 48 turbines, equaling approximately 95,680 timesteps or 1.8 years.

The standard deviation of the wind speed, needed to calculate the Turbulence Intensity (TI) of the free stream flow, is only known for measurements from the nearby meteorological masts of which most time steps do not match the ones from the SCADA data measurement campaign. To be able to estimate the TI for all time steps, the meteorological mast data is mapped on the wind direction and wind speed with a bin width of 5° and 0.5 m/s, respectively. Furthermore, a map is created for every season to capture seasonal changes. The maps are not shown for brevity, but two trends can be distinguished. The TIs are the highest in the winter and the lowest in spring. Furthermore, the TIs are generally higher for north-eastern and eastern winds, somewhat lower for southern winds, and again slightly higher for western and north-western winds. This can be explained by the geographical location of the wind farm. The land that is closest to the farm is located to the east and north-west of the farm. Wind over land generally has a higher TI than wind over sea due to the roughness length of the surface.

2.2. FLORIS Model

The FLORIS model (version 2.0.1) [33] has several engineering models available. For this study, a variant of the Gaussian wake model, called the Gauss–Curl Hybrid (GCH) model, is used. The GCH model is comprised of four main elements: the near-wake region, the 2D Gaussian shape of the velocity deficit in the far-wake region, the deflection of the wake centerline, and the secondary wake effects. The GCH model uses a single wake model for wake redirection based on the paper by Bastankhah and Porté-Agel [44]. The model further uses a turbine induced turbulence model [45] and a turbulence summation model [46]. Finally, secondary wake steering is included based on the paper of King et al. [47]. The GCH model parameters will be introduced later in this section and eventually calibrated to the measurements.

2.2.1. Near Wake

The velocity deficit does not follow a 2D Gaussian distribution in the near wake. The near-wake region is modeled as a linearly converging cone. The base of the cone is located at the rotor, and the tip is located at a distance $x_0$ downstream. This cone is visualized in Figure 2 and $x_0$ is computed as

$$\frac{x_0}{D}=\frac{cos\left(\gamma \right)\left(1+\sqrt{1-C_T}\right)}{\sqrt{2}\left[4\alpha I_{rotor}+2\beta \left(1-\sqrt{1-C_T}\right)\right]}.$$

(1)

In Equation (1), $\gamma$ is the yaw misalignment angle with the incoming flow in radians, $C_T$ is the thrust coefficient of the rotor, $I_{rotor}$ is the TI at the rotor of the turbine. $\alpha$ and $\beta$ are the first two model parameters, which governs the influence of turbulence and thrust coefficient on the location of the end of the near wake, respectively.

2.2.2. Velocity Deficit and Wake Deflection

The far-wake region starts at a distance $x_0$ downstream from the turbine. It is assumed that the wake velocity deficit has the shape of a 2D Gaussian distribution from this point onwards until the wake is fully recovered. The wake velocity deficit is expressed as

$$\frac{U(x,y,z)}{U_{\infty }}=1-\left(1-\sqrt{1-\frac{{\sigma }_{y0}{\sigma }_{z0}}{{\sigma }_y{\sigma }_z}{C}_T}\right)exp\left(\frac{{\left(y-{\delta }_f\right)}^2}{2{\sigma }_y^2}+\frac{z^2}{2{\sigma }_z^2}\right),$$

(2)

with the distance ${\delta }_f$ given by -4.6cm0cm

$${\delta }_f={\delta }_r\left(x\right)+tan\left(\theta \right)x_0+\frac{\theta }{5.2}\left(C_0^2-3e^{1/12}{C}_0+3e^{1/3}\right)\sqrt{\frac{{\sigma }_{y0}{\sigma }_{z0}}{k_y{k}_z{C}_T}}ln\left[\frac{\left(1.6+\sqrt{C_T}\right)\left(1.6S_{\sigma }-\sqrt{C_T}\right)}{\left(1.6-\sqrt{C_T}\right)\left(1.6S_{\sigma }+\sqrt{C_T}\right)}\right],$$

(3)

and the initial deflection angle defined as

$$\theta \approx \frac{0.3\gamma }{cos\left(\gamma \right)}\left(1-\sqrt{1-C_{T}cos\left(\gamma \right)}\right).$$

(4)

In Equation (2), $\sigma$ is the standard deviation in the specified direction (x, y, z) in the Euclidean space. This standard deviation defines the width of the Gaussian-shaped wake deficit.

The Gaussian-shaped wake deficit is centered around the centerline of the wake. This centerline is displaced with a distance ${\delta }_f$ in the y-direction from the x-axis as indicated in Figure 2. This displacement occurs due to yaw misalignment and wake rotation and can be calculated with Equation (3), where $\theta$ is the initial deflection angle and can be calculated using Equation (4). Note, that $\gamma$ is in radians and defined in FLORIS as positive counter-clockwise relative to the wind direction (in top view), hence Figure 2 illustrates a situation with negative yaw misalignment.

$$\begin{array}{cc}\hfill C_0& =1-\sqrt{1-C_T}\hfill \end{array}$$

(5)

$$\begin{array}{cc}\hfill S_{\sigma }& =\sqrt{\frac{{\sigma }_y{\sigma }_z}{{\sigma }_{y0}{\sigma }_{z0}}}\hfill \end{array}$$

(6)

$$\begin{array}{cc}\hfill {\sigma }_y& ={\sigma }_{y0}+(x-x_0)k_y,\mathrm{with}{\sigma }_{y0}=\frac{D}{2\sqrt{2}}cos\left(\gamma \right)\hfill \end{array}$$

(7)

$$\begin{array}{cc}\hfill {\sigma }_z& ={\sigma }_{z0}+(x-x_0)k_z,\mathrm{with}{\sigma }_{z0}=\frac{D}{2\sqrt{2}}\hfill \end{array}$$

(8)

$$\begin{array}{cc}\hfill k_y& =k_z=k_a{I}_{rotor}+k_b\hfill \end{array}$$

(9)

$$\begin{array}{cc}\hfill {\delta }_r& =a_{d}D+b_{d}x\hfill \end{array}$$

(10)

Furthermore, $C_0$ is the velocity deficit at the start of the far-wake region defined in Equation (5) and $S_{\sigma }$ is defined as Equation (6), while ${\sigma }_y$ and ${\sigma }_z$ are defined in Equations (7) and (8), respectively. The coefficients $k_y$ and $k_z$ are linear wake expansion coefficients, which governs the wake recovery, similar to the entrainment constant defined by Jensen [48]. The wake recovery is dependent on $I_{rotor}$ as shown by Pena Diaz et al. [49], and hence $k_y$ and $k_z$ are computed using two additional model parameters, namely $k_a$ and $k_b$, which gives weight to the influence of the turbulence intensity ($k_a$) and the fundamental wake recovery ($k_b$) as seen in Equation (9). Finally, ${\delta }_r$ represents an additional wake deflection induced by the rotation of the wake [37] and determined by the linear function shown in Equation (10), which depends on rotor size (D) and downstream distance (x) and the associated two model parameters $a_d$ and $b_d$, respectively.

2.2.3. Secondary Wake Effects

The Gaussian wake model is transformed into the GCH model through the secondary wake effects. The details of this sub-model can be found in the paper of [47], but the main principles and equations are described here. The model consists of two vortex systems. One is a result of the wake rotation, while the other is a result of the yaw misalignment and is called the counter-rotating vortex [44].

The circulation strength corresponding to these vortices can be used to find the spanwise V and vertical W velocity components based on the wake rotation and yaw misalignment angle anywhere in the space domain. The parameter $ϵ$ is introduced in the calculation of the velocity components, which represents the size of the vortex core. In this study, it is assumed that $ϵ=0.2D$, as it is the standard value in the FLORIS model. The additional wake recovery is defined as

$$u(x,y,z)=u_G(x,y,z)+\frac{W(x,y,z)(x-x_0)(y-y_0)}{\pi {\left({\alpha }_r(x-x_0)+\frac{D}{2}\right)}^2}.$$

(11)

Note, that $u_G$ is equal to Equation (2). Furthermore, the model parameter ${\alpha }_r$ is added, which dictates how much the vertical entrainment due to W affects the wake recovery. Here, ${\alpha }_r=0.03$ is used as in King et al. [47].

2.3. Calibration Method

The FLORIS model is an engineering model constructed based on simplified physics, and which requires calibration. However, the input power curve provided to FLORIS is crucial in getting a realistic output, and an empirical power curve is therefore derived initially before calibrating the model parameters. Once FLORIS has calculated the local wind speed for every turbine, it should also accurately calculate the power using the empirical power curve and the power-yaw loss coefficient $\alpha$ through

$$P_{\gamma }\left(U\right)=P_0\left(U\right){cos}^{\alpha }\left(\gamma \right)$$

(12)

where $\gamma$ is the yaw misalignment angle and $P_0\left(U\right)$ the power from the power curve for local wind speed U.

The empirical power curve is derived using bin averaged SCADA data from free standing turbines at the edge of the farm to exclude wake effects. It should be noted that the epistemic uncertainty in the SCADA data related to measurement uncertainty is essentially unknown and therefore not included in the subsequent calibration. An empirical power curve should reduce the effects on input uncertainty because it encompasses similar uncertainties of wind speed and wind direction measurements as will be used in the calibration. Excluding the input uncertainty from the calibration can have significant effects on the obtained model parameters, because it leads to an increased overestimation in uncertainty of the outputs, which lead to larger uncertainty of the model parameters, see Murcia Leon [50].

The power-yaw loss coefficient $\alpha$ is implemented as suggested by the paper of Liew et al. [51]. This paper describes how $\alpha$ should be adjusted based on whether a turbine is in the wake of another turbine and what the spacing is between those two turbines. Since the turbine used by Liew et al. [51] has similar characteristics as the Siemens SW-2.3-93 turbine, the values from that paper are directly used in this study. Note, that a different $\alpha$ is used per turbine spacing.

Previous model calibrations performed by for instance Fleming et al. [52], Doekemeijer et al. [14], Doekemeijer et al. [40], and Gebraad et al. [37] relied on high-fidelity flow model which was not available in the measured data. Instead, the calibration is set up to minimize the difference between the FLORIS power output and SCADA power measurements by adjusting the model parameters, similar to Göçmen and Giebel [53]. Here, the focus is only on the model parameters related to the Gaussian wake model, i.e., the following model parameters will be calibrated $\Psi =\left[\alpha \beta k_a{k}_b{a}_d{b}_d\right]$. The optimization used to find the model parameters is defined as

$${\Psi }_{opt}^i=arg\underset{{\Psi }^i}{min}\frac{1}{n}\sum \left|\frac{P_{\mathrm{SCADA}}^i-P_{FLORIS}^i(\Psi )}{P_{\mathrm{SCADA}}^i}\right|,$$

(13)

where $P^i=\left[P_1^i{P}_2^i\dots P_n^i\right]$, n indicating the number of turbines considered, and i corresponds to a combination of wind speed $U_{\infty }$ and wind direction ${\theta }_{\infty }$, and hence, also turbine spacing.

The calibration is performed using Differential Evolution [54], which yields the global minimum of a multivariate function within defined boundaries. These boundaries are obtained from Doekemeijer et al. [40]. The boundaries for $\alpha$ and $\beta$ from Doekemeijer et al. [40] are adjusted to match the Python version of FLORIS. Furthermore, the lower boundary of $k_b$ is changed from −0.01 to 0.0. The boundaries used here are listed in

Table 2. The boundaries per model parameter defined for the optimizer.

Parameter	Lower Bound	Upper Bound
$\alpha$	$0.125$	$2.5$
$\beta$	$0.015$	$0.3$
$k_a$	$0.05$	$1.5$
$k_b$	$0.0$	$0.02$
$a_d$	$-1.0$	$1.0$
$b_d$	$-0.1$	$0.1$

.

2.4. Yaw Optimization

FLORIS is calibrated for normal operating conditions and subsequently used to determine optimal yaw angles aimed at optimizing the power output and hence AEP. The potential improvements will be compared to a baseline scenario, where all turbines have a zero yaw angle and are thus aligned with the main wind direction. The yaw angles of the turbines will be optimized in two ways, namely deterministic, where input uncertainty is neglected, and robust, where input uncertainties are considered.

2.4.1. Deterministic Optimization

There are a total of 48 variables to optimize, as the yaw angle of each of the turbines is the only control parameter. The cost function used in the minimization is the negative wind farm power and the optimization tool uses a sequential least squares programming algorithm. The aim is to increase the AEP, and hence the optimization is performed for a full wind rose. The wind rose is obtained from data of meteorological mast 1, which conducted a measurement campaign of 2.5 years before the construction of the wind farm. The wind rose is discretized for $\Delta U_{\infty }=0.5$ m/s and $\Delta {\theta }_{\infty }=5°$ and is visualized in Figure 3.

The measurements from the meteorological mast include measurement uncertainty, but as previously described, the influence is reduced by using an empirical power curve, which includes the same uncertainty in wind direction and wind speed.

The deterministic yaw angle sets are indicated as ${\overrightarrow{\gamma }}_{w,0}$. Once the optimal yaw angles are determined for the full wind rose, the AEP can be calculated for the baseline and optimized case as

$$AE{P}_k=T\left[\sum _i\sum _j\left(p(U_{\infty ,i},{\theta }_{\infty ,j})P_k(U_{\infty ,i},{\theta }_{\infty ,j})\right)\right],$$

(14)

where T is the number of hours in a year, $p(U_{\infty ,i},{\theta }_{\infty ,j})$ is the probability that a specific combination of wind speed and direction occurs, and $P_k(U_{\infty ,i},{\theta }_{\infty ,j})$ is the power corresponding to that combination of wind speed and direction. Note, that the subscript k indicates whether it is the baseline or the optimized case. In order to better quantify the difference between baseline and optimized case, the AEP gain $G_f$ is used instead of the AEP alone. The AEP gain is defined as

$$G_f=100\frac{\left(AE{P}_{opt}-AE{P}_{base}\right)}{AE{P}_{base}}.$$

(15)

2.4.2. Optimization under Uncertainty

The optimization can be made more robust by taking the associated uncertainty in wind direction ${\sigma }_{\theta }$ and yaw position ${\sigma }_{\gamma }$ into account, i.e., performing Optimization Under Uncertainty (OUU) [16,17,18,55]. The present study is limited to these uncertainties, but many additional input uncertainties can be taken into account, such as uncertainties in wind speed, TI, and wind shear. Following the study of Simley et al. [18], the input uncertainties are related and derived as

$${\sigma }_w^2={\sigma }_{\theta }^2+{\sigma }_{\gamma }^2=5.25°.$$

(16)

In [18], it is also shown that ${\sigma }_{\theta }=4.95°$ and ${\sigma }_{\gamma }=1.75°$ yields the best fit for ${\sigma }_w=5.25°$. It should be mentioned that the use of these standard deviations for all inflow conditions and turbines is a significant simplification. The robust optimization, which takes the input uncertainties into account, will be denoted ${\sigma }_{w,1}$, and the deterministic input case will be denoted ${\sigma }_{w,0}$. It is assumed that two representative normal distributed Probability Density Functions (PDFs) can be obtained using those standard deviations, i.e.,

$$p\left({\Delta }_i\right)\sim \mathcal{N}(0,{\sigma }_i^2).$$

(17)

The PDFs indicate the probability of a deviation from the mean ${\Delta }_i$. The subscript i indicates either the wind direction or the yaw position in both ${\sigma }_i$ and ${\Delta }_i$. The PDFs are used to obtain a weighted average of the wind farm power for combinations of wind direction and yaw position errors. The weighted average of power as defined by Simley et al. [18] is

$$P_w({\theta }_{\infty },\gamma )={\int }_{-\pi }^{\pi }{\int }_{-\pi }^{\pi }\rho \left({\Delta }_{\theta },{\Delta }_{\gamma }\right)\sum _j^{N_t}{P}_{FLORIS}^j(\theta +{\Delta }_{\theta },\gamma +{\Delta }_{\gamma })d{\Delta }_{\theta }d{\Delta }_{\gamma },$$

(18)

where the superscript j denotes the turbine number. Note, that $\rho \left({\Delta }_{\theta },{\Delta }_{\gamma }\right)$ represents the joint probability distribution of wind direction and yaw position. The joint probability distribution derivation can be found in the paper of Simley et al. [18]. Equation (18) is implemented by approximating the integration as a double summation and by discretizing the wind direction and yaw position. Once calculated, $P_w$ can be used to obtain the robust AEP through Equation (14) by replacing P.

Determining the weighted average of the farm power takes longer than when the farm power is calculated using the deterministic approach because the model needs to be evaluated multiple times with different combinations of wind directions and yaw positions. The number of evaluations N needed for each $P_w$ depends on the number of discretization points of the wind direction PDF $N_{\theta }$ and yaw position PDF $N_{\gamma }$, through

$$N=N_{\theta }·N_{\gamma }.$$

(19)

A convergence study is done to determine the number of evenly spaced steps needed per PDF, to ensure that the weighted average calculation is as fast as possible, without losing too much accuracy. It is found that generally convergence is reached for $N_{\theta }=7$ and $N_{\gamma }=5$, which is similar to what is used by Doekemeijer et al. [14]. The robust yaw angle sets are indicated as ${\overrightarrow{\gamma }}_{w,1}$.

3. Results

3.1. Sobol Method

The current model setup of FLORIS contains a total of six model parameters, which all affect the accuracy of the model results, but the influence differ for each parameter. Therefore, it is important to examine and quantify the influence of each parameter through a sensitivity analysis. This is done using the variance-based Sobol method [56,57] which is a global and model-independent sensitivity analysis method based on variance decomposition. It can handle non-linear and non-monotonic functions and models. The method decomposes the variance of the model output into fractions that can be connected to parameters or sets of parameters in the model. Hence, the method makes it possible to distinguish highly influential parameters from non-influential parameters.

For full details on the variance-based Sobol method, the reader is encouraged to read the paper of Sobol [56,57] with examples on how to apply the Sobol method by Nossent et al. [58] and Zhang et al. [59]. The Sobol sensitivity indices are characterized by the ratio of the partial variance of the output to the total variance. The model output used for this method is the AEP, meaning that all wind speeds and wind directions are taken into account. The method depends on the choice of distribution and boundaries for each of the evaluated parameters. Uniform distributions have been used, and the boundaries for the model parameters are the same as used in Section 2.3 and shown in Table 2. Note, that the boundary values influence the Sobol index of the parameter, as a larger range for a given parameter will likely yield a higher Sobol index. In this paper, only the first-order sensitivity $S_i$ and the total-order sensitivity $S_{Ti}$ are considered, where $S_i$ is a measure of the contribution of the individual parameter ${\Psi }_i$ to the total model variance, and $S_{Ti}$ describes the main effect of ${\Psi }_i$ and all its interactions with other parameters. In other words, one can determine the impact of the interaction between ${\Psi }_i$ and other parameters by analyzing the difference between $S_{Ti}$ and $S_i$.

The reliability of the Sobol indices can be tested by using confidence intervals. Bootstrapping with re-sampling is used to obtain 100 samples from the model output distribution [60] and hence distributions of $S_i$ and $S_{Ti}$. The 95% confidence intervals can be determined from these distributions. In general, the confidence intervals of the most sensitive parameters should be as small as possible, and ideally smaller than 10% of the sensitivity indices for the most sensitive parameters [59]. If the confidence intervals are larger than 10%, it means that the Sobol indices are not fully converged yet and that the model should be evaluated more times. For the indices shown, the model has been evaluated 12,000 times.

The first and total-order sensitivity indices for a turbine spacing of 4.8D are shown in Figure 4. The confidence intervals are shown in black and the 10% intervals in red. First of all, the figure shows that the confidence intervals of the most sensitive parameters are within, or very close to, the 10% deviation. Furthermore, it can be seen that only three model parameters have a significant effect on the model output variance, namely $k_a$, $b_d$, and to some degree $\alpha$. In total these three parameters account for approximately 79% of the model output variance. This means that the three sub-models, describing the near-wake region, wake deflection, and velocity deficit, all influence the output significantly. However, the insignificant effects of some of the model parameters indicate that certain sub-models or the combined model complex might be overparameterized. The effect of $a_d$ is negligible, hence this parameter is set to zero, every time the model is used. Finally, for all parameters, the total-order sensitivity is significantly higher than the first-order sensitivity. This indicates that there are some significant second-order (and probably higher-order) sensitivities, meaning that there is significant parameter interaction.

The Sobol indices for all main directions and the associated spacings are presented in Figure 5. Figure 5 shows that the parameters are generally converged with the 95% confidence intervals being within the 10% deviation, except for $\alpha$ and $b_d$ for the largest two spacings. This again implies that the estimations of the Sobol indices are not fully converged for these two parameters. However, the results still indicate the parameters’ dependence in relation to the turbine spacing. The changing sensitivity of the model parameters for each turbine spacing highlights the need to map the model parameters for each turbine spacing during the calibration, which is done in Section 3.2. Additionally, the negative values for $a_d$ are non-physical but occur because the Sobol method is not fully converged, and since the Sobol index of $a_d$ is close to zero.

The three most sensitive parameters have different trends for varying turbine spacing. It can be seen that the sensitivity of $k_a$ increases as the turbine spacing increase. The sensitivity of $\alpha$ and $b_d$ follow the opposite trend; they are higher for relative short turbine spacings and lower for larger turbine spacings. These changes in sensitivity are expected because $\alpha$ influences the near-wake, which is less influential for larger spacings. The lateral deflection influenced by $b_d$ Equation (10) is relatively small compared to the lateral deflection described by the other parameters Equation (3) for larger spacings. The observed trends in sensitivities are similar for both the baseline and the optimized scenario.

3.2. Calibration Method

Figure 6 shows the power predictions of FLORIS before and after the calibration compared to the power measurements from the SCADA data for $U_{\infty }$ = 8 m/s and ${\theta }_{\infty }=180°$, which is a frequently occurring wind combination. The plot shows the output of FLORIS for the standard model parameters ${\Psi }_{org}$ and for the optimal model parameters ${\Psi }_{opt}$, including the Mean Absolute Percentage Error (MAPE) per turbine. It can be seen that ${\Psi }_{opt}$ performs much better than ${\Psi }_{org}$ and that all predictions of FLORIS using ${\Psi }_{opt}$ fall within the 68% confidence interval of the SCADA data. The average MAPE of all turbines shown in Figure 6 is minimized to 5.7%.

However, the average MAPE increases if ${\Psi }_{opt}^A$ is applied for other wind speeds. A new set of model parameters ${\Psi }_{opt}^B$ can be found for each wind speed, where A represents the combination $U_{\infty }$ = 8 m/s and ${\theta }_{\infty }=180°$ and B represents for example $U_{\infty }$ = 10 m/s and ${\theta }_{\infty }=180°$. Hence, the model parameters are wind speed dependent. Therefore, ${\Psi }_{opt}^i$ is found for every wind speed independently. Figure 7 shows how the model parameters change for varying wind speed. Generally, the parameter values are spread over the entire domain, except for $k_a$ and $k_b$. It can also be seen that $k_a$ is often close to its lower boundary and that $k_b$ is closer to the upper boundary. However, it should be noted that the sensitivity of $k_b$ is relatively low, as discussed in Section 3.1. Hence, $k_a$ is still the dominating parameter for the wake expansion. There is also a clear trend that the major change in parameter value occurs as the wind speed increases above rated, where the wake effects are significantly reduced. The model parameters are only obtained for wind speeds from 5 m/s up to 17.5 m/s, since all turbines operate at rated power for higher wind speeds, even at the smallest turbine spacings, and downstream waked turbines do not operate when the upstream turbine experiences wind speeds below 5 m/s.

The variation of the six model parameters shown in Figure 7 can also be determined for every turbine spacing. The corresponding average errors of the predicted power production for each wind speed and turbine spacing is shown in Figure 8. The highest MAPEs occur for small turbine spacings and the lowest MAPEs for large turbine spacings. This is expected as the FLORIS model is intended for far-wake modeling. Furthermore, the behavior of the MAPE is the same for every turbine spacing. It is generally the highest for low wind speeds and it decreases slightly with increasing wind speed before it converges to zero as all turbines in the row start to operate at rated power. The standard deviation of the MAPE follows the same trend but is not shown for brevity.

The shown MAPEs are based on the calibrated model parameters for each wind speed and each turbine spacing. The MAPE will increase significantly if the model parameters of the case with a spacing of 3.3D are applied on another spacing, e.g., 5.9D. Hence, all model parameters need to be calibrated as a function of wind speed and turbine spacing instead of a single constant value in order to achieve accurate model predictions. These dependencies might also depend on additional conditions, such as atmospheric stability.

3.3. Model Uncertainty

The model bias and uncertainty can be quantified using the SCADA data as a benchmark. The method used is proposed and fully described by Nygaard [61]. The main principle is to estimate the discrepancy distribution between measured data and the model output. The mean of this distribution gives the bias of the model and the width of the distribution gives an estimation of the uncertainty. The method relies on the power output as a function of inflow wind speed and wind direction only. The benchmark of the model is the wake loss defined as

$$L_i=1-\frac{{\langle P_{net}\rangle }_i}{{\langle P_{gross}\rangle }_i}.$$

(20)

The angle brackets 〈..〉 denote the mean value of the validation sample considered. This validation sample is derived from the available SCADA data. The subscript i indicates the case, i.e., either the observed case using measured values $L_{obs}$ or the model case $L_{model}$. The gross power $P_{gross}$ is the power estimation produced by the farm with no wake effects, and the net power $P_{net}$ is the power estimation of the wind farm including wake effects. The wake losses can be used to find the relative wake model error defined as

$$ϵ=\frac{\left(L_{obs}-L_{model}\right)}{L_{obs}}.$$

(21)

The model underestimates the wake losses for $ϵ>0$, i.e., the estimated AEP is optimistic. By contrast, the estimated AEP is conservative when $ϵ<0$.

A statistical re-sampling approach, called the circular block bootstrap method [62], is used to get a distribution of the relative wake model error. This method divides the validation sample, which essentially is a time series of SCADA data, into overlapping blocks of length k and wraps the time series to connect the beginning and the end. The created blocks are then randomly sampled with replacement until a new unique time series is created with the same size as the original validation sample. This unique time series is then one bootstrap sample of which $L_{obs}$, $L_{model}$ and $ϵ$ can be derived. This process is repeated N times to get a distribution of $ϵ$.

The method is applied with $k=798$ time steps and $N=5000$ samples. The high value of N is chosen to obtain a smooth distribution, but the sample size does not affect the resulting bias and uncertainty from N = 500 onwards. The corresponding relative wake model error distribution is shown in Figure 9. The mean value, hence the bias, is −5.8%, meaning that the model estimate tends to be conservative. The uncertainty, as mentioned before, is defined as the half-width of the 68% confidence interval. In this case, the model uncertainty is ±1.1%. Although there is a significant bias, the FLORIS model is still able to estimate the wake loss with low uncertainty. This model uncertainty means a wake loss of 30% essentially has a 68% confidence interval of (29.7%, 30.3%). However, this does not include the determined bias, which corrects the predicted wake loss with a 68% confidence interval to (27.9%, 28.6%).

The wind farm efficiency can be calculated through

$${\eta }_{f,d}=\frac{\langle P_{net,d}\rangle }{\langle P_{gross,d}\rangle }=1-\langle L_d\rangle ,$$

(22)

where the subscript d denotes the wind direction bin and the angle brackets $\langle ...\rangle$ indicate the mean value of the power measurements of all timesteps for a given wind direction measurement in d. Note, that the wind farm efficiency is the same as one minus the wake loss. The wind farm efficiency is plotted for different wind directions in Figure 10. It is seen that the FLORIS model output follows the general trends and matches most peaks and troughs of the observed data, except for two turbine spacings, namely 5.9D (75° & 255°) and 8.6D (90° & 270°). The wake loss is especially overestimated by the FLORIS model in the 10–110° sector. The error bars show the 25th and 75th quartiles of the observed and modeled values, and hence show the width of the distributions. Generally, the two distributions overlap, particularly where the model performs well, which indicates that results are statistically comparable and the variation in the model input does to some degree capture part of the variability seen in the observations. However, the modeled distribution frequently appears more skewed than the observed values, e.g., 50–70° and 130–220°, while the observed distributions appear more symmetric.

The discrepancy between the observed and modeled efficiencies can potentially be attributed to the influence of the atmospheric stability. The atmospheric stability was previously estimated for the Lillgrund wind farm from the nearby Drogden light tower for the period of 2008–2013 by Larsen et al. [63] and has been categorized as either stable, neutral, or unstable based on the Monin–Obukhov length [64,65]. However, it has not been possible to directly link all the stability measurements with the SCADA data, but the stability classification is shown as a function of wind direction in Figure 10 to visualize the correlation between the atmospheric stability and the ability of the FLORIS model to predict the wake loss accurately. The stability is predominantly unstable for eastern winds and north-western winds. This is presumably because the winds during autumn and winter are often cold winds coming from land, which can result in an unstable atmosphere when it is traveling over a relatively warmer surface as the sea, which is not fully cooled down yet.

The FLORIS model mainly overestimates the wake loss in the 10–110° wind direction sector, where the atmosphere is predominantly unstable. The analysis suggests that the FLORIS model with the current wind shear setting is, in general, better able to model stable and neutral atmospheric conditions compared to unstable atmospheric conditions. This highlights the importance of taking into account the atmospheric stability. Atmospheric stability effects can potentially be included by adjusting the wind shear accordingly and to include the separate effect on changing TI [66]. However, such adjustments are not applied in this study, because adjusting the TI and wind shear based on the stability would require additional assumptions and hence additional uncertainties.

3.4. Yaw Optimization

Section 2.4 described how the optimized yaw angle sets ${\overrightarrow{\gamma }}_{{\sigma }_{w,0}}$ and ${\overrightarrow{\gamma }}_{{\sigma }_{w,1}}$ can be derived. The performance of each yaw angle set is tested for three input cases, which are summarized in

Table 3. Gains per input condition case for the deterministic yaw angle set (${\overrightarrow{\gamma }}_{{\sigma }_{w,0}}$) and the robust yaw angle set (${\overrightarrow{\gamma }}_{{\sigma }_{w,1}}$).

Input Case	${\mathit{\sigma }}_{\mathit{i}}$ (deg)	Gain ${\overrightarrow{\mathit{\gamma }}}_{{\mathit{\sigma }}_{\mathit{w},0}}$ (%)	Gain ${\overrightarrow{\mathit{\gamma }}}_{{\mathit{\sigma }}_{\mathit{w},1}}$ (%)
${\sigma }_{w,0}$	0.0	6.1	4.3
${\sigma }_{w,1}$	5.25	2.5	3.4
${\sigma }_{w,2}$	10.5	1.1	2.2

. The first column of Table 3 shows the three input cases. Note, that an additional input case ${\sigma }_{w,2}$ is introduced which is used to test the robustness of the yaw angle sets since both yaw sets are not optimized for this input case. The second column shows the input uncertainty for each case given as a standard deviation of a normal distribution. The third column shows the potential AEP gain for each input case for the optimized yaw set determined using input case ${\sigma }_{w,0}$. Finally, the fourth column shows the potential AEP gain per input case for the optimized yaw set determined using input case ${\sigma }_{w,1}$.

The highest AEP gain of 6.1% can be observed for the deterministic yaw set ${\overrightarrow{\gamma }}_{{\sigma }_{w,0}}$ and deterministic input case ${\sigma }_{w,0}$. Assuming the same model uncertainty for the optimized case as for the baseline case yields a 68% confidence interval of the energy gain as (5.5%, 6.6%). Furthermore, it is also assumed that the model bias is the same for the baseline and optimized case, i.e., bias influences both $AE{P}_{base}$ and $AE{P}_{opt}$ in Equation (15) in the same way, and therefore, the model bias does not influence the estimated relative energy gain.

The energy gain for the more robust yaw set ${\overrightarrow{\gamma }}_{{\sigma }_{w,1}}$ and deterministic input case ${\sigma }_{w,0}$ is 2.5% with a 68% confidence interval of (1.8%, 3.1%), which is significantly lower than the 6.1% observed for ${\overrightarrow{\gamma }}_{{\sigma }_{w,0}}$. The deterministic input case assumes perfect knowledge of the wind direction and yaw position, which is usually not the case. Higher gains are observed for input cases ${\sigma }_{w,1}$ and ${\sigma }_{w,2}$ when input uncertainties are included in the yaw optimization (${\overrightarrow{\gamma }}_{{\sigma }_{w,1}}$). For example, the energy gain for input case ${\sigma }_{w,1}$ increases from 2.5% to 3.4%, which is an increase of 36%. On the other hand, a lower gain can be observed for input case ${\sigma }_{w,0}$.

The yaw set ${\overrightarrow{\gamma }}_{{\sigma }_{w,1}}$ also performs better for the larger input uncertainty case ${\sigma }_{w,2}$. Thus, ${\overrightarrow{\gamma }}_{{\sigma }_{w,1}}$ can be considered to be more robust than ${\overrightarrow{\gamma }}_{{\sigma }_{w,0}}$. The most realistic input uncertainty of the cases considered is presumably ${\sigma }_{w,1}$. The final estimation of the potential AEP gain for the Lillgrund wind farm is, therefore, 3.4% with a 68% confidence interval of (2.8%, 4.0%) using ${\overrightarrow{\gamma }}_{{\sigma }_{w,1}}$, which is in the same order as the 3.7% found by Gebraad et al. [38] for the Princess Amalia Wind Park in the Netherlands.

In order to identify for which wind directions the energy gain is achieved, and what the differences are between the gains of yaw angle sets ${\overrightarrow{\gamma }}_{{\sigma }_{w,0}}$ and ${\overrightarrow{\gamma }}_{{\sigma }_{w,1}}$, the gains are plotted as a function of wind direction in Figure 11 for input case ${\sigma }_{w,1}$. It is seen how ${\overrightarrow{\gamma }}_{{\sigma }_{w,0}}$ occasionally yields locally negative gains, while ${\overrightarrow{\gamma }}_{{\sigma }_{w,1}}$ constantly shows positive gains. The largest differences between the two curves are observed for alignment directions (90°, 155°, 180°, etc.).

Figure 12 shows the yaw angle sets for the three different scenarios for a wind direction of 80° and a wind speed of 8 m/s. The more robust yaw angle set ${\overrightarrow{\gamma }}_{{\sigma }_{w,1}}$ (lower number in red) is compared with the deterministic yaw angle set ${\overrightarrow{\gamma }}_{{\sigma }_{w,0}}$ (upper number in black). In general, it can be seen that the robust yaw angles are smaller than the deterministic ones. This is an expected consequence of including uncertainty of yaw and wind direction because the sensitivity of power is higher for larger yaw angles and the optimization does therefore result in smaller yaw angles, which are more robust against changing wind directions.

Similar behavior is observed by Rott et al. [17], who also shows that more robust yaw angles result in a lower yaw activity, which is beneficial especially for the turbine control. Furthermore, it is evident how the most downstream turbines generally have a slightly larger yaw angle for ${\overrightarrow{\gamma }}_{{\sigma }_{w,1}}$. This is again because the turbines prefer to align themselves for the combination of yaw and wind direction deviations with fewer wake effects, in order to increase the power produced in these situations [17].

4. Conclusions

Wake effects impact the aerodynamic performance of large wind farms, which results in decreased power production of wind turbines operating in clusters as opposed to free-standing turbines. This work uses the wind farm model FLORIS to estimate the potential of reducing wake losses by applying yaw-based wake steering. The Lillgrund wind farm is used as a test case. The methodology is divided into three steps: model calibration, sensitivity study, and optimization including uncertainty.

Firstly, a global sensitivity analysis was performed to assess the relative effect of the various model parameters. It was found that three model parameters account for approximately 79% of the total model output variance, while the others have a smaller influence, indicating that the model could be over-parameterized.

Secondly, the model parameters were calibrated to match the power production of the Lillgrund wind farm. The calibration clearly showed that the amount of available data and the data quality are important. In particular, the model parameters were dependent on both wind speed and turbine spacing, so their values should be considered carefully. The performance of the calibrated model was tested by comparing wake losses with measured data. Overall, the wake losses were generally overestimated by 5.8% but the model estimates were considered to be reliable with only a 1.1% model uncertainty. It was also found that the model bias correlated with the changes in predominant atmospheric conditions from different wind directions, suggesting that more data on the atmospheric stability would improve the model predictions. Thirdly, FLORIS was used to optimize the individual turbine yaw angles, in order to increase the potential energy gain of the wind farm. Three cases were considered: a deterministic case without uncertainty in the wind direction and yaw settings, and two cases with uncertainty in the wind direction and yaw settings. The potential AEP gain was larger for the deterministic input case (6.1%) compared with the deterministic yaw angle case with input uncertainty (2.5%). The potential energy gain was shown to increase further by including uncertainty in the yaw angles during the optimization procedure. In particular, the AEP gain was increased by 36% when OUU was used for the yaw angles, leading to a more robust yaw angles set. Overall, the observed AEP gain of 3.4% with OUU shows the potential of wake steering as a control strategy to increase the power production of wind farms.

This work was performed in the context of the Lillgrund wind farm, which is very dense compared to modern wind farms. Future work should assess the potential for wake steering when increasing turbine spacing. Furthermore, this work relies on two fundamental assumptions: (i) neglecting the input uncertainty during calibration, and (ii) constant bias and uncertainty for all model applications (i.e., for both yawed and non-yawed cases). These strong assumptions would need to be validated based on long term measurement campaigns. Finally, this paper showed the importance of the amount of available data and the data quality. Additional data on atmospheric stability could refine the present findings.