Assessing wind speed simulation methods

Renewable and Sustainable Energy Reviews 56 (2016) 473–483

Contents lists available at ScienceDirect

Renewable and Sustainable Energy Reviews journal homepage: www.elsevier.com/locate/rser

Assessing wind speed simulation methods Andrés Feijóo n, Daniel Villanueva Departamento de Enxeñería Eléctrica, Universidade de Vigo, EEI, Campus de Lagoas, 36310 Vigo, Spain

art ic l e i nf o

a b s t r a c t

Article history: Received 6 February 2015 Received in revised form 28 October 2015 Accepted 30 November 2015

In this paper simulation methods of wind speed series at different locations are reviewed. Over the last few decades wind farms (WF) have increasingly been introduced into electric power systems (EPS) in many countries. The strong relationship between wind speed and the power generated by a wind energy converter (WEC) has led researchers to reflect on the need to develop adequate models for simulating wind speed data. In recent years some proposals to meet this need have been published in the literature, and this paper aims to gather and discuss them. & 2015 Elsevier Ltd. All rights reserved.

Keywords: Wind speed Weibull distribution Rayleigh distribution Correlation Autoregressive model Markov process Autocorrelation

Contents 1. 2.

Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Wind speed distributions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.1. Weibull and Rayleigh distributions. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.2. Other distributions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3. Correlation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.1. Inducing Pearson parametric correlations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2. The use of evolutionary algorithms . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.3. Inducing Spearman rank correlations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.4. A different method for inducing Pearson correlations. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.5. Example . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4. Autocorrelation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5. Markov processes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.1. First-order Markov processes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.2. Second-order Markov processes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.3. Examples . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.4. Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.5. Correlation in Markov processes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6. Some comments about other methods. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7. Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . References . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

n

Corresponding author. Tel.: þ 34 986 812055. E-mail addresses: [email protected] (A. Feijóo), [email protected] (D. Villanueva).

http://dx.doi.org/10.1016/j.rser.2015.11.094 1364-0321/& 2015 Elsevier Ltd. All rights reserved.

474 474 474 475 475 476 476 477 477 478 478 478 479 479 480 480 480 481 481 482

474

A. Feijóo, D. Villanueva / Renewable and Sustainable Energy Reviews 56 (2016) 473–483

1. Introduction Integration of wind power into EPS needs to be carried out with knowledge of the behaviour of such an energy source. The possibility of fitting its cumulative behaviour to given distributions in order to make statistical analysis easier has been an issue over the years [1,2]. One of the main concerns facing those who manage wind energy in EPS is the difficulty of controlling power obtained from such a variable resource. The impossibility of storing wind directly as a primary resource together with the difficulty of predicting its behaviour has created a need to develop adequate procedures for correct simulation. It must be said that prediction tools based on meteorological models have improved over recent years and in some markets, such as the Spanish one, part of the wind energy being injected in the system is commercialised daily in the electricity pool. This is possible thanks to such models, as they provide some guarantee as to the amount that can be injected [3]. In the early years of wind energy development, some companies, particularly transmission system operators (TSO), saw the concept of simultaneousness as an interesting question that needed an answer. It is a concept that has been generally applied to loads. Conventionally, it is considered that a power flow can be run for a given area under the assumption that only a given percentage of the total load is being applied, and not the whole bulk of the installed capacity. According to this, a figure below 100% of total load amount is called the simultaneousness of the load. The question begging an answer is whether a similar concept can be defined for wind generation. In the case of loads, thousands of recorded operation hours allow operators to establish certain rules that result in that figure called simultaneousness. For conventional power stations (nuclear, coal, hydro), the possibility of programming them on the basis of the availability of their primary resources allows stakeholders to limit the degree of uncertainty in the calculations. This is not the case for wind energy, due to its variable behaviour. The way the wind behaves can determine different decisions in an EPS. Thus, the knowledge of what kind of wind speed distribution is to be handled is important and, when there are several WFs in a given network, which is the general case, the correlation between their wind speeds is also important. For instance, in a probabilistic load flow analysis, changes in these variables (distributions, correlations) can determine different solutions. There has been a steady flow of research and literature dealing with this topic over recent years, and the goal of this paper is to review that body of work and add to the discussion. The rest of the paper is devoted to discussing different issues involved in the simulation of wind speed series. In Section 2 the most used continuous statistical distributions for the description of wind speed cumulative behaviour are reviewed. In Section 3 the important concept of correlation is introduced and several methods that have been used in the literature for inducing correlation in uncorrelated series of data are explained. In Section 4 the introduction of autocorrelation in the simulation is explained. Section 5 is about Markov processes and their relationship with autocorrelation and chronological features. Finally, after very brief comments in Section 6 about other possible methods, Section 7 adds conclusions and some more discussion.

2. Wind speed distributions The main feature of wind speed is its variability. Disregarding the fact that wind speed at a given location in particular can be constant over long periods of time, it has been generally confirmed that wind speed is a very variable resource at any place. This

variability can be observed in time and in space. Wind speed is variable at a given location over time, but it can also be very different at two different locations simultaneously. Wind speed can be considered as a vector, i.e. a pair of numbers, one of them corresponding to its mean value for a fixed period of time, and the other one indicating a direction. In this paper the direction of wind speed will not be considered, and the discussion will be only about mean values. There are interesting works in the literature where the direction of wind has been considered when modelling it stochastically [4]. There seems to be general agreement on the fact that wind speed in a given location can be cumulatively represented by means of a Weibull distribution [5,6]. Saying cumulatively means that a typical probability density function (PDF) for a given period of time can be established. Wind speed can be considered as a continuous function depending on time. However, the number of wind speed values stored by any meteorological station can only be finite, due to the existence of a sampling period, which means these values can be interpreted as a discrete function. The Weibull distribution is described by means of a continuous function of time, so when it is used as a wind speed distribution it must be understood at best as a good approximation to the previously considered discrete function. Finally, when simulations are required for any process needing those values, a new discrete function will be generated, the simulated discrete function, whose values will probably be extracted from the continuous function that has been obtained from the discrete values stored. It is a complex process. The goal of the simulations is to obtain series satisfying some statistical constraints, depending on those observed in the original series. Generally, these statistical values are given in the form of a PDF, or at least given in terms of the mean, variance and correlations between series. Occasionally, autocorrelations are also considered, especially when chronological features are required. 2.1. Weibull and Rayleigh distributions The PDF of the Weibull distribution, f , has been widely described in the literature, and can be expressed as a function of three parameters, i.e. origin (γ), scale (C) and shape (k), although in the case of wind speed distributions, given that the minimum value of wind speed is 0, the origin parameter, γ, equals 0, and this is why the following 2-parameter expression (C, k) stands for it: 8 vw o 0 <0   ð1Þ f ðvw Þ ¼ k
vw k  1 exp 
vw k v :C C w Z0 C where vw is the wind speed absolute value and exp represents the exponential function. Parameters C and k are connected with the mean, μvw , and with the variance, σ 2vw , of the distribution bymeans of the Gamma



2  , function [7], μvw ¼ C Γ 1 þ 1k , and σ 2vw ¼ C 2 Γ 1 þ 2k  Γ 1 þ 1k R1 and Γ ðpÞ ¼ 0 e  x xp  1 dx. In addition, a simplified version of the Weibull distribution has been accepted as an alternative PDF for wind speed series, i.e. the Rayleigh distribution. A Weibull distribution is said to be a Rayleigh one when parameter k equals 2. The expression of the Rayleigh PDF function does not need to be written, because it consists of substituting k in (1) by 2. An advantage of this distribution is that it can be completely determined by its mean value, and this feature can be very useful in locations where no large sets of data are available, and where only an estimate of the average wind speed is known. If the Rayleigh distribution is assumed, then the pffiffiffi mean and the variance can be expressed as μvw ¼ C 2π and σ 2vw ¼ C 2 ð1  π4Þ.

A. Feijóo, D. Villanueva / Renewable and Sustainable Energy Reviews 56 (2016) 473–483

Summarizing, both Weibull and Rayleigh distributions have been widely used to represent wind speed cumulative behaviour. When both are compared, the following aspects can be highlighted. The Weibull distribution can be defined by means of three parameters in the general case, which turn into only two in the case of wind speeds, because there are no negative values for them. When one of these two parameters (k) is made equal to 2, then the particular case of the Rayleigh distribution is obtained, and it is likely that the obtained PDF will not fit the original data as accurately as in the case of the Weibull one. Instead, the advantage of the Rayleigh distribution is that it is completely defined by the mean value of the set. In many cases only this mean value is known or can be estimated to some extent, and there is no other knowledge of wind speed features at the given site. In these cases, the use of the Rayleigh distribution can be of great value because, among other reasons, there is no criterion for using a different PDF. So, in a way, it can be said that the use of Rayleigh distribution is less accurate but easier. As has been mentioned earlier, it must not be forgotten that these two functions are continuous, while in meteorological stations only finite sets of values can be stored. 2.2. Other distributions Some authors have proposed different distributions for describing cumulative wind speed behaviour. Among them, the inverse Gaussian function presented as a useful alternative to the 3-parameter Weibull distribution [8], Normal, Gamma or Lognormal distributions [9], or different combinations such as Normal-Weibull [10,11], two Weibull distributions [10–13], two Normals and Gamma-Weibull [13] and others [11,14]. According to Bardsley [8], the inverse Gaussian distribution is a possible alternative to the three-parameter Weibull one, and is appropriate for describing wind data with low frequencies of low speeds. However, the author reveals that the fit is not always good enough, but suggests that it could be more widely used than it is. Safari [9] presented an analysis of the wind potential in Rwanda by using four different distributions: Weibull, Normal, Gamma and Lognormal. An interesting conclusion of the work is that it cannot be said there is a typical distribution representing cumulative wind behaviour. For different sites, the best fit can be achieved by using different distributions among those given, assessed by means of the χ 2 goodness-of-fit test for each site. What is also interesting is the fact that Normal, Gamma or Lognormal distributions show better fit features than the Weibull one in some cases; Gamma distribution being the one that appears as the best one in many of the cases where it is not the Weibull distribution. The idea of obtaining better fits by means of the use of different distributions other than those obtained with the Weibull one is also studied by Romero-Ternero [10]. The paper studies the relative error in the calculation of annual mean power for a given set of wind turbines (WT) in the Canary Islands and compares the behaviour by using the Weibull distribution and two others: the singly truncated from below normal Weibull mixture distribution and the two-component mixture Weibull distribution. However, the author concludes that a representation such as the Weibull one seems to be sufficient because the estimation in the generated annual mean power is acceptably accurate. Carta and Ramírez also concluded that the Weibull distribution is not the best fit in the case of several stations on the Canary Islands [11]. They investigated results obtained by means of the use of other distributions, such as singly truncated from below normal distribution, singly truncated normal Weibull mixture distribution and two component mixture Weibull distribution, in contrast with the classical Weibull distribution, and concluded that the mixture distributions are a useful alternative to the Weibull one. Carta and Mentado

475

have even used different variables, such as air density, generally considered as a constant, or the correlation between air density and wind speed variability, in order to propose a bivariate model for wind power density and WT energy output estimations [12], concluding that this model is more realistic than the univariate probability ones, although they recognise their results do not differ much from those obtained through the use of univariate models due to the climatological characteristics of the area of study. Four mixture probability density functions and a function derived with the maximum entropy principle, have been compared with the conventional Weibull one by Chang [13], who also concluded that, in certain cases, the Weibull distribution is not necessarily the best fit, depending on certain considerations described in the work. In addition, the maximum entropy principle has been investigated by Akpinar and Akpinar [14], together with the Weibull distribution and two mixture distributions. The conclusion is that the proposed mixture distributions are a good alternative to the Weibull one, and some of the stations investigated can be adequately represented by the maximum entropy principle. A comparison between parametric and nonparametric methods can be read in [15], where a kernel density method (KDM) has been used for estimating wind speed probability distributions, showing better results than those given by other parametric distribution models. A KDM was also used in [16] with the same purpose of estimating wind speed probability distributions. According to these works and others existing in the literature, it seems to be clear that Weibull or Rayleigh distributions are, of course, not the only possible models for describing wind behaviour. Finally, a review of wind speed probability distributions used in wind energy analysis can be read in the work of Carta, Ramírez and Velázquez [24].

3. Correlation The correlation between two statistical variables, x and y, is a real number in the interval [  1,1] that quantifies how both variables evolve with each other. If they change in a similar way, i.e. they increase simultaneously or they decrease simultaneously, they are said to be strongly correlated and the value will be close to 1. If both variables change in the opposite way, i.e. one of them increases while the other decreases, and the trend is always similar, then they are said to be strongly correlated in a negative sense, and in this case the value will be close to  1. Both variables are said to be uncorrelated if their correlation is close to 0. In the case of wind speed series, correlation can also be considered in simulations, because it truly exists. A simulation of possible wind speed series of values at different locations will be more realistic if correlation is taken into account, because in this case not only are the statistical parameters of every distribution considered but also their joint variability. From the point of view of EPS analysis, the above explanation is relevant because, although correlation does not always mean dependency, what is unquestionable is that dependency exists, and it is quantifiable through correlation. When a cumulative situation of wind energy generation is simulated in an EPS, the results can be different in the cases when correlation is considered or not. How to induce a given correlation to sets of several data series is the object of discussion in the following sections. Before proceeding, it may be of interest to comment on two methods that have not been studied in depth by the authors, i.e., the one based on copulas and the NATAF method.

476

A. Feijóo, D. Villanueva / Renewable and Sustainable Energy Reviews 56 (2016) 473–483

Copulas have been successfully used for modelling dependence in stochastic processes with application to power system uncertainty analysis, in cases with large-scale integration of wind power [17]. A case with both photovoltaics and wind energy in a distribution network has been investigated on the basis of the use of Archimedean copulas [18], which gives better results than methods based on linear correlation matrices. Probabilistic load flow (PLF) is specifically mentioned in [18], and is the main subject of the work by Usaola [19], where correlated wind injections are considered in the PLF, and an approximation based on statistical moments and Cornish–Fisher expansion has been employed to solve the problem of dependence, obtaining accurate results assessed by means of the estimation of statistical moments. The paper highlights the importance of correlation in the accurate evaluation of distribution functions. Hagspiel et al. have used copula theory for generating increased data sample sizes for Monte Carlo simulations in stochastic modelling of wind power in Europe and the implications for the Swiss power grid [20]. Morgenstern and Nataf models have been studied by Liu and Kiureghian [21], with the conclusion that the Morgenstern one fails when the correlation of the correlated variables is outside the interval [  0.3, 0.3], although it behaves adequately when simulating PDFs. However, the Nataf procedure can describe a wider range of correlation coefficients although the simulated PDF exhibits worse behaviour than the Morgenstern one. A model for representing random vectors of arbitrary distributions and a given correlation matrix has been studied by Cario and Nelson [22], based on the idea of transforming multivariate normal random vectors into the desired random ones, a method called NORTA (Normal to anything). The impact of wind production on the marginal prices of a pool-based electricity market has been investigated by Morales et al. [23], where the Nataf transformation has been used. 3.1. Inducing Pearson parametric correlations Pearson correlations are normalised covariances, defined as σ2

ρxy ¼ σxxyσ y , where σ 2xy is the covariance between both variables x and y, and σ x and σ y , are their corresponding standard deviations. As mentioned before, correlation and dependency can be related but do not have the same meaning, as correlation involves only a linear measure, a mathematical linear relationship, and dependency involves a relationship between a cause and an effect. Having accepted that the cumulative behaviour of wind can be described, as a good approximation, by means of certain continuous distributions, such as the Weibull one, and having also accepted the fact that wind speed series at different locations are correlated, the problem to be solved consists of generating series of values that satisfy those constraints, i.e. their cumulative behaviour can be fitted to those given by the distribution parameters, and the correlations must be close to the given correlations. As a very simple numerical example for two locations, the problem can consist of obtaining a set of 100,000 samples for each location, given the parameters C ¼7.1 for the first location and 6.5 for the second one, and given parameters k¼ 2 and 1.75 respectively, and with a correlation of 0.78. In general, such an example should not be difficult to solve. The real problem generally involves a higher number of locations, and the difficulty grows with this number. In this section some of the solutions that have been proposed are discussed. Before explaining these methods, some other previous works can be mentioned. For instance, the problem of multivariate statistical simulation was applied to non-Normal distributions by Johnson [25], while Lee [26] carried out a

similar investigation for distributions with Weibull properties and Jensen [27] focused on the analysis of Rayleigh distributions. Correlation has been simulated in wind speed series in different ways [28]. The one reviewed here was presented by Feijóo, Cidrás and Dornelas [29]. Given a set of M uncorrelated Normally distributed variables, with means equal to 0 and variances equal to m mN 1, let us say, xm M m ¼ 1 , where for all m ¼1,…,M, x ¼ xn n ¼ 1 , i.e. each variable xm consists of N samples. Then, a given correlation matrix Ω can be induced over such a set, converting it into a new set m ym M m ¼ 1 , where every variable y , for all m ¼1,…,M is also Normally distributed. The means and variances are retained and the correlation matrix is as close to Ω as desired, at least theoretically [30]. The process is very simple, as it consists of operating as follows: y ¼ Lx

ð2Þ

where x is now an M  N matrix containing the n samples of each of the xm variables, y is an M  N matrix containing the n samples of each of the ym variables and L is a lower triangular matrix obtained from the correlation matrix Ω by using the Cholesky technique, i.e. Ω ¼ LLT [25,31], where Ω is an M  M matrix. Specifying it a bit more, (2) can be rewritten as in (3). 0 1 1 0 10 1 1 x1 ⋯ x1N y1 ⋯ y1N L11 ⋯ 0 B ⋮ ⋱ ⋮ C B ⋮ C B ⋱ ⋮ A@ ⋮ ⋱ ⋮ C ð3Þ @ A¼@ A M M M ⋯ x y1 ⋯ yN xM LM1 ⋯ LMM N 1 If the mean values of the desired variables are different from 0, and they are given in an M  1 matrix μ, then (2) becomes (4). y ¼ Lx þ μ

ð4Þ

The method explained above works correctly for Normal distributions. However, it does not work in a such an accurate way for non-Normal distributions. In the paper by Feijóo, Cidrás and Dornelas [29] the method was applied to Rayleigh distributions generated using the Monte Carlo method (MC) [32–34], by converting Uniform distributions into Rayleigh ones by means of the inverses of their cumulative distributions functions (CDF), although commercial software such as Matlab generates Rayleigh distributions and others directly. So xm M m ¼ 1 was not a set of Normal distributions, but it was a set of M Rayleigh distributions. The consequences of applying it were: 1. Generation of negative values for the wind speed. 2. Loss of Rayleigh nature of the distributions, which also occurs with other distributions, such as Weibull ones. Negative values are a problem because there are no negative wind speed values and so they have to be rejected. In general, fewer than 5% of the sets of simulated values contain negative wind speeds, although if the number of locations grows, this figure also grows. The solution employed, the rejection of the sets with a negative value, involves a new error in mean and covariance, as will be shown in an example later. On the other hand, from a computational point of view the method operates in a sequential way (beginning with distribution 1 and going up to distribution M, and taking into account values of the m  1 previous distributions), and the consequence of this is that the higher the index m, the more the distributions lose their Rayleigh characteristics. As mentioned, the method operates adequately for Normal distributions, but with errors for nonNormal ones. 3.2. The use of evolutionary algorithms Villanueva and Feijóo proposed a completely different method consisting of the use of genetic algorithms (GA) and local search technique (LS) to simulate correlated wind speed series [35].

A. Feijóo, D. Villanueva / Renewable and Sustainable Energy Reviews 56 (2016) 473–483

The method is based on generating series of wind speeds, considering the distribution at each location, and then to rearrange each of them in order to obtain the correlation matrix requested. The values of wind speed generated initially will not change, but their positions in the series may be different. Rearrangement of the sets means obtaining a different permutation of their values. This method uses evolutionary algorithms such as GA and LS. GA has been explained by Gen and Cheng [36] and LS by Hoos and Stützle [37]. It can be said that the objective of evolutionary techniques is to minimise the value of an error function, obtained from the comparison of a possible solution and the requested one. The drawbacks of using evolutionary techniques are the long computing time and the need for a proper initialisation. LS and GA are very common and complementary. LS is based on the fact that any possible solution has better ones around it, except the best one. Taking this premise into account, if a possible solution is slightly modified, but keeps certain conditions, the solution obtained may be better than the original one. In that case, the new solution is accepted whereas otherwise it is discarded. GA consists of using two possible solutions, combining them to obtain a pair of new ones. This technique may speed up the process to obtain the best solution if one of the new solutions is better than the former ones. Both techniques also accept, randomly and with low probability, worse solutions, in order to avoid being caught in a local minimum. When applied to obtain correlated series of wind speed values, LS consists of swapping pairs of values in the same series, with the objective of obtaining a better correlation matrix. GA is based on obtaining two new solutions that are the result of a random combination of the locations from another two already accepted as possible ones, having the same objective. The functions to be minimised, defined as the proximity to the desired results, are obtained through the maximum norm of the differences between the elements of the correlation matrix, obtained from the possible solution, and those of the desired one as in (5). n o k d ¼ max ρij  ρkij 1 r i; j r n ð5Þ When the acceptable error is obtained or a maximum number of iterations is reached, the best of the current possible solutions is considered the best one. By using the method presented in this section, and in comparison to the method explained in 3.1, the following conclusions can be given: 1. No negative values of wind speeds are generated. The resulting data series are permutations of the initial ones. 2. For the same reason the statistical values of the resulting distributions, mean and variance and all the distribution parameters, i.e. the nature of the distributions themselves, are retained. However, it must not be forgotten that the process is very timeconsuming, which seems to be its main disadvantage. 3.3. Inducing Spearman rank correlations Feijóo and Sobolewski proposed a new method [38], based on the same idea of correlation as in [29] but with a significant difference. The correlations considered are rank correlations, more specifically, Spearman rank correlations. Spearman rank correlations can be defined for two different series of data, each one

477

containing N values, as in (6). r ¼ 1

6

PN

i¼1 2

2

di

NðN  1Þ

ð6Þ

where di are the differences in rank between samples of both series that occur simultaneously, after having ranked both series previously. The main difference between this method and the one explained based on Pearson parametric correlations is that this one does not propose to calculate a new set of correlated values, mM ym M m ¼ 1 , starting with uncorrelated values x m ¼ 1 . It is based on operating permutations of the values inside the variables, based on an idea proposed by Iman and Conover [39]. In this sense, it can be said to operate in a similar way to the evolutionary methods, but consuming less time. To facilitate its understanding, in the case of M¼ 2, i.e. when there are only two sets of values, with a given Spearman rank correlation, the method proposes to operate a permutation of x2 , x1 and x2 resulting in the desired rank correlation. If x1 ¼ x11 ; x12 ; …; x1N , and x2 ¼ x21 ; x22 ; …; x2N , for a given value of N ¼ 5, let us suppose x1 ¼ 1; 2; 6; 4; 5 and x2 ¼ 3; 4; 3; 5; 8, the Spearman rank correlation between both series is 0.2050. However, if a permutation of the second one is operated, for instance π ðx2 Þ ¼ 8; 3; 4; 5; 3, then the Spearman rank correlation between x1 and π ðx2 Þ is  0.4104, but the values inside the series have not changed. The new series is a permutation of the initial one, so the mean value, the variance and the PDF are exactly the same. As the method performs permutations of the values, it operates in a similar way to that presented in 3.2, but with two differences, previously mentioned: 1. It operates by trying to fit Spearman rank correlations, while the GA based method works with Pearson parametric correlations. 2. It is much less time-consuming, as it is based on algebraic operations, and not on GA. This method, additionally, avoids problems that had appeared in the method explained in 3.1: 1. No negative values are generated because there are no negative values in the initial sets, and no new values are generated. Only permutations of initial sets are handled. 2. The nature of the initial distributions is completely retained for the same reason, with the only chronological difference derived from the permutations operated. All the parameters of the distributions are retained. If the retaining of Pearson parametric correlation is very important, then the method can be employed assuming the error committed. 3.4. A different method for inducing Pearson correlations Villanueva, Feijóo and Pazos [40] presented a modification of the method proposed in [29]. The goal of the work is to induce correlations between data series that must satisfy being part of Weibull distributions, again according to Pearson correlations. The way of operation is similar to the one presented in [29], but instead of operating directly with the Weibull distributions (in [29] it was made with Rayleigh distributions), correlations were induced over a set of Normally distributed variables and then a transformation from Normal to Weibull distributions was operated.

478

A. Feijóo, D. Villanueva / Renewable and Sustainable Energy Reviews 56 (2016) 473–483

Assuming a set of Normally distributed values xn N n ¼ 1 belong to the domain of a real function such as the following one:  111 0 0 k 1  erf pxnffiffi2 @ @ A A un ¼ C  log ð7Þ 2 where C and k are the constants of the desired Weibull distribution, log stands for natural logarithm, and erf is the error function defined as: Z x 2 2 e  t dt ð8Þ erf ðxÞ ¼ pffiffiffiffi

π

0

then it can be said that un N n ¼ 1 constitutes a set of data that are distributed according to a Weibull distribution with parameters C and k as desired. The proposal made in [40] consists of inducing the desired correlations in a set of Normal distributions, with means equal to 0 and variances equal to 1 and then transforming them by means of (7) into Weibull distributions with the desired parameters. The reason why correlations are retained to a high degree when data series are transformed from Normal to Weibull distributions has been studied by Feijóo and Villanueva [32]. But as a conclusion it has to be said that, in comparison with method explained in Section 3.1, here: 1. No negative values are generated, i.e. all the values generated are valid. 2. The nature of the resulting distributions fits better to the real ones. This will be shown with an example. What ensures the nonexistence of negative values is the use of (7), an equation that converts values of the domain of a Normal distribution into their corresponding values of a Weibull distribution, considering that the image of both distributions, i.e. what is generally considered as the CDF, is common to both. 3.5. Example The methods proposed in Sections 3.1–3.3 have been illustrated in the cited works. However, this is not the case of the method proposed in Section 3.4, presented in [40], where no results of the wind speed series simulated were given. In that work the wind speed series were used with a different purpose. The method described in Section 3.4 shows its true advantage in comparison to that reviewed in Section 3.1. In this section the first example given in [29] has been retrieved in order to compare results obtained by means of both methods. The example consisted of seven WF for which Rayleigh distributions were assumed, with the constants C given in Table 1. The seven WF are numbered from 1 to 7. In Table 1 the following information has been given for each WF: constant C of the Rayleigh distribution (C), mean wind speed in ms  1 (vwmean ), and as a result, constants C obtained by means of methods presented in 3.1 (I) and in 3.4 (II), and constants k and mean values obtained. The Table 1 Comparison of results of methods explained in Sections 3.1 and 3.4. WF

C

vwmean

CðIÞ

CðIIÞ

kðIÞ

kðIIÞ

vwmean ðIÞ

vwmean ðIIÞ

1 2 3 4 5 6 7

6.54 6.77 6.99 5.98 6.54 7.22 10.49

5.8 6.0 6.2 5.3 5.8 6.4 9.3

6.65 6.91 7.14 6.10 6.66 7.35 10.67

6.54 6.76 6.98 5.98 6.55 7.23 10.51

2.08 2.11 2.11 2.10 2.07 2.08 2.08

2.01 2.01 2.01 2.00 2.00 2.01 2.01

5.89 6.13 6.34 5.41 5.89 6.52 9.47

5.79 6.00 6.19 5.30 5.80 6.41 9.31

values of these four constants (C and k by both methods) correspond to maximum likelihood estimations with a 95% confidence. These results have been obtained with a simulation of 100,000 samples, where 3794 (3.794%) produced some negative wind speed data with the first method. From additional simulations it could be stated that the samples with some negative value of wind speeds are generally under 4% of the total. Results given in Table 1 reveal that the method explained in Section 3.4 improves results given by the method presented in Section 3.1, at least in relation to the parameters of the marginal distributions. As for retaining the correlation matrix values, which have not been included here, it has to be said that both methods achieve a good approximation to the desired one. If Ω denotes the correlation matrix of the initial series and ΩI and ΩII the correlation matrices obtained by means of both methods respectively, and if the error matrices EI ¼ Ω  ΩI and EII ¼ Ω  ΩII are defined, then the Euclidean norms of both error matrices are jEI j ¼ 0:107 and jEII j ¼ 0:106 respectively, taking into account the fact that they are 7  7 matrices, it can be considered to be very good approach in both cases.

4. Autocorrelation Autocorrelation is correlation of a series with itself after applying a given lag. Lag one autocorrelation is the correlation of a series and the same series lagged by one position, and in general lag n autocorrelation is the correlation of a series with the same series lagged by n positions. Autocorrelation in wind speed series has been studied and simulated, among others, in [1, 40–47]. Autocorrelation becomes important when chronology needs to be taken into account. Given a data series with N different data, xn N n ¼ 1 , autocorrelations can be different for all the permutations of the series. Coming back to numerical examples, let us suppose the series is given by the five values f1; 4; 3; 5; 7g: In this case lag 1 autocorrelation equals 0.1, lag 2 autocorrelation equals 0, lag 3 autocorrelation equals  0.15 and lag 4 autocorrelation equals  0.45. A permutation of this series, for example 1; 5; 3; 4; 7, presents different autocorrelation values, i.e.  0.2, 0, 0.15 and  0.45 respectively. As autocorrelation seems to be important when chronology must be considered, it will be needed in processes such as those involving the need for prediction or in dynamic processes. For example, in probabilistic load flow (PLF) autocorrelation will not be needed because the order of the samples in time is not important, but it will be needed if the injected power in an EPS needs to be scheduled, such as occurs in the electricity market in some countries, where the generation capacity has to be predicted up to 24 h in advance. In general, autocorrelation can be induced in a series of data by means of autoregressive models [48], i.e. AR, MA, ARMA, ARIMA [49,50], models that are mentioned here, although the focus will be on a different group of models, such as Markov processes, in the next section.

5. Markov processes Markov processes can be considered as autoregressive processes. Markov processes have been used for wind speed simulation in several works, in their first-order Markov chain version [51,52] or in the second-order one [53]. Nfaoui, Essiarab and Sayigh investigated the preservation of statistical characteristics of the distributions [54] when Markov processes are used, and

A. Feijóo, D. Villanueva / Renewable and Sustainable Energy Reviews 56 (2016) 473–483

479

The simulation of wind speed series with the help of a Markov process uses the MC method. The first value of the series can be randomly generated inside the total interval of wind speeds. After that, all the values can be obtained by generating a random number between 0 and 1 in a Uniform distribution, and converting it into a wind speed value with the help of the transition matrix, if every line of such a matrix is interpreted as a PDF. One of the consequences of using this method is that autocorrelations are highly retained, although they are more accurate as their lag number is lower. Lag 1 autocorrelation is more accurate than lag 2 autocorrelation, and so on, as will be seen later in an example. And another consequence is that the nature of the resulting PDF is, in general, very accurate in comparison to those of the original ones, which will also be illustrated with an example.

Brokish and Kirgley studied the pitfalls of using these methodologies [55]. In this section some results obtained from the authors' own experience will be explained. On this occasion, samples from real meteorological stations have been obtained. The sample data were published by Meteo Galicia through its web page (available at www.meteogalicia.es), a meteorological research institute located in Galicia. This region is an autonomous community in the northwest of Spain where the use of wind energy for electricity supply has seen significant development since the 1990s. 5.1. First-order Markov processes The theory relating how to use Markov chains for the simulation of wind speed series is known and has been presented in works such as those mentioned above. A brief summary of some ideas is presented here. If the interval of possible wind speeds at a location is divided into smaller intervals, for instance, if the interval [0, 25] is divided into intervals of, let us say, 1 ms  1 , [0, 1), [1, 2),…, [24,25], a matrix can be generated accounting for the probabilities that wind speed goes from one interval to a different one. These probabilities can be associated to the relative frequencies of such jumps in observed series. Such a matrix is called a Markov transition matrix. If such a matrix is at a researcher's disposal, it can be used for simulating wind speeds in a chronological way. Summarizing, a Markov matrix is a mathematical object similar to: 0 1 P 11 P 12 … P 1N B P … P 2N C B 21 P 22 C B C B ⋮ ⋮ ⋱ ⋮ C @ A P N1 P N2 … P NN

5.2. Second-order Markov processes In second-order Markov processes, as explained by Shamshad et al. [53], the matrices contain probabilities of the variable, in this case wind speed, to be in a given state, with not one, but two known past states. A second-order Markov matrix is, actually, a set of matrices that can be represented such as proposed by Shamshad et al. [53]: 0 1 … P 11N P 111 P 112 B P P 122 … P 12N C B C 121 B C B ⋮ C ⋮ ⋮ ⋱ B C B P … P 212 P 21N C B 211 C B C B P 221 P 222 P 22N C … B C B ⋮ ⋮ ⋮ C ⋱ B C @ A P NN2 P NN1 … P NNN

where a value P ij indicates the probability of going from state i to state j, and N possible states are considered. If the intervals [0, 1), [1, 2),…, [24,25] are considered, then state 1 corresponds to the interval [0, 1), …, and state 25 corresponds to [24,25]. It must be noticed that an NN matrix has to be stored. For example, in the case of 25 states, this matrix has 625 elements, although experience shows that not all these values are different from 0 and, in practice, there are many zeros in the matrix. In the Markov matrix, the values of the lines must satisfy that PN j ¼ 1 P ij ¼ 1, for all i A 1; 2; …; N. More about the formation of Markov matrices can be read in [53]. For example, if some of the Meteo Galicia stations are analysed and the first-step Markov matrix is generated, the percentage of nonzero elements in the resulting matrices are those given in Table 2 (ten meteorological stations in the province of Pontevedra, which will be denominated as stations 1–10). In all cases the initial sample size was of 52,560 and they are 10-min mean wind speed values through one year. The reason why there are meteorological stations in Table 2 with Markov matrix sizes greater than 25 is because some of them have registered extreme values during the sampling period. For example, in station 6, the maximum value of wind speed was 35.5, which makes the program developed by the authors automatically generate intervals from [1, 2) to [35,36], and to generate a 3636 matrix.

where P ijk indicates the probability of experiencing a transition from i to j and then from j to k. In a similar way to the case of first-order Markov matrices, PN PN j¼1 k ¼ 1 P ijk ¼ 1, for all i A 1; 2; …; N. One of the disadvantages that can be foreseen is the fact that these matrices are much greater than first-order ones. For instance, in the case of 25 intervals, a first order matrix has 625 elements, although not all of them are different from 0, but in the case of a second-order one, the number of elements to be stored would be 15,625. Researchers must therefore assess whether it is worth generating such a matrix instead of using the original series directly. It must be noticed that not all elements are different from 0, and also the fact that original chronological series are generally much larger. An analysis of the same ten meteorological stations of Meteo Galicia has been carried out, for the cases of first-order and second-order Markov matrices, with the results given in Table 2. The number of intervals needed for each of the stations is given along with the percentage of values different from 0 in both cases. Notice that a number of intervals of 23 means matrices of 23  23 in the first-order case and 23  23  23 in the second-order one. The use of second-order Markov processes results in an improvement in the autocorrelations and also in the retention of the nature of the PDF with respect to the original ones.

Table 2 Percentage of non-zero elements of the Markov matrix. Station

1

2

3

4

5

6

7

8

9

10

Intervals a 0ð1stÞ a 0ð2ndÞ

22 37.60 8.38

15 45.33 12.24

12 41.67 9.20

34 27.71 2.75

21 33.10 6.88

14 36.22 8.67

25 31.84 5.54

28 25.89 3.90

12 51.39 16.55

14 42.35 11.70

480

A. Feijóo, D. Villanueva / Renewable and Sustainable Energy Reviews 56 (2016) 473–483

Table 3 Comparison between real data PDF and simulated series PDF. Station

1 2 3 4 5 6 7 8 9 10

Original data

First-order

Second-order

AR(1)

C

k

C

k

C

k

C

k

5.30 2.09 1.40 5.70 4.16 1.32 4.18 6.85 2.09 3.10

1.60 1.13 1.37 1.57 1.35 1.01 1.47 1.96 1.45 1.74

5.32 2.17 1.45 5.96 4.12 1.39 4.16 6.89 2.06 3.15

1.67 1.29 1.46 1.51 1.32 1.18 1.45 2.01 1.40 1.68

5.20 2.25 1.41 5.86 4.24 1.38 4.07 6.71 2.09 3.04

1.68 1.29 1.47 1.59 1.34 1.19 1.46 2.03 1.40 1.75

5.31 2.04 1.39 5.85 4.07 1.32 4.24 6.71 2.07 3.12

1.63 1.10 1.36 1.59 1.34 1.01 1.46 1.91 1.43 1.77

1.1 Original First-order

1

Second-order Autocorrelation

0.9

AR(1)

0.8 0.7 0.6 0.5 0.4 0

1

2

3

4

5

6

7

8

9 10 11 12 13 14 15 16 17 18 19 20 lag

Fig. 1. : Autocorrelations (lag 1 to lag 20) for the original series and first-order Markov, second-order Markov and AR(1) simulated series for station 4 of Table 3.

5.3. Examples In the following example both Markov methods are used to simulate samples of wind speeds, and the results obtained in a simulation of both methods are compared with the results obtained by means of an AR(1) method. The simulations are based on the MC method, for which different runs would reveal slightly different results. This is why results given in Table 3 are only illustrative of the order of the values resulting by means of each method, in the sense that they have been collected in one of these runs. Table 3 does not show important quantitative differences between first-order and second-order processes. Both methods seem to achieve a good approximation resulting in an acceptable retention of the PDF nature. For an adequate interpretation of Table 3, notice that only the maximum likelihood estimators of parameters C and k of the Weibull distributions of the original series and both first-order and second-order simulated ones have been reflected in it. The values have been obtained for the same ten stations analyzed before and some of them seem to be rather unusual for wind speed stations. It must be remembered that they are meteorological stations, not specifically WFs. From Table 3 it cannot be clearly concluded that a second-order method is much more accurate than a first-order one. The advantage of the second-order method consists of improving the values of the different autocorrelations. This means that the series obtained by means of this method, from a chronological point of view, look more like the original series than those obtained by using the first-order method. As an example, and with the aim of not overloading the results with many numbers, in Fig. 1 the values of autocorrelations from lag 1 to lag 20 are shown for some of the stations of Table 3. In the same Table, the values for

lags 1 to 20 autocorrelations have also been included given by an AR(1) process such as that proposed in [40]. If similar graphs are generated for the different stations, similar results are obtained, and so it does not make any sense to repeat them here. Shamshad et al. [53] have produced graphs with much greater lag values. The interest of Fig. 1 here is only to show that it is the different autocorrelation lags that give a second-order method some advantage over a first-order one. Results of other simulations have always been similar to these, and reveal that autocorrelations of an AR(1) method lie between those given by first-order and second-order Markov ones. 5.4. Discussion This section contains a brief discussion about the Markov and AR (1) methods. The question here is, do some of the presented methods have any advantages or disadvantages that make them more recommendable? Glancing at the results of the simulation presented in Section 5.3., which can be seen in Table 3 and in Fig. 1, the three methods seem to be very similar. Looking at the features of the PDFs obtained, there are very small differences which perhaps can be said to be in favour of the AR(1) method. The approximation of the PDF distributions is better (both parameters seem to be a bit closer to the original ones), although it cannot be said that there is a great difference when this happens. In Fig. 1 it can be seen that the autocorrelation coefficients are also closer to those belonging to the original sample in the case of the AR (1) method. Although in Fig. 1 only the case of station 4 can be seen, the other stations are similar. According to these observations, AR (1) seems to have all the advantages with respect to Markov methods. In addition, Markov methods involve the use of Markov matrices, instead of a simple equation, which is the case for the AR (1) method. This seems to strengthen the previous comments. Relating both Markov methods, it is not clear if the improvements obtained by means of the two-step one justify the effort of programming it. Improvements can be seen in Fig. 1. All these methods can be used in the same cases, although as commented before, Markov methods need the use of Markov matrices. So if only a lag 1 autocorrelation value is known, together with the feature of the PDF, then the AR(1) method is the recommended one. 5.5. Correlation in Markov processes The main disadvantage of Markov processes, shared with other autoregressive techniques, is the difficulty of inducing desired correlation in simultaneous series. An interesting discussion could be introduced here, focussing on the possibility of simulating wind speed series while keeping statistical nature according to given parameters, and also respecting given correlation matrices and autocorrelations. The methods explained for inducing correlations work with some success, but autocorrelations are not guaranteed. And autoregressive or Markov methods operate adequately with autocorrelations but can fail when inducing correlations. However, given series of simultaneous wind speeds, all these statistical values are present, so it should be possible to simulate them. As an example of this difficulty, stations 1, 7 and 8 of Table 3 have been chosen, as they have better features for possible WF installation. The correlation matrix of the original data of these stations is the following: 0 1 1 0:6627 0:7397 B 0:6627 1 0:6579 C @ A 0:7397 0:6579 1 where station 1 is at position 1, station 7 is at position 2 and station 8 is at position 3.

A. Feijóo, D. Villanueva / Renewable and Sustainable Energy Reviews 56 (2016) 473–483

The use of a first-order Markov process for inducing autocorrelations to the final series gives the results expressed before in Table 3. However, the correlation matrix obtained is the following: 0 1 1 0:0116 0:0262 B C 1 0:0086 A @ 0:0116 0:0262 0:0086 1 A second-order Markov process gives the following correlation matrix: 0 1 1 0:0121  0:0196 B 0:0121 1  0:0028 C @ A  0:0196  0:0028 1 which does not reveal a better result than a first-order Markov method. A challenge for future work would be to generate series keeping all the statistical features, i.e. nature of the distributions, correlations and autocorrelations, and also introducing cross-correlations, not commented in this paper. Perhaps evolutionary algorithms or other different algorithms could be appropriate for this purpose.

6. Some comments about other methods Finally, and very briefly, in this section a few comments are added about other methods proposed in the literature that have not been so widely explored by the authors in this paper. Among them, Sfetsos presented a comparison between different forecasting techniques for hourly mean wind speed time series [56]. The fact that wind power generation differs from other conventional generation technologies has been made clear in [57], which reviews some methods that deal with wind power generation forecasting. As for forecast techniques, different papers can be found in the literature on Kalman Filter (KF) methods. One of them was presented by Bossanyi [58], with the idea of predicting wind speeds a short time ahead, with good results for the prediction of oneminute-average wind speeds, but not so good results for hourly wind speed predictions. These results were satisfactorily applied to increase energy yield in a WT. The use of the KF for short-term wind speed prediction was also an issue in [59], where the authors constructed a support vector regression model and formulated a nonlinear state-space framework, using KF for updating states under random uncertainty. According to their conclusions, the method can predict wind speed in the short-term with higher accuracy than other conventional methods, such as those based on artificial neural networks (ANN) or autoregressive ones. Zuluaga et al. [60] also worked on wind speed prediction with the help of KF, in an operation based on statistical indicators such as skewness and kurtosis. In [60] three different methods were experimentally compared, for making KF robust for one-step-ahead wind speed predictions, and some possible improvements to the standard KF were presented. Also with the aim of contributing to short-term wind speed forecasting, in [61] a hybrid algorithm was proposed for wind forecasting by means of the combination of time series analysis and a KF algorithm. The results showed an improvement in prediction when the proposed algorithm is used, compensating prediction errors and prediction time delay, although prediction results were affected by some errors in the case of fluctuating wind speeds. Louka et al. used a technique based on third-order polynomial functions in Kalman filtering with accurate results in wind speed numerical predictions [62], with direct application to a WF in Crete, Greece. The method was applied to data given by atmospheric numerical models with different capabilities in

481

horizontal resolution. Data filtered by the KF improved the possibility of forecast power prediction. Cassola and Burlando [63] focused on numerical weather prediction, and used KF for finding suitable configurations for wind prediction, the filter being used in that work for forecasting wind speed at a given WF. According to the conclusions of the authors, simulations show very low errors with respect to measurements. Crochet [64] presented an adaptive KF procedure applied to 10-m wind speed forecasts, developed with the goal of reducing systematic bias and improving the accuracy of local forecasts from a numerical weather prediction model. It was performed by adding two algorithms to the traditional KF procedure, with the result of removing systematic errors in the prediction. There are many prediction tools for wind speed, and Poncela et al. worked on improving the performance of statistical tools by means of KF models [65]. They estimated model parameters by quasi maximum likelihood methods, producing an automatic self-tuning of model parameters for different WFs, optimising model parameters and obtaining improvements in the results. A review on the history of short-term prediction of wind power can be read in [66]. Bayesian models and artificial networks have also been researched [67–69]. The combination of Bayesian model averaging and Markov Chain Monte Carlo sampling methods has been applied to derive wind speed distributions in a robust way [67]. A method based on Bayesian combination algorithm has been used for improving wind speed forecasting results obtained by means of three different neural network (NN) models [68], and these three different NN models, denoted as adaptive linear element, back propagation and radial basis function have been compared in [69]. Wavelet transforms have been used for wind speed simulation purposes [70,71], although not always for energy calculation [72]. The wavelet transform has been combined with the auto regressive moving average (ARMA) model to give minimum absolute prediction error (MAPE) in wind speed forecasting techniques [70], and results have been compared with those obtained by means of a combination of artificial neural network and ensemble Kalman Filter resulting in improvements in the forecasting error. Based on wavelets and also in classical time series analysis, a forecasting method has been proposed with small error in multi-step forecasting, which can be applied to wind speed and wind power forecasting [71]. There are many more papers on the interesting subject of simultaneous wind speed series simulation. In this paper an overview of some of them has been presented.

7. Conclusions Several methods for the simulation of wind speeds have been reviewed. The goal of the simulation of wind speeds in WF is to obtain series of data that satisfy certain requirements, such as fitting wind speed PDF and correlations between series. Depending on the kind of processes that must be simulated, the needs of wind speed series can include chronological features or not. When no chronological features are needed, the methods presented in Section 3 are valid. The main conclusion is that, from those methods, the one presented in Section 3.4 achieves the best results, because it retains statistical values of the initial marginal distributions and correlation with great accuracy. When chronological requirements are part of the constraints, the autoregressive methods explained in Section 5 are more advisable. A comparison between an AR(1) and two Markov based methods has been presented, where results do not reveal very important differences between them. Simultaneously inducing correlation and autocorrelation does not seem to be so easy, and this can be part of the challenge of

482

A. Feijóo, D. Villanueva / Renewable and Sustainable Energy Reviews 56 (2016) 473–483

improving methods in the future, together with the possible consideration of cross-correlations.

References [1] Brown BG, Katz RW, Murphi AH. Time series models to simulate and forecast wind speed and wind power. J Clim Appl Meteorol 1984;23:1184–95. [2] Hennessey Jr. JP. Some aspects of wind power statistics. J Appl Meteorol 1977;16:119–28. [3] REE. The Spanish Electricity System 2013. Available at 〈http://www.ree.es〉. [4] Castino F, Festa R, Ratto CF. Stochastic modelling of wind velocities time series. J Wind Eng Ind Aerodyn 1998;74–76:141–51. [5] Freris LL. Wind energy conversion systems. Cambridge, UK: Prentice Hall; 1990. [6] Arslan T, Murat Bulut Y, Yavuz Arzu Altin. Comparative study of numerical methods for determining Weibull parameters for wind energy potential. Renew Sustain Energy Rev 2014;40:820–5. [7] Polyanin AD, Manzhirov AV. Handbook of mathematics for engineers and scientists. Boca Raton, Fl, US: Taylor and Francis Group; 2007. [8] Bardsley WE. Note on the use of the inverse Gaussian distribution for wind energy applications. J Appl Meteorol 1980;19:1126–30. [9] Safari B. Modeling wind speed and wind power distributions in Rwanda. Renew Sustain Energy Rev 2011;15(2):925–35. [10] Romero-Ternero V. Influence of the fitted probability distribution type on the annual mean power generated by wind turbines. A case study at the Canary Islands. Energy Convers Manag 2008;49(8):2047–54. [11] Carta JA, Ramírez P. Use of finite misture distribution models in the analysis of wind energy in the Canarian archipelago. Energy Convers Manag 2007;48 (1):281–91. [12] Carta JA, Mentado D. A continuous bivariate model for wind power density and wind turbine energy output estimations. Energy Convers Manag 2007;48 (2):420–32. [13] Chang TP. Estimation of wind energy potential using different probability density functions. Appl Energy 2011;88(5):1848–56. [14] Akpinar S, Akpinar EK. Estimation of wind energy potential using finite misture distribution models. Energy Convers Manag 2009;50(4):877–84. [15] Qin Z, Li W, Xiong X. Estimating wind speed probability distribution using kernel density method. Electr Power Syst Res 2011;81(12):2139–46. [16] Xu X, Yan Z, Xu S. Estimating wind speed probability distribution by diffusionbased kernel density method. Electr Power Syst Res 2015;121:28–37. [17] Papaefthymiou G, Kurowicka D. Using copulas for modeling stochastic dependence in power system uncertainty analysis. IEEE Trans Power Syst 2009;24(1):40–8. [18] Haghi HV, Bina MT, Golkar MA, Moghaddas-Tafreshi SM. Using copulas for analysis of large datasets in renewable distributed generation: PV and wind power integration in Iran. Renew Energy 2010;35:1991–2000. [19] Usaola J. Probabilistic load flow with correlated wind power injections. Electr Power Syst Res 2010;80:528–36. [20] Hagspiel S, Papaemannouil A, Schmid M, Andersson G. Copula-based modeling of stochastic wind power in Europe and implications for the Swiss power grid. Appl Energy 2012;96:33–44. [21] Liu P, Kiureghian A. Multivariate distribution models with prescribed marginals and covariances. Probab Eng Mech 1986;1(2):105–12. [22] Cario MC, Nelson BL. Modeling and generating random vectors with arbitrary marginal distributions and correlation matrix. Technical report. Department of industrial engineering and management sciences, Northwestern University, Evanston, Illinois; 1997. [23] Morales JM, Conejo AJ, Pérez-Ruiz J. Simulating the impact of wind production on locational marginal prices. IEEE Trans Power Syst 2011;26(2):820–8. [24] Carta JA, Ramírez P, Velázquez S. A review of wind speed probability distributions used in wind energy analysis. Case studies in the Canary Islands. Renew Sustain Energy Rev 2009;13:933–55. [25] Johnson ME. Multivariate statistical simulation. New York, NY, US: Wiley and Sons; 1987. [26] Lee L. Multivariate distributions having Weibull properties. J Multivar Anal 1979;9:267–77. [27] Jensen DR. A generalization of the multivariate Rayleigh distribution. Shankhya: Indian J Stat Ser A 1970;32:193–208. [28] Feijóo A, Villanueva D, Pazos JL, Sobolewski R. Simulation of correlated wind speeds: a review. Renew Sustain Energy Rev 2011;15:2826–32. [29] Feijóo AE, Cidrás J, Dornelas JLG. Wind speed simulation in wind farms for steady-state security assessment of electrical power systems. IEEE Trans Energy Convers 1999;14(4):1582–8. [30] Kleijnen J, van Groendall W. Simulation. A statistical perspective. New York, NY, US: Wiley and Sons; 1992. [31] Naylor TH, Balintfy JL, Burdick DS, Kong C. Computer simulation techniques. New York, NY, US: Wiley and Sons; 1966. [32] Metropolis N, Ulam S. The Monte Carlo method. J Am Stat Assoc 1949;44 (247):335–41. [33] Rubinstein RY. Simulation and the Monte Carlo method. New York, NY, US: Wiley Interscience; 2008. [34] Gentle JE. Random number generation and Monte Carlo methods. New York, NY, US: Springer; 2005.

[35] Villanueva D, Feijóo A. A genetic algorithm for the simulation of correlated wind speeds. Int J Electr Energy Syst 2009;1(2):107–12. [36] Gen M, Cheng R. Genetic algorithms and engineering. NY: John Wiley and Sons; 2000. [37] Hoos HH, Stützle T. Stochastich local search: foundations and applications. Morgan Kauffmann/Elsevier; 2004. [38] Feijóo A, Sobolewski R. Simulation of correlated wind speeds. Int J Electr Energy Syst 2009;1(2):99–106. [39] Iman RL, Conover WJ. A distribution-free approach to inducing ranc correlation among input variables. Commun Stat Simul Comput 1982;11:311–32. [40] Villanueva D, Feijóo A, Pazos JL. Simulation of correlated wind speed data for economic dispatch evaluation. IEEE Trans Sustain Energy 2012;3(1):142–9. [41] Feijóo A, Villanueva D. Polynomial approximation of the Normal to Weibull distribution tranformation. AIMS Energy 2014;2(4):342–58. [42] Chou KC, Corotis RB. Simulation of hourly wind speed and array wind power. Solar Energy 1981;26(3):199–212. [43] Goh TN, Nathan GK. A statistical methodology for study of wind characteristics from a close array of stations. Wind Eng 1979;3:197–206. [44] Brett AC, Tuller SE. The autocorrelation of hourly wind speed observations. J Appl Meteorol Climatol 1991;30(6):823–33. [45] Mohandes MA, Rehman S, Hlawani TO. A neural networks approach for wind speed prediction. Renew Energy 1998;13(3):345–54. [46] Kamal L, Hafri YZ. Time series models to simulate and forecast hourly averaged wind speed in Quetta, Pakistan. Sol Energy 1997;61(1):23–32. [47] Kavasseri RG, Seetharaman K. Day-ahead wind speed forecasting using fARIMA models. Renew Energy 2009;34:1388–93. [48] Madsen H. Time series anaylsis. Boca Raton, Fl, US: Taylor and Francis Group; 2007. [49] Kashyap RL. Optimal choice of AR and MA parts in autoregressive moving average models. IEEE Trans Pattern Anal Mach Intell 1982;PAMI–4(2):99–104. [50] Zhang GP. Time series forecasting using a hybrid ARIMA and neural network model. Neurocomputing 2003;50:159–75. [51] Sahin AD, Sen Z. First-order Markov chain approach to wind speed modelling. J Wind Eng Ind Aerodyn 2001;89:263–9. [52] Ettoumi FY, Sauvageot H, Adane AEH. Statistical bivariate modelling of wind using first-order Markov chain and Weibull distribution. Renew Energy 2003;28:1787–802. [53] Shamshad A, Bawadi MA, Wan Hussin WMA, Majid TA, Sanusi SAM. First and second order Markov chain models for synthetic generation of wind speed time series. Energy 2005;30:693–708. [54] Nfaoui H, Essiarab H, Sayigh AAM. A stochastic Markov chain model for simulating wind speed time series at Tangiers, Morocco. Renew Energy 2004;29:1407–18. [55] Brokish K, Kirgley J. Pitfalls of modeling wind power using Makov chains. In: Proceedings of the power systems conference and exposition 2009. PSCE 09 IEEE/PES; 2009. [56] Sfetsos A. A comparison of various forecasting techniques applied to mean hourly wind speed time series. Renew Energy 2000;21:23–35. [57] Foley AM, Leahy PG, Marvuglia A, McKeogh EJ. Current methods and advances in forescasting of wind power generation. Renew Energy 2012;36:1–8. [58] Bossanyi EA. Short-term wind speed prediction using Kalman filters. Wind Eng 1985;9:1–8. [59] Chen Kuilin, Yu Jie. Short-term wind speed prediction using an unscented Kalman filter based state-space support vector regression approach. Appl Energy 2014;113:690–705. [60] Zuluaga CD, Álvarez MA, Giraldo E. Short-term wind speed prediction based on robust Kalman filtering: an experimental comparison. Appl Energy 2015;156:321–30. [61] Tian Y, Liu Q, Hu Z, Liao Y. Wind speed forecasting based on time seriesadaptive Kalman filtering algorithm. In: Proceedings of the FENDT 2014. Article number 6928287; 2014. p. 315–19. [62] Louka P, Galanis G, Siebert N, Kariniotakis G, Katsafados P, Pytharoulis I, et al. Improvements in wind speed forecasts for wind power prediction purposes using Kalman filtering. J Wind Eng Ind Aerody 2008;96:2348–62. [63] Cassola F, Burlando M. Wind speed and wind energy forecasts through Kalman filtering of numerical weather prediction model output. Appl Energy 2012;99:154–66. [64] Crochet P. Adaptive Kalman filtering of 2-m temperature and 10-m windspeed forecasts in Iceland. Meteorol Appl 2004;11:173–87. [65] Poncela M, Poncela P, Perán JR. Automatic tuning of Kalman filters by maximum likelihood methods for wind energy forecasting. Appl Energy 2013;108:349–62. [66] Costa A, Crespo A, Navarro J, Lizcano G, Madsen H, Feitosa E. A review on the young history of the wind power short-term prediction. Renew Sustain Energy Rev 2008;12:1725–44. [67] Li Gong, Shi Jing. Application of Bayesian model averaging in modeling longterm wind speed distributions. Renew Energy 2010;35:1192–202. [68] Gong Li, Jing Shi, Junyi Zhou. Bayesian adaptive combination of short-term wind speed forecasts from neural network models. Renew Energy 2011;36:352–9. [69] Gong Li Jing Shi. On comparing three artificial neural networks for wind speed forecasting. Appl Energy 2010;87:2313–20. [70] Kaur D, Lie TT, Nair NKC, Vallès B. Wind speed forecasting using hybrid Wavelet Transform-ARMA techniques. AIMS Energy 2014;3(1):13–24. [71] Liu Hui, Tian Hong-Qi, Chen Chao, Li Yan-fei. A hybrid statistical method to predict wind speed and wind power. Renew Energy 2010;35:1857–61.

A. Feijóo, D. Villanueva / Renewable and Sustainable Energy Reviews 56 (2016) 473–483

[72] Aksoy H, Toprak ZF, Aytek A, Ünal NE. Stochastic generation of hourly mean wind speed data. Renew Energy 2004;29:2111–31.

Andrés Feijóo obtained his Ph.D. degree in Electrical Engineering from the Departamento de Enxeñería Eléctrica, Universidade de Vigo, Spain, in 1998. He is now in the same department and is interested in wind energy and its impact in electrical power networks, in particular wind speed simulation and steady-state and dynamic modelling and simulation of electrical machines in wind parks.

483

Daniel Villanueva obtained his Ph.D. degree in Electrical Engineering from the Departamento de Enxeñería Eléctrica, Universidade de Vigo, Spain, in 2012, where he is currently working on the impact of wind energy in power systems, including simulation of wind speeds, analysis of wind power and probabilistic analysis of power systems including wind power generation.