# A Hybrid Forecasting Model Based on Empirical Mode Decomposition and the Cuckoo Search Algorithm: A Case Study for Power Load

Authors: Jiani Heng; Chen Wang; Xuejing Zhao; Jianzhou Wang
Journal: Mathematical Problems in Engineering
 (2016)

---

## Abstract

Power load forecasting always plays a considerable role in the management of a power system, as accurate forecasting provides a guarantee for the daily operation of the power grid. It has been widely demonstrated in forecasting that hybrid forecasts can improve forecast performance compared with individual forecasts. In this paper, a hybrid forecasting approach, comprising Empirical Mode Decomposition, CSA (Cuckoo Search Algorithm), and WNN (Wavelet Neural Network), is proposed. This approach constructs a more valid forecasting structure and more stable results than traditional ANN (Artificial Neural Network) models such as BPNN (Back Propagation Neural Network), GABPNN (Back Propagation Neural Network Optimized by Genetic Algorithm), and WNN. To evaluate the forecasting performance of the proposed model, a half-hourly power load in New South Wales of Australia is used as a case study in this paper. The experimental results demonstrate that the proposed hybrid model is not only simple but also able to satisfactorily approximate the actual power load and can be an effective tool in planning and dispatch for smart grids.


---

## Body

## 1. Introduction

In a power system, the short-term power load forecasting is very important for the stable operation of the system. Accurate forecasting is a guarantee in the development of preventive maintenance plans, which include generator safeguards, power system reliability estimation, and scheduling dispatch [1, 2]. High-accuracy power load forecasts improve the economic and social benefits of power grid management, which reduce generation costs, improve the security of power systems, and help administrators develop optimal plans. Moreover, accurate load forecasting is crucial in forecasts of the power price in power markets [3]. Therefore, developing power load forecasting techniques to achieve accurate, simple, and fast load forecasts is necessary. Thus far, many short-term power load forecasting methods have been proposed, and these methods can be mainly divided into three categories: conventional methods, modern forecasting methods, and hybrid forecasting methods. Conventional methods include multiple linear regression analysis [4, 5], time series [6, 7], state space models [8], general exponential smoothing [9], and knowledge-based methods. However, these methods cannot provide appropriate nonlinear mathematical relationships to express actual power loads. The primary modern forecasting methods are intelligent evolutionary algorithms [10, 11], expert systems [12, 13], neural networks [14–17], and fuzzy inference [18]. Intelligent algorithms and neural networks obtain good performance because of their clear patterns, easy implementation, and strong ability to address the problem. Hybrid forecasting methods, proposed to avoid the shortcomings that exist in individual forecasting methods, have become increasingly prevalent [19, 20]. A detailed introduction of the three categories is given below.The deduction processes of traditional forecasting methods are rigorous, and most of them are based on traditional mathematics theories such as statistics, calculus, and modeling by subjective data analysis [21]. The main idea of trend extrapolation technology is to look for the trend of data changes, according to the trend equation, to forecast future data. The method is simple, and, especially for smooth power load changes, it can achieve a good prediction effect. Its deficiency is that its precision is greatly influenced by the random load component [22]. The regression analysis method is often applied to short-term load forecasting [23]. This method has many advantages such as a simple principle and better quality of data which leads to better precision; however, the selection of the main factors affecting power load in the model is difficult as many factors that affect the forecasting accuracy are hard to quantify. This model is lacking any self-study capability, and the input variable and output variable cannot be revised automatically [24]. With years of development, the time series forecasting method has become a mature theory method and has been applied to power load forecasting [25]. The basic time series prediction models mainly include AR, MA, and ARMA [26]. Although the time series forecasting method has advantages such as only requiring a small volume of historical data and a small amount of calculation and the fast speed of its calculation, this method has certain limitations such as its inability to reflect the influence of meteorological factors and how its forecasting accuracy will decrease with the increase of the prediction step [27]. ANN [17] is a type of nonlinear simulation of the human brain information processing system with an intelligent processing process; for an inaccurate variation trend, this method also has a good ability to adapt, is able to grasp information and keep on learning, and has good knowledge reasoning and self-optimization [28]. An expert system is a computer system based on the knowledge of the programming approach, mainly a software system, and the main components of an expert system include the inference engine of the system, the expert knowledge base, the explain interface, and the knowledge acquisition module. An expert system is a program that has decision-making capabilities based on reasoned knowledge; however, this method is limited by whether the expert knowledge is complete [29]. The grey forecasting method is an important technique in grey theory, and it uses approximate differential equations to describe future tendencies for a time series [30]. The limitation of this method is that the greater the dispersion degree of data, the worse the forecasting accuracy. Although traditional forecasting methods and forecasting methods based on intelligent computing have their respective applications, it is difficult to achieve better results when using one of them by itself [31]. In the literature related to forecasting [31–34], the forecasting results are not quite as good with any single forecasting model. The primary reason is that single forecasting models cannot extract the complicated factors encountered in reality.Due to the limitations of the forecasting capacity of a single model, it cannot always be optimal in all cases. In this paper, a novel hybrid model was developed with the hope of obtaining more accurate power load forecasting results. The proposed hybrid wind speed forecasting model can be grouped into four steps. Firstly, we used the empirical model decomposition technique, which represents a nonstationary data analysis technique, to reconstruct the original wind speed series. Secondly, a WNN model was employed to create the power load forecasting, and the parameters in the WNN model were tuned by the CSA. The simulation results illustrate that the hybrid model is an effective method in power load forecasting. The main contributions of this paper are summarized as follows:(1)
The CSA algorithm is applied to choose the optimal initial weight in the WNN model, which always leads to unstable forecasting error.
(2)
In the field of power load forecasting, the proposed hybrid model is manifested as a valid method with efficient computation and satisfactory forecasting accuracy.
(3)
Considering the skewness and kurtosis of the forecasting accuracy distribution, the developed forecasting availability is proposed as an effective evaluation criterion for model selection in the power load forecasting field.This paper is organized as follows. First, we outline the concept of models used in this paper, including empirical model decomposition, WNN, CSA, BPNN, GABPNN, and EMD-CSAWNN. Second, the modeling processes of the methods mentioned above are introduced. Simulation results are presented and analyzed. Finally, the overall conclusion is included.


## 2. Methodology

In this section, the required individual tools will be presented concisely, including the empirical model decomposition technique, BPNN, the WNN model, and the CSA and GA algorithms. Moreover, the proposed hybrid approach will be described in detail. In addition, the structure of the feed-forward neural network will be confirmed.


### 2.1. Empirical Mode Decomposition

Empirical model decomposition was proposed by Karthikeyan and Kumar as an adaptive method for nonstationary time series analysis, and it is now widely used. It can be applied to any type of signal decomposition [35]. Thus, it has obvious advantages in processing nonstationary and nonlinear series. The foundation of this technique is to decompose a time series into a finite set of several IMFs and a residue [36].Definition 1.
All IMFs are defined to satisfy the following conditions: (1) the number of local extreme points and the number of zero crossings must be equal or at least differ by only one; (2) the mean value of the upper envelope and lower envelope is zero.Definition 2.
The stoppage criterion determined is defined as(1)
S
D
k
=
∑
t
=
0
T
h
1
k
-
1
t
-
h
1
k
t
2
∑
t
=
0
T
h
1
k
-
1
2
t
.
The sifting process stops whenS
D
k is smaller than a pregiven value. The process of decomposition is over when the value of S
D
k is between 0.2 and 0.3. Additional details of the empirical model decomposition technique are illustrated in Figure 1.Figure 1
The construction of empirical model decomposition.


### 2.2. Artificial Neural Network (ANN)

#### 2.2.1. Confirmation of the Structure of the Network

The ANN has received considerable attention as a powerful computational tool for forecasting in many fields since 1980. ANN models always outperform statistical models because of their ability to map the inputs onto outputs via simple computation [37]. We discuss the feed-forward neural network in this paper because of its strong learning ability and simple structure. The determination of the network structure is as follows [38].Definition 3.
Given an arbitrary continuous functionf
:
[
0,1
]
n
→
R
m, f can be accurately approximated using a three-layer forward neural network realization. The first layer of the network is the input layer, containing n neurons. The middle layer is the hidden layer, containing 2
n
+
1 neurons. The third layer is the output layer, which has m neurons.Definition 4.
Letϕ be bounded continuous monotone in function, let K be the bounded compact subsets, and let K
⊂
R
n, f
X
=
f
x
1
,
x
2
,
…
,
x
n be a real continuous function in f
X
∈
K. ∀
ε
>
0, ∃ integer N, real constant C
i, θ
i
i
=
1,2
,
…
,
N and W
i
j
i
=
1,2
,
…
,
N, such that f
^
i
j
x
1
,
x
2
,
…
x
n
=
∑
i
=
1
N
C
i
ϕ
∑
j
=
1
n
W
i
j
x
j
-
θ
i satisfying max
⁡
f
^
x
1
,
…
,
x
n
-
f
x
1
,
…
,
x
n
<
ε. Note that, ∀
ε
>
0, ∃ a three-layer network structure, and the output function of the hidden layer is ϕ
x, the output function of the input and output layer is linear, and the total relationship of I
/
0 isf
x
1
,
…
,
x
n, such that max
⁡
f
^
x
1
,
…
,
x
n
-
f
x
1
,
…
,
x
n
<
ε.Proof.
Becausef
X
X
=
x
1
,
…
,
x
n is a continuous function, f
X
∈
K, K
⊂
R
n, assume that f
X is a bounded C
∞ function. Based on the Paley-Wiener Theorem [39], the Fourier transform F
W
W
=
W
1
,
…
,
W
n is a real analytic function of f
X, ∀ a constant C
N, such that F
W
≤
C
N
1
+
W
-
N. We define such a function: (2)
I
A
x
1
,
…
,
x
n
=
∫
-
A
A
⋯
∫
-
A
A
ϕ
∑
i
=
1
n
x
i
W
i
-
W
0
1
2
π
n
ϕ
1
F
W
1
,
…
,
W
n
exp
⁡
i
W
0
d
W
0
⋯
d
W
n
,
I
∞
,
A
x
1
,
…
,
x
n
=
∫
-
A
A
⋯
∫
-
A
A
∫
-
∞
∞
ϕ
∑
i
=
1
n
x
i
W
i
-
W
0
1
2
π
n
ϕ
1
F
W
1
,
…
,
W
n
exp
⁡
i
W
0
d
W
0
⋯
d
W
n
,
J
A
x
1
,
…
,
x
n
=
1
2
π
n
∫
-
∞
∞
⋯
∫
-
∞
∞
F
W
1
,
…
,
W
n
exp
⁡
i
∑
i
=
1
n
x
i
W
i
d
W
1
⋯
d
W
n
.
Using an estimation ofF
·, l
i
m
A
→
∞
J
A
x
1
,
…
,
x
n
=
f
x
1
,
…
,
x
n can be proved: that is, J
A is uniform convergence. Thus, l
i
m
A
→
∞
I
∞
,
A
x
1
,
…
,
x
n
=
f
x
1
,
…
,
x
n. That is, ∀
ε
>
0, ∃
A
>
0, such that m
a
x
x
∈
R
n
I
∞
,
A
x
1
,
…
,
x
n
-
f
x
1
,
…
,
x
n
<
E
/
2. Moreover, for A
′
>
0, let I
A
′
,
A
x
1
,
…
,
x
n
=
∫
-
A
A
⋯
∫
-
A
A
∫
-
A
′
A
ϕ
∑
i
=
1
n
x
i
W
i
-
W
0
·
(
1
/
2
π
n
ϕ
1
)
F
W
1
,
…
,
W
n
·
exp
⁡
i
W
0
d
W
0
⋯
d
W
n. For ε
>
0, A
′
>
0 can be found such that m
a
x
x
∈
K
I
A
,
A
x
1
,
…
,
x
n
-
I
∞
,
A
x
1
,
…
,
x
n
<
ε
/
2. That is to say, f
X is a uniform approximation by the integral I
A
′
,
A
X, I
A
′
,
A
X
∈
K. If A
i
>
0
i
=
1,2
,
…
,
m, K is a bounded subset, K
⊂
R
n, and h
i
x
1
,
…
,
x
m
;
t
1
,
…
,
t
m is a uniformly continuous function of -
A
1
,
A
1
×
⋯
×
-
A
m
,
A
m
×
K. There are functions: (3)
H
t
=
∫
-
A
1
A
1
⋯
∫
-
A
m
A
m
h
x
1
,
…
,
x
m
;
t
1
,
…
,
t
m
d
x
1
⋯
d
x
m
,
H
N
t
=
2
A
1
⋯
2
A
m
N
m
×
∑
k
1
⋯
k
m
=
0
N
-
1
h
-
A
1
+
k
1
·
2
A
1
N
⋯
A
m
+
k
m
·
2
A
m
N
t
1
⋯
t
m
.
∀
ε
>
0, ∃
N
0, N
0 is a natural number, when N
≥
N
0, such that(4)
max
l
∈
K
⁡
H
t
-
H
N
t
<
ε
.Thus, I
A
′
,
A
X is a uniform approximate on K, (5)
ϕ
∑
i
=
0
n
x
i
W
i
-
W
0
=
ϕ
∑
i
=
1
n
W
i
x
i
δ
-
W
0
+
α
-
ϕ
∑
i
=
1
n
W
i
x
i
δ
-
W
0
-
α
.The formula above can be achieved by a three-layer neural network. Thus, f
X can be approximated by that neural network.Definition 5.
∀
ε
>
0 and f
:
[
0,1
]
n
→
R
m, there exists a three-layer structure that can approximate f in any square error precision of ε.The definition above proves that,∀
f
:
0,1
n
→
R
m, we can use a feed-forward neural network with a three-layer structure n
×
2
n
+
1
×
m to approximate it accurately. Thus, this part not only proves the existence of the mapping network but also demonstrates the network structure of the mapping. In summary, this paper adopts the three-layer neural network as the basis neural network.


#### 2.2.2. BPNN

BPNN is a type of multilayer feed-forward neural network with an error back propagation learning process. The structure of BPNN is illustrated in Figure2. Details of BPNN are introduced in [40].Figure 2
The construction of GABPNN and CSAWNN.


#### 2.2.3. WNN

WNN, a feed-forward network, is generally multilayer [41]. It is widely applied in signal processing because of its advantages of the localization property and generalization ability [42]. The structure of WNN is shown in Figure 2.


### 2.3. Neural Network Optimized by an Intelligence Algorithm

The intelligent optimization algorithm provides an efficient and powerful mathematical tool for optimizing the initial weights and thresholds of the ANN [43].


#### 2.3.1. CSA

Cuckoo search is a heuristic swarm intelligence algorithm inspired by the behavior of the obligating brood parasitism of cuckoo species [44]. CSA is utilized in this paper for its stronger capability of global optimization [43].Definition 6.
To simulate the mode of cuckoo breeding, three idealized assumptions are presented, as follows: (1) each cuckoo selects nest randomly and dumps only one egg at a time, (2) the eggs with high quality will be carried over to the next generation, and (3) the available nest numbern is fixed, and the probability of the host bird discovering the exotic egg is p
α.Definition 7.
The Lévy flight model simulates the process of the nest-seeking characteristic of cuckoo, and the update formula of the path and location is as follows:(6)
X
i
i
t
e
r
+
1
=
X
i
i
t
e
r
+
α
⊕
L
λ
.
X
i
i
t
e
r represents the location of nest i at generation i
t
e
r, α is the step length-controlled factor, ⊕ is the point-to-point multiplication, and L
e
v
y
λ obeys the Lévy distribution with a random search path of parameter λ:(7)
L
λ
=
0.01
u
v
X
i
i
t
e
r
-
X
∗
i
t
e
r
. Here u
=
t
-
λ, 1
<
λ
<
3, u and v obey the normal distribution, u
~
N
0
,
σ
u
2, v
~
N
0,1, σ
u
=
Γ
λ
sin
⁡
0.5
π
λ
-
1
/
2
λ
-
2
/
2
Γ
0.5
λ
λ
-
1
1
/
λ
-
1, X
∗
i
t
e
r denotes the location of the best nest at generation i
t
e
r, and Γ is the standard gamma function with unbounded variance and mean of the probability distribution.


#### 2.3.2. GA

GA is a population-based optimization algorithm that simulates natural genetic mechanisms and biological evolutionism. It possesses a capacity for powerful global optimization [45]. The principle of GA relies on a random process, which is constituted by the processes of selection, crossover, and mutation [46]. The implementation process is shown in Figure 2.


### 2.4. GABPNN

The primary mechanism of GABPNN is composed of three parts: GA optimization, determination of the BPNN structure, and forecasting covered by BPNN [47]. The pseudo-code for GABPNN is as shown in Algorithm 1.Algorithm 1:GABPNN.
(
1) GENERATE  i
t
e
r
=
0 /∗ Calculate the fitness of each individual in the population. ∗/
(2) DO WHILE i
t
e
r
≤
i
t
e
r
m
a
x
(3)    i
t
e
r
=
i
t
e
r
+
1 /∗ Record the best fitness values and the average fitness values. ∗/
(4)    FOR EACH  1
≤
i
≤
P  DO /∗The process of the select operation. ∗/
(5)         p
i
i
t
e
r
=
f
i
i
t
e
r
⋅
/
∑
i
=
1
P
f
i
i
t
e
r
⋅
(6)         FOR EACH  1
≤
j
≤
s  DO /∗ The process of the crossover and mutation operations. ∗/
(7)           w
k
j
i
t
e
r
=
w
k
j
i
t
e
r
1
-
b
+
w
l
j
i
t
e
r
b  w
l
j
i
t
e
r
=
w
l
j
i
t
e
r
1
-
b
+
w
k
j
i
t
e
r
b
(8)           w
i
j
i
t
e
r
=
w
i
j
i
t
e
r
+
w
i
j
i
t
e
r
-
w
m
a
x
i
t
e
r
r
2
1
-
i
t
e
r
/
G
m
a
x
r
≥
0.5
w
i
j
i
t
e
r
+
w
m
i
n
i
t
e
r
-
w
i
j
i
t
e
r
r
2
1
-
i
t
e
r
/
G
m
a
x
r
<
0.5
(9) END FOR; END WHILE
(10) W
b
e
s
t
=
ω
j
k
ω
i
j
a
j
b
k /∗Initial weights and thresholds of BPNN by the obtained best individual. ∗/
(11) x
k
=
x
k
-
x
m
i
n
/
x
m
a
x
-
x
m
i
n, k
∈
[
1
,
q
] /∗ The process of training data normalization. ∗/
(12) i
t
e
r
=
0 /∗ Adjust the weights and threshold of BPNN according to the forecast error. ∗/
(13) DO WHILE  i
t
e
r
≤
i
t
e
r
m
a
x
(14)     i
t
e
r
=
i
t
e
r
+
1
(15)     FOR EACH  1
≤
i
≤
n  DO
(16)          FOR EACH  1
≤
j
≤
n
H
i
d
d
e
n  DO
(17)              H
j
=
f
∑
i
=
1
n
ω
i
j
x
i
-
a /∗ Calculate the outputs of the hidden layer. ∗/
(18)              FOR EACH  1
≤
k
≤
m  DO
(19)                  O
k
=
∑
j
=
1
l
H
j
ω
j
k
-
b
k /∗ Calculate the outputs of the output layer. ∗/
(20)                  e
k
=
y
n
k
-
O
k /∗ Calculate the forecast error in the output layer. ∗/
(21)                  ω
i
j
=
ω
i
j
+
η
H
j
1
-
H
j
x
i
∑
k
=
1
m
ω
j
k
e
k /∗ Update the connection weights ∗/
(22)                  ω
j
k
=
ω
j
k
+
η
H
j
e
k
(23)                  a
j
=
a
j
+
η
H
j
1
-
H
j
∑
k
=
1
m
ω
j
k
e
k  b
k
=
b
k
+
e
k /∗ Update the threshold. ∗/
(24) END FOR; END WHILE
(25) RETURN  x
^
f
=
(
x
^
(
q
+
1
)
,
x
^
(
q
+
2
)
,
…
,
x
^
(
q
+
p
)
)


### 2.5. EMD-CSAWNN

In this paper, the proposed model, which incorporates the empirical model decomposition technique into the WNN model based on CSA, is adopted for short-term power load forecasts. Empirical model decomposition represents a self-adaptive decomposition technique to decompose short-term power load series into several IMFs and one residual item. WNN is adopted as a forecasting engine in the proposed approach because of its powerful approximation and high computation speed. Additionally, to avoid the deficiencies of WNN such as its unstable structure, CSA is used to initialize and determine the weights and thresholds of WNN, thereby imparting an outreach capacity to WNN. Figure3 illustrates the general structure of the hybrid power load forecasting method. The pseudo-code for the algorithm of the EMD-CSAWNN model is as shown in Algorithm 2.Algorithm 2:EMD-CSAWNN.
(
1) /∗ Initialize S
D
0
=
1, r
0
t
=
x
t, h
k
-
1
t
=
r
j
-
1
t, k
=
1,j
=
1.∗/
(2) WHILE  r
j
t is not a monotonic function DO
(3)    WHILE  S
D
k
>
0.3  DO /∗ Find all local maxima and minima of h
k
-
1
t by cubic spline.∗/
(4)         /∗ Produce the upper and lower envelopes expressed as e
m
a
x
t and e
m
i
n
t.∗/
(5)         m
k
t
=
e
m
a
x
t
+
e
m
i
n
t
/
2, h
k
t
=
h
k
-
1
t
-
m
k
-
1
t, k
=
k
+
1
(6)    END WHILE
(7)    I
M
F
j
=
h
k
t, r
j
t
=
r
j
-
1
t
-
I
M
F
j, j
=
j
+
1
(8) END WHILE
(9) x
k
=
∑
i
=
2
n
I
M
F
i
+
r
n, k
∈
[
1
,
q
]
(10) GENERATE  i
t
e
r
=
0
(11) DO WHILE  i
t
e
r
≤
i
t
e
r
m
a
x
(12)     i
t
e
r
=
i
t
e
r
+
1 /∗ Generate a cuckoo egg X
i by taking a Lévy flight from a random nest.∗/
(13)     FOR EACH  1
≤
i
≤
n  DO/
∗Calculate the fitness value f
i
⋅ and select the candidate. ∗/
(14)        X
i
i
t
e
r
+
1
=
X
i
i
t
e
r
+
α
⊕
L
λ
(15)        IF  f
i
t
e
r
⋅
≤
f
i
t
e
r
+
1
⋅  DO /
∗ Replace the worse location of the nest by the better one.∗/
(16)          X
i
t
e
r
+
1
=
X
i
t
e
r
(17)       ELSE  X
i
t
e
r
+
1
=
X
i
t
e
r
+
1
(18)       rand = randn()
(19)           IF  r
a
n
d
>
P
α  DO /∗ Discard the worst solution by P
α and produce a new solution.∗/
(20)         X
i
i
t
e
r
+
1
=
X
i
i
t
e
r
+
r
χ
j
i
t
e
r
-
χ
k
i
t
e
r
(21)         ELSE  X
i
t
e
r
+
1
=
X
i
t
e
r
+
1
(22) END IF; END FOR; END WHILE
(23) X
b
e
s
t
=
ω
j
k
ω
i
j
a
j
b
j /∗Initial weights and thresholds of WNN by obtaining the best individual. ∗/
(24) x
k
=
x
k
-
x
m
i
n
/
x
m
a
x
-
x
m
i
n, k
∈
[
1
,
q
] /∗ The process of training data normalization. ∗/
(25) i
t
e
r
=
0 /∗ Adjust the weights and threshold of WNN according to the forecast error. ∗/
(26) DO WHILE  i
t
e
r
≤
i
t
e
r
m
a
x
(27)    i
t
e
r
=
i
t
e
r
+
1
(28)    FOR EACH  1
≤
i
≤
n  DO
(29)      FOR EACH  1
≤
j
≤
n
H
i
d
d
e
n  DO
(30)         h
j
=
h
j
∑
i
=
1
k
ω
i
j
x
i
-
b
j
/
a
j /∗ Calculate the outputs of the hidden layer. ∗/
(31)         FOR EACH  1
≤
k
≤
m  DO
(32)            y
k
=
∑
i
=
1
n
H
ω
i
j
h
i /∗ Calculate the outputs of the output layer. ∗/
(33)            e
k
=
y
n
k
-
y
k /∗ Calculate the forecast error in the output layer. ∗/
(34)            ω
i
+
1
,
k
=
ω
i
,
k
+
Δ
ω
i
+
1
,
k /∗ Update the connection weights ∗/
(35)            a
k
=
a
k
+
Δ
a
k  b
k
=
b
k
+
Δ
b
k /∗ Update the threshold.∗/
(36)            Δ
ω
i
+
1
,
k
=
-
η
∂
e
k
/
∂
ω
i
,
k  Δ
a
k
=
-
η
∂
e
k
/
∂
a
k  Δ
b
k
=
-
η
∂
e
k
/
∂
b
k
(37) END FOR; END WHILE
(38) RETURN  x
^
f
=
(
x
^
(
q
+
1
)
,
x
^
(
q
+
2
)
,
…
,
x
^
(
q
+
p
)
)Figure 3
The construction of EMD-CSAWNN. (a) IMFs and the residual decomposed by empirical model decomposition. (b) The structure of the training and testing sets of the neural network. (c) The procedure of WNN combined with CSA. (d) The procedure of CSA. (e) The structure of WNN. (f) The forecasting sequence and error by using the EMD-CSAWNN model.


## 2.1. Empirical Mode Decomposition

Empirical model decomposition was proposed by Karthikeyan and Kumar as an adaptive method for nonstationary time series analysis, and it is now widely used. It can be applied to any type of signal decomposition [35]. Thus, it has obvious advantages in processing nonstationary and nonlinear series. The foundation of this technique is to decompose a time series into a finite set of several IMFs and a residue [36].Definition 1.
All IMFs are defined to satisfy the following conditions: (1) the number of local extreme points and the number of zero crossings must be equal or at least differ by only one; (2) the mean value of the upper envelope and lower envelope is zero.Definition 2.
The stoppage criterion determined is defined as(1)
S
D
k
=
∑
t
=
0
T
h
1
k
-
1
t
-
h
1
k
t
2
∑
t
=
0
T
h
1
k
-
1
2
t
.
The sifting process stops whenS
D
k is smaller than a pregiven value. The process of decomposition is over when the value of S
D
k is between 0.2 and 0.3. Additional details of the empirical model decomposition technique are illustrated in Figure 1.Figure 1
The construction of empirical model decomposition.


## 2.2. Artificial Neural Network (ANN)

### 2.2.1. Confirmation of the Structure of the Network

The ANN has received considerable attention as a powerful computational tool for forecasting in many fields since 1980. ANN models always outperform statistical models because of their ability to map the inputs onto outputs via simple computation [37]. We discuss the feed-forward neural network in this paper because of its strong learning ability and simple structure. The determination of the network structure is as follows [38].Definition 3.
Given an arbitrary continuous functionf
:
[
0,1
]
n
→
R
m, f can be accurately approximated using a three-layer forward neural network realization. The first layer of the network is the input layer, containing n neurons. The middle layer is the hidden layer, containing 2
n
+
1 neurons. The third layer is the output layer, which has m neurons.Definition 4.
Letϕ be bounded continuous monotone in function, let K be the bounded compact subsets, and let K
⊂
R
n, f
X
=
f
x
1
,
x
2
,
…
,
x
n be a real continuous function in f
X
∈
K. ∀
ε
>
0, ∃ integer N, real constant C
i, θ
i
i
=
1,2
,
…
,
N and W
i
j
i
=
1,2
,
…
,
N, such that f
^
i
j
x
1
,
x
2
,
…
x
n
=
∑
i
=
1
N
C
i
ϕ
∑
j
=
1
n
W
i
j
x
j
-
θ
i satisfying max
⁡
f
^
x
1
,
…
,
x
n
-
f
x
1
,
…
,
x
n
<
ε. Note that, ∀
ε
>
0, ∃ a three-layer network structure, and the output function of the hidden layer is ϕ
x, the output function of the input and output layer is linear, and the total relationship of I
/
0 isf
x
1
,
…
,
x
n, such that max
⁡
f
^
x
1
,
…
,
x
n
-
f
x
1
,
…
,
x
n
<
ε.Proof.
Becausef
X
X
=
x
1
,
…
,
x
n is a continuous function, f
X
∈
K, K
⊂
R
n, assume that f
X is a bounded C
∞ function. Based on the Paley-Wiener Theorem [39], the Fourier transform F
W
W
=
W
1
,
…
,
W
n is a real analytic function of f
X, ∀ a constant C
N, such that F
W
≤
C
N
1
+
W
-
N. We define such a function: (2)
I
A
x
1
,
…
,
x
n
=
∫
-
A
A
⋯
∫
-
A
A
ϕ
∑
i
=
1
n
x
i
W
i
-
W
0
1
2
π
n
ϕ
1
F
W
1
,
…
,
W
n
exp
⁡
i
W
0
d
W
0
⋯
d
W
n
,
I
∞
,
A
x
1
,
…
,
x
n
=
∫
-
A
A
⋯
∫
-
A
A
∫
-
∞
∞
ϕ
∑
i
=
1
n
x
i
W
i
-
W
0
1
2
π
n
ϕ
1
F
W
1
,
…
,
W
n
exp
⁡
i
W
0
d
W
0
⋯
d
W
n
,
J
A
x
1
,
…
,
x
n
=
1
2
π
n
∫
-
∞
∞
⋯
∫
-
∞
∞
F
W
1
,
…
,
W
n
exp
⁡
i
∑
i
=
1
n
x
i
W
i
d
W
1
⋯
d
W
n
.
Using an estimation ofF
·, l
i
m
A
→
∞
J
A
x
1
,
…
,
x
n
=
f
x
1
,
…
,
x
n can be proved: that is, J
A is uniform convergence. Thus, l
i
m
A
→
∞
I
∞
,
A
x
1
,
…
,
x
n
=
f
x
1
,
…
,
x
n. That is, ∀
ε
>
0, ∃
A
>
0, such that m
a
x
x
∈
R
n
I
∞
,
A
x
1
,
…
,
x
n
-
f
x
1
,
…
,
x
n
<
E
/
2. Moreover, for A
′
>
0, let I
A
′
,
A
x
1
,
…
,
x
n
=
∫
-
A
A
⋯
∫
-
A
A
∫
-
A
′
A
ϕ
∑
i
=
1
n
x
i
W
i
-
W
0
·
(
1
/
2
π
n
ϕ
1
)
F
W
1
,
…
,
W
n
·
exp
⁡
i
W
0
d
W
0
⋯
d
W
n. For ε
>
0, A
′
>
0 can be found such that m
a
x
x
∈
K
I
A
,
A
x
1
,
…
,
x
n
-
I
∞
,
A
x
1
,
…
,
x
n
<
ε
/
2. That is to say, f
X is a uniform approximation by the integral I
A
′
,
A
X, I
A
′
,
A
X
∈
K. If A
i
>
0
i
=
1,2
,
…
,
m, K is a bounded subset, K
⊂
R
n, and h
i
x
1
,
…
,
x
m
;
t
1
,
…
,
t
m is a uniformly continuous function of -
A
1
,
A
1
×
⋯
×
-
A
m
,
A
m
×
K. There are functions: (3)
H
t
=
∫
-
A
1
A
1
⋯
∫
-
A
m
A
m
h
x
1
,
…
,
x
m
;
t
1
,
…
,
t
m
d
x
1
⋯
d
x
m
,
H
N
t
=
2
A
1
⋯
2
A
m
N
m
×
∑
k
1
⋯
k
m
=
0
N
-
1
h
-
A
1
+
k
1
·
2
A
1
N
⋯
A
m
+
k
m
·
2
A
m
N
t
1
⋯
t
m
.
∀
ε
>
0, ∃
N
0, N
0 is a natural number, when N
≥
N
0, such that(4)
max
l
∈
K
⁡
H
t
-
H
N
t
<
ε
.Thus, I
A
′
,
A
X is a uniform approximate on K, (5)
ϕ
∑
i
=
0
n
x
i
W
i
-
W
0
=
ϕ
∑
i
=
1
n
W
i
x
i
δ
-
W
0
+
α
-
ϕ
∑
i
=
1
n
W
i
x
i
δ
-
W
0
-
α
.The formula above can be achieved by a three-layer neural network. Thus, f
X can be approximated by that neural network.Definition 5.
∀
ε
>
0 and f
:
[
0,1
]
n
→
R
m, there exists a three-layer structure that can approximate f in any square error precision of ε.The definition above proves that,∀
f
:
0,1
n
→
R
m, we can use a feed-forward neural network with a three-layer structure n
×
2
n
+
1
×
m to approximate it accurately. Thus, this part not only proves the existence of the mapping network but also demonstrates the network structure of the mapping. In summary, this paper adopts the three-layer neural network as the basis neural network.


### 2.2.2. BPNN

BPNN is a type of multilayer feed-forward neural network with an error back propagation learning process. The structure of BPNN is illustrated in Figure2. Details of BPNN are introduced in [40].Figure 2
The construction of GABPNN and CSAWNN.


### 2.2.3. WNN

WNN, a feed-forward network, is generally multilayer [41]. It is widely applied in signal processing because of its advantages of the localization property and generalization ability [42]. The structure of WNN is shown in Figure 2.


## 2.2.1. Confirmation of the Structure of the Network

The ANN has received considerable attention as a powerful computational tool for forecasting in many fields since 1980. ANN models always outperform statistical models because of their ability to map the inputs onto outputs via simple computation [37]. We discuss the feed-forward neural network in this paper because of its strong learning ability and simple structure. The determination of the network structure is as follows [38].Definition 3.
Given an arbitrary continuous functionf
:
[
0,1
]
n
→
R
m, f can be accurately approximated using a three-layer forward neural network realization. The first layer of the network is the input layer, containing n neurons. The middle layer is the hidden layer, containing 2
n
+
1 neurons. The third layer is the output layer, which has m neurons.Definition 4.
Letϕ be bounded continuous monotone in function, let K be the bounded compact subsets, and let K
⊂
R
n, f
X
=
f
x
1
,
x
2
,
…
,
x
n be a real continuous function in f
X
∈
K. ∀
ε
>
0, ∃ integer N, real constant C
i, θ
i
i
=
1,2
,
…
,
N and W
i
j
i
=
1,2
,
…
,
N, such that f
^
i
j
x
1
,
x
2
,
…
x
n
=
∑
i
=
1
N
C
i
ϕ
∑
j
=
1
n
W
i
j
x
j
-
θ
i satisfying max
⁡
f
^
x
1
,
…
,
x
n
-
f
x
1
,
…
,
x
n
<
ε. Note that, ∀
ε
>
0, ∃ a three-layer network structure, and the output function of the hidden layer is ϕ
x, the output function of the input and output layer is linear, and the total relationship of I
/
0 isf
x
1
,
…
,
x
n, such that max
⁡
f
^
x
1
,
…
,
x
n
-
f
x
1
,
…
,
x
n
<
ε.Proof.
Becausef
X
X
=
x
1
,
…
,
x
n is a continuous function, f
X
∈
K, K
⊂
R
n, assume that f
X is a bounded C
∞ function. Based on the Paley-Wiener Theorem [39], the Fourier transform F
W
W
=
W
1
,
…
,
W
n is a real analytic function of f
X, ∀ a constant C
N, such that F
W
≤
C
N
1
+
W
-
N. We define such a function: (2)
I
A
x
1
,
…
,
x
n
=
∫
-
A
A
⋯
∫
-
A
A
ϕ
∑
i
=
1
n
x
i
W
i
-
W
0
1
2
π
n
ϕ
1
F
W
1
,
…
,
W
n
exp
⁡
i
W
0
d
W
0
⋯
d
W
n
,
I
∞
,
A
x
1
,
…
,
x
n
=
∫
-
A
A
⋯
∫
-
A
A
∫
-
∞
∞
ϕ
∑
i
=
1
n
x
i
W
i
-
W
0
1
2
π
n
ϕ
1
F
W
1
,
…
,
W
n
exp
⁡
i
W
0
d
W
0
⋯
d
W
n
,
J
A
x
1
,
…
,
x
n
=
1
2
π
n
∫
-
∞
∞
⋯
∫
-
∞
∞
F
W
1
,
…
,
W
n
exp
⁡
i
∑
i
=
1
n
x
i
W
i
d
W
1
⋯
d
W
n
.
Using an estimation ofF
·, l
i
m
A
→
∞
J
A
x
1
,
…
,
x
n
=
f
x
1
,
…
,
x
n can be proved: that is, J
A is uniform convergence. Thus, l
i
m
A
→
∞
I
∞
,
A
x
1
,
…
,
x
n
=
f
x
1
,
…
,
x
n. That is, ∀
ε
>
0, ∃
A
>
0, such that m
a
x
x
∈
R
n
I
∞
,
A
x
1
,
…
,
x
n
-
f
x
1
,
…
,
x
n
<
E
/
2. Moreover, for A
′
>
0, let I
A
′
,
A
x
1
,
…
,
x
n
=
∫
-
A
A
⋯
∫
-
A
A
∫
-
A
′
A
ϕ
∑
i
=
1
n
x
i
W
i
-
W
0
·
(
1
/
2
π
n
ϕ
1
)
F
W
1
,
…
,
W
n
·
exp
⁡
i
W
0
d
W
0
⋯
d
W
n. For ε
>
0, A
′
>
0 can be found such that m
a
x
x
∈
K
I
A
,
A
x
1
,
…
,
x
n
-
I
∞
,
A
x
1
,
…
,
x
n
<
ε
/
2. That is to say, f
X is a uniform approximation by the integral I
A
′
,
A
X, I
A
′
,
A
X
∈
K. If A
i
>
0
i
=
1,2
,
…
,
m, K is a bounded subset, K
⊂
R
n, and h
i
x
1
,
…
,
x
m
;
t
1
,
…
,
t
m is a uniformly continuous function of -
A
1
,
A
1
×
⋯
×
-
A
m
,
A
m
×
K. There are functions: (3)
H
t
=
∫
-
A
1
A
1
⋯
∫
-
A
m
A
m
h
x
1
,
…
,
x
m
;
t
1
,
…
,
t
m
d
x
1
⋯
d
x
m
,
H
N
t
=
2
A
1
⋯
2
A
m
N
m
×
∑
k
1
⋯
k
m
=
0
N
-
1
h
-
A
1
+
k
1
·
2
A
1
N
⋯
A
m
+
k
m
·
2
A
m
N
t
1
⋯
t
m
.
∀
ε
>
0, ∃
N
0, N
0 is a natural number, when N
≥
N
0, such that(4)
max
l
∈
K
⁡
H
t
-
H
N
t
<
ε
.Thus, I
A
′
,
A
X is a uniform approximate on K, (5)
ϕ
∑
i
=
0
n
x
i
W
i
-
W
0
=
ϕ
∑
i
=
1
n
W
i
x
i
δ
-
W
0
+
α
-
ϕ
∑
i
=
1
n
W
i
x
i
δ
-
W
0
-
α
.The formula above can be achieved by a three-layer neural network. Thus, f
X can be approximated by that neural network.Definition 5.
∀
ε
>
0 and f
:
[
0,1
]
n
→
R
m, there exists a three-layer structure that can approximate f in any square error precision of ε.The definition above proves that,∀
f
:
0,1
n
→
R
m, we can use a feed-forward neural network with a three-layer structure n
×
2
n
+
1
×
m to approximate it accurately. Thus, this part not only proves the existence of the mapping network but also demonstrates the network structure of the mapping. In summary, this paper adopts the three-layer neural network as the basis neural network.


## 2.2.2. BPNN

BPNN is a type of multilayer feed-forward neural network with an error back propagation learning process. The structure of BPNN is illustrated in Figure2. Details of BPNN are introduced in [40].Figure 2
The construction of GABPNN and CSAWNN.


## 2.2.3. WNN

WNN, a feed-forward network, is generally multilayer [41]. It is widely applied in signal processing because of its advantages of the localization property and generalization ability [42]. The structure of WNN is shown in Figure 2.


## 2.3. Neural Network Optimized by an Intelligence Algorithm

The intelligent optimization algorithm provides an efficient and powerful mathematical tool for optimizing the initial weights and thresholds of the ANN [43].


### 2.3.1. CSA

Cuckoo search is a heuristic swarm intelligence algorithm inspired by the behavior of the obligating brood parasitism of cuckoo species [44]. CSA is utilized in this paper for its stronger capability of global optimization [43].Definition 6.
To simulate the mode of cuckoo breeding, three idealized assumptions are presented, as follows: (1) each cuckoo selects nest randomly and dumps only one egg at a time, (2) the eggs with high quality will be carried over to the next generation, and (3) the available nest numbern is fixed, and the probability of the host bird discovering the exotic egg is p
α.Definition 7.
The Lévy flight model simulates the process of the nest-seeking characteristic of cuckoo, and the update formula of the path and location is as follows:(6)
X
i
i
t
e
r
+
1
=
X
i
i
t
e
r
+
α
⊕
L
λ
.
X
i
i
t
e
r represents the location of nest i at generation i
t
e
r, α is the step length-controlled factor, ⊕ is the point-to-point multiplication, and L
e
v
y
λ obeys the Lévy distribution with a random search path of parameter λ:(7)
L
λ
=
0.01
u
v
X
i
i
t
e
r
-
X
∗
i
t
e
r
. Here u
=
t
-
λ, 1
<
λ
<
3, u and v obey the normal distribution, u
~
N
0
,
σ
u
2, v
~
N
0,1, σ
u
=
Γ
λ
sin
⁡
0.5
π
λ
-
1
/
2
λ
-
2
/
2
Γ
0.5
λ
λ
-
1
1
/
λ
-
1, X
∗
i
t
e
r denotes the location of the best nest at generation i
t
e
r, and Γ is the standard gamma function with unbounded variance and mean of the probability distribution.


### 2.3.2. GA

GA is a population-based optimization algorithm that simulates natural genetic mechanisms and biological evolutionism. It possesses a capacity for powerful global optimization [45]. The principle of GA relies on a random process, which is constituted by the processes of selection, crossover, and mutation [46]. The implementation process is shown in Figure 2.


## 2.3.1. CSA

Cuckoo search is a heuristic swarm intelligence algorithm inspired by the behavior of the obligating brood parasitism of cuckoo species [44]. CSA is utilized in this paper for its stronger capability of global optimization [43].Definition 6.
To simulate the mode of cuckoo breeding, three idealized assumptions are presented, as follows: (1) each cuckoo selects nest randomly and dumps only one egg at a time, (2) the eggs with high quality will be carried over to the next generation, and (3) the available nest numbern is fixed, and the probability of the host bird discovering the exotic egg is p
α.Definition 7.
The Lévy flight model simulates the process of the nest-seeking characteristic of cuckoo, and the update formula of the path and location is as follows:(6)
X
i
i
t
e
r
+
1
=
X
i
i
t
e
r
+
α
⊕
L
λ
.
X
i
i
t
e
r represents the location of nest i at generation i
t
e
r, α is the step length-controlled factor, ⊕ is the point-to-point multiplication, and L
e
v
y
λ obeys the Lévy distribution with a random search path of parameter λ:(7)
L
λ
=
0.01
u
v
X
i
i
t
e
r
-
X
∗
i
t
e
r
. Here u
=
t
-
λ, 1
<
λ
<
3, u and v obey the normal distribution, u
~
N
0
,
σ
u
2, v
~
N
0,1, σ
u
=
Γ
λ
sin
⁡
0.5
π
λ
-
1
/
2
λ
-
2
/
2
Γ
0.5
λ
λ
-
1
1
/
λ
-
1, X
∗
i
t
e
r denotes the location of the best nest at generation i
t
e
r, and Γ is the standard gamma function with unbounded variance and mean of the probability distribution.


## 2.3.2. GA

GA is a population-based optimization algorithm that simulates natural genetic mechanisms and biological evolutionism. It possesses a capacity for powerful global optimization [45]. The principle of GA relies on a random process, which is constituted by the processes of selection, crossover, and mutation [46]. The implementation process is shown in Figure 2.


## 2.4. GABPNN

The primary mechanism of GABPNN is composed of three parts: GA optimization, determination of the BPNN structure, and forecasting covered by BPNN [47]. The pseudo-code for GABPNN is as shown in Algorithm 1.Algorithm 1:GABPNN.
(
1) GENERATE  i
t
e
r
=
0 /∗ Calculate the fitness of each individual in the population. ∗/
(2) DO WHILE i
t
e
r
≤
i
t
e
r
m
a
x
(3)    i
t
e
r
=
i
t
e
r
+
1 /∗ Record the best fitness values and the average fitness values. ∗/
(4)    FOR EACH  1
≤
i
≤
P  DO /∗The process of the select operation. ∗/
(5)         p
i
i
t
e
r
=
f
i
i
t
e
r
⋅
/
∑
i
=
1
P
f
i
i
t
e
r
⋅
(6)         FOR EACH  1
≤
j
≤
s  DO /∗ The process of the crossover and mutation operations. ∗/
(7)           w
k
j
i
t
e
r
=
w
k
j
i
t
e
r
1
-
b
+
w
l
j
i
t
e
r
b  w
l
j
i
t
e
r
=
w
l
j
i
t
e
r
1
-
b
+
w
k
j
i
t
e
r
b
(8)           w
i
j
i
t
e
r
=
w
i
j
i
t
e
r
+
w
i
j
i
t
e
r
-
w
m
a
x
i
t
e
r
r
2
1
-
i
t
e
r
/
G
m
a
x
r
≥
0.5
w
i
j
i
t
e
r
+
w
m
i
n
i
t
e
r
-
w
i
j
i
t
e
r
r
2
1
-
i
t
e
r
/
G
m
a
x
r
<
0.5
(9) END FOR; END WHILE
(10) W
b
e
s
t
=
ω
j
k
ω
i
j
a
j
b
k /∗Initial weights and thresholds of BPNN by the obtained best individual. ∗/
(11) x
k
=
x
k
-
x
m
i
n
/
x
m
a
x
-
x
m
i
n, k
∈
[
1
,
q
] /∗ The process of training data normalization. ∗/
(12) i
t
e
r
=
0 /∗ Adjust the weights and threshold of BPNN according to the forecast error. ∗/
(13) DO WHILE  i
t
e
r
≤
i
t
e
r
m
a
x
(14)     i
t
e
r
=
i
t
e
r
+
1
(15)     FOR EACH  1
≤
i
≤
n  DO
(16)          FOR EACH  1
≤
j
≤
n
H
i
d
d
e
n  DO
(17)              H
j
=
f
∑
i
=
1
n
ω
i
j
x
i
-
a /∗ Calculate the outputs of the hidden layer. ∗/
(18)              FOR EACH  1
≤
k
≤
m  DO
(19)                  O
k
=
∑
j
=
1
l
H
j
ω
j
k
-
b
k /∗ Calculate the outputs of the output layer. ∗/
(20)                  e
k
=
y
n
k
-
O
k /∗ Calculate the forecast error in the output layer. ∗/
(21)                  ω
i
j
=
ω
i
j
+
η
H
j
1
-
H
j
x
i
∑
k
=
1
m
ω
j
k
e
k /∗ Update the connection weights ∗/
(22)                  ω
j
k
=
ω
j
k
+
η
H
j
e
k
(23)                  a
j
=
a
j
+
η
H
j
1
-
H
j
∑
k
=
1
m
ω
j
k
e
k  b
k
=
b
k
+
e
k /∗ Update the threshold. ∗/
(24) END FOR; END WHILE
(25) RETURN  x
^
f
=
(
x
^
(
q
+
1
)
,
x
^
(
q
+
2
)
,
…
,
x
^
(
q
+
p
)
)


## 2.5. EMD-CSAWNN

In this paper, the proposed model, which incorporates the empirical model decomposition technique into the WNN model based on CSA, is adopted for short-term power load forecasts. Empirical model decomposition represents a self-adaptive decomposition technique to decompose short-term power load series into several IMFs and one residual item. WNN is adopted as a forecasting engine in the proposed approach because of its powerful approximation and high computation speed. Additionally, to avoid the deficiencies of WNN such as its unstable structure, CSA is used to initialize and determine the weights and thresholds of WNN, thereby imparting an outreach capacity to WNN. Figure3 illustrates the general structure of the hybrid power load forecasting method. The pseudo-code for the algorithm of the EMD-CSAWNN model is as shown in Algorithm 2.Algorithm 2:EMD-CSAWNN.
(
1) /∗ Initialize S
D
0
=
1, r
0
t
=
x
t, h
k
-
1
t
=
r
j
-
1
t, k
=
1,j
=
1.∗/
(2) WHILE  r
j
t is not a monotonic function DO
(3)    WHILE  S
D
k
>
0.3  DO /∗ Find all local maxima and minima of h
k
-
1
t by cubic spline.∗/
(4)         /∗ Produce the upper and lower envelopes expressed as e
m
a
x
t and e
m
i
n
t.∗/
(5)         m
k
t
=
e
m
a
x
t
+
e
m
i
n
t
/
2, h
k
t
=
h
k
-
1
t
-
m
k
-
1
t, k
=
k
+
1
(6)    END WHILE
(7)    I
M
F
j
=
h
k
t, r
j
t
=
r
j
-
1
t
-
I
M
F
j, j
=
j
+
1
(8) END WHILE
(9) x
k
=
∑
i
=
2
n
I
M
F
i
+
r
n, k
∈
[
1
,
q
]
(10) GENERATE  i
t
e
r
=
0
(11) DO WHILE  i
t
e
r
≤
i
t
e
r
m
a
x
(12)     i
t
e
r
=
i
t
e
r
+
1 /∗ Generate a cuckoo egg X
i by taking a Lévy flight from a random nest.∗/
(13)     FOR EACH  1
≤
i
≤
n  DO/
∗Calculate the fitness value f
i
⋅ and select the candidate. ∗/
(14)        X
i
i
t
e
r
+
1
=
X
i
i
t
e
r
+
α
⊕
L
λ
(15)        IF  f
i
t
e
r
⋅
≤
f
i
t
e
r
+
1
⋅  DO /
∗ Replace the worse location of the nest by the better one.∗/
(16)          X
i
t
e
r
+
1
=
X
i
t
e
r
(17)       ELSE  X
i
t
e
r
+
1
=
X
i
t
e
r
+
1
(18)       rand = randn()
(19)           IF  r
a
n
d
>
P
α  DO /∗ Discard the worst solution by P
α and produce a new solution.∗/
(20)         X
i
i
t
e
r
+
1
=
X
i
i
t
e
r
+
r
χ
j
i
t
e
r
-
χ
k
i
t
e
r
(21)         ELSE  X
i
t
e
r
+
1
=
X
i
t
e
r
+
1
(22) END IF; END FOR; END WHILE
(23) X
b
e
s
t
=
ω
j
k
ω
i
j
a
j
b
j /∗Initial weights and thresholds of WNN by obtaining the best individual. ∗/
(24) x
k
=
x
k
-
x
m
i
n
/
x
m
a
x
-
x
m
i
n, k
∈
[
1
,
q
] /∗ The process of training data normalization. ∗/
(25) i
t
e
r
=
0 /∗ Adjust the weights and threshold of WNN according to the forecast error. ∗/
(26) DO WHILE  i
t
e
r
≤
i
t
e
r
m
a
x
(27)    i
t
e
r
=
i
t
e
r
+
1
(28)    FOR EACH  1
≤
i
≤
n  DO
(29)      FOR EACH  1
≤
j
≤
n
H
i
d
d
e
n  DO
(30)         h
j
=
h
j
∑
i
=
1
k
ω
i
j
x
i
-
b
j
/
a
j /∗ Calculate the outputs of the hidden layer. ∗/
(31)         FOR EACH  1
≤
k
≤
m  DO
(32)            y
k
=
∑
i
=
1
n
H
ω
i
j
h
i /∗ Calculate the outputs of the output layer. ∗/
(33)            e
k
=
y
n
k
-
y
k /∗ Calculate the forecast error in the output layer. ∗/
(34)            ω
i
+
1
,
k
=
ω
i
,
k
+
Δ
ω
i
+
1
,
k /∗ Update the connection weights ∗/
(35)            a
k
=
a
k
+
Δ
a
k  b
k
=
b
k
+
Δ
b
k /∗ Update the threshold.∗/
(36)            Δ
ω
i
+
1
,
k
=
-
η
∂
e
k
/
∂
ω
i
,
k  Δ
a
k
=
-
η
∂
e
k
/
∂
a
k  Δ
b
k
=
-
η
∂
e
k
/
∂
b
k
(37) END FOR; END WHILE
(38) RETURN  x
^
f
=
(
x
^
(
q
+
1
)
,
x
^
(
q
+
2
)
,
…
,
x
^
(
q
+
p
)
)Figure 3
The construction of EMD-CSAWNN. (a) IMFs and the residual decomposed by empirical model decomposition. (b) The structure of the training and testing sets of the neural network. (c) The procedure of WNN combined with CSA. (d) The procedure of CSA. (e) The structure of WNN. (f) The forecasting sequence and error by using the EMD-CSAWNN model.


## 3. Experiments and Evaluations

Applications of the proposed hybrid approach and five comparison models are shown in this subsection. All algorithms are operated on the given platform: 3.20 GHz CPU, 8.00 GB RAM, Windows 7, and MATLAB R2012a. Meanwhile, taking into account the randomness factors and to make sure the final results are reliable and independent from the initial weights, we carry out each ANN experiment 50 times and then take the average value.


### 3.1. Region Description and Data Collection

Australia has plentiful coal and wind resources across its coastline and land. A power load data set from NSW, which is the state with the largest population and highest levels of industrialization and urbanization in Australia, is employed to validate each model. In this paper, power load data are collected randomly from NSW from January 12, 2009, to March 8, 2009, which includes eight weeks of data. Among them, the data from January 12, 2009, to February 1, 2009, is used as a training set to obtain the appropriate model, and data from February 2, 2009, to March 1, 2009, is used as the testing set. The power load data from January 12, 2009, to February 1, 2009, together with their statistical measures, that is, minimum, maximum, mean, and standard deviations, are shown in Figure4(c); the standard deviations are all above 1800, which implies that the power load series fluctuates significantly with the minimum/maximum of week one, week two, and week three, which are 6049.94/13518.06 (MWh), 6280.27/13326.47, and 6375.84/13096.37, respectively. This can be intuitively observed from the amplitude and frequency of the series fluctuation, which can change from very high to low values and vice versa.Figure 4
Description of observations in New South Wales of China. (a) Location of the study site. (b) The original power load series from January 12 to February 1, 2009. (c) The statistical measures for the power load.
(a)
(b)
(c)Figure4(b) exhibits 1008 power load data points from January 12, 2009, to February 1, 2009, divided into three groups, with 336 data points in every group. Because the power load data of NSW are collected once every half hour, each day includes 48 data points. On different days of the week, daily life and human economic production usually have different behaviors; thus, the characteristics of the load are different on different days. To minimize forecasting error as much as possible, we forecast the load of different days in the week separately. In this paper, the cycles of data division are seven days; the first three weeks of Monday data, January 12, January 19, and January 26, are employed to forecast the next Monday load on February 2, 2009. Accordingly, the data on January 13, January 20, and January 27 are employed to forecast the power load on February 3, 2009. The rest can be conducted in the same manner. The structures of the training and testing sets are illustrated in Figure 3(b).


### 3.2. Evaluation Metrics

Forecasting accuracy is an important criterion for evaluating a forecasting model. In this paper, the basic error calculation method is as follows:(8)
M
A
E
=
1
T
∑
t
=
1
T
A
E
t
=
1
T
∑
t
=
1
T
x
t
-
x
^
t
. Here, T is the number of data points; the formula A
E represents the absolute error between the observed value and the forecasting value at time t, x is the observed value, and x
^
t is the forecasting value at time t. To avoid a positive or negative offset in forecasting error, solve the problem for which positive and negative forecasting error cannot be added, adding to the absolute value of the error, and take the average in the end. This error belongs to the comprehensive index in error analysis:(9)
R
M
S
E
=
1
T
∑
t
=
1
T
A
E
t
2
=
1
T
∑
t
=
1
T
x
t
-
x
^
t
2
.RMSE is the square root of the mean square error. It also belongs to the comprehensive index in error analysis. Take the square of the absolute error AE; thus, the role of great values in the error will be strengthened, improving the sensitivity of the indicators, which is the prevailing reason in the error analysis:(10)
M
A
P
E
=
1
T
∑
t
=
1
T
x
t
-
x
^
t
x
t
×
100
%
,where the symbolic meaning is as above. The indicator is the average of the absolute error. The index is one of the comprehensive indexes in error analysis that usually occupies a very important position in the analysis and forecasting performance of the model.To determine the degree of correlation between different forecasting model results with observed values, GRA [48] is employed in this paper.Definition 8 (determining the grey relational coefficient).
x
0
i
i
=
1
n is the reference sequence, x
j
i
i
=
1
n is the comparative sequence, j
=
1,2
,
…
,
m, and the relational coefficient of x
0 and x
j in point k is represented as(11)
ξ
j
k
=
m
i
n
j
m
i
n
k
x
0
k
-
x
j
k
+
ρ
m
a
x
j
m
a
x
k
x
0
k
-
x
j
k
x
0
k
-
x
j
k
+
ρ
m
a
x
j
m
a
x
k
x
0
k
+
x
j
k
.Here, ρ is the distinguishing coefficient, and ρ
∈
0,1, usually ρ
=
0.5.Definition 9 (grey relational degree).
By focusing the degree ofξ
j
k at utter points, the algorithm on the grey relational degree is(12)
r
i
=
1
n
∑
k
=
1
n
ξ
i
k
.
Considering the generation capacity of the proposed hybrid model, four statistical indices are employed as evaluation metrics to measure the forecasting accuracy, MAE, RMSE, MAPE, and GRA. MAPE, MAE, and RMSE measure the mean performance, and GRA illustrates how well the forecasted data points fit the trend.


## 3.1. Region Description and Data Collection

Australia has plentiful coal and wind resources across its coastline and land. A power load data set from NSW, which is the state with the largest population and highest levels of industrialization and urbanization in Australia, is employed to validate each model. In this paper, power load data are collected randomly from NSW from January 12, 2009, to March 8, 2009, which includes eight weeks of data. Among them, the data from January 12, 2009, to February 1, 2009, is used as a training set to obtain the appropriate model, and data from February 2, 2009, to March 1, 2009, is used as the testing set. The power load data from January 12, 2009, to February 1, 2009, together with their statistical measures, that is, minimum, maximum, mean, and standard deviations, are shown in Figure4(c); the standard deviations are all above 1800, which implies that the power load series fluctuates significantly with the minimum/maximum of week one, week two, and week three, which are 6049.94/13518.06 (MWh), 6280.27/13326.47, and 6375.84/13096.37, respectively. This can be intuitively observed from the amplitude and frequency of the series fluctuation, which can change from very high to low values and vice versa.Figure 4
Description of observations in New South Wales of China. (a) Location of the study site. (b) The original power load series from January 12 to February 1, 2009. (c) The statistical measures for the power load.
(a)
(b)
(c)Figure4(b) exhibits 1008 power load data points from January 12, 2009, to February 1, 2009, divided into three groups, with 336 data points in every group. Because the power load data of NSW are collected once every half hour, each day includes 48 data points. On different days of the week, daily life and human economic production usually have different behaviors; thus, the characteristics of the load are different on different days. To minimize forecasting error as much as possible, we forecast the load of different days in the week separately. In this paper, the cycles of data division are seven days; the first three weeks of Monday data, January 12, January 19, and January 26, are employed to forecast the next Monday load on February 2, 2009. Accordingly, the data on January 13, January 20, and January 27 are employed to forecast the power load on February 3, 2009. The rest can be conducted in the same manner. The structures of the training and testing sets are illustrated in Figure 3(b).


## 3.2. Evaluation Metrics

Forecasting accuracy is an important criterion for evaluating a forecasting model. In this paper, the basic error calculation method is as follows:(8)
M
A
E
=
1
T
∑
t
=
1
T
A
E
t
=
1
T
∑
t
=
1
T
x
t
-
x
^
t
. Here, T is the number of data points; the formula A
E represents the absolute error between the observed value and the forecasting value at time t, x is the observed value, and x
^
t is the forecasting value at time t. To avoid a positive or negative offset in forecasting error, solve the problem for which positive and negative forecasting error cannot be added, adding to the absolute value of the error, and take the average in the end. This error belongs to the comprehensive index in error analysis:(9)
R
M
S
E
=
1
T
∑
t
=
1
T
A
E
t
2
=
1
T
∑
t
=
1
T
x
t
-
x
^
t
2
.RMSE is the square root of the mean square error. It also belongs to the comprehensive index in error analysis. Take the square of the absolute error AE; thus, the role of great values in the error will be strengthened, improving the sensitivity of the indicators, which is the prevailing reason in the error analysis:(10)
M
A
P
E
=
1
T
∑
t
=
1
T
x
t
-
x
^
t
x
t
×
100
%
,where the symbolic meaning is as above. The indicator is the average of the absolute error. The index is one of the comprehensive indexes in error analysis that usually occupies a very important position in the analysis and forecasting performance of the model.To determine the degree of correlation between different forecasting model results with observed values, GRA [48] is employed in this paper.Definition 8 (determining the grey relational coefficient).
x
0
i
i
=
1
n is the reference sequence, x
j
i
i
=
1
n is the comparative sequence, j
=
1,2
,
…
,
m, and the relational coefficient of x
0 and x
j in point k is represented as(11)
ξ
j
k
=
m
i
n
j
m
i
n
k
x
0
k
-
x
j
k
+
ρ
m
a
x
j
m
a
x
k
x
0
k
-
x
j
k
x
0
k
-
x
j
k
+
ρ
m
a
x
j
m
a
x
k
x
0
k
+
x
j
k
.Here, ρ is the distinguishing coefficient, and ρ
∈
0,1, usually ρ
=
0.5.Definition 9 (grey relational degree).
By focusing the degree ofξ
j
k at utter points, the algorithm on the grey relational degree is(12)
r
i
=
1
n
∑
k
=
1
n
ξ
i
k
.
Considering the generation capacity of the proposed hybrid model, four statistical indices are employed as evaluation metrics to measure the forecasting accuracy, MAE, RMSE, MAPE, and GRA. MAPE, MAE, and RMSE measure the mean performance, and GRA illustrates how well the forecasted data points fit the trend.


## 4. Results and Analysis

The experiments were divided into three parts, Experiment1, Experiment 2, Experiment 3, Experiment 4, and Experiment 5. Experiment 1, the primary experiment, aims to compare the performance of different models in one day (February 2, 2009) and, meanwhile, to compare the ARIMA and ANN models. Experiment 2 aims to compare the performance of different models between weekdays and weekends. Experiment 3 aims to compare the performance of different models in 28 days (from February 2, 2009, to March 1, 2009). In order to further validate the hybrid model we proposed Experiments 4 and 5 which are studied in this paper. In addition, ENN (Elman neural network) and RBFNN (radial basis function neural network) are utilized in further comparison. Experiment 4 aims at globally testing the proposed hybrid model by using real power load data on Thursday, April 24, 2008, and Tuesday, April 29, 2008, which is selected from New South Wales of Australia randomly. Experiment 5 aims to prove the general applicability of hybrid model; thus, the power load data on Saturday, June 28, 2008, and Monday, June 30, 2008, from Victoria of Australia are chosen for forecasting.Experiment 1.
According to determination of the network structure,n
×
2
n
+
1
×
m as the structure of the three-layer neural network will achieve better forecasting results. The network structure of the ANN model is 4-9-1 because of its four nodes in the input layer and one node in the output layer. The iteration time is set to 100, the learning rate is set to 0.01, and the training requirement accuracy is set to 0.0004. To forecast the power load on Monday, February 2, the historical values from the Mondays of the first three weeks, January 12, January 19, and January 26, are chosen. Test results (MAE, RMSE, MAPE, and GRA) are presented in Tables 1 and 2, Figures 5, 6, and 7. Each individual model exhibits its best performance at a special time. For example, Figure 6 shows that BPNN provides the lowest MAPE value at 3:00 among all of the individual models, and CSAWNN yields the highest accuracy forecasting value from 1:00 to 3:00 among all of the individual models, whereas the maximum error is with the BPNN forecasting model on February 2, with a MAPE value of 10.12%. This result is due to the unstable initialization of the ANN. The result of the original data of empirical model decomposition is shown in Figure 3(a). The noise in the data is eliminated by using the empirical model decomposition technique; in this paper, the IMF1 is a high-frequency sequence with small values, which can be regarded as interference factors. As a result, the rest of the IMFs and the remainder term can be constructed as the training input of the CSAWNN model. Table 1 shows the experimental results for Monday of six types of forecasting models. The average values of MAPE for six models on February 2, 1.6% are 1.54%, 1.97%, 1.37%, 1.02%, and 1.94%, respectively; as shown in Table 2, the MAPE afforded by the hybrid model decreased by 48.22% compared with the maximum average value of MAPE. In addition, it is shown that the value of MAPE offered by the hybrid model is more stable than that of the other proposed models, and the maximum value of MAPE is 2.54%. By comparing the hybrid model with the other models, it is shown that the hybrid model can provide high and stable forecasting accuracy.Table 1
Comparison of power load forecasting result by using different methods in Feb. 2.
Time
Actual value (MWh)
BPNN
GABPNN
WNN
CSAWNN
EMD-CSAWNN
ARIMA
Forecasting
MAPE
Forecasting
MAPE
Forecasting
MAPE
Forecasting
MAPE
Forecasting
MAPE
Forecasting
MAPE
Value
(%)
Value
(%)
Value
(%)
Value
(%)
Value
(%)
Value
(%)
0:00
8086.20
7267.78
10.12
7471.59
7.60
7828.93
3.18
7708.89
4.67
8033.14
0.66
8146.33
0.74
2:00
7132.07
7381.04
3.49
7193.94
0.87
7303.83
2.41
7150.74
0.26
7013.90
1.66
7298.17
2.33
4:00
7038.50
7023.50
0.21
7010.09
0.40
6470.25
8.07
6938.94
1.41
6938.97
1.41
6899.89
1.97
6:00
8803.29
8561.04
2.75
8611.39
2.18
8159.83
7.31
8612.24
2.17
8668.60
1.53
8619.45
2.09
8:00
10646.56
10778.68
1.24
10416.06
2.16
10890.92
2.30
10474.09
1.62
10707.41
0.57
10692.11
0.43
10:00
11822.68
11877.37
0.46
11811.56
0.09
11975.11
1.29
11880.45
0.49
11787.15
0.30
12002.49
1.52
12:00
12397.22
12329.59
0.55
12453.22
0.45
12328.77
0.55
12432.34
0.28
12326.00
0.57
12597.51
1.62
14:00
12824.82
12638.87
1.45
12918.86
0.73
12708.23
0.91
12929.78
0.82
12599.61
1.76
13091.61
2.08
16:00
13088.14
12810.33
2.12
13065.94
0.17
12903.89
1.41
13171.34
0.64
12830.13
1.97
12949.91
1.06
18:00
12001.48
11857.34
1.20
12090.61
0.74
12148.48
1.22
11500.36
4.18
11987.11
0.12
11776.37
1.88
20:00
10898.97
10827.46
0.66
10896.10
0.03
10967.20
0.63
10679.01
2.02
10778.06
1.11
10705.75
1.77
22:00
9398.37
9516.59
1.26
9665.37
2.84
9415.90
0.19
9706.03
3.27
9451.51
0.57
9446.46
0.51Table 2
Comparison of power load forecasting result by using different methods in Feb. 2.
Model
MAPE (%)
Maximum
Minimum
Average
BPNN
10.12
0.05
1.6
GABPNN
7.6
0.03
1.54
WNN
8.07
0.12
1.97
CSAWNN
5.19
0.01
1.37
EMD-CSAWNN
2.54
0.03
1.02
ARIMA
8.67
0.03
1.94Figure 5
The forecasting results and actual values for February 2.Figure 6
The BPNN, GABPNN, WNN, CSAWNN, EMD-CSAWNN, and ARIMA forecasting MAPE (%) for February 2.Figure 7
The comprehensive evaluation of forecasting models for February 2. (a) The comparison of GRA by using six different models. (b) The comparison of MAE, RMSE, and MAPE by using six different models. (c) The scatterplot of forecasting versus actual levels by using six different models. (d) The evaluation results of six different models. The value in bold is the best value of each evaluation index.Experiment 2.
Figure8 and Table 3 describe the comparison of six models with values on weekdays and weekends on different evaluation metrics. On weekdays, the best performance model is the hybrid model and the value of MAPE is 0.82%; on the contrary, the worst is the ARIMA model, whose value of MAPE is 1.48%. The MAPE offered by the hybrid model is 44.59% lower than that offered by the ARIMA model. On weekends, the value of MAPE offered by the hybrid model is 0.68%, which outperforms all other models and is 49.25% lower than that of the worst model, ARIMA. Figure 8 indicates that, except for the criterion RMSE, the weekdays outperform the weekends when utilizing WNN and the criterion GRA, the weekdays outperform the weekends when utilizing the BPNN, GABPNN, WNN, and CSAWNN models, and the other performances on weekends are all better than those on weekdays when utilizing the six models. The values of MAPE offered by the six models on weekends are 1.23%, 1.05%, 1.27%, 1.13%, 0.68%, and 1.34%, which are 5.38%, 10.25%, 0.78%, 8.13%, 17.07%, and 9.46% lower than the corresponding weekday values, respectively. This illustrates that the forecasting results on weekends outperform those of weekdays.Table 3
Comparison of power load forecasting evaluation by using different methods between weekdays and weekends.
Model
Date
Evaluation criterion
MAPE (%)
RMSE
MAE
GRA
BPNN
Weekdays
1.30
164.97
118.03
0.8654
Weekends
1.23
148.35
101.40
0.8600
GABPNN
Weekdays
1.17
154.94
106.05
0.8781
Weekends
1.05
137.21
88.33
0.8779
WNN
Weekdays
1.28
162.09
114.33
0.8688
Weekends
1.27
162.83
108.65
0.8581
CSAWNN
Weekdays
1.23
157.05
109.05
0.8737
Weekends
1.13
138.86
94.96
0.8689
EMD-CSAWNN
Weekdays
0.82
101.56
72.77
0.9087
Weekends
0.68
87.00
55.19
0.9148
ARIMA
Weekdays
1.48
184.23
128.98
0.8539
Weekends
1.34
169.06
109.09
0.8548Figure 8
The BPNN, GABPNN, WNN, CSAWNN, EMD-CSAWNN, and ARIMA forecasting results between weekdays and weekends.Experiment 3.
Figures9 and 10 illustrate the forecasting results of Experiment 3. In comparison with the BPNN, WNN, GABPNN, CSWNN, and ARIMA models, the forecasting results offered by the hybrid model are more accurate. The forecasting results are shown in Figure 10, and detailed prediction results are shown in Table 4. The maximum values of MAPE offered by the six models are 1.28%, 1.14%, 1.27%, 1.2%, 0.78%, and 1.44%, respectively. The maximum value of MAPE offered by the ARIMA model is 4.4% over 28 days. Meanwhile, on all days of the test, the average values of MAE, MAPE, and RMSE offered by the hybrid model are all smaller than those of the other models. The average value of MAPE offered by the hybrid model over four weeks is 0.78%, and the highest decrease is 39.06% compared with the other ANN models. This indicates that the hybrid model is an effective power load forecasting approach. The GRA result is shown in Table 5. In addition, on March 1, the GRA of GABPNN is higher than that of the hybrid model; on the remaining 27 days, the values of GRA offered by the hybrid model are higher than those of the other five models. According to the average value of GRA over the 28 days, the forecasting effects of all six forecasting models are increasing in the following order: ARIMA, BPNN, which is the WNN, CSAWNN, GABPNN, and EMD-CSAWNN, which concludes that the effect of the hybrid forecasting model is the best model among the six forecasting models.Table 4
Comparison of power load forecasting evaluation by using different methods from Feb. 2 to Mar. 1.
Time
BPNN
GABPNN
WNN
CSAWNN
EMD-CSAWNN
ARIMA
MAE
MAPEE (%)
RMSE
MAE
MAPE (%)
RMSE
MAE
MAPE (%)
RMSE
MAE
MAPE (%)
RMSE
MAE
MAPE (%)
RMSE
MAE
MAPE (%)
RMSE
Feb. 2
334.27
4.22
211.53
304.59
3.85
196.98
278.96
3.55
233.05
280.07
3.55
186.08
190.32
2.41
130.64
344.95
4.40
225.97
Feb. 3
126.81
1.70
126.63
96.89
1.29
125.59
114.76
1.55
116.64
92.33
1.25
127.49
100.06
1.35
78.43
184.62
2.48
128.64
Feb. 4
115.99
1.67
117.88
106.40
1.53
121.43
121.62
1.75
135.78
89.46
1.30
116.67
100.59
1.46
85.54
97.21
1.40
134.84
Feb. 5
62.90
0.93
209.87
55.75
0.82
177.65
73.07
1.09
141.15
83.43
1.26
142.02
63.85
0.96
115.75
88.74
1.32
169.61
Feb. 6
71.29
1.07
205.84
57.27
0.85
165.60
70.96
1.06
139.43
81.86
1.24
185.20
57.25
0.85
98.40
88.15
1.32
184.44
Feb. 7
84.41
1.21
226.10
60.33
0.86
198.48
86.40
1.24
285.05
82.14
1.19
203.72
75.10
1.07
175.72
121.60
1.75
211.26
Feb. 8
130.54
1.71
123.10
101.31
1.33
145.93
142.19
1.85
123.60
128.98
1.70
120.03
65.08
0.87
67.65
262.37
3.47
124.33
Feb. 9
118.07
1.43
150.69
108.27
1.31
146.03
115.14
1.38
169.23
104.84
1.28
134.63
70.86
0.84
97.30
107.14
1.31
180.04
Feb. 10
107.55
1.21
161.91
96.10
1.06
139.39
128.24
1.48
156.85
123.78
1.39
142.14
61.45
0.69
96.35
145.55
1.70
182.63
Feb. 11
93.31
1.01
192.90
111.84
1.17
177.37
95.98
1.03
172.67
118.18
1.29
219.21
53.49
0.58
83.04
118.46
1.30
221.91
Feb. 12
109.95
1.14
226.38
104.93
1.08
194.03
100.88
1.04
209.32
130.04
1.34
221.86
37.33
0.39
98.76
89.20
0.91
236.75
Feb. 13
81.85
0.82
253.47
58.74
0.59
231.68
92.82
0.94
249.81
79.82
0.79
243.81
43.70
0.44
127.79
61.35
0.62
310.47
Feb. 14
64.30
0.62
141.32
52.72
0.51
124.09
73.70
0.71
152.33
69.55
0.66
124.22
35.49
0.35
66.87
58.42
0.56
176.15
Feb. 15
55.99
0.56
99.81
39.82
0.40
93.03
50.17
0.49
102.20
55.17
0.55
120.02
32.92
0.33
57.69
47.36
0.46
118.80
Feb. 16
79.41
0.77
157.56
66.18
0.63
151.21
71.83
0.66
174.28
61.84
0.60
147.61
40.62
0.39
151.27
63.28
0.60
183.95
Feb. 17
82.06
0.77
116.51
59.34
0.55
147.71
82.63
0.73
128.10
56.92
0.54
131.90
46.33
0.41
97.02
54.36
0.50
142.02
Feb. 18
77.09
0.71
134.92
60.71
0.57
123.87
72.46
0.68
123.32
46.97
0.45
129.04
57.85
0.53
91.79
66.87
0.62
148.78
Feb. 19
96.09
0.89
132.50
87.60
0.82
129.83
99.07
0.92
132.46
83.92
0.79
132.63
58.39
0.52
83.16
71.12
0.68
171.76
Feb. 20
96.07
0.88
146.61
90.49
0.84
140.51
93.47
0.87
132.77
68.44
0.66
104.60
87.30
0.77
93.94
52.59
0.49
148.44
Feb. 21
97.27
0.88
186.13
82.81
0.74
185.30
107.84
0.98
192.47
76.45
0.69
179.31
77.07
0.70
48.00
95.04
0.87
263.72
Feb. 22
131.13
1.28
111.39
109.41
1.11
105.24
128.56
1.26
117.27
93.69
0.95
93.12
61.20
0.61
53.84
138.42
1.38
105.10
Feb. 23
103.86
0.98
132.92
74.92
0.72
150.01
93.02
0.92
168.00
89.48
0.86
173.36
62.09
0.62
93.66
113.45
1.17
163.20
Feb. 24
156.24
1.59
122.25
146.41
1.49
148.16
132.83
1.34
139.12
146.53
1.45
117.79
67.67
0.70
94.62
135.60
1.37
150.05
Feb. 25
152.98
1.59
125.02
146.55
1.51
144.41
162.06
1.67
120.09
146.17
1.50
126.54
99.86
1.05
87.54
174.77
1.82
152.51
Feb. 26
177.81
1.94
128.21
153.29
1.66
118.80
183.35
2.00
157.75
188.12
2.06
141.24
69.22
0.76
82.76
188.15
2.05
192.72
Feb. 27
162.42
1.84
145.73
184.00
2.08
110.68
191.75
2.19
162.54
173.50
1.95
122.48
67.46
0.77
110.58
177.14
2.01
159.88
Feb. 28
100.06
1.18
147.69
120.10
1.41
122.11
106.96
1.25
150.10
96.13
1.12
132.03
55.77
0.64
90.29
121.94
1.43
174.16
Mar. 1
102.17
1.27
104.78
90.90
1.13
74.49
85.07
1.06
92.83
92.96
1.16
101.33
58.62
0.73
62.75
184.40
2.27
111.80
Maximum value
334.27
4.22
253.47
304.59
3.85
231.68
278.96
3.55
285.05
280.07
3.55
243.81
190.32
2.41
175.72
344.95
4.40
310.47
Minimum value
55.99
0.56
99.81
39.82
0.40
74.49
50.17
0.49
92.83
46.97
0.45
93.12
32.92
0.33
48.00
47.36
0.46
105.10
Average value
113.28
1.28
160.39
100.99
1.14
150.09
112.71
1.27
162.30
105.03
1.20
152.07
67.75
0.78
97.62
123.29
1.44
180.03Table 5
Comparison of power load forecasting GRA by using different methods from Feb. 2 to Mar. 1.
Time
BPNN
GABPNN
WNN
CSAWN
EMD-CSAWNN
ARIMA
GRA
GRA
GRA
GRA
GRA
GRA
Feb. 2
0.7962
0.8018
0.7582
0.8186
0.8470
0.7553
Feb. 3
0.7587
0.7532
0.7756
0.7626
0.8121
0.7468
Feb. 4
0.7132
0.7077
0.7216
0.7429
0.7611
0.7076
Feb. 5
0.7225
0.7381
0.7672
0.7717
0.8032
0.7305
Feb. 6
0.7660
0.7991
0.8194
0.7800
0.8598
0.8027
Feb. 7
0.8555
0.8680
0.7952
0.8431
0.8839
0.8381
Feb. 8
0.7507
0.7184
0.7562
0.7504
0.8191
0.7541
Feb. 9
0.7622
0.7682
0.7495
0.7840
0.8352
0.7554
Feb. 10
0.7644
0.8003
0.7750
0.7993
0.8376
0.7573
Feb. 11
0.8297
0.8608
0.8479
0.8035
0.9109
0.8032
Feb. 12
0.7949
0.8509
0.8439
0.8074
0.9018
0.8169
Feb. 13
0.8262
0.8606
0.8396
0.8493
0.9050
0.8156
Feb. 14
0.7628
0.8064
0.7800
0.8067
0.8856
0.8058
Feb. 15
0.7392
0.7565
0.7716
0.7674
0.8262
0.7317
Feb. 16
0.8147
0.8298
0.8093
0.8291
0.8520
0.7859
Feb. 17
0.7927
0.7751
0.7656
0.7755
0.7964
0.7577
Feb. 18
0.7908
0.8114
0.7938
0.8043
0.8400
0.7791
Feb. 19
0.8254
0.8199
0.8227
0.8234
0.8756
0.7884
Feb. 20
0.7249
0.7719
0.7641
0.7901
0.7907
0.7340
Feb. 21
0.8555
0.8756
0.8473
0.8676
0.9334
0.8174
Feb. 22
0.7267
0.7600
0.7412
0.7624
0.8356
0.7547
Feb. 23
0.8258
0.8255
0.7912
0.7951
0.8670
0.7872
Feb. 24
0.7511
0.7630
0.7404
0.7642
0.8316
0.7363
Feb. 25
0.7932
0.7916
0.8259
0.8200
0.8331
0.7758
Feb. 26
0.7628
0.7983
0.7316
0.7402
0.8058
0.7063
Feb. 27
0.7440
0.7870
0.7203
0.7615
0.8257
0.7306
Feb. 28
0.8012
0.8331
0.7970
0.8011
0.8666
0.7607
Mar. 1
0.7681
0.8433
0.7974
0.7667
0.8214
0.7775
Maximum value
0.8555
0.8756
0.8479
0.8676
0.9334
0.8381
Minimum value
0.7132
0.7077
0.7203
0.7402
0.7611
0.7063
Average value
0.8673
0.8811
0.8692
0.8756
0.9127
0.8576Figure 9
The forecasting results and actual values from February 2 to March 1.Figure 10
The comprehensive evaluation of forecasting models from February 2 to March 1. (a) The radar diagram of MAPE by using six different models. (b) Comparison of forecasting MAPE by using the EMD-CSAWNN and ARIMA models. The red line is the polynomial regression line. (c) The comparison of GRA and MAPE by using six different models. (d) The comparison of MAE and RMSE by using six different models. (e) The scatterplot of forecasting versus actual levels by using six different models. The solid line represents the perfect fit: that is,y
=
x. (f) Evaluation results of six different models. The red font is the best value of every evaluation index; the green font is the worst. (g) The comparison of the box plot by using six different models. The whiskers in the box plot indicate the primary range for the data, in which the lowest data are 1.5 times the interquartile range of the lower quartile and the highest data are 1.5 times the interquartile range of the upper quartile. The outliers, which are not included between the whiskers, are represented by the red crosshair.The higher the power load forecasting accuracy, the lower the economic cost, which has actual economic significance [49]. As is illustrated in this case, the ANN optimized by the intelligence algorithm after denoising provides a better power load forecasting effect.Experiment 4.
The power load on Thursday, April 24, 2008, and Tuesday, April 29, 2008, from New South Wales of Australia is used to globally testing the proposed hybrid model. The results are shown as in Tables6 and 7. To forecast the power load on Thursday, April 24, 2008, the historical values from the Thursdays of the first three weeks, April 3, April 10, and April 17, are chosen, respectively. To forecast the power load on Tuesday, April 29, 2008, the historical values from the Thursdays of the first three weeks, April 8, April 15, and April 22, are chosen, respectively. Test results (AE, MAE, RMSE, and MAPE) are presented in Tables 6 and 7 and part (a) of Figure 11. Table 6 indicates that EMD-CSAWNN has the highest accuracy forecasting results on Thursday, April 24, 2008; the maximum, minimum, and average MAPE values are 1.2355% at 2:00, 1.1836 at 14:00, and 1.2025%, respectively. The second-highest to sixth-highest accurate models are GABPNN, BPNN, RBFNN, CSAWNN, WNN, ENN, and ARIMA with average MAPE values of 1.6436%, 1.7795%, 1.8314%, 2.3373%, 2.4141%, 3.1381%, and 5.5872%, respectively. Table 7 indicates that EMD-CSAWNN still yields the highest accuracy forecasting value from among all of other models mentioned in this paper when forecasting power load on Tuesday, April 29, 2008; the maximum, minimum, and average MAPE values are 1.9844% at 14:00, 1.6033% at 12:00, and 1.7811%, respectively. According to the average MAPE value, CSANN is the second most accurate model, GABPNN is the third most accurate model, RBFNN is the fourth most accurate model, WNN is the fifth most accurate model, BPNN is the sixth most accurate model, ENN is the fifth most accurate model, and ARIMA is the sixth most accurate model with average MAPE values of 2.5827%, 2.8472%, 2.8473%, 3.1717%, 3.2725%, 3.5151%, and 3.7549%, respectively. As shown in Table 6, the average MAPE afforded by the hybrid model decreased by 78.48% compared with the maximum average MAPE value. In Table 7, the average MAPE afforded by the hybrid model decreased by 52.57% compared with the maximum average MAPE value. In addition, it is shown that the value of MAPE offered by the hybrid model is more stable than that of the other proposed models. By comparing GABPNN, BPNN, CSAWNN, WNN, ARIMA, ENN, and RBFNN, it is shown that the hybrid EMD-CSAWNN model is better than that of single models.Table 6
Comparison of power load forecasting result by using different methods on Thursday, Apr. 24, 2008, from New South Wales.
Time
0:00
2:00
4:00
6:00
8:00
10:00
12:00
14:00
16:00
18:00
20:00
22:00
Average
BPNN
AE
−83.6894
−49.7095
42.2092
70.1483
−1.1418
81.1176
71.0641
−6.3979
−83.6375
70.2015
−74.2667
68.5518
−10.0773
MAE
139.0489
107.2347
115.6315
153.4157
181.8889
201.3444
189.9489
186.2829
187.2521
163.6832
160.6981
127.2364
161.7255
MSE
20126.81
11597.8
13753.41
23870.75
33776.23
41545.23
37131.99
34961.75
35301.52
27696.51
25909.11
16227.94
27452.27
MAPE
1.7927
1.5953
1.6655
1.7118
1.8221
1.9574
1.8205
1.7860
1.8178
1.7268
1.7671
1.5004
1.7795
GABPNN
AE
−70.5662
−61.4181
14.0903
−5.1011
70.2694
−93.49
62.9971
−5.0262
−175.466
9.4920
87.4908
14.0313
3.7326
MAE
133.1993
118.7066
120.9694
156.3609
168.3905
175.0163
159.27
175.0882
175.4656
142.1779
169.1021
135.8504
149.0696
MSE
17861.21
14198.39
14907.39
24782.75
28896.06
30696.38
25765.8
30889.71
31099.72
20391.04
28608.41
18670.9
22988.39
MAPE
1.7178
1.7651
1.7397
1.7443
1.687962
1.7025
1.5259
1.6792
1.6996
1.4971
1.8594
1.6012
1.6436
WNN
AE
−94.9552
11.0611
86.5661
109.4117
129.0671
−120.256
−133.083
−6.5916
−114.111
−121.37
114.9819
−97.3915
1.46853
MAE
189.6161
163.6107
182.5346
226.7625
249.978
261.0484
274.6263
247.5517
241.1515
223.7553
230.5013
215.3979
219.2263
MSE
36049.23
26940.16
33456.37
51459.93
62579.29
68421.2
75544.98
61478.9
58389.42
50583.32
53412.08
46583
49476.49
MAPE
2.4456
2.4343
2.6249
2.5296
2.5013
2.5395
2.6342
2.3707
2.3368
2.3564
2.5336
2.5407
2.4141
CSAWNN
AE
−21.5793
74.9257
−92.6562
−121.433
−242.688
−4.1947
112.5843
−99.8578
108.9015
−117.485
9.2803
18.8441
0.6051
MAE
187.2812
162.3261
159.9964
230.9195
242.688
234.3061
254.4321
244.3948
207.1397
244.4175
200.3291
213.4906
212.3467
MSE
35548.74
26522.19
25927.18
53405.7
60342
55413.93
65290.47
60673.03
43149.33
60030.08
40539.89
45950.18
46910.12
MAPE
2.4171
2.4145
2.3015
2.5755
2.4308
2.2807
2.4378
2.3433
2.0082
2.5742
2.2039
2.5175
2.3373
ARIMA
AE
−215.208
165.0423
7.8473
8.6388
−561.165
316.2328
−273.573
24.7401
313.4157
212.8736
210.3343
−0.2804
11.1686
MAE
482.1711
333.4499
400.0869
495.5127
561.1652
628.2508
580.8741
607.6131
624.4892
463.577
453.6873
432.3235
507.6683
MSE
234075.7
112032
162067.4
257238.9
322221.9
396897.3
341886.4
369861.6
394541.8
216189.5
206839.9
190582.5
268746.7
MAPE
6.2181
4.9614
5.7564
5.5391
5.6139
6.1119
5.5704
5.8302
6.0641
4.8836
4.9889
5.0989
5.5872
ENN
AE
−134.3503
−116.9724
26.1762
−9.5096
134.8336
−177.9787
121.4220
−9.2200
−334.7105
17.6777
166.4174
26.0920
7.1333
MAE
254.0135
226.2179
230.6157
298.0677
321.2936
333.8439
304.7279
334.1142
334.7105
272.1785
321.8414
259.5614
284.6262
MSE
64935.4893
51544.7643
54128.8887
89996.8881
105084.4741
111677.8252
94240.0696
112436.8781
113110.3356
74693.6074
103626.7463
68117.7863
83742.8823
MAPE
3.2760
3.3639
3.3167
3.3253
3.2205
3.2476
2.9197
3.2042
3.2422
2.8661
3.5391
3.0595
3.1381
RBFNN
AE
−73.9591
−57.3006
1.4359
−115.7463
−11.4324
−103.9494
−69.3536
80.0282
145.0666
−90.0793
−80.4687
−4.8787
−31.7198
MAE
122.8559
107.0627
125.1804
199.8651
161.9459
204.8784
190.0733
161.0021
145.0666
178.5049
152.1063
180.7679
160.7758
MSE
15669.7874
11588.6074
15777.1354
40494.5143
26595.7994
43189.1618
37331.3576
26588.7852
21332.3809
31894.7300
23497.2824
33025.1939
27248.7280
MAPE
1.7706
1.5996
1.4716
2.0627
1.5745
1.9708
1.8099
1.6528
1.5475
2.0890
2.0250
2.4029
1.8314
EMD-CSAWNN
AE
−5.8199
−40.1789
3.2302
−53.9029
0.3157
56.8411
−1.4835
−2.3791
−58.6926
−55.7161
−59.25
55.1392
−9.7616
MAE
92.57402
83.04634
83.72504
106.387
120.3944
124.3859
128.0547
123.5589
124.0705
116.6259
111.4998
102.3642
109.1557
MSE
8608.774
6911.941
7025.682
11341.53
14539.97
15513.36
16423.82
15278.06
15453.49
13635.16
12451.61
10503
12183.36
MAPE
1.1938
1.2355
1.2045
1.1866
1.2037
1.2096
1.2272
1.1836
1.2016
1.2293
1.2265
1.2071
1.2025Table 7
Comparison of power load forecasting result by using different methods on Tuesday, Apr. 29, 2008, from New South Wales.
Time
0:00
2:00
4:00
6:00
8:00
10:00
12:00
14:00
16:00
18:00
20:00
22:00
Average
BPNN
AE
4.3235
104.427
1.9326
−169.582
7.9721
−163.878
−189.155
−9.9777
−133.868
−23.2861
22.6951
−8.1229
−9.6227
MAE
217.8733
215.1561
232.6085
296.2496
297.4614
329.2282
349.862
321.6289
310.8799
328.0564
307.5204
253.7871
293.1005
MSE
47552.31
46537.35
54319.34
88796.61
89366.89
109897.5
122842.7
104957.9
98292.89
108873.5
96080.46
64784.3
88590.6
MAPE
2.8131
3.2020
3.3602
3.3071
3.0174
3.2299
3.4431
3.2158
3.1271
3.4705
3.4071
3.0137
3.2725
GABPNN
AE
−128.819
116.37
−93.8387
5.9302
−25.0426
−157.356
2.0908
145.7456
−282.984
11.7947
129.7324
−11.5832
−16.1558
MAE
260.8535
181.5997
191.8484
262.899
271.7663
306.0667
293.922
285.5415
282.9845
279.7753
244.2308
206.0163
254.9086
MSE
68557.96
35322.65
39010.9
71358.5
76860.8
96621.05
91438.31
82785.86
81937.17
80203.93
60197.01
43169.73
68551.07
MAPE
3.3665
2.7102
2.7624
2.9323
2.7553
3.0142
2.8718
2.8589
2.8581
2.9660
2.7122
2.4488
2.8472
WNN
AE
−102.792
−89.2578
95.7712
−95.0562
15.7301
−8.0755
−132.3
110.6793
127.7917
−302.974
17.0601
33.8715
23.1926
MAE
237.6748
198.5173
194.0934
276.8971
362.0218
334.4987
301.0269
303.3078
290.4403
302.9738
315.3741
285.342
284.0601
MSE
57399.3
40265.36
39263.54
79362.39
131879.1
114268.9
91690.29
94464.47
85736.09
95667.19
101536.4
83146.33
85230.43
MAPE
3.0684
2.9532
2.8117
3.0867
3.6747
3.2822
2.9631
3.0161
2.9222
3.2051
3.4881
3.3871
3.1717
CSAWNN
AE
2.7198
2.9868
−88.3877
122.0205
−143.045
−12.6708
−13.6472
4.8180
110.5905
−114.276
1.9674
−127.205
0.9651
MAE
205.4517
170.8708
173.308
216.7274
244.0602
287.2786
252.2365
266.4932
268.5673
259.4214
244.6435
229.1077
231.6409
MSE
42279.49
29713.52
30554.35
47370.17
60689.6
82774.86
63979.43
72549.97
72983.45
67791.63
60267.6
52777.79
55535.78
MAPE
2.6528
2.5465
2.4986
2.4191
2.4748
2.8256
2.4824
2.6584
2.6966
2.7463
2.7141
2.7266
2.5827
ARIMA
AE
22.0044
278.4617
−1.9584
−222.701
−204.82
−201.995
−411.791
−53.4355
193.9193
−232.404
−15.6347
3.8278
−48.5181
MAE
276.2362
278.4617
255.0419
411.837
374.2884
384.3065
411.7913
355.8697
329.2014
404.3898
374.5698
299.1434
336.7285
MSE
78244.45
78470.78
65671.56
170074.2
144274.5
151190.6
173159.5
133331.2
112870.9
165704.4
141249.1
90175.37
118788.7
MAPE
3.5681
4.1448
3.6876
4.5914
3.8021
3.7762
4.0691
3.5361
3.2986
4.2747
4.1591
3.5501
3.7549
ENN
AE
−122.3804
−106.1050
112.1774
−124.5511
11.2302
−5.7063
−154.8530
138.7430
150.4820
−334.3572
13.6915
25.3581
23.7370
MAE
266.7030
226.0328
225.6456
309.9669
380.8511
365.7942
342.4249
342.0588
332.0910
334.3572
337.4444
308.6425
314.8147
MSE
71574.4601
51519.8987
51615.5316
97406.4037
145424.6816
135137.5738
117847.8733
118529.0844
111042.3847
113745.3418
115039.8651
96156.4770
102413.9259
MAPE
3.4432
3.3631
3.2647
3.4560
3.8651
3.5920
3.3700
3.4068
3.3416
3.5383
3.7352
3.6650
3.5151
RBFNN
AE
−128.8194
116.3700
−93.8387
5.9302
−25.0426
−157.3556
2.0909
145.7456
−282.9845
11.7948
129.7324
−11.5832
−16.1558
MAE
260.8535
181.5997
191.8484
262.8990
271.7663
306.0667
293.9220
285.5415
282.9845
279.7753
244.2308
206.0163
254.9086
MSE
68557.9649
35322.6478
39010.8992
71358.4965
76860.7961
96621.0544
91438.3143
82785.8632
81937.1664
80203.9338
60197.0099
43169.7279
68551.0659
MAPE
3.3665
2.7103
2.7625
2.9324
2.7553
3.0142
2.8718
2.8589
2.8582
2.9660
2.7123
2.4489
2.8473
EMD-CSAWNN
AE
71.8585
19.0263
119.4944
−0.2347
−85.8361
21.7152
162.9778
198.7374
12.5729
0.7221
7.0831
11.6529
13.8229
MAE
135.0333
124.2513
119.4944
159.9286
166.286
188.3798
162.9778
198.7374
191.4925
164.075
155.0245
145.8072
159.2127
MSE
18544.32
15861.56
14935.85
25684.91
28077.06
36387.53
26760.25
39816.23
38064.54
28373.25
24954.2
21683.46
26609.38
MAPE
1.7439
1.8477
1.7241
1.7831
1.6889
1.8577
1.6033
1.9844
1.9379
1.7341
1.7161
1.7319
1.7811Figure 11
The comprehensive evaluation of forecasting models in Experiments4 and 5. (a) Comparison of forecasting results by eight models on Thursday, April 24, 2008, and Tuesday, April 29, 2008, from New South Wales of Australia. (b) Comparison of forecasting results by eight models on Saturday, June 28, 2008, and Monday, June 30, 2008, from Victoria of Australia.
(a)
(b)Experiment 5.
The power load on Saturday, June 28, 2008, and Monday, June 30, 2008, from Victoria of Australia is used to further prove that the proposed hybrid model can improve the performance of power load forecasting in different cases. The historical values from the Saturday of the first three weeks, June 7, June 14, and June 21, are chosen in order to forecast the power load on Saturday, June 28, 2008. Moreover, the historical values from the Monday of the first three weeks, June 9, June 16, and June 23, are chosen in order to forecast the power load on Monday, June 30, 2008. The experimental results are presented in Tables8 and 9 and part (b) of Figure 11. For Saturday (June 28, 2008) data, the average MAPE values for the BPNN, GABPNN, WNN, CSAWNN, ARIMA, ENN, RBFNN, and EMD-CSAWNN models are 3.3224%, 2.8656%, 3.2301%, 2.5248%, 3.7716%, 3.3890%, 3.1754%, and 1.7699%, respectively. The maximum MAPE values for the six models are 3.5467% at 16:00, 3.3061% at 18:00, 3.5751% at 4:00, 2.9527% at 16:00, 4.2174% at 0:00, 3.7591% at 20:00, 3.3889% at 10:00, and 2.0503% at 14:00, respectively, and the minimum MAPE values are 2.9644% at 20:00, 2.2710% at 12:00, 2.7591% at 20:00, 2.1818% at 20:00, 3.4718% at 4:00, 3.1379% at 18:00, 2.8405% at 22:00, and 1.5798% at 10:00, respectively. The differences between the maximum and minimum MAPE values for the models are 0.5823%, 1.0351%, 0.816%, 0.1936%, 0.7456%, 0.6212%, 0.5484%, and 0.4705%. For Monday (June 30, 2008) data, the average MAPE values for the BPNN, GABPNN, WNN, CSAWNN, ARIMA, ENN, RBFNN, and EMD-CSAWNN model are 3.2747%, 2.8415%, 3.0806%, 2.5422%, 3.7656%, 3.4113%, 3.1982%, and 1.7312%, respectively. The maximum MAPE values for the six models are 3.9915% at 22:00, 3.4117% at 2:00, 3.6561% at 6:00, 2.8144% at 0:00, 4.1783% at 22:00, 3.6738% at 8:00, 3.5262% at 0:00, and 2.0361% at 12:00, respectively, and the minimum MAPE values are 3.0113% at 20:00, 2.3376% at 22:00, 2.4491% at 16:00, 2.1398% at 20:00, 3.2795% at 0:00, 3.2235% at 4:00, 2.9466% at 22:00, and 1.4521% at 10:00, respectively. The differences between the maximum and minimum MAPE values for the models are 0.9802%, 1.0741%, 1.207%, 0.6746%, 0.8988%, 0.4503%, 0.5796%, and 0.584%. Therefore, the proposed hybrid EMD-CSAWNN model is not only the most accurate but also the most stable of the investigated forecasting models.Table 8
Comparison of power load forecasting result by using different methods on Saturday, Jun. 28, 2008, from Victoria.
Time
0:00
2:00
4:00
6:00
8:00
10:00
12:00
14:00
16:00
18:00
20:00
22:00
Average
BPNN
AE
−74.2069
88.5026
164.3428
−94.2384
−108.445
99.3887
−131.822
−97.1293
132.6263
−9.9181
−199.556
95.9906
−25.4772
MAE
161.8614
176.3894
164.3428
175.3241
219.8255
225.0625
233.3815
229.9635
234.7046
233.2327
199.5555
182.4864
206.2306
MSE
26380.04
31503.46
27247.69
31227.59
48517.87
51693.03
54905.29
54093.61
55631.95
54814.4
39922.9
33457.23
43708.18
MAPE
2.9997
3.3342
3.4049
3.1834
3.3696
3.3340
3.5174
3.4406
3.5467
3.1652
2.9644
3.0498
3.3224
GABPNN
AE
95.3661
152.9973
70.1585
95.8643
21.4256
−10.3011
65.8767
6.8887
−16.0091
−14.2224
97.4645
19.2111
3.9264
MAE
163.6398
152.9973
142.3091
180.6853
192.6446
181.9683
151.09
198.873
168.6475
242.7444
202.9146
183.0895
178.1422
MSE
27233.94
23592.27
20427.45
33187.61
37879.35
34166.35
23242.55
40431.09
28919.45
60044.79
41793.54
34130.68
33502.97
MAPE
3.0324
2.8936
2.9494
3.2736
2.9262
2.7029
2.2710
3.0080
2.5461
3.3061
3.0207
3.0576
2.8656
WNN
AE
175.8896
27.8588
−7.9285
86.3077
139.9775
201.231
90.7231
39.9104
18.8402
−108.905
−114.363
193.1549
4.4674
MAE
175.8896
156.0152
172.6755
164.069
207.6373
201.231
193.4483
205.2855
183.8711
259.0394
186.1753
193.1549
200.6168
MSE
31845.68
25360.07
30586.03
27917.8
45875.7
42312.96
38650.67
43776.7
34700.63
69347.15
35949.54
38952.51
42677.08
MAPE
3.2604
2.9507
3.5751
2.9766
3.1492
2.9988
2.9427
3.1425
2.7626
3.5188
2.7591
3.2385
3.2301
CSAWNN
AE
−44.652
−131.152
−63.0165
11.4112
−9.5339
−0.7664
171.5844
193.7814
7.9337
106.7022
3.6645
−69.7819
12.6984
MAE
124.0692
131.1523
114.0833
146.3165
167.5085
165.3564
171.5844
193.7814
195.2323
191.7276
146.8918
141.5059
156.8164
MSE
15804.77
17589.3
13292.63
21642.06
29014.79
27346.81
30261.31
37723.03
38246.53
36979.81
21623.75
20248.95
25425.16
MAPE
2.2995
2.4781
2.3604
2.6579
2.5319
2.4803
2.5786
2.9203
2.9527
2.6113
2.1818
2.3631
2.5248
ARIMA
AE
18.5613
124.0827
−14.7831
111.8472
−20.6124
−19.7243
120.3048
26.1011
−146.313
−172.198
23.3302
145.1471
18.9775
MAE
227.3289
204.722
167.4229
212.6136
230.4365
247.6076
230.858
253.9037
267.5722
296.6621
255.3748
241.6215
233.6259
MSE
52394.21
43176.64
28355.68
45620.29
53698.66
62570.82
54191.74
65290.6
72666.23
90860.47
67238.24
59345.33
56855.39
MAPE
4.2174
3.8741
3.4718
3.8504
3.5598
3.7094
3.4799
3.8433
4.0557
4.0217
3.7627
4.0433
3.7716
ENN
AE
−98.995
−95.8081
17.1172
6.0116
90.0902
−141.101
100.5493
4.7625
−229.624
18.1665
126.2839
11.7423
6.5262
MAE
189.6312
190.0733
171.6852
196.7854
231.2569
235.0732
210.1064
229.5182
229.6243
230.8097
254.2102
198.2688
210.274
MSE
36109.8
36258.98
29924.18
39256.18
56148.35
56021.74
44222.47
53380.62
52793.31
53676.67
64753.52
39571.84
45444.62
MAPE
3.5163
3.5974
3.5538
3.5618
3.4650
3.4901
3.1873
3.4500
3.4851
3.1379
3.7591
3.3164
3.3890
RBFNN
AE
−98.2027
−73.2415
57.3594
77.1316
9.4619
96.7952
101.7499
−14.9842
−97.2306
97.9996
−110.216
86.7962
−11.4032
MAE
172.1148
156.3662
146.8716
171.0673
214.0821
228.8974
212.5963
212.2372
213.3772
228.0278
213.4688
169.8724
196.9359
MSE
30159.71
24592.22
21954
29579.24
47567.12
53850.57
45459.14
45699.94
45988.64
52566.64
45607.02
28877.88
39869.92
MAPE
3.1912
2.9543
3.0386
3.0942
3.2266
3.3889
3.2246
3.1832
3.2214
3.1122
3.1605
2.8405
3.1754
EMD-CSAWNN
AE
35.0621
51.5938
41.9313
−56.3883
−10.8655
51.1997
−42.3959
84.1310
70.1839
82.8233
46.4874
−55.1894
6.6501
MAE
91.3676
100.3773
81.2863
92.1573
108.8744
104.7747
119.9606
137.012
110.6269
141.76
109.7464
115.5855
109.712
MSE
8629.51
10214.74
6858.997
8885.95
12543.16
11063.31
14958.06
19196.88
12798.02
20368.17
12301.46
13873.76
12612.17
MAPE
1.6951
1.8968
1.6898
1.6576
1.6298
1.5798
1.8118
2.0503
1.6639
1.9325
1.6329
1.9192
1.7699Table 9
Comparison of power load forecasting result by using different methods on Monday, Jun. 30, 2008, from Victoria.
Time
0:00
2:00
4:00
6:00
8:00
10:00
12:00
14:00
16:00
18:00
20:00
22:00
Average
BPNN
AE
−97.1358
−5.79918
−77.632
−69.5865
−5.12235
181.7938
−108.862
172.1967
−11.7116
−112.91
−1.05112
−5.29961
7.270652
MAE
180.8076
174.7874
160.6262
161.5628
166.0493
181.7938
187.9272
172.1967
176.0733
215.9977
187.979
197.6079
181.9883
MSE
32946.72
30612.98
26035.71
26290.49
27633.14
33390.49
35789.18
29752.65
31224.95
46849.69
35447.84
39567.68
33681.6
MAPE
3.2431
3.2773
3.3784
3.3739
3.1979
3.2166
3.3433
3.0728
3.0901
3.2467
3.0113
3.9915
3.2747
GABPNN
AE
87.0447
182.023
62.5079
−12.5589
−92.6525
−84.4608
−25.4557
19.48
−22.7953
76.9831
78.4895
−6.6962
4.1017
MAE
148.3886
182.023
116.5217
112.0699
155.8314
181.3435
153.1462
162.2208
182.3206
182.7003
175.2353
132.1649
158.2286
MSE
23338.89
33529.04
14089.3
12900.14
24927.85
33303.2
24526.58
26867.37
34101.29
35242.19
31203.34
17935.79
26253.68
MAPE
2.6633
3.4117
2.4520
2.3381
3.005
3.2113
2.7189
2.8955
3.1883
2.7492
2.8096
2.3376
2.8415
WNN
AE
−3.4921
−74.9511
−134.269
−80.2326
−142.363
64.7736
−111.042
−89.0593
139.566
14.1565
−82.5293
66.4498
−4.1007
MAE
153.3954
185.6073
134.2689
174.9947
142.3635
174.6803
186.5582
193.4015
139.566
212.2396
156.4929
167.1148
171.6629
MSE
23561.14
35229.8
18364.93
31218.27
21459.47
33650.21
37245.82
38184.73
19503.57
47310.2
24810.93
28778.48
31298.91
MAPE
2.7532
3.4749
2.8191
3.6561
2.7458
3.0828
3.3176
3.4628
2.4491
3.1823
2.5057
2.9566
3.0806
CSAWNN
AE
−5.0573
−3.8213
1.3263
41.3898
54.9301
−66.0719
14.6333
−3.6448
149.7019
80.7116
3.7387
−73.6169
8.0382
MAE
156.8892
136.6022
125.5879
113.7244
116.9543
143.8786
140.1588
139.9815
149.7019
159.8907
133.4934
137.4735
141.3812
MSE
24786.49
18999.83
15803.74
13261.64
13725.91
20966.75
20021.99
19790.28
22684.1
25822.94
17849.43
18934.4
20481.5
MAPE
2.8144
2.5621
2.6401
2.3753
2.2519
2.5478
2.5001
2.5018
2.6261
2.4020
2.1398
2.4302
2.5422
ARIMA
AE
14.7023
24.8103
26.0758
96.5552
−103.598
67.3243
1.5328
−12.8999
229.6626
254.5464
−133.048
125.5915
−0.7485
MAE
183.057
217.205
158.1086
187.75
186.6591
198.3604
205.3784
213.5527
229.6626
254.5464
240.3672
236.3162
209.7854
MSE
34965.65
47841.28
26107.02
35347.8
35329.16
41301.58
42893.69
46719.97
54654.13
65818.94
58199.06
55984.16
45718.98
MAPE
3.2795
4.0777
3.3209
3.9201
3.5982
3.5138
3.6743
3.8162
4.0242
3.8204
3.8526
4.1783
3.7656
ENN
AE
−88.6146
−85.0312
75.6176
−69.7223
4.5320
−0.8129
−87.546
82.3113
88.4143
−227.985
3.5052
12.0972
8.3585
MAE
186.9891
175.8838
153.515
161.1809
190.7898
196.0601
185.4341
186.5996
187.2131
227.9848
223.304
199.3213
189.5443
MSE
35060.02
31090.34
23853.32
26154.44
36507.3
38584.26
34554.06
35020.38
35263.23
52336.47
50161.92
39941.99
36505.16
MAPE
3.3574
3.2973
3.2235
3.3669
3.6738
3.4690
3.3024
3.3301
3.2811
3.4287
3.5763
3.5237
3.4113
RBFNN
AE
−98.586
97.2085
−72.7937
0.0161
−10.8099
−94.714
−1.0184
92.2613
−182.302
5.0740
100.423
−5.2195
−12.4652
MAE
196.4346
166.1881
149.5233
155.7541
163.1654
186.8753
180.256
179.795
182.3024
217.8852
194.234
166.6086
177.8178
MSE
38661.36
28328.97
22697.24
24518.37
27001.22
35355.86
32968.16
32671.75
33413.99
47850.86
37890.5
27884.44
32243.32
MAPE
3.5262
3.1117
3.1447
3.2520
3.1402
3.3037
3.2137
3.2056
3.2051
3.2732
3.1130
2.9466
3.1982
EMD-CSAWNN
AE
−4.1652
−40.9803
−54.073
−28.5539
−49.4364
−46.9202
−1.3708
−51.2203
−61.0243
58.5634
44.2871
55.8282
6.1509
MAE
99.0617
81.542
83.7323
70.3801
80.9749
82.0208
114.0354
111.8617
115.3834
129.6408
110.32
105.1976
96.6039
MSE
10306.84
6696.25
7299.67
5046.59
6865.14
6815.57
13126.18
12711.91
13468.79
17077.42
12871.21
11278.86
9834.31
MAPE
1.7771
1.5295
1.7617
1.4701
1.5594
1.4521
2.0361
2.0032
2.0227
1.9479
1.7703
1.8607
1.7312


## 5. Discussion

In this section, we discuss two important evaluation metrics, convergence speed and degree of certainty [50], offered by the GABPNN and CSAWNN models to determine a more practical forecasting model by considering reality factors such as forecasting stability and calculation time. The results illustrate that the CSAWNN model is more practical than the GABPNN model in forecasting power load. In addition, we propose forecasting availability to analyze and evaluate the quality of power load forecasting.


### 5.1. Convergence Speed

The computational complexity of evolutionary algorithms and swarm intelligence still remains a challenging issue; here, we use convergence speed as one of the evaluation metrics to examine the forecasting performance of GABPNN and CSAWNN. We obtain the computation time of the best fitness by analyzing the convergence speed of GA and CSA for use in comparative evaluation of optimization algorithm performance. However, the exploration and development are always two competing goals, and the conflict would exist between the convergence speed and forecasting accuracy. We define performance less than 10−5 as the convergence criteria.We take the data from January 12, January 19, and January 26 as an example to illustrate the convergence speed of GA and CSA; Figure12 shows the results of the comparison of evolutions among GA and CSA with different population sizes. We observed that the fitness values monotonically decrease as the iterations increase. In addition, when the iterations are less than 100, the larger the population size, the faster the convergence speed. We also observed that CSA has better convergence speed than GA. At iteration 20, the convergence of CSAWNN obtained the best speed in a population of 50, and the convergence of GABPNN at the iteration between 60 and 80, at the iteration 60, obtained the best speed.Figure 12
A convergence speed comparison of the fitness values between GA and CSA with different population sizes.


### 5.2. Degree of Certainty

The forecasting results of optimization algorithm-NN are also usually different for each experiment because of the probability mechanism of the optimization algorithm. However, in the actual forecasting field, the future values are not known; thus, we cannot obtain which experiment will obtain the best result. Hence, we use these evaluation metrics to determine the randomness.We defined the degree of certainty as(13)
D
C
σ
=
∑
i
=
1
m
σ
i
-
σ
-
2
m
,where m is the number of experiments, σ
i is the value of ith forecasting experiment on σ, and σ
- is the average value of all n experiments. It is clear that a smaller DC can bring a higher degree of certainty.The scatterplot in Figure13 indicates the MAPE and GRA distributions of different results of 100 experiments for February 2 by, respectively, using GABPNN and CSAWNN. Although the minimum value of MAPE of GABPNN is smaller than that of CSAWNN and D
C
σ of CSAWNN, both are smaller than those of GABPNN. DC (MAPE) and DC (GRA) of CSAWNN are 0.0014 and 0.00016, and the DC (MAPE) and DC (GRA) of GABPNN are 0.0016 and 0.00019. Thus, CSAWNN is a better forecasting method than GABPNN in the actual forecasting field.Figure 13
The DC (MAPE) and DC (GRA) of 100 experiments by comparison of GABPNN and CSAWNN.


### 5.3. Forecasting Availability

Forecasting availability can be measured not only by the square sum of forecasting error but also by the mean and mean squared deviation of the forecasting accuracy. In certain practical circumstances, the skewness and kurtosis of the distribution of forecasting accuracy need further consideration; on that basis, this section will give a general discrete form of forecasting availability [51].Definition 10.
Let(14)
e
i
t
=
-
1
,
x
t
-
x
i
t
x
t
>
1
x
t
-
x
i
t
x
t
,
-
1
≤
x
t
-
x
i
t
x
t
1
,
x
t
-
x
i
t
x
t
>
1
.
>
1
e
i
t denotes the relative forecasting error of ith forecasting method at time t, i
=
1,2
,
…
,
m, t
=
1,2
,
…
,
N. The matrix E
=
e
i
t
m
×
N is called the relative error matrix of the forecasting model.Definition 11.
A
i
t
=
1
-
e
i
t is called the forecasting accuracy of ith forecasting method at time t, i
=
1,2
,
…
,
m, and t
=
1,2
,
…
,
N.Definition 12.
m
i
k
=
∑
t
=
1
N
Q
t
A
i
t
k is called kth-order forecasting availability unit of ith forecasting method, k is a positive integer, i
=
1,2
,
…
,
m, Q
t
,
t
=
1,2
,
…
,
N is the discrete probability distribution of m types of methods at time t, and ∑
t
=
1
N
Q
t
=
1
,
Q
t
>
0.
Especially if the priori information of the discrete probability distribution ofm types of methods is unknown, we define Q
t
=
1
/
N
,
t
=
1,2
,
…
,
N.Definition 13.
m
i
k is called the k-order forecasting availability unit of ith forecasting method, and H is a continuous function of a certain k unit. H
m
i
1
,
m
i
2
,
…
,
m
i
k is called the k-order forecasting availability of ith forecasting method.Definition 14.
WhenH
x
=
x is a continuous function of one variable, H
m
i
1
=
m
i
1 is the 1-order forecasting availability of ith forecasting method. When H
x
,
y
=
x
1
-
y
-
x
2 is a continuous function of two variables, H
m
i
1
,
m
i
2
=
m
i
1
1
-
m
i
2
-
m
i
1
2 is the 2-order forecasting availability of ith forecasting method.
Especially if the first decimal ofH
m
i
1
,
m
i
2
,
…
,
m
i
k is the same, we define H
m
i
1
,
m
i
2
,
…
,
m
i
k
=
10
H
m
i
1
,
m
i
2
,
…
,
m
i
k, where · denotes the fractional part.Definition14 illustrates that the 1st-order forecasting availability is the expectation forecasting accuracy sequence. The 2nd-order forecasting availability is the difference between the expectation and standard deviation of the forecasting accuracy sequence. We use the forecasting availability to evaluate the power load forecasting results in this paper. Through Figure 14, we obtain that the 1st-order and 2nd-order forecasting availability offered by the hybrid model are 0.9222 and 0.8366, respectively, which outperform those of the others; this evaluation result corresponds to the previous evaluation criterion. Thus, the hybrid model is a more valid model than the others.Figure 14
The forecasting availability of BPNN, GABPNN, WNN, CSAWNN, EMD-CSAWNN, and ARIMA.


## 5.1. Convergence Speed

The computational complexity of evolutionary algorithms and swarm intelligence still remains a challenging issue; here, we use convergence speed as one of the evaluation metrics to examine the forecasting performance of GABPNN and CSAWNN. We obtain the computation time of the best fitness by analyzing the convergence speed of GA and CSA for use in comparative evaluation of optimization algorithm performance. However, the exploration and development are always two competing goals, and the conflict would exist between the convergence speed and forecasting accuracy. We define performance less than 10−5 as the convergence criteria.We take the data from January 12, January 19, and January 26 as an example to illustrate the convergence speed of GA and CSA; Figure12 shows the results of the comparison of evolutions among GA and CSA with different population sizes. We observed that the fitness values monotonically decrease as the iterations increase. In addition, when the iterations are less than 100, the larger the population size, the faster the convergence speed. We also observed that CSA has better convergence speed than GA. At iteration 20, the convergence of CSAWNN obtained the best speed in a population of 50, and the convergence of GABPNN at the iteration between 60 and 80, at the iteration 60, obtained the best speed.Figure 12
A convergence speed comparison of the fitness values between GA and CSA with different population sizes.


## 5.2. Degree of Certainty

The forecasting results of optimization algorithm-NN are also usually different for each experiment because of the probability mechanism of the optimization algorithm. However, in the actual forecasting field, the future values are not known; thus, we cannot obtain which experiment will obtain the best result. Hence, we use these evaluation metrics to determine the randomness.We defined the degree of certainty as(13)
D
C
σ
=
∑
i
=
1
m
σ
i
-
σ
-
2
m
,where m is the number of experiments, σ
i is the value of ith forecasting experiment on σ, and σ
- is the average value of all n experiments. It is clear that a smaller DC can bring a higher degree of certainty.The scatterplot in Figure13 indicates the MAPE and GRA distributions of different results of 100 experiments for February 2 by, respectively, using GABPNN and CSAWNN. Although the minimum value of MAPE of GABPNN is smaller than that of CSAWNN and D
C
σ of CSAWNN, both are smaller than those of GABPNN. DC (MAPE) and DC (GRA) of CSAWNN are 0.0014 and 0.00016, and the DC (MAPE) and DC (GRA) of GABPNN are 0.0016 and 0.00019. Thus, CSAWNN is a better forecasting method than GABPNN in the actual forecasting field.Figure 13
The DC (MAPE) and DC (GRA) of 100 experiments by comparison of GABPNN and CSAWNN.


## 5.3. Forecasting Availability

Forecasting availability can be measured not only by the square sum of forecasting error but also by the mean and mean squared deviation of the forecasting accuracy. In certain practical circumstances, the skewness and kurtosis of the distribution of forecasting accuracy need further consideration; on that basis, this section will give a general discrete form of forecasting availability [51].Definition 10.
Let(14)
e
i
t
=
-
1
,
x
t
-
x
i
t
x
t
>
1
x
t
-
x
i
t
x
t
,
-
1
≤
x
t
-
x
i
t
x
t
1
,
x
t
-
x
i
t
x
t
>
1
.
>
1
e
i
t denotes the relative forecasting error of ith forecasting method at time t, i
=
1,2
,
…
,
m, t
=
1,2
,
…
,
N. The matrix E
=
e
i
t
m
×
N is called the relative error matrix of the forecasting model.Definition 11.
A
i
t
=
1
-
e
i
t is called the forecasting accuracy of ith forecasting method at time t, i
=
1,2
,
…
,
m, and t
=
1,2
,
…
,
N.Definition 12.
m
i
k
=
∑
t
=
1
N
Q
t
A
i
t
k is called kth-order forecasting availability unit of ith forecasting method, k is a positive integer, i
=
1,2
,
…
,
m, Q
t
,
t
=
1,2
,
…
,
N is the discrete probability distribution of m types of methods at time t, and ∑
t
=
1
N
Q
t
=
1
,
Q
t
>
0.
Especially if the priori information of the discrete probability distribution ofm types of methods is unknown, we define Q
t
=
1
/
N
,
t
=
1,2
,
…
,
N.Definition 13.
m
i
k is called the k-order forecasting availability unit of ith forecasting method, and H is a continuous function of a certain k unit. H
m
i
1
,
m
i
2
,
…
,
m
i
k is called the k-order forecasting availability of ith forecasting method.Definition 14.
WhenH
x
=
x is a continuous function of one variable, H
m
i
1
=
m
i
1 is the 1-order forecasting availability of ith forecasting method. When H
x
,
y
=
x
1
-
y
-
x
2 is a continuous function of two variables, H
m
i
1
,
m
i
2
=
m
i
1
1
-
m
i
2
-
m
i
1
2 is the 2-order forecasting availability of ith forecasting method.
Especially if the first decimal ofH
m
i
1
,
m
i
2
,
…
,
m
i
k is the same, we define H
m
i
1
,
m
i
2
,
…
,
m
i
k
=
10
H
m
i
1
,
m
i
2
,
…
,
m
i
k, where · denotes the fractional part.Definition14 illustrates that the 1st-order forecasting availability is the expectation forecasting accuracy sequence. The 2nd-order forecasting availability is the difference between the expectation and standard deviation of the forecasting accuracy sequence. We use the forecasting availability to evaluate the power load forecasting results in this paper. Through Figure 14, we obtain that the 1st-order and 2nd-order forecasting availability offered by the hybrid model are 0.9222 and 0.8366, respectively, which outperform those of the others; this evaluation result corresponds to the previous evaluation criterion. Thus, the hybrid model is a more valid model than the others.Figure 14
The forecasting availability of BPNN, GABPNN, WNN, CSAWNN, EMD-CSAWNN, and ARIMA.


## 6. Conclusion

The one-day-ahead power load forecasting is an extremely important problem in power load planning, secure operation, and energy expenditure economy. Assessment of the power load as accurately and quickly as possible is the primary objective in power load forecasting. However, power load is affected by various uncertain factors such as climate change and the social environment, which may lead to difficulty in obtaining accurate power load forecasts. The accuracy of traditional individual forecasting methods, which lack denoising, is not satisfactory for power load forecasting. Herein, a hybrid EMD-CSAWNN model for short-term power load forecasting is developed. The empirical model decomposition technique is applied to reduce the high-frequency items. On the basis that WNN can handle the data with nonlinear features, the ensemble forecasting method is adopted to overcome the uncertainty of the outcomes that can be attributed to the randomness of the initialization of the single WNN. Moreover, we use the CSA to optimize the parameter in the ensemble forecasting model. Experimental studies of power load forecasting in NSW demonstrate that the hybrid model has higher precision than conventional forecasting models. The proposed EMD-CSAWNN model can provide efficient computation and satisfactory forecasting accuracy for this type of data. Therefore, the developed hybrid approach is suggested for broad application in power load forecasts or even other fields such as wind speed and traffic flow forecasts.



---

*Source: 3205396-2016-12-15.xml*