Forecasting Electricity Demand Using a New Grey Prediction Model with Smoothness Operator

Abstract

A stable electricity supply is the basis for ensuring the healthy and sustained development of a regional economy. Reasonable electricity prediction is the key to guaranteeing the stability and efficiency of electricity supply. To this end, we used a reformative grey prediction model to forecast electricity demand. In order to effectively improve the smoothness of a raw modelling sequence, we employed an existing smoothing algorithm that significantly compressed the amplitude of the random oscillation sequence. Then, an improved grey forecasting model with three parameters (IGFM_TP) was deduced. In the end, a new model was used to forecast the demand for electricity of one city in the western region of China, and comparisons of simulation values and errors with those of GFM_TP, GM(1,1), DGM(1,1) and SAIGM were conducted. The findings show that the mean absolute simulation percentage error of IGFM_TP was 7.8%, and those of the other four models were 12.1%, 12.3%, 11.1%, and 10.1%, respectively. Therefore, the simulation precision of the new model achieved an optimal effect. The proposed new grey model provides is an effective method for electricity demand prediction.

1. Introduction

China is the second largest economy in the world (Figure 1a), and its manufacturing industry is one of the most developed. China consumes vast amounts of electric energy every year (Figure 1b). Hence, the adequate supply of electricity is one of the prerequisites for ensuring the sustainable development of China’s economy; reasonable electricity demand forecasting plays an important role in generation and distribution. Since the early 1970s, a variety of estimation methods have been proposed in studies related to electricity demand forecasting, including semiparametric regression [1], time series modeling [2], exponential smoothing [3], Bayesian statistics [4], time-varying splines, neural networks [5], decomposition techniques [6], transfer functions [7], grey dynamic models [8], and judgmental forecasting [9].

Electricity demand is influenced by various factors. However, a full analysis of all the factors that affect electricity demand is complicated. Electricity demand varies from season to season in the same area, which reflects its random oscillation feature. Under this circumstance, traditional methods for forecasting the electricity demand are insufficient.

In systems theory, a system can be defined by a color that demonstrates the magnitude of clear information about the system. For example, in a black box system, the internal characteristics or mathematical equations that describe the dynamics of the system are fully unknown. Conversely, a system is called a white system when its description is completely known. A system in an intermediate position is known as a grey system. In practice, due to various certainties of systems, every system can be modelled as a grey system. Due to noise both inside and outside a system and limitations of human cognitive abilities, the information we are able to acquire about a system is always limited and incomplete. Grey system theory, developed by Professor Julong Deng [10] based on “Grey Box” thinking in the 1980s, is an important and useful tool employed to study grey uncertain problems with poor information and small samples [11]. The grey forecasting model, as one of the significant constituents of grey system theory, has attracted attention and has been actively studied by many scholars. The GM(1,1) model is the most important part of grey prediction theory; however, many scholars discovered that the prediction precision of this model is unstable, and extensive research mainly on the following topics has been conducted: Converting an original sequence for improvement of smoothness [12]; enhancing the computational methods of parameters [13,14,15,16]; modifying residuals of models [17,18,19]; reforming the modelling fashion, and performing some preparations for expanding GM(1,1) models [20]; investigating the modelling conditions [21]; extending the structure of traditional grey prediction models [22]; optimizing the parameters of grey prediction models, such as initial value [23,24], background value [25,26], and the order of accumulating generation operators [27,28]; and combining the grey prediction model with other modelling methods [29,30].

The above measures and methods have significantly improved the simulation and prediction performance of grey prediction models. However, the structure of these grey prediction models is fixed in either non-homogeneous or homogeneous exponential form. In the real world, the development of the system is always influenced by many factors. Hence, the existing grey prediction models with fixed structures cannot usefully describe the development trends of complicated systems. So, we think the adaptability of the existing grey models is poor. In order to solve these difficulties, a grey forecasting model with three parameters (GFM_TP) is proposed [31]. The GFM_TP model achieves unbiased simulation and prediction for both homogenous and non-homogenous exponential and linear function sequences. The proposed GFM_TP model is compatibility with exponential functions (homogenous and non-homogenous) and linear functions, and its model structure can intelligently match and dynamically optimize through adjusting the parameters of the GFM_TP model according to the data characteristics of an original sequence. Thus, the GFM_TP model has a more flexible structure and better performance than traditional grey prediction models.

Satisfactory simulation or prediction accuracy can usually be obtained when an original sequence for building the GFM_TP model demonstrates a monotonic increase (or decrease), but the prediction or simulation accuracy of the GFM_TP model is regarded as unsatisfactory when random oscillation characteristics are shown in the modelling sequence. To this end, we applied an existing smoothing algorithm [32,33] that can compress the amplitude of a random oscillation sequence to improve the smoothness of a raw modelling sequence. After this, a new grey prediction model, SAIGM, was proposed based on the classical GM(1,1) model. However, SAIGM performs poorly, due to the simple structure of GM(1,1). Then, an improved grey forecasting model with three parameters (IGFM_TP) was deduced based on the existing smoothing algorithm and the GFM_TP model. At last, the new model was used to compare the simulation values and errors with those of GFM_TP, GM(1,1), SAIGM, and DGM(1,1), and to predict the electricity demand of one city in the western region of China. According to the results, the optimal effects of the simulation are provided in the new model. The research results of this paper significantly improve the simulation and prediction performance of grey prediction models.

In summary, the main contribution of this work is that we propose a three-parameter grey prediction model, based on the existing smoothing algorithm and the GFM_TP model, that can be used to more effectively build a grey prediction model with oscillating sequences. Therefore, the research results of this paper expand the application scope and improve the modeling ability of grey forecasting model.

The remainder of the paper is organized as follows. In Section 2, the smoothness operator is introduced. In Section 3, the grey forecasting model with three parameters is built. In Section 4, an improved GFM_TP model based on amplitude contraction is proposed. In Section 5, the IGFM_TP model is used to compare simulation values and errors with other grey models and to predict the electricity demand of one city in the western region of China. Conclusions are presented in Section 6.

2. Basic Concepts

The performance of a grey prediction model is closely related to the smoothness of the modelling sequence. The smoother the modelling sequence, the higher the accuracy of the grey prediction model. In this paper, the smoothness of the modelling sequence was improved using an oscillating sequence smoothing algorithm. Before building the grey prediction model, we first introduced a smoothness algorithm of the random oscillation sequence.

Definition 1. 

Assume that a modelling sequence is given:

$$Y^{\left(0\right)}=\left(y^{\left(0\right)}\left(1\right),y^{\left(0\right)}\left(2\right),\cdots ,y^{\left(0\right)}\left(m\right)\right),$$

Then,

(1) 

If, for k = 2, 3, …, m$y^{\left(0\right)}\left(k\right)-y^{\left(0\right)}\left(k-1\right)>0$, then$Y^{\left(0\right)}$is a monotonic increasing sequence;

(2) 

If, for$k=2,3,\cdots ,m$,$y^{\left(0\right)}\left(k\right)-y^{\left(0\right)}\left(k-1\right)<0$, then$Y^{\left(0\right)}$is a monotonic decreasing sequence;

(3) 

If, for any$k$and$k=2,3,\cdots ,m$,

$$y^{\left(0\right)}\left(k\right)-y^{\left(0\right)}\left(k-1\right)>0\mathrm{and}{y}^{\left(0\right)}\left(k\right)-y^{\left(0\right)}\left(k-1\right)<0,$$

then$Y^{\left(0\right)}$is called a random oscillation sequence [34]. Note, the upper right corner mark (0) is used to represent the original sequence, and (1) is the accumulating generation sequence.

Assume that

$$A=\max \left\{y^{\left(0\right)}\left(k\right)|k=1,2,\cdots ,m\right\},a=\min \left\{y^{\left(0\right)}\left(k\right)|k=1,2,\cdots ,m\right\},$$

then

$$T_{\mathrm{amplitude}}=A-a,$$

is called the amplitude of$Y^{\left(0\right)}$.

Definition 2. 

Assume that

$$Y=\left(y^{\left(0\right)}\left(1\right),y^{\left(0\right)}\left(2\right),\cdots ,y^{\left(0\right)}\left(m\right)\right),$$

is a random oscillation sequence, and another sequence

$$Y^{\left(0\right)}D=\left(y^{\left(0\right)}\left(1\right)d,y^{\left(0\right)}\left(2\right)d,\cdots ,y^{\left(0\right)}\left(m-1\right)d\right),$$

where,

$$y^{\left(0\right)}\left(k\right)d=0.25y^{\left(0\right)}\left(k\right)+0.25y^{\left(0\right)}\left(k+1\right)+0.5T_{\mathrm{amplitude}}$$

(1)

In Equation (1),$k=1,2,\cdots ,m-1$,$T_{\mathrm{amplitude}}$is the amplitude of$Y^{\left(0\right)}$,$D$is a sequence operator and a first-ordersmoothness operator of$Y^{\left(0\right)}$, and$Y^{\left(0\right)}D$is a smoothness sequence of$Y^{\left(0\right)}$[35].

Theorem 1. 

Assume that$T\left(Y^{\left(0\right)}\right)$and$T\left(Y^{\left(0\right)}D\right)$are theamplitudes of the random oscillation sequence$Y^{\left(0\right)}$and its smoothness sequence$Y^{\left(0\right)}D$, respectively, then$T\left(Y^{\left(0\right)}\right)>2T\left(Y^{\left(0\right)}D\right)$.

Proof. 

Assume that

$$\max \left\{y^{\left(0\right)}\left(k\right)|k=1,2,\cdots ,n\right\}=y^{\left(0\right)}\left(p\right),p=1,2,\cdots ,n,$$

$$\min \left\{y^{\left(0\right)}\left(k\right)|k=1,2,\cdots ,n\right\}=y^{\left(0\right)}\left(q\right),q=1,2,\cdots ,n.$$

Then,

$$T\left(Y^{\left(0\right)}\right)=y^{\left(0\right)}\left(p\right)-y^{\left(0\right)}\left(q\right).$$

(2)

Similarly, assume that,

$$\max \left\{y^{(0)}\left(k\right)d|k=1,2,\cdots ,n-1\right\}=y^{(0)}\left(i\right)d,i=1,2,\cdots ,n-1,$$

$$\min \left\{y^{(0)}\left(k\right)d|k=1,2,\cdots ,n-1\right\}=y^{(0)}\left(j\right)d,j=1,2,\cdots ,n-1.$$

According to Definition 2,

$$T\left(Y^{\left(0\right)}D\right)=y^{\left(0\right)}\left(i\right)d-y^{\left(0\right)}\left(j\right)d,$$

(3)

$$y^{\left(0\right)}\left(i\right)d=\frac{\left[y^{\left(0\right)}\left(i\right)+T\left(Y^{\left(0\right)}\right)\right]+\left[y^{\left(0\right)}\left(i+1\right)+T\left(Y^{\left(0\right)}\right)\right]}{4},$$

(4)

$$y^{\left(0\right)}\left(j\right)d=\frac{\left[y^{\left(0\right)}\left(j\right)+T\left(Y^{\left(0\right)}\right)\right]+\left[y^{\left(0\right)}\left(j+1\right)+T\left(Y^{\left(0\right)}\right)\right]}{4},$$

(5)

$$T\left(X^{\left(0\right)}D\right)=\frac{y^{\left(0\right)}\left(i\right)+y^{\left(0\right)}\left(i+1\right)-y^{\left(0\right)}\left(j\right)-y^{\left(0\right)}\left(j+1\right)}{4},$$

(6)

For

$$y^{\left(0\right)}\left(i\right)+y^{\left(0\right)}\left(i+1\right)-y^{\left(0\right)}\left(j\right)-y^{\left(0\right)}\left(j+1\right)\le 2T\left(Y^{\left(0\right)}\right),$$

(7)

then,

$$T\left(Y^{\left(0\right)}\right)\ge 2T\left(Y^{\left(0\right)}D\right)$$

(8)

□

The amplitude of a random oscillation sequence can be compressed using a smoothness operator. The smoothness operator D can improve the smoothness of a random oscillation sequence, which provides better modelling conditions and helps to build a more reasonable grey forecasting model.

The smoothness operator proposed in this paper does not meet the three axioms of a grey buffer operator: Axiom of fixed point, Axiom of in accordance with information, and Axiom of expressed normality. Hence, it is not a grey buffer operator. Before building a grey prediction model, if a raw sequence meets the random oscillation sequence condition in Definition 1, the algorithm to improve the smoothness of the random oscillation sequence should be applied first.

3. Grey Prediction Model with Three Parameters

Definition 3. 

Assume that a sequence is

$$Y^{\left(0\right)}=\left(y^{\left(0\right)}\left(1\right),y^{\left(0\right)}\left(2\right),\cdots ,y^{\left(0\right)}\left(m\right)\right),$$

where$y^{(0)}\left(k\right)\ge 0$,$k=1,2,\cdots ,m$. $Y^{\left(1\right)}$is the 1-accumulating generation operator (AGO) sequence of$Y^{\left(0\right)}$, i.e.,

$$Y^{\left(1\right)}=\left(y^{\left(1\right)}\left(1\right),y^{\left(1\right)}\left(2\right),\cdots ,y^{\left(1\right)}\left(m\right)\right),$$

where

$$y^{\left(1\right)}\left(k\right)={\displaystyle \sum _{i=1}^k{y}^{\left(0\right)}\left(i\right)},k=1,2,\cdots m.$$

From sequence$Y^{\left(1\right)}$, a new sequence$Z^{\left(1\right)}$can be derived:

$$Z^{\left(1\right)}=\left(z^{\left(1\right)}\left(2\right),z^{\left(1\right)}\left(3\right),\cdots ,z^{\left(1\right)}\left(m\right)\right),$$

where

$$z^{\left(1\right)}\left(k\right)=0.5\times \left[y^{\left(1\right)}\left(k\right)+y^{\left(1\right)}\left(k-1\right)\right],k=2,3\cdots ,m.$$

$Z^{\left(1\right)}$, generated by the consecutive neighbor of$Y^{\left(1\right)}$, is called the adjacent neighbor mean generation sequence.

Definition 4. 

For$Y^{(0)}$,$Y^{(1)}$, and$Z^{\left(1\right)}$given by Definition 3, and three constants u, v, and w, then the following equation

$$y^{\left(0\right)}\left(k\right)+u{z}^{\left(1\right)}\left(k\right)=vk+w,$$

(9)

is the expanded form of the GM(1,1) model.

From Definitions 3 and 4, the equations as follows can be deduced:

$$\{\begin{array}{l}{y}^{\left(0\right)}\left(k\right)+u{z}^{\left(1\right)}\left(k\right)=vk+w\\ y^{\left(0\right)}\left(k\right)=y^{\left(1\right)}\left(k\right)-y^{\left(1\right)}\left(k-1\right)\\ z^{\left(1\right)}\left(k\right)=0.5\left(y^{\left(1\right)}\left(k\right)+y^{\left(1\right)}\left(k-1\right)\right)\end{array}.$$

That is,

$$\left(1+0.5u\right)y^{\left(1\right)}\left(k\right)=vk+w+\left(1-0.5a\right)y^{\left(1\right)}\left(k-1\right).$$

Then, the following equation can be obtained:

$$y^{\left(1\right)}\left(k\right)=\frac{1-0.5u}{1+0.5u}{y}^{\left(1\right)}\left(k-1\right)+\frac{v}{1+0.5u}k+\frac{w}{1+0.5u}.$$

(10)

Letting

$${\chi }_1=\frac{1-0.5u}{1+0.5u},{\chi }_2=\frac{v}{1+0.5u},{\chi }_3=\frac{w}{1+0.5u}.$$

Equation (10) is converted into:

$$y^{\left(1\right)}\left(k\right)={\chi }_1{y}^{\left(1\right)}\left(k-1\right)+{\chi }_{2}k+{\chi }_3,k=2,3,\cdots ,m.$$

(11)

Equation (11) is the time response function of sequence $Y^{\left(1\right)}$.

In this subsection, the ordinary least square (OLS) method and Cramer’s rule are employed to estimate parameters ${\chi }_1$, ${\chi }_2$, and ${\chi }_3$ in Equation (11). After this, parameters $u$, $v$, and $w$ in Equation (11) can be calculated.

To minimize the simulation error S under the condition ${\widehat{y}}^{\left(1\right)}\left(k\right)$ being the simulation value of $y^{\left(1\right)}\left(k\right)$, the following condition must be satisfied:

$$\begin{array}{ll}S& =\min {{\displaystyle \sum _{k=2}^n\left[y^{\left(1\right)}\left(k\right)-{\widehat{y}}^{\left(1\right)}\left(k\right)\right]}}^2\\ & =\min {{\displaystyle \sum _{k=2}^n\left[y^{\left(1\right)}\left(k\right)-{\chi }_1{\widehat{y}}^{\left(1\right)}\left(k-1\right)-{\chi }_{2}k-{\chi }_3\right]}}^2\end{array}$$

Based on the OLS method, $S$ is minimized in relation to parameters ${\chi }_1,{\chi }_2,{\chi }_3$ to obtain:

$$\{\begin{array}{l}\frac{\partial S}{\partial {\chi }_1}=-2{\displaystyle \sum _{k=2}^n\left[y^{\left(1\right)}\left(k\right)-{\chi }_1{\widehat{y}}^{\left(1\right)}\left(k-1\right)-{\chi }_{2}k-{\chi }_3\right]\cdot }{\widehat{y}}^{\left(1\right)}\left(k-1\right)=0\\ \frac{\partial S}{\partial {\chi }_2}=-2{\displaystyle \sum _{k=2}^n\left[y^{\left(1\right)}\left(k\right)-{\chi }_1{\widehat{y}}^{\left(1\right)}\left(k-1\right)-{\chi }_{2}k-{\chi }_3\right]k=0}\\ \frac{\partial S}{\partial {\chi }_3}=-2{\displaystyle \sum _{k=2}^n\left[y^{\left(1\right)}\left(k\right)-{\chi }_1{\widehat{y}}^{\left(1\right)}\left(k-1\right)-{\chi }_{2}k-{\chi }_3\right]=0}\end{array}.$$

According to the formulas above, an equation set is acquired as follows:

$$\{\begin{array}{l}{\chi }_1{\displaystyle \sum _{k=2}^n{\left[{\widehat{y}}^{\left(1\right)}\left(k-1\right)\right]}^2+{\chi }_2{\displaystyle \sum _{k=2}^n\left[k{\widehat{y}}^{\left(1\right)}\left(k-1\right)\right]}}+{\chi }_3{\displaystyle \sum _{k=2}^n{\widehat{y}}^{\left(1\right)}\left(k-1\right)={\displaystyle \sum _{k=2}^n\left[y^{\left(1\right)}\left(k\right){\widehat{y}}^{\left(1\right)}\left(k-1\right)\right]}}\\ {\chi }_1{\displaystyle \sum _{k=2}^{n}k{\widehat{y}}^{\left(1\right)}\left(k-1\right)+{\chi }_2{\displaystyle \sum _{k=2}^n{k}^2+{\chi }_3{\displaystyle \sum _{k=2}^{n}k}={\displaystyle \sum _{k=2}^{n}k{y}^{\left(1\right)}\left(k\right)}}}\\ {\chi }_1{\displaystyle \sum _{k=2}^n{\widehat{y}}^{\left(1\right)}\left(k-1\right)+{\chi }_2{\displaystyle \sum _{k=2}^{n}k+{\chi }_3\left(n-1\right)={\displaystyle \sum _{k=2}^n{y}^{\left(1\right)}\left(k\right)}}}\end{array}$$

(12)

Then, the unknown parameters ${\chi }_1,{\chi }_2,{\chi }_3$ in Equation (12) are calculated as follows. Based on Cramer’s rule, the following results are obtained:

$$\begin{array}{l}\Phi =\left|\begin{array}{ccc}{\displaystyle \sum _{k=2}^n{\left[{\widehat{y}}^{\left(1\right)}\left(k-1\right)\right]}^2}& {\displaystyle \sum _{k=2}^{n}k{\widehat{y}}^{\left(1\right)}\left(k-1\right)}& {\displaystyle \sum _{k=2}^n{\widehat{y}}^{\left(1\right)}\left(k-1\right)}\\ {\displaystyle \sum _{k=2}^{n}k{\widehat{y}}^{\left(1\right)}\left(k-1\right)}& {\displaystyle \sum _{k=2}^n{k}^2}& {\displaystyle \sum _{k=2}^{n}k}\\ {\displaystyle \sum _{k=2}^n{\widehat{y}}^{\left(1\right)}\left(k-1\right)}& {\displaystyle \sum _{k=2}^{n}k}& n-1\end{array}\right|\\ {\Phi }_1=\left|\begin{array}{ccc}{\displaystyle \sum _{k=2}^n\left[y^{\left(1\right)}\left(k\right){\widehat{y}}^{\left(1\right)}\left(k-1\right)\right]}& {\displaystyle \sum _{k=2}^{n}k{\widehat{y}}^{\left(1\right)}\left(k-1\right)}& {\displaystyle \sum _{k=2}^n{\widehat{y}}^{\left(1\right)}\left(k-1\right)}\\ {\displaystyle \sum _{k=2}^{n}k{y}^{\left(1\right)}\left(k\right)}& {\displaystyle \sum _{k=2}^n{k}^2}& {\displaystyle \sum _{k=2}^{n}k}\\ {\displaystyle \sum _{k=2}^n{y}^{\left(1\right)}\left(k\right)}& {\displaystyle \sum _{k=2}^{n}k}& n-1\end{array}\right|\\ {\Phi }_2=\left|\begin{array}{ccc}{\displaystyle \sum _{k=2}^n{\left[{\widehat{y}}^{\left(1\right)}\left(k-1\right)\right]}^2}& {\displaystyle \sum _{k=2}^n\left[y^{\left(1\right)}\left(k\right){\widehat{y}}^{\left(1\right)}\left(k-1\right)\right]}& {\displaystyle \sum _{k=2}^n{\widehat{y}}^{\left(1\right)}\left(k-1\right)}\\ {\displaystyle \sum _{k=2}^{n}k{\widehat{y}}^{\left(1\right)}\left(k-1\right)}& {\displaystyle \sum _{k=2}^{n}k{\widehat{y}}^{\left(1\right)}\left(k\right)}& {\displaystyle \sum _{k=2}^{n}k}\\ {\displaystyle \sum _{k=2}^n{\widehat{y}}^{\left(1\right)}\left(k-1\right)}& {\displaystyle \sum _{k=2}^n{\widehat{y}}^{\left(1\right)}\left(k\right)}& n-1\end{array}\right|\\ {\Phi }_3=\left|\begin{array}{ccc}{\displaystyle \sum _{k=2}^n{\left[{\widehat{y}}^{\left(1\right)}\left(k-1\right)\right]}^2}& {\displaystyle \sum _{k=2}^{n}k{\widehat{y}}^{\left(1\right)}\left(k-1\right)}& {\displaystyle \sum _{k=2}^n\left[x^{\left(1\right)}\left(k\right){\widehat{y}}^{\left(1\right)}\left(k-1\right)\right]}\\ {\displaystyle \sum _{k=2}^{n}k{\widehat{y}}^{\left(1\right)}\left(k-1\right)}& {\displaystyle \sum _{k=2}^n{k}^2}& {\displaystyle \sum _{k=2}^{n}k{\widehat{y}}^{\left(1\right)}\left(k\right)}\\ {\displaystyle \sum _{k=2}^n{\widehat{y}}^{\left(1\right)}\left(k-1\right)}& {\displaystyle \sum _{k=2}^{n}k}& {\displaystyle \sum _{k=2}^n{\widehat{y}}^{\left(1\right)}\left(k\right)}\end{array}\right|\end{array}$$

According to Cramer’s rule, the parameters ${\chi }_1,{\chi }_2,{\chi }_3$ can be computed, as shown below, and the parameters $u,v,w$ can be obtained at the same time, as follows:

$$\{\begin{array}{l}{\chi }_1=\frac{{\Phi }_1}{\Phi }=\frac{1-0.5u}{1+0.5u}\\ {\chi }_2=\frac{{\Phi }_2}{\Phi }=\frac{v}{1+0.5a}\\ {\chi }_3=\frac{{\Phi }_3}{\Phi }=\frac{w}{1+0.5u}\end{array}\Rightarrow \{\begin{array}{l}u=\frac{2-2{\chi }_1}{1+{\chi }_1}\\ v=\frac{2{\chi }_2}{1+{\chi }_1}\\ w=\frac{2{\chi }_3}{1+{\chi }_1}\end{array}.$$

Then, the expression ${\widehat{y}}^{\left(1\right)}\left(k\right)$, which is the simulation value of $y^{\left(1\right)}\left(k\right)$, can be transformed into:

$${\widehat{y}}^{\left(1\right)}\left(k\right)=\frac{1-0.5u}{1+0.5u}{\widehat{y}}^{\left(1\right)}\left(k-1\right)+\frac{vk+w}{1+0.5u}.$$

(13)

Given Equation (11), we can obtain equations as follows when k equals 2 and 3, respectively:

$$y^{\left(1\right)}\left(2\right)={\chi }_1{y}^{\left(1\right)}\left(1\right)+2{\chi }_2+{\chi }_3,$$

(14)

$$y^{\left(1\right)}\left(3\right)={\chi }_1{y}^{\left(1\right)}\left(2\right)+3{\chi }_2+{\chi }_3.$$

(15)

The following equation can be obtained after the substitution of Equation (14) into Equation (15) and reorganization these two equations:

$$y^{\left(1\right)}\left(3\right)={\chi }_1\left[{\chi }_1{y}^{\left(1\right)}\left(1\right)+2{\chi }_2+{\chi }_3\right]+3{\chi }_2+{\chi }_3.$$

That is,

$$y^{\left(1\right)}\left(3\right)={\chi }_{{}_1}^2{y}^{\left(1\right)}\left(1\right)+2{\chi }_1{\chi }_2+{\chi }_1{\chi }_3+3{\chi }_2+{\chi }_3.$$

(16)

According to Equation (11), when k = 4:

$$y^{\left(1\right)}\left(4\right)={\chi }_1{y}^{\left(1\right)}\left(3\right)+4{\chi }_2+{\chi }_3.$$

(17)

The following equation can be obtained after substituting Equation (16) into (17) and reorganization of these two equations according to the method shown above:

$$y^{\left(1\right)}\left(4\right)={\chi }_{{}_1}^3{y}^{\left(1\right)}\left(1\right)+2{\chi }_{{}_1}^2{\chi }_2+{\chi }_{{}_1}^2{\chi }_3+3{\chi }_1{\chi }_2+{\chi }_1{\chi }_3+4{\chi }_2+{\chi }_3,$$

(18)

$${\widehat{y}}^{\left(1\right)}\left(j\right)=y^{\left(1\right)}\left(1\right)\cdot {\chi }_1^{\left(j-1\right)}+{\displaystyle \sum _{g=0}^{j-2}\left[\left(j-g\right)\cdot {\chi }_2+{\chi }_3\right]\cdot {\chi }_1^{\left(g\right)}},j=2,3,\cdots .$$

(19)

Equation (19) contains three parameters, ${\chi }_1,{\chi }_2,{\chi }_3$, and is therefore called the grey forecasting model with three parameters [36], or GFM_TP.

4. Improved GFM_TP Model (IGFM_TP) Based on a Smoothness Operator

Assume a random oscillation sequence:

$$Y=\left(y\left(1\right),y\left(2\right),\cdots ,y\left(n\right),y\left(n+1\right)\right).$$

Based on Definition 2, the smoothness sequence $YD$ of $Y$ is:

$$YD=H=\left(h\left(1\right),h\left(2\right),\cdots ,h\left(n\right)\right),$$

where

$$h\left(k\right)=y\left(k\right)d=0.25y\left(k\right)+0.25y\left(k+1\right)+0.5T_{\mathrm{amplitude}},$$

and $k=1,2,\cdots ,n$ and $T_{\mathrm{amplitude}}$ is the amplitude of the raw sequence $Y$.

Now, the IGFM_TP model with the smoothness sequence $H=\left(h\left(1\right),h\left(2\right),\cdots ,h\left(n\right)\right)$ is built. Given Equation (12), we obtain:

$${\widehat{h}}^{\left(1\right)}\left(j\right)=h^{\left(1\right)}\left(1\right)\cdot {\chi }_1^{\left(j-1\right)}+{\displaystyle \sum _{g=0}^{j-2}\left[\left(j-g\right)\cdot {\chi }_2+{\chi }_3\right]\cdot {\chi }_1^{\left(g\right)}}.$$

(20)

Based on Definition 3,

$${\widehat{h}}^{\left(0\right)}\left(j\right)={\widehat{h}}^{\left(1\right)}\left(j\right)-{\widehat{h}}^{\left(1\right)}\left(j-1\right),j=2,3,...,n.$$

(21)

However, the predictive data are $\widehat{y}\left(k\right)$ $\left(k=3,4,...,n+1\right)$, not ${\widehat{h}}^{\left(0\right)}\left(j\right)\left(j=2,3,...,n\right)$. Hence, we needed to deduce $\widehat{y}\left(k\right)$ from ${\widehat{h}}^{\left(0\right)}\left(j\right)$ according to Equation (1); that is,

$$y^{\left(0\right)}\left(t\right)d=0.25y^{\left(0\right)}\left(t\right)+0.25y^{\left(0\right)}\left(t+1\right)+0.5T_{\mathrm{amplitude}},$$

then,

$${\widehat{y}}^{\left(0\right)}\left(k+1\right)=4{\widehat{h}}^{\left(0\right)}\left(k\right)-{\widehat{y}}^{\left(0\right)}\left(k\right)-2T_{\mathrm{amplitude}},$$

(22)

When $k=1$, let

$${\widehat{y}}^{\left(0\right)}\left(2\right)=4{\widehat{h}}^{\left(0\right)}\left(1\right)-y^{\left(0\right)}\left(1\right)-2T_{\mathrm{amplitude}}.$$

(23)

Since ${\widehat{h}}^{\left(0\right)}\left(1\right)=h^{\left(0\right)}\left(1\right),y^{\left(0\right)}\left(1\right)$ and $T_{\mathrm{amplitude}}$ are constants, $y^{\left(0\right)}\left(1\right)$ is also a constant. Then ${\widehat{y}}^{\left(0\right)}\left(2\right)=y^{\left(0\right)}\left(2\right)$ is called the initial value of the new proposed IGFM_TP model.

When $k=2$,

$${\widehat{y}}^{\left(0\right)}\left(3\right)=4{\widehat{h}}^{\left(0\right)}\left(2\right)-{\widehat{y}}^{\left(0\right)}\left(2\right)-2T_{\mathrm{amplitude}},$$

(24)

The following equation can be obtained when Equation (16) is substituted into Equation (17) and both equations are reorganized:

$$\begin{array}{ll}{\widehat{y}}^{\left(0\right)}\left(3\right)& =4{\widehat{h}}^{\left(0\right)}\left(2\right)-\left(4{\widehat{h}}^{\left(0\right)}\left(1\right)-y^{\left(0\right)}\left(1\right)-2T_{\mathrm{amplitude}}\right)-2T_{\mathrm{amplitude}}\\ & =4{\widehat{h}}^{\left(0\right)}\left(2\right)-4{\widehat{h}}^{\left(0\right)}\left(1\right)+y^{\left(0\right)}\left(1\right)+2T_{\mathrm{amplitude}}-2T_{\mathrm{amplitude}}\end{array}$$

(25)

That is:

$${\widehat{y}}^{\left(0\right)}\left(3\right)=4{\widehat{h}}^{\left(0\right)}\left(2\right)-4{\widehat{h}}^{\left(0\right)}\left(1\right)+y^{\left(0\right)}\left(1\right).$$

(26)

When $k=3$,

$${\widehat{y}}^{\left(0\right)}\left(4\right)=4{\widehat{h}}^{\left(0\right)}\left(3\right)-{\widehat{y}}^{\left(0\right)}\left(3\right)-2T_{\mathrm{amplitude}}.$$

(27)

The following equation can be obtained when Equation (19) is substituted into Equation (20) and both equations are reorganized:

$${\widehat{y}}^{\left(0\right)}\left(4\right)=4{\widehat{h}}^{\left(0\right)}\left(3\right)-4{\widehat{h}}^{\left(0\right)}\left(2\right)+4{\widehat{h}}^{\left(0\right)}\left(1\right)-y^{\left(0\right)}\left(1\right)-2T_{\mathrm{amplitude}}.$$

(28)

Since

$${\widehat{h}}^{\left(1\right)}\left(3\right)=h^{\left(1\right)}\left(1\right)\cdot {\chi }_1^{\left(2\right)}+{\displaystyle \sum _{g=0}^1\left[\left(3-g\right)\cdot {\chi }_2+{\chi }_3\right]\cdot {\chi }_1^{\left(g\right)}},$$

(29)

$${\widehat{h}}^{\left(1\right)}\left(2\right)=h^{\left(1\right)}\left(1\right)\cdot {\chi }_1^{\left(1\right)}+{\displaystyle \sum _{g=0}^0\left[\left(2-g\right)\cdot {\chi }_2+{\chi }_3\right]\cdot {\chi }_1^{\left(g\right)}}.$$

(30)

Then

$$\begin{array}{l}{\widehat{y}}^{\left(0\right)}\left(4\right)=4h^{\left(1\right)}\left(1\right)\cdot {\chi }_1^{\left(2\right)}+4{\displaystyle \sum _{g=0}^1\left[\left(3-g\right)\cdot {\chi }_2+{\chi }_3\right]\cdot {\chi }_1^{\left(g\right)}-}4h^{\left(1\right)}\left(1\right)\cdot {\chi }_1^{\left(1\right)}+\\ 4{\displaystyle \sum _{g=0}^0\left[\left(2-g\right)\cdot {\chi }_2+{\chi }_3\right]\cdot {\chi }_1^{\left(g\right)}}-4h^{\left(1\right)}\left(1\right)\cdot {\chi }_1^{\left(1\right)}-4{\displaystyle \sum _{g=0}^0\left[\left(2-g\right)\cdot {\chi }_2+{\chi }_3\right]\cdot {\chi }_1^{\left(g\right)}+}\\ 4{\widehat{h}}^{\left(0\right)}\left(1\right)+4{\widehat{h}}^{\left(0\right)}\left(1\right)-y^{\left(0\right)}\left(1\right)-2T_{\mathrm{amplitude}}\end{array}$$

(31)

As can be seen from derivations demonstrated above, we failed to obtain regular evolution, because the time response function of the IGFM_TP model is exceedingly complicated. Actually, our main objective was to predict or simulate ${\widehat{y}}^{\left(1\right)}\left(k\right)$ and $y^{\left(0\right)}\left(k\right)$ by building the IGFM_TP model, and the time response function of the IGFM_TP model was not a vital factor for this study. A recursive program can be written to achieve the model, because the IGFM_TP model satisfies the recursive algorithm based on Equation (21).

5. Electricity Forecasting Using IGFM_TP

The electricity consumption of one city in Western China is high from January to July, as listed in

Table 1. Electricity consumption of a city in Western China.

Month	January	February	March	April	May	June	July
Electricity Consumption (Million Kilowatts)	439	320	584	481	640	635	790

. The electricity demand of the city was simulated and its simulation errors were compared with the proposed IGFM_TP model, the traditional GFM_TP model [37], the discrete grey model DGM(1,1) [38], and the classical GM(1,1) model.

Step 1. Collecting modeling data.

Based on Table 1, the original sequence $Y$ is the modelling data:

$$\begin{array}{ll}Y& =\left(y\left(1\right),y\left(2\right),y\left(3\right),y\left(4\right),y\left(5\right),y\left(6\right),y\left(7\right)\right)\\ & =\left(439,320,584,481,640,635,790\right)\end{array}$$

Step 2. Computing smoothness sequence.

Given of Definition 2, the smoothness sequence $Y^{\left(0\right)}$ of $Y$ is:

$$\begin{array}{ll}{Y}^{\left(0\right)}& =\left(y\left(1\right)d,y\left(2\right)d,y\left(3\right)d,y\left(4\right)d,y\left(5\right)d,y\left(6\right)d,y\left(7\right)d\right)\\ & =\left(424.75,461.0,501.25,515.25,553.75,591.25\right)\end{array}$$

From the sequence $Y^{\left(0\right)}$, when compared with the random oscillation sequence $Y$, the smoothness sequence is smoother, which is important for improving the simulation or prediction precision of grey prediction models.

Step 3. Computing the parameters of the IGFM_TP model.

Compute the parameters ${\chi }_1,{\chi }_2,{\chi }_3$ of the GFM_TP model as follows:

$${\chi }_1=1.099598,{\chi }_2=-19.285853,{\chi }_3=461.095615.$$

Step 4. Computing the simulated data and errors.

Firstly, the simulated data and errors of the IGFM_TP model were calculated according to Steps 1–4 and Equation (31). Four grey models were built to compare their simulated values and errors, as shown in

Table 2. Simulation values and errors with five different grey prediction models.

Month	${\mathit{y}}^{\left(0\right)}\left(\mathit{k}\right)$	IGFM_TP		Model in AUTHOR et al.		Classical GM(1,1)		Model in AUTHOR et al.		Model in AUTHOR et al.
		${\widehat{\mathit{y}}}^{\left(0\right)}\left(\mathit{k}\right)$	${\mathbf{\Delta }}_{\mathit{k}}$	${\widehat{\mathit{y}}}^{\left(0\right)}\left(\mathit{k}\right)$	${\mathbf{\Delta }}_{\mathit{k}}$	${\widehat{\mathit{y}}}^{\left(0\right)}\left(\mathit{k}\right)$	${\mathbf{\Delta }}_{\mathit{k}}$	${\widehat{\mathit{y}}}^{\left(0\right)}\left(\mathit{k}\right)$	${\mathbf{\Delta }}_{\mathit{k}}$	${\widehat{\mathit{y}}}^{\left(0\right)}\left(\mathit{k}\right)$	${\mathbf{\Delta }}_{\mathit{k}}$
February	320	320.0	0.0%	409.4	27.9%	405.2	26.6%	371.0	16.0%	334.7	4.6%
March	584	599.3	2.6%	465.1	20.4%	461.5	21.0%	464.5	20.5%	487.5	16.5%
April	481	428.1	11.0%	528.5	9.9%	525.6	9.3%	548.4	14.1%	583.2	21.2%
May	640	718.1	12.2%	600.5	6.2%	598.7	6.5%	623.5	2.6%	643.0	0.50%
June	635	558.7	12.0%	682.4	7.5%	681.9	7.4%	691.0	8.8%	680.5	7.2%
July	790	861.8	9.1%	775.3	1.9%	776.7	1.7%	751.4	4.9%	704.0	10.9%
$\Delta$			7.8%	-	12.3%	-	12.1%	-	11.1%		10.1%

.

In Table 2, $y^{\left(0\right)}\left(k\right)$ is the k-actual value in the original sequence $X$, ${\widehat{y}}^{\left(0\right)}\left(k\right)$ is the simulative or predictive value corresponding to $y^{\left(0\right)}\left(k\right)$, and ${\Delta }_k$ is the absolute simulation or prediction percentage error of ${\widehat{y}}^{\left(0\right)}\left(k\right)$.

$${\Delta }_k=\frac{\left|{\widehat{y}}^{\left(0\right)}\left(k\right)-y^{\left(0\right)}\left(k\right)\right|}{y^{\left(0\right)}\left(k\right)}\times 100\%,$$

where $\Delta$ is the mean absolute simulation or prediction percentage error of sequence $\widehat{Y}$.

$$\Delta =\frac{1}{n-1}{\displaystyle \sum _{k=2}^n{\Delta }_k}\times 100\%.$$

In accordance with statistics in Table 2, the simulated data obtained from the above five grey prediction models are depicted as dotted lines in Figure 2.

Compared with the curves of the five grey prediction models, Figure 2 shows that the curve of the simulated data of the IGFM_TP model is closest to that of the raw data, which demonstrates IGFM_TP has the best simulation precision. In other words, the proposed IGFM_TP model is effective. The curves of the simulated data with the GM(1,1), DGM(1,1), GFM_TP, and SAIGM models all monotonically increased, and they cannot be used to simulate a random oscillation sequence, because their most recent restored forms demonstrate homogeneous or non-homogeneous exponential forms.

The size of the mean absolute simulation percentage error (MASPE) is often used to test the robustness and soundness of a grey prediction model. The smaller the MASPE, the better the performance of the grey model. In this paper, the MASPE grade of the IGFM_TP model was close to Level-II, which means that the new model has good simulation performance. However, the MASPE grades of the four other models were all greater than Level-III, which proves again that the performance of the new model is better than that of the other grey prediction models.

6. Conclusions

Electricity is the foundation of a country’s economic construction. Reasonable and effective prediction of electricity demand is important for ensuring a balance between electricity supply and demand. For this, a new grey prediction model, named IGFM_TP, was proposed in this paper. A smoothness algorithm was applied in the IGFM_TP model to process a raw data sequence, and then a smoothness sequence was obtained. After this, a three-parameters grey prediction model was established based on the smoothness sequence. Hence, the IGFM_TP model was built with two existing methods and model.

The IGFM_TP model was used to simulate the electricity consumption of one city in Western China. The simulation results showed that the performance of IGFM_TP is not only better than the classical GM(1,1) model and frequently-used DGM(1,1) model, but also the improved upon the SAIGM and GFM_TP models, which confirms the effectiveness of the performance improvement. Detailed modeling steps were provided in this paper and the new model can be used to solve some actual issues in real life according to the steps. However, the proposed new model is not perfect and still requires considerable optimization, such as the initial and background values and the order of accumulation. Therefore, further optimizing and improving the simulation performance of the GFM_TP model is the next research goal of our team.