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.

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 [36].

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 by 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 into a monotonicity increased sequence . And then use the solution of a first order-ordinary differential equation, which is called the whitenization differential equation, to approximate the sequence . Finally, applying the first order-inverse accumulated generating operator (1-IAGO), one can get the fitted value of 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 [1419], 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, 2426] for examples. Modeling mechanism analysis can be found in [2729]. 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 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, 3941].

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 [4244]. 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 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. 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]). Let N and denote the set of natural numbers and the set of positive integers, respectively. Let r be a real number. Then, for , the general rising factorial function of a rational number isIt has been studied in Theorem 1 [41] that is a monotonically decreasing function and a monotonically increasing function with respect to n if and , respectively. In particular, if , then . By letting and introducing another weighted parameter used in the NIPGM(1,1) model [18], we define a new accumulation generating operator as follows.

Definition 2. Letbe the original nonnegative sequence, be two nonnegative parameters, and the weighted fractional-accumulation generating sequence of is defined aswhereAccording to Definition 2, the weighted fractional-accumulation generating operator (WF-AGO) can be simply expressed by using the following matrix-vector form:where is a unit lower triangular matrix and takes the formEvidently, if , then 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 and , then the weighted fractional-accumulation generating operator matrix reduces towhich 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 and 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) isWith the notations given above, we can use and 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 and , respectively. For the specific expression of , we have the following result.

Theorem 1. The WF-IAGO matrix satisfies

Proof. According to Definition 1 and (6), we can see that the matrix is a unit lower triangular matrix. Therefore, and its inverse matrix exits.
To prove that the matrix has form (9), we only need to show that is an identity matrix. Since both and are unit lower triangular matrices, so is . In addition, for , we haveFrom Theorem 2 [41], we can see that is an identity matrix. Thus, for we have , which means that the matrix defined as in (9) is the inverse matrix of . 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. Let and be the original sequence (2) and the WF-AGO sequence (3), respectively. We call the following linear differential equationthe whitenization equation of the WFGM(1,1) model, and the difference equationthe basic model of the WFGM(1,1) model, whereis 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 , 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 . By using the WF-IAGO defined as in (8), the predicted values of the original sequence can be obtained:or in matrix form, it satisfiesObviously, the new WFGM(1,1) model is a generalization of the classical GM(1,1) model and its improvements studied recently. In fact, if , then the WFGM(1,1) model reduces to the GM(1,1) model [3, 4]. If and , 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 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). The 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 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:

Theoretically, the MAPE is constrained by several conditions, such as the constraints on parameters r and λ and the time response function and its inverse sequence , 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. 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 [4244]. 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 formwhere denotes the position of the ith particle and denotes its coordinate. Here, and can be regarded as the choice for the parameters r and λ, respectively. At the initial iteration, i.e., , let (). Then, at the lth iteration, the optimal position for the ith particle (often called “”) is obtained by computingand the global optimal position, which is called “” and can be regarded as the global quasi-optimal solution for the nonlinear constrained optimization problem (21), is obtained:where

At the th iteration, the position for each particle is updated by the following rule:whereand 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 Appendix B and Appendix C.Step 1: determine the original data sequence according to the actual application background. Set the initial parameters r and λ, as well as the largest iteration number for the QPSO algorithm.Step 2: compute the WF-AGO sequence and the background value 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 parameters a and b into (14) to compute the approximated values of the WF-AGO sequence and compute the predicted values 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.Step 6: repeat Steps 2–5 to obtain the global quasi-optimal (“”) parameters , , and the corresponding optimal fitted and predicted values .

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 data . 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. stands for the rank of the corresponding matrix.

Lemma 1 (see [27, 45]). Suppose that x and satisfyrespectively, where with , and , . Let , , and , be the pseudoinverse of the matrix B. If and , then

Definition 5 (see [27]). Assume that is the original nonnegative sequence. Let B and Y be defined as in (19) and . We callthe perturbation bound when ε is a disturbance of .
Analogous to the analysis in [27], we have the following results for the new WFGM(1,1) model.

Theorem 2. Assume that is the original nonnegative sequence and ε is a disturbance of . Let B and Y be defined as (19). Denote byThen, we haveand is an increasing function with respect to n.

Proof. To deduce the perturbation bound for each data, we first simplify the expression of the matrix B and the vector Y. Let . Since and , we have(1)If ε is regarded as a disturbance of , i.e., . Then, we haveThus, it satisfiesBased on the definition of 2-norm, from (35) we havewhere the notations in (31) are used. Substituting (36) and (37) into Lemma 1 (30), we immediately have the disturbance boundwhen ε is regarded as a disturbance of .(2)If ε is regarded as a disturbance of , then we haveSimilar to the discussion of the first case, we obtainSubstituting (40) and (41) again into Lemma 1 (30), we have the following disturbance bound:when ε is regarded as a disturbance of .(3)Now, we consider the last case, i.e., ε is regarded as a disturbance of . Similar to the above two cases, we haveSubstituting (44) into Lemma 1 (30) again, we immediately obtainwhen ε is regarded as a disturbance of .As can be seen from (32), the disturbance bounds are increasing functions with respect to n since . 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:

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 and , 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 Table 1. Note that the official website only offers the data from 2000 to 2015 (), which have been studied in [26]. Here, the data from 2000 to 2012 () are used to build the grey models and the data from 2013 to 2015 () 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 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. 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 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 and the corresponding MAPE value is less than , 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. 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 (), 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 () are used to fit and the data from 2009 to 2012 () 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.

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 Table 5, which can be found in Table 1 [33]. These data show China’s natural gas consumption from 2003 to 2013 (). Here, we use the data from 2003 to 2009 () as the sample and the data from 2010 to 2013 () 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.

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 Table 7. The data from 2009 to 2014 () are used as the sample and the data from 2015 to 2019 () 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 . 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.

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.

Appendix

A. Matlab code for the QPSO algorithm

(1)function [gbestfit,lambda, r]=(2)WFGM_QPSO(maxgen,popsize, dimensions,populmax1, populmin1, populmax2, populmin2, x0, t, n(3) Input: maxgen-largest iteration number(4)    popsize-number of discrete particles(5)    dimensions-number of parameters(6)    populmax1-upper bound of the parameter lambda(7)    populmin1-lower bound of the parameter lambda(8)    populmax2-upper bound of the parameter r(9)    populmin2-lower bound of the parameter r(10)    x0-data sequence(11)    t-data length used for fitting(12)    n-total length of the input data(13) Output: gbestfit-best fitting error(14)    lambda-quasi-optimal weight(15)    r-quasi-optimal fractional order(16)(17) Setting bounds for parameters(18)populmax = [populmax1, populmax2];(19)populmin = [populmin1, populmin2];(20)(21) Initializing discrete particles(22)popul1 = rand(1, popsize)(populmax(1)-populmin(1)) + populmin(1);(23)popul2 = rand(1, popsize)(populmax(2)-populmin(2)) + populmin(2);(24)popul = [popul1’, popul2’]’;(25)(26) Fitting the data sequence with initial discrete particles(27)for m = 1:popsize(28)  xNFGM = WFGM(x0(1 : t), n, t, popul(1, m), popul(2, m));(29)  fitness(m) = sum((abs(xNFGM(1 : t)-x0(1 : t))./x0(1 : t))100)/t;(30)end(31)(32) Initializing parameters(33)ibestpos = popul; % Initializing individual best position(34)ibestfit = fitness; % Fittness of each individual(35)[fbestpart, g] = min(ibestfit); % The best global fit(36)gbestfit = fbestpart;(37)gbestpos = ibestpos(:,g); % The optimal parameters correspond to the best global fitness(38)(39)for count = 1:maxgen(40) for m = 1:popsize(41)  xNFGM = WFGM(x0(1 : t), n, t, popul(1, m), popul(2, m));(42)  fitness(m) = sum((abs(xNFGM(1 : t)-x0(1 : t))./x0(1 : t))100)/t;(43)  if fitness(m)  ibestfit(m);(44)   ibestfit(m) = fitness(m);(45)   ibestpos(:,m) = popul(:,m);(46)  end(47) end(48) mbest = sum(ibestpos)/popsize;(49) [fbestpart,g] = min(fitness); % Update global history best position(50) if fbestpart  gbestfit(51)   gbestfit = fbestpart;(52)   gbestpos = popul(:,g);(53) end(54) bestfitness(count) = gbestfit;(55)(56) update the position(57)W(count) = 0.5–0.5(maxgen-count)/maxgen;(58)for j = 1 : popsize(59)  for i = 1 : dimensions(60)   f = rand;(61)   pp = fibestpos(i, j)+(1-f)gbestpos(i)’;(62)   u = rand;(63)   if u 0.5(64)    popul(i, j) = pp-W(count)abs(mbest(1, j)-popul(i, j))log(1/u);(65)   else(66)    popul(i, j) = pp + W(count)abs(mbest(1, j)-popul(i, j))log(1/u);(67)   end(68)   if popul(i, j) populmax(i)(69)    popul(i,  j) = rand(populmax(i)-populmin(i)) + populmin(i);  %rand1.5; populmax(70)   elseif popul(i, j) populmin(i)(71)    popul(i,  j) = rand(populmax(i)-populmin(i)) + populmin(i);  %rand(−1.5); populmin(72)   end(73)  end(74) end(75)end(76)(77)lambda = gbestpos(1);(78)r = gbestpos(2);(79)(80)return

B. Matlab code for the WFGM(1,1) prediction model

(1)function x0_WFGM = WFGM(x0, n, t, lambda, r)(2) Input: x0-data sequence(3)   n-total length of the input data(4)   t-data length used for fitting(5)   lambda-weight(6)   r-fractional order(7) Output: x0WFGM - predicted data(8)(9) Construct weighted fractional acumulation matrix A(10)d(1) = 1;(11)for i = 2 : n(12)  d(i) = d(i − 1)lambda(r + i − 2)/(i − 1);(13)end(14)for i = 1 : n(15)  A(i : n, i) = d(1:n − i + 1);(16)end(17)(18) Data preprocess(19)x1 = A(1 : t,1 : t)x0’;(20) Compute background values(21)z1 = 0.5(x1(2 : t)+x1(1 : t − 1));(22) Compute the development coefficient and the grey input action(23)B = [-z1 ones(t − 1,1)];(24)Y = x1(2 : t)-x1(1 : t − 1);(25)p = (BB)(BY);(26)a = p(1); b = p(2);(27) Compute the constant for whitenization equation based on initial condition(28)c = exp(a)(x0(1)-b/a);(29)(30) Compute the predicted value(31)hatx1 = cexp(-a.(1 : n))+b/a;(32)x0WFGM = Ahatx1’;(33)x0WFGM = x0WFGM’;(34)(35)return

C. Matlab code for test problems

(1)clc;(2)clear;(3)format short(4)(5)Example = input(’Choose Exmaple: ’);(6)if Example = = 1(7)   x0 = [1066.9 1157.6 1286 1477.1 1695.2 1913 2180.6 2482.2(8)   2607.6 2781.7 3134.8 3497.0 3684.2 3993 4132.9 4231];(9)   t = 13;(10)   % Annual per capita electricity consumption prediction problem. Data comes from Table 3 of(11)   % “B.-L. Wei, N.-M. Xie and A. Hu, Optimal solution for novel grey polynomial prediction model,(12)   % Applied Mathematical Modeling, vol. 62, pp. 717–727, 2018”(13)elseif Example = = 2(14)   x0 = [8881 9609 11061 13652 17167 18742 20948 22232 23709 25774 27589 28850];(15)   t = 8;(16)   % Energy consumption prediction problem. Data comes from Table 1 of “W.-J. Zhou, H.-R. Zhang,(17)   % Y.-G. Dang and Z.-X. Wang, New information priority accumulated grey discrete model and(18)   % its application, China Journal of Management Science, vol. 25, no. 8, pp. 140–148, 2017.”(19)elseif Example = = 3(20)   x0 = [35 41.5 49.3 58.6 69.2 80.3 85.2 94.8 103.1 107.2 119.3];(21)   t = 7;(22)   % Natural gas output prediction problem. Data comes from Table 4 of “B. Zeng,(23)   % Forecasting the relation of supply and demand of natural gas in China during 2015–2020 using(24)   % a novel grey model, Journal of Intelligence and Fuzzy System, vol. 32, no. 1, pp. 141–155, 2017.”(25)elseif Example = = 4(26)   x0 = [76807.70 96031.10 116463.30 137217.86 160366.06 176713.40(27)   180757.47 193566.78 213943.56 235085.53];(28)   t = 6;(29)   % Total output value of construction industry prediction problem.(30)   % Data comes from the official website http://data.stats.gov.cn/easyquery.htm?cn=C01(31)End(32)(33)n = length(x0); % Compute the length of the data sequence(34)(35)fprintf(’Fitting and predicting by the WFGM(1,1) modeln’)(36)populmax1 = 1; populmin1 = 0;(37)populmax2 = 1; populmin2 = 0;(38)dimensions = 2;(39)maxgen = 500;(40)popsize = 5;(41)[gbestfitWFGM,lambdaWFGM, rWFGM] = (42)WFGM_QPSO(maxgen,popsize, dimensions,populmax1, populmin1, populmax2, populmin2, x0, t, n);(43)fprintf(’Quasi-optimal parameters for WFGM(1,1) are: lambda %1.4f, r %1.4fn’, lambdaWFGM,rWFGM);(44)xWFGM = WFGM(x0(1 : t),  n,  t, lambdaWFGM,rWFGM);(45)epsilonFitnessxWFGM = (abs(xWFGM(1 : t)-x0(1 : t))./x0(1 : t))100;(46)MAPEFitnessxWFGM = sum(epsilonFitnessxWFGM)/t;(47)fprintf(’Fitting MAPE for WFGM(1,1) is %1.4fn’, MAPEFitnessxWFGM)(48)epsilonPredictxWFGM = (abs(xWFGM(t + 1 : n)-x0(t + 1 : n))./x0(t + 1 : n))100;(49)MAPEPredictxWFGM = sum(epsilonPredictxWFGM)/(n-t);(50)fprintf(’Prediction MAPE for WFGM(1,1) is %1.4fn’, MAPEPredictxWFGM)

Data Availability

The data used to support the findings of this study are available from the corresponding author upon request.

Conflicts of Interest

The authors declare that they have no conflicts of interest.

Acknowledgments

This work was supported by the National Natural Science Foundation of China (nos. 61771265 and 11771225), the Natural Science Foundation of the Jiangsu Higher Education Institutions of China (nos. 18KJB580012 and 19KJB580003), and the Science and Technology Project of Nantong City (nos. CP12017001, GY12017006, and JC2018142).

Supplementary Materials

Matlab code for four test problems are illustrated in Section 3; WFGM: Matlab code for the weighted fractional GM(1,1) prediction model; WFGM_QPSO: Matlab code for the QPSO algorithm used in the WFGM(1,1) model. (Supplementary Materials)