Benchmarks and Custom Package for Energy Forecasting

## Abstract

Abstract Energy (load, wind, photovoltaic) forecasting is significant in the power industry as it can provide a reference for subsequent tasks such as power grid dispatch, thus bringing huge economic benefits. However, there are many differences between energy forecasting and traditional time series forecasting. On the one hand, traditional time series mainly focus on capturing characteristics like trends and cycles. In contrast, the energy series is largely influenced by many external factors, such as meteorological and calendar variables. On the other hand, energy forecasting aims to minimize the cost of subsequent tasks such as power grid dispatch, rather than simply pursuing prediction accuracy. In addition, the scale of energy data can also significantly impact the predicted results. In this paper, we collected large-scale load datasets and released a new renewable energy dataset that contains both station-level and region-level renewable generation data with meteorological data. For load data, we also included load domain-specific feature engineering and provided a method to customize the loss function and link the forecasting error to requirements related to subsequent tasks (such as power grid dispatching costs), integrating it into our forecasting framework. Based on such a situation, we conducted extensive experiments with 21 forecasting methods in these energy datasets at different levels under 11 evaluation metrics, providing a comprehensive reference for researchers to compare different energy forecasting models.

## 1 Introduction

Time series data are becoming ubiquitous in numerous real-world applications (Wen et al., ; Lai et al., ; Zhou et al., ; Wang et al., ). Among them, energy forecasting is crucial for maintaining the supply and demand balance in the power system. Thanks to the development of machine learning in recent years, various methods have been developed for energy forecasting (Yildiz et al., ; Zhang et al., ; Aslam et al., ; Benti et al., ; Yang et al., ; Zhu et al., ). To further promote the development of this field, many energy forecasting competitions like the Global Energy Forecasting (GEF) Competition have been held over the years (Hong et al., ; ; ). In addition, many competitions target specific themes, like building energy management based on electricity demand and solar PV generation (Nweye et al., ), and the impact of COVID-19 issues on the power systems (Farrokhabadi et al., ).

Although many advanced energy forecasting methods have emerged in the past decades, the winners of forecasting competitions often use conventional machine learning models (like non-deep learning models). The secret to their victory lies in targeted feature engineering, which is also the major difference between energy forecasting and general time series forecasting (Sobhani et al., ). However, no existing benchmarks focus on those parts of energy forecasting. Compared with other time series, energy data are greatly affected by external factors such as temperature and calendar variables, making it challenging to model the dynamics accurately. Therefore, exploring the impact of external factors on energy forecasting has always been an important research direction in this field (Aprillia et al., ). For load data, temperature is considered to have a significant impact and many researchers have focused on using temperature variables to assist in constructing load forecasting models (Haben et al., ; Sobhani et al., ; Liu et al., ). At present, the utilization of temperature variables can be roughly divided into two strategies. One is to make targeted transformations on temperature variables, which are often based on relatively simple statistical learning methods (Guan et al., ; Farfar & Khadir, ). The other one is to extract features by neural networks. Such models usually achieve better accuracy (Imani, ; Hafeez et al., ). However, the interpretability of this kind of model decreases due to the black-box characteristic of neural networks. Nevertheless, related feature engineering also has a guiding role for neural network-based forecasting models. Currently, no large-scale experimental results have been provided to demonstrate this. Therefore, we will provide various related feature engineering in our package and discuss the impact on load forecasting models based on temperature feature engineering. As for renewable energy data, the impact of external factors is even more significant, such as the effect of wind speed on wind power. However, unlike load data, a common problem in the community is that we often lack meteorological external factors that match the renewable energy data. As stated in (Menezes et al., ), an excellent renewable energy dataset not only requires generation data but also corresponding meteorological data, which is often difficult. For this situation, we have collected and publicly released a renewable energy output dataset from a first-level administrative region, covering different types of renewable energy data such as PV,onshore wind, and offshore wind, as well as corresponding key meteorological factors such as real-time wind speed at wind turbine hubs.

Apart from feature engineering, another difference is that the most important concern of energy forecasting models is to achieve the lowest cost instead of the best accuracy of predictions. Due to the diversity of the time series data, general time series forecasting results are rarely optimized for a specific task. However, energy forecasting results mainly be used for subsequent power grid dispatch, which inspires us to pay attention to the relationship between the prediction and subsequent decision-making cost. (Wang & Wu, ) discovered the asymmetry between cost and load forecasting error and this asymmetry comes from actual scenarios. Underestimating predictions at peak points will result in additional power purchase costs while overestimating predictions at low points will waste power generation. Therefore, bias will be introduced if we just use traditional symmetric loss functions like MSE and MAE to train the model. Then, (Zhang et al., ) proposed to use the characteristics of piecewise linearization and the Huber function to model the relationship between forecasting error and real cost. Inspired by this work, our package provides methods for modeling the relationship between forecasting error and other variables and then constructing the corresponding loss function.

To provide an accessible and extensible reference for future researchers in energy forecasting, our developed package differs from the existing time series packages (Alexandrov et al., ; Godahewa et al., ). Specifically, our package splits the entire energy forecasting process into five modules: data preprocessing, feature engineering, forecasting methods, postprocessing, and evaluation metrics.
Our package will cover both probabilistic forecasting and point forecasting, providing feature engineering methods and predictors based on traditional machine learning models and deep learning models. Users can combine any of these components and obtain their customized models. Furthermore, our package adds specific functionalities to address the characteristics of energy forecasting like load forecasting and its differences from traditional time series forecasting, greatly enhancing the user’s flexibility to construct energy forecasting models.

Lastly, we conduct extensive experiments to evaluate both point and probabilistic forecasting performance of different models on multiple load datasets at different levels. As for renewable energy, due to the influence of meteorological factors, renewable energy often has significant volatility and uncertainty. Therefore, we mainly consider the probabilistic forecasting of renewable energy. These results not only provide insights about different forecasting models under various scenarios and evaluation metrics but also show the accessibility and extensibility of our package and benchmark.

We summarize our primary contributions as follows:

- 1.
A new open-source high-quality renewable energy dataset. We release a renewable energy output dataset, covering different types of data such as PV,onshore, and offshore wind, as well as corresponding key meteorological factors. Compared with other datasets, our advantage lies in providing different types of renewable energy generation data like PV and offshore wind, at different levels within the same region, including data from specific stations and some rather big regions (with multiple stations). Meanwhile, we also provide the relative positional relationships between different stations. It is beneficial for the community to study the relationship between different renewable energy types under similar climate conditions in a certain area. To our knowledge, this is the first high-quality renewable energy dataset with these characteristics.
- 2.
Domain-specific feature engineering and self-defined loss function for load forecasting. Based on the characteristics of load, temperature, and calendar variables, we integrate the feature engineering that reflects the ternary relationship into our package for users to use in any forecasting model. At the same time, we also provide users with a function to customize the loss function. Users can define the relationship between the forecasting error and any variable (such as the dispatching cost of the power grid) and integrate it into our forecasting framework as a loss function. In our experiment, we simulate an IEEE 30-bus system, provide the relationship between simulated forecasting error and cost, and construct the corresponding loss function.
- 3.
Fully open-source platform with accessibility and extensibility. We release the relevant code publicly (^1^11[https://anonymous.4open.science/r/ProEnFo-17CC](https://anonymous.4open.science/r/ProEnFo-17CC)). Users can freely combine the components we provide to design their energy forecasting framework to cope with different energy forecasting scenarios. We provide over 20 forecasting methods, including both probabilistic forecasting and point forecasting. In addition, we also encourage users to combine their forecasting methods based on our framework. After defining the input and output of their method, users only need one command to add their forecasting model to our framework to accomplish common 24-hour-ahead electricity forecasting tasks. At the same time, we also provide various evaluation and visualization methods to facilitate users to evaluate the predictive performance of different models from multiple perspectives.
- 4.
The first large-scale benchmark for energy forecasting. Based on 12 energy datasets in different scales (e.g., building-level and aggregated-level), we conduct experiments on both probabilistic forecasting and point forecasting for load forecasting. For probabilistic forecasting, we compare 21 probabilistic forecasting methods and apply 11 metrics to comprehensively compare their performance. For point forecasting, we focus on 12 widely used deep learning methods and compare the traditional MSE loss function with our proposed loss function based on the relationship between forecasting error and cost. For the renewable energy dataset, we conduct probabilistic forecasting with 12 methods. To the best of our knowledge, this is the first work to construct comprehensive benchmarks with large-scale datasets for energy forecasting scenarios.

## 2 Data Description

In this section, we will introduce how our data is collected and the characteristics of the datasets, including electrical load datasets and the renewable energy dataset.

### 2.1 Electrical load dataset

We have collected a total of 11 load datasets for our collection, and a detailed description of each dataset is provided in the Appendix A. In summary, the load data we collect mainly comes from UCI machine learning databases (Dua & Graff, ), Kaggle data competition platforms (Nicholas, ; Yeafi, ; Shahane, ), and the famous global energy forecasting competitions (Hong et al., ; ; ). In addition, we also include a dataset reflecting the impact of the COVID-19 epidemic on the power system in our archives. Under the influence of COVID-19, an influential external factor, the power load has changed significantly, posing a challenge to the robustness of the forecasting model (Farrokhabadi et al., ). From the perspective of load hierarchy, 7 of the data we collect are aggregated-level datasets, and the remaining 4 are building-level datasets. Aggregated-level load refers to the total load that aggregates multiple independent loads (such as the power demand of various buildings in the power system) together. More specifically, we classify the load of an area greater than one building as aggregated level. Because the aggregated-level load results from multiple load aggregations, it typically exhibits more pronounced periodicity and seasonality. For this reason, calendar variables significantly impact load forecasting at this level. In contrast, the load of the building level, which can also be seen as a part of the aggregated load, change very dramatically, resulting in significant uncertainty. Therefore, many works related to building-level load forecasting often focus on probabilistic forecasting (Xu et al., ; Jeong et al., ). To provide a reference for researchers in related fields, we also collect building-level datasets from the Building Data Genome 2 (BDG2) Data-Set (Miller et al., ). In addition to different levels, the data we collect also has the characteristic of almost covering all meteorological data (actual measurement) such as temperature, which may be greatly beneficial to forecasting because of the great impact of external variables (especially temperature) on load. The number of time series contained in each dataset and their corresponding features are listed in Table 1. And all the data will be released under appropriate licenses.

**Table 1: Datasets in the load forecasting archive.**
|  | Dataset | No. of series | Length | Resolution | Missing | License | Type | External variables |
| --- | --- | --- | --- | --- | --- | --- | --- | --- |
| 1 | Covid19(Farrokhabadi et al., ) | 1 | 31912 | hourly | No | CC BY 4.0 | aggregated-level | airTemperature, Humidity, etc |
| 2 | GEF12(Hong et al., ) | 20 | 39414 | hourly | No | CC BY 4.0 | aggregated-level | airTemperature |
| 3 | GEF14(Hong et al., ) | 1 | 17520 | hourly | No | CC BY 4.0 | aggregated-level | airTemperature |
| 4 | GEF17(Hong et al., ) | 8 | 17544 | hourly | No | CC BY 4.0 | aggregated-level | airTemperature |
| 5 | PDB(Yeafi, ) | 1 | 17520 | hourly | No | CC0 1.0 | aggregated-level | airTemperature |
| 6 | Spain(Nicholas, ) | 1 | 35064 | hourly | Yes | CC0 1.0 | aggregated-level | airTemperature, seaLvlPressure, etc |
| 7 | Hog(Miller et al., ) | 24 | 17544 | hourly | Yes | CC BY 4.0 | building-level | airTemperature, wind speed, etc |
| 8 | Bull(Miller et al., ) | 41 | 17544 | hourly | Yes | CC BY 4.0 | building-level | airTemperature, wind speed, etc |
| 9 | Cockatoo(Miller et al., ) | 1 | 17544 | hourly | Yes | CC BY 4.0 | building-level | airTemperature, wind speed, etc |
| 10 | ELF(Shahane, ) | 1 | 21792 | hourly | No | CC BY 4.0 | aggregated-level | No |
| 11 | UCI(Dua & Graff, ) | 321 | 26304 | hourly | No | CC BY 4.0 | building-level | No |

### 2.2 Renewable Energy dataset

We collect and release a renewable energy dataset from a first-level administrative region. For station-level, there are 10 onshore wind series with its corresponding wind speed and wind direction at the wind turbine hubs(1 only has the wind speed), 1 offshore wind series with wind speed and wind direction, and 10 PV series with corresponding irradiance and air temperature. For the region-level, we provide 16 onshore wind series, 8 offshore wind series, and 13 PV series in different cities. In addition, 5 of onshore wind series, 4 of offshore wind series, and 5 of PV series in different regions are also provided(as summarized in Table 2 and the detailed description can be found in the Appendix A).

**Table 2: Datasets in the renewable energy forecasting archive.**
|  | Data | No. of series | Type | Level | Missing | Resolution | External variables |
| --- | --- | --- | --- | --- | --- | --- | --- |
| 1 | LW_S | 10 | Onshore Wind | Station | No | 15 min | Wind speed, Wind direction, Location |
| 2 | OW_S | 1 | Offshore Wind | Station | No | 15 min | Wind speed, Wind direction, Location |
| 3 | PV_S | 10 | PV | Station | No | 15 min | Irradiance, Temperature, Location |
| 4 | LW_C | 16 | Onshore Wind | City | No | 15 min | - |
| 5 | OW_C | 8 | Offshore Wind | City | No | 15 min | - |
| 6 | PV_C | 13 | PV | City | No | 15 min | - |
| 7 | LW_R | 5 | Onshore Wind | Region | No | 15 min | - |
| 8 | OW_R | 4 | Offshore Wind | Region | No | 15 min | - |
| 9 | PV_R | 5 | PV | Region | No | 15 min | - |

## 3 Package Functions

### 3.1 Overview of the package

Fig. 1 shows the overview of our package. As stated before, we divide the overall forecasting process into several parts to address potential issues in energy forecasting for the power industry. First of all, energy data like load is obtained by physical devices such as electricity meters. During this process, it is inevitable to encounter missing values, omissions, and other situations. In this regard, our package provides various methods such as ARIMA based on Kalman filtering (Harvey & Pierse, ), K-nearest neighbor algorithm (García-Laencina et al., ) to fill in missing values, ensuring minimum data information distortion. Secondly, our model provides a variety of feature selection strategies to meet the needs of different scenarios. For example, users can choose the corresponding data from the previous seven days for day-ahead forecasting or use Autocorrelation Function (ACF) and Partial Autocorrelation Function (PACF) metrics to help select the lagged values. In addition, our framework allows users to add external variables such as temperature and calendar variables that may impact the forecasting model. As for the forecasting methods, we provide both probabilistic forecasting and point forecasting methods. Among them, probabilistic forecasting will be based on quantile forecasting. However, quantile regression may lead to confusion about quantile, that is, the forecasting result of a larger quantile is smaller than that of a smaller quantile. To address this situation, we have provided corresponding post-processing for reordering.

After obtaining forecasting results, we provide extensive metrics for evaluation. Existing forecasting packages generally provide some metrics, such as Pinball Loss and CRPS. Although they can evaluate the quality of forecasting results, they tend to evaluate the forecasting performance as a whole and fail to provide detailed guidelines for further improvement.
For example, a model may perform poorly in a certain quantile while performing well in other quantiles. To more intuitively compare the performance of models in different quantiles, our package provides the matrix visualization function of multiple metrics in different quantiles. And the evaluation metrics we have implemented include CalibrationError (Chung et al., ), WinklerScore (Barnett, ), CoverageError, and so on (details can be seen in Appendix Section D.3).

Figure: Figure 1: Overview of the energy forecasting package.
Refer to caption: x1.png

### 3.2 Feature Engineering and Custom Loss function for Electrical load forecasting

#### 3.2.1 Temperature-calendar feature engineering strategy

The impact of temperature on load is greatly influenced by calendar variables. Inspired by the linear regression model (Hong et al., ), we apply one-hot encoding to calendar variables and then model this coupling relationship by taking their products with temperature to the first, second, and third powers as features. The specific formula is as follows:

$$ $\displaystyle\hat{y}_{t}=$ $\displaystyle\beta_{0}+\beta_{1}\underbrace{\mathrm{Trend}_{t}}_{\text{ replaced by sequence model}}+\beta_{2}M_{t}+\beta_{3}W_{t}+\beta_{4}H_{t}+ \beta_{5}W_{t}H_{t}+\beta_{6}T_{t}+\beta_{7}T_{t}^{2}$ $\displaystyle+\beta_{8}T_{t}^{3}+\beta_{9}T_{t}M_{t}+\beta_{10}T_{t}^{2}M_{t}+ \beta_{11}T_{t}^{3}M_{t}+\beta_{12}T_{t}H_{t}+\beta_{13}T_{t}^{2}H_{t}+\beta_{ 14}T_{t}^{3}H_{t},$ $$

where $M_{t}$, $W_{t}$, $H_{t}$, and $T_{t}$ represent the month, workday, and hour vectors after one-hot encoding and temperature at the corresponding time. Due to the one-hot coding, one categorical variable is changed into multiple binary categorical variables. When the corresponding variable is 0, the parameters of the linear model will not have any effect on it. Therefore, the result of doing so is constructing multiple personalized models based on calendar variables. Such a feature engineering strategy can help the forecasting model cope with situations where the temperature and load relationship shifts under different calendar variables. To preserve such characteristics and integrate existing sequence modeling methods (such as LSTM, and N-BEATS), we treat the information extracted by sequence modeling methods as trend variables and concatenate them with the previously obtained calendar temperature coupling variables. Finally, a fully connected layer is used to map the final output result. In section 5, we will compare the impact of this feature engineering on forecasting results across multiple datasets.

#### 3.2.2 Custom loss function

Based on (Zhang et al., ), our framework provides corresponding piecewise linearization functions to help users model the relationship between forecasting errors and real requirements (such as scheduling costs) and integrate it into the gradient descent training. Specifically, we need data pairs $(\epsilon_{i},C_{i})_{i=1,\dots,N}$, where $\epsilon_{i}$ is the forecasting error and $C_{i}$ is the real requirement(here we mainly refer to the dispatching cost). Here, we consider using Forecasting Error Percentage (FEP) $\epsilon_{i}=\frac{f\left(x_{i}\right)-y_{i}}{y_{i}}$ as our error metric. At the same time, we normalize $\{C\}_{i=1,\dots,N}$, making its value fall between 0 and 1. Now, our goal has become how to construct $L(\epsilon)$ to estimate $C$. To achieve high fitting accuracy, we can use a spline cubic function, denoted as $s$, to fit it. However, the disadvantage of doing so is that there will be many discontinuities, which is not convenient to integrate them into our forecasting framework as a loss function. To ensure the fitting ability of the function while making it as simple as possible, a piecewise linearization strategy is adopted here. The number of segments K can be determined by setting the upper bound of the fitting error
$\left\|s-L(\epsilon)\right\|_{2}\leq\frac{\left(\int_{\epsilon_{\min}}^{
\epsilon_{\max}}s^{\prime\prime}(\epsilon)^{\frac{2}{5}}d\epsilon\right)^{
\frac{5}{2}}}{\sqrt{120}K^{2}}$ (Berjón et al., ).
Regarding the position of the corresponding interval points, we strive to distribute the data points within the interval formed by each pair of endpoints (De Boor & De Boor, ). So far, we have obtained a piecewise linearization function. To take it as a loss function, we need to ensure its differentiability. Specifically, we use a quadratic function in a cell around each breakpoint to smooth it. Note that the quadratic function does not need to fit the data, but only needs to ensure its left and right continuity and the continuity of the corresponding first derivative to obtain the parameters.

Figure: Figure 2: Visualization of the simulated loss function.
Refer to caption: x2.png

Adapted from (Zhang et al., ), we employ a modified IEEE 30-bus test system to simulate an economic dispatch optimization problem (Hota & Naik, ). This system comprises 6 generators and 3 Battery Energy Storage Systems (BESS) connected to the network to supply power to 21 loads. We will address two dispatch optimization problems: the Day-Ahead Economic Dispatch (DAED) problem and the Intra-day Power Balancing (IPB) problem, to obtain the dispatch cost. The DAED problem primarily focuses on dispatching the power output of all online generators based on forecasted results, while the IPB problem fine-tunes the outputs of generators, charging/discharging behavior of BESS, and load shedding to balance the intra-day load deviation. The simulation process can be found in Appendix C.1. The simulated visualization result can be found in Fig. 2. From the simulation results, we can identify two key distinctions in comparison to traditional Mean Squared Error (MSE) losses. First, there is a significant asymmetry between forecasting error and dispatching cost, which is evidenced by the considerable discrepancy in dispatching costs between instances when predicted values are higher than the true value and instances when predicted values are lower than the true value. This difference stems from varying sources of dispatching costs. When the predicted value is too high, the system may over-dispatch power generation equipment, whereas when the predicted value is too low, the system may need to purchase emergency power supplies. Second, to minimize dispatching costs, we often do not need to require that predicted values perfectly match true values. Instead, we only need to meet specific accuracy requirements to minimize dispatching costs. This finding can guide the model to focus more on data points with low forecasting accuracy, thereby minimizing the final dispatching cost. We have packaged the loss function based on these simulated data for users to call.

### 3.3 Usage of the package

Our package provides multiple components for users to freely combine and construct various power forecasting scenarios. Meanwhile, we also encapsulate the default strategy, where users only need one command to achieve the default 24-hour ahead power probabilistic forecasting task and obtain detailed result files, including multiple forecasting methods and evaluation metrics for further analysis.

In addition, our package has a good extensibility and we also encourage users to implement more forecasting methods based on our framework and construct more energy forecasting scenarios. Similarly, after adding a new model, users can also add the new forecasting model to the forecasting task with a simple command.

We will provide a detailed description process in the Appendix C.2 and corresponding code documents to help users use our package to incorporate more power forecasting methods into the framework and construct different energy forecasting scenarios to adapt to different power forecasting tasks.

## 4 Benchmarking Process

### 4.1 Electricial Load Forecasting

In our archive, we will mainly discuss the results of probabilistic load forecasting. At the same time, to explain our proposed custom loss function, we will also compare the point forecasting performance of the forecasting model trained using gradient descent in Appendix D.5.

Data preprocessing and 24 hours-ahead forecasting. We first use the functions provided by our framework to fill in missing values and address the issue of zero padding. For forecasting scenarios, we chose the most common load forecasting, which is to forecast 24 hours in advance, as our main task for evaluation (our framework also supports the construction of other forecasting scenarios). To meet the needs of subsequent power grid scheduling, load forecasting needs to reserve sufficient time for subsequent tasks, which means that there is a certain gap between the available historical sequences and the forecasting range. Therefore, we adopt the widely used forecasting setting in the power industry (Qin et al., ; Wang et al., ), which uses the historical values of the previous 7 days at the same time to forecast the corresponding power load on the 8th day. And the basis for doing so mainly comes from the periodicity of the load sequence.

Feature engineering. As mentioned in Section 3.2, we apply the transformation of feature engineering based on temperature and calendar variables to our forecasting models. For sequence models like the LSTM, we concatenate the features with the output of the models and input them into a single-layer MLP to get the final output. As for the non-sequence models, we concatenate all the features and input lagged values. As a comparison, we also conduct experiments on non-transformed features simultaneously, directly inputting calendar variables and temperature as features.

Forecasting models. We introduce 21 probabilistic forecasting methods for comparison, covering multiple types. These include two simple moving quantile methods based on global historical data (BEQ) and fixed-length windows (BMQ), as well as two models that account for forecasting errors, based on the Persistence (BECP) and linear regression methods (CE).
In addition, there are 5 non-deep learning methods, and they are quantile regression methods based on the K-nearest neighbor algorithm (Hastie et al., ), quantile regression methods based on random forest and sample random forest (Meinshausen & Ridgeway, ), and quantile regression methods based on extreme random tree and sample extreme random tree (Geurts et al., ).
Finally, we introduce 12 deep learning methods, including simple forward propagation networks (Jain et al., ), LSTM networks (Hochreiter & Schmidhuber, ) for sequence modeling, convolutional neural networks (Li et al., ) (where we use one-dimensional convolutional kernels), Transformer (Vaswani et al., ) networks applying attention mechanisms. In addition, we include methods that modify the above neural network structures to make them more suitable for forecasting tasks, such as LSTNet (Lai et al., ), which is designed to simultaneously capture both long-term and short-term patterns of time series, WaveNet based on causal convolution (Oord et al., ), N-BEATS stacked into blocks using multiple linear layers (Oreshkin et al., ), Dlinear, and NLinear (Zeng et al., ) which only use linear layer to model the series. In recent years, Transformer-based forecasting models have been widely applied. Therefore, we also include related models like Informer (Zhou et al., ), Autoformer (Wu et al., ), and Fedformer (Zhou et al., ). Among them, the neural network is trained based on gradient descent. For probabilistic forecasting, we take the sum of nine quantile losses from 0.1 to 0.9 as the loss function. For point forecasting, we provide an asymmetric differentiable loss function through data fitting and integrate it into our forecasting framework as a loss function. At the same time, we also construct neural networks based on the traditional MSE Loss function for comparison(in Appendix D.5.2).

### 4.2 Renewable Energy Forecasting

Similar to load forecasting, we also conduct 24 hours ahead on renewable energy. Compared to the load, due to the great uncertainty of meteorological factors, the periodicity of the renewable energy sequence is weaker. Therefore, we will use a multi-step output approach for prediction. Firstly, we downsample the renewable energy data to a 1-hour resolution. Next, we will input the renewable energy sequence and meteorological factors of the past 24 hours, as well as the meteorological factors of the next 24 hours, into the model and predict the renewable energy sequence for the next 24 hours. Due to the need to adapt to multi-step and multi-quantile output formats. We mainly compared deep learning models based on gradient training (as described in the previous section).

## 5 Benchmark Evaluation

### 5.1 Benchmark on electrical load forecasting

With the help of the framework, we conduct extensive experiments on the collected load dataset based on the 21 probabilistic forecasting methods mentioned above. In addition, we also provide relevant point forecasting results for our proposed custom loss function. All experiments were conducted on Intel (R) Xeon (R) W-3335 CPU @ 3.40GHz and NVIDIA GeForce RTX3080Ti. Here, we primarily discuss the forecasting results of datasets with corresponding temperature data and the methods that can combine external data like temperature. The complete forecasting results as well as the running time of all the datasets are summarized in the Appendix D.4, and D.6, and our code repository.

#### 5.1.1 Comparison of temperature feature engineering

In this section, we will examine the impact of the temperature feature engineering we provide on the forecasting results from the perspective of probabilistic forecasting. We use Pinball Loss$\downarrow$ to evaluate the results of probabilistic forecasting and the point forecasting situation will be shown in the Appendix D.5.2.

Fig. 3 reports partial probabilistic and point forecasting results, where the blue part represents the results of incorporating our feature engineering, while the green one represents the results without doing so.
From the perspective of forecasting models, non-deep learning methods perform better than deep learning methods without the temperature transformation strategy. In deep learning methods, simple FFNN, LSTM, and CNN methods usually perform better than the more complicated ones. Moreover, complex deep learning models like the Wavenet and N-BEATS may even yield poor results. With the temperature transformation strategy, non-deep learning methods do not experience much improvement. The KNNR method experienced significant performance degradation on all datasets. This is because our feature engineering makes the input features very sparse, which seriously affects the performance of the K-nearest neighbor clustering algorithm, leading to a decrease in performance. However, deep learning methods have great improvements with this feature engineering.

Table 3(b) summarizes partial probabilistic forecasting results with PinBall Loss$\downarrow$ under different deep models. Among them, for the COVID-19 dataset, adding feature engineering significantly worsens the result. The characteristic of this data is that after the impact of COVID-19, the load of the power system has changed significantly, and there is a large deviation between the training set and the test set. Therefore, the decrease in forecasting performance indicates that after this feature engineering, the model tends to learn more about the relationship between temperature and load, while ignoring the influence of historical load to a certain extent. In addition, the probabilistic forecasting results in the Spanish dataset also indicate the negative effect of feature engineering. This is because the temperature data of this dataset does not exactly match the load data (see Appendix A for details). On the contrary, in datasets such as GEF12, 14, and 17, it can be seen that for relatively stable aggregated level loads with corresponding temperature data, such feature engineering can significantly improve the performance of the forecasting model.

We provide more complete results in the Appendix D.6 and the corresponding code repository.

Figure: Figure 3: Pinball Loss$\downarrow$ comparison on parts of datasets.
Refer to caption: x3.png

**Table 3: Comparison of temperature feature engineering on deep learning-based probabilistic forecasting by PinBall Loss (the lower the better, the underline indicates a performance improvement after incorporating temperature-calendar feature engineering).**
| Pinball Loss $\downarrow$ | FFNN | LSTM | CNN | Transformer | LSTNet | N-BEATS | WaveNet | Informer | DLinear | NLinear | Autoformer | Fedformer |
| --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- |
| Covid19 | 18214.73 | 17019.92 | 16278.11 | 18269.58 | 17422.96 | 19832.39 | 17772.31 | 21430.63 | 25734.82 | 25191.69 | 18859.62 | 19440.16 |
| GEF12 | 5373.21 | 7166.41 | 5360.12 | 7315.20 | 7218.51 | 7469.99 | 13636.99 | 7201.01 | 6662.59 | 6755.06 | 7054.01 | 7458.82 |
| GEF14 | 79.64 | 110.77 | 74.56 | 120.17 | 364.50 | 118.37 | 106.42 | 110.14 | 112.72 | 112.63 | 118.62 | 124.07 |
| GEF17 | 56.00 | 80.00 | 57.36 | 79.75 | 77.00 | 119.27 | 197.13 | 84.74 | 80.13 | 81.42 | 83.11 | 92.88 |
| Spain | 1270.92 | 1376.77 | 1195.38 | 1471.55 | 1338.71 | 1208.81 | 3008.40 | 1380.52 | 1678.62 | 1709.13 | 1473.18 | 1750.63 |
| Bull | 10.75 | 11.50 | 10.39 | 12.03 | 12.66 | 16.89 | 18.42 | 11.15 | 10.66 | 11.19 | 11.41 | 12.08 |
| Hog | 323.05 | 400.44 | 361.30 | 406.06 | 437.80 | 582.80 | 683.65 | 394.52 | 396.44 | 397.12 | 438.58 | 439.51 |
| Cockatoo | 41.33 | 56.91 | 30.14 | 57.50 | 54.91 | 103.64 | 103.50 | 42.38 | 45.58 | 44.92 | 46.56 | 49.56 |
| PDB | 346.48 | 552.28 | 319.26 | 531.36 | 532.23 | 1641.57 | 1645.28 | 574.33 | 572.37 | 573.86 | 569.49 | 620.75 |

### 5.2 Benchmark on renewable energy forecasting

Table 4 reports the performance of different models in renewable energy probabilistic forecasting tasks. Overall, the best-performing one is FFNN. The main reason for this is that the forecast of renewable energy is greatly influenced by external factors, namely meteorological factors. Relatively speaking, fitting time series is not as important. In this case, the advantage of advanced time series forecasting models in modeling sequence relationships is diluted. It can be seen that Fedformer and Autoformer perform poorly in most cases, while simpler Transformer and Informer models perform better compared to them. In terms of the types of renewable energy, even for the same type of renewable energy, there may be deviations in the performance of the model at different levels. Similarly, for renewable energy sources at the same level, no model can perform the best in all types. Constructing different forecasting strategies based on different situations is significant in providing a reference.

**Table 4: Pinball Loss $\downarrow$ comparison in the renewable energy forecasting.**
| Pinball Loss $\downarrow$ | Data | FFNN | LSTM | CNN | Transformer | LSTNet | WaveNet | NBEATS | Informer | DLinear | NLinear | Autoformer | Fedformer |
| --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- |
| 1 | LW_S | 2.65 | 2.74 | 2.86 | 2.77 | 2.80 | 3.58 | 3.47 | 2.76 | 2.94 | 2.96 | 3.55 | 3.37 |
| 2 | OW_S | 7.99 | 7.53 | 8.42 | 6.42 | 10.92 | 21.46 | 10.96 | 6.30 | 18.01 | 18.41 | 20.05 | 19.03 |
| 3 | PV_S | 0.55 | 0.73 | 0.67 | 0.56 | 1.05 | 1.79 | 1.01 | 0.51 | 0.82 | 1.09 | 1.36 | 1.18 |
| 4 | LW_C | 23.49 | 25.71 | 26.04 | 25.62 | 24.18 | 36.29 | 41.99 | 25.02 | 23.33 | 24.01 | 28.20 | 26.16 |
| 5 | OW_C | 27.63 | 31.85 | 30.54 | 36.84 | 30.96 | 48.64 | 57.84 | 35.06 | 27.91 | 28.53 | 33.14 | 31.20 |
| 6 | PV_C | 6.10 | 6.26 | 6.15 | 6.96 | 8.80 | 15.47 | 10.81 | 6.88 | 6.25 | 7.39 | 9.81 | 11.24 |
| 7 | LW_R | 60.04 | 64.02 | 63.28 | 65.04 | 65.38 | 99.82 | 118.34 | 63.41 | 62.82 | 63.47 | 76.78 | 70.27 |
| 8 | OW_R | 192.68 | 221.41 | 223.96 | 238.14 | 217.48 | 275.10 | 324.37 | 219.76 | 179.30 | 179.86 | 240.54 | 235.92 |
| 9 | PV_R | 25.07 | 26.26 | 26.20 | 27.66 | 46.80 | 70.51 | 62.87 | 40.09 | 25.27 | 25.52 | 39.28 | 64.21 |

## 6 Conclusions

In this paper, we construct a package and benchmark specifically designed for energy forecasting.
Our forecasting package comes with high accessibility and extensibility by dividing the entire energy forecasting process into several modules for users to freely combine and construct their forecasting frameworks. In addition, our package also provides the engineering implementation of features based on temperature and the construction method of custom loss functions by data fitting.
Meanwhile, with the help of our package, we have provided comprehensive forecasting benchmark results using multiple forecasting methods and multiple datasets as well as detailed discussion and analysis, serving as an important reference for researchers in the community.

## Appendix A Dataset Description

### A.1 Electrical Load Dataset

As shown in Table 1, we collected multiple load datasets at different levels and organized them into a user-friendly ‘.pkl’ format and users can obtain them through the URL at GitHub we provide. Now we will introduce their sources and detailed information one by one.

#### A.1.1 GEF12

The GEF12 dataset is sourced from the Global Energy Forecasting Competition 2012 (Hong et al., ). This competition has multiple tracks, and we have compiled the dataset provided by the load forecasting tracks as one of our benchmark datasets. In this dataset, there are a total of 20 aggregated-level load series data and 11 temperature series. It is worth noting that the one-to-one correspondence between these temperature data and load data has not been clearly defined. For simplicity, the strategy used in our benchmark testing is to simply use only one temperature series data as the temperature variable for all series (the randomly selected result here is the second temperature data). Each time series covers load data with a resolution of 1 hour from 0:00 on January 1, 2004, to 5:00 on June 30, 2008. Because this dataset is used for competitions and the integrity of the data is relatively good, we did not preprocess the data (such as filling in missing values).

Figure: Figure 4: Visualization of the first load series in the GEF12 dataset.
Refer to caption: x4.png

Figure: Figure 5: Visualization of the second load series in the GEF12 dataset.
Refer to caption: x5.png

Fig. 4 and Fig. 5 visualize the partial load sequence in the GEF12 dataset. From it, we can see that these load series have obvious periodicity and seasonality. And this is an important feature of aggregated-level load.

#### A.1.2 GEF14

The GEF14 dataset is from the Global Energy Forecasting Competition 2014 (Hong et al., ). This competition also has multiple tracks, and we focus on load forecasting tracks. The competition provides load data spanning up to 8 years from 2006 to 2014. Unlike the 2012 competition, We truncate the load data and only use the data from 2013 and 2014 for testing. On the one hand, it is because the impact of load data from many years ago on the current forecast is very small, and on the other hand, it is because most of the load data we collect is about 2 years in length. For relative consistency, we only took the last two years to construct our load forecast archive. Our adjusted load data covers load data with a resolution of 1 hour from 1:00 on January 1, 2013, to 0:00 on January 1, 2015. Fig. 6 shows the adjusted load data, similar to the data in GEF12, which is also aggregated level data. Therefore, the data in GEF14 also shows obvious periodicity and seasonality.

Figure: Figure 6: Visualization of the load series in the GEF14 dataset.
Refer to caption: x6.png

#### A.1.3 GEF17

The GEF17 dataset is from the Global Energy Forecasting Competition 2017 (Hong et al., ). Similar to GEF12, this dataset also provides multiple aggregated-level load data. However, the difference is that it clarifies the corresponding relationship between the temperature series and load series, providing a one-to-one temperature series corresponding to the load series. In terms of period, it provides load data from 2013 to 2017. For the reasons mentioned above, we have intercepted the load data and only used the load data from the past two years (i.e. 2016 and 2017). Finally, we used 8 aggregated-level load data from 2016 to 2017 and their corresponding temperature data to construct our load forecasting archive. Fig. 7 and Fig. 8 visualize some data in the GEF17 dataset, similarly, it also showcases the common characteristics of aggregated-level load data, namely periodicity and seasonality.

Figure: Figure 7: Visualization of the load series in the GEF17 SEMASS dataset.
Refer to caption: x7.png

Figure: Figure 8: Visualization of the load series in the GEF17 CT dataset.
Refer to caption: x8.png

#### A.1.4 Covid19

Figure: Figure 9: Visualization of the load series in the Covid19 dataset.
Refer to caption: x9.png

The Covid19 dataset is from the Day-Ahead Electricity Demand Forecasting Competition: Post-COVID Paradigm (Farrokhabadi et al., ). This dataset covers the load data from 0:00 on March 18, 2017, to 15:00 on November 5, 2020. In addition to load data and temperature, this dataset also provides other meteorological factors such as humidity and wind speed. To maintain the consistency of the forecasting archives, we did not consider such factors. Unlike the datasets mentioned above, this dataset focuses on the impact of COVID-19 on the power system. Fig. 9 shows the load data in the Covid19 dataset. The blue part indicates that the power system has not yet been impacted by COVID-19, similar to other aggregated level load data, showing periodicity. The orange section displays the load data after COVID-19. It can be seen that it is different from the blue part. The absolute value of the load rapidly decreases during the period being impacted, and then recovers smoothly after a period of time. However, compared to the same period when it was not impacted, the load value has decreased. This transformation poses a challenge to the robustness of forecasting models. As shown in the main text, our temperature-calendar variable feature engineering will make the model more inclined to remember the impact of the temperature and calendar variable on the load and ignore the historical value to a certain extent, which ultimately leads to the decline of the forecasting performance. Therefore, when encountering strong external events like this, day-ahead forecasting should focus more on historical values.

#### A.1.5 PDB

The PDB dataset is a public dataset from the Kaggle data competition platform (Yeafi, ). It covers load and temperature data from 1:00 on January 1, 2013, to 0:00 on January 1, 2015. Due to its moderate length, we did not intercept it. Fig. 10 shows its load data visualization results.

Figure: Figure 10: Visualization of the load series in the Spanish dataset.
Refer to caption: x10.png

#### A.1.6 Spanish

Figure: Figure 11: Visualization of the load series in the PDB dataset.
Refer to caption: x11.png

The Spanish dataset is also a public dataset from the Kaggle data competition platform (Nicholas, ). It provides nationwide load data for Spain from 0:00 on January 1, 2015, to 23:00 on December 31, 2018. At the same time, it also provides meteorological data (such as temperature and wind speed) corresponding to the five major cities in Spain. This situation is similar to GEF12. For the same reason, we select relevant meteorological data from the economically developed Barcelona region as the corresponding meteorological data. In addition, the load data of this dataset is partially missing (with a missing rate of 0.1$\%$). Because of the low missing rate, we used a simple Linear interpolation method to fill the data. Fig. 11 shows the corresponding load visualization results. Compared to other aggregated level loads, the periodicity and seasonality of the Spanish national load have become relatively less pronounced.

#### A.1.7 Hog

The Hog dataset comes from The Building Data Genome 2 (BDG2) Data-Set (Miller et al., ). BDG2 is an open dataset that includes building-level data collected from 3053 electricity meters, which covers 1636 buildings. From the perspective of the area where the building is located, it includes the load, cooling, and heating data of buildings in multiple areas such as Hog and Bull. From a period perspective, it covers data from 2016 and 2017. In addition, BDG2 also classifies buildings, including buildings for educational purposes, offices, and so on. Based on the characteristics of this dataset, we divide it by region, and the Hog dataset is composed of relevant load data from buildings in the Hog region in the BDG2 dataset. Because the data in this dataset is all building-level data and we often find situations such as missing values and outliers in data at this level(Jeong et al., ). Therefore, we first use the functions provided by the package to check for outliers. Specifically, we first calculate the lower quartile (Q1) and the upper quartile (Q3) and then calculate the quartile interval (IQR), that is, $IQR=Q3-Q1$. Here, the outlier is defined as the point that is lower than $Q1-q\times IQR$ or higher than $Q3+q\times IQR$. The outlier factor $q$ here is set to 1.5. We set the detected outlier as the missing value, and discard the sequence with a missing rate of more than 10$\%$. For sequences with a missing rate of less than $10\%$, we interpolate them (using linear, polynomial, etc.). Finally, we obtained 24 available load sequences and their corresponding temperature sequences for the Hog region.

Figure: Figure 12: Visualization of the Hog education Haywood in the Hog dataset.
Refer to caption: x12.png

Figure: Figure 13: Visualization of the Hog office Elizbeth series in the Hog dataset.
Refer to caption: x13.png

Fig. 12 and Fig. 13 show two load sequences in the Hog dataset. They belong to educational facilities and offices respectively. It can be seen that compared to aggregate-level datasets, building-level datasets exhibit greater uncertainty. The similarity of data for the same period in different years is also significantly lower than those aggregate-level ones. Although the data is only two years old, the building dataset also exhibits significant seasonality. Specifically, the load during summer and autumn is relatively high, while the load during winter and spring is relatively low. In addition, despite the different properties of buildings, they still maintain a relatively similar seasonality.

#### A.1.8 Bull

Similar to the Hog dataset, the Bull dataset also comes from the BDG2 dataset. Similarly, we screen and preprocess the building load data in the Bull area, resulting in 41 available sequences covering multiple building properties. Fig. 14 and Fig. 15 show the load data of two representative building types in the Bull area. Similar to other building-level load data, the manifestation of periodicity is not obvious. Meanwhile, sudden changes also occur from time to time, posing challenges for the forecasting model to accurately model and forecast.

Figure: Figure 14: Visualization of the Bull education Luke series in the Bull dataset.
Refer to caption: x14.png

Figure: Figure 15: Visualization of the Bull office Yvonne series in the Bull dataset.
Refer to caption: x15.png

#### A.1.9 Cockatoo

Cockatoo is also from the BDG2 dataset. However, after our screening, only one load sequence met our requirements, which is “Cockatoo Office Laila”. Fig. 16 shows the load characteristics of this building. It is worth noting that during the period from February to April 2016, the load data appeared relatively stable. This may be caused by the fault of the measuring meters, by human error in the reading, or it may be the real situation. These errors typically occur in building-level load data and it is difficult to avoid this situation through data cleansing unless we directly discard the relevant data. Here, we retain this data to evaluate its impact on the final forecasting performance.

Figure: Figure 16: Visualization of the Cockatoo office Laila series in the Cockatoo dataset.
Refer to caption: x16.png

#### A.1.10 ELF

The ELF dataset comes from the Kaggle data platform (Shahane, ). It is worth noting that the platform provides temperature data from multiple Panama cities as well as other meteorological data such as wind speed and humidity. However, the relationship between these meteorological data and load data has not been clarified, and unlike datasets such as Spanish, we are not clear about the detailed regions to which the load value data belongs. Therefore, we only conduct experiments on historical load series and calendar variables. However, users can also add relevant meteorological variables to the forecasting model through simple code. Fig. 17 shows the load data for this dataset. Similar to the Covid19 dataset, this dataset also shows the impact of COVID-19 on the power system (see the data after April 2020).

Figure: Figure 17: Visualization of the load series in the ELF dataset.
Refer to caption: x17.png

#### A.1.11 UCI

Figure: Figure 18: Visualization of the load series in the UCI dataset.
Refer to caption: x18.png

The UCI dataset is a power load dataset from the UCI database, which has been widely used in the field of machine learning(Dua & Graff, ). Unlike the original dataset, we select the processed version, which includes 321 load sequences with hourly resolution (Bergsma et al., ).

### A.2 Renewable Energy Dataset

We have collected and released a new dataset on renewable energy generation, which includes onshore wind power (represented by LW), offshore wind power (represented by OW), and photovoltaics (represented by PV). At the same time, we also divide the new energy output data into station level (represented by S), city level (represented by C), and regional level (represented by R) according to geographical scale from small to large. The data at the station level also includes real-time measurement of meteorological factors (see Table 5). Because these stations are all located in the same geographical zone, we have provided the relative distances between different stations in Fig. 19. People can download this dataset in our provided code repository.

Figure: Figure 19: Visualization of the distance between different station(KM).
Refer to caption: x19.png

**Table 5: Datasets in the load forecasting archive.**
| Data | Installed Capacity(MW) | External Variables |
| --- | --- | --- |
| PV_S_1 | 20 | Irradiance($W/m^{2}$), Temperature(°C) |
| PV_S_2 | 46 | Irradiance($W/m^{2}$), Temperature(°C) |
| PV_S_3 | 20 | Irradiance($W/m^{2}$), Temperature(°C) |
| PV_S_4 | 114.47 | Irradiance($W/m^{2}$), Temperature(°C) |
| PV_S_5 | 100 | Irradiance($W/m^{2}$), Temperature(°C) |
| PV_S_6 | 50 | Irradiance($W/m^{2}$), Temperature(°C) |
| PV_S_7 | 30 | Irradiance($W/m^{2}$), Temperature(°C) |
| PV_S_8 | 120 | Irradiance($W/m^{2}$), Temperature(°C) |
| PV_S_9 | 99 | Irradiance($W/m^{2}$), Temperature(°C) |
| PV_S_10 | 100 | Irradiance($W/m^{2}$), Temperature(°C) |
| LW_S_1 | 48 | Wind speed($m/s$), Wind direction |
| LW_S_2 | 40 | Wind speed($m/s$), Wind direction |
| LW_S_3 | 49.5 | Wind speed($m/s$), Wind direction |
| LW_S_4 | 36 | Wind speed($m/s$), Wind direction |
| LW_S_5 | 45.05 | Wind speed($m/s$), Wind direction |
| LW_S_6 | 66 | Wind speed($m/s$), Wind direction |
| LW_S_7 | 49.5 | Wind speed($m/s$), Wind direction |
| LW_S_8 | 49.5 | Wind speed($m/s$), Wind direction |
| LW_S_9 | 88 | Wind speed($m/s$), Wind direction |
| LW_S_10 | 39.5 | Wind speed($m/s$) |
| OW_S_1 | 300 | Wind speed($m/s$), Wind direction |

## Appendix B Feature Analysis for load forecasting

As we mentioned before, external features have a significant impact on load forecasting. And among them, temperature variables and calendar variables have the greatest impact. This is also recognized by the famous global energy forecasting competition. Based on this, the organizer developed a linear model called the HongTao vanilla model, which considers load, calendar variables, and temperature as the main variables and it serves as the benchmark for the forecasting competition (Hong et al., ). Therefore, in this section, we will visualize the relevant features in different levels of datasets to explore the relationship between load, calendar variables, and temperature. At the same time, we will also provide feature engineering based on the relationship among load, calendar variables, and temperature.

### B.1 Temperature-load analysis

Fig. 20 and 21 show scatter plots of the relationship between temperature load at two different levels, respectively. The data at the aggregated level comes from the GEF14 competition, while the building level is randomly selected from the BDG2(Bull) dataset. It is worth noting that in the BDG2 dataset, each building has its corresponding attribute usage, such as educational facilities, office space, and so on. We divided the scatter plot of load and temperature into 12 blocks by month, with the aim of expanding the relationship between temperature and load to the ternary relationship between temperature, load, and calendar variables. Among them, we can consider calendar variables as indicators of seasons, months, workdays (weekends), and hours, and explore the temperature load relationships of different seasons (months, etc.). Similarly, we will only analyze the months here and include the remaining analysis.

From figure 20, we can see that, in line with common sense, the relationship between load and temperature shows significant differences when in different months. From May to September, there is a significant positive correlation between load and temperature. Starting from October, the relationship between load and temperature gradually shifted from a significant positive correlation to an insignificant correlation. May to September is also a period of frequent high-temperature weather, indicating that when the temperature is high, there is a significant positive correlation between temperature and load. It is worth noting that this positive correlation is not always true. If we directly hand over the temperature variables to the model for modeling without processing, such changes in the relationship may cause confusion and ultimately lead to a decrease in forecasting performance.

When it comes to building-level load, figure 21 shows that the uncertainty of the detected load is significantly greater than that of the aggregated load since each month presents different temperature-load relationships. As this is the load data from educational facilities, there may be classifications such as teaching days or rest days, as shown in the figure, where there is a clear phenomenon of fragmentation within multiple months.

In summary, the two levels of load exhibit different load-temperature relationships in different months. This situation also occurs in other time scales such as the hour. Therefore, to make the forecasting model understand this relationship correctly, it is necessary to consider calendar variables and temperature together.

Figure: Figure 20: Temperature-aggregated load scatter plots for 12 months.
Refer to caption: x20.png

Figure: Figure 21: Temperature-building load scatter plots for 12 months.
Refer to caption: x21.png

## Appendix C Package usage

### C.1 Asymmetric loss function

The ultimate goal of energy is to minimize subsequent scheduling costs, which are closely related to prediction errors. Inspired by (Zhang et al., ), our package provides a piecewise linearized function to fit forecasting errors with other variables (which can be corresponding scheduling costs, etc.). At the same time, we also give an asymmetric loss function to replace the symmetric MSE loss function. Specifically, in actual power grid dispatch, the economic losses caused by forecasting values being less than the true values are often greater than the losses caused by forecasting values being greater than the true values. Therefore, we use a simple quadratic function to construct a piecewise generating function.

Here, $\epsilon$ represents the Forecasting Error Percentage(FEP) $\epsilon=\frac{f\left(x\right)-y}{y}$. We first use this function to sample and obtain many data points and then use a smoothing spline, denoted as $s$, to fit them. To avoid many breakpoints, which may make it difficult to integrate into our forecasting framework as a loss function (Perperoglou et al., ), we use piecewise linearization to approximate the spline function. The selection of breakpoints can be based on the following formula(Berjón et al., ),

$$ $\left\|s-L(\epsilon)\right\|_{2}\leq\frac{\left(\int_{\epsilon_{\min}}^{ \epsilon_{\max}}s^{\prime\prime}(\epsilon)^{\frac{2}{5}}d\epsilon\right)^{ \frac{5}{2}}}{\sqrt{120}K^{2}},$ $$

where K is the number of breakpoints, and the integration interval we choose here is $(-0.15,0.15)$. By controlling the error between piecewise linear functions and spline functions, we can obtain an appropriate number of breakpoints. Here we control the error lower than 0.005. As for the location of the breakpoints, we first calculate the cumulative breakpoint distribution function according to (De Boor & De Boor, ),

$$ $F\left(\epsilon_{k}\right)=\frac{\int_{\epsilon_{\min}}^{\epsilon_{k}}\left|s^ {\prime\prime}(x)\right|^{2/5}dx}{\int_{\epsilon_{\min}}^{\epsilon_{\min}} \left|s^{\prime\prime}(x)\right|^{2/5}dx}.$ $$

The breakpoints $\left\{\epsilon_{k}\right\}_{k=1}^{K-1}$will be placed such that each subinterval can contribute equally to the value of the cumulative breakpoint distribution function. To eventually integrate it into our forecasting framework as a loss function, we need to ensure that it is differentiable. And we can achieve this by inserting a quadratic function at each breakpoint. Specifically, we insert a quadratic function within the $0.000001$ distance before and after each breakpoint and obtain the parameters of the quadratic function by ensuring the continuity of the function and its first derivative at the two connections before and after.

As mentioned in section 3.2.2, we simulate an IEEE 30 bus system to get the relationship between forecasting error and dispatching cost. Fig. 22 shows the specific process of simulation. Here we mainly focus on two optimization problems DEAD and IPB. The mathematical definitions of these two can be found in (Zhang et al., ). To solve these two optimization problems, we have provided the corresponding [MATLAB code](https://anonymous.4open.science/r/ProEnFo-17CC). Based on the above process, we have provided the relationship between the load forecasting error and the actual dispatching cost for each hour within 24 hours of a day, as shown in Fig 23. All corresponding data is saved in the file "breakpoint$_$new.mat" we provide and we can construct the corresponding loss function through a single line of code. Note that our experiment is conducted based on hour 9.

Figure: Figure 22: Loss function generation (adapted from (Zhang et al., )).
Refer to caption: x22.png

Figure: Figure 23: Loss function visualization.
Refer to caption: x23.png

### C.2 How to use the package

#### C.2.1 how to construct forecasting scenarios

Our framework mainly constructs forecasting scenarios through the function "calculate$_$scenario". Below, we will introduce how to prepare the input for this function separately.

##### data

The load dataset, we need to input the data in Pandas format.

##### target

The column name of the variable we need to predict in the input Pandas, and other columns will be treated as external variables.

##### methods _ _ \_ _ to _ _ \_ _ train

A list containing the forecasting methods we need to include.

##### horizon

The default is to make predictions 24 hours in advance and multiple time scale forecasting can be made by adjusting this.

##### train _ _ \_ _ ratio

The division of the training and test sets.

##### feature _ _ \_ _ transformation

The strategies that used to stabilize the time series, including logarithmic transformation, data differentiation, and so on.

##### time _ _ \_ _ stationarization

The division of the training and test sets.

##### datetime _ _ \_ _ features

Define what calendar variables to consider, such as day, month, year, whether it is a holiday, etc.

##### target _ _ \_ _ lag _ _ \_ _ selection

Define how to select historical data for forecasting. In our default settings and our benchmark, we will select values from the same time point in the past seven days to forecast the corresponding values for the eighth day. In addition, we also provide a strategy for selecting highly correlated historical data based on the autocorrelation of the data.

##### external _ _ \_ _ feature _ _ \_ _ selection

We provide two strategies for selecting external variables: direct input and based on temperature calendar variable relationships.

##### post _ _ \_ _ processing _ _ \_ _ quantile

Quantile-based forecasting may sometimes result in lower quantiles being greater than higher quantiles, and the main focus here is to rearrange them.

##### post _ _ \_ _ processing _ _ \_ _ value

To limit the final output result, such as forcing the forecasting result not to exceed a certain value.

##### evaluation _ _ \_ _ metrics

We include various evaluation metrics for users to choose from, which can be referenced specifically from D.3.

#### C.2.2 How to add new models

Our framework mainly focuses on quantile-based probabilistic forecasting, and to add new models, we need to make definitions for the new models.

In addition, users also need to provide specific details of the forecasting model, that is, how to handle the output of the model and ultimately convert it into output. Note that here we must provide information on how the model handles external variables. Generally speaking, the specific form of external variables is related to the external variable processing strategy we define. Common external variables include meteorological factors such as temperature. If there are no external variables, we will use the corresponding calendar variables as external variables.

## Appendix D Benchmark Evaluation

### D.1 Electrical Load forecasting

In this section, we will introduce the hyperparameter settings in our load forecasting archive.
Table 6 shows the parameter settings for non-deep learning methods. Here, BMQ represents the moving quantity method based on a fixed number of past time points, while BEQ is based on all historical data. BECP represents that the forecasting error obtained by the persistence method on the training set is directly added to the forecasting results obtained by the persistence method as quantile forecasting. QCE is similar but replaces the persistence method with linear regression. In addition, there are quantile regression methods based on the K-nearest neighbor algorithm (Hastie et al., ), quantile regression methods based on random forest and sample random forest (Meinshausen & Ridgeway, ), and quantile regression methods based on extreme random tree and sample extreme random tree (Geurts et al., ).

**Table 6: Parameter settings for non-deep learning methods.**
| Method | Parameters |  |  |  |
| --- | --- | --- | --- | --- |
| Window size | N_neighbors | N_estimators | Quantiles |  |
| BMQ | 7 | - | - | 0.1$\sim$0.9 |
| BEQ | all | - | - | 0.1$\sim$0.9 |
| BCEP | - | - | - | 0.1$\sim$0.9 |
| CE | - | - | - | 0.1$\sim$0.9 |
| KNNR | - | 20 | - | 0.1$\sim$0.9 |
| RFR | - | - | 100 | 0.1$\sim$0.9 |
| SRFR | - | - | 100 | 0.1$\sim$0.9 |
| ERT | - | - | 100 | 0.1$\sim$0.9 |
| SERT | - | - | 100 | 0.1$\sim$0.9 |

Apart from those, we introduce several deep learning methods and they are feed-forward neural networks(FFNN) (Jain et al., ), LSTM networks (Hochreiter & Schmidhuber, ) for sequence modeling, convolutional neural networks (Li et al., ), and Transformer (Vaswani et al., ) networks applying attention mechanisms. Additionally, we also have methods that modify the above neural network structures to make them more suitable for time series forecasting, such as LSTNet (Lai et al., ), which is designed to simultaneously capture both long-term and short-term patterns of time series, Wavenet-based on causal convolution (Oord et al., ), N-BEATS stacked into blocks using multiple linear layers (Oreshkin et al., ), DLinear(NLinear) (Zeng et al., ) based on MLP layers, and some Transformer-based varients(Informer (Zhou et al., ), Autoformer (Wu et al., ), Fedformer (Zhou et al., )). Tables 7 and 8 respectively show the hyperparameter settings of the training process and the network structure and parameters of the relevant deep learning methods. We divide the entire dataset into training and test sets at a ratio of 0.2, and then divide the training set into the final training and validation sets at a ratio of 0.2. To reduce the impact of neural network overfitting, we enable the early stop mechanism. Specifically, when the loss on the validation does not decrease for 15 epochs, we will stop training.

**Table 7: Training process parameters.**
| Parameters |  |  |  |  |  |
| --- | --- | --- | --- | --- | --- |
| Loss function | Validation ratio | Epochs | Patience | Optimizer | Learning rate |
| MSE<br>custom loss function<br>PinballLoss(0.1$\sim$0.9) | 0.2 | 1000 | 15 | torch.Adam | 0.0005 |

**Table 8: Parameter settings for deep learning methods, here $e$ represents the dimension of the external variables.**
| Method | Parameters |  |  |
| --- | --- | --- | --- |
| Network structure | Network parameters | Quantiles |  |
| FFNN | 2 Dense layers | Dense_1: (1+$e$,50)<br>Dense_2: (50,quantiles) | 0.1$\sim$0.9 |
| LSTM | 1 LSTM layer<br>1 Dense layer | LSTM: (1,64,2)<br>Dense: (1+$e$,quantiles) | 0.1$\sim$0.9 |
| CNN | 3 Conv1d layers<br>2 Maxpool1d layers<br>2 Dense layers | Conv1d_1: (1+$e$, batch size, 1, 1)<br>Conv1d_2: (batch size, 128, 1, 3)<br>Conv1d_3: (128, 256, 1, 3)<br>Dense_1: (256, 128)<br>Dense_2: (128, quantiles) | 0.1$\sim$0.9 |
| Transformer | 1 Positional Encoding layer<br>1 Encoding layer<br>1 Decoding layer<br>1 Transformer layer<br>1 Dense layer | Encoding: (1, 256)<br>Decoding: (1, 256)<br>Transformer: d_model=256,<br>n_head=4,<br>dim_forward=512<br>Dense: (7$\times$256+$e$, quantiles) | 0.1$\sim$0.9 |
| LSTNet | 1 Conv2d layer<br>2 GRU layers<br>2 hidden(Dense) layers<br>1 Dense layer | Conv2d: (1, 16, 2, 1)<br>GRU_1: (16, 16)<br>GRU_2: (16, 32)<br>hidden_1: (16+1$\times$32, 1)<br>hidden_2: (16, 1)<br>Dense: (1+$e$, quantiles)<br>skip=1<br>highway=7 | 0.1$\sim$0.9 |
| WaveNet | 1 CausalConv1d layer<br>1 DilatedStack<br>2 Conv1d layers<br>1 Dense layer | CausalConv1d: (1, 16, 2, 1)<br>DilatedStack: residual size=16<br>skip size=4<br>dilation depth=2<br>Conv1d_1: (4, 1, 1, 0)<br>Conv1d_2: (1, 1, 1, 0)<br>Dense: (1+$e$, quantiles) | 0.1$\sim$0.9 |
| N-BEATS | 1 Trend Stack<br>1 Seasonal Stack<br>1 Dense layer | Trend: hidden=64,<br>theta dim=(4,8)<br>Seasonal: hidden=64,<br>theta dim=(4,8)<br>Dense: (1+$e$,quantiles) | 0.1$\sim$0.9 |
| DLinear, NLinear | 1 DLinear(NLinear) layer<br>1 Dense layer | Individual:False<br>pred_len = 1<br>Dense: (1+$e$,quantiles) | 0.1$\sim$0.9 |
| Informer,Autoformer, Fedformer | 1 Informer(Autoformer, Fedformer) layer<br>1 Dense layer | n_head=8<br>d_model = 16<br>pred_len = 1<br>Dense: (1+$e$,quantiles) | 0.1$\sim$0.9 |

### D.2 Renewable energy forecasting

In renewable energy forecasting, we conduct multi-step forecasting. For non-sequential models such as FFNN, we simply concatenate all the features into the model. For sequence models such as Informer, we have eliminated the operation of concatenating external factors with model-extracted features in load forecasting. Instead, we directly concatenate external information with the original sequence in the feature dimension and input all information into the model. Then, compared to the parameters in Table 8, we mainly changed the pred_length from 1 to 24. In addition, we adjust the total number of training epochs to 100 and the patience to 10.

### D.3 Evaluation metrics

To evaluate the forecasting performance of different methods in our set day-ahead forecasting, we will introduce many evaluation metrics, which are divided into metrics for point forecasting and metrics for probabilistic forecasting. It is worth noting that not all metrics are used to directly distinguish forecasting performance, and some of them may be used to describe the shape of probabilistic forecasting, thereby more comprehensively presenting the forecasting characteristics of different models. We will provide a detailed introduction below.

#### D.3.1 Point forecasting evaluation

Similar to (Godahewa et al., ), we adopt 4 metrics that are widely used to evaluate the results of deterministic forecasting, and they are MAPE (Mean Absolute Percentage Error), MASE (Mean Absolute Scaled Error), RMSE (Root Mean Squared Error), and MAE (Root Mean Squared Error) respectively. Their mathematical definitions are listed below, note that $\{y_{t}\}_{t=1}^{n}$ represents the actual value and $\{F_{t}\}_{t=1}^{n}$ represents the predicted one.

- •
MAPE. MAPE is a metric of forecasting accuracy that calculates the average percentage of forecasting error for all data points. The smaller the value of MAPE, the higher the forecasting accuracy. Due to its percentage error, it can be used to compare forecasting performance at different scales. However, MAPE may result in an infinite or very large error percentage for zero or near zero actual values. The formal definition of MAPE is given below
$\text{MAPE}=\frac{1}{n}\sum_{t=1}^{n}\left|\frac{y_{t}-F_{t}}{y_{t}}\right|
\times 100\%.$
- •
MASE. MASE is a scale-independent error measure that calculates errors by comparing the forecasting error with the average absolute first-order difference of the actual value sequence. The advantage of MASE is that it is not affected by the size of actual values, so it is more robust for forecasting problems of different sizes. The formal definition of MASE is given below
$\text{MASE}=\frac{1}{n}\sum_{t=1}^{n}\frac{\left|y_{t}-F_{t}\right|}{\frac{1}{
n-1}\sum_{t=1}^{n}\left|y_{t}-y_{t-1}\right|}.$
- •
RMSE. RMSE is a commonly used measure of forecasting error that calculates the square root of the average of the sum of squares of forecasting errors for all data points. The smaller the value of RMSE, the higher the forecasting accuracy. It is sensitive to outliers, which may lead to large forecasting errors. However, RMSE has good interpretability because its units are the same as the actual and predicted values. The formal definition of RMSE is given below
$\text{RMSE}=\sqrt{\frac{1}{n}\sum_{t=1}^{n}(y_{t}-F_{t})^{2}}.$
- •
MAE. MAE is also a commonly used measure of forecasting error, which calculates the average of the absolute value of forecasting errors for all data points. The smaller the value of MAE, the higher the forecasting accuracy. Compared with RMSE, MAE is less sensitive to outliers, so it may be more robust in the case of outliers. The formal definition of MAE is given below
$\text{MAE}=\frac{1}{n}\sum_{t=1}^{n}\left|y_{t}-F_{t}\right|.$

#### D.3.2 Probabilistic forecasting evaluation

Compared to point forecasting, probabilistic forecasting can provide more information. Therefore, we can evaluate the results of probabilistic forecasting from more aspects. We have summarized a total of 11 metrics to comprehensively evaluate the load probabilistic forecasting results and list them below. Note that we will also perform matrix visualization based on some metrics to help users better evaluate different prediction models.

- •
CoverageError (CE). CoverageError is a method of measuring the quality of forecasting intervals, which measures the difference between the proportion of actual observations falling within the forecasting interval and the expected coverage rate. A smaller CoverageError indicates that the forecasting interval captures actual observations more accurately. Here, $L_{t}$ and $U_{t}$ represent the lower and upper bound of the forecasting interval while $UB$ and $LB$ respectively represent the upper and lower bounds of the interval we want. It is worth noting that when we visualize it, we call it ReliabilityMatrix. Specifically, we first divide the quantiles into the upper half and the lower half with 0.5 as the boundary. And perform pairwise combinations to obtain different nominal coverage rates as the horizontal axis, while the vertical axis represents the actual coverage rate.
$CE=\frac{1}{n}\sum_{t=1}^{n}(I(L_{t}\leq y_{t}\leq U_{t})-(UB-LB)).$
- •
Winkler Score (WS). Winkler Score (WS) is a metric that measures the quality of forecasting intervals. The forecasting interval is the forecasting range for future observations, usually represented by a lower bound and an upper bound. Winkler Score is used to evaluate whether the forecasting interval accurately captures actual observations, taking into account the width of the interval. A lower Winkler Score indicates better forecasting interval quality. Here, the symbols used are the same as CE while $\delta=U_{t}-L_{t}$. Similar to CE, in the corresponding visualization matrix, the abscissa should be different nominal coverage rates, and for a central (1-$\alpha$ )$\%$ forecasting interval, it is defined as follows:
$WS_{a,t}=\begin{cases}\delta,&L_{t}\leq y_{t}\leq U_{t}.\\
\delta+\frac{2\left(y_{t}-U_{t}\right)}{\alpha},&y_{t}>U_{t}.\\
\delta+\frac{2\left(L_{t}-y_{t}\right)}{\alpha},&y_{t}<L_{t}.\end{cases}$
- •
Pinball Loss (PL). Pinball Loss considers the difference between the forecasting value and the actual observation value, and weights the error based on whether the forecasting value falls on the side of the actual observation value (above or below). This enables Pinball Loss to capture the uncertainty in probabilistic forecasting and assign different weights to symmetric errors in loss calculations. A lower Pinball Loss indicates a smaller error between probabilistic forecasting and actual observations. Here, $L_{\tau}$ represents the Pinball Loss at the quantile $\tau$ and $\hat{y}_{\tau,t}$ is the forecasting value of corresponding time and quantile. In our setting, we consider the sum of 99 quantiles from 0.1 to 0.9, and it is defined as follows:
$PL=\frac{1}{n_{\tau}\cdot n}\sum_{t=1}^{n}\sum_{i=1}^{n_{\tau}}L_{\tau}\left(
\hat{y}_{\tau,t},y_{t}\right).$
- •
RampScore (RS). RampScore measures the consistency of the slope (i.e. increasing or decreasing trend) between the forecasting sequence and the actual observation sequence. Firstly, we use the Swing Door compression algorithm (Khan et al., ) to compress the forecasting sequence and the observed sequence, and then calculate the first-order difference values of these two sequences separately. Finally, we calculate the absolute difference between the first-order difference values of the two sequences and take the average to obtain the RampScore. A lower RampScore indicates that the model is more capable of capturing trends in sequence changes. Here, we calculate RampScore for 9 quantiles from 0.1 to 0.9.
- •
CalibrationError. CalibrationError (Chung et al., ) mainly evaluates the accuracy of forecasting models in representing uncertainty. The CalibrationError represents the difference between the forecasting quantile and the actual quantile. A smaller CalibrationError means that the forecasting model has higher accuracy in representing uncertainty, while a larger calibration error means that the forecasting model has lower accuracy in representing uncertainty. In the visualization matrix, we show the proportion of the predicted value greater than the true value under different quantiles. The closer the forecasting method is to the line y=x, the better the performance will be.

In addition to the metrics mentioned above, we also provide many other metrics. Although we will not present each of them in detail here, interested users can easily visualize them with the open-source code we provide. These evaluation metrics include IntervalWidth, QuantileCrossing, BoundaryCrossing, Skewness, Kurtosis, and QuartileDispersion. Among them, IntervalWidth calculates the width of probabilistic forecasting intervals given by different methods while QuantileCrossing gives the ratio of any two quantiles in which the predicted value of the lower quantile is greater than the predicted value of the higher quantile. BoundaryCrossing calculates the probability that the true value falls outside the forecasting range. Skewness and Kurtosis are metrics that describe the shape of a probability distribution. As for QuartileDispersion, its detailed description can be found in (Bonett, ).

### D.4 Evaluation results

In this section, we will mainly demonstrate the forecasting performance of 14 out of the 16 probabilistic forecasting methods we mentioned earlier, as well as 7-point forecasting methods. The relevant results of the two methods based on moving average can be obtained in the [repository](https://anonymous.4open.science/r/ProEnFo-17CC) we provide.

#### D.4.1 Running time

In our archive, all experiments were conducted on Intel (R) Xeon (R) W-3335 CPU @ 3.40GHz and NVIDIA GeForce RTX3080Ti. Table 9(b), Table 10(b), Table 11, and Table 12 report the training and inference time of various methods separately (note that the time for calculating metrics is not included). From the perspective of deep learning probabilistic forecasting models, incorporating temperature feature engineering not only improves the forecasting performance (see section 5.1.1 and Appendix D.6) but also reduces the time spent. For non-deep learning models, incorporating temperature feature engineering greatly increases the required time. This may be because our feature engineering incorporates a large amount of sparse data, which is difficult for nondeep learning models to handle. From the perspective of point forecasting, the traditional MSE loss and asymmetric loss functions take approximately the same amount of time.

**Table 9: Comparison of running time for probabilistic forecasting (except for ELF and UCI).**
| Time(s) | FFNN | CNN | LSTM | LSTN | WaveNet | NBEATS | Transformer | Informer | DLinear | NLinear | Autoformer | Fedformer |
| --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- |
| With feature engineering | 671.898 | 668.576 | 1565.256 | 2396.536 | 3665.051 | 3915.706 | 6484.478 | 9994.297 | 2323.525 | 2339.317 | 12605.944 | 12043.003 |
| Without feature engineering | 1107.768 | 608.697 | 1738.488 | 4651.474 | 4717.738 | 5292.745 | 5097.802 | 6351.088 | 1478.425 | 1946.943 | 8697.968 | 11289.180 |

**Table 10: Comparison of running time for probabilistic forecasting (for ELF and UCI).**
| Time(s) | FFNN | CNN | LSTM | LSTN | WaveNet | NBEATS | Transformer | Informer | DLinear | NLinear | Autoformer | Fedformer |
| --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- |
| Without feature engineering | 6465.294 | 3187.983 | 8666.618 | 21245.597 | 23781.845 | 31907.537 | 26157.561 | 39601.757 | 5431.757 | 7583.688 | 46356.421 | 48030.044 |

**Table 11: Comparison of running time for different loss functions (except for ELF and UCI).**
| Time(s) | FFNN | CNN | LSTM | LSTN | WaveNet | NBEATS | Transformer |
| --- | --- | --- | --- | --- | --- | --- | --- |
| MSE without feature engineering | 1116.204 | 545.852 | 1145.189 | 2918.964 | 2869.714 | 4043.892 | 4475.159 |
| MSE with feature engineering | 585.367 | 716.680 | 1388.907 | 2369.784 | 3680.179 | 3606.027 | 6356.714 |
| custom loss function with feature engineering | 1275.936 | 1146.166 | 2106.425 | 3439.587 | 4615.414 | 4400.446 | 7074.623 |

**Table 12: Comparison of running time for different loss functions (for ELF and UCI).**
| Time(s) | FFNN | CNN | LSTM | LSTN | WaveNet | NBEATS | Transformer |
| --- | --- | --- | --- | --- | --- | --- | --- |
| MSE | 5216.818 | 2386.388 | 7562.510 | 15905.827 | 15621.957 | 23684.826 | 20519.646 |
| custom loss function | 12066.808 | 4485.706 | 14169.861 | 23616.685 | 21856.217 | 29754.229 | 23146.784 |

### D.5 Point forecasting results

#### D.5.1 Comparison of asymmetric fitting loss function

According to (Zhang et al., ), the relationship between load error and the scheduling cost it causes is not symmetrical; the cost of underestimating at peak and overestimating at low values is different. Therefore, we use the relationship between forecasting error and cost $(\epsilon_{i},C_{i})$ as shown in Fig. 2 to estimate the costs of different methods and evaluate the impact of our asymmetric loss function. Among them, because the results in the previous section show that temperature-based feature engineering significantly improves the deep learning network, we apply this feature engineering to all of the methods.

**Table 13: Comparison of different loss functions on deep learning-based forecasting by cost (see the description in section 3.2.2).**
| Cost$\downarrow$ | FFNN | LSTM | CNN | Transformer | LSTNet | N-BEATS | Wavenet |
| --- | --- | --- | --- | --- | --- | --- | --- |
| Covid19 | 0.2048 | 0.1965 | 0.2242 | 0.1749 | 0.2911 | 0.6141 | 0.6140 |
| GEF12 | 0.3259 | 0.3213 | 0.3337 | 0.3298 | 0.3314 | 0.3571 | 0.4438 |
| GEF14 | 0.0384 | 0.0417 | 0.0413 | 0.0421 | 0.0412 | 0.0648 | 0.0651 |
| GEF17 | 0.1109 | 0.0982 | 0.1028 | 0.1095 | 0.1075 | 0.1177 | 0.1614 |
| Spain | 0.2248 | 0.2192 | 0.2244 | 0.2025 | 0.2067 | 0.2038 | 0.2684 |
| Bull | 1.8616 | 1.7499 | 1.8071 | 1.7603 | 1.7768 | 1.7765 | 2.2614 |
| Hog | 1.4099 | 1.3431 | 1.2443 | 1.1334 | 1.3918 | 1.5175 | 1.8091 |
| Cockatoo | 1.7939 | 1.4710 | 1.8784 | 1.4991 | 1.3170 | 1.4124 | 1.7414 |
| PDB | 0.2487 | 0.1808 | 0.1848 | 0.1733 | 0.1906 | 0.1568 | 0.2412 |

Table 13(b) reports the point forecasting results with different loss functions.
It can be seen that asymmetric loss brings improvements in most cases.
However, in the GEF17 and Spain datasets, the performance of the asymmetric loss function is generally lower than that of the MSE loss function. Due to our estimation of dispatching costs coming from a 30-bus system, there may be some differences in scale compared to a real large-scale power system, so our loss function may perform poorly on aggregated-level data like the GEF17. On the contrary, with building-level data, our custom loss function can help the model recognize the asymmetric nature between forecasting error and cost, bringing significant positive effects on minimizing costs.
Overall, asymmetric loss functions can provide some performance improvement for most methods on most datasets, thereby minimizing the dispatching costs.

#### D.5.2 Comparison of point forecasting

Tables 14, 15, and 16 respectively report the performance comparison of forecasting models based on the MSE loss function in several different datasets. Among them, only a portion of the UCI dataset is reported. The remaining results and other evaluation metrics can be found in the Github(^2^22[https://anonymous.4open.science/r/ProEnFo-17CC](https://anonymous.4open.science/r/ProEnFo-17CC)). From the perspective of MAPE metrics, in some datasets, existing forecasting models cannot provide reasonable forecasting results. This situation is particularly severe in building-level datasets. Fortunately, in the aggregated level data, the vast majority of prediction models can provide reasonable prediction results. This is due to the stronger stationarity of aggregated level data compared to the building-level ones. From the perspective of predictive models, models based on simple structures perform better than relatively complex models such as LSTNet, WaveNet, and N-BEATS. Among them, WaveNet performs the worst in multiple datasets, indicating that it is not suitable for application in the scenario of day-ahead forecasting of the power grid.

**Table 14: MAPE comparison results based on MSE Loss function in datasets with temperature variables (I), the coloring indicates that the current forecasting model cannot obtain reasonable results on this dataset.**
|  | Methods |  |  |  |  |  |  |
| --- | --- | --- | --- | --- | --- | --- | --- |
| MAPE($\%$) | FFNN | LSTM | CNN | LSTNet | WaveNet | N-BEATS | Transformer |
| GEF12_1 | 9.21 | 9.53 | 9.16 | 9.32 | 9.77 | 11.72 | 9.73 |
| GEF12_2 | 5.58 | 5.41 | 5.53 | 5.57 | 6.99 | 7.03 | 5.49 |
| GEF12_3 | 5.65 | 5.48 | 5.59 | 5.46 | 7.02 | 6.97 | 5.51 |
| GEF12_4 | 27.36 | 26.28 | 26.67 | 26.8 | 28.4 | 27.04 | 26.96 |
| GEF12_5 | 9.3 | 9.19 | 9.17 | 9.47 | 11.05 | 9.52 | 9.62 |
| GEF12_6 | 5.53 | 5.48 | 5.65 | 5.53 | 6.89 | 5.57 | 5.57 |
| GEF12_7 | 5.65 | 5.48 | 5.59 | 5.46 | 7.02 | 6.97 | 5.51 |
| GEF12_8 | 8.48 | 8.08 | 8.12 | 8.45 | 11.16 | 8.16 | 8.02 |
| GEF12_9 | 118.52 | 107.78 | 124.49 | 122.08 | 141.54 | 140.2 | 99.19 |
| GEF12_10 | 26.21 | 31 | 27.29 | 26.81 | 39.4 | 33.38 | 31.74 |
| GEF12_11 | 6.87 | 6.97 | 6.8 | 7.01 | 10.08 | 6.6 | 6.38 |
| GEF12_12 | 7 | 7.17 | 6.46 | 7.36 | 10.03 | 7.16 | 7.1 |
| GEF12_13 | 7.67 | 7.67 | 7.96 | 7.76 | 9.67 | 7.74 | 7.88 |
| GEF12_14 | 10.57 | 10.97 | 10.58 | 10.86 | 13.29 | 11.77 | 11.08 |
| GEF12_15 | 8.32 | 8.21 | 8.29 | 8.36 | 10.61 | 8.88 | 8.67 |
| GEF12_16 | 8.46 | 8.88 | 8.54 | 8.84 | 10.26 | 9.82 | 9.21 |
| GEF12_17 | 7.26 | 7.01 | 6.79 | 7.08 | 9.76 | 7.41 | 7.01 |
| GEF12_18 | 8.21 | 8.27 | 8.14 | 8.28 | 10.35 | 8.67 | 8.32 |
| GEF12_19 | 9.85 | 9.83 | 9.56 | 9.78 | 12.66 | 10.16 | 9.98 |
| GEF12_20 | 6.8 | 6.44 | 6.64 | 6.39 | 9.36 | 9.4 | 6.5 |
| GEF14 | 2.26 | 2.48 | 2.52 | 2.36 | 3.02 | 3.01 | 2.44 |
| GEF17_CT | 3.76 | 3.75 | 3.63 | 3.78 | 5.02 | 4 | 4.29 |
| GEF17_ME | 3.08 | 3.08 | 3.21 | 3.07 | 3.2 | 3.22 | 3.06 |
| GEF17_NEMASSBOST | 3.85 | 3.51 | 3.84 | 3.69 | 4.67 | 3.98 | 3.88 |
| GEF17_NH | 3.67 | 3.7 | 3.66 | 4.76 | 4.61 | 3.85 | 3.93 |
| GEF17_RI | 3.3 | 3.38 | 3.34 | 3.36 | 4.4 | 4.48 | 3.55 |
| GEF17_SEMASS | 4.76 | 4.84 | 4.8 | 4.85 | 5.58 | 5.17 | 5.07 |
| GEF17_VT | 5.01 | 4.96 | 5.03 | 5 | 6.2 | 6.11 | 5.55 |
| GEF17_WCMASS | 3.88 | 3.8 | 3.8 | 3.7 | 4.86 | 4.11 | 3.92 |
| Covid19 | 5.19 | 5.16 | 5.55 | 6.54 | 11.71 | 11.71 | 4.9 |
| Bull_assembly_Amalia | 21.92 | 20.96 | 21.02 | 21.53 | 24.07 | 23.86 | 21.19 |
| Bull_assembly_Goldie | 20.34 | 18.49 | 18.27 | 23.78 | 24.79 | 19.79 | 20.12 |
| Bull_assembly_Lance | 26.42 | 26.6 | 26.43 | 27.26 | 27.7 | 27.38 | 27.27 |
| Bull_assembly_Maren | 86.83 | 60.86 | 89.6 | 111.2 | 124.41 | 77.62 | 53.34 |
| Bull_assembly_Nathanial | 12.12 | 11.67 | 10.88 | 11.19 | 12.68 | 13.67 | 11.42 |
| Bull_assembly_Vanessa | 21.52 | 20.12 | 21.34 | 19.4 | 23.4 | 19.88 | 20.22 |
| Bull_education_Arthur | 34.48 | 31.22 | 33.86 | 32.37 | 43.8 | 31.51 | 30.82 |
| Bull_education_Bernice | 21.11 | 20.13 | 20.99 | 20.24 | 26.06 | 22.33 | 20.74 |
| Bull_education_Brain | 29.12 | 27.88 | 27.5 | 29.52 | 31.81 | 30.89 | 27.53 |
| Bull_education_Brenda | 39.39 | 38.99 | 35.22 | 36.29 | 43.16 | 36.07 | 39.96 |
| Bull_education_Dan | 22.79 | 21.89 | 21.46 | 23.88 | 29.66 | 22.06 | 21.82 |
| Bull_education_Dora | 31.19 | 29.44 | 27.67 | 29.52 | 29.8 | 30.77 | 29.65 |
| Bull_education_Dottie | 33.92 | 21.67 | 30.64 | 20.37 | 45.3 | 19.78 | 19.93 |
| Bull_education_Elva | 37.03 | 39.39 | 34.51 | 38.09 | 43.1 | 40.74 | 36.9 |
| Bull_education_Gregory | 127.84 | 52.35 | 148.31 | 85.61 | 84.84 | 429.83 | 36.06 |
| Bull_education_Jae | 14.36 | 14.55 | 14.92 | 14.13 | 20.3 | 18.81 | 15.57 |
| Bull_education_Joseph | 22.98 | 24.14 | 23.07 | 23.27 | 31.05 | 24.45 | 26.48 |
| Bull_education_Kendra | 20.69 | 21.31 | 20.71 | 21.71 | 35.31 | 21.01 | 19.91 |
| Bull_education_Krista | 65.97 | 62.57 | 60.66 | 64.12 | 85.78 | 86.32 | 63.59 |
| Bull_education_Lenny | 34.1 | 29.99 | 31.52 | 31.53 | 47.07 | 30.59 | 30.66 |
| Bull_education_Linnie | 31.24 | 29.19 | 31.6 | 35.97 | 36.4 | 29.77 | 29.62 |

**Table 15: MAPE comparison results based on MSE Loss function in datasets with temperature variables (II), the coloring indicates that the current forecasting model cannot obtain reasonable results on this dataset.**
|  | Methods |  |  |  |  |  |  |
| --- | --- | --- | --- | --- | --- | --- | --- |
| MAPE($\%$) | FFNN | LSTM | CNN | LSTNet | WaveNet | N-BEATS | Transformer |
| Bull_education_Luke | 22.73 | 22.23 | 22.68 | 22.9 | 23.75 | 22.48 | 22.28 |
| Bull_education_Magaret | 137.24 | 108.8 | 138.13 | 109.81 | 190.14 | 109.75 | 109.24 |
| Bull_education_Mervin | 37.97 | 33.91 | 33.68 | 34.04 | 34.33 | 47.01 | 39.5 |
| Bull_education_Miranda | 21.35 | 20.8 | 21.38 | 21.35 | 25.53 | 20.36 | 20.38 |
| Bull_education_Patrina | 30.02 | 29.06 | 30 | 28.97 | 34.97 | 29.27 | 29.82 |
| Bull_education_Racheal | 56.78 | 65.22 | 58.9 | 64.3 | 55.41 | 67.88 | 68.37 |
| Bull_education_Reina | 30.04 | 28.79 | 25.6 | 26.84 | 32.55 | 28.89 | 26.5 |
| Bull_education_Roland | 378.96 | 418.75 | 356.88 | 375.33 | 381.9 | 478.95 | 467.88 |
| Bull_education_Roseann | 45.68 | 43.38 | 44.55 | 43.97 | 46.57 | 45.81 | 42.81 |
| Bull_lodging_Carie | 107.47 | 98.5 | 112.22 | 103.09 | 146.3 | 84.83 | 101.97 |
| Bull_lodging_Hugo | 55.93 | 53.55 | 70.42 | 59.78 | 65.33 | 55.56 | 61.16 |
| Bull_lodging_Jeremiah | 13.18 | 12.22 | 13.9 | 14.38 | 21.95 | 14.73 | 13.43 |
| Bull_lodging_Lettie | 25.19 | 25.65 | 24.71 | 24.7 | 29.7 | 25.14 | 27.04 |
| Bull_lodging_Melissa | 20.8 | 23.35 | 20.77 | 22.64 | 24.37 | 22.49 | 25.61 |
| Bull_lodging_Perry | 25.13 | 24.1 | 23.45 | 24.54 | 24.15 | 26.96 | 28.28 |
| Bull_lodging_Terence | 24.65 | 22.94 | 27.15 | 22.26 | 32.09 | 22.7 | 22.67 |
| Bull_lodging_Travis | 15.37 | 15.3 | 17.97 | 22.81 | 31.7 | 14.67 | 14.25 |
| Bull_office_Chantel | 78.93 | 68.54 | 79.19 | 73.42 | 96.3 | 68.03 | 68.57 |
| Bull_office_Yvonne | 29.72 | 23.36 | 34.93 | 25.04 | 48.33 | 22.8 | 23.29 |
| Bull_public_Hyun | 21.22 | 20.7 | 19.61 | 21.09 | 22.19 | 21.48 | 23.34 |
| Hog_education_Haywood | 16.44 | 16.58 | 12.97 | 20.98 | 19.96 | 16.01 | 17.47 |
| Hog_education_Jewel | 22.53 | 22.04 | 23.67 | 31.47 | 30.49 | 21.26 | 20.21 |
| Hog_education_Leandro | 120.88 | 234.08 | 211.67 | 182.21 | 241.2 | 192.01 | 208.26 |
| Hog_education_Luvenia | 21.09 | 18.64 | 16.75 | 18.26 | 20.38 | 18.4 | 19.44 |
| Hog_education_Sonia | 186.16 | 224.47 | 214.12 | 223.01 | 319.63 | 205.18 | 279.88 |
| Hog_lodging_Shanti | 106.31 | 81.9 | 101.81 | 88.25 | 161.28 | 158.54 | 75.17 |
| Hog_office_Betsy | 16.82 | 23.57 | 18 | 21.91 | 21.6 | 17.92 | 18.74 |
| Hog_office_Byron | 30.03 | 27.33 | 25.12 | 27.51 | 49.41 | 49.59 | 27.7 |
| Hog_office_Candi | 169.49 | 167.73 | 179.04 | 173.45 | 322.46 | 175.25 | 164.42 |
| Hog_office_Charla | 11.72 | 12.44 | 11.61 | 15.07 | 14.74 | 14.91 | 12.48 |
| Hog_office_Corey | 58.37 | 72.63 | 56.7 | 70.6 | 76.69 | 78.21 | 76.82 |
| Hog_office_Cornell | 22.39 | 19.87 | 25.47 | 20.82 | 20.82 | 30 | 19.27 |
| Hog_office_Elizbeth | 14.07 | 16.62 | 17.09 | 15.37 | 21.97 | 15.78 | 16.19 |
| Hog_office_Elnora | 16.56 | 16 | 14.87 | 19.62 | 19.73 | 20.04 | 16.43 |
| Hog_office_Leanne | 78.77 | 47.73 | 78.73 | 72.48 | 90.79 | 92.55 | 71.63 |
| Hog_office_Nia | 25.21 | 23.6 | 23.77 | 23.44 | 29.73 | 32.57 | 25.06 |
| Hog_office_Richelle | 100.88 | 73.73 | 96.45 | 74.58 | 74.87 | 63.8 | 74.83 |
| Hog_office_Roger | 8.05 | 7.67 | 8.47 | 6.76 | 13.67 | 9.05 | 7.54 |
| Hog_office_Rolando | 15.16 | 18.08 | 14.55 | 18.63 | 21.87 | 16.8 | 17.2 |
| Hog_office_Terry | 12.29 | 12.34 | 12.81 | 12.06 | 25.49 | 12.24 | 11.19 |
| Hog_public_Brad | 11.16 | 11.6 | 12.25 | 12.89 | 17.74 | 20.97 | 11.67 |
| Hog_public_Crystal | 450.2 | 280.87 | 426.0 | 274.48 | 291.52 | 592.19 | 281.31 |
| Hog_public_Kevin | 227.29 | 258.55 | 205.92 | 223.34 | 250.82 | 227.99 | 238.8 |
| Hog_public_Octavia | 19.01 | 17.12 | 19.94 | 17.96 | 21.71 | 21.15 | 17.67 |
| Cockatoo | 26.39 | 25.64 | 24.64 | 26.93 | 26.12 | 27.69 | 26.07 |
| PDB | 6.35 | 5.24 | 5.34 | 5.35 | 6.04 | 4.75 | 5.02 |
| Spain | 6.05 | 5.83 | 6 | 5.58 | 6.94 | 5.55 | 5.63 |

**Table 16: MAPE comparison results based on MSE Loss function in datasets without temperature variables (Partial), the coloring indicates that the current forecasting model cannot obtain reasonable results on this dataset.**
|  | Methods |  |  |  |  |  |  |
| --- | --- | --- | --- | --- | --- | --- | --- |
| MAPE($\%$) | FFNN | LSTM | CNN | LSTNet | WaveNet | N-BEATS | Transformer |
| ELF_load | 4.07 | 4.26 | 4.07 | 4.38 | 4.64 | 4.39 | 4.33 |
| UCI_0 | 85.93 | 83.47 | 68.91 | 96.69 | 90.91 | 253.25 | 89.81 |
| UCI_1 | 5.52 | 5.36 | 5.55 | 5.41 | 18.3 | 6.55 | 5.35 |
| UCI_2 | 27.81 | 15.52 | 19.39 | 186.33 | 176.38 | 17 | 18.9 |
| UCI_3 | 7.1 | 7.47 | 7.25 | 7.26 | 24.75 | 8.07 | 7.45 |
| UCI_4 | 11.5 | 11.89 | 12.08 | 38.47 | 26.91 | 12.76 | 11.74 |
| UCI_5 | 6.92 | 7.32 | 6.95 | 7.23 | 24.88 | 8.24 | 7.66 |
| UCI_6 | 36.15 | 36.92 | 40.84 | 32.01 | 72.87 | 42.37 | 41.16 |
| UCI_7 | 5.93 | 6.18 | 6.17 | 6.17 | 14.55 | 6.64 | 6.24 |
| UCI_8 | 13.51 | 14.43 | 13.6 | 15.03 | 14.94 | 28.97 | 14.02 |
| UCI_9 | 23.01 | 22.4 | 22.88 | 23.15 | 43.88 | 23.65 | 24.65 |
| UCI_10 | 8.27 | 8.48 | 8.47 | 8.51 | 21.28 | 21.32 | 8.37 |
| UCI_11 | 13.31 | 13.48 | 13.48 | 14.54 | 13.9 | 13.65 | 13.52 |
| UCI_12 | 7.64 | 7.8 | 8.07 | 23.69 | 24.29 | 10.35 | 8.31 |
| UCI_13 | 8.77 | 8.97 | 9.09 | 8.94 | 25.66 | 11.36 | 10.21 |
| UCI_14 | 5.66 | 5.91 | 5.56 | 5.81 | 16.87 | 6.35 | 5.84 |
| UCI_15 | 11.35 | 11.44 | 11.29 | 11.61 | 23.52 | 23.63 | 11.61 |
| UCI_16 | 9.79 | 9.46 | 9.81 | 19.3 | 19.31 | 10.08 | 9.47 |
| UCI_17 | 12.35 | 12.47 | 12.42 | 12.5 | 12.5 | 18.71 | 12.51 |
| UCI_18 | 5.83 | 6.02 | 5.86 | 6.09 | 21.16 | 6.66 | 6.54 |
| UCI_19 | 12.78 | 13.11 | 12.79 | 13.94 | 37.44 | 18.66 | 13.5 |
| UCI_20 | 11.32 | 11.76 | 11.23 | 26.05 | 12.08 | 26.04 | 11.79 |
| UCI_21 | 9.97 | 9.97 | 9.79 | 10.07 | 31.48 | 31.58 | 10.14 |
| UCI_22 | 8.57 | 8.58 | 8.33 | 8.67 | 25.04 | 9.25 | 8.79 |
| UCI_23 | 9.02 | 9.11 | 9.14 | 33.44 | 27.28 | 23.81 | 9.22 |
| UCI_24 | 7.25 | 7.42 | 7.56 | 7.53 | 20.24 | 20.14 | 7.51 |
| UCI_25 | 14.35 | 14.78 | 14.42 | 15.63 | 51.78 | 14.38 | 16.06 |
| UCI_26 | 6.91 | 7.12 | 6.88 | 21.33 | 21.34 | 7.76 | 7.13 |
| UCI_27 | 8.79 | 8.88 | 8.88 | 8.71 | 22.93 | 22.93 | 8.92 |
| UCI_28 | 9 | 8.97 | 9.16 | 9.18 | 18.54 | 9.94 | 9.48 |
| UCI_29 | 38.46 | 40.07 | 37.95 | 37.85 | 101.18 | 49.83 | 39.39 |
| UCI_30 | 9.94 | 9.72 | 10.49 | 10.04 | 23.17 | 10.78 | 9.9 |
| UCI_31 | 5.82 | 5.91 | 5.81 | 5.79 | 18.68 | 6.17 | 5.93 |
| UCI_32 | 41.21 | 67.08 | 23.04 | 113.77 | 362.63 | 47.05 | 89.83 |
| UCI_33 | 7.98 | 8.08 | 8.16 | 15.86 | 8.15 | 8.31 | 8.32 |
| UCI_34 | 6.71 | 6.95 | 6.77 | 7.21 | 7.03 | 7.76 | 7 |
| UCI_35 | 7.1 | 7.32 | 7.16 | 19.73 | 19.75 | 19.72 | 7.42 |
| UCI_36 | 21.95 | 22.8 | 22.11 | 23.09 | 63.61 | 24.83 | 22.97 |
| UCI_37 | 8.72 | 8.79 | 8.77 | 9.01 | 22.42 | 9.61 | 8.83 |
| UCI_38 | 9.34 | 9.71 | 9.48 | 9.73 | 21.71 | 10.4 | 9.65 |
| UCI_39 | 11.2 | 11.22 | 11.56 | 11.54 | 21.86 | 11.76 | 11.33 |
| UCI_40 | 6.28 | 6.45 | 6.25 | 6.56 | 18.2 | 18.1 | 6.62 |
| UCI_41 | 7.75 | 7.86 | 7.48 | 8.28 | 8.44 | 25.07 | 8.1 |
| UCI_42 | 10.75 | 10.8 | 10.9 | 10.87 | 22.12 | 11.37 | 10.84 |
| UCI_43 | 8.16 | 8.07 | 8.13 | 8.12 | 23.95 | 24.04 | 8.17 |
| UCI_44 | 8.98 | 9.37 | 8.88 | 9.24 | 31.73 | 9.78 | 9.32 |
| UCI_45 | 5.4 | 5.54 | 5.56 | 5.4 | 14.03 | 5.68 | 5.53 |
| UCI_46 | 15.56 | 13.98 | 15.85 | 36.2 | 33.63 | 43.56 | 13.57 |
| UCI_47 | 5.87 | 6.03 | 6.11 | 5.97 | 19.44 | 19.44 | 6 |
| UCI_48 | 14.31 | 14.27 | 14.88 | 33.58 | 33.63 | 33.55 | 14.52 |
| UCI_49 | 9.17 | 9.26 | 9.52 | 9.28 | 26.85 | 10.28 | 9.38 |
| UCI_50 | 12.4 | 12.56 | 12.41 | 13.39 | 22.93 | 15.04 | 12.58 |

### D.6 Probabilistic forecasting results

Similar to point forecasting, we present the PinballLoss results of the forecasting model on some datasets in Tables 17, 18, 19, and 20, while placing other results in [repository](https://anonymous.4open.science/r/ProEnFo-17CC) for users to read. These tables present the forecasting results for all datasets with temperature variables.

- •
From the perspective of forecasting models, the performance of non-deep learning forecasting models is generally better than that of deep learning models. Among them, ERT and SERT perform relatively well. In deep learning methods, LSTM networks perform well while Transformer, FFNN, and CNN can also achieve the best results on certain datasets. However, the relatively complex neural network methods are lagging on the vast majority of datasets, which is relatively consistent with the results of point forecasting methods.
- •
From the perspective of feature engineering, when there is no one-to-one correspondence between temperature data and load data in some data sets, such as Spain and GEF12, feature engineering on temperature may reduce forecasting accuracy. When there is a clear one-to-one correspondence, such as GEF17, feature engineering of temperature-calendar variables will greatly improve the model’s forecasting performance.

Fig. 24, 25 report the evaluation result of two data in the GEF17 dataset respectively as an illustration. These three figures show the results after adding the temperature-calendar variable feature engineering. In most of these figures, the abscissa represents different quantiles, and the ordinate is the corresponding evaluation metric. The best-performing methods on these three datasets by the PinballLoss metric are FFNN, CNN, and QCE. Take CT as an example. Although FFNN performs best in PinballLoss, it cannot maintain an advantage in other metrics. For WinklerScoreMatrix, the LSTM model performs better than the FFNN method under high nominal coverage. This indicates that the LSTM model is superior to FFNN in considering extreme scenarios. For RampScoreMatrix, although FFNN has achieved good results in low quantiles, it is not better than the simple QCE method in high quantiles. Similar results have also been observed in other datasets and it is rare for a forecasting method to overwhelm all other forecasting methods in all aspects. These examples show that different quantiles can be considered separately as well as different metrics are needed to focus on different aspects of the forecasting model so that different forecasting models can be distinguished. This is also a major contribution to our load forecasting archive.

Figure: Figure 24: Visualization evaluation metrics in the GEF17 CT dataset.
Refer to caption: x24.png

Figure: Figure 25: Visualization evaluation metrics in the GEF17 ME dataset.
Refer to caption: x25.png

Figure: Table 17: PinballLoss comparison without temperature-calendar variable feature engineering (I).

Figure: Table 18: PinballLoss comparison without temperature-calendar variable feature engineering (II).

Figure: Table 19: PinballLoss comparison with temperature-calendar variable feature engineering (I).

Figure: Table 20: PinballLoss comparison with temperature-calendar variable feature engineering (II).

## Appendix E Documentation

Long-term preserve plan: Currently, our relevant datasets and prediction results are saved in the folder in a cloud service. This is mainly because we are still updating it, and the main direction is to add more fine-grained data related to smart meters. After our dataset is fully developed, we will apply for the relevant DOI for it.

Author statement: We confirm that the relevant dataset sources comply with relevant regulations and we bear all responsibility in case of violation of rights, etc., and confirmation of the data license.

## Peer Reviews

confidence: 5
contribution: 3
initial_rating: 6
presentation: 3
questions: Please see the comments in the weakness section.
soundness: 3
strengths: An important and impactful application
- framework can handle feature engineering, missing data, common pre-processing pipeline
- can be customized to integrate more models and loss functions. designed to handle custom loss functions.
- multiple evaluation metrics, loss, public datasets have been integrated.
summary: This paper is on an open source package for benchamarking energy-forecasting time-series models. Authors have integrated a framework which can handle and showed performance of widely varying 11 energy public datasets and integrated both SOTA deep-learning models(mostly time-series forecasting) and other machine learning models. The framework can handle missing data features, and various forecasting metrics and point-based and probabilistic loss. The framework looks similar idea to open graph benchmark except on energy forecasting time-series.
weaknesses: The work is a very timely application for energy forecasting. The reviewer has some concern regarding the results and reproducability on new datasets and integration of new models.

Comments:
1. Are the model features been handled as extrinsic and intrinsic features?, e.g. Weekday vs weekend, #occupancies, static features

2. If the framework code can follow OOP format to aggregate all models, would be more convenient. I was looking for the base model in the repository, could not find it. Not sure if OoP is integrated. 
- E.g., MyQuantileModel can be used as base model and all different deep-learning models can inherit the base model, forward traning can be same as the base model or can be overridden. 
- Similar thing can be done for loss functions.

3. Is there a way to customize hyperparameters like time-lag, forecast horizon,  

4. A json config file for chyperparameters, and parameters selction for model and datasets would be a useful integration.
5. Can this repo be usable and reproducable for evaluating other time-series datasets? What things to be fixed for such cases?

6. pinball Loss are extremely high, which is unususal. 
	- How are the train and test datasets selected, is the test on unseen domain or same domain data?
	- how are the training data for building-level and aggregated level? Seems building level data is not multivariate forecasting. Are each building data trained through different models? Would be easier if authors can clarify on train and test dataset selection?
	- Can you show the forecasted vs ground-truth comaprison plot for different models?

7. The repository dont seem to include result plots for visualizing training and test forecast vs ground-truth comparison plots, etc. These plots are also missing in the paper. 

8. As the framework is designed to handle both building level and aggregated performance, showing some analysis plots on the performance of building vs aggregation would be convincing to understand the models capability.

---

confidence: 4
contribution: 2
initial_rating: 5
presentation: 2
questions: Please see the weakness section.
soundness: 2
strengths: In this work, the authors develop an open-source forecasting package for energy forecasting. They also provide a dataset that includes different types of renewable energy generation and key meteorological factors.
summary: The paper introduces an open-source package and benchmark specifically designed for energy forecasting, addressing the challenges in the field such as the significant impact of external factors like temperature and calendar variables. The authors highlight the differences between energy forecasting and general time series forecasting, emphasizing the need for targeted feature engineering and custom loss functions that consider the costs associated with power grid dispatch. They provide a modular framework that includes data pre-processing, feature engineering, forecasting methods, post-processing, and evaluation metrics. Additionally, they release a high-quality renewable energy dataset with corresponding meteorological data to facilitate research in this area.
weaknesses: While the authors claim that energy forecasting is different from general time series forecasting, they do not provide a detailed explanation of these differences. 

The feature engineering techniques are heavily centered on temperature and calendar variables. Other external factors that might affect energy forecasting, such as economic indicators or unexpected events (e.g., pandemics), are not extensively explored. 

The effectiveness of the proposed temperature-calendar feature engineering and custom loss function may be dataset-specific. 

The custom loss function is derived from simulations using an IEEE 30-bus test system. Due to its simplified structure and limited scale, this system may not fully capture the complexity and variability of real-world power networks, limiting the generalizability of the results to more complex systems.

While the feature engineering strategy is well-motivated, the paper could provide more justification for its generalization across different datasets and energy types. 

The overall clarity of writing needs improvement.

---

confidence: 4
contribution: 2
initial_rating: 3
presentation: 1
questions: ## Goal
What is the objective of this package?
- To gather datasets and models related to “Energy (load, wind, photovoltaic)”: However, why not include datasets like the London smart meter datasets or NREL Solar energy, or other existing datasets? Additionally, why not include the latest SOTA models? Such choices might prevent future users to select your repository.
- To provide datasets that have both energy and weather data: Yet why include datasets with only energy data (ELD and ELF)?
- To introduce a specific loss including subsequent tasks: But without real “tasks” in the provided dataset, how can this loss be effectively tested? Only the simulation is not sufficient.
Furthermore, what about the negative log-likelihood (NLL) loss that provides uncertainty (or probabilistic prediction) to point forecast models? Models using architecture based on diffusion or traditional models minimizing the NLL can provide the uncertainty of a given prediction (range of values for each time step). Isn’t this the solution to account for any potential task? Why would authors loss should be preferred compared to these options?

## Questions

### Q1
I agree that energy forecasting requires external factors' information to achieve better performance. However, having the information in the input is one thing; the model can learn the relation between the external factor and the target. But then comes the availability for the prediction window. Does the model have to predict the external factors along with the energy forecasting, or does the model receive the forecast of these external factors? For the latter, are the forecasts reliable? If there is a mistake in the forecast, there will be a mistake in the prediction. How is this taken into consideration in the metrics when it comes to benchmarking models?

> We will input the renewable energy sequence and meteorological factors of the past 24 hours, as well as the meteorological factors of the next 24 hours.

This sentence implies that your scenario has a perfect forecast of future weather, which is highly unfeasible in real scenarios. Noise should have been included to account for errors in weather forecasts, if authors want to provide a package that is close to real case scenarios.

### Q2
> […] ensuring minimum data information distortion.

Do we have a guarantee on that?


## Experiment Setting

### E1
> It is beneficial for the community to study the relationship between different renewable energy types under similar climate conditions in a certain area.

Why did the authors not demonstrate such an advantage in the experiment section, with a multivariate-to-multivariate scenario of different types of renewables?

### E2
There is a lack of information on the experiment setting. Are we in a Univariate-to-univariate, multivariate-to-univariate, or multivariate-to-multivariate forecasting? This uncertainty poses reproducibility problems. Figure 3 and Table 3 results would seem to be U2U, but it needs to be confirmed for each experiment. Authors should clearly state the forecasting task of each experiment.

### E3
What does MQ[model_name] stand for in Figure 3?

### E4
What does $123 456 789$ in Table 4's first columns stand for?

## Feature Engineering

### F1
How will feature engineering based on calendar info, such as workdays, benefit REF? Are there options for tailoring the engineered features according to the user's needs and datasets capabilities? Or it only provides the combination described in Section 3.2.1? Which definitely limit the usage of such features.

### F2
How are engineered features integrated into transformer-based models (informer, etc.)? Are they considered as channels? The usage of such features limits the forecasting task to univariate-to-univariate forecasting, which impedes the usage of this framework. Such features should have been better integrated into the models for easier usage.

## Dataset

### D1
> The UCI dataset is a power load dataset from the UCI database, […]

Not enough details. I assume it is the Electricity Load Diagrams, and it should be clearly mentioned.

### D2
Why is UCI ELD not in the results of Figure 3 and Table 3? It is the most well-known and used dataset in the TSF domain. Is it because there is no weather data to do your feature engineering? If yes, your platform enables feature engineering only with calendar data, right? If so, it should be possible to give results with and without feature engineering considering only calendar data.

### D3
Why, for the spatial aggregation of REF, are there no external variables? It would be better to provide min/max/std/mean of the weather indicators of the considered area to improve forecasts, as this is one of the main goals/advantages of this dataset/paper.

### D4
Where are these regions/cities located? This information might be necessary to perform some feature engineering for given tasks.

### D5
The Git repository does not contain all the datasets mentioned in the paper. Does ```load_with_weather.pkl``` correspond to REF?


## Limitations
The SOTA of transformer models is outdated. It would be more pertinent to have models like iTransformer or Pathformer.

## Revision
$L(\epsilon)$ is not introduced before being used.

## Citation Formatting Issue
 * “As stated in (Menezes et al., 2020), […]” should be “As stated by Menezes et al. (2020), […]”
 * “(Wang & Wu, 2017) discovered the asymmetry […]”
 * “Based on (Zhang et al., 2022), our framework […]”

Authors should use \citet{} as mentioned in ICLR guidelines:
> When the authors or the publication are included in the sentence, the citation should not be in parenthesis using \citet{} (as in “See Hinton et al. (2006) for more information.”). Otherwise, the citation should be in parenthesis using \citep{} (as in “Deep learning shows promise to make progress towards AI (Bengio & LeCun, 2007).”)

## Proofreading
- “[…] such as PV,onshore wind, and offshore wind, […]” missing space.
- “[…] on traditional machine learning models and deep learning models.” -> “[…] on traditional machine learning and deep learning models.”
- “[…] in data at this level(Jeong et al., 2021a).” missing space before citation.
- “[…] the real requirement(here we mainly refer […]”
- “To our knowledge, […]” -> To the best of our knowledge,
soundness: 2
strengths: - Provision of a Renewable Energy dataset.
- Important collection of Energy datasets.
- Combination of various techniques, models, and metrics.
summary: In this paper, the authors propose a novel renewable electricity generation dataset that can be utilized in various time series forecasting tasks. This dataset is embedded in a Git Repository, along with other datasets and accompanied by several pre- and post-processing techniques, metrics, point and probabilistic forecasting models, and feature engineering methods. 

The authors used this dataset to benchmark different forecasting models and test their proposed loss function, aiming to train models not only for optimal accuracy but also to comply with specific tasks.
weaknesses: - Clarity could be improved to better understand the different parts/experiments of the proposal.
- Paper flow could be enhanced.
- Lack of details for some parts, even with the appendix.
- Needs thorough proofreading.
- Missing important state-of-the-art (SOTA) Time Series Forecasting baselines.

---

confidence: 5
contribution: 1
initial_rating: 3
presentation: 2
questions: - Including datasets that mirror the complexity of operational power systems (such as the ARPA-E PERFORM dataset) would improve the package's utility and scalability. 
- It would be nice to compare with physical/hybrid models that have been well-established in the industry. The PERFORM dataset has predictions from the SAM model with outstanding accuracy for PV power prediction. 
- To align probabilistic forecasts with real-world power system applications, quantile regression should be supplemented with models that capture inter-series correlation, for example [1] and [2]. 
- Additionally, introducing metrics such as the energy score and variogram score would enable a more robust evaluation of probabilistic forecasts, especially in multivariate contexts.

[1] Hauzenberger, Niko, et al. "Gaussian process vector autoregressions and macroeconomic uncertainty." Journal of Business & Economic Statistics (2024): 1-17.
[2]Ashok, Arjun, et al. "TACTis-2: Better, faster, simpler attentional copulas for multivariate time series." arXiv preprint arXiv:2310.01327 (2023).
soundness: 2
strengths: - The paper's package allows for customized forecasting pipelines
- The paper collects many datasets and benchmarked on 21 models.
summary: The paper presents a modular package for energy forecasting that incorporates datasets for load, wind, and photovoltaic (PV) energy. It emphasizes feature engineering techniques and custom loss functions aimed at aligning forecasting error with real-world dispatch costs. The package supports both probabilistic and point forecasting methods, evaluated across a set of 21 machine learning and deep learning models. The authors highlight their contribution in terms of feature engineering for load data, customization of loss functions for energy forecasting, and extensive benchmarking results.
weaknesses: My major concern for this paper is as follows. 

- Scale of Datasets:
The datasets included are relatively small and may not fully capture the scale or complexity of modern power systems. For instance, the load and renewable datasets lack the granularity typically needed at operational levels, such as bus-level load forecasting or wind farm-level generation forecasts. A large-scale dataset, such as the ARPA-E PERFORM dataset provided by NREL, which includes hundreds of time series, would serve as a more relevant benchmark for modern power systems.
- External Factors:
The authors claim that external factors for renewables have not been extensively studied, which is inaccurate. In the energy sector, weather-informed and physics-informed models for load and renewable forecasting are well-established. Commercial software solutions, like Solcast and NREL's System Advisor Model (SAM), already rely heavily on meteorological factors and hybrid physical and machine learning models. These approaches are widely adopted and should be considered as relevant benchmarks. The software and their predictions are publicly available. 
- Limitations in Probabilistic Forecasting Approach:
While the paper extends to probabilistic forecasting through quantile regression, this approach does not account for the correlation between errors across time series, which is essential in multivariate contexts like interconnected power systems. Effective probabilistic forecasting for power systems often requires capturing correlation structures to support stochastic optimization of power dispatch problems (since the paper emphasizes power grid dispatch). The metrics used to evaluate probabilistic forecasts (e.g., Pinball Loss and CRPS) are limited. Important metrics like the energy score (I agree this is similar to CRPS, but I believe it is more general than CRPS), which handles multivariate distributions, and the variogram score, which captures spatial or temporal correlation, are missing. These metrics are commonly used in power systems to assess probabilistic forecasts and would improve the package's relevance for real-world applications.