# Parameter Optimization on the Three-Parameter Whitenization Grey Model and Its Application in Simulation and Prediction of Gross Enrollment Rate of Higher Education in China

Authors: Jihong Sun; Hui Li; Bo Zeng; Xiaoyun Zhao; Chuanhui Wang
Journal: Complexity
 (2020)

---

## Abstract

The gray prediction model, based on the GM(1,1) method, is an important branch of gray theory with the most active research and the most fruitful results, and it is the most widely used because of its small sample size, simple modeling process, and easy to use. Such advantages have been successfully applied in many fields such as transportation, agriculture, energy, medicine, and environment and have been gradually developed into a mainstream predictive modeling method. This study combines the Three-parameter Whitenization Grey Model (TWGM(1,1)), which fits the inhomogeneous exponential law sequence, and the Particle Swarm Algorithm (PSA) to optimize the order and background value coefficients under the condition of the minimum sum of squares of simulation errors, and hence, to solve the problem that the cumulative order is fixed to “1” and the background value coefficient value is fixed to “0.5.” As a result, a parameter-optimized gray system model with flexibility, adaptability, and dynamic adjustment is designed to simulate and predict China’s higher education gross enrollment rate. The application shows that the model has better overall simulation and prediction performance than others. On the one hand, the parametric optimization model significantly improves its own performance, and on the other hand, its intelligent and adjustable adaptivity improves the accuracy and further extends its application.


---

## Body

## 1. Introduction

Since Professor Ju-Long [1] established the gray system theory in 1982, the gray prediction model based on the GM(1,1) method is the most active, fruitful, and widely used branch of gray theory. With feature of the small sample size, simple modeling process, and low cost of learning, it has been successfully applied to many fields, such as transportation, agriculture, energy, medicine, and environment [2–6]. The model gradually developed into a mainstream predictive modeling method. The enrollment rate of higher education is influenced by many factors such as social, economic, cultural, and geographical factors, which are characterized by partly known and partly unknown information; thus, the GM(1,1) model, as a model to study the uncertainty system problem of “partly known and partly unknown,” has been used to predict the enrollment rate of higher education. For example, Liu et al. [7] used the GM(1, 1) prediction model to make a mathematical analysis of the gross enrollment rate of higher education. It scientifically predicted and eased the employment situation of college students. It provided reference for promoting social stability and economic development. Dong [8] used the GM(1, 1) model to predict the fluctuation of China’s medium- and long-term gross enrollment ratio in higher education, which provides a reference for government assessment and citizens’ choice of higher education.The GM(1, 1), as the first gray prediction model proposed to be applicable to univariate modeling, has shown favorable simulation and prediction performance for the sequence of unity exponential growth characteristics. However, according to the Statistical Yearbook, the gross enrollment rate of higher education in China was 51.6% in 2019, which exceeded the 50% mark for the first time, indicating that China’s higher education has moved from the mass stage to the generalized stage. The shift in development scale has undoubtedly increased the complexity and uncertainty of the higher education GER system, which makes the system behavior series more similar to nonuniform exponential growth characteristics, and the original system behavior series with similar unity exponential growth characteristics can no longer satisfy the existing system. How to develop a high-performance gray prediction model that is more flexible, adaptable, and suitable for modeling and predicting the gross enrollment ratio in higher education dynamically is the focus of this research.In recent years, researchers have carried out a lot of research on the modeling ability of approximate nonhomogeneous exponential sequences. These studies can be divided into three categories according to different modeling ideas. The first is to optimize the structure of the GM(1,1) model to ensure that the final reduction formula of the model presents an approximate homogeneous exponential form. The NHGM(1, 1,k) [9] model extends the basic form x0k+az1k=b of the GM(1,1) model to x0k+az1k=0.52k−1b+c. The SIGM model [10] expands the basic form x0k+az1k=b of the GM(1,1) model to x0k+az1k=kb+c. The KRNGM model [11] introduces the nonlinear function ft into the whitening differential equation of the GM(1,1) model. The second is to modify the adaptability of modeling objects. The conversion from a nonhomogeneous exponential sequence to a homogeneous exponential sequence can be realized by the difference between adjacent elements in the original sequence. The third is the direct modeling method. Wang [12] proposed a direct modeling method of the GM(1,1) model generated from the original data without accumulation, and it had been gradually optimized [13, 14]. Ye and Li [15] established a whitening weight gray prediction model to measure the influence of the probability of interval gray number in the prediction results. Zeng and Liu [16] established a direct modeling method of DGM(1, 1) based on the original sequence by omitting the accumulation and subtraction processes. In addition, the models mentioned above have been applied in variate fields. Xiao et al. [17–20] studied parameter optimization of Grey Riccati model and its application in the prediction on the energy consumption and carbon. A new structure of the Gray Verhulst model is proposed by Xiao et al. [21–24], which improves the ability of the gray model to model saturated S-sequences, and it is applied in China’s tight gas production forecasting.The above optimization models around approximate nonaligned exponential sequences have better properties, modeling capabilities, and a wider range of applications than the GM(1,1) model. The parameters (such as order, background value coefficients, and initial values.) of the gray prediction model are crucial to the performance of the model. However, these performance parameters are often simplified to a specific value in the above optimization model. Therefore, optimization of performance parameters is a key means to improve the stability, applicability, and flexibility of the model. A large number of studies have been launched for the optimization of performance parameters from many directions, such as initial value [25–27], order [28, 29], and background value [30, 31]. The optimal value of each performance parameter should meet the condition of the minimum sum of squared simulation errors of the model, and the optimization process requires a lot of calculations. The Particle Swarm Algorithm (PSA) provides an optional solution for the optimization of those parameters, which find the global optimum by following the currently searched optimum. A distinctive feature compared to other modern optimization methods is that the PSA requires few parameters to be adjusted, is simple and easy to implement, and converges quickly, which has become a research hotspot in the field of modern optimization methods [32–35]. Zeng and Liu [36] proposed the SAIGM_FO model and used the PSA to optimize the order of the model. Wang and Li [37] used the PSA to optimize the structural parameters (a, b) of the model. It has been proved that the PSA in the gray prediction model improves the modeling accuracy and model flexibility.In this study, it is firstly proposed that the PSA is applied in the Three-parameter Whitenization Grey Model (TWGM(1,1)) [38] to optimize the order and background value coefficients under the condition of minimum simulation error squared. The solution is not only capable to solve the problem that the cumulative order is fixed to “1” and the background value coefficient value is fixed to “0.5” but also displays a feature of flexible, adaptable, and dynamically adjustable in the application in modeling China’s higher education gross enrollment ratio (GER). The empirical analysis shows that the model has better overall simulation and prediction performance than others. On the one hand, the parametric optimization model significantly improves its own performance, and on the other hand, its intelligent and adjustable adaptivity improves the accuracy and further extends its application. The rest of the study is organized as below. Section 2 describes the original Three-parameter Whitenization Grey Model. Then, the parameter optimization process of the Three-parameter Whitenization Grey Model is followed in Section 3, and Section 4 gives an application of the optimization model on China’s higher education gross enrollment rate. The conclusion of the study is revealed in the last section.


## 2. Original TWGM(1,1) Model

According to the derivation process of the time-response function of the classic GM(1,1) model, the whitenization equation of the Three-parameter Whitenization Grey Model is established, and then the time-response function is derived by solving the differential equation. This model is called the Three-parameter Whitenization Grey Model, or TWGM (1,1) model in short.Definition 1.
Suppose sequenceX0=x01,x02,…,x0n is the original data sequence, where x0k≥0, k=1,2,…n, and X1=x11,x12,…,x1n is the 1-AGO sequence of X0, where(1)x1k=∑i=1kx0i,k=1,2,…n.
Z1=z12,z13,…,z1n is the immediate mean-generating sequence of X1, where(2)z1k=0.5×x1k+x1k−1,k=2,3…,n.Definition 2.
Assuming that the sequencesX0 and X1 are as shown in Definition 1, then(3)dx1dt+ax1=bt+c.
Equation (3) is the whitening differential equation of the TWGM(1,1) model. The equation(4)x0k+az1k=0.52k−1b+c,is the basic form of the TWGM(1,1) model.
According to the modeling ideas of the GM(1,1), the model parametersa, b, and c are estimated through the basic form of the TWGM(1,1) model. The time-response function of the TWGM(1,1) model is obtained by solving the differential equation.Theorem 1.
Suppose the sequenceX0,X1, and Z1 are as shown in Definition 1, p^=a,b,cT is the parameter list, and(5)Y=x02x03⋮x0n,B=−z12321−z13521⋮⋮⋮−z1n2n−121.
Then, the least-square estimation parameter list of the TWGM(1,1) model is satisfied.Theorem 2.
The time-response function of the TWGM(1,1) model is(6)x1t=x11−ba+ba2−cae−at−1+bat−ba2+ca.
The inverse-accumulating reduction formula of the TWGM(1,1) model is(7)x^0t=1−eax01−ba+ba2−cae−at−1+ba.


## 3. Parameter Optimization of the Three-Parameter Whitenization Gray Model

### 3.1. Background Value Optimization

In the TWGM(1,1) model modeling process, the sequence generated next to the mean is used as the background value, which is a common smoothing method to weaken the influence of extreme values or outliers in the 1-AGO sequence on the magnitude of the gray effect. Among them, the background value coefficient is the weight of adjacent elements in the process of constructing the adjacent mean sequence. The difference of its value will affect the calculation result of the adjacent mean value series and then affect the simulation and prediction effects of the model. In the above modeling process, the background value coefficient is set to 0.5. This processing method lacks flexibility and cannot guarantee that the model achieves the best simulation effect. So, it is necessary to optimize the background value of the original model. The optimization of the background value is mainly based on the optimization of the background value coefficient. Optimize the definition of the sequence immediately adjacent to the mean.Definition 3.
Suppose thatX1 is as shown in Definition 1, Zw1=zw12,zw13,…,zw1n is the optimized sequence generated by the immediate mean of X1, and w is the background value coefficient,w∈0,1, where(8)zw1k=wx1k+1−wx1k−1,k=2,3,…,n.
Then,(9)x0k+awx1k+1−wx1k−1=kb+c,is the TWGM(1,1) model with the background value coefficient w0<w<1.
Among them, the optimal value of the background value coefficientw should satisfy the minimum sum of squared simulation errors of the model, namely,(10)minfw=1n−1∑k=2nx0k−x^0k2,0<w<1.
In equation (10), x0k is the original modeling data and x^0k is the simulated data of the model.


### 3.2. Order Optimization

In the original modeling process, the cumulative order is fixed at “1,” which leads to poor flexibility and adaptability of the model. Therefore, Wu et al. [39, 40] introduced fractional order into gray modeling based on the “in between” idea. This realizes the expansion of the cumulative order of the gray prediction model from integer to fraction. Meng et al. [38, 41] used the Gamma function to give the functional expression of the fractional accumulation operator, which provided the basis for constructing the fractional gray prediction model. This section optimizes the order of the TWGM(1,1) model in order to improve model performance and model adaptability.Definition 4.
SupposeX0 is as shown in Definition 1 and Xrr∈R+ is a new sequence; then, sequence Xr=xr1,xr2,…,xrn is called the r-order cumulative generating sequence of X0, where(11)xrk=∑i=1kΓr+k−iΓk−i+1Γrx0i,k=1,2,…n.
Then, the sequenceZr=zr2,zr3,…,zrn is called the adjacent mean-value generating sequence of Xr, where(12)zrk=xrk+xrk−12,k=2,3,…,n.
Then,(13)xr−1k+azrk=0.52k−1b+c,is the TWGM(1,1) model of r order.
Among them, the optimal value of ther order should satisfy the minimum sum of squared simulation errors of the model, namely,(14)minfr=1n−1∑k=2nx0k−x^0k2,r∈R+.


### 3.3. Simultaneous Optimization of Background Value and Order

Based on the above optimization analysis, this section derives and constructs a new TWGM(1,1) model under the condition of simultaneous optimization of background value and order. The new model is called the TWGM(1,1|w,r) model.


#### 3.3.1. Definition

Definition 5.
Suppose the sequencesX0 and Xr are as shown in Definitions 1 and 4; then,(15)dxrdt+axr=bt+c,is called the whitening differential equation of the TWGM (1,1|w,r) model.
The sequenceZwr=zwr2,zwr3,…,zwrn is called the optimized sequence generated by the adjacent mean value of Xr under the condition of the background value coefficient w, where(16)zwrk=wxrk+1−wxrk−1,k=2,3,…,n.
Integrating both ends of formula (15) in interval k−1,k, we can obtain(17)∫k−1kdxrt+a∫k−1kxrtdt=∫k−1kbtdt+∫k−1kcdt.
Because(18)∫k−1kdxrt=xrt−xrt−1=xr−1t,∫k−1kbtdt=0.52k−1b,∫k−1kcdt=c,among them, the size of ∫k−1kxrtdt can be approximately replaced by the area represented by zwrt, namely,(19)∫k−1kxrtdt≈wxrt+1−wxrt−1=zwrt.
So, equation (15) can be converted to(20)xr−1k+zwrk=0.52k−1b+c.
Equation (20) is the basic form of the TWGM(1,1|w,r) model.


#### 3.3.2. Parameter Estimation

Theorem 3.
Suppose the sequencesX0,Xr, and Zwr are as shown in Definitions 1, 4, and 5, p^=a,b,cT is the parameter list, and(21)Y=xr−12xr−13⋮xr−1n=xr2−xr1xr3−xr2⋮xrn−xrn−1,B=−zwr2321−zwr3521⋮⋮⋮−zwrn2n−121=−wxr2−1−wxr1321−wxr3−1−wxr2521⋮⋮⋮−wxrn−1−wxrn−12n−121.
Then, the least-square estimation parameter list of the TWGM(1,1|w,r) model satisfies p^=a,b,cT=BTB−1BTY.


#### 3.3.3. Time-Response Derivation

According to the whitening differential equation of the TWGM(1,1|w,r) model to derive its time-response function and according to formula (15), the corresponding homogeneous equation is(22)dxrtdt+axrt=0⇒dxrtdt=−axrt.Then,(23)lnxr=−at+lnC1.The general solution of homogeneous equation (22) is(24)xrt=C1e−at.Using the constant variation method, replace formula (24) C1 with ut, and let(25)xrt=ute−at.Derivation from both ends of equation (25) with respect to t:(26)dxrtdt=u′te−at−aute−at.Substituting equation (26) into formula (15), we obtain(27)u′te−at−aute−at=bt+c−axr.According to equation (27), there is u′t=bt+ceat, so(28)ut=∫bt+ceatdt=b∫teatdt+c∫eatdt=bateat−ba2eat+caeat+C.Substituting equation (28) into (25), we obtain(29)xrt=bateat−ba2eat+caeat+Ce−at.Calculating the above equation,(30)xrt=bat−ba2+ca+Ce−at.Among them,xr1=x01 is a known item. At that time t=1, equation (30) can be obtained:(31)xr1=Ce−a+ba−ba2+ca.We obtain(32)C=xr1−b/a+b/a2−c/ae−a.Substituting equation (32) into (30), we obtain(33)xrt=xr1−ba+ba2−cae−at−1+bat−ba2+ca.Among them,t=2,3,…,n and equation (33) is the time-response formula. At this point, the time-response derivation is over.


#### 3.3.4. Derivation of Accumulative Reduction

Definition 6.
Suppose sequenceX0 is as shown in Definition 1 and if r∈R+, then sequence X−r=x−r1,x−r2,…,x−rn is called the r-order accumulative generating sequence of X0, where(34)x−rk=∑i=0k−1−1iΓr+1Γi+1Γr−i+1x0k−i,k=1,2,…,n.Theorem 4.
Suppose sequenceX0 is as defined in Definition 1, p∈R+，q∈R+, Xp is the qth order cumulative generation sequence of X0, Xq is the pth order cumulative generation sequence of X0, Xp+q is the p + q-order cumulative generation sequence of X0, Xpq is the q-order cumulative generation sequence of Xp, and Xqp is the p-order cumulative generating sequence of Xq; then, the multiple cumulative generating operator satisfies the commutative law and exponential rate, namely,(35)Xpq=Xqp=Xp+q.Corollary 1.
According to Theorem1, it can be derived as(36)X0=Xr−r=X−rr.
According to Corollary1, which is equation (36), the final reduction formula of the TWGM(1,1|w,r) model can be derived as(37)x^0k=x^rk−rk=∑i=0k−1−1iΓr+1Γi+1Γr−i+1x^rk−i.
In equation (37), when k=2,3,…,n, x^0k is called the simulated value and when k=n+1,n+2,…, x^0k is called the predicted value.


### 3.4. TWGM(1,1) Model Parameter Optimization Based on PSO Algorithm

In the TWGM(1,1|w,r) model, there are two undetermined parameters (the background value coefficient and order). In order to obtain the best performance of the TWGM(1,1|w,r) model, the optimal values of these two undetermined parameters are required. The combined solution of the two parameters is easy to cause errors due to mutual influence and mutual interference of the parameters, and at the same time, it weakens the optimization effect of the independent parameter on the model. Therefore, in the process of solving the optimal value of the parameter, a step-by-step method is adopted. Specifically, first, in 0<w<1, assuming that the order does not change, that is, r = 1, the optimal background value coefficient is solved. Secondly, on the basis of w=w∗, in r∈R+, the optimal order is solved. According to Definitions 3 and 4, the optimal background value coefficient w and the optimal order r are found within their respective value ranges and should be obtained under the condition that the sum of squared simulation errors of the model is minimized, that is, to satisfy equations (10) and (14).Obviously, the optimization process of each parameter to be determined takes a lot of time and takes up limited computer resources. Various group optimization algorithms (such as Particle Swarm Optimization and ant colony algorithm.) provide good solutions to complex distributed optimization problems. The Particle Swarm Optimization (PSO) is a swarm-intelligent global optimization algorithm that simulates the predation behavior of birds. Its basic concept comes from the study of the foraging behavior of birds. The algorithm has the advantages of simple structure, few parameters, and easy programming. At the same time, the Particle Swarm Optimization algorithm based on adaptive mutation of the population fitness variance effectively solves the phenomenon of premature convergence and can significantly improve the global convergence performance. It has been used in function optimization and neural network. It is widely used in training, engineering, and other fields [28–31]. Therefore, this paper uses PSO to optimize the background value coefficient and order step by step. The PSO algorithm solving steps are as follows:Step 1. Initialize randomly the position and velocity of the particles in the particle swarm.Step 2. Set the particle in the current position to the position of the best particle in the initial population.Step 3. Calculate the average relative simulation error of the TWGM(1,1|w,r) model when w=pBest(or r=pBest):Step 3.1. Calculate the cumulative sequence Xr of order r.Step 3.2. Calculate the sequence Zwr immediately adjacent to the background value coefficient w.Step 3.3. Construct matrix B and Y, and solve model parameter p^=a,b,cT.Step 3.4. Calculate the simulation value X^0.Step 3.5. Calculate the average relative simulation error fpBest of X^0.Step 3.6. Judge whether fpBest−fgBest is less than the given convergence value δ; if it is satisfied, go to Step 9. Otherwise, execute Step 4.Step 4. Perform the following operations for all particles in the particle swarm:Step 4.1. Update particle position and velocity:(38)V=ω×V+c1×rand×pBest−Present+c2×rand×gBest−Present,Present=Present+V,ω=ωmax−run×ωmax−ωminrunMax.Step 4.2. If the particle fitness is better than the fitness, set it to the new position.Step 4.3. If the particle fitness is better than the fitness, set it to the new position.Step 5. Calculate the variance of the group fitness σ2, and calculate fpBest:(39)σ2=∑i=1nfi−favgf2,f=maxfi−favg,maxfi−favg>1,1,others.Step 6. Calculate the probability of mutation pm.Step 7. Generate random number ε∈0,1; if ε<pm, press to execute mutation operation. Otherwise, go to Step 8.Step 8. Determine whether the algorithm convergence criterion is satisfied; if it is satisfied, execute Step 9. Otherwise, turn to Step 3.Step 9. Output the optimal value gBest of the background value coefficient w(or order r) and the simulation and prediction data of the TWGM(1,1|w,r) model at this time, and the algorithm operation ends.


## 3.1. Background Value Optimization

In the TWGM(1,1) model modeling process, the sequence generated next to the mean is used as the background value, which is a common smoothing method to weaken the influence of extreme values or outliers in the 1-AGO sequence on the magnitude of the gray effect. Among them, the background value coefficient is the weight of adjacent elements in the process of constructing the adjacent mean sequence. The difference of its value will affect the calculation result of the adjacent mean value series and then affect the simulation and prediction effects of the model. In the above modeling process, the background value coefficient is set to 0.5. This processing method lacks flexibility and cannot guarantee that the model achieves the best simulation effect. So, it is necessary to optimize the background value of the original model. The optimization of the background value is mainly based on the optimization of the background value coefficient. Optimize the definition of the sequence immediately adjacent to the mean.Definition 3.
Suppose thatX1 is as shown in Definition 1, Zw1=zw12,zw13,…,zw1n is the optimized sequence generated by the immediate mean of X1, and w is the background value coefficient,w∈0,1, where(8)zw1k=wx1k+1−wx1k−1,k=2,3,…,n.
Then,(9)x0k+awx1k+1−wx1k−1=kb+c,is the TWGM(1,1) model with the background value coefficient w0<w<1.
Among them, the optimal value of the background value coefficientw should satisfy the minimum sum of squared simulation errors of the model, namely,(10)minfw=1n−1∑k=2nx0k−x^0k2,0<w<1.
In equation (10), x0k is the original modeling data and x^0k is the simulated data of the model.


## 3.2. Order Optimization

In the original modeling process, the cumulative order is fixed at “1,” which leads to poor flexibility and adaptability of the model. Therefore, Wu et al. [39, 40] introduced fractional order into gray modeling based on the “in between” idea. This realizes the expansion of the cumulative order of the gray prediction model from integer to fraction. Meng et al. [38, 41] used the Gamma function to give the functional expression of the fractional accumulation operator, which provided the basis for constructing the fractional gray prediction model. This section optimizes the order of the TWGM(1,1) model in order to improve model performance and model adaptability.Definition 4.
SupposeX0 is as shown in Definition 1 and Xrr∈R+ is a new sequence; then, sequence Xr=xr1,xr2,…,xrn is called the r-order cumulative generating sequence of X0, where(11)xrk=∑i=1kΓr+k−iΓk−i+1Γrx0i,k=1,2,…n.
Then, the sequenceZr=zr2,zr3,…,zrn is called the adjacent mean-value generating sequence of Xr, where(12)zrk=xrk+xrk−12,k=2,3,…,n.
Then,(13)xr−1k+azrk=0.52k−1b+c,is the TWGM(1,1) model of r order.
Among them, the optimal value of ther order should satisfy the minimum sum of squared simulation errors of the model, namely,(14)minfr=1n−1∑k=2nx0k−x^0k2,r∈R+.


## 3.3. Simultaneous Optimization of Background Value and Order

Based on the above optimization analysis, this section derives and constructs a new TWGM(1,1) model under the condition of simultaneous optimization of background value and order. The new model is called the TWGM(1,1|w,r) model.


### 3.3.1. Definition

Definition 5.
Suppose the sequencesX0 and Xr are as shown in Definitions 1 and 4; then,(15)dxrdt+axr=bt+c,is called the whitening differential equation of the TWGM (1,1|w,r) model.
The sequenceZwr=zwr2,zwr3,…,zwrn is called the optimized sequence generated by the adjacent mean value of Xr under the condition of the background value coefficient w, where(16)zwrk=wxrk+1−wxrk−1,k=2,3,…,n.
Integrating both ends of formula (15) in interval k−1,k, we can obtain(17)∫k−1kdxrt+a∫k−1kxrtdt=∫k−1kbtdt+∫k−1kcdt.
Because(18)∫k−1kdxrt=xrt−xrt−1=xr−1t,∫k−1kbtdt=0.52k−1b,∫k−1kcdt=c,among them, the size of ∫k−1kxrtdt can be approximately replaced by the area represented by zwrt, namely,(19)∫k−1kxrtdt≈wxrt+1−wxrt−1=zwrt.
So, equation (15) can be converted to(20)xr−1k+zwrk=0.52k−1b+c.
Equation (20) is the basic form of the TWGM(1,1|w,r) model.


### 3.3.2. Parameter Estimation

Theorem 3.
Suppose the sequencesX0,Xr, and Zwr are as shown in Definitions 1, 4, and 5, p^=a,b,cT is the parameter list, and(21)Y=xr−12xr−13⋮xr−1n=xr2−xr1xr3−xr2⋮xrn−xrn−1,B=−zwr2321−zwr3521⋮⋮⋮−zwrn2n−121=−wxr2−1−wxr1321−wxr3−1−wxr2521⋮⋮⋮−wxrn−1−wxrn−12n−121.
Then, the least-square estimation parameter list of the TWGM(1,1|w,r) model satisfies p^=a,b,cT=BTB−1BTY.


### 3.3.3. Time-Response Derivation

According to the whitening differential equation of the TWGM(1,1|w,r) model to derive its time-response function and according to formula (15), the corresponding homogeneous equation is(22)dxrtdt+axrt=0⇒dxrtdt=−axrt.Then,(23)lnxr=−at+lnC1.The general solution of homogeneous equation (22) is(24)xrt=C1e−at.Using the constant variation method, replace formula (24) C1 with ut, and let(25)xrt=ute−at.Derivation from both ends of equation (25) with respect to t:(26)dxrtdt=u′te−at−aute−at.Substituting equation (26) into formula (15), we obtain(27)u′te−at−aute−at=bt+c−axr.According to equation (27), there is u′t=bt+ceat, so(28)ut=∫bt+ceatdt=b∫teatdt+c∫eatdt=bateat−ba2eat+caeat+C.Substituting equation (28) into (25), we obtain(29)xrt=bateat−ba2eat+caeat+Ce−at.Calculating the above equation,(30)xrt=bat−ba2+ca+Ce−at.Among them,xr1=x01 is a known item. At that time t=1, equation (30) can be obtained:(31)xr1=Ce−a+ba−ba2+ca.We obtain(32)C=xr1−b/a+b/a2−c/ae−a.Substituting equation (32) into (30), we obtain(33)xrt=xr1−ba+ba2−cae−at−1+bat−ba2+ca.Among them,t=2,3,…,n and equation (33) is the time-response formula. At this point, the time-response derivation is over.


### 3.3.4. Derivation of Accumulative Reduction

Definition 6.
Suppose sequenceX0 is as shown in Definition 1 and if r∈R+, then sequence X−r=x−r1,x−r2,…,x−rn is called the r-order accumulative generating sequence of X0, where(34)x−rk=∑i=0k−1−1iΓr+1Γi+1Γr−i+1x0k−i,k=1,2,…,n.Theorem 4.
Suppose sequenceX0 is as defined in Definition 1, p∈R+，q∈R+, Xp is the qth order cumulative generation sequence of X0, Xq is the pth order cumulative generation sequence of X0, Xp+q is the p + q-order cumulative generation sequence of X0, Xpq is the q-order cumulative generation sequence of Xp, and Xqp is the p-order cumulative generating sequence of Xq; then, the multiple cumulative generating operator satisfies the commutative law and exponential rate, namely,(35)Xpq=Xqp=Xp+q.Corollary 1.
According to Theorem1, it can be derived as(36)X0=Xr−r=X−rr.
According to Corollary1, which is equation (36), the final reduction formula of the TWGM(1,1|w,r) model can be derived as(37)x^0k=x^rk−rk=∑i=0k−1−1iΓr+1Γi+1Γr−i+1x^rk−i.
In equation (37), when k=2,3,…,n, x^0k is called the simulated value and when k=n+1,n+2,…, x^0k is called the predicted value.


## 3.3.1. Definition

Definition 5.
Suppose the sequencesX0 and Xr are as shown in Definitions 1 and 4; then,(15)dxrdt+axr=bt+c,is called the whitening differential equation of the TWGM (1,1|w,r) model.
The sequenceZwr=zwr2,zwr3,…,zwrn is called the optimized sequence generated by the adjacent mean value of Xr under the condition of the background value coefficient w, where(16)zwrk=wxrk+1−wxrk−1,k=2,3,…,n.
Integrating both ends of formula (15) in interval k−1,k, we can obtain(17)∫k−1kdxrt+a∫k−1kxrtdt=∫k−1kbtdt+∫k−1kcdt.
Because(18)∫k−1kdxrt=xrt−xrt−1=xr−1t,∫k−1kbtdt=0.52k−1b,∫k−1kcdt=c,among them, the size of ∫k−1kxrtdt can be approximately replaced by the area represented by zwrt, namely,(19)∫k−1kxrtdt≈wxrt+1−wxrt−1=zwrt.
So, equation (15) can be converted to(20)xr−1k+zwrk=0.52k−1b+c.
Equation (20) is the basic form of the TWGM(1,1|w,r) model.


## 3.3.2. Parameter Estimation

Theorem 3.
Suppose the sequencesX0,Xr, and Zwr are as shown in Definitions 1, 4, and 5, p^=a,b,cT is the parameter list, and(21)Y=xr−12xr−13⋮xr−1n=xr2−xr1xr3−xr2⋮xrn−xrn−1,B=−zwr2321−zwr3521⋮⋮⋮−zwrn2n−121=−wxr2−1−wxr1321−wxr3−1−wxr2521⋮⋮⋮−wxrn−1−wxrn−12n−121.
Then, the least-square estimation parameter list of the TWGM(1,1|w,r) model satisfies p^=a,b,cT=BTB−1BTY.


## 3.3.3. Time-Response Derivation

According to the whitening differential equation of the TWGM(1,1|w,r) model to derive its time-response function and according to formula (15), the corresponding homogeneous equation is(22)dxrtdt+axrt=0⇒dxrtdt=−axrt.Then,(23)lnxr=−at+lnC1.The general solution of homogeneous equation (22) is(24)xrt=C1e−at.Using the constant variation method, replace formula (24) C1 with ut, and let(25)xrt=ute−at.Derivation from both ends of equation (25) with respect to t:(26)dxrtdt=u′te−at−aute−at.Substituting equation (26) into formula (15), we obtain(27)u′te−at−aute−at=bt+c−axr.According to equation (27), there is u′t=bt+ceat, so(28)ut=∫bt+ceatdt=b∫teatdt+c∫eatdt=bateat−ba2eat+caeat+C.Substituting equation (28) into (25), we obtain(29)xrt=bateat−ba2eat+caeat+Ce−at.Calculating the above equation,(30)xrt=bat−ba2+ca+Ce−at.Among them,xr1=x01 is a known item. At that time t=1, equation (30) can be obtained:(31)xr1=Ce−a+ba−ba2+ca.We obtain(32)C=xr1−b/a+b/a2−c/ae−a.Substituting equation (32) into (30), we obtain(33)xrt=xr1−ba+ba2−cae−at−1+bat−ba2+ca.Among them,t=2,3,…,n and equation (33) is the time-response formula. At this point, the time-response derivation is over.


## 3.3.4. Derivation of Accumulative Reduction

Definition 6.
Suppose sequenceX0 is as shown in Definition 1 and if r∈R+, then sequence X−r=x−r1,x−r2,…,x−rn is called the r-order accumulative generating sequence of X0, where(34)x−rk=∑i=0k−1−1iΓr+1Γi+1Γr−i+1x0k−i,k=1,2,…,n.Theorem 4.
Suppose sequenceX0 is as defined in Definition 1, p∈R+，q∈R+, Xp is the qth order cumulative generation sequence of X0, Xq is the pth order cumulative generation sequence of X0, Xp+q is the p + q-order cumulative generation sequence of X0, Xpq is the q-order cumulative generation sequence of Xp, and Xqp is the p-order cumulative generating sequence of Xq; then, the multiple cumulative generating operator satisfies the commutative law and exponential rate, namely,(35)Xpq=Xqp=Xp+q.Corollary 1.
According to Theorem1, it can be derived as(36)X0=Xr−r=X−rr.
According to Corollary1, which is equation (36), the final reduction formula of the TWGM(1,1|w,r) model can be derived as(37)x^0k=x^rk−rk=∑i=0k−1−1iΓr+1Γi+1Γr−i+1x^rk−i.
In equation (37), when k=2,3,…,n, x^0k is called the simulated value and when k=n+1,n+2,…, x^0k is called the predicted value.


## 3.4. TWGM(1,1) Model Parameter Optimization Based on PSO Algorithm

In the TWGM(1,1|w,r) model, there are two undetermined parameters (the background value coefficient and order). In order to obtain the best performance of the TWGM(1,1|w,r) model, the optimal values of these two undetermined parameters are required. The combined solution of the two parameters is easy to cause errors due to mutual influence and mutual interference of the parameters, and at the same time, it weakens the optimization effect of the independent parameter on the model. Therefore, in the process of solving the optimal value of the parameter, a step-by-step method is adopted. Specifically, first, in 0<w<1, assuming that the order does not change, that is, r = 1, the optimal background value coefficient is solved. Secondly, on the basis of w=w∗, in r∈R+, the optimal order is solved. According to Definitions 3 and 4, the optimal background value coefficient w and the optimal order r are found within their respective value ranges and should be obtained under the condition that the sum of squared simulation errors of the model is minimized, that is, to satisfy equations (10) and (14).Obviously, the optimization process of each parameter to be determined takes a lot of time and takes up limited computer resources. Various group optimization algorithms (such as Particle Swarm Optimization and ant colony algorithm.) provide good solutions to complex distributed optimization problems. The Particle Swarm Optimization (PSO) is a swarm-intelligent global optimization algorithm that simulates the predation behavior of birds. Its basic concept comes from the study of the foraging behavior of birds. The algorithm has the advantages of simple structure, few parameters, and easy programming. At the same time, the Particle Swarm Optimization algorithm based on adaptive mutation of the population fitness variance effectively solves the phenomenon of premature convergence and can significantly improve the global convergence performance. It has been used in function optimization and neural network. It is widely used in training, engineering, and other fields [28–31]. Therefore, this paper uses PSO to optimize the background value coefficient and order step by step. The PSO algorithm solving steps are as follows:Step 1. Initialize randomly the position and velocity of the particles in the particle swarm.Step 2. Set the particle in the current position to the position of the best particle in the initial population.Step 3. Calculate the average relative simulation error of the TWGM(1,1|w,r) model when w=pBest(or r=pBest):Step 3.1. Calculate the cumulative sequence Xr of order r.Step 3.2. Calculate the sequence Zwr immediately adjacent to the background value coefficient w.Step 3.3. Construct matrix B and Y, and solve model parameter p^=a,b,cT.Step 3.4. Calculate the simulation value X^0.Step 3.5. Calculate the average relative simulation error fpBest of X^0.Step 3.6. Judge whether fpBest−fgBest is less than the given convergence value δ; if it is satisfied, go to Step 9. Otherwise, execute Step 4.Step 4. Perform the following operations for all particles in the particle swarm:Step 4.1. Update particle position and velocity:(38)V=ω×V+c1×rand×pBest−Present+c2×rand×gBest−Present,Present=Present+V,ω=ωmax−run×ωmax−ωminrunMax.Step 4.2. If the particle fitness is better than the fitness, set it to the new position.Step 4.3. If the particle fitness is better than the fitness, set it to the new position.Step 5. Calculate the variance of the group fitness σ2, and calculate fpBest:(39)σ2=∑i=1nfi−favgf2,f=maxfi−favg,maxfi−favg>1,1,others.Step 6. Calculate the probability of mutation pm.Step 7. Generate random number ε∈0,1; if ε<pm, press to execute mutation operation. Otherwise, go to Step 8.Step 8. Determine whether the algorithm convergence criterion is satisfied; if it is satisfied, execute Step 9. Otherwise, turn to Step 3.Step 9. Output the optimal value gBest of the background value coefficient w(or order r) and the simulation and prediction data of the TWGM(1,1|w,r) model at this time, and the algorithm operation ends.


## 4. Model Application: Forecast of China’s Higher Education Gross Enrollment Rate

To identify the optimized performance on each parameter, the variables (the background value coefficientw and order r) are simulated with four different groups. The original TWGM(1,1) model means that the model is not optimized for order and background value coefficients at all. The first group is the model with, w = 0.5 and r = 1, that is, the TWGM(1,1|0.5,1)model. The second one is the model with w = 0.5 and r = r∗, that is, the TWGM(1,1|0.5,r∗) model, which indicates that the model only optimizes the order. The third one is the model with, w = w∗ and r = 1, that is, the TWGM(1,1|w∗,1) model, which means that the model only optimizes the background value coefficient. The last group is w = w∗ and r = r∗, that is, the TWGM(1,1|w∗,r∗) model, which suggests that the model optimizes the order and background value coefficients at the same time. The simulation results with Matlab for four parameter groups are shown in Table 1, where x0k is China’s high education enrollment data from 1991 to 2019.Table 1
Simulation and prediction results of different models on China’s higher education gross enrollment rate.
Serial numberx0kTWGM(1,1|0.5,1)TWGM(1,1|0.5,r∗)TWGM(1,1|w∗,1)TWGM(1,1|w∗,r∗)x^0kΔkx^0kΔkx^0kΔkx^0kΔkSimulation datak = 23.53.1410.21883.1310.67253.1510.10193.1011.5502k = 33.94.043.50784.064.16204.043.67834.053.7288k = 454.970.62785.000.03114.980.44144.980.4818k = 565.941.01135.960.61185.950.80975.931.1729k = 67.26.953.46456.963.28066.973.25616.923.9439k = 78.38.003.56618.003.56418.023.34887.944.3312k = 89.19.100.01599.090.13249.120.24889.011.0183k = 99.810.244.537810.224.269210.274.787710.123.2693k = 1010.511.448.911311.408.558811.469.177411.287.4608k = 1112.512.681.411112.631.056212.711.663412.500.0000k = 1213.313.975.029313.924.674314.005.294513.773.5632k = 131515.322.102415.271.806315.352.363715.110.7258k = 141716.721.658616.681.863416.761.404016.512.8898k = 151918.184.319618.164.410218.234.069417.985.3818k = 162119.706.183219.716.137619.755.935719.527.0515k = 172221.293.239421.333.026421.342.982121.143.9173k = 182322.940.263923.040.153323.000.003222.840.7002k = 1923.324.665.838124.826.517024.736.123424.635.6917Forecast datak = 2024.226.459.310326.6910.281626.529.606826.509.5236k = 2126.528.326.870828.658.108828.407.162328.487.4706k = 2226.930.2712.514630.7114.146030.3512.823130.5613.5932k = 233032.297.644332.869.541532.387.940932.749.1365k = 2434.534.400.275935.131.813134.500.000035.041.5616k = 2537.536.602.388137.500.000036.712.116937.460.1167k = 264038.902.760139.990.021039.002.488940.000.0000k = 2742.741.283.318342.610.219541.403.047642.680.0541k = 2845.743.774.223445.350.765043.893.954345.490.4502k = 2948.146.363.616348.230.271446.493.344648.460.7477k = 3051.649.064.923751.250.671249.204.654851.580.0365In order to intuitively reflect the simulation effects of different models, a comparison chart of the simulation curves of different models is drawn, as shown in Figure1, where models with different parameter values show different simulation effects. The curve of the TWGM(1,1| 0.5, 1) and TWGM(1,1|w∗,1) models overlap and show a same trend. The TWGM(1,1|w∗,r∗) model initially displays better simulating and predicting effects for the gross enrollment ratio in higher education in China since it produces the most close result with raw data.In order to measure the optimization performance of different parameters, the average relative simulation percentage errorΔ¯S, average relative prediction percentage error Δ¯F, and comprehensive percentage error Δ¯ of the gray model of each model are calculated according to the results in Table 1. The results are shown in Table 2. The relevant calculation equation is as follows:(40)Δ¯S=1t−1∑k=2tΔk,Δ¯F=1n−t∑k=t+1nΔk,Δ¯=1n−1∑k=2nΔk.Table 2
Error comparison of different models.
Error typeTWGM(1,1|0.5,1) (%)TWGM(1,1|0.5,r∗) (%)TWGM(1,1|w∗,1) (%)TWGM(1,1|w∗,r∗) (%)Average relative simulation percentage errorΔ¯s3.66153.60713.64943.7155Average relative prediction percentage errorΔ¯F5.25874.16725.19463.8810Grey model comprehensive percentage errorΔ¯4.26733.81954.23553.7782Among them,Δk represents the simulation/prediction error of x^0k, n represents the number of experimental data, and t represents the number of experimental data used for simulation.Based on the results of the calculations in Table2, the relative percentage error comparison figure with different parameter values is plotted, as shown in Figure 2 in supplementary materials. From Table 2 and Figure 2 in supplementary materials, the following conclusions are drawn:(1)
The background value optimization effect: comparing the TWGM(1,1|0.5,1) model and TWGM(1,1|w∗,1) model with a w∗ = 0.4538, the graph shows that the integrated percent error of the gray model after the background value optimization is slightly lower than the original model, but the difference is not remarkable, which suggests that the background value optimization has an optimization effect on the comprehensive performance of China’s higher education gross enrollment rate, but the effect is not significant.(2)
The order optimization effect: compare the TWGM(1,1|0.5,1) model with the TWGM(1,1|0.5,r∗) model with a r∗ = 1.1127 in Figure 2 in supplementary materials, the average relative percentage prediction error and the comprehensive percentage error after the order optimization are significantly lower than the original model, which shows that the order optimization has a remarkable improvement on the model performance in simulating China’s higher education gross enrollment rate.(3)
The TWGM(1,1|w∗,r∗) model possesses the lowest error of 3.7782%, which indicates that the optimization of both the background value and the order is of the best option for improvement on the model performance. Thus, the TWGM(1,1|w∗,r∗) model provides the best fitting effect and is preferred for simulation and prediction.


## 5. Conclusions

The study proposed a new optimization model by combining the Three-parameter Whitenization Grey Model with the PSA. To demonstrate the performance of this new model, an application analysis is carried out in the predication of China’s high education enrollment ratio with Matlab simulation. The findings suggest that the new optimization model proposed in this study improves the performance of the original gray model remarkably. Specifically, compared with the original gray model, the average relative percentage prediction error and comprehensive percentage error after the optimization are significantly reduced, which indicates that the simultaneous optimization of the background value and order has a notable improvement on the overall performance in simulating and predicting China’s higher education gross enrollment ratio. Besides, it is indicated that the intelligent and adjustable adaptivity improves the accuracy of optimization and, furthermore, extends its application in the empirical study in the future.



---

*Source: 6640000-2020-12-19.xml*