MONTHLY NATURAL GAS CONSUMPTION’S MODELLING AND ITS TREND ANALYSIS FOR YOZGAT IN TURKEY

In this study, an energy consumption modelling for long term (December 2006March 2016) forecasting of monthly natural gas consumption in households and industry area for Yozgat city, Turkey was presented. In this context, it can be said that this paper has two purposes. One of them is the application and accuracy of the artificial neural networks. Estimate performances are compared with each other, and the estimates of the optimal models are evaluated with the monthly recorded natural gas consumption according to root mean square error, mean absolute error, and correlation coefficient. The other purpose of the study is to analysis trend of monthly natural gas consumption of Yozgat by using Mann-Kendall and a new method recently proposed by Şen. The results showed that the artificial neural networks gave satisfactory results in estimating monthly natural gas consumption. In the trend analysis, it was seen that both Mann-Kendall and Şen trend tests gave statistically significant increasing trend at 95% confidence level for monthly natural gas consumption of Yozgat.


INTRODUCTION
Most of Turkey's natural gas demand (roughly 57%) comes from Russia (MENR, 2008; Erdogdu, 2010a). Natural gas consumption (NGC) started with 0.5 bcm (billion cubic meters) in 1987 (Ankara, Bursa, and Istanbul provinces), reached approximately 35 bcm in 2007, and was about 50 bcm in 2015, and have been used in most regions of Turkey. Moreover, future consumption amount in Turkey is predicted as 82 bcm for 2020. When it is looked these high consumption amounts, it can be said that although Turkey has very limited natural gas reserves, NGC have been exponentially rising in households, industry, power generation, refinery, commercial and special uses (Aras and Aras, 2004;Kilic, 2006;Hacisalihoglu, 2008;MENR, 2008;Erdogdu, 2010aErdogdu, , 2010bCapik et al., 2013;Melikoğlu, 2013). But, NGC in a number of areas is decreasing according to the Energy Market Regulatory Authority (EMRA) in 2015 and 2016 (EMRA, 2017) in Turkey. With the Board decision dated 26/01/2017 and numbered 6884, the estimated amount is also decreasing for 2017 when it is compared with the previous years (http://www.epdk.org.tr/tr/anasayfa). On the other hand, Turkey is seen as an energy corridor both Europe and Asia because of construction of natural gas pipelines, distribution of natural gas networks, and establishment of natural gas market. In this concept, Turkey has a critical international advantage which is very significant for political ground, technical works and economy of the country. For example, there are some projects which have been under construction, recently. Transanatolian Pipeline (TANAP) is one of them, and is very important for natural gas transfer to the other regions especially to Europe, and is about 4000 km long, and is the biggest project in the world (http://www.tanap.com). TurkStream, from Russia to Turkey, south and south-east of Europe, is also one of the important projects in national and international platform. This project is also under the construction (http://turkstream.info). And, there are other projects such as Israel, Iraq and Turkmenistan natural gas pipelines in feasibility phases.
Studies in modelling of NGC in the literature were investigated by Soldo, 2012. This study is an extensive review paper. Moreover, some mathematical methods (Liu and Lin, 1991;Brabec et al., 2008;Vondráček et al., 2008;Vitullo et al., 2009;Sabo et al., 2011;Goncu, 2013) were also implemented in modelling NGC. In Turkey, there are some relevant studies (Gumrah et al., 2001;Ozturk and Hepbasli. 2003;Sarak and Satman, 2003;Haldenbilen and Ceylan, 2005;Ediger and Akar, 2007;Goncu, 2013) about future energy demand and modelling. In all studies' authors have explained that Turkey's energy supply and demand have increased dramatically over the last three decades as result of industrialization, urbanization, economic expansion rate, population growth, social prosperity and other factors. In this context, studies in modelling of NGC are very limited. For example, Toksari (2010) predicted the natural gas demand based on economic indicators: gross domestic product (GDP), population, import, and export variables for Turkey by using linear and quadratic forms. Demirel et. al. (2012) used the multiple regression, the autoregressive moving average with exogenous input (ARMAX) model and artificial neural networks methods to model natural gas consumption recorded at Istanbul. Melikoğlu (2013) investigated Turkey's natural gas demand for 2030, and calculated as 76.8 billion m 3 using the linear model and 83.8 billion m 3 based on the logistic model. He also compared his study with the previous studies and Official BOTAS for natural gas demand forecasts. His study was in a better agreement with the BOTAS's official forecasts than the models by Kilic (2006) and Erdogdu (2010a).
As well as modelling studies, trend analysis is one of the most important subjects in this area. In order to determine whether an increasing or descending trend exists for the incident is very important for distribution, transportation companies and government agencies related to the natural gas sector. In this context; low, medium and high values of a variable are very important issue. Moreover, these values are used to decide various design parameters based on scientific aspects and real applications everywhere in the world. Some hydrological studies (Douglas et al., 2000;Sang et al., 2014;Şen, 2015;Ay and Kişi, 2016;Ay, 2016a;Ay, 2017) can be given for this purposes. In this study, a new trend method recently proposed by Şen was firstly used for monthly NGC data recorded at Yozgat. The Mann-Kendall (MK) trend test was also applied to the same data. To my knowledge, there has not been any study about trend analysis of monthly NGC of Yozgat in Turkey; thus, this is the first trend analysis study of the monthly NGC.
The main point of this article is that estimations of natural gas demand should be done correctly in planning phase. In this study, two aims were implemented. It was firstly tried to model and trend analysis of monthly NCG for Yozgat. The results showed that the artificial neural networks gave satisfactory results in estimating monthly natural gas consumption. And, the results of the trend analysis were also compared with the modelling results.

YOZGAT, MONTHLY NATURAL GAS CONSUMPTION AND MONTHLY MEAN AIR TEMPERATURE DATA
Yozgat province is in Middle-Kizilirmak basin and is the border of the Kizilirmak River starting from Zara/Sivas province in the Central Anatolia, Bozok Plateau. The province in Figure 1 is nearly at 34º05'-36º10' eastern longitudes and 38º40'-40º18' northern latitudes. The number of people was determined as 421041 people according to the 2016 year's survey. Generally, livelihood of the province is agriculture and livestock, and industry has mediumsized structure. 14037 km 2 area of Yozgat is the 15 th largest city in Turkey in terms of geographical area. The city center that has the rugged terrain is approximately 1298 m elevation (height). The annual average rainfall is about 400-600 mm and this value is smaller than annual average rainfall (643 mm/year) of Turkey. Many irrigation, drinking water, flood protection, energy structures and dams were also constructed, and there are some projects on under construction in the province (Ay, 2016b). Sürmeli Gas Distribution Company is the first and only gas supply company in Yozgat, and has supplied natural gas to the province since November 17, 2007. There is no natural gas pipeline everywhere of Yozgat. Moreover, there are few industries using natural gas. Generally, most of the total natural gas has been used in households.
The recorded monthly NGC data in household and industry areas between December 2006 and March 2016 were used to implement this study. In Table 1, it can be seen the basic statistics of the monthly recorded NGC and monthly mean air temperature data (Station no.: 17140) between December 2006 and March 2016. The 17140-station is in Yozgat's center. Because one of the most vital factors which affect the NGC is air temperature, this variable is commented with natural gas in this study. In Table 1, it can be said that although monthly NGC data has a positive skewness, the air temperature has a negative skewness. Skewness coefficient of the monthly air temperature is too close zero. This means that air temperature value is close to normal (Gaussian) distribution. Moreover, NGC is higher variation coefficient value than the air temperature variable. It can be said that NGC is more scattered and more unsteady than the air temperature variable.

Table 1. Basic statistics of the monthly recorded NGC and monthly mean air temperature of Yozgat in Turkey
Xmax, Xmin, Xmean, Sd, Cvx, and Csx indicate the maximum, minimum, mean, standard deviation, variation coefficient, and skewness coefficient, respectively.
In Figure 2, it can be seen the time series of Yozgat's monthly NGC, and it is clearly seen that the NGC is exponentially increasing way. Relationship between the monthly NGC and monthly mean air temperature data recorded (December 2006 and March 2016) by Turkish State Meteorology Service is given in Figure 3. It can be clearly seen that although air temperature is getting higher; the consumption natural gas is getting lower. Therefore, there is an opposite relationship between two variables. The quantity, r, called the linear correlation coefficient in Figure 3, measures the strength and the direction of a linear relationship between two variables. The correlation coefficient (r) for both variables is calculated to be -0.68 (negative correlation). This means that it can be said that determination coefficient (R 2 ) of the relationship is lower than average value (0.50). It can be said that there are other variables effecting this relationship.  Table  2. Although the most increasing rate is clearly seen at 2008 with +137,96%, the least increasing rate is seen at 2013 with +1,10%. Therefore, it can be talked about that the increasing in NGC always occurs for all the years in Yozgat.

Multi-layer perceptron (MLP)
This method has the architecture of the network that is most commonly used with the backpropagation algorithm the multi-layer feedforward network. It is a massive parallel system composed of many processing elements connected by links of variable weights (Haykin 1998). An elementary neuron with R inputs is determined. Each input is weighted with an appropriate w. The sum of the weighted inputs and the bias forms the input to the transfer function f. Neurons can use any differentiable transfer function f to generate their output. The initial assigned weights are progressively corrected during the training process. In this process, the outputs predicted by MLP are compared with known outputs, and errors are back propagated (from right to left) to determine the appropriate weight adjustments necessary to minimize errors. In this study, the Levenberg-Marquardt (Marquardt, 1963) algorithm is used for adjusting the MLP weights. In backpropagation it is important to be able to calculate the derivatives of any transfer functions used. Each of the transfer functions above, logsig, tansig, and purelin, can be called to calculate its own derivative. X 1 , X 2 , X 3 ,…, X R are the input signals; W 1 , W 2 ,…, W R are the weights of neuron; b is bias value; and f activation function, R number of the elements in input vector. The linear and sigmoid are the most common used activation functions (tansig, logsig, purelin and tribas) in the form of artificial neural networks.

Radial basis neural network (RBNN)
RBNN was first introduced into the ANN literature by Broomhead and Lowe, Poggio and Girosi. The RBNN has two layers whose output nodes form a linear combination of the basis functions. RBNN is also known as a localized receptive field network because of the fact that the basis functions in the hidden layer produce a significant nonzero response to input stimulus only when the input falls within a small localized region of the input space. The RBNN has connection weights between the hidden layer and the output layer only. These weight values can be obtained by linear least-squares method, which gives an important advantage for convergence. Gaussian activation function is widely used as radial basis function. The RBNN method does not perform parameter learning as in the MLP. It performs linear adjustment of the weights for the radial bases. This characteristic gives the RBNN advantage of a very fast converging time without local minima. Because, its error function is always a convex. In this study, different numbers of hidden layer neurons are examined for the RBNN models with a simple trial-and-error method.

Mann-Kendall (MK) trend test
MK test is one of the non-parametric tests to detect trend in a time series (Mann, 1945;Kendall, 1975;Helsel and Hirsch, 2002;Onoz and Bayazit, 2003). The MK test statistic (S) is calculated in the following Equation 1 and 2: where x i and x j are the data values at times i and j, n indicates the length of the data set. While a positive value of S indicates an increasing trend, negative value indicates a decreasing trend. The following expression, as an assumption, is used for the series where the data length n>10 and data is approximately normally distributed with variance (σ=1) and mean (μ=0) value.

Şen trend test
In this test, a recorded hydrological time series is divided into two equal halves from the first date to the end date, and both sub-series are separately sorted in ascending manner. The first sub-series (X i ) is located on X-axis, and the other sub-series (X j ) is located on Y-axis ( Figure 4). If data is collected on the 1:1 (45˚) line, it can be said that there is no trend. If data is in the below triangular area, it can be said that there is a decreasing trend in time series. If data is in the upper triangular area, it can be said that there is an increasing trend in time series (Şen, 2012; 2014). The first ordered data, X i  Şen (2015) proposed to the method a new statistical process. Steps of this method are given in following by formulas 5-10.
̅̅̅̅ ̅̅̅̅ In here, ̅̅̅: mean of the first data, ̅̅̅: mean of the second data, ρ: correlation between first and second data, s: slope value, n: number of data, σ: standard deviation of all data, : slope standard deviation, and s critical denotes Z critical values in one-way hypothesis at %95 confidence level. Critical upper and lower limits values calculated by Equation 7 are established to make limits for hypothesis test. Slope value, s, of each station falls outside the lower and upper confidence limits, and thus, the alternative hypotheses, H 1 , are approved and they indicate the existence of trends (YES) as decision. The type of trend is stated depending on the slope (s) sign. Slope value (s) can be positive or negative. This means that there is an increasing (+) or a decreasing (-) trend in time series (Şen, 2015). This method was applied by Şen (2012) to annual flow data, annual total precipitation, and run-off data. Another study by Şen (2014) was applied to the long-term recorded air temperature data in Turkey. Kisi and Ay (2014) determined the trend of some water quality variables such as dissolved oxygen, chemical oxygen demand. Şen (2015) applied to air temperature, discharge and rainfall data. Ay and Kisi (2017) applied the test to streamflow data. Ay and Özyıldırım (2017) implemented a trend study for monthly total rainfall and monthly mean air temperature variables.

Modelling of monthly natural gas consumption by using artificial neural networks
In order to form all of the methods for estimating monthly NGC, the available 112 recorded data set was divided into two phases as training and testing. Roughly %60 (67) of data set was chosen for training set, %40 (45) of data was chosen for testing phase respectively. Estimations of Yozgat's monthly natural gas consumption (NGC) based on simulated MLP and RBNN approaches were applied to model NGC's a month ahead. For learning process for MLP method, the input vector (NGC (t) ) and corresponding target vector (NGC (t+1) ) were used to train the neural networks by applying LM algorithm which is both fast and has the advantage in terms of gaining the time according to other training algorithms. Because the number of hidden units directly affects the performance of the network, many experimental investigations are conducted with models. It was found that optimum number of hidden layer as "1" in all of the models. The stopped criteria for the training phase is: MSE=0.00001 or a number of epochs is 100. It was found that the optimum activation functions of MLP models were logsig for input layers, and for output layer was purelin. Unlike the MLP, the RBNN can have optimum spread and hidden node numbers were found trial-error method. After the optimum spread coefficient and hidden nodes number were found, main program was executed with these numbers.
Formulas given by equations 11-14 were used for the evaluation of the model's performances in the study. The correlation coefficient (r), determination coefficient (R 2 ), root mean square error (RMSE), and mean absolute error (MAE) criteria are expressed as: (14) where N denotes the number of observations; Y indicates the modelled and recorded NGC, and Y indicates mean of the recorded monthly NGC. Table 3 shows the results of MLP and RBNN models in terms of RMSE, MAE and r values in the training and testing phase. In this table, the MLP and RBNN models of one input parameter NGC (t) are used to model NGC (t+1) . Similarly, MLP(1,1,1) model has slightly better performance than the RBNN(1,0.3,3,1) model according to RMSE, MAE and r criteria in training and testing phases. Therefore, the RBNN(1,0.3,3,1) model has higher errors values than the MLP(1,1,1).

Table 3. RMSE, MAE and r values of the MLP and RBNN in training and testing phases for monthly natural gas consumption (NGC) of Yozgat
The models are also graphically compared in Figure 5 in the form of time series and scatter plot in which there are the recorded data in x axis and the modelled data in y axis for NGC of Yozgat by (a) MLP(1,1,1) and (b) RBNN(1,0.3,3,1) models in the testing phase. It can be seen from Table 3 and figures (a and b) the best model is MLP according to the evaluation criteria. It can be clearly seen that MLP gives satisfactory performance in estimating monthly NGC (t+1) of Yozgat. Moreover, both methods give both overestimations and underestimations in some points.   (1,1,1) and (b) RBNN(1,0.3,3,1) models in testing phase (45 data)

Trend analysis of monthly natural gas consumption in Yozgat
Assumptions such as pre-whitening process (von Storch, 1995) were not applied to the data. Original recorded data were taken into consideration in order not to lose originality of the time series in the trend methods (Douglas et Table 4 shows the results of the MK trend test for Yozgat's monthly NGC. Z value of the data was calculated and compared with normal distribution critical Z values at the 90% and 95% two-tailed confidence levels. According to the MK trend test, it can be said that Yozgat's monthly NGC has a statistically significant increasing trend at 95% and 90% confidence levels with Z=+5.86. Therefore, H 1 hypothesis is accepted. Results of Şen trend test are also given in the following Table 5, Figure 6 and Figure 7. Table 5 shows the calculation steps of Şen trend test and it can be seen the there is a statistically significant increasing trend at 95% confidence level. Figure 5 Figure 7. All of the data separate above the 1:1 straight line. In addition, the linear trend equation of NGC's time series is seen in order to describe the behavior of the time series. It is also possible to see that there is an increasing trend in low, medium and high values for Yozgat's monthly NGC.

Figure 6: Time series of Yozgat's monthly NGC for Şen trend test
The variations of the Yozgat's monthly NGC are summarized in Table 6 in terms of low, medium and high values. Şen and MK test results are particularly evaluated in this table. It is clearly seen from this table that both methods have same result which is the statistically significant increasing trend in monthly NGC for Yozgat.

CONCLUSIONS
In this study, two aims were implemented on modelling and trend analysis of natural gas consumption. Estimations of Yozgat's monthly natural gas consumption (NGC) based on simulated MLP and RBNN approaches were applied to model NGC's a month ahead. Trend analysis of NGC was investigated by using MK and a new method recently proposed by Şen trend test. Consequently, the results of both analyses were compared with each other. In this context, the following observations can be drawn from the results of the present study: 1. The developed MLP model gave slightly better estimations than the RBNN method. The MLP and RBNN approaches could be used to estimate both Yozgat and other regions in the world's natural gas consumption. Estimates of the MLP were closer values to the recorded monthly NGC than the RBNN. On the other hand, the MLP and RBNN gave both overestimates and underestimations in some data (see Figure 5a and 5b).
2. The MK and Şen trend tests gave statistically significant increasing trend at 95% confidence level for monthly NGC. The study also showed that the Şen's method compared with the MK trend test has some advantages. One advantage is that it supports to identify trends in low, medium and high values. This helps us in evaluation of the data's extreme values as a figure (see Figure 7).
3. The modelling of monthly NGC can also be investigated with other methods such as neuro-fuzzy, genetic algorithm, ant colony optimization, and wavelet analysis. The results of the different methods can be compared with each other, and it can be found the best model. This may be a subject for another study. 4. NGC may rise in point of number of subscribers and opening new industries branches in different times. On the other hand, air temperature may affect the consumption. Moreover, other some variables may affect the NGC. This issue is to be important for trend analysis and modelling. In Yozgat, available data was used in this study. In this point, it can be said that the increasing trend can be easily seen at Figure 2 and Table 2. This hypothesis was tried and proved with the methods in the study. 5. The results of this study are significant report to engineers, practitioners, policy makers and natural gas companies. Moreover, the results of the models and trend tests can also be evaluated for many purposes such as annual maximum consumption, using areas (households, industry, other places) and the optimal transportation capacity in pipelines.
6. All energy types including nuclear, hydropower, wind, geothermal, solar, biomass, coal, and electricity are in connected with each other as directly or indirectly in the nature. Furthermore, energy is the most crucial issue for all countries. In here, there is a very significant point to be considered. For example, most countries including Turkey have tried in reducing CO 2 emissions against to the increasing levels of greenhouse gases (GHGs) in the atmosphere according to Kyoto Protocol conditions. Among the fossil fuels, the natural gas is the least responsible for CO 2 emissions (Demirbas, 2006); therefore, countries must use the natural gas or other clean energy types to reduce CO 2 emissions in their countries.
7. As a main global issue, nuclear energy stations' construction projects have been thought in Kırklareli, Mersin and Sinop provinces. Other energy projects in Turkey are also in planning, designing, investment phases or already constructed to the related place, and they have been working efficiently. In the end of these energy projects, they are going to absolutely provide decreasing in NGC for some areas; therefore, this situation will reduce Turkey's dependence on imported natural gas sector. This matter has also been mentioned in National Energy Efficiency Action Plan for 2017-2023 (MENR, 2018) in term of reducing external dependence.