Stochastic gradient descent for wind farm optimization

It is important to optimize wind turbine positions to mitigate potential wake losses. To perform this optimization, atmospheric conditions, such as the inflow speed and direction, are assigned probability distributions according to measured data, which are propagated through engineering wake models to estimate the annual energy production (AEP). This study presents stochastic gradient descent (SGD) for wind farm optimization, which is an approach that estimates the gradient of the AEP using Monte Carlo simulation, allowing for the consideration of an arbitrarily large number of atmospheric conditions. SGD is demonstrated using wind farms with square and circular boundaries, considering cases with 100, 144, 225, and 325 turbines, and the results are compared to a deterministic optimization approach. It is shown that SGD finds a larger optimal AEP in substantially less time than the deterministic counterpart as the number of wind turbines is increased.

1 Introduction

Wind farms are groups of wind turbines that harness the power in the atmospheric boundary layer to provide renewable energy. When a wind turbine absorbs energy from the air, the air downstream of the wind turbine has reduced power, which often reduces the power production of downstream turbines. This is known as the wake effect (Sanderse, 2009; Hasager et al., 2013). When a new wind power plant is to be constructed, optimal turbine locations are determined using engineering wind farm models (Samorani, 2013; Ning. et al., 2020; Annoni et al., 2018). Turbine positions are optimized to exploit the benefits of the local wind resource while avoiding energy losses from turbine wakes. In the wind turbine placement problem, atmospheric conditions, such as the inflow speed and direction, are assigned probability distributions according to measured data. By propagating these probability distributions through the engineering wake model, the annual energy production (AEP) can be estimated. The AEP is often computed using rectangular quadrature, dividing the relevant speeds and directions into equal-sized bins, then computing the expected AEP as the product of the power and probability of each bin, added together, then multiplied by the number of hours per year. The cost of wind farm optimization generally increases with the number of atmospheric conditions considered during AEP computation, and this expense becomes more extreme as more complex wake models (e.g., Reynolds-averaged Navier–Stokes models) are considered. For example, there are some memory limitations when computing AEP gradients using automatic differentiation with very large wind farms. This has given rise to studies seeking convergence of the AEP, proposing methods such as polynomial chaos expansion (Padrón et al., 2019; Murcia et al., 2015) or Bayesian quadrature (King et al., 2020) to avoid discretizing the input distributions into evenly spaced intervals. In this study, we present an approach for wind farm optimization that estimates the gradient of the AEP using Monte Carlo simulation. This does not require the input to be discretized at all and allows for the consideration of an arbitrarily large number of atmospheric conditions.

Stochastic gradient descent (SGD) is an optimization algorithm commonly used in machine learning when selecting neural network weights (Ketkar, 2017). The algorithm samples the gradient of a stochastic objective, following the mean gradient by a specified distance, then repeating the process, which amounts to optimizing the expected value of the objective. The SGD algorithm is often enhanced to avoid oscillations caused by large changes in the gradient of the objective (Ruder, 2016). This includes methods to reuse previous gradient information (Qian, 1999), dampen oscillations (Riedmiller and Braun, 1993), or incorporate an estimate of the Hessian matrix (Moritz et al., 2016; Byrd et al., 2016; Liu et al., 2018; Najafabadi et al., 2017). Kingma and Ba (2014) introduced the Adam SGD algorithm, which reuses gradient evaluations and dampens oscillations, and which is the basis of the SGD method we propose in this study.

Interestingly, SGD is not often applied to problems with nonlinear constraints, although it can be fruitful to include nonlinear constraints in the context of training a machine learning algorithm. For example, when recognizing three-dimensional pictures of people, it can be useful to impose a constraint that any person's left arm should be close to the same length as their right arm (Márquez-Neila et al., 2017). Many frameworks have been proposed for constrained SGD, including the log-barrier function (Kervadec et al., 2019), penalty functions (Márquez-Neila et al., 2017), blending barrier and penalty functions (Kervadec et al., 2019), and Riemannian geometry (Roy and Harandi, 2017). In this study, we use a penalty term to transform the constrained problem into an unconstrained optimization.

The wind farm layout optimization problem presents a setting where the objective (AEP) can be formulated as being stochastic (e.g., the AEP is derived from a probability density function), while the constraints (e.g., boundaries and minimum turbine spacing) are firmly deterministic. This paper explores the potential benefits of formulating the wind farm layout optimization problem in this way. As part of this, the Adam algorithm is extended to optimize a stochastic objective with deterministic constraints. To the best of the authors' knowledge, this exact algorithm has not been published before.

This study benchmarks the performance of the proposed SGD approach when compared to conventional gradient-based optimization within the TOPFARM framework (DTU Wind Energy Systems, 2023b), considering wind farms with different shapes and sizes. We examine the open-source SLSQP algorithm (Kraft, 1988), which is employed in many engineering frameworks (King et al., 2017; Allen et al., 2020; Wu et al., 2020; Zhang et al., 2022; Zilong and Wei, 2022; Kölle et al., 2022; Clark et al., 2022; Simley et al., 2023) and has been used in previous comparisons of optimization algorithms (Lam et al., 2018; Li and Zhang, 2021; Fleming et al., 2022). The TOPFARM framework has been used with SLSQP in several wind farm optimization studies (Riva et al., 2020; Ciavarra et al., 2022; Criado Risco et al., 2023; Rodrigues et al., 2023). In future work, this approach can be extended to co-optimize layout and control strategy – the SGD framework can naturally incorporate uncertainty quantification when modeling the potential control strategies for potential layouts (similar to the work in Gebraad et al., 2017; Quick et al., 2020; Howland et al., 2022).

While there are some wind plant optimization studies that resemble our approach, we are not aware of any studies that have applied SGD to the wind farm optimization problem (although SGD has been applied to other problems in engineering optimization, e.g., De et al., 2020; Sivanantham and Gopalakrishnan, 2022). Several wind farm optimization studies have made use of gradient-based optimization techniques (Herbert-Acero et al., 2014; Guirguis et al., 2016; Graf et al., 2016; Gebraad et al., 2017; Baker et al., 2019; Riva et al., 2020; Stanley et al., 2021; Croonenbroeck and Hennecke, 2021). Feng and Shen (2015) present a random search approach, moving the wind turbines one by one using a greedy algorithm. Some studies have employed neural networks to forecast power production (Godinho and Castro, 2021), estimate local atmospheric conditions (Stengel et al., 2020), suggest control strategies (Najd et al., 2020), or optimize engineering wake models (Zhang et al., 2021; Zhang and Zhao, 2022; Hussain et al., 2022), which all use SGD algorithms to train the parameters of the neural networks.

The remainder of the paper is the following. Section 2 outlines the SGD and deterministic optimization approaches used in this study. Section 3 details the wind farm optimization application cases examined. Section 4 discusses the results of these optimization comparisons. Section 5 provides conclusions and future research directions.

2 Methods

When deciding where to put wind turbines, a typical strategy is to maximize wind farm AEP while ensuring turbines are within the prospective site and are not spaced too closely together. In this study, we examine square and circular wind farms, where the corresponding optimization problems are posed as

$$\begin{array}{}\text{(1)}& \begin{array}{ll}\underset{\mathit{x},\mathit{y}}{\mathrm{maximize}}& \mathrm{AEP}(\mathit{x},\mathit{y})\\ \mathrm{subject}\mathrm{to}& {\left(x_i-x_j\right)}^{\mathrm{2}}+{\left(y_i-y_j\right)}^{\mathrm{2}}\ge {\left(N_{\mathrm{D}}D\right)}^{\mathrm{2}},\\ & \forall i\ne j\\ & x_l\le x_i\le x_u\\ & y_l\le y_i\le y_u\end{array}\end{array}$$

and

$$\begin{array}{}\text{(2)}& \begin{array}{ll}\underset{\mathit{x},\mathit{y}}{\mathrm{maximize}}& \mathrm{AEP}(\mathit{x},\mathit{y})\\ \mathrm{subject}\mathrm{to}& {\left(x_i-x_j\right)}^{\mathrm{2}}+{\left(y_i-y_j\right)}^{\mathrm{2}}\ge {\left(N_{\mathrm{D}}D\right)}^{\mathrm{2}},\\ & \forall i\ne j\\ & \sqrt{x_i^{\mathrm{2}}+y_i^{\mathrm{2}}}\le R,\end{array}\end{array}$$

respectively, where x and y are the turbine horizontal and vertical locations, x_l and x_u are the lower and upper horizontal square wind farm boundaries, y_l and y_u are the lower and upper square wind farm vertical boundaries, R is the radius of the circular wind farm, D is the rotor diameter, and N_D is the minimum allowable spacing between turbines measured in rotor diameters. From this point forward, we will use a single variable to represent the x and y locations, $\mathit{s}=\mathit{\{}\mathit{x},\mathit{y}\mathit{\}}$. When optimizing square wind farms, x_l, x_u, y_l, and y_u are constant.

The AEP is defined as

$$\begin{array}{}\text{(3)}& \mathrm{AEP}\left(\mathit{s}\right)=\mathrm{8760}\underset{\mathrm{0}}{\overset{\mathrm{2}\mathit{\pi }}{\int }}\underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}P\left(\mathit{s},u_{\mathrm{\infty }},\mathit{\theta }\right)\mathit{\pi }\left(u_{\mathrm{\infty }},\mathit{\theta }\right)\mathrm{d}{u}_{\mathrm{\infty }}\mathrm{d}\mathit{\theta },\end{array}$$

where P is power, π is probability, u_∞ is the freestream velocity, and θ is the freestream direction. The 8760 factor reflects the number of hours per year, converting from units of power to units of energy.

The AEP is typically estimated through rectangular quadrature, where the freestream velocity and direction are discretized using evenly spaced intervals,

$$\begin{array}{}\text{(4)}& \mathrm{AEP}\left(\mathit{s}\right)\approx \mathrm{8760}\sum _{d=\mathrm{1}}^D\sum _{u=\mathrm{1}}^{U}P\left(\mathit{s},{\mathcal{U}}_{\mathrm{u}},{\mathit{\theta }}_{\mathrm{d}}\right)\mathit{\rho }\left({\mathcal{U}}_{\mathrm{u}},{\mathit{\theta }}_{\mathrm{d}}\right),\end{array}$$

where 𝓤 is a vector of evenly spaced wind speeds, θ is a vector of evenly spaced wind directions, and ρ(𝒰_u,θ_d) is a probability mass function.

The AEP can also be estimated through Monte Carlo integration,

$$\begin{array}{}\text{(5)}& \mathrm{AEP}\left(\mathit{s}\right)\approx \mathrm{8760}{\displaystyle \frac{\mathrm{1}}{K}}\sum _{k=\mathrm{1}}^{K}P\left(\mathit{s},u_{\mathrm{\infty }}^{\left(k\right)},{\mathit{\theta }}^{\left(k\right)}\right),\end{array}$$

where $u_{\mathrm{\infty }}^{\left(k\right)}$ and θ^(k) represent draw k of the probability distribution π(u_∞, θ).

The associated AEP gradient can also be approximated through Monte Carlo simulation:

$$\begin{array}{}\text{(6)}& {\displaystyle \frac{\mathrm{d}}{\mathrm{d}\mathit{s}}}\mathrm{AEP}\approx \mathrm{8760}{\displaystyle \frac{\mathrm{1}}{K}}\sum _{k=\mathrm{1}}^K{\displaystyle \frac{\mathrm{d}}{\mathrm{d}\mathit{s}}}P\left(\mathit{s},u_{\mathrm{\infty }}^{\left(k\right)},{\mathit{\theta }}^{\left(k\right)}\right).\end{array}$$

2.1 Stochastic gradient descent

SGD is built upon the steepest descent algorithm. Early SGD algorithms added a moving average term (sometimes referred to as “momentum”) to avoid spurious oscillations (Tian et al., 2023). The conventional Adam SGD algorithm uses two moving averages: one of the gradient and one of the squared gradient. The ratio of these moving averages is used to determine the search direction. SGD algorithms are often combined with a learning rate scheduler, where the step size of the gradient descent is gradually decreased, allowing the optimization algorithm to hone in on the best solution. While the Adam algorithm is already designed to dynamically change the step size, including a learning rate scheduler can further improve the performance. The conventional Adam algorithm is designed for unconstrained optimization algorithms. In the following, we extend the algorithm to allow for deterministic constraints, which is a case that is common in mechanical engineering and unusual in the context of training neural networks. The basic idea is to aggregate the constraints to a penalty term with units that are consistent with the objective. The penalty term is designed so that, initially, the penalty gradients are of similar magnitude to the AEP gradients and so that the penalty gradients overwhelm the AEP gradients as the optimization continues. The SGD algorithm is shown in Algorithm 1, where s₀ is the initial turbine position; i is the iteration number; α_i is referred to as the constraint multiplier; γ(s) is a penalty function; $P(\mathit{s},u_{\mathrm{\infty }}^{\left(k\right)},{\mathit{\theta }}^{\left(k\right)})$ is the wind farm power associated with the inflow speed and direction, $u_{\mathrm{\infty }}^{\left(k\right)}$ and θ^(k); K is the number of samples employed in each SGD iteration; β₁ and β₂ are constants; T is the number of SGD iterations; 𝒮 is the learning rate scheduler; and η_i is the learning rate. By default, the early stopping option is false.

Algorithm 1 TOPFARM stochastic gradient descent implementation.

m←0, v←0, s←s₀

for i in [0, 1, 2, …, T−1]

do

if early_stopping and ${\mathit{\eta }}_i/{\mathit{\eta }}_{\mathrm{0}}\le$ threshold:

$\mathit{j}={\mathit{\alpha }}_i\frac{\partial \mathit{\gamma }}{\partial \mathit{s}}$

if $\left|\mathit{j}\right|\equiv \mathrm{0}$:

break

else:

$\mathit{j}=-\frac{\mathrm{8760}}{K}\sum _{k=\mathrm{1}}^K\frac{\partial }{\partial \mathit{s}}P\left(\mathit{s},u_{\mathrm{\infty }}^{\left(k\right)},{\mathit{\theta }}^{\left(k\right)}\right)+{\mathit{\alpha }}_i\frac{\partial \mathit{\gamma }}{\partial \mathit{s}}$

$\mathit{m}={\mathit{\beta }}_{\mathrm{1}}\mathit{m}-\left(\mathrm{1}-{\mathit{\beta }}_{\mathrm{1}}\right)\mathit{j}$

$\mathit{v}={\mathit{\beta }}_{\mathrm{2}}\mathit{v}-\left(\mathrm{1}-{\mathit{\beta }}_{\mathrm{2}}\right){\mathit{j}}^{\mathrm{2}}$

$\widehat{\mathit{m}}=\frac{\mathit{m}}{\mathrm{1}-{\left({\mathit{\beta }}_{\mathrm{1}}\right)}^i}$

$\widehat{\mathit{v}}=\frac{\mathit{v}}{\mathrm{1}-{\left({\mathit{\beta }}_{\mathrm{2}}\right)}^i}$

$\mathit{s}=\mathit{s}-{\mathit{\eta }}_i\widehat{\mathit{m}}/\sqrt{\widehat{\mathit{v}}}$

${\mathit{\eta }}_i=\mathcal{S}\left({\mathit{\eta }}_{\mathrm{0}},\mathit{\delta },i\right)$

${\mathit{\alpha }}_i={\mathit{\alpha }}_{\mathrm{0}}\frac{{\mathit{\eta }}_{\mathrm{0}}}{{\mathit{\eta }}_i}$

The spacing between turbines is enforced using a penalty term,

$$\begin{array}{}\text{(7)}& {\mathit{\gamma }}_{\mathrm{s}}=\sum _{\forall i,j>i}\min \left[{\left(x_i-x_j\right)}^{\mathrm{2}}+{\left(y_i-y_j\right)}^{\mathrm{2}}-{\left(N_{\mathrm{D}}D\right)}^{\mathrm{2}},\mathrm{0}\right].\end{array}$$

Similarly, the distance outside of boundaries is enforced using a penalty term. When considering square wind farms this penalty term is defined as

$$\begin{array}{ll}{\displaystyle }{\mathit{\gamma }}_{\mathrm{b}}& {\displaystyle }=\sum _{i=\mathrm{1}}^{N_{\mathrm{t}}}\left[\max {\left(x_i-x_{\mathrm{ub}},\mathrm{0}\right)}^{\mathrm{2}}+\max {\left(x_{\mathrm{lb}}-x_i,\mathrm{0}\right)}^{\mathrm{2}}\right.\\ {\displaystyle }& {\displaystyle }+\left.\max {\left(y_i-y_{\mathrm{ub}},\mathrm{0}\right)}^{\mathrm{2}}+\max {\left(y_{\mathrm{lb}}-y_i,\mathrm{0}\right)}^{\mathrm{2}}\right],\end{array}$$

and, in the case of circular boundaries, it is defined as

$$\begin{array}{}\text{(8)}& {\mathit{\gamma }}_{\mathrm{b}}=\sum _{i=\mathrm{1}}^{N_{\mathrm{t}}}\max {\left(\sqrt{x_i^{\mathrm{2}}+y_i^{\mathrm{2}}}-R,\mathrm{0}\right)}^{\mathrm{2}},\end{array}$$

where N_t is the number of wind turbines.

The total penalty, γ, is defined as the sum of these two penalty terms,

$$\begin{array}{}\text{(9)}& \mathit{\gamma }\left(\mathit{s}\right)={\mathit{\gamma }}_{\mathrm{s}}\left(\mathit{s}\right)+{\mathit{\gamma }}_{\mathrm{b}}\left(\mathit{s}\right).\end{array}$$

The gradient of the penalty term, γ, is scaled before being added to the negative gradient of the AEP using the scaling factor, α_i.

In Algorithm 1, the learning rate (η_i), constraint multiplier (α_i), number of SGD iterations (T), and the samples per SGD iteration (K) are all free parameters. These parameters can be optimized to perform well for individual wind farm optimization problems. But there is no guarantee that these particular parameters will perform well for other wind farm problems – and this meta-optimization can be expensive. In the machine learning community, these parameters are sometimes optimized using evolutionary, grid search, or Bayesian optimization approaches (Alibrahim and Ludwig, 2021). In addition, it is common to schedule the learning rate to decay as the optimization proceeds (You et al., 2019; Denkowski and Neubig, 2017).

We propose a method for setting free parameters to ensure that all units are consistent. The only free parameters we manually set are the number of optimization iterations and the number of power samples per iteration. The optimization generally becomes more accurate and more expensive as these parameters increase, and users are free to balance this trade-off as they see fit. Our formulation does not guarantee that all intermediate solutions satisfy the constraints, especially in the beginning of the optimization. The constraint multiplier begins on a comparable scale to the AEP and is scheduled to increase so that the constraint gradients overwhelm the AEP gradients as the optimization progresses. The number of iterations, T, can be based on a prescribed computational budget.

We initially attempted this approach using the widely used default parameter values in the original Adam algorithm, β₁=0.9 and β₂=0.999. The parameters can be thought of as adding momentum to the moving averages of the gradient and squared gradient, m and v. We found that these default values gave too much emphasis to gradients from the penalty function, launching the turbines away from the boundaries in a dramatic fashion. Instead, we suggest the parameters β₁=0.1 and β₂=0.2, which encode a shorter memory of the presence of the penalty. With these new default parameters, and the learning rate defined below, we observed successful convergence for a wide variety of test cases.

The learning rate, η_i, can be interpreted as converting $\widehat{\mathit{m}}/\sqrt{\widehat{\mathit{v}}}$ (with unity units) to distance (units of m). In this study, the learning rate is scheduled to decay according to

$$\begin{array}{}\text{(10)}& \begin{array}{rl}& S(t=\mathrm{0})={\mathit{\eta }}_{\mathrm{0}}\\ & S(t=T-\mathrm{1})={\mathit{\eta }}_T,\end{array}\end{array}$$

where T is the number of optimization iterations, η₀ is the initial learning rate, and η_T is the scheduled final learning rate. This final learning rate can be thought of as a solution tolerance for the design variables. In this study, we set η_T=0.1 m.

The initial learning rate, η₀, is based on a length scale parameter, L, which corresponds to a reasonable initial step size for the optimization. By setting the initial learning rate according to

$$\begin{array}{}\text{(11)}& {\mathit{\eta }}_{\mathrm{0}}=L=D/\mathrm{5},\end{array}$$

where D is the turbine rotor diameter, we encourage the turbines to move at most L distance every optimization iteration.

The learning rate is scheduled to decay as

$$\begin{array}{}\text{(12)}& S\left({\mathit{\eta }}_{\mathrm{0}},\mathit{\delta },t\right)={\mathit{\eta }}_{\mathrm{0}}\prod _{i=\mathrm{0}}^t{\displaystyle \frac{\mathrm{1}}{\mathrm{1}+i\mathit{\delta }}},\end{array}$$

where δ is a parameter that controls the learning rate length, such that the final learning rate is η_T. The parameter δ is numerically set as

$$\begin{array}{}\text{(13)}& \mathit{\delta }\left({\mathit{\eta }}_{\mathrm{0}},{\mathit{\eta }}_T,T\right)=\underset{\mathit{\delta }}{\mathrm{argmin}}|{\mathit{\eta }}_T-S\left({\mathit{\eta }}_{\mathrm{0}},\mathit{\delta },T\right)|.\end{array}$$

The constraint multiplier, α_i, can be interpreted as converting the gradient of constrained square distances (in units of m) to AEP gradients. The initial constraint multiplier, α₀, is set as the mean absolute AEP gradient divided by the length scale, L, so that the separation constraint has a similar scale to AEP gradients,

$$\begin{array}{}\text{(14)}& {\mathit{\alpha }}_{\mathrm{0}}={\displaystyle \frac{\mathrm{mean}\left[\left|\mathrm{\nabla }\mathrm{AEP}\left({\mathit{s}}_{\mathrm{0}}\right)\right|\right]}{L}},\end{array}$$

where mean[|∇AEP(s₀)] is the mean of the absolute AEP gradient of the initial guess with respect to each component of the gradient. During each iteration, the constraint multiplier, α_i, is scheduled to increase based on the inverse of the learning rate,

$$\begin{array}{}\text{(15)}& {\mathit{\alpha }}_i={\mathit{\alpha }}_{\mathrm{0}}{\displaystyle \frac{{\mathit{\eta }}_{\mathrm{0}}}{{\mathit{\eta }}_i}}.\end{array}$$

The wind rose samples, $(u_{\mathrm{\infty }}^{\left(i\right)},{\mathit{\theta }}^{\left(i\right)})\sim \mathit{\pi }(u_{\mathrm{\infty }},\mathit{\theta })$, are randomly selected based on the direction frequency and direction-specific Weibull shape and scale parameters. Note that the tilde (∼) denotes a shared probability distribution. After a direction is sampled, the wind speed is sampled as a continuous Weibull-distributed random variable,

$$\begin{array}{}\text{(16)}& u_{\mathrm{\infty }}\left(\mathit{\theta }\right)\sim \mathcal{W}\left[u_{\mathrm{\infty }},A\left(\mathit{\theta }\right),k\left(\mathit{\theta }\right)\right],\end{array}$$

where the probability density of the Weibull distribution, 𝒲, is given by

$$\begin{array}{}\text{(17)}& \mathcal{W}\left(u_{\mathrm{\infty }},A,k\right)={\displaystyle \frac{k}{A}}{\left({\displaystyle \frac{u_{\mathrm{\infty }}}{A}}\right)}^{k-\mathrm{1}}\exp \left[-{\left({\displaystyle \frac{u_{\mathrm{\infty }}}{A}}\right)}^k\right].\end{array}$$

2.2 Deterministic approach

The SLSQP algorithm (Kraft, 1988) is selected to be the deterministic optimization algorithm to act as a benchmark to the SGD approach. SLSQP is a conventional deterministic optimization approach. It is employed in many open-source engineering design codes (Allen et al., 2020; Wu et al., 2020; Ciavarra et al., 2022) and has been used in previous comparisons of optimization algorithms (Lam et al., 2018; Li and Zhang, 2021; Fleming et al., 2022).

The spacing and boundary constraints are passed to the optimizer as individual inequality constraints. The spacing constraints are defined as

$$\begin{array}{}\text{(18)}& C_{ij}={\left(x_i-x_j\right)}^{\mathrm{2}}+{\left(y_i-y_j\right)}^{\mathrm{2}}-{\left(N_{\mathrm{D}}D\right)}^{\mathrm{2}}\forall i,j>i,\end{array}$$

where C is an upper triangular matrix of nonlinear inequality constraints.

Square wind farm boundaries are represented using four inequality constraints per turbine,

$$\begin{array}{}\text{(19)}& \begin{array}{rl}& D_{i\mathrm{1}}=x_i-x_u\\ & D_{i\mathrm{2}}=x_l-x_i\\ & D_{i\mathrm{3}}=y_i-y_u\\ & D_{i\mathrm{4}}=y_l-y_i,\end{array}\end{array}$$

and the circular wind farm boundaries are represented with one inequality constraint per turbine,

$$\begin{array}{}\text{(20)}& D_i=\sqrt{x_i^{\mathrm{2}}+y_i^{\mathrm{2}}}-R,\end{array}$$

where D is a matrix of boundary constraints that must be less than or equal to 0.

3 Application

We apply the optimization approaches discussed above to optimize wind power plants of various sizes using the TOPFARM framework (DTU Wind Energy Systems, 2023b). Each farm consists of turbines with 70 m hub heights, 80 m rotor diameters, and 2 MW rated powers. Power is computed using PyWake (Pedersen et al., 2019), which is an open-source wake modeling tool that has been used in several related studies (Riva et al., 2020; Rodrigues et al., 2022; Ciavarra et al., 2022; van der Laan et al., 2023; Fischereit et al., 2022). Power gradients are computed directly from PyWake using automatic differentiation. The power of each turbine is estimated by a combination of velocity deficits predicted by the Bastankhah Gaussian wake model (Bastankhah and Porté-Agel, 2014) using the default parameters in the PyWake tool and the squared sum superposition (Pedersen et al., 2019; DTU Wind Energy Systems, 2023a). We require each turbine to be spaced at minimum two rotor diameters apart (N_D=2). This is imposed as an optimization constraint. We considered wind farms with square and circular boundaries. The square wind farm boundaries are determined as

$$\begin{array}{}\text{(21)}& \begin{array}{rl}& x_l=\mathrm{0}\\ & y_l=\mathrm{0}\\ & x_u=D\left(\sqrt{N_{\mathrm{t}}}-\mathrm{1}\right)\mathrm{\Delta }\\ & y_u=D\left(\sqrt{N_{\mathrm{t}}}-\mathrm{1}\right)\mathrm{\Delta },\end{array}\end{array}$$

and, in cases with circular wind farm boundaries, the radius is determined as

$$\begin{array}{}\text{(22)}& R=D\left(\sqrt{N_{\mathrm{t}}}-\mathrm{1}\right)\mathrm{\Delta },\end{array}$$

where the Δ parameter controls the average spacing of the turbines. In this study, Δ=5.

We use the pyOptSparseDriver (Wu et al., 2020) SLSQP (Kraft, 1988) implementation (Virtanen et al., 2020) in TOPFARM. The optimizer was set to run for 300 maximum iterations with a tolerance of 10⁻¹. The TOPFARM “expected_cost” parameter is set to 10. The turbine coordinates are normalized from 0 to 1. In each optimization iteration, the AEP, and the corresponding gradient, is computed using rectangular quadrature as described in Eq. (4), using 360 wind direction bins and 23 wind speed bins, resulting in 8280 power evaluations.

The wind rose, visualized in Fig. 1, is based on PyWake's Lillgrund example site. A probability mass function is assigned to different direction bins. Each direction bin is associated with Weibull scale and shape parameters describing the distributions of wind speeds within the sector. This probability mass is derived from 7 months of measured data used in a previous study (Göçmen and Giebel, 2016). Each direction bin is 30^∘ wide. The reasoning behind this is similar to that behind the IEC 614 400 power curve standard (International Electrotechnical Commission, 2005) – it is crucial that the reference data consider a statistically significant number of data in each bin. This coarse direction discretization results in a faster convergence of the estimated probability mass function than a finer discretization would. The continuous probability density function π(u_∞,θ) is approximated as ρ(θ)π(u_∞|θ), where ρ is the previously mentioned probability mass function, linearly interpolated across 1^∘ bins, and π(u_∞|θ) is parameterized by direction-specific Weibull shape and scale parameters that are also linearly interpolated from the provided data. With this formulation, the likelihood of different wind directions is provided as a probability mass function, ρ(θ). This probability mass is used as weights passed to the Numpy “choice” function (Harris et al., 2020), allowing the wind direction to be sampled as a discrete random variable. We note that this formulation could be extended to a fully continuous formulation by drawing the direction samples from the inverse of an empirical cumulative direction density function.

Figure 1 Lillgrund wind speed and direction probability mass function with 360 direction bins and five wind speed bins, where the probability mass function is a linear interpolation of coarser measurements.

In all wind farm optimization problems considered, constraint gradients, and the associated penalty function gradients, are computed analytically. The AEP gradient is computed via automatic differentiation. The directions are discretized from 0 to 360^∘, with 1^∘ increments. In the deterministic formulation, the discretized wind speed ranges from 3–25 m s⁻¹ and is divided using increments of 1 m s⁻¹.

While each Monte Carlo estimate of AEP has significant error, the average error will be close to 0 throughout the course of the SGD optimization. We compare the accuracy of the Monte Carlo approach (Eq. 5) and the quadrature approach (Eq. 4) to estimate the AEP and the L-2 norm of the AEP gradient. The true values are estimated with a very fine discretization of speed and direction, 0.2 m s⁻¹ and 0.2^∘. These are used as reference values to assess the accuracy of the Monte Carlo and quadrature approaches by comparing the errors associated with both approaches as functions of the number of samples and the discretization level, respectively, when analyzing a 100-turbine farm with square boundaries. This convergence analysis is shown in Fig. 2. While some realizations of the Monte Carlo approach yield more accurate results than the quadrature approach, the quadrature approach is generally more accurate than Monte Carlo sampling. The Monte Carlo approach requires on average around 10 times as many power evaluations to obtain the same accuracy as the deterministic approach.

Figure 2 Convergence of AEP (a) and L-2 norm of the AEP gradient (b) with respect to the number of samples used in the quadrature and Monte Carlo techniques. The grey cloud shows the 90 % confidence interval associated with the Monte Carlo approach, using 50 samples. The dashed black line shows the average error associated with the Monte Carlo approach. The solid black line shows the error associated with the deterministic approach.

In this study, we select 50 samples for every SGD iteration. Figure 3 shows the measured computational cost of computing AEP gradients using the circular wind farm described in this study, with different wind farm sizes. The minimum measured time is reported as the minimum of 30 identical runs on the DTU Sophia supercomputer (Technical University of Denmark, 2019). The computational time generally scales logarithmically with the number of turbines. This is to be expected, as there are more interaction terms in the wake model as more turbines are considered. The computational time does not scale logarithmically with the number of wind rose samples. For small numbers of turbines, evaluating 10 wind rose samples is about as expensive as evaluating 50 samples. This scaling changes as the wind farm grows in size, and it gradually becomes more expensive to sample the wind rose. The evaluation time appears to converge to a logarithmic scaling for large numbers of wind rose samples. These scaling results are likely influenced by memory limitations.

Figure 3 Computational time associated with computing the gradient for various wind farm sizes and wind rose sampling strategies. The left panel compares the cost of computing AEP gradients when using different numbers of Monte Carlo samples with the cost of the full factorial wind rose (8280 samples) for farms with various numbers of wind turbines. The right panel compares the cost of Monte Carlo estimates of the AEP gradient for different numbers of samples of the atmospheric conditions and turbines in the wind farm.

The optimization algorithms are timed based on the time elapsed between the first and final optimization gradient evaluations. Each optimization case is run on first-generation AMD EPYC 7351 processors.

4 Results and discussion

In the following subsections, the performance of SGD and SLSQP is compared for wind farms with square and circular boundaries, and the sensitivity of the SGD algorithm is assessed.

4.1 Square wind farm

The performance of SGD is compared to the deterministic counterpart, considering wind farms with 100, 144, 225, and 324 turbines, with square boundary constraints, using 20 different initial starting conditions to obtain statistically significant results. The AEP, constraint violation (γ), and time elapsed associated with each optimization solution are plotted in Fig. 4. The SGD approach consistently yields higher AEPs than the SLSQP approach when the number of scheduled SGD iterations, T, is 2000. There is a large range of computational times associated with the SLSQP approach, though the computational expense of SLSQP generally grows much larger than SGD as the number of turbines is increased. This is largely due to the nature of the turbine spacing constraint, the size of which grows as the number of turbines squared. SLSQP takes about as much computational time as the SGD approach with 500 scheduled iterations when there are 100 turbines in the square farm. As the number of turbines grows, the average time required by SLSQP becomes more costly than SGD with 2000 scheduled iterations. The computational cost of SLSQP is a strong function of the initial layout, and the variance of the SLSQP optimization time also increases with the wind farm size. This is due to the complex interaction between the linear boundary constraints and nonlinear spacing constraints. As the proposed SGD formulation does not offer an automatic way to set the number of SGD iterations, T, results are shown for different values of T. When T is increased, the optimizer finds solutions with larger AEPs, with a computational cost that is approximately proportional to T. The SGD solution consistently improves as more optimization iterations are scheduled (larger values of T). Results associated with 1000 SGD iterations tend to yield similar AEPs to the SLSQP approach, and results with 2000 SGD iterations tend to yield higher AEPs than the SLSQP designs.

Figure 4 Optimization results associated with SGD and SLSQP for square wind farms with 100, 144, 225, and 324 turbines, using 20 random initial starting conditions. The AEP (top panels), constraint penalty (middle panels), and computational time (bottom panels) are plotted as boxplots. The SGD results are plotted for T=500, 1000, and 2000 iterations.

The final layouts associated with one of the random initial conditions used in the 324-turbine analysis, when T=2000 iterations, are shown in Fig. 5. The SGD approach generally identifies solutions with the majority of turbines packed into the side boundaries. The deterministic algorithm also packed turbines into the edges of the farm, although not as many turbines were packed into the east and west boundaries as in the SGD results. The layouts found using the SGD approach tend to have interior turbines that generally appear to be more aligned in the north–south direction than in the deterministic solutions.

Figure 5 The final layouts found using SLSQP (a) and SGD (b) using one of the random initial layouts examined in a 324-turbine wind farm with square boundaries for T=2000 iterations. Each circle has a radius of one rotor diameter.

The results of the 100-, 144-, 225-, and 324-turbine wind farm optimization cases are summarized in Table 1. The mean time, mean constraint violation, and mean and standard deviation of the AEP are reported with respect to the 20 random initial starting conditions. Constraint violation is reported as $\sqrt{\mathit{\gamma }/N_{\mathrm{t}}}$ to quantify the mean length of the constraint violations of each turbine. The final constraint violations can be reduced by lowering the η_T parameter. In all of these cases, the SLSQP optimization resulted in solutions with zero constraint violations. This is likely because of the linear formulation of the boundary constraints – when a solution satisfies the spacing constraint, any solutions that satisfy the boundary constraints can quickly be found. SGD with 2000 iterations generally yields solutions with AEP that are 0.3 %–0.5 % higher than the solutions found using SLSQP. This is likely because the SGD algorithm is able to better explore the design space by initially relaxing the constraints, allowing for some initial constraint violations.

Table 1 Results of SGD and deterministic optimizations for various square wind farm sizes. Each optimization case is run using 20 random initial starting conditions, and the mean and standard deviation are reported with respect to these 20 initial points. The SGD results are associated with T=2000 iterations.

4.2 Circular wind farms

To ensure that the previously presented results are not specific to square wind farms, we performed a similar set of analyses examining circular wind farms. The results yield similar trends to the analysis using square wind farms – SGD becomes significantly less time-consuming than SLSQP as the number of turbines increases and generally yields solutions with slightly larger AEPs.

Figure 6 Optimization results associated with SGD and SLSQP for circular wind farms with 100, 144, and 225 turbines, using 20 random initial starting conditions. The AEP (top panels), constraint penalty (middle panels), and computational time (bottom panels) are plotted as box and whisker plots. The SGD results are plotted for T=500, 1000, and 2000 iterations.

Circular winds farms were optimized using 20 random initial layouts, examining different farm sizes, using the SGD and SLSQP optimization algorithms. The results are summarized in Fig. 6. The circular wind farm optimization generally took longer than the square wind farm when using the SLSQP optimizer. This is likely due to the more complicated nature of the circular boundary when using Cartesian coordinates. These results are similar to the results in Sect. 4.1 – as the number of wind turbines and scheduled SGD iterations increases, SGD tends to find solutions with larger AEPs in less computational time.

The results of the circular wind farm optimization are compared between the SGD and SLSQP optimizers in Table 2, where SGD is scheduled to run for 2000 optimization iterations. SGD generally results in about 0.5 % more AEPs in significantly less time than SLSQP as the number of turbines is increased. The result area is also compared in Fig. 7. The SGD optimizer generally results in more turbines on the boundary edge than the SLSQP optimizer.

Table 2 Results of SGD and deterministic optimizations for various circular wind farm sizes. Each optimization case is run using 20 random initial starting conditions, and the mean and standard deviation are reported with respect to these 20 initial points. The SGD results are associated with T=2000 iterations.

Figure 7 The final layouts found using SLSQP (a) and SGD (b) using one of the random initial layouts examined in the 225-turbine wind farm with circular boundaries using T=2000 iterations. The final turbine layouts are shown as filled circles. Each circle has a radius of one rotor diameter.

4.3 Sensitivity analysis

There are several parameters in the SGD algorithm that were tuned to perform reasonably well. In this section, we investigate the sensitivity of the SGD optimization results with respect to the early stopping option in Algorithm 1, the number of Monte Carlo samples per optimization iteration, the learning rate schedule, and the initial and final learning rates.

As the optimization progresses, the constraint multiplier, α_i, becomes large (approaching 10 as the learning rate approaches 0.1), and the gradients of the AEP are overwhelmed by the gradients of the penalty, which take very little time to compute. This situation can be addressed by using the early stopping option in Algorithm 1. The solution tends to terminate quickly when the optimizer only follows the deterministic gradient (the optimization engine terminates when the constraint gradients are 0). Figure 8 shows the AEP, constraint violation, and computational time associated with five random initial layouts, using threshold parameters of 0.01, 0.05, and 0.1, as well as the SGD algorithm as applied in the previous sections, without the early stopping option activated. The use of each early stopping option results in layouts without constraint violations. As the threshold parameter is increased, the AEP is slightly reduced, and the total computational time decreases. A threshold parameter of 0.1 results in approximately 0.3 % reduction in AEP and 44 % reduction in computation time.

Figure 8 Results using different configurations of the early stopping option in TOPFARM, as well as the SGD optimization without early stopping, considering the circular wind farm with 225 turbines, with T=2000, using five random initial layouts.

The optimization results presented in this study used 50 power samples per iteration (K=50). We found this to produce high-quality results without incurring unacceptable computational expense. Figure 9 shows the behavior of the SGD approach associated with different values of K, considering 100 turbines with 2000 scheduled optimization iterations. As K increases, the optimization finds solutions with larger AEPs. There is a small increase in time elapsed and a large increase in the final AEP between the K=5 and K=50 cases, while there is a large increase in time elapsed and a small increase in the final AEP between K=50 and K=200. As K increases, we expect the maximum AEP to reach a plateau and the time and memory required to increase indefinitely. In future work, we plan to explore scheduling K to change as the optimization progresses.

Figure 9 Optimization results associated with SGD for a square 100-turbine wind farm using 20 random initial starting conditions and T=2000. The upper and lower bounds of the results are plotted as a function of the optimization iteration number. The SGD results associated with K=5, 50, and 200 iterations are shown in purple, blue, and yellow, respectively.

This study used an exotic learning rate scheduler. We tried several schedulers and observed this one to be the best at finding sufficiently large AEP solutions that reasonably satisfied the imposed constraints. Figure 10 shows the behavior of the SGD algorithm associated with the presented learning rate scheduler, referred to here as the product scheduler, as well as an exponential and a linear decay scheduler. The exponential scheduler quickly diminishes the learning rate, causing the SGD algorithm to become stuck in local minima. The linear transition from large to lower learning rates prevents the SGD algorithm from having sufficient time to follow enlarged constraint gradients. It is possible that the algorithm could be improved by using separate schedulers for the learning rate and constraint multiplier. For instance, it might be more effective to use a linear scheduler to decrease the learning rate and an exponential scheduler to increase the constraint multiplier. We leave this question for future work.

Figure 10 Influence of the learning rate scheduler on the SGD optimization, considering a 100-turbine wind farm with a square boundary. The product scheduler is shown in yellow. The exponential scheduler is shown in blue. The linear scheduler is shown in purple. The AEP (a), constraint penalty function (b), and the learning rate decay (c) are plotted as a function of the number of optimization iterations. The optimization iteration is denoted as t in the legend of panel (c).

The initial and final learning rates, η₀ and η_T, have units of distance and correspond to the initial and final step size of the optimization algorithm. The final learning rate can be interpreted as the degree to which the constraints are to be satisfied, since this will be the step size the optimization algorithm uses when α_i is large and the constraint gradients overwhelm the AEP gradients. This is illustrated in the left panel of Fig. 11, which shows the results of several SGD optimizations using different initial and final learning rates. On average, there is a linear relationship between the constraint violation of the solution and η_T, where the average final constraint violation is approximately 2 times η_T. The maximum observed constraint violation is approximately 4 times η_T. In addition, there is a trade-off between the AEP and constraint violation of the final solution. This trade-off is influenced by the initial and final learning rates. It is important to tune the initial learning rate. An initial learning rate that is too low will result in very little exploration. Initial learning rates that are too high will result in a rapid influx of penalty violations that overwhelm AEP gradients throughout the optimization. From our experiments, we found a step size of one-fifth of the rotor diameter to produce satisfactory results, as shown in the right panel of Fig. 11.

Figure 11 (a) The final constraint penalty plotted against the final learning rate. (b) The final AEP plotted against the final constraint penalty. The different colors represent different initial learning rates. The results of 20 initial starting positions are plotted as points. The average results of the 20 initial starting conditions are connected as lines. These data are associated with 100-turbine wind farms with square boundaries and T=2000 iterations.

5 Conclusions

SGD is a promising optimization tool for wind farm design. Instead of evaluating all anticipated atmospheric conditions during every optimization iteration, SGD randomly samples the defined distributions of atmospheric conditions, resulting in substantially reduced computational time required for each optimization iteration. The total optimization time can be scheduled according to a prescribed computational budget. The presented formulation allows for continuous resolution of uncertain variables, eliminating the need to choose a discretization resolution of atmospheric conditions, such as the wind speed and direction. This technique does not become exponentially more expensive as a greater number of uncertain parameters is included, allowing for consideration of other atmospheric conditions, such as turbulence intensity, air density, veer, and shear (Saint-Drenan et al., 2020; Duc et al., 2019).

The presented SGD approach was shown to become more effective than a deterministic counterpart as the number of wind turbines increased. SGD yielded slightly higher AEPs than the deterministic approach in substantially reduced computational time. The time required to optimize wind farm layouts can be a major bottleneck in corporate workflows, and the time savings associated with the SGD approach allows engineers to access optimization results sooner than a conventional approach. If the inflow conditions were discretized using extremely small bins, or if several atmospheric conditions were to be considered, we expect that the SGD approach would perform the optimization even faster and more effectively than the deterministic approach.

The SGD approach is a simple framework that is well suited to large-scale stochastic wind power plant design optimization challenges. This framework is available in the open-source TOPFARM package. Future work includes exploring separate schedulers for the constraint multiplier and learning rate and scheduling the number of Monte Carlo samples, K, to change as the optimization proceeds.