# Multidimensional Meteorological Variables for Wind Speed Forecasting in Qinghai Region of China: A Novel Approach

Authors: He Jiang; Luo Shihua; Yao Dong
Journal: Advances in Meteorology
 (2020)

---

## Abstract

The accurate, efficient, and reliable forecasting of wind speed is a hot research topic in wind power generation and integration. However, available forecasting models focus on forecasting the wind speed using historical wind speed data and ignore multidimensional meteorological variables. The objective is to develop a hybrid model with multidimensional meteorological variables for forecasting the wind speed accurately. The complementary ensemble empirical mode decomposition with adaptive noise (CEEMDAN) is applied to handle the nonlinearity of the wind speed. Then, the original wind speed will be decomposed into a series of intrinsic model functions with specified numbers of frequencies. A quadratic model that considers the two-way interactions between factors is used to pursue accurate forecasting. To reduce the model complexity, Gram–Schmidt-based feature selection (GSFS) is applied to extract the important meteorological factors. Finally, all the forecasting values of IMFs will be summed by assigning weights that are carefully determined by the whale optimization algorithm (WOA). The proposed forecasting approach has been applied on six datasets that were collected in Qinghai province and is compared with several state-of-the-art wind speed forecasting models. The forecasting results demonstrate that the proposed model can represent the nonlinearity of the wind speed and deliver better results than the competitors.


---

## Body

## 1. Introduction

The depletion of conventional fossil fuels and the deterioration of the ecological environment have been receiving increasing concern in recent years. It is inevitable that fossil fuels will be restricted or replaced by renewable energy resources, which include wind, solar, geothermal, and biomass ones [1, 2]. As a type of clean, inexhaustible, and environmentally preferable renewable energy, wind energy has been developed rapidly since the 1990s around the world and has received widespread attention [3]. Based on the annual report by the Global Wind Energy Council, the global market was 52.573 GW of the wind power installed capacity in 2017, thereby resulting in a cumulative installed capacity of 539.581 GW [4]. Figure 1 presents the top 10 countries in the world in terms of wind power installed capacity from January to December 2017. The wind power is a function of the cube of the wind speed; hence, small differences in the wind speed forecast cause large differences in the wind power forecast [5]. Moreover, the wind speed is an important variable that is related to the wind turbine output power, and the accurate forecast of wind speed facilitates the improvement of the efficiency of wind turbines and the enhancement of the stability of the power supply [6]. However, the wind speed exhibits the stochastic fluctuations and irregular intermittency due to the influence of several meteorological factors, and its distribution is not normal. Therefore, it is also critically challenging to forecast the wind speed accurately. Many scientific techniques have been developed for forecasting the wind speed. These algorithms can be divided into five types: physical-based methods, statistical-based methods, persistence methods, intelligent methods, and hybrid methods [7]. Physical-based methods, which are suitable for long-term wind speed forecasts, apply parametric models according to the physical changes in the atmosphere [8], such as the numerical weather prediction (NWP) model [9]. Statistical-based methods, such as statistical regression [10], the autoregressive moving average (ARMA) models [11], Bayesian models [12], and Markov models [13], depend on the historical wind speed data and forecast the future pattern [6]. These models are suitable for forecasting short-term wind speed but have difficulties if the wind speed series exhibit high-frequency variations. Persistence methods, which are used in short-term wind speed forecasts, are based on the assumption that the wind speed at time k+1 will not change significantly from the wind speed at time k [14]. The intelligent methods apply machine learning and artificial intelligence algorithms to forecast the wind speed and include Kalman filtering [15], artificial neural networks (ANNs) [16], support vector machine [17], extreme learning machine (ELM) [18], and high-dimensional feature selection [19]. The most popular intelligent methods for forecasting wind speed are ANNs. Zhang et al. proposed a new hybrid model that combines a shared-weight long short-term memory network (SWLSTM) and Gaussian process regression (GPR), and they applied SWLSTM-GPR to wind speed prediction cases [20]. Li et al. developed an innovative framework for multistep wind speed prediction that uses a recursive neural network that is based on the wind speed and the turbulence intensity. According to experimental results, the proposed model dramatically outperforms conventional machine learning models on multistep wind speed prediction [21].Figure 1
Distribution of the global wind power installed capacity in 2017.Currently, the main focus has shifted from single methods to advanced hybrid methods, which combine techniques such as statistical and machine learning approaches and utilize the characteristics of individual methods. Sheela and Deepa used a hybrid computing model that integrates a multilayer perceptron network and self-organizing feature maps for wind speed prediction [22–24]. Shukur and Lee proposed a hybrid model that was based on ARIMA for forecasting the daily wind speed in Iraq and Malaysia; to improve the forecasting accuracy, an ANN and a Kalman filter (KF) could be utilized to overcome uncertainty and nonlinearity problems [25]. Doucoure et al. used wavelet decomposition and ANNs with predictability analysis to forecast wind speed time series [26]. Liu et al. developed a novel wind speed prediction model that was based on the WPD (wavelet packet decomposition), a convolutional neural network and a convolutional long short-term memory network, and they demonstrated that the proposed model can perform best in wind speed 1-step to 3-step predictions [27]. Zhu et al. proposed a wind speed forecasting model by integrating Gaussian process regression with a long short-term memory network, and the results demonstrated that the proposed method improves the forecast accuracy [28]. Cai et al. proposed a novel filtering strategy that integrates the statistical predictors and the NWP model outputs into a unified framework. Based on the proposed filtering strategy, a combined predictor, namely, SVR-SDA-UKF (support vector regression, stacked denoising autoencoder, and unscented Kalman filter), is proposed and validated [29].In addition, wind speed time series exhibit nonlinearity and nonstationarity due to the weather and geographical factors. This may lead to incorrect conclusions regarding the activity of the wind speed when the techniques are used to forecast [15]. Therefore, several researchers have attempted to improve the forecasting precision via the addition of a preprocessing stage for time series decomposition and feature selection [30]. Time series decomposition approaches have been widely used to decompose the original complex-structured series into several simple series, such as the wavelet transform (WT) [31] and variational mode decomposition (VMD) [32]. Feature selection approaches facilitate the selection of suitable inputs for the forecasting models, which enhances the performances of the models substantially. Various feature selection approaches have yielded satisfactory results, such as Gram–Schmidt orthogonalization (GSO) [33], the least absolute shrinkage and selection operator (LASSO) [19], and forward regression [34].The main contribution of this paper is the developed hybrid forecasting model, which uses random forest, a Gram–Schmidt-based feature selection method, and the whale optimization algorithm. The proposed method has inherited all the advantages of the above approaches, which are embedded together. Instead of combining them in parallel, random forest and the Gram–Schmidt-based feature selection method are embedded in the whale optimization algorithm; this embedding is a challenging task in terms of computation. The objective function value of the whale optimization algorithm is determined by these two methods. The details will be presented in Section2, which describes the use of the complementary ensemble empirical mode decomposition with adaptive noise (CEEMDAN), random forest (RF), the Gram–Schmidt-based feature selection method (GSFS), and the whale optimization algorithm (WOA) to solve several chained problems.(i)
The first methodology (CEEMDAN in Section2.1) addresses the problem that the signals that are detected from the measured objects are often mixed with various levels of noise. Wavelet analysis can filter out weak or useless signals and extract the information that users need. This paper introduces a data preprocessing strategy that is designed based on the theory of decomposition and reconstruction; it filters the noise from the original wind speed efficiently and searches for the intrinsic values of the wind speed data. Finally, this preprocessing procedure utilizes wavelet packet decomposition for denoising and multiscale wavelet decomposition to complete data preprocessing, and it guarantees the accurate forecasting of the wind speed.(ii)
The second methodology (RF in Section2.2) is a machine learning model that is utilized as a forecasting model after selecting all the important variables. This forecasting model is stable since it weakens the correlations between variables to address the multicollinearity issue. Furthermore, it is designed based on bootstrapping, which theoretically guarantees its stability and accurate forecasting results. In summary, compared with traditional programming methods, the utilized machine learning algorithm simplifies the code and improves the execution performance of the code. Second, machine learning technology can be used to solve problems that cannot be solved via traditional methods. Furthermore, machine learning can adapt to new data; namely, the system can be personalized according to the environment.(iii)
The third methodology (GSFS in Section2.3) is an effective and efficient variable selection approach that is used to select the important variable after the preprocessing procedure. This approach is designed based on an orthogonalization method, which makes the selected variables orthogonal with each other. Furthermore, in contrast to other penalized variable selection methods, this novel approach does not introduce any bias into the model; hence, all nuisance variables can be eliminated effectively.(iv)
The fourth methodology (WOA in Section2.4) is a whale optimization algorithm, which is a new swarm intelligence optimization algorithm that simulates the predation behavior of humpback whales in the ocean. The algorithm has the advantages of simple structure, few parameters, strong search performance, and easy implementation. It is applied in the proposed model in an embedded way to determine the optimal weights for assignment to machine learning models. The fitness function in WOA is determined by the machine learning model and the variable selection method.The hybrid forecasting model, namely, RF-GSFS-WOA, is proposed in Section3. The experimental results and the corresponding analysis are presented in Section 4. Finally, Section 5 presents the conclusions of this study.


## 2. Materials and Methods

In this section, four techniques, namely, CEEMDAN, random forest, GSFS, and WOA, will be introduced. The proposed methodology is developed based on these important techniques.


### 2.1. Complementary Ensemble Empirical Mode Decomposition with Adaptive Noise (CEEMDAN)

The complementary ensemble empirical mode decomposition with adaptive noise(CEEMDAN) originates from the ensemble empirical mode decomposition (EEMD) algorithm. To solve the mode mixing problem, adaptive white noise that has a specified standard deviation can be added at every stage of the decomposition, and a unique residue can be calculated to generate each mode. Moreover, the decomposition is completed with a negligible error; hence, the CEEMDAN provides superior spectral separation of the modes and an accurate reconstruction of the original signal at a lower computational cost [35, 36].The process of the CEEMEDAN is as follows.Step 1: add realizations of a white Gaussian noise series and generate the corresponding signal:(1)yin=yn+ε0ωin,whereyn is the original signal, ii=1,2,…,I is the number of the trials, ε0 is the noise standard deviation, and ωin denotes the realizations of the white Gaussian noise series.Step 2: decomposeI realizations yin via EMD to obtain their first intrinsic model function (IMF):(2)IMF1n=1I∑i=1IIMF1in.Step 3: calculate the first residuer1n (at the first stage k=1):(3)r1n=yn−IMF1n.Step 4: decompose the noise-added residues until their first EMD mode and define the (k+1)-th IMF of the CEEMDAN:(4)IMFk+1n=1I∑i=1IE1rkn+εkEkωin,where the operatorEk⋅ produces the k-th mode of the EMD.Step 5: calculate thek-th residue:(5)rkn=rk−1n−IMFkn,wherek=2,…,K, K is the total number of modes.Step 6: repeat step 4 and step 5 for the nextk until the obtained residue is not suitable for decomposition. The final residue is expressed as(6)Rn=yn−∑i=1KIMFin.Hence, the original signal can be expressed as(7)yn=∑i=1KIMFin+Rn.


### 2.2. Random Forest

Random forest is an ensemble learning method that uses decision trees as base learners and measures the importance of each variable. Random forest is a variant of bagging (bootstrap aggregating), which considers allm attributes or variables when establishing a base learner. Decision trees tend to overfit the data, thereby leading to low bias but high variance. To overcome this drawback, random forest constructs decision trees and integrates the results of all individual trees to generate the final outcome [37]. In classification, the mode of the classes that are predicted by the trees is adopted as the output, while, in regression, the mean of the predictions of all individual trees is used.Since bootstrapping is used to generate the sample from which an individual tree is grown, only approximately 63.2% of the observations are used to produce each base learner, and the remainder of the training set can be used as a validation set: suppose that there arem samples. The probability that not all the samples are chosen in the training set is 1−1/mm which is 1/e≈36.8% if m⟶∞. For data z=x,y, with x representing the observations and y representing the response variable, the out-of-bag (OOB) estimate, which is denoted by HOOBx, is produced by base learners [32, 33]. Given an ensemble of classifiers h1x,h2x,...,hTx, it is defined as follows:(8)HOOBx=argminy∈Y∑t=1TIhtx=y⋅Ix∉Dt,where Y is the set that consists of all possible values of the response and Dt is the subsample that is applied to produce base learner ht. The OOB estimate of the generalization error is expressed as(9)εOOB=1D∑x,y∈DIHOOBx≠y,where D is the training set.Both bagging and random forest measure the importance of each variable via the following procedure: (1) calculate the OOB error of each observation and average them to obtain the OOB error of the forest. (2) Permute the values of thej-th feature among the training set. (3) Compute the OOB error again on this perturbed data set, and the difference in the OOB error before and after the permutation can be used to measure the importance of the j-th variable, where a larger difference corresponds to higher importance of variable j.


### 2.3. Gram–Schmidt-Based Feature Selection Method (GSFS)

Feature selection substantially facilitates the reduction of the model complexity and the enhancement of the forecasting accuracy when establishing a forecasting model. Many feature selection methods are designed based on penalized regression, which shrinks the estimate; these methods yield inaccurate results. In this study, a subset selection method, namely, the Gram–Schmidt-based feature selection method (GSFS) [38], is used for feature selection. The objective of GSFS is to identify all the relevant features that contribute to the model forecasting. The details of GSFS are presented as follows.GSFS iterations: letG0 be the empty set and initialize y^0 and k to 0. Here, y^0 denotes the initial forecasting values. Calculate the pseudoresponse Rtk=yt−y^k and select the optimal variable xj^k+1 by regressing Rtk on x:(10)j^k+1=argmin∑t=1nRtk−δjkxtj2,where δjk=∑t=1nRtkxtj/∑t=1nxtj2. Then, G1=G0∪xj^1 and let y^k+1=y^k+δj^k+1kxj^k+1. The objective in building G1 is to add the optimal variables into the empty set.In the nextMn=n/logp1/2/2 steps, we will conduct orthogonalization before selecting each new feature. If k≥1, orthogonalize Gk=xj^1,xj^2,…,xj^k to Gorthk=xj^1,xj^2⊥,…,xj^k⊥ via the Gram–Schmidt method. Similar to the first step, compute the pseudoresponse Otk=yt−y^k and consider(11)j^k+1=argmin∑t=1nOtk−ηjkxtj⊥2,where ηjk=∑t=1nOtkxtj⊥/∑t=1nxtj⊥2. Update y^k+1=y^k+ηj^k+1kxj^k+1⊥ and Gk+1=Gk∪xj^k+1. The objective of this procedure is to select variables that are orthogonal with each other sequentially and add them into an empty set.The main advantage of GSFS is that it does not involve any penalized regression. It detects the relevant features by examining each feature sequentially. Although it may take more time, it will not miss any important data information. In the second step, GSFS computes the pseudoresponse as the target vector, which can be used to find the most strongly correlated features with the pseudoresponse. This pseudoresponse differs from the ordinary response that is applied in regression or classification problems. It is updated in every iteration of GSFS to eliminate the effects of selected features. An orthogonalization procedure that is based on the Gram–Schmidt method is applied to make all the features orthogonal with each other. This procedure removes the dependencies among features and shortens the computation time for inverting the projection matrix.The determination ofMn plays an important role in guaranteeing the accuracy and the selection consistency. Let p=pn⟶∞; namely, let the number of features increase with the sample size. We must impose the following regularity condition before proving the theorem.(C1)
logpn=on. It follows that limn⟶∞logpn/n=0.Assume that the noise term is independent of each feature and letu and u0 be constants; namely,(C2)
Eexpuε<∞ for u<u0Letσj2=Exj2, which represents the noise level and let zj=xj/σj and ztj=xtj/σj. Assume that there exists u1 such that(C3)
limsupn⟶∞max1≤j≤pnEexpu1zj2<∞.(C4)
supn≥1∑j=1pnδjσj<∞.(C5)
There exists0≤γ≤1 such that nγ=on/logpn1/2 and liminfn⟶∞nγmin1≤j≤pn:δj≠0δj2σj2>0.Theorem 1.
Suppose that regularity conditions (C1)–(C5) are satisfied. LetMn/nγ⟶∞ and Mn=On/logpn1/2. Then, limn⟶∞PNn=SMn=1, where Nn=1≤j≤pn:δj∗≠0 denotes the set of important variables.Theorem1 establishes the selection consistency of the GSFS method and provides a theoretical guarantee on the choice of Mn. Due to the selection consistency, the feature selection method can detect all the important features as the sample size approaches infinity. Here, GMn contains the features that are selected in step Mn. To select the optimal model with all the important features, high-dimensional information criteria (HDIC) will be applied after GSFS. The proof of Theorem 1 is similar to [34] and is omitted here.


### 2.4. Whale Optimization Algorithm (WOA)

In 2016, Mirjalili and Lewis [39] proposed the whale optimization algorithm (WOA), which is a novel metaheuristic that mimics the bubble-net foraging behavior of humpback whales in searching for their prey. It is illustrated in Figure 2(a). The humpback whales attack small fishes or schools of krill that are close to the surface by swimming in a shrinking circle and creating distinctive bubbles along a spiral-shaped path. The algorithm includes two phases: the first phase is encircling a prey and applying the spiral bubble-net attacking method, and the second phase is searching randomly for a prey. The details of each phase are presented as follows.Figure 2
Bubble-net search mechanism that is implemented in WOA: (a) the bubble-net feeding behavior of humpback whales, (b) the shrinking encircling mechanism, and (c) the spiral position updating.


#### 2.4.1. Encircling Prey

The humpback whales dive downwards and generate bubbles in a spiral shape, which makes the prey (typically krill) feel threatened and, hence, gather together. It is assumed that there areN humpback whales in a d-dimensional search space. The position of each whale is denoted by Xi=Xi1,Xi2,...,Xidi=1,2,...,N. The position of the prey represents the globally optimal solution to the problem. The humpback whales can identify the prey’s position and encircle it. However, there is no priori information about the optimal solution. Therefore, in the WOA, the current best position is adopted as the position that corresponds to the optimal solution. Once this position has been identified, other humpback whales will gather around it and update their position as follows [36]:(12)Xt+1=X∗t−A⋅D,(13)D=C⋅X∗t−Xt,where t is the current iteration, Xt is the current position, X∗t is the current best position, and A and C are the coefficient vectors:(14)A=2a⋅rand1−a,C=2⋅rand2,where rand1 and rand2 are 2 random numbers that are in the interval (0, 1) and a is a convergence factor, which decreases linearly as the iteration proceeds. Define t as the iteration number and tmax as the maximum allowed number of iterations. Factor a is defined as(15)a=2−2ttmax.Two mechanisms are designed to describe the hunting behavior of the humpback whales, namely, the shrinking encircling mechanism and the spiral position updating mechanism, as illustrated in Figures2(b) and 2(c). Shrinking encircling is realized as the convergence factor a decreases, while, in spiral position updating, the spiral movement of the humpback whales is simulated by the following:(16)Xt+1=D′⋅ebl⋅cos2πl+X∗t,where D′=X∗t−Xt is the distance between the humpback and the prey, b is the shape parameter of the spiral, and l is a random number that is in the interval (−1, 1). Since the humpback whales perform these two behaviors simultaneously, in each iteration, each of these mechanisms has a probability of 0.5 to be utilized. The corresponding model is as follows:(17)Xt+1=X∗t−A⋅D,fp<0.5,D′⋅ebl⋅cos2πl+X∗t,fp≥0.5where p is a random number in (0, 1).


#### 2.4.2. Search for Prey

To enhance the exploration in WOA, instead of updating the positions of the humpback whales based on the best position so far, a random position is chosen to guide the search accordingly.A with random values that are greater than 1 or less than −1 is applied to force the position to move far away from the best position.The model of searching for prey is defined as [40](18)D″=C⋅Xrand−X,(19)Xt+1=Xrand−A⋅D″,where Xrand is a random position (or a random whale) that is chosen from the current population.


## 2.1. Complementary Ensemble Empirical Mode Decomposition with Adaptive Noise (CEEMDAN)

The complementary ensemble empirical mode decomposition with adaptive noise(CEEMDAN) originates from the ensemble empirical mode decomposition (EEMD) algorithm. To solve the mode mixing problem, adaptive white noise that has a specified standard deviation can be added at every stage of the decomposition, and a unique residue can be calculated to generate each mode. Moreover, the decomposition is completed with a negligible error; hence, the CEEMDAN provides superior spectral separation of the modes and an accurate reconstruction of the original signal at a lower computational cost [35, 36].The process of the CEEMEDAN is as follows.Step 1: add realizations of a white Gaussian noise series and generate the corresponding signal:(1)yin=yn+ε0ωin,whereyn is the original signal, ii=1,2,…,I is the number of the trials, ε0 is the noise standard deviation, and ωin denotes the realizations of the white Gaussian noise series.Step 2: decomposeI realizations yin via EMD to obtain their first intrinsic model function (IMF):(2)IMF1n=1I∑i=1IIMF1in.Step 3: calculate the first residuer1n (at the first stage k=1):(3)r1n=yn−IMF1n.Step 4: decompose the noise-added residues until their first EMD mode and define the (k+1)-th IMF of the CEEMDAN:(4)IMFk+1n=1I∑i=1IE1rkn+εkEkωin,where the operatorEk⋅ produces the k-th mode of the EMD.Step 5: calculate thek-th residue:(5)rkn=rk−1n−IMFkn,wherek=2,…,K, K is the total number of modes.Step 6: repeat step 4 and step 5 for the nextk until the obtained residue is not suitable for decomposition. The final residue is expressed as(6)Rn=yn−∑i=1KIMFin.Hence, the original signal can be expressed as(7)yn=∑i=1KIMFin+Rn.


## 2.2. Random Forest

Random forest is an ensemble learning method that uses decision trees as base learners and measures the importance of each variable. Random forest is a variant of bagging (bootstrap aggregating), which considers allm attributes or variables when establishing a base learner. Decision trees tend to overfit the data, thereby leading to low bias but high variance. To overcome this drawback, random forest constructs decision trees and integrates the results of all individual trees to generate the final outcome [37]. In classification, the mode of the classes that are predicted by the trees is adopted as the output, while, in regression, the mean of the predictions of all individual trees is used.Since bootstrapping is used to generate the sample from which an individual tree is grown, only approximately 63.2% of the observations are used to produce each base learner, and the remainder of the training set can be used as a validation set: suppose that there arem samples. The probability that not all the samples are chosen in the training set is 1−1/mm which is 1/e≈36.8% if m⟶∞. For data z=x,y, with x representing the observations and y representing the response variable, the out-of-bag (OOB) estimate, which is denoted by HOOBx, is produced by base learners [32, 33]. Given an ensemble of classifiers h1x,h2x,...,hTx, it is defined as follows:(8)HOOBx=argminy∈Y∑t=1TIhtx=y⋅Ix∉Dt,where Y is the set that consists of all possible values of the response and Dt is the subsample that is applied to produce base learner ht. The OOB estimate of the generalization error is expressed as(9)εOOB=1D∑x,y∈DIHOOBx≠y,where D is the training set.Both bagging and random forest measure the importance of each variable via the following procedure: (1) calculate the OOB error of each observation and average them to obtain the OOB error of the forest. (2) Permute the values of thej-th feature among the training set. (3) Compute the OOB error again on this perturbed data set, and the difference in the OOB error before and after the permutation can be used to measure the importance of the j-th variable, where a larger difference corresponds to higher importance of variable j.


## 2.3. Gram–Schmidt-Based Feature Selection Method (GSFS)

Feature selection substantially facilitates the reduction of the model complexity and the enhancement of the forecasting accuracy when establishing a forecasting model. Many feature selection methods are designed based on penalized regression, which shrinks the estimate; these methods yield inaccurate results. In this study, a subset selection method, namely, the Gram–Schmidt-based feature selection method (GSFS) [38], is used for feature selection. The objective of GSFS is to identify all the relevant features that contribute to the model forecasting. The details of GSFS are presented as follows.GSFS iterations: letG0 be the empty set and initialize y^0 and k to 0. Here, y^0 denotes the initial forecasting values. Calculate the pseudoresponse Rtk=yt−y^k and select the optimal variable xj^k+1 by regressing Rtk on x:(10)j^k+1=argmin∑t=1nRtk−δjkxtj2,where δjk=∑t=1nRtkxtj/∑t=1nxtj2. Then, G1=G0∪xj^1 and let y^k+1=y^k+δj^k+1kxj^k+1. The objective in building G1 is to add the optimal variables into the empty set.In the nextMn=n/logp1/2/2 steps, we will conduct orthogonalization before selecting each new feature. If k≥1, orthogonalize Gk=xj^1,xj^2,…,xj^k to Gorthk=xj^1,xj^2⊥,…,xj^k⊥ via the Gram–Schmidt method. Similar to the first step, compute the pseudoresponse Otk=yt−y^k and consider(11)j^k+1=argmin∑t=1nOtk−ηjkxtj⊥2,where ηjk=∑t=1nOtkxtj⊥/∑t=1nxtj⊥2. Update y^k+1=y^k+ηj^k+1kxj^k+1⊥ and Gk+1=Gk∪xj^k+1. The objective of this procedure is to select variables that are orthogonal with each other sequentially and add them into an empty set.The main advantage of GSFS is that it does not involve any penalized regression. It detects the relevant features by examining each feature sequentially. Although it may take more time, it will not miss any important data information. In the second step, GSFS computes the pseudoresponse as the target vector, which can be used to find the most strongly correlated features with the pseudoresponse. This pseudoresponse differs from the ordinary response that is applied in regression or classification problems. It is updated in every iteration of GSFS to eliminate the effects of selected features. An orthogonalization procedure that is based on the Gram–Schmidt method is applied to make all the features orthogonal with each other. This procedure removes the dependencies among features and shortens the computation time for inverting the projection matrix.The determination ofMn plays an important role in guaranteeing the accuracy and the selection consistency. Let p=pn⟶∞; namely, let the number of features increase with the sample size. We must impose the following regularity condition before proving the theorem.(C1)
logpn=on. It follows that limn⟶∞logpn/n=0.Assume that the noise term is independent of each feature and letu and u0 be constants; namely,(C2)
Eexpuε<∞ for u<u0Letσj2=Exj2, which represents the noise level and let zj=xj/σj and ztj=xtj/σj. Assume that there exists u1 such that(C3)
limsupn⟶∞max1≤j≤pnEexpu1zj2<∞.(C4)
supn≥1∑j=1pnδjσj<∞.(C5)
There exists0≤γ≤1 such that nγ=on/logpn1/2 and liminfn⟶∞nγmin1≤j≤pn:δj≠0δj2σj2>0.Theorem 1.
Suppose that regularity conditions (C1)–(C5) are satisfied. LetMn/nγ⟶∞ and Mn=On/logpn1/2. Then, limn⟶∞PNn=SMn=1, where Nn=1≤j≤pn:δj∗≠0 denotes the set of important variables.Theorem1 establishes the selection consistency of the GSFS method and provides a theoretical guarantee on the choice of Mn. Due to the selection consistency, the feature selection method can detect all the important features as the sample size approaches infinity. Here, GMn contains the features that are selected in step Mn. To select the optimal model with all the important features, high-dimensional information criteria (HDIC) will be applied after GSFS. The proof of Theorem 1 is similar to [34] and is omitted here.


## 2.4. Whale Optimization Algorithm (WOA)

In 2016, Mirjalili and Lewis [39] proposed the whale optimization algorithm (WOA), which is a novel metaheuristic that mimics the bubble-net foraging behavior of humpback whales in searching for their prey. It is illustrated in Figure 2(a). The humpback whales attack small fishes or schools of krill that are close to the surface by swimming in a shrinking circle and creating distinctive bubbles along a spiral-shaped path. The algorithm includes two phases: the first phase is encircling a prey and applying the spiral bubble-net attacking method, and the second phase is searching randomly for a prey. The details of each phase are presented as follows.Figure 2
Bubble-net search mechanism that is implemented in WOA: (a) the bubble-net feeding behavior of humpback whales, (b) the shrinking encircling mechanism, and (c) the spiral position updating.


### 2.4.1. Encircling Prey

The humpback whales dive downwards and generate bubbles in a spiral shape, which makes the prey (typically krill) feel threatened and, hence, gather together. It is assumed that there areN humpback whales in a d-dimensional search space. The position of each whale is denoted by Xi=Xi1,Xi2,...,Xidi=1,2,...,N. The position of the prey represents the globally optimal solution to the problem. The humpback whales can identify the prey’s position and encircle it. However, there is no priori information about the optimal solution. Therefore, in the WOA, the current best position is adopted as the position that corresponds to the optimal solution. Once this position has been identified, other humpback whales will gather around it and update their position as follows [36]:(12)Xt+1=X∗t−A⋅D,(13)D=C⋅X∗t−Xt,where t is the current iteration, Xt is the current position, X∗t is the current best position, and A and C are the coefficient vectors:(14)A=2a⋅rand1−a,C=2⋅rand2,where rand1 and rand2 are 2 random numbers that are in the interval (0, 1) and a is a convergence factor, which decreases linearly as the iteration proceeds. Define t as the iteration number and tmax as the maximum allowed number of iterations. Factor a is defined as(15)a=2−2ttmax.Two mechanisms are designed to describe the hunting behavior of the humpback whales, namely, the shrinking encircling mechanism and the spiral position updating mechanism, as illustrated in Figures2(b) and 2(c). Shrinking encircling is realized as the convergence factor a decreases, while, in spiral position updating, the spiral movement of the humpback whales is simulated by the following:(16)Xt+1=D′⋅ebl⋅cos2πl+X∗t,where D′=X∗t−Xt is the distance between the humpback and the prey, b is the shape parameter of the spiral, and l is a random number that is in the interval (−1, 1). Since the humpback whales perform these two behaviors simultaneously, in each iteration, each of these mechanisms has a probability of 0.5 to be utilized. The corresponding model is as follows:(17)Xt+1=X∗t−A⋅D,fp<0.5,D′⋅ebl⋅cos2πl+X∗t,fp≥0.5where p is a random number in (0, 1).


### 2.4.2. Search for Prey

To enhance the exploration in WOA, instead of updating the positions of the humpback whales based on the best position so far, a random position is chosen to guide the search accordingly.A with random values that are greater than 1 or less than −1 is applied to force the position to move far away from the best position.The model of searching for prey is defined as [40](18)D″=C⋅Xrand−X,(19)Xt+1=Xrand−A⋅D″,where Xrand is a random position (or a random whale) that is chosen from the current population.


## 2.4.1. Encircling Prey

The humpback whales dive downwards and generate bubbles in a spiral shape, which makes the prey (typically krill) feel threatened and, hence, gather together. It is assumed that there areN humpback whales in a d-dimensional search space. The position of each whale is denoted by Xi=Xi1,Xi2,...,Xidi=1,2,...,N. The position of the prey represents the globally optimal solution to the problem. The humpback whales can identify the prey’s position and encircle it. However, there is no priori information about the optimal solution. Therefore, in the WOA, the current best position is adopted as the position that corresponds to the optimal solution. Once this position has been identified, other humpback whales will gather around it and update their position as follows [36]:(12)Xt+1=X∗t−A⋅D,(13)D=C⋅X∗t−Xt,where t is the current iteration, Xt is the current position, X∗t is the current best position, and A and C are the coefficient vectors:(14)A=2a⋅rand1−a,C=2⋅rand2,where rand1 and rand2 are 2 random numbers that are in the interval (0, 1) and a is a convergence factor, which decreases linearly as the iteration proceeds. Define t as the iteration number and tmax as the maximum allowed number of iterations. Factor a is defined as(15)a=2−2ttmax.Two mechanisms are designed to describe the hunting behavior of the humpback whales, namely, the shrinking encircling mechanism and the spiral position updating mechanism, as illustrated in Figures2(b) and 2(c). Shrinking encircling is realized as the convergence factor a decreases, while, in spiral position updating, the spiral movement of the humpback whales is simulated by the following:(16)Xt+1=D′⋅ebl⋅cos2πl+X∗t,where D′=X∗t−Xt is the distance between the humpback and the prey, b is the shape parameter of the spiral, and l is a random number that is in the interval (−1, 1). Since the humpback whales perform these two behaviors simultaneously, in each iteration, each of these mechanisms has a probability of 0.5 to be utilized. The corresponding model is as follows:(17)Xt+1=X∗t−A⋅D,fp<0.5,D′⋅ebl⋅cos2πl+X∗t,fp≥0.5where p is a random number in (0, 1).


## 2.4.2. Search for Prey

To enhance the exploration in WOA, instead of updating the positions of the humpback whales based on the best position so far, a random position is chosen to guide the search accordingly.A with random values that are greater than 1 or less than −1 is applied to force the position to move far away from the best position.The model of searching for prey is defined as [40](18)D″=C⋅Xrand−X,(19)Xt+1=Xrand−A⋅D″,where Xrand is a random position (or a random whale) that is chosen from the current population.


## 3. Proposed Hybrid Forecasting Model: RF-GSFS-WOA

In this paper, the proposed RF-GSFS-WOA is defined by the following procedure.Step 1: the original wind speed time seriesyn is decomposed into I IMFs (let IMFI = R (n)) by following Step 1–Step 6 in Section 2.1.Step 2: a quadratic model (QM) is established for eachIMFi:(20)IMF1=a11x1+a21x2+…+ap1xp+a111x12+a121x1x2+…+app1xp2+ε1,⋮IMFi=a1ix1+a2ix2+…+apixp+a11ix12+a12ix1x2+…+appixp2+ε2,⋮IMFI=a1Ix1+a2Ix2+…+apIxp+a11Ix12+a12Ix1x2+…+appIxp2+εn,wherex1,…,xp are variables (first order) in the original dataset; xjxk are quadratic variables (second order) that are constructed by calculating the component-wise product of xj and xk; a1,a2,…,ap,…,a11,…,app are the corresponding coefficients of the variables and quadratic variables; and ε1,…,εn are independent Gaussian noise errors that follow N0,σ2. I quadratic models will be obtained, which are denoted by QM1,…,QMI. As p2+p variables are included in the QM, the computational cost is high even if the value of p is moderate. Therefore, a variable selection method is urgently needed to reduce the computational cost.Step 3: to reduce the computational cost of QM in Step 2, GSFS is used to complete the variable selection task. With slight abuse of notation,IMFi1≤i≤I is denoted by IMF. Let IMF∧0=0,k=0, where IMF∧0 represents the initial forecasting value, Rtk=IMFt−IMF∧k, and the optimal variables are selected by regressing Rtk on Xo:=x1,…,xp,x12,x1x2,…,xp2=xjo1≤j≤p2+p, which is expressed as(21)j^k+1=argmin∑t=1nRtk−δjkxtjo2,whereδjk=∑t=1nRtkxtjo/∑t=1nxtjo2. Then, G1=G0∪xj^1o and let IMF∧k+1=IMF∧k+δj^k+1kxj^1o. In the next Mn=n/logp1/2/2 steps, the data will be orthogonalized. If k≥1, orthogonalize Gk=xj^1o,xj^2o,…,xj^ko to Gorthk=xj^1o,xj^2o⊥,…,xj^ko⊥ via the Gram–Schmidt method. Similar to the first step, compute the pseudoresponse Otk=IMFt−IMF∧k and consider(22)j^k+1=argmin∑t=1nOtk−γjkxtjo⊥2,whereγjk=∑t=1nOtkxtjo⊥/∑t=1nxtjo⊥2. Update IMF∧k+1=IMF∧k+γj^k+1kxj^k+1o⊥ and Gk+1=Gk∪xj^k+1o.Step 4: random forest modelRFi is built over the variables that were selected in Step 3.Step 5: assign a weightωi to each RFi, and the final forecast model is provided by the ensemble model (EM):(23)EM=∑i=1IωiRFi,with ∑i=1Iωi=1, where ωi is obtained via the WOA algorithm.The RF-GSFS-WOA algorithm is presented in detail as Algorithm1. The proposed combination forecasting method, namely, RF-GSFS-WOA, has the following advantages. Instead of applying a single decision tree, RF grows many decision trees and, hence, can handle the correlation among the features. It bootstraps samples, which promotes the diversity of the base learners and, consequently, reduces the number of variables of the model. Furthermore, RF outperforms bagging as the number of base learners increases. It trains fast since it applies a fraction of the features instead of considering all possible features. The time that is required for training RF is often shorter than that for bagging. The overfitting problem of RF can be controlled well via the selection of a suitable number of features in each decision tree. GSFS is a greedy feature selection method, which is a type of forward selection method. In contrast to penalized feature selection methods such as LASSO and SCAD, it does not introduce any bias into the model. WOA is a metaheuristic method that can obtain the globally optimal solutions of nonlinear optimization problems. A flowchart of the proposed RF-GSFS-WOA is presented in Figure 3. The RF-GSFS-WOA algorithm is presented in the Appendix.Algorithm 1: RF-GSFS-WOA algorithm.
(i)
INPUTS: D=X,y, the wind speed dataset.
(ii)
PARAMETERS:
Nstd: noise standard deviation in CEEMDAN
NR: number of realizations in CEEMDAN
MaxIter: maximum number of sifting iterations that are allowed in CEEMDAN
q: variable that is selected in GSFS
MaxIterRF: maximum number of iterations in random forest
NTress: number of trees in random forest
N: number of whales in WOA
MaxIterWOA: maximum number of iterations in WOA
(iii)
OUTPUTS : the forecasting error
(1)
The original datasetD is divided into training dataset Dtrn=Xtrn,ytrn and test dataset Dtst=Xtst,ytst.
(2)
Decomposition using CEEMDAN, which results in I IMFs:IMF1,IMF2,...,IMFI.
(3)
Establish a quadratic model ofXtrn and Xtst, which results in Xtrno and Xtsto
(4)
INITIALIZATION: setG0=∅, IMF∧0=0,k=0.
(5)
IFk=0 DO
(6)
Calculate the pseudoresponse
(7)
Perform regression on the data with a pseudoresponse and select the most strongly correlated variable
(8)
Update the responseIMF∧k and Gk+1=Gk∪xj^k+1o.
(9)
END IF
(10)
FORk=1 to Mn DO
(11)
OrthogonalizeGk via the Gram–Schmidt method.
(12)
Perform regression on the data with a pseudoresponse and select the most strongly correlated variable
(13)
Update the responseIMF∧k and Gk+1=Gk∪xj^k+1o.
(14)
END FOR
(15)
Build the solution pathG=Gk:1≤k≤Mn and use HDIC to select the optimal solution, including q important variables.
(16)
Grow a random forest with NTress = 3000 trees using the selected variables with MaxIterRF iterations.
(17)
Initialize the whale population for weight matrixW∈RN×I in which the columns are weights for established RF models.
(18)
Calculate the fitness function usingW and RF models
(19)
LetW∗ be the current best weights
(20)
WHILEt<MaxIterWOA DO
(21)
FOR each weight DO
(22)
Updatea,A,C,l and p
(23)
IFp<0.5 THEN
(24)
IFA|<1 THEN
(25)
Update the positions of the current weight via (12) and (13)
(26)
ELSE IFA|≥1 THEN
(27)
Randomly select a search agent
(28)
Update the position of the current search agent via (17)
(29)
END IF
(30)
ELSE IFp≥0.5 THEN
(31)
Update the position of the current weight via (16)
(32)
END FOR
(33)
Check if any search agent goes beyond the search space
(34)
Calculate the fitness function of each weight
(35)
UpdateW∗ if there is a better solution
(36)
t
⟵
t
+
1
(37)
END WHILE
(38)
Obtain the optimal weights
(39)
Calculate the forecasting error using the optimal weights and the RF model that is generated byXtsto and ytst.Figure 3
Flowchart of the developed model: RF-GSFS-WOA.


## 4. Case Studies

This section describes the data collection process, establishes the forecasting assessment criteria, and presents, discusses, and analyzes the experiment results. The main objective of the experiments is to evaluate the performance of the proposed hybrid forecasting model, namely, RF-GSFS-WOA.


### 4.1. Data Collection

Qinghai province is located in the northeast of the Qinghai-Tibetan Plateau, and most of the area is rich in wind energy resources; hence, it is suitable for the installation of large-scale wind farm. Thus, it is meaningful to investigate sites in this region. The hourly data were collected between 0 : 00 and 24 : 00 from 1/1/2014 to 12/31/2014 at six sites in the Qinghai region. The dataset that is used in this study was collected by National Renewable Energy Laboratory (NREL) and is available athttps://www.nrel.gov/gis/tools.html. To simplify the description, the six sites are denoted by Site 1, Site 2, Site 3, Site 4, Site 5, and Site 6. The longitudes and latitudes of these locations are specified in Figure 4. Five meteorological factors, namely, the temperature (T), pressure (P), relative humidity (RH), wind precipitable water (PW), and wind direction (WD), are considered in the forecasting of the wind speed for the Qinghai region, which facilitates the design of wind power plants to efficiently generate electricity for the local population [41]. There are four seasons in each year: spring (Mar., Apr., and May), summer (Jun., Jul., and Aug.), autumn (Sep., Oct., and Nov.), and winter (Dec., Jan., and Feb.). The data are divided into training data and forecasting data as follows: 20 days in each season are randomly chosen as the training data (20 days × 24 hours/days = 480) for building the model, and the forecasting data include 5 days (5 days × 24 hours/day = 120) that are randomly chosen from the remaining days in each season. The performances of various forecasting models are evaluated in the four seasons. The training sets and forecasting sets for the six sites are described in Table 1. The forecasting sets account for approximately 20% of the total data, which is reasonable.Figure 4
Distribution and location information of the observation sites.Table 1
Training sets and forecasting sets for hourly wind speed data at six sites.
Season
Training set (size: 20 × 24 = 480)
Forecasting set (size: 5 × 24 = 120)
Relative meteorological factors
Spring (Mar., Apr., and May)
Mar. 9th, Mar. 12th, Mar. 19th, and Mar. 21st
May 16th and May 20th
Temperature, pressure, precipitable water, and wind direction
Mar. 26th, Mar. 29th, Apr. 5th, and Apr. 6th
May 25th, May 28th, and May 29th
Apr. 13th, Apr. 15th, Apr. 16th, and Apr. 21st
Apr. 22nd, Apr. 23th, Apr. 24th, and Apr. 25th
Apr. 28th, May 1st, May 8th, and May 9th
Summer (Jun., Jul., and Aug.)
Jun. 1st, Jun. 3rd, Jun. 4th, and Jun. 6th
Aug. 12th and Aug. 20th
Jun. 13th, Jun. 16th, Jun. 22nd, and Jun. 28th
Aug. 22nd, Aug. 24th, and Aug. 25th
Jun. 30th, Jul. 2nd, Jul. 7th, and Jul. 13th
Jul. 15th, Jul. 18th, Jul. 22nd, and Jul. 28th
Jul. 31st, Aug. 1st, Aug. 3rd, and Aug. 11th
Autumn (Sep., Oct., and Nov.)
Sep. 3rd, Sep. 6th, Sep. 8th, and Sep. 11th
Nov. 18th and Nov. 23rd
Sep. 14th, Sep. 15th, Sep. 21st, and Sep. 24th
Nov. 24th, Nov. 27th, and Nov. 29th
Sep. 30th, Oct. 4th, Oct. 6th, and Oct. 13th
Oct. 18th, Oct. 31st, Nov. 2nd, and Nov. 3rd
Nov. 4th, Nov. 8th, Nov. 12th, and Nov. 14th
Winter (Dec., Jan., and Feb.)
Dec. 2nd, Dec. 3rd, Dec. 4th, and Dec 8th
Feb. 7th and Feb. 8th
Dec. 9th, Dec. 13th, Dec. 14th, and Dec. 18th
Feb. 15th, Feb. 24th, and Feb. 26th
Dec. 20th, Dec. 25th, Dec. 26th, and Dec. 29th
Dec. 30th, Jan. 5th, Jan. 8th, and Jan 9th
Jan. 16th, Jan. 18th, Jan. 28th, and Feb. 5th


### 4.2. Forecasting Assessment Criteria

To evaluate the results of the forecasting models quantitatively, four assessment measures, namely, the mean absolute error (MAE), root mean square error (RMSE), mean absolute percent error (MAPE), and Theil inequality coefficient (TIC), are applied. All are calculated between the actual data and the forecasted values, and the smaller their values, the better the forecasting results. These criteria are expressed as follows:(24)MAE=1N∑i=1Nyi−y^i,RMSE=1N∑i=1Nyi−y^i2,MAPE=1N∑i=1Nyi−y^iyi×100%,TIC=1/N∑i=1Nyi−y^i21/N∑i=1Ny^i2+1/N∑i=1Nyi2×100%,where N is the size of the forecasting sample and yi and y^i are the actual value and the forecasted value, respectively, for time period i.


### 4.3. Experiment Results and Corresponding Analysis

The experimental settings are described in detail in Table2. The parameters in the complementary ensemble empirical mode decomposition with adaptive noise (CEEMDAN) are set as follows: the noise standard deviation (Nstd) ε0 is set as 0.2, the number of realizations (Nr) is set as 500, and the maximum allowed number of sifting iterations (MaxIter) is set as 500. The regularization parameters in LASSO and SCAD are chosen via 10-fold cross-validation. In the experiments on the neural networks, the maximum number of iterations (MaxIterNN) is set as 2000 for both BP and Elman. Each activation function from the hidden layer to the output layer is set as the Tansig function (Func). The number of hidden neurons is crucial since it mainly determines the network complexity. To determine the optimal number of hidden neurons, values are selected from the following set: 5,8,10,15,20,25. Furthermore, Elman is trained using backpropagation through time, and the corresponding weights are selected based on a uniform distribution. Based on the experimental results, the gradient descent with momentum is set as 0.82, and the adaptive learning rate is set as 0.02. In GSFS, we set q as 30; hence, 30 variables will be chosen after the GSFS procedure, and, thus, there are 1200 variables in the quadratic model. The number of trees in RF is set as 3000, and the maximum number of iterations (MaxIterRF) is 500. RF is implemented using MATLAB function TreeBagger. The fitness function of WOA is selected as “F1,” which is the square error loss function, and the range is (0, 1). The weights are determined after normalization such that the sum of all the weights is 1.Table 2
Important parameters of seven forecasting models.
Nstd
Nr
MaxIter
λ
MaxIterNN
Func
h
q
MaxIterRF
Elman
0.2
500
500
–
2000
Tansig
10
–
–
Elman-LASSO
0.2
500
500
0.875
2000
Tansig
10
–
–
Elman-SCAD
0.2
500
500
0.625
2000
Tansig
10
–
–
BP
0.2
500
500
–
2000
Tansig
10
–
–
BP-LASSO
0.2
500
500
0.875
2000
Tansig
10
–
–
BP-SCAD
0.2
500
500
0.625
2000
Tansig
10
–
–
RF-GSFS-WOA
0.2
500
500
–
–
Tansig
-
30
500To evaluate the forecasting performance of the proposed hybrid model for forecasting hourly wind speed, experiments for six sites and four seasons are designed in this paper. Four types of assessment criteria are used in all experiments to evaluate the performance. The proposed hybrid forecasting model is RF-GSFS-WOA, and six forecasting models, namely, an Elman neural network, Elman-LASSO, Elman-SCAD, a BP neural network (BP), BP-LASSO, and BP-SCAD, are applied as comparison models to evaluate the performance of the developed hybrid model. Wind speed data exhibit volatility or instability due to changes in the weather; hence, this study conducts a preprocessing step before building all the models; namely, the complementary ensemble empirical mode decomposition with adaptive noise(CEEMDAN) is used to eliminate the noise and smooth the wind speed time series data before the developed hybrid model and six comparison models are established for hourly wind speed forecasting in four seasons at six sites.In each season, 20 days of hourly wind speed time series are selected randomly as the training phase for the construction of the models, and 5 days of wind speed time series are selected randomly as the forecasting phase. Thus, 80% of the data will be used as training data; the remaining data (20%) will be used as test data. The reason for this splitting scheme is that the model can be trained sufficiently with more observations (80% of the data) and overfitting can be avoided. However, the number of observations in the test data is also important: the test accuracy cannot be evaluated accurately if the number of observations is too small (less than 20%). The multiple output strategy (5 days are used as outputs), which is one of the multistep forecasting techniques, is applied in the experiments.The wind speed is influenced by several meteorological variables. Consequently, seven forecasting models consider multidimensional meteorological factors, such as the temperature, pressure, relative humidity, precipitable water, and wind direction, as dependent variables in multivariate models. The use of multivariable inputs increases the difficulty of determining the model parameters.Four performance criteria, namely, MAPE, RMSE, MSE, and TIC, for each forecasting model in the four seasons at the six sites, are listed in Tables3–5. Figures 5–7 present the forecasting results of the seven models for observation sites 1–6, which successfully illustrate that the proposed hybrid model RF-GSFS-WOA yields closer results to the real wind speed data than any other forecasting model in all instances. Furthermore, the developing hybrid framework realizes the best MAPE, RMSE, MSE, and TIC values. For hourly wind speed forecasting with multidimensional meteorological variables, the forecasting performance is described in detail in Tables 3–5 and Figures 5–8, from which the following findings are obtained.(I)
From Table3 and Figure 5, it is found that Elman-based models, namely, Elman, Elman-LASSO, and Elman-SCAD, outperform the BP-based models, namely, BP, BP-LASSO, and BP-SCAD, and that Elman and BP yield the worst forecasting results. For instance, at Site 1, the developed hybrid model, namely, RF-GSFS-WOA, outperforms the other six forecasting models. The RF-GSFS-WOA model realizes the smallest MAPE of 2.44% in summer, and similar results are obtained in terms of the MAE, RMSE, and TIC criteria.(II)
According to Table4 and Figure 6, the proposed RF-GSFS-WOA still outperforms other forecasting models by realizing lower values of MAE, RMSE, MAPE, and TIC. The forecasting models that were established based on Elman substantially outperform the models that were constructed based on BP. Elman-SCAD outperforms BP-SCAD. For example, the MAPE of Elman-SCAD is 7.61, which is critically lower than that of BP-SCAD, which is 10.37 in spring at Site 3. In autumn at Site 4, the RMSE of Elman-SCAD is 0.27, which is approximately 34.15% lower than that of BP-SCAD. Furthermore, similar to the results from Table 3 and Figure 5, the SCAD-related forecasting models outperform the LASSO-related forecasting models. This is mainly because SCAD can better handle the collinearity of variables than LASSO. Similar phenomena are observed in Table 5 and Figure 7. RF-GSFS-WOA realizes remarkable forecasting performances.(III)
According to Figure8, which presents the average MAPE and RMSE values of the compared forecasting methods in the four seasons, RF-GSFS-WOA realizes the lowest MAPE values over all seasons, which are 4.04%, 3.73%, 3.02%, 3.08%, 2.81%, and 4.13% at the six sites. Elman-SCAD still outperforms Elman-LASSO and BP-SCAD. The worst value is obtained by BP, which is approximately 20%. BP-LASSO and BP-SCAD slightly outperform BP, which does not perform as well as Elman.Table 3
Forecasting performances of the proposed model and other models in four seasons at Site 1 and Site 2.
Site
Seasons
Criteria
Forecasting models
Elman
Elman-LASSO
Elman-SCAD
BP
BP-LASSO
BP-SCAD
RF-GSFS-WOA
Site 1
Spring
MAE (m/s)
0.38
0.35
0.31
0.47
0.40
0.35
0.08
RMSE (m/s)
0.49
0.43
0.37
0.59
0.49
0.43
0.10
MAPE (%)
17.27
13.32
12.27
20.71
17.15
13.78
3.03
TIC (%)
6.35
5.51
4.79
7.67
6.38
5.60
1.29
Summer
MAE (m/s)
0.38
0.33
0.29
0.49
0.39
0.36
0.08
RMSE (m/s)
0.49
0.42
0.36
0.59
0.49
0.44
0.10
MAPE (%)
12.76
10.30
8.99
15.38
11.86
11.85
2.44
TIC (%)
6.11
5.24
4.49
7.25
6.15
5.52
1.24
Autumn
MAE (m/s)
0.39
0.28
0.22
0.44
0.40
0.34
0.08
RMSE (m/s)
0.48
0.35
0.27
0.56
0.49
0.43
0.10
MAPE (%)
26.46
19.85
13.76
30.44
27.60
21.55
5.90
TIC (%)
8.57
6.31
4.90
10.17
8.75
7.72
1.81
Winter
MAE (m/s)
0.38
0.33
0.29
0.41
0.40
0.31
0.08
RMSE (m/s)
0.49
0.42
0.36
0.51
0.51
0.40
0.10
MAPE (%)
18.39
17.22
15.68
22.21
18.24
16.64
4.78
TIC (%)
7.15
6.13
5.27
7.43
7.38
5.99
1.46
Site
Seasons
Criteria
Forecasting models
Elman
Elman-LASSO
Elman-SCAD
BP
BP-LASSO
BP-SCAD
RF-GSFS-WOA
Site 2
Spring
MAE (m/s)
0.38
0.33
0.27
0.52
0.39
0.37
0.11
RMSE (m/s)
0.49
0.42
0.34
0.64
0.48
0.48
0.14
MAPE (%)
14.34
14.06
10.36
22.23
15.38
13.12
4.87
TIC (%)
6.93
5.95
4.77
9.07
6.81
6.81
1.98
Summer
MAE (m/s)
0.38
0.33
0.29
0.49
0.39
0.36
0.08
RMSE (m/s)
0.49
0.42
0.36
0.59
0.49
0.44
0.10
MAPE (%)
19.23
17.71
14.15
25.68
20.56
19.24
4.16
TIC (%)
9.10
7.81
6.69
10.78
9.16
8.24
1.86
Autumn
MAE (m/s)
0.39
0.29
0.22
0.42
0.37
0.37
0.08
RMSE (m/s)
0.48
0.34
0.28
0.52
0.46
0.48
0.10
MAPE (%)
10.92
8.28
6.35
11.88
10.21
10.21
2.17
TIC (%)
5.08
3.67
2.97
5.60
4.94
5.15
1.07
Winter
MAE (m/s)
0.38
0.33
0.29
0.43
0.38
0.30
0.08
RMSE (m/s)
0.49
0.42
0.36
0.52
0.47
0.37
0.10
MAPE (%)
8.14
6.73
6.30
9.02
7.90
6.28
1.63
TIC (%)
4.44
3.80
3.27
4.68
4.29
3.39
0.90Table 4
Forecasting performances of the proposed model and other models in four seasons at Site 3 and Site 4.
Site
Seasons
Criteria
Forecasting models
Elman
Elman-LASSO
Elman-SCAD
BP
BP-LASSO
BP-SCAD
RF-GSFS-WOA
Site 3
Spring
MAE (m/s)
0.38
0.33
0.29
0.53
0.39
0.36
0.08
RMSE (m/s)
0.49
0.42
0.36
0.64
0.49
0.44
0.10
MAPE (%)
10.62
9.03
7.61
13.97
10.82
10.37
2.09
TIC (%)
6.00
5.13
4.41
7.71
6.04
5.43
1.22
Summer
MAE (m/s)
0.38
0.33
0.29
0.49
0.39
0.36
0.08
RMSE (m/s)
0.49
0.42
0.36
0.59
0.49
0.44
0.10
MAPE (%)
12.37
10.59
8.83
15.97
12.75
11.65
2.54
TIC (%)
6.51
5.58
4.78
7.74
6.58
5.88
1.32
Autumn
MAE (m/s)
0.38
0.26
0.22
0.43
0.38
0.36
0.08
RMSE (m/s)
0.49
0.33
0.27
0.52
0.47
0.44
0.10
MAPE (%)
16.01
12.03
9.66
19.39
17.64
15.32
3.44
TIC (%)
7.73
5.14
4.21
8.08
7.43
7.01
1.57
Winter
MAE (m/s)
0.38
0.29
0.28
0.42
0.39
0.31
0.08
RMSE (m/s)
0.49
0.36
0.33
0.52
0.51
0.39
0.10
MAPE (%)
16.97
13.78
14.81
18.16
17.39
16.42
3.99
TIC (%)
7.26
5.32
4.92
7.70
7.58
5.73
1.48
Site
Seasons
Criteria
Forecasting models
Elman
Elman-LASSO
Elman-SCAD
BP
BP-LASSO
BP-SCAD
RF-GSFS-WOA
Site 4
Spring
MAE (m/s)
0.29
0.26
0.22
0.37
0.29
0.23
0.08
RMSE (m/s)
0.37
0.31
0.28
0.46
0.36
0.29
0.10
MAPE (%)
14.17
12.38
9.78
16.45
14.09
9.97
3.51
TIC (%)
6.16
5.15
4.55
7.51
5.84
4.71
1.64
Summer
MAE (m/s)
0.38
0.33
0.25
0.42
0.35
0.29
0.08
RMSE (m/s)
0.49
0.39
0.32
0.52
0.43
0.37
0.10
MAPE (%)
15.78
13.77
10.70
18.66
14.62
11.82
3.21
TIC (%)
8.24
6.56
5.35
8.74
7.16
6.22
1.68
Autumn
MAE (m/s)
0.38
0.26
0.22
0.47
0.38
0.33
0.08
RMSE (m/s)
0.49
0.33
0.27
0.59
0.47
0.41
0.10
MAPE (%)
13.58
7.87
7.75
17.97
14.40
10.09
2.63
TIC (%)
5.64
3.75
3.06
6.78
5.44
4.68
1.14
Winter
MAE (m/s)
0.38
0.33
0.29
0.43
0.38
0.30
0.08
RMSE (m/s)
0.49
0.42
0.36
0.52
0.47
0.37
0.10
MAPE (%)
16.10
13.29
11.24
17.53
14.17
11.19
2.96
TIC (%)
6.06
5.19
4.46
6.38
5.86
4.63
1.23Table 5
Forecasting performances of the proposed model and other models in four seasons at Site 5 and Site 6.
Site
Seasons
Criteria
Forecasting models
Elman
Elman-LASSO
Elman-SCAD
BP
BP-LASSO
BP-SCAD
RF-GSFS-WOA
Site 5
Spring
MAE (m/s)
0.38
0.33
0.29
0.53
0.39
0.36
0.08
RMSE (m/s)
0.49
0.42
0.36
0.64
0.49
0.44
0.10
MAPE (%)
10.72
9.69
8.43
15.08
11.16
10.13
2.46
TIC (%)
5.39
4.60
3.96
6.93
5.43
4.89
1.10
Summer
MAE (m/s)
0.38
0.33
0.29
0.49
0.39
0.36
0.08
RMSE (m/s)
0.49
0.42
0.36
0.59
0.49
0.44
0.10
MAPE (%)
13.01
10.32
8.42
15.75
12.30
11.49
2.39
TIC (%)
5.94
5.09
4.36
7.07
5.99
5.37
1.21
Autumn
MAE (m/s)
0.38
0.26
0.22
0.43
0.38
0.36
0.08
RMSE (m/s)
0.49
0.33
0.27
0.52
0.48
0.44
0.10
MAPE (%)
12.59
8.55
7.03
14.66
12.60
11.98
2.58
TIC (%)
6.37
4.25
3.47
6.70
6.15
5.77
1.30
Winter
MAE (m/s)
0.38
0.31
0.25
0.42
0.38
0.29
0.08
RMSE (m/s)
0.49
0.38
0.29
0.52
0.47
0.36
0.10
MAPE (%)
16.24
14.03
11.94
18.39
16.97
12.04
3.82
TIC (%)
6.45
4.98
3.87
6.87
6.20
4.81
1.31
Site
Seasons
Criteria
Forecasting models
Elman
Elman-LASSO
Elman-SCAD
BP
BP-LASSO
BP-SCAD
RF-GSFS-WOA
Site 6
Spring
MAE (m/s)
0.29
0.26
0.22
0.41
0.33
0.29
0.08
RMSE (m/s)
0.37
0.31
0.28
0.53
0.41
0.36
0.10
MAPE (%)
13.95
11.40
9.65
18.60
14.59
12.41
3.45
TIC (%)
6.19
5.16
4.56
8.79
6.70
5.84
1.65
Summer
MAE (m/s)
0.26
0.22
0.18
0.30
0.23
0.18
0.08
RMSE (m/s)
0.31
0.28
0.23
0.37
0.28
0.23
0.10
MAPE (%)
16.73
12.62
10.06
18.37
13.56
11.60
4.99
TIC (%)
5.60
4.97
4.14
6.57
4.98
4.17
1.79
Autumn
MAE (m/s)
0.38
0.27
0.21
0.45
0.32
0.28
0.08
RMSE (m/s)
0.49
0.34
0.26
0.54
0.39
0.37
0.10
MAPE (%)
16.37
11.46
9.36
19.20
13.45
11.73
3.26
TIC (%)
7.18
5.03
3.82
7.76
5.77
5.44
1.45
Winter
MAE (m/s)
0.27
0.20
0.15
0.26
0.25
0.21
0.08
RMSE (m/s)
0.33
0.27
0.19
0.34
0.31
0.27
0.10
MAPE (%)
14.38
11.93
9.19
15.63
14.07
13.06
4.82
TIC (%)
6.27
5.06
3.61
6.37
5.88
5.02
1.89Figure 5
The forecasting results of seven models at Site 1 and Site 2.Figure 6
Forecasting results of seven models at Site 3 and Site 4.Figure 7
Forecasting results of seven models at Site 5 and Site 6.Figure 8
Average values of MAPE and RMSE for four seasons at six sites.
(a)
(b)Therefore, from the discussions above, the proposed hybrid forecasting model, namely, RF-GSFS-WOA, performs well in the wind speed forecasting task by combining the advantages of machine learning module RF, feature selection approach GSFS, and metaheuristic method WOA.


### 4.4. Statistical Test

To further evaluate the performance of the proposed forecasting model, a nonparametric statistical test, namely, the Friedman test, is conducted. The forecasting performances of the compared methods differ among the six sites [42]. The following hypothesis, namely, the null hypothesis H0, along with the alternative hypothesis H1, is considered in the Friedman test:(25)H0:L1=L2=…=Lk,where Lj=1/B∑iaij, with aij representing the j-th rank of k algorithms on a dataset of size B. The Friedman statistic is defined as follows [43]:(26)χF2=12Bkk+1∑jLj2−kk+124.This statistic follows chi-square distributionχF2k−1. An F-distribution statistic, namely, Fk−1,k−1B−1, is calculated as follows:(27)FF=B−1χF2Bk−1−χF2.Suppose that there are significant differences among the comparing forecasting methods. Then, the null hypothesisH0 will be rejected. If this occurs, critical differences (CDs) will be applied for further comparisons, which is known as a post hoc test. A Bonferroni-Dunn test is used with CD values q0.1kk+1/6B, where q0.1 is determined as 2.291 in [44]. The Friedman test is described in detail in Table 6. All the FF values are much larger than the critical value F6,30 = 2.42. This suggests that the null hypothesis H0 should be rejected and a post hoc test is needed. Using the above-described method, the CD value is calculated as q0.1kk+1/6B=2.291×7×7+1/6×6=2.8574, which is used as a benchmark for comparison. The proposed RF-GSFS-WOA significantly outperforms BP-LASSO, BP, and Elman since the differences between the average ranks of these methods exceed 2.85. For instance, the difference between RF-GSFS-WOA and BP-LASSO is 4.67 for MAE, which is larger than 2.85. Compared to other methods, the proposed method does not show strong advantages in forecasting. However, based on the above analysis, the proposed method delivers remarkable forecasting results.Table 6
Results of Friedman test for all forecasting methods.
Models
MAE
Site 1
Site 2
Site 3
Site 4
Site 5
Site 6
Average
χ
F
2
F
F
Elman
5
5
5
6
5
6
5.33
35.07
188.84
Elman-LASSO
3
3
3
4
3
3
3.17
Elman-SCAD
2
2
2
2
2
2
2.00
BP
7
7
7
7
7
7
7.00
BP-LASSO
6
6
6
5
6
5
5.67
BP-SCAD
4
4
4
3
4
4
3.83
RF-GSFS-WOA
1
1
1
1
1
1
1.00
Models
RMSE
Site 1
Site 2
Site 3
Site 4
Site 5
Site 6
Average
χ
F
2
F
F
Elman
5
6
5
6
6
6
5.67
35.07
188.84
Elman-LASSO
3
3
3
4
3
3
3.17
Elman-SCAD
2
2
2
2
2
2
2.00
BP
7
7
7
7
7
7
7.00
BP-LASSO
6
5
6
5
5
5
5.33
BP-SCAD
4
4
4
3
4
4
3.83
RF-GSFS-WOA
1
1
1
1
1
1
1.00
Models
MAPE
Site 1
Site 2
Site 3
Site 4
Site 5
Site 6
Average
χ
F
2
F
F
Elman
5
6
5
6
6
6
5.50
35
175
Elman-LASSO
3
3
3
4
3
3
3.17
Elman-SCAD
2
2
2
2
2
2
2.00
BP
7
7
7
7
7
7
7.00
BP-LASSO
6
5
6
5
5
5
5.50
BP-SCAD
4
4
4
3
4
4
3.83
RF-GSFS-WOA
1
1
1
1
1
1
1.00
Models
TIC
Site 1
Site 2
Site 3
Site 4
Site 5
Site 6
Average
χ
F
2
F
F
Elman
5
6
5
6
6
6
5.67
35.07
188.84
Elman-LASSO
3
3
3
4
3
3
3.17
Elman-SCAD
2
2
2
2
2
2
2.00
BP
7
7
7
7
7
7
7.00
BP-LASSO
6
5
6
5
5
5
5.33
BP-SCAD
4
4
4
3
4
4
3.83
RF-GSFS-WOA
1
1
1
1
1
1
1.00


## 4.1. Data Collection

Qinghai province is located in the northeast of the Qinghai-Tibetan Plateau, and most of the area is rich in wind energy resources; hence, it is suitable for the installation of large-scale wind farm. Thus, it is meaningful to investigate sites in this region. The hourly data were collected between 0 : 00 and 24 : 00 from 1/1/2014 to 12/31/2014 at six sites in the Qinghai region. The dataset that is used in this study was collected by National Renewable Energy Laboratory (NREL) and is available athttps://www.nrel.gov/gis/tools.html. To simplify the description, the six sites are denoted by Site 1, Site 2, Site 3, Site 4, Site 5, and Site 6. The longitudes and latitudes of these locations are specified in Figure 4. Five meteorological factors, namely, the temperature (T), pressure (P), relative humidity (RH), wind precipitable water (PW), and wind direction (WD), are considered in the forecasting of the wind speed for the Qinghai region, which facilitates the design of wind power plants to efficiently generate electricity for the local population [41]. There are four seasons in each year: spring (Mar., Apr., and May), summer (Jun., Jul., and Aug.), autumn (Sep., Oct., and Nov.), and winter (Dec., Jan., and Feb.). The data are divided into training data and forecasting data as follows: 20 days in each season are randomly chosen as the training data (20 days × 24 hours/days = 480) for building the model, and the forecasting data include 5 days (5 days × 24 hours/day = 120) that are randomly chosen from the remaining days in each season. The performances of various forecasting models are evaluated in the four seasons. The training sets and forecasting sets for the six sites are described in Table 1. The forecasting sets account for approximately 20% of the total data, which is reasonable.Figure 4
Distribution and location information of the observation sites.Table 1
Training sets and forecasting sets for hourly wind speed data at six sites.
Season
Training set (size: 20 × 24 = 480)
Forecasting set (size: 5 × 24 = 120)
Relative meteorological factors
Spring (Mar., Apr., and May)
Mar. 9th, Mar. 12th, Mar. 19th, and Mar. 21st
May 16th and May 20th
Temperature, pressure, precipitable water, and wind direction
Mar. 26th, Mar. 29th, Apr. 5th, and Apr. 6th
May 25th, May 28th, and May 29th
Apr. 13th, Apr. 15th, Apr. 16th, and Apr. 21st
Apr. 22nd, Apr. 23th, Apr. 24th, and Apr. 25th
Apr. 28th, May 1st, May 8th, and May 9th
Summer (Jun., Jul., and Aug.)
Jun. 1st, Jun. 3rd, Jun. 4th, and Jun. 6th
Aug. 12th and Aug. 20th
Jun. 13th, Jun. 16th, Jun. 22nd, and Jun. 28th
Aug. 22nd, Aug. 24th, and Aug. 25th
Jun. 30th, Jul. 2nd, Jul. 7th, and Jul. 13th
Jul. 15th, Jul. 18th, Jul. 22nd, and Jul. 28th
Jul. 31st, Aug. 1st, Aug. 3rd, and Aug. 11th
Autumn (Sep., Oct., and Nov.)
Sep. 3rd, Sep. 6th, Sep. 8th, and Sep. 11th
Nov. 18th and Nov. 23rd
Sep. 14th, Sep. 15th, Sep. 21st, and Sep. 24th
Nov. 24th, Nov. 27th, and Nov. 29th
Sep. 30th, Oct. 4th, Oct. 6th, and Oct. 13th
Oct. 18th, Oct. 31st, Nov. 2nd, and Nov. 3rd
Nov. 4th, Nov. 8th, Nov. 12th, and Nov. 14th
Winter (Dec., Jan., and Feb.)
Dec. 2nd, Dec. 3rd, Dec. 4th, and Dec 8th
Feb. 7th and Feb. 8th
Dec. 9th, Dec. 13th, Dec. 14th, and Dec. 18th
Feb. 15th, Feb. 24th, and Feb. 26th
Dec. 20th, Dec. 25th, Dec. 26th, and Dec. 29th
Dec. 30th, Jan. 5th, Jan. 8th, and Jan 9th
Jan. 16th, Jan. 18th, Jan. 28th, and Feb. 5th


## 4.2. Forecasting Assessment Criteria

To evaluate the results of the forecasting models quantitatively, four assessment measures, namely, the mean absolute error (MAE), root mean square error (RMSE), mean absolute percent error (MAPE), and Theil inequality coefficient (TIC), are applied. All are calculated between the actual data and the forecasted values, and the smaller their values, the better the forecasting results. These criteria are expressed as follows:(24)MAE=1N∑i=1Nyi−y^i,RMSE=1N∑i=1Nyi−y^i2,MAPE=1N∑i=1Nyi−y^iyi×100%,TIC=1/N∑i=1Nyi−y^i21/N∑i=1Ny^i2+1/N∑i=1Nyi2×100%,where N is the size of the forecasting sample and yi and y^i are the actual value and the forecasted value, respectively, for time period i.


## 4.3. Experiment Results and Corresponding Analysis

The experimental settings are described in detail in Table2. The parameters in the complementary ensemble empirical mode decomposition with adaptive noise (CEEMDAN) are set as follows: the noise standard deviation (Nstd) ε0 is set as 0.2, the number of realizations (Nr) is set as 500, and the maximum allowed number of sifting iterations (MaxIter) is set as 500. The regularization parameters in LASSO and SCAD are chosen via 10-fold cross-validation. In the experiments on the neural networks, the maximum number of iterations (MaxIterNN) is set as 2000 for both BP and Elman. Each activation function from the hidden layer to the output layer is set as the Tansig function (Func). The number of hidden neurons is crucial since it mainly determines the network complexity. To determine the optimal number of hidden neurons, values are selected from the following set: 5,8,10,15,20,25. Furthermore, Elman is trained using backpropagation through time, and the corresponding weights are selected based on a uniform distribution. Based on the experimental results, the gradient descent with momentum is set as 0.82, and the adaptive learning rate is set as 0.02. In GSFS, we set q as 30; hence, 30 variables will be chosen after the GSFS procedure, and, thus, there are 1200 variables in the quadratic model. The number of trees in RF is set as 3000, and the maximum number of iterations (MaxIterRF) is 500. RF is implemented using MATLAB function TreeBagger. The fitness function of WOA is selected as “F1,” which is the square error loss function, and the range is (0, 1). The weights are determined after normalization such that the sum of all the weights is 1.Table 2
Important parameters of seven forecasting models.
Nstd
Nr
MaxIter
λ
MaxIterNN
Func
h
q
MaxIterRF
Elman
0.2
500
500
–
2000
Tansig
10
–
–
Elman-LASSO
0.2
500
500
0.875
2000
Tansig
10
–
–
Elman-SCAD
0.2
500
500
0.625
2000
Tansig
10
–
–
BP
0.2
500
500
–
2000
Tansig
10
–
–
BP-LASSO
0.2
500
500
0.875
2000
Tansig
10
–
–
BP-SCAD
0.2
500
500
0.625
2000
Tansig
10
–
–
RF-GSFS-WOA
0.2
500
500
–
–
Tansig
-
30
500To evaluate the forecasting performance of the proposed hybrid model for forecasting hourly wind speed, experiments for six sites and four seasons are designed in this paper. Four types of assessment criteria are used in all experiments to evaluate the performance. The proposed hybrid forecasting model is RF-GSFS-WOA, and six forecasting models, namely, an Elman neural network, Elman-LASSO, Elman-SCAD, a BP neural network (BP), BP-LASSO, and BP-SCAD, are applied as comparison models to evaluate the performance of the developed hybrid model. Wind speed data exhibit volatility or instability due to changes in the weather; hence, this study conducts a preprocessing step before building all the models; namely, the complementary ensemble empirical mode decomposition with adaptive noise(CEEMDAN) is used to eliminate the noise and smooth the wind speed time series data before the developed hybrid model and six comparison models are established for hourly wind speed forecasting in four seasons at six sites.In each season, 20 days of hourly wind speed time series are selected randomly as the training phase for the construction of the models, and 5 days of wind speed time series are selected randomly as the forecasting phase. Thus, 80% of the data will be used as training data; the remaining data (20%) will be used as test data. The reason for this splitting scheme is that the model can be trained sufficiently with more observations (80% of the data) and overfitting can be avoided. However, the number of observations in the test data is also important: the test accuracy cannot be evaluated accurately if the number of observations is too small (less than 20%). The multiple output strategy (5 days are used as outputs), which is one of the multistep forecasting techniques, is applied in the experiments.The wind speed is influenced by several meteorological variables. Consequently, seven forecasting models consider multidimensional meteorological factors, such as the temperature, pressure, relative humidity, precipitable water, and wind direction, as dependent variables in multivariate models. The use of multivariable inputs increases the difficulty of determining the model parameters.Four performance criteria, namely, MAPE, RMSE, MSE, and TIC, for each forecasting model in the four seasons at the six sites, are listed in Tables3–5. Figures 5–7 present the forecasting results of the seven models for observation sites 1–6, which successfully illustrate that the proposed hybrid model RF-GSFS-WOA yields closer results to the real wind speed data than any other forecasting model in all instances. Furthermore, the developing hybrid framework realizes the best MAPE, RMSE, MSE, and TIC values. For hourly wind speed forecasting with multidimensional meteorological variables, the forecasting performance is described in detail in Tables 3–5 and Figures 5–8, from which the following findings are obtained.(I)
From Table3 and Figure 5, it is found that Elman-based models, namely, Elman, Elman-LASSO, and Elman-SCAD, outperform the BP-based models, namely, BP, BP-LASSO, and BP-SCAD, and that Elman and BP yield the worst forecasting results. For instance, at Site 1, the developed hybrid model, namely, RF-GSFS-WOA, outperforms the other six forecasting models. The RF-GSFS-WOA model realizes the smallest MAPE of 2.44% in summer, and similar results are obtained in terms of the MAE, RMSE, and TIC criteria.(II)
According to Table4 and Figure 6, the proposed RF-GSFS-WOA still outperforms other forecasting models by realizing lower values of MAE, RMSE, MAPE, and TIC. The forecasting models that were established based on Elman substantially outperform the models that were constructed based on BP. Elman-SCAD outperforms BP-SCAD. For example, the MAPE of Elman-SCAD is 7.61, which is critically lower than that of BP-SCAD, which is 10.37 in spring at Site 3. In autumn at Site 4, the RMSE of Elman-SCAD is 0.27, which is approximately 34.15% lower than that of BP-SCAD. Furthermore, similar to the results from Table 3 and Figure 5, the SCAD-related forecasting models outperform the LASSO-related forecasting models. This is mainly because SCAD can better handle the collinearity of variables than LASSO. Similar phenomena are observed in Table 5 and Figure 7. RF-GSFS-WOA realizes remarkable forecasting performances.(III)
According to Figure8, which presents the average MAPE and RMSE values of the compared forecasting methods in the four seasons, RF-GSFS-WOA realizes the lowest MAPE values over all seasons, which are 4.04%, 3.73%, 3.02%, 3.08%, 2.81%, and 4.13% at the six sites. Elman-SCAD still outperforms Elman-LASSO and BP-SCAD. The worst value is obtained by BP, which is approximately 20%. BP-LASSO and BP-SCAD slightly outperform BP, which does not perform as well as Elman.Table 3
Forecasting performances of the proposed model and other models in four seasons at Site 1 and Site 2.
Site
Seasons
Criteria
Forecasting models
Elman
Elman-LASSO
Elman-SCAD
BP
BP-LASSO
BP-SCAD
RF-GSFS-WOA
Site 1
Spring
MAE (m/s)
0.38
0.35
0.31
0.47
0.40
0.35
0.08
RMSE (m/s)
0.49
0.43
0.37
0.59
0.49
0.43
0.10
MAPE (%)
17.27
13.32
12.27
20.71
17.15
13.78
3.03
TIC (%)
6.35
5.51
4.79
7.67
6.38
5.60
1.29
Summer
MAE (m/s)
0.38
0.33
0.29
0.49
0.39
0.36
0.08
RMSE (m/s)
0.49
0.42
0.36
0.59
0.49
0.44
0.10
MAPE (%)
12.76
10.30
8.99
15.38
11.86
11.85
2.44
TIC (%)
6.11
5.24
4.49
7.25
6.15
5.52
1.24
Autumn
MAE (m/s)
0.39
0.28
0.22
0.44
0.40
0.34
0.08
RMSE (m/s)
0.48
0.35
0.27
0.56
0.49
0.43
0.10
MAPE (%)
26.46
19.85
13.76
30.44
27.60
21.55
5.90
TIC (%)
8.57
6.31
4.90
10.17
8.75
7.72
1.81
Winter
MAE (m/s)
0.38
0.33
0.29
0.41
0.40
0.31
0.08
RMSE (m/s)
0.49
0.42
0.36
0.51
0.51
0.40
0.10
MAPE (%)
18.39
17.22
15.68
22.21
18.24
16.64
4.78
TIC (%)
7.15
6.13
5.27
7.43
7.38
5.99
1.46
Site
Seasons
Criteria
Forecasting models
Elman
Elman-LASSO
Elman-SCAD
BP
BP-LASSO
BP-SCAD
RF-GSFS-WOA
Site 2
Spring
MAE (m/s)
0.38
0.33
0.27
0.52
0.39
0.37
0.11
RMSE (m/s)
0.49
0.42
0.34
0.64
0.48
0.48
0.14
MAPE (%)
14.34
14.06
10.36
22.23
15.38
13.12
4.87
TIC (%)
6.93
5.95
4.77
9.07
6.81
6.81
1.98
Summer
MAE (m/s)
0.38
0.33
0.29
0.49
0.39
0.36
0.08
RMSE (m/s)
0.49
0.42
0.36
0.59
0.49
0.44
0.10
MAPE (%)
19.23
17.71
14.15
25.68
20.56
19.24
4.16
TIC (%)
9.10
7.81
6.69
10.78
9.16
8.24
1.86
Autumn
MAE (m/s)
0.39
0.29
0.22
0.42
0.37
0.37
0.08
RMSE (m/s)
0.48
0.34
0.28
0.52
0.46
0.48
0.10
MAPE (%)
10.92
8.28
6.35
11.88
10.21
10.21
2.17
TIC (%)
5.08
3.67
2.97
5.60
4.94
5.15
1.07
Winter
MAE (m/s)
0.38
0.33
0.29
0.43
0.38
0.30
0.08
RMSE (m/s)
0.49
0.42
0.36
0.52
0.47
0.37
0.10
MAPE (%)
8.14
6.73
6.30
9.02
7.90
6.28
1.63
TIC (%)
4.44
3.80
3.27
4.68
4.29
3.39
0.90Table 4
Forecasting performances of the proposed model and other models in four seasons at Site 3 and Site 4.
Site
Seasons
Criteria
Forecasting models
Elman
Elman-LASSO
Elman-SCAD
BP
BP-LASSO
BP-SCAD
RF-GSFS-WOA
Site 3
Spring
MAE (m/s)
0.38
0.33
0.29
0.53
0.39
0.36
0.08
RMSE (m/s)
0.49
0.42
0.36
0.64
0.49
0.44
0.10
MAPE (%)
10.62
9.03
7.61
13.97
10.82
10.37
2.09
TIC (%)
6.00
5.13
4.41
7.71
6.04
5.43
1.22
Summer
MAE (m/s)
0.38
0.33
0.29
0.49
0.39
0.36
0.08
RMSE (m/s)
0.49
0.42
0.36
0.59
0.49
0.44
0.10
MAPE (%)
12.37
10.59
8.83
15.97
12.75
11.65
2.54
TIC (%)
6.51
5.58
4.78
7.74
6.58
5.88
1.32
Autumn
MAE (m/s)
0.38
0.26
0.22
0.43
0.38
0.36
0.08
RMSE (m/s)
0.49
0.33
0.27
0.52
0.47
0.44
0.10
MAPE (%)
16.01
12.03
9.66
19.39
17.64
15.32
3.44
TIC (%)
7.73
5.14
4.21
8.08
7.43
7.01
1.57
Winter
MAE (m/s)
0.38
0.29
0.28
0.42
0.39
0.31
0.08
RMSE (m/s)
0.49
0.36
0.33
0.52
0.51
0.39
0.10
MAPE (%)
16.97
13.78
14.81
18.16
17.39
16.42
3.99
TIC (%)
7.26
5.32
4.92
7.70
7.58
5.73
1.48
Site
Seasons
Criteria
Forecasting models
Elman
Elman-LASSO
Elman-SCAD
BP
BP-LASSO
BP-SCAD
RF-GSFS-WOA
Site 4
Spring
MAE (m/s)
0.29
0.26
0.22
0.37
0.29
0.23
0.08
RMSE (m/s)
0.37
0.31
0.28
0.46
0.36
0.29
0.10
MAPE (%)
14.17
12.38
9.78
16.45
14.09
9.97
3.51
TIC (%)
6.16
5.15
4.55
7.51
5.84
4.71
1.64
Summer
MAE (m/s)
0.38
0.33
0.25
0.42
0.35
0.29
0.08
RMSE (m/s)
0.49
0.39
0.32
0.52
0.43
0.37
0.10
MAPE (%)
15.78
13.77
10.70
18.66
14.62
11.82
3.21
TIC (%)
8.24
6.56
5.35
8.74
7.16
6.22
1.68
Autumn
MAE (m/s)
0.38
0.26
0.22
0.47
0.38
0.33
0.08
RMSE (m/s)
0.49
0.33
0.27
0.59
0.47
0.41
0.10
MAPE (%)
13.58
7.87
7.75
17.97
14.40
10.09
2.63
TIC (%)
5.64
3.75
3.06
6.78
5.44
4.68
1.14
Winter
MAE (m/s)
0.38
0.33
0.29
0.43
0.38
0.30
0.08
RMSE (m/s)
0.49
0.42
0.36
0.52
0.47
0.37
0.10
MAPE (%)
16.10
13.29
11.24
17.53
14.17
11.19
2.96
TIC (%)
6.06
5.19
4.46
6.38
5.86
4.63
1.23Table 5
Forecasting performances of the proposed model and other models in four seasons at Site 5 and Site 6.
Site
Seasons
Criteria
Forecasting models
Elman
Elman-LASSO
Elman-SCAD
BP
BP-LASSO
BP-SCAD
RF-GSFS-WOA
Site 5
Spring
MAE (m/s)
0.38
0.33
0.29
0.53
0.39
0.36
0.08
RMSE (m/s)
0.49
0.42
0.36
0.64
0.49
0.44
0.10
MAPE (%)
10.72
9.69
8.43
15.08
11.16
10.13
2.46
TIC (%)
5.39
4.60
3.96
6.93
5.43
4.89
1.10
Summer
MAE (m/s)
0.38
0.33
0.29
0.49
0.39
0.36
0.08
RMSE (m/s)
0.49
0.42
0.36
0.59
0.49
0.44
0.10
MAPE (%)
13.01
10.32
8.42
15.75
12.30
11.49
2.39
TIC (%)
5.94
5.09
4.36
7.07
5.99
5.37
1.21
Autumn
MAE (m/s)
0.38
0.26
0.22
0.43
0.38
0.36
0.08
RMSE (m/s)
0.49
0.33
0.27
0.52
0.48
0.44
0.10
MAPE (%)
12.59
8.55
7.03
14.66
12.60
11.98
2.58
TIC (%)
6.37
4.25
3.47
6.70
6.15
5.77
1.30
Winter
MAE (m/s)
0.38
0.31
0.25
0.42
0.38
0.29
0.08
RMSE (m/s)
0.49
0.38
0.29
0.52
0.47
0.36
0.10
MAPE (%)
16.24
14.03
11.94
18.39
16.97
12.04
3.82
TIC (%)
6.45
4.98
3.87
6.87
6.20
4.81
1.31
Site
Seasons
Criteria
Forecasting models
Elman
Elman-LASSO
Elman-SCAD
BP
BP-LASSO
BP-SCAD
RF-GSFS-WOA
Site 6
Spring
MAE (m/s)
0.29
0.26
0.22
0.41
0.33
0.29
0.08
RMSE (m/s)
0.37
0.31
0.28
0.53
0.41
0.36
0.10
MAPE (%)
13.95
11.40
9.65
18.60
14.59
12.41
3.45
TIC (%)
6.19
5.16
4.56
8.79
6.70
5.84
1.65
Summer
MAE (m/s)
0.26
0.22
0.18
0.30
0.23
0.18
0.08
RMSE (m/s)
0.31
0.28
0.23
0.37
0.28
0.23
0.10
MAPE (%)
16.73
12.62
10.06
18.37
13.56
11.60
4.99
TIC (%)
5.60
4.97
4.14
6.57
4.98
4.17
1.79
Autumn
MAE (m/s)
0.38
0.27
0.21
0.45
0.32
0.28
0.08
RMSE (m/s)
0.49
0.34
0.26
0.54
0.39
0.37
0.10
MAPE (%)
16.37
11.46
9.36
19.20
13.45
11.73
3.26
TIC (%)
7.18
5.03
3.82
7.76
5.77
5.44
1.45
Winter
MAE (m/s)
0.27
0.20
0.15
0.26
0.25
0.21
0.08
RMSE (m/s)
0.33
0.27
0.19
0.34
0.31
0.27
0.10
MAPE (%)
14.38
11.93
9.19
15.63
14.07
13.06
4.82
TIC (%)
6.27
5.06
3.61
6.37
5.88
5.02
1.89Figure 5
The forecasting results of seven models at Site 1 and Site 2.Figure 6
Forecasting results of seven models at Site 3 and Site 4.Figure 7
Forecasting results of seven models at Site 5 and Site 6.Figure 8
Average values of MAPE and RMSE for four seasons at six sites.
(a)
(b)Therefore, from the discussions above, the proposed hybrid forecasting model, namely, RF-GSFS-WOA, performs well in the wind speed forecasting task by combining the advantages of machine learning module RF, feature selection approach GSFS, and metaheuristic method WOA.


## 4.4. Statistical Test

To further evaluate the performance of the proposed forecasting model, a nonparametric statistical test, namely, the Friedman test, is conducted. The forecasting performances of the compared methods differ among the six sites [42]. The following hypothesis, namely, the null hypothesis H0, along with the alternative hypothesis H1, is considered in the Friedman test:(25)H0:L1=L2=…=Lk,where Lj=1/B∑iaij, with aij representing the j-th rank of k algorithms on a dataset of size B. The Friedman statistic is defined as follows [43]:(26)χF2=12Bkk+1∑jLj2−kk+124.This statistic follows chi-square distributionχF2k−1. An F-distribution statistic, namely, Fk−1,k−1B−1, is calculated as follows:(27)FF=B−1χF2Bk−1−χF2.Suppose that there are significant differences among the comparing forecasting methods. Then, the null hypothesisH0 will be rejected. If this occurs, critical differences (CDs) will be applied for further comparisons, which is known as a post hoc test. A Bonferroni-Dunn test is used with CD values q0.1kk+1/6B, where q0.1 is determined as 2.291 in [44]. The Friedman test is described in detail in Table 6. All the FF values are much larger than the critical value F6,30 = 2.42. This suggests that the null hypothesis H0 should be rejected and a post hoc test is needed. Using the above-described method, the CD value is calculated as q0.1kk+1/6B=2.291×7×7+1/6×6=2.8574, which is used as a benchmark for comparison. The proposed RF-GSFS-WOA significantly outperforms BP-LASSO, BP, and Elman since the differences between the average ranks of these methods exceed 2.85. For instance, the difference between RF-GSFS-WOA and BP-LASSO is 4.67 for MAE, which is larger than 2.85. Compared to other methods, the proposed method does not show strong advantages in forecasting. However, based on the above analysis, the proposed method delivers remarkable forecasting results.Table 6
Results of Friedman test for all forecasting methods.
Models
MAE
Site 1
Site 2
Site 3
Site 4
Site 5
Site 6
Average
χ
F
2
F
F
Elman
5
5
5
6
5
6
5.33
35.07
188.84
Elman-LASSO
3
3
3
4
3
3
3.17
Elman-SCAD
2
2
2
2
2
2
2.00
BP
7
7
7
7
7
7
7.00
BP-LASSO
6
6
6
5
6
5
5.67
BP-SCAD
4
4
4
3
4
4
3.83
RF-GSFS-WOA
1
1
1
1
1
1
1.00
Models
RMSE
Site 1
Site 2
Site 3
Site 4
Site 5
Site 6
Average
χ
F
2
F
F
Elman
5
6
5
6
6
6
5.67
35.07
188.84
Elman-LASSO
3
3
3
4
3
3
3.17
Elman-SCAD
2
2
2
2
2
2
2.00
BP
7
7
7
7
7
7
7.00
BP-LASSO
6
5
6
5
5
5
5.33
BP-SCAD
4
4
4
3
4
4
3.83
RF-GSFS-WOA
1
1
1
1
1
1
1.00
Models
MAPE
Site 1
Site 2
Site 3
Site 4
Site 5
Site 6
Average
χ
F
2
F
F
Elman
5
6
5
6
6
6
5.50
35
175
Elman-LASSO
3
3
3
4
3
3
3.17
Elman-SCAD
2
2
2
2
2
2
2.00
BP
7
7
7
7
7
7
7.00
BP-LASSO
6
5
6
5
5
5
5.50
BP-SCAD
4
4
4
3
4
4
3.83
RF-GSFS-WOA
1
1
1
1
1
1
1.00
Models
TIC
Site 1
Site 2
Site 3
Site 4
Site 5
Site 6
Average
χ
F
2
F
F
Elman
5
6
5
6
6
6
5.67
35.07
188.84
Elman-LASSO
3
3
3
4
3
3
3.17
Elman-SCAD
2
2
2
2
2
2
2.00
BP
7
7
7
7
7
7
7.00
BP-LASSO
6
5
6
5
5
5
5.33
BP-SCAD
4
4
4
3
4
4
3.83
RF-GSFS-WOA
1
1
1
1
1
1
1.00


## 5. Conclusions

Accurate wind speed forecasting has attracted the attention of researchers of wind energy for decades. This paper proposed and investigated a novel hybrid wind speed forecasting model that combines the advantages of CEEMDAN, GSFS, RF, and WOA. The proposed model has been applied to data from six sites in Qinghai province to forecast the wind speed. The forecasting performances are evaluated using four statistical measures: MAE, RMSE, MAPE, and TIC. The forecasting results demonstrate that the proposed hybrid model outperforms other compared forecasting models, namely, Elman, Elman-LASSO, Elman-SCAD, BP, BP-LASSO, and BP-SCAD, which are famous wind speed forecasting models and variable selection methods. Moreover, the Elman-based forecasting models deliver better results than the BP-based models, and the SCAD-based models outperform the LASSO-based models. Furthermore, by using variable selection methods such as GSFS, LASSO, and SCAD, the model performances are boosted substantially. Thus, variable selection is necessary since it can be used to select the meteorological factors that make the most important contributions to the forecasting model while discarding useless factors. To further evaluate the performance of the proposed model, a nonparametric test, namely, the Friedman test, is conducted to evaluate the forecasting results. The test results are consistent with the findings of the analysis of the results. Consequently, according to the experimental results, the proposed hybrid model forecasts the wind speed accurately and efficiently [38, 45].



---

*Source: 5396473-2020-05-06.xml*