A VVWBO-BVO-based GM (1,1) and its parameter optimization by GRA-IGSA integration algorithm for annual power load forecasting

Annual power load forecasting is not only the premise of formulating reasonable macro power planning, but also an important guarantee for the safety and economic operation of power system. In view of the characteristics of annual power load forecasting, the grey model of GM (1,1) are widely applied. Introducing buffer operator into GM (1,1) to pre-process the historical annual power load data is an approach to improve the forecasting accuracy. To solve the problem of nonadjustable action intensity of traditional weakening buffer operator, variable-weight weakening buffer operator (VWWBO) and background value optimization (BVO) are used to dynamically pre-process the historical annual power load data and a VWWBO-BVO-based GM (1,1) is proposed. To find the optimal value of variable-weight buffer coefficient and background value weight generating coefficient of the proposed model, grey relational analysis (GRA) and improved gravitational search algorithm (IGSA) are integrated and a GRA-IGSA integration algorithm is constructed aiming to maximize the grey relativity between simulating value sequence and actual value sequence. By the adjustable action intensity of buffer operator, the proposed model optimized by GRA-IGSA integration algorithm can obtain a better forecasting accuracy which is demonstrated by the case studies and can provide an optimized solution for annual power load forecasting.


Motivation
Accurate power load forecasting is the key to realize the sustainable development of the power industry and ensure the reliable and safe operation of the power grid. In the past few years, many scholars have given a lot of attention to short-term power load forecasting based on different factors and have reached a higher level. Annual power load forecasting can provide reliable guidance for power grid operation and power construction planning [1][2][3][4]. Affected by PLOS  economic, climate and other factors, an annual load curve often shows a non-linear characteristic and has a certain degree of mutation. As a result, the difficulty of annual power load forecasting is obviously larger than short-term power load prediction.

Literature review
In recent years, different kinds of power load forecasting technology have been developed to improve the forecasting accuracy. The traditional forecasting methods, such as time series analysis [5,6], regression model [7] and auto-regressive moving average model [8], have good forecasting accuracy, but the nonlinear simulating ability is poor. With the development of intelligent technology, many new intelligent methods have been used for annual power load forecasting. Bozkurt et al. [9] introduced two different models for current Turkish Market using Seasonal Autoregressive Integrated Moving Average (SARIMA) and Artificial Neural Network (ANN) and presented an Artificial neural network and SARIMA based models for power load forecasting. Lang et al. [10] proposed a new forecasting model based on the improved neural networks with random weights (INNRW). Fu et al. [11] combined the fruit fly optimization algorithm with the grey neural network and finally improved the prediction accuracy, but this method was easy to fall into local optimal or premature convergence. On the whole, annual power load of an area is growing year by year. However, due to the impact of various external factors, such as economic policy and weather condition, there are a variety of random fluctuations in the growth process of annual power load. According to the grey theory proposed by Deng [12], annual power load forecasting can be regarded as a typical grey system. Therefore, the GM (1,1) model and its optimization model have been widely applied in annual power load prediction and other similar fields [13][14][15][16][17][18][19][20]. Representative researches are as follows. Zhao et al. [13] put forward a new hybrid electricity consumption forecasting method based on GM (1,1) optimized by moth-flame optimization (MFO) algorithm with rolling mechanism. Karadede et al. [18] showed a general natural gas demand forecasting model of the breeder hybrid algorithm based nonlinear regression and applied it to Turkey natural gas demand forecasting to show its superiority and applicability. Boroojeni et al. [19] presented a generalized technique to model historical load data in the form of timeseries with different cycles of seasonality for electric power demand forecasting. In day-ahead natural gas demand predictions, Panapakidis et al. [20] tested the robustness of a novel hybrid computational intelligence model, which combined the Wavelet Transform (WT), Genetic Algorithm (GA), Adaptive Neuro-Fuzzy Inference System (ANFIS) and Feed-Forward Neural Network (FFNN).
If the forecasting model is directly constructed based on original annual power load data, it may lead to a larger forecasting error, or the simulating effect is good but the forecasting effect is poor. For this reason, Liu [21] proposed the concept of buffer operator and established the GM (1,1) model after using the buffer operator to pre-process the historical annual power load data. So the smoothness of the modeling sequence is effectively improved, the randomness of the historical annual power load data is weakened and the overall trend of load development is strengthened, and the accuracy of GM (1,1) model is improved. However, the buffer operator still has problems in the practical load forecasting application. The traditional buffer operator is not adjustable in intensity. According to the qualitative analysis and the subjective judgment, users can only select a fixed buffer operator which is less adaptable, so the buffer effect is sometimes too weak and sometimes too strong. If the buffer effect is too weak, the pre-process effect will be not obvious and the forecasting accuracy of the model will be limited. If the buffer effect is too strong, the inherent variation of the historical annual power load sequence will be changed. Although high accuracy simulating of the pre-process load sequence can be achieved, the forecasting result will actually deviate from the initial annual power load forecasting problem.

Our contribution
Aiming at the above practical application problems of traditional GM (1,1) forecasting model in annual power load forecasting, variable-weight weakening buffer operator (VWWBO) and background value optimization (BVO) are introduced into GM (1,1) and a VWWBO-BVObased GM (1,1) is proposed. Next, grey relational analysis (GRA) and improved gravitational search algorithm (IGSA) are integrated for the parameter optimization of proposed model. By the adjustable action intensity of weakening buffer operator, the dynamic pre-process of historical annual power load data is realized and the forecasting accuracy is improved. Meanwhile, the inherent trend and change rule of the original power load sequence are maintained to the maximum extent in the forecasting and the stability of simulating and forecasting are improved.

Organization of paper
The rest of this paper is organized as follows. Details of the methods, which we used to create our model, are given in Methods. Case studies and detailed discussion on experimental results are presented in Experimental results. Finally, we conclude the paper and provide guidelines for future work in Conclusions.

Basic knowledge of the traditional GM (1,1) forecasting model
Based on the known or uncertain information of past and present, grey theory [12,21] is aiming to build a grey model from past to future, so as to determine the trend of future development of the system. Any random process is considered to be a grey amount that changes within a certain range in grey theory. By accumulating the irregular original power data sequence, a new sequence of randomness weakening and regularity strengthening will be generated. The grey model of GM (1,1) is a differential equation model based on the generated sequence. The specific implementation steps of GM (1,1) are given below.
Step 2. For x (1) (t), its one-order differential equation is constructed as follows: where a and u are undetermined coefficients and called the development coefficient and the grey amount of action, respectively. The effective interval of a is (-2, 2).
The matrix consisting of a and u is defined as the grey parameter as follows: x (1) (t) can be calculated just by calculating the value of a and u, and then the forecasted value of x (0) can be calculated. For x (1) , the mean generation vector B and the constant term vector Y are constructed as follows: The least square method is used to solve the grey parameterÂ, that is: whereâ andû are the approximate value of a and u, respectively.
Step 5. The established GM (1,1) model is used for forecasting as follows:

The proposal of VWWBO-BVO-based GM (1,1)
Generally, the grey model of GM (1,1) has many advantages, such as less load data, no need to consider distribution rules, easy operation and test. However, there are some limitations in its application of annual power load forecasting. In the original power data sequence of load forecasting, the former part is growing faster, but the latter part is growing slowly. It's not a strict regular sequence. The original power data sequence cannot correctly reflect the real change rule because of the impact disturbance. By the effect of the weakening buffer operator on the original power data sequence, the effect of the exception value is weakened, and the change trend of the original data is strengthened. In addition, the weakening buffer operator satisfies the principle of full utilization of information and the principle of "new information priority", that is, the new data sequence is generated by the buffer operator on the basis of the information in the existing data sequence, and the latest information remains unchanged under the action of the buffer operator.
Therefore, using the weakening buffer operator on the original power data sequence can weaken the impact of interference data and strengthen the function of the latest information. As a result, the problem of inconsistent quantitative forecasting result and qualitative analysis conclusion can be effectively solved [21].
The action intensity of the traditional buffer operator is nonadjustable and the fixed buffer operator is less adaptable. This will lead to an uncertain buffer effect, so we introduce VWWBO into the traditional GM (1,1) forecasting model for the dynamic pre-process of historical annual power data.
The original power load data sequence is x = (x(1), x(2),. . ., x(q)) and the variable-weight weakening buffer operator is D. After the pre-processing by D, x is converted to y as follows: where Whenever x is a monotone increasing sequence, a monotone reducing sequence or an oscillating sequence, D is always a weakening buffer operator [15]. Here, λ is the variable-weight buffer coefficient.
The average change rate from x(k) to x(q) in x is represented as r(k): The adjustment degree of D at k is represented as δ(k): The adjustment degree of D at each point is a constant which equals to the variable-weight weakening buffer coefficient λ (i.e., δ(k) = λ) [15]. The adjustment degree δ(k) reflects the buffer operator's effect on the original power load data sequence. By dynamically adjusting the variable-weight buffer coefficient λ, it can realize the fine-tuning of the buffer operator's effect and effectively solve the problem of traditional fixed buffer operator. As a result, the flexibility, controllability and adaptability of the buffer operator in the pre-processing of power load data sequence have been enhanced.
In order to lighten the fluctuation of original power load data sequence, strengthen the trend of modeling sequence and improve the forecasting accuracy, VWWBO and BVO [17] are used to improve the traditional GM (1,1) forecasting model. So, the VWWBO-BVO-based GM (1,1) is built as follows.

GRA-IGSA integration algorithm for model parameter optimization
GSA algorithm. The key factors that affect the forecasting accuracy of the above VWWBO-BVO-based GM (1,1) are the variable-weight weakening buffer coefficient λ and the background-value weight generating coefficient η. A GRA-IGSA integration algorithm is proposed to find the optimal value of λ and η.
In 2009, Rashedi et al. [22] proposed the gravitation search algorithm (GSA). GSA is a new meta-heuristic algorithm and has been used to solve nonlinear benchmark function problem [23] and nonlinear filtering modeling problem [24]. In GSA [25][26][27], each object has the position and velocity which represent the solution. The mass of an object is determined by its position and fitness value. An object with a better fitness value has a greater mass. From the gravitational motion mechanism shown in Fig 1, an object with a greater mass will generate a greater gravitation. This will lead to the motion of small mass objects toward large mass objects, so the information interaction among objects and population evolution is realized. In  Fig 1, F, M and a are the symbols of force, mass and acceleration respectively.
For a system consisting of N objects, in gravitational search, the position P i and the velocity V i of the object i (i = 1,2,. . .,N) are defined as: where L is the dimension number of search space, p l i and v l i are the position and velocity of the object i in the dimension l and l = 1,2,. . .,L.
On a specific time T, the hypotheses are as follows: • The mass of the object i and the object j are M i (T) and M j (T), respectively.
• The gravitational constant is G(T) and the Euclidean distance between the object i and the object j is R ij (T).
• The positions of the object i and the object j are P i (T) and P j (T), respectively.
• The positions of the object i and the object j in dimension l are p l i ðTÞ and p l j ðTÞ, respectively.
Therefore, the gravitation acting on the object i by the object j on the time T can be defined as: where the component of F ij (T) in the dimension l (l = 1,2,. . .,L) is F l ij ðTÞ ¼ GðTÞ ðp l j ðTÞ À p l i ðTÞÞ. Here, ξ is a small constant and R ij (T) = kP i (T),P j (T)k 2 . To increase the random characteristics of the algorithm, the resultant of forces acting on the object i in the dimension l is defined as follows: where rand j is a random variable in the interval [0, 1]. The acceleration of the object i in the dimension l on the time T is as follows: Therefore, the position and the velocity of the object i in the dimension l on the next time T+1 are as follows: where rand i is a random variable in the interval [0, 1]. (20), and then normalized by Eq (21):

The mass of each object is computed first by the fitness function values shown in Eq
where f P i ðTÞ is the fitness value of the object i on the time T, f best (T) and f worst (T) are the best and worst fitness value on the time T, respectively. Gravitational constant, which plays a very important role in GSA, is defined as a function of the time T: where G 0 is the initial value, s is a constant, I is the number of iterations. From Eq (22), we can see that the value of gravitational constant decreases as the number of iterations increases. Thus the search stride is reduced and the search accuracy is improved. IGSA algorithm based on local search. GSA has a fast convergence rate, but it is easy to fall into local optimum. In order to develop the local search capability, a local search algorithm is introduced to increase the chance of better solution and improve GSA. In the local search algorithm, single exchange operator and double exchange operator are used to generate a new neighborhood structure and then new solution is obtained. The exchange operation is shown in Fig 2. In single exchange operation shown in Fig 2(a), the positions of object i in two random dimensions are chosen and exchanged. In double exchange operation shown in Fig 2(b), single exchange operator is done twice.
The local search algorithm that incorporates the exchange operators follows below steps: Step 1. Begin with i = 1, P i ¼ ðp 1 i ; p 2 i ; . . . ; p L i Þ represents the solution of the object i. Step Step 6. If i reaches N, finish the local search. Otherwise, i = i+1 and turn Step 2.

GRA-IGSA integration algorithm.
When constructing the fitness function, the usual method is to minimize the average simulating error. This method is easy to fall into local optimum and leads to higher simulating accuracy and lower forecasting accuracy. In the proposed VWWBO-BVO-based GM (1,1), the weakened data sequence after pre-processing is used for forecasting. If minimizing the average simulating error of weakened data sequence is taken as the optimization target, the forecasting result will be highly fitted to the weakening data sequence and deviate from the historical annual power load forecasting problem. To avoid these problems, we build the fitness function with the optimization target of maximizing the grey relativity between simulating value and actual value based on GRA.
The original power load sequence is x (0) = (x (0) (1),x (0) (2),. . .,x (0) (q)) and the simulating power load sequence based on x isŷ ð0Þ ¼ ðŷ ð0Þ ð1Þ;ŷ ð0Þ ð2Þ; . . . ;ŷ ð0Þ ðqÞÞ. The grey relativity coefficient between x (0) (k) andŷ ð0Þ ðkÞ is as follows: Grey relativity reflects the degree of similarity of sequence curves. A bigger grey relativity between simulating power load sequence and original power load sequence shows that simulating power load sequence can better maintain the intrinsic variation of original power load sequence better. The comprehensive simulating and forecasting effect of the proposed model is better. Aiming to maximize the grey relativity between simulating power load sequence and original power load sequence, the fitness function is defined as follows: gðx ð0Þ ðkÞ;ŷ ð0Þ ðkÞÞ q 8 > > > < > > > : Therefore, a GRA-IGSA integration algorithm is proposed. Combining the global optimization ability of GSA and the space expansion ability of local search operator, IGSA algorithm based on local search is formed. Its pseudo codes are as follows: While T<I Calculate the fitness function value of each object; Calculate the mass of each object by Eqs (20) and (21); Calculate the acceleration of each object by Eq (18); Update the velocity of each object by Eq (19); Update the mass of each object by Eq (19); Do the local search; T+1;

Example analysis
To evaluate the effectiveness of the proposed forecasting model in annual power load forecasting, load data from Literature [28,29] are used for validation. Traditional GM (1,1) model (model 1), GM (1,1) model optimized by buffer operator and time response function (model 2) [20] and the proposed forecasting model (proposed model) are compared and analyzed.
As shown in Table 1, the power consumption data of an area of China from 1998 to 2005 are taken as historical annual power load data for the load forecasting from 2006 to 2008. From Table 1, we can see that the annual power load data have the trend of increasing year by year and fluctuate with economic, climatic and other factors.
The above three models are applied for annual power load forecasting, and percentage error (PE), mean absolute percentage error (MAPE) and grey relativity between simulating value and actual value are used as the indicators to evaluate the above three models. The results are as follows.
The simulating results (1998)(1999)(2000)(2001)(2002)(2003)(2004)(2005) of three models are shown in Table 2 and Fig 3. According to Table 2, MAPE of simulating value relative to actual value of three models are obtained as 3.0218%, 28.1078% and 7.3237%, and the grey relativity of simulating value relative to actual value of three models are obtained as 0.6438, 0.5159 and 0.8662. Therefore, we get  Table 3 and Fig 7. According to Table 3, MAPE of forecasting value relative to actual value of three models are obtained as 12.4555%, 3.3665% and 0.6657%. Therefore, we get the comparisons of PE and MAPE of forecasting value relative to actual value of three models in the stage of 2006-2008 as shown in Figs 8 and 9, respectively. Based on Figs 3-9, model 1, model 2 and proposed model are compared and analyzed as follows:  • For model 1, MAPE of simulating value relative to actual value of three models (1998-2005) is 3.0218% (in Fig 5), so the simulating effect is good. However, it is a failure to control and weaken the fluctuations in historical annual power load sequence (actual value in Table 1). The forecasting effect is poor and MAPE of forecasting value relative to actual value (2006-2008) reaches 12.4555% (in Fig 9).
• For model 2, mean weakening buffer operator is used to pre-process the historical annual power load sequence and GM (1,1) model optimized by time response function is built. High accuracy simulating to weakened value is realized as shown in Fig 3, but MAPE of simulating value relative to actual value (1998-2005) reaches 28.1078% (in Fig 5) and the grey relativity is only 0.5159 (in Fig 6). This illustrates that the buffering effect of mean weakening buffer operator is too strong. The main trend and fluctuation law of weakened value sequence and forecasting value sequence will be completely different from historical annual power load sequence. It can be considered that the high accuracy simulating and forecasting object of model 2 is just the weakened value sequence after pre-processing. In fact, it deviates from the intrinsic annual power load forecasting problem.
• For proposed model, IGSA is used to maximize the grey relativity between simulating value sequence and actual value sequence. The optimal value of variable-weight weakening buffer coefficient λ and background-value weight generating coefficient η are found to do the pre-processing of the historical annual power load data, and then the VWWBO-BVObased GM (1,1) is built for annual power load forecasting. The basic parameters of IGSA are set as: the initial value of gravitational constant G 0 = 10, the constant l = 1.0, the number of objects N = 40, the number of iterations T = 300. The optimal value are found as λ = 0.1136 and η = 0.7659 by GRA-IGSA integration algorithm. By proposed model, MAPE of simulating value relative to actual value (1998-2005) is 7.3237% (in Fig 5) and the grey relativity of simulating value relative to actual value reaches 0.8662 (in Fig 6). MAPE of forecasting value relative to actual value (2006-2008) is 0.6657%, which is the lowest in three models (in Fig 9).
Compared with model 1, proposed model uses buffer operator for historical annual power load sequence pre-processing, weakens the random fluctuation of historical annual power load sequence, and strengthens the main development trend. Although the simulating accuracy is slightly worse, a high forecasting accuracy is achieved. Compared with model 2, proposed model aims to maximize the grey relativity of simulating value relative to actual value and uses GRA-IGSA integration algorithm to optimize the selection of weakening buffer operator. The problem of nonadjustable buffering effect has been solved. While using the buffer operator to pre-process the historical annual power load sequence to improve the forecasting accuracy, the inherent trend and the fluctuation law of historical annual power load sequence are maintained to the maximum extent. The simulating accuracy is improved by an order of magnitude (in Figs 4 and 5). Therefore, the proposed model is more consistent with the historical annual power load forecasting problem and has more practical significance.

Practical application
Using the nationwide power consumption of China (2003)(2004)(2005)(2006)(2007)(2008) [29,30] as the historical annual power load data, we forecast the nationwide power consumption in 2009-2011. In  2003-2007, China shook off the impact of the 1998 financial crisis. The economy was booming and the nationwide power consumption increased by an average of 13.16%. In 2008, the emergence of major natural disasters and the 2008 financial crisis had a great impact on power consumption. Although power consumption continued to grow, the growth rate has dropped considerably. Overall, annual power consumption data in 2003-2008 years shows the main trend of sustained growth. Due to economic factors, natural disasters and other effects, it also shows some fluctuations. Therefore, using the proposed model for annual power load forecasting is suitable. The basic parameters of IGSA are set as same as Section Example analysis. The simulating results and forecasting of model 1 and proposed model are shown in Table 4 and Fig 10. According to Table 3, MAPE of forecasting value relative to actual value of three models are obtained as 12.4555%, 3.3665% and 0.6657%. Therefore, we get the comparisons of PE and MAPE of forecasting value relative to actual value of three models in the stage of 2006-2008 as shown in Figs 8 and 9, respectively.
According to Table 4  and forecasting results are stable. Therefore, the proposed model can meet the requirements of annual power load forecasting.

Conclusions
The main findings of this study are summarized as follows.
1. The proposed model introduces variable-weight weakening buffer operator to dynamically pre-process the original power load data and realizes the action intensity adjustment of buffer operator by variable weight buffer coefficient. It effectively solves the problem of nonadjustable action intensity of traditional buffer operator, enhances the adaptability and improves the forecasting precision.
2. Based on grey relational analysis and improved gravitational search algorithm, the proposed model aims to maximize the grey relativity between simulation value sequence and actual value sequence for parameter optimization.
3. The intrinsic variation law of the historical annual power load data is maintained to the maximum extent. It realizes the organic unification of data pre-processing and historical annual power load forecasting, so the simulation and forecast effect are more stable and the forecast accuracy is higher. Author Contributions