Polynomial chaos to efficiently compute the annual energy production in wind farm layout optimization

In this paper, we develop computationally-efficient techniques to calculate statistics used in wind farm optimization with the goal of enabling the use of higher-fidelity models and larger wind farm optimization problems. We apply these techniques to maximize the Annual Energy Production (AEP) of a wind farm by optimizing the position of the individual wind turbines. The AEP (a statistic) is the expected power produced by the wind farm over a period of one year subject to uncer5 tainties in the wind conditions (wind direction and wind speed) that are described with empirically-determined probability distributions. To compute the AEP of the wind farm, we use a wake model to simulate the power at different input conditions composed of wind direction and wind speed pairs. We use polynomial chaos (PC), an uncertainty quantification method, to construct a polynomial approximation of the power over the entire stochastic space and to efficiently (using as few simulations as possible) compute the expected power (AEP). We explore both regression and quadrature approaches to compute the 10 PC coefficients. PC based on regression is significantly more efficient than the rectangle rule (the method most commonly used to compute the expected power). With PC based on regression, we have reduced on average by a factor of five the number of simulations required to accurately compute the AEP , thus enabling the :::: when ::::::::: compared :: to ::: the :::::::: rectangle :::: rule ::: for ::: the ::::::: different :::: wind ::::: farm :::::: layouts ::::::::: considered. ::: In :: the ::::: wind :::: farm :::::: layout :::::::::: optimization :::::::: problem, :::: each ::::::::::: optimization ::: step ::::::: requires :: an ::::: AEP :::::::::: computation. :::::: Thus, ::: the :::::: ability :: to ::::::: compute ::: the ::::: AEP ::::::::: accurately :::: with ::::: fewer :::::::::: simulations :: is ::::::::: beneficial :: as :: it :::::: reduces :::: the ::: cost ::: to 15 :::::: perform ::: an ::::::::::: optimization, ::::: which ::::::: enables ::: the use of more expensive, ::::::::::::: computationally :::::::: expensive : higher-fidelity models or larger wind farm optimizations :: the :::::::::::: consideration :: of :::::: larger :: or ::::::: multiple ::::: wind :::: farm ::::::::::: optimization ::::::::: problems. ::: We ::::::: perform :: a :::: large ::::: suite :: of :::::::::::: gradient-based :::::::::::: optimizations :: to :::::::: compare ::: the :::::: optimal ::::::: layouts :::::::: obtained ::::: when ::::::::: computing ::: the :::: AEP ::::: with ::::::::: polynomial :::::: chaos ::::: based :: on ::::::::: regression ::: and ::: the :::::::: rectangle :::: rule. :::: We ::::::: consider :::: three :::::::: different :::::: starting ::::::: layouts ::::: (Grid, ::::::: Amalia, :::::::: Random) ::: and :::: find :::: that :: the ::::::::::: optimization ::: has ::::: many ::::: local :::::: optima ::: and :: is :::::::: sensitive :: to ::: the :::::: starting :::::: layout :: of ::: the :::::::: turbines. ::: We ::::::: observe ::: that ::::::: starting :::: from :: a 20 :::: good :::::: layout ::::: (Grid, ::::::: Amalia) :::: will, ::: in ::::::: general, ::: find ::::: better ::::::: optima :::: than :::::: starting ::::: from : a :::: bad ::::: layout ::::::::: (Random) :::::::::: independent ::: of ::: the :::::: method :::: used :: to ::::::: compute ::: the ::::: AEP.We show how the PC method can be used to efficiently compute the gradients of the AEPand we perform a large suite of gradient-based optimizations with different initial turbine locations and with different numbers of samples ::: For :::: both ::: PC ::::: based ::: on :::::::: regression :::: and :: the :::::::: rectangle :::: rule, ::: we ::::::: consider :::: both :: a :::::: coarse :::::: (⇠ 225) :::: and : a ::: fine :::::: (⇠ 625) ::::::: number :: of :::::::::: simulations to compute the AEP. We find that :: for ::::::: roughly :::::::: one-third :: of :::: the :::::::::::: computational :::: cost, ::: the ::::::::::: optimizations ::::: with the 25 optimizations with ::::: coarse PC based on regression result in optimized layouts that produce comparable AEP as the optimized

1 Introduction 20 In 2015, wind energy growth accounted for almost half of the global electricity supply growth. In the United States, it accounted for 41 % of new power capacity, raising the wind energy supply to 4.7 % of the total electricity generated in 2015 and on target to reach 10 % by 2020 (U.S. Department of Energy, 2015;AWEA, 2016;GWEC, 2016). Most of the current and upcoming wind energy comes from large turbines (greater than 1 MW) situated in clusters-wind farms. A problem with putting turbines together in confined spaces is that they operate in the wakes of other turbines, i.e., in regions of reduced speed and increased turbulence. This leads to an underproduction of power and decreased (10-20 %) energy output for the farm (Barthelmie et al., 2007(Barthelmie et al., , 2009Briggs, 2013) when compared to ideal conditions. This loss in energy capture results in millions of dollars of loss for operators and investors and increased economic uncertainty for new installations. Many current wind farms have grid-like layouts, where wind turbines are aligned in rows, which further exacerbate the wake losses. By optimizing the layout of the 5 wind farm, the wake losses can be minimized, with a corresponding increase in energy production and revenue.
Wind farm optimization is a complex, multi-disciplinary and high-dimensional problem. The wind farm may contain dozens or even hundreds of wind turbines, where each turbine may be parameterically described using several design variables. Furthermore, the wind conditions (wind direction, wind speed, wind turbulence, etc.) are stochastic (uncertain), and thus we need a statistic to evaluate the performance of the wind farm. A common statistic is the expected power or the Annual Energy 10 We see three approaches to improving wind farm optimization capabilities. Each approach focuses on the different blocks of the optimization under uncertainty problem (Fig. 1c). The first approach is to improve the modeling quality of entire wind farms, i.e., improve the models in all the disciplines (aerodynamics, structures, controls, electrical, acoustics, atmospheric physics, policy, economics, etc.) that are relevant to building and operating a wind farm, as well as the interaction between the different turbines. The second approach is to improve the optimization problem formulation and the algorithms to solve the 5 optimization. And the third approach is to improve the treatment of the stochastic nature of the problem, i.e., develop better uncertainty quantification methods to efficiently compute the relevant wind farm statistics (and their gradients with respect to the design variables of the problem) used in the optimization under uncertainty problem.
The first approach increases the fidelity of the model whereas the second (optimization) and third (uncertainty quantification) approaches seek to reduce the number of model evaluations, as this enables the study of larger and more realistic problems. 10 Here, we focus on the uncertainty quantification approach (the third approach), as it has not been considered in detail before.
The most recent and thorough review of the wind farm optimization literature (Herbert-Acero et al., 2014) does not mention it. It only mentions the first two approaches. In the existing work in the literature, the third approach typically focuses on simple integration methods to compute the statistics, which quantify the effect of the stochastic inputs (Kusiak and Song, 2010;Kwong et al., 2012;Chowdhury et al., 2013;Gebraad et al., 2017). Simple integration methods, 15 such as the rectangle rule, are inefficient in the number of samples (simulations of the model) needed to accurately estimate a statistic, such as the AEP. They are especially inefficient if multiple stochastic inputs are considered simultaneously. Normally only the wind direction and/or the wind speed are considered as stochastic input variables. Recent work (Padrón et al., 2016;Murcia et al., 2015) is starting to move beyond these simple integration techniques to compute the AEP and instead using the uncertainty quantification method of polynomial chaos to compute the AEP. 20 In this paper, which is meant as a comprehensive introduction to uncertainty quantification methods applied to wind-farm simulations, we describe in detail the polynomial chaos (PC) method and show that, for the efficient (small number of model simulations) computation of the AEP, the PC method based on regression should be used. An additional benefit of the PC method is that it makes it feasible to consider multiple uncertain variables (e.g., wind direction, wind speed, wind turbulence, wake model parameter) that impact the computation of the AEP (Padrón, 2017). 25 In addition, we show how to compute gradients of the ::: the ::: PC :::::: method :::: can :: be ::::: used :: to ::::::::: efficiently ::::::: compute ::: the ::::::: gradient ::: of statistics, such as the AEP, with a modified version of the PC method. : . :::::::::::::::::::: Eldred (2011) describes ::: how ::: to ::::::: compute ::: the :::::::: gradients D w,1 D w,2 D w,3 Note: D w,q := f (k e , m e,q ), q 2 {1, 2, 3} Figure 2. Schematic of the Floris wake model. The model has three zones with varying diameters, Dw,q, that depend on tuning parameters ke and me,q. The effective hub velocity is computed using the overlap ratio, A OL q , of the part of the rotor-swept area overlapping each wake zone respectively to the total rotor-swept area.
2 Computing the power and the annual energy production of a wind farm We first describe the aerodynamic wake model we use (Sect. 2.1). The wake model gives an estimate of the hub-height velocity at each wind turbine, from which we can compute the power produced by the wind farm (Sect. 2.2). Then, to obtain the Annual Energy Production (AEP) we need to integrate the power over all wind conditions that occur in a year (Sect. 2.3) and weigh the results proportionally to the frequency with which such wind conditions manifest themselves.

Floris
The Floris (FLow Redirection and Induction in Steady-state) 1  wake model is an enhancement of the Jensen wake model (Jensen, 1983) and the wake deflection model presented in (Jiménez et al., 2009). The Floris model builds on the Jensen model by defining three separate wake zones with differing expansion and decay rates (controlled by tunable coefficients) to more accurately describe the velocity deficit across the wake region. A simple overlap ratio of the area of 10 the rotor in each zone of each shadowing wake to the full rotor area is used to determine the effective hub velocity of a given turbine. A simple overview of the Floris model, showing the zones and overlap areas, is shown in Fig. 2. We use the Floris wake model with changes to remove discontinuities and add curvature to regions of non-physical zero gradient to make the model more suitable for gradient-based optimization . In this work, we use the parameter values recommended in  and  and set the yaw-offset angle of each turbine to zero. Design x x -The x location of each turbine.
y -The y location of each turbine.
Parameters θ yaw angles, turbine characteristics, and wake model parameters.

Computing the power of a wind farm
We will consider the power of the wind farm to be a function of three classes of variables: uncertain variables ξ, design variables x, and parameters θ, Uncertain variables are variables that follow a probability distribution, design variables are variables that an optimizer can vary, 5 and parameters are important constants that govern the behavior of the system. The classification of the variables is problem dependent. For instance, the rotor yaw could be considered a design variable, a parameter, or an uncertain variable to account for yaw-offset measurement error. A tunable parameter of a wake model, such as the wake expansion coefficient, could be considered as a parameter or as an uncertain variable given by a particular distribution.
For the problems considered in this work, Table 1 lists in which category we place each variable that influences the power 10 computation. The uncertain variables are the wind direction and wind speed with probability distributions described in Sect. 5.1.
The power of the wind farm for a given wind direction and wind speed is equal to the sum of the power produced by each turbine P = n turb i=1 P i . The power of each turbine is calculated from where ρ is the air density, A is the rotor swept area, C P is the power coefficient, and U i is the effective hub velocity for each 15 turbine, which is calculated by the wake model and is a function of the three types of variables described above The power coefficient captures both the aerodynamic and electromechanical properties of the wind turbine. It is a complex function of many variables (Herbert-Acero et al., 2014) and is usually reported by wind turbine manufacturers as a function of the tip-speed ratio, which depends on wind speed at hub height U i . A simple expression for C p can also be computed using the 20 classical actuator disk theory (Sanderse, 2009).

Computing the Annual Energy Production (AEP) of a wind farm
The Annual Energy Production (AEP) is an important metric used to describe a wind farm. The AEP is a statistic. Specifically, it is a mean, as it is a function of the expected power multiplied by the number of hours in a year, The expected power, E[P (ξ)], or the mean of the power, µ P , is defined as where ξ = (ξ 1 , ξ 2 , . . . , ξ n ) is a vector of random variables, which we refer to as the uncertain variables, ρ(ξ) is the joint probability density function of the uncertain variables, Ω is the domain of the uncertain variables, and P is the power produced by the wind farm (Eq. (1)). Common uncertain variables are the wind direction and the freestream wind speed.
The expected power, and hence the AEP, is normally computed as a weighted average, which amounts to the rectangle rule of 10 integration (Sect. 3.1.1). Other uncertainty quantification methods (Sect. 3) can be used to compute the expected value (AEP).
Specifically, we can compute the AEP efficiently by polynomial chaos (Sect. 4).

Uncertainty quantification
Uncertainty quantification (UQ) is the process of (1) characterizing input uncertainties, and then (2) propagating these input uncertainties through a computational model with the goal of quantifying their effect on the model's output. There are 15 many sources of uncertainty in the modeling of a problem, and different classifications of the uncertainties have been proposed (Kennedy and O'Hagan, 2001;Oberkampf et al., 2001;Beyer and Sendhoff, 2007). A common classification is to divide the uncertainty into aleatory and epistemic uncertainties (Oberkampf et al., 2001).
In this work, we consider aleatory uncertainties that arise from the variability in the inputs to our model caused by changing environmental conditions. We describe this input variability as random variables with associated probability distributions. Thus, 20 the first step of characterizing the input uncertainties is concerned with finding the probability distributions that describe the model's inputs. This process is known as statistical inference, model calibration and inverse uncertainty quantification (Smith, 2014). Here, we assume that this step has been completed, i.e., we have distributions that characterize the uncertain inputs (Sect. 5.1). We focus on the second step of propagating the input uncertainties to find the statistics that describe the output.

25
The goal of uncertainty propagation methods is to compute the statistics that describe the effect of the uncertain inputs on the model output. There are several methods to propagate the uncertainties and compute statistics (Le Maître and Knio, 2010;Smith, 2014), where each method has its advantages/disadvantages depending on the type and size of the problem. The most common methods are sampling or Monte Carlo methods (Caflisch, 1998). Other methods include direct integration methods and stochastic expansion methods. Direct integration methods are the currently used method to compute the AEP of a wind farm; we briefly describe them below (Sect. 3.1.1). We describe the stochastic expansion method of polynomial chaos in detail in Sect. 4, and compare the different methods to propagate the uncertainties to compute the AEP in Sect. 6.2.

Direct numerical integration (Rectangle rule)
As the name implies, this method numerically evaluates the integrals in the definition of the statistics. The integrals to evaluate 5 for the mean (or expected value) and the variance are 10 where R(ξ) is the model output and ξ the uncertain input variables. The random vector ξ = (ξ 1 , ξ 2 , . . . , ξ n ) with joint probability distribution ρ(ξ) describes the input variability over the domain Ω. Each random input can follow a particular distribution ρ i (ξ i ). For the case of independent random variables, the joint distribution is the product of each univariate distribution In this work, we will assume that the random input variables are independent. For a description of methods for dealing with cases when the variables are dependent see Padrón (2017). 15 There are many quadrature methods to evaluate integrals (Ascher and Greif, 2011). We describe the rectangle rule, as this is what is currently used in the wind farm community to compute the AEP.
Rectangle rule. The rectangle rule, or mid-point rule, is the simplest and most straightforward quadrature method. To approximate the mean or expected value, A simple improvement is to integrate the density exactly within each subinterval This is easily done as the density is known. This modification is helpful for a small number of evaluations m. We will use this modified rectangle rule and simply refer to it as the rectangle rule.
4 Polynomial chaos 5 Polynomial chaos (PC) is the name of an uncertainty quantification (UQ) method that approximates a function with a polynomial expansion made up of orthogonal polynomials. This function has random variables as inputs, and we are interested in the effects of the random (uncertain) inputs on the output of this function. Statistics of the output can describe the effects of the inputs. We use the polynomial chaos method to efficiently compute these statistics and the gradients of these statistics.

Polynomial chaos expansion
Let R(ξ) be a function of interest that depends on the uncertain variable ξ. We can approximate the function by a polynomial The approximate responseR(ξ) is a polynomial of order p. Usually, the larger the polynomial order the closer the approximation is to the true response R(ξ).
For the case of multiple uncertain variables ξ = (ξ 1 , ξ 2 , . . . , ξ n ) and using a multi-index i = (i 1 , i 2 , . . . , i n ), we write the 5 multi-dimensional polynomial approximation as The multi-dimensional basis functions Φ i (ξ) are given by products of the 1-dimensional orthogonal polynomials When the uncertain variables are independent, the multi-dimensional basis functions are also orthogonal. The values of the 10 elements i j of the multi-index depend on how the expansion is truncated, i.e., on how the index set I p is defined. There are two common ways in which to define the index set: total-order expansion and tensor-product expansion.
In total-order expansion a total polynomial order bound p is enforced: And in tensor-product expansion a per-dimension polynomial order bound p j is enforced An example showing the multi-dimensional basis polynomials Eq. (14) for both the total-order expansion and tensor-product expansion can be found in Padrón (2017). The tensor-product expansion is the preferred approach when the coefficients are computed with quadrature (Sect. 4.2.1) because of increased monomial coverage and accuracy (Eldred and Burkardt, 2009).
The total-order expansion is the preferred approach when the coefficients are computed with regression (Sect. (4.2.2)) because 20 it keeps the sampling requirements lower (Eldred and Burkardt, 2009).

Mean and variance from the polynomial chaos expansion
The mean and variance of the function of interest R(ξ) are a function of the coefficients α i of the polynomial chaos expansion.
The statistics are obtained by substituting the polynomial chaos expansion Eq. (13) into the definitions of the mean Eq. (6) and variance Eq. (9), and by integrating the expansion and simplifying using the orthogonality of the polynomials. 25 The mean is the zeroth coefficient and the variance is the sum of the product of the square of the coefficients-excluding the zeroth coefficient-with the inner where the inner product is defined as Φ 2 i (ξ) = Ω Φ i (ξ)Φ i (ξ)ρ(ξ) dξ and 0 is the first multi-index-the one with all zero elements.

Calculating polynomial chaos coefficients
The coefficients of the polynomial chaos expansion Eq. (13) can be calculated via quadrature or by linear regression.

Quadrature
To obtain the coefficients of the polynomial chaos expansion 10 via quadrature, we take the inner product of both sides of Eq.
Making use of the orthogonality of the polynomials and solving for the coefficients in Eq. (20), we obtain where the domain Ω is the Cartesian product of 1D domains Ω j for each dimension, Ω = Ω 1 ×· · ·×Ω n , and ρ(ξ) = n j=1 ρ j (ξ j ) 15 is the joint probability density of the stochastic parameters. The inner product Φ 2 i (ξ) is known analytically or inexpensively computed. Thus, most of the computational expense in solving for the coefficients resides in evaluating the model R(ξ) in the This integral is solved with quadrature (numerical integration). Note that the zero coefficient in Eq. (21) reduces to the definition of the mean 20 which the direct numerical integration methods attempt to compute directly (Sect. 3.1.1).

Regression
To obtain the coefficients of the polynomial chaos expansion Eq. (13) via regression, we construct a linear system and solve for the coefficients α that best represent a set of responses R. The set of responses is generated by evaluating the model at m realizations of the uncertain vector ξ. The m uncertain vectors are most commonly obtained by sampling the density of the uncertain variables (Hosder et al., 2007).
Each row of the matrix Φ contains the orthogonal polynomials Φ j evaluated at a sample ξ i The size of the m × n matrix is determined by the number of samples m and by how the polynomial chaos expansion is truncated (Sect. 4.1) which results in n terms. It is common to specify a total order expansion along with a collocation ratio cr = m/n to determine the number of samples m. The collocation ratio determines if the system is overdetermined cr > 1 or underdetermined cr < 1.
For overdetermined systems, the most popular method (and the one we use) to estimate the coefficients is least squares, in 10 which we pick coefficients α = (α 0 , α 1 , . . . , α n−1 ) that minimize the residual sum of squares For underdetermined systems, solving a regularized least squares problem is preferred (Doostan and Owhadi, 2011).
For a given number of samples m in the linear system Eq. (24) we can use cross-validation (Hastie et al., 2009) to pick the best polynomial order n to approximate the response.

Gradients of statistics with polynomial chaos
Let R(ξ, x) be a function of interest that depends on uncertain variables ξ and also on design variables x. We assume independence between the design and uncertain variables 3 . Now the polynomial chaos expansion-over the uncertain variablesbecomes 20 This expansion is only valid for a particular design vector-the coefficients α i (x) are a function of the design variables.
Therefore, the statistics are also a function of the design variables. Specifically the mean and the variance are

Gradients of the statistics with polynomial chaos
We want to know the gradients of the statistics with respect to the design variables, and we proceed to derive them below. For simplicity, we drop the subscript from the statistics µ R = µ, the explicit variable dependence R(ξ, x) = R, the bolded notation, and we use the following notation for the gradient df dx ≡ ∇f . The gradient of the mean from Eq. (27) is and the gradient of the variance from Eq. (28) is Both, the mean and the variance gradients, depend on the gradient of the coefficients dα i dx .

Gradients of the coefficients 10
The gradient of the coefficients can be computed with quadrature or regression, similarly to how the coefficients can be calcu- Quadrature. We start from the equation for the coefficients Eq. (21) and take the gradient to obtain Replacing this equation into the gradient of the mean Eq. (29) we obtain And replacing Eq. (31) into the gradient of the variance Eq. (30) we obtain To obtain the gradients of the statistics with respect to each design variable we need to evaluate the multi-dimensional integral containing dR dx . The integral is evaluated with quadrature (Sect. 4.2.1) and requires the computation of the gradient of Regression. We start from the linear system Eq. (23) and take the gradient to obtain To solve for the gradient of the coefficients, we solve the linear system one column at a time of the dα dx matrix with the 5 corresponding column of the matrix of the gradients dR dx . The linear system for the multiple right hand sides can be solved with the methods described in Sect. 4.2.2.
Again, the gradient of the mean Eq. (29) and the gradient of the variance Eq. (30) are a function of the gradient of the coefficients. Thus, to obtain the gradient of the mean take the first row of the dα dx matrix; and for the gradient of the variance use the gradient of the coefficients from all the other rows.

Gradients of the statistics by direct numerical integration
Similarly to computing the mean and variance with direct numerical integration (Sect. 3.1.1), we can also compute the gradients of the mean and variance directly with numerical integration by differentiating the definitions of the mean, Eq. (6), and the variance, Eq. (9). 15 Here we describe the details of computing AEP in the wind farm. We first describe and discuss the inputs, which are the wind direction and wind speed (uncertain variables) (Sect. 5.1) and the wind turbine positions (design variables) (Sect. 5.2).

The wind farm layouts
To showcase the results, we will focus on four representative layouts: Grid, Amalia, Optimized and Random (Fig. 4). The where the diameter represents the rotor swept area. The Grid layout fits in a box of equivalent area to that of the Amalia layout.
The Amalia layout, which is grid-like, is that of the Princess Amalia wind farm located 23 km off the coast of the Netherlands.
Export as pdf to maintain quality. Then open in Preview and crop to only get the figure.   Maybe make the text weight different, from medium to regular. In reality, the :: 60 : turbines in the Princess Amalia wind farm are the Vestas V80 model. For each of the layouts in our study, we use the NREL 5-MW reference turbine (Jonkman et al., 2009), as this turbine is similar the the Vestas V80, and has an open source design.

Convergence metric -the average AEP error
We use an ensemble of 10 AEP results to compute the average AEP error. The average AEP error allows us to better illustrate the 5 differences between the different methods used to compute the AEP and to avoid drawing conclusions from one-off solutions.

Baseline AEP
We take as the baseline or true AEP the AEP computed with 200,000 Monte Carlo samples. We picked 200,000 MC samples 15 to ensure the 99 % confidence interval for the true AEP was smaller than 1 % of the computed AEP value for all layouts.
We consider an AEP within 1 % of the baseline AEP to be accurate, and we will use it as a reference for the results. AEP predictions of real wind farms usually have an error of 10-20 % ( Barthelmie et al., 2007;Briggs, 2013) due to uncertainty in wind conditions and to the errors of wake models, thus resolving the AEP to less than 1 % is unnecessary.

Methods to compute the AEP 20
Here, we provide the details of the methods used to compute the AEP, as well as, the abbreviations for the methods:  (Caflisch, 1998) For the quadrature-based methods (rect, PC-Q), we use tensor product quadrature to compute multi-dimensional integrals. We use the same number of points for each dimension because we did not see any benefit in favoring a particular dimension. For PC-Q, we use Gaussian quadrature. For Monte Carlo, we use the traditional Monte Carlo method, i.e., random samples. For

25
PC-R, we use Latin-Hypercube sampling (McKay et al., 1979) to generate the samples needed to construct the linear system.
We solve the linear system with least-squares. For a given number of samples, given by the total polynomial order p, we use 10-fold cross-validation to find the least-squares best fit from polynomials of total order 1 up to total order p. The sampling methods and the polynomial chaos methods we use are implemented in the open-source DAKOTA toolkit (Adams et al., 2017). For all layouts as the speed increases the power increases until the wind farm reaches its rated power 300 MW.

Results
The Annual Energy Production (AEP) of a wind farm is obtained by integrating the power output of the farm over the different wind conditions. Thus, we first characterize the power output of the wake model for different input conditions (Sect. 6.1). Next, we focus on the convergence of the AEP (Sect. 6.2), and then on the wind farm layout optimization problem to maximize the AEP (Sect. 6.3).

Power response as a function of the uncertain variables
The power production, computed with the Floris wake model, for the four wind farm layouts (Sect. 5.2) as a function of both the wind direction and the wind speed is shown in Fig. 5. The wind direction is measured from North and increases clockwise.
6.2 AEP convergence: polynomial chaos vs. rectangle rule 15 We consider the AEP as a function of two uncertain variables: the wind speed and the wind direction. The AEP is usually considered a function of these two variables. We compare the convergence of the AEP for the different methods to compute the AEP: polynomial chaos (based on quadrature and regression), the rectangle rule and Monte Carlo (Fig. 6). The polynomial chaos based on regression (PC-R) performs the best for all layouts. It is followed, by the polynomial chaos based on quadrature (PC-Q) and the rectangle rule, which perform similarly. The worst performer is Monte Carlo. The slow convergence of statistics 20 with Monte Carlo is well known. Monte Carlo will start to outperform the other methods when the AEP is a function of a large number (5)(6)(7)(8)(9)(10) of uncertain variables, as it does not suffer from the curse of dimensionality.
The superior performance of the polynomial chaos based on regression, especially for the grid like-layouts (Grid and Amalia), is due to the following: the polynomial chaos fit based on regression does not chase all the high-frequency oscillations in the power response (Fig. 5). : , ::: i.e., :: it :::::::: smooths ::: out ::: the :::::::: response. : The PC-R fit is usually not higher than an eight 25 total-order polynomial (Sect. 4.1). Whereas, the PC-Q order fit is higher, as it is directly proportional to the number of samples per dimension 7 . A downside of the PC-R being able to predict the mean (AEP) accurately, is that it can underpredict the true variance (standard deviation) of the response (Fig. 7). Usually, the standard deviation of the power response (energy) over a year is not considered as a function to optimize. A common objective in wind farm optimization is to maximize the total amount of energy produced over a year independent of the variability in the power production over the year. For a wind farm,  In what follows we will compare the PC-R (the best performing method) with the rectangle rule (the method currently used in practice) to quantify the reduction of samples needed to compute the AEP accurately. Also, we will sometimes use polynomial chaos to refer to polynomial chaos based on regression. Figure 8 only keeps PC-R and the rectangle rule results from Fig. 6, and in addition, the figure shows the average AEP error computed with 10 and 100 sets of samples for each method. For the rectangle rule, there is hardly any difference between the average AEP error computed with 10 or 100 sets. 5 For the PC-R method, the average error with 100 sets shows a smoother convergence. In general, averaging the AEP error over 10 sets of samples is enough to ::::::: minimize ::: the :::::: AEP's :::::: sample ::::::: location ::::::::: sensitivity ::: and :: to : clearly see the differences between the methods used for computing the AEP (Sect. 5.3).
In Fig. 8, we see that the PC-R convergence curve is consistently below the rectangle rule curve, i.e., PC-R has a smaller error for the same number of samples or the same error for a smaller number of samples. Using the 1 % average AEP error wind speed. Note that the PC-R response is biased for the Grid and Amalia layouts. The PC-R underpredicts the true STD at the expense of computing the mean (AEP) more accurately (see Fig. 6).
6.3 Wind farm layout optimization

Optimization problem
The objective of the wind farm layout optimization is to maximize the AEP (Sect. 2.3) by changing the position of the wind turbines. We assume a fixed number of turbines, 60, of the same type (NREL 5-MW (Jonkman et al., 2009)) and constrain the turbines to stay within a given area and with a minimum separation between them. This objective and constraints result in a 5 nonlinear optimization under uncertainty problem with deterministic constraints: where S i,j is the distance between each pair of turbines i and j, and D is the turbine diameter. The normal distance, N i,b , from each turbine i to each boundary b is defined as positive when a turbine is inside the boundary, and negative when it is outside of the boundary. The boundary is the convex hull of the Princess Amalia layout-a fourteen-sided convex polygon (dashed-line We solve the optimization problem with the gradient-based sequential quadratic programming optimizer SNOPT (Gill et al., 2005). We use OpenMDAO (Gray et al., 2010) and its wrapper for pyOptSparse (Perez et al., 2012) to call SNOPT from 15 Python. We scale the variables, constraints and the objective to make them of order one, and set the tolerances to 1 × 10 −4 per the function (AEP) precision.

Optimization results
The AEP in the optimization objective is a function of two uncertain variables the wind direction and wind speed, with probability distributions given in Sect. 5.1. We solve the optimization under uncertainty problem, Eq. (38), for two different starting 20 layouts: Amalia and Random. We compute the AEP (objective) with different precision (number of samples) and with different methods: the rectangle rule, and polynomial chaos based on quadrature and regression. For each method, we run 10 optimizations, where each optimization uses different sample points to compute the AEP (see Sect. 5.3). The 10 optimizations enable us to get a better understanding of which method is better at finding layouts with high AEP and to avoid drawing conclusions from one-off local optima. 25 The results of the optimizations are reported in Table 2. The table contains the statistics of the AEP of the optimal layouts for the 10 optimizations of each method. The values of the AEP reported in Table 2 are computed with 200,000 Monte Carlo samples. The same 200,000 samples-wind direction, wind speed pairs-are used to test the optimal layouts obtained by each method. Table 2. The AEP statistics of the optimized layouts. We generate the statistics for each method from a set of 10 different optimizations.
Polynomial chaos based on regression on average produces the best optimums when the different methods to compute the AEP use roughly the same number of samples (∼ 225). Furthermore, PC-R optimums are comparable with the optimums obtained with the rectangle rule that used 625 samples to compute the AEP. Starting the optimizations from an already good layout (Amalia), will in general, find layouts that are better than starting from a bad layout (Random Polynomial chaos based on regression on average produces the best optimums when the different methods to compute the AEP use roughly the same number of samples (∼ 225). Furthermore, PC-R optimums are comparable with the optimums obtained with the rectangle rule that used 625 samples to compute the AEP.
Starting the optimizations from an already good layout (Amalia) will, in general, find layouts that are better than starting from a bad layout (Random). When considering all the methods, on average, the optimal layouts starting from the Amalia have 5 an AEP that is 0.62 % higher than those starting from the Random layout. However, starting from random layouts gives the turbines more freedom to move around the design space with the potential of finding novel and better layouts.
The optimal layouts with the maximum AEP for each method and starting layout are shown in Fig. 10. For the optimums obtained starting from the Amalia layout, we see that the distribution of the turbines is similar to the one of the starting layout.
And for the optimums starting from the Random layout, we see that the optimal layouts obtained are less grid-like and produce 10 more novel layouts. For each method, as shown in Fig. 10, the optimal layouts starting from the Amalia have a higher AEP than those from the Random. But a random start has the potential to find better optimums as the turbines explore more of the design space (Fig. 11).
A typical optimization history as a number function calls is shown in Fig. 12. A function call is the AEP computation which requires hundreds (the number of samples specified for each method in Table 2) of calls to the wake model for the wind farm 15 resulting in tens of thousands of calls to the wind farm wake model per optimization. If the gradients of the AEP are computed with a first-order finite difference, 120 (the number of design variables) times as many wind farm wake model calls would be needed per optimization. The optimization starting from the Amalia layout (Fig. 12a) converges faster than that starting from the Random layout (Fig. 12b).  Table 2. The first row shows the optimums obtained starting from the Amalia layout; we see that the distribution of the turbines is similar to the one of the starting layout. The second row shows the optimums starting from the Random layout; we see that the optimal layouts obtained are less grid-like and produce more novel layouts. For each method, the optimal layouts starting from the Amalia have a higher AEP than those from the Random.
Export as pdf to maintain quality. Then open in Preview and crop to only get the figure.   displacement of all turbines measured in turbine diameters is 315 and for the Random case 2790. In general, most optimizations we performed behaved similarly to those shown in Fig. 12. The optimizations starting from the Amalia layout run for fewer iterations (usually on the order of 100) with less design variable change than those starting from the Random layout. Most of the optimizations we ran either successfully converged or were on their way to satisfying the convergence criteria before they reached the time limit we set to save computational cost.

5
To properly compare the results obtained by the different methods used to compute the AEP in the optimization, we should compute the AEP of the optimal layout with a method that was not used for the optimization. We use 200,000 Monte Carlo samples-wind direction, wind speed pairs-to test the AEP of the optimal layouts. For the optimal layouts in Fig. 10, with the Amalia starting layout, the average improvement in AEP over the starting layout for all the methods is 0.14 % when measured with the Monte Carlo samples and 1.74 % when measured by the method used in the optimization. For the optimizations used in the optimization (see Table 3). The reporting of the AEP computed by Monte Carlo explains why some of the results reported in Table 2 and Fig. 10 have a smaller AEP than its starting layout.
The improvements starting from the Amalia layout are similar to those found by Gebraad et al. (2017) for turbine position optimization. The large improvements in AEP for the layouts starting from the Random layout show that the optimizer can find good layouts from a bad starting position. Thus, to search for the best layout, we would perform many optimizations with 5 random starts.
From all the optimizations runs, we conclude that computing the AEP with polynomial chaos based on regression produces the best layouts-around 1 % higher AEP when starting from a random start-for the same number of samples. This is because PC-R computes the AEP more accurately than the other methods for the same number of samples (Sect. 6.2). Furthermore, the optimization with PC-R finds optimums comparable to those found with the rectangle rule that used roughly three times 10 as many simulations 8 . We found that starting an optimization from a good layout will, in general, find better optimums than starting from a random layout, but starting from the random layout can lead to novel layouts and possibly better layouts as the turbines explore more of the design space. Finally, to properly compare between methods an ensemble of optimizations should be used and evaluated in the same way.
The layout of the wind farm influences the convergence of the AEP because the layout has a significant effect on the power output of the farm as the wind conditions vary . :::: (Fig. :: 5). : We considered four representative layouts: Grid, Amalia, Optimized and Random. The power response of the grid-like layouts (Grid and Amalia) has large oscillations caused by the large drops in power that occur when rows of wind turbines are aligned with particular wind directions. Because of the larger variability in 20 the power response, the grid-like layouts require more simulations than the non-grid-like layouts (Optimized and Random) to converge the AEP, independent of the method used to compute the AEP..
Competing interests. The authors declare that they have no conflict of interest.