Multi-objective regression modeling for natural gas prediction with ridge regression and CMARS

Article history: Received: 17 February 2021 Accepted: 23 November 2021 Available Online: 1 January 2022 Residential customers are the main users generally need a great quantity of natural gas in distribution systems, especially, in the wintry weather season since it is particularly consumed for cooking and space heating. Hence, it ought to be noninterruptible. Since distribution systems have a restricted ability for supply, reasonable planning and prediction through the whole year, especially in winter seasons, have emerged as vital. The Ridge Regression (RR) is formulated mainly to decrease collinearity results through shrinking the regression coefficients and reducing the impact in the model of variables. Conic multivariate adaptive regression splines ((C)MARS) model is constructed as an effective choice for MARS by using inverse problems, statistical learning, and multi-objective optimization theories. In this approach, the model complexity is penalized in the structure of RR and it is constructed a relaxation by utilizing continuous optimization, called Conic Quadratic Programming (CQP). In this study, CMARS and RR are applied to obtain forecasts of residential natural gas demand for local distribution companies (LDCs) that require short-term forecasts, and the model performances are compared by using some criteria. Here, our analysis shows that CMARS models outperform RR models. For one-day-ahead forecasts, CMARS yields a MAPE of about 4.8%, while the same value under RR reaches 8.5%. As the forecast horizon increases, it can be seen that the performance of the methods becomes worse, and for a forecast one week ahead, the MAPE values for CMARS and RR are 9.9% and 18.3%, respectively.


Introduction
Local Distribution Companies (LDCs) are the service carriers of private natural gas purchasers within a particular transportation structure. When LDCs target unbalanced gas consumption from Transmission System Operator (TSO) pipelines, regulations, and related policies impose fairly high costs on them in the form of penalties. So correct forecasting is crucial here, as LCD demand is met through spot markets and much of the total demand is organized with global supply contracts through pipelines or liquefied natural gas. The operation must be flexible to compensate for fluctuations in demand. The fluctuation adjustment of energy demand under the constraints of system operation should be realized by predictive models organized for individual types of customers. For LCDs, a unique one-day and one-week forecast gives a discount to cost operations and elimination of penalties that occur due to unbalanced demand-supply quantities [1,2]. A problem is described as an ill-posed problem if a solution is not unique, present, or stable under perturbation on data. This means that if a small perturbation of the data may bring about a large perturbation of the solution. Ridge Regression is one of the most well-known structures to make these problems regular and stable [3,4]. It is also known as Tikhonov Regularization. In the statistical literature, there exist some approaches such as principal components regression, partial least squares, least absolute shrinkage and selection operator (LASSO), and ridge regression (RR) to prevent collinearity in traditional linear regression models. Here, RR and LASSO are regularization/penalization strategies that impose a constraint on the regression coefficients while principal components and partial least squares regression are variable subset decision strategies that employ linear combinations of the independent variables in the regression model. RR is formulated mainly to decrease collinearity results by shrinking the regression coefficients and reducing the impact in the model of variables. Here, it seems that shrinking the coefficient estimates may extensively decrease their variance [5]. MARS creates flexible models by employing piecewise linear functions. This continuous model gives a highquality way to model nonlinearities [6]. Nowadays, as a nonparametric model, MARS is efficiently applied to many areas of technology and science. Here, Conic Multivariate Adaptive Regression Splines (CMARS) model [7][8][9][10] is evolved for the backward stage of MARS. It is obtained as a model-based alternative and an effective choice to MARS. In this algorithm, a Penalized Residual Sum of Squares (PRSS) is applied for MARS as RR problem and it is worked out by Conic Quadratic Programming (CQP). So, Interior Point Methods [11,12] and their codes, e.g., MOSEK [13] can be applied by the technique of CMARS. In this paper, we represent natural gas forecasts for one day and one week in advance for a transmission system operator by using RR and CMARS. In this study, two multi-objective regression models are developed for short-term natural gas demand prediction using RR and CMARS for one day, and one week ahead with daily forecasting intervals. Here, the minimum temperature, the maximum temperature, and heating degree days of the daily average temperature are taken into consideration as different input variables. The MAPE values of CMARS reach 4.8 % and 8.5 for one day, and one week ahead forecasts, respectively. On the other hand, RR gives MAPE values of around 9.9% and 18.3%. We reveal that CMARS performs better than RR in terms of the main performance criteria. The rest of the paper is organized into four parts. In Sections 2 and 3, we provide a brief review of the models applied in this study. In Section 3, we present our mathematical models. The numerical results of the model are presented and discussed in Section 4. We conclude our study with a discussion of the results and giving future research in Section 5.

Ridge regression
The subset determination strategies include applying the least-squares to fit a linear model which involves a subset of the predictors. As an alternative to subset selection, we may fit a model including all k predictors utilizing a procedure that contracts the coefficient estimates towards zero, or equivalently, that compels or regularizes the coefficient estimates [4,14]. The shrinkage strategies become smaller the regression coefficients by implying a penalty on their size. These techniques bias the estimator of the regression coefficients to decrease the variance in addition to the mean squared error of the estimator and to forestall the model from overfitting [15]. The least-squares fitting procedure estimates 0 The coefficients of RR are the values that minimize the quantity  [4,15]. The penalty term is also known as 2 penalty and the 2 norm of a coefficient vector  is provided by

Conic multivariate adaptive regression splines
In general, with the CMARS algorithm, implementation of spline function acquires extreme and very important benefits especially in the modeling of dynamics. For both MARS and CMARS, considering a one-dimensional case (input variable), splines are piecewise polynomials. MARS obeys the following general model assumed to exist between the variables [6,15]: Here, each function is piecewise linear with a knot value,  . Therefore, in MARS, the Basis Functions (BFs) are determined as [15]   Here, (2) can be closer represented by a successively constructed linear combination of 58 A. Ozmen / IJOCTA, Vol. 12, No.1, pp.56-65 (2022) functions obtained from the set C and Y takes the following form:  + are coefficients that are estimated by least-squares with all other M+1 coefficients. Then the "winning" products are added into the model and then, for example, the possible candidates for BFs as follows [17]: x is in the model already, The model obtained by forward stage overfits the data. Therefore, a backwards-pruning procedure is applied to find BFs that contribute least to model fit, and then these BFs are progressively deleted. This iterative process continues until an optimal number of terms are presented in the final model [18,19]. Here, MARS used a modified recursive partitioning strategy to simplify high-dimensional problems. The sequence of models generated from this procedure is evaluated by using generalized cross-validation (GCV), and the model with the best predictive fit is selected. Consequently, an estimated best model f  of each number of terms,  , is found at the end of that process. In MARS model, generalized cross-validation is applied to define the optimal number of terms, ,  and it also shows the lack of fit. The GCV criterion defined by Friedman [4] is given by . (1 1)-M + parameter vector to be estimated by using the data points [16]. Furthermore, m Q is some appropriate integration domains in m K -dimensional parallel-pipe. After a discretization is employed to approximate the multi-dimensional integral Here, L is to be assigned an ( Then, the PRRS problem looks like a classical RR problem with 0   , 2  = for some ,   and this problem in (6) is expressed as a CQP [21]. So, based on a suitable selection of the bound , K the optimization problem in (6) can be rearranged [16]: At this point, we state that a careful learning process has to be followed for the choice of K .

Natural gas demand modeling by RR and CMARS
In this study, we consider two multi-objective models which use 2 norm regularization in linear and nonlinear modeling for the prediction of natural gas demand. Here, the advantage of regularization is to decrease the risk of overfitting, which usually occurs in high-dimensional learning. The primary goal of the regularization technique is to make the machine learning algorithm "learn" but "not memorize" by penalizing the algorithm to decrease its generalization error to avoid the risk of overfitting. As a consequence, the variance of the model may be considerably declined, without losing any important properties in the data. Moreover, since regularization is a kind of robustification, these kinds of models can also be called robust models.

Data
In this study, the data set comprising of daily natural gas demand data from 2004 to 2013 are provided by Baskentgaz, the LDC of Ankara. To check the performance of our models, we draw on the validation technique. We divide the dataset into two subsets as training and testing sets. Here, the dataset is not divided randomly since it includes a time series of natural gas and meteorological variables. Instead, the first six years (from 2004 to 2009) of each variable under consideration are selected as the training dataset, while the last four years of the series are selected as the test dataset. In this study, demand for residential customers of Baskentgaz is modeled by applying RR and CMARS algorithms. Our dataset contains some meteorological variables namely heating degree day (average temperature), relative humidity, wind speed, daily maximum and minimum temperatures. in this study, the meteorological variables used are proved by the Turkish State Meteorological Service. Indeed, energy consumption, especially natural gas, is highly dependent on weather conditions. If the temperature drops below a certain value of the heating threshold, households use more energy owing to the excessive need for space heating [1,2]. Here, firstly, the model details are given and then their performances are discussed for each forecasting time horizon. In the application of RR, the MATLAB regularization toolbox is utilized. In applications of CMARS, Salford MARS [19] is utilized to obtain BFs for the large model provided by the forward MARS algorithm. Afterward, MOSEK optimization software and MATLAB are used to solve the CQP problem and estimate unknown parameters.

Criteria for performance evaluations
Here, we tested the prediction accuracy of two specific prediction methods on real-time data. The main performance indicator for checking the accuracy of the models is the mean absolute percentage error (MAPE). In addition to MAPE, we also evaluate the multiple coefficient of determination ( 2 R ), correlation coefficient (r), average absolute error (AAE), and root mean square error (RMSE) to check the performances of proposed models. These measures and their formulas are presented in Table 1.

RR models
In the RR algorithm, first, many different models were obtained based on different penalty parameters,  , by using the MATLAB regularization toolbox. Here, the penalty parameter, , controlled the relative effect of both criteria (bias and variance) on the estimation of regression coefficient in (1). Then, the penalty parameter that tries to minimize two criteria in a balanced manner was selected. So, the following RR model based on selected penalty value provided the best solution for complexity and accuracy in (1).

One-day-ahead prediction
Day-ahead forecast is generally used to reduce operational costs and eliminate the drawbacks that can occur owing to the imbalance between supply and demand quantities. The system operator requires to identify the supply security issues that may exist the next day. The RR model for one-day-ahead forecasting can be represented as follows: Here, for RR and CMARS models based on one-day ahead and one-week ahead, X1, X2, X3, X4, X5, X6, X7, and X8 represents the first-order lagged, second-order lagged, third-order lagged, fourth-order lagged, fifthorder lagged, sixth-order lagged, seventh-order lagged and fourteenth-order lagged natural gas consumption, respectively. Hence, the RR and CMARS models obtained are "lag" models since they involve lagged dependent variables. Furthermore, X9, X10, X11, and X12 are the heating degree days, maximum temperature, minimum temperature, and wind speed, respectively.

One-week-ahead prediction
One-week-ahead forecast is usually used for generating unit production schedules for the next week. This gives companies insight into generation and consumption trends as well as assisting them to plan their weekly generation schedules. The RR model for a one-week forecast can be stated as follows: Here, RR models find all variables to be significant as expected.

CMARS Models
In the CMARS algorithm, first of all, using original data, the MARS algorithm is run via Salford MARS [22] and many MARS models were developed changing the number of BF and interaction. After selecting the optimal model for MARS among obtained models, the set of BFs in the forward part of MARS in the optimal MARS model was taken for the CMARS model. After the greatest models are constructed with the selected BFs and the L matrices in (7) are obtained, the PRSS in (6) is reformulated as a CQP problem. Here, based on different K values in (7), several different CQP models in (7) were solved individually via MOSEK [11]. Here, MOSEK applies an interior point algorithm [12,13] to treat the CQP problems CMARS yields many solutions. Finally, the CQP model that has the minimum value of approximate PRSS was selected and the unknown parameters were estimated for the CMARS model. Here, this point tries to minimize the criteria, Here, we should note that the knot values of BFs in the CMARS model are selected differently but extremely close to the corresponding input data to avoid nondifferentiability in the optimization problem.

One-day ahead prediction
In the CMARS algorithm, using the MARS software [22], the highest degree of interactions and For the CMARS model based on one-day ahead forecasting, the first-order lagged natural gas consumption ( 1 X ), the sixth-order lagged natural gas consumption ( 6 X ), the heating degree days ( 9 X ), and minimum temperature ( 11 X ) are significant.

One-week ahead prediction
For this natural gas data set based on one-week ahead forecasting, max M is 27 and the highest degree of interaction is 2. So, the largest model constructed by the forward MARS algorithm includes the following BFs.
For the CMARS model based on one-week ahead forecasting, the first-order lagged natural gas consumption ( 1 X ), the third-order lagged natural gas consumption ( 3 X ), the sixth-order lagged natural gas consumption ( 6 X ), the seventh-order lagged natural gas consumption ( 7 X ), the fourteenth-order lagged natural gas consumption ( 8 X ), the heating degree days ( 9 X ), maximum temperature ( 10 X ), and wind speed ( 12 X ) are significant.

Results
In this study, we propose to model seasonal patterns and trends directly without using any transformations. As shown in Figures 1-8, the presented models capture the seasonality and trend of natural gas consumption. Moreover, with an additional plot, Figures 1-8 display the daily error (residual) time series for each model and let us analyze the residual dynamics. For natural gas demand prediction, the performance of the RR and CMARS models was evaluated using the performance measures presented in Table 1. For the training and test cases of the RR and CMARS models, Table 2 and Table  3 show the performance matrices for the forecasting of one-day-ahead and one-week-ahead, respectively.     For the one-day-ahead forecast, Table 2 shows the performances of the RR and CMARS models based on training and test. Depending on the results given in Table 2, the CMARS model outperforms the RR model in terms of all performance measures for training and testing data. For the one-week-ahead forecast, the performances of RR and CMARS models are compared in Table 3 in terms of training and test cases. Based on the results shown in Table 3, the CMARS model outperforms the RR model in terms of all performance measures for training and test cases. Here, the MAPE performance becomes worse as the prediction horizon increases. Based on residual, RR and CMARS have similar plots for one-day-ahead forecast whereas CMARS has better results than RR for the one-week-ahead forecast as in terms of training and test cases, as seen in Figures 1-8.
Here, the proposed models have very small residuals for one-day-ahead forecasting. However, when the forecasting horizon increases, the residuals of the models become worse as expected.
In this application, the most essential variables in the analyzed residential demand models are based on temperature since a huge volume of natural gas is needed for space heating in Ankara, Therefore, in cold weather, natural gas usage reduces virtually linearly as the temperature rises. As you can see in Figures 1-7, this event also constructs a large difference between the winter and summer periods of private natural gas demand. Among the presented models, the CMARS model significantly outperforms the RR model in terms of all the training and test case performance criteria presented in Table 1 based on the short-term forecast. Therefore, CMARS should be the preferred model for this particular problem based on the short-term natural gas forecast.

Conclusion
In this study, we present two innovative models to short-term natural gas forecasting problems in the energy market. For residential users of LDCs, we evaluate daily and weekly forecasts of natural gas demand with daily intervals. We produce out-ofsample forecasts and compare them to observed data in terms of test datasets to assess the models' prediction accuracy. The prediction accuracy is assessed for each forecast horizon using the performance criteria listed in Table 1.
The proposed models, as illustrated in Figures 1-8, can capture the natural gas demand trend and seasonal pattern. However, it is found that CMARS performs better for one-day-ahead and one-week-ahead forecasts with MAPE values of 4.8% and 9.9%, respectively. Moreover, using CMARS as a substitute for MARS provides an integrated representation of all parameter identification tasks as a model-based optimization model instead of a model-free problem. As a result, in the multicriteria regression models offered, CMARS models should be favored for this particular situation. However, since the knot selection does not require for the linear parts in CMARS, as further study, Conic Generalized Partial Linear Model (CGPLM) [23] may also be added to the CMARS algorithms so that these semi-parametric algorithms may provide to decrease the complexity of CMARS. Moreover, the same study may be repeated with the robust counterparts being R(C)PLM [24,25] and R(C)MARS [26][27][28][29][30][31], and then the comparative results may be provided. In future researches, proposed models can be compared to time series models and many other traditional and recent methods.

Acknowledgments
The authors would like to thank the anonymous reviewers for their valuable comments. This work was supported by TUBITAK under the 2219-International Postdoctoral Research Scholarship Programme.