Constrained optimization of offshore wind turbine positions under uncertain wind conditions with correlated data

Offshore wind power production is subject to substantial uncertainty due to variability in wind conditions. Planning of new wind parks should take these uncertainties into account by means of stochastic modeling and uncertainty quantification. Wind speed and wind direction exhibit dependence that needs to be properly modeled to yield reliable uncertainty estimates. In particular, the dependence is strong when data are averaged over a short time-horizon, e.g., on a monthly rather than an annual basis. In this work we introduce a stochastic model for wind speed and wind direction using Rosenblatt transformation and an empirical model. This allows efficient numerical quadrature that replaces thousands of Monte Carlos samples by a few tens of model evaluations at quadrature points in parametric space. A robust design formulation for maximizing power production honoring dependent data is proposed and demonstrated on data from the North Sea site Sørlige Nordsjø II.


Introduction
The European Union has experienced a 40-fold increase in the offshore wind power capacity, from 620 MW in 2004, to almost 25,000 MW in 2020 [19].With the increasing importance of offshore wind energy resources, it is essential to achieve efficient utilization of sites by optimal location of the wind turbines, as available area could be limited, due to, e.g., licensing agreements, integration to off-shore power-grids, and subsurface topography.In particular, efficient optimization under uncertainty is a prerequisite for real world applications where wind speed and directions are dynamically changing.
Wind farm layout optimization aims at determining those wind turbine positions that maximizes energy production from a wind field while accounting for the reduction in wind speed downstream of every turbine via wake models [7,9].The wake model and its interference from nearby wind turbines makes the problem non-convex and challenging to solve, and it could also create a non-smooth objective function.Common solution strategies include heuristic techniques such as genetic algorithm and particle swarm optimization [17].However, formulation of continuous wake models [14] render gradient-based solution strategies possible.While none of these strategies may lead to the globally optimal solution, both can yield relevant optimized solution estimates.
Optimization that considers only the most likely wind direction, may lead to a suboptimal solution that does not attain the maximum possible total annual energy production (AEP) available if all scenarios are considered.Hence, the total distribution of wind over time needs to be taken into account.This, in combination with wind turbines of different sizes and designs, leads to non-trivial optimization problems.Numerical optimization of offshore wind turbine (OWT) farm layouts under uncertainty is a challenging task, as the computational cost of brute-force approaches is simply the product of the cost of optimization and the cost of uncertainty quantification.Both optimization and uncertainty quantification rely on repeated evaluation of a computational OWT model, with a cost that increases quickly with the number of turbines.To target large OWT configurations, it is therefore essential to use efficient methods that do not require an excessive amount of model evaluations.Optimization under uncertainty in general assumes independent uncertain variables, but wind data are typically highly variable, and wind speed and direction often exhibit strong dependency [4].Ignoring this dependence can result in erroneous results.Uncertain parameters are often modeled as independent, wellknown random variables in wind simulation, but should rather be estimated directly from data, including dependencies between parameters.
Early works on optimization under uncertainty of wind turbine layouts include Monte Carlo simulation on a discretized physical domain [10].Uncertainty quantification and correction of the wake model as a function of yaw misalignment was investigated in [6].Monte Carlo simulation with dependent variables for wind farms was carried out in [5].Optimization under uncertainty was considered in [13], where generalized polynomial chaos was introduced to represent uncertainty assuming independence in wind speed and direction.Sensitivity analysis of energy production and levelized cost of energy for a number of uncertain input parameters was performed for the three wind farms Dan Tysk, Horns Rev 1, and Sandbank using sparse grid stochastic collocation and Quasi-Monte Carlo methods in [15].
In this work, we consider a non-convex and smooth objective function in the optimization formulation using the continuous wake model of Park and Law [14], as described in Sect. 2. In Sect.3, we suggest methods for accelerated uncertainty quantification for dependent data based on transformations and numerical quadrature.As a test case for the OWT layout optimization, we use wind data and wind turbine data relevant for the Sørlige Nordsjø II site.A robust design formulation for OWT layout optimization under uncertainty is presented in Sect.4, and numerical results in Sect.5, followed by concluding remarks in Sect.6.

Wind farm and wake models
Consider a wind turbine site with air density ρ, ambient wind speed v and wind direction θ.The site is populated with N T turbines, each of which experiences a local wind speed v n and produces a time-averaged power P n (n = 1, . . ., N T ), given by where A is the swept rotor area of the turbine and C p is the power coefficient with maximum C * p .The maximum power coefficient is upward limited by the Betz's limit, that is, C * p ≤ 16/27.The turbine only operates above the cut-in wind speed v cut-in and below the cut-out speed v cut-out .For ambient wind speeds above the rated speed v rated and less than v cut-out , the induction is modified to maintain a constant power production, resulting in an effective speed equal to v rated .Here, we will relate the induction factor, α, to the power coefficient, C p , through the relation [2] C p (α) = 4α(1 − α) 2 , and the wind turbine power model in Eq. (1).For v n < v cut−in and v n > v cut−out we have that α = 0.For v cut−in ≤ v n ≤ v cut−out the power coefficient is a constant, C * p , and In this work we use parameters given by the IEA 15-MW offshore reference wind turbine [1] for all turbines, as specified in Table 1.
The local wind speed v n can be smaller than the ambient wind speed v as a result of wake effects from other wind turbines.In general, where the wind speed deficit factor, δv n , determines the reduction of the freestream wind speed due to effects from other wind turbines.The deficit factor on turbine n, resulting from a single upstream turbine m, is typically a function of the downstream-wake distance d mn and downstream-wake radial distance r mn as well as the (upstream turbine) induction factor α. The exact functional form is determined by the wake model.There are numerous analytical wake models, e.g., Jensen's model [7] and Larsen's model [9].In this work, we consider the continuous wake model from Park and Law [14].Here, δv n has the form of a scaled Gaussian where d ≥ 0 is the downstream distance and r is the radial distance from an OWT, and R is the radial extension of the wake.The radial extension function is defined as where R 0 is the radius of the actuator disk and κ is the wake expansion rate and is set to 0.05 in this work.We use the wake interference from Park and Law [14], where the velocity used in the power production function is defined using the wind speed deficit factors.For a given turbine n the influence from all upstream turbines is calculated through the wind speed deficit factor, averaged over turbine n's actuator disc with area A. This average is calculated using a Gaussian quadrature approach where the integral over the actuator disc is approximated as a sum over a set of N A q quadrature points, r q n , and weights w q where q = 1, 2, ..., N A q .Given an upstream wind turbine m at position r m , we can calculate the downstream and radial distances from that turbine to each of the integration points d q nm and r q nm at the downstream turbine's actuator disc quadrature points.The mean wind speed deficit factor is now approximated as w q δv(d q nm , r q nm , α m ).
Finally, the total wind speed deficit factor for turbine n is calculated by the sum of squares of the deficit factors for all upstream turbines [8] so that When the total wind speed deficit factor is found, we know from Eq. ( 2) the velocity and we obtain the total power produced on the site as

Stochastic wind models and uncertainty quantification
Relevant wind data are often measured over time and can be interpreted as a time series.Despite the obvious temporal dependence in most wind data, in the current work it is convenient to consider wind data as homogeneous, i.e., wind speed and directions are assumed to follow some unknown distributions that are independent of time.This is motivated by the fact that we wish to study the overall variability in wind conditions, but not the temporal changes per se.Hence, we follow the common practice of estimating a single joint distribution for wind speed and wind direction based on data collected over time at a given site.

Model for dependent wind data
Let θ and v be dependent random variables representing wind direction and speed with cumulative distribution functions F v and F θ , respectively, and let U = (U 1 , U 2 ) be independent uniform random variables on the unit interval.Using the Rosenblatt transformation [16], we have the distributional equalities or, equivalently Hence, we can express quantities of interest, e.g., P , as functions of independent uniform random variables, Mean values, variances, and failure probabilities can all be expressed as integrals over some function g of the quantity of interest, i.e.,

E(g(P
where g(P ) = P corresponds to the expectation, g(P ) = P 2 the second moment (used for variance), and g(P ) = 1 P ∈E the probability of occurrence of some sets of events E, for instance a failure probability.We note that the second integral in Eq. ( 3) can be evaluated by, e.g., standard numerical quadrature due to the independence of the variables U, in contrast to the first integral, where the dependence between v and θ excludes direct application of standard methods that assume independence.

Numerical quadrature
A numerical approximation of Eq. ( 3) by means of standard Monte Carlo integration is computationally expensive, in particular if many evaluations are needed as part of an optimization procedure.
Using numerical quadrature can significantly alleviate the computational cost, in particular if the integrand is a smooth function of the variables, so that, e.g., Gauss-type quadrature can be employed.In the current work, the power is a piecewise smooth function of the input variables, and a compound quadrature rule will be introduced, based on the following considerations.If the wind speed is below or above the interval [v cut-in , v rated ] with some non-negligible probability, then the power (Eq.( 1)) is a non-smooth function of v n .With this a priori information, in the situation where v n = v (i.e., no wake effects), we may assign probability point masses to the events v < v cut-in , v rated ≤ v < v cut-out and v > v cut-out , and use a quadrature rule, e.g., Gauss-type quadrature for [v cut-in , v rated ].However, when wake effects are present, the effective wind speed an individual turbine experiences can be less than v rated even when v > v rated .In this case, the conditional expectation on [v rated , v cut-out ] should be approximated by numerical quadrature and not a point mass.To accommodate for that, we will employ a compound quadrature rule: Gauss quadrature for v cut-in ≤ v < v rated and v rated ≤ v < v cut-out , and a discrete probability for all other events.The integral in Eq. ( 3) can then be approximated by where P q = P (U q ), with quadrature nodes U q and weights w q for q = 1, . . ., N q .

Wind data: Sørlige Nordsjø II
Hourly wind data for the entire period 2016-2018 are retrieved from the NORA3 dataset [18,3] for a location within the field Sørlige Nordsjø II (SN II), located in the North Sea west of Denmark.A random model as outlined in Sect.3.1 is fitted to the data using the Rosenblatt transformation method from the software pyvinecopulib [11].The wind data are visualized in Fig. 1, and compared to a wind rose generated by sampling from the fitted model, and a subset (every 20th observation) of all wind data where speed and direction are treated as independent variables.The dependence in data is hard to assess from the windroses only.As an illustration of the dependence of the current data set, consider the function g = v 3 θ which is somewhat reminiscent of Eq. ( 1).The relative difference in statistics of that function with dependent data compared to assuming independence of data gives an indication of the effect of dependence between wind speed and direction on, e.g., power production.In particular, we investigate the effect on the mean value µ = E(g) (first moment), variance σ 2 = E(g − µ) 2 (second central moment), skewness E(g − µ) 3 /σ 3 (standardized third moment), and kurtosis E(g − µ) 4 /σ 4 (standardized fourth moment).The relative errors in mean, variance, skewness, and kurtosis of the function g, assuming independence compared to accounting for the dependence displayed in actual data, are -1%, 1.7%, -8.5 %, and -7.3% respectively (see all data in Fig. 2).While these errors are not very large, it is interesting to consider the errors for a model that assumes stochastic homogeneity on a shorter time-span, e.g., a model that assumes data-driven distributions for each month.Figure 2 shows the relative errors in mean, variance, skewness, and kurtosis for the same test function g = v 3 θ.Clearly, the dependence is more pronounced on a monthly basis compared to averaging over years.The relative error in the mean of the function is modest, but the errors displayed in higher-order moments are often significant and suggest that the true distributions are quite far from independent.

Robust design of OWT
A robust design of an OWT layout should maximize the power production while ensuring that the risk for undesired performance is as small as possible.The task of robust design optimization under uncertainty is to formulate and solve a problem that takes both these aims into account.Consider an objective function g( P (x; U)), where x are design variables confined to some set x ∈ X and U are random variables.In this work, x will always denote spatial locations of the OWTs.With an appropriate choice of g and by applying a function I U that involves expectations with respect to the random variables U (i.e., integration in random space), we obtain the robust optimization problem min x I U (g( P (x; U))), subject to constraints on x.Several choices of g and I U may be relevant, for instance letting and g = 1 {P ≤q} that leads to minimization of the probability that the power is lower than some threshold value q.Here, we will consider a mean-value penalty optimization with mean, µ, and higher-order central moments m k (k = 2, . . . ) of the power, A negative β 1 and positive β k (k = 2, . . ., M ) leads to a formulation that favors large expected power with small variability.

Optimization method
The optimization problem Eq. ( 6) is solved using Sequential Least Squares Quadratic Programming (SLSQP).This method for nonlinear optimization is suitable as it can handle complex constraints that are twice differentiable, which is the case for OWT layout optimization with, e.g., minimum inter-distance constraints.SLSQP relies on both the objective function and its gradient with respect to the design variables, i.e., the physical coordinates of all turbines.Numerical Gauss-Legendre compound quadrature as described in Sect.3.2 is used both to compute the objective function and its gradient.

Gradient calculations
For the gradients in robust design, let P q denote power quadrature point, and w q the corresponding quadrature weight, q = 1, . . ., N q .With µ = E(P ) and and 2P q ∇P q w q − 2µ∇µ, By the chain rule, we obtain ∇σ from ∇σ = ∇σ 2 /(2σ).This expression is clearly undefined if σ = 0.For this case, since there is no variability, we have that P q = µ for all q = 1, . . ., N q , and we have that Thus, if the standard deviation is 0, so is its gradient.

Numerical optimization of OWT layouts
Consider a rectangular OWT site of size 2400 m x 2000 m subject to the uncertain wind conditions determined by the stochastic model and data from Sørlige Nordsjø II.We solve the robust optimization problem Eq. ( 6) on the rectangular domain where turbines are subject to the constraints x m − x n 2 ≥ 2R 0 for m, n = 1, . . ., N T , m = n, i.e., the turbines centers must be at least two rotor radii apart.We have chosen this to be the minimal distance to avoid collisions between turbines.In practice, a larger inter-turbine distance may be required due to constraints given by other considerations, e.g., turbulence effects.We use M = 2, The integrals in Eq. ( 6) are approximated with the compound quadrature rule described in Sect.3.2 with N q = 81 quadrature points (one point mass for regimes of zero Fig. 3 shows the robust optimized OWT locations for increasing number of turbines.For reference, the expected wind field due to wind condition variability and wake effects, is also shown.As a comparison, Fig. 4 shows the optimized OWT layout under the mean wind speed and mean direction, which may be different from the most common speed and direction.Note, that in this case only one wind speed and wind direction is used in the wind farm modelling.The wind field also shows this, as the wake from each OWT is clearly visible in Fig. 4. The mean wind field in the robust optimized case is the spatial manifestation of the wind rose in Fig. 1, and has a more isotropic structure than the simple wake from just one wind direction (although the wind rose is far from being perfectly isotropic).For a perfectly isotropic system we would have that the mean wind speed would only depend on the radial distance from all turbines.Intuitively, if we only consider energy production optimization in a limited area, we would assume that for a perfectly isotropic system a good strategy would simply be to place each turbine as far as possible from all other turbines, which would suggest that we end up with a regular pattern of wind turbines.On the other hand, if we only have one wind direction each turbine should stay out of each others wake which would suggest that they, ideally, should be placed on a line normal to the wind direction.We see these trends in Figs. 3 and 4. In Fig. 3    pattern, but still a much more homogeneous density of OWTs than for the ones in Fig. 4. We must also keep in mind that the standard deviation is also part of the objective function used in the minimization in the case for the robust design.In Fig. 4 we see that the OWTs are aligned with the diagonal, normal to the mean wind direction, for N T ≤ 12.We note a topological change in OWT positions for higher turbine densities where we go from the diagonal of the domain to the boundaries.On the diagonal of the domain we can place 13 turbines so this transition can be related to the maximum density of turbines along the most ideal position.
Figure 5 shows how wind farms optimized according to the robust design and farms optimized with respect to mean wind conditions operate under conditions given by the wind data presented in Sect.3.3.All measures are calculated using the full wind data set.From Fig. 5(a), we see that the relative expected efficiency E[P (N T )]/(N T •E[P (1)]) is generally lower in farms only optimized with regards to the mean wind conditions at this site.In this figure, we also observe that the relative expected efficiency in the mean-wind optimized farms does not decrease monotonically as a function of turbines N T in the farm.We see for instance that this quantity increases from N T = 12 and N T = 15.This is simply because, in the actual wind conditions experienced at the site of operations, a larger number of the turbines ends up in the wake of other turbines in the farm layout presented in Fig. 4(a) compared to the layout presented in Fig. 4(b).The coefficient of variation, i.e., the ratio between standard deviation and mean value, of the power production is high for both robust design and optimization under mean conditions, and increases with the number of turbines, see Fig. 5(b).The robust design results are not sensitive to a higher penalty on the standard deviation, since the effect of uncertainty in the weather conditions can hardly be reduced.The locations that minimize the objective functions are typically not unique.The physical space that corresponds to locations satisfying minimum objective function values would shrink by increasing the minimum distance constraints.This would affect the numerical results and could possibly lead to local minima, but the numerical cost of the method would not be significantly affected.
For development phase one of the SN II site, the Norwegian Ministry of Petroleum and Energy gives as prerequisite a maximum power production of 1500 MW with a suggested minimum area efficiency of 5 MW/km 2 [12].The areal extent of SN II is 605 km 2 and with efficiency of 5 MW/km 2 about half of the area must be developed.Figure 5(c) shows the area efficiency from our simulations.Based on our model, a robust design simulation with 30 OWTs yields an area efficiency of 58.7 MW/km 2 .This requires development of only 4.2 % of the SN II to meet the requirements from the authorities, which allows for a substantial room for other usages of the marine space in this area.

Conclusions
An efficient method for optimization of wind turbine layout under uncertain wind conditions has been presented.As investigated for the site Sørlige Nordsjø II, wind speed and direction are dependent.The dependence is particularly strong if one is interested in including time dependence in stochastic wind models, i.e., when data are averaged over a shorter time-frame instead of assuming that data are stochastically homogeneous in time.Robust mean-value penalty optimization has been performed, where statistics evaluated by means of numerical quadrature replaces Monte Carlo simulation to reduce the computational cost.
We have demonstrated that this procedure can be used to optimize OWT layouts with use of actual site specific wind condition data.The results, exemplified by Figs. 4 and 3, demonstrate the importance of optimizing the wind turbine positions with respect to the variations in the actual wind conditions at the targeted production site.The simulation examples performed, although using simplified turbine and wake models, indicate the power production, suggested by the Norwegian Ministry of Petroleum and Energy, in the stage-1 development of SN II may be reached by utilizing only a bit above 4 % of the area in SN II.
In this work, we have only considered maximisation of the power production.Investigation of additional effects requires augmentation of the objective function with extra terms.For instance, if one wants to include additional fatigue loads one will need to add the cost of repair to the value of the produced power or the added risk of turbine damage with appropriately assigned weights.

Figure 1 .
Figure 1.Wind roses for Sørlige Nordsjø II for hourly data 2016-2018 (left), wind rose based on equally many (26304) samples of the fitted model (center), and wind data treated as independent observations of direction and speed (right).

A l l d-Figure 2 .
Figure 2. Relative error in statistics for the function g = v 3 θ when assuming independence in wind speed and direction relative to actual dependent data.Hourly data 2016-2018 month-bymonth from Sørlige Nordsjø II.

Figure 4 .
Figure 4. Layout with respect to mean wind conditions.Optimized turbine locations are indicated by circles and initial location by crosses.The colors represent the local wind speed field under mean conditions with the related wake effects.

Figure 5 .
Figure 5. Plots showing measures relating to the overall power production of the wind farm.

Table 1 .
[1]ameters for the IEA 15-MW offshore reference wind turbine[1]used in the wind farm simulations.Using the actuator disc turbine power model presented in Eq. (1), these parameter values result in a maximal power production of 16.09 MW for one turbine., by solving the equation C p (α * ) = C * p .When v rated < v n ≤ v cut−out the induction factor, α, becomes a function of the wind speed, implicitly defined by the equation C p * we see that for N T ≤ 15 shows almost regular patterns, while for N T ≥ 18 we note a more irregular