# Interpretable Short-Term Electrical Load Forecasting Scheme Using Cubist

Authors: Jihoon Moon; Sungwoo Park; Seungmin Rho; Eenjun Hwang
Journal: Computational Intelligence and Neuroscience
 (2022)

---

## Abstract

Daily peak load forecasting (DPLF) and total daily load forecasting (TDLF) are essential for optimal power system operation from one day to one week later. This study develops a Cubist-based incremental learning model to perform accurate and interpretable DPLF and TDLF. To this end, we employ time-series cross-validation to effectively reflect recent electrical load trends and patterns when constructing the model. We also analyze variable importance to identify the most crucial factors in the Cubist model. In the experiments, we used two publicly available building datasets and three educational building cluster datasets. The results showed that the proposed model yielded averages of 7.77 and 10.06 in mean absolute percentage error and coefficient of variation of the root mean square error, respectively. We also confirmed that temperature and holiday information are significant external factors, and electrical loads one day and one week ago are significant internal factors.


---

## Body

## 1. Introduction

The increase in urban population has caused various problems, such as resource depletion, traffic congestion, and environmental pollution [1]. For the effective operation of complex urban systems, many municipalities and governments have been trying to transform existing cities into smart cities [2]. The smart city concept aims to improve the efficiency and security of urban infrastructure as well as the quality of life for its citizens [3]. For instance, smart cities can reduce GHG emissions by reducing traffic congestion and energy consumption and introducing technologies such as electric vehicles, energy storage systems (ESSs), and renewable energy (RE) [4]. In particular, improving the energy efficiency of buildings with ESSs and RE is an important issue for cities because building energy consumption is one of the main sources of GHG emissions [5]. Most smart city systems use recent technologies such as the Internet of Things and big data to implement various city services [6]. For instance, the energy efficiency of existing buildings can be improved using building energy management systems (BEMSs) [7].A BEMS is a computer-aided tool that improves energy efficiency between the grid operator and consumers through bidirectional interaction [7, 8]. It collects and analyzes data related to electrical energy consumption to establish operational plans for building energy use [8]. On the demand side, a BEMS provides ways for customers to reduce or shift peak energy consumption and trade the remaining energy [9]. On the supply side, it serves as a tool for optimal allocation of RE, ESS, and demand response in electric utility grids [10]. Here, short-term load forecasting (STLF) has been widely used to determine the amount of power necessary for the reliable operation of electric utility grids from the next hour to the next week [9, 11]. It includes daily peak load forecasting (DPLF), total daily load forecasting (TDLF), hourly electrical load forecasting, and very short-term load forecasting (VSTLF) [12]. DPLF and TDLF are used to predict from one day to one week later as an essential procedure for unit commitment, energy trading, security analysis, and tight scheduling of outages and fuel supplies in power systems [13, 14].Accurate STLF is challenging because typical electrical energy consumption has various patterns accompanied by uncertainties due to unforeseeable external factors [15]. Furthermore, when predicting electrical loads, it is necessary to adequately consider complex correlations between historical and current data [9]. Many studies have been conducted to achieve accurate STLF based on machine learning (ML) methods because they can properly extract implicit nonlinear relationships between input and output variables [16–19]. Table 1 briefly summarizes recent STLF models based on ML methods. For instance, Lee and Han [20] developed a day-ahead (DA) DPLF model using multiple linear regression (MLR). They collected daily peak electrical load data in South Korea from 2012 to 2016 through the Korea Power Exchange (KPX) and constructed an MLR model using the days of the week, seasons, average temperature, and historical loads from one day to one year before the prediction time point. The model achieved better predictive performance than the forecasting model of KPX, extreme learning machine, and autoregressive moving average (ARMA). Fan et al. [21] proposed a STLF model based on weighted K-nearest neighbor (KNN), called W-KNN. When constructing the KNN model, they considered the inverse of the Euclidean distance to give appropriate weights to each data point after selecting a value for K. In experiments using electrical load data from the National Electricity Market (Australia), their model outperformed KNN, ARMA, and artificial neural network (ANN) in predicting performance.Table 1
Summary of recent STLF studies based on ML techniques.
Author (Year)DatasetGranularityML methodRolling procedureModel interpretabilityLee and Han [20] (2017)South Korea provided by Korea Power Exchange (KPX)Daily peak loadMLRYesYesFan et al. [21] (2019)Australian Energy Market Operator (AEMO)8 hKNNNoNoDong et al. [22] (2021)Qingdao City in China1 hBaggingNoNoSun et al. [23] (2021)Tai’an City, Shandong Province in China1 hANNNoNoTruong et al. [24] (2021)Residential building with a renewable energy system1 hAANNNoNoFan et al. [25] (2020)New South Wales (NSW) in Australia30 minEMDSVRPSOAR-GARCHNoYesZhang et al. [26] (2020)Queensland (QLD) in Australia30 minVMDSRSVRCBCSYesNoBouktif et al. [27] (2018)Metropolitan France30 minGALSTM-RNNYesYesWang et al. [28] (2018)University campus in Florida1 hRFNoYesRuiz-Abellón et al. [29] (2018)University campus in Spain1 hBaggingRFCRFXGBNoYesAbbasi et al. [30] (2019)AEMO30 minXGBNoYesZhang et al. [31] (2020)More than 1,400 enterprises in Yangzhong High-Tech Zone, ChinaDailyK-means clusteringCatBoostNoYesOn the other hand, Dong et al. [22] developed an hourly electrical load forecasting model based on bootstrap aggregating (Bagging) for Chinese special days. They collected three years of hourly electrical load in Qingdao, China. To construct their STLF model, they defined a holiday variable, with 2 representing statutory days, 1 representing working days, and 0 representing bridging days, proximity days, and weekends. Their Bagging model showed better prediction performance than ANN and Bagging, which trained a holiday variable that included 1 for working days and 0 for holidays. Sun et al. [23] proposed an hourly electrical load forecasting model based on ANNs. They first collected hourly electrical load data of Tai’an City, Shandong Province, China, from 2016 to 2018. They then generated various input variables, considering timestamp, temperature, and historical electricity consumption to construct their STLF model. Their model showed a mean absolute percentage error (MAPE) of 4.11 and a mean absolute error of 33.88. Truong et al. [24] developed an additive ANN (AANN) model to predict building electrical energy consumption. They collected hourly electrical load data for one year from a residential building with an RE system and configured five input variables, such as days of the week, hours of the day, isolation, temperature, and historical electricity consumption, to train the AANN model. Their concept was based on a gradient boosting machine (GBM). Unlike GBM, which generally uses decision trees (DTs) as weak learners, the AANN trains iteratively by estimating an ANN as one weak learner and passing the remaining residuals back to other ANNs. The AANN model outperformed MLR, DT, ANN, and support vector regression (SVR) in predictive performance.Recently, various hybrid STLF models based on two or more ML techniques have been proposed to derive better prediction performance than single ML-based STLF models. For instance, Fan et al. [25] proposed an SVR-based STLF model, namely, EMD-SVR-PSO-AR-GARCH, by hybridizing with empirical mode decomposition (EMD), particle swarm optimization (PSO), and autoregressive-generalized autoregressive conditional heteroscedasticity (AR-GARCH). They firstly decomposed an original electrical load data sequence into conventional intrinsic mode functions (IMFs) and residuals using the EMD. Then, they used SVR-optimized PSO and AR-GARCH to fit IMF1 and other IMF components and residuals, respectively. Finally, they obtained the final prediction value by integrating and fitting the prediction values of the models. Their model outperformed ARMA, AR-GARCH, SVR, and others in predictive performance. Zhang et al. [26] developed a hybrid STLF model, called VMD-SR-SVRCBCS, using variational mode decomposition (VMD) and self-recurrent (SR)-SVR by optimizing the parameters through the cuckoo bird search process of the cuckoo search algorithm (CBCS). They performed data preprocessing using the VMD to obtain more accurate IMFs and applied SR-SVRCBCS to model each decomposed IMF for more accurate forecasting results. Here, the SR mechanism, inspired from the combination of Jordan’s and Elman’s recurrent neural network, was used to learn more recurrent information from the hidden layer values at the previous time point concerning the outcomes of the SVRCBCS models. Their model outperformed single ML models such as autoregressive integrated moving average (ARIMA), seasonal ARIMA, ANN, and SVR in predictive performance.However, because the decision-making process inside most of these models mentioned above is opaque (i.e., a black box), forecasting results derived from these models cannot be entirely accepted and utilized. Therefore, their interpretation has been another challenging task [32]. Recently, interpretability methods in ML have attracted increasing attention for constructing accurate and interpretable forecasting models [33, 34]. Here, “interpretable” means that the user can understand how the model employs the input variables to make predictions [33]. The variable importance (also referred to as feature importance) measure is the basis for enhancing the interpretability of a model [34]. Several studies have used the variable importance measure to confirm the most significant factors in STLF. Bouktif et al. [27] proposed a long short-term memory (LSTM) network-based STLF model using Metropolitan France’s electrical load data. They used a genetic algorithm to optimize the hyperparameter tuning of the LSTM model. They also confirmed that the historical load was the most significant input variable for model training using the variable importance of extra trees (ETs). Their LSTM model outperformed ridge regression, KNN, RF, GBM, ANN, and ET. Wang et al. [28] proposed an hourly electrical load forecasting model based on RF. They configured ten input variables representing weather, occupancy, and time-related data and constructed a forecasting model for each academic semester. They predicted the electrical load of two educational buildings and identified the most influential variables that differed by semester. The RF model performed better than the DT and SVR models.Ruiz-Abellón et al. [29] developed four 48-hour-ahead electrical load forecasting models using hourly electrical load data collected from the Technical University of Cartagena in Spain. They used Bagging, RF, conditional RF (CRF), and XGB as tree-based ensemble methods to construct the forecasting models. They also described the essential input variables for training each forecasting model through variable importance. In the experiments, the XGB model achieved better prediction performance than Bagging, RF, and CRF. Abbasi et al. [30] proposed a 30-minute-interval electrical load forecasting model based on XGB. They used variable importance to extract input variables from the historical load during a week and confirmed that the historical loads close to the prediction time point and from a week before the prediction time point had high importance for the model construction. Subsequently, they constructed the XGB-based forecasting model using the extracted input variables. The XGB model exhibited a MAPE of 10% and an accuracy of 97%. Zhang et al. [31] proposed a TDLF model based on K-means clustering and categorical boosting (CatBoost). They collected total daily electrical load data from Yangzhong High-Tech Zone in China and executed K-means clustering to group industrial customers with similar load features into the same cluster. They then constructed a CatBoost-based TDLF model for each cluster and showed the variable importance of the CatBoost model. Their proposed model outperformed ARMA, LSTM, GBM, and CatBoost in predictive performance.Despite these efforts, there are still limitations. For instance, when a forecasting model based on these methods is constructed using conventional time-series forecasting model evaluation [35], we can determine the essential input variables by analyzing the variable importance in the model constructed from the training dataset. However, conventional time-series forecasting model evaluation performs unsatisfactorily when there is a significant gap between the training set period and the test set period [14]. This makes it difficult to ensure confidence in the decision-making process of the model. In addition, most of the studies were conducted mainly on STLF with high time resolution, such as hourly or subhourly intervals, for DA energy planning. Therefore, further studies are needed on the quantitative DPLF and TDLF for DA and weak-ahead (WA) energy planning. Furthermore, although a Cubist regression model has shown excellent performance in time-series forecasting [36–38], its use for STLF has rarely been reported [38].To address these issues, this study proposes a robust interpretable short-term electrical load forecasting model for accurate DPLF and TDLF. To this end, we collected five electrical load datasets from two commercial buildings and three educational building clusters. We configured various input variables that highly correlate with DPLF and TDLF. Then, we constructed STLF models using Cubist and time-series cross-validation (TSCV) to achieve high accuracy. The main contributions of this paper are as follows:(1)
We use electrical load data from two public buildings and three building clusters with different purposes to predict building-level electrical energy consumption, which has more complex patterns.(2)
We configure different input variables to predict both DPLF and TDLF by considering both DA and WA energy planning for effective BEMS operation.(3)
When a forecasting model is constructed for WA forecasting, we obtain multistep-ahead forecasting of all prediction time points (from one day to seven days) to compensate for uncertainty.(4)
To address the data shortage problem and reflect current electrical load trends and patterns, we use the cross-validation procedure based on a rolling forecasting origin.(5)
We compare our proposed model with other popular statistical and machine learning methods in terms of four forecasting types: DA-DPLF, DA-TDLF, WA-DPLF, and WA-TDLF.(6)
We perform an in-depth analysis of which input variables are the most important factors in electrical load forecasting for each dataset by using the variable importance of the proposed model.The remainder of this paper is organized as follows. Section2 describes the data preprocessing used to configure different input variables for forecasting types and presents the process of constructing the proposed model. In Section 3, we analyze the experimental results to demonstrate the superiority of the proposed model. Finally, we conclude our study and present the directions for future research in Section 4.


## 2. Materials and Methods

In this section, we describe in detail the data preprocessing and forecasting model construction. Figure1 illustrates the overall flowchart of the proposed method.Figure 1
Architecture of interpretable short-term electrical load forecasting model (DA: day-ahead, WA: week-ahead, DPLF: daily peak load forecasting, and TDLF: total daily load forecasting).


### 2.1. Data Collection and Preprocessing

In this section, we first present a data resampling process to construct data suitable for forecasting purposes. Then, we describe the process of configuring the input variables according to the purpose (i.e., DA and WA forecasting) using timestamps, weather, and historical load information. Figure2 illustrates the framework of the data preprocessing for constructing a Cubist model.Figure 2
Data preprocessing for Cubist modeling (DA: day-ahead, WA: week-ahead, DPLF: daily peak load forecasting, and TDLF: total daily load forecasting).We collected electrical load datasets from five different types of buildings or building clusters, as summarized in Table2. We first collected publicly available datasets from two buildings in Richland, WA, USA [39, 40]. The datasets consist of three years’ worth of information, including the hourly electrical load, the hourly outdoor temperature, and the corresponding timestamps. We filled in the missing values (i.e., daylight saving time in North America) in both datasets using linear interpolation. Because the dataset contains the hourly electrical load data for a day (24 rows), we calculated the maximum load value and the sum of all load values for each day for DPLF and TDLF, respectively.Table 2
Building information.
Dataset #Number of buildingsBuilding type (description)LocationDataset periodPublic accessBuilding 11Commercial (office)Richland, WashingtonJan. 2, 2009–Dec. 31, 2011YesBuilding 21Commercial (office)Richland, WashingtonJan. 2, 2009–Dec. 31, 2011YesCluster 116Educational (dormitory)Seoul, South KoreaJan. 1, 2016–Dec. 31, 2018NoCluster 232Educational (humanities bldg.)Seoul, South KoreaJan. 1, 2016–Dec. 31, 2018NoCluster 35Educational (engineering bldg.)Seoul, South KoreaJan. 1, 2016–Dec. 31, 2018NoIn addition, we collected typical 15-minute-interval electrical load datasets from three clusters of buildings at a private university in Seoul, South Korea [41]. The dataset collection period was three years. The first cluster comprised 16 residential buildings, and the electrical load had a residential pattern. The second cluster consisted of 32 academic buildings, including the main hall, library, classrooms, and offices. The third cluster contained five science and engineering buildings. This cluster exhibited much higher electrical loads per building than the other clusters, mainly due to the various experimental equipment and devices in the laboratories. In South Korea, the daily peak load of the building was calculated by multiplying the highest value among the 15-minute-interval electrical load used per day by 4. We took this into account when calculating the daily peak electrical load for DPLF. Likewise, we took the sum of 96 values per day to perform TDLF.Tables3 and 4 provide some statistics on daily peak loads and total daily loads for the five datasets, respectively. Weather conditions and holiday information are closely related to electrical load [42]. Therefore, we used these data to configure the input variables for Cubist modeling. Because we performed two types of forecasting, namely, DA and WA forecasting, we used different input variable configurations for each type of forecasting.Table 3
Statistics on daily peak electrical load data (unit: kW).
StatisticsBuilding 1Building 2Cluster 1Cluster 2Cluster 3Number of valid cases10941094109610961096Mean49.1654.461575.944132.022606.01Standard deviation21.6921.15308.911327.07451.57Trimmed mean50.4056.231552.324325.762623.04Median48.5954.521561.974176.242670.00Median absolute deviation19.8718.52321.311537.16437.66Minimum8.8610.97878.401426.561579.20Maximum141.11135.002623.686900.483549.60Range132.25124.031745.285473.921970.40Skew0.340.050.42–0.25–0.31Kurtosis0.430.09–0.15–1.03–0.71Standard error0.660.649.3340.0913.66Table 4
Statistics on total daily electrical load data (unit: kW).
StatisticsBuilding 1Building 2Cluster 1Cluster 2Cluster 3Number of valid cases10941094109610961096Mean719.31852.6329802.4162563.1049440.52Standard deviation250.53290.955350.3917167.956409.23Trimmed mean723.27850.1529390.7662514.4849506.54Median714.88844.9729569.9262872.2349771.80Median absolute deviation229.40244.275407.2820269.876802.47Minimum198.22242.0919013.0427961.4432546.40Maximum1527.302130.5049235.7698475.8464403.70Range1329.081888.4130222.7270514.4031857.30Skew0.160.440.45–0.13–0.11Kurtosis–0.320.98–0.03–1.06–0.64Standard error7.578.80161.61518.58193.86Time is another important factor for electrical loads. We considered various temporal variables, such as months, days, and days of the week. However, these variables cannot reflect periodic information when applied to forecasting models. For instance, 31 December and 1 January are temporally contiguous, yet, the range of both values in sequence form is 30 and 11 as day and month, respectively. Hence, we represented them as continuous data in the two-dimensional (2D) space to reflect their periodicity using equations (1)–(6) [41, 42]. Here, the variable for seven days of the week is defined as 1 to 7 from Monday to Sunday, according to international standard ISO 8601. LDMmonth represents the last day of the month to which the day belongs (e.g., January: 31, February: 28 or 29, March: 31, and so on). Consequently, we used six input variables to describe the date and time of the prediction time points.(1)Monthx=sin36012×Month,(2)Monthy=cos36012×Month,(3)Dayx=sin360LDMmonth×Day,(4)Dayy=cos360LDMmonth×Day,(5)Day_of_the_Weekx=sin3607×Day_of_the_week,(6)Day_of_the_Weeky=cos3607×Day_of_the_week.To verify the validity and applicability of the 2D representation, we computed several regression statistics on the electrical loads in one-dimensional (1D) space, consisting of three temporal variables (i.e., months, days, and days of the week), and 2D space, as shown in Tables5 and 6. Here, the residual standard error (RSE), multiple R-squared, and adjusted R-squared were calculated using equations (7)–(9).(7)RSE=1n−2∑t=1nyt−yt^2,(8)R2=1−∑yt−yt^2∑yt−y¯2,(9)Adjusted_R2=1−1−R2n−1n−p−1,where yt and yt^ represent the actual and estimated values, respectively, at the time t. y¯ represents the mean of the actual values, and n and p indicate the number of observations and variables, respectively. From the tables, we can observe that the 2D representation exhibits correlations more effectively than the 1D representation.Table 5
Residual standard error and R-squared statistics for daily peak electrical load data.
StatisticsBuilding 1Building 2Cluster 1Cluster 2Cluster 31D2D1D2D1D2D1D2D1D2DResidual standard error (unit: kW)17.717.117.416.5285.8271.411281121365.9351.9Multiple R-squared (unit: %)33.938.632.639.714.623.227.929.136.441.3Adjusted R-squared (unit: %)33.738.232.439.414.422.827.728.736.241.0Table 6
Residual standard error and R-squared statistics for total daily electrical load data.
StatisticsBuilding 1Building 2Cluster 1Cluster 2Cluster 31D2D1D2D1D2D1D2D1D2DResidual standard error (unit: kW)21120325622849584698146801467056995417Multiple R-squared (unit: %)29.334.922.939.014.423.327.127.427.434.6Adjusted R-squared (unit: %)29.134.622.738.714.122.926.927.027.234.2Typical electrical load patterns are different on weekdays and holidays, depending on the type of building [16, 40]. To reflect this, holiday information was also used as an input variable of the forecasting model. We collected holiday information by country from https://www.timeanddate.com [43]. The holidays included Saturdays, Sundays, and public holidays. We used one-hot encoding (i.e., “1” for the relevant data and “0” otherwise) as a nominal scale. Hence, we used seven time factors as input variables at the prediction time point.In general, the electrical load increases in the summer and winter due to the heavy use of air-conditioning and electrical heating appliances, respectively [13, 44]. Here, the most influential factor in the electrical load is temperature [13]. In this study, we focused on temperature-related variables, such as the daily maximum, average, and minimum temperatures [44]. In the case of Buildings 1 and 2, because the datasets we collected also included the outdoor temperature along with the electrical load, we calculated the daily maximum, average, and minimum temperatures (in the Fahrenheit scale). For Clusters 1 to 3, we first collected hourly temperature data using regional synoptic meteorological data provided by the Korea Meteorological Administration (KMA) and collected by the Seoul Meteorological Observatory located about 6 km from the university campus. As in Buildings 1 and 2, we calculated the daily maximum, average, and minimum temperature in Celsius. In South Korea, the KMA provides weather forecasts, including short and mid-term forecasts, for the most significant regions. The KMA mid-term forecast service provides the maximum and minimum temperatures from day 3 to day 10, as shown in Figure 3.Figure 3
Example of mid-term forecast provided by the KMA.Historical electrical load data are particularly important in electrical load forecasting because they exhibit the trend of recent electrical loads [17, 45]. To reflect the recent trend in the prediction, we used the daily peak and total daily loads of the previous seven days as input variables for the DPLF and TDLF models, respectively. Because the electrical load trends of weekdays and holidays can differ, we added a holiday indicator to indicate whether the day was a holiday [46]. We used a total of 24 input variables to build the DA-DPLF and TDLF models.Now, we describe the configuration of the input variables for the WA-DPLF and TDLF models. First, we used the same time and temperature variables as in the DA forecasting models. However, the daily peak and total daily loads from the previous day to the sixth day before the forecast date are unknown when considering WA forecasting. To compensate for this, we configured the input variables using the daily peak and total daily loads and holiday indicators of the same day of the previous four weeks for the WA-DPLF and TDLF, respectively [46]. Therefore, we used a total of 18 input variables to construct our WA electrical load forecasting models. Table 7 presents all the input variables that we considered for DA and WA forecasts.Table 7
Input variables for day-ahead and week-ahead forecasts.
IV #Input variables for day-ahead forecastingInput variables for week-ahead forecastingVariable type01MonthxMonthxContinuous [–1, 1]02MonthyMonthyContinuous [–1, 1]03DayxDayxContinuous [–1, 1]04DayyDayyContinuous [–1, 1]05Day of the weekxDay of the weekxContinuous [–1, 1]06Day of the weekyDay of the weekyContinuous [–1, 1]07HolidayHolidayBinary08Minimum temperatureMinimum temperatureContinuous09Average temperatureAverage temperatureContinuous10Maximum temperatureMaximum temperatureContinuous11Holiday (the day before seven days)Holiday (the day before four weeks)Binary12Electrical load (the day before seven days)Electrical load (the day before four weeks)Continuous13Holiday (the day before six days days)Holiday (the day before six three weeks)Binary14Electrical load (the day before six days)Electrical load (the day before three weeks)Continuous15Holiday (the day before five days)Holiday (the day before five two weeks)Binary16Electrical load (the day before five days)Electrical load (the day before two weeks)Continuous17Holiday (the day before four days)Holiday (the day before one week)Binary18Electrical load (the day before four days)Electrical load (the day before one week)Continuous19Holiday (the day before three days)—Binary20Electrical load (the day before three days)—Continuous21Holiday (the day before two days)—Binary21Electrical load (the day before two days)—Continuous23Holiday (the day before one day)—Binary24Electrical load (the day before one day)—Continuous


### 2.2. Forecasting Model Construction

Variable importance is a technique that assigns scores to input variables based on the relative importance of each variable for accurate prediction [34]. Variable importance scores can be evaluated for both classification and regression problems. Variable importance scores play an essential role in constructing forecasting models because they can provide insight into the dataset or the model supported by the dataset [47]. For instance, relative scores can highlight which input variables have the greatest or least effect on the output variable. These scores can be interpreted by domain experts and used as a reference for collecting more or different data. Because most variable importance scores are calculated by forecasting models that fit the dataset, these scores can be calculated for the model interpretation. In addition, analyzing variable importance can offer suggestions to improve the efficiency and effectiveness of the forecasting model through dimensionality reduction and feature selection [48].To date, the R language has been widely used for data cleansing, preparation, and analysis [49]. To facilitate accessibility, we also adopted multiple R packages, including variable importance evaluation functions. To calculate variable importance measures in the R environment, we used Cubist, a regression tree-based model, because it provides a balance between interpretability and predictive power [36]. Figure 4 exhibits the flowchart of interpretable electrical load forecasting based on Cubist modeling. Cubist was developed based on Quinlan’s M5 model tree [48, 50]. The Cubist method generates a series of “if-after-after” rules. Each rule holds a linked multivariate linear model. As long as the covariate set meets the rule conditions, the corresponding model is used to compute the predicted value. The Cubist output includes variable usage statistics and provides the percentage of times each variable was adopted in a condition and/or a linear model.Figure 4
Flowchart of interpretable electrical load forecasting based on Cubist modeling (DA: day-ahead, WA: week-ahead, DPLF: daily peak load forecasting, and TDLF: total daily load forecasting).The general concept of a Cubist regression model can be explained as follows.(1)
While the tree is growing, many leaves and branches grow.(2)
Branches can be considered a series of “if-then” rules, whereas terminal leaves remain a connected multivariate linear model.(3)
Assuming that the covariate set complies with the rule conditions, the relevant model is used to measure the predicted value.The Cubist model sequentially develops a series of trees with adjusted weights and strengthens them with training committees (usually one or more), similar to the “boosting” method. The number of neighbors in the Cubist model is used to correct rule-based predicted values, and the final predicted value denotes a function of all the linear models from the initial node to the terminal node. The percentages displayed in the Cubist output reflect all the models related in the predicted value rather than the terminal models displayed in the output. The variable importance used here is a linear combination of the rule condition usage and model [50].When constructing a forecasting model, datasets are usually divided into a training set and a test set, and then the forecasting model is built using the training set and verified using the test set [35]. However, conventional time-series forecasting model evaluation could exhibit unsatisfactory prediction performance when there is a significant gap between the training set period and the test set period. Additionally, when the dataset is insufficient, it is challenging to obtain satisfactory prediction performance using a small amount of training data [19].To solve these problems, we utilized TSCV based on a rolling forecasting origin [51]. TSCV focuses on a single or several forecast horizons for each test set. In this study, we used several different training sets, each containing one or more observations not used in the previous training set, depending on the scheduling period. To perform DPLF and TDLF, we took the test sets one day after the current time and a week after the current time, as shown in Figure 5. We evaluated the forecasting model performance by calculating the prediction accuracy at each time point and then calculated their average value.Figure 5
. Time-series cross-validation for day-ahead forecasting and week-ahead forecasting. (a) Day-ahead forecasting. (b) Week-ahead forecasting.
(a)(b)Therefore, it is possible to solve the data shortage problem because more data can be used over time than in the conventional time-series forecasting model evaluation. We can also expect satisfactory prediction performance because it can adequately reflect recent electrical load patterns and adjust the weights of the input variables in the forecasting model. Here, we presented interpretable electrical load forecasting results by calculating the importance of the variables for each training set in the model.


## 2.1. Data Collection and Preprocessing

In this section, we first present a data resampling process to construct data suitable for forecasting purposes. Then, we describe the process of configuring the input variables according to the purpose (i.e., DA and WA forecasting) using timestamps, weather, and historical load information. Figure2 illustrates the framework of the data preprocessing for constructing a Cubist model.Figure 2
Data preprocessing for Cubist modeling (DA: day-ahead, WA: week-ahead, DPLF: daily peak load forecasting, and TDLF: total daily load forecasting).We collected electrical load datasets from five different types of buildings or building clusters, as summarized in Table2. We first collected publicly available datasets from two buildings in Richland, WA, USA [39, 40]. The datasets consist of three years’ worth of information, including the hourly electrical load, the hourly outdoor temperature, and the corresponding timestamps. We filled in the missing values (i.e., daylight saving time in North America) in both datasets using linear interpolation. Because the dataset contains the hourly electrical load data for a day (24 rows), we calculated the maximum load value and the sum of all load values for each day for DPLF and TDLF, respectively.Table 2
Building information.
Dataset #Number of buildingsBuilding type (description)LocationDataset periodPublic accessBuilding 11Commercial (office)Richland, WashingtonJan. 2, 2009–Dec. 31, 2011YesBuilding 21Commercial (office)Richland, WashingtonJan. 2, 2009–Dec. 31, 2011YesCluster 116Educational (dormitory)Seoul, South KoreaJan. 1, 2016–Dec. 31, 2018NoCluster 232Educational (humanities bldg.)Seoul, South KoreaJan. 1, 2016–Dec. 31, 2018NoCluster 35Educational (engineering bldg.)Seoul, South KoreaJan. 1, 2016–Dec. 31, 2018NoIn addition, we collected typical 15-minute-interval electrical load datasets from three clusters of buildings at a private university in Seoul, South Korea [41]. The dataset collection period was three years. The first cluster comprised 16 residential buildings, and the electrical load had a residential pattern. The second cluster consisted of 32 academic buildings, including the main hall, library, classrooms, and offices. The third cluster contained five science and engineering buildings. This cluster exhibited much higher electrical loads per building than the other clusters, mainly due to the various experimental equipment and devices in the laboratories. In South Korea, the daily peak load of the building was calculated by multiplying the highest value among the 15-minute-interval electrical load used per day by 4. We took this into account when calculating the daily peak electrical load for DPLF. Likewise, we took the sum of 96 values per day to perform TDLF.Tables3 and 4 provide some statistics on daily peak loads and total daily loads for the five datasets, respectively. Weather conditions and holiday information are closely related to electrical load [42]. Therefore, we used these data to configure the input variables for Cubist modeling. Because we performed two types of forecasting, namely, DA and WA forecasting, we used different input variable configurations for each type of forecasting.Table 3
Statistics on daily peak electrical load data (unit: kW).
StatisticsBuilding 1Building 2Cluster 1Cluster 2Cluster 3Number of valid cases10941094109610961096Mean49.1654.461575.944132.022606.01Standard deviation21.6921.15308.911327.07451.57Trimmed mean50.4056.231552.324325.762623.04Median48.5954.521561.974176.242670.00Median absolute deviation19.8718.52321.311537.16437.66Minimum8.8610.97878.401426.561579.20Maximum141.11135.002623.686900.483549.60Range132.25124.031745.285473.921970.40Skew0.340.050.42–0.25–0.31Kurtosis0.430.09–0.15–1.03–0.71Standard error0.660.649.3340.0913.66Table 4
Statistics on total daily electrical load data (unit: kW).
StatisticsBuilding 1Building 2Cluster 1Cluster 2Cluster 3Number of valid cases10941094109610961096Mean719.31852.6329802.4162563.1049440.52Standard deviation250.53290.955350.3917167.956409.23Trimmed mean723.27850.1529390.7662514.4849506.54Median714.88844.9729569.9262872.2349771.80Median absolute deviation229.40244.275407.2820269.876802.47Minimum198.22242.0919013.0427961.4432546.40Maximum1527.302130.5049235.7698475.8464403.70Range1329.081888.4130222.7270514.4031857.30Skew0.160.440.45–0.13–0.11Kurtosis–0.320.98–0.03–1.06–0.64Standard error7.578.80161.61518.58193.86Time is another important factor for electrical loads. We considered various temporal variables, such as months, days, and days of the week. However, these variables cannot reflect periodic information when applied to forecasting models. For instance, 31 December and 1 January are temporally contiguous, yet, the range of both values in sequence form is 30 and 11 as day and month, respectively. Hence, we represented them as continuous data in the two-dimensional (2D) space to reflect their periodicity using equations (1)–(6) [41, 42]. Here, the variable for seven days of the week is defined as 1 to 7 from Monday to Sunday, according to international standard ISO 8601. LDMmonth represents the last day of the month to which the day belongs (e.g., January: 31, February: 28 or 29, March: 31, and so on). Consequently, we used six input variables to describe the date and time of the prediction time points.(1)Monthx=sin36012×Month,(2)Monthy=cos36012×Month,(3)Dayx=sin360LDMmonth×Day,(4)Dayy=cos360LDMmonth×Day,(5)Day_of_the_Weekx=sin3607×Day_of_the_week,(6)Day_of_the_Weeky=cos3607×Day_of_the_week.To verify the validity and applicability of the 2D representation, we computed several regression statistics on the electrical loads in one-dimensional (1D) space, consisting of three temporal variables (i.e., months, days, and days of the week), and 2D space, as shown in Tables5 and 6. Here, the residual standard error (RSE), multiple R-squared, and adjusted R-squared were calculated using equations (7)–(9).(7)RSE=1n−2∑t=1nyt−yt^2,(8)R2=1−∑yt−yt^2∑yt−y¯2,(9)Adjusted_R2=1−1−R2n−1n−p−1,where yt and yt^ represent the actual and estimated values, respectively, at the time t. y¯ represents the mean of the actual values, and n and p indicate the number of observations and variables, respectively. From the tables, we can observe that the 2D representation exhibits correlations more effectively than the 1D representation.Table 5
Residual standard error and R-squared statistics for daily peak electrical load data.
StatisticsBuilding 1Building 2Cluster 1Cluster 2Cluster 31D2D1D2D1D2D1D2D1D2DResidual standard error (unit: kW)17.717.117.416.5285.8271.411281121365.9351.9Multiple R-squared (unit: %)33.938.632.639.714.623.227.929.136.441.3Adjusted R-squared (unit: %)33.738.232.439.414.422.827.728.736.241.0Table 6
Residual standard error and R-squared statistics for total daily electrical load data.
StatisticsBuilding 1Building 2Cluster 1Cluster 2Cluster 31D2D1D2D1D2D1D2D1D2DResidual standard error (unit: kW)21120325622849584698146801467056995417Multiple R-squared (unit: %)29.334.922.939.014.423.327.127.427.434.6Adjusted R-squared (unit: %)29.134.622.738.714.122.926.927.027.234.2Typical electrical load patterns are different on weekdays and holidays, depending on the type of building [16, 40]. To reflect this, holiday information was also used as an input variable of the forecasting model. We collected holiday information by country from https://www.timeanddate.com [43]. The holidays included Saturdays, Sundays, and public holidays. We used one-hot encoding (i.e., “1” for the relevant data and “0” otherwise) as a nominal scale. Hence, we used seven time factors as input variables at the prediction time point.In general, the electrical load increases in the summer and winter due to the heavy use of air-conditioning and electrical heating appliances, respectively [13, 44]. Here, the most influential factor in the electrical load is temperature [13]. In this study, we focused on temperature-related variables, such as the daily maximum, average, and minimum temperatures [44]. In the case of Buildings 1 and 2, because the datasets we collected also included the outdoor temperature along with the electrical load, we calculated the daily maximum, average, and minimum temperatures (in the Fahrenheit scale). For Clusters 1 to 3, we first collected hourly temperature data using regional synoptic meteorological data provided by the Korea Meteorological Administration (KMA) and collected by the Seoul Meteorological Observatory located about 6 km from the university campus. As in Buildings 1 and 2, we calculated the daily maximum, average, and minimum temperature in Celsius. In South Korea, the KMA provides weather forecasts, including short and mid-term forecasts, for the most significant regions. The KMA mid-term forecast service provides the maximum and minimum temperatures from day 3 to day 10, as shown in Figure 3.Figure 3
Example of mid-term forecast provided by the KMA.Historical electrical load data are particularly important in electrical load forecasting because they exhibit the trend of recent electrical loads [17, 45]. To reflect the recent trend in the prediction, we used the daily peak and total daily loads of the previous seven days as input variables for the DPLF and TDLF models, respectively. Because the electrical load trends of weekdays and holidays can differ, we added a holiday indicator to indicate whether the day was a holiday [46]. We used a total of 24 input variables to build the DA-DPLF and TDLF models.Now, we describe the configuration of the input variables for the WA-DPLF and TDLF models. First, we used the same time and temperature variables as in the DA forecasting models. However, the daily peak and total daily loads from the previous day to the sixth day before the forecast date are unknown when considering WA forecasting. To compensate for this, we configured the input variables using the daily peak and total daily loads and holiday indicators of the same day of the previous four weeks for the WA-DPLF and TDLF, respectively [46]. Therefore, we used a total of 18 input variables to construct our WA electrical load forecasting models. Table 7 presents all the input variables that we considered for DA and WA forecasts.Table 7
Input variables for day-ahead and week-ahead forecasts.
IV #Input variables for day-ahead forecastingInput variables for week-ahead forecastingVariable type01MonthxMonthxContinuous [–1, 1]02MonthyMonthyContinuous [–1, 1]03DayxDayxContinuous [–1, 1]04DayyDayyContinuous [–1, 1]05Day of the weekxDay of the weekxContinuous [–1, 1]06Day of the weekyDay of the weekyContinuous [–1, 1]07HolidayHolidayBinary08Minimum temperatureMinimum temperatureContinuous09Average temperatureAverage temperatureContinuous10Maximum temperatureMaximum temperatureContinuous11Holiday (the day before seven days)Holiday (the day before four weeks)Binary12Electrical load (the day before seven days)Electrical load (the day before four weeks)Continuous13Holiday (the day before six days days)Holiday (the day before six three weeks)Binary14Electrical load (the day before six days)Electrical load (the day before three weeks)Continuous15Holiday (the day before five days)Holiday (the day before five two weeks)Binary16Electrical load (the day before five days)Electrical load (the day before two weeks)Continuous17Holiday (the day before four days)Holiday (the day before one week)Binary18Electrical load (the day before four days)Electrical load (the day before one week)Continuous19Holiday (the day before three days)—Binary20Electrical load (the day before three days)—Continuous21Holiday (the day before two days)—Binary21Electrical load (the day before two days)—Continuous23Holiday (the day before one day)—Binary24Electrical load (the day before one day)—Continuous


## 2.2. Forecasting Model Construction

Variable importance is a technique that assigns scores to input variables based on the relative importance of each variable for accurate prediction [34]. Variable importance scores can be evaluated for both classification and regression problems. Variable importance scores play an essential role in constructing forecasting models because they can provide insight into the dataset or the model supported by the dataset [47]. For instance, relative scores can highlight which input variables have the greatest or least effect on the output variable. These scores can be interpreted by domain experts and used as a reference for collecting more or different data. Because most variable importance scores are calculated by forecasting models that fit the dataset, these scores can be calculated for the model interpretation. In addition, analyzing variable importance can offer suggestions to improve the efficiency and effectiveness of the forecasting model through dimensionality reduction and feature selection [48].To date, the R language has been widely used for data cleansing, preparation, and analysis [49]. To facilitate accessibility, we also adopted multiple R packages, including variable importance evaluation functions. To calculate variable importance measures in the R environment, we used Cubist, a regression tree-based model, because it provides a balance between interpretability and predictive power [36]. Figure 4 exhibits the flowchart of interpretable electrical load forecasting based on Cubist modeling. Cubist was developed based on Quinlan’s M5 model tree [48, 50]. The Cubist method generates a series of “if-after-after” rules. Each rule holds a linked multivariate linear model. As long as the covariate set meets the rule conditions, the corresponding model is used to compute the predicted value. The Cubist output includes variable usage statistics and provides the percentage of times each variable was adopted in a condition and/or a linear model.Figure 4
Flowchart of interpretable electrical load forecasting based on Cubist modeling (DA: day-ahead, WA: week-ahead, DPLF: daily peak load forecasting, and TDLF: total daily load forecasting).The general concept of a Cubist regression model can be explained as follows.(1)
While the tree is growing, many leaves and branches grow.(2)
Branches can be considered a series of “if-then” rules, whereas terminal leaves remain a connected multivariate linear model.(3)
Assuming that the covariate set complies with the rule conditions, the relevant model is used to measure the predicted value.The Cubist model sequentially develops a series of trees with adjusted weights and strengthens them with training committees (usually one or more), similar to the “boosting” method. The number of neighbors in the Cubist model is used to correct rule-based predicted values, and the final predicted value denotes a function of all the linear models from the initial node to the terminal node. The percentages displayed in the Cubist output reflect all the models related in the predicted value rather than the terminal models displayed in the output. The variable importance used here is a linear combination of the rule condition usage and model [50].When constructing a forecasting model, datasets are usually divided into a training set and a test set, and then the forecasting model is built using the training set and verified using the test set [35]. However, conventional time-series forecasting model evaluation could exhibit unsatisfactory prediction performance when there is a significant gap between the training set period and the test set period. Additionally, when the dataset is insufficient, it is challenging to obtain satisfactory prediction performance using a small amount of training data [19].To solve these problems, we utilized TSCV based on a rolling forecasting origin [51]. TSCV focuses on a single or several forecast horizons for each test set. In this study, we used several different training sets, each containing one or more observations not used in the previous training set, depending on the scheduling period. To perform DPLF and TDLF, we took the test sets one day after the current time and a week after the current time, as shown in Figure 5. We evaluated the forecasting model performance by calculating the prediction accuracy at each time point and then calculated their average value.Figure 5
. Time-series cross-validation for day-ahead forecasting and week-ahead forecasting. (a) Day-ahead forecasting. (b) Week-ahead forecasting.
(a)(b)Therefore, it is possible to solve the data shortage problem because more data can be used over time than in the conventional time-series forecasting model evaluation. We can also expect satisfactory prediction performance because it can adequately reflect recent electrical load patterns and adjust the weights of the input variables in the forecasting model. Here, we presented interpretable electrical load forecasting results by calculating the importance of the variables for each training set in the model.


## 3. Results and Discussion

In this section, we first introduce metrics to compare the prediction performance of the forecasting models and describe the experimental design and results in detail. We also present the results of several statistical tests to prove the validity of our experiment. We exhibit the interpretable short-term electrical load forecasting using the proposed model. Finally, we discuss experimental procedures.


### 3.1. Experimental Design

The quantitative experiments were conducted with an Intel® CoreTM i7-8700k CPU with 32 GB DDR4 RAM. We performed the input variable configuration and forecasting model construction in RStudio (v. 1.1.453) with R (v. 3.5.1). We used three years of electrical load data from 2009 to 2011 for Buildings 1 and 2 and from 2016 to 2018 for Clusters 1 to 3, respectively. We divided the dataset into training (in-sample) and test (out-of-sample) sets in an approximate proportion of 2 : 1. For Buildings 1 and 2, we confirmed that the electrical load from November to December 2011, the out-of-sample period, was higher than that for the remaining days. The R random number generator seed was set to 1234 for all methods.To evaluate the predictive performance of forecasting models, we used the MAPE and coefficient of variation of the root mean square error (CVRMSE). MAPE and CVRMSE values show the accuracy as a relative percentage error. Hence, they are easier to understand than other well-known metrics, such as the mean absolute error, mean square error, and root mean square error [40, 52]. The MAPE measures the prediction accuracy for constructing fitted time-series values in statistics, specifically in trend estimation. The CVRMSE is used to aggregate the residuals into a single measure of predictive ability and is more useful when significant errors are particularly undesirable. The lower the MAPE and CVRMSE values, the better the forecasting model’s predictive performance. However, it is known that the MAPE and CVRMSE increase significantly when the actual value tends to zero [35, 52]. The MAPE and CVRMSE are calculated using (10) and (11), respectively, where yt and yt^ are the actual and forecasted values at time t, respectively, y¯ is an average of the actual values, and n is the number of observations.(10)MAPE=100n∑t=1nyt−yt^yt,(11)CVRMSE=100y¯∑t=1nyt−yt^2n.To demonstrate the validity of the proposed model, we considered a total of 12 machine learning methods including MLR, partial least squares (PLS), multivariate adaptive regression splines (MARS), KNN, SVR, DT, Bagging, RF, GBM, XGB, and CatBoost. Most machine learning methods have hyperparameters that can influence model performance. To further improve their performance, it is critical to tune the hyperparameters effectively. To find the optimal hyperparameters, we performed 10-fold cross-validation, all on the training set. Table8 describes the hyperparameters for each method, the R package used, the references where details about their optimal value can be found, and the hyperparameter ranges.Table 8
List of the hyperparameters used to build optimal forecasting models.
MethodsPackageHyperparameters and their rangeMLR [53]lmNone (automatic identification)PLS [53]pls, caretncomp (vector of positive integers): 1 : 1 less than the number of input variablesMARS [49]earth, caretdegree (maximum degree of interactions): 1 : 3nprune (number of terms retained in the final model): 2, 13, 24, 35, 46, 56, 67, 78, 89, 100KNN [21]caretk (number of neighbors): 2SVR [54]kernlab, caretsigma (sigma): 0.35, 0.4, 0.1C (cost): 1, 3, 5, 8, 10, 12DT [53]rpartmaxdepth (maximum depth of any node of the final tree): automatic identificationBagging [53]IpredNone (automatic identification)RF [40, 48]randomForestmtry (number of variables randomly chosen at each split): number of input variables divided by 3ntree (number of trees to grow): 128GBM [55]gbmn.trees (number of trees to grow): 3000, 4000, 5000, 6000, 7000, 8000, 9000, 10000interaction.depth (maximum depth of each tree): 5shrinkage (shrinkage or learning rate parameter): 0.001bag.fraction (subsampling rate): 0.5XGB [56]xgboost, caretnrounds (number of trees to grow): 50, 100, 250, 500eta (shrinkage or learning rate parameter): 0.01, 0.1, 1lambda (L2 regularization term on weights): 0.1, 0.5, 1alpha (L1 regularization term on weights): 0.1, 0.5, 1CatBoost [57]catboost, caretlearning rate (shrinkage or learning rate parameter): 0.03, 0.1depth (maximum depth of each tree): 4, 6, 10l2_leaf_reg (coefficient at the L2 regularization term of the cost function): 1, 3, 5, 7, 9Cubist [50]Cubistcommittees (sequence generation of rule-based models (similar to boosting)): 1, 10, 50, 100neighbors (single integer value to adjust the rule-based predictions from the training set): 0, 1, 5, 9For DA forecasting, we also considered the MLR model in [20], which achieved better prediction performance than KPX’s forecasting model, as a baseline model in the evaluation of the proposed model. The MLR model adapts a rolling procedure using a dataset from one day to one year before the prediction time point and is specified in the following equation:(12)YD=β0+β1YD−1+β2YD−2+β3YD−7+β4WD+β5SD+β6TD+β7SDTD+ϵ,where YD is the expected daily peak or total daily electrical load at day D; YD−1, YD−2, and YD−7 are historical loads before one day, two days, and one week, respectively; WD, which represents the day of the week at day D, is a categorical variable and consists of seven categories from 1 (Monday) to 7 (Sunday); SD is a categorical variable for season and consists of six categories from 1 (Jan. and Feb.) to 6 (Nov. and Dec.); TD is the average temperature at day D; SDTD is a quadratic variable to reflect appropriate weather information according to the season; and β0, βi, and ϵ are the constant term, the slope coefficient of the ith independent variable, and the error term, respectively.


### 3.2. Experimental Results

Tables9–12 present the MAPE and CVRMSE results of DA-DPLF and TDLF, respectively, through conventional time-series forecasting model evaluation, also known as holdout, and TSCV. Tables 13–16 present the MAPE and CVRMSE forecasting results for DPLF and TDLF, respectively. In Tables 9–16, we demonstrate that the Cubist method produced lower MAPE and CVRMSE values than any other forecasting method that we considered. We also exhibit the prediction performance for each time point through holdout and TSCV. In most cases, TSCV exhibited better prediction performance than holdout, except for PLS, where the input variable weight was adjusted in the forecasting model to reflect the recent electrical load pattern.Table 9
MAPE comparison of DA-DPLF (%).
MethodsBuilding 1Building 2Cluster 1Cluster 2Cluster 3HoldoutTSCVHoldoutTSCVHoldoutTSCVHoldoutTSCVHoldoutTSCVMLR21.1520.8818.2717.955.555.558.918.953.683.81PLS21.1033.3018.4528.215.5110.778.9522.903.709.15MARS20.3218.8318.0816.284.964.705.905.952.832.87KNN22.3522.1519.9719.496.656.029.498.803.983.87SVR22.1720.8217.3715.955.944.667.426.562.812.29DT22.1123.2026.0423.108.718.0512.0711.565.225.05Bagging23.0822.1122.1320.567.306.5810.9110.394.674.43RF19.7217.7818.4215.425.905.126.455.992.982.63GBM18.1716.9917.9015.385.424.876.215.822.992.71XGB19.9317.1718.5015.086.235.326.255.702.942.68CatBoost22.6120.3119.3116.886.495.497.486.653.583.05Cubist18.6016.9814.9713.514.904.685.095.033.152.78Values in bold indicate the lowest values for the respective datasets.Table 10
CVRMSE comparison of DA-DPLF (%).
MethodsBuilding 1Building 2Cluster 1Cluster 2Cluster 3HoldoutTSCVHoldoutTSCVHoldoutTSCVHoldoutTSCVHoldoutTSCVMLR24.0323.8519.3519.187.367.3510.0310.005.505.85PLS24.0734.3319.4727.737.3112.8310.0422.385.5011.38MARS28.6623.1118.4717.886.946.456.836.623.413.71KNN29.7527.7022.2822.159.718.7012.0511.065.425.30SVR30.2826.6117.6517.4411.067.668.787.664.533.99DT26.3426.8624.7322.6712.2810.6612.6112.495.815.62Bagging26.9025.1820.5419.5910.819.0111.0410.715.275.00RF27.9323.3318.2617.038.817.197.096.553.573.24GBM26.6822.8718.0816.928.076.756.916.443.533.28XGB29.3723.1318.4317.609.327.277.286.653.563.32CatBoost31.2026.0819.4818.079.457.568.097.314.353.90Cubist26.9023.2217.3416.166.726.295.635.674.453.62Values in bold indicate the lowest values for the respective datasets.Table 11
MAPE comparison of DA-TDLF (%).
MethodsBuilding 1Building 2Cluster 1Cluster 2Cluster 3HoldoutTSCVHoldoutTSCVHoldoutTSCVHoldoutTSCVHoldoutTSCVMLR11.1811.1210.9810.663.613.637.067.112.552.73PLS11.1826.2711.1027.343.608.627.0318.132.577.49MARS9.398.6411.599.973.283.364.304.262.082.09KNN13.0312.0518.8217.676.005.057.426.894.794.63SVR9.848.7417.1410.743.972.615.194.462.892.15DT17.4017.2923.6819.716.966.169.128.864.724.22Bagging15.8915.2418.9016.115.805.018.338.004.183.48RF10.929.9515.019.584.443.974.514.292.462.05GBM10.749.7714.4410.044.083.594.614.422.502.08XGB10.989.8613.669.474.543.844.544.382.632.18CatBoost11.4910.2014.9410.074.923.815.274.702.952.23Cubist8.808.8911.208.393.263.243.553.602.361.98Values in bold indicate the lowest values for the respective datasets.Table 12
CVRMSE comparison of DA-TDLF (%).
MethodsBuilding 1Building 2Cluster 1Cluster 2Cluster 3HoldoutTSCVHoldoutTSCVHoldoutTSCVHoldoutTSCVHoldoutTSCVMLR14.2314.1712.9812.794.844.838.318.314.434.45PLS14.1326.5912.9525.954.8210.198.3119.074.459.14MARS12.5911.0212.8511.604.444.394.995.003.443.82KNN17.2316.7319.3618.618.416.9710.379.656.576.45SVR12.4511.8016.6712.568.674.847.126.055.545.04DT22.0722.3823.0821.579.618.3810.4110.465.414.95Bagging19.9119.9618.1917.538.377.059.579.194.734.04RF14.4813.5914.8911.686.865.655.435.162.952.66GBM13.4712.8614.5412.176.255.105.565.312.942.79XGB14.8613.7714.7112.457.005.455.655.483.123.66CatBoost14.3413.5215.0712.347.345.366.745.863.572.85Cubist12.0911.7213.1610.554.464.364.154.243.273.34Values in bold indicate the lowest values for the respective datasets.Table 13
MAPE comparison of WA-DPLF (%).
DatasetsEvaluationForecasting methodsMLRPLSMARSKNNSVRDTBaggingRFGBMXGBCatBoostCubistBuilding 1Holdout25.8926.2522.1524.1120.0326.9724.6419.9219.4419.4521.1520.68TSCV (avg.)26.0828.6221.0623.5419.3625.5823.5919.4819.1619.0020.1519.13Building 2Holdout21.5221.6019.2821.9718.9224.8222.0418.6018.3117.4418.7616.92TSCV (avg.)21.2623.8917.1021.1716.9923.5820.9715.8715.5615.6516.8714.55Cluster 1Holdout8.008.016.828.537.919.608.897.897.558.257.896.77TSCV (avg.)8.1610.477.018.076.809.478.717.266.997.777.266.53Cluster 2Holdout11.2711.388.759.798.1112.4811.497.757.607.878.487.04TSCV (avg.)11.1818.277.329.887.5611.9211.077.467.187.418.246.87Cluster 3Holdout4.774.853.824.863.145.765.243.503.533.573.933.40TSCV (avg.)4.937.383.574.632.835.564.943.193.223.433.573.25Table 14
CVRMSE comparison of WA-DPLF (%).
DatasetsEvaluationForecasting methodsMLRPLSMARSKNNSVRDTBaggingRFGBMXGBCatBoostCubistBuilding 1Holdout31.4931.5548.7733.5732.3834.5432.3231.4429.7432.2032.9030.09TSCV (avg.)31.4035.8339.1532.8829.0332.1629.1327.8227.4026.7329.6127.44Building 2Holdout23.3123.3619.8323.6219.7723.0720.5419.0918.9319.5120.1719.02TSCV (avg.)23.0427.0318.9922.7919.4323.8419.9818.2317.5018.4219.3118.07Cluster 1Holdout10.8910.899.2111.9711.9413.4012.4511.0510.5711.6910.989.17TSCV (avg.)10.8213.669.2811.399.8112.7611.7710.029.5710.599.878.92Cluster 2Holdout12.1312.199.6011.949.1113.2911.958.478.579.369.418.08TSCV (avg.)12.1820.308.3512.108.5712.9011.608.298.248.709.008.00Cluster 3Holdout6.616.474.585.934.216.335.804.124.124.254.924.14TSCV (avg.)7.1510.074.445.693.896.135.503.943.874.214.434.01Table 15
MAPE comparison of WA-TDLF (%).
DatasetsEvaluationForecasting methodsMLRPLSMARSKNNSVRDTBaggingRFGBMXGBCatBoostCubistBuilding 1Holdout17.1817.6112.6816.0212.5120.4116.6912.5612.4513.3713.3111.62TSCV (avg.)17.2318.3512.0515.1111.9618.4015.9511.8011.8712.0512.5411.46Building 2Holdout17.4717.8615.7721.7817.0124.4219.6416.2816.4515.1316.9815.24TSCV (avg.)16.1916.8012.9520.3711.2817.8215.5910.8211.8810.8512.5710.94Cluster 1Holdout6.766.876.118.586.709.037.737.046.617.087.145.90TSCV (avg.)6.909.055.718.015.498.767.826.466.056.546.255.51Cluster 2Holdout9.739.925.987.906.189.799.135.895.575.856.474.91TSCV (avg.)9.8415.075.637.765.639.819.185.865.535.676.155.11Cluster 3Holdout4.514.543.455.973.885.404.793.203.283.293.632.96TSCV (avg.)4.765.973.055.563.624.874.242.752.862.913.062.88Table 16
CVRMSE comparison of WA-TDLF (%).
DatasetsEvaluationForecasting methodsMLRPLSMARSKNNSVRDTBaggingRFGBMXGBCatBoostCubistBuilding 1Holdout21.7122.2515.1619.5815.7124.9120.6516.3515.1016.5717.4915.04TSCV (avg.)21.5925.7314.1818.2615.3322.9720.2615.4314.8215.7316.3314.82Building 2Holdout19.7320.4316.7522.7717.0023.6719.6916.1516.7715.8817.2016.33TSCV (avg.)19.2722.7714.3921.4413.9520.8317.6813.7313.9814.2115.4213.41Cluster 1Holdout9.219.277.8211.2410.2412.0610.609.578.939.849.677.78TSCV (avg.)9.1711.857.3910.607.9411.4710.528.628.088.738.307.39Cluster 2Holdout11.3511.437.4710.808.0411.4610.166.936.727.377.826.08TSCV (avg.)11.5517.956.6310.627.1811.4810.456.946.747.237.476.27Cluster 3Holdout6.486.244.197.786.555.965.323.823.823.944.353.68TSCV (avg.)7.437.833.957.496.405.694.893.413.443.723.823.65We made the following five observations from the experiments:(1)
The day-ahead electrical load forecasting models exhibited better prediction performance than the week-ahead electrical load forecasting models.(2)
The week-ahead electrical load forecasting models exhibited a lower prediction performance as the prediction time moved further from the current time.(3)
Despite the same time point, the day-ahead electrical load forecasting models showed better prediction performance than the first prediction time point of the week-ahead electrical load forecasting models.(4)
The total electrical load forecasting models displayed higher prediction performance than the peak electrical load forecasting models.(5)
The prediction accuracy for Clusters 1–3 was higher than the prediction accuracy for Buildings 1 and 2.In general, the electrical load pattern of buildings or building clusters could change from a variety of causes, and hence the further apart the current and forecast times, the higher the uncertainty [9]. Therefore, both the DA and the first prediction time point of the WA electrical load forecasting models yielded more accurate predictions because they modeled electrical load patterns of the day before using TSCV. In addition, because the DA electrical load forecasting models utilized the electrical load of the day before as an input variable, they exhibited more accurate prediction performance than the first prediction time point of the WA electrical load forecasting models. Hence, we confirmed that the electrical loads during the week were more significant input variables than the electrical loads of the same days of the week. Overall, we also confirmed that the forecasting models displayed low prediction accuracy when the electrical load was close to zero. From November to December 2011 (out-of-sample), Buildings 1 and 2 showed a sudden higher electrical load than that on other days. Thus, the electrical load forecasting models had a large prediction error because they were not trained on these electrical loads.To find the best forecasting method, we ranked them by considering all the TSCV performance metrics for each building or cluster. We then calculated the average rank using the rank for all buildings and clusters for the DA-DPLF and TDLF and WA-DPLF and TDLF for each method. Table17 presents the ranks for each method for each performance metric and the average ranks. We confirmed that the Cubist method exhibited the best rank in the table. In addition, we demonstrated that the proposed Cubist model outperformed the MLR model in all aspects of DA forecasting, as shown in Table 18.Table 17
Ranks of each model based on the performance metrics and average rank.
MethodsMAPECVRMSEAverage rankDA-DPLFDA-TDLFWA-DPLFWA-TDLFDA-DPLFDA-TDLFWA-DPLFWA-TDLFMLR8.27.49.89.47.879.49.68.6PLS12121211.4121211.81211.9MARS4.82.85.44.63.82.663.24.2KNN99.68.29.89.89.89.29.89.4SVR44.43.83.876.24.64.84.8DT1110.810.810.810.610.610.41010.6Bagging9.89.49.28.68.6988.28.9RF44.44.23.23.24.23.83.63.8GBM3.24.22.63.223.81.833.0XGB3.24.844.64.264.65.64.6CatBoost6.86.65.86.26.856.26.26.2Cubist21.61.822.21.82.21.61.9Values in bold indicate the lowest values for the respective electrical load forecasting types (DA: day-ahead; WA: week-ahead).Table 18
Prediction performance of Lee and Han’s MLR and Cubist.
DatasetsMetricsDPLFTDLFMLR [20]CubistMLR [20]CubistBuilding 1MAPE28.9816.9817.198.89CVRMSE29.5823.2219.5711.72Building 2MAPE24.4413.5116.058.39CVRMSE24.5716.1617.8210.55Cluster 1MAPE7.224.685.003.24CVRMSE9.176.296.494.36Cluster 2MAPE13.625.0316.593.60CVRMSE10.225.6713.664.24Cluster 3MAPE6.202.785.121.98CVRMSE8.833.627.713.34


### 3.3. Statistical Verification

To demonstrate the validity of the proposed method, we conducted three statistical tests. The paired samplet-test was used to demonstrate the effectiveness of TSCV, and the Wilcoxon signed-rank and Friedman tests were used to confirm the difference in prediction performance between the proposed model and the other models. Here, the p value gives the probability of observing test results under the null hypothesis. The cutoff value for determining statistical significance is usually a value of less than 0.05, which corresponds to a 5% or lower chance of obtaining a result like the one observed if the null hypothesis is correct.The paired samplet-test compares two means of the same individual, object, or related units. The test objective is to discover whether there is statistical evidence that the mean difference between a pair of observations for a critical result is significantly different from zero. Typical applications of the paired sample t-test include case-control studies or repeated-measurement designs. Because PLS performed poorly for all forecasting tasks in TSCV, we performed the paired sample t-test considering all the holdout and TSCV values for each method except PLS for MAPE and CVRMSE. Table 19 shows that the MAPE and CVRMSE p values were both less than 0.05. Therefore, we can confirm the validity of the TSCV used in this study.Table 19
Paired samplet-test of holdout and TSCV.
StatisticsMAPECVRMSET-test statistic value (t)11.13611.167Degrees of freedom (df)219219Significance level of thet-test (p value)2.2 × 10−162.2 × 10−16Confidence interval (conf.int) of the mean differences at 95%[0.787, 1.125][0.821, 1.173]Mean differences between pairs (sample estimates)0.9560.997The Wilcoxon signed-rank test [58, 59] is used to confirm the null hypothesis to determine a significant difference between the two models. In contrast, the Friedman test [58] is a multiple-comparison test that aims to identify significant differences between three or more forecasting models. To verify these test results, we used all the MAPE and CVRMSE values (the DA-DPLF and TDLF and the WA-DPLF and TDLF) for each forecasting method. The Wilcoxon signed-rank test results with a significance level of 0.5 and the Friedman test results are shown in Table 20. Because the p value in all cases is below the significance level, the proposed model is superior to the other models.Table 20
Results of Wilcoxon signed-rank and Friedman tests.
MethodsWilcoxon signed-rank testFriedman testMAPECVRMSEMAPECVRMSEMLR1.907 × 10−61.907 × 10−62.2 × 10−162.2 × 10−16PLS1.907 × 10−61.907 × 10−6MARS3.624 × 10−50.005841KNN1.907 × 10−61.907 × 10−6SVR0.0094635.722 × 10−6DT1.907 × 10−61.907 × 10−6Bagging1.907 × 10−61.907 × 10−6RF0.0001680.001432GBM8.202 × 10−50.019580XGB6.294 × 10−50.000210CatBoost1.907 × 10−60.000175


### 3.4. Model Interpretation

We determined the variable importance of the proposed model at each test point (one day or one week), according to the TSCV cycle. Figures6 and 7 present heat map graphs for the DA forecast, revealing the importance of the input variables listed in Table 7. Figures 8 and 9 present heat map graphs for the WA forecast, revealing the importance of the input variables listed in Table 7. The analytical results obtained, as shown in the figures, can be summarized as follows. In the calendar data, we confirmed that the holiday was the most significant variable for the forecasting model, and the variables for the days of the week were highly important. The variables for the day did not significantly affect the model performance. We also confirmed that, overall, the temperature data were essential for the model performance, and the average temperature was a key input variable for the proposed model. In the historical load data, the electrical load from the previous day was the most important variable, and the electrical load one-week prior was also crucial for the DA forecasting model. In the WA forecasting models, the electrical load one-week prior was the most significant input variable in the historical load data. Here, we can see that the adjacent autocorrelation variables (e.g., historical load) are essential for the model performance. We also presume that the DA forecasting models performed better than the WA forecasting models because they could reflect the crucial historical load variables, both the day before and the week before.Figure 6
Example of variable importance for DA-DPLF (Building 1).Figure 7
Example of variable importance for DA-TDLF (Building 1).Figure 8
Example of variable importance for WA-DPLF (Building 1).Figure 9
Example of variable importance for WA-TDLF (Building 1).


### 3.5. Discussion

The experimental results showed that PLS exhibited poor prediction performance in TSCV. PLS is a popular method to deal with multicollinear relationships between output and input variables [48]. Even though we predicted the electrical load by setting the PLS hyperparameter automatically, it performed poorly. Although the historical load was highly correlated with the actual electrical load, PLS did not adequately reflect the historical load in TSCV, and hence we conclude that this caused the prediction accuracy to be low. XGB and CatBoost are state-of-the-art technologies. XGB performed satisfactorily, but CatBoost did not. Because most input variables are continuous, CatBoost could not use its advantages, such as ordered target statistics. Therefore, we concluded that XGB is better suited for DPLF and TDLF than CatBoost. The main advantage of the Cubist method is the addition of multiple training committees and “reinforcement” to balance the weights better. Therefore, we presume that Cubist can achieve satisfactory prediction performance because it predicts the next time point by adjusting the weights of the input variables better through TSCV.The light GBM (LightGBM) and neural network methods also performed well in STLF but were not considered here. These methods require a sufficient dataset for model training [18, 60]; however, it takes a significant amount of time to collect enough data because only one dataset for daily peak load and total daily load is generated per day. Moreover, these methods are better configured for the Python environment [60]. In particular, neural network methods require high-performance computer specifications [9]. In this paper, we only considered several datasets collected over three years, and a little over 700 tuples were used for the first model training. Therefore, we did not expect these methods to formulate a robust forecasting model with such small datasets. In the future, we will apply these methods to perform interpretable hourly electrical load forecasting or VSTLF.


## 3.1. Experimental Design

The quantitative experiments were conducted with an Intel® CoreTM i7-8700k CPU with 32 GB DDR4 RAM. We performed the input variable configuration and forecasting model construction in RStudio (v. 1.1.453) with R (v. 3.5.1). We used three years of electrical load data from 2009 to 2011 for Buildings 1 and 2 and from 2016 to 2018 for Clusters 1 to 3, respectively. We divided the dataset into training (in-sample) and test (out-of-sample) sets in an approximate proportion of 2 : 1. For Buildings 1 and 2, we confirmed that the electrical load from November to December 2011, the out-of-sample period, was higher than that for the remaining days. The R random number generator seed was set to 1234 for all methods.To evaluate the predictive performance of forecasting models, we used the MAPE and coefficient of variation of the root mean square error (CVRMSE). MAPE and CVRMSE values show the accuracy as a relative percentage error. Hence, they are easier to understand than other well-known metrics, such as the mean absolute error, mean square error, and root mean square error [40, 52]. The MAPE measures the prediction accuracy for constructing fitted time-series values in statistics, specifically in trend estimation. The CVRMSE is used to aggregate the residuals into a single measure of predictive ability and is more useful when significant errors are particularly undesirable. The lower the MAPE and CVRMSE values, the better the forecasting model’s predictive performance. However, it is known that the MAPE and CVRMSE increase significantly when the actual value tends to zero [35, 52]. The MAPE and CVRMSE are calculated using (10) and (11), respectively, where yt and yt^ are the actual and forecasted values at time t, respectively, y¯ is an average of the actual values, and n is the number of observations.(10)MAPE=100n∑t=1nyt−yt^yt,(11)CVRMSE=100y¯∑t=1nyt−yt^2n.To demonstrate the validity of the proposed model, we considered a total of 12 machine learning methods including MLR, partial least squares (PLS), multivariate adaptive regression splines (MARS), KNN, SVR, DT, Bagging, RF, GBM, XGB, and CatBoost. Most machine learning methods have hyperparameters that can influence model performance. To further improve their performance, it is critical to tune the hyperparameters effectively. To find the optimal hyperparameters, we performed 10-fold cross-validation, all on the training set. Table8 describes the hyperparameters for each method, the R package used, the references where details about their optimal value can be found, and the hyperparameter ranges.Table 8
List of the hyperparameters used to build optimal forecasting models.
MethodsPackageHyperparameters and their rangeMLR [53]lmNone (automatic identification)PLS [53]pls, caretncomp (vector of positive integers): 1 : 1 less than the number of input variablesMARS [49]earth, caretdegree (maximum degree of interactions): 1 : 3nprune (number of terms retained in the final model): 2, 13, 24, 35, 46, 56, 67, 78, 89, 100KNN [21]caretk (number of neighbors): 2SVR [54]kernlab, caretsigma (sigma): 0.35, 0.4, 0.1C (cost): 1, 3, 5, 8, 10, 12DT [53]rpartmaxdepth (maximum depth of any node of the final tree): automatic identificationBagging [53]IpredNone (automatic identification)RF [40, 48]randomForestmtry (number of variables randomly chosen at each split): number of input variables divided by 3ntree (number of trees to grow): 128GBM [55]gbmn.trees (number of trees to grow): 3000, 4000, 5000, 6000, 7000, 8000, 9000, 10000interaction.depth (maximum depth of each tree): 5shrinkage (shrinkage or learning rate parameter): 0.001bag.fraction (subsampling rate): 0.5XGB [56]xgboost, caretnrounds (number of trees to grow): 50, 100, 250, 500eta (shrinkage or learning rate parameter): 0.01, 0.1, 1lambda (L2 regularization term on weights): 0.1, 0.5, 1alpha (L1 regularization term on weights): 0.1, 0.5, 1CatBoost [57]catboost, caretlearning rate (shrinkage or learning rate parameter): 0.03, 0.1depth (maximum depth of each tree): 4, 6, 10l2_leaf_reg (coefficient at the L2 regularization term of the cost function): 1, 3, 5, 7, 9Cubist [50]Cubistcommittees (sequence generation of rule-based models (similar to boosting)): 1, 10, 50, 100neighbors (single integer value to adjust the rule-based predictions from the training set): 0, 1, 5, 9For DA forecasting, we also considered the MLR model in [20], which achieved better prediction performance than KPX’s forecasting model, as a baseline model in the evaluation of the proposed model. The MLR model adapts a rolling procedure using a dataset from one day to one year before the prediction time point and is specified in the following equation:(12)YD=β0+β1YD−1+β2YD−2+β3YD−7+β4WD+β5SD+β6TD+β7SDTD+ϵ,where YD is the expected daily peak or total daily electrical load at day D; YD−1, YD−2, and YD−7 are historical loads before one day, two days, and one week, respectively; WD, which represents the day of the week at day D, is a categorical variable and consists of seven categories from 1 (Monday) to 7 (Sunday); SD is a categorical variable for season and consists of six categories from 1 (Jan. and Feb.) to 6 (Nov. and Dec.); TD is the average temperature at day D; SDTD is a quadratic variable to reflect appropriate weather information according to the season; and β0, βi, and ϵ are the constant term, the slope coefficient of the ith independent variable, and the error term, respectively.


## 3.2. Experimental Results

Tables9–12 present the MAPE and CVRMSE results of DA-DPLF and TDLF, respectively, through conventional time-series forecasting model evaluation, also known as holdout, and TSCV. Tables 13–16 present the MAPE and CVRMSE forecasting results for DPLF and TDLF, respectively. In Tables 9–16, we demonstrate that the Cubist method produced lower MAPE and CVRMSE values than any other forecasting method that we considered. We also exhibit the prediction performance for each time point through holdout and TSCV. In most cases, TSCV exhibited better prediction performance than holdout, except for PLS, where the input variable weight was adjusted in the forecasting model to reflect the recent electrical load pattern.Table 9
MAPE comparison of DA-DPLF (%).
MethodsBuilding 1Building 2Cluster 1Cluster 2Cluster 3HoldoutTSCVHoldoutTSCVHoldoutTSCVHoldoutTSCVHoldoutTSCVMLR21.1520.8818.2717.955.555.558.918.953.683.81PLS21.1033.3018.4528.215.5110.778.9522.903.709.15MARS20.3218.8318.0816.284.964.705.905.952.832.87KNN22.3522.1519.9719.496.656.029.498.803.983.87SVR22.1720.8217.3715.955.944.667.426.562.812.29DT22.1123.2026.0423.108.718.0512.0711.565.225.05Bagging23.0822.1122.1320.567.306.5810.9110.394.674.43RF19.7217.7818.4215.425.905.126.455.992.982.63GBM18.1716.9917.9015.385.424.876.215.822.992.71XGB19.9317.1718.5015.086.235.326.255.702.942.68CatBoost22.6120.3119.3116.886.495.497.486.653.583.05Cubist18.6016.9814.9713.514.904.685.095.033.152.78Values in bold indicate the lowest values for the respective datasets.Table 10
CVRMSE comparison of DA-DPLF (%).
MethodsBuilding 1Building 2Cluster 1Cluster 2Cluster 3HoldoutTSCVHoldoutTSCVHoldoutTSCVHoldoutTSCVHoldoutTSCVMLR24.0323.8519.3519.187.367.3510.0310.005.505.85PLS24.0734.3319.4727.737.3112.8310.0422.385.5011.38MARS28.6623.1118.4717.886.946.456.836.623.413.71KNN29.7527.7022.2822.159.718.7012.0511.065.425.30SVR30.2826.6117.6517.4411.067.668.787.664.533.99DT26.3426.8624.7322.6712.2810.6612.6112.495.815.62Bagging26.9025.1820.5419.5910.819.0111.0410.715.275.00RF27.9323.3318.2617.038.817.197.096.553.573.24GBM26.6822.8718.0816.928.076.756.916.443.533.28XGB29.3723.1318.4317.609.327.277.286.653.563.32CatBoost31.2026.0819.4818.079.457.568.097.314.353.90Cubist26.9023.2217.3416.166.726.295.635.674.453.62Values in bold indicate the lowest values for the respective datasets.Table 11
MAPE comparison of DA-TDLF (%).
MethodsBuilding 1Building 2Cluster 1Cluster 2Cluster 3HoldoutTSCVHoldoutTSCVHoldoutTSCVHoldoutTSCVHoldoutTSCVMLR11.1811.1210.9810.663.613.637.067.112.552.73PLS11.1826.2711.1027.343.608.627.0318.132.577.49MARS9.398.6411.599.973.283.364.304.262.082.09KNN13.0312.0518.8217.676.005.057.426.894.794.63SVR9.848.7417.1410.743.972.615.194.462.892.15DT17.4017.2923.6819.716.966.169.128.864.724.22Bagging15.8915.2418.9016.115.805.018.338.004.183.48RF10.929.9515.019.584.443.974.514.292.462.05GBM10.749.7714.4410.044.083.594.614.422.502.08XGB10.989.8613.669.474.543.844.544.382.632.18CatBoost11.4910.2014.9410.074.923.815.274.702.952.23Cubist8.808.8911.208.393.263.243.553.602.361.98Values in bold indicate the lowest values for the respective datasets.Table 12
CVRMSE comparison of DA-TDLF (%).
MethodsBuilding 1Building 2Cluster 1Cluster 2Cluster 3HoldoutTSCVHoldoutTSCVHoldoutTSCVHoldoutTSCVHoldoutTSCVMLR14.2314.1712.9812.794.844.838.318.314.434.45PLS14.1326.5912.9525.954.8210.198.3119.074.459.14MARS12.5911.0212.8511.604.444.394.995.003.443.82KNN17.2316.7319.3618.618.416.9710.379.656.576.45SVR12.4511.8016.6712.568.674.847.126.055.545.04DT22.0722.3823.0821.579.618.3810.4110.465.414.95Bagging19.9119.9618.1917.538.377.059.579.194.734.04RF14.4813.5914.8911.686.865.655.435.162.952.66GBM13.4712.8614.5412.176.255.105.565.312.942.79XGB14.8613.7714.7112.457.005.455.655.483.123.66CatBoost14.3413.5215.0712.347.345.366.745.863.572.85Cubist12.0911.7213.1610.554.464.364.154.243.273.34Values in bold indicate the lowest values for the respective datasets.Table 13
MAPE comparison of WA-DPLF (%).
DatasetsEvaluationForecasting methodsMLRPLSMARSKNNSVRDTBaggingRFGBMXGBCatBoostCubistBuilding 1Holdout25.8926.2522.1524.1120.0326.9724.6419.9219.4419.4521.1520.68TSCV (avg.)26.0828.6221.0623.5419.3625.5823.5919.4819.1619.0020.1519.13Building 2Holdout21.5221.6019.2821.9718.9224.8222.0418.6018.3117.4418.7616.92TSCV (avg.)21.2623.8917.1021.1716.9923.5820.9715.8715.5615.6516.8714.55Cluster 1Holdout8.008.016.828.537.919.608.897.897.558.257.896.77TSCV (avg.)8.1610.477.018.076.809.478.717.266.997.777.266.53Cluster 2Holdout11.2711.388.759.798.1112.4811.497.757.607.878.487.04TSCV (avg.)11.1818.277.329.887.5611.9211.077.467.187.418.246.87Cluster 3Holdout4.774.853.824.863.145.765.243.503.533.573.933.40TSCV (avg.)4.937.383.574.632.835.564.943.193.223.433.573.25Table 14
CVRMSE comparison of WA-DPLF (%).
DatasetsEvaluationForecasting methodsMLRPLSMARSKNNSVRDTBaggingRFGBMXGBCatBoostCubistBuilding 1Holdout31.4931.5548.7733.5732.3834.5432.3231.4429.7432.2032.9030.09TSCV (avg.)31.4035.8339.1532.8829.0332.1629.1327.8227.4026.7329.6127.44Building 2Holdout23.3123.3619.8323.6219.7723.0720.5419.0918.9319.5120.1719.02TSCV (avg.)23.0427.0318.9922.7919.4323.8419.9818.2317.5018.4219.3118.07Cluster 1Holdout10.8910.899.2111.9711.9413.4012.4511.0510.5711.6910.989.17TSCV (avg.)10.8213.669.2811.399.8112.7611.7710.029.5710.599.878.92Cluster 2Holdout12.1312.199.6011.949.1113.2911.958.478.579.369.418.08TSCV (avg.)12.1820.308.3512.108.5712.9011.608.298.248.709.008.00Cluster 3Holdout6.616.474.585.934.216.335.804.124.124.254.924.14TSCV (avg.)7.1510.074.445.693.896.135.503.943.874.214.434.01Table 15
MAPE comparison of WA-TDLF (%).
DatasetsEvaluationForecasting methodsMLRPLSMARSKNNSVRDTBaggingRFGBMXGBCatBoostCubistBuilding 1Holdout17.1817.6112.6816.0212.5120.4116.6912.5612.4513.3713.3111.62TSCV (avg.)17.2318.3512.0515.1111.9618.4015.9511.8011.8712.0512.5411.46Building 2Holdout17.4717.8615.7721.7817.0124.4219.6416.2816.4515.1316.9815.24TSCV (avg.)16.1916.8012.9520.3711.2817.8215.5910.8211.8810.8512.5710.94Cluster 1Holdout6.766.876.118.586.709.037.737.046.617.087.145.90TSCV (avg.)6.909.055.718.015.498.767.826.466.056.546.255.51Cluster 2Holdout9.739.925.987.906.189.799.135.895.575.856.474.91TSCV (avg.)9.8415.075.637.765.639.819.185.865.535.676.155.11Cluster 3Holdout4.514.543.455.973.885.404.793.203.283.293.632.96TSCV (avg.)4.765.973.055.563.624.874.242.752.862.913.062.88Table 16
CVRMSE comparison of WA-TDLF (%).
DatasetsEvaluationForecasting methodsMLRPLSMARSKNNSVRDTBaggingRFGBMXGBCatBoostCubistBuilding 1Holdout21.7122.2515.1619.5815.7124.9120.6516.3515.1016.5717.4915.04TSCV (avg.)21.5925.7314.1818.2615.3322.9720.2615.4314.8215.7316.3314.82Building 2Holdout19.7320.4316.7522.7717.0023.6719.6916.1516.7715.8817.2016.33TSCV (avg.)19.2722.7714.3921.4413.9520.8317.6813.7313.9814.2115.4213.41Cluster 1Holdout9.219.277.8211.2410.2412.0610.609.578.939.849.677.78TSCV (avg.)9.1711.857.3910.607.9411.4710.528.628.088.738.307.39Cluster 2Holdout11.3511.437.4710.808.0411.4610.166.936.727.377.826.08TSCV (avg.)11.5517.956.6310.627.1811.4810.456.946.747.237.476.27Cluster 3Holdout6.486.244.197.786.555.965.323.823.823.944.353.68TSCV (avg.)7.437.833.957.496.405.694.893.413.443.723.823.65We made the following five observations from the experiments:(1)
The day-ahead electrical load forecasting models exhibited better prediction performance than the week-ahead electrical load forecasting models.(2)
The week-ahead electrical load forecasting models exhibited a lower prediction performance as the prediction time moved further from the current time.(3)
Despite the same time point, the day-ahead electrical load forecasting models showed better prediction performance than the first prediction time point of the week-ahead electrical load forecasting models.(4)
The total electrical load forecasting models displayed higher prediction performance than the peak electrical load forecasting models.(5)
The prediction accuracy for Clusters 1–3 was higher than the prediction accuracy for Buildings 1 and 2.In general, the electrical load pattern of buildings or building clusters could change from a variety of causes, and hence the further apart the current and forecast times, the higher the uncertainty [9]. Therefore, both the DA and the first prediction time point of the WA electrical load forecasting models yielded more accurate predictions because they modeled electrical load patterns of the day before using TSCV. In addition, because the DA electrical load forecasting models utilized the electrical load of the day before as an input variable, they exhibited more accurate prediction performance than the first prediction time point of the WA electrical load forecasting models. Hence, we confirmed that the electrical loads during the week were more significant input variables than the electrical loads of the same days of the week. Overall, we also confirmed that the forecasting models displayed low prediction accuracy when the electrical load was close to zero. From November to December 2011 (out-of-sample), Buildings 1 and 2 showed a sudden higher electrical load than that on other days. Thus, the electrical load forecasting models had a large prediction error because they were not trained on these electrical loads.To find the best forecasting method, we ranked them by considering all the TSCV performance metrics for each building or cluster. We then calculated the average rank using the rank for all buildings and clusters for the DA-DPLF and TDLF and WA-DPLF and TDLF for each method. Table17 presents the ranks for each method for each performance metric and the average ranks. We confirmed that the Cubist method exhibited the best rank in the table. In addition, we demonstrated that the proposed Cubist model outperformed the MLR model in all aspects of DA forecasting, as shown in Table 18.Table 17
Ranks of each model based on the performance metrics and average rank.
MethodsMAPECVRMSEAverage rankDA-DPLFDA-TDLFWA-DPLFWA-TDLFDA-DPLFDA-TDLFWA-DPLFWA-TDLFMLR8.27.49.89.47.879.49.68.6PLS12121211.4121211.81211.9MARS4.82.85.44.63.82.663.24.2KNN99.68.29.89.89.89.29.89.4SVR44.43.83.876.24.64.84.8DT1110.810.810.810.610.610.41010.6Bagging9.89.49.28.68.6988.28.9RF44.44.23.23.24.23.83.63.8GBM3.24.22.63.223.81.833.0XGB3.24.844.64.264.65.64.6CatBoost6.86.65.86.26.856.26.26.2Cubist21.61.822.21.82.21.61.9Values in bold indicate the lowest values for the respective electrical load forecasting types (DA: day-ahead; WA: week-ahead).Table 18
Prediction performance of Lee and Han’s MLR and Cubist.
DatasetsMetricsDPLFTDLFMLR [20]CubistMLR [20]CubistBuilding 1MAPE28.9816.9817.198.89CVRMSE29.5823.2219.5711.72Building 2MAPE24.4413.5116.058.39CVRMSE24.5716.1617.8210.55Cluster 1MAPE7.224.685.003.24CVRMSE9.176.296.494.36Cluster 2MAPE13.625.0316.593.60CVRMSE10.225.6713.664.24Cluster 3MAPE6.202.785.121.98CVRMSE8.833.627.713.34


## 3.3. Statistical Verification

To demonstrate the validity of the proposed method, we conducted three statistical tests. The paired samplet-test was used to demonstrate the effectiveness of TSCV, and the Wilcoxon signed-rank and Friedman tests were used to confirm the difference in prediction performance between the proposed model and the other models. Here, the p value gives the probability of observing test results under the null hypothesis. The cutoff value for determining statistical significance is usually a value of less than 0.05, which corresponds to a 5% or lower chance of obtaining a result like the one observed if the null hypothesis is correct.The paired samplet-test compares two means of the same individual, object, or related units. The test objective is to discover whether there is statistical evidence that the mean difference between a pair of observations for a critical result is significantly different from zero. Typical applications of the paired sample t-test include case-control studies or repeated-measurement designs. Because PLS performed poorly for all forecasting tasks in TSCV, we performed the paired sample t-test considering all the holdout and TSCV values for each method except PLS for MAPE and CVRMSE. Table 19 shows that the MAPE and CVRMSE p values were both less than 0.05. Therefore, we can confirm the validity of the TSCV used in this study.Table 19
Paired samplet-test of holdout and TSCV.
StatisticsMAPECVRMSET-test statistic value (t)11.13611.167Degrees of freedom (df)219219Significance level of thet-test (p value)2.2 × 10−162.2 × 10−16Confidence interval (conf.int) of the mean differences at 95%[0.787, 1.125][0.821, 1.173]Mean differences between pairs (sample estimates)0.9560.997The Wilcoxon signed-rank test [58, 59] is used to confirm the null hypothesis to determine a significant difference between the two models. In contrast, the Friedman test [58] is a multiple-comparison test that aims to identify significant differences between three or more forecasting models. To verify these test results, we used all the MAPE and CVRMSE values (the DA-DPLF and TDLF and the WA-DPLF and TDLF) for each forecasting method. The Wilcoxon signed-rank test results with a significance level of 0.5 and the Friedman test results are shown in Table 20. Because the p value in all cases is below the significance level, the proposed model is superior to the other models.Table 20
Results of Wilcoxon signed-rank and Friedman tests.
MethodsWilcoxon signed-rank testFriedman testMAPECVRMSEMAPECVRMSEMLR1.907 × 10−61.907 × 10−62.2 × 10−162.2 × 10−16PLS1.907 × 10−61.907 × 10−6MARS3.624 × 10−50.005841KNN1.907 × 10−61.907 × 10−6SVR0.0094635.722 × 10−6DT1.907 × 10−61.907 × 10−6Bagging1.907 × 10−61.907 × 10−6RF0.0001680.001432GBM8.202 × 10−50.019580XGB6.294 × 10−50.000210CatBoost1.907 × 10−60.000175


## 3.4. Model Interpretation

We determined the variable importance of the proposed model at each test point (one day or one week), according to the TSCV cycle. Figures6 and 7 present heat map graphs for the DA forecast, revealing the importance of the input variables listed in Table 7. Figures 8 and 9 present heat map graphs for the WA forecast, revealing the importance of the input variables listed in Table 7. The analytical results obtained, as shown in the figures, can be summarized as follows. In the calendar data, we confirmed that the holiday was the most significant variable for the forecasting model, and the variables for the days of the week were highly important. The variables for the day did not significantly affect the model performance. We also confirmed that, overall, the temperature data were essential for the model performance, and the average temperature was a key input variable for the proposed model. In the historical load data, the electrical load from the previous day was the most important variable, and the electrical load one-week prior was also crucial for the DA forecasting model. In the WA forecasting models, the electrical load one-week prior was the most significant input variable in the historical load data. Here, we can see that the adjacent autocorrelation variables (e.g., historical load) are essential for the model performance. We also presume that the DA forecasting models performed better than the WA forecasting models because they could reflect the crucial historical load variables, both the day before and the week before.Figure 6
Example of variable importance for DA-DPLF (Building 1).Figure 7
Example of variable importance for DA-TDLF (Building 1).Figure 8
Example of variable importance for WA-DPLF (Building 1).Figure 9
Example of variable importance for WA-TDLF (Building 1).


## 3.5. Discussion

The experimental results showed that PLS exhibited poor prediction performance in TSCV. PLS is a popular method to deal with multicollinear relationships between output and input variables [48]. Even though we predicted the electrical load by setting the PLS hyperparameter automatically, it performed poorly. Although the historical load was highly correlated with the actual electrical load, PLS did not adequately reflect the historical load in TSCV, and hence we conclude that this caused the prediction accuracy to be low. XGB and CatBoost are state-of-the-art technologies. XGB performed satisfactorily, but CatBoost did not. Because most input variables are continuous, CatBoost could not use its advantages, such as ordered target statistics. Therefore, we concluded that XGB is better suited for DPLF and TDLF than CatBoost. The main advantage of the Cubist method is the addition of multiple training committees and “reinforcement” to balance the weights better. Therefore, we presume that Cubist can achieve satisfactory prediction performance because it predicts the next time point by adjusting the weights of the input variables better through TSCV.The light GBM (LightGBM) and neural network methods also performed well in STLF but were not considered here. These methods require a sufficient dataset for model training [18, 60]; however, it takes a significant amount of time to collect enough data because only one dataset for daily peak load and total daily load is generated per day. Moreover, these methods are better configured for the Python environment [60]. In particular, neural network methods require high-performance computer specifications [9]. In this paper, we only considered several datasets collected over three years, and a little over 700 tuples were used for the first model training. Therefore, we did not expect these methods to formulate a robust forecasting model with such small datasets. In the future, we will apply these methods to perform interpretable hourly electrical load forecasting or VSTLF.


## 4. Conclusions

In this paper, we developed a novel forecasting model for interpretable short-term electrical load forecasting. To do this, we collected five different electrical load datasets with temperature and holiday information. We constructed different input variables by considering four forecasting types: day-ahead DPLF and TDLF and week-ahead DPLF and TDLF. We built the proposed model based on Cubist, a rule-based model, and applied TSCV to address the lack of data and reflect recent electrical load trends. The experimental results demonstrated that the proposed model showed excellent prediction performance. In addition, we conducted interpretable electrical load forecasting for each building or building cluster using the variable importance produced by the proposed model.We found that applying the TSCV method can improve prediction performance, except for PLS, and that the Cubist method performed satisfactorily using a small dataset. It was challenging for CatBoost, a state-of-the-art technique, to produce excellent prediction performance because almost all input variables were configured as continuous. Overall, we confirmed that the higher the electrical load, the higher the prediction accuracy. TDLF and the day-ahead forecasting model had a better prediction performance than DPLF and the weak-ahead forecasting model. However, it was difficult to adequately train the forecasting models on sudden electrical load fluctuations because the amount of data was smaller than the amount of data for hourly load forecasting or VSTLF.We plan to perform interpretable VSTLF, such as 10-minute or 15-minute-interval load forecasting, using neural network methods such as activation maps or an attention mechanism. In addition, we will make an effort to develop various methodologies for explainable forecasting in interpretable forecasting. We also plan to find variables that can reflect building characteristics and include them in the forecasting model.



---

*Source: 6892995-2022-02-08.xml*