Ensemble Machine Learning Approaches for Prediction of Türkiye's Energy Demand

1. Introduction

Energy plays a vital role in supporting social and economic development of a country from past to present. It drives economic growth, improves living standards, supports social services, enhances national security, and contributes to environmental sustainability. Ensuring reliable, affordable, and sustainable energy access is essential for a country's overall development and progress. Therefore, the development of energy policies and the estimation of energy demand are a most important priority for developed and developing countries. Türkiye is the 19th largest economy in the world, with a gross domestic product (GDP) of roughly $906 billion [1]. Therefore, as a fast-growing economy, Türkiye's energy consumption has undergone significant growth and diversification over the years. As a rapidly developing country, Türkiye has experienced a substantial increase in energy demand due to population growth, urbanization, industrialization, and economic expansion. In terms of energy sources, Türkiye has a mix of fossil fuels, renewable energy, and imports. Electricity consumption constitutes a significant portion of Türkiye's energy consumption. However, the energy production in Türkiye is rather low in spite of the considerable increase in energy consumption.

In 2021, Türkiye produced approximately 41.3 MTOE (million tons of oil equivalent), mostly based on coal and lignite (57.5 %). In contrast, the country's energy consumption reached approximately 109.8 MTOE in the same year. The substantial gap between energy production and consumption led to Türkiye becoming one of the major energy importing countries in Europe [2]. As a report of Ministry of Energy and Natural Resources in September 2021, Türkiye's energy dependency rate was reported to be around 74%. This means that Türkiye relied on imports to meet approximately 74% of its total energy consumption. The country heavily depended on imports of natural gas and oil to bridge the energy gap between domestic production and demand [3].

To provide energy either by importing or by producing, forecasting energy consumption, analyzing the relationship between energy demand and supply are crucial issues in short and long term energy planning. Managing energy demand also involves identifying and prioritizing energy resources, optimizing energy utilization, improving energy efficiency, shaping policy decisions, and devising strategies to reduce emissions.

Furthermore, it is important to emphasize that the United Nations' SDGs provide a comprehensive blueprint for addressing global challenges and promoting sustainability by 2030 [4]. This study aligns closely with several of these goals: including Goal 7, Goal 8, and Goal 13, making a significant contribution to the broader aims of sustainable development. Accurate energy demand forecasting plays a central role in achieving 'Goal 7: Affordable and Clean Energy.' By optimizing energy production, distribution, and consumption, this study facilitates the provision of affordable, reliable, and clean energy. This, in turn, supports economic growth, enhances energy access, and reduces environmental impacts. Moreover, this study directly addresses 'Goal 13: Climate Action' by mitigating the effects of climate change. Through precise energy demand predictions, it empowers Türkiye to make informed decisions that reduce greenhouse gas emissions, promote renewable energy adoption, and foster a low-carbon, sustainable energy sector. Energy efficiency serves as a catalyst for economic growth and job creation, aligning with 'Goal 8: Decent Work and Economic Growth.' This study enables Türkiye to implement energy-efficient measures, leading to cost savings for industries and households.

Researchers have developed various statistical techniques, meta-heuristic algorithms, artificial intelligence techniques in energy modeling. Artificial Neural Networks (ANNs), have gained significant interest in energy planning due to their ability to handle complex nonlinear relationships between input and output data [5]. ANNs have been applied in various energy forecasting applications, including gas consumption [6], energy demand [7], electricity consumption [8], transportation energy demand [9,10,11] energy source analysis [12] and energy dependency [7]. Apart from ANNs, other prediction methods have emerged, such as Fuzzy Logic, Adaptive Network-based Fuzzy Inference System (ANFIS), and general Machine Learning algorithms [13,14,15]. It's important to recognize that Artificial Neural Network (ANN) is a subfield of Machine Learning (ML), which, in itself, is a subset of Artificial Intelligence (AI). Frequently, the terms AI, ML, and Deep Learning (DL) are used interchangeably to refer to intelligent systems or software. DL, specifically, extends the concept of ANNs by incorporating extra hidden layers and employing specialized activation functions that aren't typically found in traditional ANN models.

The AI-based prediction models have received considerable interest to solve a variety of problems in energy planning recently. These models leverage the power of artificial intelligence techniques to forecast future outcomes based on historical data patterns. These models use advanced algorithms and machine learning methods to analyze large datasets, identify patterns, and make predictions.

Research on predicting Türkiye's energy requirements began in the 1960s, with the State Planning Organization (SPO) employing basic regression techniques for energy forecasting. In the late 1970s, The Ministry of Energy and Natural Resources (MENR) and the Turkish Statistical Institute (TSK) started preparing energy demand projections [16], but the estimated values provided by MENR were found to be higher than the actual energy demand [17]. Numerous econometric modeling techniques were applied to forecasting energy consumption after 1984. the Model for Analysis of Energy Demand (MAED) is most frequently used approach developed by MENR [18]. Nevertheless, the energy demand predictions generated by MAED continued to overstate actual demand, rendering them unreliable [19,20]. Utgikar and Scott [21] conducted an inquiry to understand the reasons behind unsuccessful energy predictions. They found that although statistical models are favored by researchers due to their simplicity, they tend to deliver acceptable results only for short-term periods, while becoming increasingly unstable for longer-term forecasts. While statistical models have contributed valuable insights to energy forecasting studies, El-Telbany and El-Karmi [22] argue that these models, which rely on statistical methods, perform well under normal conditions but struggle to account for sudden changes in environmental or sociological variables. Due to the above discussed reasons, researchers have garnered the interest in the field of energy demand modeling.

AI-based prediction methods utilize advanced algorithms and models that learn patterns, relationships, and trends from historical data to make accurate predictions about future events or outcomes. AI-based prediction methods have gained significant popularity due to their ability to handle complex, non-linear relationships and adapt to changing data patterns. Machine Learning (ML) is an AI-based prediction method, which encompasses a range of algorithms that automatically learn and improve from data without being explicitly programmed. Deep Learning (DL), a subset of ML, has emerged as a powerful AI-based prediction method.

AI-based prediction methods have demonstrated their effectiveness in various fields, including finance, healthcare, weather forecasting, sales forecasting, demand prediction, and fraud detection [23,24,25,26]. However, it's important to note that when applying these methods to real-world scenarios, Factors such as data quality, feature selection, model complexity, and interpretability should be carefully considered. Machine learning approaches can be broadly categorized into three main types: supervised learning, unsupervised learning, and reinforcement learning. Supervised learning algorithms include Classification and Regression. The regression is used for energy demand in this work. In machine learning, ensemble learning and deep learning methods outperform traditional algorithms. Ensemble methods are learning algorithms that build a set of classifiers and then classify new data points by taking (weighted) votes of their predictions. The effectiveness of an ensemble method depends on several factors, including how the underlying models are trained and how they are combined. In the literature, there are common approaches to building ensemble models that have been successfully demonstrated in various domains [27,28,29].

Ensemble machine learning is a powerful technique in the field of machine learning that involves combining the predictions of multiple individual models (base models) to create a more accurate and robust predictive model. This approach, preferred over single methods, offers several key advantages. Firstly, ensembles enhance prediction accuracy by aggregating multiple models, reducing errors, and biases. Secondly, they mitigate overfitting, a common issue in machine learning, by balancing out individual model weaknesses. Moreover, ensembles prove their robustness by effectively handling noisy data and outliers, making them suitable for real-world applications [30]. In automated decision-making applications, especially in engineering, ensemble methods have demonstrated superior performance compared to individual learners. This is attributed to their ability to capture diverse patterns, reduce bias and variance, and improve generalization. Ensemble methods are particularly effective when there is a large amount of data, complex relationships, and a need for high predictive accuracy.

Common ensemble method strategies comprise bagging, boosting, blending, and stacking. Bagging, exemplified by the Random Forest algorithm, enhances model robustness by reducing overfitting and improving prediction accuracy through the wisdom of the crowd [31]. Boosting is another ensemble technique that iteratively builds a strong predictive model by giving more weight to the data points that previous models misclassified [32]. Bagging reduces variance by averaging over multiple models, boosting focuses on reducing bias through weighted data points. Blending, sometimes referred to as model stacking or meta-ensembling, involves training multiple diverse base models on the same dataset and then combining their predictions using a separate model trained on the validation set. Stacking, similar to blending, combines multiple base models to form a meta-model but differs in its approach. In stacking, the predictions of the base models serve as input features for a meta-model, which learns to make the final predictions [33]. Blending combines diverse models with a separate meta-model, and stacking uses base models to create a meta-model for predictions. The choice among these ensemble methods depends on the specific problem, dataset, and the trade-off between bias and variance in the model [34].

The main objective of this study is to use ensemble machine learning methodologies, which have not received much attention in prior research on energy, to assess energy demand in Türkiye. In this paper, several significant contributions were presented, including empirical evidence substantiating the efficacy of ensemble methods in forecasting Türkiye's energy demand, comprehensive hyperparameter tuning, and an exhaustive examination of 19 distinct machine learning algorithms along with the application of 5 performance metrics to evaluate their predictive capabilities. To the best of our knowledge, this paper is the first to investigate Türkiye's energy demand using ensemble machine learning models.

This paper is organized as follows: Section 2 will provide an overview of the scope and definition of energy demand studies. Section 3 will present the primary methods and approaches employed in energy demand study, along with the main data sources and challenges associated with energy. In Section 4, the principal findings and trends will be discussed using various ensemble learning algorithms. Section 5 will summarize the main implications and recommendations from the results of energy demand study. Finally, our paper will conclude with a discussion of limitations and directions for future research.

2. Literature Review

Energy demand forecasting is an important task for planning and managing energy systems. It involves predicting the future energy consumption of different sectors, regions, or appliances based on various factors such as weather, economic activity, population, lifestyle, etc. This literature review aims to summarize some of the key findings and trends from recent articles on this topic.

2.1. Review of Energy Demand Forecasting in the World

The global energy crisis of 2022 has brought renewed attention to the challenges and opportunities of energy demand in a rapidly changing world. According to various sources, such as the World Energy Outlook 2022 by the International Energy Agency (IEA), the Outlook for Energy Demand by McKinsey, and several articles from reputable media outlets, energy demand is undergoing a major transformation due to a combination of factors, such as economic slowdown, high energy prices, energy security concerns, and climate policies.

The global energy demand has been affected by the Covid-19 pandemic and the economic recovery in 2021. According to the Global Energy Review 2021 by the International Energy Agency (IEA), the global energy demand is expected to grow by 4.6% in 2021. The IEA projects that global energy demand will grow by 0.8% per year on average between 2021 and 2030 in its Stated Policies Scenario (STEPS), which reflects current and announced policies and targets [35]. This literature review has provided a brief overview of some of the main findings and trends from recent articles on energy demand.

There are many methods for energy demand forecasting, spanning from conventional approaches like econometric and time series models to contemporary soft computing techniques, including artificial intelligence methods and heuristics.

Zhang et al. [36] aimed to predict the transport energy demand in China for the next two decades using a statistical method called partial least square regression (PLSR). The study estimated that by 2020, China would need about 433 to 468 million tons of coal equivalent (Mtce) of energy for its transportation sector, mainly from petroleum products. The study demonstrated the usefulness of PLSR as a simple and reliable technique for forecasting transport energy demand.

Kumar and Jain [37] presented three time series models for forecasting the consumption of conventional energy in India. The study assessed the accuracy of these forecasts by comparing them with both actual data and projections made by the Planning Commission of India.

Chaturvedi et al. [38] evaluated four time-series models for predicting total and peak monthly energy demand in India. The study compared the performance of the Seasonal Auto-Regressive Integrated Moving Average (SARIMA), Long Short Term Memory Recurrent Neural Network (LSTM RNN), and Facebook Prophet (Fb Prophet) models in forecasting Indian energy demand. Among these models, Fb Prophet performed well for both total and peak demand, exhibiting lower prediction errors compared to the others.

Sahraei et al. [39] presented the transportation energy demand in Türkiye using a statistical method called Multivariate Adaptive Regression Splines (MARS), which is a nonparametric regression technique. They established five models and compared with real data which are collected from MENR. The approach employed in that study illustrates the utility of the MARS method in examining and forecasting energy consumption within the transportation sector.

Machine learning algorithms have been applied to predict energy consumption in various studies. Javanmard and Ghaderi [40] studied about predicting energy demand in Iran until 2040 using machine learning and optimization algorithms. The authors collect data from different energy sectors and apply six machine learning algorithms (ARIMA, SARIMA, SARIMAX, ANN, AR, and LSTM) and mathematical programming to forecast the energy demand. Then, they use an integrated model with Particle Swarm and Grey-Wolf optimization algorithms to improve the prediction accuracy. They compare the results of the proposed method with the machine learning algorithms and report a high accuracy for the proposed method. The study also estimated the growth rate of energy demand in each sector.

Ye et al. [41] employed the discrete grey model (DGM) to enhance the accuracy of energy consumption predictions. Meanwhile, Zhang and Wang [42] introduced a fuzzy wavelet neural network (FWNN) method to forecast annual electricity consumption in a city known for its high energy consumption in China.

El-Telbany and El-Karmi [22] presented and discussed the results of applying ANN to forecast the energy demand using PSO method, which was a novel adaptive algorithm based on a social-psychological metaphor. They compared the performance of this method with the performance of BP algorithm and autoregressive moving average technique.

Mason et al. [43] utilized the covariance matrix adaptation evolutionary strategy, a powerful evolutionary optimization algorithm, to train neural networks for predicting Ireland's energy needs. The neural network trained with the covariance matrix adaptation evolutionary strategy performed well compared to other state-of-the-art prediction methods.

Muratlıharan et al. [44] introduced a novel approach for energy demand prediction by utilizing neural networks in an optimization framework. Initially, the Conventional Neural Network (CNN) approach was employed to predict energy demand at the consumer level. Then, two variations, namely the Neural Network based Genetic Algorithm (NNGA) and Neural Network based Particle Swarm Optimization (NNPSO), were proposed to automatically adjust the neural network weights. The results indicated that the NNGA approach outperforms in short-term load forecasting, while the NNPSO approach is more suitable for long-term energy prediction.

Yu et al. [45] introduced a hybrid algorithm known as the Particle Swarm Optimization and Genetic Algorithm Optimal Energy Demand Estimating (PSO–GA EDE) model to estimate China's energy demand. Their study demonstrated the superiority of this model when compared to alternative approaches.

In general, statistical technique and artificial intelligence can be used for solving the real-world problems. Statistical techniques involve mathematical models, including regression, multiple regression models, Autoregressive Moving Average, and Kalman filters [39,41]. Evolutionary algorithms, such as genetic algorithms and particle swarm optimization, are also employed. Artificial intelligence approaches involve a combination of fuzzy logic, artificial neural networks, genetic algorithms, support vector machines, and machine learning methods [46].

A systematic literature review encompassing 419 articles on energy demand modeling, covering the period between 2015 and 2020, was conducted by Verwiebe et al. [46]. They analyzed the methodologies, prediction accuracy, input variables, energy sources, sectors, temporal scopes, and spatial resolutions employed in these models. They found that machine learning techniques were the most used, followed by engineering-based models, metaheuristic and uncertainty techniques, and statistical techniques. They also discussed the drawbacks and countermeasures of each technique. Another systematic literature review of energy demand forecasting methods published in 2005–2015 was conducted by Ghalehkhondabi et al. [47]. They focused on the methods that are used to predict energy consumption and compared their performance and applicability. They reported that neural networks were the most cited technique and had notable performance, but also high computation time. They suggested that hybrid methods could be a promising field of future research.

2.2. Energy Demand Forecasting in Türkiye

A summary of studies for the Türkiye's energy demand forecasting is tabulated in Table 1. However, to the best of our knowledge, there is no research paper that employs ensemble machine learning methods and compares them with each other to forecast Türkiye's energy demand.

3. Materials and Methods

In this section, proposed methodology is introduced in detail. Ensemble methods refer to algorithms that combine multiple machine learning models into a unified framework. These methods have gained significant attention and recognition in the machine learning community due to their ability to enhance prediction accuracy and robustness [41,42]. By combining the predictions of multiple models, ensemble methods can mitigate the limitations of individual models and provide more accurate and reliable results. Several types of ensemble methods commonly used in machine learning are Bagging, Boosting, Blending, Random Forest and Stacking [31,32,33,34]. This paper proposes and analyzes different ensemble combination models that can be achieved by using diverse base models, varying model architectures, or training on different subsets of the data.

3.1. ML Algorithms

In the context of forecasting Türkiye's energy demand, 19 machine learning algorithms have been selected that have demonstrated excellent performance in energy consumption forecasting problems. These algorithms are briefly described below.

• Light Gradient Boosting Machine (LightGBM) [72]

LightGBM is a popular machine learning algorithm used for both regression and classification tasks. It is designed to efficiently handle large-scale datasets with high-dimensional features. LightGBM is known for its speed, accuracy, and ability to handle complex problems. LightGBM is based on the gradient boosting framework, similar to other boosting algorithms. However, it incorporates several optimizations that make it more efficient and scalable. One key optimization is the implementation of a novel histogram-based algorithm for binning and feature discretization. This approach reduces the memory usage and speeds up the training process by avoiding the need to sort data during each iteration. Another important feature of LightGBM is the implementation of a leaf-wise tree growth strategy. Unlike traditional depth-wise growth, where trees are expanded level by level, LightGBM grows trees’ leaf by leaf. This approach reduces the number of levels in the trees, resulting in faster training and improved accuracy.

• XGBoost [73]

XGBoost Regressor is a powerful machine learning algorithm used for regression tasks. XGBoost Regressor is known for its efficiency, accuracy, and ability to handle complex datasets. The algorithm minimizes a loss function by iteratively adding decision trees to the ensemble. Each tree is trained to predict the residuals (the differences between the actual and predicted values) of the previous ensemble. The process continues until a specified number of trees is reached or the desired level of performance is achieved. To further improve the performance, XGBoost Regressor incorporates several advanced techniques. It employs a regularization term in the objective function to control the complexity of the model and prevent overfitting. It uses a technique called gradient boosting with shrinkage, which adds a weighting factor to each tree's contribution, allowing for more conservative updates.

XGBoost Regressor offers several advantages. It is highly customizable, allowing users to fine-tune various hyperparameters to achieve optimal performance. It can handle large-scale datasets with a high number of features and is computationally efficient due to parallelization and tree pruning techniques.

• Extra Tree Regression [74]

Extra Trees Regression is a machine learning algorithm used for regression tasks. It belongs to the ensemble learning family and is an extension of the popular Random Forest algorithm. Extra Trees Regression combines multiple decision trees to make predictions by aggregating their outputs. The algorithm builds a user-defined number of decision trees using random subsets of the training data and random subsets of features. Each tree is grown to its maximum depth, resulting in a collection of unpruned and highly random trees. The final prediction is then obtained by averaging or aggregating the predictions made by each tree.

Extra Trees Regression offers several advantages. First, it reduces overfitting since the randomization helps to capture different aspects of the data. It can handle high-dimensional datasets with a large number of features and works well with both continuous and categorical variables. Additionally, it is computationally efficient as the trees can be grown in parallel.

• Passive Aggressive Regressor (PAR) [75]

(PAR) is a machine learning algorithm used for regression tasks. In Passive Aggressive Regressor, the algorithm updates the regression model incrementally, making predictions on new instances as they arrive. It adapts to new data points by adjusting the model's parameters without revisiting the entire training set. This property makes it suitable for handling large-scale datasets or scenarios where data arrives in a streaming fashion.

Passive Aggressive Regressor is particularly useful when dealing with dynamic or non-stationary data, as it can update the model in an online fashion. It is also memory-efficient since it does not require storing the entire training set.

• Elastic Net [76]

Elastic Net is a regression method that combines the Lasso and Ridge regression techniques. It is used for feature selection and regularization in linear models, providing a balance between the two methods. In Elastic Net, the algorithm aims to minimize the sum of squared residuals between the predicted and actual values, similar to ordinary least squares (OLS) regression. However, Elastic Net adds a penalty term to the cost function that combines both the Lasso (L1) and Ridge (L2) penalties. The penalty term is a linear combination of the absolute values of the coefficients (L1) and the squared values of the coefficients (L2), each multiplied by their respective tuning parameters. The combination of L1 and L2 penalties in Elastic Net allows for both feature selection and regularization.

• Least Angle Regression (LARS) [77]

LARS is a regression method used for feature selection and model building. It is an algorithm that selects predictors in a direction that minimizes the prediction error and correlates strongly with the target variable. LARS starts with an empty set of selected features and gradually adds features in a way that balances their correlations and coefficients. In each iteration, LARS identifies the most correlated feature with the residual and increases the coefficient of that feature along its equiangular direction. The algorithm continues this process until it reaches the desired number of selected features or the maximum number of available features.

• Lasso Least Angle Regression [78]

Lasso Least Angle Regression is a regression method that combines the features of the Lasso regularization and the Least Angle Regression algorithm. It is used for feature selection and regularization in linear regression tasks. Lasso Least Angle Regression aims to estimate the coefficients of a linear regression model while simultaneously performing feature selection by encouraging sparsity in the solution. It starts with an empty set of selected features and gradually adds features in a way that balances their correlations and coefficients.

This algorithm provides a continuous path of solutions for different regularization levels, allowing for more flexibility in choosing the amount of regularization. It is particularly useful when dealing with high-dimensional datasets where feature selection and regularization are crucial to prevent overfitting and improve interpretability.

• Orthogonal Matching Pursuit (OMP) [79]

Orthogonal Matching Pursuit (OMP) is an algorithm used for sparse signal recovery and feature selection tasks. It aims to find the most relevant features or components of a signal by iteratively selecting and reconstructing the signal based on a small subset of measurements or features. In OMP, the algorithm starts with an initial guess of the solution and iteratively improves it by selecting the most correlated features with the residual error. In each iteration, the algorithm identifies the feature that has the highest correlation with the residual and adds it to the solution. It then updates the residual by subtracting the contribution of the selected feature from the original signal. This process continues until a predefined stopping criterion is met, such as a desired level of approximation or a fixed number of selected features.

OMP leverages the orthogonality property to efficiently select features and estimate the signal. At each iteration, the algorithm ensures that the selected features are orthogonal or nearly orthogonal to each other, which helps in accurate signal reconstruction and efficient convergence.

• Random Forest Regressor [80]

Random Forest Regressor is a popular machine learning algorithm used for regression tasks. It belongs to the ensemble learning family and is built upon the concept of decision trees. In Random Forest Regressor, a user-defined number of decision trees are constructed. Each tree is built using a random subset of the training data and a random subset of features. The process of constructing each tree involves recursively splitting the data based on different features and their respective splitting points. The splitting is done in a way that minimizes the variance of the target variable within each resulting subset. During prediction, each tree in the random forest independently generates a prediction. The final prediction is obtained by averaging or taking the majority vote (in the case of classification) of the predictions from all the trees.

Random Forest Regressor offers several advantages. It can handle high-dimensional datasets with a large number of features and works well with both continuous and categorical variables. Additionally, it provides a measure of feature importance, allowing for better understanding of the underlying relationships in the data.

• Gradient Boosting Regressor [81]

Gradient Boosting Regressor is a powerful machine learning algorithm used for regression tasks. It belongs to the boosting family of algorithms and is designed to create a strong learner by iteratively combining weak learners. Gradient Boosting Regressor works by minimizing a loss function through an additive approach, where each new model is built to correct the errors made by the previous models. In Gradient Boosting Regressor, the algorithm starts by creating an initial model, often a simple one like a decision tree. This model makes initial predictions on the training data. The subsequent models are then built to minimize the differences between the predicted and actual values of the previous models. The new models are trained on the negative gradients of the loss function with respect to the previous ensemble's predictions. This means the new models learn to predict the residuals or errors made by the existing ensemble. The final prediction is obtained by aggregating the predictions from all the models in the ensemble.

Gradient Boosting Regressor offers several advantages. It can handle complex datasets and capture intricate relationships between features. It is also capable of handling various loss functions, making it adaptable to different regression problems. The algorithm is robust against overfitting due to regularization techniques employed during the training process.

• AdaBoost Regressor [82]

AdaBoost Regressor, short for Adaptive Boosting Regressor, is a machine learning algorithm used for regression tasks. The algorithm iteratively trains a series of weak regressors, each focusing on the instances that were mis predicted by the previous regressors, to improve the overall prediction accuracy. In AdaBoost Regressor, each weak regressor is trained on a subset of the training data. During training, the algorithm assigns weights to each instance, with initially equal weights for all instances. In each iteration, a weak regressor is trained on the weighted training data, and the weights are adjusted based on the performance of the regressor. Instances that are mispredicted by the current regressor are given higher weights, while correctly predicted instances are given lower weights. This emphasizes the importance of the mispredicted instances in subsequent iterations. During prediction, the weak regressors are combined by assigning weights to their predictions. The final prediction is obtained by aggregating the weighted predictions of all the weak regressors. The weights are determined based on the performance of each regressor during training. The regressors with better performance are given higher weights, indicating their greater influence on the final prediction. AdaBoost Regressor has several advantages. It can handle complex relationships in the data and is resistant to overfitting. By focusing on mispredicted instances, it iteratively improves the model's ability to capture difficult patterns in the data.

• Linear Regression [83]

Linear regression is a statistical method that models the relationship between a dependent variable (y) and one or more independent variables (x). It can be used to estimate how the dependent variable changes as the independent variables change, and to test hypotheses about the strength and direction of the relationship. There are different types of linear regression, such as simple linear regression (one independent variable), multiple linear regression (more than one independent variable), and multivariate linear regression (more than one dependent variable). Linear Regression is a widely used statistical modeling method for predicting a continuous target variable based on one or more predictor variables. It assumes a linear relationship between the predictors and the target variable, where the goal is to find the best-fitting line or hyperplane that minimizes the difference between the predicted and actual values.

• Lasso Regression [84]

Lasso Regression (LASSO) is a method of regression analysis that performs both variable selection and regularization. It aims to improve the prediction accuracy and interpretability of the regression model by shrinking the coefficients of some predictor variables to zero and reducing the magnitude of others.

In Lasso Regression, the algorithm aims to minimize the sum of squared residuals between the predicted and actual values, similar to ordinary least squares (OLS) regression. However, Lasso Regression adds a penalty term to the cost function, known as the Lasso or L1 penalty. This penalty term is the sum of the absolute values of the coefficients multiplied by a tuning parameter called the regularization parameter or lambda (λ).

• K Neighbors Regressor [85]

K Neighbors Regressor is a machine learning algorithm used for regression tasks. It is a non-parametric method that predicts the target value of an instance by considering the average or weighted average of the target values of its k nearest neighbors in the training data. In K Neighbors Regressor, the algorithm identifies the k nearest neighbors of a given instance based on a distance metric, such as Euclidean distance. The target values of these neighbors are then used to calculate the predicted value for the instance. In the case of average-based prediction, the predicted value is the average of the target values of the k neighbors. In the case of weighted average prediction, the predicted value is the weighted average, where the weights are determined based on the proximity of the neighbors to the instance. The value of k in K Neighbors Regressor determines the number of neighbors to consider in the prediction.

K Neighbors Regressor offers simplicity and flexibility. It can handle both numerical and categorical features and can capture complex relationships in the data.

• Bayesian Ridge Regression [86]

Bayesian Ridge Regression is a regression method that incorporates Bayesian principles into the linear regression framework. In Bayesian Ridge Regression, the algorithm places a prior distribution on the regression coefficients, typically assuming a Gaussian distribution. This prior distribution represents the initial belief about the likely values of the coefficients before observing the data. As data is observed, the prior distribution is updated to form the posterior distribution, which represents the updated belief about the coefficient values given the observed data.

The posterior distribution is obtained using Bayes' theorem, which combines the prior distribution with the likelihood function of the data. The likelihood function expresses the probability of observing the data given the regression coefficients. The posterior distribution captures the uncertainty in the coefficient estimates and can be used to compute credible intervals or conduct hypothesis tests.

• Decision Tree Regressor [87]

Decision Tree Regressor is a machine learning algorithm used for regression tasks. It is based on the concept of a decision tree, which partitions the input space into regions and predicts the target value based on the average or majority value of the training instances within each region. In Decision Tree Regressor, the algorithm recursively splits the data based on different features and their respective splitting points to create a tree-like structure. The splitting is done in a way that minimizes the variance or mean squared error of the target variable within each resulting subset. Each internal node of the tree represents a splitting point, and each leaf node represents a prediction for a subset of data.

Decision Tree Regressor offers several advantages. It is intuitive and easy to interpret, as the resulting tree structure can be visualized and understood. It can handle both numerical and categorical features and can capture complex nonlinear relationships in the data. Decision Tree Regressor is also robust to outliers and can handle missing values in the data.

• Ridge Regression [88]

Ridge Regression is a linear regression method used for modeling and prediction tasks. It is an extension of ordinary least squares (OLS) regression that introduces a regularization term to handle multicollinearity and prevent overfitting. In Ridge Regression, the algorithm seeks to minimize the sum of squared residuals between the predicted and actual values, similar to OLS regression. However, Ridge Regression adds a penalty term, known as the Ridge or L2 penalty, to the cost function. This penalty term is the squared sum of the coefficients multiplied by a tuning parameter called the regularization parameter or lambda (λ). The regularization parameter (λ) controls the degree of shrinkage applied to the coefficients. A larger λ value increases the amount of shrinkage, resulting in a more heavily regularized model with smaller coefficient values. Conversely, a smaller λ value reduces the amount of shrinkage, making the model behave more like ordinary least squares regression.

Ridge Regression offers several advantages. It is particularly useful when dealing with multicollinearity, where the predictor variables are highly correlated. It also provides a balance between bias and variance, helping to prevent overfitting in situations where the number of predictors is larger than the number of observations.

• Huber Regressor [89]

Huber Regressor is a robust regression method that combines the benefits of both least squares’ regression and robust regression techniques. Huber Regressor addresses these issues by introducing a hybrid loss function that behaves like least squares for small residuals and like a scaled absolute loss for large residuals. The loss function is controlled by a tuning parameter called the Huber loss threshold. Residuals below the threshold are squared, while residuals above the threshold are linearly scaled. This combination of squared and absolute loss enables Huber Regressor to provide a balance between robustness to outliers and accuracy for non-outlying data points. During training, Huber Regressor estimates the coefficients of the linear regression model by minimizing the modified loss function. This estimation process reduces the influence of outliers while still accurately capturing the patterns in the majority of the data.

• Dummy Regressor

Dummy Regressor is a simple baseline model used for regression tasks. It provides a straightforward way to establish a baseline performance against which other regression models can be compared. Dummy Regressor makes predictions based on simple rules or heuristics rather than learning patterns from the data. In Dummy Regressor, the predictions are determined based on a predefined strategy.

The Materials and Methods should be described with sufficient details to allow others to replicate and build on the published results. Please note that the publication of your manuscript implicates that you must make all materials, data, computer code, and protocols associated with the publication available to readers. Please disclose at the submission stage any restrictions on the availability of materials or information. New methods and protocols should be described in detail while well-established methods can be briefly described and appropriately cited.

Research manuscripts reporting large datasets that are deposited in a publicly available database should specify where the data have been deposited and provide the relevant accession numbers. If the accession numbers have not yet been obtained at the time of submission, please state that they will be provided during review. They must be provided prior to publication.

Interventionary studies involving animals or humans, and other studies that require ethical approval, must list the authority that provided approval and the corresponding ethical approval code.

3.2. Structure of the Proposed Methods

This section presents the structure and abstract overview of the study, with a basic conceptual flow shown in Figure 1. The methodology consists of several key steps to improve the accuracy of our models. The first step entails data preparation, encompassing data preprocessing, normalization, and transformation to ensure the dataset is primed for analysis. Following this, the performance of 19 different machine learning algorithms have been assesed. This evaluation forms the basis for creating various ensemble combinations that leverage the strengths of individual models. In the third step, four different ensemble techniques like bagging, boosting, blending, and stacking models have been used to to create powerful ensemble models that can capture complex patterns and relationships in the data. Finally, our study concludes with the fourth step, where it has carefully evaluated and compared these ensemble models, gaining insights into their strengths and weaknesses. This comprehensive evaluation process has informed our final prediction, guiding us towards data-driven decisions that hold the potential to advance the field of artificial intelligence and machine learning.

3.2.1. K-Fold Cross-Validation

The next step in the processing block involves selecting machine learning algorithms that exhibit superior performance and diverse learning capabilities when applied to our energy dataset. K-fold cross-validation is a technique used in machine learning algorithms to evaluate the performance of a model by dividing the available data into k equally sized subsets or "folds." The process entails iteratively training the model on k-1 folds and then evaluating it on the remaining fold. This cycle is repeated k times, with each fold being used as the test set exactly once. The final evaluation is obtained by averaging the performance results from each iteration. The value of k is usually chosen as k=5 or k=10. A 5-fold cross-validation is used to obtain groups of performance measures in this study. By using different subsets of data for training and validation, k-fold cross-validation helps to assess the model's ability to generalize well to unseen data. This process helps to mitigate issues such as overfitting and provides a more reliable estimate of our model's performance by averaging the results obtained from each iteration.

3.2.2. Model Hyperparameters Tuning

Hyperparameter tuning is the process of finding the optimal values for the hyperparameters of a machine learning algorithm. Hyperparameters are parameters that are set before the learning process begins and determine how the algorithm learns and generalizes from the training data. By selecting the right combination of hyperparameter values, we can improve the model's accuracy, generalization, and robustness [90].

In this study, hyperparameter tuning was employed to improve model performance and prevent overfitting before proceeding to the next stage of our framework. The commonly used approaches for hyperparameter tuning are Grid Search, Random Search, and Bayesian Optimization. Grid search involves defining a grid of hyperparameter values and searching through a predefined set of hyperparameter values for a given model. In random search, hyperparameters are randomly sampled from predefined distributions. Bayesian optimization uses probabilistic models to estimate the performance of different hyperparameter configurations. The choice of the technique ultimately depends on the specific problem, available computational resources, and the characteristics of the hyperparameter search space. In this work Grid Search approach has used.

3.2.3. Performance Metrics

The performance of ML algorithms in our energy demand problem was estimated using powerful validation techniques. Five validation methods, represented by equations (1)-(5), were employed to evaluate the models. Additionally, the regularization variables were updated on the approval set to optimize the forecasting performance. The data was divided into training and testing sets to assess the models' functionality.

In order to evaluate and compare the models, various standard predictive performance metrics are [10]:

Mean Square Error (MSE) calculates the average of the squared differences between predicted and actual values, which emphasizes larger errors and makes it sensitive to outliers. This metric is commonly employed in regression problems, and lower MSE values indicate better model performance.

$$\mathrm{M}\mathrm{S}\mathrm{E}=\sum _{\mathrm{i}=1}^{\mathrm{n}}\frac{{\left({\breve{\mathrm{y}}}_{\mathrm{i}}-{\mathrm{O}}_{\mathrm{i}}\right)}^2}{\mathrm{n}}$$

(1)

Root Mean Square Error (RMSE) is calculated as the square root of MSE and serves as a measure of the average magnitude of residuals, which represent the differences between predicted and actual values. It offers an overall assessment of prediction accuracy, where lower RMSE values correspond to better model performance.

$$\mathrm{R}\mathrm{M}\mathrm{S}\mathrm{E}=\sqrt{\sum _{\mathrm{i}=1}^{\mathrm{n}}\frac{{\left({\breve{\mathrm{y}}}_{\mathrm{i}}-{\mathrm{O}}_{\mathrm{i}}\right)}^2}{\mathrm{n}}}$$

(2)

Mean Absolute Error (MAE) calculates the average absolute difference between the predicted values and the actual values. It offers a measure of the average error magnitude, irrespective of their direction.

$$\mathrm{M}\mathrm{A}\mathrm{E}=\frac{\sum _{\mathrm{i}=1}^{\mathrm{n}}\left|{\breve{\mathrm{y}}}_{\mathrm{i}}-{\mathrm{O}}_{\mathrm{i}}\right|}{\mathrm{n}}$$

(3)

R-squared (R2) measures the proportion of the variance in the dependent variable that's predictable from the independent variables in a regression model. It ranges from 0 to 1, with higher values indicating a better fit. However, R² can be misleading when overfitting occurs.

$${\mathrm{R}}^2=1-\frac{{\sum }_{\mathrm{i}}{\left({\breve{\mathrm{y}}}_{\mathrm{i}}-{\mathrm{O}}_{\mathrm{i}}\right)}^2}{{\sum }_{\mathrm{i}}\left({\mathrm{O}}_{\mathrm{i}}-{\overline{\mathrm{O}}}_{\mathrm{i}}\right)2}$$

(4)

Mean Absolute Percentage Error (MAPE) calculates the percentage difference between predicted and actual values, averaging these percentages across all data points. MAPE can help assess the model's accuracy in terms of percentage errors.

$$MAPE=\frac{1}{m}\sum _{i=1}^m\left|\frac{{\breve{\mathrm{y}}}_{\mathrm{i}}-{\mathrm{O}}_{\mathrm{i}}}{{\breve{\mathrm{y}}}_{\mathrm{i}}}\right|$$

(5)

Here, $\left({\overline{\mathbf{O}}}_{\mathbf{i}}\right)$ represents the magnitude of the actual values, $\left({\breve{\mathbf{y}}}_{\mathbf{i}}\right)$ is the model's predicted value, (Oi) stands for the real data, and (n) indicates the number of observed data points.

3.3. Data Collection

The dataset used in this paper includes independent variables such as population (in millions), gross domestic product (GDP), import, and export, which were selected based on a comprehensive literature review. This dataset spans the years 1979-2021 and was sourced from various government agencies, including the Turkish Statistical Institute [91], the Turkish Ministry of Energy and Natural Resources (MENR) [92], and the World Bank [1,2]. Additionally, the energy consumption data (measured in million tons of oil equivalents, MTOE) was obtained from the MENR. The details of the variables are given in Table 2.

The predictors mentioned above have been commonly utilized in numerous energy forecasting studies, as seen in Table 1. Considering the data collection period from 1979 to 2021, the population grew from 43.19 million to 84.78 million, while the GDP increased from $82 billion to $819.04 billion, indicating a roughly 2 times and 10 times increase, respectively, by 2021. Import and export volumes also saw significant growth, rising from 5.07 and 2.26 to 271.42 and 225.29, respectively, marking approximately 55 times and 100 times increases by 2021. Furthermore, the demand for transportation energy surged nearly fivefold, from 26.37 Mtoe in 1979 to 123.86 Mtoe in 2021. Detailed historical data for these parameters from 1979 to 2021 can be found in Table 3.

The dataset used for predicting of energy is divided into training and test subsets, comprising approximately 75% and 25% of the total observations, respectively. The training set consists of 32 observations, while the test set has 11 samples. Data normalization is performed as the first step in the data preprocessing block. Normalization ensures that all features have a similar scale, preventing any particular feature from dominating the model's learning process. MinMax scale method with the intervals [0, 1] is used to normalize each feature individually.

4. Results and Discussion

4.1. Implementation Setup

A detailed overview of our implementation setup is given in this section. Python, a popular and general-purpose programming language that allows to work quickly and integrate systems more effectively was used. PyCaret, a Python-based open-source machine learning library, was employed for automating machine learning workflows. PyCaret serves as a wrapper around multiple machine learning libraries and frameworks, simplifying the process. Python version 3.11.0 and PyCaret version 3.0.4 were used, which were the latest versions as of October 2022 and July 2023, respectively. All Codes were implemented in Google Colab, a cloud-based platform that provides free access to GPUs and TPUs for running machine learning experiments. All experiments were conducted on a system equipped with an Intel i7 3.40 GHz processor and 8 GB of memory.

4.2. Feature Selection

In machine learning, a correlation matrix is a table that shows how different features in a dataset are related to each other and how they affect the outcome of a model. Figure 2 presents the correlation matrix of the dataset, which includes the target variable as one of the features.

The values within the matrix describe both the intensity and direction of the correlation between pairs of features. Each element represents the correlation between two specific features. In Figure 2, the maximum correlation value is 1, while the minimum is 0.92, observed between the 'import-population' and 'population-GDP' features. A positive correlation between two features implies that as one property's value increases, the other feature's value also tends to increase. It's worth noting that all features exhibit correlations with each other. Additionally, population exhibits the strongest correlation with the target variable 'Energy,' while GDP shows the weakest correlation.

The selection and combination of features is important in machine learning because it can affect the performance and complexity of the model. Different features may have different levels of relevance, redundancy, and noise for a given problem and a given algorithm. By selecting and combining the most appropriate features, the dimensionality of the data, the computational cost and the risk of overfitting can reduce. All combinations of predictor variables (i.e., population, GDP, import and export) are outlined within Table 4. For instance, Model 1 (M1) comprises of two independent variables, i.e. GDP and population, Model 7 (M7): GDP, population, import. The ML performance results with 5-fold cross-validation on the training set, using the created models that include different combinations of features, are given in Table 4. The best three results, sorted by the highest R-square, are provided in the Table 5.

During the model building process, all possible combinations were explored, and finally, the configuration with four inputs, displaying the highest R-squared and lowest error terms, was chosen for application in the next part of the study.

4.3. Performance Evaluation

Firstly we compared the performance of 19 ML individual algorithms using five different metrics with all features (GDP, Population, Import, Export). Table 6 presents the performance results achieved by training 19 ML algorithms using 5-fold cross-validation. The first column in Table 6 lists the base ML algorithms. The subsequent columns, numbered 2nd through 6th, display the best values for various training-phase metrics, including MAE, MSE, RMSE, R2, and MAPE. The ML algorithms results are organized in descending order of R-squared values, from the highest to the lowest.

Table 6. Results of ML algorithms based on train set.

ML Algorithm	MAE	MSE	RMSE	R2	MAPE
Extra Trees Regressor	2296.86	8756864.57	2932.96	0.9788	0.0464
Random Forest Regressor	3186.05	14777499.11	3817.37	0.9684	0.0658
Ridge Regression	3676.12	21641675.00	4466.14	0.9655	0.0736
Linear Regression	3780.00	23739669.80	4668.86	0.9635	0.0825
Lasso Regression	3779.85	23736214.80	4668.54	0.9635	0.0825
Least Angle Regression	3780.00	23739655.90	4668.86	0.9635	0.0825
Lasso Least Angle Regression	3779.85	23736205.30	4668.54	0.9635	0.0825
Orthogonal Matching Pursuit	3780.00	23739655.90	4668.86	0.9635	0.0825
Huber Regressor	3828.58	22823167.86	4595.40	0.9634	0.0785
AdaBoost Regressor	3575.50	15583915.43	3934.88	0.9611	0.0691
Gradient Boosting Regressor	3772.78	16872570.61	4096.12	0.9556	0.0707
Decision Tree Regressor	3768.74	16856622.84	4094.38	0.9554	0.0707
Extreme Gradient Boosting	3768.71	16856282.40	4094.34	0.9554	0.0706
K Neighbors Regressor	3987.62	29417635.00	5274.57	0.9493	0.0848
Elastic Net	7402.15	81150325.20	8721.88	0.8666	0.1248
Bayesian Ridge	22303.79	705872003.20	25684.29	-0.0970	0.4306
Light Gradient Boosting Machine	22303.79	705872041.79	25684.29	-0.0970	0.4306
Dummy Regressor	22303.79	705872041.60	25684.29	-0.0970	0.4306
Passive Aggressive Regressor	40863.75	2361836689.9	48178.33	-3.9041	0.5635

Table 6. Results of ML algorithms based on train set.

ML Algorithm	MAE	MSE	RMSE	R2	MAPE
Extra Trees Regressor	2296.86	8756864.57	2932.96	0.9788	0.0464
Random Forest Regressor	3186.05	14777499.11	3817.37	0.9684	0.0658
Ridge Regression	3676.12	21641675.00	4466.14	0.9655	0.0736
Linear Regression	3780.00	23739669.80	4668.86	0.9635	0.0825
Lasso Regression	3779.85	23736214.80	4668.54	0.9635	0.0825
Least Angle Regression	3780.00	23739655.90	4668.86	0.9635	0.0825
Lasso Least Angle Regression	3779.85	23736205.30	4668.54	0.9635	0.0825
Orthogonal Matching Pursuit	3780.00	23739655.90	4668.86	0.9635	0.0825
Huber Regressor	3828.58	22823167.86	4595.40	0.9634	0.0785
AdaBoost Regressor	3575.50	15583915.43	3934.88	0.9611	0.0691
Gradient Boosting Regressor	3772.78	16872570.61	4096.12	0.9556	0.0707
Decision Tree Regressor	3768.74	16856622.84	4094.38	0.9554	0.0707
Extreme Gradient Boosting	3768.71	16856282.40	4094.34	0.9554	0.0706
K Neighbors Regressor	3987.62	29417635.00	5274.57	0.9493	0.0848
Elastic Net	7402.15	81150325.20	8721.88	0.8666	0.1248
Bayesian Ridge	22303.79	705872003.20	25684.29	-0.0970	0.4306
Light Gradient Boosting Machine	22303.79	705872041.79	25684.29	-0.0970	0.4306
Dummy Regressor	22303.79	705872041.60	25684.29	-0.0970	0.4306
Passive Aggressive Regressor	40863.75	2361836689.9	48178.33	-3.9041	0.5635

Table 7. Results of the Extra Trees Regressor algorithm based on the test set.

ML Algorithm	MAE	MSE	RMSE	R2	MAPE
Extra Trees Regressor	2989.27	17145375.48	4140.6975	0.9811	0.0406

Table 7. Results of the Extra Trees Regressor algorithm based on the test set.

ML Algorithm	MAE	MSE	RMSE	R2	MAPE
Extra Trees Regressor	2989.27	17145375.48	4140.6975	0.9811	0.0406

The hyperparameters were tuned via grid searches because it is a critical step in the machine learning model development process. When ML performance degraded, this step was skipped, and the model was applied to the successive stages without tuning the hyperparameters. After this stage, ensemble methods were suggested and applied for predicting Türkiye’s energy demand. The prediction performance of the ensemble methods (bagging, boosting, blending and stacking) in both training and test sets are shown in Table 8. Compared to the preliminary results of data from the training of 19 ML algorithms shown in Table 6, the mean R-squared values are as follows: 0.9801 with bagging, 0.9809 with boosting, 0.9874 with blending, and 0.9882 with stacking methods. Among these, the stacking ensemble model yielded the highest R-squared value, indicating its superior performance. Additionally, when considering other evaluation metrics such as MAE, MSE, RMSE, and MAPE, the stacking ensemble model consistently outperforms the others, further confirming its superiority in predictive accuracy.

Bagging and Boosting techniques are used to improve the accuracy and robustness of our individual machine learning model. The Extra Tree Regressor (ET) algorithm is trained on different subsets, which are created through a process called bootstrapping, of the training data by the Bagging ensemble method. In boosting, the focus is on correcting the errors made by previous models. The base model ET was trained until a certain level of accuracy is achieved.

In our blending approach, we harnessed the predictive power of three distinct machine learning algorithms: Extra Trees Regressor (ET), Random Forest Regressor (RF), and Ridge Regression (Ridge). To execute blending, each of these algorithms were initially trained separately on a portion of our training dataset, generating individual predictions for our target variable. Subsequently, we combined the predictions from ET, RF, and Ridge using a straightforward averaging technique. By averaging these predictions, we effectively created an ensemble prediction that capitalizes on the strengths of each algorithm.

Unlike traditional ensemble methods like bagging and boosting, stacking takes a more sophisticated approach by using the predictions of base models as input features to train a higher-level model that makes the final predictions. The 14 ML algorithms with R-squared values higher than 0.90 in Table 6 were separately combined to select a set of diverse base models. After trying many combined base models and conducting trial and error experiments, we leveraged the capabilities of three diverse base machine learning algorithms: ET, RF and Ridge. To implement stacking, each of these base models were initially trained separately on a portion of our training dataset, obtaining individual predictions for our target variable. Next, a new dataset was created where each data point consisted of these base model predictions. This dataset served as the input for the meta-learner. Based on our empirical investigation, it was determined that employing a linear regression algorithm as the meta learner for the second level of the stacking regressor was strong, as it consistently demonstrated superior performance in terms of R-squared compared to alternative machine learning algorithms. Ridge Regression, as the meta-learner on the second level of our stacking ensemble, was selected and trained to learn how to best combine the predictions from ET, RF, and Ridge. To prevent overfitting during the training phase, we employed 5-fold cross-validation.

Stacking's flexibility in utilizing both weak and strong learners make it a powerful technique for enhancing predictive performance in various machine learning tasks. In practice, researchers often employ a mixture of weak and strong learners to construct a versatile ensemble that performs effectively across diverse datasets and problem domains. The choice of whether the base models are weak or strong is flexible and depends on the problem and the effectiveness of the ensemble. In this study, we experimented with different combinations of both weak and strong base models to create a diverse ensemble. It was observed that the combinations composed of strong base learners consistently delivered superior results in forecasting Türkiye's energy demand.

In order to evaluate the efficacy of the developed ensemble methods, the prediction performance on the test set is presented in Table 8. Notably, the Stacking Model achieved a remarkable R-squared value of 0.9826. When compared with the findings in Table 7 and other ensemble models, the Stacking Model's metrics consistently reveal a significantly enhancement. These observations lead us to assert, in accordance with our scientific paper, that our proposed stacking ensemble model does not exhibit signs of either overfitting or underfitting. The utilization of features (population, GDP, import, export) enhances the model's accessibility and interpretability, facilitating its utility for generating accurate and reliable forecasts of Türkiye's energy demand, as detailed in the paper.

Table 9 provides detailed descriptions of the model's predicted outcomes in the 'Prediction' column, alongside the corresponding ground truth values for 'Energy.' It presents the prediction performance of the stacking ensemble model for each of the five folds, utilizing all dataset features and evaluating it using five different metrics.

5. Conclusions

This paper has presented a comprehensive methodology for applying ensemble techniques and machine learning algorithms to the crucial task of forecasting Türkiye's energy demand. The primary objectives of this methodology were threefold: Firstly, to enhance the accuracy of energy demand predictions in Türkiye. Secondly, to provide authorities and institutions with an interpretable model that facilitates informed decision-making and policy development. Lastly, this study aligns closely with the United Nations' Sustainable Development Goals (SDGs), contributing to the broader aims of sustainable development by addressing global challenges.

Accordingly, the following key findings can be derived based on the current research.

• The GDP, population, import, export, energy data taken between 1979 and 2021 were used and it is observed that there has a strong correlation among them.
• Five statistical metrics are discussed to evaluate the performance of the algorithms in the forecast.
• A total of 19 machine learning algorithms were constructed and analyzed to select models for diverse ensemble combinations.
• Considering all metrics collectively, the stacking ensemble model utilizing Ridge Regressor as a meta-learner outperforms single ML algorithms as well as other bagging, boosting, and blending models.
• The results of predicted values demonstrate that stacking ensemble model has presented very satisfied results when comparing the truth energy-demand outputs.
• These ensemble models can readily be adapted and recommended for future energy demand forecasts in other countries. Notably, the stacking ensemble model demonstrates statistically superior results compared to other models, making it a more suitable choice for accurate forecasting.

It is anticipated that the outcomes of this study will make a significant contribution to the field of energy forecasting, laying the groundwork for Türkiye's sustainable energy future. Furthermore, this research represents a meaningful step toward a more equitable, prosperous, and sustainable world for all. As future work, further improvements can be explored through the use of different hybrid techniques for optimizing hyperparameter tuning, feature selection, and more.