Modelling Wind for Wind Farm Layout Optimization Using Joint Distribution of Wind Speed and Wind Direction

: Reliable wind modelling is of crucial importance for wind farm development. The common practice of using sector-wise Weibull distributions has been found inappropriate for wind farm layout optimization. In this study, we propose a simple and easily implementable method to construct joint distributions of wind speed and wind direction, which is based on the parameters of sector-wise Weibull distributions and interpolations between direction sectors. It is applied to the wind measurement data at Horns Rev and three different joint distributions are obtained, which all fit the measurement data quite well in terms of the coefficient of determination (cid:1844) (cid:2870) . Then, the best of these joint distributions is used in the layout optimization of the Horns Rev 1 wind farm and the choice of bin sizes for wind speed and wind direction is also investigated. It is found that the choice of bin size for wind direction is especially critical for layout optimization and the recommended choice of bin sizes for wind speed and wind direction is finally presented.


Introduction
In the last two decades, the progress of technologies, together with the increased experience of wind farm (WF) construction and operation, has enabled the development of modern WFs, which typically consist of tens or hundreds of utility-scale (multi-MW) wind turbines (WTs) and can have a total capacity of hundreds of MW. In parallel with this trend, the efforts to increase the percentage of wind power in the energy mix have led to the proliferation of modern WFs [1].
One of the critical problems for modern WF development is layout optimization, which seeks to determine the optimal positions of WTs inside the WF by maximizing and/or minimizing a single objective or multiple objectives, while satisfying certain constraints. This problem has been investigated in many studies with different problem formulations and using various optimization algorithms [2]. Most of these studies focused on developing and improving the optimization methods, with relatively less attention paid to the appropriate modelling of wind.
In their seminal work published in 1994, Mosetti et al. constructed an ideal test problem of WF layout optimization and solved it with genetic algorithm (GA) [3]. They used three wind cases: (1) uniform north wind with a speed of 12 m/s; (2) equally distributed (36 directions) wind with a speed of 12 m/s; (3) non-uniformally distributed (36 directions) wind with speeds of 8, 12 and 17 m/s. Clearly these ideal wind cases are quite simple compared with the real wind, and they do not need much consideration in wind modelling. The same ideal wind cases were used in many following studies which mainly aimed at developing various optimization methods, including GA [4], particle swarm optimization [5], extended pattern search [6], and so on.
As the energy source for wind power, wind needs to be characterized appropriately in order to obtain a profitable utilization of wind energy. For the wind speed variation, different statistical models were proposed [7], among which Weibull distribution is the most widely used one. For the wind direction variation, the frequencies of occurrence of different direction sectors are usually used, typically presented in a wind rose. Combining these two, the measured wind data can be fitted into sector-wise Weibull distributions, which is the common practice in wind modelling for wind resource assessment and annual energy production (AEP) calculation. This more realistic wind modelling method has been adopted in several studies on WF layout optimization. Some studies used 12 sectors for wind direction [8,9], while some studies used 24 sectors [10,11].
Modelling wind with sector-wise Weibull distributions has been proved to be suitable for wind resource assessment or AEP calculation for a given WF in the industry's past experience. However, in our previous study [12], we found that this modelling method might cause problematic results when optimizing the layout of a large offshore WF, with 12 or even 72 sectors for wind direction. The reason lies in the discreteness and discontinuity of wind direction in the wind modelling. Additionally, in a numerical study of the same WF [13], Porté-Agel et al. demonstrated that the power production of a large WF is highly sensitive to wind direction due to the complex wake effects between WTs. These studies suggest that more advanced wind modelling methods that can better characterize wind direction variations are quite needed, especially for application in WF layout optimization.
Considering the fact that the wind speed and wind direction are not independent random variables, we can expect that a bivariate distribution of both wind speed and wind direction would characterize the wind better and also be beneficial for WF layout optimization based on its continuous nature in both dimensions of speed and direction. Such distributions have been proposed by different studies, which mainly aimed at wind modelling. Carta et al. [14] proposed to construct a joint distribution from two marginal distributions, i.e., truncated Normal-Weibull mixture distribution for wind speed and a finite mixture of von Mises distributions for wind direction. Erdem and Shi [15] constructed seven different bivariate joint distributions using three construction approaches, namely, angular-linear, Farlie-Gumbel-Morgenstern and anisotropic lognormal approaches. Zhang et al. [16] presented a multivariate distribution of wind speed, wind direction and air density by using a non-parametric approach, multivariate kernel density estimation. However, these joint distributions have not yet been widely known or used in the wind energy field. One possible reason might be that the method of sector-wise Weibull distributions has been proved to be adequate for the tasks such as wind resource assessment and AEP calculation, and has become the widely adopted common practice.
In this study, with the WF layout optimization as the targeted application, we propose a simple and easily implementable method to construct joint distributions of wind speed and direction based on the parameters of sector-wise Weibull distributions. Since in the common practice these parameters can usually be obtained from a wind measurement campaign and are often presented in the wind measurement report, this method is quite convenient for a WF developer to use. Three types of joint distributions are constructed using this method. We apply this method to the 3 years' measurement data at Horn Rev 1 and use one of the obtained distributions, namely, spline joint distribution, in the layout optimization of the Horns Rev 1 WF. We also investigate the choice of bin sizes for wind speed and wind direction in numerical calculation, both for WF power calculation and layout optimization, and find that the choice of bin size for wind direction is especially critical for layout optimization. The recommended choice of bin sizes for wind speed and wind direction is presented after the case study for Horns Rev 1 WF.

Wind Farm Layout Optimization
Wind farm layout optimization can be formulated as a general optimization problem, which seeks to find the optimal layout by minimizing a cost function subject to certain constraints. Considering a WF with WTs located at , , … , , , , … , , we can denote the layout as , and write the layout optimization problem in the following form: where , denotes the cost function of layout optimization, , represent the constraint functions and is the number of constraints.
The cost function , can be derived for different optimization objectives, such as maximizing the power, AEP, profit, net present value, or minimizing the cost of energy (CoE), levelized production cost (LPC) [2]. In all these cost functions, one of the essential components is the expected power output of the WF. Note that for a given WF site, the wind resource is given; if we assume the type of WTs is also given, then the expected power output is a function of the layout , , i.e., , . For maximizing the power of a WF with WTs, the cost function can be defined as , , .
The constraints in WF layout optimization may come from various considerations, such as technical, environmental and economic ones. The boundary of WF and the minimal distance between WTs are the two commonly considered constraints.
In order to calculate , the wake effects between WTs need to be modelled appropriately. The most widely used wake model is the Jensen wake model [17], which was developed by assuming conservation of momentum within the wake and linear expansion of wake region. In this study, we use the same problem formulation and modelling methods as in our previous study [12], in which the detailed calculation methods for power and constraints can be found.

Wind Modelling
Wind resource assessment is the starting point of a WF development project, since wind resource will mostly determine the amount of power production and therefore the profit of the WF in its lifetime. To get a reliable assessment, it is typical to first carry out a wind measurement campaign at the planned WF site and collect a large quantity of wind data. The measured wind data might be obtained at a reference height and then used to predict the power production of WTs with hub height . It is thus necessary to first convert the measured wind speed into the inflow wind speed at hub height , usually using the logarithmic law: where denotes the wind speed at the reference height and is the surface roughness length.
The converted wind data can be processed using the method of bins [18]. Considering all wind directions and the interested speed range , for power production, where is the cut-in wind speed, is the cut-out wind speed, we can first discretize the interested wind conditions into two dimensional bins by using bin size for wind speed and for wind direction. Then the middle point of each bin can be denoted as: where , are the numbers of bins for wind speed and wind direction respectively, which are determined by the bin sizes and . Note that in Equation (3) the range of interest for wind direction is deliberately chosen as 0.5 , 2 0.5 , so that the middle point of the first bin is always at wind direction 0°, which coincides with the common practice in wind measurement. Then the converted wind data can be summarized in matrix form as: (4) where , ) denotes the frequency of occurrence of the ijth wind bin ( , . Based on the frequency of occurrence, we can also calculate the bivariate probability density function (PDF) of wind speed and wind direction at the discretized points , as: In the common practice, the wind measurement data can also be processed using tools like WAsP [19] to obtain the wind rose, which describes the probability distribution of wind direction variation, and a fitted Weibull distribution for every direction sector (usually 12 sectors with width of 30° each), which describes the wind speed variations.
In both methods, it is necessary to choose appropriate bin sizes and , as discretization is always needed when calculating the expected power output of the WF by using numerical integration. Usually in AEP assessment, 1 m/s and 30° are used. But in the context of WF layout optimization, we have shown in Ref. [12] that a large number of direction sectors, i.e., much smaller than 30°, have to be used in order to get consistent and reliable results.
To better demonstrate the importance of choosing appropriate , we can construct an ideal WF that is composed of three Vestas V80 WTs. The three turbines are located on the edge of a circle with radius of 400 m and placed with equal distances between any two ones. Two possible layouts of the WF and the characteristics of the WTs are shown in Figure 1. Note that in Figure 1a, the small colored circles represent the covering area of the turbine rotor and Layout 2 is obtained by rotating Layout 1 with 15°. Assume the wind condition here is represented by a constant wind speed 10 m/s with uniformly distributed wind directions, i.e., the wind blows from any direction with the same possibility. Observing the characteristics of these two layouts and the wind condition, we can conclude that Layout 1 and 2 will have the same expected power output.
If we choose 30°, then we are assuming 12 wind directions as shown in Figure 1a. Using the Jensen wake model and assuming the wake decay coefficient 0.04, we can plot the wake zone developments of WT 1 along certain wind directions, which are shown in Figure 2 for both layouts. We can easily see from Figure 2 that for Layout 1, the wake effects will reduce the expected power output, whereas for Layout 2, the expected power output will not be reduced. If we use this kind of wind modeling ( 30°) and calculate the expected power production, we will come to the conclusion that Layout 2 will produce more power than Layout 1, which contradicts with the observation we have made before, i.e., these two layouts are identical in terms of power production.
We can put this constructed ideal case in the context of WF layout optimization. If we use the common practice, i.e., 12 sectors for wind direction, the optimization algorithm might give an 'improved' layout (Layout 2) over the original one (Layout 1). This might be achieved by moving the locations of each WT to escape the wake effects of other WTs, since there are some big 'gaps' between different wind directions. But in reality wind comes from all possible directions instead of only the discretized 12 directions, so the achieved 'improvement' might actually be artificial, as in this constructed case; or even deterioration, as shown in Ref. [12]. Therefore, it is of crucial importance to choose small enough value for , so that we can obtain consistent and reliable results for layout optimization.

Data Source
The wind farm investigated in this study is the Horns Rev 1 offshore wind farm, which is located about 15 km off the westernmost point of Denmark and has a rated capacity of 160 MW. It consists of 80 Vestas V80 WTs, which have a rotor diameter 80 m and a hub height 70 m . Its construction started in February 2002 and its commissioning occurred in December 2002. Before the construction, meteorological measurements had been conducted at Horns Rev since May 1999, using primarily a wind mast. Three years' (01/06/1999-31/05/2002) wind measurement data from this wind mast are used, which is imposed of statistical 10-min mean values for wind speed and wind direction and has 145,048 data points in total. Since the used measurement data is from a period when no WT has been in operation, it can be treated as clean data, i.e., wind data without wake effects involved. With the assumption that the three years' wind characteristics are representative for the following 20 years, this set of data is suitable for evaluating the power production of the WF in Horns Rev.
Note that the wind speed is measured at 62 m and will first be converted into wind speed at the hub height using Equation (2). The converted wind data and the corresponding total power of the WF, calculated using the same method and code as in Ref. [12], are shown as time series in Figure 3. The expected power output can be calculated as 78.79 MW by averaging the simulated time series of power. This can be viewed as a direct way of using the wind measurement data for power calculation, but it is much more time consuming than using statistical wind modelling method. If "one WF simulation" represents simulating the power production of a given WF for a given wind speed and wind direction once, then using a time series of wind measurement composed of data points to calculate the expected power requires WF simulations. For the Horns Rev 1 site, calculating the expected power output of a given layout once will need 145,048 WF simulations, which in turn puts a quite high computation cost for the optimization process. From this perspective of view, the statistical modelling of wind is not only needed for characterizing the long term wind conditions, but also beneficial for reducing computation cost in layout optimization.
The converted wind data at hub height can also be processed with wind rose and Weibull distribution, which is a widely used probability distribution for wind speed modelling and governed by: where is the scale parameter and is the shape parameter. Assuming that the mean value ̅ and the standard deviation of wind speed have been calculated from the data, we can estimate the two parameters as [18]: where is the gamma function. Choosing 30°, i.e., using 12 direction sectors, we can obtain the wind rose and Weibull distribution of the 3 years' wind data as shown in Figure 4.  For each of the 12 direction sectors, the same procedure can be used to derive 12 sector-wise Weibull distributions, the obtained parameters and the frequency of occurrence for each sector are listed in Table 1.

Construction of Joint Distributions
In this section, we present a method for constructing joint distributions of wind speed and wind direction, which gives the bivariate probability density functions (PDFs) for three wind distribution models. This method is based on the parameters obtained by fitting the measurement data into sector-wise Weibull distributions, as shown in Table 1, and it can be implemented easily. Three variants of this method are presented, in which one obtains a piecewise bivariate PDF, while the other two obtain continuous bivariate PDFs by using parameter interpolation.

Piecewise Bivariate PDF
Denoting the number of direction sectors as 12 , we can write the processed wind resource data shown in Table 1 as: , , , with 1, 2, … , . Note that when modelling wind with sector-wise Weibull distributions as described above, we are assuming that the wind speed satisfies the same probability distribution inside a direction sector.
Based on the sector-wise parameters , , , , the joint distribution of wind speed and wind direction can be constructed as: , , , 360/N , when is in the th direction sector, which describes the piecewise bivariate PDF of the joint distribution. Using the converted wind data and choosing a discretization of 1 m/s and 10°, we can calculate the bivariate PDF at the discretized points using Equation (5) as , . The bivariate PDF at the same points can also be obtained by the constructed joint distribution using Equation (8).
In order to evaluate the goodness-of-fit of the proposed piecewise joint distribution and the recorded data, the coefficient of determination , can be used.
is calculated from two data series: one is the observed values , , … , , and the other is the corresponding values , , … , that are estimated by the model. Suppose is the mean of the observed data, then can be calculated as: In this paper, is calculated for the bivariate PDF values at discretized points, i.e., , from the measurement data and , from the piecewise joint distribution. For discretization with 1 m/s, 10°, we get 0.9116, which suggests the model fits the data quite well. Two joint distributions of wind speed and wind direction, calculated from measurement data and modelled by the proposed piecewise joint distribution respectively, are shown in Figure 5. It can be seen in Figure 5 that the overall shape of the observed joint distribution of wind speed and wind direction has been quite well captured by the proposed piecewise joint distribution, and also some of the sharp variations shown in measurement have been smoothed. It is also worthy to note that the piecewise joint distribution shows significant discontinuity between different direction sectors, which is an artificial phenomenon introduced by the way this model is constructed.

Continuous Joint Distributions
In order to avoid the discontinuity in wind direction introduced by the piecewise model, we use interpolation to better utilize the sector-wise parameters, i.e., , , , . Two interpolation methods are applied here: one is linear interpolation, and the other is spline interpolation. In this study, the spline interpolation is done by using the cubic spline data interpolation (spline function in Matlab). From interpolations, we can obtain Weibull parameters and frequency of occurrence as continuous functions of wind direction: , and , while in the piecewise joint distribution, they are all assumed to be piecewise functions. These different functions are shown against each other in Figure 6. which describes the continuous bivariate PDF constructed by using linear interpolation. Similarly, the continuous bivariate PDF constructed by using spline interpolation is governed by where , and represent the functions obtained from spline interpolation.
Using the same procedure as in Section 4.1, we can also calculate the coefficient of determination for these two continuous PDFs, which is 0.9208 for the model using linear interpolation, and 0.9224 for the model using spline interpolation. These bivariate PDFs are shown in Figure 7. Note that we have proposed three joint distributions of wind speed and wind direction, which are described by Equation (8), Equation (10) and Equation (11), and they can be called as "piecewise joint distribution", "linear joint distribution" and "spline joint distribution", respectively. Comparing these three distributions, we can conclude that the spline joint distribution fits the measurement data best.

Application in Wind Farm Layout Optimization
Considering the Horns Rev 1 WF, whose layout is shown in Figure 8, we can formulate the layout optimization problem as described in Equation (1), i.e., maximizing the expected total power output , while satisfying the constraints of WF boundary and minimal distance between any two WTs (set as five rotor diameters). The problem can be solved by using the random search algorithm, which was developed in our previous study [12]. Since this study is focused on wind modelling, the details of the wake modelling and the optimization algorithm are referred to Ref. [12]. In this section, the important task of choosing appropriate bin sizes and for wind modelling in numerical calculation is investigated, both for power calculation and layout optimization.

Assessment of Bin Sizes for Power Calculation
Since the expected power output of a WF is calculated by numerical integration over all interested wind conditions, the result of calculation can depend on the choice of bin sizes. Thus, we first make a sensitivity study of the power calculation for the Horns Rev 1 WF regarding different combinations of and . The results are shown in Table 2 for using the three proposed joint distributions. From the results obtained by using each distribution, and comparing the variations of the calculated along rows and along columns, we can see that power calculation is relatively more sensitive to the choice of than that of , which is consistent with what we have shown through the ideal WF in Section 2 and also in our previous study [12]. We can also see that the three joint distributions proposed in the last section give similar results when using the same combination of and , which suggests that they can all be applied for wind modelling in power calculation or layout optimization. In the rest of this study, the spline joint distribution will be used, since it fits the measurement data best.
Note that the expected power output, when calculated using the 3 years' time series directly, is 78.79 MW , while the value obtained by using joint distributions with 1 m/s and 30° is 76.86 MW. Comparing these two values, we can see the relative difference is within 2.45%, and we may conclude that the common practice of choosing 1 m/s and 30° is appropriate, if the task is to assess the power production of a given large WF with a layout of regular shape, such as that of Horns Rev 1 WF. Besides, we can see that using 1 m/s and 10° is a better choice, since it obtains a result closer with the value from using measurement data directly, but without increasing the computation cost dramatically.
However, as we have demonstrated through the ideal WF in Section 2, the choice of appropriate bin sizes, especially bin size of wind direction, has to be made through more careful investigations, if the task is to carry out layout optimization.

Choice of Bin Size for Layout Optimization
The effect of bin size for layout optimization is investigated here, by solving the layout optimization problem with the RS algorithm using a fixed value of and four different values of . The spline joint distribution is applied in wind modelling and the optimized layouts using these four bin size choices, i.e., bs-1: 1 m/s, 30°; bs-2: 1 m/s, 5°; bs-3: 1 m/s, 3°; bs-4: 1 m/s, 1°, are shown in Figure 9. Note that 'bs-i' denotes the ith kind of bin size choice. It is obvious that the optimized layout in Figure 9a shows the largest deviation from the original layout. In order to better assess the superiority of the optimized layouts over the original layout, the power improvements in the four cases, obtained with bin size choices: bs-1, bs-2, bs-3 and bs-4, are re-evaluated with different combinations of bin sizes of wind speed and wind direction, and also with measured time-series directly. The results are shown in Figure 10. It can be seen that, when assessing the superiority of an optimized layout over the original layout, the calculated power improvement is more sensitive to in re-evaluation than to , and the re-evaluated improvement with 1 m/s and 1° is more consistent with that obtained by using measured time-series directly.
In this figure, the optimized layout with bin size choice bs-1, i.e., 1 m/s , 30°, shows an impressive improvement of up to 6.5% when evaluated with the same choice of bin sizes, but the improvement turns into actually a negative value when evaluated with smaller or measured time-series. Comparing the four optimized layouts, we can see that the optimized layout with bin size choice bs-4, i.e., 1 m/s, 1°, obtains the largest improvement when re-evaluated with smaller ( 1°) or measured time-series, and it also shows the most consistent improvement when evaluated with different combinations of bin sizes.
From these observations, we can see that the choice of bin size is of crucial importance for layout optimization. In order to obtain reliable and consistent results, using 1° in the optimization process is necessary, and the common practice of using 30° tends to obtain artificial improvement, which may looks impressive but actually leads to decreased power. This is consistent with our previous finding [12].

Choice of Bin Size for Layout Optimization
Similarly, the layout optimization is carried out here by using a fixed value of and four different values of . The optimized layouts, obtained using these four combinations of bin sizes, i.e., bs-5: 2 m/s, 1°; bs-6: 1 m/s, 1°; bs-7: 0.5 m/s, 1°; bs-8: 0.1 m/s, 1°, are shown in Figure 11. The four optimized layouts, obtained with bin size choices: bs-5, bs-6, bs-7 and bs-8, are then assessed by using the same re-evaluation procedure as in Section 5.2, and the obtained results are shown in Figure 12.
It can be seen that the optimized layouts with different choices of bin sizes bs-5, bs-6, bs-7 and bs-8, i.e., 1° and different , show nearly the same level of improvements when re-evaluated with smaller ( 1°) or measured time-series. From this figure, we can see that the choice of bin size is not so important when 1° is chosen. Thus, the common practice of using 1 m/s is valid for application in layout optimization. Based on the above investigations, we can conclude:  The proposed continuous joint distributions of wind speed and wind direction can be applied in both wind farm power calculation and layout optimization;  The common practice of using 1 m/s and 30° might be appropriate to assess the power production of a given wind farm with a layout of regular shape, but it is not suitable to be applied in layout optimization;  The choice of using 1 m/s and 1° is recommended in layout optimization, in order to obtain reliable and consistent optimization results.

Conclusions
In this study, a simple and easily implementable method for constructing joint distributions of wind speed and wind direction is proposed. First, the measurement data is fitted into sector-wise Weibull distributions. The obtained parameters are used to construct three types of joint distributions, namely, piecewise joint distribution, linear joint distribution and spline joint distribution. For the three years' wind measurement data at Horns Rev, reasonable levels of good fitness have been obtained by the three proposed distributions, among which the spline joint distribution shows the best fitness.
Second, the choice of bin sizes for wind modelling in numerical calculation is investigated. The importance of choosing small enough value for bin size of wind direction is addressed by considering a constructed ideal wind farm. Furthermore, the wind modelling problem is investigated in a realistic setting, i.e., for the layout optimization problem of Horns Rev 1 wind farm. The proposed spline joint distribution is applied for wind modelling, and the problem of choosing appropriate bin sizes for wind speed and wind direction is solved by carrying out a detailed sensitivity study. It has been found that the choice of bin size for wind direction is crucial for layout optimization and the choice of using 1 m/s and 1° is recommended, so that reliable and consistent optimization results can be obtained.
Future study will consider the uncertainties in wind modelling for wind farm layout optimization and their impacts on optimization results. The uncertainties may come from various sources, such as the wind measurement, the statistical models, the representative limitations (considering the fact that we are using the measured data in several years to predict the wind condition in the future life time of the wind farm, which is typically up to 20 years).