Evaluation of Wind Energy Potential Using an Optimum Approach based on Maximum Distance Metric

The integration of wind power as an alternative energy source has gotten much attention globally. In this paper, the Weibull distribution model based on different artificial intelligent algorithms and numerical methods is used to evaluate the wind profile. The application of Weibull distribution in wind data assessment can be extensively found, but the methods applied for estimating the parameters still need improvement. Three artificial intelligent algorithms are presented as an alternative method for estimation of Weibull parameters, and an objective function is proposed through the concept of maximum distance metric. Its convergence was proven mathematically through its boundedness for all wind data types. The optimization methods based on the proposed objective function are compared with the conventional numerical approaches for Weibull parameter estimation. Two-year wind data from the site in the southern area of Pakistan has been used to conduct this analysis. Furthermore, this work provides an eloquent way for the selection of a suitable area, evaluation of parameters, and appropriate wind turbine models through real-time data for power production.


Introduction
The renewable energy sector has been encouraged to develop due to its endless resource, lower index of greenhouse gas emissions, and better economics. The World Wind Energy Association highlights the increased development of wind energy investments in recent years [1]. According to their report, the world's total wind installed capacity reached a staggering 597 GW in the year 2018, with a high rate of progress about 13% in the last three years [1,2].
Wind resource analyses are performed to predict the performance and output of potential power projects. This allows investors and developers to gather more accurate feasibility of the project under consideration. The process consists of many steps and methods, but one of the most critical steps is wind characterization through distribution functions. Several distribution models have been used for this purpose. Some of the most commonly used are the Rayleigh distribution, the Gaussian distribution, and the two-parameter Weibull (Wb) distribution [3][4][5][6][7][8]. Luna and Church [9] determined that the application of probability distribution on wind data depends more on the considered data. According to their analysis, lognormal distribution was the best representation of wind data. Several distribution functions have been utilized for wind data representation and can be found in the literature extensively [10][11][12][13].
The proposed site is located near Sujawal city in the south-east coastal area of Pakistan and about 160 km from provincial capital (Sindh) Karachi. The geographical location of the site is latitude 24.515563° N, longitude 68.18865° E, and an elevation of 17 m. It is an open area site and easily accessible for all kinds of vehicles. The land profile of the site is mostly flat with no obstacle in the radius of 300 m from the Met-mast. The surface roughness length was found to be 0.054 m for the open field with low buildings and trees. A vertical triangular lattice tower, usually named metrological mast (Met-mast), was designed and installed at the proposed location according to the up-to-date version of the IEC Standard 61400-12-1 [30]. The data is calculated for two years (June 2017-May 2018) at 10 minutes average of time. The geographical map of Sujawal city is shown in Figure 1.

Data Analysis Method
The wind data analysis method is shown in Figure 2. After the collection of wind data, it is analyzed through two approaches. The first approach is through the data descriptors from the real data. This includes the wind mean, kurtosis, skewness, standard deviation, and coefficient of variance. The second approach is the Weibull model for the given wind data. The Weibull parameter estimation is performed through two approaches. Four numerical methods and three AI optimization methods based on the maximum-distance metric are used and compared to determine the better approach for Weibull parameter estimation, and the mentioned aspects are discussed below.

Weibull Distribution Function and Estimation Methods
The Weibull probability density function and Weibull cumulative distribution function are expressed as follows [11][12][13][14][15][16]: Weibull (Wb) parameters have been computed through several numerical methods. This study compares four of the most commonly used methods with artificial intelligence optimization algorithms. The four methods selected here are given as follows: i. Empirical method of Justus (EMJ) [31,16], ii.

Artificial Optimization Algorithm for Weibull Parameters Estimation
Based on natural mechanisms, artificial intelligence optimization algorithms can be used to search the optimal parameter values. These algorithms have vast applications in many fields. The basic structure of these algorithms is defined as follows.
 Initialize the search agents and their positions. In order for an optimization algorithm to converge, the most important criteria is the convergence of the objective function for the given data. Therefore, the selection of these objective functions can influence the final outcome. The objective functions are based on the error formulas. Some of the conventional errors that have been used in wind data analysis are the mean bias error, root mean square error and many more.

Explanation:
Assume that 1 and * , , then the metric becomes mean bias error. Similarly, if 1/2 and * , , then becomes the root mean square error. This is the definition of mean-based errors that are conventionally used with artificial intelligence optimization techniques. However, it has been observed that these objective functions are divergent in many cases of Weibull parameter estimation. The above theorem states the property that makes these objective functions divergent.

Explanation:
Now define a metric as follows: Definition: Let ℝ be the set of length vectors with real components. Assume that is a metric on ℝ ℝ , based on 1-D maximum distance, defined as , maximum | |.
The properties defined above allow to always converge for Weibull parameter estimation.

Particle Swarm Optimization (PSO) Algorithm
Eberhart and Kennedy [33] proposed the idea of an artificial intelligence optimization algorithm based on the behavior of particle swarms. In the PSO algorithm, each particle is considered to be a part of multidimensional space, adjusting its position for the optimum results in the swarm on each iteration. The position of the swarm is also considered in the process. It optimizes itself due to the continuous motion of the particles. Therefore, the whole swarm moves towards the optimum solution and ultimately finishes the search with the best value. Assume that for each particle , the best position is , . This position is termed as the local best. Then the best result termed as global best position among the all particles is recorded , . The movement of each particle depends on the particle's initial position and the weighted average of the particle's best position and swarm's best position. It can be mathematically defined as [33] , , , . Here, where, ώ stands for particle inertia weight, and refer to acceleration coefficients, µ µ , , and are positive random constants within the range of (0,1), and ώ is expressed as Here, represents the maximum number of iterations. Further, the PSO algorithm can be summarized as follows:  Initialization of particle space and positive random constants.  Evaluate and calculate the objective function for each particle.  Evaluate the best position (by updating local best) of each particle and global best within all of the particles.  Update the particles according to the initial step.  Termination criterion: stop if the output is the global best, otherwise, repeat the iteration process and return to the objective function calculation step for the best output.

Cuckoo Search optimization (CSO) Algorithm
Many fruit flies and bird species have Lévy flight behavior. It was combined with cuckoo species' brood parasitic behavior to form an artificial intelligence optimization algorithm. Cuckoo has a mechanism of searching the optimum nest for the hatching of eggs. The bird considers an initial set of nests for the possible solution. These nests are updated based on Lévy distribution. In this context, the Mantegna′s algorithm [34] is utilized with the Lévy flight of Cuckoo to find the optimal solutions. Mathematically, the new position is stated as û where, , is a constant, usually 0.01. û and ṽ are obtained through normal distribution, is the new solution generated from . ά is the scale factor with a value of 1.5, and is calculated according to the following equation: The Cuckoo discards the weak solutions and keeps better solutions for the next iteration, which is achieved through the probability of detection, represented as . The new solution is obtained from the following expression [35,36]: In Equation (13), stands for a uniformly distributed random number of intervals (0, 1), whereas , and , represent two random solutions for the ℎ generation.
Cuckoo Search Optimization algorithm can be summarized in the following steps:  Parameters initialization: generate the initial population of nests (q) and assume the objective function (Fq).  Compare each nest and record the best solution, including its objective function.  Generate new solution (i) by the Levy flight and calculate the new objective function (Fi)  If then exchange "q" for the new solution generated by "i", otherwise keep the solution "q", while the nest fraction smaller than "Pro" is abandoned, and new nets are generated.  If the required iteration is reached with an optimal solution, continue to the next step; otherwise, go to the objective function (Fi) calculation step.  Finally, the output is the best solution or best nest.

Grey Wolf Optimization (GWO) Algorithm
Another similar artificial intelligence optimization algorithm is the Grey Wolf optimization algorithm. This algorithm was developed in 2014 and has been used in various applications ever since [37]. Based on the leadership hierarchy of gray wolves, this method initializes a set of wolves and their positions. It follows the gray wolf leadership hierarchy and selects the leaders based on the objective function. If the current location is not acceptable, the distance vector is computed. The new position is determined based on the distance-vector, and the process is repeated. The distance vector is determined as follows: .
Here, and represent the pray position and wolf position, respectively, denotes iteration number, whereas and are the vectors of coefficients that are computed as 2 , 2 .
Here shows a linear dropping from 2 to 0 throughout the iterations and , 0,1 . The new position is determined based on the following formula: GWO is based on the following steps:  Parameters initialization: generate the initial population of search agents, leaders, iteration count, and search domain.  Compare the agents' fitness based on the objective function.  If the location is acceptable, then terminate. If not, find a new distance through Equation (14) and find new leaders through Equation (17).  Repeat from the second step.  Finally, the output is the best location.

Wind Shear and Wind Indicators
Wind shear is a variation of the speed of wind with altitudes, and due to the air viscosity and drag friction of surface (land or sea), the air flows faster at higher altitudes. Usually, daytime variations follow the one-seventh rule of thumb (1/7 power law), which predicts a proportional increase in wind speed to the 7th root of altitude. The power-law model is generally used to characterize the effect of roughness of the earth's surface on wind speed, which is expressed as [38,39] , where and represent the wind velocities at the respective heights and , respectively, and α is the friction coefficient, called wind shear exponent, which is given as [38,40] .
The critical indicators for wind data are the power density (PD), which depicts the power generated by numerous wind speeds at a specific location. Both real wind data collected from an area and the Wb function could be utilized to calculate wind PD and energy density (ED). The formulations are provided as [11,20] 1 2 , Further relevant numerical descriptors for wind data are the speed average , speed variance , speed standard deviation , maximum energy-carrying wind speed (VmaxE), and the most probable wind speed (Vmp) are defined in [16,19,41,42].

Fitness Test Analysis
Many statistical methods have been used to compute the errors and then be analyzed to find the better fit technique for these calculations. The performance test analysis methods are considered in this study are root mean square error (RMSE), the coefficient of determination (R 2 ), and mean bias error (MBE) and expressed as [16,28,42] ∑ , 1 .
where m, y , x and , respectively, show average wind, actual i wind, predicted i wind, and the observation count of wind speeds. The average accuracy percentage is computed as follows: 2 . 100.
Here, is the average coefficient of determination, is the average root mean square and is the mean bias errors over the four heights.

Results and Discussion
Data for this work comprises real-time wind data for a period of two years with the measurement of ten minutes intervals at different heights. Wind shear exponent is measured using the power-law model to analyze the vertical deviation of wind speed. The Weibull distribution is employed to appraise the wind power potential and wind characteristics of the given site. Three optimization algorithms and four numerical methods are compared to elicit the optimal Weibull parameters. Furthermore, turbulence intensity, extreme wind speed, wind directions, and annual energy production are examined to give an optimal assessment of the respective area.

Wind Speed Characteristics Analysis
The statistical descriptors such as mean wind speed standard deviation (S.D), power density (PD), energy density (ED) skewness (Sk.), kurtosis (Kur.), and coefficient of variance (CV) defining the behavior of wind data are presented in Table 1. It is evident that the two-year mean wind speeds of the site are 5.3, 6.2, 6.9, and 7.4 m/s at mast heights 20, 40, 60, and 80 m, respectively. The standard deviation value changes with an increase in the Met-mast height. However, this change is considerably small when compared to the change in mean wind speeds. So, it can be seen that skewness for both years is slightly negative at the height 80 m, but it stays positive for all other heights with a decreasing trend as Met-mast height increases. It is because of which increase in the mast height, the mean wind speed changes, which indicates a more stable wind at higher mast heights. The kurtosis of the wind also follows a similar trend where it decreases as the height is increased. Figure 3 depicts monthly mean wind speed (MWS) variations at different heights for the years 2016-2017 and 2017-2018. It can be analyzed from the monthly statistics; March to September is found to be the best wind speed period in both years at all mast heights. The months of May and June show the highest mean wind speeds, while the lowest mean wind speed occurred in the months of January and February for the two years at four-mast heights. It was observed that mean wind speed for the mast height above 60 m in all months was above 5 m/s.

Wind Shear and Turbulence Intensity Analysis
Figure 4a-c shows the wind shear profile, monthly mean wind shear, and diurnal shear exponent (α) of the investigated site. The exponent "α" of power-law was evaluated for the atmospheric boundary layer, which is calculated by using 10 minutes measured mean wind speed data at heights of 40, 60, and 80 m. The mean value of "α" was 0.228 at the selected site. The monthly highest values of "α" were in January (0.307) and February (0.319), while the lowest was in June (0.140) and July (0.144). It is apparent in Figure 3c that the wind shear exponent reaches its minimum value during the day time due to the higher temperature. The shear exponent starts increasing at 1400 hrs as the temperature starts decreasing and reaches its maximum value in the midnight.   T.I is calculated using 10 minutes time interval measured wind speed and standard deviation data of the candidate site. It was observed that T.I reduces as the height increases. The variation in T.I level at 80 and 60 m heights was found to be less as compared to T.I variation at 40 and 20 m because the wind speed at higher heights becomes smoother and steadier. The International Electro-Technical Commission (IEC) designed standards IEC 61400-1 for wind turbines, and it states that the allowable T.I level is 18% for 15 m/s wind speed [43,44]. The variation in T.I level at the selected site was found in the allowable range of IEC, as exhibited in Figure 5. Moreover, the surface roughness length at the site was found 0.054 m for the open field [30,45]; thus, the site has an acceptable range of wind shear and T.I level.

Extreme Wind Speed and Direction Analysis
Wind turbine optimization is essential for the maximum output of a wind system. In this context, the analysis of extreme wind speed values and wind direction is essential. Figure 6 represents the annually extreme wind speed and the return period of extreme wind speed values at different heights. Gumbel distribution function is applied to compute the value of extreme wind speed (EWS) [46][47][48]   In the following year, 2017-2018, there is no significant change in the direction of the wind; around 55 percent wind lies in between 240 and 270 degrees, and almost 17 percent lies between 30 and 45 degrees clockwise. The maximum percentage of the wind lies in the south-west direction, which is ideal for wind turbines (WTs). Therefore, WTs are able to capture the maximum optimum power from the wind and work effectively, instead of changing the direction frequently. Table 2 illustrated the Wb parameters "k" and "c" at different heights for both years (2016-2018) calculated using the seven methods. It is apparent from the table that the "k" ranged from 2.37 to 3.24, whereas the "c" ranged from 5.79 to 8.34 for the heights of 20 to 80 m. The shape parameter value increase with the increase in height, which indicates steady wind at the candidate site. The scale parameter for the two years was 7.1 m/s, which is almost the same as that for the two years independently. Fitness tests are used to compute the accuracy of the used methods. These tests also help us to compare different methods with each other. In this context, MBE, R 2 , and RMSE errors are applied to the Weibull distribution curves and real data. Table A1 (see appendix) illustrates the fitness test results for different approaches. It shall be noted that all the used methods show acceptable fitness, suggesting the good fit of wind data with Weibull distribution. The accuracy percentage is also illustrated in Table A1. It shall be noted that although all the used methods show an acceptable accuracy, the accuracy percentage of all the AI methods is over 99%. The good fit of these methods can also be seen in Figure 8a-d, where the probability density function and cumulative distribution for Weibull are presented at all mast heights for the year 2017-18. Moreover, the cumulative frequency for wind speeds more than 4 m/s was found to be 71%, 79%, 82%, and 85% at 20, 40, 60, and 80 m mast heights, respectively. It also exhibited that the probability for wind speed higher than 20 m/s is zero. The PDFs and CDFs for the year 2017-2018 at considered heights are depicted in Figure A1a-d. It is challenging to analyze the fitness of distribution with multiple error functions. Therefore, a new fitness criterion is proposed based on the three error formulas. It is defined as

Weibull Parameters Test Analysis
where n is the total number of error entries. Figure 9 presents the comparison of net fitness of the methods used for the estimation of Weibull parameters. It is apparent that all methods at all mast heights show acceptable performance. The error for AI algorithms ranges between 0.004276 and 0.005583, while that of the numerical methods ranges between 0.056606 to 0.067614. Figure 10 compares the convergence for all AI methods at the height of 80 m with 100 iterations. In both years, the grey wolf algorithm converges the fastest. However, all the methods perform similarly for the higher number of iterations.

Wind Energy Density and Specific Wind Characteristics
For the year 2016-2017, Weibull MWS, specific wind speeds, power density, and energy density at four heights are summarized in Table 3. The average wind speed ranged from 5.3 to 7.42 m/s, while Vmp and VmaxE are found between 4.7 to 7.2 m/s and 7.7 to 10.1 m/s, respectively. The annual Weibull PD and ED ranged between 147 to 351 W/m 2 and 1993 to 3079 kWh/m 2 , respectively. The annual Weibull wind PD and ED, mean wind speed, and specific wind speeds for the year 2017-2018 are given in Table A2. Table 3. Weibull mean wind speed, specific wind speeds, power, and energy densities for the year 2016-2017. 1323.781 According to wind power classes (WPC), regions that fall in WPC-III or above are suitable for most wind power conversion systems [19,49]. Results of measured wind speed and power density revealed that the Sujawal area stands in WPC-III and is suitable to harness the wind power.

Energy Production and Cost analysis
Total energy production and capacity factor are fundamental aspects of a wind power project. Furthermore, the capacity factor (CF) can be used to characterize the annual energy production (AEP) of the WTs. From [50], the random nature of wind usually affects the value of CF, which in result affects the efficiency and serviceability of the WTs. In practice, CF generally varies from 20% to 70%. The technical specifications of selected WTs are illustrated in Table 4. The energy produced is calculated using the WTs power curves data and the number of hours when the possibility of wind speed stayed in different bins at the turbine hub height; energy production, CF, and cost analysis are defined in [19,42,51]. The wind speed for different hub heights of the WTs is calculated using the power-law extrapolation model. The capacity factor (CF) and annual energy production (AEP) and of WTs with different power ratings ranged from 0.33 to 2.75 MW are estimated, respectively, and shown in Table 5. The calculations are made by assuming constant losses of 16.005%, i.e., turbine performance losses, wake effect losses, availability losses, and electrical losses. Table 5. WTs power production, annual energy production (AEP), capacity factor (CF), and cost/MWh. The CF is ranged between 0.20 to 0.51 for the selected WTs at their corresponding hub heights. The maximum CF is 0.51 for Gamesa G114-2.5, while minimum CF is 0.20 for Nordex N60/1300 wind turbine.

WTs Model
The estimated cost of energy per unit for WTs is given in Table 5. The cost of energy per unit is calculated by taking the assumptions (for instance, initial investment cost, operation and maintenance cost, real interest rate, life span, and cost of wind turbines) from the literature [19,42,[51][52][53]. The lowest cost of energy (COE) per unit is 21.31 US Dollar $/MWh for Windtec DF 2000 followed by Gamesa G114-2.5 and GE 1.7-100 with COE/unit 21.63 $/MWh and 22.72 $/MWh, respectively, while the highest COE is 54.47 $/MWh for Nordex N60/1300 followed by Enercon E-48/800 and Enercon E-33/330 with COE/unit 47.75 $/MWh and 43.32 $/MWh, respectively. Moreover, the cost of energy per unit ranges from 21.31 to 54.47 $/MWh. The capacity factor and cost/unit of all WTs are found in the appropriate range at the proposed site. Moreover, the capacity factor and cost/unit of selected WTs show that the Gamesa G114-2.5, Windtec DF 2000-WT55, and GE 1.7-100 wind turbines with the highest CF and lowest cost/unit are the most suitable for the selected site. Therefore, the WT models with corresponding specifications and power curves can be suitable for power production in the proposed area.

Conclusion
The Weibull distribution is considered as the most capable model for wind resource assessment in any region worldwide. Nevertheless, the precise estimation of its parameters remains a challenge for wind community researchers and developers. In this paper, a new criterion is introduced to estimate the Weibull parameters, which is based on the maximum distance metric. The proposed convergence objection function results are verified through mathematical expression. Furthermore, three artificial intelligence (AI) optimization algorithms are compared with four numerical methods using 2-year wind data (acquired from the Sujawal site). The goodness-fit tests show that the optimization algorithms have adequate accuracy for Weibull model parameter estimation in comparison to numerical methods.
The net fitness results for AI algorithms are ranged between 0.004276 and 0.005583, while that of the numerical methods ranged between 0.056606 to 0.067614. Furthermore, the average accuracy percentage of all the considered AI methods is over 99%. According to these findings, optimization algorithms present supremacy over numerical approaches. Also, the proposed convergent objective functions based on the maximum distance metric give a better approach for the optimization algorithms. The outcomes of the data analysis showed that the site has an energy density of about 3014.477 kWh/m 2 , and a power density of 344.118 W/m 2 at the 80 m height.
The Windtec DF 2000-WT55 and Gamesa G114-2.5 wind turbines showed the highest annual energy production with good capacity factor and lowest cost of energy per unit, whereas, the capacity factor was found between 0.20 and 0.51 for the selected wind turbines. It is determined that the wind turbine models with corresponding specifications and power curves can be suitable for the investigated area. Besides, artificial intelligence optimization algorithms can be an alternative approach to the conventional numerical methods for Weibull model parameter estimation.