# A Compound Structure for Wind Speed Forecasting Using MKLSSVM with Feature Selection and Parameter Optimization

Authors: Sizhou Sun; Jingqi Fu; Feng Zhu; Nan Xiong
Journal: Mathematical Problems in Engineering
 (2018)

---

## Abstract

The aims of this study contribute to a new hybrid model by combining ensemble empirical mode decomposition (EEMD) with multikernel function least square support vector machine (MKLSSVM) optimized by hybrid gravitation search algorithm (HGSA) for short-term wind speed prediction. In the forecasting process, EEMD is adopted to make the original wind speed data decomposed into intrinsic mode functions (IMFs) and one residual firstly. Then, partial autocorrelation function (PACF) is applied to identify the correlation between the corresponding decomposed components. Subsequently, the MKLSSVM using multikernel function of radial basis function (RBF) and polynomial (Poly) kernel function by weight coefficient is exploited as core forecasting engine to make the short-term wind speed prediction. To improve the regression performance, the binary-value GSA (BGSA) in HGSA is utilized as feature selection approach to remove the ineffective candidates and reconstruct the most relevant feature input-matrix for the forecasting engine, while real-value GSA (RGSA) makes the parameter combination optimization of MKLSSVM model. In the end, these respective decomposed subseries forecasting results are combined into the final forecasting values by aggregate calculation. Numerical results and comparable analysis illustrate the excellent performance of the EEMD-HGSA-MKLSSVM model when applied in the short-term wind speed forecasting.


---

## Body

## 1. Introduction

Owing to the abundant, renewable, and economical characteristics, the exploitation and utilization technique of renewable wind energy have attracted extensive attention of the scientific researchers. Wind energy has been considered as an effective way to address the global energy demands and eliminate green-house gas emissions [1]. In the past few years, wind energy has experienced fast growth worldwide.World Wind Energy Association reports that the total installed wind turbine capacity of the top 10 countries by the end of 2016 has approximately amounted to 410.613 GW and all the wind turbines worldwide by mid-2016 can generate about 4.7% of the global electricity demand [2]. However, for the high fluctuation and nonlinear and uncontrollable nature of wind speed, the integration of power system with large capacity of wind power has brought new challenges to the operation security and reliability of power system and the management of wind farms. Accurate short-term wind power output forecasting has been considered as one of the most economical and effective approaches to eliminate these problems; therefore, wind speed forecasting is a fundamental task in the routine operation management of wind farms [2–4].Over the past decades, many methods and models, mainly including physical model, statistical model, and artificial intelligent method, are widely applied to predict the short-term wind speed [5, 6]. Physical model is generally applied in the large-term wind speed forecasting by usage of detailed meteorological data and environmental information, while statistical models are constructed commonly for short-term wind speed forecasting by revealing explicitly the linear relationship among the wind speed time series [7, 8]. Different from physical models and statistical methods, the artificial intelligent methods can tackle nonlinear problems better, thus, they are the most popular and extensive approaches to apply in the short-term wind speed forecasting. To name a few here, backpropagation neural network (BPNN) [9], artificial neural networks (ANN) [10], fuzzy neural network (FNN) [11], support vector machine (SVM) [12–14], extreme learning machine (ELM) [2, 15], and least square support vector machine (LSSVM) [5] are mainly artificial intelligent methods for wind speed prediction.As stated in [5], the single artificial intelligent model cannot work well when applied in wind speed forecasting in that wind speed exhibits high nonlinearity. Wind speed forecasting by the single model using directly the raw wind speed data without disposal is easily subjected from large errors; hence, multiscale decomposition or denoising processing techniques are utilized to preprocess wind data, and intelligent algorithms are used to tune the parameters in the forecasting engine. For example, Liu at al. [13] developed a hybrid forecasting model combining Wavelet Transform (WT) with SVM tuned by general algorithm (GA). Meng et al. [16] applied the signal analysis method WPD to realize the decomposition of the original wind speed data into several different subseries; then, each decomposed component with different frequency was submitted to ANN tuned by crisscross optimization algorithm for the multistep wind speed forecasting. Wang et al. [17] developed a hybrid prediction method using EEMD and BPNN. Abdoos [18] took advantages of the combination of VMD with ELM for short-term wind power prediction. These hybrid forecasting models discussed above improve the prediction performance mainly by integration of the individual advantages of signal preprocessing technique and optimization algorithm and artificial intelligent model.Among these data preprocessing-based techniques as discussed and analyzed above, WT has sensitivity in the choice of threshold and the figuration of its wavelet basis should be determined beforehand, while EMD is sensitive to noise and suffers from mode mixing problems [6, 19]. EEMD method can eliminate the drawbacks of the decomposition approaches to some extent. EEMD is an empirical and self-adaptive signal processing approach which is widely used to analyze the nonlinear and nonstationary signal so that we use EEMD to decompose and analyze the original wind speed data in this study. LSSVM, proposed by Suykens [20], is an improved version of SVM, which lowers calculation complexity by translating convex quadratic programming problems into solving linear equations [21]. LSSVM can exhibit some advantages in solving small samples, nonlinearity, and pattern recognition with excellent generalization ability [22] and has been successfully applied in time series-based wind speed forecasting [10, 21, 23], and therefore LSSVM algorithm is adopted as the core forecasting engine for short-term wind speed forecasting.Even though these signal decomposition based models have obtained good forecasting results, Wang et al. [24] pointed out that not all decomposed subseries are a benefit for the final wind speed forecasting. To address this problem, the feature selection method is utilized widely [25, 26]. In [8], Kullback-Leibler divergence-based and energy-based feature selections were exploited to identify the illusive components caused by the decomposed method EEMD. In [27], Salcedo-Sanz developed a hybrid model of physical model and ELM, where coral reefs optimization algorithm (CRO) was utilized as feature selection to select the useful meteorological predictive information from the output of the physical approach. In the hybrid model, removing the ineffective input candidate and decreasing the dimension of the input-matrix by feature selection, ELM model can better train and regress, thus improving the forecasting performance. In [28], a hybrid GSA integrating the binary-value GSA and real-value GSA is introduced for feature selection and optimization of the weights and biases in the ELM to diagnose fault of rolling element bearings. In these hybrid algorithms, feature selection removes the ineffective variables and determines the useful input candidate to construct the input-matrix for the forecasting engine while the optimization algorithm tunes the parameters in the forecasting engine other than random initialization, and therefore the hybrid forecasting model can obtain better forecasting results.Inspired by these forecasting mechanisms, a novel hybrid model EEMD-HGSA-MKLSSVM is proposed for short-term wind speed prediction. Four sets of actual historical wind speed data from a wind farm located in Anhui of China are utilized as training and test samples to evaluate the proposed forecasting model. Considering the previous studies in the same research fields, the main works and contributions of this paper are summarized as(i)
A novel proposed forecasting model takes individual advantage of EEMD and HGSA and MKLSSVM to enhance wind speed prediction accuracy.(ii)
To improve the regression performance, radial basis function (RBF) with local exploitation capacity and polynomial (poly) kernel function with global exploration capacity are employed to construct multikernel function by weight coefficient for LSSVM, namely, MKLSSVM.(iii)
The forecasting accuracy and stability of the forecasting engine are enhanced by identification of the useful input variables and determination of optimal parameters through HGSA algorithm simultaneously.(iv)
To examine the performance of the proposed combined architecture, a number of simulation experiments are carried out and compared using four sets of wind speed data.The remaining of this study are structured as follows. In Section2, the individual models, including EEMD and HGSA and MKLSSVM, are introduced. Section 3 illustrates the detailed working principle of the proposed EEMD-HGSA-MKLSSVM model and performance evaluation indices. Case studies are implemented to evaluate the proposed hybrid model for short-term wind speed forecasting in Section 4. Conclusions are drawn in the final section.


## 2. Methodology

### 2.1. Wind Speed Decomposition Method

EMD, a self-adaptive signal analysis technique, is developed to analyze and decompose nonlinear signals by sifting process [6]. However, EMD easily suffers from the mode mixing problem that defined as a single IMF containing signals with dramatically disparate scales or a component of a similar scale residing in different IMFs, which causes easily intermittency in analyzing signals [6]. To eliminate the mode mixing problems caused by the EMD, a novel nonlinear signal analysis method EEMD was developed by adding white noises with finite amplitude to the original signals and offsetting themselves through ensemble averaging [17]. The analysis process of signals by EEMD algorithm can be described as in the following steps:(i)
Step1. Add Gaussian distribution white noise n(t) with finite amplitude to the original signal data x(t) to obtain a new signal s(t) expressed as (1)st=xt+nt.(ii)
Step2. Decompose s(t) into a series of IMF components ci and residual component rn using a standard EMD method. After decomposition, s(t) can be mathematically expressed as (2)st=∑i=1Nci+rN.(iii)
Step3. Repeating step from 1 to 2 and add random white noise each time, these N groups of different white noise have characteristics of uncorrelated relationship and its statistical mean is zero.(iv)
Step4. Offsetting the impact of the Gaussian white noise by the final mean of the corresponding IMFs as (3), by the same way, the final residue can be obtained as(3)cj=1N∑i=1Ncij.(v)
Step5. In the end, the noise-free signal data are obtained by reconstruction of the IMF components cj and the residual component rN.The effects of the added Gaussian distribution white noise are ensured by the statistical rule proved by Wu and Huang [29],(4)εne=εNe,whereNe is the ensemble members, ε is the amplitude of the added noise, and εne is defined as the difference between the input signal and the corresponding IMFs.To better illustrate the effectiveness of EEMD in overcoming the mode mixing problem, a given synthetic test signaly(t) expressed as (5) consisting of a sinusoid signal y1(t) expressed as (6) and an intermittent signal y2(t) expressed as (7), which are displayed in Figure 1, is utilized to test EMD and EEMD.(5)y=y1+y2,(6)y1=sin⁡20πt0≤t≤0.4,(7)y2=0.3sin⁡100πt0.05≤t≤0.15,0.2sin⁡300πt0.325≤t≤0.375,0.3sin⁡500πt0.225≤t≤0.275,0.Figure 1
Synthesized signal y.The decomposed results of the synthesized signaly by EMD and EEMD are shown in the Figures 2 and 3, respectively. It is obviously seen from Figures 2 and 3 that the mode mixing problems exist in the different components decomposed by EMD, while the mode mixing problems have been eliminated by EEMD and the intermittent signals embedded in the synthesized signal have been extracted successfully.Figure 2
Decomposed results of the synthesized signaly by EMD.Figure 3
Decomposed results of the synthesized signaly by EEMD.


### 2.2. Multikernel LSSVM(MKLSSVM)

In this study, LSSVM is adopted in that it is expert in addressing small sample problems [21] and has high generalization performance [30]. LSSVM is a type of powerful artificial intelligence technology based on the structural risk minimization. The basic principle of single output LSSVM regression is shown as follows.Assume that the training sample set{xi,yi}, where i=1⋯N, N is the total number of training samples; xi is input training sample and yi is its corresponding output. The regression function can be expressed as(8)y=wTφx+b,wherew and b denote the weight vector and the bias term, respectively. φ(·) is a nonlinear mapping function which maps the training samples into a high-dimension feature space where regression is carried out. The regression can be calculated by minimizing a cost function expressed as(9)min⁡C=min⁡12ωTω+12γ∑i=1Nei2,subject to the equality constraints as(10)yi=ωTφxi+b+ei,i=1⋯N.The first part of the cost function in (9) is utilized to regularize the weight sizes and penalize weights and γ is a regularization parameter which is optimized by user to control the trade-off between the bias and variance of LSSVM. The Lagrange function is constructed as follows to solve the convex optimization problem:(11)Lω,b,e,α=C-∑i=1NβiωTφx+b-yi+ei,whereβi are Lagrange multipliers. By partially differentiating with respect to ω, b, e, and βi, eliminating ω and e, the solution of (11) can be obtained by the following linear regression function:(12)y=∑i=1Nβikx,xi+b,wherek(x,xi) is a positive definite kernel function which meets Mercer’s condition. The different type of kernel function influences the regressive performance of LSSVM model. In this paper, to improve the generalization capacity, a weighted multikernel LSSVM combining RBF kernel function and Poly kernel function is constructed and expressed as (13)kmixxi,xj=μkRBFxi,xj+1-μkPolyxi,xj,whereμ represents the weight coefficient within [01] and RBF kernel function and Poly kernel function are expressed as (14) and (15), respectively.(14)kRBFxi,xj=exp⁡-xi-xj22δ2,(15)kpolyxi,xj=xi·xj+rd.Whenμ=1, the weighted multiple-kernel function exhibits the characteristic of RBF function and μ=0; this shows the characteristic of Poly function. By adjusting the μ value, the weighted multiple-kernel function can be suitable for different input samples.


### 2.3. Hybrid Gravitational Search Algorithm (HGSA)

GSA, a novel heuristic optimization algorithm, was firstly developed by Rashedi [31], and the standard GSA has been successfully used to solve the engineering optimization problems in real-value parameters domain. However, in the actual engineering application, there exist many binary-encoded optimization problems, such as feature selection, which need to be solved. In this study, the optimization objectives including input variables binary feature selection as well as the real-value kernel parameters and weighted coefficients in MKLSSVM are dealt simultaneously; thus, we develop a hybrid GSA combing BGSA for feature selection with standard GSA for parameter optimization.


#### 2.3.1. Gravitational Search Algorithm (GSA)

In GSA algorithm, the performance of the agent is evaluated by their masses, namely, the heavier agent represents the optimal solution; thus the agent’s mass is considered as the objective. All the agents move towards the agents with heavier masses by the gravitational force between them. As a result, after many iterations, the heavier masses with the higher fitness are obtained.Assume that there areN random agents in the GSA. Firstly, the speed and position of each agent should be initialized and the position of ith agent is {(Xi=xi1,xi2,⋯,xid,…,xiD)∣i=1,2,…,N} and D denotes the dimension of the search space. At the tth iteration, the mass of the ith agent is mathematically expressed as (16)mit=fitit-worsttbestt-worstt,Mit=mit∑j=1Nmjt,wherefit(t), M(t), worst(t), and best(t) are fitness value, mass, worst, and best value, respectively. For the wind speed forecasting error minimization problem, the worst(t) and best(t) are expressed as (17)bestt=minj∈1,⋯N⁡fitjt,worstt=maxj∈1,⋯N⁡fitjt,The gravitational force that theith mass xid(t) acts on the jth mass xjd(t) is expressed as follows according to the Newton gravitation theory:(18)Fijdt=GtMit×MjtXit,Xjt2+εxjdt-xidt,Gt=G0×exp⁡-αtitermax,whereG0, α, k, and itermax are the initial gravitational constant, attenuation factor, the current iteration number, and the maximum iteration number, respectively.The resultant force that other masses act on theith mass is calculated as (19) by a randomly weighted sum of Fijd(t). Afterwards, the acceleration aid(t), velocity vid(t), and position xid(t) of an agent are calculated as (20):(19)Fidt=∑i≠jrandj×Fijdt(20)vidt+1=randi×vidt+aidtxidt+1=xidt+vidt+1whererandi and randj are random value within [01] and aid(t)=Fid(t)/Mi(t).


#### 2.3.2. Binary GSA(BGSA)

The real-value GSA is exploited to solve the optimization problems in the continuous space, which can be not directly employed to solve the binary problems. The binary GSA (BGSA), developed by Rashedi et al. [32], is used to deal with discrete binary problems, and its working mechanism is that the velocity of each agent is converted by the Hyperbolic tangent function into a probability value which is expressed as (21)Svidt=tanh⁡vidt,wheretanh⁡(·) denotes the Hyperbolic tangent function. Based on the above (21), each dimension in the discrete binary space takes on “0” or “1” binary value by (22)ifrand<Svidtxidt+1=complementxidtelsexidt+1=xidtwherecomplement(·) stands for the logical negation.


## 2.1. Wind Speed Decomposition Method

EMD, a self-adaptive signal analysis technique, is developed to analyze and decompose nonlinear signals by sifting process [6]. However, EMD easily suffers from the mode mixing problem that defined as a single IMF containing signals with dramatically disparate scales or a component of a similar scale residing in different IMFs, which causes easily intermittency in analyzing signals [6]. To eliminate the mode mixing problems caused by the EMD, a novel nonlinear signal analysis method EEMD was developed by adding white noises with finite amplitude to the original signals and offsetting themselves through ensemble averaging [17]. The analysis process of signals by EEMD algorithm can be described as in the following steps:(i)
Step1. Add Gaussian distribution white noise n(t) with finite amplitude to the original signal data x(t) to obtain a new signal s(t) expressed as (1)st=xt+nt.(ii)
Step2. Decompose s(t) into a series of IMF components ci and residual component rn using a standard EMD method. After decomposition, s(t) can be mathematically expressed as (2)st=∑i=1Nci+rN.(iii)
Step3. Repeating step from 1 to 2 and add random white noise each time, these N groups of different white noise have characteristics of uncorrelated relationship and its statistical mean is zero.(iv)
Step4. Offsetting the impact of the Gaussian white noise by the final mean of the corresponding IMFs as (3), by the same way, the final residue can be obtained as(3)cj=1N∑i=1Ncij.(v)
Step5. In the end, the noise-free signal data are obtained by reconstruction of the IMF components cj and the residual component rN.The effects of the added Gaussian distribution white noise are ensured by the statistical rule proved by Wu and Huang [29],(4)εne=εNe,whereNe is the ensemble members, ε is the amplitude of the added noise, and εne is defined as the difference between the input signal and the corresponding IMFs.To better illustrate the effectiveness of EEMD in overcoming the mode mixing problem, a given synthetic test signaly(t) expressed as (5) consisting of a sinusoid signal y1(t) expressed as (6) and an intermittent signal y2(t) expressed as (7), which are displayed in Figure 1, is utilized to test EMD and EEMD.(5)y=y1+y2,(6)y1=sin⁡20πt0≤t≤0.4,(7)y2=0.3sin⁡100πt0.05≤t≤0.15,0.2sin⁡300πt0.325≤t≤0.375,0.3sin⁡500πt0.225≤t≤0.275,0.Figure 1
Synthesized signal y.The decomposed results of the synthesized signaly by EMD and EEMD are shown in the Figures 2 and 3, respectively. It is obviously seen from Figures 2 and 3 that the mode mixing problems exist in the different components decomposed by EMD, while the mode mixing problems have been eliminated by EEMD and the intermittent signals embedded in the synthesized signal have been extracted successfully.Figure 2
Decomposed results of the synthesized signaly by EMD.Figure 3
Decomposed results of the synthesized signaly by EEMD.


## 2.2. Multikernel LSSVM(MKLSSVM)

In this study, LSSVM is adopted in that it is expert in addressing small sample problems [21] and has high generalization performance [30]. LSSVM is a type of powerful artificial intelligence technology based on the structural risk minimization. The basic principle of single output LSSVM regression is shown as follows.Assume that the training sample set{xi,yi}, where i=1⋯N, N is the total number of training samples; xi is input training sample and yi is its corresponding output. The regression function can be expressed as(8)y=wTφx+b,wherew and b denote the weight vector and the bias term, respectively. φ(·) is a nonlinear mapping function which maps the training samples into a high-dimension feature space where regression is carried out. The regression can be calculated by minimizing a cost function expressed as(9)min⁡C=min⁡12ωTω+12γ∑i=1Nei2,subject to the equality constraints as(10)yi=ωTφxi+b+ei,i=1⋯N.The first part of the cost function in (9) is utilized to regularize the weight sizes and penalize weights and γ is a regularization parameter which is optimized by user to control the trade-off between the bias and variance of LSSVM. The Lagrange function is constructed as follows to solve the convex optimization problem:(11)Lω,b,e,α=C-∑i=1NβiωTφx+b-yi+ei,whereβi are Lagrange multipliers. By partially differentiating with respect to ω, b, e, and βi, eliminating ω and e, the solution of (11) can be obtained by the following linear regression function:(12)y=∑i=1Nβikx,xi+b,wherek(x,xi) is a positive definite kernel function which meets Mercer’s condition. The different type of kernel function influences the regressive performance of LSSVM model. In this paper, to improve the generalization capacity, a weighted multikernel LSSVM combining RBF kernel function and Poly kernel function is constructed and expressed as (13)kmixxi,xj=μkRBFxi,xj+1-μkPolyxi,xj,whereμ represents the weight coefficient within [01] and RBF kernel function and Poly kernel function are expressed as (14) and (15), respectively.(14)kRBFxi,xj=exp⁡-xi-xj22δ2,(15)kpolyxi,xj=xi·xj+rd.Whenμ=1, the weighted multiple-kernel function exhibits the characteristic of RBF function and μ=0; this shows the characteristic of Poly function. By adjusting the μ value, the weighted multiple-kernel function can be suitable for different input samples.


## 2.3. Hybrid Gravitational Search Algorithm (HGSA)

GSA, a novel heuristic optimization algorithm, was firstly developed by Rashedi [31], and the standard GSA has been successfully used to solve the engineering optimization problems in real-value parameters domain. However, in the actual engineering application, there exist many binary-encoded optimization problems, such as feature selection, which need to be solved. In this study, the optimization objectives including input variables binary feature selection as well as the real-value kernel parameters and weighted coefficients in MKLSSVM are dealt simultaneously; thus, we develop a hybrid GSA combing BGSA for feature selection with standard GSA for parameter optimization.


### 2.3.1. Gravitational Search Algorithm (GSA)

In GSA algorithm, the performance of the agent is evaluated by their masses, namely, the heavier agent represents the optimal solution; thus the agent’s mass is considered as the objective. All the agents move towards the agents with heavier masses by the gravitational force between them. As a result, after many iterations, the heavier masses with the higher fitness are obtained.Assume that there areN random agents in the GSA. Firstly, the speed and position of each agent should be initialized and the position of ith agent is {(Xi=xi1,xi2,⋯,xid,…,xiD)∣i=1,2,…,N} and D denotes the dimension of the search space. At the tth iteration, the mass of the ith agent is mathematically expressed as (16)mit=fitit-worsttbestt-worstt,Mit=mit∑j=1Nmjt,wherefit(t), M(t), worst(t), and best(t) are fitness value, mass, worst, and best value, respectively. For the wind speed forecasting error minimization problem, the worst(t) and best(t) are expressed as (17)bestt=minj∈1,⋯N⁡fitjt,worstt=maxj∈1,⋯N⁡fitjt,The gravitational force that theith mass xid(t) acts on the jth mass xjd(t) is expressed as follows according to the Newton gravitation theory:(18)Fijdt=GtMit×MjtXit,Xjt2+εxjdt-xidt,Gt=G0×exp⁡-αtitermax,whereG0, α, k, and itermax are the initial gravitational constant, attenuation factor, the current iteration number, and the maximum iteration number, respectively.The resultant force that other masses act on theith mass is calculated as (19) by a randomly weighted sum of Fijd(t). Afterwards, the acceleration aid(t), velocity vid(t), and position xid(t) of an agent are calculated as (20):(19)Fidt=∑i≠jrandj×Fijdt(20)vidt+1=randi×vidt+aidtxidt+1=xidt+vidt+1whererandi and randj are random value within [01] and aid(t)=Fid(t)/Mi(t).


### 2.3.2. Binary GSA(BGSA)

The real-value GSA is exploited to solve the optimization problems in the continuous space, which can be not directly employed to solve the binary problems. The binary GSA (BGSA), developed by Rashedi et al. [32], is used to deal with discrete binary problems, and its working mechanism is that the velocity of each agent is converted by the Hyperbolic tangent function into a probability value which is expressed as (21)Svidt=tanh⁡vidt,wheretanh⁡(·) denotes the Hyperbolic tangent function. Based on the above (21), each dimension in the discrete binary space takes on “0” or “1” binary value by (22)ifrand<Svidtxidt+1=complementxidtelsexidt+1=xidtwherecomplement(·) stands for the logical negation.


## 2.3.1. Gravitational Search Algorithm (GSA)

In GSA algorithm, the performance of the agent is evaluated by their masses, namely, the heavier agent represents the optimal solution; thus the agent’s mass is considered as the objective. All the agents move towards the agents with heavier masses by the gravitational force between them. As a result, after many iterations, the heavier masses with the higher fitness are obtained.Assume that there areN random agents in the GSA. Firstly, the speed and position of each agent should be initialized and the position of ith agent is {(Xi=xi1,xi2,⋯,xid,…,xiD)∣i=1,2,…,N} and D denotes the dimension of the search space. At the tth iteration, the mass of the ith agent is mathematically expressed as (16)mit=fitit-worsttbestt-worstt,Mit=mit∑j=1Nmjt,wherefit(t), M(t), worst(t), and best(t) are fitness value, mass, worst, and best value, respectively. For the wind speed forecasting error minimization problem, the worst(t) and best(t) are expressed as (17)bestt=minj∈1,⋯N⁡fitjt,worstt=maxj∈1,⋯N⁡fitjt,The gravitational force that theith mass xid(t) acts on the jth mass xjd(t) is expressed as follows according to the Newton gravitation theory:(18)Fijdt=GtMit×MjtXit,Xjt2+εxjdt-xidt,Gt=G0×exp⁡-αtitermax,whereG0, α, k, and itermax are the initial gravitational constant, attenuation factor, the current iteration number, and the maximum iteration number, respectively.The resultant force that other masses act on theith mass is calculated as (19) by a randomly weighted sum of Fijd(t). Afterwards, the acceleration aid(t), velocity vid(t), and position xid(t) of an agent are calculated as (20):(19)Fidt=∑i≠jrandj×Fijdt(20)vidt+1=randi×vidt+aidtxidt+1=xidt+vidt+1whererandi and randj are random value within [01] and aid(t)=Fid(t)/Mi(t).


## 2.3.2. Binary GSA(BGSA)

The real-value GSA is exploited to solve the optimization problems in the continuous space, which can be not directly employed to solve the binary problems. The binary GSA (BGSA), developed by Rashedi et al. [32], is used to deal with discrete binary problems, and its working mechanism is that the velocity of each agent is converted by the Hyperbolic tangent function into a probability value which is expressed as (21)Svidt=tanh⁡vidt,wheretanh⁡(·) denotes the Hyperbolic tangent function. Based on the above (21), each dimension in the discrete binary space takes on “0” or “1” binary value by (22)ifrand<Svidtxidt+1=complementxidtelsexidt+1=xidtwherecomplement(·) stands for the logical negation.


## 3. The Proposed Forecasting Strategy

### 3.1. Feature Selection and Parameter Optimization by HGSA

To eliminate the uncorrelated and ineffective variables in the input samples and avoid the drawbacks of trapping in local optima or overfitting of MKLSSVM, HGSA is used as feature selection and parameter optimization to solve this problems simultaneously. As seen from Figure4, a hybrid optimization algorithm HGSA is exploited for the optimization of the kernel parameters and weighted coefficient in MKLSSVM by GSA and feature selection of input variables by BGSA to improve the forecasting performance of the forecasting strategy. During the optimization procedure, the root mean square error (RMSE), expressed in (23), between the training results of MKLSSVM and the measured wind speed data is used as fitness function.(23)RMSE=1N∑t=1Nst-s^t2,Figure 4
Feature selection and parameter optimization using HGSA.wheres(t) and s^(t) denote the measured wind speed and the forecasting wind speed value, respectively. N stands for the total number of samples.With this in mind, the hybrid optimization problem can be illustrated in Figure5, and there are n real-valued kernel parameters and weighted coefficient and m binary-valued parameters. The binary value “1” in BGSA algorithm means that the input variable is selected while “0” means that the input variable is not considered. Thus, the initialization encoding dimension of an agent in HGSA is set as (n+m)-length vector.Figure 5
Solution representation.


### 3.2. Specific Steps of EEMD-HGSA-MKLSSVM Model

In this study, LSSVM, based on weighted multikernel function, is utilized as the core forecasting engine in the forecasting strategy. The single LSSVM model has many advantages in solving small sample and nonlinearity. However, owing to the irregularity and randomness of wind speed time series, it cannot obtain the favorable prediction performance when independently applied in wind speed forecasting. To better catch the characteristics of wind speed, a hybrid approach combining EEMD and MKLSSVM with HGSA together is proposed as shown in Figure6 and its working process is as follows:(i)
Use EEMD method to break down the empirical wind speed into IMFs andRes with different frequency.(ii)
Prior to forecasting by MKLSSVM model, the PACF that is widely used as a lag identification approach in Auto Regression (AR) (p) is applied to determine the correlation coefficients of the inputs. When the PACF values at lags bigger than p are approximately independent N(0,1/n) random variables, the lag of the sample can be determined as p.(iii)
To lower the forecasting difficulties of the MKLSSVM model, the inputs are normalized linearly to interval[0,1].(iv)
Train the MKLSSVMs optimized by the HGSA algorithm using the different frequency subseries (IMFs andRes). The 1st-480th wind speed series in Figure 7 are adopted as the training dataset.(v)
Apply the well-trained HGSA-MKLSSVM models to do the multistep ahead wind speed forecasting using each subseries. The subsequent 481th-576th sampling points in Figure 7 are used as the test dataset. Additionally, a rolling forecasting mechanism is adopted in the prediction processes.(vi)
Obtain the final forecasting results by aggregating the calculation after denormalization. In the end, compare the forecasting performance between EEMD-HGSA-MKLSSVM and other wind speed forecasting methods.Figure 6
The flowchart of the hybrid EEMD-HGSA-MLLSSVM model.Figure 7
Four datasets of empirical wind speed time series (red and blue lines stand for testing and training data, respectively).
(a)
Wind speed of different seasons (b)
Empirical wind speed time series for training and testing


### 3.3. Forecasting Performance Evaluation Indices

To verify the prediction performance of the hybrid EEMD-HGSA-MKLSSVM model, three statical indices, namely RMSE, MAE, and MAPE, are utilized to measure the prediction accuracy, and these indices are expressed as (23), (24), and (25).(24)MAE=1N∑t=1Nst-s^t,(25)MAPE=1N∑t=1Nst-s^tst×100%,wheres(i), s^(t), and N are the same meaning as that in (23). It is well known that small statistical index values indicate high forecast accuracy. Among the statistical indices, the MAE reveals the similarity between the predicted and observed wind speed value, MAPE indicates forecasting percent error at every observation spot, whereas the RMSE measures the overall deviation; therefore, RMSE is adopted as the fitness function in this paper.To further illustrate the improvement of modelI over model II, three improved percentage indices based on the RMSE, MAPE and MAE, are introduced to describe the improvement degree, and these three indices PRMSE, PMAE and PMAPE are expressed as follows [2]:(26)PRMSE=RMSE1-RMSE2RMSE1×100%,(27)PMAE=MAE1-MAE2MAE1×100%,(28)PMAPE=MAPE1-MAPE2MAPE1×100%.


## 3.1. Feature Selection and Parameter Optimization by HGSA

To eliminate the uncorrelated and ineffective variables in the input samples and avoid the drawbacks of trapping in local optima or overfitting of MKLSSVM, HGSA is used as feature selection and parameter optimization to solve this problems simultaneously. As seen from Figure4, a hybrid optimization algorithm HGSA is exploited for the optimization of the kernel parameters and weighted coefficient in MKLSSVM by GSA and feature selection of input variables by BGSA to improve the forecasting performance of the forecasting strategy. During the optimization procedure, the root mean square error (RMSE), expressed in (23), between the training results of MKLSSVM and the measured wind speed data is used as fitness function.(23)RMSE=1N∑t=1Nst-s^t2,Figure 4
Feature selection and parameter optimization using HGSA.wheres(t) and s^(t) denote the measured wind speed and the forecasting wind speed value, respectively. N stands for the total number of samples.With this in mind, the hybrid optimization problem can be illustrated in Figure5, and there are n real-valued kernel parameters and weighted coefficient and m binary-valued parameters. The binary value “1” in BGSA algorithm means that the input variable is selected while “0” means that the input variable is not considered. Thus, the initialization encoding dimension of an agent in HGSA is set as (n+m)-length vector.Figure 5
Solution representation.


## 3.2. Specific Steps of EEMD-HGSA-MKLSSVM Model

In this study, LSSVM, based on weighted multikernel function, is utilized as the core forecasting engine in the forecasting strategy. The single LSSVM model has many advantages in solving small sample and nonlinearity. However, owing to the irregularity and randomness of wind speed time series, it cannot obtain the favorable prediction performance when independently applied in wind speed forecasting. To better catch the characteristics of wind speed, a hybrid approach combining EEMD and MKLSSVM with HGSA together is proposed as shown in Figure6 and its working process is as follows:(i)
Use EEMD method to break down the empirical wind speed into IMFs andRes with different frequency.(ii)
Prior to forecasting by MKLSSVM model, the PACF that is widely used as a lag identification approach in Auto Regression (AR) (p) is applied to determine the correlation coefficients of the inputs. When the PACF values at lags bigger than p are approximately independent N(0,1/n) random variables, the lag of the sample can be determined as p.(iii)
To lower the forecasting difficulties of the MKLSSVM model, the inputs are normalized linearly to interval[0,1].(iv)
Train the MKLSSVMs optimized by the HGSA algorithm using the different frequency subseries (IMFs andRes). The 1st-480th wind speed series in Figure 7 are adopted as the training dataset.(v)
Apply the well-trained HGSA-MKLSSVM models to do the multistep ahead wind speed forecasting using each subseries. The subsequent 481th-576th sampling points in Figure 7 are used as the test dataset. Additionally, a rolling forecasting mechanism is adopted in the prediction processes.(vi)
Obtain the final forecasting results by aggregating the calculation after denormalization. In the end, compare the forecasting performance between EEMD-HGSA-MKLSSVM and other wind speed forecasting methods.Figure 6
The flowchart of the hybrid EEMD-HGSA-MLLSSVM model.Figure 7
Four datasets of empirical wind speed time series (red and blue lines stand for testing and training data, respectively).
(a)
Wind speed of different seasons (b)
Empirical wind speed time series for training and testing


## 3.3. Forecasting Performance Evaluation Indices

To verify the prediction performance of the hybrid EEMD-HGSA-MKLSSVM model, three statical indices, namely RMSE, MAE, and MAPE, are utilized to measure the prediction accuracy, and these indices are expressed as (23), (24), and (25).(24)MAE=1N∑t=1Nst-s^t,(25)MAPE=1N∑t=1Nst-s^tst×100%,wheres(i), s^(t), and N are the same meaning as that in (23). It is well known that small statistical index values indicate high forecast accuracy. Among the statistical indices, the MAE reveals the similarity between the predicted and observed wind speed value, MAPE indicates forecasting percent error at every observation spot, whereas the RMSE measures the overall deviation; therefore, RMSE is adopted as the fitness function in this paper.To further illustrate the improvement of modelI over model II, three improved percentage indices based on the RMSE, MAPE and MAE, are introduced to describe the improvement degree, and these three indices PRMSE, PMAE and PMAPE are expressed as follows [2]:(26)PRMSE=RMSE1-RMSE2RMSE1×100%,(27)PMAE=MAE1-MAE2MAE1×100%,(28)PMAPE=MAPE1-MAPE2MAPE1×100%.


## 4. Case Study and Forecasting Results

### 4.1. Wind Speed Data Statistical Description

In this study, the empirical wind speed data were collected and stored in a wind farm located in Anhui province of China to verify the proposed model. The wind speed time series were measured every 10 min from the anemometer at the top of a nacelle, and the wind turbines with 88m are installed at the top of the 300 m mountain (32°28′N, 118°26′E). Accurate half-an-hour to several days ahead wind speed forecasting are very beneficial to routine management of wind farms, and thus the original 10 m wind speed time series are transformed to 30 min by averaging. The developed models in the most studies are tested and evaluated with no more than two sets of wind speed data, and the conclusion that the models are better than the other approaches is generally not convincing enough [26]. Therefore, four sets of wind speed data in 2015 are randomly selected to evaluate the proposed model. The historical wind speed time series of four seasons in 2015 are displayed in the Figure 7(a). The red mark variables in Figure 7(a), which are clearly shown in Figure 7(b), are utilized to train and test the proposed models. The blue variables and red ones in Figure 7(b) are employed to train and test the models, respectively.The quantity of input samples of LSSVM model affects the wind speed forecasting performance. However, there are not any unified standard and clear definition for section of the input sample quantity [5]. As for the LSSVM which is expert in solving small samples and the suggestion made by Wang et al. [6], four sets of 576 half-hourly wind speed samples in 2015 are selected randomly from Spring, Summer, Autumn, and Winter, respectively. Thereinto, the 1st-480th sampling points are employed to train the proposed model while the subsequent 481th-576th data are applied to test the well-trained model.The descriptive statistics of the original wind speed times series list in Table1. It can be obviously seen that the empirical wind speed time series presents the characters of fluctuations and stochastic volatility, fluctuating around 7.80 m/s, 5.48 m/s, 8.47 m/s, and 5.55 m/s, respectively. There is no apparent regularity in the wind speed series.Table 1
Statistical description of empirical wind speed (m/s).
data set NO. Max Min Mean St.dev. Median Data set A All samples 576 14.59 1.22 7.80 2.79 7.91 Training data 480 14.59 1.22 8.16 2.68 7.91 Testing data 96 10.68 1.32 6.01 2.64 6.01  Data set B All samples 576 12.44 0.94 5.48 2.46 6.69 Training data 480 12.44 0.94 5.58 2.46 6.69 Testing data 96 9.91 1.34 4.99 2.38 5.62  Data set C All samples 576 15.34 0.87 8.47 3.08 8.10 Training data 480 15.34 0.87 8.56 3.28 8.10 Testing data 96 11.57 4.08 7.99 1.69 7.82  Data set D All samples 576 12.87 0.41 5.55 2.40 6.64 Training data 480 12.87 0.41 5.42 2.52 6.64 Testing data 96 9.91 2.57 6.17 1.56 6.25


### 4.2. Wind Speed Preprocessing Using EEMD

Prior to submitting the original wind speed to the forecasting engine, the signal decomposition technique EEMD is applied to break down the empirical wind speed into several relatively stable IMFs and oneRes for reducing the forecasting difficulties of HGSA-MKLSSVM. According to [8], the amplitude of the added white noise and the ensemble number in this study are set 0.2 and 100, respectively. As shown in Figure 8, the empirical wind speed time series are broken down into six IMFs and one Res. It is obviously observed from the figures that all the components present distinct respective characteristic, frequencies decrease gradually from IMF1 to Res, and high-frequency IMF1 and IMF2 reflect the stochastic characteristic of wind speed data, while the signals IMF3 to IMF6 with periodic trend features represent the periodic components of the wind speed data and Res are named the trend components.Figure 8
Wind speed decomposition using EEMD.
(a)
Wind speed dataA (b)
Wind speed dataB (c)
Wind speed dataC (d)
Wind speed dataD


### 4.3. Input-Matrix Construction by PACF and BGSA

Before the forecasting processes are carried out by the HGSA-MKLSSVM model, the input variables matrix requires to be determined. The PACF technique is employed to identify the correlations among each corresponding decomposed component for determination of the input combination. The PACF of four sets of the original wind speed and their decomposed components are illustrated in Figure9. As seen in Figure 9, the PACF values of all wind speed and decomposed components from lags 1 to 30 are calculated. For original wind speed data A, PACF values at lags bigger than 7 are between the 95% confidence lines (the upper and lower lines); thus, the time lag is determined as 7 which can be used to identify the input variables dimension, and the 7 previous continuous wind speed time series contribute the most correlative information to forecast the subsequent wind speed value. The lag order of the original wind speed data and their decomposed components are shown in Table 2. For the original wind speed data A, the input vector for the MKLSSVM model can be described as X=xt-1, xt-2, ⋯, xt-7 and the corresponding forecasting value is Y=xt-k, where k(≥7) is the forecasting horizontal. The other input variables combination for the HGSA-MKLSSVM model determined by PACF values is expressed in Table 3. The feature selection results obtained by BGSA for different subseries are shown in Table 4.Table 2
Lag order of the autoregressive process.
Samples Original IMF1 IMF2 IMF3 IMF4 IMF5 IMF6 Res A 7 11 10 15 8 12 7 8 B 8 10 13 14 11 8 7 9 C 7 11 9 10 7 13 9 6 D 6 11 10 8 8 10 9 10Table 3
Inputs combinations of wind speed and their decomposed components for the HGSA-MKLSSVM model.
Data set Time series input variable combination A Original x t - 1, xt-2, xt-3, xt-4, xt-5, xt-6, xt-7 IMF1 x t - 1, xt-2, xt-3, xt-4,xt-5, xt-6, xt-7, xt-8, xt-9, xt-10, xt-11 IMF2 x t - 1, xt-2, xt-3, xt-4, xt-5, xt-6, xt-7, xt-8, xt-9, xt-10 IMF3 x t - 1, xt-2, xt-3, xt-4, xt-5, xt-6, xt-7, xt-8, xt-9, xt-10, xt-11, xt-12, xt-13, xt-14, xt-15 IMF4 x t - 1, xt-2, xt-3, xt-4, xt-5, xt-6, xt-7, xt-8, IMF5 x t - 1, xt-2, xt-3, xt-4, xt-5, xt-6, xt-7, xt-8, xt-9, xt-10, xt-11, xt-12 IMF6 x t - 1, xt-2, xt-3, xt-4, xt-5, xt-6, xt-7 Res x t - 1, xt-2, xt-3, xt-4, xt-5, xt-6, xt-7, xt-8  B Original x t - 1, xt-2, xt-3, xt-4, xt-5, xt-6, xt-7, xt-8 IMF1 x t - 1, xt-2, xt-3, xt-4, xt-5, xt-6, xt-7, xt-8, xt-9, xt-10 IMF2 x t - 1, xt-2, xt-3, xt-4, xt-5, xt-6, xt-7, xt-8, xt-9, xt-10, xt-11, xt-12, xt-13 IMF3 x t - 1, xt-2, xt-3, xt-4, xt-5, xt-6, xt-7, xt-8, xt-9, xt-10, xt-11, xt-12, xt-13, xt-14 IMF4 x t - 1, xt-2, xt-3, xt-4, xt-5, xt-6, xt-7, xt-8, xt-9, xt-10, xt-11 IMF5 x t - 1, xt-2, xt-3, xt-4, xt-5, xt-6, xt-7, xt-8 IMF6 x t - 1, xt-2, xt-3, xt-4, xt-5, xt-6, xt-7 Res x t - 1, xt-2, xt-3, xt-4, xt-5, xt-6, xt-7, xt-8, xt-9  C Original x t - 1, xt-2, xt-3, xt-4, xt-5, xt-6, xt-7 IMF1 x t - 1, xt-2, xt-3, xt-4, xt-5, xt-6, xt-7, xt-8, xt-9, xt-10, xt-11 IMF2 x t - 1, xt-2, xt-3, xt-4, xt-5, xt-6, xt-7, xt-8, xt-9 IMF3 x t - 1, xt-2, xt-3, xt-4, xt-5, xt-6, xt-7, xt-8, xt-9, xt-10 IMF4 x t - 1, xt-2, xt-3, xt-4, xt-5, xt-6, xt-7 IMF5 x t - 1, xt-2, xt-3, xt-4, xt-5, xt-6, xt-7, xt-8, xt-9, xt-10, xt-11, xt-12, xt-13 IMF6 x t - 1, xt-2, xt-3, xt-4, xt-5, xt-6, xt-7, xt-8, xt-9, xt-8, xt-9 Res x t - 1, xt-2, xt-3, xt-4, xt-5, xt-6  D Original x t - 1, xt-2, xt-3, xt-4, xt-5, xt-6 IMF1 x t - 1, xt-2, xt-3, xt-4, xt-5, xt-6, xt-7, xt-8, xt-9, xt-10, xt-11 IMF2 x t - 1, xt-2, xt-3, xt-4, xt-5, xt-6, xt-7, xt-8, xt-9, xt-10 IMF3 x t - 1, xt-2, xt-3, xt-4, xt-5, xt-6, xt-7, xt-8 IMF4 x t - 1, xt-2, xt-3, xt-4, xt-5, xt-6, xt-7, xt-8 IMF5 x t - 1, xt-2, xt-3, xt-4, xt-5, xt-6, xt-7, xt-8, xt-9, xt-10 IMF6 x t - 1, xt-2, xt-3, xt-4, xt-5, xt-6, xt-7, xt-8, xt-9 Res x t - 1, xt-2, xt-3, xt-4, xt-5, xt-6, xt-7, xt-8, xt-9, xt-10Table 4
Feature selection results obtained by HGSA algorithm.
Data set Time series 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 A Original 1 1 1 0 1 0 1 IMF1 1 1 0 1 1 1 0 1 1 0 1 IMF2 1 1 1 1 0 0 1 1 1 1 IMF3 1 1 0 1 1 0 1 1 1 1 1 0 1 1 1 IMF4 1 1 1 0 1 1 0 1 IMF5 1 0 1 0 1 1 1 0 1 1 1 1 IMF6 1 1 0 1 1 0 1 Res 1 1 0 1 1 1 1 1  B Original 1 0 1 1 1 1 1 1 IMF1 1 0 1 1 1 1 1 1 0 1 IMF2 1 1 0 1 0 1 1 1 1 1 1 1 0 IMF3 1 0 1 1 1 1 1 1 0 0 1 1 1 1 IMF4 1 1 0 1 1 0 1 1 1 1 1 IMF5 1 1 1 0 1 1 0 1 IMF6 1 1 0 1 1 1 1 Res 1 0 1 1 1 1 0 1 1  C Original 1 1 1 0 1 1 1 IMF1 1 0 1 1 0 1 1 0 1 1 1 IMF2 1 0 1 1 1 1 0 1 1 IMF3 1 1 1 0 1 1 0 1 1 1 IMF4 1 1 0 1 1 1 1 IMF5 1 1 1 1 0 1 1 1 1 1 0 0 1 IMF6 1 1 0 1 1 1 1 0 1 Res 1 0 1 1 1 1  D Original 1 1 1 0 1 1 IMF1 1 0 1 1 1 1 0 1 1 0 1 IMF2 1 0 1 1 0 1 1 0 1 1 IMF3 1 0 1 1 1 0 1 1 IMF4 1 1 1 0 1 1 1 1 IMF5 1 1 0 1 1 1 1 1 1 1 IMF6 1 1 0 1 0 1 1 1 1 Res 1 0 1 1 1 1 1 1 0 1 1 1 stands for selection, 0 represents discard.Figure 9
PACF values of original wind speed and decomposed components (95% confidence interval line).
(a)
Wind speed data setA (b)
Wind speed data setB (c)
Wind speed data setC (d)
Wind speed data setDIn order to provide better training conditions for MKLSSVM model to enhance wind speed forecasting accuracy, the inputs of each MKLSSVM model are normalized linearly according to (29) into the interval [01]. The prediction outcomes of each MKLSSVM model are denormalized through (30) by the contrary process of the corresponding normalization approach.(29)si′=si-sminsmax-smin,(30)s^i=s^i′smax-smin+smin,wheresi, smax, and smin are wind speed at time i and the maximum and the minimum wind speed, respectively. s^i′ is the output of MKLSSVM.


### 4.4. Parameter Settings

After the lag values of each input sample are determined by PACF technique, the LSSVM model based on the weighted multikernel function is ready for wind speed forecasting by the parameter optimization and feature selection. As described in the aforementioned section, the regression performance of MKLSSVM model is greatly affected by the weighted coefficient, penalty parameter, and kernel parameters. These parameters are tuned by the GSA algorithm according to the fitness function with the training wind speed data and they are set in Table5.Table 5
Parameters setup in HGSA-MKLSSVM model.
Parameters Value Maximum iteration number 100  Particle number 30  Dimension number 5  Gravitational constantG0 2  Kernel parameterγ ( 0,200 )  Kernel parameterδ ( 0,20 )  Kernel parameterr ( 0,20 )  Kernel parameterd ( 0,20 )  Weighted coefficientμ [ 0,1 ]


### 4.5. Forecasting Results, Comparison, and Analysis

#### 4.5.1. Forecasting Results and Comparison for LSSVM-Based Model

In an effort to evaluate comprehensively the proposed combined model, the comparisons and analysis of EEMD-HGSA-LSSVM based on RBF, Poly and multikernel function, EMD-HGSA-MKLSSVM, WT-HGSA-MKLSSVM, HGSA-MKLSSVM, and EEMD-MKLSSVM are given and these comparisons are divided into three parts, namely, experiments I, II, and III. The tests with the empirical wind speed samples are carried out in Matlab 2014a environment on windows 7 with 2.4GHz Intel Core i5-4210U and 64bit 8G RAM. The statistical indices and the improved percentage indices are applied to evaluate the forecasting performance.Experiment I. In this part, the proposed EEMD-HGSA-LSSVM model based on multikernel function is compared with the model based on RBF and Poly kernel functions, respectively, to show the superiority of the multikernel function. Table 6 illustrates the forecasting results obtained by EEMD-HGSA-LSSVM based on multikernel, RBF and Poly functions and the detailed improved percentage indices between the proposed model and the compared models are listed in Table 7 and the forecasting curves are shown in Figure 10. As seen from the tables, the RMSE errors of the model based on the weighted multikernel function are 0.2787 m/s, 0.2856 m/s, 0.2691 m/s, and 0.2903 m/s for wind speed data A, B, C, and D, respectively, which are smallest. The improved percentage PRMSE between the models based on the weighted multikernel function vs. based on the RBF function and the Poly function are 2.48% and 5.75% for data A, 2.79% and 5.59% for data B, 3.93% and 7.43% for data C, and 3.04% and 6.59% for data D, respectively. In the scatter diagram of the forecasting values versus the empirical measured wind speed time series displayed in Figure 10, the dashed straight red lines represent the actual original wind speed data which are the same as the forecasting results, which means that the forecasting errors are bigger if the forecasting points get further away from the line. The whole forecasting points obtained by EEMD-HGSA-MKLSSVM model are much closer to the line. From the statistical indices in Table 6 and the improved percentage in Table 7, the EEMD-HGSA-LSSVM based on the weighted multikernel function performs best although the models based on RBF function and Poly function perform slightly worse than that based on the weighted multikernel function in terms of performance indices for all wind speed data sets.Table 6
Forecasting results by EEMD-HGSA-LSSVM based on different kernel function.
Data set Models RMSE ( m / s ) MAPE ( % ) MAE ( m / s ) A Model 1 0.2858 5.6439 0.2597 Model 2 0.2957 5.8311 0.2673 Model 3 0.2787 5.457 0.2558  B Model 1 0.2938 6.8854 0.274 Model 2 0.3025 7.4624 0.2786 Model 3 0.2856 6.4077 0.2538  C Model 1 0.2801 3.433 0.2607 Model 2 0.2907 3.382 0.2546 Model 3 0.2691 3.232 0.2392  D Model 1 0.2994 4.9216 0.2819 Model 2 0.3108 4.9586 0.2865 Model 3 0.2903 4.5088 0.2607 1 Models 1, 2, and 3 denote EEMD-HGSA-LSSVM based on RBF, Poly, and multikernel function, respectively.Table 7
Improved percentage obtained by EEMD-HGSA-MKLSSVM over that based on RBF or Poly kernel function (%).
Data set Contrast P R M S E P M A P E P M A E A Model 1Vs. Model 3 2.48 3.31 1.51 Model 2Vs. Model 3 5.75 6.42 4.30  B Model 1Vs. Model 3 2.79 6.94 7.37 Model 2Vs. Model 3 5.59 14.13 8.91  C Model 1Vs. Model 3 3.93 5.85 8.25 Model 2Vs. Model 3 7.43 4.43 6.05  D Model 1Vs. Model 3 3.04 8.39 7.52 Model 2Vs. Model 3 6.59 9.07 9.01 1 Models 1, 2, and 3 denote EEMD-HGSA-LSSVM based on RBF, Poly and multikernel function, respectively.Figure 10
Wind speed forecasting results obtained by EEMD-based HGSA-LSSVM with different kernel function.
(a)
Forecasting diagram for dataA (b)
Scatter diagram and histogram for dataA (c)
Forecasting diagram for dataB (d)
Scatter diagram and histogram for dataB (e)
Forecasting diagram and histogram for dataC (f)
Scatter diagram for dataC (g)
Forecasting diagram for dataD (h)
Scatter diagram and histogram for dataDRemark. Overall, LSSVM with multikernel function exhibits the highest forecasting accuracy whereas LSSVM with RBF kernel function obtains the worst forecasting performance in that there are high fluctuation and randomness in wind speed. The combination of individual advantages of the RBF kernel function and Poly kernel function by weight coefficient can catch the nonlinear character in the wind speed.Experiment II. In this part, comparisons involving EEMD-HGSA-MKLSSVM, EMD-HGSA-MKLSSVM, WT-HGSA-MKLSSVM, and HGSA-MKLSSVM without signal preprocessing method are employed to illustrate the indispensability of the wind speed preprocessing technique in the forecasting model and superiority of EEMD approach. The optimal parameters in the HGSA-MKLSSVM without signal processing are illustrated in Table 8. The forecasting curves, forecasting indices and the detailed improved percentage indices, are shown in Figure 11 and Tables 9 and 10, respectively. As seen from the figures and tables, the proposed EEMD-based model has better forecasting performance than the corresponding WT-HGSA-MKLSSVM and EMD-HGSA-MKLSSVM models and HGSA-MKLSSVM model. In these models, HGSA-MKLSSVM without wind speed decomposition performs the worst from the statistical indices viewpoint. For example, compared with WT-based, EMD-based, and without wind speed data decomposition based HGSA-MKLSSVM, RMSE of the proposed EEMD-based model are cut by 0.1192 m/s, 0.1056 m/s, and 0.2332 m/s for data set A, 0.1312 m/s, 0.1173 m/s, and 0.255 m/s for data set B, 0.1399 m/s, 0.1191 m/s, and 0.249 m/s for data set C, 0.1122 m/s, 0.1005 m/s, and 0.2503 m/s for data set D.Table 8
Tuned parameters of MKLSSVM by HGSA directly for original wind speed.
Data set γ δ μ r d A 8.0192 0.6309 0.2935 0.9763 4.0871 B 123.7445 3.4637 0.5853 1.1434 8.3528 C 42.4292 2.3582 0.7329 1.8739 9.3567 D 33.6749 1.1526 0.4323 2.1457 7.4737Table 9
Forecasting results by EEMD-HGSA-LSSVM based on different decomposition method.
Data set Models RMSE (m/s) MAPE ( % ) MAE (m/s) A Model 4 0.5119 9.9989 0.4785 Model 5 0.3979 7.8679 0.3514 Model 6 0.3843 7.883 0.3596 Model 7 0.2787 5.457 0.2558  B Model 4 0.5406 13.3218 0.5138 Model 5 0.4168 9.9403 0.3927 Model 6 0.4029 9.658 0.3824 Model 7 0.2856 6.4077 0.2538  C Model 4 0.5181 6.5663 0.4985 Model 5 0.409 4.5731 0.3456 Model 6 0.3882 4.8105 0.3642 Model 7 0.2691 3.232 0.2392  D Model 4 0.5406 8.7821 0.5064 Model 5 0.4025 6.3819 0.3701 Model 6 0.3908 6.3498 0.3687 Model 7 0.2903 4.5088 0.2607 1 Model 4 denotes HGSA-MKLSSVM. 2 Models 5, 6, and 7 denote HGSA-MKLSSVM with WT, EMD, and EEMD, respectively.Table 10
Improved percentage obtained by EEMD-HGSA-MKLSSVM over other models with different decomposition method (%).
Data set Contrast P R M S E PMAPE PMAE A Model 4Vs. Model 7 45.56 45.42 46.54 Model 5Vs. Model 7 29.96 30.64 27.21 Model 6Vs. Model 7 27.48 30.78 28.87  B Model 4Vs. Model 7 47.17 51.90 50.60 Model 5Vs. Model 7 31.48 35.54 35.37 Model 6Vs. Model 7 29.11 33.65 33.63  C Model 4Vs. Model 7 48.06 50.78 52.02 Model 5Vs. Model 7 34.21 29.33 30.79 Model 6Vs. Model 7 30.68 32.81 34.32  D Model 4Vs. Model 7 46.30 48.66 48.52 Model 5Vs. Model 7 27.88 29.35 29.56 Model 6Vs. Model 7 25.72 28.99 29.29 1 Model 4 denotes HGSA-LSSVM. 2 Models 5, 6, and 7 denote HGSA-MKLSSVM with WT, EMD, and EEMD, respectively.Figure 11
Wind speed forecasting results obtained by HGSA-LSSVM with different signal preprocessing.
(a)
Forecasting diagram for dataA (b)
Scatter diagram for dataA (c)
Forecasting diagram for dataB (d)
Scatter diagram for dataB (e)
Forecasting diagram for dataC (f)
Scatter diagram for dataC (g)
Forecasting diagram for dataD (h)
Scatter diagram for dataDRemark. It can be obviously seen from the histograms that the forecasting accuracy of these models from the highest to the lowest are HGSA-MKLSSVM with EEMD, EMD, and WT and without signal decomposition. Compared with other models, the scatter points of the forecasting approach EEMD-HGSA-MKLSSVM distribute closest to the regression line. The preprocessing technique EEMD is more effective than WT, EMD signal decomposition when applied in the wind speed decomposition. Besides, the forecasting accuracy of the HGSA-MKLSSVM model can be improved greatly through the signal decomposition technique for all wind speed data sets; therefore, wind speed decomposition technique is indispensable in the application wind speed forecasting in that wind speed time series exhibit random and highly fluctuant, and wind speed decomposition technique makes wind speed data decomposed into relatively stable components which reduce the prediction difficulties of the forecasting engine.Experiment III. In this part, comparisons between the proposed model with EEMD-GSA-MKLSSVM and EEMD-MKLSSVM are carried out to illustrate the necessity of feature selection and parameter optimization in the forecasting model. In the EEMD-MKLSSVM model, μ=0.5 and the other kernel parameters are set according to [21]. The forecasting curves, forecasting indices and the detailed improved percentage indices, are shown in Figure 12 and Tables 11 and 12, respectively. As seen from the figures and tables, compared with EEMD-GSA-MKLSSVM and EEMD-MKLSSVM, the RMSE errors of EEMD-HGSA-MKLSSVM are reduced by 0.1155 m/s and 0.2012 m/s for data set A, 0.1113 m/s and 0.2295 m/s for data set B, 0.1227 m/s and 0.2163 m/s for data set C, 0.1026 m/s and 0.2192 m/s for data set D, respectively. From the histogram in the subplot of Figure 12, the forecasting accuracy of the models ranks from low to high as EEMD-MKLSSVM, EEMD-GSA-MKLSSVM, and EEMD-HGSA-MKLSSVM for all wind speed data sets. The most forecasting points of EEMD-MKLSSVM locate farthest from the regression line, while those of EEMD-HGSA-MKLSSVM are closest to the line.Table 11
Forecasting results by EEMD-MKLSSVM based model.
Data set Models RMSE (m/s) MAPE ( % ) MAE (m/s) A Model 8 0.4799 9.9012 0.4557 Model 9 0.3942 7.0633 0.3397 Model 10 0.2787 5.457 0.2558  B Model 8 0.5151 12.4774 0.4844 Model 9 0.3969 9.4711 0.3527 Model 10 0.2856 6.4077 0.2538  C Model 8 0.4854 5.9519 0.4591 Model 9 0.3918 4.365 0.329 Model 10 0.2691 3.232 0.2392  D Model 8 0.5095 8.5347 0.4884 Model 9 0.3929 6.479 0.3673 Model 10 0.2903 4.5088 0.2607 1 Model 8 denotes EEMD-MKLSSVM. 2 Model 9 denotes EEMD-GSA-MKLSSVM. 3 Model 10 denotes EEMD-HGSA-MKLSSVM.Table 12
Improved percentage obtained by EEMD-HGSA-MKLSSVM over other EEMD-MKLSSVM based models (%).
Data set Contrast P R M S E PMAPE PMAE A Model 8Vs. Model 10 41.93 44.89 43.87 Model 9Vs. Model 10 29.29 22.74 24.69  B Model 8Vs. Model 10 44.55 48.64 47.60 Model 9Vs. Model 10 28.04 32.34 28.04  C Model 8Vs. Model 10 44.56 45.69 47.89 Model 9Vs. Model 10 31.32 25.96 27.29  D Model 8Vs. Model 10 43.05 47.17 46.62 Model 9Vs. Model 10 26.11 30.41 29.02 1 Model 8 denotes EEMD-MKLSSVM. 2 Model 9 denotes EEMD-GSA-MKLSSVM. 3 Model 10 denotes EEMD-HGSA-MKLSSVM.Figure 12
Wind speed forecasting results obtained by EEMD-based models.
(a)
Forecasting diagram for dataA (b)
Scatter diagram for dataA (c)
Forecasting diagram for dataB (d)
Scatter diagram for dataB (e)
Forecasting diagram for dataC (f)
Scatter diagram for dataC (g)
Forecasting diagram for dataD (h)
Scatter diagram for dataDRemark. It illustrates that the application of feature selection and parameter optimization in the forecasting model is helpful to wind speed prediction accuracy; the EEMD-MKLSSVM model works worse than both EEMD-GSA-MKLSSVM model and EEMD-HGSA-MKLSSVM model, while EEMD-GSA-MKLSSVM model without feature selection performs also worse than EEMD-HGSA-MKLSSVM model in that the direct application of MKLSSVM in wind speed forecasting without parameters optimization by GSA algorithm may result in overfitting or trapping into local optima, and the redundant and illusive components within each subseries are not identified by feature selection through BGSA algorithm.


#### 4.5.2. Compared with Other Forecasting Models

In this section, Persistence method, generally adopted as a benchmark approach to validate a new developed forecasting method [1], is also taken as the benchmark approach to see how much the EEMD-HGSA-MKLSSVM approach improves the prediction performance. Moreover, the forecasting performances of the EEMD-based ELM and SVM with parameter optimization and feature selection by HGSA algorithm are compared to further evaluate the effectiveness of the proposed forecasting model in terms of the three statistical indices. In the test, the hidden nodes number in ELM is determined by the grid search(GS) method whose searching range and grid step are set as [1200] and 1, respectively. Tables 13 and 14 list forecasting statistical indices and improved percentage, respectively. As seen from the tables, compared with Persistence, EEMD-HGSA-ELM, and EEMD-HGSA-SVM, the RMSE values of the proposed model are cut by 0.4517 m/s, 0.0712 m/s, and 0.0827 m/s for data set A, 0.4664 m/s, 0.1042 m/s, and 0.1238 m/s for data set B, 0.4925 m/s, 0.0984 m/s, and 0.1148 m/s for data set C, 0.4506 m/s, 0.0844 m/s, and 0.0929 m/s for data set D, respectively.Table 13
Forecasting results by different forecasting model.
Data set Models RMSE (m/s) MAPE ( % ) MAE (m/s) A Model 11 0.7304 15.3973 0.717 Model 12 0.3499 6.2955 0.2929 Model 13 0.3614 7.5089 0.3353 Model 14 0.2787 5.457 0.2558  B Model 11 0.752 18.7638 0.7344 Model 12 0.3898 9.2387 0.3638 Model 13 0.4094 10.4421 0.3869 Model 14 0.2856 6.4077 0.2538  C Model 11 0.7616 9.9466 0.7445 Model 12 0.3675 4.642 0.3472 Model 13 0.3839 4.6887 0.355 Model 14 0.2691 3.232 0.2392  D Model 11 0.7409 12.5176 0.7193 Model 12 0.3747 6.1382 0.353 Model 13 0.3832 6.3973 0.3623 Model 14 0.2903 4.5088 0.2607 1 Model 11 denotes Persistence. 2 Model 12 denotes EEMD-HGSA-ELM. 3 Model 13 denotes EEMD-HGSA-SVM. 4 Model 14 denotes EEMD-HGSA-MKLSSVM.Table 14
Improved percentage obtained by EEMD-HGSA-MKLSSVM over other models (%).
Data set Contrast P R M S E PMAPE PMAE A Model 11Vs. Model 14 61.84 64.56 64.33 Model 12Vs. Model 14 20.34 13.32 12.67 Model 13Vs. Model 14 22.88 27.33 23.71  B Model 11Vs. Model 14 62.02 65.85 65.44 Model 12Vs. Model 14 26.73 30.64 30.24 Model 13Vs. Model 14 30.24 38.64 34.40  C Model 11Vs. Model 14 64.67 67.51 67.87 Model 12Vs. Model 14 26.78 30.37 31.11 Model 13Vs. Model 14 29.91 31.07 32.62  D Model 11Vs. Model 14 60.82 63.98 63.76 Model 12Vs. Model 14 22.52 26.55 26.15 Model 13Vs. Model 14 24.24 29.52 28.04 1 Model 11 denotes Persistence. 2 Model 12 denotes EEMD-HGSA-ELM. 3 Model 13 denotes EEMD-HGSA-SVM. 4 Model 14 denotes EEMD-HGSA-MKLSSVM.Remark. The proposed EEMD-HGSA-MKLSSVM model not only performs better than EEMD-HSA-ELM and EEMD-HGSA-SVM models, but also it obtains remarkably higher forecasting accuracy than benchmark method Persistence. In Persistence approach, the current sample at time t is utilized to predict the future time t+Δt, then the t+Δt value as the current observation is employed to forecast the next data; therefore, it can be easily established in the wind speed forecasting. Compared with Persistence approach, the developed EEMD-HGSA-MKLSSVM model is more complicated to accomplish the total prediction process; however, thanks to the advance of computer technology, this is acceptable. Moreover, the superiority of the proposed model over EEMD-HGSA-SVM can explain that the basic idea of SVM is to map the input samples into high-dimensional space through nonlinear function, while the basic working mechanism of LSSVM model is that the quadratic programming problems are converted to solve linear equations, thus enhancing its regressive performance. Therefore, the discussion and analysis by comparison with other forecasting models can provide sufficient evidence that the proposed hybrid model with feature selection and parameter optimization is an excellent approach for short-term wind speed prediction.


## 4.1. Wind Speed Data Statistical Description

In this study, the empirical wind speed data were collected and stored in a wind farm located in Anhui province of China to verify the proposed model. The wind speed time series were measured every 10 min from the anemometer at the top of a nacelle, and the wind turbines with 88m are installed at the top of the 300 m mountain (32°28′N, 118°26′E). Accurate half-an-hour to several days ahead wind speed forecasting are very beneficial to routine management of wind farms, and thus the original 10 m wind speed time series are transformed to 30 min by averaging. The developed models in the most studies are tested and evaluated with no more than two sets of wind speed data, and the conclusion that the models are better than the other approaches is generally not convincing enough [26]. Therefore, four sets of wind speed data in 2015 are randomly selected to evaluate the proposed model. The historical wind speed time series of four seasons in 2015 are displayed in the Figure 7(a). The red mark variables in Figure 7(a), which are clearly shown in Figure 7(b), are utilized to train and test the proposed models. The blue variables and red ones in Figure 7(b) are employed to train and test the models, respectively.The quantity of input samples of LSSVM model affects the wind speed forecasting performance. However, there are not any unified standard and clear definition for section of the input sample quantity [5]. As for the LSSVM which is expert in solving small samples and the suggestion made by Wang et al. [6], four sets of 576 half-hourly wind speed samples in 2015 are selected randomly from Spring, Summer, Autumn, and Winter, respectively. Thereinto, the 1st-480th sampling points are employed to train the proposed model while the subsequent 481th-576th data are applied to test the well-trained model.The descriptive statistics of the original wind speed times series list in Table1. It can be obviously seen that the empirical wind speed time series presents the characters of fluctuations and stochastic volatility, fluctuating around 7.80 m/s, 5.48 m/s, 8.47 m/s, and 5.55 m/s, respectively. There is no apparent regularity in the wind speed series.Table 1
Statistical description of empirical wind speed (m/s).
data set NO. Max Min Mean St.dev. Median Data set A All samples 576 14.59 1.22 7.80 2.79 7.91 Training data 480 14.59 1.22 8.16 2.68 7.91 Testing data 96 10.68 1.32 6.01 2.64 6.01  Data set B All samples 576 12.44 0.94 5.48 2.46 6.69 Training data 480 12.44 0.94 5.58 2.46 6.69 Testing data 96 9.91 1.34 4.99 2.38 5.62  Data set C All samples 576 15.34 0.87 8.47 3.08 8.10 Training data 480 15.34 0.87 8.56 3.28 8.10 Testing data 96 11.57 4.08 7.99 1.69 7.82  Data set D All samples 576 12.87 0.41 5.55 2.40 6.64 Training data 480 12.87 0.41 5.42 2.52 6.64 Testing data 96 9.91 2.57 6.17 1.56 6.25


## 4.2. Wind Speed Preprocessing Using EEMD

Prior to submitting the original wind speed to the forecasting engine, the signal decomposition technique EEMD is applied to break down the empirical wind speed into several relatively stable IMFs and oneRes for reducing the forecasting difficulties of HGSA-MKLSSVM. According to [8], the amplitude of the added white noise and the ensemble number in this study are set 0.2 and 100, respectively. As shown in Figure 8, the empirical wind speed time series are broken down into six IMFs and one Res. It is obviously observed from the figures that all the components present distinct respective characteristic, frequencies decrease gradually from IMF1 to Res, and high-frequency IMF1 and IMF2 reflect the stochastic characteristic of wind speed data, while the signals IMF3 to IMF6 with periodic trend features represent the periodic components of the wind speed data and Res are named the trend components.Figure 8
Wind speed decomposition using EEMD.
(a)
Wind speed dataA (b)
Wind speed dataB (c)
Wind speed dataC (d)
Wind speed dataD


## 4.3. Input-Matrix Construction by PACF and BGSA

Before the forecasting processes are carried out by the HGSA-MKLSSVM model, the input variables matrix requires to be determined. The PACF technique is employed to identify the correlations among each corresponding decomposed component for determination of the input combination. The PACF of four sets of the original wind speed and their decomposed components are illustrated in Figure9. As seen in Figure 9, the PACF values of all wind speed and decomposed components from lags 1 to 30 are calculated. For original wind speed data A, PACF values at lags bigger than 7 are between the 95% confidence lines (the upper and lower lines); thus, the time lag is determined as 7 which can be used to identify the input variables dimension, and the 7 previous continuous wind speed time series contribute the most correlative information to forecast the subsequent wind speed value. The lag order of the original wind speed data and their decomposed components are shown in Table 2. For the original wind speed data A, the input vector for the MKLSSVM model can be described as X=xt-1, xt-2, ⋯, xt-7 and the corresponding forecasting value is Y=xt-k, where k(≥7) is the forecasting horizontal. The other input variables combination for the HGSA-MKLSSVM model determined by PACF values is expressed in Table 3. The feature selection results obtained by BGSA for different subseries are shown in Table 4.Table 2
Lag order of the autoregressive process.
Samples Original IMF1 IMF2 IMF3 IMF4 IMF5 IMF6 Res A 7 11 10 15 8 12 7 8 B 8 10 13 14 11 8 7 9 C 7 11 9 10 7 13 9 6 D 6 11 10 8 8 10 9 10Table 3
Inputs combinations of wind speed and their decomposed components for the HGSA-MKLSSVM model.
Data set Time series input variable combination A Original x t - 1, xt-2, xt-3, xt-4, xt-5, xt-6, xt-7 IMF1 x t - 1, xt-2, xt-3, xt-4,xt-5, xt-6, xt-7, xt-8, xt-9, xt-10, xt-11 IMF2 x t - 1, xt-2, xt-3, xt-4, xt-5, xt-6, xt-7, xt-8, xt-9, xt-10 IMF3 x t - 1, xt-2, xt-3, xt-4, xt-5, xt-6, xt-7, xt-8, xt-9, xt-10, xt-11, xt-12, xt-13, xt-14, xt-15 IMF4 x t - 1, xt-2, xt-3, xt-4, xt-5, xt-6, xt-7, xt-8, IMF5 x t - 1, xt-2, xt-3, xt-4, xt-5, xt-6, xt-7, xt-8, xt-9, xt-10, xt-11, xt-12 IMF6 x t - 1, xt-2, xt-3, xt-4, xt-5, xt-6, xt-7 Res x t - 1, xt-2, xt-3, xt-4, xt-5, xt-6, xt-7, xt-8  B Original x t - 1, xt-2, xt-3, xt-4, xt-5, xt-6, xt-7, xt-8 IMF1 x t - 1, xt-2, xt-3, xt-4, xt-5, xt-6, xt-7, xt-8, xt-9, xt-10 IMF2 x t - 1, xt-2, xt-3, xt-4, xt-5, xt-6, xt-7, xt-8, xt-9, xt-10, xt-11, xt-12, xt-13 IMF3 x t - 1, xt-2, xt-3, xt-4, xt-5, xt-6, xt-7, xt-8, xt-9, xt-10, xt-11, xt-12, xt-13, xt-14 IMF4 x t - 1, xt-2, xt-3, xt-4, xt-5, xt-6, xt-7, xt-8, xt-9, xt-10, xt-11 IMF5 x t - 1, xt-2, xt-3, xt-4, xt-5, xt-6, xt-7, xt-8 IMF6 x t - 1, xt-2, xt-3, xt-4, xt-5, xt-6, xt-7 Res x t - 1, xt-2, xt-3, xt-4, xt-5, xt-6, xt-7, xt-8, xt-9  C Original x t - 1, xt-2, xt-3, xt-4, xt-5, xt-6, xt-7 IMF1 x t - 1, xt-2, xt-3, xt-4, xt-5, xt-6, xt-7, xt-8, xt-9, xt-10, xt-11 IMF2 x t - 1, xt-2, xt-3, xt-4, xt-5, xt-6, xt-7, xt-8, xt-9 IMF3 x t - 1, xt-2, xt-3, xt-4, xt-5, xt-6, xt-7, xt-8, xt-9, xt-10 IMF4 x t - 1, xt-2, xt-3, xt-4, xt-5, xt-6, xt-7 IMF5 x t - 1, xt-2, xt-3, xt-4, xt-5, xt-6, xt-7, xt-8, xt-9, xt-10, xt-11, xt-12, xt-13 IMF6 x t - 1, xt-2, xt-3, xt-4, xt-5, xt-6, xt-7, xt-8, xt-9, xt-8, xt-9 Res x t - 1, xt-2, xt-3, xt-4, xt-5, xt-6  D Original x t - 1, xt-2, xt-3, xt-4, xt-5, xt-6 IMF1 x t - 1, xt-2, xt-3, xt-4, xt-5, xt-6, xt-7, xt-8, xt-9, xt-10, xt-11 IMF2 x t - 1, xt-2, xt-3, xt-4, xt-5, xt-6, xt-7, xt-8, xt-9, xt-10 IMF3 x t - 1, xt-2, xt-3, xt-4, xt-5, xt-6, xt-7, xt-8 IMF4 x t - 1, xt-2, xt-3, xt-4, xt-5, xt-6, xt-7, xt-8 IMF5 x t - 1, xt-2, xt-3, xt-4, xt-5, xt-6, xt-7, xt-8, xt-9, xt-10 IMF6 x t - 1, xt-2, xt-3, xt-4, xt-5, xt-6, xt-7, xt-8, xt-9 Res x t - 1, xt-2, xt-3, xt-4, xt-5, xt-6, xt-7, xt-8, xt-9, xt-10Table 4
Feature selection results obtained by HGSA algorithm.
Data set Time series 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 A Original 1 1 1 0 1 0 1 IMF1 1 1 0 1 1 1 0 1 1 0 1 IMF2 1 1 1 1 0 0 1 1 1 1 IMF3 1 1 0 1 1 0 1 1 1 1 1 0 1 1 1 IMF4 1 1 1 0 1 1 0 1 IMF5 1 0 1 0 1 1 1 0 1 1 1 1 IMF6 1 1 0 1 1 0 1 Res 1 1 0 1 1 1 1 1  B Original 1 0 1 1 1 1 1 1 IMF1 1 0 1 1 1 1 1 1 0 1 IMF2 1 1 0 1 0 1 1 1 1 1 1 1 0 IMF3 1 0 1 1 1 1 1 1 0 0 1 1 1 1 IMF4 1 1 0 1 1 0 1 1 1 1 1 IMF5 1 1 1 0 1 1 0 1 IMF6 1 1 0 1 1 1 1 Res 1 0 1 1 1 1 0 1 1  C Original 1 1 1 0 1 1 1 IMF1 1 0 1 1 0 1 1 0 1 1 1 IMF2 1 0 1 1 1 1 0 1 1 IMF3 1 1 1 0 1 1 0 1 1 1 IMF4 1 1 0 1 1 1 1 IMF5 1 1 1 1 0 1 1 1 1 1 0 0 1 IMF6 1 1 0 1 1 1 1 0 1 Res 1 0 1 1 1 1  D Original 1 1 1 0 1 1 IMF1 1 0 1 1 1 1 0 1 1 0 1 IMF2 1 0 1 1 0 1 1 0 1 1 IMF3 1 0 1 1 1 0 1 1 IMF4 1 1 1 0 1 1 1 1 IMF5 1 1 0 1 1 1 1 1 1 1 IMF6 1 1 0 1 0 1 1 1 1 Res 1 0 1 1 1 1 1 1 0 1 1 1 stands for selection, 0 represents discard.Figure 9
PACF values of original wind speed and decomposed components (95% confidence interval line).
(a)
Wind speed data setA (b)
Wind speed data setB (c)
Wind speed data setC (d)
Wind speed data setDIn order to provide better training conditions for MKLSSVM model to enhance wind speed forecasting accuracy, the inputs of each MKLSSVM model are normalized linearly according to (29) into the interval [01]. The prediction outcomes of each MKLSSVM model are denormalized through (30) by the contrary process of the corresponding normalization approach.(29)si′=si-sminsmax-smin,(30)s^i=s^i′smax-smin+smin,wheresi, smax, and smin are wind speed at time i and the maximum and the minimum wind speed, respectively. s^i′ is the output of MKLSSVM.


## 4.4. Parameter Settings

After the lag values of each input sample are determined by PACF technique, the LSSVM model based on the weighted multikernel function is ready for wind speed forecasting by the parameter optimization and feature selection. As described in the aforementioned section, the regression performance of MKLSSVM model is greatly affected by the weighted coefficient, penalty parameter, and kernel parameters. These parameters are tuned by the GSA algorithm according to the fitness function with the training wind speed data and they are set in Table5.Table 5
Parameters setup in HGSA-MKLSSVM model.
Parameters Value Maximum iteration number 100  Particle number 30  Dimension number 5  Gravitational constantG0 2  Kernel parameterγ ( 0,200 )  Kernel parameterδ ( 0,20 )  Kernel parameterr ( 0,20 )  Kernel parameterd ( 0,20 )  Weighted coefficientμ [ 0,1 ]


## 4.5. Forecasting Results, Comparison, and Analysis

### 4.5.1. Forecasting Results and Comparison for LSSVM-Based Model

In an effort to evaluate comprehensively the proposed combined model, the comparisons and analysis of EEMD-HGSA-LSSVM based on RBF, Poly and multikernel function, EMD-HGSA-MKLSSVM, WT-HGSA-MKLSSVM, HGSA-MKLSSVM, and EEMD-MKLSSVM are given and these comparisons are divided into three parts, namely, experiments I, II, and III. The tests with the empirical wind speed samples are carried out in Matlab 2014a environment on windows 7 with 2.4GHz Intel Core i5-4210U and 64bit 8G RAM. The statistical indices and the improved percentage indices are applied to evaluate the forecasting performance.Experiment I. In this part, the proposed EEMD-HGSA-LSSVM model based on multikernel function is compared with the model based on RBF and Poly kernel functions, respectively, to show the superiority of the multikernel function. Table 6 illustrates the forecasting results obtained by EEMD-HGSA-LSSVM based on multikernel, RBF and Poly functions and the detailed improved percentage indices between the proposed model and the compared models are listed in Table 7 and the forecasting curves are shown in Figure 10. As seen from the tables, the RMSE errors of the model based on the weighted multikernel function are 0.2787 m/s, 0.2856 m/s, 0.2691 m/s, and 0.2903 m/s for wind speed data A, B, C, and D, respectively, which are smallest. The improved percentage PRMSE between the models based on the weighted multikernel function vs. based on the RBF function and the Poly function are 2.48% and 5.75% for data A, 2.79% and 5.59% for data B, 3.93% and 7.43% for data C, and 3.04% and 6.59% for data D, respectively. In the scatter diagram of the forecasting values versus the empirical measured wind speed time series displayed in Figure 10, the dashed straight red lines represent the actual original wind speed data which are the same as the forecasting results, which means that the forecasting errors are bigger if the forecasting points get further away from the line. The whole forecasting points obtained by EEMD-HGSA-MKLSSVM model are much closer to the line. From the statistical indices in Table 6 and the improved percentage in Table 7, the EEMD-HGSA-LSSVM based on the weighted multikernel function performs best although the models based on RBF function and Poly function perform slightly worse than that based on the weighted multikernel function in terms of performance indices for all wind speed data sets.Table 6
Forecasting results by EEMD-HGSA-LSSVM based on different kernel function.
Data set Models RMSE ( m / s ) MAPE ( % ) MAE ( m / s ) A Model 1 0.2858 5.6439 0.2597 Model 2 0.2957 5.8311 0.2673 Model 3 0.2787 5.457 0.2558  B Model 1 0.2938 6.8854 0.274 Model 2 0.3025 7.4624 0.2786 Model 3 0.2856 6.4077 0.2538  C Model 1 0.2801 3.433 0.2607 Model 2 0.2907 3.382 0.2546 Model 3 0.2691 3.232 0.2392  D Model 1 0.2994 4.9216 0.2819 Model 2 0.3108 4.9586 0.2865 Model 3 0.2903 4.5088 0.2607 1 Models 1, 2, and 3 denote EEMD-HGSA-LSSVM based on RBF, Poly, and multikernel function, respectively.Table 7
Improved percentage obtained by EEMD-HGSA-MKLSSVM over that based on RBF or Poly kernel function (%).
Data set Contrast P R M S E P M A P E P M A E A Model 1Vs. Model 3 2.48 3.31 1.51 Model 2Vs. Model 3 5.75 6.42 4.30  B Model 1Vs. Model 3 2.79 6.94 7.37 Model 2Vs. Model 3 5.59 14.13 8.91  C Model 1Vs. Model 3 3.93 5.85 8.25 Model 2Vs. Model 3 7.43 4.43 6.05  D Model 1Vs. Model 3 3.04 8.39 7.52 Model 2Vs. Model 3 6.59 9.07 9.01 1 Models 1, 2, and 3 denote EEMD-HGSA-LSSVM based on RBF, Poly and multikernel function, respectively.Figure 10
Wind speed forecasting results obtained by EEMD-based HGSA-LSSVM with different kernel function.
(a)
Forecasting diagram for dataA (b)
Scatter diagram and histogram for dataA (c)
Forecasting diagram for dataB (d)
Scatter diagram and histogram for dataB (e)
Forecasting diagram and histogram for dataC (f)
Scatter diagram for dataC (g)
Forecasting diagram for dataD (h)
Scatter diagram and histogram for dataDRemark. Overall, LSSVM with multikernel function exhibits the highest forecasting accuracy whereas LSSVM with RBF kernel function obtains the worst forecasting performance in that there are high fluctuation and randomness in wind speed. The combination of individual advantages of the RBF kernel function and Poly kernel function by weight coefficient can catch the nonlinear character in the wind speed.Experiment II. In this part, comparisons involving EEMD-HGSA-MKLSSVM, EMD-HGSA-MKLSSVM, WT-HGSA-MKLSSVM, and HGSA-MKLSSVM without signal preprocessing method are employed to illustrate the indispensability of the wind speed preprocessing technique in the forecasting model and superiority of EEMD approach. The optimal parameters in the HGSA-MKLSSVM without signal processing are illustrated in Table 8. The forecasting curves, forecasting indices and the detailed improved percentage indices, are shown in Figure 11 and Tables 9 and 10, respectively. As seen from the figures and tables, the proposed EEMD-based model has better forecasting performance than the corresponding WT-HGSA-MKLSSVM and EMD-HGSA-MKLSSVM models and HGSA-MKLSSVM model. In these models, HGSA-MKLSSVM without wind speed decomposition performs the worst from the statistical indices viewpoint. For example, compared with WT-based, EMD-based, and without wind speed data decomposition based HGSA-MKLSSVM, RMSE of the proposed EEMD-based model are cut by 0.1192 m/s, 0.1056 m/s, and 0.2332 m/s for data set A, 0.1312 m/s, 0.1173 m/s, and 0.255 m/s for data set B, 0.1399 m/s, 0.1191 m/s, and 0.249 m/s for data set C, 0.1122 m/s, 0.1005 m/s, and 0.2503 m/s for data set D.Table 8
Tuned parameters of MKLSSVM by HGSA directly for original wind speed.
Data set γ δ μ r d A 8.0192 0.6309 0.2935 0.9763 4.0871 B 123.7445 3.4637 0.5853 1.1434 8.3528 C 42.4292 2.3582 0.7329 1.8739 9.3567 D 33.6749 1.1526 0.4323 2.1457 7.4737Table 9
Forecasting results by EEMD-HGSA-LSSVM based on different decomposition method.
Data set Models RMSE (m/s) MAPE ( % ) MAE (m/s) A Model 4 0.5119 9.9989 0.4785 Model 5 0.3979 7.8679 0.3514 Model 6 0.3843 7.883 0.3596 Model 7 0.2787 5.457 0.2558  B Model 4 0.5406 13.3218 0.5138 Model 5 0.4168 9.9403 0.3927 Model 6 0.4029 9.658 0.3824 Model 7 0.2856 6.4077 0.2538  C Model 4 0.5181 6.5663 0.4985 Model 5 0.409 4.5731 0.3456 Model 6 0.3882 4.8105 0.3642 Model 7 0.2691 3.232 0.2392  D Model 4 0.5406 8.7821 0.5064 Model 5 0.4025 6.3819 0.3701 Model 6 0.3908 6.3498 0.3687 Model 7 0.2903 4.5088 0.2607 1 Model 4 denotes HGSA-MKLSSVM. 2 Models 5, 6, and 7 denote HGSA-MKLSSVM with WT, EMD, and EEMD, respectively.Table 10
Improved percentage obtained by EEMD-HGSA-MKLSSVM over other models with different decomposition method (%).
Data set Contrast P R M S E PMAPE PMAE A Model 4Vs. Model 7 45.56 45.42 46.54 Model 5Vs. Model 7 29.96 30.64 27.21 Model 6Vs. Model 7 27.48 30.78 28.87  B Model 4Vs. Model 7 47.17 51.90 50.60 Model 5Vs. Model 7 31.48 35.54 35.37 Model 6Vs. Model 7 29.11 33.65 33.63  C Model 4Vs. Model 7 48.06 50.78 52.02 Model 5Vs. Model 7 34.21 29.33 30.79 Model 6Vs. Model 7 30.68 32.81 34.32  D Model 4Vs. Model 7 46.30 48.66 48.52 Model 5Vs. Model 7 27.88 29.35 29.56 Model 6Vs. Model 7 25.72 28.99 29.29 1 Model 4 denotes HGSA-LSSVM. 2 Models 5, 6, and 7 denote HGSA-MKLSSVM with WT, EMD, and EEMD, respectively.Figure 11
Wind speed forecasting results obtained by HGSA-LSSVM with different signal preprocessing.
(a)
Forecasting diagram for dataA (b)
Scatter diagram for dataA (c)
Forecasting diagram for dataB (d)
Scatter diagram for dataB (e)
Forecasting diagram for dataC (f)
Scatter diagram for dataC (g)
Forecasting diagram for dataD (h)
Scatter diagram for dataDRemark. It can be obviously seen from the histograms that the forecasting accuracy of these models from the highest to the lowest are HGSA-MKLSSVM with EEMD, EMD, and WT and without signal decomposition. Compared with other models, the scatter points of the forecasting approach EEMD-HGSA-MKLSSVM distribute closest to the regression line. The preprocessing technique EEMD is more effective than WT, EMD signal decomposition when applied in the wind speed decomposition. Besides, the forecasting accuracy of the HGSA-MKLSSVM model can be improved greatly through the signal decomposition technique for all wind speed data sets; therefore, wind speed decomposition technique is indispensable in the application wind speed forecasting in that wind speed time series exhibit random and highly fluctuant, and wind speed decomposition technique makes wind speed data decomposed into relatively stable components which reduce the prediction difficulties of the forecasting engine.Experiment III. In this part, comparisons between the proposed model with EEMD-GSA-MKLSSVM and EEMD-MKLSSVM are carried out to illustrate the necessity of feature selection and parameter optimization in the forecasting model. In the EEMD-MKLSSVM model, μ=0.5 and the other kernel parameters are set according to [21]. The forecasting curves, forecasting indices and the detailed improved percentage indices, are shown in Figure 12 and Tables 11 and 12, respectively. As seen from the figures and tables, compared with EEMD-GSA-MKLSSVM and EEMD-MKLSSVM, the RMSE errors of EEMD-HGSA-MKLSSVM are reduced by 0.1155 m/s and 0.2012 m/s for data set A, 0.1113 m/s and 0.2295 m/s for data set B, 0.1227 m/s and 0.2163 m/s for data set C, 0.1026 m/s and 0.2192 m/s for data set D, respectively. From the histogram in the subplot of Figure 12, the forecasting accuracy of the models ranks from low to high as EEMD-MKLSSVM, EEMD-GSA-MKLSSVM, and EEMD-HGSA-MKLSSVM for all wind speed data sets. The most forecasting points of EEMD-MKLSSVM locate farthest from the regression line, while those of EEMD-HGSA-MKLSSVM are closest to the line.Table 11
Forecasting results by EEMD-MKLSSVM based model.
Data set Models RMSE (m/s) MAPE ( % ) MAE (m/s) A Model 8 0.4799 9.9012 0.4557 Model 9 0.3942 7.0633 0.3397 Model 10 0.2787 5.457 0.2558  B Model 8 0.5151 12.4774 0.4844 Model 9 0.3969 9.4711 0.3527 Model 10 0.2856 6.4077 0.2538  C Model 8 0.4854 5.9519 0.4591 Model 9 0.3918 4.365 0.329 Model 10 0.2691 3.232 0.2392  D Model 8 0.5095 8.5347 0.4884 Model 9 0.3929 6.479 0.3673 Model 10 0.2903 4.5088 0.2607 1 Model 8 denotes EEMD-MKLSSVM. 2 Model 9 denotes EEMD-GSA-MKLSSVM. 3 Model 10 denotes EEMD-HGSA-MKLSSVM.Table 12
Improved percentage obtained by EEMD-HGSA-MKLSSVM over other EEMD-MKLSSVM based models (%).
Data set Contrast P R M S E PMAPE PMAE A Model 8Vs. Model 10 41.93 44.89 43.87 Model 9Vs. Model 10 29.29 22.74 24.69  B Model 8Vs. Model 10 44.55 48.64 47.60 Model 9Vs. Model 10 28.04 32.34 28.04  C Model 8Vs. Model 10 44.56 45.69 47.89 Model 9Vs. Model 10 31.32 25.96 27.29  D Model 8Vs. Model 10 43.05 47.17 46.62 Model 9Vs. Model 10 26.11 30.41 29.02 1 Model 8 denotes EEMD-MKLSSVM. 2 Model 9 denotes EEMD-GSA-MKLSSVM. 3 Model 10 denotes EEMD-HGSA-MKLSSVM.Figure 12
Wind speed forecasting results obtained by EEMD-based models.
(a)
Forecasting diagram for dataA (b)
Scatter diagram for dataA (c)
Forecasting diagram for dataB (d)
Scatter diagram for dataB (e)
Forecasting diagram for dataC (f)
Scatter diagram for dataC (g)
Forecasting diagram for dataD (h)
Scatter diagram for dataDRemark. It illustrates that the application of feature selection and parameter optimization in the forecasting model is helpful to wind speed prediction accuracy; the EEMD-MKLSSVM model works worse than both EEMD-GSA-MKLSSVM model and EEMD-HGSA-MKLSSVM model, while EEMD-GSA-MKLSSVM model without feature selection performs also worse than EEMD-HGSA-MKLSSVM model in that the direct application of MKLSSVM in wind speed forecasting without parameters optimization by GSA algorithm may result in overfitting or trapping into local optima, and the redundant and illusive components within each subseries are not identified by feature selection through BGSA algorithm.


### 4.5.2. Compared with Other Forecasting Models

In this section, Persistence method, generally adopted as a benchmark approach to validate a new developed forecasting method [1], is also taken as the benchmark approach to see how much the EEMD-HGSA-MKLSSVM approach improves the prediction performance. Moreover, the forecasting performances of the EEMD-based ELM and SVM with parameter optimization and feature selection by HGSA algorithm are compared to further evaluate the effectiveness of the proposed forecasting model in terms of the three statistical indices. In the test, the hidden nodes number in ELM is determined by the grid search(GS) method whose searching range and grid step are set as [1200] and 1, respectively. Tables 13 and 14 list forecasting statistical indices and improved percentage, respectively. As seen from the tables, compared with Persistence, EEMD-HGSA-ELM, and EEMD-HGSA-SVM, the RMSE values of the proposed model are cut by 0.4517 m/s, 0.0712 m/s, and 0.0827 m/s for data set A, 0.4664 m/s, 0.1042 m/s, and 0.1238 m/s for data set B, 0.4925 m/s, 0.0984 m/s, and 0.1148 m/s for data set C, 0.4506 m/s, 0.0844 m/s, and 0.0929 m/s for data set D, respectively.Table 13
Forecasting results by different forecasting model.
Data set Models RMSE (m/s) MAPE ( % ) MAE (m/s) A Model 11 0.7304 15.3973 0.717 Model 12 0.3499 6.2955 0.2929 Model 13 0.3614 7.5089 0.3353 Model 14 0.2787 5.457 0.2558  B Model 11 0.752 18.7638 0.7344 Model 12 0.3898 9.2387 0.3638 Model 13 0.4094 10.4421 0.3869 Model 14 0.2856 6.4077 0.2538  C Model 11 0.7616 9.9466 0.7445 Model 12 0.3675 4.642 0.3472 Model 13 0.3839 4.6887 0.355 Model 14 0.2691 3.232 0.2392  D Model 11 0.7409 12.5176 0.7193 Model 12 0.3747 6.1382 0.353 Model 13 0.3832 6.3973 0.3623 Model 14 0.2903 4.5088 0.2607 1 Model 11 denotes Persistence. 2 Model 12 denotes EEMD-HGSA-ELM. 3 Model 13 denotes EEMD-HGSA-SVM. 4 Model 14 denotes EEMD-HGSA-MKLSSVM.Table 14
Improved percentage obtained by EEMD-HGSA-MKLSSVM over other models (%).
Data set Contrast P R M S E PMAPE PMAE A Model 11Vs. Model 14 61.84 64.56 64.33 Model 12Vs. Model 14 20.34 13.32 12.67 Model 13Vs. Model 14 22.88 27.33 23.71  B Model 11Vs. Model 14 62.02 65.85 65.44 Model 12Vs. Model 14 26.73 30.64 30.24 Model 13Vs. Model 14 30.24 38.64 34.40  C Model 11Vs. Model 14 64.67 67.51 67.87 Model 12Vs. Model 14 26.78 30.37 31.11 Model 13Vs. Model 14 29.91 31.07 32.62  D Model 11Vs. Model 14 60.82 63.98 63.76 Model 12Vs. Model 14 22.52 26.55 26.15 Model 13Vs. Model 14 24.24 29.52 28.04 1 Model 11 denotes Persistence. 2 Model 12 denotes EEMD-HGSA-ELM. 3 Model 13 denotes EEMD-HGSA-SVM. 4 Model 14 denotes EEMD-HGSA-MKLSSVM.Remark. The proposed EEMD-HGSA-MKLSSVM model not only performs better than EEMD-HSA-ELM and EEMD-HGSA-SVM models, but also it obtains remarkably higher forecasting accuracy than benchmark method Persistence. In Persistence approach, the current sample at time t is utilized to predict the future time t+Δt, then the t+Δt value as the current observation is employed to forecast the next data; therefore, it can be easily established in the wind speed forecasting. Compared with Persistence approach, the developed EEMD-HGSA-MKLSSVM model is more complicated to accomplish the total prediction process; however, thanks to the advance of computer technology, this is acceptable. Moreover, the superiority of the proposed model over EEMD-HGSA-SVM can explain that the basic idea of SVM is to map the input samples into high-dimensional space through nonlinear function, while the basic working mechanism of LSSVM model is that the quadratic programming problems are converted to solve linear equations, thus enhancing its regressive performance. Therefore, the discussion and analysis by comparison with other forecasting models can provide sufficient evidence that the proposed hybrid model with feature selection and parameter optimization is an excellent approach for short-term wind speed prediction.


## 4.5.1. Forecasting Results and Comparison for LSSVM-Based Model

In an effort to evaluate comprehensively the proposed combined model, the comparisons and analysis of EEMD-HGSA-LSSVM based on RBF, Poly and multikernel function, EMD-HGSA-MKLSSVM, WT-HGSA-MKLSSVM, HGSA-MKLSSVM, and EEMD-MKLSSVM are given and these comparisons are divided into three parts, namely, experiments I, II, and III. The tests with the empirical wind speed samples are carried out in Matlab 2014a environment on windows 7 with 2.4GHz Intel Core i5-4210U and 64bit 8G RAM. The statistical indices and the improved percentage indices are applied to evaluate the forecasting performance.Experiment I. In this part, the proposed EEMD-HGSA-LSSVM model based on multikernel function is compared with the model based on RBF and Poly kernel functions, respectively, to show the superiority of the multikernel function. Table 6 illustrates the forecasting results obtained by EEMD-HGSA-LSSVM based on multikernel, RBF and Poly functions and the detailed improved percentage indices between the proposed model and the compared models are listed in Table 7 and the forecasting curves are shown in Figure 10. As seen from the tables, the RMSE errors of the model based on the weighted multikernel function are 0.2787 m/s, 0.2856 m/s, 0.2691 m/s, and 0.2903 m/s for wind speed data A, B, C, and D, respectively, which are smallest. The improved percentage PRMSE between the models based on the weighted multikernel function vs. based on the RBF function and the Poly function are 2.48% and 5.75% for data A, 2.79% and 5.59% for data B, 3.93% and 7.43% for data C, and 3.04% and 6.59% for data D, respectively. In the scatter diagram of the forecasting values versus the empirical measured wind speed time series displayed in Figure 10, the dashed straight red lines represent the actual original wind speed data which are the same as the forecasting results, which means that the forecasting errors are bigger if the forecasting points get further away from the line. The whole forecasting points obtained by EEMD-HGSA-MKLSSVM model are much closer to the line. From the statistical indices in Table 6 and the improved percentage in Table 7, the EEMD-HGSA-LSSVM based on the weighted multikernel function performs best although the models based on RBF function and Poly function perform slightly worse than that based on the weighted multikernel function in terms of performance indices for all wind speed data sets.Table 6
Forecasting results by EEMD-HGSA-LSSVM based on different kernel function.
Data set Models RMSE ( m / s ) MAPE ( % ) MAE ( m / s ) A Model 1 0.2858 5.6439 0.2597 Model 2 0.2957 5.8311 0.2673 Model 3 0.2787 5.457 0.2558  B Model 1 0.2938 6.8854 0.274 Model 2 0.3025 7.4624 0.2786 Model 3 0.2856 6.4077 0.2538  C Model 1 0.2801 3.433 0.2607 Model 2 0.2907 3.382 0.2546 Model 3 0.2691 3.232 0.2392  D Model 1 0.2994 4.9216 0.2819 Model 2 0.3108 4.9586 0.2865 Model 3 0.2903 4.5088 0.2607 1 Models 1, 2, and 3 denote EEMD-HGSA-LSSVM based on RBF, Poly, and multikernel function, respectively.Table 7
Improved percentage obtained by EEMD-HGSA-MKLSSVM over that based on RBF or Poly kernel function (%).
Data set Contrast P R M S E P M A P E P M A E A Model 1Vs. Model 3 2.48 3.31 1.51 Model 2Vs. Model 3 5.75 6.42 4.30  B Model 1Vs. Model 3 2.79 6.94 7.37 Model 2Vs. Model 3 5.59 14.13 8.91  C Model 1Vs. Model 3 3.93 5.85 8.25 Model 2Vs. Model 3 7.43 4.43 6.05  D Model 1Vs. Model 3 3.04 8.39 7.52 Model 2Vs. Model 3 6.59 9.07 9.01 1 Models 1, 2, and 3 denote EEMD-HGSA-LSSVM based on RBF, Poly and multikernel function, respectively.Figure 10
Wind speed forecasting results obtained by EEMD-based HGSA-LSSVM with different kernel function.
(a)
Forecasting diagram for dataA (b)
Scatter diagram and histogram for dataA (c)
Forecasting diagram for dataB (d)
Scatter diagram and histogram for dataB (e)
Forecasting diagram and histogram for dataC (f)
Scatter diagram for dataC (g)
Forecasting diagram for dataD (h)
Scatter diagram and histogram for dataDRemark. Overall, LSSVM with multikernel function exhibits the highest forecasting accuracy whereas LSSVM with RBF kernel function obtains the worst forecasting performance in that there are high fluctuation and randomness in wind speed. The combination of individual advantages of the RBF kernel function and Poly kernel function by weight coefficient can catch the nonlinear character in the wind speed.Experiment II. In this part, comparisons involving EEMD-HGSA-MKLSSVM, EMD-HGSA-MKLSSVM, WT-HGSA-MKLSSVM, and HGSA-MKLSSVM without signal preprocessing method are employed to illustrate the indispensability of the wind speed preprocessing technique in the forecasting model and superiority of EEMD approach. The optimal parameters in the HGSA-MKLSSVM without signal processing are illustrated in Table 8. The forecasting curves, forecasting indices and the detailed improved percentage indices, are shown in Figure 11 and Tables 9 and 10, respectively. As seen from the figures and tables, the proposed EEMD-based model has better forecasting performance than the corresponding WT-HGSA-MKLSSVM and EMD-HGSA-MKLSSVM models and HGSA-MKLSSVM model. In these models, HGSA-MKLSSVM without wind speed decomposition performs the worst from the statistical indices viewpoint. For example, compared with WT-based, EMD-based, and without wind speed data decomposition based HGSA-MKLSSVM, RMSE of the proposed EEMD-based model are cut by 0.1192 m/s, 0.1056 m/s, and 0.2332 m/s for data set A, 0.1312 m/s, 0.1173 m/s, and 0.255 m/s for data set B, 0.1399 m/s, 0.1191 m/s, and 0.249 m/s for data set C, 0.1122 m/s, 0.1005 m/s, and 0.2503 m/s for data set D.Table 8
Tuned parameters of MKLSSVM by HGSA directly for original wind speed.
Data set γ δ μ r d A 8.0192 0.6309 0.2935 0.9763 4.0871 B 123.7445 3.4637 0.5853 1.1434 8.3528 C 42.4292 2.3582 0.7329 1.8739 9.3567 D 33.6749 1.1526 0.4323 2.1457 7.4737Table 9
Forecasting results by EEMD-HGSA-LSSVM based on different decomposition method.
Data set Models RMSE (m/s) MAPE ( % ) MAE (m/s) A Model 4 0.5119 9.9989 0.4785 Model 5 0.3979 7.8679 0.3514 Model 6 0.3843 7.883 0.3596 Model 7 0.2787 5.457 0.2558  B Model 4 0.5406 13.3218 0.5138 Model 5 0.4168 9.9403 0.3927 Model 6 0.4029 9.658 0.3824 Model 7 0.2856 6.4077 0.2538  C Model 4 0.5181 6.5663 0.4985 Model 5 0.409 4.5731 0.3456 Model 6 0.3882 4.8105 0.3642 Model 7 0.2691 3.232 0.2392  D Model 4 0.5406 8.7821 0.5064 Model 5 0.4025 6.3819 0.3701 Model 6 0.3908 6.3498 0.3687 Model 7 0.2903 4.5088 0.2607 1 Model 4 denotes HGSA-MKLSSVM. 2 Models 5, 6, and 7 denote HGSA-MKLSSVM with WT, EMD, and EEMD, respectively.Table 10
Improved percentage obtained by EEMD-HGSA-MKLSSVM over other models with different decomposition method (%).
Data set Contrast P R M S E PMAPE PMAE A Model 4Vs. Model 7 45.56 45.42 46.54 Model 5Vs. Model 7 29.96 30.64 27.21 Model 6Vs. Model 7 27.48 30.78 28.87  B Model 4Vs. Model 7 47.17 51.90 50.60 Model 5Vs. Model 7 31.48 35.54 35.37 Model 6Vs. Model 7 29.11 33.65 33.63  C Model 4Vs. Model 7 48.06 50.78 52.02 Model 5Vs. Model 7 34.21 29.33 30.79 Model 6Vs. Model 7 30.68 32.81 34.32  D Model 4Vs. Model 7 46.30 48.66 48.52 Model 5Vs. Model 7 27.88 29.35 29.56 Model 6Vs. Model 7 25.72 28.99 29.29 1 Model 4 denotes HGSA-LSSVM. 2 Models 5, 6, and 7 denote HGSA-MKLSSVM with WT, EMD, and EEMD, respectively.Figure 11
Wind speed forecasting results obtained by HGSA-LSSVM with different signal preprocessing.
(a)
Forecasting diagram for dataA (b)
Scatter diagram for dataA (c)
Forecasting diagram for dataB (d)
Scatter diagram for dataB (e)
Forecasting diagram for dataC (f)
Scatter diagram for dataC (g)
Forecasting diagram for dataD (h)
Scatter diagram for dataDRemark. It can be obviously seen from the histograms that the forecasting accuracy of these models from the highest to the lowest are HGSA-MKLSSVM with EEMD, EMD, and WT and without signal decomposition. Compared with other models, the scatter points of the forecasting approach EEMD-HGSA-MKLSSVM distribute closest to the regression line. The preprocessing technique EEMD is more effective than WT, EMD signal decomposition when applied in the wind speed decomposition. Besides, the forecasting accuracy of the HGSA-MKLSSVM model can be improved greatly through the signal decomposition technique for all wind speed data sets; therefore, wind speed decomposition technique is indispensable in the application wind speed forecasting in that wind speed time series exhibit random and highly fluctuant, and wind speed decomposition technique makes wind speed data decomposed into relatively stable components which reduce the prediction difficulties of the forecasting engine.Experiment III. In this part, comparisons between the proposed model with EEMD-GSA-MKLSSVM and EEMD-MKLSSVM are carried out to illustrate the necessity of feature selection and parameter optimization in the forecasting model. In the EEMD-MKLSSVM model, μ=0.5 and the other kernel parameters are set according to [21]. The forecasting curves, forecasting indices and the detailed improved percentage indices, are shown in Figure 12 and Tables 11 and 12, respectively. As seen from the figures and tables, compared with EEMD-GSA-MKLSSVM and EEMD-MKLSSVM, the RMSE errors of EEMD-HGSA-MKLSSVM are reduced by 0.1155 m/s and 0.2012 m/s for data set A, 0.1113 m/s and 0.2295 m/s for data set B, 0.1227 m/s and 0.2163 m/s for data set C, 0.1026 m/s and 0.2192 m/s for data set D, respectively. From the histogram in the subplot of Figure 12, the forecasting accuracy of the models ranks from low to high as EEMD-MKLSSVM, EEMD-GSA-MKLSSVM, and EEMD-HGSA-MKLSSVM for all wind speed data sets. The most forecasting points of EEMD-MKLSSVM locate farthest from the regression line, while those of EEMD-HGSA-MKLSSVM are closest to the line.Table 11
Forecasting results by EEMD-MKLSSVM based model.
Data set Models RMSE (m/s) MAPE ( % ) MAE (m/s) A Model 8 0.4799 9.9012 0.4557 Model 9 0.3942 7.0633 0.3397 Model 10 0.2787 5.457 0.2558  B Model 8 0.5151 12.4774 0.4844 Model 9 0.3969 9.4711 0.3527 Model 10 0.2856 6.4077 0.2538  C Model 8 0.4854 5.9519 0.4591 Model 9 0.3918 4.365 0.329 Model 10 0.2691 3.232 0.2392  D Model 8 0.5095 8.5347 0.4884 Model 9 0.3929 6.479 0.3673 Model 10 0.2903 4.5088 0.2607 1 Model 8 denotes EEMD-MKLSSVM. 2 Model 9 denotes EEMD-GSA-MKLSSVM. 3 Model 10 denotes EEMD-HGSA-MKLSSVM.Table 12
Improved percentage obtained by EEMD-HGSA-MKLSSVM over other EEMD-MKLSSVM based models (%).
Data set Contrast P R M S E PMAPE PMAE A Model 8Vs. Model 10 41.93 44.89 43.87 Model 9Vs. Model 10 29.29 22.74 24.69  B Model 8Vs. Model 10 44.55 48.64 47.60 Model 9Vs. Model 10 28.04 32.34 28.04  C Model 8Vs. Model 10 44.56 45.69 47.89 Model 9Vs. Model 10 31.32 25.96 27.29  D Model 8Vs. Model 10 43.05 47.17 46.62 Model 9Vs. Model 10 26.11 30.41 29.02 1 Model 8 denotes EEMD-MKLSSVM. 2 Model 9 denotes EEMD-GSA-MKLSSVM. 3 Model 10 denotes EEMD-HGSA-MKLSSVM.Figure 12
Wind speed forecasting results obtained by EEMD-based models.
(a)
Forecasting diagram for dataA (b)
Scatter diagram for dataA (c)
Forecasting diagram for dataB (d)
Scatter diagram for dataB (e)
Forecasting diagram for dataC (f)
Scatter diagram for dataC (g)
Forecasting diagram for dataD (h)
Scatter diagram for dataDRemark. It illustrates that the application of feature selection and parameter optimization in the forecasting model is helpful to wind speed prediction accuracy; the EEMD-MKLSSVM model works worse than both EEMD-GSA-MKLSSVM model and EEMD-HGSA-MKLSSVM model, while EEMD-GSA-MKLSSVM model without feature selection performs also worse than EEMD-HGSA-MKLSSVM model in that the direct application of MKLSSVM in wind speed forecasting without parameters optimization by GSA algorithm may result in overfitting or trapping into local optima, and the redundant and illusive components within each subseries are not identified by feature selection through BGSA algorithm.


## 4.5.2. Compared with Other Forecasting Models

In this section, Persistence method, generally adopted as a benchmark approach to validate a new developed forecasting method [1], is also taken as the benchmark approach to see how much the EEMD-HGSA-MKLSSVM approach improves the prediction performance. Moreover, the forecasting performances of the EEMD-based ELM and SVM with parameter optimization and feature selection by HGSA algorithm are compared to further evaluate the effectiveness of the proposed forecasting model in terms of the three statistical indices. In the test, the hidden nodes number in ELM is determined by the grid search(GS) method whose searching range and grid step are set as [1200] and 1, respectively. Tables 13 and 14 list forecasting statistical indices and improved percentage, respectively. As seen from the tables, compared with Persistence, EEMD-HGSA-ELM, and EEMD-HGSA-SVM, the RMSE values of the proposed model are cut by 0.4517 m/s, 0.0712 m/s, and 0.0827 m/s for data set A, 0.4664 m/s, 0.1042 m/s, and 0.1238 m/s for data set B, 0.4925 m/s, 0.0984 m/s, and 0.1148 m/s for data set C, 0.4506 m/s, 0.0844 m/s, and 0.0929 m/s for data set D, respectively.Table 13
Forecasting results by different forecasting model.
Data set Models RMSE (m/s) MAPE ( % ) MAE (m/s) A Model 11 0.7304 15.3973 0.717 Model 12 0.3499 6.2955 0.2929 Model 13 0.3614 7.5089 0.3353 Model 14 0.2787 5.457 0.2558  B Model 11 0.752 18.7638 0.7344 Model 12 0.3898 9.2387 0.3638 Model 13 0.4094 10.4421 0.3869 Model 14 0.2856 6.4077 0.2538  C Model 11 0.7616 9.9466 0.7445 Model 12 0.3675 4.642 0.3472 Model 13 0.3839 4.6887 0.355 Model 14 0.2691 3.232 0.2392  D Model 11 0.7409 12.5176 0.7193 Model 12 0.3747 6.1382 0.353 Model 13 0.3832 6.3973 0.3623 Model 14 0.2903 4.5088 0.2607 1 Model 11 denotes Persistence. 2 Model 12 denotes EEMD-HGSA-ELM. 3 Model 13 denotes EEMD-HGSA-SVM. 4 Model 14 denotes EEMD-HGSA-MKLSSVM.Table 14
Improved percentage obtained by EEMD-HGSA-MKLSSVM over other models (%).
Data set Contrast P R M S E PMAPE PMAE A Model 11Vs. Model 14 61.84 64.56 64.33 Model 12Vs. Model 14 20.34 13.32 12.67 Model 13Vs. Model 14 22.88 27.33 23.71  B Model 11Vs. Model 14 62.02 65.85 65.44 Model 12Vs. Model 14 26.73 30.64 30.24 Model 13Vs. Model 14 30.24 38.64 34.40  C Model 11Vs. Model 14 64.67 67.51 67.87 Model 12Vs. Model 14 26.78 30.37 31.11 Model 13Vs. Model 14 29.91 31.07 32.62  D Model 11Vs. Model 14 60.82 63.98 63.76 Model 12Vs. Model 14 22.52 26.55 26.15 Model 13Vs. Model 14 24.24 29.52 28.04 1 Model 11 denotes Persistence. 2 Model 12 denotes EEMD-HGSA-ELM. 3 Model 13 denotes EEMD-HGSA-SVM. 4 Model 14 denotes EEMD-HGSA-MKLSSVM.Remark. The proposed EEMD-HGSA-MKLSSVM model not only performs better than EEMD-HSA-ELM and EEMD-HGSA-SVM models, but also it obtains remarkably higher forecasting accuracy than benchmark method Persistence. In Persistence approach, the current sample at time t is utilized to predict the future time t+Δt, then the t+Δt value as the current observation is employed to forecast the next data; therefore, it can be easily established in the wind speed forecasting. Compared with Persistence approach, the developed EEMD-HGSA-MKLSSVM model is more complicated to accomplish the total prediction process; however, thanks to the advance of computer technology, this is acceptable. Moreover, the superiority of the proposed model over EEMD-HGSA-SVM can explain that the basic idea of SVM is to map the input samples into high-dimensional space through nonlinear function, while the basic working mechanism of LSSVM model is that the quadratic programming problems are converted to solve linear equations, thus enhancing its regressive performance. Therefore, the discussion and analysis by comparison with other forecasting models can provide sufficient evidence that the proposed hybrid model with feature selection and parameter optimization is an excellent approach for short-term wind speed prediction.


## 5. Conclusion

In this article, a compound MKLSSVM model optimized by HGSA algorithm integrated with signal decomposition technique EEMD, namely, EEMD-HGSA-MKLSSVM, is proposed for short-term wind speed forecasting. Four sets of mean half-hour wind speed selected randomly from the historical wind speed data in 2015 collected from a wind farm located in Anhui of China are utilized as case studies to evaluate the forecasting performance of EEMD-HGSA-MKLSSVM model. From the comparison and analysis carried out in the previous sections, some conclusions can be drawn as follows:(i)
Considering that EEMD is an effective approach to decompose and analyze the nonlinear and nonstationary signal, we adopt it as the wind speed data preprocessing tool in the hybrid models. Correspondingly, the EEMD-based forecasting model is trained and tested with the four sets of empirical wind speed data. Compared with HGSA-MKLSSVM model without signal preprocessing, the EEMD-based HGSA-MKLSSVM model has obvious improvement in the forecasting results, thus, wind speed decomposition is indispensable in the wind speed forecasting. The forecasting results show that the EEMD-based HGSA-MKLSSVM model yields better forecasting accuracy than the corresponding EMD-based and WT-based models; thus signal decomposition technique EEMD is suitable in this hybrid forecasting model.(ii)
The EEMD-HGSA-LSSVM based on multikernel function has better forecasting results than that based on RBF kernel function or Poly kernel function in that the multikernel function takes advantages of individual merits of RBF and Poly kernel functions by optimal weighted coefficient.(iii)
The hybrid algorithm HGSA utilizes the respective advantages of BGSA and RGSA to realize the feature selection and parameter optimization simultaneously. Compared with the EEMD-GSA-MKLSSVM without feature selection, the proposed model obtains smaller RMSE for all the data sets, which means that BGSA selects the useful candidates for the forecasting engine.(iv)
The proposed model outperforms the EEMD-HGSA-ELM and EEMD-HGSA-SVM. Especially, the improvements obtained by EEMD-based HGSA-MKLSSVM model over the benchmark Persistence model in terms of improved percentage RMSE are about 61.84%, 62.02%, 64.67%, and 60.82% for the corresponding wind speed data setsA, B, C, and D, respectively.Therefore, the proposed model EEMD-HGSA-MKLSSVM is an effective short-term wind speed prediction approach. For further studies, this hybrid model will be utilized for other wind farms, and some environmental and climate information should be taken into consideration as potential input samples.



---

*Source: 9287097-2018-11-14.xml*