A Neighborhood Search Integer Programming Approach for Wind Farm Layout Optimization

. Two models and a heuristic algorithm to address the wind farm layout optimization problem are presented. The models are linear integer programming formulations where candidate locations of wind turbines are described by binary variables. One formulation considers an approximation of the power curve by means of a step-wise constant function. The other model is based on a power-curve-free model where minimization of a measure closely related to total wind speed deficit is aimed. A special-purpose neighborhood search heuristic wraps the formulations in order to increase tractability and effective-5 ness compared to the full model. The heuristic iteratively searches neighborhoods around the incumbent using a branch-and-cut algorithm. The number of candidate locations and neighborhood sizes are adjusted adaptively. Numerical results on a set of publicly available benchmark problems indicate that a proxy for total velocity deficit as objective is a functional approach, since high-quality solutions of an annual energy production metric are found. Furthermore, the proposed heuristic is able to match and in some cases improve the results obtained when considering the turbine positions as continuous variables. 10


Introduction
Cost reductions for renewable energy generation is on the top of political agendas, with the objective of supporting the worldwide proliferation of clean energy production systems.Subsidy-free tendering processes become more frequent, as is the case for offshore wind auctions in Germany since 2017 and in Netherlands since 2018, or in China for onshore wind from 2021 (GWEC, 2020a).The fast evolution of offshore wind in the last decade, with an sharp growth of global installed capacity (GWEC, 2020b), is yet another clear indicator of the maturity of the industry.As wind energy is unleashing its potential to contribute for a successful green energy transition, the need to further economically refine wind farm designs turns into a crucial task.
The basic Wind Farm Layout Optimization (WFLO) problem consists in deciding the positioning of Wind Turbines (WTs) within a given project area to maximize the Annual Energy Production (AEP) while respecting a minimum separation distance between the generators.The classic problem definition is to place a fixed number n T of homogeneous (single type) WTs.This problem has been studied broadly and intensively since at least three decades (Herbert-Acero et al., 2014).The first effort in the topic was the pioneering work of Mosetti et al. (Mosetti et al., 1994), where the Katic-Jensen wake decay model (Katic  et al., 1986), implemented to compute wake losses, is coupled with a genetic algorithm (Deb, 2013) as optimizer to iteratively improve the layout.
In this perspective, the main components when building an optimization workflow for the WFLO problem are the wake and the optimization models, and the associated numerical algorithms.For the first, in order to be able to formulate tractable frameworks, the designer needs to rely on the so-called engineering models.These are essentially mathematical representations which can be expressed in terms of analytical equations after greatly simplifying complex physics modelling, while still capturing at a good extent the underlying nature of the phenomenon under analysis.Recently, scientific articles in this field have proposed and validated engineering wake models with smooth and differentiable velocity deficit shape, as the Bastankhah's Gaussian (Bastankhah and Porté-Agel, 2016) or its simplified version (Thomas and Ning, 2018), the Niayifar and Porté-Agel (Niayifar and Porté-Agel, 2015) or the Jensen cosine model (Thomas et al., 2022b).Likewise, the aggregation of individual wake velocity deficits can be done through linear superposition or root sum squares (Porté-Agel et al., 2020).
Optimization techniques for the WFLO problem can be classified, depending on the choice of variables, into continuous and discrete optimization.In the field of continuous optimization, the location p i of a WT i, in terms of the abscissa (x i ) and ordinate variables (y i ) in the Cartesian plane, p i = (x i , y i ), can take any real values, while ensuring that the point is within the project area F, and simultaneously satisfying the minimum distance constraints.Several gradient-free algorithms have been applied to this problem, including metaheuristics, as genetic algorithm (Réthoré et al., 2014) or particle swarm optimization (Wan et al., 2010).Likewise, gradient-based methods can be utilized for this problem, as for example the Sparse Nonlinear OPTimizer (SNOPT), that uses a Sequential Quadratic Programming (SQP) approach (Thomas et al., 2022a), or interior-point solvers (Pérez et al., 2013).In general, gradient-free algorithms, although highly flexible for modelling aspects, have considerably poorer scalability for larger problem sizes than gradient-based approaches.Re-parametrization approaches aiming to reduce the number of variables through simplified geometrical representations of the problem, such as row and column spacing or inclination angle, are also emerging within this context bringing along various enhancements (Stanley and Ning, 2019).Additionally, multi-start strategies are frequently implemented as a workaround for the intrinsic multi-modal nature of the WFLO problem.Finally, hybrid methods combining gradient-free and gradient-based algorithms are also an alternative able to come up with competitive results (Mittal and Mitra, 2017).
Discrete optimization models can be formulated for this problem by means of sampling the available project area in form of N candidate location points.Thus, only a set of finite options from the continuous search space are considered, where the n T WTs are installable, with in principle N ≫ n T .In contrast to continuous optimization, a candidate point i is then represented by a binary variable ξ i , that gets a value of one if a WT is installed at that location, or zero otherwise.The vast majority of articles in the literature implement gradient-free algorithms for this technique, as the works of Mosetti et al. (Mosetti et al., 1994) and Grady et al. (Grady et al., 2005), where both use genetic algorithms.This modelling technique, however, fits very well in the well-studied general framework of integer programming.The main advantage of this approach is the possibility to utilize exact solvers based on branch-and-cut method; theoretically, being able to solve a problem to optimality while supporting common engineering constraints (Wolsey, 2020).Nevertheless, the low tractability and poor scalability of this method in function of the size of N and the number of state variables is well-known.tional experiments are deployed in Sect.5, followed up by discussions in Sect.6, and lastly the manuscript is finalized with the conclusions in Sect.7.

Physics Modelling
The proposed MILP models and general optimization framework in this article can be easily applied to many wake deficit models.No particular properties on smoothness or differentiability, or specific demands on mathematical structure are required in connection with controlling wake diameter and deficit (Thomas et al., 2022b).Since the computational results in the article are obtained after solving open access case studies from the IEA37 Wind Task 37 (Baker et al., 2019;Dykes et al., 2015), the implemented wake model from that source is presented in Sect.2.1, along with the used superposition techniques in Sect.2.2, WT power curve in Sect.2.3, and the AEP calculation procedure in Sect.2.4.Variations on ways of computing the absolute velocity deficits and linear wakes superposition under the framework of MILP are introduced as well.

Wake Model
A simplified version of the Bastankhah's Gaussian is considered (Thomas and Ning, 2018).The relative velocity deficit behind a single WT located at ℓ, and evaluated at point i, is described using the model and notation from (IEA Wind Task 37, 2019). (1) where u ∞ is the inflow wind speed, C T is the thrust coefficient, d ∥ ij = xi − xℓ is the stream-wise distance from the hub generating wake (x ℓ ) to hub of interest (x i ) along freestream, d ⊥ ij = ȳi − ȳℓ is the span-wise distance from the hub generating wake to hub of interest perpendicular to freestream, σ y is the standard deviation of the wake deficit, k y is a variable based on a turbulence intensity, and D is the WT diameter.

Wake velocity deficit superposition
The absolute velocity deficit ∆u iℓ (θ j , k) at wind direction θ j and wind speed k, ∆u iℓ (θ j , k), can be estimated in two ways.
Either it is based on the inflow wind speed through or it is based on the wind speed u ℓjk at WT ℓ creating the wake at point i for wind direction θ j and speed k, here δu iℓ (θ j , k) is the relative velocity deficit of i over j at same operation condition after Eq.( 1) and Eq.( 2).Note that Eq. ( 3) leads to a greater value and therefore is considered a conservative approach compared to (the potentially more realistic) Eq. ( 4) (Niayifar and Porté-Agel, 2015).Nonetheless, implementing Eq. ( 3) greatly simplifies the resultant system of equations and allow for preprocessing calculations.
Let the set U θ j i collects the WTs creating wake over WT at point i for wind direction θ j and speed k as per The wake velocity deficit superposition, to calculate the total velocity deficit at WT i, ∆u i (θ j , k), can be obtained through two mechanisms.Either it is based on linear superposition model through or it is based on the root sum squares superposition model

WT Power Curve
For computing AEP suitable power curves are required for the available turbine types.These power curves are often not perfectly suitable for optimization, due to the non-differentiability at rated wind speed.The general characteristics are that the power curve is zero below cut-in wind speed, zero above the cut-out wind speed, constant between the rated wind speed and the cut-out wind speed.In between the cut-in and rated wind speeds the curve is assumed to be smooth, convex and monotonically increasing.The simplified power curve for a generic turbine as a function of wind speed u is modelled through where p rated is the nominal power at (and above) rated wind speed u rated .The other turbine characteristics are the cut-in wind speed u cut-in , and the cut-out wind speed u cut-out .

Annual Energy Production, AEP
The AEP is obtained after the application of the following expression https://doi.org/10.
where w jk is the joint probability of wind direction j and wind speed k, and 8760 is the number of hours of a standard year.

Optimization Models
The MILP model with explicit modelling of the WT power curve, wake model and wakes superposition is first introduced in Sect.3.1.Then, the power-curve-free formulation, based on total wind speed within the farm, is deployed in Sect.3.2.
For both models, the main type of variables ξ i ∈ {0, 1} represent presence or absence of turbines at the candidate locations as mentioned before.Given N points, i.e. candidate locations for turbine positions, with positions p i inside the domain F (i.e. p i ∈ F for all i WT candidate locations), binary variables ξ i ∈ {0, 1} are associated with the following interpretation , if a turbine is located at point i with position p i , and Let the index sets N i storing the candidate locations violating the minimum distance constraints for a WT i be defined as where d min > 0 is the minimum required distance between two turbines (2D in this study).If ξ i = 1 then all binary variables in the set N i should be forced to zero, whereas if ξ i = 0 these variables should be free to take any value in {0,1}.
All relevant distances can be preprocessed for all combinations of points i and j.These pertinent parameters are then defined as function of the Cartesian plane positions p and wind direction θ j , as the Euclidean distances the stream-wise distances d ∥ ij (p; θ j ) and the span-wise distances d ⊥ ij (p; θ j ), extending the concept introduced in Sect.2.1.

Power-curve-based model
For the purpose of wake modelling and power computation continuous state variables u ijk are used.A variable u ijk represents the wind speed at WT location i, for wind direction j, and wind speed k.
The general characteristics of the power curve outlined in Sect.2.3 is used to approximate the curve by a step-wise function.
The cubic domain of the power curve is first partitioned into m intervals, plus one interval from a negative point (−u ini ) to the cut-in speed, and a final one to cover the range from rated to cut-out speed.Each isometric interval of length ∆u is approximated with a constant power value, see Fig. 1.
An interval l is characterized by three parameters u l s , u l m , and u l h with the next properties 6 https://doi.org/10.5194/wes-2022-82Preprint.Discussion started: 20 September 2022 c Author(s) 2022.CC BY 4.0 License.Equation (12) defines the lower and upper limits for the extreme intervals l = 1 and l = m + 2, Eq. ( 13) formalizes them for a interval l in the cubic part of the curve (2 ≤ l ≤ m + 1), and Eq. ( 14) presents how to determine the extracted wind speed 170 associated to the power value in l, which is the average value of u l s and u l h .Let define the binary state variables η l ijk ∈ {0, 1} for l = 1, . . ., m + 2 with the interpretation subject to: The objective function in Eq. ( 16a) is an approximation of the AEP computation presented in Eq. ( 9).Equation (16b) models the minimum distance constraints as explained in the introduction of Sect. 3. If a binary variable ξ i is active, then all candidate points closer than d min should be excluded, i.e. set to zero.If a binary variable ξ i is inactive then the other candidates are still eligible.The definition of the set N i is provided in Eq. (11).Equation (16c) models the situation that the designer requires at least n min and at most n max WTs to be located in the domain.Note that for the classic problem definition n min = n max =n T .
Equation (16d) connects both state variables u and η as explained in Eq. (15).Eq. (16e) forces that there is one operation case active for each WT candidate at each wind direction and speed.The last constraint in Eq. ( 16f) is for the wake velocity deficit and wakes superposition modelling to calculate wind speed for each candidate location at each wind direction and inflow speed The presented model supports a conservative velocity deficit approach (Eq.( 3)) with linear superposition (Eq.( 6)).The definition of set U θ j i is provided in Eq. ( 5).Note that an extension, consisting in creating extra continuous state variables and associated constraints, could allow for considering the more realistic approach in Eq. ( 4).It is still unknown if the root sum squares model of Eq. ( 7) could be implemented in the framework of MILP.Finally, Eq. ( 16f) defines the domain of the required variables.A value for u ini of u cut-out is set up.

Power-curve-free model
Albeit the formulation of Sect.3.1 succeeds at representing at a very good extent the physics ruling the problem, it has a considerable number of variables and constraints that may hinder the capacity to tackle interestingly sized problems.The next model that neglects power curve and AEP calculation is deployed, is intended to simplify it.
The model presented in this section deploys a strategy to account for the combination of Eq. ( 3) and Eq. ( 7) to calculate velocities, since the case studies from the IEA37 Wind Task 37 follow this outlook.It would be possible though to consider the linear superposition model if necessitated.However, this model does not support the application of Eq. ( 4).Combining Eq. ( 3) and Eq. ( 7) and extending the summation range in Eq. ( 7) to all candidate locations, the total wind speed in the farm, U , can be modelled through where new binary variables z iℓ are introduced.The variable z iℓ is equal to one if both WTs i and ℓ are active (i.e. if ) and zero otherwise.Nevertheless, the previous expression is not linear for variable z iℓ due to the presence of the square root in each total relative velocity deficit term.Dropping the square roots, the following expression is obtained: Combining Eq. ( 18) and Eq. ( 19) results in

Ũ =
Total inflow wind speed Eq. ( 20) defines the objective function to be maximized of the power-curve-free model.In comparison to the objective function Eq. (16a), no power curve or continuous state variables are demanded.
Nonetheless, the presence of variables z iℓ can be troublesome.For the complete model, in addition to having these variables of combinatorial nature, constraints of the same kind must be incorporated: results show the heavy computational burden incurred when solving this formulation, impeding to solve large-scale problems (Fischetti et al., 2016).A big-M trick is then incorporated to come up with a model which is exactly equivalent, as reflected in formulation of Eq. ( 21).
subject to: The new objective function in Eq. ( 21a) modifies the component linked to the total wind speed deficit proxy by creating variables τ i ; this variable means total wind speed deficit proxy for WT in candidate location i. Equation (21b) defines τ i , if a WT candidate location is inactive ξ i = 0, then there is no deficit at this location, therefore τ i = 0, because of M i = N ℓ=1:i̸ =ℓ b iℓ , and the minimization nature of the problem for wind speed deficits.Oppositely, if ξ i = 1, then τ i is forced to be equal to Note that for the classic problem definition n min = n max =n T , the first part of the objective function becomes For this situation, the objective function is thus equivalent to minimize 4 Neighborhood Search Heuristic In order to be able to address large-scale problems, a heuristic is proposed to wrap the previously presented MILP formulations in Sect.3. It is based on neighborhood search and local branching theory (Fischetti and Lodi, 2003).The algorithm solves a sequence of MILPs, with different candidates number N or/and neighborhood search size K, taking advantage of robust and efficient implementations of branch-and-cut methods for MILP.
The heuristic relies on the observation that for a fixed layout described by ξ ∈ {0, 1} N , the other design and state variables are straight forward to determine.This observations is valid for all presented problem formulations in Sect.The pseudo code of the Neighborhood Search Heuristic (NSH) is described in detail in Algorithm 1.

9:
Add Hamming distance constraint centered around the incumbent ξ, i:ξi=0 Solve opt.model from algorithm lines (8) to (9) until optimality or computing time T with ξ as warm-starter 11: Get the solution pool S, where ξ ∈ S represents the activation binary variables for WTs of an individual solution Break 25: end if

26: end for
The first three lines are for the main inputs of the algorithm: the candidates set C, the times set T , and neighborhood sizes set V .The first set contains the sizes N of the meshes to be considered, the second one is for the maximum computing time T for the MILP solver for each size N , the last one is for the search size defined as the maximum number of changes K allowed The true objective function refers to the real equation that represents the ultimate aim to be optimized.For example, if this is the AEP, then it is the product of the power calculation process, applying the considered wake and superposition models and the original power curve, and not the objective function of the implemented formulation, as in Eq. ( 16a), which in any case is always an approximation.
The next step is to start the iterative process in line 6.Values for N , T , and K are fetched in line 7, followed by the formulation of the MILP model for candidates N accounting for the points of ξ.The Hamming distance, see e.g.(Fischetti and Lodi, 2003), centered around the incumbent point ξ, is added to the optimization model in line 9; this constraint reduces the search space as the number of changes to ξ are limited to K. The complete model is sent to the MILP solver with ξ as warm-starter, stopped until it reaches optimality or the assigned maximum computing time T .
After solver termination, the solution pool S is retrieved in line 11.It is very important to emphasize the aim of getting the whole pool instead of the best solution.This is done because of the imperfect correspondence between the true objective function and the objective function of the applied MILP model.For example, a solution which may have worse objective value, may actually have a better AEP according to the real model.In this order of ideas, the whole pool is examined, and the best solution indexed by i t with AEP of o t is obtained in line 13.If o t is actually greater than o b , then the whole algorithm is recentered around the new ξ (lines 14 to 16), and in the next iteration κ, the same values of N and K are maintained.Otherwise, the next value of K is taken (line 18), unless all of them have been traversed.In the last situation, the next candidates size N is considered given by countern, restarting the neighborhood set counter counterv to one (lines 20 to 22).The NSH algorithm is terminated when all candidates set C have been processed (line 23 to 25).

Computational Experiments
As mentioned earlier in the manuscript, in order to have a transparent benchmark of the proposed methods, the open access case studies from the IEA37 Wind Task 37 (Baker et al., 2019;Dykes et al., 2015)  The results of the statistical correlation between the proxy function given by the argument in Eq. ( 22), and AEP of the problem definition (Baker et al., 2019;Dykes et al., 2015) for each case are presented in Sect.5.1.After this, the performances of the proposed models in the case studies are deployed in Sects.5.2 (Case I), 5.3 (Case II), 5.4 (Case III).The power-curve-free model of Eq. ( 21) is implemented with the Eq. ( 22 The main parameters of the wake model in Sect.2.1 are fixed to C T = 8/9 and k y = 0.0324555.Furthermore, a wind rose approach for modelling wind resource is utilized in this work.In a wind rose, the wind resource is binned in J directions, and for a specific direction j (θ j ), wind speeds are correspondingly discretized in V sectors.For the case studies, the wind rose is composed of 16 directions and a single wind speed k of 9.8 ms −1 .The considered wind rose is displayed in Fig. 2.
Lastly, for the purpose of replicability of the numerical results, the power curve from Eq. ( 8) modelling the IEA37 3.35 MW reference turbine (with diameter of D = 130 m) is used in this manuscript (IEA Wind Task 37, 2019;Baker et al., 2019).The main parameters are p rated = 3.35 MW, u rated = 9.8 ms −1 , u cut-in = 4 ms −1 , and u cut-out = 25 ms −1 .The power curve for these specific parameter choices is plotted in Fig. 1.
The experiments for Sects.5.2, 5.3, and 5.4 have been carried out on an Intel Core i7-6600U CPU running at 2.80 GHz with four logical processors and 16 GB of RAM.For the Sect.5.5 a larger resource is used for being able to exploit the powercurve-based model capabilities, an Intel Xeon Gold 6226R CPU running at 2.90 GHz with 32 virtual cores and 640 GB of RAM.
The selected MILP solver is the commercial branch-and-cut algorithm implemented in IBM ILOG CPLEX Optimization Studio V20.1 (IBM, 2022).Some parameters are set to different values compared to the default choices.For example, the parameter emphasizing the finding of high-quality feasible solutions early in the process (CPX_MIPEMPHASIS_HEURISTIC) is activated.This is intended to generate more feasible layouts compared to the default setting which is important for the neighborhood search algorithm.Additionally, strong branching is used for variable selection given the large size of the models (CPX_VARSEL_STRONG).This setting is intended to reduce the size of the search tree and thus the memory requirements compared to default settings.

Parameters
The wind turbine is the IEA37 3.35 MW onshore reference turbine [1] with the following characteristics: Rotor All turbine data is also contained in the enclosed iea37-335mw.yaml.The power curve is defined as:

P (MW)
The farm wind speed for all scenarios is constant at 9.8 m/s.The +y axis is coincident with 0 • , and the CW wind rose is defined by 16 discrete bins tabulated in iea37-windrose.yaml,depicted pictorially below:

Case Study 1: Optimization Only
This problem defines three different wind farm sizes, and corresponding number of turbines, intended to test scalability of your optimization approach.The three scenarios are: 1. 16 turbines, boundary radius of 1,300 m.
For this Case Study the user is only free to choose the optimization approach.The wake model is fixed and is a simplified version of Bastankhah's Gaussian wake model [2,3,4].A Python implementation is supplied for convenience (iea37-aepcalc.py).Alterations to this implementation are permitted, as long as the In the above figure is noticeable that the boundary is densely sampled, as a candidate point is defined every natural angle from 0 • to 359 • , i.e. 360 candidate points are provided.This is done as it is intuitively expected that a good portion of the WTs will be placed in the boundaries, as WTs will be spread out in the available area as much as possible to decrease wake losses.
For the interior, a set of finite parallel line segments are generated in between of the extreme abscissa and ordinate points of the available project area.The candidates points are then taken along those segments, according to sampling parameters, such as, distance between the points along a line, and vertical distance between lines.In the example of Fig. 3, the slope of the line segments is zero, and the distance between points and lines is equal to 1.7D.Although in this example a circle with radius of 1300 m is displayed, the underlying principles apply for any project area shape.

Correlations
To validate the approach modelled by the MILP formulation of Eq. ( 21) (i.e. the power-curve-free model), 5000 random feasible WT layouts are created.For each of them, the AEP (Baker et al., 2019;Dykes et al., 2015), the total wind speed theoretical, U , given by Eq. ( 17), the total wind speed proxy, Ũ , defined by Eq. ( 20), and total wind speed deficit proxy, N i=1 τ i , argument of Eq. ( 22), are calculated.Although the random way of generating the layouts biases the obtention of high-quality solutions, the general trend is the aspect of interest and can be assumed to be fairly representative for the whole domain.
In all cases Pearson product-moment linear correlation coefficients (Pearson, 1895) are used to extract information from the data, the results are collected in Table 1 for every single pair of the aforementioned measures.This coefficient illustrates the degree to which the movement of pairs of variables is associated in a linear fashion.The correlation plots of Fig. 4 present the graphical representation of the relations for Case I.The correlation between AEP and the total theoretical wind speed is shown in Fig. 4a for the Case I.The main observation is the very strong linear relation between these two variables as illustrated by the correlation coefficient of 0.97.Interestingly, this reflects the rather low influence of the WT power curve to determine high-quality feasible points.The next observation is present in Fig. 4b.In this plot, the relation between U and Ũ is represented, resulting in almost an identical linear connection between them, as in the previous graph.When one looks into AEP vs N i=1 τ i , however, it is noticeable that the Pearson coefficient worsens, decreasing to −0.88.There is a wider area in the formed body of points that causes this behaviour.Note that in contrast to the previous two figures, there is a negative correlation because the comparison is done in terms of wind speed deficit instead of total wind speed.In spite of this deterioration, the linear correlation is still strong enough.These results motivate an approach where the minimization of a proxy total wind speed deficit can lead to high-quality AEP solutions.The NSH Algorithm 1 helps correcting the imperfect correspondence between these two variables during the optimization routine as reflected in Sect. 4.
For Case II the general trends of the correlation plots are very similar to the observed in the previous case.Correlations between AEP versus total wind speed theoretical (0.97), and total wind speed theoretical versus total wind speed proxy (0.95) are still very strong.Nonetheless, there is a slight worsening between the relation of AEP vs total wind speed proxy (it moved from −0.88 in the previous case to −0.85), as the spread for middle velocity values is enlarged.The linear relation is deemed as satisfactory enough to carry on with the application of model of Eq. ( 21).
Lastly, for Case III yet another very strong the linear relation between AEP and the total wind speed theoretical (0.96) is observed, just as for the two previous situations.This is a very interesting outcome of this manuscript.Almost all research in the WFLO problem strictly focuses in power modelling, which brings a great deal in complexity due to the non-linear and non-differentiable properties of a typical WT power curve.Using an exact model for determining total wind speed as objective function alleviates the computational complexity, while being able to find high-quality solutions in terms of AEP for the classic WFLO problem.
Nevertheless, deterioration in the correlations stemming from the proxy to calculate total wind speed deficit is present in this case.Both with the total wind speed theoretical (0.88) and the AEP (−0.72).Keep in mind that the reason to formulate such approximation is to fit in the context of integer programming to leverage theory and state-of-the-art algorithms of this mature field by having a compact formulation.However, the price to pay is to lose fidelity to represent the real (true) target to optimize.
The deterioration in the correlation of these pairs of variables may also suggest the need to resort to the power-curve-based model for some applications.Whether the price is too high or not is reflected in the reachable solution quality.Sects.5.2, 5.3, and 5.4 present the optimization results for the cases of fixed number of WTs that will ultimately help to elaborate a final evaluation regarding the adopted modelling technique.

Case I: 16 WTs
The performance evolution of two of the proposed optimization frameworks in this article is depicted in Fig. 5.The green line of the full model is obtained after solving the model of Eq. ( 21) with objective function as in Eq. ( 22) for N = 1014 without implementing the NSH.In order to plot this line, the CPLEX's log is postprocessed by finding out the incumbent solution in 360 terms of AEP and not for total wind speed deficit proxy.The blue line results after applying the NSH with the model of Eq. ( 21) plus objective Eq. ( 22 they start increasing to refine the search.The red line is for establishing a reference of AEP value, this comes from the best performing method in the survey of IEA37 Wind Task 37 (Baker et al., 2019), the SNOPT plus Wake Expansion Continuity (WEC) (Thomas and Ning, 2018;Thomas et al., 2022b).Time evolution for the SNOPT+WEC is not reflected in this graph, as this information is unavailable.Results for the benchmark against a testbed of different algorithms are available in Table 2  Table 2. Results for all three benchmark cases from other algorithms (G, gradient-based and GF, gradient-free) obtained while allowing turbine locations to vary continuously.The AEP values for each case/algorithm in GWh are reproduced from the IEA37 study (Baker et al., 2019).The difference column shows how the proposed heuristic with the power-curve-free model performs in comparison in %.Negative percentages means that the proposed method performs better than the corresponding algorithm.The benefit of the proposed neighborhood search strategy is shown in Fig. 5. Solving the full model is significantly slower, leading actually to a worse solution (3.31% lower).The capacity of the NSH to iterate over different values of candidate points N and search size K brings alone not only improvements in terms of solution time and solution quality, but also in less computational resources as the RAM memory generally escalates faster when solving the single model.

Method
The initial and final solution layouts for this case study are illustrated in Fig. 6.The importance of finely sampling the boundaries of the available area is evident in Fig. 6b, because 7 out of the 16 WTs are placed in that subdomain.
Finally, Table 2 compares the proposed method to a large number of different approaches from the IEA37 reference study (Baker et al., 2019).The results for all case studies are presented in this table, where I, II, and III make reference to cases from this section, Sect.5.3, and Sect.5.4, respectively.For details of each of the benchmark algorithms, the reader is referred to the mentioned reference.
The third column of Table 2 reports the difference of the AEP with respect to the proposed method for the smallest case study.The obtained AEP in this case is better than almost all the other alternatives, except to the SNOPT+WEC, where a nearly identical objective value is achieved.When directly comparing to the gradient-free (GF) methods, the best found solution by them (full pseudo-gradient with 402318.75MWh) is determined in around 2 h by the proposed method, which is typically way faster than the overall performing computing time of these kind of algorithms.Discrete optimization approaches, as the MILP ones presented in this article, could be formulated to cope with problem definitions with required functionalities that in theory continuous optimization methods can not support (or at least the implementation becomes strenuous).
The power-curve-based model of Eq. ( 16) within the NSH using the same AEP formulation as true objective function, provides a solution 1.18% lower in objective value in around 36 h using the computer system with 32 virtual cores.Although the quality of the layout is very close to the one schematized in Fig. 6b, the need for larger computational resources tips the scales for implementing the power-curve-free model for these type of problems with fixed number of WTs.For this reason, in Sects.5.3 and 5.4 are presented only the results reached after the application of the power-curve-free model embedded into the NSH. -

Case II: 36 WTs
The performance evolution of the proposed methods, and the initial and final WT layouts are plotted in Fig. 7     Opposed to Fig. 5, in Fig. 7 is distinguishable that the full model (i.e.without implement the NSH algorithm) initially provides better solutions within the first 3 h, but then it remains in a considerably worst solution than when applying the NSH algorithm in the long run (lower 3.05%).
In this case, the proposed method reaches the best solution, as given in the fifth column of Table 2.The SNOPT+WEC is again the closest contender.When uniquely comparing to GF methods, the proposed method matches the best solution from those algorithms in around 3 h, which is generally a very fast computing time when contrasted to other methods where gradients are not explicitly utilized in the optimization process.

Case III: 64 WTs
The performance evolution of the proposed methods, and the initial and final WT layouts are displayed in Fig. 9    Seventh column of Table 2 shows that the SNOPT+WEC and the preconditioned SQP provide slightly better layouts than the proposed method for this case.However, the algorithm provides feasible layouts that improve the objective compared to all the listed gradient-free approaches.

Case IV: 10-50 WTs
Although in most projects today the total capacity for grid connection is decided already in the early planning phases, in the future one can envisage situations where flexibility in optimizing the number of wind turbines in a project would yield benefits.
As much as the power-curve-free model presented in Sect.3.2 exhibits a quite good performance both in terms of AEP and computing time for the scenario of fixed number of WTs, it is not very well suited for optimizations aiming at improving economic indices, like NPV.For such an optimization, the power-curve-based mathematical program of Sect.3.1 may be handy as the number of generators is allowed to vary between a lower and upper bound, n min and n max , respectively.For illustration, a domain defined by a circle with radius 1300 m, and variable number of WTs between 10 and 50 are utilized.
These parameters are set relatively arbitrarily but with sufficient distance to reasonably expect that the limits are not reached.Three runs launching the model of Eq. ( 16) with modified objective function Eq. ( 23), embedded in the NSH Algorithm 1 with NPV as target function are executed.For the first example the number of turbines is fixed to n min = n max = 10.
In the second example the number of turbines remains fixed but is increased to n min = n max = 50.When the number of turbines is fixed to 10, the NPV evolution (Fig. 11b) is driven by the AEP (Fig. 11a Whether there is a trade-off point in between of both extreme cases is the main question to answer.The evolution of the WTs number in Fig. 13a and the AEP in Fig. 13b exhibits a perfect correspondence.The more WTs the larger AEP, in spite of the increased wake losses.The curves increase in time, up to a point where the model estimates that further increase of WTs would not lead to a better NPV.The final number of WTs is 34.The NPV evolution in Fig. 13c naturally only improves with time, resulting in a final value of 795.86 Mill.Eur.Note that the NPV in this case is greater than when a larger number of WTs (i.e.50) was considered and of course when only 10 were considered.
This result shows the benefit of having optimization models that support variable number of WTs accounting for metrics beyond AEP.The advantages may become even more pronounced for more complex situations, as for instance, if the WT investment costs are dependent on the exact installation area or different WT sizes are considered.

Discussions and future work
The two models proposed in this article have many of the common characteristics of mixed integer linear programming models.
They require significant resources like computational time and memory.Furthermore, they exhibit rather low tractability and scalability for global optimization algorithms.
The power-curve-based model, albeit requiring large computational resources as number of cores and RAM memory, manages to provide reasonably good solutions for small-sized problem instances, being only 1.18% worse than its power-curve-free counterpart for the 16 WTs case and 4.41% for the 36 WTs case.This diminishing effectivity is to be expected, given the large number of variables and constraints.The power-curve-free model on the other hand along with the heuristic is much faster due to its associated more compact formulation.This translates into the ability to be highly competitive compared to a large set of benchmark algorithms.In situations where there is an interest for optimizing metrics beyond AEP, such as the NPV, the power-curve-based model becomes very useful given its intrinsic capacity to support this kind of objective functions, as shown in the final test case.
It is relevant to mention that there are limitations in the used wake models compared to recent models from the literature (Thomas et al., 2022b).For example, the used wake model does not consider the modifications of the turbulence intensity or thrust coefficient variations from wind speed inside the wind farm.It is uncertain if these modifications would still allow an integer linear programming formulation or approximation of the WFLO problem.It is also uncertain what impact on the final solution quality these detailed modelling aspects imply.These questions are left for future work.
Notwithstanding the listed shortcomings, it is very enthralling that these models in combination with the neighborhood search heuristic are able to match and in some cases improve the results obtained when considering the turbine positions as continuous variables (see Table 2).This opens the door to trying out new case studies with functionalities easily adaptable to discrete parametrization techniques, which can be very challenging for continuous variable modelling approaches.Future work can include the application of the proposed methods to real-world problem instances with for example, forbidden areas, complex turbine cost functions or integrated optimization with electrical systems.Furthermore, the assessment of the robustness of the method for different initial layouts can be assessed in following studies.

Conclusions
This manuscript contributes both methodologically and empirically to address the WFLO problem.It is very difficult to cope with extended versions of the classic problem definition.Relevant extensions include variable number of WTs or complex forbidden areas.For these situations there is a need to develop new modeling and algorithmic approaches in this direction.
A neighborhood search heuristic embedding integer programming formulations is therefore proposed to satisfy these new demands.For both presented formulations, the step-wise power curve and power-curve-free, the heuristic notably improves a single execution of full models when calling a state-of-the-art branch-and-cut solver in terms of solution quality.An improvement of up to 3.42% in the AEP is achieved by applying the neighborhood search strategy for cases where the WTs number is fixed compared to solving the full model.
Another important takeaway is the satisfactory performance of the power-curve-free model, which uses an approximation of the total wind speed deficit, when (implicitly) optimizing for AEP.This is due to the good correlation between the two measures, and the correction capability of the heuristic.For the classic WFLO problem definition, the proposed model is able to considerably enhance (from 1% to around 10%) the AEP compared to benchmark layouts obtained by several different gradient-based and gradient-free algorithms.Even when directly compared to methods implementing a continuous variables technique, the proposed heuristic provides very competitive results, being able to come up with layouts with similar objective function values, or even improving it.These are very promising results that would enable to get high-quality solutions for problem instances where continuous variables modelling approaches may not be able to run or come up with good incumbents.
Finally, the model with explicit representation of the power curve embedded within the neighborhood search heuristic is able to propose non-trivial solutions when implementing an objective function beyond AEP, such as NPV.For these cases the trade-off between energy revenues and investment costs is inherently studied.For the case study, the model suggests that is not worth to install the maximum number of allowed wind turbines, but a lower value that reflects a better NPV value.

Figure 1 .
Figure1.Piece-wise constant approximation of a wind turbine power curve through sampling with m = 10 intervals between the cut-in and rated wind speeds.
3. Given ξ ∈ {0, 1} N , for the power-curve-based model, the continuous state variables u can be determined through classical wake analysis, and the 10 https://doi.org/10.5194/wes-2022-82Preprint.Discussion started: 20 September 2022 c Author(s) 2022.CC BY 4.0 License.binary state variables η are directly determined by inspection of the velocities.Similarly, for the power-curve-free model, the τ variables are trivially computed.
250 to the incumbent solution.11 https://doi.org/10.5194/wes-2022-82Preprint.Discussion started: 20 September 2022 c Author(s) 2022.CC BY 4.0 License.If the incumbent is improved, then the candidates set C, and neighborhood size K are kept, otherwise at least one of them is increased.The first step (line 5) is to obtain an initial incumbent binary variables, with the set ξ storing the acquired value (0 or 1) for each variable ξ i : i ≤ N .The incumbent has an objective value of o b calculated after the true objective function.
are used for comparison.Those cases have circular project areas with three different radius (1300 m, 2000 m, and 3000 m) and number of WTs (16, 36, and 64), n T .Thus, Case I has a radius of 1300 m and n T = 16 WTs, whereas Case II has radius 2000 m and n T = 36, and Case III has radius 3000 m and n T = 64, correspondingly.
) as objective function.The true objective function in the NSH Algorithm 1 for these cases is the AEP of the problem definition.In the end, to prove the capabilities of power-curve-based model of Eq. https://doi.org/10.5194/wes-2022-82Preprint.Discussion started: 20 September 2022 c Author(s) 2022.CC BY 4.0 License.(16), Sect.5.5 displays results after applying this formulation with a modified objective function to express a metric similar to NPV.

Figure 2 .Figure 3 .
Figure 2. Wind rose used in the computational experiments.Taken from open access source (IEA Wind Task 37, 2019).
Correlation between the theoretical and approximated wind speeds.
Correlation between AEP and the approximated wind speed deficit.

Figure 4 .
Figure 4. Correlation plots for 5000 randomly generated wind turbine layouts for Case I.
), and AEP as true objective function in Algorithm 1.The main inputs are C = {467, 590, 1014}, 16 https://doi.org/10.5194/wes-2022-82Preprint.Discussion started: 20 September 2022 c Author(s) 2022.CC BY 4.0 License.T = {1, 1.5, 2} h, V = {2, 4, 6, 16}.These inputs are tuned after evaluating the performance of the method using different values.In general, the first two elements of C consists of moderately big values, relatively close to each other, while the last element is sizeably greater in the seek of the best possible solution.Each element N ∈ C has associated a computing time T .Finally, the first elements of V are relatively low values to favour termination of the solver due to optimality, and then

Figure 5 .
Figure 5. Performance of two different optimization approaches for Case I and comparison with existing best benchmark results.
Final wind farm layout obtained by the heuristic.

Figure 6 .
Figure 6.Generated wind farm layouts for the benchmark Case I with 16 turbines.
and Fig. 8, respectively.Main inputs are C = {477, 684, 1907}, T = {1, 1.5, 2} h, V = {2, 4, 8, 16, 36}.The blue line (model of Eq. (21) with objective function Eq. (22) plus NSH Algorithm 1) has clearly three sectors stemming from each value of N ∈ C. The initial WT layout (Fig. 8a) -also determined by choosing roughly equidistant candidate locations in the boundary -has an AEP of 796514.61MWh.After seven iterations of the NSH (point 8) in 41 s, the incumbent is improved by 1.84%, when N = 477 and 2 ≤ K ≤ 4, being able to solve each model instantation to optimality.As shown in Fig.7, after a three-hours-plateau linked to 8 ≤ K ≤ 36 (four iterations), N is raised to 684, resulting in the sharpest AEP enhancement.The energy production increases with 4.51% after only 23 min in point 27.This noticeable improvement comes after solving to optimality models with rather small neighborhood search sizes 2 ≤ K ≤ 4. The convenience of allowing large neighborhood search sizes as K = 16 or K = 36 is reflected from this moment.From point 30 to 33 (6 h) with K = 16 the incumbent is slowly boosted by nearly 1%.Again, after a three-hours-plateau, N becomes equal to 1907, and https://doi.org/10.5194/wes-2022-82Preprint.Discussion started: 20 September 2022 c Author(s) 2022.CC BY 4.0 License. in around 32 min for 2 ≤ K ≤ 4, the AEP is augmented by 0.41%.Then, the large neighborhood search starts for K = 16 and K = 36, and after a total of 16 h, the final solution of 865327.78MWh (increment of 0.61%) is achieved (Fig.8b).
Final wind farm layout obtained by the heuristic.

Figure 8 .
Figure 8. Generated wind farm layouts for the benchmark Case II with 36 wind turbines.
and Fig.10, respectively.Main inputs are C = {625, 1017, 2741}, T = {1, 1.5, 2} h, V = {2, 4, 8, 16, 32, 64}.Note that in comparison the number of elements of V has been increased by one after each study case.This has been done taking into account the number of WTs.Likewise, the values of N ∈ C are greater to cover for the wider project areas.

Figure 9 .
Figure 9. Performance of two different optimization approaches for Case III and comparison with existing best benchmark results.

Figure 10 .
Figure 10.Generated wind farm layouts for benchmark Case III with 64 wind turbines.
22 https://doi.org/10.5194/wes-2022-82Preprint.Discussion started: 20 September 2022 c Author(s) 2022.CC BY 4.0 License.The aim is to illustrate the ability of the method in reaching non-trivial solutions, resulting in a an optimized design with an an intermediate number of wind turbines.Keep in mind that for this case, a linear superposition model for the AEP component in the NPV calculation is considered.In this sense, the original WT power curve as deployed in Fig.1is supported.NPV is the true objective function when applying the NSH Algorithm 1.The modified objective function of MILP model of Eq. (16) for this case has the form (Investopedia, c wt is the cost per WT in Mill.Eur, c e the energy price in Mill.EurMWh −1 , r is the discount rate in %, and Y is the number of years of lifetime of the project.For this case study, values of c wt = 6.7 Mill.Eur(Nord Pool, 2022), c e = 0.00015 Mill.EurMWh −1 (Mishnaevsky Jr and Thomsen, 2020), r = 5%, and Y = 20 are assumed.

Figure 11 .
Figure 11.Evolution of the AEP and NPV with number of wind turbines fixed to 10 applying the power-curve-based model.

Figure 12 .
Figure 12.Evolution of the AEP and NPV with number of wind turbines fixed to 50 applying the power-curve-based model.

Figure 13 .
Figure13.Evolution of the AEP and NPV with variable number of wind turbines between 10 and 50 applying the power-curve-based model.
).Both curves are monotonically increasing, reaching a final value of NPV of = 456.40Mill.Eur.The same behaviour is visible for n T = 50, 24 https://doi.org/10.5194/wes-2022-82Preprint.Discussion started: 20 September 2022 c Author(s) 2022.CC BY 4.0 License.although the final NPV is greater (683.53Mill.Eur), see Figure 12.The difference in discounted cash flow revenues surpasses the associated extra investment costs from the additional number of wind turbines.

Table 1 .
Pearson product-moment linear correlation coefficients for all case studies.
CaseAEP vs Theoretical wind speed Theoretical vs Proxy wind speed AEP vs Proxy wind speed deficit .
Performance of two different optimization approaches for Case II and comparison with existing best benchmark results.
(a) Initial wind farm layout provided to the heuristic.