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

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 rst 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 tting accuracy. Finally, four numerical examples from dierent practical applications are present. Numerical results show that the new proposed prediction model is very ecient and has both the best tting 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.


Introduction
Given some historic data and the current data, tting 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 e ective prediction models has been proposed in the past few decades. Among numerous existing prediction models, the grey theorybased models, which were rst proposed by Deng in 1982 [2], have attracted many researchers' attention due to its simple form, high computing e ciency, and high prediction accuracy especially for limited data and insu cient information [3][4][5][6]. e classical GM(1,1) model, in which the rst one denotes the rst order and the second one denotes the single variable, is the foundation of grey theory. Denote by X (0) the original sequence. e basic idea of the GM(1,1) model is rst to use the rst order-accumulated generating operator (1-AGO) to transform the original sequence X (0) into a monotonicity increased sequence X (1) . And then use the solution X (1) of a rst orderordinary di erential equation, which is called the whitenization di erential equation, to approximate the sequence X (1) . Finally, applying the rst order-inverse accumulated generating operator (1-IAGO), one can get the tted value X (0) of X (0) 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 tra c 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 di erent aspects, such as nding new accumulation generating operators [14][15][16][17][18][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][25][26] for examples. Modeling mechanism analysis can be found in [27][28][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]. en, 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. is 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. ese 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 CO 2 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][40][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).
e 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. en, 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][43][44].
eoretically, 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. e remainder of this paper is organized as follows. In Section 2, 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.
eoretical 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.1.
e 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]). Let N and N + denote the set of natural numbers and the set of positive integers, respectively. Let r be a real number. en, for n ∈ N, the general rising factorial function of a rational number is

Definition 2. Let
be the original nonnegative sequence, r, λ ∈ (0, 1] be two nonnegative parameters, and the weighted fractional-accumulation generating sequence of X (0) is defined as where According to Definition 2, the weighted fractionalaccumulation generating operator (WF-AGO) can be simply expressed by using the following matrix-vector form: where A rλ is a unit lower triangular matrix and takes the form Evidently, if r � λ � 1, then A rλ 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 A rλ reduces to 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. Let X (0) and X (rλ) be the original sequence and the weighted fractional-accumulation generating sequence defined as in (2) and (3), respectively. e weighted fractional-inverse accumulation generating operator (WF-IAGO) is With the notations given above, we can use D r and D λ to denote the fractional-inverse accumulation generating operator (F-IAGO) [17] and the new information priorityinverse accumulation generating operator (NIP-IAGO) [18] matrices, which are the inverse matrices of A r and A λ , respectively. For the specific expression of D rλ , we have the following result.

Complexity
Proof. According to Definition 1 and (6), we can see that the matrix A rλ is a unit lower triangular matrix. erefore, de t(A rλ ) � 1 and its inverse matrix D rλ exits.
To prove that the matrix D rλ has form (9), we only need to show that A rλ D rλ is an identity matrix. Since both A rλ and D rλ are unit lower triangular matrices, so is A rλ D rλ . In addition, for i > j, we have From eorem 2 [41], we can see that A r D r is an identity matrix. us, for i > j we have (A rλ D rλ ) ij � 0, which means that the matrix D rλ defined as in (9) is the inverse matrix of A rλ . erefore, the proof is completed.  (2) and the WF-AGO sequence (3), respectively. We call the following linear differential equation the whitenization equation of the WFGM(1,1) model, and the difference equation the basic model of the WFGM(1,1) model, where 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 guess x (rλ) (t) |t�1 � x (0) (1), we can easily obtain the solution of the whitenization equation (11): 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 X (rλ) � x (rλ) (1), . . . , x (rλ) (n) . By using the WF-IAGO defined as in (8), the predicted values X (0) � x (0) (1), . . . , x (0) (n) of the original sequence can be obtained: (15) or in matrix form, it satisfies Obviously, the new WFGM(1,1) model is a generalization of the classical GM(1,1) model and its improvements studied recently. In fact, if r � λ � 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.

Parameters Estimation and Computational
Steps. As can be seen from Section 2.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 parameters a and b can be obtained by minimizing the following unconstrained optimization problem: which is constructed by difference equation (12). e solution of the above optimization problem can be simply written as the following linear system: where Next, we turn to study how to determine the parameters r and λ, which play very important roles to make the WF-AGO sequence X (rλ) better satisfy the grey exponential law and the principle of new information priority, and thus get better prediction accuracy. at is to say, the parameters must be chosen to get both the best fitting accuracy and the best prediction accuracy. ere 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. erefore, we use the MAPE as an objective function in this paper: eoretically, the MAPE is constrained by several conditions, such as the constraints on parameters r and λ and the time response function x (rλ) (k) and its inverse sequence x (0) (k), 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 λ: 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 parameters a and b in advance. is 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. e 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][43][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 where E i (l) (i � 1, . . . , m) denotes the position of the ith particle and (E i1 (l), E i2 (l)) T denotes its coordinate. Here, E i1 (l) and E i2 (l) can be regarded as the choice for the parameters r and λ, respectively. At the initial iteration, i.e., m). en, at the lth iteration, the optimal position for the ith particle (often called "pbest") is obtained by computing 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: where At the (l + 1)th iteration, the position (E i1 (l + 1), E i2 (l + 1)) T for each particle is updated by the following rule: where and l max 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. e corresponding Matlab code can be found in Appendix B and Appendix C.
Step 1: determine the original data sequence X (0) according to the actual application background. Set the initial parameters r and λ, as well as the largest iteration number l max for the QPSO algorithm.
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 parameters a and b into (14) to compute the approximated values X (rλ) of the WF-AGO sequence X (rλ) and compute the predicted values X (0) of the WFGM(1,1) model according to (16).
Step 5: determine the local optimal parameters r and λ to minimize the MAPE (21) by using the QPSO algorithm.

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. ey defined a perturbation bound to show how the solution is affected by the change of the original data X (0) . 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.
Analogous to the analysis in [27], we have the following results for the new WFGM(1,1) model.
Proof. To deduce the perturbation bound for each data, we first simplify the expression of the matrix B and the vector Y.
As can be seen from (32), the disturbance bounds L[x (0) (k)] (k � 1, 2, · · · , n − 1) are increasing functions with respect to n since θ i > 0 and φ i > 0. at is to say the larger the sample size of the original sequence is, the larger the perturbation bound of the solution will be. erefore, 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.

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: 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 Section 2.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 l max � 500, i.e., 5 discrete particles are distributed in the problem domain and the largest iteration number is 500. e 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]. e 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 Table 1. 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 quasioptimal 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. e fitted data and the predicted data for each model are listed in Table 1. e 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. eir relative errors are also plotted.
From Table 1 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. e 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    Complexity classical GM(1,1) model is nice, its prediction error is over 10% and seems unacceptable. Both the tting accuracy and the prediction accuracy show that the new WFGM(1,1) model is the best one. Most important, the tting 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.

Example 2.
Consider the energy consumption prediction problem for Jiangsu Province, China. e 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 t and the data from 2009 to 2012 (k 9, . . . , 12) are used to predict. e 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. e tted 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 tting 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 tting accuracy and the best prediction accuracy, which con rms our original intention of this work. More specically, from the aspect of the tting 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.

Example 3.
Consider the China's natural gas consumption prediction problem (unit: billion cubic meters). is problem has been deeply studied in [13] and was used to be a test problem in [33] to show the e ciency of the new information priority discrete grey model. e raw data of the third example are listed in Table 5, which can be found in Table 1 [33].
ese 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 di erent 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 di erence, we plot the tted data and the predicted data for each prediction model in Figure 3. e relative errors for the tted 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 tting 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. Di erent from the second example, the FGM(1,1) model has better prediction accuracy than the NIPGM(1,1) model. is 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 tting 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 7. e 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 di erent prediction models. e 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 tted data, the predicted data, and the corresponding values of error indexes for each model are listed.  Table 7 and Figure 4, we can see that all prediction models can e ectively simulate the sample data. Although the GM(1,1) model has the worst tness MAPE, its value is only 1.97%. e tness 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 e ectively simulate the near further trend and greatly reduce the prediction errors. Among them, the WFGM(1,1) model preforms the best. is further con rms that the new proposed WFGM(1,1) model is a valuable improvement of the grey models.

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. e 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 di erent applications, we have shown that the new WFGM(1,1) prediction model has much better tting 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.