2.1 Weibull distribution

The statistical distribution of wind speed varies from region to region worldwide. These changes result from local climatic conditions, soil structure, terrain and surface topography. They can be statistically represented by a probability density and cumulative function. A probability density function may be described as a probability distribution which corresponds to the relative possibility for a random variable $v$ to take on a given value^{(Reference Kiran and Kiran51)}. The Weibull distribution, which was created by the Swedish physicist Waloddi Weibull in the 1930s, provides information regarding the wind speed distribution resulting from these changes. Since the Weibull distribution is a good fit to represent the wind speed distribution at a location, it is widely used in literature to estimate, model and analyse both the potential of wind energy and the wind speed at a specific location^{(Reference Akdağ and Dinler15)}. The probability density function of the Weibull distribution can be described as Equation (1):

(1) \begin{equation}f(v) = {\rm{\;}}\left( {\frac{k}{c}} \right){\left( {\frac{v}{c}} \right)^{k - 1}}\exp \left( { - {{\left( {\frac{v}{c}} \right)}^k}} \right),{\rm{\;}}\end{equation}

where $f(v)$ is the probability of occurrence of a wind speed $v$, $k$ represents the shape parameter and $c$ denotes the scale parameter. The cumulative distribution function corresponding to this probability density can be easily derived from Equation (2):

(2) \begin{equation}F(v) = \int\limits_{ - \infty }^v f(x){\rm{\;d}}x.{\rm{\;\;}}\end{equation}

Therefore, the cumulative distribution function of the Weibull distribution can be written as Equation (3):

(3) \begin{equation}F{\rm{\;}}(v) = {\rm{\;\;}}1 - \exp \left( { - {{\left( {\frac{v}{c}} \right)}^k}} \right),{\rm{\;\;\;}}\end{equation}

where $F(v)$ denotes the cumulative distribution function of the wind speed $v$.

The shape and scale parameters of the Weibull distribution play a crucial role in the wind distribution of a region. The shape parameter, which is an important parameter for the Weibull distribution, helps to calculate the frequency of the difference between the measured speed and wind speed^{(Reference Siddiqui52)}. The value of the shape parameter describes the variation in mean wind speed for a given sample^{(Reference Siddiqui52,Reference Rinne53)} . When the value of the shape parameter is high, the stability of the wind speed is also high. If the wind speed in a region varies little, the $k$ parameter is high, indicating that the wind speed is approximately constant (which may be low or high). High values of the shape parameter signify that the Weibull distribution is intensified around a given value and that the peak value of the distribution is higher. Meanwhile, low values of the shape parameter imply that a broad range of values are included in the distribution and that the peak value of the distribution is lower^{(Reference Chang17,Reference Dagdougui, Ouammi and Sacile54)} . The scale parameter directly indicates the wind potential at a location^{(Reference Siddiqui52,Reference Rinne53)} . The main function of the scale parameter is to fix the position of the curve. For a location with strong wind, it takes higher values, whereas it has lower values for locations with weak wind^{(Reference Dagdougui, Ouammi and Sacile54)}. This means that the $c$ parameter varies depending on the mean wind speed. If the mean wind speed is high, then the $c$ parameter is also high.

2.2 Geographic information

Eskişehir is located in the north-western part of Turkey. It has an area of 13,960km², representing approximately 1.75% of the territory of Turkey, and measures about 80km from north to south. Its population is 887,475 according to the Turkish Statistical Institute⁽⁵⁵⁾. It has its own climatic features. Like other cities in the Central Anatolian region, it has a terrestrial climate. The annual average temperature at Eskişehir, which is one of the coldest regions in Turkey, is 10.9°C. The mean wind speed at Eskişehir is approximately 3.355m/s. The proposed methods were applied to HPA, which is 5km from Eskişehir. The wind speed data were sampled hourly in the period between 1 January 2012 and 31 December 2014, at a height of 10m above ground level owing to operational conditions. A total of over 27,000 measurements were taken from the Turkish State Meteorological Service.

The geographical location of the site is shown in Fig. 1. The diamond symbol indicates the exact location of the meteorological station at the site. In addition to its geographical position, the location of HPA on the Wind Atlas is shown in Fig. 2. A detailed description of the site, including latitude and longitude as geographical coordinates, the elevation corresponding to the meteorological station and the period of the measured data, is presented in Table 1.

Figure 1. The geographical location of the study site on the map (revised from Google Earth).

Table 1 Description of the meteorological station and wind speed data

Figure 2. The location of HPA at Wind Atlas⁽⁵⁶⁾.

3.1 Graphical method

The graphical method, introduced by Deaves et al. in 1997, is directly based on the cumulative distribution function^{(Reference Deaves and Lines57)}. The Weibull parameters can be calculated by using the cumulative distribution function^{(Reference Akdağ and Dinler15,Reference Costa Rocha, de Sousa, de Andrade and da Silva19)} . Applying a logarithmic transform twice, the cumulative distribution function of the Weibull distribution, given in Equation (3), becomes Equation (4):

(4) \begin{equation}\ln\!\left( { - {\rm{ln}}\!\left( {1 - F(v)} \right)} \right) = k\ln(v) - k\ln(c).{\rm{\;}}\end{equation}

Note that there is a special relationship between $\ln \left( { - {\rm{ln}}\left( {1 - F(v)} \right)} \right)$ and $k\ln (v) - k{\rm{ln}}\left( c \right)$. By plotting $\ln \left( { - {\rm{ln}}\left( {1 - F\left( {{v_i}} \right)} \right)} \right)$ on the y-axis versus ${\rm{ln}}\left( {{v_i}} \right)$ on the x-axis for $i = 1,\;2, \ldots N$, an approximately straight line is obtained^{(Reference Mathew58)}. The slope of the straight line gives us the $k$ parameter, while the $c$ parameter is simply determined from the intersection with the y-ordinate^{(Reference Akdağ and Dinler15,Reference Kang, Ko and Huh59)} .

A second approach to determine the $k$ and $c$ parameters is the least-squares method. As can be seen, Equation (4) adopts the linear form $y = ax + b$ with respect to $\ln \left( { - \ln \left( {1 - F(v)} \right)} \right)$ and${\rm{\;ln}}(v)$. If the straight line fits best to $\ln \left( { - {\rm{ln}}\left( {1 - F\left( {{v_i}} \right)} \right)} \right)$ against ${\rm{ln}}\left( {{v_i}} \right)$ for $i = 1,2, \ldots N$ which is obtained using the least-squares method, the coefficients of the straight line, $a$ and $b$, can easily be calculated^{(Reference Akdağ and Dinler15,Reference Chang17)} . In such a case,

(5) \begin{equation}k = a \& c={\rm{exp}}\left( { - \frac{b}{a}} \right)\!.\end{equation}

It can be clearly seen that the slope of the straight line fitted to the data gives us the $k$ parameter while the $c$ parameter is simply determined from the y-intercept of the straight line in Equation (5)^{(Reference Akdağ and Dinler15,Reference Chang17,Reference Kang, Ko and Huh59)} .

3.2 Empirical method

The empirical method was proposed by Justus et al. in 1978^{(Reference Justus, Hargraves, Mikhail and Graber60)}. It is a simple approach that can be used to determine the Weibull parameters. For this purpose, it makes use of the standard deviation and mean wind speed. In this method, the shape and scale parameters can be obtained using Equations (6) and (7):

(6) \begin{equation}k = {\left( {\frac{\sigma }{{\mathop v\limits}}} \right)^{ - 1.086}},\end{equation}

(7) \begin{equation}c = \frac{{\bar{v}}}{{{\rm{\Gamma }}\!\left( {1 + \frac{1}{k}} \right)}},\end{equation}

where $\sigma $ represents the standard deviation, $\bar{v}$ denotes the mean wind speed and ${\rm{\Gamma }}$ indicates the gamma function given in Equation (8):

(8) \begin{equation}{\rm{\Gamma }}({\rm{z}}) = {\rm{\;}}\int\limits_0^\infty {t^{z - 1}}{e^{ - t}}dt\end{equation}

3.3 Maximum-likelihood method

The maximum-likelihood method, proposed by Stevens and Smulders in 1979, has a high computational load since the parameters of the Weibull distribution are determined iteratively^{(Reference Akdağ and Dinler15,Reference Chang17,Reference Stevens and Smulders61)} . For this reason, it is more difficult and complicated to apply this method compared with other methods. The shape and scale parameters can be calculated from Equations (9) and (10):

(9) \begin{equation}k = {\left( {\frac{{\mathop \sum \nolimits_{j = 1}^N {v_{}}^k\ln\!\left({{v_j}} \right)}}{{\mathop \sum \nolimits_{j = 1}^N {v_j}^k}} - \frac{{\mathop \sum \nolimits_{j = 1}^N \ln\!\left( {{v_j}} \right)}}{N}} \right)^{ - 1}},\end{equation}

(10) \begin{equation}c = {\left( {\frac{{\mathop \sum \nolimits_{j = 1}^N {v_j}^k}}{N}} \right)^{\frac{1}{k}}},\end{equation}

Iterative methods, such as the Newton–Raphson algorithm, can be used to find the $k$ and $c$ parameters. However, it should not be forgotten that zero velocities are excluded from the wind data due to the fact that there is no solution for zero velocity values.

3.4 Energy pattern factor method

The energy pattern factor method, also commonly known as the power density method, was proposed by Akdağ and Dinler in 2009^{(Reference Akdağ and Dinler15)}. This method is based on the energy pattern factor, a ratio that is often used in the aerodynamic design of turbine blades^{(Reference Singh, Bule, Khan and Ahmed62)}. Edward William Golding described it as the total power of the wind divided by the power computed by taking the cube of the mean wind speed, which is determined from sample measurements as shown in Equations (11) and (12)^{(Reference Haack63,Reference Mani and Mooley64)} :

(11) \begin{equation}{E_{pf}} = \frac{{{\rm{Total\;power\;of\;wind}}}}{{{\rm{Power\;computed\;by\;cubing\;the\;mean\;wind\;speed}}}},\end{equation}

(12) \begin{equation}{E_{pf}} = \frac{{\left( {\frac{1}{N}\mathop \sum \nolimits_{j = 1}^N {v_j}^3} \right)}}{{{{\left( {\frac{1}{N}\mathop \sum \nolimits_{j = 1}^N {v_j}} \right)}^3}}} = \frac{{\overline{{v^3}}}}{{{{( \bar{v} )}^3}}} = \frac{{{\rm{\Gamma }}\left( {1 + \frac{3}{k}} \right)}}{{{{\rm{\Gamma }}^3}\left( {1 + \frac{1}{k}} \right)}},\end{equation}

where ${\left( {\bar{v}} \right)^3}$ is the cube of the mean wind speed, $\mathop {{\overline{v^3}}}\limits$ is the mean of the wind speed cubed and $\frac{{\mathop {{\overline{v^3}}}\limits}}{{{{\left( \bar{v} \right)}^3}}}$ is defined as the energy pattern factor. ${E_{pf}}$ is always greater than 1 and, generally, its range is from 1.44 to 4.4 for most wind distributions^{(Reference Akdağ and Dinler15)}. After the energy pattern factor ${E_{pf}}$ has been calculated, the $k$ and $c$ parameters can be easily computed using Equations (13) and (14):

(13) \begin{equation}k = 1 + \frac{{3.69}}{{{{\left( {{E_{pf}}} \right)}^2}}},{\rm{\;\;\;}}\left( {13} \right)\end{equation}

(14) \begin{equation}c = \frac{\bar{v}}{{{\rm{\Gamma }}\left( {1 + \frac{1}{k}} \right)}}.\end{equation}

Figure 3. The framework of the study with numerical and metaheuristic methods.

4.1 Genetic algorithm-based data fitting method

Genetic algorithms are population-based, stochastic optimisation methods that are inspired by the phenomenon of natural selection. Genetic algorithms do not need any gradients or Hessians to operate in the search space, which most numerical methods would need. The crossover, mutation and selection operations are important phases of genetic algorithms. While crossover remains the main variation operator, mutation serves as a random perturbation. The selection of the next generation is based on the survival-of-the-fittest rule; that is, individuals with the highest probabilities will be transferred into offspring^{(Reference Ahn65,Reference Chong and Żak66)} .

The parameters of the genetic algorithms are presented in Table 2. The lower and upper bounds for the shape and scale parameters are chosen as [0, 100], which generously extends the search space. The initial values of the population are chosen as [0.1, 0.1], which are close to the lower bound to prevent any bias in the results. The population size is chosen as ten, which reduces the time demands of the proposed method at the cost of decreasing the chance of converging to the global optimal solution. The algorithm is run three times, and the mean value is taken as the final solution.

Table 2 Parameters of genetic algorithm

The pseudo-code of the proposed genetic algorithm-based data fitting method is given in Algorithm 1. The proposed algorithm starts with the initialisation of the parameters given in Table 2. Next, the counter i is reset and the for loop is started. After this, the while loop is created to evaluate the fitness of the population and also to apply the genetic operators for selection, crossover, mutation and replacement. When predefined conditions are met, for example on the generation number or a stall limit, the optimised variables are saved and the algorithm is restarted with the increased counter value until the run number n is reached. Next, the mean values of the optimised parameters are calculated and returned to the main program. Two separate fitness steps are evaluated, since two different objective functions are proposed. The roulette wheel selection method is used in the proposed algorithm, in which individuals with the highest probabilities will likely survive in the next generation. In the crossover phase, the weighted arithmetic means of the parents are used to create the offspring. In the mutation phase, randomly generated direction and steps are applied to the population to ensure exploration of the search space.

Algorithm 1. Genetic algorithm-based data fitting method

4.2 Particle swarm optimisation-based data fitting method

Particle swarm optimisation is an intelligent method which was introduced in Ref. Reference Kennedy and Eberhart26 as a tool to optimise continuous nonlinear functions. The algorithm is inspired by the social and cognitive behaviour of flocks of birds, which is computationally inexpensive due to the primitive mathematical operators used^{(Reference Kennedy and Eberhart26)}. Each individual in the swarm is called a particle. A particle is a candidate for a global solution which has:

• • a tendency for randomness in the search space

• • a tendency to return to the best location visited by itself and

• • a tendency to go where most particles in the neighborhood go^{(Reference Gümüşboğa and İftar67)}.

The equations to update the position and velocity of the particles are given in Equations (15) and (16).

(15) \begin{equation}{v_{t + 1}} = {\rm{\;}}w\,{\rm{*}}\,{v_t} + {\rm{\;}}{c_1}\,{\rm{*}}\,{u_1}\left( {{p_t} - {\rm{\;}}{x_t}} \right) + {\rm{\;}}{c_2}\,{\rm{*}}\,{u_2}\,{\rm{*}}\,\left( {{g_t} - {x_t}} \right),\end{equation}

(16) \begin{equation}{x_{t + 1}} = {\rm{\;}}{x_t} + {\rm{\;}}{v_{t + 1}},\end{equation}

where ${v_t}$ is the current velocity, ${v_{t + 1}}$ is the updated velocity, $w$ is the inertial weight, ${c_1}$ is the cognitive weight, ${u_1}$ and ${u_2}$ are uniform random numbers $\; \in \left[ {0,1} \right]$, ${p_t}$ is the best position of the particle, ${c_2}$ is the social weight, ${g_t}$ is the best global position of the swarm, ${x_t}$ is the current position and ${x_{t + 1}}$ is the updated position^{(Reference Carneiro, Melo, Carvalho and Braga31)}.

The parameters for the particle swarm optimisation method are presented in Table 3. It can be seen from Tables 2 and 3 that all the common parameters are chosen the same to achieve a fair comparison between the metaheuristic algorithms. These parameters are swarm size, iteration number, stall limit, lower and upper bounds, initialisation values and the total number of runs of the proposed algorithms. The inertial weight in the particle swarm algorithm lies in the range of 0.1–1.1, and its value is assigned adaptively. The cognitive and social weights are chosen as 1.49, whereas in literature they are generally taken in the range of 1–2^{(Reference Li, Cui, Li and Cui68)}.

Table 3 Parameters of particle swarm optimisation

The pseudo-code for the proposed particle swarm optimisation-based data fitting method is given in Algorithm 2. The initialisation step is conducted according to Table 3. The counter j is set to zero, and the for loop is started. Next, a while loop is created to assess the fitness measures of the particles and also implement the exploration as well as cognitive and social behaviour of the particles. When the maximum number of iterations is reached or a stall occurs, the gbest value is stored and the proposed method is restarted until the run number n is reached. Finally, the average values of the gbest values are calculated and returned to the main program as optimal solutions. The algorithm is evaluated twice, since two different fitness functions are proposed.

Algorithm 2. PSO algorithm-based data fitting method

4.3 Fitness functions

The objective function is an important part of the metaheuristic methods, since their optimisation performance is directly related to the fitness measure. The definition of the objective function clearly determines how close an individual is to the global minimum. Therefore, two different objective functions are proposed in this paper to achieve a higher performance standard according to the different figures of merit. While the first objective function deals with minimisation of the probability density error, the second deals with minimisation of the cumulative distribution error. The proposed objective functions are given in Equations (17) and (18), respectively.

(17) \begin{equation}fi{t_{pdf}} = {\rm{\;}}\sqrt {\frac{1}{{n}}\mathop \sum \limits_{i = 0}^n {{\left( {{\rm{\;}}{f_{measPdf}}\!\left( {{v_i}} \right) - {\rm{\;\;}}{f_{calcPdf}}\!\left( {{v_i}} \right)} \right)}^2}},\end{equation}

(18) \begin{equation}fi{t_{cdf}} = {\rm{\;}}\sqrt {\frac{1}{n}\sum\limits_{i = 0}^n {{\left( {{\rm{\;}}{f_{measCdf}}\left( {{v_i}} \right) - {\rm{\;\;}}{f_{calcCdf}}\left( {{v_i}} \right){\rm{\;}}} \right)}^2}{\rm{\;}}} {\rm{\;}},\end{equation}

where $fi{t_{pdf}}$ is the first proposed objective function dependent on the RMSE of the probability density, $fi{t_{cdf}}$ is the second proposed objective function dependent on the RMSE of the cumulative distribution, $n$ is the number of unique wind speed values, $\;{f_{measPdf}}$ and $\;{f_{measCdf}}$ are the measured probability and cumulative functions and, lastly, $\;{f_{calcPdf}}$ and $\;{f_{calcCdf}}$ are the calculated probability and cumulative functions based on the $k$ and $c$ values obtained from the proposed methods, respectively.

5.1 Root-mean-square error

The Root-Mean-Square Error (RMSE) is a statistic that is often used for evaluating the performance of a model. The RMSE of a model, with regard to a given dataset, is the square root of the mean of the square of each estimation error on all the samples in the dataset. Each individual estimation error is equal to the difference between the actual observation and the estimated observation for the sample^{(Reference Fürnkranz72)}. The RMSE can be expressed mathematically as

(19) \begin{equation}RMSE{\rm =}\sqrt{\frac{{\rm 1}}{N}\sum^N_{k{\rm =1}}{{\left(y_k{\rm -}{\hat{y}}_k\right)}^{{\rm 2}}.}}\end{equation}

It can be seen that, the smaller the RMSE value is, the better the accuracy of the model.

5.2 R squared

R squared (R ²), which is also defined as the coefficient of determination, is another way of measuring the error between the data and model. R ² values range from 0 to 1. As the value of R ² converges to 1, the accuracy of the model increases. It can be defined as follows:

(20) \begin{equation}R^{{\rm 2}}{\rm =1-}\frac{\sum^N_{k{\rm =1}}{{\left(y_k{\rm -}{\hat{y}}_k\right)}^{{\rm 2}}}}{\sum^N_{k{\rm =1}}{{\left(y_k{\rm -}\overline{y}\right)}^{{\rm 2}}}}\end{equation}

where $\bar{y}$ represents the mean wind speed.

5.3 Chi squared

Chi squared is used to determine the magnitude of the difference between the actual and estimated observations. If the actual observation exactly matches with the estimated observation, then the value of chi squared will be zero. As the difference between the actual and estimated observations grows, the value of chi squared will increase as well. It can be expressed as follows:

(21) \begin{equation}X^{{\rm 2}}{\rm =}\frac{\sum^N_{k{\rm =1}}{{\left(y_k{\rm -}{\hat{y}}_k\right)}^{{\rm 2}}}}{N{\rm -}n}\end{equation}

where $n$ is the number of parameters used in the model.

5.4 Mean absolute error

The Mean Absolute Error (MAE) of a model, with regard to a given dataset, is the mean of the absolute values of each estimation error on all the samples in the dataset. Each individual estimation error is equal to the difference between the actual and estimated observation for the sample^{(Reference Fürnkranz72)}. The MAE can be expressed mathematically as

(22) \begin{equation}MAE{\rm =}\frac{{\rm 1}}{N}\sum^N_{k{\rm =1}}{{\rm |}y_k{\rm -}{\hat{y}}_k{\rm |}},\end{equation}

where ${y_k}$ represents the actual observation, ${\hat{y}_k}$ depicts the estimated observation, $N$ denotes the number of observations and $\left| \cdot \right|$ is the absolute value operator.

Table 4 Monthly and yearly mean and standard deviation values of the measured wind speed of HPA