# A Novel Hybrid Method for Short-Term Power Load Forecasting

Authors: Huang Yuansheng; Huang Shenhai; Song Jiayin
Journal: Journal of Electrical and Computer Engineering
 (2016)

---

## Abstract

Influenced by many uncertain and random factors, nonstationary, nonlinearity, and time-variety appear in power load series, which is difficult to forecast accurately. Aiming at locating these issues of power load forecasting, an innovative hybrid method is proposed to forecast power load in this paper. Firstly, ensemble empirical mode decomposition (EEMD) is used to decompose the power load series into a series of independent intrinsic mode functions (IMFs) and a residual term. Secondly, genetic algorithm (GA) is then applied to determine the best weights of each IMF and the residual term named ensemble empirical mode decomposition based on weight (WEEMD). Thirdly, least square support vector machine (LSSVM) and nonparametric generalized autoregressive conditional heteroscedasticity (NPGARCH) are employed to forecast the subseries, respectively, based on the characteristics of power load series. Finally, the forecasted power load of each component is summed as the final forecasted result of power load. Compared with other methods, the forecasting results of this proposed model applied to the electricity market of Pennsylvania-New Jersey-Maryland (PJM) indicate that the proposed model outperforms other models.


---

## Body

## 1. Introduction

In the operation of power system, power load forecasting is not only an important part of power system planning but also one of the most influential factors for the improvement of social economy developing, which has a significant impact on generation, transmission, and distribution. Therefore, it is crucial to design an efficient method to forecasting power load accurately by improving both economic efficiency and power supply quality and enhancing reliability of power system operation [1, 2].However, the power load forecasting results are influenced by many uncertain random factors, including changes of nature environment and electricity price and factors influencing power load [3, 4], causing the complex volatility characteristics of power load, which lead the time-series of power load series to inherent nonstationary and discrete in practice, even in the short-term dynamics. Notably, due to the factors, influencing the short-term power load, previous studies [5–7] mainly focus on the establishment of scientific forecasting method to enhance the adaptability for these multifactors. Therefore, the time-varying components of power load are hardly separated to extract the information in them.Under this circumstance, a novel hybrid method for power load forecasting has been proposed in this paper. Based on artificial algorithm, the advantages of econometric models and signal processing are emphasized by this research. Specifically, this hybrid method for load power forecasting includes ensemble empirical mode decomposition based on variable weights (WEEMD), genetic algorithm (GA), least square support vector machine (LSSVM), and nonparametric generalized autoregressive conditional heteroscedasticity (NPGARCH). Notably, in order to test the validity and feasibility of the model, amount of historical data of power load in American Electric Power (AEP) has been adopted to apply this new hybrid method and compare it with previously well-known methods for the power load forecasting accuracy. Results indicate that this hybrid method outperforms the compared methods with the forecasting accuracy.The rest of the paper is organized as follows: Section2 summaries previous literature over methodology for power load forecasting; Section 3 briefly describes EEMD, GA, nonparametric GARCH, and LSSVM and proposes the improved WEEMD, as well as presenting the procedure of the newly hybrid method; Section 4 evaluates the forecasting accuracy and tests the robustness of the proposed method, by comparing with other forecasting methods. Section 5 draws conclusions.


## 2. Related Literature Review

Current power load forecasting methods can be broadly divided into two categories: statistical analysis method [8–11] and artificial algorithm [12–15]. The traditional statistical analysis methods are always under the framework of time-series analysis [16]. By using ARMA and ARMA with weather inputs, Bolzern and Fronza [8] propose two short-term predictors of the winter power load. Respectively, the power load forecast results are satisfied with both predictors and similar in these two cases. In order to model certain behavior of energy consumption, Wang [9] employs conventional fuzzy systems based on the integrated algorithm, and results show that the monthly electricity consumption of Iran is accurately forecasted with this proposed method. In [10], considering the nonlinearity and volatility of power load series, analyzing the threshold characteristics, Wang et al. adopt a novel double-threshold generalized autoregressive conditional heteroscedasticity (DT-GARCH) model for short-term power load forecasting. To forecast the power load efficiently and scientifically, Liu [11] synthesizes the time-series modeling with the regression modeling. According to the residual error from regression forecasting, the accuracy of power load forecasting can be improved by correcting the error of GARCH model.As the artificial algorithm applied in the power load forecasting is prevalent in recent decades, there are plenty of researches related to that. Taking the characteristics of randomness, tendency, and periodicity of short-term power load into account, Sun and Ye [12] propose a model based on LSSVM and fruit fly algorithm (FOA) for short-term load forecasting. In order to reduce the nonlinear power load sequence and improve the accuracy of forecasting result, Hu and Chang [13] decompose the time-series into several components by local wave method and use optimal parameters to establish the DEMD-LSSVM model for component forecasting. Lauret et al. [14] utilize Bayesian techniques to design an optimal neural network (NN) based model for power load forecasting, where Bayesian technique offers great advantages on traditional neural network learning methods. Based on an artificial neural network (ANN), Shao et al. [15] proposed an approach combined with a fuzzy system for short-term power load forecasting. The load error is obtained from the historical information and past forecasted load errors are caused by fuzzy systems, and the final forecasted load can be obtained by adding the load error to preliminary load forecasted by ANN. In [16], to improve the forecasting performance by searching a suitable parameters combination, the paper presents an SVR-based electric load forecasting model by applying a novel algorithm named chaotic ant swarm optimization (CAS). Combined with the chaotic behavior of single ant and self-organization behavior of ant colony in the foraging process, the CAS is proposed to overcome premature local optimum. In [17], a combination of the wavelet transform (WT) and gray model is proposed for short-term power load forecasting, which is improved by particle swarm optimization (PSO). With this proposed method taking mean temperature, mean relative humidity, mean wind speed, and previous days load data into consideration and eliminating the high frequency of historical data built by the WT, the accuracy can be largely improved.With the improvement of the forecasting accuracy, Ghelardoni et al. [18] decompose time power load series with empirical mode decomposition (EMD) into two sets of components, respectively, describing the trend of energy consumption values. LSSVM is built to forecast these two components. The hybrid method proposed in [19] enhances the capability for forecasting. This proposed method is modeled based on supporting vector regression (SVR), EMD, and regression (AR), which can simultaneously provide forecasting with high accuracy and interpretability. In order to solve the problem of mode mixing and high frequency random components, Liu et al. [19] proposed an optimized method based on ensemble empirical mode decomposition (EEMD) and subsection particle swarm optimization (SS-PSO). By extracting and reconstructing intrinsic mode function (IMF), the power load series movement is well forecasted. In [20], the proposed method based on complementary ensemble empirical mode decomposition (CEEMD) fuzzy entropy and echo state network (ESN) with leaky integrator neurons (LiESN) enhances the forecasting accuracy of power load.In summary, although the characteristics of power load series are complex, such as being rich multidimensional, nonlinear, time-varying, and nonstationary, the previous literature has accumulated a great deal of experience about power load forecasting. Particularly, the forecasting methods have been constantly proposed and the forecasting performance has been continuously improved, all of which provide a significant foundation for the present study. However, we can still find that the method of power load forecasting can be problematical based on the analysis and research of previous literature, while there are also some uncertain factors of power load forecasting. For example, both EEMD and EMD decompose the original power load series into several components, and other forecasting methods are directly built based on the decomposed components in most of the literature, without analysis of the weights of each component in original power load series. Besides, as the power load is affected by great uncertainties, the power load series movement is difficult to capture, and there will be nonlinear and time-varying feature in power load series caused by no matter the external factors or the internal factors. But the existing methods for power load forecasting are usually not effective to separate nonlinear and time-varying components of power load and unable to extract their inherent moving mechanisms, which consequently affects the forecasting accuracy.Based on the potential faultiness mentioned above, a hybrid method for power load forecasting is developed and applied to PJM by this study. Specially, the ensemble empirical mode decomposition (EEMD) model decomposes power load into series of intrinsic mode functions (IMFs) and residual. And genetic algorithm (GA) is used to determine the weights of each component, which is based on deviation between fitted values and actual values. For weighted components with the heteroscedasticity character, GARCH is more favorable for forecasting. However, the LSSVM method is used in the forecasting of the subseries with the characteristics of nonstationary and nonlinearity subseries. On the one hand, power load forecasting can be more convenient and accurate, after the decomposition of the power load components because of the separation of other influential factors. On the other hand, with other neural network models, which cannot fundamentally solve problems of the local minimum, the difficulty in determination of hidden layer, and the slow training rate, the LSSVM can not only get over these disadvantages, but also can improve the accuracy of forecasting. Therefore, based on the components weighted by GA method, LSSVM method can be more in line with the power load study in this research.


## 3. Power Load Forecasting Methodology

### 3.1. Genetic Algorithms (GA)

Genetic algorithms code the candidate solutions of an optimization algorithm as a string of characters which are usually binary digits [21]. In accordance with the terminology that is borrowed from the field of genetics, this bit string is usually named as chromosome. The solution represented by its chromosome is considered as an individual. The algorithm starts with the initial generation of the population. The fitness of the individuals within the population is assessed, and new individuals are generated for the further generation. A number of genetic operators, containing selector operator, crossover, and mutation, are available for this purpose. By a number of fixed generations which is the termination condition, the best individual will be obtained with the max fitness, which is the global optimal solution of the issues.


### 3.2. The WEEMD Method

EMD, an effective method for signal processing, is gradually replaced by EEMD which overcomes the mixing model problem [22–25]. The essence of EEMD is to decompose a time-series into a set of independent intrinsic mode functions (IMFs) and the residue obtained by adding a random Gaussian white noise sequence, which is different from EMD. While the time-series is decomposed by EEMD, the IMF is a function, satisfying the following two conditions:(1)
In the whole data set, the number of extreme and the number of zero crossings must either differ or differ at most by one.(2)
At any point, mean values of the envelope, defined by both the local maxima and minima, are zero.For an arbitrary time-series,X(t), procedures of EEMD method can be described as follows:(1)
Add the white noise to power load seriesX(t), with i=1, and set the number of ensemble (M):(1)Xit=Xt+nit,whereni(t) denotes the ith added white noises series and Xi(t) represents the noise-added power load of the ith trial.(2)
Decompose the noise-added seriesXi(t) into J IMFs cij(j=1,2,…,J) by using EMD, where cij is the jth IMF of the ith trial and J is the number of IMFs.(3)
Repeat Steps(1) and (2) until i=M.(4)
Calculate the ensemble meancj(t) of N trails for each IMF, then cj(t)=(1/M)∑i=1Mcij(t), where cj(t)(j=1,…,J) is the jth IMF component by using EEMD.In this paper, weights assignment, for each IMF, is proposed based on EEMD, and the rationality of this weights assignment proved by this study is as follows: for an arbitrary time-series,X(t), setting any value of each IMF satisfies the function f(t), and by the conditions mentioned above, the number of extreme and zero crossings of each IMF is q, which are obtained by the following equations:(2)fjt=0,dfjtdt=0.Assigning weight to each IMF, then(3)ωjfjt=0,ωjdfjtdt=0,where ωj, a constant, is the weights of jth IMF.Obviously, the weights have no effect on the values, which satisfies the first conditions. Specially, mean values of the envelope will not be changed by weight assignment, and the deduction can be expressed as follows:(4)mmax+mmin=0.Assigning weight to it,(5)ωjmmax+mmin=0,where mmax is the mean value of the envelope of the local maxima and mmin is the mean values of the envelope of the local minima.GA method is built to determine the weights of each IMF, and individualsW=(ω1,ω2,…,ωJ) are first randomly generated as initial population W0. By using the genetic operators, selector operator, crossover, and mutation, the best individual W=(ω1,ω2,…,ωJ) can be obtained. Specially, in the selector operator, the fitness function is defined as follows:(6)εkt=Wk·C-Xt=ωk1,ωk2,…,ωkJ,ωr·c1c2⋮cJ-Xt=∑j=1Jωkjcjt-Xt,k=1,2,…,m,where εk(t) is error and m is the size of initial population.As the GA operator is designed to maximize the fitness function, the above minimization problem can be solved by using the following transformation:(7)fk=1εkt,where fk is the fitness of kth individual.The election probability of each individual ispk:(8)pk=fk∑k=1mfk.As described above, the process of EEMD is indeed like sifting, which has an effect to eliminate riding waves. The IMFs are extracted from the power load series and contain “information” about the time-series. This paper uses weights as the contribution of IMF toX(t); the greater the weight is, the larger the amount of “information” IMF contains.


### 3.3. The Nonparametric GARCH (1, 1) Method

Nonparametric GARCH (1, 1) model for the error fluctuation [26], which requires less assumptions, is defined as(9)z^t=σtφt,σt2=gσt-12,z^t-1,where z^t is estimation of random processes; φt is specific transformation sequence of zero mean and unit variance. σt2 is the conditional variance. σt is the volatility.Nonparametric model can be used to estimate the conditional variance, and (9) can be rewritten as follows:(10)z^t2=gσt-12,z^t-1+Vt,Vt=gσt-12,z^t-1φt2-1,where Vt is the martingale difference sequence.According to the equations above, the functionz^t2 is regressed between lagged variables z^t-1 and σt-12g can be estimated by using nonparametric smoothing method and the autoregression function. The calculation of nonparametric GARCH [27–30] is as follows:(1)
Using parameters GARCH (1, 1) model fits the volatilityσ^t,0, with estimation using maximum likelihood estimate parameter being employed as parameter estimation, where a=1.(2)
σt,m2 is as weight, and z^t2 is estimated by using zt-1 and σ^t-1,m-12. Smooth nonparametric estimation method is applied to obtain the estimation g^a of autoregression function g.(3)
Standard deviation is obtained by using(11)σ^t,a2=gaσ^t-1,m-12,z^t-1,2≤t≤n.(4)
Incrementm and return to step (2) if a>A, where M is a prespecified maximum number of iterations.The nonparametric estimation is an improvement over the parametric GARCH estimation of volatility. By means of step(1)–(4) continuous iteration, there is little to pick and choose between volatility estimates σ^t,a2 for various values of a. However, the algorithm can often be improved by averaging over the final k estimates to obtain(12)σ^t,∗=∑a=A-k+1Aσ^t,ak.


### 3.4. The LSSVM Method

LSSVM is proposed as an improved algorithm based on support vector machine (SVM) [31–35], with the given training data set {zi,yi}(i=1,2,…,N), with the input zi∈Rn and the output yi∈Rn. The following regression model is constructed by using nonlinear mapping function ϕ(zi):(13)fzi=WTϕzi+b.With the given training data set{zi,yi}(i=1,2,…,N), the optimization problem of LSSVM is defined as follows:(14)min⁡JW,ξi=12WTW+C2∑i=1Mξi,s.t.yi=WTφzi+b+ξi,i=1,2,…,n.According to the Kuhn-Tucker conditions, the LSSVM regression model can be expressed as(15)yzi=∑i=1naikz,zi+b.k ( z , z i ) is the kernel function, which can map variables to the feature space and avoid high dimensional complex difficulties. This paper applies RBF as the kernel function, which is defined as(16)kx,xi=exp⁡-x-xi22σ2,where σ is the kernel function parameter. The LSSVM method can be used by establishing the parameters σ and c.


### 3.5. The Hybrid Method for Power Load Forecasting

Considering the complex volatility characteristics of power load series, much more scientific forecasting models are required to address the nonlinearity and time variations. Under this circumstance, the WEEMD, based on EEMD, is proposed to extract different components of power load series and assign weights to each component (IMF) according to its standard deviation, where the LSSVM is presented to forecast the subseries with the characteristics of nonstationary and nonlinearity, and the nonparametric GARCH (1, 1) is used to forecast the subseries with heteroscedasticity. With this hybrid model, the power load movement can be well forecasted. The procedures of the improved model can be described as Figure1 and the concrete steps are given as follows:(1)
The power load series is first decomposed by EEMD inton intrinsic mode functions (IMFs) and one residual series, and then(17)Xt=∑j=1J-1Cjt+Rt,whereX(t) is the original power load series and Cj(t) and R(t) are decomposed from the series.(2)
Each IMF series and the residual series are assigned to be weighted by GA, which can be represented as follows:(18)Xt=∑j=1J-1ωjCjt+ωrRt.(3)
To verify the existence of ARCH effect and heteroscedasticity, ARCH-LM is used to test the subseries which is related to the stochastic error, and then(19)Xt=∑j=1mωjNjt+∑i=m+1JωiSit,whereNj(t) is the subseries with heteroscedasticity and Si(t) represents the subseries without heteroscedasticity.(4)
The LSSVM is built to forecast the future values ofSi(t); meanwhile, the nonparametric GARCH model is presented to forecast the future values of Nj(t), and their forecasted results are S^i(t) and N^j(t), respectively, which can be represented as follows:(20)X^t=∑j=1mωjN^jt+∑i=m+1JωiS^it,whereX^(t) is the forecasted value of power load series.Figure 1
The procedures of power load forecasting using the novel hybrid model.To examine the proposed hybrid method performance, three criteria—mean absolute error (MAE), mean absolute percentage error (MAPE), and root mean square error (RMSE)—are represented as follows:(21)MAE=1N∑t=1NOt-O^t,MAPE=1N∑t=1NOt-O^tOt,RMSE=1N∑t=1NOt-O^t2,where Ot and O^t represent the real and forecasted values, respectively, and N is the number of the forecasting performance evaluations.


## 3.1. Genetic Algorithms (GA)

Genetic algorithms code the candidate solutions of an optimization algorithm as a string of characters which are usually binary digits [21]. In accordance with the terminology that is borrowed from the field of genetics, this bit string is usually named as chromosome. The solution represented by its chromosome is considered as an individual. The algorithm starts with the initial generation of the population. The fitness of the individuals within the population is assessed, and new individuals are generated for the further generation. A number of genetic operators, containing selector operator, crossover, and mutation, are available for this purpose. By a number of fixed generations which is the termination condition, the best individual will be obtained with the max fitness, which is the global optimal solution of the issues.


## 3.2. The WEEMD Method

EMD, an effective method for signal processing, is gradually replaced by EEMD which overcomes the mixing model problem [22–25]. The essence of EEMD is to decompose a time-series into a set of independent intrinsic mode functions (IMFs) and the residue obtained by adding a random Gaussian white noise sequence, which is different from EMD. While the time-series is decomposed by EEMD, the IMF is a function, satisfying the following two conditions:(1)
In the whole data set, the number of extreme and the number of zero crossings must either differ or differ at most by one.(2)
At any point, mean values of the envelope, defined by both the local maxima and minima, are zero.For an arbitrary time-series,X(t), procedures of EEMD method can be described as follows:(1)
Add the white noise to power load seriesX(t), with i=1, and set the number of ensemble (M):(1)Xit=Xt+nit,whereni(t) denotes the ith added white noises series and Xi(t) represents the noise-added power load of the ith trial.(2)
Decompose the noise-added seriesXi(t) into J IMFs cij(j=1,2,…,J) by using EMD, where cij is the jth IMF of the ith trial and J is the number of IMFs.(3)
Repeat Steps(1) and (2) until i=M.(4)
Calculate the ensemble meancj(t) of N trails for each IMF, then cj(t)=(1/M)∑i=1Mcij(t), where cj(t)(j=1,…,J) is the jth IMF component by using EEMD.In this paper, weights assignment, for each IMF, is proposed based on EEMD, and the rationality of this weights assignment proved by this study is as follows: for an arbitrary time-series,X(t), setting any value of each IMF satisfies the function f(t), and by the conditions mentioned above, the number of extreme and zero crossings of each IMF is q, which are obtained by the following equations:(2)fjt=0,dfjtdt=0.Assigning weight to each IMF, then(3)ωjfjt=0,ωjdfjtdt=0,where ωj, a constant, is the weights of jth IMF.Obviously, the weights have no effect on the values, which satisfies the first conditions. Specially, mean values of the envelope will not be changed by weight assignment, and the deduction can be expressed as follows:(4)mmax+mmin=0.Assigning weight to it,(5)ωjmmax+mmin=0,where mmax is the mean value of the envelope of the local maxima and mmin is the mean values of the envelope of the local minima.GA method is built to determine the weights of each IMF, and individualsW=(ω1,ω2,…,ωJ) are first randomly generated as initial population W0. By using the genetic operators, selector operator, crossover, and mutation, the best individual W=(ω1,ω2,…,ωJ) can be obtained. Specially, in the selector operator, the fitness function is defined as follows:(6)εkt=Wk·C-Xt=ωk1,ωk2,…,ωkJ,ωr·c1c2⋮cJ-Xt=∑j=1Jωkjcjt-Xt,k=1,2,…,m,where εk(t) is error and m is the size of initial population.As the GA operator is designed to maximize the fitness function, the above minimization problem can be solved by using the following transformation:(7)fk=1εkt,where fk is the fitness of kth individual.The election probability of each individual ispk:(8)pk=fk∑k=1mfk.As described above, the process of EEMD is indeed like sifting, which has an effect to eliminate riding waves. The IMFs are extracted from the power load series and contain “information” about the time-series. This paper uses weights as the contribution of IMF toX(t); the greater the weight is, the larger the amount of “information” IMF contains.


## 3.3. The Nonparametric GARCH (1, 1) Method

Nonparametric GARCH (1, 1) model for the error fluctuation [26], which requires less assumptions, is defined as(9)z^t=σtφt,σt2=gσt-12,z^t-1,where z^t is estimation of random processes; φt is specific transformation sequence of zero mean and unit variance. σt2 is the conditional variance. σt is the volatility.Nonparametric model can be used to estimate the conditional variance, and (9) can be rewritten as follows:(10)z^t2=gσt-12,z^t-1+Vt,Vt=gσt-12,z^t-1φt2-1,where Vt is the martingale difference sequence.According to the equations above, the functionz^t2 is regressed between lagged variables z^t-1 and σt-12g can be estimated by using nonparametric smoothing method and the autoregression function. The calculation of nonparametric GARCH [27–30] is as follows:(1)
Using parameters GARCH (1, 1) model fits the volatilityσ^t,0, with estimation using maximum likelihood estimate parameter being employed as parameter estimation, where a=1.(2)
σt,m2 is as weight, and z^t2 is estimated by using zt-1 and σ^t-1,m-12. Smooth nonparametric estimation method is applied to obtain the estimation g^a of autoregression function g.(3)
Standard deviation is obtained by using(11)σ^t,a2=gaσ^t-1,m-12,z^t-1,2≤t≤n.(4)
Incrementm and return to step (2) if a>A, where M is a prespecified maximum number of iterations.The nonparametric estimation is an improvement over the parametric GARCH estimation of volatility. By means of step(1)–(4) continuous iteration, there is little to pick and choose between volatility estimates σ^t,a2 for various values of a. However, the algorithm can often be improved by averaging over the final k estimates to obtain(12)σ^t,∗=∑a=A-k+1Aσ^t,ak.


## 3.4. The LSSVM Method

LSSVM is proposed as an improved algorithm based on support vector machine (SVM) [31–35], with the given training data set {zi,yi}(i=1,2,…,N), with the input zi∈Rn and the output yi∈Rn. The following regression model is constructed by using nonlinear mapping function ϕ(zi):(13)fzi=WTϕzi+b.With the given training data set{zi,yi}(i=1,2,…,N), the optimization problem of LSSVM is defined as follows:(14)min⁡JW,ξi=12WTW+C2∑i=1Mξi,s.t.yi=WTφzi+b+ξi,i=1,2,…,n.According to the Kuhn-Tucker conditions, the LSSVM regression model can be expressed as(15)yzi=∑i=1naikz,zi+b.k ( z , z i ) is the kernel function, which can map variables to the feature space and avoid high dimensional complex difficulties. This paper applies RBF as the kernel function, which is defined as(16)kx,xi=exp⁡-x-xi22σ2,where σ is the kernel function parameter. The LSSVM method can be used by establishing the parameters σ and c.


## 3.5. The Hybrid Method for Power Load Forecasting

Considering the complex volatility characteristics of power load series, much more scientific forecasting models are required to address the nonlinearity and time variations. Under this circumstance, the WEEMD, based on EEMD, is proposed to extract different components of power load series and assign weights to each component (IMF) according to its standard deviation, where the LSSVM is presented to forecast the subseries with the characteristics of nonstationary and nonlinearity, and the nonparametric GARCH (1, 1) is used to forecast the subseries with heteroscedasticity. With this hybrid model, the power load movement can be well forecasted. The procedures of the improved model can be described as Figure1 and the concrete steps are given as follows:(1)
The power load series is first decomposed by EEMD inton intrinsic mode functions (IMFs) and one residual series, and then(17)Xt=∑j=1J-1Cjt+Rt,whereX(t) is the original power load series and Cj(t) and R(t) are decomposed from the series.(2)
Each IMF series and the residual series are assigned to be weighted by GA, which can be represented as follows:(18)Xt=∑j=1J-1ωjCjt+ωrRt.(3)
To verify the existence of ARCH effect and heteroscedasticity, ARCH-LM is used to test the subseries which is related to the stochastic error, and then(19)Xt=∑j=1mωjNjt+∑i=m+1JωiSit,whereNj(t) is the subseries with heteroscedasticity and Si(t) represents the subseries without heteroscedasticity.(4)
The LSSVM is built to forecast the future values ofSi(t); meanwhile, the nonparametric GARCH model is presented to forecast the future values of Nj(t), and their forecasted results are S^i(t) and N^j(t), respectively, which can be represented as follows:(20)X^t=∑j=1mωjN^jt+∑i=m+1JωiS^it,whereX^(t) is the forecasted value of power load series.Figure 1
The procedures of power load forecasting using the novel hybrid model.To examine the proposed hybrid method performance, three criteria—mean absolute error (MAE), mean absolute percentage error (MAPE), and root mean square error (RMSE)—are represented as follows:(21)MAE=1N∑t=1NOt-O^t,MAPE=1N∑t=1NOt-O^tOt,RMSE=1N∑t=1NOt-O^t2,where Ot and O^t represent the real and forecasted values, respectively, and N is the number of the forecasting performance evaluations.


## 4. Case Studies

Pennsylvania-New Jersey-Maryland (PJM), reliable operations and efficient wholesale market, is a fair and efficient electricity market, which provides information on electrical field. For our simulation, the short-term load power data from American Electric Power (AEP) are obtained hourly from 1/5/2015 to 10/8/2015. Data points from 1/5/2015 to 31/7/2015 are selected as training samples, and the data from 1/8/2015 to 10/8/2015 are selected as the test sample. Figure2 presents the power load of the training sample, which shows that the power load severely fluctuates periodically.Figure 2
Power load changes from 1/5/2015 to 31/7/2015.


### 4.1. Data Processing for IMF

According to procedure proposed above, first, the original power load series is decomposed by the EEMD into seven independent intrinsic mode functions and one residual, which are shown in Figure3.Figure 3
The decomposition of the power load series.In order to assign weights for each IMF, GA method is used to select the best individual. In Figure4, the results show that the average fitness is 0.545, while the algorithm iterates to 745 times, which do not change in later iterations, and this individual is the best one. Besides, the decomposition of the power load series with assigned weights is shown in Figure 5. Obviously, the assigned weights will not change the movement of the power load series but the domain which will be significant to the predicted values.Figure 4
Fitness curve.Figure 5
The decomposition of the power load series with assigned weights.Considering the impact caused by the closer prices on the further data, three training samples near the forecasting point are selected as input variable:(22)X=xd,t-3xd,t-2xd,t-1,where x(d,t) is the price on dth at time tth. Actually, the power load series is generally nonstationary, which will affect the forecasting results of LSSVM without handling. So, ARCH-LM test is proposed to resolve this problem. 1-hour return of the power load is calculated as R=ln⁡(xt/xt-1). Figure 6 provides the distribution of 1-hour return series from 1/5/2015 to 31/7/2015. Obviously, volatility cluster appeared in the residual, and the variance in the area of a is larger than b.Figure 6
The residual of the power load series with assigned weights.ARCH-LM test based on AR(3) is used to quantitatively test the heteroscedasticity, and WIMF1 is taken as an example to test the heteroscedasticity.In Table1, obviously, the significance of coefficients is close to zero, which indicates that significant autocorrelation appears in WIMF1 with Lag1, Lag2, and Lag3 and displays autocorrelation characteristic of WIMF1. Therefore, ARCH-LM test is applied to exam the conditional heteroscedasticity of returns series, and the result of test is shown in Table 2.Table 1
Autocorrelation test.
t-statistic Prob. Lag1 −22.61395 0.0000 Lag2 −30.69547 0.0000 Lag3 −18.65080 0.0000Table 2
Test of conditional heteroscedasticity of returns series.
F-statistic 4.571244 Prob. 0.0034 t-statistic Prob.  Lag1 3.281471 0.0010 Lag2 0.723646 0.4694 Lag3 −1.499284 0.1339TheF-statistic and the Lag1 of residual in Table 2 is under the significance level of 0.005. So the assumption of ARCH is not accepted. It is the WIMF1 which illustrates the heteroscedasticity.Each WIMF is tested by the same method; generally, the subseries of WIMF1, WIMF2, and WIMF5 have heteroscedasticity characteristic; and without the notable heteroscedasticity, there are WIMF3, WIMF4, WIMF6, WIMF7, and WIMF8. Then nonparametric GARCH model is established to forecast the power load. A large number of domestic and international demonstration analyses have shown that. GARCH(1,1) can accurately describe the fluctuation characteristics of the model, and therefore this paper adopts nonparametric GARCH (1,1) to imitate the power load series.


### 4.2. Result Analysis of the Hybrid Forecasting Method

The performance of LSSVM relays on the parameters ofγ and σ2, which, respectively, represent the regulation parameter and kernel parameter. A large number of domestic and international demonstration analyses indicate that two parameters are experientially determined. Based on a number of tests, this paper adopts 100 for γ and 0.1 for σ2. The main parameters in GA are listed in Table 3.Table 3
Parameter of GA.
Main parameter Value Initial population size 10 Evolution iteration 1000 Crossover probability 0.4 Mutation probability 0.05In Figure7, the error caused by the hybrid forecasting method is clearly shown. The error movement stably changes, the maximum of relative error should be no more than 10%, and meanwhile, the MAE, MAPE, and RMSE are 244.469, 1.56, and 407.38, respectively, which mean that the forecasting results are acceptable.Figure 7
The error of the hybrid method.


### 4.3. Comparative Analysis

#### 4.3.1. Hourly Power Load Forecasting Analysis

To demonstrate the forecasting performance of the novel hybrid method, LSSVM, BPNN, EEMD plus LSSVM, EEMD plus BPNN, and GA-LSSVM are employed as the comparative methods, which are shown in Figure8. Table 4 summarizes the values of the three error criteria, including MAE, RMSE, and MAPE, and the forecasted results of the six methods show that, using the proposed hybrid method, the power load series forecasted errors can be accepted. Notably, the MAE is less than 2%, and meanwhile, it is evident that MAE and RMSE are lower than the other methods, which implies that the forecasting accuracy of the proposed method appears better than the comparative methods. Compared with the MAPE of the hybrid method 1.56%, the suboptimum with the MAE 4.05% is worse than the proposed method, and the result indicates that the hybrid method, by using nonparametric GARCH (1, 1) to forecast the subseries with heteroscedasticity, has well captured the time-varying volatility features of the power load series. Meanwhile, by the results of EEMD-LSSVM and WEEMD-LSSVM, it is obvious that assigning weights to each IMF improves the forecasting accuracy. Besides, it verifies that EEMD method decomposes power load series to constitutive subseries forecasted more accurately than original series by directly comparing LSSVM with EEMD. Generally, the forecasted results of the proposed method are reasonable and much more accurate than the other method based on the hourly observations.Table 4
The error of power load forecasting using different methods based on hourly observation.
Hybrid method LSSVM BPNN WEEMD-LSSVM GA-LSSVM EEMD-LSSVM MAE (MW) 244.469 1313.7 1812.2 622.200 1083.4 909.020 MAPE (%) 1.56 8.44 11.81 4.05 7.12 6.08 RMSE (MW) 407.308 1552.5 2089.4 720.567 1220.7 1031.3Figure 8
The comparison between the prediction results based on hourly observation.


#### 4.3.2. Daily Power Load Forecasting Analysis

As the data frequency is a significant factor for the sensitivity of the time-series forecasting, to examine the robustness of the hybrid method this study adopts a daily observation method to forecast the power load. And the forecasted power load is decomposed to ten parts shown in Figure9.Figure 9
The comparison between the prediction results based on daily observation.Table5 shows the errors of forecasted results among different methods, and the MAE, MAPE, and RMSE of the hybrid method can be accepted with smaller MAE and RSME values and the MAPE is less than 1% comparing with the other methods, which indicates that the hybrid method has a better performance than other five methods. Besides, due to the differences in their characteristics, the forecasting accuracy can be improved and clustered by using WEEMD method. Hence, this newly proposed hybrid method for power load forecasting in this paper has relatively reliable robustness with respect to the data frequency.Table 5
The error of power load forecasting using different methods based on daily observation.
Hybrid method LSSVM BPNN WEEMD-LSSVM GA-LSSVM EEMD-LSSVM MAE (MW) 1441.263 13343.236 16027.205 5023.101 10623.505 10729.767 MAPE (%) 0.389 3.488 4.146 1.313 2.749 2.854 RMSE (MW) 1845.988 17460.604 22325.684 6561.656 14855.986 13914.859


## 4.1. Data Processing for IMF

According to procedure proposed above, first, the original power load series is decomposed by the EEMD into seven independent intrinsic mode functions and one residual, which are shown in Figure3.Figure 3
The decomposition of the power load series.In order to assign weights for each IMF, GA method is used to select the best individual. In Figure4, the results show that the average fitness is 0.545, while the algorithm iterates to 745 times, which do not change in later iterations, and this individual is the best one. Besides, the decomposition of the power load series with assigned weights is shown in Figure 5. Obviously, the assigned weights will not change the movement of the power load series but the domain which will be significant to the predicted values.Figure 4
Fitness curve.Figure 5
The decomposition of the power load series with assigned weights.Considering the impact caused by the closer prices on the further data, three training samples near the forecasting point are selected as input variable:(22)X=xd,t-3xd,t-2xd,t-1,where x(d,t) is the price on dth at time tth. Actually, the power load series is generally nonstationary, which will affect the forecasting results of LSSVM without handling. So, ARCH-LM test is proposed to resolve this problem. 1-hour return of the power load is calculated as R=ln⁡(xt/xt-1). Figure 6 provides the distribution of 1-hour return series from 1/5/2015 to 31/7/2015. Obviously, volatility cluster appeared in the residual, and the variance in the area of a is larger than b.Figure 6
The residual of the power load series with assigned weights.ARCH-LM test based on AR(3) is used to quantitatively test the heteroscedasticity, and WIMF1 is taken as an example to test the heteroscedasticity.In Table1, obviously, the significance of coefficients is close to zero, which indicates that significant autocorrelation appears in WIMF1 with Lag1, Lag2, and Lag3 and displays autocorrelation characteristic of WIMF1. Therefore, ARCH-LM test is applied to exam the conditional heteroscedasticity of returns series, and the result of test is shown in Table 2.Table 1
Autocorrelation test.
t-statistic Prob. Lag1 −22.61395 0.0000 Lag2 −30.69547 0.0000 Lag3 −18.65080 0.0000Table 2
Test of conditional heteroscedasticity of returns series.
F-statistic 4.571244 Prob. 0.0034 t-statistic Prob.  Lag1 3.281471 0.0010 Lag2 0.723646 0.4694 Lag3 −1.499284 0.1339TheF-statistic and the Lag1 of residual in Table 2 is under the significance level of 0.005. So the assumption of ARCH is not accepted. It is the WIMF1 which illustrates the heteroscedasticity.Each WIMF is tested by the same method; generally, the subseries of WIMF1, WIMF2, and WIMF5 have heteroscedasticity characteristic; and without the notable heteroscedasticity, there are WIMF3, WIMF4, WIMF6, WIMF7, and WIMF8. Then nonparametric GARCH model is established to forecast the power load. A large number of domestic and international demonstration analyses have shown that. GARCH(1,1) can accurately describe the fluctuation characteristics of the model, and therefore this paper adopts nonparametric GARCH (1,1) to imitate the power load series.


## 4.2. Result Analysis of the Hybrid Forecasting Method

The performance of LSSVM relays on the parameters ofγ and σ2, which, respectively, represent the regulation parameter and kernel parameter. A large number of domestic and international demonstration analyses indicate that two parameters are experientially determined. Based on a number of tests, this paper adopts 100 for γ and 0.1 for σ2. The main parameters in GA are listed in Table 3.Table 3
Parameter of GA.
Main parameter Value Initial population size 10 Evolution iteration 1000 Crossover probability 0.4 Mutation probability 0.05In Figure7, the error caused by the hybrid forecasting method is clearly shown. The error movement stably changes, the maximum of relative error should be no more than 10%, and meanwhile, the MAE, MAPE, and RMSE are 244.469, 1.56, and 407.38, respectively, which mean that the forecasting results are acceptable.Figure 7
The error of the hybrid method.


## 4.3. Comparative Analysis

### 4.3.1. Hourly Power Load Forecasting Analysis

To demonstrate the forecasting performance of the novel hybrid method, LSSVM, BPNN, EEMD plus LSSVM, EEMD plus BPNN, and GA-LSSVM are employed as the comparative methods, which are shown in Figure8. Table 4 summarizes the values of the three error criteria, including MAE, RMSE, and MAPE, and the forecasted results of the six methods show that, using the proposed hybrid method, the power load series forecasted errors can be accepted. Notably, the MAE is less than 2%, and meanwhile, it is evident that MAE and RMSE are lower than the other methods, which implies that the forecasting accuracy of the proposed method appears better than the comparative methods. Compared with the MAPE of the hybrid method 1.56%, the suboptimum with the MAE 4.05% is worse than the proposed method, and the result indicates that the hybrid method, by using nonparametric GARCH (1, 1) to forecast the subseries with heteroscedasticity, has well captured the time-varying volatility features of the power load series. Meanwhile, by the results of EEMD-LSSVM and WEEMD-LSSVM, it is obvious that assigning weights to each IMF improves the forecasting accuracy. Besides, it verifies that EEMD method decomposes power load series to constitutive subseries forecasted more accurately than original series by directly comparing LSSVM with EEMD. Generally, the forecasted results of the proposed method are reasonable and much more accurate than the other method based on the hourly observations.Table 4
The error of power load forecasting using different methods based on hourly observation.
Hybrid method LSSVM BPNN WEEMD-LSSVM GA-LSSVM EEMD-LSSVM MAE (MW) 244.469 1313.7 1812.2 622.200 1083.4 909.020 MAPE (%) 1.56 8.44 11.81 4.05 7.12 6.08 RMSE (MW) 407.308 1552.5 2089.4 720.567 1220.7 1031.3Figure 8
The comparison between the prediction results based on hourly observation.


### 4.3.2. Daily Power Load Forecasting Analysis

As the data frequency is a significant factor for the sensitivity of the time-series forecasting, to examine the robustness of the hybrid method this study adopts a daily observation method to forecast the power load. And the forecasted power load is decomposed to ten parts shown in Figure9.Figure 9
The comparison between the prediction results based on daily observation.Table5 shows the errors of forecasted results among different methods, and the MAE, MAPE, and RMSE of the hybrid method can be accepted with smaller MAE and RSME values and the MAPE is less than 1% comparing with the other methods, which indicates that the hybrid method has a better performance than other five methods. Besides, due to the differences in their characteristics, the forecasting accuracy can be improved and clustered by using WEEMD method. Hence, this newly proposed hybrid method for power load forecasting in this paper has relatively reliable robustness with respect to the data frequency.Table 5
The error of power load forecasting using different methods based on daily observation.
Hybrid method LSSVM BPNN WEEMD-LSSVM GA-LSSVM EEMD-LSSVM MAE (MW) 1441.263 13343.236 16027.205 5023.101 10623.505 10729.767 MAPE (%) 0.389 3.488 4.146 1.313 2.749 2.854 RMSE (MW) 1845.988 17460.604 22325.684 6561.656 14855.986 13914.859


## 4.3.1. Hourly Power Load Forecasting Analysis

To demonstrate the forecasting performance of the novel hybrid method, LSSVM, BPNN, EEMD plus LSSVM, EEMD plus BPNN, and GA-LSSVM are employed as the comparative methods, which are shown in Figure8. Table 4 summarizes the values of the three error criteria, including MAE, RMSE, and MAPE, and the forecasted results of the six methods show that, using the proposed hybrid method, the power load series forecasted errors can be accepted. Notably, the MAE is less than 2%, and meanwhile, it is evident that MAE and RMSE are lower than the other methods, which implies that the forecasting accuracy of the proposed method appears better than the comparative methods. Compared with the MAPE of the hybrid method 1.56%, the suboptimum with the MAE 4.05% is worse than the proposed method, and the result indicates that the hybrid method, by using nonparametric GARCH (1, 1) to forecast the subseries with heteroscedasticity, has well captured the time-varying volatility features of the power load series. Meanwhile, by the results of EEMD-LSSVM and WEEMD-LSSVM, it is obvious that assigning weights to each IMF improves the forecasting accuracy. Besides, it verifies that EEMD method decomposes power load series to constitutive subseries forecasted more accurately than original series by directly comparing LSSVM with EEMD. Generally, the forecasted results of the proposed method are reasonable and much more accurate than the other method based on the hourly observations.Table 4
The error of power load forecasting using different methods based on hourly observation.
Hybrid method LSSVM BPNN WEEMD-LSSVM GA-LSSVM EEMD-LSSVM MAE (MW) 244.469 1313.7 1812.2 622.200 1083.4 909.020 MAPE (%) 1.56 8.44 11.81 4.05 7.12 6.08 RMSE (MW) 407.308 1552.5 2089.4 720.567 1220.7 1031.3Figure 8
The comparison between the prediction results based on hourly observation.


## 4.3.2. Daily Power Load Forecasting Analysis

As the data frequency is a significant factor for the sensitivity of the time-series forecasting, to examine the robustness of the hybrid method this study adopts a daily observation method to forecast the power load. And the forecasted power load is decomposed to ten parts shown in Figure9.Figure 9
The comparison between the prediction results based on daily observation.Table5 shows the errors of forecasted results among different methods, and the MAE, MAPE, and RMSE of the hybrid method can be accepted with smaller MAE and RSME values and the MAPE is less than 1% comparing with the other methods, which indicates that the hybrid method has a better performance than other five methods. Besides, due to the differences in their characteristics, the forecasting accuracy can be improved and clustered by using WEEMD method. Hence, this newly proposed hybrid method for power load forecasting in this paper has relatively reliable robustness with respect to the data frequency.Table 5
The error of power load forecasting using different methods based on daily observation.
Hybrid method LSSVM BPNN WEEMD-LSSVM GA-LSSVM EEMD-LSSVM MAE (MW) 1441.263 13343.236 16027.205 5023.101 10623.505 10729.767 MAPE (%) 0.389 3.488 4.146 1.313 2.749 2.854 RMSE (MW) 1845.988 17460.604 22325.684 6561.656 14855.986 13914.859


## 5. Conclusions

To address the problem of power load forecasting with the characteristic of nonstationary, nonlinearity, and time-varying, this paper proposes a novel hybrid method for power load forecasting. The data frequency has been changed to test the robustness of the proposed method. Besides, other five methods presented by this study are compared with the proposed one to verify the accuracy of hybrid method by different criteria presented above. In the end, several conclusions are drawn as follows.(a) The newly proposed decomposition algorithm named WEEMD has a better performance than EEMD method. (b) Due to the differences in their characteristics improvement of the forecasting accuracy, the components are clustered. (c) Regardless of the influence of data frequency or the fluctuation of time-series, the proposed hybrid method has excellent forecasting performance for power load.



---

*Source: 2165324-2016-08-11.xml*