Five different distributions and metaheuristics to model wind speed distribution

This paper presents a comprehensive empirical study of five distribution functions to analyze wind energy potential: Rayleigh, Weibull, Gamma, Burr Type XII, and Generalized Extreme Value. In addition, two metaheuristics optimization methods, Grey Wolf optimization and Whale optimization algorithm, are utilized to determine the optimal parameter values of each distribution. Five error measures are investigated and compared to test the accuracy of the introduced distributions and optimization methods, such as mean absolute error, root mean square error, regression coefficient, correlation coefficient, and net fitness. The Catalca site in Istanbul, Turkey, was selected to be the case study to conduct this analysis. The obtained results confirm that all introduced distributions based on optimization methods efficiently model wind speed distribution in the selected site. Although Gamma distribution based on GWO and WOA outperformed other distributions for all datasets at all heights, it was the worst in terms of computation complexity. Rayleigh distribution occupied the latest rank, but it was the best in terms of computation complexity. MATLAB 2020b and Excel 365 were used to perform this study. Cite this article as: Mohammed W. Five different distributions and metaheuristics to model wind speed distribution. J Ther Eng 2021;7,Supp14:1898–1920. Journal of Thermal Engineering Web page info: https://jten.yildiz.edu.tr DOI: 10.18186/thermal.1051262


INTRODUCTION
Wind energy is one of the most rapidly growing renewable energies worldwide due to its many advantages. Wind speed distribution for a specific site shows the available wind potential. Once the distribution of wind patterns is established, the wind energy potential can be determined. Wind energy availability allows investors and developers to collect more accurate feasibility of the underlying site. The feasibility study process consists of several steps: however, wind regime characterization via distribution functions is one of the most critical steps. During the last decades, many Probability Density Functions (PDFs) have been appeared in literature to represent wind speed distribution, such as Weibull [1], Rayleigh [2], Gamma [3], Normal [4], Lognormal [5], Logistical [6], Beta [7], Nakagami [8], Burr [9] distribution functions, and others. However, Rayleigh, Gamma, and Weibull functions are the most widely used among these functions [10], [11], [12].
Wadi et al. [12] presented a statistical study based on the two-parameter Weibull distribution function for wind data for three years at Catalca in Istanbul, Turkey. Three estimation methods for Weibull parameters: approximation, graphical, and Energy Pattern Factor (EPF) methods, are studied and compared. Wind data for two years from the southern area in Pakistan at four heights were used to evaluate the wind potential based on Weibull distribution [13]. The Weibull parameters were modeled based on three different algorithms: Grey Wolf Optimization (GWO), Particle Swarm Optimization (PSO), and Cuckoo Search Optimization (CSO). Moreover, four numerical methods; EPF, Method of Moments (MOM), Empirical Method of Justus (EMJ), and Modified Maximum Likelihood (MML) estimation methods were used. Junk and Schindler [14] presented 24 distributions to assess their Goodness-Of-Fit (GOF) based on four years of wind speed data at different sites worldwide. Parameters of distributions were varied between one and five parameters. MOM, L-Moment (LM), Maximum Likelihood (ML), and least-squares (LS) methods are used to assess these parameters.
Two-parameter Weibull distribution function based on four estimation methods: graphical, empirical, EPF, and MML to estimate the capacity factor of wind turbines in Jaisalmer district of western Rajasthan in India was presented in [15]. MML method provided the best matching while the graphical method had the least matching. Twoparameter Weibull, three-parameter Weibull, two-parameter Gamma, and two-parameter Lognormal were presented to model wind speed at the airport site in Dolny Hricov [16]. ML method was applied to estimate the parameters of the distributions. The three-parameter Weibull and the two-parameter Weibull achieved the best first and second fitness, respectively.
Due to the inability of the Weibull distribution function to achieve the required matching in some wind regimes and the calculations complexity of mixture distribution functions, many researchers introduced different distribution functions like Birnbaum Saunders [17], [18]. Mohammadi et al. [19] presented the Birnbaum Saunders distribution to model the wind speed frequency distribution for ten sites distributed in the Ontario province of Canada. The GOF of the two-parameter Birnbaum Saunders distribution was compared with earlier nine one-component distributions, and the obtained results have confirmed the best fitness of the Birnbaum Saunders distribution. Several types of research based on Burr and inverse Burr distribution functions have appeared in the literature [9], [20], [21], [22]. A statistical analysis of wind speed data based on a four-parameter Burr distribution function in Antakya, Turkey, was proposed in [9]. GOF proved that the Burr distribution was more convenient than the Weibull and generalized-Gamma distributions for describing the real distribution of wind speed data. Chiodo and Falco in [21] presented the two-parameter inverse Burr distribution function to represent the extreme values of wind speed. Three estimation methods: moment, ML, and quantile methods are used. The inverse Burr distribution was provided the best matching.
Wind speed and frequency are basically based on the site. Hence, some distribution functions that can successfully describe wind patterns at a specific site could fail at another. Therefore, this paper presents an extensive study to compare the performance of five distributions, namely, Rayleigh, Gamma, Weibull, Burr XII, and Generalized Extreme Value. Two metaheuristics optimization methods: GWO and Whale Optimization Algorithm (WOA) are utilized to evaluate the optimal parameter values of each distribution. As various GOF measures can lead to varied fit evaluation results, five measures such like Mean Absolute Error (MAE), Root Mean Square Error (RMSE), Regression Coefficient (R 2 ), Correlation Coefficient (R), and net fitness are applied to assess the accuracy of introduced distributions. Wind data for three years from 2008 to 2010 at three different heights, i.e., 30, 60, and 80 m, from the Marmara area (Catalca) are utilized to carry out this study.
The rest of this paper is organized as follows: Section II is dedicated to exploring the statistical distribution functions, including the formulas of PDF, Cumulative Density Function (CDF), and Inverse CDF (ICDF) for each distribution. Furthermore, Section III introduces the methodology of selecting the optimal parameters of the distributions using GWO and WOA. Section IV elaborates the wind speed extrapolation procedure and wind energy potential. Section V explains the accuracy measures that are exploited to examine the performance of each distribution. Besides, Section VI discusses the obtained results. Finally, Section VII concludes the paper.

STATISTICAL DISTRIBUTION
As said, determination of the wind speed distribution is an essential step in assessing wind energy potential at a particular location. Once the wind speed distribution pattern is accurately determined, the technical and economic features belonging to the location can be appropriately specified. To describe the distributions of wind speed data, five different distributions have been studied. Hereinafter, a concise explanation of these distributions is given.

RAYLEIGH DISTRIBUTION
Rayleigh distribution (RD) was derived by British physicist Lord Rayleigh. The one-parameter Rayleigh distribution is a particular case of Weibull distribution in which the shape parameter (k W ) is assumed to equal 2 [23]. Due to its simplicity and the ability to accurately describe wind regimes, many researches used Rayleigh distribution to evaluate the potential of wind at various locations worldwide [10], [24], [25]. One-parameter Rayleigh PDF is defined as follows [26]. (1) where f(v) is the probability of observing wind speed v, v = 0, 1, 2, . . ., N, N is the size of wind speed vector, and b R > 0. The Rayleigh CDF, F(v) and ICDF, G(p) are defined according to equation (2) and equation (3), respectively.

WEIBULL DISTRIBUTION
Weibull distribution (WD) is one of the commonly used functions to represent wind speed measurements due to its high flexibility and accuracy [2], [15], [25]. Bi-parameter Weibull distribution depends on the shape (k W ) and the scale (c W ) parameters. The PDF of Weibull distribution is given by: where k W and c W are the shape and scale parameters. Weibull distribution has several vital characteristics that its parameters are evaluated at a particular height, and it is conceivable to extrapolate them for different heights [27]. The Weibull CDF, F(v) and ICDF, G(p) [28], [29] are defined as equations (5) and (6), respectively

Gamma Distribution
Gamma distribution (GD) is also a widely used distribution; it is largely related to Exponential and Normal distributions. Besides, many distributions like Exponential, Chi-squared, and Erlang distributions are special cases of the Gamma distribution. The Gamma PDF is expressed as [10]: where k G and c G > 0, and Γ is the gamma function and can be computed by the following formula: The Gamma distribution formulas for the CDF, F(v) and ICDF, G(p) [10] are according to equations (9) and (11), respectively.
where Γ v (γ) can be computed using the following formula.

Burr XII Distribution
Burr distribution (BD), also known as Burr Type XII distribution, is a continuous distribution which Irving W. Burr originally introduces. The PDF of the three-parameter Burr distribution is given by the following equation [30], [31], [32]: The Burr distribution CDF, F(v) [33] and ICDF, G(p) [34] can be computed from equations (13) and (14), respectively.
Generalized Extreme Value Distribution Generalized Extreme Value distribution (GEVD) is a continuous probability function developed by combining three simpler distributions: Extreme Value, Frechet, and Weibull. Generalized Extreme Value distribution is often used as an approximation to model the maxima of long (finite) sequences of random variables. The following formula gives the PDF of the three-parameter Generalized Extreme Value distribution [35].
where k GE , c GE , and λ GE are the shape, scale, and location parameters, respectively. Furthermore, Generalized Extreme Value distribution is pivotal in modeling the extreme wind events [36]. Generalized Extreme Value distribution formulas for the CDF, F(v) and ICDF, G(p) are as follows [35]:  Table 1 summarizes all the used distributions along with the name and notation of their parameters.

MATERIALS AND METHODOLOGY
Evolutionary Algorithms (EAs) are prevalent in many research areas. This is because such algorithms are simple, flexible, and capable of avoiding local optima. On the other hand, EAs may suffer from various drawbacks such as long computation time, no guarantee to converge, and having several operating parameters to be adjusted before starting [37]. In order to overcome the limitations mentioned above, GWO and WOA are utilized. In this section, GWO and WOA are introduced briefly. Then, the methodology of exploiting the metaheuristics to select the optimal parameters of the introduced distributions is explained.

Grey Wolf Optimization
GWO algorithm imitates the leadership hierarchy and hunting mechanism of grey wolves in nature. Four types of grey wolves, ordered from the top of the pack hierarchy as alpha, beta, delta, and omega, are used for simulating this mechanism. Besides, three fundamental steps of hunting, namely, searching for prey, encircling prey, and attacking prey, are implemented to perform optimization [37].
The process of GWO can be summarized as follows. Firstly, all grey wolves (search agents) in the pack (population) are initialized randomly in [LB, UB], where LB and UB are the lower and upper bounds of the problem variables, respectively. Afterward, the fitness score for each search agent is evaluated using the accompanying objective function. The fittest solution so far is alpha, followed by beta and delta, respectively. Meanwhile, the rest of the wolves are grouped under the omega. The optimization in GWO is carried out by alpha, beta, and delta, whereas omega follows them. Then, search agents update their position based on prey position using the following equations: X t X X X + where t, X(t), X (t + 1), X α , X β , X δ denote the current iteration, the current position of the prey, the position of the prey in the next iteration, the position of alpha, the position of beta, the position of the delta, respectively. The coefficient vectors A and C are computed according to the following equations.
where r 1 and r 2 are random vectors uniformly distributed in [0, 1], and the control vector a is linearly decreased from 2 to 0 for better exploration and exploitation of candidate solutions. GWO continues in the same procedure until the maximum number of iterations is reached. Figure 1 (a) shows the flowchart of GWO.

Whale Optimization Algorithm
WOA is a new metaheuristic for optimization problems proposed in 2016 [38]. It has been inspired by the unique hunting behavior of humpback whales. In WOA, a population of whales (search agents) evolves to find the global optima after a specified number of iterations. WOA begins with the initialization of search agents randomly upon the interval of LB and UB of the problem variables. After that, WOA evaluates the fitness score for each search agent by using the fitness function. The best solution is saved for further processing later.
Then, WOA updates the position of each search agent depending on the following cases. If a random number (z) and the A vector are less than 1, then the particular search agent applies the Encircling method by updating its position for the next iteration using the following formulas: where X*(t) is the position of the best solution so far. Furthermore, A and C are coefficient vectors. The a vector is linearly decreased from 1 to 0. The r vector is a random vector in [0, 1]. Else if a random number (z) is less than 0.5, but A vector is bigger than 1, then the particular search agent applies the Exploration method by updating its position for the next iteration using the following formulas: where X rand is the random whale in the current iteration, otherwise, the particular search agent applies the Spiral method by updating its position for the next iteration using the following equation: .exp .cos *, . π (29) where Dʹ = |X* -X(t)| represents the distance between the search agent to the prey, b is a constant, and l is a random number [-1,1]. Finally, WOA continues in the same process until the maximum number of iterations is reached. Figure  1 (b) depicts the flowchart of WOA.

Methodology
The parameters selection of a distribution may be defined as a non-linear optimization task that minimizes the mean absolute error between the collected and predicted wind speed vectors. Mathematically as shown in Equation (30).
where Vm is the collected wind speed vector, and Vd is the predicted wind speed vector. Furthermore, Vd can be calculated based on the ICDF of the distribution, as expressed in Section II.
To solve the optimization as mentioned earlier problem, EAs can be used. Many algorithms were appeared in the literature to solve non-linear optimization problems and applied in many engineering fields [39], [40], [41], [42], [43]. In GWO and WOA, a group of search agents designates candidate solutions to the optimization problem. Concerning the parameter selection problem of the distribution, every search agent consists of integer values representing the values of the distribution parameters. The first group of search agents is generated randomly within specified boundaries of the parameters.
Then, within the initial iteration, the fitness score values per search agent are calculated via Equation (30). Next, the population increases by seeking the optimum solution applying the specified operation of the used EA. It continues in the same manner until the condition of the maximum number of iterations is satisfied. Table 2 gives the simulation parameters of GWO and WOA decided after performing a comprehensive exploration for all parameters over their recommended range. Table 3 displays the simulation time of GWO and WOA in seconds. It can be perceived that GWO is faster than WOA in convergence regarding all used distributions and datasets. Table 4 presents the resulting parameter values. Notably, GWO and WOA selected almost the same parameter values in most cases.

WIND SPEED EXTRAPOLATION AND WIND ENERGY POTENTIAL
Wind data are generally collected at 10 m hub height; then, extrapolation can be done for any different height at the same site. The extrapolated wind speed at any height can be calculated as follows [44]: where V 1 and V 2 are the real and extrapolated wind speed values at the heights of h 1 and h 2 , respectively, the roughness  correlated to the shape of wind distribution. Wind power expectation is mainly based on power curves matching with wind distribution. The wind energy stored in the actual wind data is computed as follows [41]: where v 3 is the average of wind speed cubes, ρ is the air density (ρ = 1.225 kg/m 3 ), and T is the hourly time equals 730 and 8760 hours for monthly and annual periods, respectively.
The annual average wind speed is a preliminary sign of the suitability of wind energy generation at the selected location. Wind locations with average wind speed greater than 7.9 m/s at 30 m height are classified as excellent locations for generating wind energy [51]. Table 6 presents the annual average wind speeds at the three different heights at the studied site. These values reflect a good to a very good implication of wind potential compared to corresponding values in Table 7 [44]. Table 8 shows the annual average power density classes.

Evaluation metrics
Many statistical metrics have appeared in the literature to determine the fittest method to represent the actual wind data. In this study, five error metrics are utilized. The list of the applied evaluation metrics is briefly explained here: coefficient α varies due to many factors such as the nature of the site surface, steadiness of atmospheric, wind speed shape, and height interval [45]. Table 5 provides different values of α based on different surface characteristics [44]. Besides, the value of α can be calculated for any site when the wind speed values are available at any two different hub heights as follows [46], [47], [48]: where v 1 and v 2 are the wind speed values at the two different hub heights of h 1 and h 2 , respectively. There is another way to calculate α if the two wind speed values are not available using the under-analyzed site's roughness length (zo), as shown in Equation (33) In this study, z o is 0.0095 m [50] while h 1 equals 10 m. Moreover, α is computed by Equation (33) and found to be 0.143. Finally, α is checked at two different 10 m and 30 m heights using Equation (32).
Wind turbines transform kinetic wind energy into electrical energy. Wind turbines performance is fundamentally  between -1 (opposite correlation) and 1 (a perfect positive correlation). However, its zero value indicates that the two datasets are completely different (no correlation). The correlation coefficient is given by the following equation [56].
where (x -,y -) and (σ x , σ y ) denote the mean and the standard deviation of the measured wind speed vector and the predicted wind speed vector, respectively.
• Net Fitness is an outstanding metric that is used to average a group of metrics in a fairly way. It is used to make a decision over a group of metrics and simplify performing the ranks for comparisons goals. In this paper, the four metrics mentioned above are utilized to calculate net fitness. To find the net fitness, • Mean Absolute Error (MAE) is the mathematical mean between the collected (x) and the predicted (y) wind speed vectors, which created using any of the distributions as in Equation (35) where N is the vector length.
• Root Mean Square Error (RMSE) is the square root of the differences between the collected and predicted wind speeds [54]. It can be calculated as follows: • Regression Coefficient (R 2 ) shows the degree of linearity between the collected and predicted wind speed data, as given in Equation (37). If R 2 equals one, the relationship between the collected and predicted data can be linearly drawn [55].
where μ i is the i th mean of collected wind speed data.
• Correlation Coefficient (R) expresses the degree of correlation between two datasets. Its values are where n represents the total number of error metrics.

RESULTS AND DISCUSSION
The hourly data gathered from the Catalca site in the Marmara area for three years from 2008 to 2010 at different heights have been utilized to test the performance of the introduced five distributions. Table 9 provides the location information of the studied site. The wind speed data at this site were gathered at 10 m height then extrapolated at three different heights, namely, 30, 60, and 80 m, respectively.
Two optimization algorithms, GWO and WOA, were exploited to model the parameters for each distribution. The statistical descriptors such as mean, standard deviation, variance, minimum, maximum, skewness, kurtosis, and average power describe the characteristics of wind data for all datasets are shown in Tables 10, 11, and 12. It is obvious that the mean wind speeds of the collected data are (6.906, 7.646, 7.967), (7.251, 7.663, 8.000), and (7.875, 8.738, 9.124)  Table 7, it can be observed that the wind potential at the analyzed site occupies a good to very good rank. The standard deviation value slightly increases with an increase in the wind tower height. Variance, the square of standard deviation, is a measurement of the spread between wind speed values and their mean value. Also, variance slightly increases with an increase in the wind tower height. The minimum real wind speed was zero, whereas the maximum was varied between 25 and 33.5 m/s. The skewness indicates the level of asymmetry from the average wind speed. The skewness values of all data sets demonstrate that the actual data conform to the positive-skewness shape. Kurtosis indicates the tops of a recurrence distribution. There are three types of kurtosis: zero, positive, and negative. Zero kurtosis distribution generally follows Normal distribution; positive kurtosis has heavier tails and a higher peak than Normal, whereas Generalized Extreme Value the negative kurtosis has lighter tails and is flatter [58].
The real data of all datasets tend to the positive kurtosis (Leptokurtic  [59], [60], the underlying location can be classified in the fifth and sixth classes, which prove the location appropriateness for the installation of large-scale wind power turbines [61], [62]. Tables  10, 11, and 12, it can be noticed that Gamma distribution achieved the best matching, whereas Rayleigh distribution presented the worst matching.

Regarding the statistical descriptor values in
The distribution function reaches the optimal matching when the difference between the collected and the predicted wind speed values approaches zero [63], [64], [65]. Four evaluation metrics are used in this paper, namely, MAE, RMSE, R 2 , and R. Tables 13, 14, and 15 summary the GOF of the introduced distribution functions based on GWO and WOA. The bold values in these tables point to the best at each height, whilst the underlined values point to the best between GWO and WOA. In most cases, Gamma distribution achieved the best GOF. On the other hand, Rayleigh was the worst. In some cases, the 2008 dataset at 80 m height, Generalized Extreme Value distribution based on WOA provided the best matching in MAE, RMSE, and R 2 measures.
The rank evaluation based on the net fitness metric is required to specify the accuracy of the best distribution. Distributions are ranked according to four GOF criteria. The rankings are performed by a maximum of R 2 and R, whereas a minimum of MAE and RMSE.
Based on net fitness, the top-down rank of the five distributions based on GWO and WOA is Gamma, Generalized Extreme Value, Burr Type XII, Weibull, and Rayleigh. Table 16 shows the ranking of the five distributions. In addition, based on WOA, it can be observed that there is a slight difference between Weibull and Generalized Extreme Value to occupy the second rank. Hence, distributions based on different optimization methods provide different degrees of matching. Hence, it can be concluded that there is no unique, fully trusted, and best optimization method to assess the parameters of any distribution. Although Gamma distribution based on GWO and WOA outperformed the other distributions in terms of matching, it was the worst in computations. The computation time of Gamma was about sixteenth times of Generalized Extreme Value. Figures 2 and 3 depict the fitted PDFs and CDFs, respectively, to explicate the obtained results visually for all datasets. For the PDF and CDF plots, the horizontal axis represents the wind speed data in m/s. Regarding PDF plots, two different-scale vertical axes are used; the left axis is for the histogram of the measured wind data, whereas the right axis is for the introduced distributions. These vertical axes represent the probability density which varies between the lowest and highest possible values. The vertical axis for CDF plots represents the cumulative density. It is observed from Figures 2 and 3 that all the introduced distributions achieved good matching. Rayleigh distribution occupied the last rank in terms of matching, whereas it was the best in computation complexity.

Many valuable inferences can be deduced from this study as follows:
• One of the most paramount deductions is the distribution pattern of wind regime. This study confirms that wind distribution patterns vary from one location to another; therefore, different distribution functions should be applied to describe its pattern accurately. • The second important inference is that wind regime varies from one location to another and from height to another at the same location. Therefore, both utilized distributions and estimation methods need to be tested at different heights. • The third important deduction is the selection of the optimization method. The success of the optimization method mainly depends on the characteristics of wind speed data. For instance, both GWO and WOA with the Gamma distribution achieved the best matching; however, Rayleigh distribution did not. Therefore, various optimization methods need to be tested and compared. • The fourth deduction is the applied error metric. For instance, a specific error metric may demonstrate a particular distribution function as the best, but for another the worst. Hence, it is crucial to utilize many error metrics. Next, to decide the rank of the estimation technique accurately, the net fitness calculation is indispensable.

CONCLUSION
Selecting convenient distribution functions for representing wind speed frequency distribution is a critical for feasibility studies, wind turbine design, and long-term investment decisions. Moreover, selecting the suitable estimation method is also crucial since one method can achieve best GOF with specific distribution but may fail with others. This paper presents five different distributions to estimate the wind energy potential. The optimal parameter values for each distribution were selected based on GWO and WOA methods. The statistical characteristics of the analyzed site and the GOF of each distribution function were studied and compared via many statistical descriptors and different error measures. In most cases, Gamma distribution based on both GWO and WOA outperformed the other distributions for all datasets. On the contrary, it was the worst in respect of computation complexity. Conversely, Rayleigh distribution was the worst in matching, but the best in computation time. In most cases, GWO was more robust and faster than WOA. The skewness and kurtosis statistical descriptors are also crucial to describe the wind regime since they can display the whole wind distribution pattern. In this study, skewness and kurtosis values are positive, therefore, the wind pattern takes the shape of positively skewed and leptokurtic. Accordingly, the selection of the convenient distribution can be perceived. Besides, it can be observed that particular distribution with particular error measures can achieve the best matching, for example, in terms of RMSE and R2, Weibull distribution based on WOA at 30 m height for 2009 dataset achieved the best fitness, and in terms of MAE, Burr distribution was the best. Therefore, to accurately specify the best GOF, net fitness computation is required. Ultimately, the used distribution function, optimization method and error measure are essential factors to make an accurate decision of which the best GOF for wind pattern at any site. This study can be extended to include other distribution functions, optimization algorithms, and error measures.

PDF
Probability

DATA AVAILABILITY STATEMENT
The authors confirm that the data that supports the findings of this study are available within the article. Raw data that support the finding of this study are available from the corresponding author, upon reasonable request.

CONFLICT OF INTEREST
The author declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.