Enhancing the Predictability of Wintertime Energy Demand in The Netherlands Using Ensemble Model Prophet-LSTM

Abstract

Energy demand forecasting is crucial for maintaining stable and affordable energy supplies, especially for vulnerable populations most affected by shortages and high costs. In the Netherlands, transmission system operator TenneT has raised concerns about potential electricity shortages by 2030. Rising energy prices and the impact of climate change on the energy demand further complicate today’s energy market. Policymakers lack clear insights into demand patterns, which complicates the optimization of energy use and the protection of at-risk communities. Accurate and timely forecasts are essential for addressing these issues and supporting sustainable energy management. This research focuses on enhancing the accuracy and lead time of wintertime energy demand forecasts in the Netherlands using advanced machine learning. The ensemble model Prophet-LSTM is trained on hourly load consumption data combined with climate change-related and energy price predictors. The results demonstrate significant improvements over baseline models, achieving a Pearson correlation coefficient of $r=0.93$ compared to $r=0.50$ in prior studies, as well as accurate forecasts up to 180 days ahead, compared to 2 months. Incorporating climate change-related predictors is challenging due to multicollinearity, highlighting the importance of careful predictor selection. Including energy price predictors yielded modest yet hopeful results, suggesting their ability to optimize energy demand forecasting.

1. Introduction

Wintertime energy demand forecasting (EDF) plays a crucial role in ensuring a reliable energy supply and managing available resources. It also helps mitigate the impact of extreme weather events caused by climate change on vulnerable populations, particularly in the Netherlands. In 2023, Dutch transmission system operator TenneT expressed concerns about potential electricity shortages in the Netherlands by 2030 [1]. From a global perspective, the current security of the electricity supply in the Netherlands is extremely high at $99.99963\%$. However, the number of controllable coal and gas power stations in the Netherlands is set to decrease, while the country becomes increasingly dependent on sustainable electricity, which is highly dependent on weather conditions. This makes it uncertain whether sufficient electricity can be generated at any time to meet demands [2,3]. Therefore, accurate and timely EDF is crucial.

Important influences on the Dutch energy market are energy prices and climate change. Energy prices have been quickly increasing in the Netherlands, resulting in very high costs for Dutch households compared to other European Union (EU) countries [4]. The average gas price in the Netherlands is EUR $0.182$ per kWh, significantly higher than the EU average of EUR $0.110$ per kWh [5]. In 2022, high energy prices caused inflation to reach $5.2\%$, decreasing purchasing power by an average of $2.7\%$. The Dutch government has taken measures to help low- and middle-income households facing rising costs [6], but many continued to struggle financially [4]. In addition, gas consumption decreased by $15\%$ due to significantly higher gas prices in 2022, with people saving energy by turning down heating, taking shorter showers, or better insulating their homes [7]. Climate-wise, it is notable that climate change is causing shifts in the weather in the Netherlands, leading to more rain and heavier storms [8]. The wintertime energy demand fluctuates yearly and is closely tied to the variability of weather conditions [9]. Thus, changing weather patterns can lead to shifts in the energy demand, highlighting the importance of capturing these changes to enable people to meet their energy demands. Moreover, the share of renewable energy of the total energy consumption in the Netherlands was $15\%$ in 2022 [10]. The Climate Agreement aims for $70\%$ renewable electricity by 2030, and the EU aims for at least $32\%$ sustainable energy by 2030 [11]. However, climate change is causing less certainty in weather patterns in Europe, which has a negative long-term effect on the production of renewable energy [12]. For people to still meet their energy demands, less available renewable energy requires them to use more nonrenewable energy.

Ultimately, enhancing the accuracy and lead time of EDF in the Netherlands is crucial to protect the energy supply against evolving challenges such as increasing energy prices and the impact of climate change. In particular, TenneTs’ concerns highlight the importance of accurately forecasting the energy demand at any time to mitigate the effects of potential shortages and thereby protect vulnerable populations.

This study addresses the research gap in the unique dynamics of the Dutch energy market by introducing an advanced forecasting approach tailored to this specific context. While EDF has been extensively studied, examining the Dutch context provides novel insights. High energy prices negatively impact the energy demand [7], while the Netherlands is shifting toward renewable energy due to national climate goals and the phasing out of coal and gas power stations [1,11]. However, both nonrenewable and renewable energy consumption are reliant on weather conditions, which are impacted by climate change [8]. These factors highlight the complexity of the Dutch energy market and the need for an appropriate approach to EDF in the Netherlands.

The primary contribution of this study lies in extending current EDF research by applying advanced machine learning (AML) techniques specifically to Dutch energy data, building upon the research of [9]. In emphasizing the critical interaction between environmental, economic, and energy dynamics, this research aimed to improve forecasting accuracy and lead time—the length of time into the future for which predictions are made. By fitting this approach to Dutch energy data and addressing climate change impacts and energy price fluctuations, this study improves academic understanding and offers practical insights for energy planning strategies. This promises a more adaptive and resilient forecasting framework for the unique dynamics of the Netherlands. A novel approach is introduced that combines predictors closely related to the Dutch energy market to demonstrate its effectiveness in improving Dutch EDF performance.

2. Related Work

EDF is an extensively researched topic with research directions varying in both predictor and model selection. This section discusses the relevant predictors and models for EDF.

2.1. Baseline Approach

Challenges in EDF are two-fold: (1) rising energy prices are influencing the energy demand [7] and (2) climate change is shifting weather patterns, which influences energy demand patterns [8].

One of the most dominant approaches is proposed in [9], using a multiple linear regression (MLR) model to predict winter seasonal energy consumption two months in advance in the United Kingdom (UK), using climate change-related predictors. Leave-one-out cross-validation was used to evaluate the model using the Pearson correlation coefficient (PCC), yielding a significant correlation of $r=0.50$ between the predicted and observed time series. The study explored the relationship between surface climate conditions and daily UK wintertime gas and electricity consumption, revealing that the prediction skill of energy consumption arises from the predictability of the North Atlantic oscillation (NAO) variability and the joint estimation of surface climate condition variables temperature, dew point depression, and wind speed. One main limitation of this work was the lack of validation of the model’s generalizability.

However, recent developments regarding rising energy prices and climate change, which affect the energy demand, require a more advanced approach to EDF.

2.2. Impact of Climate Change on Energy Demand

Currently, the EDF challenge is addressed by focusing mostly on climate change. Based on recent papers, it is notable that climate change-related features are frequently used as predictors, emphasizing the importance of these predictors in today’s energy market.

In 2022, the residential gas consumption in the Netherlands decreased by $22\%$ compared to 2021, depicting the lowest energy usage since 1990. Warmer weather conditions contributed to this reduction, demonstrating how changing weather conditions influence the energy demand. In the agricultural sector, a $31\%$ decrease in gas consumption was observed, mostly driven by high gas prices, but also influenced by warmer weather [13]. These trends emphasize the relevance of climate change-related predictors in EDF, such as temperature, which play an important role in understanding fluctuations in the energy demand.

Firstly, relevant climate change-related predictors are the sea ice concentration (SIC), stratospheric circulation (SC), sea surface temperature (SST), and the NAO variability index. The NAO is the dominant mode of winter variability over the North Atlantic and strongly influences surface weather conditions of the winter period such as average temperature, wind speed, and storminess, over much of the European continent and the United States [9,14,15]. Autumn conditions in the Arctic SIC, SC (especially the geopotential at 70 pHa), and SST influence the NAO variability [9,16]. SIC measures the amount of sea ice in a given area, typically expressed as a percentage [17]. SC is the air movement in the stratosphere [18]. SST indicates the surface temperatures in a broad swath of the tropical North Atlantic Ocean [19]. These are all important predictors as they influence NAO variability, which in turn affects energy consumption.

Secondly, relevant climate change-related predictors are weather-related such as surface air temperature, dew point depression, and wind speed, all showing significant correlations with gas consumption [9]. In the UK, winter means of temperature, the NAO, a UK-centered north–south pressure difference, and an additional predictor based on the frequency of high-demand weather types over the winter period are used to predict the gas demand and the number of extreme gas demand days over the winter period [15], demonstrating that seasonal weather forecasts can skillfully predict the weather-driven component of the winter gas demand and the number of extreme gas demand days in the UK. Similar research has also been performed outside of the UK, namely in Hanoi where temperature, absolute humidity, and wind speed are used to forecast the short-term electricity load [20].

2.3. Impact of Energy Price and Renewable Energy on Energy Demand

Furthermore, considering energy prices and renewable energy in EDF is also crucial. The relationship between the renewable and nonrenewable energy demand is complex. Intuitively, an increase in renewable energy demand results in a decrease in nonrenewable energy demand, given that only so much is needed. However, [21] suggests that it is unlikely that renewable and nonrenewable energy are perfect substitutes, suggesting some level of complementarity between these two forms of energy. The relationship between the nonrenewable and renewable energy demand was further investigated in [22], studying the elasticity of renewable and nonrenewable energy in countries that are part of the Organisation for Economic Co-operation and Development (OECD), which includes the Netherlands. The study showed that nonrenewable and renewable energy are substitutes in some industries and complements in others.

Additionally, the impact of nonrenewable energy prices on US renewable energy consumption was studied, showing a positive relationship [21]. This indicates that increases in nonrenewable energy prices cause increases in the renewable energy demand, resulting in decreases in the nonrenewable energy demand. This finding is supported by the authors of [23], who investigated the response of the electricity demand to price changes in Portugal, showing that increases in price lead to decreases in electricity consumption, with lower-income groups being most affected.

As briefly mentioned in Section 2.2, a $31\%$ decrease in gas consumption was observed in the Dutch agricultural sector, primarily driven by high gas prices [13]. This shows how significantly price can influence the energy demand, illustrating the importance of incorporating energy prices into EDF.

2.4. Energy Demand Forecasting Models

The consideration of modeling approaches is necessary to tackle the challenging problem of forecasting energy demand accurately.

2.4.1. More Traditional Models

Traditional statistical models are foundational in EDF. For instance, in [15], an LR model was utilized to predict Britain’s gas demand and the number of extreme gas demand days over the winter period, showing skillful predictions. Furthermore, an LR model and a long short-term memory (LSTM) model—a type of recurrent neural network that is effective for tasks such as time series forecasting due to its ability to learn and model long-term dependencies in sequential data [24]—were compared for predicting the daily and monthly El Niño–Southern Oscillation [25], showing that monthly performances from LSTM and LR are similar, whereas daily results show that the LSTM has some advantage over LR regarding correlation coefficient. This suggests that LSTM is superior at capturing nonlinear relationships compared to LR, although at a computational cost.

2.4.2. AML Models

State-of-the-art (SOTA) EDF research focuses on utilizing AML models, often comparing them with more traditional models. For example, a combined model based on LSTM and LightGBM was used to forecast power load in [26]. This model decomposes historical power load data with Empirical Mode Decomposition and combines this with historical weather data to establish LSTM and LightGBM prediction models, showing improved forecasting accuracy and application prospects for power load forecasting compared to traditional methods and standard LSTM and LightGBM methods. Moreover, the ensemble model NeuralProphet-LightGBM was used in [27] to forecast power load. By first decoupling the original sequence using NeuralProphet and then using LightGBM for further feature extraction, their model shows significantly better performance in comparison with comparative models.

Overall, the finding that LSTM models are better at capturing nonlinear relationships compared to LR models [25], in combination with the success of ensemble models in EDF [26,27], motivates the exploration of similar approaches for tackling the challenging task of EDF, combined with the use of the important predictors of climate change and energy price.

3. Methodology

This research concerns a multivariate time series forecasting task. This necessitates five key steps: data collection, exploratory data analysis (EDA), data cleaning and preparation, and model implementation and evaluation.

3.1. Data Collection

Multiple open source data sets were integrated to create a data set comprising predictors and a target variable, intended for model training and evaluation. The data set included energy consumption as the target variable, together with predictors related to climate change and energy prices.

3.1.1. Target Variable

Energy consumption in the Netherlands served as the target variable. For this purpose, load consumption data from the European Network of Transmission System Operators for Electricity (ENTSO-E) were utilized. ENTSO-E defines load as “an end-use device or customer that receives power from the electric system” [28]. The data set, ranging from 2009 to 2019, included hourly load information in megawatts from the Netherlands.

3.1.2. Predictors

The predictors in this study were subdivided into two categories: climate change-related and energy prices.

Climate change-related: Daily temperature, precipitation, and wind speed were used—previously utilized in studies such as [9,15,20]—obtained from the Royal Dutch Meteorological Institute [29]. The data, ranging from 1901 to 2024, originated from a weather station located in De Bilt, the Netherlands, which is representative of the average Dutch weather due to its central location [30].

Moreover, the NAO variability index was included, similar to [9,15]. Monthly data from 1950 to 2023 were sourced from the National Oceanic and Atmospheric Administration (NOAA) [31].

Furthermore, predictors that influence the NAO index—stratospheric circulation (geopotential height at 70 pHa), SIC and SST—were included [9,16]. The data were accessed from the Copernicus Climate Data Store and the NOAA on daily, every-three-day, and monthly levels [17,18,19], ranging from 2009 to 2019. The geopotential height and SIC data have high dimensions due to the numerous longitude and latitude values. Principal Component Analysis was applied to reduce the dimensionality to a single daily value, as previously discussed in [16].

In addition, the yearly ratio of renewable energy consumption relative to total energy consumption was incorporated as a predictor, sourced from the Central Bureau of Statistics (CBS) [32], ranging from 1990 to 2023.

Energy price: Finally, Dutch energy prices were included. Half-yearly average energy prices for consumers in the Netherlands were used, obtained from CBS and ranging from 2009 to 2023 [33].

3.2. EDA

EDA was conducted on all data sets individually and collectively, revealing daily, monthly, and yearly seasonality patterns in the load consumption data, as illustrated in Figure 1 and Figure 2. Daily patterns show higher load consumption on weekdays compared to weekends, with peaks typically at 7:00 and 18:00, and lower load consumption during nighttime and afternoon hours, as depicted in Figure 3. Figure 4 shows that load consumption peaks in the winter months (December, January, and February) and decreases in the summer.

3.2.1. Seasonal Decomposition

Seasonal decomposition was used to improve understanding of the underlying patterns in the load consumption data by decomposing it into its seasonal, trend-cycle, and remainder components [34,35], as displayed in Figure 5. This revealed an upward trend over the years and a seasonal pattern showing that consumption increases in the winter and decreases in the summer.

3.2.2. Autocorrelation Analysis

The seasonality in the data was further investigated using the autocorrelation function (ACF) to measure the linear relationship between lagged values of a time series [36]. The ACF was used to identify the appropriate lags of the load consumption target variable to use as predictors in the model.

Figure 6, Figure 7 and Figure 8 depict the ACF of hourly, daily, and monthly load consumption data, respectively. The plots reveal significant peaks at 24 h, 7 days, and 12 months. In Figure 6, spikes at both 24 h and 48 h indicate daily seasonality. Figure 7 shows spikes every 7 days, indicating weekly seasonality. Figure 8 displays a spike at 12 months, indicating yearly seasonality. The daily, weekly, and yearly seasonality in the data advocates for the use of lagged features of 24 h, 7 days, and 12 months, corresponding to these periods.

3.3. Data Cleaning and Preparation

Preprocessing steps are applied to each individual data set, after which they are combined into a unified data set that includes all predictors and the target variable.

3.3.1. Handling Outliers

Based on the EDA findings, the data exhibit outliers. This is addressed using the interquartile range ($IQR$) method as a fitting outlier detection method [37,38]. The $IQR$ is robust against outliers because it relies on the middle $50\%$ of the data distribution. The IQR method considers a data point an outlier if it is more than $1.5$ × $IQR$ below the first quartile $q1$ or above the third quartile $q3$, where $IQR=q3-q1$. Outliers are only cleared from the train folds to prevent overfitting.

3.3.2. Shifting Data

The aim is to forecast the load consumption (at least) two months in advance. Hence, at the point of forecasting, only the reality of two months ago is known. The data are shifted for this purpose, such that day t has the predictors of day $t-60$ and the lags of load consumption of 24 h, 7 days, and 12 months of day $t-60$. The 12-month lag results in the loss of the first year of the data set, 2009, since this requires a year of data available before a day t. Hence, with a data set that starts 1 January 2009, the first moment possible to take a lag of 60 days and an additional 12 months is 2 March 2010.

3.3.3. Handling Missing Values

Outlier removal resulted in some missing values in the data set. The missing values were imputed using Multivariate Imputation by Chained Equations (MICE), similar to in [39]. MICE involves a series of linear regressions run in sequence and repeated until the estimated values for missing values have sufficiently converged. To prevent data leakage, which is the case when data outside of the training data set is used to train the model, MICE is solely fitted to the training fold and is exclusively used to transform the test fold.

3.3.4. Predictor Selection

The model has 10 available predictors. The forecasting model was tested both including and excluding predictors. The impact of energy price and climate change-related predictors was tested by either including or excluding them from the model. Furthermore, the seasonal lags (as mentioned in Section 3.2.2) of the target variable load consumption were included in each of the predictor selections.

3.3.5. MinMaxScaling

MinMaxScaling was used to scale all of the features of the LSTM model, similar to in [40]. An advantage of this scaler is that it provides a uniform $[0,1]$ scale across all features [41]. The scaling was carried out using the following formula:

$$x_{scaled}={\displaystyle \frac{x-x_{min}}{x_{max}-x_{min}}}$$

where x is the original value of a feature, and $x_{min}$ and $x_{max}$ are the minimum value and the maximum value, respectively, of that feature in the data set. To prevent data leakage, the MinMaxScaler was solely fitted to the training data, and the fitted scaler was used to transform both the training and testing data sets.

3.4. Model Implementation: Baseline Model

For a baseline model, the traditional statistical forecasting model MLR was used. MLR has been shown to perform competitively with sophisticated methods [9,25] and offers several advantages, including its ease of implementation and interoperability. MLR allows for the use of multiple predictors, which makes the model suitable for capturing complex relationships between the predictors and target variables, aligning well with the overall goal of this research.

3.5. Model Implementation: AML Model

3.5.1. Ensemble Model

The ensemble model of Prophet and LSTM was used as the AML model, using Prophet for capturing seasonal patterns and holiday effects in the time series data and using LSTM to model the more complex and nonlinear relationships in the data, which might not be captured yet by Prophet. The use of an LSTM model for EDF was suggested by [9], supported by the performances of LSTM models in other studies [25,26,42]. The use of ensemble models is based on [26,27,42].

3.5.2. Prophet

Prophet is an additive time series forecasting tool that was developed by Facebook [43]. It is a nonlinear regression model and is widely recognized for its effectiveness in estimating seasonal fluctuations, performing well for time series with strong seasonality and several seasons of historical data [44]. Prophet automatically fits weekly and yearly seasonalities for time series more than two cycles long [44]. Considering that the data also show daily seasonality, automatic detection for daily seasonalities was enabled. Furthermore, predictors were added based on the predictor selection. Prophet was trained on the entire training data set, and this model was used for predictions on the entire training data set and test set. Prophet’s predictions $y_{hat}$ and its trends for daily, weekly, and yearly seasonal components were combined with the selected predictors as input for the LSTM model [27].

3.5.3. LSTM

LSTM, a type of recurrent neural network, was first proposed in [24]. It is designed to capture and utilize long-term dependencies in sequential data, making it effective for tasks such as time series forecasting. Two LSTM setups were deployed: one including dropout layers (Model 1) and one excluding dropout layers (Model 2), visualized in Figure 9 and Figure 10.

Model 1 comprises an input layer, followed by two LSTM layers, and concludes with a dense layer. The first LSTM layer contains 30 nodes, and the second has 90 nodes. This architecture is similar to that introduced in [25]. However, the referenced study used 50 nodes in the first LSTM layer, and 150 nodes in the second. Due to the smaller data set in this research compared to in the study, the number of nodes was reduced while maintaining the same node ratio. Model 2 comprises two LSTM layers with dropout layers in between, and concludes with a dense layer. The first LSTM layer contains 30 nodes, and the second has 90 nodes. This architecture is similar to that introduced in [42], the only difference being that the study used 50 nodes in the first LSTM layer, and 100 nodes in the second. Again, node numbers were reduced while maintaining the node ratio. Both models were trained using 50 epochs, determined by the number of epochs needed for convergence, and a batch size of 128. Early stopping was implemented to prevent overfitting, allowing the model to stop training if the validation loss does not decrease by $0.0001$ for 3 epochs, starting from epoch 10 [45].

3.6. Model Evaluation

3.6.1. Rolling Origin Cross-Validation

The model was trained and evaluated using rolling origin cross-validation (ROCV). In ROCV, the model is estimated on the training set, and forecasts are created for the test set. For every iteration, another year is added to the end of the training set, and the test set is advanced by the same period. The process stops when there are no more data left [46]. The minimum length of the training set was set at two years. To illustrate, the observations of 2010–2011 were used to train the model, which was then used to predict the load consumption in 2012. Next, the observations of 2010–2012 were used to train the model, which was used to predict for 2013, etcetera. This was repeated every year from 2010 to 2019 to ensure that the model was evaluated on 8 training–test folds.

ROCV was chosen over other methods, such as k-fold cross-validation (KCV), because of its suitability for time series data, where the order of observations matters. Where KCV would randomize the order of the observations, ROCV keeps their temporal order. This ensures that the model is only trained on past observations and predictions are only made on future observations.

To ensure the model’s generalizability, the test fold was subdivided into a wintertime test set, which includes the months of December, January, and February. Similar performance in the winter months and the general test set should give insights into the model’s ability to generalize and robustness.

3.6.2. PCC

The PCC measures the strength of the relationship between two variables, returning a value (r) between $-1$ and 1, where 1 indicates a strong positive relationship; $-1$, a strong negative relationship; and 0, no relationship. It requires the following assumptions to be met: (1) the two variables must be measured at the interval or ratio scale; (2) there is a linear relationship between the two variables; (3) there should be no significant outliers; and (4) the data should be approximately normally distributed [47]. PCC was used to measure the strength of the relationship between the target variable and forecasted value and was appropriate for this work as the load consumption was measured at an interval scale, there was a linear relationship between the forecasted and actual value of the load consumption, outliers were removed, and the data were approximately normal distributed. Furthermore, the relevance of using the PCC is being able to compare this work to reference works.

3.6.3. RMSE

The RMSE is the distance of the prediction from the actual value and is calculated as $RMSE=\sqrt{mean\left(e_t^2\right)}$, where $e_t=y_t-{\widehat{y}}_t$. A disadvantage of the RMSE is its sensitivity to outliers due to its nature of being scale-dependent. Given that this metric takes the square of errors, the RMSE value has the potential to become enormous for extreme errors, as then the forecasted value and actual value are no longer on the same scale [48]. Especially considering that the average RMSE over the ROCV folds was calculated in this work, this has the potential to return extreme values if just a single fold shows bad performance, even though the rest might not.

3.6.4. MAPE

The MAPE is the accuracy of the prediction relative to the actual value and is calculated as $MAPE={\displaystyle \frac{1}{n}}{\sum }_{i=1}^n\left|{\displaystyle \frac{Y_i-{\widehat{Y}}_i}{Y_i}}\right|$. It gives the percentage error and is commonly used in forecasting to compare performances. A limitation of the MAPE is that it can be infinite or undefined if $Y_i=0$ and that it can have extreme values if any $Y_i$ is close to zero [48]. In this research, the values of $Y_i$ were significantly larger than zero, so this disadvantage was neglectable.

4. Results

This section presents the results of training and evaluating the baseline MLR model and the Prophet-LSTM model with different predictor selections. The impacts of including energy price and climate change-related predictors and using a longer forecasting horizon are also investigated.

4.1. Performance of Baseline MLR

To investigate how the wintertime performance differs between the baseline model and the AML model (SRQ1), first, a baseline MLR model was trained, including both including and excluding predictors. The model was tested against multiple test folds consisting of wintertime months.

The results are displayed in

Table 1. Performance of baseline MLR model on wintertime months across different folds in RMSE, MAPE, and PCC.

Pred.	Fold	RMSE	MAPE	PCC
excl.	1	1910.53	11.45	0.71
	2	2207.39	11.85	0.72
	3	1856.90	10.55	0.72
	4	1887.11	9.47	0.71
	5	1743.45	9.29	0.70
	6	1826.48	9.51	0.69
	7	1959.85	9.97	0.63
	8	1878.80	9.29	0.62
	Avg.	1908.81	10.17	0.69
incl.	1	2199.54	12.80	0.71
	2	1958.13	12.62	0.66
	3	1896.96	12.64	0.73
	4	1967.44	10.42	0.73
	5	1753.80	10.04	0.70
	6	1693.44	10.26	0.70
	7	1757.51	9.43	0.65
	8	1777.61	9.11	0.65
	Avg.	1874.75	11.11	0.70

. The results show similar performances for the MLR models, where the inclusion of the predictors ($RMS{E}_{avg}=1874.75$, $MAP{E}_{avg}=11.11$, $PC{C}_{avg}=0.70$) has a slightly improved RMSE, yet a slightly worse MAPE and PCC compared to the exclusion of predictors ($RMS{E}_{avg}=1908.81$, $MAP{E}_{avg}=10.17$, $PC{C}_{avg}=0.69$). Compared to the PCC baseline of $r=0.50$ as stated in the RQ, the MLR shows promising improvement.

4.2. Wintertime Performance of Prophet-LSTM

To further continue the investigation of SRQ1, the Prophet-LSTM model was trained using ROCV and evaluated in the wintertime months.

Table 2. Performance of two Prophet-LSTM model setups with different predictor selections on wintertime months across different folds in RMSE, MAPE, and PCC.

Pred.	Fold	Model 1			Model 2
		RMSE	MAPE	PCC	RMSE	MAPE	PCC
excl.	1	824.56	4.91	0.95	1158.00	7.70	0.92
	2	1291.12	7.47	0.92	1340.38	8.05	0.89
	3	847.10	4.50	0.95	920.69	4.99	0.93
	4	853.98	4.64	0.95	989.01	5.68	0.91
	5	956.52	5.50	0.92	1045.81	6.05	0.89
	6	838.51	4.97	0.94	930.09	4.99	0.92
	7	938.36	4.96	0.91	1113.21	5.89	0.90
	8	922.93	4.94	0.91	990.38	5.24	0.90
	Avg.	934.13	5.23	0.93	1060.95	6.07	0.91
incl.	1	1054.20	7.13	0.95	1217.58	7.45	0.92
	2	2211.39	13.59	0.90	2026.11	11.99	0.85
	3	1269.56	7.23	0.94	935.21	4.90	0.94
	4	1501.63	9.28	0.93	1550.03	9.20	0.88
	5	1109.64	5.89	0.90	1525.65	8.16	0.87
	6	948.57	5.08	0.92	1138.56	5.89	0.91
	7	971.95	5.30	0.90	1413.26	7.59	0.90
	8	922.14	5.05	0.91	1125.58	5.87	0.90
	Avg.	1248.63	7.32	0.92	1366.50	7.63	0.90

shows the performance of the Prophet-LSTM models during the wintertime period. In comparing the predictor selections for Model 1, it is notable that the exclusion of predictors ($RMS{E}_{avg}=934.13$, $MAP{E}_{avg}=5.23$, $PC{C}_{avg}=0.93$) outperforms the inclusion of predictors ($RMS{E}_{avg}=1248.63$, $MAP{E}_{avg}=7.32$, $PC{C}_{avg}=0.92$), especially showing improvements in the RMSE and MAPE. Excluding predictors gives more stable results across the different folds than including predictors.

In comparing the predictor selections for Model 2, it is notable that the exclusion of predictors ($RMS{E}_{avg}=1060.95$, $MAP{E}_{avg}=6.07$, $PC{C}_{avg}=0.91$) outperforms the inclusion of predictors ($RMS{E}_{avg}=1366.50$, $MAP{E}_{avg}=7.63$, $PC{C}_{avg}=0.90$), again especially showing large improvements in terms of the RMSE and MAPE. Again, performance across the different folds is more stable when excluding predictors.

When comparing the performances of Model 1 and Model 2, the results show that Model 1 consistently outperforms Model 2 across all average metrics. Moreover, Model 1 shows more stable results across the different folds when excluding predictors in comparison to Model 2. In contrast, both models show variability in performance metrics when including predictors. Overall, Model 1 outperforms Model 2 in terms of average metrics and demonstrates better stability across the different folds, especially when predictors are excluded.

The results show that for the models excluding predictors, Model 1 ($RMS{E}_{avg}=934.13$, $MAP{E}_{avg}=5.23$, $PC{C}_{avg}=0.93$) outperforms Model 2 ($RMS{E}_{avg}=1060.95$, $MAP{E}_{avg}=6.07$, $PC{C}_{avg}=0.91$) in terms of all average metrics. For the models including predictors, the performances of Model 1 ($RMS{E}_{avg}=1441.66$, $MAP{E}_{avg}=8.56$, $PC{C}_{avg}=0.91$) and Model 2 ($RMS{E}_{avg}=1452.50$, $MAP{E}_{avg}=8.06$, $PC{C}_{avg}=0.90$) are very similar. Model 1 has a slightly improved RMSE and PCC, yet a lower MAPE than Model 2.

MLR vs. Wintertime Prophet-LSTM

To answer SRQ1, a comparison is made between the obtained results from the MLR and Prophet-LSTM models on the wintertime test set. Observable from the results is that the Prophet-LSTM model significantly outperforms the MLR model with both model setups and with all predictor selections, showing improved RMSEs, MAPEs, and PCCs compared to the best-performing MLR model ($RMS{E}_{avg}=1874.75$, $MAP{E}_{avg}=11.11$, $PC{C}_{avg}=0.70$). The Prophet-LSTM model performs best when excluding predictors, whereas the MLR model performs slightly better when including them. Again, compared to the PCC baseline of $r=0.50$, as stated in the research question, both the MLR and especially the Prophet-LSTM model show significant improvement in PCC.

4.3. General Performance of Prophet-LSTM

To further assess the predictive power and generalization ability of the Prophet-LSTM models, their performances were examined on the entire test set.

The results in

Table 3. Performance of two Prophet-LSTM model setups with different predictor selections on entire test fold across different folds in RMSE, MAPE, and PCC.

Pred.	Fold	Model 1			Model 2
		RMSE	MAPE	PCC	RMSE	MAPE	PCC
excl.	1	701.18	4.41	0.96	1604.70	12.01	0.87
	2	1077.49	6.83	0.95	1452.22	9.94	0.86
	3	674.92	3.97	0.96	1094.25	7.12	0.93
	4	704.01	4.31	0.96	1450.01	9.64	0.86
	5	1215.83	7.78	0.86	1584.58	10.39	0.76
	6	962.01	5.81	0.89	1105.99	6.84	0.85
	7	957.63	5.54	0.88	1063.14	6.25	0.84
	8	988.72	5.75	0.91	1021.10	5.81	0.88
	Avg.	910.22	5.55	0.92	1297.00	8.50	0.86
incl.	1	1478.74	10.98	0.93	1168.60	8.15	0.89
	2	1878.18	12.37	0.93	1596.42	9.91	0.90
	3	1258.06	7.82	0.95	770.80	4.55	0.95
	4	2044.24	13.94	0.90	2041.10	14.02	0.88
	5	1266.63	7.97	0.82	1722.65	11.26	0.71
	6	1045.57	6.32	0.86	1126.89	6.32	0.85
	7	1192.92	7.59	0.84	1342.90	7.52	0.85
	8	1172.86	7.21	0.86	1262.91	6.88	0.87
	Avg.	1417.15	9.27	0.89	1379.03	8.58	0.86

show the performances of the Prophet-LSTM models for both predictor selections on the entire test set. Firstly, Model 1 performs better when excluding predictors ($RMS{E}_{avg}=910.22$, $MAP{E}_{avg}=5.55$, $PC{C}_{avg}=0.92$) rather than including them ($RMS{E}_{avg}=1417.15$, $MAP{E}_{avg}=9.27$, $PC{C}_{avg}=0.89$). A similar pattern is seen for Model 2, showing improved performance when excluding predictors ($RMS{E}_{avg}=1297.00$, $MAP{E}_{avg}=8.50$, $PC{C}_{avg}=0.86$) contrary to including them ($RMS{E}_{avg}=1379.03$, $MAP{E}_{avg}=8.58$, $PC{C}_{avg}=0.86$). Moreover, the results show that Model 1 outperforms Model 2 when excluding predictors. However, the evaluation of the entire test set shows that Model 2 performs slightly better than Model 1 when including predictors.

Best Model Selection

For computational reasons, the decision was made to continue the rest of the evaluation with a single model, which was the one considered the best Prophet-LSTM model. This model was chosen based on its performance and its ability to generalize.

Across all the Prophet-LSTM models, Model 1 excluding predictors outperformed all the other models on both the general ($RMS{E}_{avg}=910.22$, $MAP{E}_{avg}=5.55$, $PC{C}_{avg}=0.92$) and the wintertime test sets ($RMS{E}_{avg}=934.13$, $MAP{E}_{avg}=5.23$, $PC{C}_{avg}=0.93$). Additionally, it had the best ability to generalize, showing very stable performances across the different folds as well as very similar performances on the wintertime and general test sets.

4.4. Impact of Predictors

4.4.1. Climate Change-Related Predictors

The impact of climate change-related predictors (SRQ2) was investigated by integrating them into the model, after which the performance was compared to the model excluding them.

Figure 11 shows the performance of including the climate change-related predictors, on both the entire test set and wintertime subset. The detailed performances are noted in

Table 4. Performance of Prophet-LSTM model trained with climate-change related predictors in RMSE, MAPE, and PCC.

Fold	General			Wintertime
	RMSE	MAPE	PCC	RMSE	MAPE	PCC
1	1819.07	10.91	0.92	2560.01	15.66	0.94
2	1569.33	9.57	0.93	1448.72	8.11	0.91
3	851.30	4.93	0.93	1069.90	5.90	0.95
4	919.16	5.35	0.91	1181.79	6.23	0.93
5	1029.53	5.85	0.87	958.65	5.23	0.91
6	1165.02	7.35	0.87	858.01	4.84	0.93
7	1077.90	5.98	0.86	1265.33	6.65	0.89
8	1032.16	5.94	0.90	1096.53	5.97	0.89
Avg.	1182.93	6.98	0.90	1304.87	7.32	0.92

.

The results show that including climate change-related predictors negatively impacts the forecasting ability of the model, showing lower performance across all evaluation metrics on both the entire ($RMSE=1182.93$, $MAPE=6.98$, $PCC=0.90$) and wintertime test data sets ($RMSE=1304.87$, $MAPE=7.32$, $PCC=0.92$).

Subsequently, further research was conducted on the impact of the inclusion of individual climate change-related predictors on the performance of the model. Due to computational considerations, the data set was split into seven years of training data and two years of testing data (roughly a training–test split of 75–25%), rather than applying cross-validation. Of these two years of testing data, about 25% were considered wintertime.

Table 5. Performance of Prophet-LSTM trained with individual climate-change predictors in RMSE, MAPE, and PCC.

Predictor	Test Set	RMSE	MAPE	PCC
dailyPrecipitation	Entire	1020.0	5.92	0.86
	Winter	947.45	5.09	0.90
dailyMeanTemperature	Entire	1035.79	5.74	0.87
	Winter	991.35	5.12	0.90
dailyMeanWindspeed	Entire	988.06	5.84	0.87
	Winter	938.96	5.13	0.91
NAO	Entire	1079.08	6.01	0.86
	Winter	993.43	5.38	0.89
SIC	Entire	1053.74	6.69	0.88
	Winter	1009.51	5.89	0.91
SST	Entire	999.13	5.75	0.87
	Winter	902.04	4.82	0.91
geopotential	Entire	1113.94	7.02	0.87
	Winter	1042.45	5.94	0.91
renewableEnergyRatio	Entire	1215.77	7.96	0.86
	Winter	1069.89	6.30	0.90

shows the impact of evaluating the Prophet-LSTM model trained with the individual climate change-related predictors. In comparing the performance metrics obtained during evaluation from the model including all climate change-related features ($RMS{E}_{general}=1182.93$, $MAP{E}_{general}=6.98$, $PC{C}_{general}=0.90$, $RMS{E}_{winter}=1304.87$, $MAP{E}_{winter}=7.32$, $PC{C}_{winter}=0.92$), it is notable that the performances of the models with the individual predictors are now much closer to the model excluding any climate change-related predictors ($RMS{E}_{general}=934.13$, $MAP{E}_{general}=5.23$, $PC{C}_{general}=0.93$, $RMS{E}_{winter}=910.22$, $MAP{E}_{winter}=5.55$, $PC{C}_{winter}=0.92$).

To understand why using individual predictors improves performance, multicollinearity amongst the climate change-related variables was considered. Multicollinearity occurs when strongly correlated predictors are collectively used in a model, making it difficult for the model to determine the individual effects of the predictors when incorporating them collectively. This leads to overfitting of the model due to possibly inaccurate estimates of the predictor’s coefficients [49]. The Variance Inflation Factor (VIF) was used to measure the multicollinearity between the climate change-related predictors. The results indicate that the predictors dailyPrecipitation, NAO, and SIC exhibit nearly no multicollinearity with VIF scores of $VIF=\pm 1$, whereas the predictors dailyMeanWindspeed and dailyMeanTemperature display moderate to high multicollinearity with VIF scores of $VIF=\pm 10$. The predictors SST, geopotential, and renewableEnergyRatio show extremely high multicollinearity with VIF scores of $VIF>42$ [50]. The VIF scores are depicted in Figure 12.

4.4.2. Energy Price Predictors

To investigate the impact of using energy price predictors (SRQ3), the best model was tested by including the energy price predictors.

Figure 13 depicts the performance of including the energy price predictors in the Prophet-LSTM model, on both the entire test set and the wintertime subset. The detailed performances are listed in

Table 6. Performance of Prophet-LSTM model trained with energy price predictors in RMSE, MAPE, and PCC.

Fold	General			Wintertime
	RMSE	MAPE	PCC	RMSE	MAPE	PCC
1	704.70	4.56	0.96	812.67	4.53	0.95
2	1114.92	7.13	0.95	1248.88	7.28	0.92
3	697.76	3.97	0.96	777.21	3.81	0.95
4	710.90	4.11	0.95	873.96	4.79	0.95
5	1065.19	6.23	0.86	987.40	5.79	0.91
6	976.60	6.06	0.89	873.03	5.31	0.94
7	947.21	5.60	0.88	927.27	5.00	0.91
8	948.67	5.54	0.90	960.23	5.17	0.91
Avg.	895.74	5.40	0.92	932.58	5.21	0.93

. The results show that the addition of energy price predictors is very similar to that of excluding them. A modest improvement in the RMSE is seen for the entire test set including price predictors ($RMS{E}_{avg}=895.74$, $MAP{E}_{avg}=5.40$, $PC{C}_{avg}=0.92$) compared to that excluding them ($RMS{E}_{avg}=934.13$, $MAP{E}_{avg}=5.23$, $PC{C}_{avg}=0.93$). Alternatively, the wintertime forecasts including price predictors ($RMS{E}_{avg}=932.58$, $MAP{E}_{avg}=5.21$, $PC{C}_{avg}=0.93$) show slightly worse results compared to excluding the price predictors ($RMS{E}_{avg}=910.22$, $MAP{E}_{avg}=5.55$, $PC{C}_{avg}=0.92$), but the differences are minor.

4.5. Longer Horizon Forecasts

The feasibility of forecasting the energy demand for more than two months in advance (RQ) was investigated using an iterative process, iterating through lags of 60 to 180 days (with 15-day increments) to evaluate its impact on all performance metrics. Due to computational considerations, the same trainingtest split of 75–25% was used as in Section 4.4.1 to conduct further research.

Figure 14 illustrates the trends in the metrics across different lags. The detailed performance values of each lag can be found in

Table 7. Performance of Prophet-LSTM model with varying lags in evaluation metrics RMSE, MAPE and PCC.

Days Ahead	General			Wintertime
	RMSE	MAPE	PCC	RMSE	MAPE	PCC
60	986.69	5.72	0.87	926.11	5.06	0.91
75	980.67	5.76	0.88	895.73	4.90	0.91
90	898.11	5.18	0.90	696.31	3.73	0.95
105	939.55	5.52	0.89	771.91	4.18	0.94
120	1529.96	8.76	0.82	1260.11	6.54	0.86
135	1301.76	7.29	0.84	1182.48	6.27	0.86
150	1112.05	6.25	0.85	1078.16	5.91	0.87
165	1053.46	6.06	0.85	1043.45	5.74	0.88
180	1033.69	5.98	0.86	949.45	5.10	0.90

. The results show that overall, the performance generally decreases with increasing forecast horizons, exhibiting a slight increase again for lag 180. The wintertime forecasts perform competitively with and even slightly outperform the general forecasts. The best-performing lag is observed at the 90-day lag ($RMS{E}_{general}=898.11$, $MAP{E}_{general}=5.18$, $PC{C}_{general}=0.90$, $RMS{E}_{general}=696.31$, $MAP{E}_{general}=3.73$, $PC{C}_{general}=0.95$), showing improved performance metrics compared to all other lags.

5. Discussion

This section reflects on the obtained results, compares them to SOTA approaches, and provides limitations for this study.

5.1. Baseline MLR vs. Prophet-LSTM

The results demonstrate that Prophet-LSTM consistently outperforms the MLR model across all predictor selections and model setups. Interestingly, the Prophet-LSTM model performs best when predictors are excluded, whereas the MLR model benefits from their inclusion. This difference arises from the models’ natures: Prophet-LSTM combines Prophet’s ability to capture seasonality and trends with LSTM’s ability to model complex patterns. The additional predictors lead Prophet-LSTM to overfit, capturing noise rather than underlying patterns. In contrast, MLR benefits from the predictors since it is not able to fit such complex patterns and is thus less likely to overfit. Instead, including the predictors improves MLR’s forecast accuracy. This finding aligns with [25], which emphasized LSTMs’ strength in capturing nonlinearities compared to LR models.

Another important finding is that Prophet-LSTM Model 1 outperforms Model 2’s forecasts for both predictor selections in the wintertime forecasts. Contrarily, for general forecasts, Model 2 shows better performance compared to Model 1 when predictors are included. This difference can be explained by the use of dropout layers in Model 2. These likely mitigate overfitting caused by the large dimensionality of the data due to the numerous predictors. Likely, the overfitting is not an issue when forecasting on the wintertime test set given its lower data variability.

Unlike the model in [9], the Prophet-LSTM model in this work was extensively validated for its ability to generalize, showing stable performance across the different folds, especially when predictors are excluded. The general and wintertime forecasts of Prophet-LSTM are quite similar, indicating that the models can generalize well given their similar performances on different test sets.

Comparing the wintertime performances of the MLR and Prophet-LSTM to SOTA approaches using a similar MLR setup in [9] shows that both models outperform the paper’s approach consistently, comparing PCCs of $0.69$ and $0.70$ for the MLR model and PCCs ranging from 0.90 to 0.93 for Prophet-LSTM, compared to the obtained PCC of $0.50$ in the reference paper.

5.2. Impact of Adding Climate Change-Related and Energy Price Predictors

An important finding is that incorporating energy price predictors in the Prophet-LSTM model results in a slight improvement in performance. This finding aligns with past research [51], and considering that energy prices have been increasing over the past years, the expectation is that incorporating energy prices when modeling with more recent data will show even more impact from adding this predictor.

Furthermore, it was observed that climate change-related predictors negatively impact the performance of the Prophet-LSTM model, as shown in Table 4 in Section 4.4.1. Upon further investigation, it was found that including the predictors collectively results in poorer performance compared to including the predictors individually. Further analysis using the VIF revealed multicollinearity amongst some of the climate change-related predictors. This makes it difficult for the model to determine the individual effects of the predictors when incorporating them collectively. In turn, this leads to overfitting of the model due to possibly inaccurate estimates of the predictor’s coefficients [49]. Overall, this illustrates the importance of careful predictor selection.

5.3. Lead Time Improvement

Noteworthy is the performance of the 60-day lag in the lead time experiment conducted on a 75–25% training–test split, which compares well with previously obtained results through ROCV. This indicates that the model generalizes well, as forecasting on a larger test set with more years and data variance yields similar results to forecasting on a smaller test set. Moreover, the results reveal that longer forecasting horizons exhibit lower performance, except for the 90-day lag, which demonstrates increased performance compared to the 60-day lag as consistently used in this paper. Although performance tends to decrease for longer lags, a slight performance improvement is observed at the 180-day lag. The performance at the 90-day lag is noteworthy given the insignificant autocorrelation at the 3-month lag (roughly 90 days) in the monthly ACF in Figure 8 and might be caused by uncaptured seasonalities in the data. In contrast, the improvement at the 180-day lag corresponds with the significant autocorrelation observed at the 6-month lag in the monthly ACF in Figure 8. Overall, the findings show improved performance on the 60-day lag compared to similar works [9].

5.4. Limitations

Some limitations of this work are important to consider. Firstly, one limitation arises from the different data granularities. While the target variable was at an hourly level, not all data were available at the same granularity. During ROCV, some training folds contained predictors with minimal variance, such as half-yearly values, which might not fully represent reality. However, performance did not vary much when comparing ROCV results and 75–25% test split results, in which case the training set resembles reality more closely. This weakens the impact of this limitation to some extent, though it should still be considered. Furthermore, when comparing the findings to baseline work [9], which uses hourly energy consumption data as opposed to the hourly data used in this paper, it may not be entirely fair to directly compare correlation factors. While the differences in granularity could account for some variability, the magnitude of the observed differences suggests that other factors very likely contribute as well. Secondly, the primary focus on wintertime forecasts limited the amount of available data for training and testing. Since a year only contains three winter months, it is hard to form a test set large enough while also retaining sufficient training data. The limited amount of data remains a concern but was mitigated using ROCV to validate the model on different winters across different years, as well as validation on a general test set to assess how well the model generalizes. Thirdly, multicollinearity among predictors is a potential limitation. Model performance may have varied if the predictors had not displayed this behaviour, potentially showing enhanced performance. Finally, while the performance of the Prophet-LSTM model has been extensively researched within the context of the Dutch energy market, its generalizability may be limited. Despite testing the model’s performance on both wintertime and a general test set, its applicability to other countries remains untested. This limitation arises from the model’s specific tailoring to the unique characteristics of the Dutch energy market, in addition to differing energy demand factors and challenges faced in other countries.

6. Conclusions

In conclusion, this study sheds light on accurate long-horizon EDF in the Netherlands. EDF is crucial for maintaining stable and affordable energy supplies, especially for vulnerable populations most affected by shortages and high costs, as emphasized by multiple studies [9,15]. With potential energy shortages such as those Tennet is concerned with, rising energy prices, and the impact of climate change on the energy demand complicating today’s energy market, the importance of forecasts being very accurate and readily available upfront to mitigate the effects of electricity shortages, rising energy prices, and climate change is highlighted. The results demonstrate that the accuracy and lead time of energy demand forecasts in the Netherlands can be enhanced to a PCC of $r=0.93$ and lead time of 180 days, compared to the PCC of $r=0.50$ and lead time of 2 months observed in reference work [9]. Additionally, the Prophet-LSTM model shows significant improvement compared to a baseline MLR model in terms of RMSE, MAPE, and PCC. It consistently outperforms the MLR model, demonstrating its ability to capture complex patterns in the data. Furthermore, collectively incorporating climate change-related features illustrated poorer performance as opposed to incorporating them individually. This finding highlights the importance of careful prediction selection, trying to avoid including predictors that exhibit multicollinearity. Conversely, including energy price predictors yielded slightly improved energy demand forecasts. Aligning with existing research [23], this finding gives valuable insight for enhancing forecasting models in the future.

Despite the aforementioned limitations, this paper provides valuable insights into enhancing the accuracy and lead time of energy demand forecasts, especially in the context of the Netherlands. The societal impact of these findings is important to note. Accurately forecasting energy demand during the winter months is crucial for ensuring reliable energy supply, managing available resources, and mitigating the impact of extreme weather events caused by climate change on vulnerable populations. This is especially important with the prospect of potential future energy shortages and in times when climate change and energy prices play such an important role.

Future Work

For future work, a suggestion is to explore increasing the forecasting horizon, taking bigger lags such as 180 days, where it is especially important that lags are in line with the seasonal patterns observed in the data. Moreover, it would be interesting to see what the impact of higher granularity predictors would be on the forecasting ability of the models. Moreover, combining predictors that show multicollinearity into a single variable could enhance the forecasting ability of the model, as this allows the model to accurately estimate the predictors’ coefficients. Finally, the model’s generalizability could be improved by adapting it to fit a wider range of countries. This would allow for a better evaluation of the model’s ability to generalize across different energy environments.