Forecasting of Day-Ahead Natural Gas Consumption Demand in Greece Using Adaptive Neuro-Fuzzy Inference System

Abstract

(1) Background: Forecasting of energy consumption demand is a crucial task linked directly with the economy of every country all over the world. Accurate natural gas consumption forecasting allows policy makers to formulate natural gas supply planning and apply the right strategic policies in this direction. In order to develop a real accurate natural gas (NG) prediction model for Greece, we examine the application of neuro-fuzzy models, which have recently shown significant contribution in the energy domain. (2) Methods: The adaptive neuro-fuzzy inference system (ANFIS) is a flexible and easy to use modeling method in the area of soft computing, integrating both neural networks and fuzzy logic principles. The present study aims to develop a proper ANFIS architecture for time series modeling and prediction of day-ahead natural gas demand. (3) Results: An efficient and fast ANFIS architecture is built based on neuro-fuzzy exploration performance for energy demand prediction using historical data of natural gas consumption, achieving a high prediction accuracy. The best performing ANFIS method is also compared with other well-known artificial neural networks (ANNs), soft computing methods such as fuzzy cognitive map (FCM) and their hybrid combination architectures for natural gas prediction, reported in the literature, to further assess its prediction performance. The conducted analysis reveals that the mean absolute percentage error (MAPE) of the proposed ANFIS architecture results is less than 20% in almost all the examined Greek cities, outperforming ANNs, FCMs and their hybrid combination; and (4) Conclusions: The produced results reveal an improved prediction efficacy of the proposed ANFIS-based approach for the examined natural gas case study in Greece, thus providing a fast and efficient tool for utterly accurate predictions of future short-term natural gas demand.

1. Introduction

The increasing technological advancements and the rapid global population growth have led to a remarkable increase in energy consumption all over the world and especially in the developed and developing countries. After all, energy consumption is an index of a society’s economical welfare and represents the economic development of a city or country [1]. Due to this unexpected boost in energy consumption over the past few decades, the need for energy demand management became crucial for achieving economic success that will result in self-sufficiency and economic development [2]. Thus, energy consumption forecasting is essential, as it predicates an energy efficient policy, the optimization of usage [3] and energy supplies optimum management. Even though a variety of methods have been investigated for energy demand forecasting, this is not an easy task, as it is affected by uncertain exogenous factors such as weather, technological development and government policies [4].

Among all energy resources, natural gas (NG) has received the largest increase in consumption lately [5], mostly due to its popularity as a clean energy source, with respect to environmental concerns. In particular, this energy source is characterized by low-level emissions of greenhouse gases in comparison with other non-renewable energy sources [6,7] and is considered as the cleanest-burning fossil fuel [8,9,10,11]. One fifth of the world’s primary energy demand is covered by NG [7] and is linked to industrial production, transportation, health, agricultural output and household use.

According to the sharp increase in NG consumption, the forecasting of NG consumption has attracted great attention since it is essential for project planning, gas imports, tariff design, optimal scheduling of the NG supply system [12], indigenous production, infrastructures planning and cost reduction at different levels [7]. Moreover, some noteworthy factors like the need for distribution planning, especially in residential areas, the increasing demand of NG and the restricted NG network in many countries, make consumption forecasting on an hourly, daily, weekly, monthly or yearly basis highly important [13].

Especially in national energy strategy, NG demand forecasting is of high importance and can help policymakers all over the world to choose certain strategies in this direction. The development of systems that model NG consumption could assist in good government policymaking. Over the past decade, there has also been a significant increase in NG consumption in Greece and so the prediction of demand has become crucial accordingly. In particular, there is a need for further distribution planning in the Greek territory, especially in residential areas and in cases of high demand, when the accumulation ability of the network itself is decreased [13].

Selecting appropriate load forecasting methods is the most important step. Various models have been proposed in the field of energy forecasting so far, and can be divided into three groups: (i) traditional statistical models such as regression analysis, time series methods and ARIMA [14]; (ii) artificial intelligence (AI)-based methods such as wavelet analysis, artificial neural network (ANN) methods, neuro-fuzzy, machine learning, gray theory prediction; and (iii) hybrid models, which are a combination of various forecasting methods [15,16,17,18,19,20,21,22,23]. As regards the forecasting period horizon, prediction of demand is classified into three categories: short-term, medium-term and long-term forecasting. Short-term forecasting is used for hourly, daily and weekly demand predictions and in regard to NG, it is oriented mainly in system management and balancing, storage capacities optimization, as well as in making optimum purchasing and operating decisions [24,25]. Medium-term forecasting deals with seasonal demand prediction (one to several months) and mainly plans the fuel purchases, while long-term forecasts (more than a year ahead) aim to develop the power supply and delivery system [26,27].

ANN, genetic algorithm (GA) and fuzzy inference systems (FIS), as artificial intelligence (AI) techniques, are among those methods that are often used in the energy domain and specifically for energy demand forecasting, due to their high flexibility and reasonable estimation and prediction ability. From the relevant literature, there have been several attempts in energy forecasting by ANNs, like those in [28,29,30,31,32]. ANNs have been applied to forecast electric energy consumption in Saudi Arabia [33], energy consumption of a passive solar building [34], energy consumption of the Canadian residential sector [35], the peak load of Taiwan [36], while ANNs have been further explored in [37,38,39,40,41] for short-term load forecasting. An abductive network machine learning for predicting monthly electric energy consumption in domestic sector of Eastern Saudi Arabia was proposed in [42]. Moreover, support vector machines (SVMs) and genetic algorithms (GA) were explored in [29,30,43] to predict electricity load. In [44], SVMs coupled with empirical mode decomposition were used to perform long term load forecasting. GAs have been used in [45] and [46] to estimate Turkey’s energy demand and electricity demand in the industrial sector, respectively. In addition, Azadeh et al. have proposed the integration of GAs and ANNs to estimate and predict electrical energy consumption [47]. A number of literature reviews have also been performed regarding various energy fields. For example, in [48], a review of the conventional methods and AI methods for electricity consumption forecasting was provided, while in [49] the strengths, shortcomings, and purpose of numerous AI-based approaches in the energy consumption forecasting of urban and rural-level buildings were discussed. In [50], a review about conventional models, including time series models, regression models and gray models was conducted with respect to energy consumption forecasting. Moreover, an overview of AI methods in short term electric load forecasting area was discussed in [51].

1.1. Indicative Related Work on AI Applied in Natural Gas Consumption Forecasting

On the other hand, there are many studies where different AI methods like neural networks, neuro-fuzzy and other ANN topologies have been investigated and applied in NG demand forecasting [52,53,54,55,56,57,58,59,60]. In this context, ANNs were extensively used in [21,54,57,60,61,62,63,64,65,66,67,68,69] to investigate short-term NG forecasts, while in [56] different types of ANN algorithm were explored to forecast gas consumption for residential and commercial consumers in Istanbul, Turkey. ANNs were also used in [70] for the daily and weekly prediction of NG consumption of Siberia, using historical temperature and NG consumption data, in [71] for NG output prediction of USA until 2020, as well as in [72] for the prediction of NG consumption and production in China from 2008 to 2015, applying the grey theory along with NNs. Furthermore, a combination of ANNs was applied in [54,55] for the prediction of NG consumption at a citywide distribution level.

More recent studies presented numerous techniques for NG demand forecasting, including computational intelligence-based models (ANNs), fuzzy logic and support vector machines [12,73,74,75,76]. In this sense, a combination of recurrent neural network and linear regression model was used in [77] to generate forecasts for future gas demand, whereas a multilayered perceptron (MLP) neural network was deployed in [78] to estimate the next day gas consumption. A day-ahead forecast was also examined in [79] by developing a functional autoregressive model with exogenous variables (FARX). Moreover, machine learning tools such as multiple linear regression (MLR), ANN and support vector regression (SVR) were devised in [80] to project NG consumption in the province of Istanbul, as well as in [77] to forecast the residential NG demand in the city of Ljubljana, Slovenia. When considering AI methods, self-adapting intelligent grey models were also deployed for forecasting NG demand, as in [81,82].

Regarding neural network algorithms, the multilayer perceptron and the radial basis function network with different activation functions were trained and tested in [21,63,65], while the authors of [66] used a multilayer perceptron algorithm for neural network and compared this model with two time series models. In addition, Taspinar et al. explored daily gas consumption forecasting through different methods including the seasonal autoregressive integrated moving average model with exogenous inputs (SARIMAX), multi-layer perceptron ANN (ANN-MLP), ANN with radial basis functions (ANN-RBF), and multivariate ordinary least squares (OLS) [63]. Different sets of AI methods were also implemented in the following studies regarding NG consumption. ANNs with linear regression models were used in [83] for daily prediction, while ANN and Fuzzy ANN models were investigated in [84] regarding consumption in a certain region of Poland. Finally, it is worth mentioning that there has been similar research in [12], which proposed the hybrid wavelet-ANFIS/NN model to compute day-ahead forecasts for 40 distribution nodes in the national NG system of Greece.

1.2. Related Work on ANFIS in Energy Consumption Forecasting

There are plenty of works in the relevant literature regarding the application of the ANFIS model in forecasting energy consumption. As regards the electricity domain, ANFIS model was applied to forecast annual regional load in Taiwan [85] and annual demand in Turkey [86], showing in both cases that the results are good and the ANFIS model performed better than regression, neural network and fuzzy hybrid systems. ANFIS was also used in [87] for short-term electricity demand forecasting, using weekly electricity load data, as well as in [88], to estimate possible improvement of electricity consumption. Also, for electricity load forecasting, ANFIS was used in [89] to highlight its superiority to the ANN model, while it was furthermore applied in the field of transportation, forecasting the corresponding energy demand for the years 2010 to 2030, in the country of Jordan, revealing the efficiency of the examined model. Another study regarding the energy domain, where the ANFIS model was applied, is that of [90]. A long-term prediction of oil consumption was studied, further examining the interrelationship between oil consumption and economic growth in Turkey, for the years 2012 to 2030.

1.3. Related Work on ANFIS in Natural Gas Consumption Forecasting

Casting a view on the literature that refers to NG consumption forecast, the authors came across only one study that devised solely an ANFIS model. Specifically, ANFIS was used in [91] in order to estimate the daily NG demand in Iran, which actually used an extremely small dataset of historical data for both testing and training (December 2007–June 2008). Models trained on a small dataset tend to overfit, which results in high variance and very high error on a test set, producing inaccurate results. In this case, the predicting error decreases monotonically with the size of training set [92]. The rest of the studies dealt with approaches that combine ANFIS with other methods. For example, in [21], statistical time series analysis along with ANN and ANFIS methods were applied in order to predict weekly NG consumption in Turkey. Moreover, an ANFIS-fuzzy data envelopment analysis (FDEA) was developed in [93] for long-term NG consumption forecasting and analysis. In this study, 104 ANFIS were constructed and tested and six models were proposed to forecast annual NG consumption. The same approach was proposed in [94] for accurate gas consumption estimation in South America with noisy inputs. An ANFIS-stochastic frontier analysis (ANFIS-SFA) approach was formulated in [95] for long-term NG consumption prediction and analysis. Three patterns of the hybrid ARIMA–ANFIS model were tested in [2] to predict the annual energy consumption in Iran, using a set of data like population, GDP, export and import. Finally, a hybrid model of adaptive neuro fuzzy inference system and computer simulation for the prediction of NG consumption was developed in [96].

1.4. Related Work on Fuzzy Cognitive Maps (FCMs) in Energy and Natural Gas Consumption Forecasting

Moreover, other soft computing techniques, like evolutionary fuzzy cognitive maps (FCMs) have been applied for the modeling and prediction of time series problems. The dynamic modeling structure of FCMs inheriting the learning capabilities of recurrent neural networks works properly for modeling and time series prediction. Salmeron and Froelich further investigated the applicability of FCMs in univariate time series prediction by proposing an FCM simplification approach with the removal of nodes and weights [97]. Regarding the task of multivariate time series prediction, Froelich and Salmeron proposed a nonlinear predictive model based on an evolutionary algorithm for learning fuzzy grey cognitive maps [98], while Papageorgiou et al. [99] and Poczeta et al. [100] applied a new type of evolutionary FCM enhanced with the Structure Optimization Genetic Algorithm (SOGA) in energy for electricity load forecasting. Through the SOGA algorithm, an FCM model can automatically be constructed by taking into consideration any available historical data. A two-stage prediction model for multivariate time series prediction, based on the efficient capabilities of evolutionary FCMs and enhanced by structure optimization algorithms and ANNs, was introduced in [101]. In the first stage of the prediction model, SOGA-FCM was applied for selecting the most significant concepts and defining the relationships between them. Next, that model was fed into the second stage to define the initial features and weights of the training ANN. This generic prediction approach was applied in four common prediction problems, one of which dealt with electric power consumption.

In [102], Poczeta and Papageorgiou conducted a preliminary study on implementing FCMs with ANNs for NG prediction, showing for the first time the capabilities of evolutionary FCMs in this domain. Furthermore, the research team in [13] recently contacted a study for time series analysis devoted to NG demand prediction in three Greek cities, implementing an efficient ensemble forecasting approach through combining ANN, real coded genetic algorithm (RCGA)-FCM, SOGA-FCM, and hybrid FCM-ANN. In this research study, the advantageous features of intelligent methods through an ensemble to multivariate time series prediction in NG demand forecasting are explored.

1.5. Research Gap and the Novelty of This Study

Based on the reported literature survey and reviews [7,15], regarding the application of ANN-based and hybrid forecasting methods, a research gap has been identified in the field of NG day-ahead demand prediction. The observed gap mainly refers to the lack of model simplicity and flexibility, the insufficient exploration of certain modelling aspects, and inadequacy to cope with the inherent fuzziness in data handling. Most of these forecasting methods need a large dataset to be trained and a relatively large number of features to make accurate predictions. Furthermore, they are complex in their structure, time consuming and difficult to be used by non-experienced AI users. There has been hardly any research on successfully applying the ANFIS technique on the field of NG demand prediction, having performed a deep exploration process for determining the best model configuration, thus producing a highly accurate model with generalization capabilities.

Considering the aforementioned limitations, this work aims to fill the observed research gap and seeks to develop an easy to use, robust and flexible ANFIS model, which is at the same time fast, simple in structure and able to cope with fuzziness. More specifically, the proposed ANFIS architecture uses as model’s inputs the most important and commonly used input variables according to the literature [7,15], such as day, month and daily average temperature, along with past NG consumption data. Moreover, a relatively large dataset was used for both testing and training the model, resulting in a systematic improvement of the model’s predictive accuracy [92]. The current work pays great attention to the generalization of the proposed method and tries to properly evaluate the model’s generalization capabilities in two ways: (i) by applying the model on a city level, where 10 different cities were properly examined, and (ii) by carrying out an exploration process, where 94 different models’ configuration sets were examined for each one of the cities that participated in this research.

To sum up, the innovations offered in this paper are as follows, highlighting the contribution of this work to the research community:

• The creation and demonstration of a simple, fast, robust ANFIS prediction tool to forecast NG demand using historical time series data. The proposed model is characterized by high flexibility, especially in large datasets, easiness of use and low execution time requirements.
• The rigorous ANFIS fine-tuning for determining the most appropriate architecture for an enhanced prediction performance.

1.6. Aim of This Research Work

The motivation of this work is to propose an ANFIS-based forecasting approach with generalization capabilities for short-term (day-ahead) city-level NG prediction in Greek areas. Also, a comparative analysis is conducted, applying ANNs, evolutionary FCMs and hybrid combinations of them on the same dataset to show the capabilities of the proposed ANFIS architecture.

The objectives of the present paper regarding NG demand forecasting, are briefly summarized in the following:

(a)

To develop a robust ANFIS model to provide accurate short-term forecasts for a number of cities in Greece, using a relatively large dataset. At the same time, the authors perform model fine-tuning that can lead to high accuracy in most distribution points. The proposed model is characterized by high flexibility, easiness of use and low execution time requirements.

(b)

To apply FCMs, ANNs and hybrid combinations of them to forecast NG demand in the same dataset, since these approaches have been proved as efficient techniques for NG demand forecasting according to the relevant literature.

(c)

To assess the performance of these soft computing methods in terms of prediction accuracy using well-known evaluation metrics.

(d)

To compare forecasting accuracy results of the proposed approach with those of the other soft computing and ANN methods that were examined, and finally decide on which model offers the best forecasting accuracy.

The outline of this paper is as follows. Section 2 describes the datasets of NG demand, collected for 10 Greek distribution points, as well as the proposed methodology of ANFIS for NG demand prediction using a well-defined set of evaluation performance metrics. Section 3 presents the results of the investigated ANFIS architectures. In the same section, a comparative analysis with other traditional neural networks and soft computing methods was performed for the same dataset. The discussion of results, which is also included in Section 3, presents the main outcomes of a meticulous ANFIS exploration analysis, along with ANFIS advantageous features. These are compared with ANNs and FCMs, and their overall contribution in NG forecasting is presented. Section 4 summarizes the paper, presenting future challenges in energy demand forecasting and highlighting further research directions.

2. Materials and Methods

This study aims to develop an ANFIS architecture capable of forecasting short-term NG consumption demand of the 10 main cities in Greece using the dataset that was provided by the Hellenic Gas Transmission System Operator S.A. (DESFA) [103]. The developed ANFIS approach deployed the aforementioned dataset along with other variables like the average daily temperature data for all the examined cities to accomplish forecasting. The results produced were further compared with those calculated by ANN and other soft computing techniques like FCM and hybrid-ANN to prove the prediction performance of the ANFIS prediction tool. Details on the dataset and its features, as well as the proposed methodology, are provided below. MATLAB M-file environment version 9.3.0.71 (R2017b) was used to program ANFIS networks and develop ANFIS models.

2.1. Dataset

The dataset covers ten different prediction datasets of historical data referring to ten cities all over Greece (Alexandroupoli, Athens, Drama, Karditsa, Larissa, Markopoulo, Serres, Thessaloniki, Trikala and Volos) and was linked to the values of gas demand for eight (8) previous years, in total. It should be mentioned that the time period for each dataset (city) was not the same in duration and did not correspond to the same years of data with all the other datasets collected.

Table 1. Time period referred to in each time-series dataset for all cities.

City	Time Period of the Examined Data	City	Time Period of the Examined Data
Alexandroupoli	2/2013–10/2018	Markopoulo	3/2010–10/2018
Athens	3/2010–10/2018	Serres	6/2013–10/2018
Drama	9/2011–10/2018	Thessaloniki	3/2012–10/2018
Karditsa	5/2014–10/2018	Trikala	9/2012–10/2018
Larissa	3/2010–10/2018	Volos	3/2010–10/2018

depicts the duration in years that is linked to each dataset collected and used in this case study. The historical datasets for 15 Greek cities were initially provided by the NG grid company of Greece, DESFA, which is responsible for the operation, management, exploitation and development of the Greek NG system and its interconnections. However, the authors, after thoroughly reviewing the available datasets, decided to include only 10 out of 15 cities in their case study, since these datasets contained less outliers and missing values than the other 5 datasets that were finally rejected, for data consistency purposes. For the datasets that were finally included in this work, a preliminary preprocessing phase was performed, where the insignificant outliers were removed, and any missing values were substituted with the average real value of the previous two days demand. The real data that were used for ANFIS modeling, performance evaluation and comparison with other popular forecasting methods were then split into training and testing samples. For all cities, the last year of each dataset (from November 2017 up to October 2018) was devoted to testing, whereas the rest of the years were used for training the developed ANFIS model.

In order to properly forecast day-ahead NG consumption demand of Greece, the proper number and type of input parameters should be selected. So, five factors were carefully considered as input parameters and the amount of one-day-ahead NG consumption demand of each distribution point was the output parameter. The prediction model was based on observations of past NG consumption, weather data, and calendar indicators, which are all among the most important input variables for prediction of NG consumption [15]. In particular, the dataset contains historical data of NG consumption of each city’s distribution point, the daily average temperature of the area in Celsius degrees, a month indicator and a day indicator. As regards the previous NG consumption data, these are linked to two different input variables: demand of a day before and current day demand. The temperature data are obtained by the nearest to the distribution gas point meteorological station. Concerning the calendar indicators (month and day), they need to undergo certain data form preprocessing before their use. Specifically, two different input indicators need to be considered for each one of the two variables. We define k = 1,2, …, 12 as the month index (1 January, 2 February, …, 12 December) and l = 1,2, …, 7 as the day index (1 Monday, 2 Tuesday, …7 Sunday). Following the coding procedure as presented in [104], the index for the month is scaled to the range [1/12, 1] in which the months of the year from January to December take successive values of the scaled index. That is, January has the value of 1/12 and December the value of 1. Similarly, the days of the week take successive values in the scaled range [1/7, 1], in which Monday and Sunday take the values of 1/7 and 1, respectively. All these parameters constituting the actual recorded data are briefly presented in

Table 2. Input and output parameters.

Type	Parameter	Unit
Input	Demand of a day before	MWh
Input	Current day demand	MWh
Input	Daily average temperature	Celsius degrees
Input	Month indicator	K = 1/12, 2/12, …, 1
Input	Day indicator	l = 1/7, 2/7, …, 1
Output	A day ahead NG demand	MWh

.

All data that compose the investigated dataset underwent a normalization process. This was necessary because all entries needed to have the same limited range of values so the model produces meaningful results [105].

The algorithm that was used for data normalization is the Min-Max, which scales the values of the dataset linearly over a specific range. As described in previous works [13,105], each variable was normalized to [0,1] before the forecasting model was applied. The normalized variable took its original value when the testing phase was implemented. Data normalization was carried out mathematically, as follows:

$$x_i^{\left(new\right)}=\frac{x_i-x^{\left(min\right)}}{x^{\left(max\right)}-x^{\left(min\right)}},\forall i=1,2,\dots ,N$$

(1)

where $x^{\left(new\right)}$ is the normalized value of the variable x, and $x^{\left(min\right)}$ and $x^{\left(max\right)}$ are, respectively, the minimum and maximum values of the concerned variable x.

2.2. Methods

2.2.1. Adaptive Neuro-Fuzzy Inference System (ANFIS)

The adaptive neuro-fuzzy inference system (ANFIS) uses an architecture that is based on both ANN and fuzzy logic principles and takes advantage of the benefits of both in a single framework. It can be described by the fuzzy “IF-THEN” rules from the Takagi and Sugeno (TS) type [106] as follows:

$$\begin{array}{c}{R}_i:if x_1=A_{i,1} and\dots and x_k=A_{i,k}\\ then y_i=b_{i,0}+b_{i,1}{x}_1+b_{i,2}{x}_2+\dots +b_{i,k}{x}_k\end{array}$$

(2)

where $A_{i,k}$ is the membership function associated with input variables $x_k$ and n is the number of inputs.

A typical ANFIS network is a five-layer structure consisting of the fuzzy layer, the product layer, the normalized layer, the de-fuzzy layer and the total output layer [3,107,108], as depicted in Figure 1.

In the first layer, every node i represents a linguistic label and is described by the following membership function, as given in Equation (3).

$$A_{i,k}\left(x_r\right)=e^{-{(\frac{\left(x_r-v_{i,k}\right)}{{\sigma }_{\iota ,\kappa }})}^2},\hspace{1em}for r=1, 2,\dots ,i$$

(3)

where $A_{i,k}$ is the membership function which is considered to be Gaussian and is described by the center ν and the spread σ.

In the second layer, the firing strength of the rule is computed using multiplicative operator, as presented in Equation (4). Firing strength is the weight degree of the IF-THEN rule and determines the shape of the output function for that rule.

$$w_i={\displaystyle \prod }_{k=1}^n{A}_{i,k}\left(x_k\right)$$

(4)

In the third layer, the i-th node calculates the ratio of the i-th rule’s firing strength to the sum of the firing strength of all rules. This is the normalization layer which normalizes the strength of all rules and the output of each node is given by Equation (5).

$${\overline{w}}_i=\frac{w_i}{{{\displaystyle \sum }}_{i=1}{w}_i}$$

(5)

In the fourth layer, each node is an adaptive node with a function given by Equation (6). In this layer, each node calculates a linear function where its coefficients are adapted by using the error function of the multilayer feed-forward neural network.

$${\overline{y}}_i={\overline{w}}_i{y}_i$$

(6)

In the fifth layer, there is only a fixed node indicated as the sum of the net outputs of the nodes in Layer 4. It computes the overall output as the sum of all incoming inputs and is expressed by Equation (7).

$$y={\displaystyle \sum }_i{\overline{y}}_i$$

(7)

ANFIS uses a hybrid learning algorithm to train the model. The back-propagation algorithm is used to train the parameters in Layer 1, whereas a variation of least-squares approximation or back-propagation algorithm is used for training the parameters of the fourth layer [108,109].

2.2.2. Proposed ANFIS Architecture Applied in Natural Gas Consumption Forecasting

In order to develop an efficient ANFIS model for NG demand forecasting, the authors needed to follow a certain process regarding the design of model’s architecture as well as an exploration process that will properly configure the input and training parameters of the examined model. Priority was given to the definition of the FIS architecture before the training of the network [110]. Among various fuzzy inference system (FIS) models, the Sugeno fuzzy model is the most widely used because of its higher interpretability and computational ability, that includes embedded optimal and adaptive techniques [111]. In order to create a fuzzy rule, the input space needs first to be divided. Two methods are used to divide space, comprised by input variables: the grid partitioning method and the subtractive clustering method. The main difference between these two functions refers to the way the partition of the input space is created.

In grid partitioning [109], the input space is divided into a grid-like structure without overlapping parts. Grid partitioning performs partitioning of the input space using all possible combinations of membership functions of each variable. This method is used when the number of input variables is small. For example, for 10 input variables and two membership functions for each input variable, then the input space is divided into 2¹⁰ = 1024 specific areas, representing one rule for each specific area, and the total number of rules is 1024, which is a very complicated structure. Therefore, the grid partitioning method is mainly used when the number of input variables is small.

On the other hand, the subtractive clustering method divides the input space into appropriate clusters, even if the user does not specify their number. If the size of the cluster becomes small, then the number of clusters increases, thus increasing the number of fuzzy rules. A rule is created for each cluster, whereas different values for parameters, like range of influence, squash factor, accept ratio and reject ratio, need to be explored for determining an efficient architecture, which will keep the balance between the total number of ANFIS parameters and the total number of rules.

Considering the above specifications, the authors used the Grid partition option to define the FIS architecture due to its simplicity, less time-consuming performance as well as it can easily explore the number and type of membership function (MF). In this stage, the number and type of membership functions of each input variable, along with the rules and values of parameters that belong to these functions, were determined using the option of Grid partition.

When implementing an ANFIS architecture, researchers should have in mind that there is one main restriction: the number of input variables. When these are more than five, then the number of the IF-THEN rules and the computational time also increase, hindering ANFIS to model output with respect to inputs [110]. Thus, in this study, five variables were chosen as input parameters, i.e., month, day, temperature, demand of a day before and demand of current day. As described above, a day-ahead consumption demand was selected as the output variable whose value can be produced by choosing between the option of linear or constant type of MF.

Finding the most efficient ANFIS architecture is a demanding task and entails a rigorous exploration process. Since our concern focuses on the increment of network’s accuracy and decrement of the errors, five necessary configurations should be considered in this direction: (i) the number of membership functions, (ii) types of MF (triangular, trapezoidal, bell-shaped, Gaussian and sigmoid), (iii) types of output MF (constant or linear), (iv) optimization methods (hybrid or back propagation) and (v) the number of epochs [112]. For the convenience of readers, these steps are visually represented in the flowchart in Figure 2.

The aforementioned set of configurations needs to be deployed in order to generate FIS and next to train the ANFIS model. Accordingly, the dataset that included the five input variables (i.e., month, day, temperature, demand of a day before, demand of current day) was selected to determine the only output (day-ahead demand). Initially, the training dataset was loaded in the ANFIS tool, as shown in Figure 3a. The next step was the design of the neuro-fuzzy model using the option “Generate FIS”. The Grid partition option was selected according to the description above (see Figure 3a). These two settings, concerning the fuzzy input variables along with their membership functions, are the most important parts to design the ANFIS. An example of selecting the number and type of MFs is illustrated in Figure 3b.

The number and type of membership function were assigned to the input parameters following the trial-and-error approach. The different types of MF that are offered by the MATLAB ANFIS editor include the triangular, trapezoidal, generalized bell (Gbell), Gaussian curve, Gaussian combination, difference between two sigmoid functions and product of two sigmoid functions (see Figure 3b). Regarding the type of output MFs, in the Sugeno-type fuzzy system, there are two options: a constant-type conclusion or a linear-type conclusion function. In the case of linear function, the output y is defined as:

$$y=k_0+k_1\ast x_1+k_2\ast x_2+\dots +k_n\ast x_n$$

(8)

where x₁, x₂, …, x_n are the n inputs. In this case, ANFIS needs to define k₀, k₁, k₂ up to k_n, and it is very time consuming to efficiently calculate the outputs when a large number of parameters are considered. On the other hand, when a constant MF is selected, the algorithm needs to define only one parameter to provide a reliable forecasted value. Thus, the computational time is really low.

The selected configuration also includes the hybrid optimization method, while the number of epochs selected to train the model was between 10 and 50. The hybrid optimization method uses the back propagation learning algorithm for parameters associated with input MF and the least-square estimation algorithm for parameters associated with output MF; thus, it was selected as the most proper one [113]. Various sets of ANFIS configurations are presented in

Table 3. Different configurations of the selected ANFIS architectures regarding constant output membership function (MF).

ANFIS Run	Type of Input MF	Number of MFs	Type of Output MF	Number of Epochs	Learning Method
1	trimf	2-2-2-2-2	Constant	10	Hybrid
2	trapmf	2-2-2-2-2	Constant	10	Hybrid
3	gbellmf	2-2-2-2-2	Constant	10	Hybrid
4	Gaussmf	2-2-2-2-2	Constant	10	Hybrid
5	Gauss2mf	2-2-2-2-2	Constant	10	Hybrid
6	pimf	2-2-2-2-2	Constant	10	Hybrid
7	dsigmf	2-2-2-2-2	Constant	10	Hybrid
8	psigmf	2-2-2-2-2	Constant	10	Hybrid
9	trimf	2-2-3-3-3	Constant	10	Hybrid
10	trapmf	2-2-3-3-3	Constant	10	Hybrid
11	gbellmf	2-2-3-3-3	Constant	10	Hybrid
12	Gaussmf	2-2-3-3-3	Constant	10	Hybrid
13	Gauss2mf	2-2-3-3-3	Constant	10	Hybrid
14	pimf	2-2-3-3-3	Constant	10	Hybrid
15	dsigmf	2-2-3-3-3	Constant	10	Hybrid
16	psigmf	2-2-3-3-3	Constant	10	Hybrid
17	trimf	3-3-3-2-2	Constant	10	Hybrid
18	trapmf	3-3-3-2-2	Constant	10	Hybrid
19	gbellmf	3-3-3-2-2	Constant	10	Hybrid
20	Gaussmf	3-3-3-2-2	Constant	10	Hybrid
21	trimf	3-3-3-3-3	Constant	10	hybrid
22	trimf	3-3-3-3-3	Constant	10	backpropa
23	trapmf	3-3-3-3-3	Constant	10	hybrid
24	trapmf	3-3-3-3-3	Constant	10	backpropa
25	gbellmf	3-3-3-3-3	Constant	10	hybrid
26	gbellmf	3-3-3-3-3	Constant	10	backpropa
27	trimf	3-3-3-3-3	Constant	30	hybrid
28	trimf	3-3-3-3-3	Constant	50	hybrid
29	trapmf	3-3-3-3-3	Constant	30	hybrid
30	trapmf	3-3-3-3-3	Constant	50	hybrid
31	gbellmf	3-3-3-3-3	Constant	30	hybrid
32	gbellmf	3-3-3-3-3	Constant	50	hybrid
33	trimf	3-3-4-4-4	Constant	10	hybrid
34	trimf	3-3-5-5-5	Constant	10	hybrid
35	trapmf	3-3-4-4-4	Constant	10	hybrid
36	trapmf	3-3-5-5-5	Constant	10	hybrid
37	gbellmf	3-3-4-4-4	Constant	10	hybrid
38	gbellmf	3-3-5-5-5	Constant	10	hybrid
39	gaussmf	3-3-3-3-3	Constant	10	hybrid
40	gaussmf	3-3-4-4-4	Constant	10	hybrid
41	gaussmf	3-3-5-5-5	Constant	10	hybrid
42	gauss2mf	3-3-3-3-3	Constant	10	hybrid
43	gauss2mf	3-3-4-4-4	Constant	10	hybrid
44	gauss2mf	3-3-5-5-5	Constant	10	hybrid
45	pimf	3-3-3-3-3	Constant	10	hybrid
46	pimf	3-3-4-4-4	Constant	10	hybrid
47	pimf	3-3-5-5-5	Constant	10	hybrid

, regarding different sets of number for MFs, as considered by the authors of this work.

Regarding the output MFs, constant and linear MF were accordingly investigated after certain numbers of experiments conducted. From these experiments and for the linear output, it was observed that the number of rules increases significantly, as well as the computational time, even in the case of problems with a small number of inputs (see

Table A1. Different configurations of the selected ANFIS architectures regarding linear output MF.

Type of Input MF	Number of MFs	Type of Output MF	Number of Rules	MSE	RMSE	MAE	MAPE	R2	Time (s)
trimf	2-2-2-2-2	Linear	32	0.001195	0.034572	0.019426	11.72180	0.982121	148
trapmf	2-2-2-2-2	Linear	32	0.001358	0.036859	0.020861	12.12378	0.979559	148
gbellmf	2-2-2-2-2	Linear	32	0.001267	0.035603	0.019921	11.46446	0.980963	148
Gaussmf	2-2-2-2-2	Linear	32	0.001298	0.036038	0.020259	11.97794	0.980468	148
Gauss2mf	2-2-2-2-2	Linear	32	0.001406	0.037496	0.020878	11.26382	0.978860	148
pimf	2-2-2-2-2	Linear	32	0.001635	0.040442	0.022176	12.08298	0.975405	148
dsigmf	2-2-2-2-2	Linear	32	0.001423	0.037733	0.021062	11.26721	0.978592	148
psigmf	2-2-2-2-2	Linear	32	0.001423	0.037733	0.021062	11.26722	0.978592	148
trimf	2-2-3-3-3	Linear	108	0.001476	0.038430	0.020941	11.17862	0.977773	328
Gaussmf	2-2-3-3-3	Linear	108	0.002038	0.045149	0.023286	12.71720	0.969241	328

in Appendix A). Thus, the linear type was not considered as an appropriate parameter of output MF since it is extremely time consuming. In this context, the trial-and-error approach was followed for the selection of the input-output type of MFs. Figure 4 illustrates an indicative ANFIS model, which was constructed with the following configuration set: 3-3-3-2-2, gbell MF, constant output MF, 10 epochs, hybrid.

2.2.3. Testing and Evaluation

The testing process for the ANFIS model was accomplished by using the testing data, which were completely unknown to the model. The predictor makes predictions on each day and finally compares the calculated predicted value with the real value. For example, considering the city of Volos, the predicted values that are illustrated in red in Figure 5 are compared with the real values (in blue color).

In order to evaluate the prediction of NG demand, five well known and commonly used statistical indicators were introduced, i.e., mean square error (MSE), root mean square error (RMSE), mean absolute error (MAE), mean absolute percentage error (MAPE) and coefficient of determination (R²). The mathematical equations of the statistical indicators are described below.

• Mean squared error:

  $$\mathrm{MSE}=\frac{1}{T}{\displaystyle \sum }_{t=1}^T{\left(Z\left(t\right)-X\left(t\right)\right)}^2$$

  (9)
• Root mean squared error:

  $$\mathrm{R}\mathrm{MSE}=\sqrt{\mathrm{MSE}}$$

  (10)
• Mean absolute error:

  $$\mathrm{MAE}=\frac{1}{T}{\displaystyle \sum }_{t=1}^T\left|Z\left(t\right)-X\left(t\right)\right|$$

  (11)
• Mean absolute percentage error:

  $$\mathrm{MAPE}=\frac{1}{T}{\displaystyle \sum }_{t=1}^T\left|\frac{Z\left(t\right)-X\left(t\right)}{Z\left(t\right)}\right|$$

  (12)
• Coefficient of determination:

  $$\mathrm{R}=\frac{T{{\displaystyle \sum }}_{t=1}^{T}Z\left(t\right)·X\left(t\right)-\left({{\displaystyle \sum }}_{t=1}^{T}Z\left(t\right)\right)\left({{\displaystyle \sum }}_{t=1}^{T}X\left(t\right)\right)}{\sqrt{T{{\displaystyle \sum }}_{t=1}^T{\left(Z\left(t\right)\right)}^2-{\left({{\displaystyle \sum }}_{t=1}^{T}Z\left(t\right)\right)}^2}·\sqrt{T{{\displaystyle \sum }}_{t=1}^T{\left(X\left(t\right)\right)}^2-{\left({{\displaystyle \sum }}_{t=1}^{T}X\left(t\right)\right)}^2}}$$

  (13)

  where X(t) is the forecasted value of the NG at the t-th iteration, and Z(t) is the actual value of the NG at the t-th iteration, t = 1, …, T, where T is the number of testing records.

Higher values of R², i.e., closer to 1, mean better model performance and the regression line fits the data well. A coefficient of determination value of 1.0 points out that the regression curve fits the data perfectly.

3. Results

This section presents the exploration analysis results for the various ANFIS architectures as proposed in Section 2.2.2. Considering the steps proposed in Section 2.2, the initial dataset is split into training and testing. During the training process, the ANFIS model is designed for each one of the suggested configurations. After the training process of the ANFIS finishes, NG consumption demands for the next day (one day ahead prediction) are calculated from the generated FIS. The NG consumption results only for the city of Athens (as an indicative example) regarding all configurations tested, are presented in

Table 4. Testing results for the examined ANFIS architectures concerning the city of Athens.

Anfis Run	Type of Input MF	Number of MFs	Type of Output MF	Number of Epochs	Optimization	MSE	RMSE	MAE	MAPE	R2
1	trimf	2-2-2-2-2	Constant	10	Hybrid	0.0010	0.0320	0.0192	12.6882	0.9849
2	trapmf	2-2-2-2-2	Constant	10	Hybrid	0.0013	0.0366	0.0245	19.8878	0.9806
3	gbellmf	2-2-2-2-2	Constant	10	Hybrid	0.0011	0.0335	0.0209	14.7498	0.9834
4	Gaussmf	2-2-2-2-2	Constant	10	Hybrid	0.0011	0.0326	0.0201	13.8422	0.9842
5	Gauss2mf	2-2-2-2-2	Constant	10	Hybrid	0.0011	0.0324	0.0197	13.5785	0.9845
6	pimf	2-2-2-2-2	Constant	10	Hybrid	0.0015	0.0389	0.0254	19.5486	0.9782
7	dsigmf	2-2-2-2-2	Constant	10	Hybrid	0.0014	0.0378	0.0244	18.7851	0.9794
8	psigmf	2-2-2-2-2	Constant	10	Hybrid	0.0014	0.0378	0.0244	18.7851	0.9794
9	trimf	2-2-3-3-3	Constant	10	Hybrid	0.0015	0.0388	0.0232	15.8840	0.9774
10	trapmf	2-2-3-3-3	Constant	10	Hybrid	0.0020	0.0448	0.0269	19.2727	0.9698
11	gbellmf	2-2-3-3-3	Constant	10	Hybrid	0.0014	0.0379	0.0227	15.5056	0.9785
12	Gaussmf	2-2-3-3-3	Constant	10	Hybrid	0.0014	0.0379	0.0226	15.7640	0.9784
13	Gauss2mf	2-2-3-3-3	Constant	10	Hybrid	0.0017	0.0410	0.0241	15.8227	0.9747
14	pimf	2-2-3-3-3	Constant	10	Hybrid	0.0130	0.1141	0.0347	21.6717	0.8552
15	dsigmf	2-2-3-3-3	Constant	10	Hybrid	0.0020	0.0448	0.0254	16.7809	0.9698
16	psigmf	2-2-3-3-3	Constant	10	Hybrid	0.0020	0.0448	0.0254	16.7809	0.9698
17	trimf	3-3-3-2-2	Constant	10	Hybrid	0.0012	0.0348	0.0210	14.6116	0.9819
18	trapmf	3-3-3-2-2	Constant	10	Hybrid	0.0018	0.0430	0.0297	27.0255	0.9723
19	gbellmf	3-3-3-2-2	Constant	10	Hybrid	0.0013	0.0355	0.0212	14.4247	0.9810
20	Gaussmf	3-3-3-2-2	Constant	10	Hybrid	0.0011	0.0337	0.0198	12.8988	0.9829
21	trimf	3-3-3-3-3	Constant	10	hybrid	0.0021	0.0455	0.0242	15.4964	0.9698
22	trimf	3-3-3-3-3	Constant	10	backpropa	0.0559	0.2365	0.1610	74.9654	0.7447
23	trapmf	3-3-3-3-3	Constant	10	hybrid	0.0031	0.0556	0.0281	21.5814	0.9562
24	trapmf	3-3-3-3-3	Constant	10	backpropa	0.0501	0.2238	0.1527	72.2629	0.7404
25	gbellmf	3-3-3-3-3	Constant	10	hybrid	0.0014	0.0374	0.0217	14.6538	0.9791
26	gbellmf	3-3-3-3-3	Constant	10	backpropa	0.0015	0.0392	0.0265	25.1796	0.9793
27	trimf	3-3-3-3-3	Constant	30	hybrid	0.0016	0.0403	0.0224	13.3194	0.9759
28	trimf	3-3-3-3-3	Constant	50	hybrid	0.0017	0.0417	0.0224	13.2276	0.9745
29	trapmf	3-3-3-3-3	Constant	30	hybrid	0.0029	0.0539	0.0245	17.2238	0.9590
30	trapmf	3-3-3-3-3	Constant	50	hybrid	0.0017	0.0416	0.0233	16.4612	0.9745
31	gbellmf	3-3-3-3-3	Constant	30	hybrid	0.0013	0.0366	0.0213	13.2077	0.9799
32	gbellmf	3-3-3-3-3	Constant	50	hybrid	0.0019	0.0432	0.0236	13.4445	0.9724
33	trimf	3-3-4-4-4	Constant	10	hybrid	0.0023	0.0479	0.0251	15.4225	0.9662
34	trimf	3-3-5-5-5	Constant	10	hybrid	0.0078	0.0884	0.0320	17.3158	0.9006
35	trapmf	3-3-4-4-4	Constant	10	hybrid	0.0021	0.0454	0.0275	23.1769	0.9695
36	trapmf	3-3-5-5-5	Constant	10	hybrid	0.0098	0.1084	0.0450	19.3158	0.8806
37	gbellmf	3-3-4-4-4	Constant	10	hybrid	0.0022	0.0472	0.0256	16.1637	0.9669
38	gbellmf	3-3-5-5-5	Constant	10	hybrid	0.0044	0.0660	0.0307	18.0977	0.9376
39	gaussmf	3-3-3-3-3	Constant	10	hybrid	0.0013	0.0365	0.0212	13.8235	0.9800
40	gaussmf	3-3-4-4-4	Constant	10	hybrid	0.0019	0.0431	0.0241	14.7715	0.9720
41	gaussmf	3-3-5-5-5	Constant	10	hybrid	0.0056	0.0746	0.0314	17.5307	0.9185
42	gauss2mf	3-3-3-3-3	Constant	10	hybrid	0.0017	0.0409	0.0235	16.2626	0.9755
43	gauss2mf	3-3-4-4-4	Constant	10	hybrid	0.0040	0.0632	0.0260	17.3863	0.9407
44	gauss2mf	3-3-5-5-5	Constant	10	hybrid	0.0072	0.0847	0.0290	17.7331	0.9048
45	pimf	3-3-3-3-3	Constant	10	hybrid	0.1224	0.3499	0.0482	26.0901	0.3608
46	pimf	3-3-4-4-4	Constant	10	hybrid	0.0026	0.0510	0.0307	24.7626	0.9615
47	pimf	3-3-5-5-5	Constant	10	hybrid	0.0022	0.0466	0.0285	22.3553	0.9678

, whereas

Table 5. Testing results for the best three ANFIS architectures of each city based on MAPE value.

City	Anfis Run	Type of Input MF	Number of MFs	Number of Rules	Time (s)	MSE	RMSE	MAE	MAPE	R2
Alexandroupoli	17	trimf	3-3-3-2-2	72	5	0.0024	0.0494	0.0351	10.5278	0.9638
	39	gaussmf	3-3-3-3-3	243	47	0.0031	0.0557	0.0355	10.1556	0.9538
	20	gaussmf	3-3-3-2-2	72	5	0.0023	0.0480	0.0341	10.1123	0.9659
Athens	1	trimf	2-2-2-2-2	32	7	0.0021	0.0457	0.0295	20.1799	0.9825
	17	trimf	3-3-3-2-2	108	19	0.0026	0.0511	0.0315	19.7972	0.9786
	20	gaussmf	3-3-3-2-2	108	19	0.0022	0.0467	0.0306	21.2929	0.9818
Drama	17	trimf	3-3-3-2-2	108	19	0.0026	0.0511	0.0363	6.2547	0.8997
	1	trimf	2-2-2-2-2	32	5	0.0026	0.0513	0.0361	6.2235	0.8975
	20	gaussmf	3-3-3-2-2	108	13	0.0026	0.0508	0.0371	6.4071	0.8995
Karditsa	17	trimf	3-3-3-2-2	108	12	0.0019	0.0434	0.0242	13.8394	0.9789
	1	trimf	2-2-2-2-2	32	4	0.0018	0.0421	0.0236	11.6196	0.9801
	4	gaussmf	2-2-2-2-2	32	4	0.0019	0.0431	0.0248	13.3841	0.9792
Larissa	1	trimf	2-2-2-2-2	32	4	0.0012	0.0352	0.0203	10.9568	0.9817
	4	gaussmf	2-2-2-2-2	32	4	0.0012	0.0352	0.0204	10.9833	0.9817
	20	gaussmf	3-3-3-2-2	108	19	0.0010	0.0314	0.0184	10.5236	0.9858
Markopoulo	1	trimf	2-2-2-2-2	32	5	0.0091	0.0956	0.0728	25.0887	0.6593
	4	gaussmf	2-2-2-2-2	32	5	0.0096	0.0980	0.0755	26.7510	0.6364
	17	trimf	3-3-3-2-2	108	19	0.0259	0.1609	0.1087	36.7174	0.5126
Serres	1	trimf	2-2-2-2-2	32	5	0.0007	0.0271	0.0176	10.4721	0.9839
	4	gaussmf	2-2-2-2-2	32	5	0.0008	0.0279	0.0185	11.2421	0.9831
	39	gaussmf	3-3-3-3-3	243	45	0.0008	0.0285	0.0194	12.1163	0.9824
Thessaloniki	17	trimf	3-3-3-2-2	108	13	0.0015	0.0382	0.0229	16.1046	0.9773
	20	gaussmf	3-3-3-2-2	108	13	0.0013	0.0363	0.0219	14.1944	0.9795
	39	gaussmf	3-3-3-3-3	243	45	0.0021	0.0459	0.0256	15.2032	0.9672
Trikala	1	trimf	2-2-2-2-2	32	4	0.0019	0.0433	0.0232	10.5817	0.9815
	4	gaussmf	2-2-2-2-2	32	4	0.0020	0.0450	0.0245	11.1412	0.9800
	20	gaussmf	3-3-3-2-2	108	13	0.0028	0.0530	0.0271	11.7631	0.9708
Volos	1	trimf	2-2-2-2-2	32	4	0.0021	0.0459	0.0317	13.2520	0.9564
	4	gaussmf	2-2-2-2-2	32	4	0.0021	0.0460	0.0314	13.1629	0.9563
	20	gaussmf	3-3-3-2-2	108	12	0.0020	0.0445	0.0323	13.9710	0.9588

gathers the best three results of NG consumption demand obtained from ANFIS for each of the 10 cities. To evaluate the performance of the models, the statistical indicator mean absolute percent error (MAPE) was used [91]. MAPE is a relative measurement, independent of scale, and it is the most common performance metric in time series forecasting, due to being reliable and valid [21].

The corresponding graphical representation of the results regarding the best three out of 47 total ANFIS architectures for the city of Athens is illustrated in Figure 6. Also, the best ANFIS model for each city can be found in

Table 6. The best ANFIS models for all the cities under investigation (10 epochs and hybrid optimization).

Title 1	Anfis Run	Type of Input MF	Number of MFs	Type of Output MF	Optimization	MSE	RMSE	MAE	MAPE	R2
Alexandroupoli	20	gaussmf	3-3-3-2-2	Constant	Hybrid	0.0023	0.0480	0.0341	10.1123	0.9659
Athens	17	trimf	3-3-3-2-2	Constant	Hybrid	0.0026	0.0511	0.0315	19.7972	0.9786
Drama	1	trimf	2-2-2-2-2	Constant	Hybrid	0.0026	0.0513	0.0361	6.2235	0.8975
Karditsa	1	trimf	2-2-2-2-2	Constant	Hybrid	0.0018	0.0421	0.0236	11.6196	0.9801
Larissa	20	gaussmf	3-3-3-2-2	Constant	Hybrid	0.0010	0.0314	0.0184	10.5236	0.9858
Markopoulo	1	trimf	2-2-2-2-2	Constant	Hybrid	0.0091	0.0956	0.0728	25.0887	0.6593
Serres	4	gaussmf	2-2-2-2-2	Constant	Hybrid	0.0008	0.0279	0.0185	11.2421	0.9831
Thessaloniki	20	gaussmf	3-3-3-2-2	Constant	Hybrid	0.0013	0.0363	0.0219	14.1944	0.9795
Trikala	4	gaussmf	2-2-2-2-2	Constant	Hybrid	0.0020	0.0450	0.0245	11.1412	0.9800
Volos	4	gaussmf	2-2-2-2-2	Constant	Hybrid	0.0021	0.0460	0.0314	13.1629	0.9563

, which provides the most reliable ANFIS architecture results for each city. All the results have been previously ranked, based on the minimum value of MAPE and subsequently the minimum values of MSE, RMSE and MAE. The priority was given to MAPE as one of the most crucial evaluation metrics, according to the literature [91,96], which was used in this study to compare various models obtained from ANFIS and other soft computing and neural networks methods. As a relative and easy to interpret measurement, MAPE is reliable, valid and independent of scale. The smaller the values of MAPE are, the closer the forecasted values are to the actual values.

As illustrated in Table 6, ANFIS models appear to perform best mostly when triangular MFs are used for the input variables: three MFs for the first three input variables (month, day of week and mean temperature) and two or three MFs for the other two input variables (daily demand for current day and one day before). Also, constant MFs are selected for the output variable and hybrid optimization method. The graphical representation of the best ANFIS models for each city is illustrated in Figure 7.

It is worth mentioning that all three most efficient ANFIS architectures with respect to MAPE values have triangular or gaussian MFs and 2-2-2-2-2 or 3-3-3-2-2 number of input MFs, whereas the output MF is constant and the learning algorithm is hybrid. In addition, the application of other MFs combinations does not seem to give results that could be on top of the list. Due to the limitation of the total number of parameters that should not exceed the number of training data pairs, the number of MFs was chosen based on the number of input parameters. Figure 8 shows the exponential increase in the number of rules when the number of MFs increases, whereas

Table A2. Running time and number of rules for all the proposed ANFIS configurations.

Type of Input MF	Number of MFs	Type of Output MF	Number of Epochs	Optimization	Number of Rules	Time Run
trimf, trapmf, gbell, gauss, pim, sigm	2-2-2-2-2	Constant	10	Hybrid	32	7 s
trimf, trapmf, gbell, gauss, pim, sigm	2-2-3-3-3	Constant	10	Hybrid	108	11 s
trimf, trapmf, gbell, gauss, pim, sigm	3-3-3-2-2	Constant	10	Hybrid	108	19 s
trimf, trapmf, gbell	3-3-3-3-3	Constant	10	Hybrid	243	68 s
trimf	3-3-4-4-4	Constant	10	Hybrid	576	10 min 10 s
trimf	3-3-5-5-5	Constant	10	Hybrid	1125	40 min
trapmf	3-3-4-4-4	Constant	10	Hybrid	576	12min
trapmf	3-3-5-5-5	Constant	10	Hybrid	1125	70 min
gbellmf	3-3-4-4-4	Constant	10	Hybrid	576	12 min 35 s
gbellmf	3-3-5-5-5	Constant	10	Hybrid	1125	50 min
gaussmf	3-3-3-3-3	Constant	10	Hybrid	243	4 min
gaussmf	3-3-4-4-4	Constant	10	Hybrid	576	25 min
gaussmf	3-3-5-5-5	Constant	10	Hybrid	1125	47 min
gauss2mf	3-3-3-3-3	Constant	10	Hybrid	243	4 min
gauss2mf	3-3-4-4-4	Constant	10	Hybrid	576	25 min
gauss2mf	3-3-5-5-5	Constant	10	Hybrid	1125	47 min
pimf	3-3-3-3-3	Constant	10	Hybrid	243	3.5 min
pimf	3-3-4-4-4	Constant	10	hybrid	576	20 min
pimf	3-3-5-5-5	Constant	10	hybrid	1125	42 min

in Appendix A gathers the time and number of rules for all the proposed ANFIS configurations.

3.1. Comparison with ANNs, FCMs and Hybrid FCM-ANN

To further investigate the performance of the proposed ANFIS architectures, an extensive comparative analysis between the state of the art ANNs, soft computing methods of FCMs and their hybrid combination of FCMs with ANNs was performed.

The architecture of the analyzed ANN was a multilayer feed forward network with an input layer containing five inputs (month, day, temperature, demand of a day before, current demand), a hidden layer with 10 neurons, and an output layer with one output (a day-ahead demand prediction). The authors used the sigmoidal activation function in all layers and implemented Levenberg–Marquardt algorithm to train the network.

The soft computing method of fuzzy cognitive map with evolutionary learning capabilities, such as the real-coded genetic algorithm (RCGA-FCM) and structure optimization genetic algorithm (SOGA-FCM) [100], were used for time series modeling and prediction of day-ahead NG energy demand. For FCM learning, we implemented RCGA-FCM and SOGA-FCM. A short description on the applied evolutionary-based FCM approaches is given in Appendix B. The implementations of FCMs differ from ANFIS, even though they both belong to the soft computing family.

In this research study, we used the dynamic model type (Equation (B1)) which is found in Appendix B, with sigmoidal transformation function. FCMs learned with the use of RCGA and SOGA algorithm contain five concepts (month, day, temperature, demand of a day before, current demand) [114].

The applied hybrid approach for time series prediction is based on FCMs and ANNs and was previously proposed in [18,19]. It allows us to select the most significant concepts for FCM using SOGA. These concepts are used as the inputs for ANN. In the hybrid approach, we used artificial neural networks with an input layer with five inputs selected by the SOGA-FCM approach, a hidden layer with 10 neurons and an output layer with one output (one day-ahead demand prediction). Sigmoidal activation function and Levenberg–Marquardt learning algorithm were used. All the simulations for FCMs and hybrid FCM-ANN configurations were performed with the software tool ISEMK [115] which has been developed for time series forecasting purposes. An analytical description of FCM-based models and hybrid FCM-ANN can be found in [13,101,102,116], whereas they are used in this work only for comparison purposes.

In what follows,

Table 7. Comparison results among the artificial neural network (ANN), fuzzy cognitive map (FCM), hybrid FCM-ANN and best ANFIS architectures of each city.

City	Method	MSE	RMSE	MAE	MAPE	R2
Alexandroupoli	RCGA-FCM	0.0047	0.0684	0.0538	17.6233	0.9450
	SOGA-FCM	0.0045	0.0672	0.0526	17.1707	0.9484
	ANN	0.0042	0.0645	0.0505	16.1131	0.9439
	Hybrid FCM-ANN	0.0034	0.0579	0.0427	14.3034	0.9498
	Best ANFIS	0.0023	0.0480	0.0341	10.1123	0.9659
Athens	RCGA-FCM	0.0022	0.0473	0.0303	23.5985	0.9676
	SOGA-FCM	0.0029	0.0539	0.0337	22.7453	0.9646
	ANN	0.0010	0.0323	0.0198	14.2464	0.9844
	Hybrid FCM-ANN	0.0014	0.0374	0.0230	17.5418	0.9790
	Best ANFIS	0.0026	0.0511	0.0315	19.7972	0.9786
Drama	RCGA-FCM	0.0080	0.0894	0.0749	12.9942	0.8691
	SOGA-FCM	0.0056	0.0748	0.0600	10.1766	0.8796
	ANN	0.0025	0.0501	0.0357	6.1657	0.9025
	Hybrid FCM-ANN	0.0028	0.0526	0.0363	6.2502	0.8941
	Best ANFIS	0.0026	0.0513	0.0361	6.2235	0.8975
Karditsa	RCGA-FCM	0.0039	0.0624	0.0379	27.5914	0.9591
	SOGA-FCM	0.0488	0.2210	0.1397	50.2112	0.9711
	ANN	0.0016	0.0405	0.0245	17.4579	0.9819
	Hybrid FCM-ANN	0.0017	0.0407	0.0245	18.4095	0.9817
	Best ANFIS	0.0018	0.0421	0.0236	11.6196	0.9801
Larissa	RCGA-FCM	0.0027	0.0515	0.0331	22.2481	0.9638
	SOGA-FCM	0.0025	0.0505	0.0328	22.9579	0.9649
	ANN	0.0013	0.0355	0.0209	13.2479	0.9812
	Hybrid FCM-ANN	0.0013	0.0356	0.0215	13.1974	0.9811
	Best ANFIS	0.0010	0.0314	0.0184	10.5236	0.9858
Markopoulo	RCGA-FCM	0.0075	0.0868	0.0726	26.0003	0.6975
	SOGA-FCM	0.0078	0.0883	0.0739	26.3345	0.6955
	ANN	0.0172	0.1310	0.1048	34.8594	0.4765
	Hybrid FCM-ANN	0.0070	0.0836	0.0667	23.7166	0.7094
	Best ANFIS	0.0091	0.0956	0.0728	25.0887	0.6593
Serres	RCGA-FCM	0.0017	0.0409	0.0274	16.5199	0.9648
	SOGA-FCM	0.0495	0.2225	0.1632	72.9785	0.9772
	ANN	0.0008	0.0275	0.0179	10.9948	0.9842
	Hybrid FCM-ANN	0.0008	0.0289	0.0190	11.5000	0.9821
	Best ANFIS	0.0008	0.0279	0.0185	11.2421	0.9831
Thessaloniki	RCGA-FCM	0.0029	0.0541	0.0339	29.9713	0.9565
	SOGA-FCM	0.0029	0.0539	0.0340	30.1471	0.9568
	ANN	0.0017	0.0412	0.0262	23.8748	0.9735
	Hybrid FCM-ANN	0.0019	0.0441	0.0266	23.8835	0.9696
	Best ANFIS	0.0013	0.0363	0.0219	14.1944	0.9795
Trikala	RCGA-FCM	0.0059	0.0770	0.0453	21.9722	0.9528
	SOGA-FCM	0.0433	0.2082	0.1287	42.7427	0.9715
	ANN	0.0020	0.0443	0.0258	14.1183	0.9804
	Hybrid FCM-ANN	0.0019	0.0432	0.0251	13.9034	0.9815
	Best ANFIS	0.0020	0.0450	0.0245	11.1412	0.9800
Volos	RCGA-FCM	0.0028	0.0526	0.0397	17.8195	0.9436
	SOGA-FCM	0.0027	0.0520	0.0395	17.8988	0.9445
	ANN	0.0020	0.0444	0.0319	13.2504	0.9588
	Hybrid FCM-ANN	0.0020	0.0446	0.0307	12.7881	0.9587
	Best ANFIS	0.0021	0.0460	0.0314	13.1629	0.9563

gathers the results of the explored ANN and soft computing models, which are straightforward compared with our best performed ANFIS configuration, for each one out of the 10 cities, suggested in this research work. In Figure 9, three indicative graphs of the cities Alexandroupoli, Athens, and Drama are illustrated regarding the predicted values of NG demand for all the best proposed architectures. In Appendix A, the corresponding graphs of the rest of the cities are presented (see Figure A1).

For a deeper analysis of the examined architectures (ANFIS, ANN, FCM, hybrid FCM-ANN), the authors report on further details regarding the parameters of each model used in this study. The ANN and FCM models were previously applied for NG demand prediction in several research works, such as those in [13,101,102,116]. The models were sufficiently described and the hyperparameters were properly configured to offer optimum performance of the investigated FCM models.

Table 8. Parameters and average running time for each architecture.

Architectures	Parameters for Athens City	Average Running Time
ANN	Multilayer feed forward network, six inputs, 10 neurons, one output, sigmoidal activation function, Levenberg-Marquardt learning, epochs = 20	16–20 s
RCGA-FCM	Uniform crossover with probability 0.4, Mühlenbein’s mutation with probability 0.4, ranking selection, elite strategy, population size 200, maximum number of generations 200	808 s
SOGA-FCM	Uniform crossover with probability 0.4, Mühlenbein’s mutation with probability 0.4, ranking selection, elite strategy, population size 200, maximum number of generations 200, learning parameters b1 = b2 = 0.01	799 s
Hybrid FCM-ANN	Multilayer feed forward network, four inputs selected by SOGA-FCM (month, temperature, demand of a day before, current demand), one hidden layer with 10 neurons, one output, sigmoidal activation function, Levenberg-Marquardt learning, epochs = 20	811 s
Best ANFIS	Triangular mf, 2-2-2-2-2 or 3-3-3-2-2, Constant output, epochs = 10, Hybrid optimization	4–19 s

depicts the optimum parameters for all cities considering the neural and FCM evolutionary methods (ANN, RCGA-FCM, SOGA-FCM, Hybrid), compared with the proposed best performed ANFIS. The average running time is also presented in Table 8, which was calculated for each soft computing architecture for all models. It is worth mentioning that we have conducted a rigorous exploratory analysis for all the investigated neuro-fuzzy, soft computing techniques and ANNs, with different parameters, for training and model optimization, to reach the highest prediction accuracy with respect to the evaluation metrics.

3.2. Discussion of Results

In this work, several ANFIS architectures were investigated, with respect to all the variables that were carefully determined in the developed model as reported in Section 2.1 and after different sets of model configurations were tested. However, only one ANFIS architecture reached the optimum performance, in terms of forecasting accuracy, considering the minimum value of MAPE and subsequently the minimum values of MSE, RMSE and MAE values produced. In particular, it emerged that the optimum ANFIS configuration is 2-2-2-2-2 with triangular MFs for input variables, which produces the most simple (concerning the number of rules), fast (see Table A2 in Appendix A) and accurate model for this energy forecasting problem. In general, it is observed that the best results are produced from the combination of triangular or gaussian MFs regarding the input variables, and the constant MFs regarding the output layers.

To further discuss the results produced and to show the effectiveness of the proposed forecasting methodology of ANFIS, the authors conducted a comparative analysis regarding the forecasting performance between the proposed technique, and other ANN and soft computing methods too, such as FCM, which were reported in the literature and have already been applied in the specific domain. The MAPE criterion was used to compare various models from ANN, evolutionary FCM, Hybrid FCM-ANN and ANFIS. The results are given in Table 7. For example, the MAPE values of RCGA-FCM, SOGA-FCM, ANN, hybrid and ANFIS models for the city of Alexandroupoli were calculated as 17.62%, 17.17%, 16.11%, 14.30% and 10.11%, respectively. The smaller the values of MAPE are, the closer to the actual values the forecasted values are. The best result was obtained from the ANFIS model. The respective figures in the text and in Appendix A have been updated with the new prediction values.

Considering the same dataset linked to only three cities (Athens, Thessaloniki, and Larissa) out of the ten that participated in our study, a day-ahead NG consumption prediction was investigated in [117], applying ANN and LSTM approaches and in [102], implementing the SOGA-FCM method and a hybrid combination of it. Furthermore, an ensemble FCM prediction methodology concerning the same dataset was presented in [13], in which a recent soft computing technique for time series forecasting, using evolutionary fuzzy cognitive maps and their ensemble combination was compared to ANNs, as benchmark forecasting methods. These methods and their results in terms of MSE and MAE values for three benchmark cities are all gathered in the following table and certain results can be concluded. The main reason for selecting the statistical indicators MSE and MAE in the following figure is to accomplish a straightforward comparison with the results published in previous works.

In Figure 10, it can be noted that all methods achieve high accuracy in NG consumption predictions, using the same dataset. The best ANFIS approach seems to excel over the ensemble and hybrid methods. Consequently, the proposed ANFIS architecture, which handles the fuzziness of data more efficiently, outperforms all the other examined methods in most cases, with a rather remarkable difference. ANFIS is less time consuming and more flexible than ANN, and as it employs fuzzy rules and membership functions incorporating with real-world systems, it can be used as alternate method to ANN forecasting.

It is observed that the proposed method exhibits better or similar performance to other well-known ANN, FCM or hybrid FCM-ANN architectures for the ten cities under investigation. The produced results highlight the significance and superiority of neuro-fuzzy methods over the other examined methods in terms of prediction accuracy, when they deal with time series forecasting problems in energy. This is in accordance with the main advantageous features of ANFIS models, which are their ability to capture the nonlinear structure of a process, their adaptation capability, and fast training characteristics. As reported in the literature, ANFIS models are able to cope with the uncertainty and fuzziness that characterize the energy domain [118,119] when other intelligent methods cannot tackle them.

The main outcomes of this study can be summarized as follows:

• The proposed ANFIS method exhibits the best performance when certain configuration settings are selected for the examined datasets which are linked to ten cities of Greece. The authors concluded that a certain configuration is best for the examined ANFIS model, after having conducted a number of experiments and following a trial-and error approach. The best ANFIS model is based on a distinct architecture that features a 2-2-2-2-2 triangular or gaussian MF.
• The proposed ANFIS architecture is superior to the four benchmark and well-known ANN and FCM methods (ANN, SOGA-FCM, RCGA-FCM, Hybrid FCM-ANN), which have been efficiently used in NG consumption forecasting. The results presented in Table 7, which gathers various error indicators and the R², as prediction accuracy indices for all five architectures, show that the best ANFIS model holds the best prediction accuracy among all the methods that were included in this comparative analysis.
• The proposed ANFIS model shows significant capacity when applied to forecasting NG demand, since it exhibits better performance (see Table 7) with less running time (see Table 8) and more flexibility to handle fuzziness than other well-known ANN and FCM architectures.

4. Conclusions

This study proposes the ANFIS method to predict short-term demand of NG consumption. This approach is applied on the Greek territory and uses 10 different datasets provided by DESFA, that regard previous energy consumption historical data, for ten main cities. To decide the model’s proper architecture, the authors follow an exploration process regarding the best configuration of input and training parameters. The best ANFIS model is then compared to other well-known ANN and soft computing models that are commonly used for energy demand prediction purposes. The ANFIS method demonstrates significant performance in the field of energy demand prediction, outweighing the traditional ANN and FCM architectures. In addition, the running time of the proposed architecture is much less than those of other examined models, making it the right decision for day-ahead demand forecasting of NG. The findings of this study reveal that the highest forecasting accuracy emerged when the same model configuration was used for most of the cities, highlighting the generalization capabilities of the proposed architecture.

This work can be widely used in short-term demand forecasting for other countries too, with the same or similar input parameters, and can be also useful especially for distribution operators, providing them with the ability to make long-term planning decisions and apply the correct strategic policies in this direction. Following the literature, the ANFIS approach can be applied in various other domains such as medicine, environmental modelling, various energy systems, like solar and wind, as well as other engineering applications. As can be seen, the ANFIS application area is wide, and as regards the energy sector, this method finds great applicability due to its high prediction accuracy, robustness, and easiness to use.

The results show that the proposed algorithm, which was proven to be efficient, fast and robust, can be adopted by regulatory authorities and decision makers to perform rigorous forecasting of natural gas demand for the respective case cities and other cities in Greece too. The investigated approach is an accurate estimation method as it makes efficient short-term predictions in natural gas demand, showing minor deviations between the real and the predicting values. Since short-term natural gas forecasting is mostly used for the timely reservation of transport, storage capacity optimization, timely purchase of natural gas deliveries and capacity allocation, this method becomes critical to determine the energy policy for Greece and the wider area too, having overall a positive impact in natural gas consumption.

Future work is oriented in developing more advanced neuro-fuzzy models providing explainability and transparency in prediction tasks in diverse research domains, in order to evaluate the generalization capabilities of this approach. Furthermore, new forecast combination architectures of efficient deep learning and regularized recurrent neural networks for time series modelling and prediction in the energy sector will be investigated.