Forecasting Thermal Demand in Citizen Energy Communitites Using Machine Learning: Application to the Urberoa Case Study

Forecasting Thermal Demand in Citizen Energy Communities Using Machine Learning: Application to the Urberoa Case Study Raul Saenz-Herreros1,2[0009−0005−7997−6624], Alejandro Muro-Belloso1[0009−0009−7044−5180], Ignacio Guisasola-Iparraguirre1[0009−0007−0412−7059], Ignacio Lázaro-Llorente1[0000−0003−3128−8150], and J. David Nuñez-Gonzalez[0000−0002−5047−1033] 1Fundación Tekniker, Basque Research and Technology Alliance (BRTA), Iñaki Goeanaga, 5, 20600 Eibar, Spain 2Applied Mathematics Department Engineering School of Gipuzkoa – University of the Basque Country - UPV/EHU Avda. Otaola, 29, 20600 Eibar, Spain Abstract. Accurate thermal energy demand forecasting is vital for optimizing district heating networks and ensuring efficient energy management in Citizen Energy Communities, particularly those incorporating decentralized renewable sources. This work introduces a data-driven methodology for shortand mid-term thermal demand prediction based on insights from the FEDECOM project. Historical consumption records and meteorological forecasts are used to segment seasonal patterns and apply targeted preprocessing techniques, including outlier detection and temporally coherent imputation of missing data. An extensive set of features is engineered to represent meteorological conditions, calendar cycles, and domain-specific variables. Feature relevance is evaluated, and predictive models are trained using ensemble learning techniques. Forecasts are generated at 15-minute resolution across 24-hour horizons using recursive and direct multi-step prediction strategies validated through cross-validation. The proposed method consistently achieves strong performance, with LightGBM outperforming baseline models in terms of mean absolute error, root mean square error, coefficient of determination, and relative error metrics. This approach offers a scalable, generalizable solution to improve forecasting accuracy in urban energy systems, supporting proactive planning, better renewable energy integration, and enhanced operational efficiency in sustainable environments. Keywords: Thermal Energy Demand ·Forecasting ·Machine Learning ·Feature Engineering ·FEDECOM Project. 1 Introduction The transition toward a sustainable energy system requires innovative strategies to reduce greenhouse gas emissions, improve energy efficiency, and integrate 2 R. Saenz-Herreros et al. renewable energy sources. In this context, District Heating (DH) networks and Citizen Energy Communities (CECs) have emerged as key enablers of the decarbonization objectives outlined in the European Green Deal [8]. CECs, formally recognized by the European Union, empower citizens to actively participate in the production, consumption, and management of energy, often incorporating distributed renewable generation [9]. A major operational challenge in these systems lies in the accurate forecasting of thermal energy demand, particularly for DH networks serving diverse buildings with fluctuating consumption profiles. Thermal demand is highly sensitive to several dynamic factors, including ambient temperature, solar radiation, wind speed, building insulation quality, and user behavior [10, 15]. These factors introduce significant temporal variability, complicating energy planning and control. Machine learning and data mining methods have shown promise in this domain due to their adaptability and scalability [1]. However, their success hinges on the quality and resolution of the input data. To address this, the present study proposes a hybrid, data-driven methodology for shortand mid-term thermal demand forecasting, developed within the framework of the FEDECOM project. The approach combines advanced preprocessing techniques, comprehensive feature engineering, and modern ensemble learning algorithms, with a focus on the Light Gradient Boosting Machine (LightGBM) model. The dataset comprises high-resolution (15-minute) thermal demand records from the URBEROA community [23] and weather forecasts from the Meteomatics API [16]. The proposed methodology is evaluated using real-world data across different seasons and the forecasts are generated over a 24-hour horizon using recursive, direct, and hybrid strategies, and validated through walk-forward crossvalidation. Results indicate strong predictive performance, demonstrating the method’s practical applicability in improving energy planning, optimizing resource efficiency, and supporting renewable integration within sustainable urban communities. The remainder of this paper is structured as follows: Section 2 presents the thermal demand dataset and preprocessing techniques. Section 3 describes the forecasting models and strategies, including the previous baseline based on Random Forest and LightGBM under recursive and direct schemes. Section 4 reports the experimental results and performance comparison. Finally, Section 5 summarizes the main findings and outlines future research directions. 2 Data Processing 2.1 Thermal Demand Data The dataset used in this study consists of historical thermal energy demand records (in kW) sampled every 15 minutes from January 2021 to November 2023. It captures the aggregate demand for both domestic hot water (DHW) and space heating during the winter season, and for DHW only during the summer season. Figure 1 illustrates the raw time series data segmented by winter and summer periods. Thermal Demand Forecasting in Citizen Energy Communities 3 Fig. 1: Thermal demand power data expressed in kilowatts ranging from January 2021 to November 2023. Winter and summer periods are marked in blue and red, respectively. Outliers and Missing Data A critical step in data preprocessing is the identification and treatment of outliers. In this work, a conservative approach based on the Interquartile Range (IQR) was applied, due to its robustness for time series data [24]. An upper threshold was defined as Q3+k·IQR (with k= 2 for winter and k= 4 for summer), while the lower threshold was fixed at zero as negative thermal demand is not physically possible. Very low values close to zero may be considered anomalous but not neccessarily incorrect. Outlier thresholds were computed using season-specific windows to avoid distortion from inter-seasonal variability. Although the number of flagged outliers was very low—well below 1% of the entire dataset—their removal was important to prevent distortion in model training and to preserve the temporal integrity of the time series. Figures 2 and 3 illustrate the outlier detection process for winter and summer periods, respectively. Missing values were identified at specific intervals throughout the dataset. Although they represented less than 3% of the total data, their presence could still negatively impact model performance, particularly in time series forecasting tasks. Therefore, these gaps were imputed using a multistage strategy that respects the temporal and seasonal structure of the data [17]. Isolated gaps were filled using moving average smoothing to preserve local continuity. For longer gaps, trend-based imputation was applied using values from the same hour on previous days. A small random variation was added to reflect natural variability. This methodology ensured a consistent and realistic reconstruction of both summer and winter time series. Temporal Structure To uncover underlying seasonalities and temporal correlations, a detailed exploratory time series analysis was conducted. Figure 4 presents seasonal decompositions by month and hour using boxplots. 4 R. Saenz-Herreros et al. Fig. 2: Detection of thermal demand outliers during the winter season using the Interquartile Range (IQR) method. Fig. 3: Detection of thermal demand outliers during the summer season using the Interquartile Range (IQR) method. Thermal Demand Forecasting in Citizen Energy Communities 5 (a) Thermal demand distribution by month. (b) Thermal demand distribution by hour of the day. Fig. 4: Seasonal decomposition of thermal demand highlighting intra-day and annual patterns relevant for feature engineering. The monthly analysis (Figure 4a) shows a clear annual pattern with peak demand in the colder months (November to February) and reduced consumption in the summer (June to September), aligned with the trends of heating use. This seasonal pattern reflects the heating-oriented nature of thermal demand and justifies the use of month-based cyclical features as well as interactions with meteorological inputs such as temperature and solar radiation. The hourly analysis (Figure 4b) reveals a pronounced intra-day seasonality with clear peaks in the morning (7–8 a.m.) and evening (8–9 p.m.), consistent with residential routines. These patterns clearly demonstrate that the thermal demand data exhibit multiple seasonalities and time-dependent correlations. Therefore, the integration of these structures is essential for capturing the temporal context. 2.2 Exogenous Variables Following the preprocessing of thermal demand data, a comprehensive set of exogenous variables was constructed to enhance the model’s ability to capture the complex interactions between external drivers and demand dynamics. These features capture environmental conditions, human behavioral patterns, and temporal structures that influence thermal energy usage. Meteorological Variables from Meteomatics A key component of the exogenous feature set consists of meteorological variables obtained via the Meteomatics API [16]. This platform provides access to high-resolution, spatiotemporally consistent atmospheric data at 15-minute intervals, perfectly aligned with the sampling frequency of the thermal demand data. The robustness of Meteomatics in terms of preprocessing, data continuity, and harmonization with official weather stations ensures high data reliability and usability for forecasting tasks. Table 1 summarizes the selected variables. 6 R. Saenz-Herreros et al. Table 1: Meteorological variables used as model inputs. Variable Units Description t_2m_C ◦C Outdoor temperature relative_humidity_2m % Relative humidity direct_rad_W W/m2Direct solar radiation diffuse_rad_W W/m2Diffuse solar radiation global_rad_W W/m2Global solar radiation dew_point_2m_C ◦C Dew point temperature Calendar and Cyclical Time Features To capture calendar-based seasonality and behavioral routines, features such as hour of day, day of week, week of year, and month were extracted. As these variables are cyclical, sine and cosine transformations were applied to preserve their periodic nature and ensure continuity for models sensitive to input geometry (e.g., hour 23 is just as close to 0 hour as the hour 22). These features align with the temporal dependencies identified in the seasonal analysis (Figure 4) and improve generalization across different times of day or year. Astronomical and Holiday Effects Solar-related variables including sunrise and sunset times, daylight duration, and daylight presence (binary flag) were computed based on the site’s geographic coordinates. These features reflect solar exposure, which influences heating dynamics and occupancy behavior. Additionally, a binary holiday flag and indicators for adjacent days were also included to model irregular usage patterns associated with holidays and long weekends. Sliding Window and Lag Features Short-term effects and autocorrelation were captured using lag-based features and rolling window aggregates (mean, max, min) applied to outdoor temperature. These features are critical to modeling delayed thermal responses, such as residual heat from prior periods or thermal inertia in buildings. Polynomial Interactions and Feature Expansion To model nonlinear dependencies and interaction effects, second-order polynomial features were engineered between key meteorological, temporal, and astronomical variables. However, to limit redundancy and reduce overfitting risk a post-processing step involving correlation filtering and feature selection was applied. This comprehensive feature engineering strategy resulted in a rich, interpretable, and domain-aligned set of predictors for thermal demand forecasting. Thermal Demand Forecasting in Citizen Energy Communities 7 3 Energy Demand Modeling 3.1 Baseline Approach: Random Forest In earlier work on thermal demand forecasting for the Urberoa community [21], a modeling strategy based on seasonal segmentation was adopted. The dataset was divided into summer and winter subsets, and separate Random Forest regressors were trained for each season. This approach aimed to better capture seasonal consumption behaviors and distinct relationships between input features and target values. Random Forest [12] is an ensemble learning method that constructs multiple decision trees during training and outputs their average prediction, improving accuracy and reducing overfitting [20]. The feature set included both temporal and meteorological variables, as well as lagged demand values at 24 and 168 hours, to capture daily and weekly autocorrelation patterns commonly observed in thermal energy consumption. However, the feature engineering step discussed previously was not taken into account. 3.2 Proposed Approach: Light Gradient Boosting Machine LightGBM (Light Gradient Boosting Machine) [19] is a state-of-the-art gradient boosting framework known for its efficiency and strong predictive performance. It builds additive models in a forward stage-wise fashion and trains decision trees sequentially, where each tree is optimized to reduce the residual errors of the previous one. In this work, LightGBM was chosen not only for its computational efficiency and modeling capabilities but also to enable a fair performance comparison with the Random Forest approach used in the previous approach. Furthermore, this method utilizes a single unified model trained on both summer and winter data and it learns seasonal effects directly from the input features. Two multi-step forecasting strategies were implemented using the Skforecast library [2] -recursive and directand the dataset was partitioned as follows to ensure that both winter and summer periods are present in each set: – Training: 2021-01-02 00:00:00 to 2023-01-31 23:45:00 – Validation: 2023-02-01 00:00:00 to 2023-08-01 23:45:00 – Test: 2023-08-02 00:00:00 to 2023-11-19 23:45:00 Multi-Step Recursive Forecasting In the recursive strategy (Figure 5), a single model is trained to predict the next 15-minute time step using the previous 672 lags (equivalent to one week of historical demand). Once trained, the model is recursively applied: the prediction for t+ 1 is fed back as input to forecast t+ 2, and so on, until completing the 24-hour horizon (96 steps) [25]. Model performance was evaluated using a backtesting strategy [3] that simulates daily operation by generating 24-hour forecasts at 23:59 for each day in the test period. Hyperparameter tuning was conducted via Bayesian optimization [11], using time series cross-validation [4] to maintain temporal integrity. Feature selection based on recursive elimination and model-based importance measures [22] further improved model generalization and reduced redundancy. 8 R. Saenz-Herreros et al. (a) Recursive predictions during winter. (b) Recursive predictions during summer. Fig. 5: Multi-step recursive forecasting: actual vs. predicted thermal demand for selected test intervals in both winter and summer. Multi-Step Direct Forecasting As an alternative, the direct multi-step strategy (Figure 6) was implemented to mitigate error accumulation over the prediction horizon. In this approach, a separate model is trained for each forecast step. This ensures that each prediction is based on observed values during training, avoiding the propagation of prior forecast errors. The implementation used the same input features and lag structures as the recursive model. Model evaluation was performed under the same backtesting and time series cross-validation conditions to allow a fair comparison. While this method improves prediction robustness [6], it requires training 96 distinct models (one per time step), significantly increasing computational cost. This trade-off must be considered when assessing practical deployment. 4 Results To ensure a consistent and robust evaluation of model performance, all error metrics were computed over the full set of actual and predicted values across the entire test horizon. This strategy preserves the original 15-minute resolution and allows accurate comparison of forecast quality across different timeframes. Thermal Demand Forecasting in Citizen Energy Communities 9 (a) Direct predictions during winter. (b) Direct predictions during summer. Fig. 6: Multi-step direct forecasting: actual vs. predicted thermal demand for selected test intervals in both winter and summer. Five widely accepted metrics were used to quantify predictive performance: Mean Absolute Error (MAE) [5], Mean Absolute Percentage Error (MAPE) [7], Coefficient of Determination (R2) [18], Root Mean Square Error (RMSE) [13], and the Coefficient of Variation of RMSE (CVRMSE) [14]. The comparative results show that both recursive and direct implementations of the LightGBM model significantly outperform the Random Forest baseline, while also eliminating the need for seasonal model segmentation (i.e., separate models for summer and winter). Winter Forecasting Results As shown in Table 2, Random Forest delivered the weakest performance in winter, with a MAE of 150.28 kW, RMSE of 204.63 kW, and a relatively low R2 of 0.58. In contrast, both LightGBM configurations achieved lower errors and improved explanatory power. LightGBM Recursive and Direct models produced MAEs around 113 kW, RMSEs near 162 kW, and R2values of 0.69, reflecting a substantial improvement over the baseline. The CVRMSE was also reduced by 4 percentage points.