A New Hybrid Forecasting Strategy Applied to Mean Hourly Wind Speed Time Series

Abstract

An alternative electric power source, such as wind power, has to be both reliable and autonomous. An accurate wind speed forecasting method plays the key role in achieving the aforementioned properties and also is a valuable tool in overcoming a variety of economic and technical problems connected to wind power production. The method proposed is based on the reformulation of the problem in the standard state space form and on implementing a bank of Kalman filters (KF), each fitting an ARMA model of different order. The proposed method is to be applied to a greenhouse unit which incorporates an automatized use of renewable energy sources including wind speed power.

1. Introduction

Energy is considered amongst the most significant factors that are closely related to both economic and social developments. It is also a fact that nowadays the majority of the electrical energy production is based on the fossil fuels, which on one hand are, without any doubt, highly efficient but on the other are responsible for the emission of greenhouse gases and their reserves are limited.

Consequently renewable sources of energy, such as wind, biomass, solar power, and wave power, have been already adopted for electric power production. It is well known that the wind power generation raises issues of reliability due to the fact that the wind speed is significantly and directly affected by various factors such as the type of the terrain, the height, season of the year, atmospheric conditions, obstacles present, and many more. This leads to the conclusion that unless the reliability of the wind power generation is at an acceptable level, wind power is not eligible for constant electrical energy supply to the power system [1, 2].

Recent studies have shown that combined forecasting methods can offer robust solutions and can be efficiently implemented to various real-life problems in diverging fields such as chemical processes, economics, load forecasting, tourism demand, environmental issues, medicine, and many more [3–7].

In this study a hybrid model is presented that reveals the advantages of an ARMA and SVM model in wind speed modelling and prediction problem. Initially successful model identification and parameter estimation have to be performed in order to choose the most appropriate ARMA models. For tackling this task the well-established MMPA was used. This approach was introduced by Lainiotis [8, 9] and summarizes the parametric model uncertainty into an unknown, finite dimensional parameter vector whose values are assumed to lie within a known set of finite cardinality. A nonexhaustive list of the reformulation, extension, and application of the MMPA approach as well as its application to a variety of problems can be found in [10–19].

In this research real data were used, provided by Vestas Hellas; the simulation results appear to be very promising.

2. Hybrid Model Presentation

2.1. The ARMA Model

The problem of fitting an ARMA model in a given time series is present for more than half a century and is still appearing in many different fields such as in remote monitoring of civil infrastructure [20], predicting the demand for auto spare parts in China due to the fierce market competition [21], forecasting equipment failures in order to adjust maintenance policies in manufacturing plants [22], estimating retail sales volumes [23], predicting the outbreak and development of animal infectious diseases [24], and many more.

Considering the general case an 𝑚-variate (i.e., multivariate) ARMA model of order (𝑝,𝑞) [ARMA(𝑝,𝑞)] for a stationary time series of vectors 𝑦 observed at equally spaced instants 𝑘=1,2,…,𝑛 is defined as 𝐲_𝑘=^𝑝_𝑖=1𝚨_𝑖𝐲_𝑘−𝑖+^𝑞_𝑗=1𝐁_𝑗𝐯_𝑘−𝑗+𝐯_𝑘,,𝐸𝐯_𝑘𝐯^𝑇_𝑘=𝐑, (1) where the 𝑚-dimensional vector 𝐯_𝑘 is uncorrelated random noise, not necessarily Gaussian, with zero mean; covariance matrix 𝐑, 𝜃=(𝑝,𝑞) is the order of the predictor; and 𝐀₁,…,𝐀_𝑝, 𝐁₁,…,𝐁_𝑞 are the 𝑚×𝑚 coefficient matrices of the multivariate (MV) ARMA model.

It is obvious that the problem requires both the predictor’s order 𝜃=(𝑝,𝑞) determination and computation of the predictor’s matrix coefficients {𝐀_𝑖,𝐁_𝑗}.

The major disadvantage of the ARMA models is that their performance can be limited by any significant data nonlinearities.

2.2. Multimodel Partition Algorithm (MMPA)

Due to the fact that the wind speed does not have a constant or periodic behaviour, it was noted, by trial and error, that not a single ARMA model that was able to describe the whole data set satisfactory. It is actually the combination of various ARMA models, each one used for different time intervals and applied for different time durations that describes in the best manner the existing data. So instead of having various ARMA models of different order 𝜃 running in parallel with the SVM it was decided to load all the data to an adaptive filter programmed with the MMPA, and it will be the job of that filter to decide which ARMA model will be used each time.

If we assume that the model order fitting the data is known and is equal to 𝜃=(𝑝,𝑞), we can rewrite (1) in standard state-space form as 𝐱()(𝑘),𝐲(𝑘)(𝑘)𝐱(𝑘)(𝑘).𝑘+1=𝐱=𝐇+𝐯 (2)

Now assign a new variable 𝜆 such as 𝜆=max(𝑝,𝑞). Then 𝐱(𝑘) is an 𝑚²  (𝜆+𝜆)×1 vector made up from the coefficients of the matrices {𝐀₁,…,𝚨_𝜆,𝐁₁,…,𝐁_𝜆}; 𝐇(𝑘) is an 𝑚×𝑚²(𝜆+𝜆) observation history matrix of the process {𝐲(𝑘)} up to time 𝑘−(𝜆+𝜆).

If the general form of the matrices 𝐀_𝜆 and 𝐁_𝜆 is, respectively, 𝑎^𝜆₁₁…𝑎^𝜆_1𝑚𝑎⋮⋱⋮_𝜆𝑚1⋯𝑎^𝜆_𝑚𝑚,𝑏^𝜆₁₁…𝑏^𝜆_1𝑚𝑏⋮⋱⋮_𝜆𝑚1⋯𝑏^𝜆_𝑚𝑚, (3) then 𝐱(𝑘)≜𝛼¹₁₁𝛼¹₂₁⋯𝛼¹_𝑚1⋮𝛼¹₁₂𝛼¹₂₂⋯𝛼¹_𝑚2⋮⋯𝛼¹_𝑚𝑚⋮⋯𝛼^𝜆_𝑚𝑚⋮𝑏¹₁₁𝑏¹₂₁⋯𝑏¹_𝑚1⋮𝑏¹₁₂𝑏¹₂₂⋯𝑏¹_𝑚2⋮⋯𝑏¹_𝑚𝑚⋮⋯𝑏^𝜆_𝑚𝑚^𝑇. (4) If 𝑝>𝑞 then 𝜆=𝑝; the last (𝑚²(𝑝−𝑞)) MA coefficients are zero. If 𝑞>𝑝 then 𝜆=𝑞; the last (𝑚²(𝑞−𝑝)) AR coefficients are zero. Consider 𝐇(𝑘)≜𝑦₁(𝑘−1)𝐈⋯𝑦_𝑚(𝑘−1)𝐈⋮⋯⋮𝑦₁(𝑘−𝜆)𝐈⋯𝑦_𝑚(𝑘−𝜆)𝐈⋮𝑣₁(𝑘−1)𝐈⋯𝑣_𝑚(𝑘−1)𝐈⋮⋯⋮𝑣₁(𝑘−𝜆)𝐈⋯𝑣_𝑚(𝑘−𝜆)𝐈, (5) where 𝐈 is the 𝑚×𝑚 identity matrix.

Assuming that the system model and its statistics were completely known, the Kalman filter (KF) in its various forms would be the optimal estimator in the minimum variance sense.

However, if the system model is not completely known the MMPA, introduced by Lainiotis [8, 9], is one of the most widely used approaches for similar problems [10–19, 25–28].

In the case under consideration assume that the model uncertainty is the lack of knowledge of the model order 𝜃. Let us further assume that the model order 𝜃 lies within a known sample space of finite cardinality; that is, 1≤𝜃≤𝑀, 𝜃∈ℑ, where ℑ denotes the set of integers. The MMPA operates on the following discrete-time model: 𝐱()𝑘𝑘+1=𝐅𝑘+1,𝜃𝐱(𝑘)(𝑘),𝐲(𝑘)𝑘+𝐰=𝐇𝜃𝐱(𝑘)+𝐯(𝑘), (6) where 𝜃 is the unknown parameter, the model order in this case; 𝐅 is the state transition matrix; and 𝐰(𝑘) is independent, zero mean, white noise not necessarily Gaussian with covariance 𝐐 which is usually set to a small positive nonzero constant. The optimal MMSE (minimum mean square error) estimate of 𝐱(𝑘) is given by 𝐱𝑘𝑘=^𝑀_𝑗=1𝐱𝑘𝑘;𝜃_𝑗𝑝𝜃_𝑗𝑘. (7) A set of 𝑀 models is designed, each matching one value of the parameter vector, {(1,1),(2,2),…(𝑀,𝑀)}. The probabilities 𝑝(𝜃_𝑗/𝑘) for each model are set to 1/𝑀, where 𝑀 is the cardinality of the model set.

A bank of Kalman filters is then applied, one for each model, which can be run in parallel, thus saving enormous computational time. At each iteration, the MMPA selects the model that corresponds to the maximum a posteriori probability as the correct one. This probability tends (asymptotically) to one, while the remaining probabilities tend to zero. The overall optimal estimate can be taken either to be the individual estimate of the elemental filter exhibiting the maximum posterior probability, for example, a value of 0.9 or higher [14], or the weighted average of the estimates produced by the elemental ARMA filters, as described in (7), which is the case used in this paper.

The probabilities are calculated on-line in a recursive manner as it is shown by 𝑝𝜃_𝑗𝑘=𝐿𝑘/𝑘;𝜃_𝑗∑^𝑀_𝑗=1𝐿𝑘/𝑘;𝜃_𝑗𝑝𝜃_𝑗𝑝𝜃/𝑘−1𝑘^𝑗,𝐿𝑘−1𝑘;𝜃_𝑗=𝐏_𝐲𝑘𝑘−1;𝜃_𝑗^−1/2−1⋅exp2𝐲^𝑇𝑘𝑘−1;𝜃_𝑗×𝐏_𝐲⁻¹𝑘𝑘−1;𝜃_𝑗𝐲𝑘𝑘−1;𝜃_𝑗, (8) where the innovation process 𝐲𝑘𝑘−1;𝜃_𝑗=𝐲(𝑘)−𝐇𝑘;𝜃_𝑗𝐱𝑘𝑘−1;𝜃_𝑗 (9) is a zero mean white process with covariance matrix 𝐏_𝐲𝑘𝑘−1;𝜃_𝑗=𝐇𝑘;𝜃_𝑗𝐏𝑘𝑘;𝜃_𝑗𝐇^𝐓𝑘;𝜃_𝑗+𝐑. (10)

For equations (8)–(10) 𝑗=1,2,…,𝑀.

2.3. Support Vector Machines (SVM)

In support vector machines, as they were proposed in [29], the training data set 𝑥_𝑖∈𝑅^𝑑, 𝑖=1,…,𝑁, is mapped into a higher dimensional feature space via an operator Φ.

A mathematical representation of the SVM function is 𝑦=𝜔⋅Φ(𝑥)+𝑏, (11) where 𝜔 and 𝑏 can be found by the minimization of the following equations: 𝑆(𝐶)1=𝐶𝑀^𝑀_𝑗=1𝐿𝑓_𝑒ℎ_𝑗−𝑦_𝑗+12‖𝜔‖²,𝐿𝑓ℎ,𝑦=ℎ−𝑦−𝑒ℎ−𝑦≥𝑒,0others, (12)(13) where parameters 𝐶 and 𝑒 are user-defined. The term ℎ_𝑖 is the actual wind speed at the time instant 𝑗 and term 𝐿𝑓_𝑒(ℎ_𝑖−𝑦_𝑖) is the loss function. By looking at (13) it is obvious that there is any penalty for errors below 𝑒. The width of the function is given by the term (1/2)‖𝜔‖² and finally the training error term is given by 𝐶(1/𝑀)∑^𝑀_𝑗=1𝐿𝑓_𝑒(ℎ_𝑖−𝑦_𝑖), where 𝐶 is the tradeoff between the width of the function and the minimum training error. For dealing with nonlinear cases, like wind speed data, one may introduce slack variables 𝜉 and 𝜉^∗ into (11) such that 𝜔⋅Φ𝑗+𝑏_𝑗−ℎ_𝑗≤𝑒+𝜉^∗_𝑗,−𝜔⋅Φ𝑗+𝑏_𝑗+ℎ_𝑗≤𝑒+𝜉𝑗, (14) where 𝜉_𝑗,𝜉^∗_𝑗≥0, and 𝑗=1,2,…𝑀.

By considering the above slack variables and in order to include any extra cost of the training errors, (12), which represents the objective function to be minimized, is rearranged to 𝑆𝜔,𝜉,𝜉^∗=𝐶^∗𝑀_𝑗=1𝜉_𝑗+𝜉^∗_𝑗+12𝜔𝜔^𝑇, (15) where again 𝐶^∗ is user-defined and is the tradeoff between the maximum margin defined by ‖𝜔‖ and the minimum training error as defined by ∑^𝑀_𝑗=1(𝜉_𝑗+𝜉^∗_𝑗).

Finally by introducing positive Lagrangian multipliers and maximizing (15) the latter equation is reformed to 𝑆𝑎_𝑗−𝑎^∗_𝑗=^𝑀_𝑗=1ℎ_𝑖𝑎_𝑗−𝑎^∗_𝑗−𝑒^𝑀_𝑗=1𝑎_𝑗−𝑎^∗_𝑗−12^𝑀𝑀_𝑗=1𝑖=1𝑎_𝑗−𝑎^∗_𝑗×𝑎_𝑖−𝑎^∗_𝑖𝐾𝑥_𝑗−𝑥_𝑖 (16) subject to ^𝑁_𝑗=1𝑎_𝑗−𝑎^∗_𝑗=0,0≤𝑎_𝑗,𝑎^∗_𝑗≤𝐶, (17) where 𝑗=1,2,…𝑀.

The Lagrangian multipliers 𝑎_𝑗, 𝑎^∗_𝑗 satisfy 𝑎_𝑗∗𝑎^∗_𝑗=0 and also 𝑓𝑥,𝑎,𝑎^∗=^𝑙_𝑗=1𝑎_𝑗−𝑎^∗_𝑗𝐾𝑥−𝑥_𝑖+𝑏. (18) The Kernel function 𝐾(𝑥−𝑥_𝑖) introduced in (18) is defined such that 𝐾(𝑥_𝑗−𝑥_𝑖)=Φ(𝑥_𝑖)⋅Φ(𝑥_𝑗), meaning that its value is equal to the inner product of the vectors 𝑥_𝑖 and 𝑥_𝑗, included in the featured space Φ(𝑥_𝑖) and Φ(𝑥_𝑗).

In this study the Gaussian kernel function (19), also known as radial basis function (RBF), is used. Consider 𝐾𝑥_𝑗−𝑥_𝑖=𝑒^{(−‖𝑥𝑗−𝑥𝑖‖2/2𝜎2)}. (19) The most significant feature of the SVM compared to other similar algorithms is that they manage to achieve optimum performance by restricting the decision’s function complexity so that is the most suitable according to the quantity of the data present.

2.4. The Hybrid Model

The wind speed behaviour is unpredictable and it is difficult to be represented. This is the reason for combining two different techniques for modelling the linear and the nonlinear parts of the series. The hybrid model proposed is based on a linear pattern, 𝐿(𝑡) produced by the MMPA, and a nonlinear one, 𝑁𝐿(𝑡) produced by the SVM. It can be represented as 𝑄(𝑡)=𝐿(𝑡)+𝑁𝐿(𝑡). (20) Both parts are directly calculated from the wind speed time series.

If 𝑒(𝑡) is the MMPA estimation error at any time instant 𝑡, then 𝑒(𝑡)=𝑍(𝑡)−𝐿(𝑡). (21) It is now the SVM that models these residuals as 𝑒(𝑡)=𝑓𝑒_(𝑡−1),𝑒_(𝑡−2),…,𝑒_{(𝑡−𝑛)}+Δ𝑡, (22) where 𝑓 is nonlinear and Δ𝑡 is random error. It is obvious that 𝑁𝐿(𝑡) is the forecast of (22).

Consequently the forecast of the hybrid model is 𝑄(𝑡)=𝐿(𝑡)+𝑁𝐿(𝑡). (23) The proposed method is schematically represented in Figure 1.

Figure 1 

Proposed method.

3. Results and Discussion

3.1. Results

In this method the weighted average of the estimates produced by the elemental ARMA filters was used as a data preprocessor in order to detect the data’s linearities. This was succeeded using a bank of 10 Kalman filters of order (1,1),(2,2),(3,3),…,(10,10) programmed with the MMPA. Then the MMPA’s estimation error was applied as input to the SVM that were able to achieve a further error reduction and come up with a better forecasting outcome. As far as the SVM are concerned the three parameters (𝐶, 𝜎, and 𝑒) had to be carefully adjusted. Unsuitable values for these parameters may lead to either overfitting or underfitting of the training data. The values used in this research were 𝐶=35, 𝜎=2.5, and 𝑒=0.6.

This research was conducted based on the hour average of daily wind speed recorded by the Vestas Hellas from November 2010 up to February 2011. The obtained time series did not follow any periodic pattern and it was also presenting irregular amplitudes, making it hard to both model and predict (Figures 2, 3, 4, and 5; “raw” data).

Figure 2 

November 2010: predicted and observed (raw) time series.

Figure 3 

December 2010: predicted and observed (raw) time series.

Figure 4 

January 2011: predicted and observed (raw) time series.

Figure 5 

February 2011: forecasted predicted and observed (raw) time series.

The aim of this work is to generate a single-step prediction based on past observations. The data were normalized to take values from zero to one, before using them as input data to the hybrid model.

From the 2725 available data points, 720 were for November, 744 for December, 744 for January, and 517 for February. For each month 20% of the available data was used for training, 20% for validation, and 60% for testing.

The performance of the hybrid method is judged by (a) comparing the predicted and the observed (raw) wind time series (Figures 2, 3, 4, and 5), (b) drawing scatter diagrams of the predicted and the observed sequences (Figures 6, 7, 8, and 9), and (c) computing the mean percentage absolute error (MAPE) for the testing data set, using the mathematical formula given in (22). Table 1 summarizes the results. Consider 1MAPE=𝑀^𝑀_𝑗=1P_real,𝑗−P_{predicted,𝑗}P_real,𝑗×100. (24)

Table 1 

Performance of the hybrid model.

Figure 6 

November 2010: scatter of predicted and observed wind speed time series.

Figure 7 

December 2010: scatter of predicted and observed wind speed time series.

Figure 8 

January 2011: scatter of predicted and observed wind speed time series.

Figure 9 

February 2011: scatter of predicted and observed wind speed time series.

3.2. Discussion

Figures 2–5 indicate that the predictions are very close to the real values of the wind speed time series. Additionally the hybrid model is capable of following satisfactorily the irregularities of the observed series. The value of MAPE between the real and the predicted series can be considered quite low since it has an average value below 3.01% and the individual value for each month is not above 3.3%. Furthermore the value of 𝑅² is another indication of close is the predicted to real series. A value of 𝑅²=1 indicates a straight line of the form 𝑦=𝑎𝑥 meaning that the predicted values are exactly the same as the real ones. In this study 𝑅² is at an average of 0.8685 which means that the predicted series is indeed very close to the real one.

For the sake of completeness of this research the results of forecasting using only the adaptive combination of MMPA with the ARMA models as well as using only the ANN architecture with the SVM are also presented. The first method was applied successfully in [22] for electric load modeling and forecasting. The difference is that electric load data is seasonal and after careful manipulation it can be converted to Gaussian, something that cannot be done with a wind speed time series.

Using the data for November, Figures 10 and 11 compare the predicted and the observed (raw) wind time series and Figures 12 and 13 show the scatter diagram of the predicted and observed sequences for MMPA-ARMA and SVM, respectively. Table 2 summarizes the performance of the aforementioned methods and the hybrid method proposed for the whole data set.

Table 2 

Summary of results.

Figure 10 

November 2010: predicted and observed (raw) time series, MMPA-ARMA.

Figure 11 

November 2010: predicted and observed (raw) time series, SVM.

Figure 12 

November 2010: scatter of predicted and observed wind speed time series, MMPA-ARMA.

Figure 13 

November 2010: scatter of predicted and observed wind speed time series, SVM.

At this point it should be mentioned that when dealing with the wind speed forecasting problem it is quite difficult to attain very high prediction accuracy. The hybrid method proposed gave some very satisfactory results and the accuracy level reached can be considered sufficient for decision making as far as electric power production is concerned.

Last but not least it should be also mentioned that a designer’s crucial task when using MMPA is to assign a proper value for the cardinality 𝑀. Trying to find an answer to that problem one has to take in mind the following considerations:(i)the quality of the overall MMPA estimate (in terms of MAPE) increases with the number of the elemental filters applied, so a small value of 𝑀 leads to estimates of a poor quality;(ii)the computational load of the MMPA is proportional to the number of filters implemented; so a large value 𝑀 increases the computational burden, a factor that makes the implementation of the MMPA for real time problems difficult.

A way of assigning the proper value of 𝑀 is to record the behaviour of MMPA during its operation. It is already mentioned that, at each iteration, the MMPA selects the model that corresponds to the maximum a posteriori probability as the correct one. This probability tends (asymptotically) to one, while the remaining probabilities tend to zero. Figure 14 shows the behaviour of the MMPA for the whole data set concerning December.

(a)

(b)

Figure 14 

December 2010: MMPA probability sequence (ARMA model order selection) for the whole data set provided.

It can be seen that the highest order 𝜃 of the appropriate ARMA model for the time series under investigation is (4,4) which makes the choice of 𝑀=10 acceptable. Another also significant conclusion is that if a model change is required, the algorithm senses the variation and increases its corresponding a posteriori probability while decreasing the remaining ones. Thus, the algorithm is adaptive in the sense of being able to track model changes in real time. This procedure incorporates the algorithm’s intelligence.

4. Conclusions

The area of forecasting is very demanding and it is over 50 years that ARMA models were exclusively used for tackling real-life problems. Recently ANN were applied in difficult prediction problems showing very satisfactory results especially due to their ability of manipulating the nonlinearities of the dataset. The aim of this work was not to just add yet another technique of wind speed prediction but to actually validate the fact that different forecasting methods fulfil each other and lead to accurate results. As it was shown the two individual forecasting methods, adaptive MMPA-ARMA and SVM, cannot match the performance of the hybrid method proposed. The results of the proposed method are to be applied to a greenhouse unit which incorporates an automatized use of renewable energy sources including wind speed power. Future work can include adjustments for wind speed prediction for time intervals smaller than 1 hour, say ten-minute intervals, and also for on-line wind speed prediction.

Conflict of Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.

Acknowledgments

This work was supported in part by the European Union (European Social Fund-ESF) and Greek national funds through the Operational Program “Education and Lifelong Learning” of the National Strategic Reference Framework (NSRF), Research Funding Program: ARCHIMEDES III. Investing in knowledge society through the European Social Fund.