High resolution wind speed forecasting based on wavelet decomposed phase space reconstruction and self-organizing map

Renewable Energy 140 (2019) 17e31

Contents lists available at ScienceDirect

Renewable Energy journal homepage: www.elsevier.com/locate/renene

High resolution wind speed forecasting based on wavelet decomposed phase space reconstruction and self-organizing map € kmen a, Pengfei Li a, Qi Huang b, Zhe Chen a Rui Hu a, Weihao Hu b, *, Nuri Go a b

Department of Energy Technology, Aalborg University, Denmark School of Mechanical and Electrical Engineering, University of Electronic Science and Technology of China, China

a r t i c l e i n f o

a b s t r a c t

Article history: Received 28 May 2018 Received in revised form 16 November 2018 Accepted 9 March 2019 Available online 13 March 2019

Wind power has already brought lots of benefits for energy market and environment. To promote the utilization of wind energy further and achieve the ambitious goals of applying sustainable energy, the improvement of the performance of wind speed forecasting will play a significant role. The whole system will benefit from an accurate wind forecasting system with high resolution. This paper focused on shortterm wind speed forecasting and proposed a method providing the wind speed forecasting with shortterm horizon and high resolution. In this paper, a novel method based on the wavelet decomposition, the phase space reconstruction and the self-organizing map (SOM) artificial neural network was proposed. Study case was provided to demonstrate the performance of the proposed forecasting method. Comparative results from several classic forecasting methods were provided as well. Based on the outcomes of study cases and comparison results, it can be concluded that the proposed method performed well in short term forecasting horizon. © 2019 Elsevier Ltd. All rights reserved.

Keywords: Wind speed forecasting Short-term High-resolution Phase space reconstruction Paralleling computing Self-organizing map

1. Introduction People are seeking for various renewable energy sources to substitute for traditional fossil ones. A new era of sustainable energy is dawning and many researches and practical applications have already been done [1]. Compared with other renewable energy such as solar power, marine power, etc., the application of wind power is being deployed widely in these years [1,2]. Countries like Denmark, Germany, US and so on are leaders of wind power technologies [3e5]. To promote wind power further and increase the wind power penetration, several paths could be utilized. In addition to methods such as improving the Energy Storage Systems (ESS) and wind turbine technologies, enlarging the capacity of ESS at a higher cost, or increasing the reserve of the whole system, enhancing the performance of wind forecasting will also play a significant role as well. The wind forecasting system with higher resolution and larger time horizon is always in need. An accurate wind forecasting could benefit all aspects of the utilization of wind power and overcome many problems caused by high wind power

* Corresponding author. E-mail addresses: [email protected] (R. Hu), [email protected] (W. Hu), ngo@et. €kmen), [email protected] (P. Li), [email protected] (Q. Huang), zch@ aau.dk (N. Go et.aau.dk (Z. Chen). https://doi.org/10.1016/j.renene.2019.03.041 0960-1481/© 2019 Elsevier Ltd. All rights reserved.

penetration [6]. It can provide references for the maintenance of wind turbines, reservation regulation of synchronous machines, and participation of energy market, etc. Therefore, the wind power will be more economic and competitive. Unfortunately, the movement of atmosphere is a high dimension nonlinear system with chaotic properties [7], which makes it very hard or even impossible to predict in some specified scale and accuracy. This opinion has been gradually accepted by the meteorologists since the 1960s when Edward Lorenz developed a simplified mathematical model for atmospheric convection [8]. This kind of system has several adverse characteristics for accurate forecasting. For example, it's highly sensitive to the initial values, and the orbits of the system may fold and unfold in the phase space. These orbits may have very complicated bifurcations with multi-parameters and high dimensions [9e11]. During hundreds of years, people keep trying forecasting weather accurately and long-termly. Not only the knowledge of meteorology itself, but also diverse methods from other subjects are introduced to predict weather. As the prevailing of wind energy, many latest achievements of science and technologies are applied to forecast wind status for wind farms quickly. These forecasting methods can be classified to 4 categories [12] according to different time horizons, namely very short-term wind forecasting, shortterm wind forecasting, middle-term wind forecasting, and long-

18

R. Hu et al. / Renewable Energy 140 (2019) 17e31

term wind forecasting. From technical view [13e15], these methods could be classified as numerical weather prediction (NWP) method, statistical methods, heuristic methods, artificial intelligence methods, and hybrid methods which have their applicable occasions respectively. For example, the simplest persistent method has a good performance for very short-term wind forecast, while it causes quite large errors in middle-term use cases [13]. Artificial Intelligence (AI) based methods such as BP network and SVM performs well in short-term and middle-term wind power prediction. Other classic methods such as ARMA and ARIMA also work well in some specific situations [15]. Normally the ARIMA forecasting model will perform better than the ARMA one due to its integral part in the algorithm, which enables the whole model to deal with the nonstationary of the time series [15,16]. Recently, the hybrid application of AI algorithms and statistical algorithms is a hot topic, in Ref. [17], a hybrid ANN-ARIMA model was proposed and tested. Reference [18] gave a comparison of the performance between ARIMA-ANN and ARIMA-Kalman. This paper mainly focused on improving short-term wind speed prediction whose performance has a great influence on daily operation of the power system integrated with large scale wind farms. A novel method, which is able to provide much higher resolution and keep accuracy and time horizon at the same time, is proposed and verified. As aforementioned, for the studied time horizon, forecasting models could be extremely diverse based on historic data, physical equations, A.I. technologies, or their combinations. Although some of them can perform very well in certain conditions, the resolution of their results is still limited. With power electronic interface-based devices penetrating into the traditional power system continuously, the resolution of wind forecasting results is becoming more and more meaningful. As shown in Fig. 1, low resolution results contain much less information of wind variations and large forecasting errors can happen during the intervals. Meanwhile, power electronic devices response much faster than traditional ones and their capacity is limited by the switches and topologies. Therefore, these errors could force some devices or loads to be tripped off in microgrids or higher reservation at higher costs. It can be directly seen from this viewpoint that a high-resolution forecasting system can help the whole grid operating safely and efficiently. For the methods purely based on statistics and historical data, such as ARMA, ARIMA, etc. the difficulty to increase resolution can be summarized as how to reflect and predict the nonlinearity variation of wind correctly. The heuristic and A.I. based methods are introduced in wind forecasting area later to overcome this

obstacle, but the tuning of these algorithms usually needs lots of attempts and experience, otherwise the overfitting or under fitting problem could limit their application. The obstacle of applying physical based model such as CFD or NWP in short-term wind forecasting is the time validity, economic and time costs brought by oversize computation loads. The huger challenge for the physical way is brought by the atmosphere's chaos essence. As the development of mathematics and physics, in 1960s, chaotic phenomenon was revealed by scientists gradually. It has been proven that even a very simple non-linear equation such as Logistic mapping [19,20] can have chaotic behavior. Lots of typical chaotic features have been found in nature so far. They exist in the atmospheric movement, air turbulence and weather variation. But different from Logistic mapping, they are nonlinear dynamic systems with high dimension and multiple parameters which are tricky and difficult to have closed-form expressions such as Navier-Strokers equations [21]. Once this kind of system enters into the chaotic state, even a small change in initial states may cause a totally different orbit in phase space. Therefore, it's difficult to predict in which orbit will the dynamic system stay in the future and many negative effects will be introduced to the simplified physical forecasting models. In meteorology, accurate long-term weather forecasting even was demonstrated to be impossible [22] even at current technology level. Despite of these frustrating facts, people have found the order existed in the chaos attractors in high dimension could be a useful tool to give a prediction of the nonlinear dynamic system's behavior within certain accuracy tolerance. As a raising subject and hot topic, chaos theory has already had lots of applications. Besides the time horizon and accuracy problems, resolution of the prediction also should be taken into consideration when evaluating the performance of a forecasting system. However, the exploration and research in providing high resolution wind forecasting is not enough. Most of the conventional short-term methods are giving a resolution above 10 min. Some methods could have a very high resolution such as the simplest persistence method, but their time horizons are seriously limited in time domain or geography [6,23]. Therefore, short-term and high-resolution wind speed forecasting is needed, but little information can be found in the literature. In this paper, a forecasting method, which can give the acceptable forecasting results of short-term time horizon and high resolution at the same time, is proposed and verified. The method increases the resolution of prediction results significantly than conventional methods and this contribution could be one of the effective solutions to the aforementioned problems. Additionally,

Fig. 1. The concept and benefit of high-resolution forecasting.

R. Hu et al. / Renewable Energy 140 (2019) 17e31

the proposed method uses only historical time series as inputs, which is very simple to achieve in practical application. As a result, complicated nonlinear mathematic models are avoided. Then the method reconstructs them into large-scale phase space orbits in different time scales to give an extended depiction of wind time series of the observation spot. The reconstruction algorithm is coded based on Compute Unified Device Architecture (CUDA) technology, which promises the time validity of the whole method. At last, SOM based nets were trained to extract the characteristics of wind variation and to predict the future wind speed. The proposed algorithm in this paper reserves the complex non-linearity of the wind and gives a relatively clear physical meaning. Comparison of the outcomes between classic ARIMA methods is also provided. This paper was organizing as following: Chapter 1 was the introduction and background. Chapter 2 described the main theories and necessary principles employed by the proposed method. Chapter 3 gave both the architecture and the design details of the

Fig. 2. The principle of proposed method.

19

proposed model. Chapter 4 was the study case and comparison results which demonstrated the feasibility and performance of the method. Chapter 5 concluded the whole work and contributions. 2. Modelling theories related to proposed forecasting method It could be briefly described that the principle of proposed method is to find approximate historical wind speed data as predicting results. According to the investigation on the adopted data set of 5 years in this paper, it is verified that the approximation of future high-resolution wind speed time series can be found in historical data with the same resolution, as shown in Fig. 2. This discovery in data set is the support of the proposed method. Thus, the proposed method is designed to find the acceptable approximate evolvement of wind speed time series in history as forecasting output. Based on this, the general composition of proposed method is shown in Fig. 3. The architecture is composed of three main parts. Firstly, the historical data is decomposed by wavelets in order to find the most approximate historical periods with different resolutions. Searching the best approximate period in history with different resolutions will be helpful for the accuracy of prediction, for instance, the lower resolution time series can reflect the average wind speed in larger time scale and remove the details of faster variations that caused by turbulences or other factors, which makes the predicting results fit the general trend of wind speed. Higher resolution time series will give more information about the turbulence level and other factors. Thus, it's helpful for the results to cover more details of the variation of wind speed. The second part is phase space reconstruction, which extend 1-dimension time series to a higher dimension time series. The wind speed time series could be treated as the projection of one variable in time axis of a high dimension system, therefore the reconstructed phase space have better potential when forecasting the movement of the system compared with models that only use one-dimension data. At last, the approximate states to forecasting moment in history are searched and SOM nets are employed to clustering these states to obtain a

Fig. 3. The composition of proposed method.

20

R. Hu et al. / Renewable Energy 140 (2019) 17e31

high-resolution short-term forecasting result. Several concepts and theories that related to the proposed method are introduced in detail as following. 2.1. Phase space reconstruction from time series As aforementioned, a single dimension time series of wind speed is insufficient to reflect the variation trend of the chaotic system that it belongs to. In mathematics and physics, the phase space of a dynamic system is a space in which all possible states of a system are represented with a unique point [24]. The characteristic analysis in phase space is significant for the study of dynamic systems. However, in practical situation, the observed data usually is in form of time series. Theoretically, the orbit of a dynamic system in phase space can be analyzed by the partial differential equations (PDEs) using numeric methods like Lorenz system and Logistic system. However, because of the complexity and high dimensions of the real atmosphere system, it is very difficult to get a set of accurate partial differential equations (PDEs) at arbitrary time. The most complicated earth weather methods based on governing PDEs could have hundreds of variables and states [25], which make it a tough task for computers to solve in an available time span, and it will also suffer from the chaos effect [26]. Therefore, transforming the time series into phase space to get the features of strange attractor, which called phase space reconstruction technology, is another important direction of constructing the phase space of dynamic system. The reconstructed phase space will preserve the properties of the dynamical system and can reflect the characteristics of the system's dynamics. The feasibility of time space reconstruction has been proved by Floris Takens, and been concluded as the Takens' theorem [27]. It is significant for promoting the practical application of chaos theory. The Takens' theorem elaborates the relation between embedding dimension obtained from time series and the internal dimension of the dynamic system. For instance, it has been applied in the study of fluids' states by reconstructing the system dynamics from time series only. The conclusions and outcomes have demonstrated the feasibility of the theorem in most situations [28e30]. The crucial task of phase space reconstruction is to determine the time delay and embedding dimension of the time series data appropriately. Diverse methods for determining them have been proposed. The time delay and embedding dimension can be settled separately as two independent parameters or together as coupled parameters. For the previous situation, time delay is usually determined at first using such as autocorrelation method, average displacement method, mutual information method and so on. Then the appropriate embedding dimension will be calculated according to the selected time delay. For the latter one, some researchers thought that the time delay was associated with embedding dimension and proposed many assumptions of the relation between them. Thus, there may be some relation between them and they can be determined at the same time. Some famous methods like C-C method [31,32] belong to the latter strategy. In this paper, the time delay and embedding dimension are determined separately because there are no clear clues that indicate the time delay and embedding dimension of wind speed time series have a deterministic relation. The autocorrelation method can only analyze the linear correlation of 2 time-series. Similarly, the results obtained by average displacement method also cannot ensure that the correlation of data series in each embedded dimension is the lowest [33]. The noise level of the time series may cause a big problem for these two methods as well. In this paper, mutual information (M.I.) based algorithm is chosen to determine appropriate time delay. Different from the above two methods, mutual information method proposed by Fraser and Swinney [34] can determine the nonlinear

correlation between different time series. It employed the information entropy concept from information theory and provide the measurement of statistical independence [35]. For the determination of embedding dimension, several methods could be used such as Lyapunov exponent method [36], False Nearest Neighbors (FNN) method [31], and CAO method [32], etc. FNN is an effective method in this area, but it is sensitive to noises and the selection of the parameters can be subjective sometimes [37]. In this paper, the embedding dimension was settled by using CAO method, which is the improved method of FNN method and overcomes the flaws of FNN. The structure of reconstructed phase space from time series is shown in Fig. 4. The function of aforementioned two parameters can be clearly seen from this figure. Time delay corresponds to M and embedding dimension corresponds to N. 2.2. Extracting characteristics of reconstructed phase space of wind speed Once the phase space of the dynamic system is successfully reconstructed, the trajectories in it can be utilized to forecast the next movement of the system. The forecasting methods in phase space can also be classified into several categories according to different viewpoints. From the number of steps which will be given by the forecasting algorithms, they can be divided into single step forecasting and multi-step forecasting. Regard to single step forecasting methods, as the point on the trajectories in the phase space can determine the system's state uniquely and it is usually considered that the dynamic system will move continuously along the orbit, several points in the phase space that have minimum Euclidean distances with the current point in phase space could be used to predict the next movement after being handled by some averaging method coarsely. Furthermore, to improve these basic methods which could be easily affected by the quality and accuracy of reconstructed phase space [37], other classical methods such as ARMA, ARIMA or A.I. based methods in time series are also applicable in phase space [38]. However, to give a multi-step high-resolution prediction of a high dimension nonlinear system like the air movement in the atmosphere, situations in reconstructed phase space still remain complicated if one wants to find accurate mathematic models. Regard to the multi-step forecasting, the forecasting methods can be classified into direct multi-step forecasting and indirect multi-step forecasting basically. Direct multi-

Fig. 4. The structure of reconstructed phase space.

R. Hu et al. / Renewable Energy 140 (2019) 17e31

step forecasting directly gives all the predicted value of demanded steps in the future at one time, thus the model of it will probably be very complicated and should be in high accuracy. In contrast, indirect multi-step forecasting gives the multi-step values through making prediction for each step in the future using the single step forecasting model [18]. In the application of wind speed forecasting which has aforementioned chaotic features, the errors of indirect forecasting could become unacceptable with increasing of steps. The forecasting methods could also be categorized according to which aforementioned mathematic model that they are using, such as ANN based methods, ARIMA based methods, etc. Considering the application in wind energy, once the forecasting results are given, it will directly affect a series of later operations of wind farms, TSOs, and energy markets, etc. The high-performance multi-step forecasting system will be more attractive and meaningful. 2.3. High resolution wind speed forecasting As narrated above, the resolution is a very important index for a high-performance wind forecasting system. As the prevailing of power electronic devices and facilities in power system nowadays, researchers have come up with innovative interfaces and applications constantly [39e42]. They do provide lots of flexibility for energy integration and energy conversion, but some of their characteristics also bring challenges and troubles for traditional power system such as fast response and low capacity. Consequently, from the view of forecasting, the errors and resolution of the energy forecasting system will be more influential. In the power system with high wind power penetration, the unbalancing caused by forecasting errors have to be successfully compensated by traditional synchronous machines and energy storage systems etc. as well, otherwise the errors can probably cause other severe issues [43] such as tripping off converters or cutting off loads when the system trying to keep stable especially in high penetration system or microgrids. Thus, the high-resolution forecasting could be helpful for the system to scheduling, dispatching, and adjust reservation in time as the errors caused by turbulences or fast wind variation during the forecasting interval are reduced. The more the predicting steps can be provided, the better the whole system will benefit. To choose an expedient step size, the following aspects are considered. Firstly, the output of the wind turbines currently has been required to be able to vary in very fast rate by many grid codes [3,44] in order to deal with contingencies. The modern wind turbines based on pitch control, active stall control combined with PMSG or DFIG can easily output from 0 to rated power in tens of seconds. Forecasting steps around this scale won't lose too much details and can help wind farms to arrange their generation plan. Secondly, as the wind turbine themselves have rotational components with certain inertial, wind speed variations within roughly dozens of seconds could be filtered in some extend. As a result, too small step size will be less meaningful for wind forecasting applied in wind energy, and it will increase the computing loads on the contrary. In this paper, the appropriate resolution will be chosen according to the results of wavelet analysis. For the study case in this paper, the determined resolution is around 1 min, namely the proposed method will give a 2-h short-term time horizon prediction with 120 wind speed forecasting points. 2.4. Chaos identification in time series Generally, the behavior of nonlinear system could be considered as the result of the competition between the driving force and dissipating force. Some relations of these two forces may cause a

21

chaotic behavior, but other relations between them could break the chaotic state and bring the dynamic system to other diverse orbits. For instance, when the driving force of the logistic system is dominant force, the system will stay at a well-determined state. When the dissipating force increase and is able to compete with the driving force, the system's behavior will become quite complicated and finally fall into chaos. For high dimension non-linear systems such as atmosphere system, things will get much more complicated. Hence, it's very difficult to forecast its movement in chaos state using only one-dimension time series and this is also the reason that phase space reconstruction is needed. In this paper, the existence of chaos attractor in adopted time series of wind speed is proven by using Lyapunov exponent method. The flow of wind can be considered as a complex fluid dynamic system with relatively high Reynolds number. This kind of system may have another important feature that should be considered when conduct forecasting, namely multi-scale characteristic. Turbulence is a typical phenomenon of chaos and is very common in nature, but there hasn't been a universal theoretical model for it [45]. It still remains a huge challenge for CFD to simulate the turbulence accurately especially in a large and complex space. Some researches [46e48] based on fractal theory have indicated that the turbulences of air flow can be treated as composed by multi-scale eddies. These eddies will interact with each other in multi-scales and finally form the shapes of turbulences. Based on this idea, the time series used in this paper was decomposed into multi-scale by wavelets, and then the phase space of each scale was reconstructed. The multi-scale phase space reflected the variation of wind speed in different time scale. After the decomposition, Lyapunov exponent method was employed to identify the chaos behavior, and the scale which contains chaotic behavior will be proceeded by the SOM network. 3. Design and implement of proposed method using GPU computing In this part, the design and implement details of the proposed method is narrated. As discussed in Section 2, the input of the proposed method in training stage will be the wind speed time series. The time series will be decomposed into multi-scale time series by ‘Coiflet’ wavelets at first, then the data handling in the following procedures of each level is identical. After the decomposition, phase space reconstruction of each scale will be conducted respectively. To achieve a large-scale high-resolution reconstructed phase space, GPU computing is employed for accelerating. With the help of CUDA libraries all the single thread programs are redesigned and recoded into paralleled ones which run on GPU. As elaborated above, the existence of chaotic behavior in each scale will be examined by Lyapunov exponent method. After above steps, the algorithm is ready for forecasting. The input for forecasting stage is the decomposed time series before the forecasting time point. These segments of input time series data are generated from the same wavelet decomposed method in the training process. Then, many historical segments of time series that approximate to current states in each scale will be found from the large-scale phase space that reconstructed in training stage. The SOM nets will cluster the segments that have same variation characteristics and the clustered historical segment related to winner neuron can reflect the evolvement of current states. The implement of the architecture of proposed method was shown in the following Fig. 5. 3.1. Multi-resolution decomposition based on wavelets Nowadays, wavelets are powerful tools which are widely used in

22

R. Hu et al. / Renewable Energy 140 (2019) 17e31

Fig. 5. The architecture of the proposed method.

signal analysis, noise cancelling and data compression, etc. [49]. In this paper, wavelet is used for decomposing time series into multiscales. Compared to traditional methods such as Fourier transformation, or short-time Fourier transformation, wavelet transformation has a lot of advantages [49,50] which could benefit the application in wind speed time series analysis. For instance, the wavelet functions have compact support which could give the information of the signal both in time domain and frequency domain simultaneously by using shifting and scaling. Thus, the results not only contain the frequency components of wind speed, but also the variation of each frequency component with time. Different from Fourier transformation that have a single fixed resolution, some wavelets are able to conduct multi-resolution analysis (MRA), which means if they are used in wind speed analysis, both the low frequency changes in average speed and high frequency changes could be observed in the results. In this paper, wavelets are used as filters to get multi-resolution time series which can reflect wind speed variations in different frequencies. These decomposed time series are the inputs of paralleled phase space reconstruction program. There are lots of wavelet mother functions and scale functions that can be used for wind speed data decomposition. In this paper, the “Coiflets” [51] which designed by Ingrid Daubechies were chosen to be the mother wavelet function. This wavelet set has a lot of advantages when applied in wind speed analysis. Compared to the simplest Daubechies wavelet, namely Haar wavelet, it can have much higher vanishing moment that makes it able to keep more details of smooth signals. Compared to the other Daubechies wavelets, it is a near symmetric function which has wavelet and scale function vanishing moment at the same time. The near symmetric feature makes it able to have a linear phase characteristic. Meanwhile, compared to the wavelets represented by Mexhat or Gaussian wavelet, it has a compact support that will reduce the computation loads significantly.

3.2. Phase space reconstruction using mutual information based KDE and CAO method After decomposition, the data series in each scale will be reconstructed into phase space. There are two key parameters need to be settled in this stage, namely time delay and embedding dimension. In this paper, the appropriate time delay was selected by analyzing mutual information calculated by kernel density estimation (KDE) [52] based probability density computing method. According to the Takens' theorem, when the embedding dimension satisfies the following condition:

edim > 2d þ 1

(1)

where d is the dimension of set A and edim is the embedding dimension. The phase space orbits composed by delayed time series will be the diffeomorphism of the strange attractor which could provide orderly features for forecasting system's behavior. To settle the value of embedding dimension, the mutual information of the delayed time series and the original time series can be calculated by the following equations [34]:

HðSÞ ¼ 

n X

Ps ðsi Þlog2 Ps ðsi Þ

(2)

i¼1

Hd ðSÞ ¼ 

n X

    Pq sj log2 Pq sj

(3)

j¼1

Where H(S) and Hd(S) is the entropy of the system composed by the sampled data set S and delayed one, respectively. Ps(si) and Pq(sj) are the probability of event si and sj, respectively. Then the mutual

R. Hu et al. / Renewable Energy 140 (2019) 17e31

information I(Q,S) can be defined by the following equations:

IðQ ; SÞ ¼ HðQ Þ  HðQ jSÞ

HðQ jsi Þ ¼ 

n h  . h  . i i X Ps si ; qj Ps ðsi Þ log2 Ps si ; qj Ps ðsi Þ

(4)

(5)

EðmÞ ¼

Nm Xt 1 aði; mÞ N  mt i¼1

E1 ðmÞ ¼ Eðm þ 1Þ=EðmÞ

i¼1

E* ðmÞ ¼

Substitute equation (5) into equation (4), then we can get:

  3 2  Psq si ; qj XX   5 IðQ ; SÞ ¼ Psq si ; qj log2 4 Ps ðsi ÞPq qj i j

(6)

The probabilities used by mutual information computation can be calculated by aforementioned KDE method proposed in Ref. [52]. The delayed time corresponding to the first minimum point of the calculated mutual information then could be chosen as the time delay. The above algorithm is realized on GPU and its diagram is shown in Fig. 6. The determination of embedding dimension comes after the time delay determination. In this paper, CAO method was employed to give an index to choose embedding dimension. The principle of CAO method is quite simple, it is based on identifying the variation of the distances among nearest neighbors in phase space. The method can be depicted by the following equations [32]:

  xdþ1 ðiÞ  xdþ1 NN ðiÞ  aði; dÞ ¼    xd ðiÞ  xd NN ðiÞ

 Nm Xt 1  xði þ mtÞ  xNN ði þ mtÞ N  mt i¼1

 E2 ðmÞ ¼ E* ðm þ 1Þ E* ðmÞ

23

(8)

(9)

(10)

(11)

Where a(i,d) is the parameter to determine the false nearest neighbor (FNN). xd(i) and xNN d (i) is the vector and its nearest neighbor at d dimension. N is the number of data in time series and m is the embedding dimension. E1 and E2 are the indices of choosing appropriate embedding dimension. If the time series is deterministic, then the embedding dimension will exist and E1 will be stable when the embedding dimension m exceeds a certain value. E2 is another index to identify whether the E1 is stable because in practical cases, E1 could be fluctuated at last due to noises and disturbances. As a consequence, it is hard to tell whether E1 is stable directly. The introduced E2 will change with the embedding dimension, then after E1 is stable, E2 will be equal to 1. So, the appropriate embedding dimension will be the value when E1 become stable and E2 approximates 1. This algorithm is also accelerated by GPU, and its diagram is shown in Fig. 7.

(7) 3.3. SOM based forecasting in phase space A self-organizing map (SOM) or self-organizing feature map (SOFM) is a type of competitive artificial neural networks which is first introduced by the Finish professor Teuvo Kohonen in the 1980s [53]. This kind of network is capable to be trained with unsupervised learning to produce a low-dimensional, normally 2dimension discretized representation of the input space of the training samples. This algorithm based on competitive learning is fast and efficient. But different from the “winner-takes-all” principle of the basic competitive network, the weights of nerves in winner nerve's neighborhood will also be adjusted. The winner's neighborhood is determined by the distance function, which could be Mexican hat function, Gaussian function and so on. After adjusting the weights of neural in output layer continuously during each training iteration, the position of output neural in one or twodimension plane will be in order. The update formula for a neuron v with weight vector Wv(s) is

Wv ðs þ 1Þ ¼ Wv ðsÞ þ qðu; v; sÞ  aðsÞ  ðDðtÞ  Wv ðsÞÞ

Fig. 6. The diagram of Paralleled Mutual Information Method.

(12)

Where s is the step index, t is the index of the training sample, u is the index of the best matching unit for D(t), a(s) is a monotonically decreasing learning coefficient and D(t) is the input vector; q(u,v,s) is the neighborhood function which gives the distance between the neuron u and the neuron v in step s. The coordinates of these output nerves reflect the distribution of different input patterns, in other words, the coordinates represent the topographic features of the input samples. Therefore, the statistic features of input patterns could be found. Given these advantages, SOM based networks are widely used in data clustering applications. In this paper, SOM net is employed to cluster the trajectory segments from reconstructed phase space. For the phase space in each decomposed scale, there will be a SOM net responsible for clustering the approximate trajectories that found in reconstructed phase space. These nets use the data which have

24

R. Hu et al. / Renewable Energy 140 (2019) 17e31

Fig. 7. The diagram of Paralleled CAO Method.

been processed by the methods introduced in this chapter to train themselves. Each sample was an orbit segment in phase space. After the net is trained, each cluster in the output layer will represent for the similar system movement occurred in history. There will be enough nerves in output layer to make the distance variation in each cluster stay in acceptable range.

reachability  distε;MinPts ðo; pÞ ¼ ( UNDEFINED   Max core  distε;MinPts ðpÞ; distðp; oÞ

if jNε ðpÞj < MinPts Otherwise (14) 0

3.4. OPTICS based forecasting in phase space In this part, another clustering algorithm, namely OPTICS (Ordering Points to Identify the Clustering Structure) [54] is employed for providing comparative outcome. The clustering algorithm can be described as:

core  distε;MinPts ¼ UNDEFINED MinPts  th smallest distance to Nε ðpÞ

if jNε ðpÞj < MinPts Otherwise (13)

A point p is a core point if at least MinPts points are found within its ε-neighborhood Nε(p). ε describes the maximum distance (radius) to consider, and MinPts, describing the number of points required to form a cluster. The reachability-distance of another point o from a point p is either the distance between o and p, or the core distance of p, whichever is bigger:

If point p and point o are nearest neighbors, this is the ε <ε that need to assume in order to have p and o belong to the same cluster. The whole structure of the proposed framework is shown in following Fig. 8. The inputs and outputs of OPTICS in each layer are identical with previous method based on SOM. 4. Study cases In this chapter, the data collected from National Renewable Energy Laboratory [55] was used to demonstrate the proposed method. The anemometers and sensors are located in National Wind Technology Center (NWTC) whose latitude and longitude are 39.91 N, 105.235 W, respectively. The wind speed time series sampled every minute from year 2011 to year 2016 at 80 m height were employed to be representative for the hub height of mainstream wind turbines. 4.1. Study case of proposed method The method proposed in this paper was demonstrated by the

R. Hu et al. / Renewable Energy 140 (2019) 17e31

25

Fig. 8. The proposed structure based on OPTICS.

sampled wind speed data with resolution of 1 min. The wind speed time series from year 2013 to year 2015 was decomposed and reconstructed into phase space orbits, the segments of the orbits were employed to train the SOM nets. Then the wind speed series of year 2016 was chosen to be the demonstrating set to estimate the performance of the proposed method. As narrated above, the first step is wavelet decomposing. In this paper, all the data was decomposed into 8 levels at first. Different decomposed level contains different information of wind speed

Table 1 Choosing of time delay for phase space reconstruction in 2015. Levels

Delayed times

First Minimum Mutual Information

A8 A7 A6 A5 A4 A3 A2 A1

4 9 14 19 23 40 47 53

0.0176 0.0321 0.0233 0.0338 0.0683 0.0063 0.0488 0.1082

Fig. 9. The mutual information results of different decomposed levels of 2015 wind speed data.

Fig. 10. E1 and E2 for choosing the embedding dimensions in 2015. Table 2 Reconstructed phase space information of each decomposed level from 2011 to 2015. Years Level A5

2011 2012 2013 2014 2015

A4

A3

A2

A1

Time Delay

Embedding Dimension

Time Delay

Embedding Dimension

Time Delay

Embedding Dimension

Time Delay

Embedding Dimension

Time Delay

Embedding Dimension

35 35 32 28 19

27 20 27 29 27

61 25 41 31 23

20 30 28 27 30

58 42 54 30 40

24 27 25 24 25

38 79 66 67 47

25 29 19 35 22

44 66 56 62 53

21 35 28 26 25

R. Hu et al. / Renewable Energy 140 (2019) 17e31

27

Then the mutual information of each scale was calculated, and the time delay corresponding to first minimum value of mutual information was selected as the optimal time delay for phase space reconstruction. The variation of the mutual information and the optimal time delay were shown in Fig. 9 and Table 1. After the time delay was determined, the appropriate embedding dimension could be settled by the CAO method. The calculated E1 and E2 of each scale are shown in Fig. 10. According to the principle descripted in Chapter 2, the embedding dimension where the E1 becomes stable and E2 approximate to 1 can be selected as the optimal embedding dimension. The final results of choosing the embedding dimensions were shown in Table 2. So far, the phase space of the wind can be reconstructed by these two parameters in each scale. For wavelet levels with D6, D7 and D8, the parameter E2 has already trend to random variation, the predictability of these two scales are limited and the reconstruction will be imprecise. Therefore, the forecasting results will be composed using D1-D5 and A1. After the phase spaces of each scale were reconstructed, the existence of chaos attractors could be examined by apply maximum Lyapunov Exponent methods. These outcomes are shown in Fig. 11. All the maximum Lyapunov exponents of each scale from year 2011e2015 were positive, which proved the existence of chaos attractor and using only 1-dimension time series to conduct highresolution forecasting may not obtain an ideal result. The phase space reconstruction of massive data will even be impractical if no paralleled computation technology is applied. In this paper, the huge data set caused this problem inevitably. To overcome this obstacle, NVIDIA CUDA technology was applied in the implementation of proposed method. It effectively reduces the time cost of the whole method and the reconstruction task can be finished in one or two weeks. 4.2. Comparisons with conventional methods In this paper, the comparison with other popular conventional forecasting methods was also conducted. As the simplest method, persistence model is not suitable in a multi-step and massive-step occasion and it will be excluded in the comparison. The ARMA model and ARIMA model was compared because they are already very matured in time series forecasting. The seasonal and nonseasonal “AR”, Difference lags and “MA” were obtained by automatic evaluation tools in R language, which has a standard time series analyzing library. These parameters were determined by the classic autocorrelation function (ACF) and partial autocorrelation function (PACF) methods. The model information is shown in the following Table 3. Besides, the forecast method based on ELM (Extreme Learning Machine) [56] is also employed for comparison. To achieve multi steps ahead forecast results, the ELM is invoked iteratively. To demonstrate the proposed method, the forecasting results produced by the proposed method, iterative ELM and conventional ARIMA method are shown in Fig. 12. Their RMSE are plotted in Fig. 13. It is calculated by the following equation: Fig. 11. The maximum lyapunov exponents of the reconstructed phase space in each decomposed level from 2011 to 2015. Table 3 Parameters of conventional seasonal ARIMA method.

variation. The decomposition can also reduce the requirement of computation resources. It may cause insufficient memory or unacceptable calculation time when treating the original sample points as a whole data set without decomposition.

Model

Parameters AR

D

MA

SAR

SMA

Seasonality

Seasonal ARIAM model

4

1

3

8

41

1440

28

R. Hu et al. / Renewable Energy 140 (2019) 17e31

Fig. 12. The prediction results and comparison of proposed method.

R. Hu et al. / Renewable Energy 140 (2019) 17e31

Fig. 12. (continued).

29

30

R. Hu et al. / Renewable Energy 140 (2019) 17e31

Fig. 13. The RMSE of Proposed Method and Comparative methods.

RMSE ¼

vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi u T uP 2 u b tt¼1 ð w t  wt Þ T

based paralleled computation. The outcomes of comparison also demonstrated the positive effects of the wavelet decomposition and phase space reconstruction in wind speed forecasting.

(15)

b t is the predicted wind speed value, wt is the actual wind Where w speed value and T is the predicted period at a time. It can be seen that the proposed method can perform well in high resolution situation. And the two clustering methods performs identically in proposed framework. The RMSE of them is smaller than ARIMA and ELM based method. The iterative ELM can give good results in limited steps ahead situation. However, the conventional structure based on ELM may suffer from error accumulation or other aforementioned features when forecast steps increase largely. The conventional ARIMA method could give acceptable results in much less steps, thus it may not be suitable for high-resolution wind speed forecasting application.

5. Conclusions The performance of wind speed forecasting methods plays a significant role for the whole grid. It is the input of wind power forecasting system and has a direct relation with the efficiency of wind farm planning and daily operation. Many issues that caused by high wind power penetration could also be solved with the help of an accurate forecasting system. As a consequence, the pursuit of constantly improving the accuracy, time horizon and resolution of the wind forecasting system will never be stopped. This paper gave a brief review of the conventional wind forecasting methods and proposed an innovative hybrid method aiming at increasing both the time horizon and forecasting resolution. The goals were achieved by rationally using basic physical and mathematical principles, wavelet decomposition, phase space reconstruction and clustering methods synthetically. By the proposed configuration, the chaotic characteristics of wind were identified from phase space and utilized to forecast the wind time series. The study case of the proposed method was provided and gave a more satisfying result compared to the other methods based on ARIMA, and ELM. The obstacles brought by heavy computation loads introduced by phase spaced reconstruction were successfully overcome by CUDA

Acknowledgment This work was supported by the National Natural Science Foundation of China (51707029). References [1] C.C. Mitigation, IPCC Special Report on Renewable Energy Sources and Climate Change mitigation, 2011. [2] S. Chu, A. Majumdar, Opportunities and challenges for a sustainable energy future, Nature 488 (7411) (2012) 294e303. [3] V. Akhmatov, P.B. Eriksen, A large wind power system in almost island operationda Danish case study, IEEE Trans. Power Syst. 22 (3) (2007) 937e943. [4] I. Erlich, W. Winter, A. Dittrich, Advanced grid requirements for the integration of wind turbines into the German transmission system, in: Power Engineering Society General Meeting, 2006. IEEE, IEEE, 2006, p. 7. [5] P. Meibom, K.B. Hilger, H. Madsen, et al., Energy comes together in Denmark: the key to a future fossil-free Danish power system, IEEE Power Energy Mag. 11 (5) (2013) 46e55. [6] S.S. Soman, H. Zareipour, O. Malik, et al., A review of wind power and wind speed forecasting methods with different time horizons, in: North American Power Symposium (NAPS), 2010, IEEE, 2010, pp. 1e8. [7] J. Shukla, Predictability in the midst of chaos: a scientific basis for climate forecasting, Science 282 (5389) (1998) 728e731. [8] E.N. Lorenz, Deterministic nonperiodic flow, J. Atmos. Sci. 20 (2) (1963) 130e141. [9] H. Tong, K.S. Chan, D.R. Cox, et al., A personal overview of non-linear time series analysis from a chaos perspective [with discussion and rejoinder], Scand. J. Stat. (1995) 399e445. [10] R.C. Hilborn, J.C. Sprott, Chaos and nonlinear dynamics: an introduction for scientists and engineers, Am. J. Phys. 62 (9) (1994) 861e862. [11] R.C. Hilborn, J.C. Sprott, Chaos and nonlinear dynamics: an introduction for scientists and engineers, Am. J. Phys. 62 (9) (1994) 861e862. [12] L. Jones, R. Zavadil, W. Grant, The future of wind forecasting and utility operations, IEEE Power Energy Mag. 3 (6) (2005) 57e64. [13] M. Lydia, S.S. Kumar, A comprehensive overview on wind power forecasting, in: IPEC, 2010 Conference Proceedings, IEEE, 2010, pp. 268e273. [14] X. Wang, P. Guo, X. Huang, A review of wind power forecasting models, Energy Procedia 12 (2011) 770e778. [15] Y.K. Wu, J.S. Hong, A literature review of wind forecasting technology in the world, in: Power Tech, 2007 IEEE Lausanne, IEEE, 2007, pp. 504e509. [16] M. Lei, L. Shiyan, J. Chuanwen, et al., A review on the forecasting of wind speed and generated power, Renew. Sustain. Energy Rev. 13 (4) (2009) 915e920. [17] E. Cadenas, W. Rivera, Wind speed forecasting in three different regions of Mexico, using a hybrid ARIMAeANN model, Renew. Energy 35 (12) (2010)

R. Hu et al. / Renewable Energy 140 (2019) 17e31 2732e2738. [18] H. Liu, H. Tian, Y. Li, Comparison of two new ARIMA-ANN and ARIMA-Kalman hybrid methods for wind speed prediction, Appl. Energy 98 (2012) 415e424. [19] R.M. May, Simple mathematical models with very complicated dynamics, Nature 261 (5560) (1976) 459e467. [20] J.D. Farmer, J.J. Sidorowich, Predicting chaotic time series, Phys. Rev. Lett. 59 (8) (1987) 845. [21] H.K. Versteeg, W. Malalasekera, An Introduction to Computational Fluid Dynamics: the Finite Volume method, Pearson Education, 2007. [22] E. Lorenz, Some aspects of atmospheric predictability, in: Problems And Prospects in Long and Medium Range Weather Forecasting, Springer Berlin Heidelberg, 1984, pp. 1e20. rez, L.D. Aponte-Bermúdez, J.M. Morell, et al., CariCOOS: Real-Time [23] J.M. Pe Data Validation of High-resolution wind forecast, OCEANS'15 MTS/IEEE Washington. IEEE, 2015, pp. 1e7. [24] D.D. Nolte, The tangled tale of phase space, Phys. Today 63 (4) (2010) 33e38. [25] J. Smagorinsky, General circulation experiments with the primitive equations: I. The basic experiment, Mon. Weather Rev. 91 (3) (1963) 99e164. [26] E.N. Lorenz, A study of the predictability of a 28-variable atmospheric model, Tellus 17 (3) (1965) 321e333. [27] J.P. Huke, Embedding Nonlinear Dynamical Systems: A Guide to Takens' theorem, 2006. [28] Y. Wang, J. Wang, X. Wei, A hybrid wind speed forecasting model based on phase space reconstruction theory and Markov model: a case study of wind farms in northwest China, Energy 91 (2015) 556e572. [29] D. Lei, W. Lijie, H. Shi, et al., Prediction of wind power generation based on chaotic phase space reconstruction models, in: Power Electronics and Drive Systems, 2007. PEDS'07. 7th International Conference on, IEEE, 2007, pp. 744e748. [30] B. Sivakumar, A.W. Jayawardena, T. Fernando, River flow forecasting: use of phase-space reconstruction and artificial neural networks approaches, J. Hydrol. 265 (1) (2002) 225e245. [31] M.B. Kennel, R. Brown, H.D.I. Abarbanel, Determining embedding dimension for phase-space reconstruction using a geometrical construction, Phys. Rev. A 45 (6) (1992) 3403. [32] L. Cao, Practical method for determining the minimum embedding dimension of a scalar time series, Phys. D Nonlinear Phenom. 110 (1e2) (1997) 43e50. [33] H. Ma, C. Han, Selection of embedding dimension and delay time in phase space reconstruction, Front. Electr. Electron. Eng. China 1 (1) (2006) 111e114. [34] A.M. Fraser, H.L. Swinney, Independent coordinates for strange attractors from mutual information, Phys. Rev. A 33 (2) (1986) 1134. [35] S.S. Haykin, S.S. Haykin, S.S. Haykin, et al., Neural Networks and Learning machines[M], Pearson, Upper Saddle River, NJ, USA, 2009. [36] H.D.I. Abarbanel, R. Brown, M.B. Kennel, Local Lyapunov exponents computed from observed data, J. Nonlinear Sci. 2 (3) (1992) 343e365.

31

[37] C.J. Cellucci, A.M. Albano, P.E. Rapp, Comparative study of embedding methods, Phys. Rev. E 67 (6) (2003), 066210. [38] Z. Xue-Qing, L. Jun, Chaotic Time Series Prediction Model of Wind Power Based on Ensemble Empirical Mode Decomposition-Approximate Entropy and reservoir, 2013. [39] D. Liu, F. Deng, Z. Chen, Five-level active-neutral-point-clamped DC/DC converter for medium voltage DC grids, IEEE Trans. Power Electron. (2016). [40] W. Wang, L. Yan, X. Zeng, et al., Principle and design of a single-phase inverter based grounding system for neutral-to-ground voltage compensation in distribution networks, IEEE Trans. Ind. Electron. (2016). [41] M. Barnes, J. Kondoh, H. Asano, et al., Real-world microgrids-an overview, in: System of systems engineering, 2007. SoSE'07. IEEE International Conference on, IEEE, 2007, pp. 1e8. [42] N. Wang, J. Li, W. Hu, B. Zhang, Q. Huang, Z. Chen, Optimal reactive power dispatch of a full-scale converter based wind farm considering loss minimization, Renew. Energy 139 (2019) 292e301. [43] R. Piwko, D. Osborn, R. Gramlich, et al., Wind energy delivery issues [transmission planning and competitive electricity market operation], IEEE Power Energy Mag. 3 (6) (2005) 47e56. [44] I. Erlich, W. Winter, A. Dittrich, Advanced grid requirements for the integration of wind turbines into the German transmission system, in: Power Engineering Society General Meeting, 2006. IEEE, IEEE, 2006, p. 7. [45] E.C. Waymire, Probability & incompressible Navier-Stokes equations: an overview of some recent developments, Probab. Surv. 2 (2005) 1e32. -Agel, Synthetic turbulence, fractal [46] S. Basu, E. Foufoula-Georgiou, F. Porte interpolation, and large-eddy simulation, Phys. Rev. E 70 (2) (2004), 026310. [47] New Approaches in Modeling Multiphase Flows and Dispersion in Turbulence, Fractal Methods and Synthetic Turbulence, Springer Science & Business Media, 2011. [48] K.R. Sreenivasan, C. Meneveau, The fractal facets of turbulence, J. Fluid Mech. 173 (1986) 357e386. [49] B. Jawerth, W. Sweldens, An overview of wavelet based multiresolution analyses, SIAM Rev. 36 (3) (1994) 377e412. [50] R.M.C. Fern andes, H.N.D. Rojas, An overview of wavelet transform applications in power systems, in: 14 th PSCC, vol. 1, 2002, pp. 1e8. Servilla, Session. [51] I. Daubechies, B.J. Bates, Ten lectures on wavelets, J. Acoust. Soc. Am. 93 (3) (1993), 1671-1671. [52] Y.I. Moon, B. Rajagopalan, U. Lall, Estimation of mutual information using kernel density estimators, Phys. Rev. E 52 (3) (1995) 2318. [53] T. Kohonen, The self-organizing map, Proc. IEEE 78 (9) (1990) 1464e1480. [54] Hans-Peter Kriegel, et al., Density-based clustering, Wiley Interdiscip. Rev.: Data Min. Knowl. Discov. 1 (2011) 231e240, 3. [55] http://www.nrel.gov/midc/. [56] Guang-Bin Huang, Qin-Yu Zhu, Chee-Kheong Siew, Extreme learning machine: theory and applications, Neurocomputing 70 (1e3) (2006) 489e501.