# An Advanced Hybrid Forecasting System for Wind Speed Point Forecasting and Interval Forecasting

Authors: Haipeng Zhang; Hua Luo
Journal: Complexity
 (2020)

---

## Abstract

Ultra-short-term wind speed prediction can assist the operation and scheduling of wind turbines in the short term and further reduce the adverse effects of wind power integration. However, as wind is irregular, nonlinear, and nonstationary, to accurately predict wind speed is a difficult task. To this end, researchers have made many attempts; however, they often use only point forecasting or interval forecasting, resulting in imperfect prediction results. Therefore, in this paper, we developed a prediction system integrating an advanced data preprocessing strategy, a novel optimization model, and multiple prediction algorithms. This combined forecasting system can overcome the inherent disadvantages of the traditional forecasting methods and further improve the prediction performance. To test the effectiveness of the forecasting system, the 10-min and one-hour wind speed sequences from the Sotavento wind farm in Spain were applied for conducting comparison experiments. The results of both the interval forecasting and point forecasting indicated that, in terms of the forecasting capability and stability, the proposed system was better than the compared models. Therefore, because of the minimum prediction error and excellent generalization ability, we consider this forecasting system to be an effective tool to assist smart grid programming.


---

## Body

## 1. Introduction

Along with the rapid development of modern industry, the issues of environmental deterioration and fossil fuel depletion are becoming increasingly serious. Seeking clean and renewable energy is a hotspot topic in the world. Across the globe, the scale of wind power generation is growing rapidly and wind energy has attracted attention [1].However, when wind power is connected to the grid, the intermittence and volatility of wind power will affect the balance of power demand and supply, the security of the power system, and the power quality. Therefore, in order to ensure the stability of the wind power system and promote the large-scale integration of wind power, it is necessary to predict the short-term wind power with high accuracy [2]. For obtaining more accurate prediction results, a great deal of efforts has been made and many studies have been performed. On the basis of applied data processing models, we summarized the previous research into four classes: (a) physical class, (b) statistical class, (c) spatial correlation class, and (d) artificial intelligence class [3].Generally speaking, the physical strategy uses the atmospheric data and numerical weather simulation to predict the wind speed. Cassola and Burlando [4] used an improved numerical weather prediction (NWP) model to provide 2-year-long wind speed dataset of two wind farms in the eastern Liguria (Italy). Zhao et al. [5] put forward a novel multistep forecasting method according to the weather research and forecasting for predicting the one-day ahead wind speed and power. Dong et al. [6] proposed a new clustering analysis model for NWP information, which proved that the clustering analysis of NWP can be effectively used in the field of wind speed forecasting. Yang and Tsai [7] developed two linear regression models and a microgenetic algorithm (MGA) model to predict the gust wind speed during a typhoon. NWP needs to observe not only the conventional data of temperature, humidity, wind pressure, and visibility but also the cloud top conditions, ozone, solar radiation, sea surface temperature, and so on, and it also requires the use of supercomputers. Therefore, it cannot deal with the short-term prediction problem.Compared with the physical strategy, the statistical strategy has clear advantages in addressing short-term forecasting problems due to the full use of historical data. From the perspective of existing statistical methods, exponential smoothing (ES), ARMA, GARCH, and other typical models are universally employed and can achieve reliable prediction effects.Kavasseri and Seetharaman [8] applied ARIMA for forecasting wind speeds on the day-ahead (24 h) and two-day-ahead (48 h) horizons, which resulted in a better forecasting accuracy. Qian and Wang [9] proposed an improved seasonal gray prediction model (GM (1, 1)) to forecast the wind generation of 2020 and 2021 in China for further development of the wind power industry. The spatial correlation algorithm takes the spatial relationship between different locations’ wind speeds into consideration. For instance, Chen et al. [10] examined the application of a multifactor spatial-temporal correlation model and achieved an improvement from other benchmark models. Nevertheless, these models had to measure the wind speed from multiple spatial correlation points. Because of the strict measurement requirements and time delay, the requirements are difficult to meet.Apart from the abovementioned models, artificial intelligence arithmetic has attracted a great deal of attention since its development. Different from the previously mentioned strategies, artificial neural networks are completely dependent on historical data, which means that the models have good fault tolerance, strong robustness, and adaptability of wind speed prediction. In addition, artificial neural networks have a good ability in dealing with nonlinear data, self-adaptive adjustment, determination of prediction mode, fuzzy relationship judgments, and so on. Scholars have used artificial intelligence algorithms to forecast the speed of wind, using SVM [11, 12], ANN [13–15], methods of fuzzy logic (FL) [16, 17], and so on.Some utilizations of artificial intelligence algorithms from previous works are shown below. Guo et al. [18] proposed a new method according to BPNN mixture with seasonal exponential regulation, and this greatly improved the wind speed forecasting accuracy. Liu et al. [19] designed a novel predicting method according to the neural network of Elman (ENN) and the quadratic decomposition method, which demonstrated excellent abilities in multistep wind speed predicting.At present, artificial neural networks are always combined with other intelligent algorithms. Yang et al. [20] proposed a new forecasting system for electricity market management with multistep electricity price forecasting. A deep learning framework with two-stage feature selection is developed by Niu et al. [21] for multivariate financial time series prediction. Although the abovementioned models achieved satisfactory results, there are still some inherent defects in the neural network model, such as the high degree of fit, ease of falling into local optimization, and the slow convergence speed.However, the majority of the methods only obtain the average point of forecasting results and ignore the reflecting forecasting volatility, which makes it difficult to meet the needs and the risk management of a wind system. Therefore, many researchers have studied wind speed interval forecasting models. Uncertainty in the forecast process is an accompanying factor, which has an important impact on the decision-making process of the energy system [22]. Therefore, some papers use the interval forecast to predict this uncertainty. Based on the existing literature, the interval forecasting model can be summed up in four types of models: point forecasting error models, bootstrap models, nonparametric theory models, and quantile regression models.The quantile regression method focuses on quantum regression, and within this, each component must be modelled [23]. Bootstrap is a statistical model that applies data sampling and substitution to evaluate the robustness of the statistics, including the parameter of the confidence interval, the standard error, the coefficient correlation institution, and the regression coefficient [24]. In addition, an interval forecasting model, which applies both the ANNs and lower upper bound estimation (LUBE), was established to structure the predicting intervals for generating high-quality forecasting intervals in a short time [25].By reviewing the previous literature, we concluded that there is no single model that can currently be regarded as the best model for wind speed prediction [26]. A more intuitive comparison of all models is drawn in Table 1.Table 1
Comparison on the advantages and disadvantages of existing wind speed forecasting methods.
CategoryAdvantageDisadvantageMethod sampleSample advantageSample disadvantageNWP modelsSuitable for long-term forecastHard to deal with short-term horizons and costs a lot of computing time and resourcesComputational fluid dynamics (FLUENT, CFD)High parallel efficiencyHigh requirement of dataset length, limited simulation rangeMesoscale numerical model (MM5, WRF)Accurate simulation and a large number of alternative parameterization schemesSuitable for weather forecast, not specifically optimized for wind predictionPhysical schema (predictor)Can predict spatial propertiesExpensive calculation costStatistical strategiesWide application, less running timeAssumption difficult to achieveTime series (ARIMA, ARCH, GARCH, GM (1, 1))Simple model assumptions, good self-fitnessPoor extrapolation, limited prediction rangeClusterStrong fault tolerance, simple human-computer interactionLow accuracy, lack of systematicnessMachine learningExcellent robustness, little model risk, wide applicationHigh complexity, high requirements of knowledgeNeural network (FNN, BPNN, GRNN, ELMAN)Excellent fitting effect nonlinear propertyThe model is unstable, high degree of data dependenceArtificial intelligence (PSO, GA, BA, DA)Easy to be understood and combined with other methodsEasy to fall into local extremumData denoising (EMD, EEMD, CEEMD, SSA)Eliminate the negative impact of noise dataThe mode mixing problem in EMD and the residual noise in EEMDDeep learning (CNN, RNN, LSTM)Strong fault tolerance, simple human-computer interactionLow accuracy when the amount of data is small, lack of systematicnessReview of the previous literature shows that there is no single forecasting model that can be considered the best and applied in all cases for wind speed forecasting, as there are considerable differences among the wind speed time series.From the above analysis, a conclusion can be drawn that each method has inherent disadvantages due to the particular characteristics. Therefore, in 1969, Bates and Granger [27] first proposed the combined forecasting theory, which significantly improved the prediction performance of a single model and has been widely used up to now. Stated thus, a novel combined system is put forward in this paper for conducting wind speed prediction. The system integrates an advanced data pretreatment theory, a fresh optimization technique, and several prediction methods, namely, ARIMA, GRNN, BPNN, and BiLSTM, respectively.Specifically speaking, in contrast to the previous subjective judgment of the number of denoising layers, we performed a weighted reconstruction of each sequence obtained through decomposition to effectively control the interference of noise data. For capturing the linear and nonlinear characteristics hidden in collected wind speed sequence, BPNN, GRNN, BiLSTM, and ARIMA were selected as subsequences to forecast the pretreatment series. Then, an advanced weight-deciding technology, nondominated neighbor immune algorithm (NNIA), was effectively introduced to determine the output results.This paper makes up for the deficiencies in the previous research in the following aspects:(1)
In contrast to the previous processing ideology, this paper proposes a more effective data preprocessing theory for controlling the interference of noise data. According to the disintegration and reconstruction theory, the collected time series was disintegrated into a set of series with different frequencies and then these series were reconstructed with combination weights determined by NNIA. Through such preprocessing, the randomness and volatility of a wind speed sequence are successfully alleviated.(2)
For fully capturing all data characteristics and making up for the defects of individual methods, one linear algorithm and three neural networks were employed as submodels to build up a combination forecasting system. To expand, ARIMA performs well in seizing linear information, while BPNN, ELMNN, and BiLSTM are effective at grasping nonlinear characteristics. The combination system not only can integrate the advantages of these submodels but is also more suitable for wind speed time series with complex characteristics.(3)
To comprehensively compare the capability of different models, we established an efficient assessment system, mining the point prediction module and interval prediction module. This system contained eight assessment indicators, five of which were used to indicate point forecasting errors, while three were employed for the accuracy of interval prediction.(4)
The proposed combination system with high accuracy and stability is a powerful tool for the smart grid. On the basis of the dataset from the Sotavento wind farm in Spain, the simulation results demonstrated that relative to other benchmark models, the system was able to achieve reliable improvements.The structure of this article is as follows. Section2 provides the modelling structure and the related methodologies. In Section 3, we introduce the implementation process of the three experiments and the obtained results. And, Section 4 offers our main conclusions. The abbreviations that appear in the text are listed in Table 2.Table 2
The abbreviations that appear in this article.
NomenclatureNNIANondominated neighbor immune algorithmVMDVariational mode decompositionBi-LSTMBidirectional long short-term memoryLSTMLong short-term memoryRNNRecurrent neural networkIMFsIntrinsic mode functionsGRNNGeneralized regression neural networkBPNNBackpropagation neural networkARMAAutoregressive moving average modelESExponential smoothingEMDEmpirical mode decompositionCEEMDComplementary ensemble empirical mode decompositionSSASingular spectrum analysisENNEvolutionary neural networkELMExtreme learning machine neural networkSVMSupport vector machineNWPNumerical weather predictionGARCHGeneralized autoregressive conditional heteroskedasticityLUBELower upper bound estimationANNArtificial neural networkPSOParticle swarm optimizationGAThe genetic algorithmMAPEThe average ofN absolute percentage errorMAEThe mean absolute error ofN forecasting resultsRMSEThe square root of the average of the error squaresMSEMean squared errorDCDirectional changeFICPForecast interval coverage probabilityFINAWForecast interval normalized average widthAWDAccumulated width deviation of testing datasetWACDWeight adaptive combined denoising based models


## 2. Design of the Combined Forecasting System

This section provides the detailed introduction of the modelling structure and the related methodology.


### 2.1. Modelling Structure

Figure1 illustrates the whole process of building the combination system in this paper.Figure 1
The flowchart of the proposed forecasting system. The proposed model comprises four stages: the first stage is the self-adaption data preprocessing; the second stage is the construction of the contrastive model; the third stage is the construction of the combined system for the interval forecasting and point forecasting; and the fourth stage is the simulation of experiments.In the 1960s, Bates et al. found that by combining two or more prediction models, the prediction accuracy could be greatly improved. Since then, combinatorial theory has been widely used. We propose a novel combination wind speed prediction system, and the specific steps are expressed as follows.Step 1. Data preprocessing ideology:The collected wind speed sequence was first decomposed into a set of filters with the different frequencies through the variable mode decomposition (VMD) technique. Then, an advanced optimization model, NNIA, was employed to determine the combination weights between these filters. This weight adaptive combined denoising ideology successfully removed the adverse impacts of high frequency noise. The objective function of the adaptive denoising algorithm is as follows:(1)objective:minMAE,function:x^i=ω1∗IMFi1+ω2∗IMFi2+⋯+ωn∗IMFin,MAE=1m∗∑m=1ixi−x^i,wherex^i is the i-th data after denoising, xi is the i-th data before denoising, and IMFin is the n-th intrinsic mode function of the i-th data. The constraint for the weighting parameters ω is (0, 1).Step 2. Design of contrast models:Verifying the validity of the model is on the basis of a certain benchmark. In order to demonstrate the advantages of the developed forecasting system from multiple perspectives, we constructed three different types of contrast models, which were single prediction models, nonoptimized models with different data preprocessing methods, and hybrid models with the same data preprocessing and internal structure optimized by the same optimization algorithm.Step 3. Point forecasting and interval forecasting:For point forecasting, the proposed combined system combined four different types of forecasting models to provide high-precision point-prediction results. The objective function is as follows:(2)objective1:minMAE,objective2:minMAPE,function:y^i=ω1∗yi1+ω2∗yi2+⋯+ω4∗yi4,MAE=1n∗∑n=1iyi−y^i,MAPE=100%n∗∑n=1iyi−y^iyi,wherey^i is the ith forecasting data, yij is the ith forecasting data of jth model, and yi is the ith actual data. The constraint for the weighting parameters ω is 0,1.On the basis of the point forecasting, we utilized the fuzzy c-means clustering (FCM) algorithm and optimization algorithm to build an appropriate forecast interval. As the value of the wind speed will affect the construction of the wind speed prediction interval, FCM was used to cluster the wind speed data into 10 categories and, then, to construct the prediction interval for each category. After that, NNIA was employed to search for the global optimal ofω, which is the coefficient of the interval, and the same data have the same ω. From the optimization point of view, a forecast interval coverage probability (FICP) higher than AWD and a lower FINRW are the two objectives for high-quality interval forecasting. Therefore, the primary problem for NNIA is modelled as a constrained single-objective problem:(3)objective:minFINAWω,constraints:FINAWω>0,μ≤FICPω≤100%,whereμ is the nominal confidence level.Finally, the complete prediction interval is constructed.Step 4. Establishment of developed forecasting system:The preprocessed wind speed sequence was divided into two parts: the training set and test set. The best ratio, 4 : 1, was determined from many experiments; therefore, in the dataset with a total length of 1440, the first 1152 data were used to train the model’s structure and the remaining 288 data were input as the test set. After obtaining the fitting results of a single model, a new weight determination technique, NNIA, was employed to integrate these individual models and obtain the final prediction results.


### 2.2. The Related Methods

This part introduces the implementation process of several algorithms related to this paper.


#### 2.2.1. Data Pretreatment Methodology

For each signal, the correlation analysis signal was calculated using the Hilbert transform to obtain a unilateral frequency spectrum. The Hilbert transform process can be expressed as the Cauchy principal value (demotedp.v.) of the convolution integral [28]:(4)Hilbertft=1πp.v.∫ℝfvt−vdv.Next, the bandwidth was calculated according to the Gaussian smoothness of the demodulated signal, i.e., the square gradient. Then, the decomposition problem was transformed into a constrained variational problem:(5)min,∑K∂tδt+jπt∗μkte−jφkt22,s.t.∑Kμk=f0,where μk=μ1,…,μK and φk=φ1,…,φK stand for the abbreviations for all pattern sets and their center frequencies, respectively. And, ∑k:=∑k=1K is regarded as the sum of all modes. Through leading into the Lagrangian multiplierλt and penalty factorα, the unconstrained variational problems are transformed by constrained variational problems, and it can be expressed as [29](6)Lμk,φk,λ=α∑k∂tδt+jπt∗μke−jφkt22+ft−∑kμkt22+λt,ft−∑kμkt.By alternately updatingμk, φk, and λ, the minimum point of the extended Lagrangian expression can be obtained. The equation of μk, φk, and λ are as follows:(7)μ^kn+1φ=fφ−∑i<kμ^kn+1φ−∑i>kμ^kn+1φ+λ^nφ/21+2αφ−φkn2,(8)φkn+1=∫0∞φμ^kn+1φ2dφ∫0∞μ^kn+1φ2dφ,(9)λ^n+1φ=λ^nφ+τf^φ−∑kμ^kn+1φ.Finally, based on the idea of sequence reconstruction, the processed sequence is expressed in the following form:(10)signalRE=∑kωk∗μk=ω1∗μ1+ω2∗μ2+⋯+ωk∗μk,where the weights, w1,w2,…,wk, are determined by a multiobjective optimization algorithm, NNIA.


#### 2.2.2. Optimization Theory

To build up a complete combination system, we introduced a novel multiobjective optimization algorithm, nondominated neighbor immune algorithm (NNIA). The algorithm is used to optimize the weights in a combined system, including a data preprocessing module and different single forecasting models. The main steps of NNIA are as follows [30]:Step 1. Initialize the parameters in the model: create one initial antibody population B0 with a length of ND. When t=0, D0=ϕ, A0=ϕ, and C0==ϕ.Step 2. Update the dominated populations: determine the dominated antibodies in Bt. Copy the whole dominated antibodies to form the temporary dominated population (expressed as DTt+1). When the size of DTt+1 is smaller or equal to ND, set Dt+1=DTt+1. If not, after calculating all individuals’ crowding distances, according to these values sort the individuals in DTt+1 in descending order, and then the first ND individuals are selected to form Dt+1.Step 3. Determine whether the termination conditions are met: if t ≥ Gmax (maximum number of generations), Dt+1 is derived as the output and stops; otherwise, t = t + 1.Step 4. The selection of the nondominant neighbor: if the size of Dt is smaller or equal to NA (maximum size of active population), set At = Dt. If not, after calculating all individuals’ crowding distance, according to these values sort the individuals in Dt in descending order, and the first NA individuals are chosen for forming At.Step 5. Proportional cloning: obtain the clone population Ct through making proportional cloning to At. The proportional cloning TC on the active population A = {a1, a2, ..., a|A|} can be expressed as(11)TCa1+a2+⋯+aiqi=TCa1+TCa2+⋯+TCaiqi=a1+a12+⋯+a1qi+a21+a22+⋯+a2qi+⋯+aA1+aA2+⋯+aAqi.In this formula,TCai=ai1+ai2+⋯+aiqi, where aij=ai, i=1,2,…,A,j=1,2,…,qi, and qi stands for a self-adaptive parameter. The symbol + not only represents the arithmetical operator but also separates the antibodies here. qi = 1 means no cloning on antibody ai.The greater the crowding-distance value of the individual, the larger of theqi of the individual is. The calculation of qi is as follows:(12)qi=nC×ζai,A∑j=1Aζai,A,whereζaj,A is the crowding-distance value of the active antibodies aj and nC stands for an expectant value of the size of the clone population.Step 6. Hypermutation and recombination: carry out hypermutation and recombination on Ct and let Ct be the generated population.Step 7. Combine Ct and Dt to obtain the antibody population Bt; return to repeat Step 2.


#### 2.2.3. The Basic Theory about Submodels in a Combined System

The core forecasting algorithms involved in the combination system are ARIMA, BPNN, GRNN, and BiLSTM. In this section, their basic theories are described briefly [31].(1) ARIMA. The model is also known as ARIMA (p, d, q). The ARIMA model has three parameters: p, d, and q.p is the lag number of the data of ARIMAd is the orders that the data are differential d times in order to make the data stableq is the lag number of forecasting errorARIMA can be expressed as follows:(13)yt=∑i=1pφiyt−i+∑j=1qθjet−j+εt,where εt is the error term at time t, θj is the j-th coefficient of AR, and φiis the i-th coefficient of MA.(2) BPNN. BPNN simulates the activation and transmission of human neurons. Taking a three-layer neural network as an example, the BP neural network consists of three layers: input, hidden, and output layers. The information of each layer of neurons is transmitted to the next layer of neurons through the activation function:(14)finput=11−e−input.The BP algorithm is used to optimize the weights of the network to keep the mapping between input and output consistent. Then, BP minimizes the objective function by changing the weight of each layer:(15)Fobjective=12×Ntraining∑i=1n∑j=1doutputyij−y^ij2,where N is the number of training samples, y^ij is the forecasting value, and yij is the actual value.Sigmoid function is used as neuron transfer function. No matter how complex the network structure is, each parameter can be evaluated by calculating the gradient. This is the basic idea of multilayer perceptron backpropagation algorithm. The sigmoid function is as follows:(3) GRNN. The learning speed and nonlinear fitting capability of GRNN are strong. Even if the sample size is small, the forecasting effect of GRNN is accurate. Furthermore, the instability data can be forecasted by GRNN. Generally, the GRNN can be established by linear neurons and radial basis function neurons.GRNN consists of four layers: output layer, mode layer, sum layer, and input layer. The activation function of GRNN can be summarized as(16)pi=exp−input−inputiTinput−inputj2σ2.(4) BiLSTM. The LSTM was firstly proposed by Hochreiter and Schmidhuber [32], and they have been refined and promoted by many others. BiLSTM is a kind of special RNN network, which can learn the dependence between long sequences. BiLSTM has been widely used in function fitting.(17)Zt=gϑzXt+ζzYt−1+bz,It=σϑiXt+ζiYt−1+ψiCt−1+bi,Ft=σϑfXt+ζfYt−1+ψfCt−1+bf,Ct=Zt⋅It+Ct−1⋅Ft,ot=σϑoXt+ζoYt−1+ψoCt−1+bo,where Zt, It, and Ft are the input block, input gate, and forget gate activation, respectively. Ct is the cell state, Yt is the output of the cell at time t, ot is the activation of the output gate, and gx, σx, and hx are activation functions.


## 2.1. Modelling Structure

Figure1 illustrates the whole process of building the combination system in this paper.Figure 1
The flowchart of the proposed forecasting system. The proposed model comprises four stages: the first stage is the self-adaption data preprocessing; the second stage is the construction of the contrastive model; the third stage is the construction of the combined system for the interval forecasting and point forecasting; and the fourth stage is the simulation of experiments.In the 1960s, Bates et al. found that by combining two or more prediction models, the prediction accuracy could be greatly improved. Since then, combinatorial theory has been widely used. We propose a novel combination wind speed prediction system, and the specific steps are expressed as follows.Step 1. Data preprocessing ideology:The collected wind speed sequence was first decomposed into a set of filters with the different frequencies through the variable mode decomposition (VMD) technique. Then, an advanced optimization model, NNIA, was employed to determine the combination weights between these filters. This weight adaptive combined denoising ideology successfully removed the adverse impacts of high frequency noise. The objective function of the adaptive denoising algorithm is as follows:(1)objective:minMAE,function:x^i=ω1∗IMFi1+ω2∗IMFi2+⋯+ωn∗IMFin,MAE=1m∗∑m=1ixi−x^i,wherex^i is the i-th data after denoising, xi is the i-th data before denoising, and IMFin is the n-th intrinsic mode function of the i-th data. The constraint for the weighting parameters ω is (0, 1).Step 2. Design of contrast models:Verifying the validity of the model is on the basis of a certain benchmark. In order to demonstrate the advantages of the developed forecasting system from multiple perspectives, we constructed three different types of contrast models, which were single prediction models, nonoptimized models with different data preprocessing methods, and hybrid models with the same data preprocessing and internal structure optimized by the same optimization algorithm.Step 3. Point forecasting and interval forecasting:For point forecasting, the proposed combined system combined four different types of forecasting models to provide high-precision point-prediction results. The objective function is as follows:(2)objective1:minMAE,objective2:minMAPE,function:y^i=ω1∗yi1+ω2∗yi2+⋯+ω4∗yi4,MAE=1n∗∑n=1iyi−y^i,MAPE=100%n∗∑n=1iyi−y^iyi,wherey^i is the ith forecasting data, yij is the ith forecasting data of jth model, and yi is the ith actual data. The constraint for the weighting parameters ω is 0,1.On the basis of the point forecasting, we utilized the fuzzy c-means clustering (FCM) algorithm and optimization algorithm to build an appropriate forecast interval. As the value of the wind speed will affect the construction of the wind speed prediction interval, FCM was used to cluster the wind speed data into 10 categories and, then, to construct the prediction interval for each category. After that, NNIA was employed to search for the global optimal ofω, which is the coefficient of the interval, and the same data have the same ω. From the optimization point of view, a forecast interval coverage probability (FICP) higher than AWD and a lower FINRW are the two objectives for high-quality interval forecasting. Therefore, the primary problem for NNIA is modelled as a constrained single-objective problem:(3)objective:minFINAWω,constraints:FINAWω>0,μ≤FICPω≤100%,whereμ is the nominal confidence level.Finally, the complete prediction interval is constructed.Step 4. Establishment of developed forecasting system:The preprocessed wind speed sequence was divided into two parts: the training set and test set. The best ratio, 4 : 1, was determined from many experiments; therefore, in the dataset with a total length of 1440, the first 1152 data were used to train the model’s structure and the remaining 288 data were input as the test set. After obtaining the fitting results of a single model, a new weight determination technique, NNIA, was employed to integrate these individual models and obtain the final prediction results.


## 2.2. The Related Methods

This part introduces the implementation process of several algorithms related to this paper.


### 2.2.1. Data Pretreatment Methodology

For each signal, the correlation analysis signal was calculated using the Hilbert transform to obtain a unilateral frequency spectrum. The Hilbert transform process can be expressed as the Cauchy principal value (demotedp.v.) of the convolution integral [28]:(4)Hilbertft=1πp.v.∫ℝfvt−vdv.Next, the bandwidth was calculated according to the Gaussian smoothness of the demodulated signal, i.e., the square gradient. Then, the decomposition problem was transformed into a constrained variational problem:(5)min,∑K∂tδt+jπt∗μkte−jφkt22,s.t.∑Kμk=f0,where μk=μ1,…,μK and φk=φ1,…,φK stand for the abbreviations for all pattern sets and their center frequencies, respectively. And, ∑k:=∑k=1K is regarded as the sum of all modes. Through leading into the Lagrangian multiplierλt and penalty factorα, the unconstrained variational problems are transformed by constrained variational problems, and it can be expressed as [29](6)Lμk,φk,λ=α∑k∂tδt+jπt∗μke−jφkt22+ft−∑kμkt22+λt,ft−∑kμkt.By alternately updatingμk, φk, and λ, the minimum point of the extended Lagrangian expression can be obtained. The equation of μk, φk, and λ are as follows:(7)μ^kn+1φ=fφ−∑i<kμ^kn+1φ−∑i>kμ^kn+1φ+λ^nφ/21+2αφ−φkn2,(8)φkn+1=∫0∞φμ^kn+1φ2dφ∫0∞μ^kn+1φ2dφ,(9)λ^n+1φ=λ^nφ+τf^φ−∑kμ^kn+1φ.Finally, based on the idea of sequence reconstruction, the processed sequence is expressed in the following form:(10)signalRE=∑kωk∗μk=ω1∗μ1+ω2∗μ2+⋯+ωk∗μk,where the weights, w1,w2,…,wk, are determined by a multiobjective optimization algorithm, NNIA.


### 2.2.2. Optimization Theory

To build up a complete combination system, we introduced a novel multiobjective optimization algorithm, nondominated neighbor immune algorithm (NNIA). The algorithm is used to optimize the weights in a combined system, including a data preprocessing module and different single forecasting models. The main steps of NNIA are as follows [30]:Step 1. Initialize the parameters in the model: create one initial antibody population B0 with a length of ND. When t=0, D0=ϕ, A0=ϕ, and C0==ϕ.Step 2. Update the dominated populations: determine the dominated antibodies in Bt. Copy the whole dominated antibodies to form the temporary dominated population (expressed as DTt+1). When the size of DTt+1 is smaller or equal to ND, set Dt+1=DTt+1. If not, after calculating all individuals’ crowding distances, according to these values sort the individuals in DTt+1 in descending order, and then the first ND individuals are selected to form Dt+1.Step 3. Determine whether the termination conditions are met: if t ≥ Gmax (maximum number of generations), Dt+1 is derived as the output and stops; otherwise, t = t + 1.Step 4. The selection of the nondominant neighbor: if the size of Dt is smaller or equal to NA (maximum size of active population), set At = Dt. If not, after calculating all individuals’ crowding distance, according to these values sort the individuals in Dt in descending order, and the first NA individuals are chosen for forming At.Step 5. Proportional cloning: obtain the clone population Ct through making proportional cloning to At. The proportional cloning TC on the active population A = {a1, a2, ..., a|A|} can be expressed as(11)TCa1+a2+⋯+aiqi=TCa1+TCa2+⋯+TCaiqi=a1+a12+⋯+a1qi+a21+a22+⋯+a2qi+⋯+aA1+aA2+⋯+aAqi.In this formula,TCai=ai1+ai2+⋯+aiqi, where aij=ai, i=1,2,…,A,j=1,2,…,qi, and qi stands for a self-adaptive parameter. The symbol + not only represents the arithmetical operator but also separates the antibodies here. qi = 1 means no cloning on antibody ai.The greater the crowding-distance value of the individual, the larger of theqi of the individual is. The calculation of qi is as follows:(12)qi=nC×ζai,A∑j=1Aζai,A,whereζaj,A is the crowding-distance value of the active antibodies aj and nC stands for an expectant value of the size of the clone population.Step 6. Hypermutation and recombination: carry out hypermutation and recombination on Ct and let Ct be the generated population.Step 7. Combine Ct and Dt to obtain the antibody population Bt; return to repeat Step 2.


### 2.2.3. The Basic Theory about Submodels in a Combined System

The core forecasting algorithms involved in the combination system are ARIMA, BPNN, GRNN, and BiLSTM. In this section, their basic theories are described briefly [31].(1) ARIMA. The model is also known as ARIMA (p, d, q). The ARIMA model has three parameters: p, d, and q.p is the lag number of the data of ARIMAd is the orders that the data are differential d times in order to make the data stableq is the lag number of forecasting errorARIMA can be expressed as follows:(13)yt=∑i=1pφiyt−i+∑j=1qθjet−j+εt,where εt is the error term at time t, θj is the j-th coefficient of AR, and φiis the i-th coefficient of MA.(2) BPNN. BPNN simulates the activation and transmission of human neurons. Taking a three-layer neural network as an example, the BP neural network consists of three layers: input, hidden, and output layers. The information of each layer of neurons is transmitted to the next layer of neurons through the activation function:(14)finput=11−e−input.The BP algorithm is used to optimize the weights of the network to keep the mapping between input and output consistent. Then, BP minimizes the objective function by changing the weight of each layer:(15)Fobjective=12×Ntraining∑i=1n∑j=1doutputyij−y^ij2,where N is the number of training samples, y^ij is the forecasting value, and yij is the actual value.Sigmoid function is used as neuron transfer function. No matter how complex the network structure is, each parameter can be evaluated by calculating the gradient. This is the basic idea of multilayer perceptron backpropagation algorithm. The sigmoid function is as follows:(3) GRNN. The learning speed and nonlinear fitting capability of GRNN are strong. Even if the sample size is small, the forecasting effect of GRNN is accurate. Furthermore, the instability data can be forecasted by GRNN. Generally, the GRNN can be established by linear neurons and radial basis function neurons.GRNN consists of four layers: output layer, mode layer, sum layer, and input layer. The activation function of GRNN can be summarized as(16)pi=exp−input−inputiTinput−inputj2σ2.(4) BiLSTM. The LSTM was firstly proposed by Hochreiter and Schmidhuber [32], and they have been refined and promoted by many others. BiLSTM is a kind of special RNN network, which can learn the dependence between long sequences. BiLSTM has been widely used in function fitting.(17)Zt=gϑzXt+ζzYt−1+bz,It=σϑiXt+ζiYt−1+ψiCt−1+bi,Ft=σϑfXt+ζfYt−1+ψfCt−1+bf,Ct=Zt⋅It+Ct−1⋅Ft,ot=σϑoXt+ζoYt−1+ψoCt−1+bo,where Zt, It, and Ft are the input block, input gate, and forget gate activation, respectively. Ct is the cell state, Yt is the output of the cell at time t, ot is the activation of the output gate, and gx, σx, and hx are activation functions.


## 2.2.1. Data Pretreatment Methodology

For each signal, the correlation analysis signal was calculated using the Hilbert transform to obtain a unilateral frequency spectrum. The Hilbert transform process can be expressed as the Cauchy principal value (demotedp.v.) of the convolution integral [28]:(4)Hilbertft=1πp.v.∫ℝfvt−vdv.Next, the bandwidth was calculated according to the Gaussian smoothness of the demodulated signal, i.e., the square gradient. Then, the decomposition problem was transformed into a constrained variational problem:(5)min,∑K∂tδt+jπt∗μkte−jφkt22,s.t.∑Kμk=f0,where μk=μ1,…,μK and φk=φ1,…,φK stand for the abbreviations for all pattern sets and their center frequencies, respectively. And, ∑k:=∑k=1K is regarded as the sum of all modes. Through leading into the Lagrangian multiplierλt and penalty factorα, the unconstrained variational problems are transformed by constrained variational problems, and it can be expressed as [29](6)Lμk,φk,λ=α∑k∂tδt+jπt∗μke−jφkt22+ft−∑kμkt22+λt,ft−∑kμkt.By alternately updatingμk, φk, and λ, the minimum point of the extended Lagrangian expression can be obtained. The equation of μk, φk, and λ are as follows:(7)μ^kn+1φ=fφ−∑i<kμ^kn+1φ−∑i>kμ^kn+1φ+λ^nφ/21+2αφ−φkn2,(8)φkn+1=∫0∞φμ^kn+1φ2dφ∫0∞μ^kn+1φ2dφ,(9)λ^n+1φ=λ^nφ+τf^φ−∑kμ^kn+1φ.Finally, based on the idea of sequence reconstruction, the processed sequence is expressed in the following form:(10)signalRE=∑kωk∗μk=ω1∗μ1+ω2∗μ2+⋯+ωk∗μk,where the weights, w1,w2,…,wk, are determined by a multiobjective optimization algorithm, NNIA.


## 2.2.2. Optimization Theory

To build up a complete combination system, we introduced a novel multiobjective optimization algorithm, nondominated neighbor immune algorithm (NNIA). The algorithm is used to optimize the weights in a combined system, including a data preprocessing module and different single forecasting models. The main steps of NNIA are as follows [30]:Step 1. Initialize the parameters in the model: create one initial antibody population B0 with a length of ND. When t=0, D0=ϕ, A0=ϕ, and C0==ϕ.Step 2. Update the dominated populations: determine the dominated antibodies in Bt. Copy the whole dominated antibodies to form the temporary dominated population (expressed as DTt+1). When the size of DTt+1 is smaller or equal to ND, set Dt+1=DTt+1. If not, after calculating all individuals’ crowding distances, according to these values sort the individuals in DTt+1 in descending order, and then the first ND individuals are selected to form Dt+1.Step 3. Determine whether the termination conditions are met: if t ≥ Gmax (maximum number of generations), Dt+1 is derived as the output and stops; otherwise, t = t + 1.Step 4. The selection of the nondominant neighbor: if the size of Dt is smaller or equal to NA (maximum size of active population), set At = Dt. If not, after calculating all individuals’ crowding distance, according to these values sort the individuals in Dt in descending order, and the first NA individuals are chosen for forming At.Step 5. Proportional cloning: obtain the clone population Ct through making proportional cloning to At. The proportional cloning TC on the active population A = {a1, a2, ..., a|A|} can be expressed as(11)TCa1+a2+⋯+aiqi=TCa1+TCa2+⋯+TCaiqi=a1+a12+⋯+a1qi+a21+a22+⋯+a2qi+⋯+aA1+aA2+⋯+aAqi.In this formula,TCai=ai1+ai2+⋯+aiqi, where aij=ai, i=1,2,…,A,j=1,2,…,qi, and qi stands for a self-adaptive parameter. The symbol + not only represents the arithmetical operator but also separates the antibodies here. qi = 1 means no cloning on antibody ai.The greater the crowding-distance value of the individual, the larger of theqi of the individual is. The calculation of qi is as follows:(12)qi=nC×ζai,A∑j=1Aζai,A,whereζaj,A is the crowding-distance value of the active antibodies aj and nC stands for an expectant value of the size of the clone population.Step 6. Hypermutation and recombination: carry out hypermutation and recombination on Ct and let Ct be the generated population.Step 7. Combine Ct and Dt to obtain the antibody population Bt; return to repeat Step 2.


## 2.2.3. The Basic Theory about Submodels in a Combined System

The core forecasting algorithms involved in the combination system are ARIMA, BPNN, GRNN, and BiLSTM. In this section, their basic theories are described briefly [31].(1) ARIMA. The model is also known as ARIMA (p, d, q). The ARIMA model has three parameters: p, d, and q.p is the lag number of the data of ARIMAd is the orders that the data are differential d times in order to make the data stableq is the lag number of forecasting errorARIMA can be expressed as follows:(13)yt=∑i=1pφiyt−i+∑j=1qθjet−j+εt,where εt is the error term at time t, θj is the j-th coefficient of AR, and φiis the i-th coefficient of MA.(2) BPNN. BPNN simulates the activation and transmission of human neurons. Taking a three-layer neural network as an example, the BP neural network consists of three layers: input, hidden, and output layers. The information of each layer of neurons is transmitted to the next layer of neurons through the activation function:(14)finput=11−e−input.The BP algorithm is used to optimize the weights of the network to keep the mapping between input and output consistent. Then, BP minimizes the objective function by changing the weight of each layer:(15)Fobjective=12×Ntraining∑i=1n∑j=1doutputyij−y^ij2,where N is the number of training samples, y^ij is the forecasting value, and yij is the actual value.Sigmoid function is used as neuron transfer function. No matter how complex the network structure is, each parameter can be evaluated by calculating the gradient. This is the basic idea of multilayer perceptron backpropagation algorithm. The sigmoid function is as follows:(3) GRNN. The learning speed and nonlinear fitting capability of GRNN are strong. Even if the sample size is small, the forecasting effect of GRNN is accurate. Furthermore, the instability data can be forecasted by GRNN. Generally, the GRNN can be established by linear neurons and radial basis function neurons.GRNN consists of four layers: output layer, mode layer, sum layer, and input layer. The activation function of GRNN can be summarized as(16)pi=exp−input−inputiTinput−inputj2σ2.(4) BiLSTM. The LSTM was firstly proposed by Hochreiter and Schmidhuber [32], and they have been refined and promoted by many others. BiLSTM is a kind of special RNN network, which can learn the dependence between long sequences. BiLSTM has been widely used in function fitting.(17)Zt=gϑzXt+ζzYt−1+bz,It=σϑiXt+ζiYt−1+ψiCt−1+bi,Ft=σϑfXt+ζfYt−1+ψfCt−1+bf,Ct=Zt⋅It+Ct−1⋅Ft,ot=σϑoXt+ζoYt−1+ψoCt−1+bo,where Zt, It, and Ft are the input block, input gate, and forget gate activation, respectively. Ct is the cell state, Yt is the output of the cell at time t, ot is the activation of the output gate, and gx, σx, and hx are activation functions.


## 3. Experiments and Analysis

We adopted three experiments and the associated analyses to verify the forecasting capacity of the proposed model. The data sources, evaluation metrics, and the experiments are also introduced in this section.


### 3.1. Study Area Description and Datasets

Original wind speed sequences were collected from Sotavento, which is located in the southwest of Europe, in Galicia, Spain. Sotavento is highlighted as an example of good practice by UNESCO in 2014. The time intervals were 10 minutes and one hour. The preprocessed wind speed sequence was divided into two parts: the training set and test set. The best ratio, 4 : 1, was determined from many experiments; therefore, in the dataset with a total length of 1440, the first 1152 data were used to train the model’s structure and the remaining 288 data were input as the test set. Dataset 1 and dataset 2 contained 10-minute intervals; dataset 3 contained one-hour intervals. The data structure of the three sites is listed in Table3.Table 3
The data structure of three sites.
SitesTrain testSample sizeData structureMaxMinMeanSth.Site 1All148815.70.55.41352.3068Train112515.70.55.38252.3377Test28813.31.215.41352.3068Site 2All148823.20.47.39214.2494Train112523.21.08.26384.3627Test2889.280.44.69012.3223Site 3All148821.10.35.30962.8843Train112521.10.55.60303.1751Test28812.10.94.41372.0628


### 3.2. Evaluation Metrics

After obtaining the prediction results, there must be a clear and unified standard to evaluate the prediction ability of the model. Therefore, we established a comprehensive evaluation system composed of seven indexes and evaluated the prediction ability from two aspects of point prediction performance and interval precision. In this paper, MAE, MAPE, MSE, RMSE, and DC were selected to reflect the prediction error. The FICP, FINAW, and AWD represent the interval prediction capability. Table4 shows the specific information regarding these indexes.Table 4
The indexes for evaluating the prediction ability.
MetricDefinitionEquationMAPEThe average ofN absolute percentage errorMSE=∑i=1NAi−Fi/A×100%/NMAPE=∑i=1NAi−Fi/Ai×100%/NMAEThe mean absolute error ofN forecasting resultsMSE=∑i=1NAi−Fi/NRMSEThe square root of the average of the error squaresRMSE=∑i=1NAi−Fi2/NRMSE=∑i=1NAi−Fi2/NMSEMean squared errorMSE=∑i=1Nyi−y^i2/NDCDirectional changeDC=100/N−1∑t=1N−1at,at=0,otherwise,1,ifyt+1−yty^t+1−yt>0,FICPForecast interval coverage probabilityFICP=∑i=1Nci×100%/NFINAWForecast interval normalized average widthFINAW=∑i=1NUi−Li×100%/N∗maxdata−mindataAWDAccumulated width deviation of testing datasetAWD=∑i=1NAWDi/NIn the formulas,Li and Ui represent the lower and the upper limit of the production interval, respectively. And, ci is the number of the truth value contained in constructed interval. N is the number of the testing set [16]. Data are the actual data. Besides, AWDi is the width deviation of the construction interval of each sample.(18)AWDi=forecast−UiUi−LiorLi−forecastUi−Li.


### 3.3. Experiment Design

According to the obtained wind speed time series, three experiments were conducted to make clear the comparison between the proposed forecasting system and other models. The input was the wind speed value of the first four days, and the output was the wind speed value of the fifth day.


#### 3.3.1. The Related Parameters of the MATLAB Operation of the Proposed Forecasting Model

In this part, the relevant parameters of each model used in this paper are discussed to achieve the best prediction effect.In the VMD method, there are six parameters that must be set to achieve a better denoising effect. Alpha is the balancing parameter of the data fidelity constraint, and we set it to 2000; we chose 0 for the noise slack (tau); and the number of modes K to be recovered was 10 in this paper. We did not require the DC part; therefore, the DC was 0. All omegas started at 0; therefore, init = 0. The tolerance of the convergence criterion (tol) was 1e − 7.As for NNIA, the population size, analog binary cross parameters, polynomial variation parameter, frequency minimum, frequency maximum, and iterations were the parameters that were required to be set. A small range of changes in one of the parameters has little effect on the accuracy; however, the program run time is quite different. Thus, considering the two factors, the default parameters were set up as follows: the population size was 10, analog binary cross parameter was 0.25, polynomial variation parameter was 0.5, frequency minimum was 0, frequency maximum was 2, and iterations were 100. The weights of IMFs of dataset 1 calculated by NNIA of the Sotavento wind farm were 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 0.79, 0.53, 0.08, and 0.00.The parameters of BiLSTM were as follows: the learning method was deep learning; the maximum number of training epochs was 300; the minimum batch size was 120; the count of hidden layers in BiLSTM unit one was 250; the count of hidden layers in BiLSTM unit two was 250; the number of hidden layers in the feedforward network was 20–50; and the execution environment was CPU.


#### 3.3.2. Experiment I: Comparison with Other Weight Adaptive Combined Denoising (WACD)-Based Models

In this experiment, a group of hybrid models that were also preprocessed by weight adaptive combined denoising (WACD) was selected as benchmark models. The experimental results are exhibited in Figure2 and Table 5. The following results can be drawn.Figure 2
Experiment I: comparison with other WACD-based models.Table 5
Prediction performance comparison table of the proposed model and hybrid models using the same data pretreatment technique.
DatasetModelOne stepTwo stepsMAEMAPERMSEMSEDCFICPFINAWAWDMAEMAPERMSEMSEDCFICPFINAWAWDSet 1WACD-ARIMA0.16623.94330.22960.052790.785089.45580.10170.03120.546912.99430.73210.535948.630153.58360.10740.3382WACD-BPNN0.19594.73820.27160.073785.665587.07480.10700.04120.40739.87690.55090.303556.164461.77470.10490.2207WACD-GRNN0.39309.52450.51610.266363.481255.10200.08990.26750.528812.80150.69250.479659.246642.66210.08360.5223WACD-BILSTM0.31657.04250.41500.172269.965965.98640.09460.14630.33127.99950.40880.167160.068366.66670.10600.0923Proposed model0.16483.83540.22840.049892.013993.40280.09840.01930.24355.63170.33290.110985.958983.27650.12840.0494Set 2WACD-ARIMA0.22737.38850.31260.097788.395990.13610.13460.02550.753224.17221.07961.165546.917851.19450.13730.2955WACD-BPNN0.314912.36390.41250.170280.546180.95240.12470.04800.652926.44950.89180.795456.164452.90100.12550.3144WACD-GRNN0.511217.35990.75410.568753.242369.72790.14180.14900.793932.03721.07431.154254.452131.74060.08491.4376WACD-BILSTM0.427212.85020.58000.336370.648572.10880.12760.09580.433514.73000.59620.355565.870372.10880.12820.1073Proposed model0.21596.949870.30650.093790.277894.21770.12600.01280.397411.67540.54770.299980.479581.57000.19870.0533Set 3WACD-ARIMA0.30218.77530.38280.146587.713384.69390.19040.02510.891125.18291.14311.306750.342549.14680.22160.2982WACD-BPNN0.31119.19700.39620.157087.713384.69390.19510.01960.549515.84000.73130.534775.684970.30720.21830.1445WACD-GRNN0.705918.95320.89340.798162.798646.93880.16500.27380.921525.72011.14581.312957.534233.10580.15150.5364WACD-BILSTM0.408710.85510.51950.269882.252674.48980.17580.06650.453311.83730.57210.327276.791873.12920.18720.0900Proposed model0.29468.26570.37940.145689.930594.55780.18600.00720.36559.68880.47360.224385.616488.73720.23490.0188The most notable results are bolded.WACD is weight adaptive combined denoising. MAPE is the average of absolute percentage error. MAE is the mean absolute error of forecasting results. RMSE is the square root of the average of the error squares. MSE is mean squared error. DC is directional change. FICP is forecast interval coverage probability. FINAW is forecast interval normalized average width. AWD is accumulated width deviation of testing dataset.In both one-step and two-step forecasting for dataset 1, the developed system obtained a remarkable performance. For example, the developed system achieved the smallest MAE values at 0.1648 and 0.2435, which were better than the other hybrid models. As for the DC, the correct directional change was 92.0139%, which proved that the forecasting of trend of the proposed model was very accurate. For interval prediction, FICP and FINAW cannot be analyzed separately. If the FICP is large but the interval normalized average width is bigger at the same time, the model is bad. According to Table5, when the FINAW is similar to other models, the FICP of the proposed model is bigger than the WACD-ARIMA, WACD-BPNN, WACD-GRNN, and WACD-BILSTM. In addition, the accumulated width deviation of the proposed model was 0.0193, which is smaller than the other models. Dataset 2 was similar to dataset 1. Therefore, we concluded that the forecasting effect of the developed model was better than other models in both point forecasting and interval forecasting for 10-minute wind speed data.For the 1-hour wind speed data (dataset 3), the WACD-NNIA-EBGB system still provided more distinguished prediction performance, not only in one step but also in two steps. For instance, for two-step prediction, our proposed system with the smallest errors was the most satisfied and it could precisely fit the changing trend of 89.9305% points, whereas the most nonideal error indicator values (MAPE at 18.9532%, MAE at 0.7059, RMSE at 0.8934, and MSE at 0.7981) were obtained by the WACD-GRNN model in one-step forecasting.Remark 1.
According to Figure2, compared with the mixed model, there is clear evidence that the developed prediction system with the lowest statistical error values in both point and interval prediction had the optimal prediction performance. In addition, other assessment indicator values were also the most ideal, which fully proves the feasibility of this combination system in short-term wind speed prediction.


#### 3.3.3. Experiment II: Comparison with Models with Different Data Preprocessing Technologies

The aim of experiment II was to make clear the comparison situation between several common data preprocessing technologies, namely, EMD, CEEMD, SSA, and VMD and the developed pretreatment strategy. The assessment indicator values are listed in Table6 and Figure 3. The prediction results are shown as Figure 4.Table 6
Prediction performance comparison table of the proposed model and combined models using different data pretreatment techniques.
DatasetModelOne stepTwo stepsMAEMAPERMSEMSEDCFICPFINAWAWDMAEMAPERMSEMSEDCFICPFINAWAWDSet 1EMD-based model0.23885.47280.31470.099180.204880.95240.11640.05380.33517.67730.44860.201279.794569.96590.11700.1303CEEMD-based model0.23965.73450.31340.098280.887481.63270.11200.04720.40779.62140.53110.282175.342560.75090.11330.1796SSA-based model0.39339.10760.50150.251560.750960.54420.13310.14250.39589.17050.50240.252460.958960.40960.13360.1449VMD-based model0.23795.50630.32250.104083.276582.65310.12170.05480.513212.06690.68440.468468.150752.90100.13420.2736Proposed model0.16483.83540.22840.049892.013993.40280.09840.01930.24355.63170.33290.110985.958983.27650.12840.0494Set 2EMD-based model0.359410.35800.55600.309176.109283.67350.19770.04860.520915.21060.73100.534371.575369.96590.23420.0990CEEMD-based model0.329510.42100.46970.220678.839688.09520.18100.03350.596918.60630.84200.709068.493267.57680.18750.1524SSA-based model0.503014.57620.71260.507860.409670.40820.22080.09830.512014.96080.74700.558071.232969.28330.22160.1105VMD-based model0.371511.05450.52390.274579.522284.01360.18060.05500.771523.41331.07001.145062.671251.87710.20050.2378Proposed model0.21596.949870.30650.093790.277894.21770.12600.01280.397411.67540.54770.299980.479581.57000.19870.0533Set 3EMD-based model0.446212.37840.56370.317776.450577.55100.20710.04270.596715.97070.74180.550374.315168.94200.20250.0716CEEMD-based model0.400410.94720.51950.269981.570083.33330.20140.02970.688518.80340.86870.754676.369960.75090.20840.1217SSA-based model1.002826.96961.24041.538756.314044.21770.33720.25721.011526.85551.25341.571156.849344.36860.35500.2524VMD-based model0.32658.88090.41500.172287.713391.15650.21080.01410.934925.24131.18511.404465.068548.12290.25130.2248Proposed model0.29468.26570.37940.145689.930594.55780.18600.00720.36559.68880.47360.224385.616488.73720.23490.0188The most notable results are bolded.MAPE is the average of absolute percentage error. MAE is the mean absolute error of forecasting results. RMSE is the square root of the average of the error squares. MSE is mean squared error. DC is directional change. FICP is forecast interval coverage probability. FINAW is forecast interval normalized average width. AWD is accumulated width deviation of testing dataset.Figure 3
Experiment II: comparison with models employing diverse data pretreatment techniques.Figure 4
Experiment III: comparison with classic individual models.For Set 1, the WACD-based combination system was more successful in processing noise data, as observed from all indicator values. In addition, this strategy obtained higher DC values, meaning that after removing the interference, the forecasting models better fit the changes of wind speed. The denoising effect of EMD technology followed the developed system, with MAPE values of 5.4278%, and 7.6773%, respectively. Moreover, when the number of prediction steps varied, the prediction ability of CEEMD-based and VMD-based combined models undulated greatly. When the FINAW of all models was almost the same, the forecast interval coverage probability of the proposed model was approximately 10% more than the other models in one-step forecasting. The AWD of the proposed model was 0.0345 smaller than the EMD-based model.For Set 2, the results showed that the prediction errors obtained by the WACD-based combined system under 7% were smaller than the other data that were pretreatment-based. Taking one-step prediction as example, the developed system’s RMSE was 0.3065, improving 30.40% and 54.12%, respectively, relative to the values of the CEEMD-based (0.4697) and SSA-based prediction systems (0.7126).In Set 3, for both one-step and multistep prediction, any evaluation index of the prediction system put forward in this paper was at the optimal level. More specifically, under the one-step forecasting situation, the WACD-based system had the best MAE, MAPE, RMSE, MSE, and DC at 0.2946, 8.2657%, 0.3794, 0.1456, and 89.9305, respectively. Nevertheless, the prediction performance of the SSA-based systems was the worst with the biggest MAPE at 26.9696%.Remark 2.
Different from the recombination idea of the denoising layer number based on subjective judgment, each subsequence obtained though decomposition was weighted and reconstructed to eliminate the adverse influence of high frequency noise. The results of experiment II show that the combination system based on WACD was always better than the other preprocessing technologies, which indicates that the idea of adaptive combination reconstruction of weight has great potential in the future.


#### 3.3.4. Experiment III: Comparison with Traditional Forecasting Models

This experiment highlighted the contrast that exists between the prediction capability of the developed system and several traditional prediction methods, including ARIMA, BPNN, ELMNN, GRNN, and ELMAN. The detailed numerical values are shown in Table7 and Figure 4.Table 7
Prediction performance comparison table of the proposed model and classic individual models.
DatasetModelOne stepTwo stepsMAEMAPERMSEMSEDCFICPFINAWAWDMAEMAPERMSEMSEDCFICPFINAWAWDSet 1ARIMA0.40119.55130.53440.285653.924963.94560.10420.19690.646215.73360.82020.672753.082241.97950.10420.4227BPNN0.40659.73950.54670.298953.924956.12240.10420.24330.539513.40910.71080.505255.137047.09900.09780.4435ELMAN0.537812.97080.68280.466251.535844.21770.10020.41160.612514.82760.77400.599051.712341.97950.10520.4647GRNN0.455911.44200.59190.350455.290146.93880.09410.35360.569814.40800.73850.545456.849338.22530.08870.6553ELMNN0.39769.62030.53140.282459.385759.52380.10150.21350.539113.38730.69520.483351.027444.36860.09670.4563Proposed model0.16483.83540.22840.049892.013993.40280.09840.01930.24355.63170.33290.110985.958983.27650.12840.0494Set 2ARIMA0.510416.00710.74490.554953.924971.76870.13480.15300.836327.62481.19281.422852.054844.02730.11690.4541BPNN0.636026.65980.84680.717154.607548.63950.12210.32090.745027.58241.07281.150953.082247.44030.12650.4329ELMAN0.697922.72320.99210.984353.583650.00000.11750.32130.799426.13531.13221.281953.424746.07510.11780.4054GRNN0.672027.83850.92800.861254.607540.81630.08960.69600.786430.35581.09221.192855.137034.12970.09120.9722ELMNN0.494816.08980.70330.494668.600761.22450.10690.21640.751628.84631.06301.129954.794542.66210.12070.4272Proposed model0.21596.949870.30650.093790.277894.21770.12600.01280.397411.67540.54770.299980.479581.57000.19870.0533Set 3ARIMA0.707620.49920.91310.833859.385759.52380.21550.21661.153733.80701.46142.135853.767137.88400.19780.4784BPNN0.714920.34670.91880.844256.314057.48300.20290.22350.987028.49731.24631.553253.767139.24910.18310.4745ELMAN0.961127.70131.22501.500556.655342.85710.18530.41711.301937.02471.64732.713757.191840.61430.24490.4282GRNN0.765823.33170.98230.964858.703148.97960.17730.31371.004431.02811.25581.577157.191833.78840.16290.6429ELMNN0.691220.53300.90170.813161.092260.20410.19830.22720.973129.53521.23521.525754.109640.27300.18630.4719Proposed model0.29468.26570.37940.145689.930594.55780.18600.00720.36559.68880.47360.224385.616488.73720.23490.0188The most notable results are bolded.MAPE is the average of absolute percentage error. MAE is the mean absolute error of forecasting results. RMSE is the square root of the average of the error squares. MSE is mean squared error. DC is directional change. FICP is forecast interval coverage probability. FINAW is forecast interval normalized average width. AWD is accumulated width deviation of testing dataset.For Set 1, although the assessment index values were changed with the number of prediction steps, the point prediction and interval prediction performance of the proposed combination system was always optimal among all models. This phenomenon reveals that the combined system put forward in this study not only had the smallest prediction error but also achieved the highest stability.For Set 2, the RMSE values of ARIMA, BPNN, ELMAN, GRNN, and ELM were 1.1928, 1.0728, 1.1324, 1.0922, and 1.0630, respectively, in two-step forecasting. Apart from the worst ELMAN model, there was no significant difference among the prediction results obtained from the other single models; however, they were significantly weaker than the proposed composite system, whose RMSE was 0.5477, even smaller than 50% of the above model.In Set 3, when comparing the interval forecasting capability, the developed system was still the best relative to the other comparison methods in terms of the values of the predictive performance indexes. The values of FICP and FINAW should be studied together. For instance, with the two-step prediction, the developed system had the smallest FINAW and the biggest FICP, which were almost two-times greater than the single models.Remark 3.
Compared with the single models, the advantages of the combined strategy were outstanding. This conclusion has nothing to do with the selected evaluation index and the number of prediction steps. Therefore, these traditional individual prediction models are no longer suitable for current wind speed prediction.


## 3.1. Study Area Description and Datasets

Original wind speed sequences were collected from Sotavento, which is located in the southwest of Europe, in Galicia, Spain. Sotavento is highlighted as an example of good practice by UNESCO in 2014. The time intervals were 10 minutes and one hour. The preprocessed wind speed sequence was divided into two parts: the training set and test set. The best ratio, 4 : 1, was determined from many experiments; therefore, in the dataset with a total length of 1440, the first 1152 data were used to train the model’s structure and the remaining 288 data were input as the test set. Dataset 1 and dataset 2 contained 10-minute intervals; dataset 3 contained one-hour intervals. The data structure of the three sites is listed in Table3.Table 3
The data structure of three sites.
SitesTrain testSample sizeData structureMaxMinMeanSth.Site 1All148815.70.55.41352.3068Train112515.70.55.38252.3377Test28813.31.215.41352.3068Site 2All148823.20.47.39214.2494Train112523.21.08.26384.3627Test2889.280.44.69012.3223Site 3All148821.10.35.30962.8843Train112521.10.55.60303.1751Test28812.10.94.41372.0628


## 3.2. Evaluation Metrics

After obtaining the prediction results, there must be a clear and unified standard to evaluate the prediction ability of the model. Therefore, we established a comprehensive evaluation system composed of seven indexes and evaluated the prediction ability from two aspects of point prediction performance and interval precision. In this paper, MAE, MAPE, MSE, RMSE, and DC were selected to reflect the prediction error. The FICP, FINAW, and AWD represent the interval prediction capability. Table4 shows the specific information regarding these indexes.Table 4
The indexes for evaluating the prediction ability.
MetricDefinitionEquationMAPEThe average ofN absolute percentage errorMSE=∑i=1NAi−Fi/A×100%/NMAPE=∑i=1NAi−Fi/Ai×100%/NMAEThe mean absolute error ofN forecasting resultsMSE=∑i=1NAi−Fi/NRMSEThe square root of the average of the error squaresRMSE=∑i=1NAi−Fi2/NRMSE=∑i=1NAi−Fi2/NMSEMean squared errorMSE=∑i=1Nyi−y^i2/NDCDirectional changeDC=100/N−1∑t=1N−1at,at=0,otherwise,1,ifyt+1−yty^t+1−yt>0,FICPForecast interval coverage probabilityFICP=∑i=1Nci×100%/NFINAWForecast interval normalized average widthFINAW=∑i=1NUi−Li×100%/N∗maxdata−mindataAWDAccumulated width deviation of testing datasetAWD=∑i=1NAWDi/NIn the formulas,Li and Ui represent the lower and the upper limit of the production interval, respectively. And, ci is the number of the truth value contained in constructed interval. N is the number of the testing set [16]. Data are the actual data. Besides, AWDi is the width deviation of the construction interval of each sample.(18)AWDi=forecast−UiUi−LiorLi−forecastUi−Li.


## 3.3. Experiment Design

According to the obtained wind speed time series, three experiments were conducted to make clear the comparison between the proposed forecasting system and other models. The input was the wind speed value of the first four days, and the output was the wind speed value of the fifth day.


### 3.3.1. The Related Parameters of the MATLAB Operation of the Proposed Forecasting Model

In this part, the relevant parameters of each model used in this paper are discussed to achieve the best prediction effect.In the VMD method, there are six parameters that must be set to achieve a better denoising effect. Alpha is the balancing parameter of the data fidelity constraint, and we set it to 2000; we chose 0 for the noise slack (tau); and the number of modes K to be recovered was 10 in this paper. We did not require the DC part; therefore, the DC was 0. All omegas started at 0; therefore, init = 0. The tolerance of the convergence criterion (tol) was 1e − 7.As for NNIA, the population size, analog binary cross parameters, polynomial variation parameter, frequency minimum, frequency maximum, and iterations were the parameters that were required to be set. A small range of changes in one of the parameters has little effect on the accuracy; however, the program run time is quite different. Thus, considering the two factors, the default parameters were set up as follows: the population size was 10, analog binary cross parameter was 0.25, polynomial variation parameter was 0.5, frequency minimum was 0, frequency maximum was 2, and iterations were 100. The weights of IMFs of dataset 1 calculated by NNIA of the Sotavento wind farm were 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 0.79, 0.53, 0.08, and 0.00.The parameters of BiLSTM were as follows: the learning method was deep learning; the maximum number of training epochs was 300; the minimum batch size was 120; the count of hidden layers in BiLSTM unit one was 250; the count of hidden layers in BiLSTM unit two was 250; the number of hidden layers in the feedforward network was 20–50; and the execution environment was CPU.


### 3.3.2. Experiment I: Comparison with Other Weight Adaptive Combined Denoising (WACD)-Based Models

In this experiment, a group of hybrid models that were also preprocessed by weight adaptive combined denoising (WACD) was selected as benchmark models. The experimental results are exhibited in Figure2 and Table 5. The following results can be drawn.Figure 2
Experiment I: comparison with other WACD-based models.Table 5
Prediction performance comparison table of the proposed model and hybrid models using the same data pretreatment technique.
DatasetModelOne stepTwo stepsMAEMAPERMSEMSEDCFICPFINAWAWDMAEMAPERMSEMSEDCFICPFINAWAWDSet 1WACD-ARIMA0.16623.94330.22960.052790.785089.45580.10170.03120.546912.99430.73210.535948.630153.58360.10740.3382WACD-BPNN0.19594.73820.27160.073785.665587.07480.10700.04120.40739.87690.55090.303556.164461.77470.10490.2207WACD-GRNN0.39309.52450.51610.266363.481255.10200.08990.26750.528812.80150.69250.479659.246642.66210.08360.5223WACD-BILSTM0.31657.04250.41500.172269.965965.98640.09460.14630.33127.99950.40880.167160.068366.66670.10600.0923Proposed model0.16483.83540.22840.049892.013993.40280.09840.01930.24355.63170.33290.110985.958983.27650.12840.0494Set 2WACD-ARIMA0.22737.38850.31260.097788.395990.13610.13460.02550.753224.17221.07961.165546.917851.19450.13730.2955WACD-BPNN0.314912.36390.41250.170280.546180.95240.12470.04800.652926.44950.89180.795456.164452.90100.12550.3144WACD-GRNN0.511217.35990.75410.568753.242369.72790.14180.14900.793932.03721.07431.154254.452131.74060.08491.4376WACD-BILSTM0.427212.85020.58000.336370.648572.10880.12760.09580.433514.73000.59620.355565.870372.10880.12820.1073Proposed model0.21596.949870.30650.093790.277894.21770.12600.01280.397411.67540.54770.299980.479581.57000.19870.0533Set 3WACD-ARIMA0.30218.77530.38280.146587.713384.69390.19040.02510.891125.18291.14311.306750.342549.14680.22160.2982WACD-BPNN0.31119.19700.39620.157087.713384.69390.19510.01960.549515.84000.73130.534775.684970.30720.21830.1445WACD-GRNN0.705918.95320.89340.798162.798646.93880.16500.27380.921525.72011.14581.312957.534233.10580.15150.5364WACD-BILSTM0.408710.85510.51950.269882.252674.48980.17580.06650.453311.83730.57210.327276.791873.12920.18720.0900Proposed model0.29468.26570.37940.145689.930594.55780.18600.00720.36559.68880.47360.224385.616488.73720.23490.0188The most notable results are bolded.WACD is weight adaptive combined denoising. MAPE is the average of absolute percentage error. MAE is the mean absolute error of forecasting results. RMSE is the square root of the average of the error squares. MSE is mean squared error. DC is directional change. FICP is forecast interval coverage probability. FINAW is forecast interval normalized average width. AWD is accumulated width deviation of testing dataset.In both one-step and two-step forecasting for dataset 1, the developed system obtained a remarkable performance. For example, the developed system achieved the smallest MAE values at 0.1648 and 0.2435, which were better than the other hybrid models. As for the DC, the correct directional change was 92.0139%, which proved that the forecasting of trend of the proposed model was very accurate. For interval prediction, FICP and FINAW cannot be analyzed separately. If the FICP is large but the interval normalized average width is bigger at the same time, the model is bad. According to Table5, when the FINAW is similar to other models, the FICP of the proposed model is bigger than the WACD-ARIMA, WACD-BPNN, WACD-GRNN, and WACD-BILSTM. In addition, the accumulated width deviation of the proposed model was 0.0193, which is smaller than the other models. Dataset 2 was similar to dataset 1. Therefore, we concluded that the forecasting effect of the developed model was better than other models in both point forecasting and interval forecasting for 10-minute wind speed data.For the 1-hour wind speed data (dataset 3), the WACD-NNIA-EBGB system still provided more distinguished prediction performance, not only in one step but also in two steps. For instance, for two-step prediction, our proposed system with the smallest errors was the most satisfied and it could precisely fit the changing trend of 89.9305% points, whereas the most nonideal error indicator values (MAPE at 18.9532%, MAE at 0.7059, RMSE at 0.8934, and MSE at 0.7981) were obtained by the WACD-GRNN model in one-step forecasting.Remark 1.
According to Figure2, compared with the mixed model, there is clear evidence that the developed prediction system with the lowest statistical error values in both point and interval prediction had the optimal prediction performance. In addition, other assessment indicator values were also the most ideal, which fully proves the feasibility of this combination system in short-term wind speed prediction.


### 3.3.3. Experiment II: Comparison with Models with Different Data Preprocessing Technologies

The aim of experiment II was to make clear the comparison situation between several common data preprocessing technologies, namely, EMD, CEEMD, SSA, and VMD and the developed pretreatment strategy. The assessment indicator values are listed in Table6 and Figure 3. The prediction results are shown as Figure 4.Table 6
Prediction performance comparison table of the proposed model and combined models using different data pretreatment techniques.
DatasetModelOne stepTwo stepsMAEMAPERMSEMSEDCFICPFINAWAWDMAEMAPERMSEMSEDCFICPFINAWAWDSet 1EMD-based model0.23885.47280.31470.099180.204880.95240.11640.05380.33517.67730.44860.201279.794569.96590.11700.1303CEEMD-based model0.23965.73450.31340.098280.887481.63270.11200.04720.40779.62140.53110.282175.342560.75090.11330.1796SSA-based model0.39339.10760.50150.251560.750960.54420.13310.14250.39589.17050.50240.252460.958960.40960.13360.1449VMD-based model0.23795.50630.32250.104083.276582.65310.12170.05480.513212.06690.68440.468468.150752.90100.13420.2736Proposed model0.16483.83540.22840.049892.013993.40280.09840.01930.24355.63170.33290.110985.958983.27650.12840.0494Set 2EMD-based model0.359410.35800.55600.309176.109283.67350.19770.04860.520915.21060.73100.534371.575369.96590.23420.0990CEEMD-based model0.329510.42100.46970.220678.839688.09520.18100.03350.596918.60630.84200.709068.493267.57680.18750.1524SSA-based model0.503014.57620.71260.507860.409670.40820.22080.09830.512014.96080.74700.558071.232969.28330.22160.1105VMD-based model0.371511.05450.52390.274579.522284.01360.18060.05500.771523.41331.07001.145062.671251.87710.20050.2378Proposed model0.21596.949870.30650.093790.277894.21770.12600.01280.397411.67540.54770.299980.479581.57000.19870.0533Set 3EMD-based model0.446212.37840.56370.317776.450577.55100.20710.04270.596715.97070.74180.550374.315168.94200.20250.0716CEEMD-based model0.400410.94720.51950.269981.570083.33330.20140.02970.688518.80340.86870.754676.369960.75090.20840.1217SSA-based model1.002826.96961.24041.538756.314044.21770.33720.25721.011526.85551.25341.571156.849344.36860.35500.2524VMD-based model0.32658.88090.41500.172287.713391.15650.21080.01410.934925.24131.18511.404465.068548.12290.25130.2248Proposed model0.29468.26570.37940.145689.930594.55780.18600.00720.36559.68880.47360.224385.616488.73720.23490.0188The most notable results are bolded.MAPE is the average of absolute percentage error. MAE is the mean absolute error of forecasting results. RMSE is the square root of the average of the error squares. MSE is mean squared error. DC is directional change. FICP is forecast interval coverage probability. FINAW is forecast interval normalized average width. AWD is accumulated width deviation of testing dataset.Figure 3
Experiment II: comparison with models employing diverse data pretreatment techniques.Figure 4
Experiment III: comparison with classic individual models.For Set 1, the WACD-based combination system was more successful in processing noise data, as observed from all indicator values. In addition, this strategy obtained higher DC values, meaning that after removing the interference, the forecasting models better fit the changes of wind speed. The denoising effect of EMD technology followed the developed system, with MAPE values of 5.4278%, and 7.6773%, respectively. Moreover, when the number of prediction steps varied, the prediction ability of CEEMD-based and VMD-based combined models undulated greatly. When the FINAW of all models was almost the same, the forecast interval coverage probability of the proposed model was approximately 10% more than the other models in one-step forecasting. The AWD of the proposed model was 0.0345 smaller than the EMD-based model.For Set 2, the results showed that the prediction errors obtained by the WACD-based combined system under 7% were smaller than the other data that were pretreatment-based. Taking one-step prediction as example, the developed system’s RMSE was 0.3065, improving 30.40% and 54.12%, respectively, relative to the values of the CEEMD-based (0.4697) and SSA-based prediction systems (0.7126).In Set 3, for both one-step and multistep prediction, any evaluation index of the prediction system put forward in this paper was at the optimal level. More specifically, under the one-step forecasting situation, the WACD-based system had the best MAE, MAPE, RMSE, MSE, and DC at 0.2946, 8.2657%, 0.3794, 0.1456, and 89.9305, respectively. Nevertheless, the prediction performance of the SSA-based systems was the worst with the biggest MAPE at 26.9696%.Remark 2.
Different from the recombination idea of the denoising layer number based on subjective judgment, each subsequence obtained though decomposition was weighted and reconstructed to eliminate the adverse influence of high frequency noise. The results of experiment II show that the combination system based on WACD was always better than the other preprocessing technologies, which indicates that the idea of adaptive combination reconstruction of weight has great potential in the future.


### 3.3.4. Experiment III: Comparison with Traditional Forecasting Models

This experiment highlighted the contrast that exists between the prediction capability of the developed system and several traditional prediction methods, including ARIMA, BPNN, ELMNN, GRNN, and ELMAN. The detailed numerical values are shown in Table7 and Figure 4.Table 7
Prediction performance comparison table of the proposed model and classic individual models.
DatasetModelOne stepTwo stepsMAEMAPERMSEMSEDCFICPFINAWAWDMAEMAPERMSEMSEDCFICPFINAWAWDSet 1ARIMA0.40119.55130.53440.285653.924963.94560.10420.19690.646215.73360.82020.672753.082241.97950.10420.4227BPNN0.40659.73950.54670.298953.924956.12240.10420.24330.539513.40910.71080.505255.137047.09900.09780.4435ELMAN0.537812.97080.68280.466251.535844.21770.10020.41160.612514.82760.77400.599051.712341.97950.10520.4647GRNN0.455911.44200.59190.350455.290146.93880.09410.35360.569814.40800.73850.545456.849338.22530.08870.6553ELMNN0.39769.62030.53140.282459.385759.52380.10150.21350.539113.38730.69520.483351.027444.36860.09670.4563Proposed model0.16483.83540.22840.049892.013993.40280.09840.01930.24355.63170.33290.110985.958983.27650.12840.0494Set 2ARIMA0.510416.00710.74490.554953.924971.76870.13480.15300.836327.62481.19281.422852.054844.02730.11690.4541BPNN0.636026.65980.84680.717154.607548.63950.12210.32090.745027.58241.07281.150953.082247.44030.12650.4329ELMAN0.697922.72320.99210.984353.583650.00000.11750.32130.799426.13531.13221.281953.424746.07510.11780.4054GRNN0.672027.83850.92800.861254.607540.81630.08960.69600.786430.35581.09221.192855.137034.12970.09120.9722ELMNN0.494816.08980.70330.494668.600761.22450.10690.21640.751628.84631.06301.129954.794542.66210.12070.4272Proposed model0.21596.949870.30650.093790.277894.21770.12600.01280.397411.67540.54770.299980.479581.57000.19870.0533Set 3ARIMA0.707620.49920.91310.833859.385759.52380.21550.21661.153733.80701.46142.135853.767137.88400.19780.4784BPNN0.714920.34670.91880.844256.314057.48300.20290.22350.987028.49731.24631.553253.767139.24910.18310.4745ELMAN0.961127.70131.22501.500556.655342.85710.18530.41711.301937.02471.64732.713757.191840.61430.24490.4282GRNN0.765823.33170.98230.964858.703148.97960.17730.31371.004431.02811.25581.577157.191833.78840.16290.6429ELMNN0.691220.53300.90170.813161.092260.20410.19830.22720.973129.53521.23521.525754.109640.27300.18630.4719Proposed model0.29468.26570.37940.145689.930594.55780.18600.00720.36559.68880.47360.224385.616488.73720.23490.0188The most notable results are bolded.MAPE is the average of absolute percentage error. MAE is the mean absolute error of forecasting results. RMSE is the square root of the average of the error squares. MSE is mean squared error. DC is directional change. FICP is forecast interval coverage probability. FINAW is forecast interval normalized average width. AWD is accumulated width deviation of testing dataset.For Set 1, although the assessment index values were changed with the number of prediction steps, the point prediction and interval prediction performance of the proposed combination system was always optimal among all models. This phenomenon reveals that the combined system put forward in this study not only had the smallest prediction error but also achieved the highest stability.For Set 2, the RMSE values of ARIMA, BPNN, ELMAN, GRNN, and ELM were 1.1928, 1.0728, 1.1324, 1.0922, and 1.0630, respectively, in two-step forecasting. Apart from the worst ELMAN model, there was no significant difference among the prediction results obtained from the other single models; however, they were significantly weaker than the proposed composite system, whose RMSE was 0.5477, even smaller than 50% of the above model.In Set 3, when comparing the interval forecasting capability, the developed system was still the best relative to the other comparison methods in terms of the values of the predictive performance indexes. The values of FICP and FINAW should be studied together. For instance, with the two-step prediction, the developed system had the smallest FINAW and the biggest FICP, which were almost two-times greater than the single models.Remark 3.
Compared with the single models, the advantages of the combined strategy were outstanding. This conclusion has nothing to do with the selected evaluation index and the number of prediction steps. Therefore, these traditional individual prediction models are no longer suitable for current wind speed prediction.


## 3.3.1. The Related Parameters of the MATLAB Operation of the Proposed Forecasting Model

In this part, the relevant parameters of each model used in this paper are discussed to achieve the best prediction effect.In the VMD method, there are six parameters that must be set to achieve a better denoising effect. Alpha is the balancing parameter of the data fidelity constraint, and we set it to 2000; we chose 0 for the noise slack (tau); and the number of modes K to be recovered was 10 in this paper. We did not require the DC part; therefore, the DC was 0. All omegas started at 0; therefore, init = 0. The tolerance of the convergence criterion (tol) was 1e − 7.As for NNIA, the population size, analog binary cross parameters, polynomial variation parameter, frequency minimum, frequency maximum, and iterations were the parameters that were required to be set. A small range of changes in one of the parameters has little effect on the accuracy; however, the program run time is quite different. Thus, considering the two factors, the default parameters were set up as follows: the population size was 10, analog binary cross parameter was 0.25, polynomial variation parameter was 0.5, frequency minimum was 0, frequency maximum was 2, and iterations were 100. The weights of IMFs of dataset 1 calculated by NNIA of the Sotavento wind farm were 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 0.79, 0.53, 0.08, and 0.00.The parameters of BiLSTM were as follows: the learning method was deep learning; the maximum number of training epochs was 300; the minimum batch size was 120; the count of hidden layers in BiLSTM unit one was 250; the count of hidden layers in BiLSTM unit two was 250; the number of hidden layers in the feedforward network was 20–50; and the execution environment was CPU.


## 3.3.2. Experiment I: Comparison with Other Weight Adaptive Combined Denoising (WACD)-Based Models

In this experiment, a group of hybrid models that were also preprocessed by weight adaptive combined denoising (WACD) was selected as benchmark models. The experimental results are exhibited in Figure2 and Table 5. The following results can be drawn.Figure 2
Experiment I: comparison with other WACD-based models.Table 5
Prediction performance comparison table of the proposed model and hybrid models using the same data pretreatment technique.
DatasetModelOne stepTwo stepsMAEMAPERMSEMSEDCFICPFINAWAWDMAEMAPERMSEMSEDCFICPFINAWAWDSet 1WACD-ARIMA0.16623.94330.22960.052790.785089.45580.10170.03120.546912.99430.73210.535948.630153.58360.10740.3382WACD-BPNN0.19594.73820.27160.073785.665587.07480.10700.04120.40739.87690.55090.303556.164461.77470.10490.2207WACD-GRNN0.39309.52450.51610.266363.481255.10200.08990.26750.528812.80150.69250.479659.246642.66210.08360.5223WACD-BILSTM0.31657.04250.41500.172269.965965.98640.09460.14630.33127.99950.40880.167160.068366.66670.10600.0923Proposed model0.16483.83540.22840.049892.013993.40280.09840.01930.24355.63170.33290.110985.958983.27650.12840.0494Set 2WACD-ARIMA0.22737.38850.31260.097788.395990.13610.13460.02550.753224.17221.07961.165546.917851.19450.13730.2955WACD-BPNN0.314912.36390.41250.170280.546180.95240.12470.04800.652926.44950.89180.795456.164452.90100.12550.3144WACD-GRNN0.511217.35990.75410.568753.242369.72790.14180.14900.793932.03721.07431.154254.452131.74060.08491.4376WACD-BILSTM0.427212.85020.58000.336370.648572.10880.12760.09580.433514.73000.59620.355565.870372.10880.12820.1073Proposed model0.21596.949870.30650.093790.277894.21770.12600.01280.397411.67540.54770.299980.479581.57000.19870.0533Set 3WACD-ARIMA0.30218.77530.38280.146587.713384.69390.19040.02510.891125.18291.14311.306750.342549.14680.22160.2982WACD-BPNN0.31119.19700.39620.157087.713384.69390.19510.01960.549515.84000.73130.534775.684970.30720.21830.1445WACD-GRNN0.705918.95320.89340.798162.798646.93880.16500.27380.921525.72011.14581.312957.534233.10580.15150.5364WACD-BILSTM0.408710.85510.51950.269882.252674.48980.17580.06650.453311.83730.57210.327276.791873.12920.18720.0900Proposed model0.29468.26570.37940.145689.930594.55780.18600.00720.36559.68880.47360.224385.616488.73720.23490.0188The most notable results are bolded.WACD is weight adaptive combined denoising. MAPE is the average of absolute percentage error. MAE is the mean absolute error of forecasting results. RMSE is the square root of the average of the error squares. MSE is mean squared error. DC is directional change. FICP is forecast interval coverage probability. FINAW is forecast interval normalized average width. AWD is accumulated width deviation of testing dataset.In both one-step and two-step forecasting for dataset 1, the developed system obtained a remarkable performance. For example, the developed system achieved the smallest MAE values at 0.1648 and 0.2435, which were better than the other hybrid models. As for the DC, the correct directional change was 92.0139%, which proved that the forecasting of trend of the proposed model was very accurate. For interval prediction, FICP and FINAW cannot be analyzed separately. If the FICP is large but the interval normalized average width is bigger at the same time, the model is bad. According to Table5, when the FINAW is similar to other models, the FICP of the proposed model is bigger than the WACD-ARIMA, WACD-BPNN, WACD-GRNN, and WACD-BILSTM. In addition, the accumulated width deviation of the proposed model was 0.0193, which is smaller than the other models. Dataset 2 was similar to dataset 1. Therefore, we concluded that the forecasting effect of the developed model was better than other models in both point forecasting and interval forecasting for 10-minute wind speed data.For the 1-hour wind speed data (dataset 3), the WACD-NNIA-EBGB system still provided more distinguished prediction performance, not only in one step but also in two steps. For instance, for two-step prediction, our proposed system with the smallest errors was the most satisfied and it could precisely fit the changing trend of 89.9305% points, whereas the most nonideal error indicator values (MAPE at 18.9532%, MAE at 0.7059, RMSE at 0.8934, and MSE at 0.7981) were obtained by the WACD-GRNN model in one-step forecasting.Remark 1.
According to Figure2, compared with the mixed model, there is clear evidence that the developed prediction system with the lowest statistical error values in both point and interval prediction had the optimal prediction performance. In addition, other assessment indicator values were also the most ideal, which fully proves the feasibility of this combination system in short-term wind speed prediction.


## 3.3.3. Experiment II: Comparison with Models with Different Data Preprocessing Technologies

The aim of experiment II was to make clear the comparison situation between several common data preprocessing technologies, namely, EMD, CEEMD, SSA, and VMD and the developed pretreatment strategy. The assessment indicator values are listed in Table6 and Figure 3. The prediction results are shown as Figure 4.Table 6
Prediction performance comparison table of the proposed model and combined models using different data pretreatment techniques.
DatasetModelOne stepTwo stepsMAEMAPERMSEMSEDCFICPFINAWAWDMAEMAPERMSEMSEDCFICPFINAWAWDSet 1EMD-based model0.23885.47280.31470.099180.204880.95240.11640.05380.33517.67730.44860.201279.794569.96590.11700.1303CEEMD-based model0.23965.73450.31340.098280.887481.63270.11200.04720.40779.62140.53110.282175.342560.75090.11330.1796SSA-based model0.39339.10760.50150.251560.750960.54420.13310.14250.39589.17050.50240.252460.958960.40960.13360.1449VMD-based model0.23795.50630.32250.104083.276582.65310.12170.05480.513212.06690.68440.468468.150752.90100.13420.2736Proposed model0.16483.83540.22840.049892.013993.40280.09840.01930.24355.63170.33290.110985.958983.27650.12840.0494Set 2EMD-based model0.359410.35800.55600.309176.109283.67350.19770.04860.520915.21060.73100.534371.575369.96590.23420.0990CEEMD-based model0.329510.42100.46970.220678.839688.09520.18100.03350.596918.60630.84200.709068.493267.57680.18750.1524SSA-based model0.503014.57620.71260.507860.409670.40820.22080.09830.512014.96080.74700.558071.232969.28330.22160.1105VMD-based model0.371511.05450.52390.274579.522284.01360.18060.05500.771523.41331.07001.145062.671251.87710.20050.2378Proposed model0.21596.949870.30650.093790.277894.21770.12600.01280.397411.67540.54770.299980.479581.57000.19870.0533Set 3EMD-based model0.446212.37840.56370.317776.450577.55100.20710.04270.596715.97070.74180.550374.315168.94200.20250.0716CEEMD-based model0.400410.94720.51950.269981.570083.33330.20140.02970.688518.80340.86870.754676.369960.75090.20840.1217SSA-based model1.002826.96961.24041.538756.314044.21770.33720.25721.011526.85551.25341.571156.849344.36860.35500.2524VMD-based model0.32658.88090.41500.172287.713391.15650.21080.01410.934925.24131.18511.404465.068548.12290.25130.2248Proposed model0.29468.26570.37940.145689.930594.55780.18600.00720.36559.68880.47360.224385.616488.73720.23490.0188The most notable results are bolded.MAPE is the average of absolute percentage error. MAE is the mean absolute error of forecasting results. RMSE is the square root of the average of the error squares. MSE is mean squared error. DC is directional change. FICP is forecast interval coverage probability. FINAW is forecast interval normalized average width. AWD is accumulated width deviation of testing dataset.Figure 3
Experiment II: comparison with models employing diverse data pretreatment techniques.Figure 4
Experiment III: comparison with classic individual models.For Set 1, the WACD-based combination system was more successful in processing noise data, as observed from all indicator values. In addition, this strategy obtained higher DC values, meaning that after removing the interference, the forecasting models better fit the changes of wind speed. The denoising effect of EMD technology followed the developed system, with MAPE values of 5.4278%, and 7.6773%, respectively. Moreover, when the number of prediction steps varied, the prediction ability of CEEMD-based and VMD-based combined models undulated greatly. When the FINAW of all models was almost the same, the forecast interval coverage probability of the proposed model was approximately 10% more than the other models in one-step forecasting. The AWD of the proposed model was 0.0345 smaller than the EMD-based model.For Set 2, the results showed that the prediction errors obtained by the WACD-based combined system under 7% were smaller than the other data that were pretreatment-based. Taking one-step prediction as example, the developed system’s RMSE was 0.3065, improving 30.40% and 54.12%, respectively, relative to the values of the CEEMD-based (0.4697) and SSA-based prediction systems (0.7126).In Set 3, for both one-step and multistep prediction, any evaluation index of the prediction system put forward in this paper was at the optimal level. More specifically, under the one-step forecasting situation, the WACD-based system had the best MAE, MAPE, RMSE, MSE, and DC at 0.2946, 8.2657%, 0.3794, 0.1456, and 89.9305, respectively. Nevertheless, the prediction performance of the SSA-based systems was the worst with the biggest MAPE at 26.9696%.Remark 2.
Different from the recombination idea of the denoising layer number based on subjective judgment, each subsequence obtained though decomposition was weighted and reconstructed to eliminate the adverse influence of high frequency noise. The results of experiment II show that the combination system based on WACD was always better than the other preprocessing technologies, which indicates that the idea of adaptive combination reconstruction of weight has great potential in the future.


## 3.3.4. Experiment III: Comparison with Traditional Forecasting Models

This experiment highlighted the contrast that exists between the prediction capability of the developed system and several traditional prediction methods, including ARIMA, BPNN, ELMNN, GRNN, and ELMAN. The detailed numerical values are shown in Table7 and Figure 4.Table 7
Prediction performance comparison table of the proposed model and classic individual models.
DatasetModelOne stepTwo stepsMAEMAPERMSEMSEDCFICPFINAWAWDMAEMAPERMSEMSEDCFICPFINAWAWDSet 1ARIMA0.40119.55130.53440.285653.924963.94560.10420.19690.646215.73360.82020.672753.082241.97950.10420.4227BPNN0.40659.73950.54670.298953.924956.12240.10420.24330.539513.40910.71080.505255.137047.09900.09780.4435ELMAN0.537812.97080.68280.466251.535844.21770.10020.41160.612514.82760.77400.599051.712341.97950.10520.4647GRNN0.455911.44200.59190.350455.290146.93880.09410.35360.569814.40800.73850.545456.849338.22530.08870.6553ELMNN0.39769.62030.53140.282459.385759.52380.10150.21350.539113.38730.69520.483351.027444.36860.09670.4563Proposed model0.16483.83540.22840.049892.013993.40280.09840.01930.24355.63170.33290.110985.958983.27650.12840.0494Set 2ARIMA0.510416.00710.74490.554953.924971.76870.13480.15300.836327.62481.19281.422852.054844.02730.11690.4541BPNN0.636026.65980.84680.717154.607548.63950.12210.32090.745027.58241.07281.150953.082247.44030.12650.4329ELMAN0.697922.72320.99210.984353.583650.00000.11750.32130.799426.13531.13221.281953.424746.07510.11780.4054GRNN0.672027.83850.92800.861254.607540.81630.08960.69600.786430.35581.09221.192855.137034.12970.09120.9722ELMNN0.494816.08980.70330.494668.600761.22450.10690.21640.751628.84631.06301.129954.794542.66210.12070.4272Proposed model0.21596.949870.30650.093790.277894.21770.12600.01280.397411.67540.54770.299980.479581.57000.19870.0533Set 3ARIMA0.707620.49920.91310.833859.385759.52380.21550.21661.153733.80701.46142.135853.767137.88400.19780.4784BPNN0.714920.34670.91880.844256.314057.48300.20290.22350.987028.49731.24631.553253.767139.24910.18310.4745ELMAN0.961127.70131.22501.500556.655342.85710.18530.41711.301937.02471.64732.713757.191840.61430.24490.4282GRNN0.765823.33170.98230.964858.703148.97960.17730.31371.004431.02811.25581.577157.191833.78840.16290.6429ELMNN0.691220.53300.90170.813161.092260.20410.19830.22720.973129.53521.23521.525754.109640.27300.18630.4719Proposed model0.29468.26570.37940.145689.930594.55780.18600.00720.36559.68880.47360.224385.616488.73720.23490.0188The most notable results are bolded.MAPE is the average of absolute percentage error. MAE is the mean absolute error of forecasting results. RMSE is the square root of the average of the error squares. MSE is mean squared error. DC is directional change. FICP is forecast interval coverage probability. FINAW is forecast interval normalized average width. AWD is accumulated width deviation of testing dataset.For Set 1, although the assessment index values were changed with the number of prediction steps, the point prediction and interval prediction performance of the proposed combination system was always optimal among all models. This phenomenon reveals that the combined system put forward in this study not only had the smallest prediction error but also achieved the highest stability.For Set 2, the RMSE values of ARIMA, BPNN, ELMAN, GRNN, and ELM were 1.1928, 1.0728, 1.1324, 1.0922, and 1.0630, respectively, in two-step forecasting. Apart from the worst ELMAN model, there was no significant difference among the prediction results obtained from the other single models; however, they were significantly weaker than the proposed composite system, whose RMSE was 0.5477, even smaller than 50% of the above model.In Set 3, when comparing the interval forecasting capability, the developed system was still the best relative to the other comparison methods in terms of the values of the predictive performance indexes. The values of FICP and FINAW should be studied together. For instance, with the two-step prediction, the developed system had the smallest FINAW and the biggest FICP, which were almost two-times greater than the single models.Remark 3.
Compared with the single models, the advantages of the combined strategy were outstanding. This conclusion has nothing to do with the selected evaluation index and the number of prediction steps. Therefore, these traditional individual prediction models are no longer suitable for current wind speed prediction.


## 4. Conclusions

Wind energy, as one of the fastest growing natural resources, is gradually occupying the leading position in energy structures. Wind speed forecasting, as the basis of wind energy utilization, has become a research focus [21].In this paper, an advanced combination system, which integrated the fortes of data preprocessing techniques and multiobjective optimization theory, was successfully proposed. The system achieved great improvements regarding both the prediction accuracy and stationarity of the previous models. To expand, the proposed forecasting system overcame the defects exposed in previous studies in the following aspects. (1) In contrast to the recombination idea of the denoising layer number based on subjective judgment, each subsequence obtained though decomposition was weighted and reconstructed to eliminate the adverse influence of high frequency noise. (2) One widespread applied statistical method, two traditional neural networks, and one deep learning model were used as base models for fully capturing all important characteristics in wind speed series. (3) An advanced multiobjective optimization algorithm, NNIA, was successfully applied in the forecasting system to decide the weights in the data preprocessing link and combined forecasting link.Empirical studies on the basis of realistic wind speed data indicated that the developed WACD-NNIA-ABGB system was clearly superior to other benchmark models. (1) Experiment I made a comparison between the proposed system and hybrid models with the same data pretreatment technique, and the results showed that the developed system obtained the lowest prediction error, with the MAPE mean values at 4.7405%, 9.2661%, and 8.9389%, which were reduced by more than 50%, compared with other hybrid models. (2) When contrasted with the models with diverse data preprocessing techniques, similar results indicated that the forecasting ability assessment values were still better than the other comparative models, revealing that the proposed decomposition and weighted reconstruction strategy surpassed the previous strategies and was a more effective tool. (3) In experiment III, the proposed forecasting system was compared with several traditional algorithms. According to the experiment results, we observed that regardless of the sites, the prediction performance of our system always ranked first, indicating that the combination theory was effective and successful. (4) Further discussions indicated that the developed system realized higher precision and stability relative to the other models. In conclusion, it is evident that the prediction system proposed in this paper shows great application potential and will be a reliable tool in smart grids.



---

*Source: 7854286-2020-11-23.xml*