# Research on the Stability of NDGM Model with the Fractional Order Accumulation and Its Optimization

Authors: Huiming Duan; Kailiang Shao; Xinping Xiao; Jinwei Yang
Journal: Mathematical Problems in Engineering
 (2017)

---

## Abstract

The grey forecasting model has been successfully applied in numerous fields since it was proposed. The nonhomogeneous discrete grey model (NDGM) was approximately constructed based on the nonhomogeneous index trend; it increased the applicability of discrete grey model. However, the NDGM required accurate data and better effect when the original data did not meet the conditions and fitting and prediction errors were larger. For this, the NDGM with the fractional order accumulating operator (abbreviated asNDGMp/q) has higher performance. In this paper, the matrix perturbation bound of the parameters was used to analyze the stability of NDGMp/q and the NDGMp/q can decrease the disturbance bound. Subsequently, the parameter estimation method of NDGMp/q was studied and the Particle Swarm Optimization algorithm was employed to optimize the order number of NDGMp/q and some steps were provided. In addition, the results of two practical examples demonstrated that the perturbation of NDGMp/q was smaller than that of NDGM and provided remarkable predication performance compared with the traditional NDGM model and DGM model.


---

## Body

## 1. Introduction

Forecasting the future values of time series data plays a very important role in our research; thus, many forecasting methods have been developed for many years, such as the Rough sets theory proposed by Pawlak (see [1, 2]) and fuzzy mathematics proposed by Zadeh (see [3]). However, because of limited knowledge and information, only part of system structure could be fully known. To address this problem, professor Deng proposed grey forecasting models to catch the system development tendency [4, 5].As a core model of grey prediction, GM (1, 1) model has been widely used in some fields such as transportation, agriculture, economy, and management [6–10]. Meanwhile, many scholars have improved the GM (1, 1) model a lot [11–19], thus enhancing its simulative accuracy and predictive accuracy. However, during the processes of both practical application and theoretical research of the GM (1, 1) model, GM (1, 1) model directly jumped from discrete form to continuous form, which resulted in failing in completely fitting homogeneous exponential sequence in simulation and prediction. Then the discrete grey model [20] was put forward to solve the transformation from discrete to continuous GM (1, 1) model. However, GM (1, 1) model and DGM model were constructed based on a hypothesis that the original data sequence was a homogeneous index sequence. However, the fact was not consistent, and the most original data sequence was the nonhomogeneous index sequence.Hence, Xie and Liu [21] come up with nonhomogeneous discrete grey model (NDGM), and the model was constructed based on the approximate nonhomogeneous index trend. The results indicate there is no error between original value and simulative value based on pure nonhomogeneous index sequence. And NDGM model was the extension of DGM model while the latter was the special case of the former. NDGM model increased the applicability of discrete grey model. According to improve simulation and prediction accuracies, there were some the results [22, 23] of nonhomogeneous discrete grey model. However, these models had a higher requirement for data; when the data did not meet the requirements, the errors of both model-fitting and prediction were larger. The actual value cannot always meet the definition of monotonic increasing (decreasing) concave (convex) function. So the original data were accumulated to increase the exponent law according to the classical grey modeling mechanism. However, the accumulation of integer order was not suitable for some data, as modeling, and the effects of simulation and prediction were poor. Nevertheless, the current fraction grey model which used factional order accumulation [24–27] has the important significance for improving the performance of grey model.Therefore, Wu et al. [28] defined the actual data with fractional order accumulation and defined the NDGM model with fractional order accumulation (NDGMp/q). But they only presented the algorithm of model and some random fractional order values without considering the initial value of the model and presenting the nature of NDGMp/q model. When a system was studied, the stability must be considered. They did not study the stability of the NDGMp/q model; hence, the parameters of matrix perturbation bound were used to analyze the stability of NDGMp/q model; then it was concluded that the solution of the NDGMp/q model perturbation bound of solution was smaller than that of NDGM model. When 0<p/q<1, the solution of NDGMp/q model perturbation bound of solution was smaller and could decrease the disturbance bound. And they did not conclude which order number was the best as well as how to obtain the best order number and other properties. In order to solve these issues, the parameter estimation method of NDGMp/q was studied. Furthermore, the Particle Swarm Optimization algorithm was employed to optimize the order number of the NDGMp/q model and obtained the better simulation and prediction results. Finally, the results from previous works demonstrated that the perturbation bound solution of that NDGMp/q model was smaller than traditional NDGM model in Case 1, and the results the Case 2 practical demonstrated that NDGMp/q provided better predication performance than the traditional NDGM model and DGM model.This paper is organized as follows. In Section2, we introduced the NDGM model with integer order accumulating operator and the NDGMp/q model. In Section 3, the stability of the NDGMp/q model was discussed. In Section 4, the Particle Swarm Optimization algorithm was employed to optimize the order number of NDGMp/q model; some steps were provided. Finally, in Section 5, the paper was concluded.


## 2. NDGM Model with Fractional Order Accumulating Operator

Definition 1 (see [29]).
Assume that the sequence(1)X0=x01,x02,…,x0nis an original data sequence, and the symbol D represents a kind of mathematical operational method. When D is applied once on the sequence X0, we have that(2)X0D1=X1=x11,x12,…,x1nis the accumulated generation sequence of X0, where (3)x1k=∑k=1nx0i,k=1,2,…,n.Then D is called the first-order accumulating generation operator of X0, denoted by 1-AGO. If D is applied r times on X0, we obtain (4)X0Dr=Xr=xr1,xr2,…,xrn,where (5)xrk=∑k=1nxr-1i,k=1,2,…,n,r∈Z+.
ThenD is called the r-order accumulating generation operator of X0, denoted by r-AGO.
Accordingly, the inverse accumulating generation operator is the inverse operation of the accumulating generation process and plays a role in recovery from the acts of accumulating operators; the inverse accumulating operator is defined as follows:(6)X0Dr=αrX0=αrx01,αrx02,…,αrx0n,where (7)αrx0k=αr-1x0k-αr-1x0k-1,k=2,3,…,n.In Definition 1, the order r∈Z+, and we call the two operators accumulating and inverse accumulating generation operators with integer order.
The sequence(8)Z1=z12,z13,…,z1nis the mean sequence of X1, where (9)z1k=12x1k+x1k-1,k=2,3,…,n.Definition 2 (see [29]).
The equation(10)x0k+az1k=b,x10=x01is called a GM1,1 model, where z1k=1/2x1k+x1k-1, k=2,3,…,n. The whitenization equation dx1/dt+ax1=b of GM1,1 model is solved to obtain (11)x^1t+1=x1-bae-at+ba.Definition 3 (see [20]).
The sequencesX0 and X1 are defined as (1) and (2). Then the equation (12)x1k+1=β1x1k+β2is called discrete grey model (DGM). Similar with GM (1, 1) model, β1 and β2 are the parameters of DGM model. The equation (13)x^1k+1=β1kx01-β21-β1+β21-β1,k=1,2,…,n-1,is called the recursive function of DGM model.Definition 4 (see [21]).
The sequencesX0 and X1 are defined as (1) and (2). Then the equation (14)x^1k+1=β1x^1k+β2k+β3x^11=x11+β4is called nonhomogenous discrete grey model (NDGM). x^1k is the simulative value of x1k and x11 is the iterative value of the NDGM model. β1, β2, β3, and β4 are parameters of NDGM model. The equation (15)x^1k+1=x^11β1k+β2∑j=1kjβ1k-j-1-β1k1-β1×β3is called the recursive function of DGM model.Definition 5.
Let thep/q0<p/q<1 order accumulated generating operator of the original nonhomogeneous index sequence X0 be Xp/q. Let Cp/q0=1, Ck-1k=0; then xp/qk=∑i=1kCk-i+p/q-1k-ix0i‍, k=1,2,…,n, where (16)Ck-i+p/q-1k-i=p/q+k-i-1p/q+k-i-2⋯r+1p/qk-i!.p/q order inverse accumulated generating operator of X0 is(17)αp/qX0=α1X1-p/q=α1X1-p/q1,α1X1-p/q2,…,α1X1-p/qn.Definition 6 (see [28]).
Assume that(18)x^p/qk+1=β1x^p/qk+β2k+β3x^p/q1=x11+β4is established.
The sequencesX0 and Xp/q are defined as (2) and (18). Then (18) is called nonhomogenous discrete grey model with fractional order accumulation (abbreviated as NDGMp/q model). x^p/qk is the simulative value of x1k and x^p/q1 is the iterative value of the NDGM model. β1, β2, β3, and β4 are parameters of NDGMp/q model.That is least square method. So we can get the expressions of parameters in Proposition7.Proposition 7.
Based on the least square method the first level parametersβ1, β2, and β3 satisfy the matrix equation(19)β1β2β3=BTB-1BTY,where (20)B=xp/q111xp/q221⋮⋮⋮xp/qn-1n-11,Y=xp/q2xp/q3⋮xp/qn.Proposition 8.
The recursive function ofNDGMp/q model is (21)x^p/qk+1=β1kx^p/q1+β2∑j=1kjβ1k-j+1-β1k1-β1β3,k=1,2,…,n-1.Proof.
By Definition6, we have (22)x^p/qk+1=β1x^p/qk+β2k+β3=β1β1x^p/qk-1+β2k-1+β3+β2k+β3=β12x^p/qk-1+β2β1k-1+k+β31+β1=β13x^p/qk-2+β2β12k-2+β1k-1+k+β31+β1+β12=⋯=β1kx^p/q1+β2∑j=1kjβ1k-j+1-β1k1-β1β3,k=1,2,…,n-1.
We use the least square method to calculate the value of parameterβ4. By minimizing the error of x^p/qk and xp/qk, minβ4=∑k=1n[x^(1)(k)-x(1)(k)]. We can construct a nonrestraint optimized model; then (23)β4=x^p/qk+1-β1kx11-β2∑j=1kjβ1k-j-β31-β1k/1-β1β1k1+∑k=1n-1β1k2,where k=1,2,…,n-1.


## 3. The Stability of the Solution of NDGM Model with Fractional Order Accumulation

### 3.1. The Stability ofNDGMp/q Model

In order to illustrate the stability of theNDGMp/q model, the theorem related to the matrix perturbation analysis was introduced as follows.Lemma 9 (see [29]).
LetA∈Cn×n, b∈Cn, M=A+E, and c=b+k. Suppose x+h and x satisfy Mx-c2=min and Ax-b2=min. If rankA=rank(M) and A†2E2<1, where A† is the pseudo-inverse of matrix A, then (24)h≤κ†γ†E2Ax+kA+κ†γ†E2ArxA,where A, E, and rx are the tolerance norm and rx=b-Ax, κ†=A†2A, and γ†=1-A†2E2.Theorem 10.
For theNDGMp/q model of original data x01,x02,…,x0n based on nonhomogeneous index sequence, if the rth data is disturbed, that is, x^0r=x0r+ε, r=1,2,…,n-1, based on the least square method BX-Y=min, let x be a solution of the Lemma 9; then the perturbation bound is L1x0r and (25)L1x01=εκ†γ†∑k=1n-1Ck+p/q-2k-12Bx+∑k=2nCk+p/q-2k-12B+κ†γ†∑k=1n-1Ck+p/q-2k-12BrxB,L1x0r=εκ†γ†∑k=1n-rCk+p/q-2k-12Bx+∑k=1n-r+1Ck+p/q-2k-12B+κ†γ†∑k=1n-rCk+p/q-2k-12BrxB,where r=2,3,…,n-1, and when x^0n=x0n+ε, then (26)L1x0n=κ†γ†εB,where mathematical notation is the same as Lemma 9.Proof.
Ifε is regarded as a disturbance of x01, then (27)B^=B+ΔB=xp/q111xp/q221⋮⋮⋮xp/qn-1n-11+ε00pqε00⋮⋮⋮Cn-3+p/qn-2ε00,Y^=Y+ΔYxp/q2xp/q3⋮xp/qn+pqεC1+p/q2ε⋮Cn-2+p/qn-1ε;
thus, we have(28)ΔBTΔB=1+pq2+Cp/q+122+⋯+Cn-2+p/qn-12ε200000000.Because 1+p/q2+Cp/q+122+⋯+Cn-2+p/qn-12ε2 is ΔBTΔB maximum of eigenvalue, then (29)ΔB2=λmaxΔBTΔB=∑k=1n-1Ck+p/q-2k-12ε,ΔY2=pq2+Cp/q+122+⋯+Cn-2+p/qn-12ε=∑k=2nCk+p/q-2k-12ε.
By Lemma9 and based on the least square method Bx-Y2=min, Δx is the disturbed solution of equation Bx=Y, and we have (30)Δx≤κ†γ†ΔB2Bx+ΔYB+κ†γ†ΔB2BrxB,L1x01=εκ†γ†∑k=1n-1Ck+p/q-2k-12Bx+∑k=2nCk+p/q-2k-12B+κ†γ†∑k=1n-1Ck+p/q-2k-12BrxB.
Whenx^02=x02+ε, then (31)ΔB=000ε00pqε00⋮⋮⋮Cn-4+p/qn-3ε00,ΔY=εpqεC1+p/q2ε⋮Cn-3+p/qn-2ε.
According to the proof ofL2x01, we have (32)L1x02=εκ†γ†∑k=1n-2Ck+p/q-2k-12Bx+∑k=1n-1Ck+p/q-2k-12B+κ†γ†∑k=1n-2Ck+p/q-2k-12BrxB.
According to this rule, ifε is regarded as a disturbance of x0r, r=3,4,…,n-1, both ΔB and ΔY changed; thus, (33)L1x0r=εκ†γ†∑k=1n-rCk+p/q-2k-12Bx+∑k=1n-r+1Ck+p/q-2k-12B+κ†γ†∑k=1n-rCk+p/q-2k-12BrxB.
Ifε is regarded as a disturbance of x0n, we have (34)ΔBTΔB=000000000,Y=00⋮ε,L1x0n=κ†γ†εB.
Whenp/q=1, let L2(x(0)(r) be the perturbation bound of the NDGM model. By Theorem 10, we have (35)L2x0r=εκ†γ†n-rxB1+n-r+1B1+κ†γ†n-rrxB1,r=1,2,…,n-1,L2x0n=κ†γ†εB1.Obviously,κ†/γ† and ΔB are all positive and increasing function of r; other variables are not changed. Therefore, perturbation bounds L1x0r and L2x0r are the increasing function of the effect of sample size n; that is, L1x0r and L2x0r will change large when n→+∞. When r is a constant value, p/q is lager; then L1x0r is lager too. When 0<p/q≤1, we have L1x0r<L2x0r, r=1,2,…,n-1. When the perturbation bound of solution was large, the perturbation was not always large. This was because the perturbation cannot exceed the perturbation bound. Hence, from the perspective of the size of perturbation bound, when x^0r=x0+ε, r=1,2,…,n-1, the disturbed solution of the NDGMp/q model is lower. The NDGMp/q model is more stable than NDGM model.


### 3.2. Verification of the StabilityNDGMp/q Model

In this section, the effectiveness of theNDGMp/q model is verified by two real cases study. Mean absolute percentage error (MAPE) compares the real and forecasted values to evaluate the precision. MAPE is defined as (MAPE=100%1/n∑k=1n|x(k)-x^(k)/x(k)|), where x0k is the actual value at time k, and x^0k is the predicted value for time k.Case 1 (syphilis incidence predicted in China [30]).
The trends of syphilis incidence in China forecasting example [31] compare the precision; the historical incidence of China from 2000 to 2008 was employed as the model-fitting. Then the actual values of 2009 and 2010 were predicted. The NDGMp/q model, NDGM model, and a classical GM(1, 1) model were built, respectively, to simulate the trends of syphilis incidence in China. The simulated/forecasted values and absolute percentage error with different grey models were shown in Table 1.Table 1
The fitted values and MAPE of four grey models (unit: millionth).
Year
Actual value
GM(1, 1)
NDGM
N
D
G
M
0.5
N
D
G
M
0.01
2000
5.08
5.08
5.08
5.08
5.08
2001
4.80
3.55
3.65
3.74
4.23
2002
4.67
4.52
4.61
4.41
4.43
2003
4.50
5.76
5.84
5.71
5.50
2004
7.12
7.34
7.42
7.47
7.29
2005
9.67
9.34
9.45
9.70
9.67
2006
12.80
11.90
12.04
12.42
12.54
2007
15.88
15.15
15.37
15.66
15.84
2008
19.49
19.30
19.65
19.51
19.48
MAPE
10.01
8.93
8.01
5.49
2009
23.07
24.57
25.13
24.02
23.42
2010
26.86
31.29
32.16
29.28
27.59
MAPE
11.50
14.32
6.55
2.11As can be seen from Table1, from a short-term forecasting viewpoint, NDGM0.01 obtained lower MAPE than the GM(1, 1) model, which implied that the NDGM0.01 could significantly enhance the precision of grey forecasting model. Moreover, NDGM0.01 obtained lower MAPE than the NDGM0.5; therefore, it was also shown that the perturbation bound was smaller and the prediction accuracy was higher. According to Table 1, a scatter broken-line figure regarding the simulated and forecasted data was shown in Figure 1.Figure 1
Simulates and forecasts obtained using different grey models.Figure1 showed that the whole performances of the four models for simulating and forecasting the trends of syphilis incidence in China. Obviously, the prediction effect of GM(1, 1) model was better than that of NDGM model, but the simulation and prediction precisions of other NDGM models were better than those of GM(1, 1) model. For NDGM model, the simulation and prediction precisions of NDGM0.01 were the best, followed by NDGM0.5, but those of NDGM model of integer order were the worst. Moreover, for the NDGM model of fractional order, the simulation, and prediction precisions of smaller fraction were higher than those of bigger one. In addition, Theorem 10 can also be proved to be correct from the stability of system. Namely, the smaller the perturbation order was, the higher the system stability was, and the higher the precisions of simulation and prediction were.


## 3.1. The Stability ofNDGMp/q Model

In order to illustrate the stability of theNDGMp/q model, the theorem related to the matrix perturbation analysis was introduced as follows.Lemma 9 (see [29]).
LetA∈Cn×n, b∈Cn, M=A+E, and c=b+k. Suppose x+h and x satisfy Mx-c2=min and Ax-b2=min. If rankA=rank(M) and A†2E2<1, where A† is the pseudo-inverse of matrix A, then (24)h≤κ†γ†E2Ax+kA+κ†γ†E2ArxA,where A, E, and rx are the tolerance norm and rx=b-Ax, κ†=A†2A, and γ†=1-A†2E2.Theorem 10.
For theNDGMp/q model of original data x01,x02,…,x0n based on nonhomogeneous index sequence, if the rth data is disturbed, that is, x^0r=x0r+ε, r=1,2,…,n-1, based on the least square method BX-Y=min, let x be a solution of the Lemma 9; then the perturbation bound is L1x0r and (25)L1x01=εκ†γ†∑k=1n-1Ck+p/q-2k-12Bx+∑k=2nCk+p/q-2k-12B+κ†γ†∑k=1n-1Ck+p/q-2k-12BrxB,L1x0r=εκ†γ†∑k=1n-rCk+p/q-2k-12Bx+∑k=1n-r+1Ck+p/q-2k-12B+κ†γ†∑k=1n-rCk+p/q-2k-12BrxB,where r=2,3,…,n-1, and when x^0n=x0n+ε, then (26)L1x0n=κ†γ†εB,where mathematical notation is the same as Lemma 9.Proof.
Ifε is regarded as a disturbance of x01, then (27)B^=B+ΔB=xp/q111xp/q221⋮⋮⋮xp/qn-1n-11+ε00pqε00⋮⋮⋮Cn-3+p/qn-2ε00,Y^=Y+ΔYxp/q2xp/q3⋮xp/qn+pqεC1+p/q2ε⋮Cn-2+p/qn-1ε;
thus, we have(28)ΔBTΔB=1+pq2+Cp/q+122+⋯+Cn-2+p/qn-12ε200000000.Because 1+p/q2+Cp/q+122+⋯+Cn-2+p/qn-12ε2 is ΔBTΔB maximum of eigenvalue, then (29)ΔB2=λmaxΔBTΔB=∑k=1n-1Ck+p/q-2k-12ε,ΔY2=pq2+Cp/q+122+⋯+Cn-2+p/qn-12ε=∑k=2nCk+p/q-2k-12ε.
By Lemma9 and based on the least square method Bx-Y2=min, Δx is the disturbed solution of equation Bx=Y, and we have (30)Δx≤κ†γ†ΔB2Bx+ΔYB+κ†γ†ΔB2BrxB,L1x01=εκ†γ†∑k=1n-1Ck+p/q-2k-12Bx+∑k=2nCk+p/q-2k-12B+κ†γ†∑k=1n-1Ck+p/q-2k-12BrxB.
Whenx^02=x02+ε, then (31)ΔB=000ε00pqε00⋮⋮⋮Cn-4+p/qn-3ε00,ΔY=εpqεC1+p/q2ε⋮Cn-3+p/qn-2ε.
According to the proof ofL2x01, we have (32)L1x02=εκ†γ†∑k=1n-2Ck+p/q-2k-12Bx+∑k=1n-1Ck+p/q-2k-12B+κ†γ†∑k=1n-2Ck+p/q-2k-12BrxB.
According to this rule, ifε is regarded as a disturbance of x0r, r=3,4,…,n-1, both ΔB and ΔY changed; thus, (33)L1x0r=εκ†γ†∑k=1n-rCk+p/q-2k-12Bx+∑k=1n-r+1Ck+p/q-2k-12B+κ†γ†∑k=1n-rCk+p/q-2k-12BrxB.
Ifε is regarded as a disturbance of x0n, we have (34)ΔBTΔB=000000000,Y=00⋮ε,L1x0n=κ†γ†εB.
Whenp/q=1, let L2(x(0)(r) be the perturbation bound of the NDGM model. By Theorem 10, we have (35)L2x0r=εκ†γ†n-rxB1+n-r+1B1+κ†γ†n-rrxB1,r=1,2,…,n-1,L2x0n=κ†γ†εB1.Obviously,κ†/γ† and ΔB are all positive and increasing function of r; other variables are not changed. Therefore, perturbation bounds L1x0r and L2x0r are the increasing function of the effect of sample size n; that is, L1x0r and L2x0r will change large when n→+∞. When r is a constant value, p/q is lager; then L1x0r is lager too. When 0<p/q≤1, we have L1x0r<L2x0r, r=1,2,…,n-1. When the perturbation bound of solution was large, the perturbation was not always large. This was because the perturbation cannot exceed the perturbation bound. Hence, from the perspective of the size of perturbation bound, when x^0r=x0+ε, r=1,2,…,n-1, the disturbed solution of the NDGMp/q model is lower. The NDGMp/q model is more stable than NDGM model.


## 3.2. Verification of the StabilityNDGMp/q Model

In this section, the effectiveness of theNDGMp/q model is verified by two real cases study. Mean absolute percentage error (MAPE) compares the real and forecasted values to evaluate the precision. MAPE is defined as (MAPE=100%1/n∑k=1n|x(k)-x^(k)/x(k)|), where x0k is the actual value at time k, and x^0k is the predicted value for time k.Case 1 (syphilis incidence predicted in China [30]).
The trends of syphilis incidence in China forecasting example [31] compare the precision; the historical incidence of China from 2000 to 2008 was employed as the model-fitting. Then the actual values of 2009 and 2010 were predicted. The NDGMp/q model, NDGM model, and a classical GM(1, 1) model were built, respectively, to simulate the trends of syphilis incidence in China. The simulated/forecasted values and absolute percentage error with different grey models were shown in Table 1.Table 1
The fitted values and MAPE of four grey models (unit: millionth).
Year
Actual value
GM(1, 1)
NDGM
N
D
G
M
0.5
N
D
G
M
0.01
2000
5.08
5.08
5.08
5.08
5.08
2001
4.80
3.55
3.65
3.74
4.23
2002
4.67
4.52
4.61
4.41
4.43
2003
4.50
5.76
5.84
5.71
5.50
2004
7.12
7.34
7.42
7.47
7.29
2005
9.67
9.34
9.45
9.70
9.67
2006
12.80
11.90
12.04
12.42
12.54
2007
15.88
15.15
15.37
15.66
15.84
2008
19.49
19.30
19.65
19.51
19.48
MAPE
10.01
8.93
8.01
5.49
2009
23.07
24.57
25.13
24.02
23.42
2010
26.86
31.29
32.16
29.28
27.59
MAPE
11.50
14.32
6.55
2.11As can be seen from Table1, from a short-term forecasting viewpoint, NDGM0.01 obtained lower MAPE than the GM(1, 1) model, which implied that the NDGM0.01 could significantly enhance the precision of grey forecasting model. Moreover, NDGM0.01 obtained lower MAPE than the NDGM0.5; therefore, it was also shown that the perturbation bound was smaller and the prediction accuracy was higher. According to Table 1, a scatter broken-line figure regarding the simulated and forecasted data was shown in Figure 1.Figure 1
Simulates and forecasts obtained using different grey models.Figure1 showed that the whole performances of the four models for simulating and forecasting the trends of syphilis incidence in China. Obviously, the prediction effect of GM(1, 1) model was better than that of NDGM model, but the simulation and prediction precisions of other NDGM models were better than those of GM(1, 1) model. For NDGM model, the simulation and prediction precisions of NDGM0.01 were the best, followed by NDGM0.5, but those of NDGM model of integer order were the worst. Moreover, for the NDGM model of fractional order, the simulation, and prediction precisions of smaller fraction were higher than those of bigger one. In addition, Theorem 10 can also be proved to be correct from the stability of system. Namely, the smaller the perturbation order was, the higher the system stability was, and the higher the precisions of simulation and prediction were.


## 4. Optimizing the Order Number of theNDGMp/q Model by Particle Swarm Optimization (PSO) Algorithm

### 4.1.NDGMp/q Model Steps

We have discussed theNDGMp/q model and the main steps of modeling; the steps were as follows.Step 1.
Calculation of the orderp/q accumulation sequence (36)Xp/q=xp/q1,xp/q2,…,xp/qn.Step 2.
Putxp/qk, k=1,2,…,n into (19); then use the least squares estimate parameter of β1, β2, and β3.Step 3.
Using(37)x^p/qk+1=β1kx^p/q1+β2∑j=1kjβ1k-j+1-β1k1-β1β3,k=1,2,…,n-1,then we will predict xp/q1,x^p/q2,….Step 4.
Assume thatXp/q=x^p/q1,x^p/q2,…,x^p/qn, is the fractional order accumulation sequence and the equation(38)αp/qX0=α1x^1-p/q1,α1x^1-p/q2,…,α1x^1-p/qn+1,…is established.


### 4.2. Optimizing the Order Number of theNDGMp/q Model by Particle Swarm Optimization (PSO) Algorithm

The value of the mean absolute percentage error (MAPE) is often used to judge the merits of modeling. The order of theNDGMp/q model is sought under the condition of the least mean absolute percent as follows: (39)minfr=1n-1∑k=2nx0k-x^0kx0k,r∈R+.In this section, we will use the Particle Swarm Optimization (PSO) algorithm to optimize the order number of the NDGMp/q model. The searing process of the optimization order of the NDGMp/q model is as follows.Step 1.
Defining structure and initializing randomly the position (pBest) and speed for each particle, letpBest=1.Step 2.
Set the firstpBest=1 as the current position and gBest=1 as the optimal particle position in initial swarm.Step 3.
All the particles in the particle swarm were operated according to the following order:(1)
Update the locations and speeds of particles; letC1=1, C2=1, and (40)P·V=V+C1×rand×pBest-Present+C2×rand×gBest-PresentPresent=Present+V.(2)
Judge whether the particle location Present meets the scope; if beyond the scope, a new location should be set for Present. If within the scope, then the calculation continued according to the following steps. Whenr=pBest, the average relative error of fractional order operator model, the concrete steps were as follows:
(a)
ComputingXr, Xr is the rth order accumulating generating sequence of X0.(b)
Compute the mean generated sequence with consecutive neighboursZr of Xr.(c)
Compute ther-order inverse accumulating generation sequence X-r of Xr.(d)
Solve the parametersa^=β1,β2,β3.(e)
Deduce the time response expression ofx^rk.(f)
Compute the simulative valuesx^rk.(g)
Compute the restored simulative valuesx^0k of x^rk.(h)
Compute the MAPEf(pBest).(i)
Judge whether the value offpBest is less than f(gBest); if fpBest<fgBest, the new position is set to pBest. If the f(pBest) is not superior to gBest, the new position is set to gBest.Step 4.
Step3 was carried out according to the order from the second particle to the last one.Step 5.
Judging whether the algorithm meets the convergence rule, if it meets then go to Step6; else go to Step 3.Step 6.
Output Best “g,” which is the optimal value of the order outputting the simulated or forecasted values of the NDGMp/q model, when r=Best.According to the modeling mechanism of the proposedNDGMp/q model and the Particle Swarm Optimization (PSO) algorithm of the order number, the flowchart of the new model is summarized in Figure 2.Figure 2
The flowchart of theNDGMp/q model.


### 4.3. Verification ofNDGMp/q Model

Case 2 (logistics demand forecasting in Jiangsu province [30]).
We consider an example from paper [30] which provides the sample data. We will build four models that simulate Jiangsu province’s logistics demand from 2002 to 2006 and to forecast it from 2007 to 2009. We will compare the simulation and prediction mean absolute percentage error of the NDGMp/q model with a classical NDGM model and DGM model.
According to MATLAB programs of theNDGMp/q model and Particle Swarm Optimization (PSO) algorithm, the optimal order of the NDGMp/q is p/q=0.3298. The simulated/forecasted values and mean absolute percentage error with the four different grey models are shown in Table 2.Table 2
The fitted values and MAPE of four grey models (unit:108 ton-km).
Year
Actual value
DGM
NDGM
NDGM0.3298
2002
1549.12
1549.12
1549.12
1549.12
2003
1817.44
1899.10
1807.62
1811.13
2004
2398.13
2374.15
2426.34
2418.20
2005
3068.30
2968.18
3041.35
3048.45
2006
3644.14
3710.85
3652.70
3648.96
MAPE
2.65
0.71
0.49
2007
4098.42
4639.35
4260.39
4212.31
2008
4707.50
5800.16
4864.46
4741.86
2009
5154.46
7251.42
5464.90
5242.83
MAPE
25.70
4.44
1.74According to Table2, in order to show the performances of four models clearly, we can draw Figures 3(a), 3(b), and 3(c).Figure 3
(a)
TheNDGMp/q model
(b)
The NDGM model
(c)
The DGM modelAs can be seen from Figures3(a)–3(c), obviously, both of the simulation and prediction performance of the grey model with the optimal fractional order were superior to those of the corresponding grey model with an integer order. This showed that the fractional order accumulating generation could really improve the performance of a grey prediction model. It was shown that the improved effect of the fractional order on grey prediction model was also influenced by the raw model. Moreover, we can rank the effect of three models from the worst to the best in terms of their abilities to fit the data, namely, DGM model, NDGM model, and NDGMp/q model.


## 4.1.NDGMp/q Model Steps

We have discussed theNDGMp/q model and the main steps of modeling; the steps were as follows.Step 1.
Calculation of the orderp/q accumulation sequence (36)Xp/q=xp/q1,xp/q2,…,xp/qn.Step 2.
Putxp/qk, k=1,2,…,n into (19); then use the least squares estimate parameter of β1, β2, and β3.Step 3.
Using(37)x^p/qk+1=β1kx^p/q1+β2∑j=1kjβ1k-j+1-β1k1-β1β3,k=1,2,…,n-1,then we will predict xp/q1,x^p/q2,….Step 4.
Assume thatXp/q=x^p/q1,x^p/q2,…,x^p/qn, is the fractional order accumulation sequence and the equation(38)αp/qX0=α1x^1-p/q1,α1x^1-p/q2,…,α1x^1-p/qn+1,…is established.


## 4.2. Optimizing the Order Number of theNDGMp/q Model by Particle Swarm Optimization (PSO) Algorithm

The value of the mean absolute percentage error (MAPE) is often used to judge the merits of modeling. The order of theNDGMp/q model is sought under the condition of the least mean absolute percent as follows: (39)minfr=1n-1∑k=2nx0k-x^0kx0k,r∈R+.In this section, we will use the Particle Swarm Optimization (PSO) algorithm to optimize the order number of the NDGMp/q model. The searing process of the optimization order of the NDGMp/q model is as follows.Step 1.
Defining structure and initializing randomly the position (pBest) and speed for each particle, letpBest=1.Step 2.
Set the firstpBest=1 as the current position and gBest=1 as the optimal particle position in initial swarm.Step 3.
All the particles in the particle swarm were operated according to the following order:(1)
Update the locations and speeds of particles; letC1=1, C2=1, and (40)P·V=V+C1×rand×pBest-Present+C2×rand×gBest-PresentPresent=Present+V.(2)
Judge whether the particle location Present meets the scope; if beyond the scope, a new location should be set for Present. If within the scope, then the calculation continued according to the following steps. Whenr=pBest, the average relative error of fractional order operator model, the concrete steps were as follows:
(a)
ComputingXr, Xr is the rth order accumulating generating sequence of X0.(b)
Compute the mean generated sequence with consecutive neighboursZr of Xr.(c)
Compute ther-order inverse accumulating generation sequence X-r of Xr.(d)
Solve the parametersa^=β1,β2,β3.(e)
Deduce the time response expression ofx^rk.(f)
Compute the simulative valuesx^rk.(g)
Compute the restored simulative valuesx^0k of x^rk.(h)
Compute the MAPEf(pBest).(i)
Judge whether the value offpBest is less than f(gBest); if fpBest<fgBest, the new position is set to pBest. If the f(pBest) is not superior to gBest, the new position is set to gBest.Step 4.
Step3 was carried out according to the order from the second particle to the last one.Step 5.
Judging whether the algorithm meets the convergence rule, if it meets then go to Step6; else go to Step 3.Step 6.
Output Best “g,” which is the optimal value of the order outputting the simulated or forecasted values of the NDGMp/q model, when r=Best.According to the modeling mechanism of the proposedNDGMp/q model and the Particle Swarm Optimization (PSO) algorithm of the order number, the flowchart of the new model is summarized in Figure 2.Figure 2
The flowchart of theNDGMp/q model.


## 4.3. Verification ofNDGMp/q Model

Case 2 (logistics demand forecasting in Jiangsu province [30]).
We consider an example from paper [30] which provides the sample data. We will build four models that simulate Jiangsu province’s logistics demand from 2002 to 2006 and to forecast it from 2007 to 2009. We will compare the simulation and prediction mean absolute percentage error of the NDGMp/q model with a classical NDGM model and DGM model.
According to MATLAB programs of theNDGMp/q model and Particle Swarm Optimization (PSO) algorithm, the optimal order of the NDGMp/q is p/q=0.3298. The simulated/forecasted values and mean absolute percentage error with the four different grey models are shown in Table 2.Table 2
The fitted values and MAPE of four grey models (unit:108 ton-km).
Year
Actual value
DGM
NDGM
NDGM0.3298
2002
1549.12
1549.12
1549.12
1549.12
2003
1817.44
1899.10
1807.62
1811.13
2004
2398.13
2374.15
2426.34
2418.20
2005
3068.30
2968.18
3041.35
3048.45
2006
3644.14
3710.85
3652.70
3648.96
MAPE
2.65
0.71
0.49
2007
4098.42
4639.35
4260.39
4212.31
2008
4707.50
5800.16
4864.46
4741.86
2009
5154.46
7251.42
5464.90
5242.83
MAPE
25.70
4.44
1.74According to Table2, in order to show the performances of four models clearly, we can draw Figures 3(a), 3(b), and 3(c).Figure 3
(a)
TheNDGMp/q model
(b)
The NDGM model
(c)
The DGM modelAs can be seen from Figures3(a)–3(c), obviously, both of the simulation and prediction performance of the grey model with the optimal fractional order were superior to those of the corresponding grey model with an integer order. This showed that the fractional order accumulating generation could really improve the performance of a grey prediction model. It was shown that the improved effect of the fractional order on grey prediction model was also influenced by the raw model. Moreover, we can rank the effect of three models from the worst to the best in terms of their abilities to fit the data, namely, DGM model, NDGM model, and NDGMp/q model.


## 5. Conclusion

When some systems were studied, the stability must be considered. The grey theory model mainly deals with the incomplete information cases and the uncertainty of the system models. From the perspective of model stability,NDGMp/q model is more stable than traditional NDGM, and it was the result of a case to explain the stability. Because GM(1, 1) model, DGM(1, 1) model, and NDGM model had a higher requirement for data and when the data did not meet the requirements, the errors of both simulation and prediction were larger, and the NDGMp/q model has the important significance for improving the performance of grey model. Furthermore, we provided the Particle Swarm Optimization algorithm which was employed to optimize the order number of NDGMp/q model. Using a case studied in Section 3, it has been shown that the MAPE of the NDGMp/q model was always the lowest among the four grey models.



---

*Source: 9728587-2017-06-27.xml*