GM(1,1;λ) with Constrained Linear Least Squares

: The only parameters of the original GM(1,1) that are generally estimated by the ordinary least squares method are the development coefficient a and the grey input b . However, the weight of the background value, denoted as λ , cannot be obtained simultaneously by such a method. This study, therefore, proposes two simple transformation formulations such that the unknown parameters � , � and � can be simultaneously estimated by the least squares method. Therefore, such a grey model is termed the GM(1,1;λ). On the other hand, because the permission zone of the development coefficient is bounded, the parameter estimation of the GM(1,1) could be regarded as a bound-constrained least squares problem. Since constrained linear least squares problems generally can be solved by an iterative approach, this study applies the Matlab function lsqlin to solve such constrained problems. Numerical results show that the proposed GM(1,1;λ) performs better than the GM(1,1) in terms of its model fitting accuracy and its forecasting precision.


Introduction
Uncertain systems with small samples and poor information exist commonly in the real world.Among them, a system with partially known (white) and partially unknown (black) information is generally regarded as a grey system.Grey system theory, introduced by Deng in the early 1980s, is a methodology that focuses on the study of uncertain systems through generating, excavating, and extracting useful information from what is available [1,2].The research areas of the grey system theory can generally be categorized into the following six fields: grey generating, grey relational analysis, grey models, grey prediction, grey decision making and grey control [3].The grey model is the kernel of grey prediction, where the latter is one of the most important and widely used fields in the grey system theory.In the class of grey models, the GM (1,1), namely the first-order and single-variable grey model, is the most basic and widely used model.It requires a relatively small amount of data (typically four or more) to construct a mathematical model, and then employs a simple calculation process to estimate the behavior of unknown systems.Generally speaking, the GM(1,1) has a high forecasting accuracy and has been successfully applied to many fields, such as prediction control [4], power load forecasting [5], fuel production prediction [6,7], agricultural output [8], energy demand forecasting [9][10][11], etc.
The response of the classic GM (1,1) is essentially an exponential model with two internal parameters: the development coefficient a and the grey input b.Generally speaking, these internal parameters are estimated by the least squares method with the background values being part of the observed data.In the grey model, the background values are obtained by a mean generating operation, which is generally defined as the equal-weighted average of two selected data.Such an equal-weighted average formula also gives rise to the motivation to optimize the weight of the background value, hereafter denoted as λ, to improve the precision of the fitting and prediction of GM(1,1) [12][13][14].Until now, many researchers have proposed improvements to GM (1,1) to overcome these problems.The improved methods can generally be divided into the following three categories: (i) compensation with a residual GM(1,1) model [2], (ii) a combination with different expert schemes, e.g., fuzzy control [15], (iii) optimization of internal parameters a, b, and/or λ by evolutionary algorithms, such as the genetic algorithm (GA) [8,12,14] and differential evolution (DE) [13].Rather than using the classic GM(1,1), the adaptive GA-based optimization approach proposed in [7] has been employed to simultaneously optimize the parameters a, b and λ of GM(1,1;λ) to improve accuracy compared with the traditional GM (1,1).
Linear regression models are often fitted using the least squares approach.Thus, the internal parameters of a classic GM(1,1), i.e., the development coefficient a and the grey input b, are generally estimated by the ordinary least squares method.Since both a and b are bounded [16], the parameter estimation could be regarded as a bound-constrained least squares problem.Linear least squares problems with bound constraints are commonly solved to find model parameters within bounds based on physical considerations.Common algorithms include bounded-variable least squares (BVLS) and the Matlab function lsqlin [17].While the parameter estimation of the GM(1,1) also takes the weight of the background value λ into consideration, the issue will become a nonlinear regression problem.Since some nonlinear regression problems can be moved to a linear domain by a suitable transformation of the model formulation [15], this study, therefore, proposes two simple transformation formulations to solve this problem.By the proposed transformation, the unknown parameters, a, b and λ, can also be estimated by the ordinary least squares method.In addition, this study also attempts to improve the precision of the fitting and prediction of the GM(1,1) by the proposed parameter estimation approach.
The remainder of this paper is organized as follows.Section 2 briefly describes the GM(1,1) and linear least squares problems.The proposed GM(1,1;λ) with constrained linear least squares is given in Section 3. Section 4 presents the simulation results on two numerical problems.Finally, Section 5 contains some conclusions of this study.

GM(1,1)
Assume that the raw data sequence is x (0) = x (0) (1), x (0) (2), . . . ,x (0) (n) , where x (0) (k) ≥ 0, ∀k.Performing the accumulated generating operation (AGO) on x (0) yields a new generated sequence x (1) , where the AGO is defined by and x (1) is termed the first-order accumulated generating sequence of x (0) .While x (1) is obtained, the whitenization function of the traditional GM(1,1) can be defined as: where a and b represent the development coefficient and grey input (control variable), respectively [2,3].In the GM(1,1), the forecasting value of x (1) (k) can be determined by the time response function of ( 2) with an initial condition x (1) (1) = x (0) (1), that is, In addition, the grey differential equation corresponding to the whitenization function given in (2) is where Herein, z (1) (k) is termed the background value and λ ∈ [0, 1].The matrix form of the grey differential Equation ( 4) is where Then, a and b could be determined by using the ordinary least squares method: The forecasting value of x (0) (k) can be obtained by applying the inverse AGO (IAGO) to x(1) (k) as follows.

Linear Least Squares Problems
As far as the linear system given in ( 6) is concerned, a classic approach to choosing p is to minimize the least squares (LS) error between Y and Bp: where α and β are given constant vectors and • denotes the Euclidean distance.The previous linear least squares problem with an additional constraint on the solution is termed the bound-constrained least squares problem.Generally speaking, the solution methods for the bound-constrained least squares problems of the form (13) can be categorized as active set or interior point methods [17].In the active set methods, a sequence of equality constrained problems is solved with efficient solution methods.The interior point methods use variants of Newton's method to solve the Karush-Kuhn-Tucker (KKT) equality conditions for (13).
Without considering the constraints α ≤ p ≤ β, the ordinary least squares method given in (10) represents a simple method for solving the LS problem of (13).However, while the parameter estimation of the GM(1,1) takes the weight of the background value λ into consideration, the parameter estimation will become a nonlinear regression problem.Even so, some nonlinear regression problems can be moved to a linear domain by a suitable transformation of the model formulation [18][19][20].With the help of suitable transformation, this study could apply the Matlab function lsqlin to solve the constrained problem as given in (13).Note that the Matlab function lsqlin is a linear least squares solver with bounds or linear constraints.

GM(1,1;λ) with Constrained Linear Least Squares
Upon inspection of (4), the grey differential equation of the GM(1,1) only contains two undetermined parameters (a and b), but does not contain the weight of the background value λ.This is also the main reason that the weighting factor λ cannot be determined by the least squares method given in (10).Substituting ( 5) into (4) yields As can be seen, two nonlinear terms, aλ and a(1 − λ), appear in the new equation ( 14), hereafter termed the grey differential equation of the GM(1,1;λ).Unfortunately, a and λ are not independent variables and they are inseparable.The current grey differential equation of the GM(1,1;λ), therefore, is not in the form of a nonlinear regression model.Hence, it is difficult to apply the method of least squares to estimate the values of the unknown parameters in (14).Thus, the main issue to be solved in this study is how to simultaneously estimate the unknown parameters a, b and λ by the least squares method.

Parameters Estimation of GM(1,1;λ)
A nonlinear model can sometimes be transformed into a linear one.For example, an exponential model, say y = αe βx , where α and β are unknown parameters, can be transformed into a linear model by taking logarithms, i.e., ln y = ln α + βx.Then, the unknown parameters can be estimated by the ordinary least squares method [18][19][20].The main purpose of this section is how to find a suitable transformation of the model formulation (14) such that the nonlinear regression problems can be moved to a linear domain.
Let a 1 = aλ and a 2 = a(1 − λ), where a 1 and a 2 are the components of the development coefficient a.Then, (14) becomes Obviously, ( 15) is a linear model.The corresponding matrix form is where (1) (2) −x (1) (1) 1 −x (1) (3) −x (1) (2) According to the ordinary least squares method, a 1 , a 2 and b can be estimated by Once the components a 1 and a 2 are determined, the development coefficient a and the weighting factor λ are obtained by With the proposed simple transformation, the development coefficient a, the grey input b, and the weight of background value λ can be simultaneously obtained from the ordinary least squares method (19) and the transformation formulations given in (20) and (21).

Boundary Constraint on Estimated Parameters
Generally speaking, the parameter estimation of the classic GM(1,1) is solving a linear least squares problem without any additional constraint on the solution, i.e., a simple unconstrained linear least squares problem.Let λ = 0.5, and then substituting (1) and ( 5) into (10) will yield [16] Assume that the length of the raw data sequence x (0) used to construct a classic GM(1,1) is n and the sequence x (0) is bounded, i.e., 0 max , where x (0) max = max k x (0) (k) and k = 1, 2, . . ., n.Then, it can be derived from ( 22), ( 23) and ( 24) that [16] − 2/(n + 1) < a < 2/(n + 1), (25) It is obvious that (25) and ( 26) are the boundary constraints on a and b, respectively.However, these two constraints are neglected while estimating the parameters a and b by the least squares method (10).Moreover, while taking the weight of the background value λ into consideration, the parameter estimation of the GM(1,1;λ) will become a linear least squares problem with the newly added constraints on the parameters a 1 and a 2 .The main issue of this section is how to properly determine the newly added constraints on the estimated parameters.
It can be inferred from the boundary constraint (25), the condition 0 ≤ λ ≤ 1, and the proposed transformations, a 1 = aλ and a 2 = a(1 − λ), that −2/(n + 1) < a i < 2/(n + 1), i = 1, 2. In addition, the condition 0 ≤ λ ≤ 1 also implies that a 1 , a 2 and a must have the same sign.Then, the parameter estimation of the GM(1,1;λ) can be obtained by one of the following two bound-constrained least squares problems: where C = 2. Our experiences show that the development coefficient a is generally bounded in (−1/(n + 1), +1/(n + 1)), i.e., C = 1 [16].This fact reveals that the constant C could take a value from the interval [1,2].Generally, C = 1.5 or 2.0 in the simulation.Note that this study also applies the Matlab function lsqlin to solve such constrained problems.

Simulation Results
The study selects two numerical examples to verify the effectiveness of the proposed GM(1,1;λ) against the original GM(1,1) and GA-based GM(1,1).In addition, the MAPE (mean absolute percentage error) was employed to measure fitting/forecasting performance as follows [21]: The lower the MAPE value, the more accurate the grey model.Table 1 lists the MAPE criteria for evaluating a forecasting model [22].

Example 1
Example 1 takes inbound tourists' arrivals in Taiwan from 2003 to 2014 as an experimental example [12].Table 2 lists the corresponding real values.The fitting values obtained by the original GM(1,1), GA-based GM(1,1) [12] and the proposed GM(1,1;λ) are also given in the same table.Figure 1 also depicts the real and forecasting values of different forecasting models.As can be seen, the original GM(1,1), i.e., λ = 0.5, has a MAPE of 6.67%.In [12], the parameter settings of the GA-based GM(1,1) are given as follows: population size = 70, number of generations = 100, reproduction probability = 0.85, and mutation probability = 0.005.With the previous parameter settings, the optimal weight of the background value obtained by GA is λ = 0.47323 with a minimum MAPE of 6.58%.While applying the proposed GM(1,1;λ) with C = 1.5 to solve the same problem, the estimated parameters are a 1 = −0.0769, a 2 = −0.0385and b = 2, 334, 485.66.Then, it can be obtained from ( 20) and ( 21) that a = −0.1154and λ = 0.6667.The corresponding overall fitting result (MAPE) is 6.33%.Although all the grey models could provide highly accurate forecasting according to MAPE criteria, the proposed GM(1,1;λ) has the best fitting ability.This fact also reveals that the proposed approach actually improves the model fitting precision of the original GM(1,1).

Example 2
Example 2 takes the processing volumes of crude oil given in [7] as an illustration.In the example, the processing volumes from 1983 to 1992 are selected as the fitting (in-sample) data to examine model fitting ability, whereas the volumes from 1993 to 1994 are selected as the predicting (out-of-sample) data for the ex post testing.The real values and the corresponding fitting and forecasting results obtained by the GM(1,1), GA-based GM(1,1) and the proposed GM(1,1;λ) with C = 1.5 are depicted in Table 3 and Figure 2. It can be seen from Table 3 that the MAPEs of the GM(1,1), GA-based GM(1,1) and the proposed GM(1,1;λ) for model fitting were 4.89%, 3.60% and 4.18%, respectively, whereas the MAPEs were 5.04%, 4.84% and 3.30%, respectively, for the ex post testing.According to the MAPE criteria, all the grey models can attain high forecasting ability in both in-sample and out-of-sample data.The GA-based GM(1,1) has the best fitting accuracy, whereas the proposed GM(1,1;λ) performs the best forecasting accuracy.In addition, the proposed GM(1,1;λ) also performs better than the original GM(1,1) in terms of both fitting and forecasting abilities.That is to say, the proposed constrained least squares method can effectively enhance the fitting and forecasting accuracies of GM(1,1).

Conclusions
This study proposes two simple transformation formulations such that the development coefficient a, the grey input b, and the weight of the background value λ can be simultaneously obtained from the ordinary least squares method.The corresponding boundary constraints are also derived in the study.With the help of the estimation of λ, the proposed approach could also be applied to improve the precision of the fitting and prediction of the original GM(1,1).Two real cases, the inbound arrivals in Taiwan and the processing volumes of crude oil, were used to evaluate the model fitting and forecasting performances of the proposed GM(1,1;λ).Numerical results showed that the GM(1,1), GA-based GM(1,1) and the proposed GM(1,1;λ) could provide highly accurate forecasting according to the MAPE criteria.In particular, the proposed GM(1,1;λ) performs better than the traditional GM(1,1) in terms of its accuracy in performing model fitting and forecasting.That is to say, the proposed approach actually improves the model fitting and forecasting precision of the original GM(1,1).

Figure 1 .
Figure 1.Real and forecasting values of different forecasting models for Example 1.

Figure 2 .
Figure 2. Real and forecasting values of different forecasting models for Example 2.

Table 3 .
Fitting and forecasting results of Example 2.