- Email: [email protected]
Wind speed forecasting approach based on Singular Spectrum Analysis and Adaptive Neuro Fuzzy Inference System
Accepted Manuscript Wind speed forecasting approach based on Singular Spectrum Analysis and Adaptive Neuro Fuzzy Inference System Sinvaldo Rodrigues Moreno, Leandro dos Santos Coelho PII:
S0960-1481(17)31190-4
DOI:
10.1016/j.renene.2017.11.089
Reference:
RENE 9491
To appear in:
Renewable Energy
Received Date: 6 August 2016 Revised Date:
13 October 2017
Accepted Date: 29 November 2017
Please cite this article as: Moreno SR, dos Santos Coelho L, Wind speed forecasting approach based on Singular Spectrum Analysis and Adaptive Neuro Fuzzy Inference System, Renewable Energy (2018), doi: 10.1016/j.renene.2017.11.089. This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.
ACCEPTED MANUSCRIPT
RI PT
Sinvaldo Rodrigues Moreno1 and Leandro dos Santos Coelho1,2*
1 Department of Electrical Engineering, Federal University of Parana Zip code 81531-980, Curitiba, PR, Brazil, 2 Industrial and Systems Engineering Graduate Program (PPGEPS)
SC
Pontifical Catholic University of Parana (PUCPR) Imaculada Conceicao, 1155, Zip code 80215-901, Curitiba, PR, Brazil * Corresponding author Tel/fax numbers: (55) 41 32711345 Email: [email protected]
M AN U
3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31
Abstract: As a promising renewable energy source, wind power has environmental benefits, as well as economic and social ones. Due these characteristics, wind farm has grown fast in the last five years, and in some countries, it has already surpassed conventional sources, such as hydro and coal plants. However, owing to the uncertainty of wind speed, it is essential to build an accurate forecasting model for large-scale wind power penetration. This study proposes a hybrid approach that combines the Singular Spectrum Analysis (SSA), which rarely presents application in literature on wind speed forecasting, and a Computing Natural paradigm called Adaptive Neuro Fuzzy Inference System (ANFIS). The SSA decomposes the original wind speed into various components, so these components are pre-processed regarding to the level of original wind series information remained. The main components selected to reconstruct the original series have in their structure the information about trend and harmonic components. The final remaining components grouped are labeled as noise. The ANFIS model uses these two information to construct the model applied to forecasting the next wind speed value. In this way, the hybrid model can learn the trend and the harmonic structure of the wind time series. Experimental results show that prediction errors are significantly reduced using the proposed technique to perform 10min one-step-ahead and -step-ahead wind speed forecast.
TE D
2
Wind Speed Forecasting Approach Based on Singular Spectrum Analysis and Adaptive Neuro Fuzzy Inference System
EP
1
Keywords: wind speed, hybrid forecasting approach, wind farm, machine learning, singular spectrum
33
analysis, neuro fuzzy inference system, time series forecasting.
34
1 Introduction
35
The renewable energy sources, particularly wind power is becoming more popular because of
36
ecological concerns and continuously rising prices of rapidly depleting conventional energy
37
sources. However, several drawbacks complicate the large-scale penetration into large fully
38
connected systems, as the intermittent nature of wind power. Mostly due to wind speed related
39
uncertainties affecting the reliability of energy supply, reliable forecasts methods can effectively
AC C
32
1
ACCEPTED MANUSCRIPT
reduce the risk to the power system. Machine Learning methods have attracted more attention than
41
traditional methods (such as numerical weather predictions and statistical approaches) when the
42
forecasting horizon are comprised between 1 to 3 hours or steps ahead [1-3].
RI PT
40
43
Many approaches applying Natural Computing and Machine Learning concepts have been proposed
45
recently to predicted wind speed. The approaches include promising mathematical techniques such
46
as Support Vector Machine (SVM), Wavelet Transform Method, Artificial Neural Networks,
47
Adaptive Neuro Fuzzy Inference System (ANFIS), Extreme Learning Machine (ELM), among
48
others. Despite the quality of the results, each technique has its own advantages and limitations.
49
Recent work conducted by [4] presents a combination of ARIMA (Autoregressive Integrated
50
Moving Average) model, ELM, SVM and LSSVM (Least Square SVM) to build a model of inputs,
51
and applied it a GPR (Gaussian Process Regression) to perform short-term wind speed forecasting.
52
In [5] a novel bidirectional mechanism and a backward forecasting model based on ELM for wind
53
power time series forecasting have been introduced. This model, as reported by the authors, can
54
achieve good results when applied to forecast wind farm power production for six time horizons as
55
1 up to 6 h.
M AN U
TE D
56
SC
44
A hybrid intelligent forecasting model based on Least Square Support Vector Machine (LSSVM)
58
and the Markov model is proposed in [6]. The LSSVM get the parameters optimized by PSOGA
59
algorithm (Particle Swarm Optimization combined with Gravitational Search Algorithm).
60
Moreover, the reported results shown that the approach outperformed other techniques in terms of
61
forecasting quality. Other approach, which applies Support Vector Regression (SVR) with
62
parameters tuned by Genetic Algorithm is proposed in [7] to forecast wind speed and outperformed
63
the persistence model and well-known autoregressive models: Autoregressive (AR), Autoregressive
64
Moving Average (ARMA), and Autoregressive Integrated Moving Average (ARIMA). An
65
ensemble of Mixture Density Neural Networks developed by [8] to forecasting wind speed and
66
wind power with 72h ahead shown promising results. The model was applied to one wind farm
67
located in Taiwan and achieved good results when compared to Nonlinear Autoregressive Neural
68
Network with eXogenous inputs (NARX), ARMA and GARCH approaches.
AC C
EP
57
69 2
ACCEPTED MANUSCRIPT
In the context of time series forecasting, Artificial Neural Networks (ANNs) is the approach, which
71
have been numerous types of new models proposed [9-13], such as Radial Basis Function neural
72
network (RBF), the Multi-Layer Perceptron network (MLP) and the Support Vector Machine
73
(SVM). In addition, approaches applying hybrid models also have been grown fast due to best
74
results presented. As example can be cited [14], which applies a hybrid model that combines the
75
Extreme Learning Machine (ELM), the Ljung-Box Q-test (LBQ) and the Seasonal Autoregressive
76
Integrated Moving Average (SARIMA) on wind forecasting. Other example presented by [15],
77
applies ANNs, and Kalman filter (KF) based on the ARIMA model.
SC
RI PT
70
78
Various authors recently have pointed out the efficiency and quality of machine learning approach
80
to perform wind speed forecasting. In the present work, we intend to show a new hybrid forecasting
81
approach applied to wind speed forecast. The accuracy of the method is analyzed in terms of Mean
82
Squared Error (MSE), Root Mean Squared Error (RMSE), Mean Absolute Error (MAE) and
83
Standard Deviation (STD). According to the analysis of the obtained results, the proposed hybrid
84
model show better results on forecasting wind speed than ANFIS-FCM and Group Method Data
85
Handling (GMDH).
86
TE D
M AN U
79
The remainder of this paper is organized as follows. Section 2 describes Singular Spectrum
88
Analysis (SSA) theoretical approach and common uses for trend extraction and time series
89
decomposing. In addition, this section also presents the ANFIS and GMDH theoretical foundation.
90
Section 3 presents the proposed method. Section 4 presents the data description, the experimental
91
setup, and our numerical results. After, Section 5 contains the conclusions of this work.
92
2 Background
93
Singular Spectrum Analysis (SSA) is a technique of time series analysis and forecasting [16]. SSA
94
aims at decomposing the original series into a sum of a small number of interpretable components
95
such as a slowly varying trend, oscillatory components and a ‘structure less’ noise [17]. It is based
96
on the singular value decomposition (SVD) of a specific matrix constructed upon the time series.
AC C
EP
87
3
ACCEPTED MANUSCRIPT
Neither a parametric model nor stationary-type conditions have to be assumed for the time series.
98
This makes SSA a model-free method and hence enables SSA to have a wide range of applicability
99
[18]. In this section, we present the foundation of SSA and ANFIS techniques, with detailed
100
presentation, explore SSA foundations to facilitate the clear understanding. Regarding the ANFIS
101
technique, due the popularity and numerous works found in literature applying it, only a general
102
idea of the ANFIS framework is presented in this work. A briefly introduction of the GMDH
103
concepts also is provided.
104
2.1 Singular Spectrum Analysis - SSA
105
Singular Spectrum Analysis was published for the first time in 1986, but this approach became
106
popular only after the book published in 2000 by [17]. The SSA performs four steps:
110 111 112 113
SC M AN U
109
1. Embedding step: the one-dimensional time series is represented as a multidimensional series whose dimension is called the window length. The multidimensional time series forms the trajectory matrix. The only parameter of this step is the window length . Usually has value as < , where is the time series length;
TE D
107 108
RI PT
97
2. SVD step: is the singular value decomposition of the trajectory matrix into a sum of rankone bi-orthogonal matrices;
3. Grouping step: corresponds to splitting the matrices, computed at the SVD step, into
115
several groups and summing the matrices within each group. The result of the step is a
116
representation of the trajectory matrix as a sum of several resultant matrices;
AC C
EP
114
117
4. Diagonal averaging step: It is a linear operation and maps the trajectory matrix of the
118
initial series into the initial series itself. In this way, we obtain a decomposition of the initial
119
series into several additive components.
120 121
The Embedding and SDV steps form the Decomposition Stage and the last two steps form the
122
Reconstruction Stage. The power of the SSA technique lays on the analysis, as example, if we are
123
not interested in certain aspects of the components of the time series; we subsume them under the
124
noise components. The mathematical procedure to perform the two stages described below, 4
ACCEPTED MANUSCRIPT
125 126
considering a real-valued time series = = ( , … , ) of length . The full detailed explanation can be found in [17], please refers to it for detailed mathematical foundation.
RI PT
127 128 129
2.1.1 Decomposition Stage
130 131
The Decomposition Stage starts with the Embedding Step, which maps the original time series to a
132
= ( , … , ) , which have dimension and where 1 ≤ ≤ . The matrix formed by these vectors is the trajectory matrix presented in Equation 1
135
139 140 141 142 143 144 145 146 147 148 149 150 151 152
$
% ⋮
$
%
& ⋮
… … + … * ⋱ ⋮ … )
(1)
The matrix is a Hankel matrix, it means, has equal elements on the diagonals + - =
TE D
138
# = " $ ⋮ !
./012302. The next step is the singular value decomposition (SVD) of the trajectory matrix into a sum of biorthogonal elementary matrices of rank one. Take 4 = and the eigenvalues of 4 in decreasing order of magnitude (5 ≥ 5 ≥ ⋯ ≥ 5 ≥ 0), also the orthonormal system of the
eigenvectors denoted by 9 , … , 9 corresponding these eigenvalues. Let : = max>, 1?.ℎ 5 >
EP
137
= : … : = ( ), ,
0B = C30 , so if we denote D = 9 /F5 and = F5 9 D , where = 1, … , :, then the SVD of the trajectory matrix can be write as the sum
AC C
136
SC
134
M AN U
133
sequence of multidimensional lagged vectors of size by forming = − + 1 lagged vectors
= + + ⋯ + G
(2)
This work applies the useful property of Equation 2 to evaluate the contribution of each component
in the expansion of the original time series. It can be done as follows: among all the matrices
(H) of rank C < : (: is the number of 5 > 0), the matrix ∑H provides the best approximation to the trajectory matrix , it means that the Frobenius norm is minimum ‖ − H ‖ℳ . In practical 5
ACCEPTED MANUSCRIPT
terms, it means that the components contribution to form the original time series can be measure by
154
Equation 3. Therefore, we can decide which components select to perform the Reconstruction Stage
155
by selecting the first values sorted in decreasing order of their contribution
156
RI PT
153
%./02CM = ∑P
NO
OQR NO
157 158
2.1.2 Reconstruction Stage
SC
159
(3)
160 161
The Reconstruction stage is the following step after the expansion described in Equation 2. It starts
162
Now, let T = V , … , W X, so the matrix Y corresponding to the group T is defined as Y = +
164 165
⋯ + W , and these matrices are computed for T = T , … , TU . By applying the expansion described in Equation 2 leads to the decomposition
= Y + Y + ⋯ + YU
170 171 172 173 174
175
176
(4)
In this final step, the matrices Y , … , YU are converted into a new series of length , through a process called diagonal averaging. For step in this procedure, let Z be a matrix with elements
[ , where 1 ≤ ≤ and 1 ≤ - ≤ . Then set ∗ = min (, ), ∗ = max (, ) and = + −
EP
169
∗ ∗ 1. Finally, let [ = [ if < otherwise [ = [ . The diagonal averaging procedure transforms
the matrix Z into the time series [ , … , [ using the formula
AC C
168
TE D
166 167
M AN U
163
with the grouping procedure partitions the set of indices 1, … , : into S disjoint subsets T , … , TU .
_ ∗ ∗ c _ ∑U [U,_U d/C 0 ≤ < , a ∗ ∗ d/C ∗ ≤ < ∗ , [_ = ∗ ∑U [U,_U b ∗ a ∑ [ ∗ d/C ∗ ≤ < ` _ U _ ∗ U,_U
(5)
6
ACCEPTED MANUSCRIPT
179
equivalent to the decomposition of the initial series , … , into the sum of reconstructed S series
180
g = ∑U f (0 = 1,2, … , ) _
(_)
181 182
RI PT
178
(_) (_) Equation 5 applied a resultant matrix Y_ , produces the series e (_) = ( f , … , f ). This is
(6)
SC
177
2.2 Adaptive Neuro Fuzzy Inference System - ANFIS
184
Adaptive Neuro Fuzzy Inference System (ANFIS) is able to combine a fuzzy system’s ability to
185
model a reasoning process and to handle uncertainty, with the learning ability and adaptivity of a
186
neural network [19]. The ANFIS is a type of artificial neural network that is based on Takagi–
187
Sugeno fuzzy and it uses a hybrid learning algorithm. The main structure consists of five layers,
188
which can be called inputs, if part, rules and normalization, then part and output. For simplification, it is assumed that the framework of ANFIS has two inputs , [ and one output i as presented in
TE D
189
M AN U
183
190
Figure 1. Then, the corresponding rule set with two fuzzy if–then rules for a first-order Sugeno
191
consequent parameters. l , m and l , m are the linguistic labels. Each layer is described below.
193
AC C
194
EP
192
fuzzy model can be expressed as shown in Equation 7, where j , k , C and j , k , C are the
195 196
Figure 1 – ANFIS structure for a two-input , [ and one output i [20] 7
ACCEPTED MANUSCRIPT
nopq r: Td 1 l 30: [ 1 m 2ℎs0 i = j + k [ + C nopq t: Td 1 l 30: [ 1 m 2ℎs0 i = j + k [ + C
197 198
201
Layer 1: This is the fuzzification layer. Every node in this layer is an adaptive node with a node output defined by u , = vwO ( ) d/C = 1,2
202 203
u , = vxOyz ([)
204
206 207 208
output is the product of all the incoming signals. The output of this layer is given by u, = ~ = vwO ( )vxO ([) d/C = 1,2
210
213
the normalized firing strength as
217 218 219 220 221 222 223
O
R z
d/C = 1,2
(11)
Layer 4: This layer is the defuzzification layer. Every node in this layer is an adaptive node with a
AC C
216
u$, = ~ =
EP
214 215
(10)
Layer 3: This layer is the normalization layer. The 2ℎ node of this layer, labeled as , calculates
TE D
212
(9)
Layer 2: This layer is the rule layer. The nodes in this layer are fixed nodes, labeled as Π, whose
209
211
(8)
d/C = 3,4
M AN U
205
RI PT
200
SC
199
(7)
node function
u%, = ~ d = ~ (j + k [ + C )
(12)
where j , k , C is the parameter set of this node. Layer 5: This is the output layer. The single node in this layer is a fixed node labeled as Σ, which computes the overall output as the summation of all incoming signals. Overall output is given by
224 8
ACCEPTED MANUSCRIPT
u&, = ∑ ~ d =
225
∑O O O ∑O O
(13)
2.3 GMDH
227
The GMDH was introduced by [21] in the 1966’s. The main idea of the GMDH algorithm is
228
RI PT
226
searching for optimal structure within the space of multipolynomial functions d: g → , which it
realizes as a multilayered polynomial network. An important characteristic of GMDH it is to be a
230
self-organizing neural network model. This means that in these networks, the most important input
231
variables, number of layers, neurons in hidden layers and optimal model structure are determined
232
automatically through an iterative process, in which the model structure is modified in order to find
233
the best data prediction [22-24].
234
M AN U
SC
229
The GMDH algorithms start by creating simple polynomials that roughly approximate the aimed
236
systems and incrementally expand the complexity of the polynomials in order to create models that
237
are more accurate [25]. The connections between input and output variables can be expressed by
238
Volterra functional series, discrete analog of which is Kolmogorov–Gabor polynomial [25], as
239
presented in Equation 14
242 243 244 245
[ = 3 + ∑ 3 + ∑ ∑ 3 + ∑ ∑ ∑_ 3 _ _
EP
241
(14)
where ( , , … , ) is the input vector of variables, (3 , 3 , … , 3 ) is the vector of coefficients or weights and is the number of input variables.
AC C
240
TE D
235
246
3 Proposed hybrid SSA-ANFIS approach
247
The proposed hybrid model combines SSA and ANFIS. In literature can be found one hybrid model
248
proposed by [26], which resembles our approach, but the cited model differ by the hybridization
249
form applied in ANFIS and SSA. In the cited work, the authors apply SSA to decompose a daily 9
ACCEPTED MANUSCRIPT
250
average wind power time series, obtaining (0 − 9) components. Each component has the upper and
251
lower bound determined by Firefly Algorithm combined to Subtractive Cluster, and the final step
252
bounds on -step-ahead horizon by applying (0 − 9) ANFIS structures. Finally, the time series
254
RI PT
253
comprises to do the forecasting of the (0 − 9) time series components, also the upper and lower reconstruction is performed by summing the (0 − 9) components. The cited work was developed to forecast the uncertainty of the average wind power in daily basis, besides the complexity of the
256
model, the authors reported promising results.
257
The main context of SSA application is slightly different in our approach; where the main purpose
258
of SSA is, decompose the wind speed series into additive components. After this, we evaluate the
259
analysis of the contribution of each component, discarding those components that have no
260
significant contribution to the original signal recovery, and then the reconstruction of the original
261
series is made by grouping the components in two major components: trend plus harmonic and
262
noise components. These two components are the input for two different ANFIS model, which
263
performs one-step-ahead forecasting for each one.
M AN U
SC
255
264
266
Figure 2 shows the main process of the proposed model, where the trend plus harmonic component
TE D
265
is grouped as ( ) , and the noise component as () , both are applied as input for the two instances of ANFIS algorithm. We always implicitly assume that the wind speed series is a sum of several
268
simpler series. The main purpose is training two ANFIS models, one for each component.
269
Therefore, the goal is at the end of the training process, have a final model able to forecasting each
270
one of these components and their structure. Based on aforementioned definition, the original time
271
series can be obtained by summing these two components, the final stage is the summation of the
272
forecasted trend plus harmonic and noise components to get one-step-ahead forecasted value for the
273
original time series.
AC C
274
EP
267
275
Direct Recursive forecasting procedure, as presented in [24], is the approach applied to perform k-
276
step-ahead forecasting in this work. It means that we use in the recursive way the one-step-ahead
277
forecasting method, and each forecasted value is added to the original time series, so the new time
278 279
series length become + 1. After forecasting, the decomposition stage is performed again, also all
the algorithm steps, as presented in Figure 2. The entire procedure need be repeated after one-step10
ACCEPTED MANUSCRIPT
281 282 283
ahead value be forecasted, until perform all the k-step-ahead forecasting. At the end of the forecasting process, the length of the time series resulting is equals to + .
The SSA has only two parameters that need to be set, the first parameter is the window length ,
RI PT
280
and the second parameter is the way of grouping components. In [22] can be found explanation
285
about how to do the parameters choice by different approaches. A simple way to do this selection
286
can be done by taking into account the length, seasonality and periodicity of the time series. In this
287
work, we will take in consideration this former procedure, and in the next section provides detailed
288
explanation.
SC
284
TE D
Load the Time Series ; Get the Time Series length ; Plot the Time Series to visualize possible trend and seasonality; Set window length L by applying the Equation (15); Apply SSA-Decomposition Stage (Embedding and SVD step); Evaluate the contribution of each component by Eq.(3); Apply SSA-Reconstruction Stage (Grouping and Diagonal Averaging step) to get the two components: ( ) and () ; Apply the desired Lag (usually all lags containing autocorrelation values above a given threshold (0.2 as example) is selected); Split the ( ) in Training and Test Set and apply ANFIS algorithm; Split the () in Training and Test Set and apply ANFIS algorithm; Evaluated the model quality, if necessary improve the ANFIS parameters to achieve the minimal error desired and retrain the ANFIS Do the one-step-ahead forecast applying all ( ) data into ANFIS model; Do the one-step-ahead forecast applying all () data into ANFIS model; ( ) () Reconstruct the Time Series ( ) = ( ) + ( ) Complete the Time Series with the forecasted value; Update the new Time Series length = g If k-step-forecasting then While (t < k) Apply SSA-Decomposition Stage (Embedding and SVD step) to the new Time Series with length g ; Apply SSA-Reconstruction Stage (Grouping and Diagonal Averaging step) to get the two components: ( ) and () ;
EP
291
The pseudo-code presented below introduces the main idea of the proposed approach.
AC C
290
M AN U
289
11
ACCEPTED MANUSCRIPT
model;
Apply the Lag; Do the one-step-ahead forecast applying all ( ) data into ANFIS
Do the one-step-ahead forecast applying all () data into ANFIS
Reconstruct the Time Series ( ) = ( ) + ( ) Complete the Time Series with the forecasted value; Update the new Time Series length = g End While End if ()
RI PT
( )
SC
model;
Pseudo Code 1 – Hybrid SSA-ANFIS
AC C
EP
TE D
M AN U
292
12
AC C
EP
TE D
M AN U
SC
RI PT
ACCEPTED MANUSCRIPT
13
293
ACCEPTED MANUSCRIPT
Figure 2 – Proposed approach hybrid SSA-ANFIS
294
4 Numerical Experiments
296
As previously mentioned, the aim of this paper is to develop an approach to perform a wind speed
297
forecasting by applying a hybrid strategy, which allows the hybrid model to learn the structure of
298
trend and harmonic behavior of the original wind time series. Therefore, the main goal is allow
299
better energy production by uncertainty wind speed reducing. In this section, the potential benefit of
300
developing this strategy is analyzed from a quantitative point of view.
301
4.1 SSA - Parameters choice and Wind Time Series
305 306 307 308 309 310 311 312 313 314
SC
M AN U
which help to capture the periodic components with period = ~ℎsCs
TE D
304
grouping components. The choice of considering wind time series consider the equation (15), = 240 d 2ℎs 0?SMsC /d 13Sjs1 jsC ℎ/?C 0 > 1 = 24 d 2ℎs 0?SMsC /d 13Sjs1 jsC ℎ/?C 0 = 1
(15)
where is the length of the time series, 0 is the number of samples observed on interval within one
EP
303
The two parameters of SSA as discussed in Section 3 are the window length and the way of
hour. The result of this equation must result on that need be an integer, if it is not, a rounded
AC C
302
RI PT
295
process must be applied. The only restriction that must be observed, as presented in [17], must lays into a close interval given by 2 ≤ ≤
. As example, if the time series has length = 1000
and the sampling occurs at each 10 minutes, then we have 6 samples per hour; it means that 0 = 6 , in addition =
%
=
%%
= 6.94 ≅ 7, which can help capture the seasonality for multiple
periods, if the series has multiple periods.
315 316
The way of grouping, the last parameter that need be choose to perform the SSA algorithm, is done
317
by applying the Equation 3. The result of this equation is the contribution rate of each component 14
ACCEPTED MANUSCRIPT
into the original time series formation. The Table 1 presents an example of this result, when a time
319
series of 1000 data points of wind speed measured in a Wind Farm located in Rio Grande do Norte,
320
Northeast of Brazil, sampled by each 10 min on the first week of May 2016, is decomposed in 7
321
components.
RI PT
318
322 323 324
Table 1 – Component Contribution for original Time Series
326 327
Contribution [%] 99.36685 0.32204 0.13847 0.06714
Component #5 #6 #7
M AN U
Component #1 #2 #3 #4
SC
325
Contribution [%] 0.04531 0.03554 0.02465
Per the results of Table 1, the components 1 to 4 are the best choice to group them into ( ) component, mostly due to the level of contribution on the reconstruction of the original time series
329
as presented in Figure 3.
AC C
EP
330
TE D
328
331 332
Figure 3 – Contribution of each component to recovery the original time series 15
ACCEPTED MANUSCRIPT
333
To select which components can be grouped or discarded during the reconstruction step, the main
335
procedure described in [17] consists by analyzing the decreasing behavior of the singular values for
336
each component; a pure noise series produces a slowly decreasing sequence of singular values.
337
Because the axis scale, this slow variation cannot be identified after the second component in Figure
338
3, but the same information of Figure 3 are presented on numerical values form in Table 1.
339
Therefore, after the forth component, we can identify that moving from the fifth to the next
340
component; the decreasing behavior is very slow (decaying less from 0.01).
341
Considering the aforementioned procedure, the four initial components were grouped as ( ) ,
343 344
SC
M AN U
342
RI PT
334
representing all information related to the trend and harmonic time series behavior. After the forth component, all remaining information were grouped as noise and represented as () . In practical
terms, if the percentage of contribution is less than 0.1%, considering the insights found in [17],
346
which states that after this value a slowly decreasing sequence of singular values is observed; the
347
remaining components can be grouped as noise without loss information.
348
In this work, the wind time series have length = 1000 and the sampling occurs at each 10
351 352 353 354
minutes, then ≅ 7. For the grouping step, the components 1 to 4 was choose to form the
component ( ) and the remaining components will form the noise component () . Also, the 3 = 6 was adopted. The Figure 4 presents the original wind series and the Figure 5 the
EP
350
decomposed time series.
AC C
349
TE D
345
16
M AN U
SC
RI PT
ACCEPTED MANUSCRIPT
355
Figure 4 – Wind Time Series to be forecasted.
357
Figure 5 – Time Series decomposed: Components 1 to 4 reconstruced as ( ) and Noise as () .
AC C
EP
TE D
356
358
17
ACCEPTED MANUSCRIPT
4.2 ANFIS - Parameters choice
360
The ANFIS applied is the variant called Sugeno ANFIS-FCM (FCM - Fuzzy c-means for short).
361
Fuzzy c-means is a data clustering technique wherein each data point belongs to a cluster to some
362
degree that is specified by a membership grade. It provides a method that shows how to group data
363
points that populate some multidimensional space into a specific number of different clusters. The
364
number of clusters determine the number of rules and membership functions in the generated FIS.
365
The parameters was chosen by trial and error, also the initial set of values followed the
366
recommendation found in [19].
SC
M AN U
367 368
RI PT
359
Table 2 – Sugeno ANFIS – FCM parameters Value 20 4 200 1e-5
Parameter Max. # Epochs Error Goal Initial Step Size Step Decrease Rate Step Increase Rate
Value 200 0 0.01 0.9 1.1
TE D
Parameter Number of Cluster Partition Matrix Exponent Max. Number of Iterations Minimal Improvement
4.3 Training, Validation and Test Sets
370
The time series presented in Figure 4 must be split in three subsets before performing the model
371
fitting. The subset of training was composed of 70% of original length and 15% composes the
372
validation set of ANFIS one-step-ahead model. The remaining 15% was applied as test set to
compare the quality of the -steps-ahead forecasted values. Also, after the training and validation
AC C
373
EP
369
374
process, the 85% of the original time series is presented to the ANFIS model to evaluate the model
375
performance. The quality indicators to measure the trained and the validation for the one-step-ahead
376
forecasting model are the Mean Squared Error (MSE), Root Mean Squared Error (RMSE), Mean
377 378 379 380
Absolute Error (MAE), Standard Deviation (StD) and Coefficient of Determination ( ). Hgg = ∑ − H ) (
(16)
381 18
ACCEPTED MANUSCRIPT
Hgg = ∑ − H ) (
383
385
R
¡¢O£O£¤ ¡¥¢¦ )z ∑ O OQR(O
= 1 − ∑
386
RI PT
Hgg l = ∑ − H
384
OQR(O
387
¡¢O£O£¤ ¡¥¢¦ )z O
§
(17)
(18)
(19)
SC
382
4.4 Results for Training and Validation Set
389
Figure 6, 7, and 8 present the results for the training, validation set and for both data, resulting in
390
M AN U
388
85% of the all data related to the component ( ) , when applied the hybrid SSA-ANFIS-FCM
391
model to forecast one-step-ahead.
392 393
The quality of the model was evaluated regarding the correlation factor for training, validation
397 398 399 400
TE D
396
Figure 9 presents these results. The hybrid SSA-ANFIS-FCM model learned on the training set and the correlation coefficient for the one-step-ahead forecasting achieved = 0.99472. The
EP
395
and for both data set merged as a single one (training and validation) related to Component ( ) .
correlation achieved on Forecasting model related to the validation set was lower than the correlation for the training set, and its value = 0.86626. The last quality measure was applied to evaluate the result when both data set as presented as one time series to the hybrid model, with 85% length, the value achieved was very good: = 0.97458.
AC C
394
401 402
The same procedure was applied to the component () . It worth remember this component has in
403
its structure only noise and few information about the original wind time series. The Figures 10 to
404
12 present the results for one-step-ahead forecasting.
405 406
The quality indices of the model was evaluated regarding the correlation factor for the training,
407
validation and for the combination of both data set related to Component () . The Figure 13 19
ACCEPTED MANUSCRIPT
408
presents these results. The hybrid SSA-ANFIS-FCM model learned the structure of the training set
409
set = 0.47734. The last performance analysis was done by presenting all data of both set to the
411
model, and the = 0.72985.
RI PT
410
and the correlation coefficient for the forecasted values achieved = 0.97158, for the validation
Figure 6 – Results for Training Set - 70% length component ( )
EP
414 415
AC C
413
TE D
M AN U
SC
412
20
TE D
Figure 7 – Results for Validation Set – 15% length component ( )
AC C
EP
416 417 418 419 420 421
M AN U
SC
RI PT
ACCEPTED MANUSCRIPT
422 423
21
ACCEPTED MANUSCRIPT
430 431 432
Figure 9 – coefficient for training, validation and for both data set merged as single one for
AC C
429
EP
TE D
M AN U
SC
425 426 427 428
Figure 8 – Results for both data set (training and validation merged as single one) for Component ( )
RI PT
424
Component ( )
22
433
Figure 10 – Results for Training Set – 70% length Component ()
AC C
EP
TE D
434
M AN U
SC
RI PT
ACCEPTED MANUSCRIPT
435 436 437 438
Figure 11 – Results for Validation Set – 15% length Component ()
23
M AN U
SC
RI PT
ACCEPTED MANUSCRIPT
439
TE D EP
441 442 443
Figure 12 – Results for both data set (training and validation merged as single one) for Component ()
AC C
440
24
444 445
TE D
M AN U
SC
RI PT
ACCEPTED MANUSCRIPT
Figure 13 – coefficient for training, validation and for both data set merged as single one for
EP
446 447
449
As presented in Figure 2 and in Pseudo-Code 1, after obtain the one-step-ahead forecast for each
component ( ) and () , we need to perform the summation of these two components to
AC C
448
Component ()
450
reconstruct the original time series and the continuation of the series with the forecasted values, it
451
means that
452
() ( )
= ( ) + ( ). The Figure 14 presents the result for the original time ( )
()
series reconstructed after forecasting one-step-ahead process. Therefore, Figure 15 presents the
453
coefficient for all data after performs one-step-ahead forecasting.
454
The results for -steps-ahead forecasting are presented in Figure 16. Where the forecasting was
455 456
performed in recursive way, it means the first value forecasted became an input for the next value in 25
ACCEPTED MANUSCRIPT
460 461 462
463
hours forecast horizon. The 15% remaining of the original time series was applied to test the quality of the forecasted values, and the = 0.97261 was the coefficient resulting when the forecasted time series was compared to the test data set.
RI PT
459
850 points, 3 = 6, sampled by each 10min, were applied to obtain 144 steps resulting in 24
ª(2 + 1) = d( (2), (2 − 1), (2 − 2), (2 − 3), (2 − 4), (2 − 5))
ª(2 + 2) = d( ª(2 + 1), (2), (2 − 1), (2 − 2), (2 − 3), (2 − 4))
ª(2 + 3) = d« ª(2 + 2), ª(2 + 1), (2), (2 − 1), (2 − 2), (2 − 3)¬ ⋮
ª(2 + 6) = d« ª(2 + 5), ª(2 + 4), ª(2 + 3), ª(2 + 2), ª(2 + 1), (2)¬
SC
458
the forecasting process as Equation 20. As input for the forecasting model, a wind time series with
(20)
M AN U
457
AC C
EP
TE D
464
465 466 467 468
Figure 14 – Final Time Series reconstructed after forecasting one-step-ahead each component ( ) 30: () by SSA-ANFIS-FCM model
26
EP
TE D
Figure 15 – for the reconstructed series after forecasting one-step-ahead by SSA-ANFIS-FCM model
AC C
469 470 471 472 473
M AN U
SC
RI PT
ACCEPTED MANUSCRIPT
27
EP
TE D
Figure 16 – Wind Time Series forecasting 144-step-ahead
AC C
474 475 476 477 478
M AN U
SC
RI PT
ACCEPTED MANUSCRIPT
28
Figure 17 – after forecasting 144-step-ahead by SSA-ANFIS-FCM model
TE D
479 480
M AN U
SC
RI PT
ACCEPTED MANUSCRIPT
4.5 Comparisons Results
482
In this section, we present comparisons related to the proposed model with the ANFIS-FCM and
483
GMDH. For keeping the validity of comparisons between hybrid SSA-ANFIS and the ANFISFCM, the same parameters were used as presented in Table 2. In addition, 3 = 6 for the Time
AC C
484
EP
481
485
Series as input for all models also was kept. For GMDH algorithm, the Table 3 presents the
486
parameters applied, which the values was obtained considering the recommendation found in [24-
487
and Figure 19 the results for one-step-ahead forecasting. Figures 20 and 21 the results for -steps-
488 489
25, 27] . Results related to ANFIS-FCM model are depicted in Figure 18, which presents values ahead forecasting by ANFIS-FCM model.
490 491 492
The results for GMDH model are presented in Figures 22 to 25. For the one-step-ahead forecasting
the = 0.93664, decreasing to 0.89661 for -steps-ahead forecasting. In Table 4 the quality 29
ACCEPTED MANUSCRIPT
indices for all approaches, is presented and through it can be confirmed that the Hybrid SSA-ANFIS
494
outperformed the other two approaches.
495
On the results for -steps-ahead forecasting related to ANFIS-FCM and GMDH depicted on Figure
496
RI PT
493
20 and 24, respectively, a lag is observed between the target values and those that were forecasted.
498
We believe the main cause is due to recursive process applied on forecasting model. As presented in
499
Equation 20, after the maximum lag value, all inputs for the model are previous outputs forecasted,
500
causing some instability. For the proposed hybrid model, due to the filtering process into SSA step,
501
the outputs became more stable achieving better results on the forecasting process.
SC
497
M AN U
502
Table 3 – GMDH parameters
503
EP
505
AC C
504
Value 25 5 0 0.7
TE D
Parameter Max. Number of Neurons in a Layer Maximum Number of Layers Selection Pressure () Train Ratio
30
508
Figure 18 – for ANFIS-FCM model one-step-ahead forecasting
EP
507
AC C
506
TE D
M AN U
SC
RI PT
ACCEPTED MANUSCRIPT
31
M AN U
SC
RI PT
ACCEPTED MANUSCRIPT
509 510
Figure 19 – One-step-forecasting applying ANFIS-FCM
AC C
EP
TE D
511
512 513
Figure 20 – 144-step-forecasting applying ANFIS-FCM 32
ACCEPTED MANUSCRIPT
M AN U
SC
RI PT
514
TE D
518
Figure 21 - after forecasting 144-steps-ahead by ANFIS-FCM model
EP
516 517
AC C
515
33
521
Figure 22 – for one-step-ahead forecasting GMDH model
EP
520
AC C
519
TE D
M AN U
SC
RI PT
ACCEPTED MANUSCRIPT
34
M AN U
SC
RI PT
ACCEPTED MANUSCRIPT
522 523
Figure 23 – One-step-forecasting applying GMDH
AC C
EP
TE D
524
525 526
Figure 24 – 144-steps-forecasting applying GMDH 35
ACCEPTED MANUSCRIPT
531 532
Figure 25 - after forecasting 144-steps-ahead by GMDH model
EP
529 530
Table 4 – Quality Indices for one and k-step-ahead forecasting Technique
AC C
528
TE D
M AN U
SC
RI PT
527
SSA-ANFIS-FCM
ANFIS-FCM
Quality Indices MAE MSE RMSE St.D R2 MAE MSE RMSE St.D
Values for One-step-ahead 0.24053 0.20449 0.45221 0.45238 0.96886 0.48044 0.78746 0.88739 0.88770
Values for -steps-ahead 0.20583 0.08776 0.29625 0.29646 0.97261 0.51360 0.46828 0.68431 0.66950 36
ACCEPTED MANUSCRIPT
GMDH
533
0.88125 0.46025 0.38117 0.61739 0.61774 0.93664
0.85195 0.44280 0.31985 0.56555 0.55554 0.89661
RI PT
R2 MAE MSE RMSE St.D R2
The results provided by SSA-ANFIS-FCM are very close to the real values, we can infer through it
535
that the model learned the structure and of the peaks and trend of the time series. Otherwise, the
536
ANFIS-FCM forecasted quality indices presented in Table 4 were worse than GMDH approach, but
537
when the SSA is combined to ANFIS-FCM, the improvement are really promising.
M AN U
SC
534
538
5 Conclusion
540
In this study, a novel hybrid modeling method is proposed for short-term wind speed forecasting by
541
integrating SSA to ANFIS-FCM approach. This hybrid model combines the power of components
542
extraction and reconstruction provided by SSA algorithm. Moreover, the ANFIS characteristics of
543
learning a nonlinear system behavior introduce the precision and stability of the proposed approach
544
in forecasting wind speed time series. Since the dynamic behavior is recursively updated in the
545
proposed SSA–ANFIS approach, when more than one-step-ahead forecast were performed, the
546
stochastic uncertainty and fluctuations of wind speed can be better accounted for by the presented
547
method than by applying only the ANFIS-FCM and GMDH strategies. The better performance of
EP
hybrid SSA-ANFIS-FCM in both one-step-ahead and -steps-ahead wind speed forecasting
AC C
548
TE D
539
549
indicates the robustness of the proposed approach. Future researches will comprise to derive a
550
hybrid model by applying SSA into GMDH algorithm; also, the SSA-ANFIS-FCM model will be
551
compared, aiming to achieve the high accuracy on wind power forecasting.
552
37
ACCEPTED MANUSCRIPT
References
554
[1] O. Kramer, F. Gieseke, B. Satzger. Wind energy prediction and monitoring with neural
555
computation. Neurocomputing, 109 (2013), pp. 84–93
556
RI PT
553
[2] S. Salcedo-Sanz, J. Rojo-Álvarez, M. Martínez-Ramón, G. Camps-Valls. Support vector
558
machines in engineering: an overview. Wiley Interdiscip. Rev. Data Min. Knowl. Discov. 4
559
(2014), pp. 234–267
SC
557
560
562
[3] J. Heinermann, O. Kramer, Machine learning ensembles for wind power prediction, Renewable
M AN U
561
Energy, Volume 89, April 2016, Pages 671-679, ISSN 0960-1481.
563
[4] J. Wang, J. Hu, A robust combination approach for short-term wind speed forecasting and
565
analysis – Combination of the ARIMA (Autoregressive Integrated Moving Average), ELM
566
(Extreme Learning Machine), SVM (Support Vector Machine) and LSSVM (Least Square SVM)
567
forecasts using a GPR (Gaussian Process Regression) model, Energy, Volume 93, Part 1, 15
568
December 2015, Pages 41-56, ISSN 0360-5442.
569
TE D
564
[5] Y. Zhao, L. Ye, Z. Li, X. Song, Y. Lang, J. Su, A novel bidirectional mechanism based on time
571
series model for wind power forecasting, Applied Energy, Volume 177, 1 September 2016,
572
Pages 793-803, ISSN 0306-2619.
AC C
573
EP
570
574
[6] Y. Wang, J. Wang, X. Wei, A hybrid wind speed-forecasting model based on phase space
575
reconstruction theory and Markov model: A case study of wind farms in northwest China,
576
Energy, Volume 91, November 2015, Pages 556-572, ISSN 0360-5442.
577 578
[7] G. Santamaría-Bonfil, A. Reyes-Ballesteros, C. Gershenson, Wind speed forecasting for wind
579
farms: A method based on support vector regression, Renewable Energy, Volume 85, January
580
2016, Pages 790-809, ISSN 0960-1481.
581 38
ACCEPTED MANUSCRIPT
[8] Z. Men, E. Yee, F.-S. Lien, D. Wen, Y. Chen, Short-term wind speed and power forecasting
583
using an ensemble of mixture density neural networks, Renewable Energy, Volume 87, Part 1,
584
March 2016, Pages 203-211, ISSN 0960-1481.
RI PT
582
585
[9] Y-L. Tu, T.-J. Chang, C.-L. Chen, Y.-J. Chang, Estimation of monthly wind power outputs of
587
WECS with limited record period using artificial neural networks, Energy Conversion and
588
Management, Volume 59, July 2012, Pages 114-121, ISSN 0196-8904.
SC
586
589
[10] K. G. Sheela, S.N. Deepa, Neural network based hybrid computing model for wind speed
591
prediction, Neurocomputing, Volume 122, 25 December 2013, Pages 425-429, ISSN 0925-2312.
592
M AN U
590
593
[11] Q. Cao, B. T. Ewing, Mark A. Thompson, Forecasting wind speed with recurrent neural
594
networks, European Journal of Operational Research, Volume 221, Issue 1, 16 August 2012,
595
Pages 148-154, ISSN 0377-2217.
596
[12] D.A. Fadare, The application of artificial neural networks to mapping of wind speed profile for
598
energy application in Nigeria, Applied Energy, Volume 87, Issue 3, March 2010, Pages 934-942,
599
ISSN 0306-2619.
TE D
597
EP
600
[13] J. Zhou, J. Shi, G. Li, Fine tuning support vector machines for short-term wind speed
602
forecasting, Energy Conversion and Management, Volume 52, Issue 4, April 2011, Pages 1990-
603
1998, ISSN 0196-8904.
604 605 606
AC C
601
[14] J. Wang, J. Hu, K. Ma, Y. Zhang, A self-adaptive hybrid approach for wind speed forecasting, Renewable Energy, Volume 78, June 2015, Pages 374-385, ISSN 0960-1481.
607 608 609
[15] O. B. Shukur, M. H. Lee, Daily wind speed forecasting through hybrid KF-ANN model based on ARIMA, Renewable Energy, Volume 76, April 2015, Pages 637-647, ISSN 0960-1481.
610
39
ACCEPTED MANUSCRIPT
611 612
[16] A.D.D. Zhigljavsky (eds) (1997) Principal components of time series: the “Caterpillar” method. St.Petersburg Press, St. Petersburg (in Russian).
613
615
[17] N. Golyandina, V. Nekrutkin, A. Zhigljavsky, Analysis of time series structure: SSA and related techniques, (2001) Chapman & Hall/CRC, 305 pgs.
616
618
[18] J.B. Elsner, A.A.Tsonis, (1996) Singular spectrum analysis: a new tool in time series analysis. Plenum, New York.
SC
617
619
621
[19] J. S. R. Jang, “ANFIS: adaptive-network-based fuzzy inference system,” IEEE Transactions on
M AN U
620
Systems, Man and Cybernetics, vol. 23, pp. 665-685, 1993.
622 623 624
[20] E. Akkaya, ANFIS based prediction model for biomass heating value using proximate analysis components, Fuel, Volume 180, 15 September 2016, Pages 687-693, ISSN 0016-2361.
625
627
[21] A. G. Ivakhnenko, Polynomial Theory of Complex Systems. IEEE Transactions on Systems,
TE D
626
Man, and Cybernetics, Vol. SMC-1, No. 4, pp. 364-378 (1971).
628
630
[22] A.G. Ivakhnenko, A.J. Müller, Self-organization of neuron nets with active neurons. Pattern Recognition Image Analysis 1994;4(2):177e88.
EP
629
631
633 634
[23] A.G. Ivakhnenko, D. Wunsch, Inductive sorting/out GMDH algorithms with polynomial
AC C
632
RI PT
614
complexity for active neurons of neural network. IEEE 1999;6 (99): 1169e73.
635
[24] M.G. De Giorgi, M. Malvoni, P.M. Congedo, Comparison of strategies for multi-step ahead
636
photovoltaic power forecasting models based on hybrid group method of data handling networks
637
and least square support vector machine, Energy, Volume 107, 15 July 2016, Pages 360-373,
638
ISSN 0360-5442.
639 640
[25] M. Dorn, A.L.S. Braga, C.H. Llanos, L.S. Coelho, A GMDH polynomial neural network-based
641
method to predict approximate three-dimensional structures of polypeptides, Expert Systems 40
ACCEPTED MANUSCRIPT
642
with Applications, Volume 39, Issue 15, 1 November 2012, Pages 12268-12279, ISSN 0957-
643
4174.
644
[26] Z. Zhang, Y. Song, F. Liu, J. Liu, Daily Average Wind Power Interval Forecasts Based on an
646
Optimal Adaptive-Network-Based Fuzzy Inference System and Singular Spectrum Analysis.
647
Sustainability 2016, 8, 125.
RI PT
645
648
[27] S.-Kwun Oh, W. Pedrycz, B.-Jun Park, Polynomial neural networks architecture: analysis and
650
design, In Computers & Electrical Engineering, Volume 29, Issue 6, 2003, Pages 703-725, ISSN
651
0045-7906, https://doi.org/10.1016/S0045-7906(02)00045-9.
M AN U
SC
649
AC C
EP
TE D
652
41
ACCEPTED MANUSCRIPT Research highlights •
An hybrid SSA-ANFIS-FCM approach for wind speed forecasting is proposed
•
The hybridization of SSA and ANFIS-FCM is an original idea on context of Machine Learning applied to wind speed forecast The model presents excellent results with stable performance
AC C
EP
TE D
M AN U
SC
RI PT
•
S0960-1481(17)31190-4
DOI:
10.1016/j.renene.2017.11.089
Reference:
RENE 9491
To appear in:
Renewable Energy
Received Date: 6 August 2016 Revised Date:
13 October 2017
Accepted Date: 29 November 2017
Please cite this article as: Moreno SR, dos Santos Coelho L, Wind speed forecasting approach based on Singular Spectrum Analysis and Adaptive Neuro Fuzzy Inference System, Renewable Energy (2018), doi: 10.1016/j.renene.2017.11.089. This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.
ACCEPTED MANUSCRIPT
RI PT
Sinvaldo Rodrigues Moreno1 and Leandro dos Santos Coelho1,2*
1 Department of Electrical Engineering, Federal University of Parana Zip code 81531-980, Curitiba, PR, Brazil, 2 Industrial and Systems Engineering Graduate Program (PPGEPS)
SC
Pontifical Catholic University of Parana (PUCPR) Imaculada Conceicao, 1155, Zip code 80215-901, Curitiba, PR, Brazil * Corresponding author Tel/fax numbers: (55) 41 32711345 Email: [email protected]
M AN U
3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31
Abstract: As a promising renewable energy source, wind power has environmental benefits, as well as economic and social ones. Due these characteristics, wind farm has grown fast in the last five years, and in some countries, it has already surpassed conventional sources, such as hydro and coal plants. However, owing to the uncertainty of wind speed, it is essential to build an accurate forecasting model for large-scale wind power penetration. This study proposes a hybrid approach that combines the Singular Spectrum Analysis (SSA), which rarely presents application in literature on wind speed forecasting, and a Computing Natural paradigm called Adaptive Neuro Fuzzy Inference System (ANFIS). The SSA decomposes the original wind speed into various components, so these components are pre-processed regarding to the level of original wind series information remained. The main components selected to reconstruct the original series have in their structure the information about trend and harmonic components. The final remaining components grouped are labeled as noise. The ANFIS model uses these two information to construct the model applied to forecasting the next wind speed value. In this way, the hybrid model can learn the trend and the harmonic structure of the wind time series. Experimental results show that prediction errors are significantly reduced using the proposed technique to perform 10min one-step-ahead and -step-ahead wind speed forecast.
TE D
2
Wind Speed Forecasting Approach Based on Singular Spectrum Analysis and Adaptive Neuro Fuzzy Inference System
EP
1
Keywords: wind speed, hybrid forecasting approach, wind farm, machine learning, singular spectrum
33
analysis, neuro fuzzy inference system, time series forecasting.
34
1 Introduction
35
The renewable energy sources, particularly wind power is becoming more popular because of
36
ecological concerns and continuously rising prices of rapidly depleting conventional energy
37
sources. However, several drawbacks complicate the large-scale penetration into large fully
38
connected systems, as the intermittent nature of wind power. Mostly due to wind speed related
39
uncertainties affecting the reliability of energy supply, reliable forecasts methods can effectively
AC C
32
1
ACCEPTED MANUSCRIPT
reduce the risk to the power system. Machine Learning methods have attracted more attention than
41
traditional methods (such as numerical weather predictions and statistical approaches) when the
42
forecasting horizon are comprised between 1 to 3 hours or steps ahead [1-3].
RI PT
40
43
Many approaches applying Natural Computing and Machine Learning concepts have been proposed
45
recently to predicted wind speed. The approaches include promising mathematical techniques such
46
as Support Vector Machine (SVM), Wavelet Transform Method, Artificial Neural Networks,
47
Adaptive Neuro Fuzzy Inference System (ANFIS), Extreme Learning Machine (ELM), among
48
others. Despite the quality of the results, each technique has its own advantages and limitations.
49
Recent work conducted by [4] presents a combination of ARIMA (Autoregressive Integrated
50
Moving Average) model, ELM, SVM and LSSVM (Least Square SVM) to build a model of inputs,
51
and applied it a GPR (Gaussian Process Regression) to perform short-term wind speed forecasting.
52
In [5] a novel bidirectional mechanism and a backward forecasting model based on ELM for wind
53
power time series forecasting have been introduced. This model, as reported by the authors, can
54
achieve good results when applied to forecast wind farm power production for six time horizons as
55
1 up to 6 h.
M AN U
TE D
56
SC
44
A hybrid intelligent forecasting model based on Least Square Support Vector Machine (LSSVM)
58
and the Markov model is proposed in [6]. The LSSVM get the parameters optimized by PSOGA
59
algorithm (Particle Swarm Optimization combined with Gravitational Search Algorithm).
60
Moreover, the reported results shown that the approach outperformed other techniques in terms of
61
forecasting quality. Other approach, which applies Support Vector Regression (SVR) with
62
parameters tuned by Genetic Algorithm is proposed in [7] to forecast wind speed and outperformed
63
the persistence model and well-known autoregressive models: Autoregressive (AR), Autoregressive
64
Moving Average (ARMA), and Autoregressive Integrated Moving Average (ARIMA). An
65
ensemble of Mixture Density Neural Networks developed by [8] to forecasting wind speed and
66
wind power with 72h ahead shown promising results. The model was applied to one wind farm
67
located in Taiwan and achieved good results when compared to Nonlinear Autoregressive Neural
68
Network with eXogenous inputs (NARX), ARMA and GARCH approaches.
AC C
EP
57
69 2
ACCEPTED MANUSCRIPT
In the context of time series forecasting, Artificial Neural Networks (ANNs) is the approach, which
71
have been numerous types of new models proposed [9-13], such as Radial Basis Function neural
72
network (RBF), the Multi-Layer Perceptron network (MLP) and the Support Vector Machine
73
(SVM). In addition, approaches applying hybrid models also have been grown fast due to best
74
results presented. As example can be cited [14], which applies a hybrid model that combines the
75
Extreme Learning Machine (ELM), the Ljung-Box Q-test (LBQ) and the Seasonal Autoregressive
76
Integrated Moving Average (SARIMA) on wind forecasting. Other example presented by [15],
77
applies ANNs, and Kalman filter (KF) based on the ARIMA model.
SC
RI PT
70
78
Various authors recently have pointed out the efficiency and quality of machine learning approach
80
to perform wind speed forecasting. In the present work, we intend to show a new hybrid forecasting
81
approach applied to wind speed forecast. The accuracy of the method is analyzed in terms of Mean
82
Squared Error (MSE), Root Mean Squared Error (RMSE), Mean Absolute Error (MAE) and
83
Standard Deviation (STD). According to the analysis of the obtained results, the proposed hybrid
84
model show better results on forecasting wind speed than ANFIS-FCM and Group Method Data
85
Handling (GMDH).
86
TE D
M AN U
79
The remainder of this paper is organized as follows. Section 2 describes Singular Spectrum
88
Analysis (SSA) theoretical approach and common uses for trend extraction and time series
89
decomposing. In addition, this section also presents the ANFIS and GMDH theoretical foundation.
90
Section 3 presents the proposed method. Section 4 presents the data description, the experimental
91
setup, and our numerical results. After, Section 5 contains the conclusions of this work.
92
2 Background
93
Singular Spectrum Analysis (SSA) is a technique of time series analysis and forecasting [16]. SSA
94
aims at decomposing the original series into a sum of a small number of interpretable components
95
such as a slowly varying trend, oscillatory components and a ‘structure less’ noise [17]. It is based
96
on the singular value decomposition (SVD) of a specific matrix constructed upon the time series.
AC C
EP
87
3
ACCEPTED MANUSCRIPT
Neither a parametric model nor stationary-type conditions have to be assumed for the time series.
98
This makes SSA a model-free method and hence enables SSA to have a wide range of applicability
99
[18]. In this section, we present the foundation of SSA and ANFIS techniques, with detailed
100
presentation, explore SSA foundations to facilitate the clear understanding. Regarding the ANFIS
101
technique, due the popularity and numerous works found in literature applying it, only a general
102
idea of the ANFIS framework is presented in this work. A briefly introduction of the GMDH
103
concepts also is provided.
104
2.1 Singular Spectrum Analysis - SSA
105
Singular Spectrum Analysis was published for the first time in 1986, but this approach became
106
popular only after the book published in 2000 by [17]. The SSA performs four steps:
110 111 112 113
SC M AN U
109
1. Embedding step: the one-dimensional time series is represented as a multidimensional series whose dimension is called the window length. The multidimensional time series forms the trajectory matrix. The only parameter of this step is the window length . Usually has value as < , where is the time series length;
TE D
107 108
RI PT
97
2. SVD step: is the singular value decomposition of the trajectory matrix into a sum of rankone bi-orthogonal matrices;
3. Grouping step: corresponds to splitting the matrices, computed at the SVD step, into
115
several groups and summing the matrices within each group. The result of the step is a
116
representation of the trajectory matrix as a sum of several resultant matrices;
AC C
EP
114
117
4. Diagonal averaging step: It is a linear operation and maps the trajectory matrix of the
118
initial series into the initial series itself. In this way, we obtain a decomposition of the initial
119
series into several additive components.
120 121
The Embedding and SDV steps form the Decomposition Stage and the last two steps form the
122
Reconstruction Stage. The power of the SSA technique lays on the analysis, as example, if we are
123
not interested in certain aspects of the components of the time series; we subsume them under the
124
noise components. The mathematical procedure to perform the two stages described below, 4
ACCEPTED MANUSCRIPT
125 126
considering a real-valued time series = = ( , … , ) of length . The full detailed explanation can be found in [17], please refers to it for detailed mathematical foundation.
RI PT
127 128 129
2.1.1 Decomposition Stage
130 131
The Decomposition Stage starts with the Embedding Step, which maps the original time series to a
132
= ( , … , ) , which have dimension and where 1 ≤ ≤ . The matrix formed by these vectors is the trajectory matrix presented in Equation 1
135
139 140 141 142 143 144 145 146 147 148 149 150 151 152
$
% ⋮
$
%
& ⋮
… … + … * ⋱ ⋮ … )
(1)
The matrix is a Hankel matrix, it means, has equal elements on the diagonals + - =
TE D
138
# = " $ ⋮ !
./012302. The next step is the singular value decomposition (SVD) of the trajectory matrix into a sum of biorthogonal elementary matrices of rank one. Take 4 = and the eigenvalues of 4 in decreasing order of magnitude (5 ≥ 5 ≥ ⋯ ≥ 5 ≥ 0), also the orthonormal system of the
eigenvectors denoted by 9 , … , 9 corresponding these eigenvalues. Let : = max>, 1?.ℎ 5 >
EP
137
= : … : = ( ), ,
0B = C30 , so if we denote D = 9 /F5 and = F5 9 D , where = 1, … , :, then the SVD of the trajectory matrix can be write as the sum
AC C
136
SC
134
M AN U
133
sequence of multidimensional lagged vectors of size by forming = − + 1 lagged vectors
= + + ⋯ + G
(2)
This work applies the useful property of Equation 2 to evaluate the contribution of each component
in the expansion of the original time series. It can be done as follows: among all the matrices
(H) of rank C < : (: is the number of 5 > 0), the matrix ∑H provides the best approximation to the trajectory matrix , it means that the Frobenius norm is minimum ‖ − H ‖ℳ . In practical 5
ACCEPTED MANUSCRIPT
terms, it means that the components contribution to form the original time series can be measure by
154
Equation 3. Therefore, we can decide which components select to perform the Reconstruction Stage
155
by selecting the first values sorted in decreasing order of their contribution
156
RI PT
153
%./02CM = ∑P
NO
OQR NO
157 158
2.1.2 Reconstruction Stage
SC
159
(3)
160 161
The Reconstruction stage is the following step after the expansion described in Equation 2. It starts
162
Now, let T = V , … , W X, so the matrix Y corresponding to the group T is defined as Y = +
164 165
⋯ + W , and these matrices are computed for T = T , … , TU . By applying the expansion described in Equation 2 leads to the decomposition
= Y + Y + ⋯ + YU
170 171 172 173 174
175
176
(4)
In this final step, the matrices Y , … , YU are converted into a new series of length , through a process called diagonal averaging. For step in this procedure, let Z be a matrix with elements
[ , where 1 ≤ ≤ and 1 ≤ - ≤ . Then set ∗ = min (, ), ∗ = max (, ) and = + −
EP
169
∗ ∗ 1. Finally, let [ = [ if < otherwise [ = [ . The diagonal averaging procedure transforms
the matrix Z into the time series [ , … , [ using the formula
AC C
168
TE D
166 167
M AN U
163
with the grouping procedure partitions the set of indices 1, … , : into S disjoint subsets T , … , TU .
_ ∗ ∗ c _ ∑U [U,_U d/C 0 ≤ < , a ∗ ∗ d/C ∗ ≤ < ∗ , [_ = ∗ ∑U [U,_U b ∗ a ∑ [ ∗ d/C ∗ ≤ < ` _ U _ ∗ U,_U
(5)
6
ACCEPTED MANUSCRIPT
179
equivalent to the decomposition of the initial series , … , into the sum of reconstructed S series
180
g = ∑U f (0 = 1,2, … , ) _
(_)
181 182
RI PT
178
(_) (_) Equation 5 applied a resultant matrix Y_ , produces the series e (_) = ( f , … , f ). This is
(6)
SC
177
2.2 Adaptive Neuro Fuzzy Inference System - ANFIS
184
Adaptive Neuro Fuzzy Inference System (ANFIS) is able to combine a fuzzy system’s ability to
185
model a reasoning process and to handle uncertainty, with the learning ability and adaptivity of a
186
neural network [19]. The ANFIS is a type of artificial neural network that is based on Takagi–
187
Sugeno fuzzy and it uses a hybrid learning algorithm. The main structure consists of five layers,
188
which can be called inputs, if part, rules and normalization, then part and output. For simplification, it is assumed that the framework of ANFIS has two inputs , [ and one output i as presented in
TE D
189
M AN U
183
190
Figure 1. Then, the corresponding rule set with two fuzzy if–then rules for a first-order Sugeno
191
consequent parameters. l , m and l , m are the linguistic labels. Each layer is described below.
193
AC C
194
EP
192
fuzzy model can be expressed as shown in Equation 7, where j , k , C and j , k , C are the
195 196
Figure 1 – ANFIS structure for a two-input , [ and one output i [20] 7
ACCEPTED MANUSCRIPT
nopq r: Td 1 l 30: [ 1 m 2ℎs0 i = j + k [ + C nopq t: Td 1 l 30: [ 1 m 2ℎs0 i = j + k [ + C
197 198
201
Layer 1: This is the fuzzification layer. Every node in this layer is an adaptive node with a node output defined by u , = vwO ( ) d/C = 1,2
202 203
u , = vxOyz ([)
204
206 207 208
output is the product of all the incoming signals. The output of this layer is given by u, = ~ = vwO ( )vxO ([) d/C = 1,2
210
213
the normalized firing strength as
217 218 219 220 221 222 223
O
R z
d/C = 1,2
(11)
Layer 4: This layer is the defuzzification layer. Every node in this layer is an adaptive node with a
AC C
216
u$, = ~ =
EP
214 215
(10)
Layer 3: This layer is the normalization layer. The 2ℎ node of this layer, labeled as , calculates
TE D
212
(9)
Layer 2: This layer is the rule layer. The nodes in this layer are fixed nodes, labeled as Π, whose
209
211
(8)
d/C = 3,4
M AN U
205
RI PT
200
SC
199
(7)
node function
u%, = ~ d = ~ (j + k [ + C )
(12)
where j , k , C is the parameter set of this node. Layer 5: This is the output layer. The single node in this layer is a fixed node labeled as Σ, which computes the overall output as the summation of all incoming signals. Overall output is given by
224 8
ACCEPTED MANUSCRIPT
u&, = ∑ ~ d =
225
∑O O O ∑O O
(13)
2.3 GMDH
227
The GMDH was introduced by [21] in the 1966’s. The main idea of the GMDH algorithm is
228
RI PT
226
searching for optimal structure within the space of multipolynomial functions d: g → , which it
realizes as a multilayered polynomial network. An important characteristic of GMDH it is to be a
230
self-organizing neural network model. This means that in these networks, the most important input
231
variables, number of layers, neurons in hidden layers and optimal model structure are determined
232
automatically through an iterative process, in which the model structure is modified in order to find
233
the best data prediction [22-24].
234
M AN U
SC
229
The GMDH algorithms start by creating simple polynomials that roughly approximate the aimed
236
systems and incrementally expand the complexity of the polynomials in order to create models that
237
are more accurate [25]. The connections between input and output variables can be expressed by
238
Volterra functional series, discrete analog of which is Kolmogorov–Gabor polynomial [25], as
239
presented in Equation 14
242 243 244 245
[ = 3 + ∑ 3 + ∑ ∑ 3 + ∑ ∑ ∑_ 3 _ _
EP
241
(14)
where ( , , … , ) is the input vector of variables, (3 , 3 , … , 3 ) is the vector of coefficients or weights and is the number of input variables.
AC C
240
TE D
235
246
3 Proposed hybrid SSA-ANFIS approach
247
The proposed hybrid model combines SSA and ANFIS. In literature can be found one hybrid model
248
proposed by [26], which resembles our approach, but the cited model differ by the hybridization
249
form applied in ANFIS and SSA. In the cited work, the authors apply SSA to decompose a daily 9
ACCEPTED MANUSCRIPT
250
average wind power time series, obtaining (0 − 9) components. Each component has the upper and
251
lower bound determined by Firefly Algorithm combined to Subtractive Cluster, and the final step
252
bounds on -step-ahead horizon by applying (0 − 9) ANFIS structures. Finally, the time series
254
RI PT
253
comprises to do the forecasting of the (0 − 9) time series components, also the upper and lower reconstruction is performed by summing the (0 − 9) components. The cited work was developed to forecast the uncertainty of the average wind power in daily basis, besides the complexity of the
256
model, the authors reported promising results.
257
The main context of SSA application is slightly different in our approach; where the main purpose
258
of SSA is, decompose the wind speed series into additive components. After this, we evaluate the
259
analysis of the contribution of each component, discarding those components that have no
260
significant contribution to the original signal recovery, and then the reconstruction of the original
261
series is made by grouping the components in two major components: trend plus harmonic and
262
noise components. These two components are the input for two different ANFIS model, which
263
performs one-step-ahead forecasting for each one.
M AN U
SC
255
264
266
Figure 2 shows the main process of the proposed model, where the trend plus harmonic component
TE D
265
is grouped as ( ) , and the noise component as () , both are applied as input for the two instances of ANFIS algorithm. We always implicitly assume that the wind speed series is a sum of several
268
simpler series. The main purpose is training two ANFIS models, one for each component.
269
Therefore, the goal is at the end of the training process, have a final model able to forecasting each
270
one of these components and their structure. Based on aforementioned definition, the original time
271
series can be obtained by summing these two components, the final stage is the summation of the
272
forecasted trend plus harmonic and noise components to get one-step-ahead forecasted value for the
273
original time series.
AC C
274
EP
267
275
Direct Recursive forecasting procedure, as presented in [24], is the approach applied to perform k-
276
step-ahead forecasting in this work. It means that we use in the recursive way the one-step-ahead
277
forecasting method, and each forecasted value is added to the original time series, so the new time
278 279
series length become + 1. After forecasting, the decomposition stage is performed again, also all
the algorithm steps, as presented in Figure 2. The entire procedure need be repeated after one-step10
ACCEPTED MANUSCRIPT
281 282 283
ahead value be forecasted, until perform all the k-step-ahead forecasting. At the end of the forecasting process, the length of the time series resulting is equals to + .
The SSA has only two parameters that need to be set, the first parameter is the window length ,
RI PT
280
and the second parameter is the way of grouping components. In [22] can be found explanation
285
about how to do the parameters choice by different approaches. A simple way to do this selection
286
can be done by taking into account the length, seasonality and periodicity of the time series. In this
287
work, we will take in consideration this former procedure, and in the next section provides detailed
288
explanation.
SC
284
TE D
Load the Time Series ; Get the Time Series length ; Plot the Time Series to visualize possible trend and seasonality; Set window length L by applying the Equation (15); Apply SSA-Decomposition Stage (Embedding and SVD step); Evaluate the contribution of each component by Eq.(3); Apply SSA-Reconstruction Stage (Grouping and Diagonal Averaging step) to get the two components: ( ) and () ; Apply the desired Lag (usually all lags containing autocorrelation values above a given threshold (0.2 as example) is selected); Split the ( ) in Training and Test Set and apply ANFIS algorithm; Split the () in Training and Test Set and apply ANFIS algorithm; Evaluated the model quality, if necessary improve the ANFIS parameters to achieve the minimal error desired and retrain the ANFIS Do the one-step-ahead forecast applying all ( ) data into ANFIS model; Do the one-step-ahead forecast applying all () data into ANFIS model; ( ) () Reconstruct the Time Series ( ) = ( ) + ( ) Complete the Time Series with the forecasted value; Update the new Time Series length = g If k-step-forecasting then While (t < k) Apply SSA-Decomposition Stage (Embedding and SVD step) to the new Time Series with length g ; Apply SSA-Reconstruction Stage (Grouping and Diagonal Averaging step) to get the two components: ( ) and () ;
EP
291
The pseudo-code presented below introduces the main idea of the proposed approach.
AC C
290
M AN U
289
11
ACCEPTED MANUSCRIPT
model;
Apply the Lag; Do the one-step-ahead forecast applying all ( ) data into ANFIS
Do the one-step-ahead forecast applying all () data into ANFIS
Reconstruct the Time Series ( ) = ( ) + ( ) Complete the Time Series with the forecasted value; Update the new Time Series length = g End While End if ()
RI PT
( )
SC
model;
Pseudo Code 1 – Hybrid SSA-ANFIS
AC C
EP
TE D
M AN U
292
12
AC C
EP
TE D
M AN U
SC
RI PT
ACCEPTED MANUSCRIPT
13
293
ACCEPTED MANUSCRIPT
Figure 2 – Proposed approach hybrid SSA-ANFIS
294
4 Numerical Experiments
296
As previously mentioned, the aim of this paper is to develop an approach to perform a wind speed
297
forecasting by applying a hybrid strategy, which allows the hybrid model to learn the structure of
298
trend and harmonic behavior of the original wind time series. Therefore, the main goal is allow
299
better energy production by uncertainty wind speed reducing. In this section, the potential benefit of
300
developing this strategy is analyzed from a quantitative point of view.
301
4.1 SSA - Parameters choice and Wind Time Series
305 306 307 308 309 310 311 312 313 314
SC
M AN U
which help to capture the periodic components with period = ~ℎsCs
TE D
304
grouping components. The choice of considering wind time series consider the equation (15), = 240 d 2ℎs 0?SMsC /d 13Sjs1 jsC ℎ/?C 0 > 1 = 24 d 2ℎs 0?SMsC /d 13Sjs1 jsC ℎ/?C 0 = 1
(15)
where is the length of the time series, 0 is the number of samples observed on interval within one
EP
303
The two parameters of SSA as discussed in Section 3 are the window length and the way of
hour. The result of this equation must result on that need be an integer, if it is not, a rounded
AC C
302
RI PT
295
process must be applied. The only restriction that must be observed, as presented in [17], must lays into a close interval given by 2 ≤ ≤
. As example, if the time series has length = 1000
and the sampling occurs at each 10 minutes, then we have 6 samples per hour; it means that 0 = 6 , in addition =
%
=
%%
= 6.94 ≅ 7, which can help capture the seasonality for multiple
periods, if the series has multiple periods.
315 316
The way of grouping, the last parameter that need be choose to perform the SSA algorithm, is done
317
by applying the Equation 3. The result of this equation is the contribution rate of each component 14
ACCEPTED MANUSCRIPT
into the original time series formation. The Table 1 presents an example of this result, when a time
319
series of 1000 data points of wind speed measured in a Wind Farm located in Rio Grande do Norte,
320
Northeast of Brazil, sampled by each 10 min on the first week of May 2016, is decomposed in 7
321
components.
RI PT
318
322 323 324
Table 1 – Component Contribution for original Time Series
326 327
Contribution [%] 99.36685 0.32204 0.13847 0.06714
Component #5 #6 #7
M AN U
Component #1 #2 #3 #4
SC
325
Contribution [%] 0.04531 0.03554 0.02465
Per the results of Table 1, the components 1 to 4 are the best choice to group them into ( ) component, mostly due to the level of contribution on the reconstruction of the original time series
329
as presented in Figure 3.
AC C
EP
330
TE D
328
331 332
Figure 3 – Contribution of each component to recovery the original time series 15
ACCEPTED MANUSCRIPT
333
To select which components can be grouped or discarded during the reconstruction step, the main
335
procedure described in [17] consists by analyzing the decreasing behavior of the singular values for
336
each component; a pure noise series produces a slowly decreasing sequence of singular values.
337
Because the axis scale, this slow variation cannot be identified after the second component in Figure
338
3, but the same information of Figure 3 are presented on numerical values form in Table 1.
339
Therefore, after the forth component, we can identify that moving from the fifth to the next
340
component; the decreasing behavior is very slow (decaying less from 0.01).
341
Considering the aforementioned procedure, the four initial components were grouped as ( ) ,
343 344
SC
M AN U
342
RI PT
334
representing all information related to the trend and harmonic time series behavior. After the forth component, all remaining information were grouped as noise and represented as () . In practical
terms, if the percentage of contribution is less than 0.1%, considering the insights found in [17],
346
which states that after this value a slowly decreasing sequence of singular values is observed; the
347
remaining components can be grouped as noise without loss information.
348
In this work, the wind time series have length = 1000 and the sampling occurs at each 10
351 352 353 354
minutes, then ≅ 7. For the grouping step, the components 1 to 4 was choose to form the
component ( ) and the remaining components will form the noise component () . Also, the 3 = 6 was adopted. The Figure 4 presents the original wind series and the Figure 5 the
EP
350
decomposed time series.
AC C
349
TE D
345
16
M AN U
SC
RI PT
ACCEPTED MANUSCRIPT
355
Figure 4 – Wind Time Series to be forecasted.
357
Figure 5 – Time Series decomposed: Components 1 to 4 reconstruced as ( ) and Noise as () .
AC C
EP
TE D
356
358
17
ACCEPTED MANUSCRIPT
4.2 ANFIS - Parameters choice
360
The ANFIS applied is the variant called Sugeno ANFIS-FCM (FCM - Fuzzy c-means for short).
361
Fuzzy c-means is a data clustering technique wherein each data point belongs to a cluster to some
362
degree that is specified by a membership grade. It provides a method that shows how to group data
363
points that populate some multidimensional space into a specific number of different clusters. The
364
number of clusters determine the number of rules and membership functions in the generated FIS.
365
The parameters was chosen by trial and error, also the initial set of values followed the
366
recommendation found in [19].
SC
M AN U
367 368
RI PT
359
Table 2 – Sugeno ANFIS – FCM parameters Value 20 4 200 1e-5
Parameter Max. # Epochs Error Goal Initial Step Size Step Decrease Rate Step Increase Rate
Value 200 0 0.01 0.9 1.1
TE D
Parameter Number of Cluster Partition Matrix Exponent Max. Number of Iterations Minimal Improvement
4.3 Training, Validation and Test Sets
370
The time series presented in Figure 4 must be split in three subsets before performing the model
371
fitting. The subset of training was composed of 70% of original length and 15% composes the
372
validation set of ANFIS one-step-ahead model. The remaining 15% was applied as test set to
compare the quality of the -steps-ahead forecasted values. Also, after the training and validation
AC C
373
EP
369
374
process, the 85% of the original time series is presented to the ANFIS model to evaluate the model
375
performance. The quality indicators to measure the trained and the validation for the one-step-ahead
376
forecasting model are the Mean Squared Error (MSE), Root Mean Squared Error (RMSE), Mean
377 378 379 380
Absolute Error (MAE), Standard Deviation (StD) and Coefficient of Determination ( ). Hgg = ∑ − H ) (
(16)
381 18
ACCEPTED MANUSCRIPT
Hgg = ∑ − H ) (
383
385
R
¡¢O£O£¤ ¡¥¢¦ )z ∑ O OQR(O
= 1 − ∑
386
RI PT
Hgg l = ∑ − H
384
OQR(O
387
¡¢O£O£¤ ¡¥¢¦ )z O
§
(17)
(18)
(19)
SC
382
4.4 Results for Training and Validation Set
389
Figure 6, 7, and 8 present the results for the training, validation set and for both data, resulting in
390
M AN U
388
85% of the all data related to the component ( ) , when applied the hybrid SSA-ANFIS-FCM
391
model to forecast one-step-ahead.
392 393
The quality of the model was evaluated regarding the correlation factor for training, validation
397 398 399 400
TE D
396
Figure 9 presents these results. The hybrid SSA-ANFIS-FCM model learned on the training set and the correlation coefficient for the one-step-ahead forecasting achieved = 0.99472. The
EP
395
and for both data set merged as a single one (training and validation) related to Component ( ) .
correlation achieved on Forecasting model related to the validation set was lower than the correlation for the training set, and its value = 0.86626. The last quality measure was applied to evaluate the result when both data set as presented as one time series to the hybrid model, with 85% length, the value achieved was very good: = 0.97458.
AC C
394
401 402
The same procedure was applied to the component () . It worth remember this component has in
403
its structure only noise and few information about the original wind time series. The Figures 10 to
404
12 present the results for one-step-ahead forecasting.
405 406
The quality indices of the model was evaluated regarding the correlation factor for the training,
407
validation and for the combination of both data set related to Component () . The Figure 13 19
ACCEPTED MANUSCRIPT
408
presents these results. The hybrid SSA-ANFIS-FCM model learned the structure of the training set
409
set = 0.47734. The last performance analysis was done by presenting all data of both set to the
411
model, and the = 0.72985.
RI PT
410
and the correlation coefficient for the forecasted values achieved = 0.97158, for the validation
Figure 6 – Results for Training Set - 70% length component ( )
EP
414 415
AC C
413
TE D
M AN U
SC
412
20
TE D
Figure 7 – Results for Validation Set – 15% length component ( )
AC C
EP
416 417 418 419 420 421
M AN U
SC
RI PT
ACCEPTED MANUSCRIPT
422 423
21
ACCEPTED MANUSCRIPT
430 431 432
Figure 9 – coefficient for training, validation and for both data set merged as single one for
AC C
429
EP
TE D
M AN U
SC
425 426 427 428
Figure 8 – Results for both data set (training and validation merged as single one) for Component ( )
RI PT
424
Component ( )
22
433
Figure 10 – Results for Training Set – 70% length Component ()
AC C
EP
TE D
434
M AN U
SC
RI PT
ACCEPTED MANUSCRIPT
435 436 437 438
Figure 11 – Results for Validation Set – 15% length Component ()
23
M AN U
SC
RI PT
ACCEPTED MANUSCRIPT
439
TE D EP
441 442 443
Figure 12 – Results for both data set (training and validation merged as single one) for Component ()
AC C
440
24
444 445
TE D
M AN U
SC
RI PT
ACCEPTED MANUSCRIPT
Figure 13 – coefficient for training, validation and for both data set merged as single one for
EP
446 447
449
As presented in Figure 2 and in Pseudo-Code 1, after obtain the one-step-ahead forecast for each
component ( ) and () , we need to perform the summation of these two components to
AC C
448
Component ()
450
reconstruct the original time series and the continuation of the series with the forecasted values, it
451
means that
452
() ( )
= ( ) + ( ). The Figure 14 presents the result for the original time ( )
()
series reconstructed after forecasting one-step-ahead process. Therefore, Figure 15 presents the
453
coefficient for all data after performs one-step-ahead forecasting.
454
The results for -steps-ahead forecasting are presented in Figure 16. Where the forecasting was
455 456
performed in recursive way, it means the first value forecasted became an input for the next value in 25
ACCEPTED MANUSCRIPT
460 461 462
463
hours forecast horizon. The 15% remaining of the original time series was applied to test the quality of the forecasted values, and the = 0.97261 was the coefficient resulting when the forecasted time series was compared to the test data set.
RI PT
459
850 points, 3 = 6, sampled by each 10min, were applied to obtain 144 steps resulting in 24
ª(2 + 1) = d( (2), (2 − 1), (2 − 2), (2 − 3), (2 − 4), (2 − 5))
ª(2 + 2) = d( ª(2 + 1), (2), (2 − 1), (2 − 2), (2 − 3), (2 − 4))
ª(2 + 3) = d« ª(2 + 2), ª(2 + 1), (2), (2 − 1), (2 − 2), (2 − 3)¬ ⋮
ª(2 + 6) = d« ª(2 + 5), ª(2 + 4), ª(2 + 3), ª(2 + 2), ª(2 + 1), (2)¬
SC
458
the forecasting process as Equation 20. As input for the forecasting model, a wind time series with
(20)
M AN U
457
AC C
EP
TE D
464
465 466 467 468
Figure 14 – Final Time Series reconstructed after forecasting one-step-ahead each component ( ) 30: () by SSA-ANFIS-FCM model
26
EP
TE D
Figure 15 – for the reconstructed series after forecasting one-step-ahead by SSA-ANFIS-FCM model
AC C
469 470 471 472 473
M AN U
SC
RI PT
ACCEPTED MANUSCRIPT
27
EP
TE D
Figure 16 – Wind Time Series forecasting 144-step-ahead
AC C
474 475 476 477 478
M AN U
SC
RI PT
ACCEPTED MANUSCRIPT
28
Figure 17 – after forecasting 144-step-ahead by SSA-ANFIS-FCM model
TE D
479 480
M AN U
SC
RI PT
ACCEPTED MANUSCRIPT
4.5 Comparisons Results
482
In this section, we present comparisons related to the proposed model with the ANFIS-FCM and
483
GMDH. For keeping the validity of comparisons between hybrid SSA-ANFIS and the ANFISFCM, the same parameters were used as presented in Table 2. In addition, 3 = 6 for the Time
AC C
484
EP
481
485
Series as input for all models also was kept. For GMDH algorithm, the Table 3 presents the
486
parameters applied, which the values was obtained considering the recommendation found in [24-
487
and Figure 19 the results for one-step-ahead forecasting. Figures 20 and 21 the results for -steps-
488 489
25, 27] . Results related to ANFIS-FCM model are depicted in Figure 18, which presents values ahead forecasting by ANFIS-FCM model.
490 491 492
The results for GMDH model are presented in Figures 22 to 25. For the one-step-ahead forecasting
the = 0.93664, decreasing to 0.89661 for -steps-ahead forecasting. In Table 4 the quality 29
ACCEPTED MANUSCRIPT
indices for all approaches, is presented and through it can be confirmed that the Hybrid SSA-ANFIS
494
outperformed the other two approaches.
495
On the results for -steps-ahead forecasting related to ANFIS-FCM and GMDH depicted on Figure
496
RI PT
493
20 and 24, respectively, a lag is observed between the target values and those that were forecasted.
498
We believe the main cause is due to recursive process applied on forecasting model. As presented in
499
Equation 20, after the maximum lag value, all inputs for the model are previous outputs forecasted,
500
causing some instability. For the proposed hybrid model, due to the filtering process into SSA step,
501
the outputs became more stable achieving better results on the forecasting process.
SC
497
M AN U
502
Table 3 – GMDH parameters
503
EP
505
AC C
504
Value 25 5 0 0.7
TE D
Parameter Max. Number of Neurons in a Layer Maximum Number of Layers Selection Pressure () Train Ratio
30
508
Figure 18 – for ANFIS-FCM model one-step-ahead forecasting
EP
507
AC C
506
TE D
M AN U
SC
RI PT
ACCEPTED MANUSCRIPT
31
M AN U
SC
RI PT
ACCEPTED MANUSCRIPT
509 510
Figure 19 – One-step-forecasting applying ANFIS-FCM
AC C
EP
TE D
511
512 513
Figure 20 – 144-step-forecasting applying ANFIS-FCM 32
ACCEPTED MANUSCRIPT
M AN U
SC
RI PT
514
TE D
518
Figure 21 - after forecasting 144-steps-ahead by ANFIS-FCM model
EP
516 517
AC C
515
33
521
Figure 22 – for one-step-ahead forecasting GMDH model
EP
520
AC C
519
TE D
M AN U
SC
RI PT
ACCEPTED MANUSCRIPT
34
M AN U
SC
RI PT
ACCEPTED MANUSCRIPT
522 523
Figure 23 – One-step-forecasting applying GMDH
AC C
EP
TE D
524
525 526
Figure 24 – 144-steps-forecasting applying GMDH 35
ACCEPTED MANUSCRIPT
531 532
Figure 25 - after forecasting 144-steps-ahead by GMDH model
EP
529 530
Table 4 – Quality Indices for one and k-step-ahead forecasting Technique
AC C
528
TE D
M AN U
SC
RI PT
527
SSA-ANFIS-FCM
ANFIS-FCM
Quality Indices MAE MSE RMSE St.D R2 MAE MSE RMSE St.D
Values for One-step-ahead 0.24053 0.20449 0.45221 0.45238 0.96886 0.48044 0.78746 0.88739 0.88770
Values for -steps-ahead 0.20583 0.08776 0.29625 0.29646 0.97261 0.51360 0.46828 0.68431 0.66950 36
ACCEPTED MANUSCRIPT
GMDH
533
0.88125 0.46025 0.38117 0.61739 0.61774 0.93664
0.85195 0.44280 0.31985 0.56555 0.55554 0.89661
RI PT
R2 MAE MSE RMSE St.D R2
The results provided by SSA-ANFIS-FCM are very close to the real values, we can infer through it
535
that the model learned the structure and of the peaks and trend of the time series. Otherwise, the
536
ANFIS-FCM forecasted quality indices presented in Table 4 were worse than GMDH approach, but
537
when the SSA is combined to ANFIS-FCM, the improvement are really promising.
M AN U
SC
534
538
5 Conclusion
540
In this study, a novel hybrid modeling method is proposed for short-term wind speed forecasting by
541
integrating SSA to ANFIS-FCM approach. This hybrid model combines the power of components
542
extraction and reconstruction provided by SSA algorithm. Moreover, the ANFIS characteristics of
543
learning a nonlinear system behavior introduce the precision and stability of the proposed approach
544
in forecasting wind speed time series. Since the dynamic behavior is recursively updated in the
545
proposed SSA–ANFIS approach, when more than one-step-ahead forecast were performed, the
546
stochastic uncertainty and fluctuations of wind speed can be better accounted for by the presented
547
method than by applying only the ANFIS-FCM and GMDH strategies. The better performance of
EP
hybrid SSA-ANFIS-FCM in both one-step-ahead and -steps-ahead wind speed forecasting
AC C
548
TE D
539
549
indicates the robustness of the proposed approach. Future researches will comprise to derive a
550
hybrid model by applying SSA into GMDH algorithm; also, the SSA-ANFIS-FCM model will be
551
compared, aiming to achieve the high accuracy on wind power forecasting.
552
37
ACCEPTED MANUSCRIPT
References
554
[1] O. Kramer, F. Gieseke, B. Satzger. Wind energy prediction and monitoring with neural
555
computation. Neurocomputing, 109 (2013), pp. 84–93
556
RI PT
553
[2] S. Salcedo-Sanz, J. Rojo-Álvarez, M. Martínez-Ramón, G. Camps-Valls. Support vector
558
machines in engineering: an overview. Wiley Interdiscip. Rev. Data Min. Knowl. Discov. 4
559
(2014), pp. 234–267
SC
557
560
562
[3] J. Heinermann, O. Kramer, Machine learning ensembles for wind power prediction, Renewable
M AN U
561
Energy, Volume 89, April 2016, Pages 671-679, ISSN 0960-1481.
563
[4] J. Wang, J. Hu, A robust combination approach for short-term wind speed forecasting and
565
analysis – Combination of the ARIMA (Autoregressive Integrated Moving Average), ELM
566
(Extreme Learning Machine), SVM (Support Vector Machine) and LSSVM (Least Square SVM)
567
forecasts using a GPR (Gaussian Process Regression) model, Energy, Volume 93, Part 1, 15
568
December 2015, Pages 41-56, ISSN 0360-5442.
569
TE D
564
[5] Y. Zhao, L. Ye, Z. Li, X. Song, Y. Lang, J. Su, A novel bidirectional mechanism based on time
571
series model for wind power forecasting, Applied Energy, Volume 177, 1 September 2016,
572
Pages 793-803, ISSN 0306-2619.
AC C
573
EP
570
574
[6] Y. Wang, J. Wang, X. Wei, A hybrid wind speed-forecasting model based on phase space
575
reconstruction theory and Markov model: A case study of wind farms in northwest China,
576
Energy, Volume 91, November 2015, Pages 556-572, ISSN 0360-5442.
577 578
[7] G. Santamaría-Bonfil, A. Reyes-Ballesteros, C. Gershenson, Wind speed forecasting for wind
579
farms: A method based on support vector regression, Renewable Energy, Volume 85, January
580
2016, Pages 790-809, ISSN 0960-1481.
581 38
ACCEPTED MANUSCRIPT
[8] Z. Men, E. Yee, F.-S. Lien, D. Wen, Y. Chen, Short-term wind speed and power forecasting
583
using an ensemble of mixture density neural networks, Renewable Energy, Volume 87, Part 1,
584
March 2016, Pages 203-211, ISSN 0960-1481.
RI PT
582
585
[9] Y-L. Tu, T.-J. Chang, C.-L. Chen, Y.-J. Chang, Estimation of monthly wind power outputs of
587
WECS with limited record period using artificial neural networks, Energy Conversion and
588
Management, Volume 59, July 2012, Pages 114-121, ISSN 0196-8904.
SC
586
589
[10] K. G. Sheela, S.N. Deepa, Neural network based hybrid computing model for wind speed
591
prediction, Neurocomputing, Volume 122, 25 December 2013, Pages 425-429, ISSN 0925-2312.
592
M AN U
590
593
[11] Q. Cao, B. T. Ewing, Mark A. Thompson, Forecasting wind speed with recurrent neural
594
networks, European Journal of Operational Research, Volume 221, Issue 1, 16 August 2012,
595
Pages 148-154, ISSN 0377-2217.
596
[12] D.A. Fadare, The application of artificial neural networks to mapping of wind speed profile for
598
energy application in Nigeria, Applied Energy, Volume 87, Issue 3, March 2010, Pages 934-942,
599
ISSN 0306-2619.
TE D
597
EP
600
[13] J. Zhou, J. Shi, G. Li, Fine tuning support vector machines for short-term wind speed
602
forecasting, Energy Conversion and Management, Volume 52, Issue 4, April 2011, Pages 1990-
603
1998, ISSN 0196-8904.
604 605 606
AC C
601
[14] J. Wang, J. Hu, K. Ma, Y. Zhang, A self-adaptive hybrid approach for wind speed forecasting, Renewable Energy, Volume 78, June 2015, Pages 374-385, ISSN 0960-1481.
607 608 609
[15] O. B. Shukur, M. H. Lee, Daily wind speed forecasting through hybrid KF-ANN model based on ARIMA, Renewable Energy, Volume 76, April 2015, Pages 637-647, ISSN 0960-1481.
610
39
ACCEPTED MANUSCRIPT
611 612
[16] A.D.D. Zhigljavsky (eds) (1997) Principal components of time series: the “Caterpillar” method. St.Petersburg Press, St. Petersburg (in Russian).
613
615
[17] N. Golyandina, V. Nekrutkin, A. Zhigljavsky, Analysis of time series structure: SSA and related techniques, (2001) Chapman & Hall/CRC, 305 pgs.
616
618
[18] J.B. Elsner, A.A.Tsonis, (1996) Singular spectrum analysis: a new tool in time series analysis. Plenum, New York.
SC
617
619
621
[19] J. S. R. Jang, “ANFIS: adaptive-network-based fuzzy inference system,” IEEE Transactions on
M AN U
620
Systems, Man and Cybernetics, vol. 23, pp. 665-685, 1993.
622 623 624
[20] E. Akkaya, ANFIS based prediction model for biomass heating value using proximate analysis components, Fuel, Volume 180, 15 September 2016, Pages 687-693, ISSN 0016-2361.
625
627
[21] A. G. Ivakhnenko, Polynomial Theory of Complex Systems. IEEE Transactions on Systems,
TE D
626
Man, and Cybernetics, Vol. SMC-1, No. 4, pp. 364-378 (1971).
628
630
[22] A.G. Ivakhnenko, A.J. Müller, Self-organization of neuron nets with active neurons. Pattern Recognition Image Analysis 1994;4(2):177e88.
EP
629
631
633 634
[23] A.G. Ivakhnenko, D. Wunsch, Inductive sorting/out GMDH algorithms with polynomial
AC C
632
RI PT
614
complexity for active neurons of neural network. IEEE 1999;6 (99): 1169e73.
635
[24] M.G. De Giorgi, M. Malvoni, P.M. Congedo, Comparison of strategies for multi-step ahead
636
photovoltaic power forecasting models based on hybrid group method of data handling networks
637
and least square support vector machine, Energy, Volume 107, 15 July 2016, Pages 360-373,
638
ISSN 0360-5442.
639 640
[25] M. Dorn, A.L.S. Braga, C.H. Llanos, L.S. Coelho, A GMDH polynomial neural network-based
641
method to predict approximate three-dimensional structures of polypeptides, Expert Systems 40
ACCEPTED MANUSCRIPT
642
with Applications, Volume 39, Issue 15, 1 November 2012, Pages 12268-12279, ISSN 0957-
643
4174.
644
[26] Z. Zhang, Y. Song, F. Liu, J. Liu, Daily Average Wind Power Interval Forecasts Based on an
646
Optimal Adaptive-Network-Based Fuzzy Inference System and Singular Spectrum Analysis.
647
Sustainability 2016, 8, 125.
RI PT
645
648
[27] S.-Kwun Oh, W. Pedrycz, B.-Jun Park, Polynomial neural networks architecture: analysis and
650
design, In Computers & Electrical Engineering, Volume 29, Issue 6, 2003, Pages 703-725, ISSN
651
0045-7906, https://doi.org/10.1016/S0045-7906(02)00045-9.
M AN U
SC
649
AC C
EP
TE D
652
41
ACCEPTED MANUSCRIPT Research highlights •
An hybrid SSA-ANFIS-FCM approach for wind speed forecasting is proposed
•
The hybridization of SSA and ANFIS-FCM is an original idea on context of Machine Learning applied to wind speed forecast The model presents excellent results with stable performance
AC C
EP
TE D
M AN U
SC
RI PT
•