# A Novel Weighted Fractional GM(1,1) Model and Its Applications

Authors: Qin-Qin Shen; Quan Shi; Tian-Pei Tang; Lin-Quan Yao
Journal: Complexity
 (2020)

---

## Abstract

Based on the ideas of the new information priority principle and the fractional-accumulation generating operator, in this paper we propose a novel weighted fractional GM(1,1) (WFGM(1,1)) prediction model. In the new model, the original sequence is first transformed by using the weighted fractional-accumulation generating operator, which involves two parameters. With special choices of these parameters, the proposed WFGM(1,1) model reduces to the classical GM(1,1) model and the fractional GM(1,1) (FGM(1,1)) model, as well as the new information priority GM(1,1) (NIPGM(1,1)) model studied recently. Stability property of the WFGM(1,1) model is studied in detail. In practice, the quantum particle swarm optimization algorithm is adopted to choose the quasi-optimal parameters for the new model so as to get the best fitting accuracy. Finally, four numerical examples from different practical applications are present. Numerical results show that the new proposed prediction model is very efficient and has both the best fitting accuracy and the best prediction accuracy compared with the GM(1,1) and the FGM(1,1) as well as the NIPGM(1,1) prediction models.


---

## Body

## 1. Introduction

Given some historic data and the current data, fitting the data trend and predicting or forecasting what will happen in the near future is of great interest and has been an ongoing hot topic in many engineering managements and applications [1]. A great deal of effective prediction models has been proposed in the past few decades. Among numerous existing prediction models, the grey theory-based models, which were first proposed by Deng in 1982 [2], have attracted many researchers’ attention due to its simple form, high computing efficiency, and high prediction accuracy especially for limited data and insufficient information [3–6].The classical GM(1,1) model, in which the first one denotes the first order and the second one denotes the single variable, is the foundation of grey theory. Denote byX0 the original sequence. The basic idea of the GM(1,1) model is first to use the first order-accumulated generating operator (1-AGO) to transform the original sequence X0 into a monotonicity increased sequence X1. And then use the solution X^1 of a first order-ordinary differential equation, which is called the whitenization differential equation, to approximate the sequence X1. Finally, applying the first order-inverse accumulated generating operator (1-IAGO), one can get the fitted value X^0 of X0 and obtain the future trend. In the past few decades, it has been widely utilized in a great deal of engineering applications, such as the global integrated circuit industry prediction [7], the vehicle fatality risk estimation [8], the short-term freeway traffic parameter prediction [9], the fashion retailing prediction problem [10], the electricity consumption forecasting [11], and the natural gas consumption forecasting [12, 13]. To improve the prediction accuracy of the GM(1,1) model, many researchers have carried out a lot of works from different aspects, such as finding new accumulation generating operators [14–19], constructing more accurate background value formula [20, 21], choosing parameter optimization methods [22], improving initial guess [23], and reducing residuals based on Fourier analysis and Markov chain [9, 20]. Recently, some nonhomogeneous, nonlinear, hybrid, and multivariable grey models are proposed, see [6, 24–26] for examples. Modeling mechanism analysis can be found in [27–29]. For detailed descriptions of the GM(1,1) model and its improvements, we refer the readers to [4] and a very good survey published recently [5] and also references therein.According to the modeling process of the GM(1,1) prediction model, we can see that the 1-AGO, which generally reduces stochastic of the original data, is the first and the most important step [4, 5]. So, in this paper we mainly focus on the accumulation generating operator. To improve the precision of the GM(1,1) prediction model, many works have been carried out in the past few decades to propose new accumulation generating operators according to data sequences’ different characteristics so as to meet with the grey exponential law. For instance, for decreasing sequences, a reverse-AGO and a reciprocal-AGO were proposed, respectively, in [14, 15]. Based on random fluctuation characteristics of the data, Li and Xie designed a varied weighted-AGO with the arctangent function, which not only takes the overall growth trend but also the local fluctuation characteristics into account [19]. For fashion retailing prediction problem, in which the data sequence exhibits obvious periodicity, a cycle truncation accumulation generating operator was put forward to obtain a discrete GM(1,1) model [10]. Then, this idea was further applied to the short-term traffic flow prediction problem [30, 31]. According to the grey system theory, the effect of new information is greater than that of old information. This is called the principle of new information priority. However, the above accumulation generating operators do not actually reflect this principle. To remedy this, many researchers have proposed many new accumulation generating operators based on the idea of new information priority. These operators can be divided into two categories. On one hand, by setting different weights on different elements, a weighted-AGO (W-AGO) [16] and a new information priority-AGO (NIP-AGO) [18] were proposed to construct new GM(1,1) prediction models. In fact, the NIP-AGO can be regarded as a special case of the W-AGO and has been used in constructing multivariable prediction model [32] and discrete grey model [33]. On the other hand, by analyzing the modeling mechanism of r−AGO and setting the integer r to be a fractional number, a fractional-AGO (F-AGO) and the corresponding fractional GM(1,1) (FGM(1,1)) model were proposed by Wu et al. [17]. Based on the idea of F-AGO, several other grey models such as fractional GM(q,1) model [34], discrete fractional grey model [35, 36], multivariable fractional grey model [37], and self-adaptive intelligence fractional grey prediction model [38] have been put forward and successfully applied in many practical applications, including the fuel production [35], the CO2 emission [36], and the electricity consumption [37]. Numerical results show that the fractional grey prediction models have much higher precision than the traditional GM(1,1) model. Modeling mechanism of the fractional grey models has been studied, see [17, 29, 35, 39–41].In this paper, based on the principle of the new information priority and the F-AGO, we first try to construct a new accumulation generating operator and call it the weighted fractional-accumulation generating operator (WF-AGO). The new WF-AGO contains two parameters, which not only can reflect the new information priority but also adjust the order of the summation according to different data sequences. Then, a new weighted fractional GM(1,1) (WFGM(1,1)) model, which involves the classical GM(1,1), the NIPGM(1,1), and the FGM(1,1) prediction models as special cases, is constructed. In addition, by defining a nonlinear constrained optimization problem for fitted values, the quantum particle swarm optimization (QPSO) method is adopted to find the best parameters [42–44]. Theoretically, stability property of the WFGM(1,1) model is analyzed in detail. With four numerical examples from different practical problems, we show that both the fitting accuracy and the prediction accuracy of the proposed model are better than those of the original GM(1,1) model [3], the FGM(1,1) model [17], and the NIPGM(1,1) model [18] proposed recently.The remainder of this paper is organized as follows. In Section2, a new weighted fractional-accumulation generating operator and the weighted fractional GM(1,1) model are established. Parameter estimation methods are studied. In addition, the quantum particle swarm optimization algorithm is adopted to choose the optimal parameters. Computational steps of the new model are also presented. Stability property of the new WFGM(1,1) prediction model is analyzed. Theoretical results show that the WFGM(1,1) model is also suitable for small sample data. In Section 3, four cases from different applications are studied to verify the efficiency of the new proposed WFGM(1,1) model. In addition, comparisons of the new WFGM(1,1) model with the classical GM(1,1) model, the FGM(1,1) model, and the NIPGM(1,1) model are discussed in detail. Finally, we end this paper with some conclusions and outlooks in Section 4.


## 2. The Weighted Fractional GM(1,1) Model

### 2.1. The Weighted Fractional-Accumulation Generating Operator

In this section, a new weighted fractional-accumulation generating operator (WF-AGO) and the corresponding inverse accumulation generating operator (WF-IAGO) are defined.Definition 1 (see [34]).
LetN and N+ denote the set of natural numbers and the set of positive integers, respectively. Let r be a real number. Then, for n∈N, the general rising factorial function of a rational number is(1)rn=1,ifn=0,rr+1⋯r+n−1n!,ifn∈N+.
It has been studied in Theorem1 [41] that rn is a monotonically decreasing function and a monotonically increasing function with respect to n if r∈0,1 and r∈1,+∞, respectively. In particular, if r=1, then rn≡1. By letting r∈0,1 and introducing another weighted parameter λ∈0,1 used in the NIPGM(1,1) model [18], we define a new accumulation generating operator as follows.Definition 2.
Let(2)X0=x01,x02,…,x0nT,be the original nonnegative sequence, r,λ∈0,1 be two nonnegative parameters, and the weighted fractional-accumulation generating sequence of X0 is defined as(3)Xrλ=xrλ1,xrλ2,…,xrλnT,where(4)xrλk=∑i=1krk−iλk−ix0i,k=1,…,n.
According to Definition2, the weighted fractional-accumulation generating operator (WF-AGO) can be simply expressed by using the following matrix-vector form:(5)Xrλ=ArλX0,where Arλ is a unit lower triangular matrix and takes the form(6)Arλ=r000⋯0r1λr00⋯0r2λ2r1λr0⋯0⋮⋮⋮⋱⋮rn−1λn−1rn−2λn−2rn−3λn−3⋯r0.
Evidently, ifr=λ=1, then Arλ is a lower triangular matrix with all nonzero elements being one, which reduces to the well-known first order-accumulation generating operator (1-AGO) matrix. If r=1 and λ=1, then the weighted fractional-accumulation generating operator matrix Arλ reduces to(7)Aλ=100⋯0λ10⋯0λ2λ1⋯0⋮⋮⋮⋱⋮λn−1λn−2λn−3⋯1,Ar=r000⋯0r1r00⋯0r2r1r0⋯0⋮⋮⋮⋱⋮rn−1rn−2rn−3⋯r0,which are called the new information priority-accumulation generating operator (NIP-AGO) matrix [18] and the fractional-accumulation generating operator (F-AGO) matrix [17], respectively.Definition 3.
LetX0 and Xrλ be the original sequence and the weighted fractional-accumulation generating sequence defined as in (2) and (3), respectively. The weighted fractional-inverse accumulation generating operator (WF-IAGO) is(8)X0=Arλ−1Xrλ≜DrλXrλ.
With the notations given above, we can useDr and Dλ to denote the fractional-inverse accumulation generating operator (F-IAGO) [17] and the new information priority-inverse accumulation generating operator (NIP-IAGO) [18] matrices, which are the inverse matrices of Ar and Aλ, respectively. For the specific expression of Drλ, we have the following result.Theorem 1.
The WF-IAGO matrixDrλ satisfies(9)Drλ=−r000⋯0−r1λ−r00⋯0−r2λ2−r1λ−r0⋯0⋮⋮⋮⋱⋮−rn−1λn−1−rn−2λn−2−rn−3λn−3⋯−r0.Proof.
According to Definition1 and (6), we can see that the matrix Arλ is a unit lower triangular matrix. Therefore, detArλ=1 and its inverse matrix Drλ exits.
To prove that the matrixDrλ has form (9), we only need to show that ArλDrλ is an identity matrix. Since both Arλ and Drλ are unit lower triangular matrices, so is ArλDrλ. In addition, for i>j, we have(10)ArλDrλij=λi−jri−jλi−j−r0+ri−j−1λi−j−1−r1λ+⋯+r1λ−ri−j−1λi−j−1+r0−ri−jλi−j=λi−jArDrij.
From Theorem2 [41], we can see that ArDr is an identity matrix. Thus, for i>j we have ArλDrλij=0, which means that the matrix Drλ defined as in (9) is the inverse matrix of Arλ. Therefore, the proof is completed.


### 2.2. The Representation of the WFGM(1,1) Model

Based on the idea of the grey theory, the new WFGM(1,1) model is constructed in this section.Definition 4.
LetX0=x01,x02,…,x0nT and Xrλ=xrλ1,xrλ2,…,xrλnT be the original sequence (2) and the WF-AGO sequence (3), respectively. We call the following linear differential equation(11)dxrλtdt+axrλt=b,the whitenization equation of the WFGM(1,1) model, and the difference equation(12)xrλk−xrλk−1+azrλk=b,the basic model of the WFGM(1,1) model, where(13)zrλk=12xrλk+xrλk−1,k=2,3,…,n,is called the background value, and a and b are called the development coefficient and the grey input action, respectively.
From the definition of the WFGM(1,1) model and given the initial guessxrλtt=1=x01, we can easily obtain the solution of the whitenization equation (11):(14)x^rλk=x01−bae−ak−1+ba,k=1,2,…,n,which is called the time response function of the WFGM(1,1) model. Note that the time response function (14) is often used to approximate the values of the WF-AGO sequence Xrλ=xrλ1,…,xrλn. By using the WF-IAGO defined as in (8), the predicted values X^0=x^01,…,x^0n of the original sequence can be obtained:(15)x^0k=∑i=1k−rk−iλk−ix^rλi,k=1,2,…,n,or in matrix form, it satisfies(16)X^0=DrλX^rλ.
Obviously, the new WFGM(1,1) model is a generalization of the classical GM(1,1) model and its improvements studied recently. In fact, ifr=λ=1, then the WFGM(1,1) model reduces to the GM(1,1) model [3, 4]. If λ=1 and r=1, then the new WFGM(1,1) model reduces to the FGM(1,1) model [17] and the NIPGM(1,1) model [18], respectively. With the aid of two parameters r and λ, we can choose suitable values to not only meet with the principle of new information priority but also improve the prediction accuracy.


### 2.3. Parameters Estimation and Computational Steps

As can be seen from Section2.2, four parameters need to be determined in the WFGM(1,1) model. To obtain the WF-AGO sequence, we have to determine the parameters r and λ in advance. To get the solution of the WFGM(1,1) model, the development coefficient a and the grey input action b need to be fixed in the linear differential equation (11). Fortunately, all these parameters can be obtained by solving two optimization problems.Similar to the idea of the GM(1,1) model and its improvements, the parametersa and b can be obtained by minimizing the following unconstrained optimization problem:(17)mina,b∑k=2nxrλk−xrλk−1+azrλk−b2,which is constructed by difference equation (12). The solution of the above optimization problem can be simply written as the following linear system:(18)ab=BTB−1BTY,where(19)B=−zrλ21−zrλ31⋮⋮−zrλn1,Y=xrλ2−xrλ1xrλ3−xrλ2⋮xrλn−xrλn−1.Next, we turn to study how to determine the parametersr and λ, which play very important roles to make the WF-AGO sequence Xrλ better satisfy the grey exponential law and the principle of new information priority, and thus get better prediction accuracy. That is to say, the parameters must be chosen to get both the best fitting accuracy and the best prediction accuracy. There are many ways to define the model prediction accuracy, such as the mean absolute percentage error (MAPE), the mean absolute error (MAE), and the root mean square error (RMSE). Among these, the MAPE is one of the most common indexes to judge the merits of prediction models. Therefore, we use the MAPE as an objective function in this paper:(20)MAPEr,λ=1n∑k=1nx^0k−x0kx0k×100%.Theoretically, the MAPE is constrained by several conditions, such as the constraints on parametersr and λ and the time response function x^rλk and its inverse sequence x^0k, as well as the parameters a and b. In summary, the following constrained optimization problem can be defined to obtain the optimal parameters r and λ:(21)minMAPEr,λ,s.t0<λ≤1,0<r≤1,ab=BTB−1BTY,x^rλk=x01−bae−ak−1+ba,x^0k=∑i=1k−rk−iλk−ix^rλi,k=1,2,…,n.In general, the above optimization problem is nondifferentiable, nonlinear, and is very difficult to solve. To find its optimal solution, we need to determine the parametersa and b in advance. This is a contradiction. Fortunately, we can use the quantum particle swarm optimization (QPSO) algorithm to fast get a quasi-optimal solution of (21). Note that the QPSO algorithm is a valuable development and quality improvement of the classical PSO algorithm. The greatest advantage of the QPSO algorithm over the PSO algorithm is that the QPSO ameliorates the bug of PSO preferably and improves the search efficiency [42]. Now, it has been used to solve such nonlinear constrained optimization problems [42–44]. In the following, we present a brief introduction to the QPSO algorithm to solve the nonlinear constrained optimization problem (21). Matlab code for the QPSO algorithm can be found in Appendix A.Since there are two variables in (21), only a two-dimensional problem needs to be considered. Assume that there are m discrete particles distributed in the two-dimensional space. At the lth iteration, they have the form(22)El=E1l,E2l,…,Eml=E11lE21l⋯Em1lE12lE22l⋯Em2l,where Eil i=1,…,m denotes the position of the ith particle and Ei1l,Ei2lT denotes its coordinate. Here, Ei1l and Ei2l can be regarded as the choice for the parameters r and λ, respectively. At the initial iteration, i.e., l=1, let Pi1=Ei1 (i=1,…,m). Then, at the lth iteration, the optimal position for the ith particle (often called “pbest”) is obtained by computing(23)pbest=Pil=Eil,if MAPEEil<MAPEPil−1,Pil−1,if MAPEEil≥MAPEPil−1,i=1,…,m,and the global optimal position, which is called “gbest” and can be regarded as the global quasi-optimal solution for the nonlinear constrained optimization problem (21), is obtained:(24)gbest=Pgl=Pg1lPg2l≜r∗lλ∗l,where(25)g=argmin1≤i≤mMAPEPil.At thel+1th iteration, the position Ei1l+1,Ei2l+1T for each particle is updated by the following rule:(26)Eidl+1=Qidl±αCdl−Eidl⋅ln1uidl,uidl∼U0,1,d=1,2,where(27)Qidl=ϕdlPidl+1−ϕdlPgdl,ϕd∼U0,1,Cdl=1m∑i=1mPidl,α=0.5−0.5lmax−ll,and lmax is the given largest iteration number.According to the above discussion, the computational steps of the new WFGM(1,1) model can be summarized as follows. The corresponding Matlab code can be found in AppendixB and Appendix C.Step 1: determine the original data sequenceX0 according to the actual application background. Set the initial parameters r and λ, as well as the largest iteration number lmax for the QPSO algorithm.Step 2: compute the WF-AGO sequenceXrλ=ArλX0 and the background value zrλk according to (3) and (13), respectively.Step 3: solve the unconstrained optimization problem (17) or its equivalent form, i.e., the least squares problem (18), to obtain the parameters a and b.Step 4: substitute the parametersa and b into (14) to compute the approximated values X^rλ of the WF-AGO sequence Xrλ and compute the predicted values X^0 of the WFGM(1,1) model according to (16).Step 5: determine the local optimal parametersr and λ to minimize the MAPE (21) by using the QPSO algorithm.Step 6: repeat Steps 2–5 to obtain the global quasi-optimal (“gbest”) parameters r∗, λ∗, and the corresponding optimal fitted and predicted values X^0.


### 2.4. Stability Analysis

Stability is a very important property for a prediction model since some data may be not accurate in many engineering applications due to several reasons. For the original GM(1,1) model, Wu et al. first studied its stability by using the tools of matrix perturbation analysis. They defined a perturbation bound to show how the solution is affected by the change of the original dataX0. In this section, we extend this idea to study the stability property for the new WFGM(1,1) model.To this end, we first present a useful lemma and give the definition of the perturbation bound, which can be found in [27] and some recent works. Here and in the sequel, ⋅ denotes the 2-norm of either a vector or a matrix. rank⋅ stands for the rank of the corresponding matrix.Lemma 1 (see [27, 45]).
Suppose thatx and x^ satisfy(28)Bx−Y=min,B+ΔBx^−Y+ΔY=min,respectively, where B,ΔB∈ℂm×n with m≤n, and Y≠0, ΔY∈ℂm. Let κ†=B†⋅ΔB, γ†=1−B†⋅ΔB, and rx=Y−Bx, B† be the pseudoinverse of the matrix B. If rankB=rankB+ΔB and B†⋅ΔB<1, then(29)x^−x≤κ†γ†ΔBBx+ΔYB+κ†γ†ΔBBrxB.Definition 5 (see [27]).
Assume thatX0=x01,x02,…,x0n is the original nonnegative sequence. Let B and Y be defined as in (19) and B†⋅ΔB<1. We call(30)Lx0k=κ†γ†ΔBBx+ΔYB+κ†γ†ΔBBrxB,k=1,2,…,n,the perturbation bound when ε is a disturbance of x0k k=1,2,…,n.
Analogous to the analysis in [27], we have the following results for the new WFGM(1,1) model.Theorem 2.
Assume thatX0=x01,x02,…,x0nT is the original nonnegative sequence and ε is a disturbance of x0k k=1,2,…,n. Let B and Y be defined as (19). Denote by(31)θi=ri−1λi−1+riλi2,φi=riλi−ri−1λi−12.
Then, we have(32)Lx0k=κ†γ†ε∑i=1n−1θi2Bx+∑i=1n−1φiB+κ†γ†∑i=1n−1θi2BrxB,k=1,κ†γ†ε1+∑i=1n−kθi2Bx+1+∑i=1n−kφiB+κ†γ†1+∑i=1n−kθi2BrxB,k=2,⋯,n−1,κ†γ†εx2B+1B+κ†γ+12BrxB,k=n,and Lx0kk=1,2,…,n−1 is an increasing function with respect to n.Proof.
To deduce the perturbation bound for each data, we first simplify the expression of the matrixB and the vector Y. Let e=1,1,…,1T∈ℝn. Since Xrλ=ArλX0 and e=ArλDrλe, we have(33)B=−zrλ21−zrλ31⋮⋮−zrλn1n−1×2=−12−120⋯000−12−12⋯00⋮⋮⋮⋱⋮⋮000⋯−12−12n−1×nxrλ1−1xrλ2−1⋮⋮xrλn−1n×2≜B1ArλX0−Drλe,Y=xrλ2−xrλ1xrλ3−xrλ2⋮xrλn−xrλn−1n−1×1=−110⋯000−11⋯00⋮⋮⋮⋱⋮⋮000⋯−11n−1×nxrλ1xrλ2⋮xrλnn×1≜B2ArλX0.(1)
Ifε is regarded as a disturbance of x01, i.e., X10ε=x01+ε,x02,…,x0nT. Then, we have(34)B+ΔB=B1ArλX10ε−Drλe=B+B1Arλε000⋮⋮00,ΔB=−12r0+r1λε0−12r1λ+r2λ2ε0⋮⋮−12rn−2λn−2+rn−1λn−1ε0,Y+ΔY=B2ArλX10ε=Y+B2Arλε0⋮0,ΔY=r1λ−r0εr2λ2−r1λε⋮rn−1λn−1−rn−2λn−2ε.Thus, it satisfies(35)ΔBTΔB=ε24∑i=1n−1ri−1λi−1+riλi2000,ΔYTΔY=ε2∑i=1n−1riλi−ri−1λi−12.Based on the definition of 2-norm, from (35) we have(36)ΔB=ε2∑i=1n−1ri−1λi−1+riλi2≜ε2∑i=1n−1θi,(37)ΔY=ε∑i=1n−1riλi−ri−1λi−12≜ε∑i=1n−1φi,where the notations in (31) are used. Substituting (36) and (37) into Lemma 1 (30), we immediately have the disturbance bound(38)Lx01=κ†γ†ε∑i=1n−1θi2Bx+∑i=1n−1φiB+κ†γ†∑i=1n−1θi2BrxB,whenε is regarded as a disturbance of x01.(2)
Ifε is regarded as a disturbance of x0k k=2,3,…,n−1, then we have(39)B+ΔB=B+B1Arλ00⋮⋮00ε000⋮⋮00,ΔB=00⋮⋮00−12r0ε0⋮⋮−12rn−kλn−k+rn−k−1λn−k−1ε0,Y+ΔY=Y+B2Arλ0⋮0ε0⋮0,ΔY=0⋮0r0ε⋮rn−kλn−k−rn−k−1λn−k−1ε.Similar to the discussion of the first case, we obtain(40)ΔB=ε21+∑i=1n−kri−1λi−1+riλi2≜ε21+∑i=1n−kθi,(41)ΔY≜ε1+∑i=1n−kriλi−ri−1λi−12=ε1+∑i=1n−kφi.Substituting (40) and (41) again into Lemma 1 (30), we have the following disturbance bound:(42)Lx0k=κ†γ†ε1+∑i=1n−kθi2Bx+1+∑i=1n−kφiB+κ†γ†1+∑i=1n−kθi2BrxB,whenε is regarded as a disturbance of x0k k=2,…,n−1.(3)
Now, we consider the last case, i.e.,ε is regarded as a disturbance of x0n. Similar to the above two cases, we have(43)B+ΔB=B+B1Arλ00⋮⋮00ε0,ΔB=00⋮⋮00−12r0ε0,Y+ΔY=Y+GArλ0⋮0ε,ΔY=0⋮0r0ε,(44)ΔB=ε2r02=ε2,ΔY=εr02=ε.Substituting (44) into Lemma 1 (30) again, we immediately obtain(45)Lx0n=κ†γ†ε12Bx+1B+κ†γ†12BrxB,whenε is regarded as a disturbance of x0n.
As can be seen from (32), the disturbance bounds Lx0k k=1,2,⋯,n−1 are increasing functions with respect to n since θi>0andφi>0. That is to say the larger the sample size of the original sequence is, the larger the perturbation bound of the solution will be. Therefore, similar to the traditional GM(1,1), the FGM(1,1) model, and the NIPGM(1,1) model, the new WFGM(1,1) model is also suitable for prediction problem with small samples.


## 2.1. The Weighted Fractional-Accumulation Generating Operator

In this section, a new weighted fractional-accumulation generating operator (WF-AGO) and the corresponding inverse accumulation generating operator (WF-IAGO) are defined.Definition 1 (see [34]).
LetN and N+ denote the set of natural numbers and the set of positive integers, respectively. Let r be a real number. Then, for n∈N, the general rising factorial function of a rational number is(1)rn=1,ifn=0,rr+1⋯r+n−1n!,ifn∈N+.
It has been studied in Theorem1 [41] that rn is a monotonically decreasing function and a monotonically increasing function with respect to n if r∈0,1 and r∈1,+∞, respectively. In particular, if r=1, then rn≡1. By letting r∈0,1 and introducing another weighted parameter λ∈0,1 used in the NIPGM(1,1) model [18], we define a new accumulation generating operator as follows.Definition 2.
Let(2)X0=x01,x02,…,x0nT,be the original nonnegative sequence, r,λ∈0,1 be two nonnegative parameters, and the weighted fractional-accumulation generating sequence of X0 is defined as(3)Xrλ=xrλ1,xrλ2,…,xrλnT,where(4)xrλk=∑i=1krk−iλk−ix0i,k=1,…,n.
According to Definition2, the weighted fractional-accumulation generating operator (WF-AGO) can be simply expressed by using the following matrix-vector form:(5)Xrλ=ArλX0,where Arλ is a unit lower triangular matrix and takes the form(6)Arλ=r000⋯0r1λr00⋯0r2λ2r1λr0⋯0⋮⋮⋮⋱⋮rn−1λn−1rn−2λn−2rn−3λn−3⋯r0.
Evidently, ifr=λ=1, then Arλ is a lower triangular matrix with all nonzero elements being one, which reduces to the well-known first order-accumulation generating operator (1-AGO) matrix. If r=1 and λ=1, then the weighted fractional-accumulation generating operator matrix Arλ reduces to(7)Aλ=100⋯0λ10⋯0λ2λ1⋯0⋮⋮⋮⋱⋮λn−1λn−2λn−3⋯1,Ar=r000⋯0r1r00⋯0r2r1r0⋯0⋮⋮⋮⋱⋮rn−1rn−2rn−3⋯r0,which are called the new information priority-accumulation generating operator (NIP-AGO) matrix [18] and the fractional-accumulation generating operator (F-AGO) matrix [17], respectively.Definition 3.
LetX0 and Xrλ be the original sequence and the weighted fractional-accumulation generating sequence defined as in (2) and (3), respectively. The weighted fractional-inverse accumulation generating operator (WF-IAGO) is(8)X0=Arλ−1Xrλ≜DrλXrλ.
With the notations given above, we can useDr and Dλ to denote the fractional-inverse accumulation generating operator (F-IAGO) [17] and the new information priority-inverse accumulation generating operator (NIP-IAGO) [18] matrices, which are the inverse matrices of Ar and Aλ, respectively. For the specific expression of Drλ, we have the following result.Theorem 1.
The WF-IAGO matrixDrλ satisfies(9)Drλ=−r000⋯0−r1λ−r00⋯0−r2λ2−r1λ−r0⋯0⋮⋮⋮⋱⋮−rn−1λn−1−rn−2λn−2−rn−3λn−3⋯−r0.Proof.
According to Definition1 and (6), we can see that the matrix Arλ is a unit lower triangular matrix. Therefore, detArλ=1 and its inverse matrix Drλ exits.
To prove that the matrixDrλ has form (9), we only need to show that ArλDrλ is an identity matrix. Since both Arλ and Drλ are unit lower triangular matrices, so is ArλDrλ. In addition, for i>j, we have(10)ArλDrλij=λi−jri−jλi−j−r0+ri−j−1λi−j−1−r1λ+⋯+r1λ−ri−j−1λi−j−1+r0−ri−jλi−j=λi−jArDrij.
From Theorem2 [41], we can see that ArDr is an identity matrix. Thus, for i>j we have ArλDrλij=0, which means that the matrix Drλ defined as in (9) is the inverse matrix of Arλ. Therefore, the proof is completed.


## 2.2. The Representation of the WFGM(1,1) Model

Based on the idea of the grey theory, the new WFGM(1,1) model is constructed in this section.Definition 4.
LetX0=x01,x02,…,x0nT and Xrλ=xrλ1,xrλ2,…,xrλnT be the original sequence (2) and the WF-AGO sequence (3), respectively. We call the following linear differential equation(11)dxrλtdt+axrλt=b,the whitenization equation of the WFGM(1,1) model, and the difference equation(12)xrλk−xrλk−1+azrλk=b,the basic model of the WFGM(1,1) model, where(13)zrλk=12xrλk+xrλk−1,k=2,3,…,n,is called the background value, and a and b are called the development coefficient and the grey input action, respectively.
From the definition of the WFGM(1,1) model and given the initial guessxrλtt=1=x01, we can easily obtain the solution of the whitenization equation (11):(14)x^rλk=x01−bae−ak−1+ba,k=1,2,…,n,which is called the time response function of the WFGM(1,1) model. Note that the time response function (14) is often used to approximate the values of the WF-AGO sequence Xrλ=xrλ1,…,xrλn. By using the WF-IAGO defined as in (8), the predicted values X^0=x^01,…,x^0n of the original sequence can be obtained:(15)x^0k=∑i=1k−rk−iλk−ix^rλi,k=1,2,…,n,or in matrix form, it satisfies(16)X^0=DrλX^rλ.
Obviously, the new WFGM(1,1) model is a generalization of the classical GM(1,1) model and its improvements studied recently. In fact, ifr=λ=1, then the WFGM(1,1) model reduces to the GM(1,1) model [3, 4]. If λ=1 and r=1, then the new WFGM(1,1) model reduces to the FGM(1,1) model [17] and the NIPGM(1,1) model [18], respectively. With the aid of two parameters r and λ, we can choose suitable values to not only meet with the principle of new information priority but also improve the prediction accuracy.


## 2.3. Parameters Estimation and Computational Steps

As can be seen from Section2.2, four parameters need to be determined in the WFGM(1,1) model. To obtain the WF-AGO sequence, we have to determine the parameters r and λ in advance. To get the solution of the WFGM(1,1) model, the development coefficient a and the grey input action b need to be fixed in the linear differential equation (11). Fortunately, all these parameters can be obtained by solving two optimization problems.Similar to the idea of the GM(1,1) model and its improvements, the parametersa and b can be obtained by minimizing the following unconstrained optimization problem:(17)mina,b∑k=2nxrλk−xrλk−1+azrλk−b2,which is constructed by difference equation (12). The solution of the above optimization problem can be simply written as the following linear system:(18)ab=BTB−1BTY,where(19)B=−zrλ21−zrλ31⋮⋮−zrλn1,Y=xrλ2−xrλ1xrλ3−xrλ2⋮xrλn−xrλn−1.Next, we turn to study how to determine the parametersr and λ, which play very important roles to make the WF-AGO sequence Xrλ better satisfy the grey exponential law and the principle of new information priority, and thus get better prediction accuracy. That is to say, the parameters must be chosen to get both the best fitting accuracy and the best prediction accuracy. There are many ways to define the model prediction accuracy, such as the mean absolute percentage error (MAPE), the mean absolute error (MAE), and the root mean square error (RMSE). Among these, the MAPE is one of the most common indexes to judge the merits of prediction models. Therefore, we use the MAPE as an objective function in this paper:(20)MAPEr,λ=1n∑k=1nx^0k−x0kx0k×100%.Theoretically, the MAPE is constrained by several conditions, such as the constraints on parametersr and λ and the time response function x^rλk and its inverse sequence x^0k, as well as the parameters a and b. In summary, the following constrained optimization problem can be defined to obtain the optimal parameters r and λ:(21)minMAPEr,λ,s.t0<λ≤1,0<r≤1,ab=BTB−1BTY,x^rλk=x01−bae−ak−1+ba,x^0k=∑i=1k−rk−iλk−ix^rλi,k=1,2,…,n.In general, the above optimization problem is nondifferentiable, nonlinear, and is very difficult to solve. To find its optimal solution, we need to determine the parametersa and b in advance. This is a contradiction. Fortunately, we can use the quantum particle swarm optimization (QPSO) algorithm to fast get a quasi-optimal solution of (21). Note that the QPSO algorithm is a valuable development and quality improvement of the classical PSO algorithm. The greatest advantage of the QPSO algorithm over the PSO algorithm is that the QPSO ameliorates the bug of PSO preferably and improves the search efficiency [42]. Now, it has been used to solve such nonlinear constrained optimization problems [42–44]. In the following, we present a brief introduction to the QPSO algorithm to solve the nonlinear constrained optimization problem (21). Matlab code for the QPSO algorithm can be found in Appendix A.Since there are two variables in (21), only a two-dimensional problem needs to be considered. Assume that there are m discrete particles distributed in the two-dimensional space. At the lth iteration, they have the form(22)El=E1l,E2l,…,Eml=E11lE21l⋯Em1lE12lE22l⋯Em2l,where Eil i=1,…,m denotes the position of the ith particle and Ei1l,Ei2lT denotes its coordinate. Here, Ei1l and Ei2l can be regarded as the choice for the parameters r and λ, respectively. At the initial iteration, i.e., l=1, let Pi1=Ei1 (i=1,…,m). Then, at the lth iteration, the optimal position for the ith particle (often called “pbest”) is obtained by computing(23)pbest=Pil=Eil,if MAPEEil<MAPEPil−1,Pil−1,if MAPEEil≥MAPEPil−1,i=1,…,m,and the global optimal position, which is called “gbest” and can be regarded as the global quasi-optimal solution for the nonlinear constrained optimization problem (21), is obtained:(24)gbest=Pgl=Pg1lPg2l≜r∗lλ∗l,where(25)g=argmin1≤i≤mMAPEPil.At thel+1th iteration, the position Ei1l+1,Ei2l+1T for each particle is updated by the following rule:(26)Eidl+1=Qidl±αCdl−Eidl⋅ln1uidl,uidl∼U0,1,d=1,2,where(27)Qidl=ϕdlPidl+1−ϕdlPgdl,ϕd∼U0,1,Cdl=1m∑i=1mPidl,α=0.5−0.5lmax−ll,and lmax is the given largest iteration number.According to the above discussion, the computational steps of the new WFGM(1,1) model can be summarized as follows. The corresponding Matlab code can be found in AppendixB and Appendix C.Step 1: determine the original data sequenceX0 according to the actual application background. Set the initial parameters r and λ, as well as the largest iteration number lmax for the QPSO algorithm.Step 2: compute the WF-AGO sequenceXrλ=ArλX0 and the background value zrλk according to (3) and (13), respectively.Step 3: solve the unconstrained optimization problem (17) or its equivalent form, i.e., the least squares problem (18), to obtain the parameters a and b.Step 4: substitute the parametersa and b into (14) to compute the approximated values X^rλ of the WF-AGO sequence Xrλ and compute the predicted values X^0 of the WFGM(1,1) model according to (16).Step 5: determine the local optimal parametersr and λ to minimize the MAPE (21) by using the QPSO algorithm.Step 6: repeat Steps 2–5 to obtain the global quasi-optimal (“gbest”) parameters r∗, λ∗, and the corresponding optimal fitted and predicted values X^0.


## 2.4. Stability Analysis

Stability is a very important property for a prediction model since some data may be not accurate in many engineering applications due to several reasons. For the original GM(1,1) model, Wu et al. first studied its stability by using the tools of matrix perturbation analysis. They defined a perturbation bound to show how the solution is affected by the change of the original dataX0. In this section, we extend this idea to study the stability property for the new WFGM(1,1) model.To this end, we first present a useful lemma and give the definition of the perturbation bound, which can be found in [27] and some recent works. Here and in the sequel, ⋅ denotes the 2-norm of either a vector or a matrix. rank⋅ stands for the rank of the corresponding matrix.Lemma 1 (see [27, 45]).
Suppose thatx and x^ satisfy(28)Bx−Y=min,B+ΔBx^−Y+ΔY=min,respectively, where B,ΔB∈ℂm×n with m≤n, and Y≠0, ΔY∈ℂm. Let κ†=B†⋅ΔB, γ†=1−B†⋅ΔB, and rx=Y−Bx, B† be the pseudoinverse of the matrix B. If rankB=rankB+ΔB and B†⋅ΔB<1, then(29)x^−x≤κ†γ†ΔBBx+ΔYB+κ†γ†ΔBBrxB.Definition 5 (see [27]).
Assume thatX0=x01,x02,…,x0n is the original nonnegative sequence. Let B and Y be defined as in (19) and B†⋅ΔB<1. We call(30)Lx0k=κ†γ†ΔBBx+ΔYB+κ†γ†ΔBBrxB,k=1,2,…,n,the perturbation bound when ε is a disturbance of x0k k=1,2,…,n.
Analogous to the analysis in [27], we have the following results for the new WFGM(1,1) model.Theorem 2.
Assume thatX0=x01,x02,…,x0nT is the original nonnegative sequence and ε is a disturbance of x0k k=1,2,…,n. Let B and Y be defined as (19). Denote by(31)θi=ri−1λi−1+riλi2,φi=riλi−ri−1λi−12.
Then, we have(32)Lx0k=κ†γ†ε∑i=1n−1θi2Bx+∑i=1n−1φiB+κ†γ†∑i=1n−1θi2BrxB,k=1,κ†γ†ε1+∑i=1n−kθi2Bx+1+∑i=1n−kφiB+κ†γ†1+∑i=1n−kθi2BrxB,k=2,⋯,n−1,κ†γ†εx2B+1B+κ†γ+12BrxB,k=n,and Lx0kk=1,2,…,n−1 is an increasing function with respect to n.Proof.
To deduce the perturbation bound for each data, we first simplify the expression of the matrixB and the vector Y. Let e=1,1,…,1T∈ℝn. Since Xrλ=ArλX0 and e=ArλDrλe, we have(33)B=−zrλ21−zrλ31⋮⋮−zrλn1n−1×2=−12−120⋯000−12−12⋯00⋮⋮⋮⋱⋮⋮000⋯−12−12n−1×nxrλ1−1xrλ2−1⋮⋮xrλn−1n×2≜B1ArλX0−Drλe,Y=xrλ2−xrλ1xrλ3−xrλ2⋮xrλn−xrλn−1n−1×1=−110⋯000−11⋯00⋮⋮⋮⋱⋮⋮000⋯−11n−1×nxrλ1xrλ2⋮xrλnn×1≜B2ArλX0.(1)
Ifε is regarded as a disturbance of x01, i.e., X10ε=x01+ε,x02,…,x0nT. Then, we have(34)B+ΔB=B1ArλX10ε−Drλe=B+B1Arλε000⋮⋮00,ΔB=−12r0+r1λε0−12r1λ+r2λ2ε0⋮⋮−12rn−2λn−2+rn−1λn−1ε0,Y+ΔY=B2ArλX10ε=Y+B2Arλε0⋮0,ΔY=r1λ−r0εr2λ2−r1λε⋮rn−1λn−1−rn−2λn−2ε.Thus, it satisfies(35)ΔBTΔB=ε24∑i=1n−1ri−1λi−1+riλi2000,ΔYTΔY=ε2∑i=1n−1riλi−ri−1λi−12.Based on the definition of 2-norm, from (35) we have(36)ΔB=ε2∑i=1n−1ri−1λi−1+riλi2≜ε2∑i=1n−1θi,(37)ΔY=ε∑i=1n−1riλi−ri−1λi−12≜ε∑i=1n−1φi,where the notations in (31) are used. Substituting (36) and (37) into Lemma 1 (30), we immediately have the disturbance bound(38)Lx01=κ†γ†ε∑i=1n−1θi2Bx+∑i=1n−1φiB+κ†γ†∑i=1n−1θi2BrxB,whenε is regarded as a disturbance of x01.(2)
Ifε is regarded as a disturbance of x0k k=2,3,…,n−1, then we have(39)B+ΔB=B+B1Arλ00⋮⋮00ε000⋮⋮00,ΔB=00⋮⋮00−12r0ε0⋮⋮−12rn−kλn−k+rn−k−1λn−k−1ε0,Y+ΔY=Y+B2Arλ0⋮0ε0⋮0,ΔY=0⋮0r0ε⋮rn−kλn−k−rn−k−1λn−k−1ε.Similar to the discussion of the first case, we obtain(40)ΔB=ε21+∑i=1n−kri−1λi−1+riλi2≜ε21+∑i=1n−kθi,(41)ΔY≜ε1+∑i=1n−kriλi−ri−1λi−12=ε1+∑i=1n−kφi.Substituting (40) and (41) again into Lemma 1 (30), we have the following disturbance bound:(42)Lx0k=κ†γ†ε1+∑i=1n−kθi2Bx+1+∑i=1n−kφiB+κ†γ†1+∑i=1n−kθi2BrxB,whenε is regarded as a disturbance of x0k k=2,…,n−1.(3)
Now, we consider the last case, i.e.,ε is regarded as a disturbance of x0n. Similar to the above two cases, we have(43)B+ΔB=B+B1Arλ00⋮⋮00ε0,ΔB=00⋮⋮00−12r0ε0,Y+ΔY=Y+GArλ0⋮0ε,ΔY=0⋮0r0ε,(44)ΔB=ε2r02=ε2,ΔY=εr02=ε.Substituting (44) into Lemma 1 (30) again, we immediately obtain(45)Lx0n=κ†γ†ε12Bx+1B+κ†γ†12BrxB,whenε is regarded as a disturbance of x0n.
As can be seen from (32), the disturbance bounds Lx0k k=1,2,⋯,n−1 are increasing functions with respect to n since θi>0andφi>0. That is to say the larger the sample size of the original sequence is, the larger the perturbation bound of the solution will be. Therefore, similar to the traditional GM(1,1), the FGM(1,1) model, and the NIPGM(1,1) model, the new WFGM(1,1) model is also suitable for prediction problem with small samples.


## 3. Numerical Experiments

In this section, four numerical examples, which come from the annual per capita electricity consumption prediction problem [26], the energy consumption prediction problem [18], the natural gas consumption prediction problem [13, 33], and the total output value of construction industry prediction problem, respectively, are adopted to show the high precision of the new WFGM(1,1) model. To better show its advantages, numerical results of the WFGM(1,1) model are compared with those of the original GM(1,1) model, the fractional GM(1,1) model (FGM(1,1)) [17], and the new information priority GM(1,1) model (NIPGM(1,1)) [18] studied recently. MAPE (20) is used as the evaluation index. In addition, we also compare the relative percentage error of each data for the four discussed grey models:(46)Δ^k=x^0k−x0kx0k×100%.Note that all the FGM(1,1) model, the NIPGM(1,1) model, and the new WFGM(1,1) model can be regarded as the accumulation generating operator improvements of the traditional GM(1,1) model. So, it is fair to compare their performance. Besides, the QPSO algorithm, which is derived for the WFGM(1,1) model in Section2.3, is also suitable for the FGM(1,1) model and the NIPGM(1,1) model. In actual implementations of the QPSO algorithm, we take m=5 and lmax=500, i.e., 5 discrete particles are distributed in the problem domain and the largest iteration number is 500. The detailed Matlab codes for these test problems are listed in appendices.Example 1.
Consider the annual per capita electricity consumption prediction problem in China [26]. The raw data (unit: kilowatt hour) are collected from the official website (http://data.stats.gov.cn/english/easyquery.htm?_cn=C01) of the National Bureau of Statistics of China.
For the first example, we list the raw data in Table1. Note that the official website only offers the data from 2000 to 2015 (k=1,…,16), which have been studied in [26]. Here, the data from 2000 to 2012 (k=1,…,13) are used to build the grey models and the data from 2013 to 2015 (k=14,15,16) are used to verify the prediction accuracy. By using the QPSO algorithm studied in Section 2.3, the quasi-optimal parameters for the FGM(1,1), the NIPGM(1,1), and the new WFGM(1,1) prediction models can be obtained. We list them in Table 2. The fitted data and the predicted data for each model are listed in Table 1. The relative error for each data and the error index MAPE are also listed in Table 1. To better show their comparison, we plot the raw data, the fitted data, and the predicted data for each prediction model in Figure 1. Their relative errors are also plotted.
From Table1 and Figure 1, we can see that the accumulation generating operator has great effects on both the fitting accuracy and the prediction accuracy of the grey model. The fitting accuracy of the FGM(1,1) model is slightly better than that of the NIPGM(1,1) model, while the NIPGM(1,1) model has slightly better prediction accuracy than the FGM(1,1) model. Although the fitting error of the classical GM(1,1) model is nice, its prediction error is over 10% and seems unacceptable. Both the fitting accuracy and the prediction accuracy show that the new WFGM(1,1) model is the best one. Most important, the fitting error of the WFGM(1,1) model for each data is less than 5% and the corresponding MAPE value is less than 4%, which show very high precision of the new proposed WFGM(1,1) model.Table 1
Fitted and predicted values of annual per capita electricity consumption of China (Example 1).
Year
k
Actual values
GM(1,1)
FGM(1,1)
NIPGM(1,1)
WFGM(1,1)
x
0
k
x
^
0
k
Δ
^
k
%
x
^
0
k
Δ
^
k
%
x
^
0
k
Δ
^
k
%
x
^
0
k
Δ
^
k
%
2000
1
1066.9
1066.9
0.00
1066.9
0.00
1066.9
0.00
1066.9
0.00
2001
2
1157.6
1266.0
9.36
1143.0
1.26
1088.8
5.95
1111.9
3.95
2002
3
1286.0
1399.9
8.85
1314.2
2.19
1289.4
0.27
1288.0
0.15
2003
4
1477.1
1548.0
4.80
1508.2
2.10
1496.7
1.33
1490.3
0.89
2004
5
1695.2
1711.7
0.97
1714.1
1.12
1710.8
0.92
1704.9
0.57
2005
6
1913.0
1892.8
1.06
1929.4
0.86
1932.0
0.99
1928.2
0.79
2006
7
2180.6
2093.0
4.02
2153.4
1.25
2160.4
0.93
2158.9
0.99
2007
8
2482.2
2314.4
6.76
2386.1
3.87
2396.3
3.46
2396.8
3.44
2008
9
2607.6
2559.2
1.86
2627.9
0.78
2639.9
1.24
2641.7
1.31
2009
10
2781.7
2829.9
1.73
2879.3
3.51
2891.6
3.95
2893.6
4.02
2010
11
3134.8
3129.3
0.18
3140.8
0.19
3151.5
0.53
3152.8
0.57
2011
12
3497.0
3460.3
1.05
3413.0
2.40
3420.0
2.20
3419.4
2.22
2012
13
3684.2
3826.3
3.86
3696.5
0.33
3697.3
0.35
3693.5
0.25
MAPE
3.42
1.53
1.70
1.47
2013
14
3993.0
4231.1
5.96
3992.1
0.20
3983.7
0.23
3975.4
0.44
2014
15
4132.9
4678.7
13.21
4300.3
4.05
4279.5
3.55
4265.3
3.20
2015
16
4321.0
5173.6
22.28
4621.8
9.24
4585.0
8.37
4563.4
7.86
MAPE
13.82
4.44
4.05
3.83Table 2
Parametersr and λ used in grey prediction models for Example 1.
GM(1,1)
FGM(1,1)
NIPGM(1,1)
WFGM(1,1)
r
1
0.2401
1
0.3337
λ
1
1
0.2074
0.7307Figure 1
Numerical results for Example 1: fitted and predicted values (a) and relative errors (b).
(a)
(b)Example 2.
Consider the energy consumption prediction problem for Jiangsu Province, China. The raw data (unit: ten thousand tons standard coal) come from ([18], Table 2) and have been studied to show high precision of the NIPGM(1,1) model.
For the second example, the raw data, which show the energy consumption of Jiangsu Province, China from 2001 to 2012 (k=1,…,12), are listed in Table 3. Here, in order to verify the accuracy of the WFGM(1,1), the data of energy consumption from 2001 to 2008 (k=1,…,8) are used to fit and the data from 2009 to 2012 (k=9,…,12) are used to predict. The quasi-optimal parameters of the FGM(1,1), the NIPGM(1,1) and the new WFGM(1,1) prediction models for the second example are listed in Table 4. The fitted data and the predicted data of the four discussed prediction models are listed in Table 3 and are plotted in Figure 2. In Figure 2, we also plot the relative error of each data for the four discussed prediction models. Both the fitting accuracy and the prediction accuracy show that all the FGM(1,1) model, the NIPGM(1,1) model, and the WFGM(1,1) model greatly improve the classical GM(1,1) model. From Table 3 and Figure 2, we can see again that the new WFGM(1,1) model has both the best fitting accuracy and the best prediction accuracy, which confirms our original intention of this work. More specifically, from the aspect of the fitting accuracy, the NIPGM(1,1) model and the proposed WFGM(1,1) model have almost the same result and both are better than the FGM(1,1) model. From the aspect of the prediction accuracy, the proposed WFGM(1,1) model is better than the NIPGM(1,1) model and the later one is better than the FGM(1,1) model.Table 3
Fitted and predicted values of energy consumption of Jiangsu Province, China (Example 2).
Year
k
Actual values
GM(1,1)
FGM(1,1)
NIPGM(1,1)
WFGM(1,1)
x
0
k
x
^
0
k
Δ
^
k
%
x
^
0
k
Δ
^
k
%
x
^
0
k
Δ
^
k
%
x
^
0
k
Δ
^
k
%
2001
1
8881
8881
0.00
8881
0.00
8881
0.00
8881
0.00
2002
2
9609
10483
9.09
9998
4.05
9609
0.00
9621
1.26
2003
3
11061
11988
8.38
11990
8.40
11826
6.92
11850
7.14
2004
4
13652
13710
0.42
14048
2.90
14031
2.78
14087
3.18
2005
5
17167
15678
8.67
16129
6.05
16224
5.50
16276
5.19
2006
6
18742
17930
4.33
18239
2.69
18404
1.80
18411
1.76
2007
7
20948
20505
2.12
20391
2.66
20572
1.80
20493
2.17
2008
8
22232
23449
5.47
22599
1.65
22728
2.23
22522
1.31
MAPE
4.81
3.55
2.63
2.61
2009
9
23709
26816
13.12
24875
4.92
24872
4.90
24499
3.33
2010
10
25774
30667
18.99
27228
5.64
27004
4.77
26426
2.53
2011
11
27589
35071
27.12
29669
7.54
29124
5.56
28304
2.59
2012
12
28850
40108
39.02
32207
11.64
31232
8.26
30134
4.45
MAPE
24.56
7.43
5.87
3.23Table 4
Parametersr and λ used in grey prediction models for Example 2.
GM(1,1)
FGM(1,1)
NIPGM(1,1)
WFGM(1,1)
r
1
0.4481
1
0.7499
λ
1
1
0.6562
0.5325Figure 2
Numerical results for Example 2: fitted and predicted values (a) and relative errors (b).
(a)
(b)Example 3.
Consider the China’s natural gas consumption prediction problem (unit: billion cubic meters). This problem has been deeply studied in [13] and was used to be a test problem in [33] to show the efficiency of the new information priority discrete grey model.
The raw data of the third example are listed in Table5, which can be found in Table 1 [33]. These data show China’s natural gas consumption from 2003 to 2013 (k=1,…,11). Here, we use the data from 2003 to 2009 (k=1,…,7) as the sample and the data from 2010 to 2013 (k=8,…,11) to test the performance of different prediction models. In Table 6, we list the quasi-optimal parameters of the FGM(1,1), the NIPGM(1,1), and the new WFGM(1,1) prediction models computed by the QPSO algorithm for the third example. In Table 5, we list the numerical results for the GM(1,1), the FGM(1,1), the NIPGM(1,1), and the new WFGM(1,1) prediction models. To better show their difference, we plot the fitted data and the predicted data for each prediction model in Figure 3. The relative errors for the fitted data and the predicted data are also plotted. From these numerical results, we can see again that the new WFGM(1,1) model has both the best fitting accuracy and the best prediction accuracy compared with the classical GM(1,1) model, the FGM(1,1) model, and the recent developed NIPGM(1,1) model. Different from the second example, the FGM(1,1) model has better prediction accuracy than the NIPGM(1,1) model. This shows that for some problems the FGM(1,1) model has advantages over the NIPGM(1,1) model. Our new model can take full use of advantages of both models. In addition, the fitting error and prediction error for each data are less than 5%, which shows that the new proposed WFGM(1,1) model is a very powerful prediction model.Table 5
Fitted and predicted values of natural gas consumption of China (Example 3).
Year
k
Actual values
GM(1,1)
FGM(1,1)
NIPGM(1,1)
WFGM(1,1)
x
0
k
x
^
0
k
Δ
^
k
%
x
^
0
k
Δ
^
k
%
x
^
0
k
Δ
^
k
%
x
^
0
k
Δ
^
k
%
2003
1
35
35.0
0.00
35.0
0.00
35.0
0.00
35.0
0.00
2004
2
41.5
43.7
5.22
41.6
0.31
41.5
0.00
41.5
0.00
2005
3
49.3
50.3
2.03
50.6
2.70
50.3
2.05
50.6
2.70
2006
4
58.6
57.9
1.12
59.6
1.71
59.3
1.14
59.8
2.00
2007
5
69.2
66.7
3.55
68.5
1.06
68.4
1.19
68.7
0.74
2008
6
80.3
76.9
4.25
77.3
3.74
77.6
3.12
77.3
3.69
2009
7
85.2
88.6
3.95
86.2
1.13
87.1
2.17
85.8
0.60
MAPE
2.87
1.52
1.41
1.39
2010
8
94.8
102.2
7.61
95.1
0.35
96.6
1.92
93.8
1.03
2011
9
103.1
117.5
13.98
104.2
1.10
106.3
3.15
101.7
1.38
2012
10
107.2
135.4
26.28
113.5
5.90
116.2
8.43
109.2
1.94
2013
11
119.3
155.9
30.71
123.0
3.12
126.3
5.86
116.6
2.23
MAPE
19.65
2.62
4.84
1.65Table 6
Parametersr and λ used in grey prediction models for Example 3.
GM(1,1)
FGM(1,1)
NIPGM(1,1)
WFGM(1,1)
r
1
0.4318
1
0.3842
λ
1
1
0.7464
0.4446Figure 3
Numerical results for Example 3: fitted and predicted values (a) and relative errors (b).
(a)
(b)Example 4.
Consider the China’s total output value of construction industry prediction problem. The raw data (unit: ¥100 million) are collected from 2008 to 2018 and can be found in the official website (http://data.stats.gov.cn/easyquery.htm?cn=C01) of the National Bureau of Statistics of China.
For the fourth example, the raw data are listed in Table7. The data from 2009 to 2014 (k=1,…,6) are used as the sample and the data from 2015 to 2019 (k=7,…,10) are used to test the performance of four different prediction models. The quasi-optimal parameters of the FGM(1,1), the NIPGM(1,1) model, and the WFGM(1,1) models are computed by the QPSO algorithm and listed in Table 8. In Table 7, the fitted data, the predicted data, and the corresponding values of error indexes for each model are listed. Figure 4 clearly shows these numerical results. From Table 7 and Figure 4, we can see that all prediction models can effectively simulate the sample data. Although the GM(1,1) model has the worst fitness MAPE, its value is only 1.97%. The fitness MAPE values of other three models are almost the same and only a quarter of the GM(1,1) model. However, the predicted data again show that the original GM(1,1) model is unacceptable. All other improved GM(1,1) models can effectively simulate the near further trend and greatly reduce the prediction errors. Among them, the WFGM(1,1) model preforms the best. This further confirms that the new proposed WFGM(1,1) model is a valuable improvement of the grey models.Table 7
Fitted and predicted values of China’s total output value of construction industry (Example 4).
Year
k
Actual values
GM(1,1)
FGM(1,1)
NIPGM(1,1)
WFGM(1,1)
x
0
k
x
^
0
k
Δ
^
k
%
x
^
0
k
Δ
^
k
%
x
^
0
k
Δ
^
k
%
x
^
0
k
Δ
^
k
%
2009
1
76807.70
76807.70
0.00
76807.70
0.00
76807.70
0.00
76807.70
0.00
2010
2
96031.10
99920.93
4.05
96031.28
0.00
96030.98
0.00
96031.21
0.00
2011
3
116463.30
115815.32
0.56
117392.32
0.80
117091.24
0.54
117352.95
0.76
2012
4
137217.86
134238.02
2.17
138102.74
0.64
137746.77
0.39
138112.61
0.65
2013
5
160366.06
155591.22
2.98
157967.27
1.50
158005.36
1.47
157985.41
1.48
2014
6
176713.40
180341.07
2.05
177039.71
0.18
177874.64
0.66
176929.72
0.12
MAPE
1.97
0.52
0.51
0.50
2015
7
180757.47
209027.86
15.64
195402.49
8.10
197362.09
9.19
194964.79
7.86
2016
8
193566.78
242277.86
25.16
213131.48
10.11
216475.04
11.83
212125.83
9.59
2017
9
213943.56
280816.92
31.26
230290.99
7.64
235220.70
9.95
228451.99
6.78
2018
10
235085.53
325486.37
38.45
246934.64
5.04
253606.11
7.88
243982.59
3.78
MAPE
27.63
7.72
9.71
7.00Table 8
Parametersr and λ used in grey prediction models for Example 4.
GM(1,1)
FGM(1,1)
NIPGM(1,1)
WFGM(1,1)
r
1
0.2282
1
0.2061
λ
1
1
0.6888
0.4505Figure 4
Numerical results for Example 4: fitted and predicted values (a) and relative errors (b).
(a)
(b)


## 4. Concluding Remarks and Outlooks

In this paper, we have proposed a novel weighted fractional GM(1,1) model, which is a valuable development and a quality improvement of the classical GM(1,1) model from the aspect of the accumulation generating process. The new model involves some existing improved GM(1,1) models as special cases. We have derived the modeling process of the new model and studied its stability property. In addition, the QPSO algorithm is used to obtain two optimal parameters involved in the new accumulation generating operator. With four numerical examples from different applications, we have shown that the new WFGM(1,1) prediction model has much better fitting accuracy and prediction accuracy than the classical GM(1,1) model and two improved GM(1,1) models studied recently, i.e., the FGM(1,1) model [17] and the NIPGM(1,1) model [18]. Future works should focus on testing the WFGM(1,1) model on other applications where data exhibit strong nonlinear or periodicity characteristics, constructing new weighted fractional self-adaptive nonlinear and nonhomogeneous grey models, error estimates, and so on.



---

*Source: 6570683-2020-01-20.xml*