A combined model for short-term wind speed forecasting based on empirical mode decomposition, feature selection, support vector regression and cross-validated lasso

The whole process of the proposed model

The architecture of our proposed model is shown in Fig. 1. The whole process is as follows:

1. Use EMD to decompose wind speed into a series of IMFs. EMD algorithm is introduced in ‘Empirical model decomposition’

2. Combine FS and SVR to predict the high-frequency IMF obtained by EMD. FS and SVR algorithms are provided in ‘Feature selection’ and ‘Support vector regression’, respectively.

3. Use LassoCV to complete the prediction of the low-frequency IMF and trend. LassoCV algorithm is listed in ‘Cross-validated lasso’.

4. Performance evaluation. The performance indicators are introduced in ‘Prediction performance criteria’, and the experimental results and analysis are given in ‘Results’ and ‘Discussion’.

Empirical model decomposition

Due to the non-stationarity, intermittent and inherent nature of wind speed, it is difficult to directly predict the future wind speed. One possible solution is to decompose different frequencies from chaotic wind data (Bokde et al., 2019) and use models to predict them separately. Based on this idea, the study introduces signal processing technology to decompose wind speed. Common signal decomposition algorithms include Wavelet transform, morphology filters, EMD and many others. Wavelet transform is not adaptive and follows the prior knowledge of its mother wavelet, so somewhat limits its ability to extract nonlinear and non-stationary components from the data. Similarly, the morphology filters have to select the shape and the length of the structural element. There is no uniform standard and depends on human experience, whereas EMD has received great attention from researchers because of its superior performance and easy-to-understand. Therefore, in this study, we used EMD for preprocessing the wind speed.

EMD is essentially a non-linear signal analysis method that can handle non-linear and non-stationary time series (Huang et al., 1998). EMD uses the time-scale characteristics of the data to decompose the signal, and does not need to set any basis functions in advance. In theory, EMD can be applied to any type of signal. Since EMD was proposed, it has been rapidly applied to many different engineering fields such as marine and atmospheric research, seismic record analysis and mechanical fault diagnosis (Gao & Liu, 2021).

The basic idea of EMD is to decompose non-stationary time series signals into a series of IMFs along with a residue (Huang et al., 1998). The IMF should meet two principles: (1) the number of extreme and zero values must be equal or differ by at most one; (2) the average value of upper envelop and lower envelope must be zero (Ziqiang & Puthusserypady, 2007). Let $s\left(t\right),$t =1 , 2, …, l be a time series. EMD decomposition steps are as follows:

Step 1: Identify the local minima and maxima of the time series.

Step 2: Use cubic splines to interpolate local minima and maxima values to generate lower $s_l\left(t\right)$ and upper $s_u\left(t\right)$.

Step 3: Computer the average envelope of the upper and lower envelopes ${\displaystyle m_t=\frac{s_u\left(t\right)+s_l\left(t\right)}{2}}$

Step 4: Subtract the average envelope from the original time series $h\left(t\right)=s\left(t\right)-m_t$

Step 5: Check $h\left(t\right)$ if meets the two principles of IMF. If so, treat $h\left(t\right)$ as the new IMF $c\left(t\right)$ and calculate the residual signal $r\left(t\right)=s\left(t\right)-h\left(t\right)$. Otherwise, replace $h\left(t\right)$ with $s\left(t\right)$, and then repeat steps 1 to 5.

Step 6: Set $r\left(t\right)$ as new $s\left(t\right)$ and repeat steps 1 to 5 until all IMFs are obtained.

Through the whole process, a set of IMFs from high to low frequency can be extracted from the time series. Therefore, the original time series can be expressed as: ${\displaystyle s\left(t\right)=\sum _{i=1}^n{c}_i\left(t\right)+r_n\left(t\right)}$ where n is the number of IMFs. $c_i\left(t\right)$ refers to the IMF, which is periodic and almost orthogonal to each other (Li et al., 2018b). $r_n\left(t\right)$ is the final residual representing the trend of $s\left(t\right)$.

Feature selection

After obtaining the IMF components of wind speed, we need to predict it. In the study, we use the observed and lag of the IMF components as the raw features, respectively forecast each IMF component, and add all the predicted IMF components to get the final wind speed. Despite, the raw features contain sufficient information for forecasting, some irrelevant or partially relevant features in the raw features may have a negative impact on the model. To avoid the impact, a common strategy is to use feature selection to remove irrelevant features. Commonly used feature selection algorithms include filter method, wrapper method, heuristic search algorithm, embedded method (Chandrashekar & Sahin, 2014). In this study, we use the filter method. In order to obtain scores of different variables, we use the univariate linear regression test to calculate the correlation between features and output (Liu et al., 2019b), which is defined as: ${\displaystyle Co{r}_i=\frac{\left(X\left[:,i\right]-mean\left(X\left[:,i\right]\right)\right)\mathrm{\ast }\left(y-mean\left(y\right)\right)}{std\left(X\left[:,i\right]\right)\mathrm{\ast }std\left(y\right)}}$ where X is an N × M matrix, each column is a feature. y is the N × 1 vector of the output we are interested in. Based on the rank of correlation, the irrelevant or partially relevant features are removed.

Support vector regression

The support vector machine (SVM) is a learning method based on structural risk minimization criteria, which can minimize the expected risk and obtain better generalization performance on unknown data. The support vector regression (SVR) is an extension of SVM for regression problems (Drucker et al., 1997). Due to the nonlinear and non-stationary nature of wind speed, SVR is widely used in short-term wind speed forecasting (Khosravi et al., 2018; Liu et al., 2019a; Santamaría-Bonfil, Reyes-Ballesteros & Gershenson, 2016). In the research, we use EMD to decompose the IMF components of wind speed, and the high-frequency IMF component contains the nonlinear and non-stationary part of wind speed. In order to obtain better generalization performance, we refer to existing research and use SVR to predict it.

The main idea of SVR is to implement linear regression in the high-dimensional feature space obtained by mapping the original input through a predefined function $\varnothing \left(x\right)$, and to minimize structure risks (Chen et al., 2018). Given a set of samples $\left\{x_i,y_i\right\},$i =1 , 2, …, N, y_i is the output and x_i is the input. The objective is: ${\displaystyle \begin{array}{c}\hfill f\left(x\right)=W^T\varnothing \left(x\right)+b\hfill \\ \hfill R\left[f\right]=\frac{1}{2}{∥W∥}^2+C\sum _{i=1}^{N}L\left(x_i,y_i,f\left(x_i\right)\right)\hfill \end{array}}$ where W and b are the regression coefficient and bias, respectively. C is the penalty coefficient. $L\left(x_i,y_i,f\left(x_i\right)\right)$ represents the loss function, and $R\left[f\right]$ is the structure risk. The corresponding constrained optimization problem can be expressed as: ${\displaystyle \begin{array}{c}\hfill \mathit{min}\frac{1}{2}{∥W∥}^2+C\sum _{i=1}^N\left({\xi }_i+{\xi }_i^{\mathrm{\ast }}\right)\hfill \\ \hfill \begin{array}{cc}\hfill s.t.\hfill & \hfill y_i-W^T\varphi \left(x\right)-b\le \varepsilon +{\xi }_i\hfill \\ \hfill \hfill & \hfill W^T\varphi \left(x\right)+b-y_i\le \varepsilon +{\xi }_i^{\mathrm{\ast }}\hfill \\ \hfill \hfill & \hfill {\xi }_i,{\xi }_i^{\mathrm{\ast }}\ge 0,i=1,2,\dots ,n\hfill \end{array}\hfill \end{array}}$ where ξ_i and ${\xi }_i^{\ast }$ refer to the slack variables. By introducing the Lagrange multiplier, the regression can be expressed as: ${\displaystyle f\left(x\right)=\sum _{i=1}^N\left({\alpha }_i-{\alpha }_i^{\mathrm{\ast }}\right)K\left(x_i,x\right)+b}$ where α_i and ${\alpha }_i^{\ast }$ are the Lagrange multipliers that satisfy the conditions ${\alpha }_i\ge 0,{\alpha }_i^{\ast }\ge 0$ and ${\sum }_{i=1}^N\left({\alpha }_i-{\alpha }_i^{\ast }\right)=0.K\left(x_i,x\right)$ is the kernel function conforming to Mercer’s theorem.

Cross-validated lasso

The Lasso algorithm is a regression model that can perform feature selection and regularization at the same time. It was originally proposed by Robert Tibshirani of Stanford University, with better prediction accuracy and interpretability (Tibshirani, 1996). Normally, in regression, we want to find a coefficient $\beta =\left({\beta }_1,\dots ,{\beta }_p\right)$ that satisfies the following: ${\displaystyle Y=X\beta +\varepsilon ,E\left[\varepsilon |X\right]=0}$ where Y is the dependent variable, $X=\left(X_1,\dots ,X_N\right)$ is the covariate, and ɛ is the unobserved noise. Lasso tries to minimize the objective function while forcing the sum of the absolute values of the coefficients to be less than a fixed value t (Hung, Yen & Li, 2016): ${\displaystyle {min}_{{\beta }_0,\beta }\left\{\frac{1}{N}\sum _{i=1}^N{\left(y_i-{\beta }_0-x_i^T\beta \right)}^2\right\}}$ ${\displaystyle s.t.\sum _{j=1}^p\left|{\beta }_j\right|\le t.}$

Rewritten in the Lagrangian form: ${\displaystyle {\hat{\beta }}_{lasso}=\underset{\beta \in R^p}{\mathit{argmin}}\left\{\frac{1}{N}{∥y-X\beta ∥}_2^2+\lambda {∥\beta ∥}_1\right\}}$

The L₁-norm is used instead of the L₂-norm in Lasso. Since the constraint region is diamond-shaped, it is more likely to pick the solution that lies at the corner of the region. As a result, the solution of the lasso is sparse, with some coefficients set to exactly equal to zero, that is, Lasso performs a straightforward feature selection.

To estimate ${\hat{\beta }}_{lasso}$, the value of the penalty parameter λ is critically important. However, the optimal λ is not given automatically. If λ is chosen appropriately, Lasso achieves the fast convergence under fairly general conditions; On the other hand (chosen inappropriately), Lasso may be inconsistent or have a slower convergence. In the paper, we adopt the cross-validated Lasso algorithm, in which the penalty parameter λ is chosen based on cross-validation, and this is also the leading recommendation way in the theoretical literature (Park & Casella, 2008).

Prediction performance criteria

In the study the mean absolute percentage error (MAPE) , mean absolute error (MAE) and RMSE are used as performance indicators to evaluate the proposed wind forecasting model, which are defined as follows:

${\displaystyle MAPE=\frac{1}{N}\sum _{i=1}^N\left|\left(Y_i-{\hat{Y}}_i\right)/Y_i\right|}$ ${\displaystyle MAE=\frac{1}{N}\sum _{i=1}^N\left|Y_i-{\hat{Y}}_i\right|}$ ${\displaystyle RMSE=\sqrt{\frac{1}{N-1}\sum _{i=1}^N{\left(Y_i-{\hat{Y}}_i\right)}^2}}$ where Y_i and ${\hat{Y}}_i$ refer to the observed and predicted wind speed of data point i, respectively. For MAPE, MAE, RMSE, the smaller value, the better the performance.

Wind speed data

The wind speed data used in the study is gathered from two wind stations in Michigan, USA from September 2019 to October 2019. The number of data is 1,464. The initial 50 days from September 1, 2019 to October 20, 2019 are employed as input for model training, and the remaining days, i.e., from October 21, 2019 to October 31, 2019 are used to test. Figure 2 shows these two wind speed time series, and the corresponding statistics are listed in Table 1.

Experiments and result analysis

To verify the effectiveness of the proposed model, we compare it with five classic individual models, including Persistence, ELM, SVR and ANN, ARIMA. The 1- to 3-step forecasting results of these models under time series #1 and #2 are displayed in Figs. 3–4, and the corresponding error estimated results are listed in Tables 2–5. It is worth noting that for a fair comparison, the parameters of the involved models are selected based on cross-validation. Based on the experimental results, we can get the following conclusions:

1. In the 1-step forecasting, for wind station #1, the proposed model obtains the best accuracy: RMSE, MAE, and MAPE are 0.5859, 0.4426, and 21.11%, respectively. The classic individual models from low to high based on RMSE are ELM, ANN, Persistence, SVR, and ARIMA, with MAPE values of 36.20%, 36.24%, 36.20%, 34.87%, and 34.25%, respectively. Likely, in wind station #2, compared with the classic individual models, the proposed model still obtains the best performance, and the MAPE value is 17.10%.

2. In the 2-step forecasting, when wind station #1 is used, the proposed model has the lowest performance criteria, i.e., the values of RMSE, MAE, and MAPE are 0.7531, 0.5848, and 24.78%, respectively. In addition, for wind station #2, the proposed model still achieves the lowest performance criteria value. Take MAPE as an example, the value of MAPE is 22.99%, which is significantly lower than other models.

3. In the 3-step forecasting, the proposed model is still the model with the highest prediction accuracy, and the MAPE of wind stations #1 and #2 are 27.55% and 24.59%, respectively. Persistence has the worst RMSE value among these models, with MAPE of 57.64% and 47.99%, respectively.

In general, under 1- to 3-step forecasting, the proposed model can obtain the best prediction performance compared with the classic individual models.

Compared with traditional EMD methods

As a nonlinear signal analysis method for processing nonlinear and non-stationary time series, EMD has been widely used in time series. To further verify the effectiveness of our EMD model, we compare it with four widely used EMD models, namely EMD-ELM, EMD-SVR, EMD-SP-SVR, and EMD-ANN. It is worth noting that in this study, these methods used the same way as our proposed model, using EMD to decompose the wind speed, using a single classifier to predict each IMF component separately, and adding all the prediction results to get the final prediction wind speed. The prediction results and the error estimated results of these four EMD-based methods and the proposed method are displayed in Figs. 5–6 and Tables 6–9. Based on Figs. 5–6 and Tables 6–9, it can be observed that:

1. Compared with the above-mentioned classic individual models, the performance of the EMD-based method is significantly improved. Take wind station #1 as an example, in the 1-step forecasting, the value of RMSE of the EMD-based methods is around 0.60, while the classic individual model is around 1.20. After the wind speed is decomposed by EMD, the value of RMSE is reduced almost doubled.

2. For wind station #1, except for the MAE in the 3-step forecasting, the performance indicators obtained from the proposed model are significantly better than those EMD-based combined models. For the 3-step forecasting, the performance of EMD-SVR and EMD-SVR-SP in MAE is slightly better than the proposed combined model, but in other evaluation indicators, the proposed combined model achieves a significantly better performance. Furthermore, EMD-ANN is always worse in MAPE as compared with the other three combined models, with MAPE of 23.55%, 27.67%, and 29.31% for 1- to 3-step forecasting.

3. For wind station #2, in 1- to 3-step wind speed forecasting, the proposed combined model obtains the best prediction results. The RMSE, MAE and MAPE in the 1-step forecasting are 0.5593, 0.419, and 17.10%, respectively. In comparison, among the other four EMD-based combined models, the EMD-ELM and EMD-ANN models have similar prediction performance in 1- to 3-step forecasting, with MAPE values of 21.59%, 27.49%, 27.65% and 21.83%, 25.3%, 27.86%, respectively.

In total, the EMD-based method has obvious advantages over traditional methods, and the proposed method that using EMD, FS, SVR and LassoCV can achieve better performance.

Performance of SVR-SP and LassoCV on different IMFs

According to the EMD principle, the frequency of the IMF components is from high to low. The non-linear and non-stationary information of wind speed data is mainly concentrated in the high-frequency IMF, and the low-frequency IMF presents a Sin-like function curve. Based on its characteristics, in this study we use SVR-SP and LassoCV to predict IMFs of different frequencies. In order to verify the effectiveness of this hybrid EMD model, in this section, we take wind station #2 as an example to analyze the performance of the two methods on different IMF components. Table 10 lists the RMSE of SVR-SP and LassoCV on different IMF components. It is worth mentioning that in multi-step prediction, the prediction accuracy of the first step is more important than the other steps, which is of great significance for the accurate estimation of wind power. It can be seen from Table 10 that SVR-SP can obtain significantly better performance than LassoCV at high frequency (IMF1), while LassoCV can obtain better performance at low frequencies (IMF2∼IMF7, Trend), and its RMSE is already close to zero at IMF4. Moreover, SVR-SP has a risk of overfitting when predicting low frequencies, resulting in poor performance. In total, the proposed model that combines the EMD decomposition characteristics and the advantages of the algorithm can achieve better performance than the traditional EMD model.

Comparison of different signal decomposition techniques

Besides EMD, Variational Mode Decomposition (VMD) and Ensemble Empirical Mode Decomposition (EEMD) are also widely used in short-term wind speed forecasting. Here, we analyze the impact of different signal decomposition techniques on the performance of our proposed method. Table 11 shows the prediction performance of the three signal decomposition techniques on two wind stations. For wind station #1, it can be found that compared with VMD and EEMD, EMD obtains the best RMSE value in the 1-step forecasting. The performance obtained by VMD in the 1-step and 2-step forecasting is relatively close, but it drops significantly in the 3-step forecasting. EEMD inherits from EMD, similar to EMD, as the step size increases, the performance will decrease significantly. For wind station #2, EMD also obtained the best predictive performance. VMD has a similar conclusion on wind station #1, and the performance of the 1-step and 2-step forecasting is relatively close. It should be pointed out that in multi-step forecasting, the 1-step forecasting is usually used for wind energy estimation, and other steps are used to assist decision-making, so more attention is paid to the performance of the 1-step forecasting.

The impact of the number of selected features on performance

Feature selection is used to remove redundant features in the study. However, the number of selected significant features will more or less affect the short-term wind speed forecasting. In order to ensure the stability in the complicated industrial system, we analyzed the performance of our proposed method under the different number of selected features. Figure 7 shows the RMSE value between the number of selected features and the performance of our proposed method. It should be pointed out that in the study based on the characteristics of EMD decomposition we use FS and SVR to predict high-frequency component (i.e., IMF₁), and use LassoCV to predict low-frequency components. Feature selection is mainly used in the prediction of IMF₁ component. From Fig. 7, we can be seen that feature selection can slightly improve the performance of 1-step forecasting, but has little effect on 1-step and 2-step forecasting. Overall, as the number of selected features decreases, the generalization performance of the method will improve, but when the selected features are too scarce, the performance will drop sharply due to the deletion of useful features. In order to determine the appropriate number of features, by following (Bradley, Mangasarian & Street, 1998; Chizi, Rokach & Maimon, 2009) , this study uses cross-validation to select.

Performance under different signal-to-noise ratios

In the process of collecting wind speed, it is often affected by the environment and the anemometer itself, resulting in a certain amount of noise in the data. In order to verify the reliability of the method, we analyzed the prediction performance under different signal-to-noise ratios (SNRs). Figure 8 shows the 1-step to 3-step prediction performance of the method from 30∼60db SNR. Take wind station #1 as an example, it can be seen from Fig. 8 that the performance of the proposed method is relatively stable under different signal-to-noise ratios. The RMSE value of 1-step forecasting is about 0.6, the RMSE value of 2-step forecasting is about 0.75, and the RMSE value of 3-step forecasting is about 0.85. In general, as the signal-to-noise ratio increases, the prediction performance of the proposed method will be improved. Similar performance also exists on site #2. These experimental results show that the proposed method can accurately predict wind speed under certain noise.