Comparison of Optimum Spline-Based Probability Density Functions to Parametric Distributions for the Wind Speed Data in Terms of Annual Energy Production

The common approach to wind energy feasibility studies is to use Weibull distribution for wind speed data to estimate the annual energy production (AEP). However, if the wind speed data has more than one mode in the probability density, the conventional distributions including Weibull fail to fit the wind speed data. This highly affects the technical and economic assessment of a wind energy project by causing crucial errors. This paper presents a novel way to define the probability density for wind speed data using splines. The splines are determined as a solution of constrained optimization problems. The constraints are the characteristics of probability density functions. The proposed method is implemented for different wind speed distributions including multimodal data and compared with Weibull, Weibull and Weibull and Beta Exponentiated Power Lindley (BEPL) distributions. It is also compared with two other nonparametric distributions. The results show that the spline-based probability density functions produce a minimum fitting error for all the analyzed cases. The AEP calculated based on this method is considered to have high fidelity, which will decrease the investment risk.


Introduction
Population and industry increase in the world which means that the demand for energy increases too. Fossil-fuel based energy sources have negative effects on the environment and they are running out. The greenhouse gases, such as carbon dioxide (CO 2 ), Methane (CH 4 ) and Nitrous Oxide (N 2 O), are the main factors of increasing global warming [1,2]. Those gases arise from burning fossil fuels. Alternative renewable sources must be developed to decrease the emission of greenhouse gases. Wind energy is one of the renewable energy sources which is rapidly growing in most of the countries due to its relatively low cost [3].
Before the installation of a wind farm, detailed analyses of the wind energy potential of the site including the estimation of the annual energy production (AEP), capacity factor calculation and cost analysis should be performed to minimize the investment risk [4]. The available wind power in a certain site is evaluated using wind speed data. Finding accurate and reliable methods to represent the wind speed data is a crucial issue in wind energy assessment. Most of the methods used in the literature are based on the probability distribution functions (PDF). Due to its simplicity, Weibull function [5] is the most commonly used method to fit the wind speed data [6][7][8][9][10][11][12][13][14][15]. The Weibull based

Spline-Based Probability Density Function
The wind speed data is originated from meteorological observations. It provides the frequency of the observed wind speed along a certain time, which is usually one year. A typical wind speed data is shown in Figure 1 as a histogram plot. This frequency data of the wind speed gives natural information about the probability density distribution of the wind speed and shown by OD w j , j = 1, 2, . . . , M where M is the number of the distinct values of the observed wind speed. A common approach is to fit this data using a parametric probability density function such as Weibull or Rayleigh distributions.

Spline-Based Probability Density Function
The wind speed data is originated from meteorological observations. It provides the frequency of the observed wind speed along a certain time, which is usually one year. A typical wind speed data is shown in Figure 1 as a histogram plot. This frequency data of the wind speed gives natural information about the probability density distribution of the wind speed and shown by ( ), = 1,2, … , where is the number of the distinct values of the observed wind speed. A common approach is to fit this data using a parametric probability density function such as Weibull or Rayleigh distributions. In this study, the wind speed data is fitted using splines to construct piecewise polynomial approximations on a number of intervals out of the observed wind speed range. In statistics, a probability density function obtained using such an approach is known as nonparametric distribution [28].
The procedure of implementing the proposed method is as follow. The range of the wind speed from 0 m/s ( ) to the maximum observed wind speed ( ) is first evenly divided into intervals that lead to polynoms which pass through + 1 nodes. It should be noted that is an input by the user and its value is less than , which means that the spline to be fitted is not required to pass through all the observed wind speed data. This condition is needed to satisfy some constraints on the probability density function such as providing the unity probability for the whole range. Also, the fitted probability density function will be smooth similar to the conventional distributions. In the present paper, cubic polynomials are used for constructing the spline. The piecewise cubic polynomial ( ), between each successive pair of nodes has the following form: ( ) is defined between nodes − 1 and . The spline, ( ), is the union of all ( ) and is continuous through each node. ( ) is zero when < and > .The unknown values of the function, ( ), at each node and its first derivative values in the first and last nodes are determined together with the spline coefficients, , , and by solving the constrained optimization problem given below: In this study, the wind speed data is fitted using splines to construct piecewise polynomial approximations on a number of intervals out of the observed wind speed range. In statistics, a probability density function obtained using such an approach is known as nonparametric distribution [28].
The procedure of implementing the proposed method is as follow. The range of the wind speed from 0 m/s (w 0 ) to the maximum observed wind speed (w M ) is first evenly divided into N intervals that lead to N polynoms which pass through N + 1 nodes. It should be noted that N is an input by the user and its value is less than M, which means that the spline to be fitted is not required to pass through all the observed wind speed data. This condition is needed to satisfy some constraints on the probability density function such as providing the unity probability for the whole range. Also, the fitted probability density function will be smooth similar to the conventional distributions. In the present paper, cubic polynomials are used for constructing the spline. The piecewise cubic polynomial f i (w), between each successive pair of nodes has the following form: f i (w) is defined between nodes i − 1 and i. The spline, f (w), is the union of all f i (w) and is C 2 continuous through each node. f (w) is zero when w < w 0 and w > w N .The unknown values of the function, f i (w), at each node and its first derivative values in the first and last nodes are determined together with the spline coefficients, a i , b i , c i and d i by solving the constrained optimization problem given below: where, w is the mean wind speed of the observed data while w 2 is the summation of variance of the observed wind speed and square of the mean wind speed.
The optimization variables are f (w i ), (i = 0, 1, . . . , N), f (w 0 ) and f (w M ). Prime sign ( ) denotes the first derivative with respect to the wind speed, w. Therefore, the total number of the unknowns is N + 3. Note that, a i , b i , c i and d i , (i = 1, 2, . . . , N) are not optimization variables since they are determined while constructing the spline f (w) based on the node values f (w i ). For example, we have 9 variables in the 7-node optimum spline, 7 variables in 5-node optimum spline and 5 variables in the 3-node optimum spline. However, since there are 5 constraints, three of which are equality constraints, the effective number of variables in an optimum spline decreases at least by 3. Therefore, it is quite fair to compare a 7-node optimum spline with Weibull and Weibull and BEPL methods which have 5 parameters each.
RMSE is the root mean square error which is a measure of the residuals between the fitted distribution and the observed distribution. Constraint 1 guarantees that the cumulative probability over the observed wind speed range is unity because the first and the third terms in the right hand side of Equation 3 vanish since f (w) is defined as zero beyond the observed wind speed range: Constraint 2 is due to the nature of the probability density functions which are always positive in value. Constraint 3 is required to avoid statistically possible but physically meaningless solutions. f max is selected to be 1.0 because there is no site (to be considered for wind energy project) that has a wind speed PDF which is greater than 1.0. It should be noted that inequalities (Constraints 2 and 3) in the optimization problem are given in a discrete way while solving the problem numerically.
In Constraints 4 and 5, the first and second moments of the probability density function are used. This is proposed in order to force the same mean and standard deviation values as the observed data. This approach is known as the method of moments in fitting the conventional distributions. However, the conventional distributions obtained using methods of moments do not minimize the RMSE and might result in large errors. On the other hand, the proposed method in this study does not only satisfy this empirical condition but also provides minimum RMSE. Therefore, it is considered that the fitting spline will provide the optimum probability density function for the related wind speed data.
In order to solve the problem given in Equation 2, the constrained optimization is transformed into the unconstrained optimization by adding the inequality constraints as exterior penalty functions and the equality constraints with Lagrange multipliers to the objective function [29]. The final form of the objective problem is given in Equation 3. It is solved using a gradient-based optimization algorithm, the conjugate gradient method [30].
The penalty multipliers, r p2 and r p3 , in Equation (4) are chosen as a small value at the beginning of the optimization process. Along the optimization iterations, the penalty multipliers are increased by a factor. In this study the initial value of the multipliers is selected as 1, and the factor as 3. It should be noted that the number of the optimization variables is increased in this unconstrained form. The new added optimization variables are the Lagrange multipliers λ 1 , λ 4 and λ 5 [29].
At each optimization iteration, l, the optimization variables are updated as follow [30]: where g 1 , g 2 , g 3 , g 4 , g 5 and g 6 are the components of the conjugate gradient direction for the six optimization variables f (w i ), f (w 0 ), f (w M ), λ 1 , λ 4 and λ 5 , respectively. β is the step size along the conjugate gradient direction and is calculated using a simple line search based on trial and error method.

Weibull Distribution Based Probability Density Function
The Weibull distribution is a two-parametric probability density function. These two parameters are the shape factor, k, and the scale factor, c, as shown in Equation (6) [31]: There are different ways to estimate the Weibull parameters, k and c. The method used in this study is based on the minimization of the root mean square error, RMSE.

Weibull and Weibull Distribution Based Probability Density Function
The Weibull and Weibull (W and W) distribution is defined in the convex hull of two Weibull distributions containing 5 parameters as given in Equation (7) as [17]: These five parameters, p, k 1 , k 2 , c 1 and c 2 are estimated based on the minimization of RMSE.

Beta Exponentiated Power Lindley Distribution Based Probability Density Function
The recently proposed BEPL distribution is defined in Equation (8).
where, B(a, b) is the Beta function. The five parameters (k, c, ω, a and b) used in the distribution are again estimated by minimizing the RMSE.

Empirical Probability Density Function
Assigning the normalized frequency as the probability for each observation is known as the empirical probability density function. The normalization has been done by the total number of observations throughout the recorded time period. It should be noted that the empirical method is a nonparametric distribution.

Kernel Probability Density Function Based on 3rd Order Polynomial
A kernel distribution produces a nonparametric probability density fit which adapts itself to the data. In this study, we defined this distribution by a kernel density estimator as a 3rd order polynomial. Moreover, we implemented two more conditions while defining the distribution: first and second moments of the distribution are set same as the observation data.

Annual Energy Production
The annual energy production calculations are very vital in the evaluation of any wind energy project. The wind speed distribution is combined with the power curve of a specific wind turbine to give the energy generated at each wind speed and hence the expected total energy generated overall the year. The annual energy production (AEP) can be expressed mathematically as follow: (9) where, P(w), is the power output of a wind turbine at wind speed, w, and, 8760, is the number of hours in the year. f (w) is the distribution function which has a closed form.

Analyzed Cases
The wind speed data cases analyzed in this study are shown in Figure 2. Those cases are selected from Yürüşen and Melero [20] as illustrative examples to show the effectiveness of the proposed method over the parametric distributions. It should be noted that those data are synthetic wind speed data. Visual inspection of the histograms hints that Case 1, 2 and Case 4 correspond to multimodal distributions. As for Case 3, the distribution is a simple unimodal one. The statistical properties of the observation data are summarized in Table 1. It should be noted that the mean and the standard deviation values of the distributions are calculated within the range of the observation data.

Results
Three different optimum splines are fitted for the cases analyzed. The splines are obtained using 3, 5 and 7 nodes. The mean speed, standard deviation and the RMS error of the different fitted distributions for all the analyzed cases are summarized in Table 2.
The mean speeds and standard deviations of the optimum splines are the same as the observation data thanks to the constraints specified in the optimization problem given in Section 2.1. As seen from Table 2, the RMS errors of Weibull distribution are much higher than the errors of the optimum splines in Cases 1, 2 and 4 (around 4 times more). In Case 3, the RMS error of the Weibull distribution is about twice of the errors of the optimum spline method. Although the number of parameters in the BEPL distribution is 5 and therefore a good performance was expected, the results for all the analyzed cases were observed to be poor in terms of the RMS error. As for the Weibull and Weibull distribution, the RMS errors for all cases were as low as the error of the distribution based on the 5-node optimum spline. However, the distribution based on the 7-node optimum spline is superior to the Weibull and Weibull distribution as seen from the RMS error results in the table. The kernel distribution shows almost the same RMS error performance as the 3-node optimum spline for all the cases.

Results
Three different optimum splines are fitted for the cases analyzed. The splines are obtained using 3, 5 and 7 nodes. The mean speed, standard deviation and the RMS error of the different fitted distributions for all the analyzed cases are summarized in Table 2.
The mean speeds and standard deviations of the optimum splines are the same as the observation data thanks to the constraints specified in the optimization problem given in Section 2.1. As seen from Table 2, the RMS errors of Weibull distribution are much higher than the errors of the optimum splines in Cases 1, 2 and 4 (around 4 times more). In Case 3, the RMS error of the Weibull distribution is about twice of the errors of the optimum spline method. Although the number of parameters in the BEPL distribution is 5 and therefore a good performance was expected, the results for all the analyzed cases were observed to be poor in terms of the RMS error. As for the Weibull and Weibull distribution, the RMS errors for all cases were as low as the error of the distribution based on the 5-node optimum spline. However, the distribution based on the 7-node optimum spline is superior to the Weibull and Weibull distribution as seen from the RMS error results in the table. The kernel distribution shows almost the same RMS error performance as the 3-node optimum spline for all the cases. A comparison of the 7-node optimum spline based distribution fit with the Weibull distribution for the analyzed cases is given in Figure 3. One may notice that the optimum spline based distributions fits very well against the observation data for all the analyzed cases. On the other hand, Weibull based distribution fits well only for Case 3 in which the observation data has a simple unimodal distribution. However, for the other more complicated cases (multimodal distributions), the distributions given by the Weibull method significantly deviate from the observation data. This will obviously affect the accuracy of the expected value of the annual energy production.
Fitting the observation data using the Weibull and Weibull method is compared to the 7-node optimum spline fit in Figure 4. In Cases 2 and 4, the distributions based on the optimum spline and the Weibull and Weibull method are observed to be very close to each other. An eye inspection of Case 1 states that there are more than two modes. Therefore, the proposed method based on the 7-node optimum spline captures the observed data better than the Weibull and Weibull method which is suitable for dual modes. Since there is only one mode in Case 3, there is no superiority of the Weibull and Weibull distribution over the original Weibull distribution as expected. It is also noticed that the optimum spline captures the peak point more accurately for this case. Fitting the observation data using the Weibull and Weibull method is compared to the 7-node optimum spline fit in Figure 4. In Cases 2 and 4, the distributions based on the optimum spline and the Weibull and Weibull method are observed to be very close to each other. An eye inspection of Case 1 states that there are more than two modes. Therefore, the proposed method based on the 7-node optimum spline captures the observed data better than the Weibull and Weibull method which is suitable for dual modes. Since there is only one mode in Case 3, there is no superiority of the Weibull and Weibull distribution over the original Weibull distribution as expected. It is also noticed that the optimum spline captures the peak point more accurately for this case.   Fitting the observation data using the Weibull and Weibull method is compared to the 7-node optimum spline fit in Figure 4. In Cases 2 and 4, the distributions based on the optimum spline and the Weibull and Weibull method are observed to be very close to each other. An eye inspection of Case 1 states that there are more than two modes. Therefore, the proposed method based on the 7-node optimum spline captures the observed data better than the Weibull and Weibull method which is suitable for dual modes. Since there is only one mode in Case 3, there is no superiority of the Weibull and Weibull distribution over the original Weibull distribution as expected. It is also noticed that the optimum spline captures the peak point more accurately for this case.   The BEPL fitting is plotted and compared in Figure 5. As seen from the figure, the BEPL distribution does not show a satisfactory performance in capturing the observation data. The fitting results of the optimum splines with different node numbers are compared in Figure  6. It is noticed that the optimum splines based on 5 nodes and 7 nodes are very close to each other and to the observation data. This result was already seen in Table 2 with very small RMS error The BEPL fitting is plotted and compared in Figure 5. As seen from the figure, the BEPL distribution does not show a satisfactory performance in capturing the observation data.  The BEPL fitting is plotted and compared in Figure 5. As seen from the figure, the BEPL distribution does not show a satisfactory performance in capturing the observation data. The fitting results of the optimum splines with different node numbers are compared in Figure  6. It is noticed that the optimum splines based on 5 nodes and 7 nodes are very close to each other and to the observation data. This result was already seen in Table 2 with very small RMS error The fitting results of the optimum splines with different node numbers are compared in Figure 6. It is noticed that the optimum splines based on 5 nodes and 7 nodes are very close to each other and Energies 2018, 11, 3190 11 of 15 to the observation data. This result was already seen in Table 2 with very small RMS error values. The optimum spline based on 3 nodes fit only well for Case 3. However, it is still better than the Weibull distribution in terms of RMS error values as seen in Table 2.
Energies 2018, 11, x FOR PEER REVIEW 11 of 15 values. The optimum spline based on 3 nodes fit only well for Case 3. However, it is still better than the Weibull distribution in terms of RMS error values as seen in Table 2. In this study, the selected wind turbine is Enercon E53-800 kW. This turbine, like the majority of large size wind turbines, is pitch regulated. This leads to a fixed power output at speeds higher than the rated wind speed. The technical specifications and the power curve of this wind turbine are given in Table 3 and Figure 7, respectively.  In this study, the selected wind turbine is Enercon E53-800 kW. This turbine, like the majority of large size wind turbines, is pitch regulated. This leads to a fixed power output at speeds higher than the rated wind speed. The technical specifications and the power curve of this wind turbine are given in Table 3 and Figure 7, respectively. In this study, the selected wind turbine is Enercon E53-800 kW. This turbine, like the majority of large size wind turbines, is pitch regulated. This leads to a fixed power output at speeds higher than the rated wind speed. The technical specifications and the power curve of this wind turbine are given in Table 3 and Figure 7, respectively.  Using Equation (9), the annual energy production is calculated in all the cases for spline and parametric methods. The results are shown in Table 4. The baseline for comparison with all the methods is selected as the empirical probability density function, that is, the observation data itself. The AEP results based on 7-node optimum spline are almost the same as the baseline with a difference of less than 0.4%. 5-node optimum spline gives also similar AEP values. The Weibull and Weibull method gives AEP estimations with less than 1% deviation. The Weibull distribution resulted in about 5 to 8% differences from the baseline AEP for Cases 2 and 4. The AEP results based on the BEPL method differ within the range of 4 to 22%. The 4% difference belongs to Case 3, which is a unimodal case. The kernel distribution deviates between 1 to 2% from the baseline AEP.
Inaccurate estimations in the AEP will obviously mislead the assessment project and highly increase the investment risk. For example, a difference of 4% leads to about 0.2 MWh in AEP which should not be neglected.

Conclusions
In this work, a new method is proposed for the probability distribution function of the observed wind speed data. The method is based on fitting splines for a given wind speed versus probability data. The splines obtained are optimized for minimum root mean square error subjected to a number of constraints such as satisfying given mean and standard deviations. The method is compared to the Weibull, Weibull and Weibull and BEPL methods. Among those three distributions, the Weibull method is the most widely used in wind assessment applications. The method is also compared with two nonparametric probability distribution functions: empirical and 3rd order polynomial based kernel distributions.
Four different wind speed distributions, three of which involve multi-modality, are analyzed. The results show that the Weibull and BEPL methods fail in fitting the multi-modal distributions with high RMS errors. The disagreement in fitting causes miscalculation in annual energy production. On the other hand, the optimum spline and the Weibull and Weibull methods produce very small RMS errors with a slight superiority of the distribution based on the 7-node optimum spline. Accurate prediction of AEP is vital in wind energy project assessment to minimize the investment risk. Since the investment cost of wind energy projects is relatively high, even slight prediction error in AEP calculations leads to a great difference in the payback period of the project. Therefore, one may conclude that the AEP calculated based on the optimum splines distributions are considered to be more accurate and reliable.
However, since obtaining the optimum splines requires the solution of a constrained optimization problem, a computation involving a lot of mathematical operations is necessary. This method is recommended when multimodality is present in the wind speed data.