Use of Available Daylight to Improve Short-Term Load Forecasting Accuracy

: This paper introduces a new methodology to include daylight information in short-term load forecasting (STLF) models. The relation between daylight and power consumption is obvious due to the use of electricity in lighting in general. Nevertheless, very few STLF systems include this variable as an input. In addition, an analysis of one of the current STLF models at the Spanish Transmission System Operator (TSO), shows two humps in its error proﬁle, occurring at sunrise and sunset times. The new methodology includes properly treated daylight information in STLF models in order to reduce the forecasting error during sunrise and sunset, especially when daylight savings time (DST) one-hour time shifts occur. This paper describes the raw information and the linearization method needed. The forecasting model used as the benchmark is currently used at the TSO’s headquarters and it uses both autoregressive (AR) and neural network (NN) components. The method has been designed with data from the Spanish electric system from 2011 to 2017 and tested over 2018 data. The results include a justiﬁcation to use the proposed linearization over other techniques as well as a thorough analysis of the forecast results yielding an error reduction in sunset hours from 1.56% to 1.38% for the AR model and from 1.37% to 1.30% for the combined forecast. In addition, during the weeks in which DST shifts are implemented, sunset error drops from 2.53% to 2.09%.


Introduction
Short-term load forecasting (STLF) is one of the key functions of a Transmission System Operator (TSO) in maintaining technical viability of the system with low operational costs. STLF includes forecasts made from one hour to several days ahead. It provides useful information for system operators to guarantee the reliability of the system, for generators to optimize schedules [1] and for market participants to generate market biddings on both sides of the market.
STLF has been a subject of research for several decades. However, its relevancy has not decayed [2][3][4][5]. Several techniques are used as forecasting engines: linear regressive models [6][7][8], neural networks [6,9,10] and other sorts of artificial intelligence like fuzzy logic [9,11], or evolutionary algorithms [11][12][13][14], and recurrent neural networks [15,16]. Recent advances in technology have freed the use of neural networks with more than the typical three layers [17]. These deep neural networks (DNNs) have been applied to load forecasting [18][19][20] and they have allowed researchers to reduce the workload of hand-designed feature inputs, which makes them good candidates for smart grid applications [21]. However, when forecasting large systems, the most relevant factor or innovation may not be the forecasting engine but rather the selection and treatment of the input information fed to these engines. Therefore, for a 40 GW system like Spain's inland network, it is still worth designing these stages specifically for the system.
There are several stages in input treatment and selection. The following paragraphs review different techniques used to address this issue. - The relation between sunrise and sunset times and hourly load is modeled. -A valid linearization for this relation is presented so that it can be included in both linear or non-linear models is presented. -This new input is included in both types of models to reduce the forecasting error on sunrise and sunset times, especially when sunrise and sunset times vary faster and when DST time-shifts are implemented.
This paper is organized as follows. Section 2 describes the mathematical models tested, the data used and a detailed description and justification of the linearization proposed. It also includes a description of the different tests carried out in the design and assessment phases. In Section 3, the results of each test are detailed, to guide the reader through the design decisions and the assessment conclusions. Section 4 includes the conclusions and a recommendation to roll-out this methodology into the current forecasting system at the Spanish TSO.

Materials and Methods
The main objective of this paper is to reduce forecasting error on hours near sunrise and sunset times. To this end, the starting point for this research will be the forecasting model currently at use by the Spanish TSO, which was designed by our team and is thoroughly described in [6]. The most relevant characteristics of the data used by this model and its mathematical features are explained in the following section to provide a proper background to the reader. Nevertheless, the section will focus on the new additions to this model: characteristics of the daylight variables, treatment and linearization of these variables and modifications to the existing model. Further details about the input data treatment, configuration of forecasting engines, combination of output and general data flow in the original model, can be found in [6].

Model Structure
The model employs a hybrid technique in which two forecasting engines are used to provide hourly forecasts for the day in course and the next 9 days. This paper will focus on the forecast for the next day, as most papers do in order to provide a more comparable result. One of the forecasting engines uses an autoregressive neural network (NARX) while the other one uses a linear autoregressive model with errors (AR). The input analysis carried out in this paper will focus on the autoregressive (AR) model because it allows interpreting the effect of each variable. Still, the effect of the new input will be tested on both engines.
The models consider each hour of the day individually. This means that to forecast a full day, 24 models are used, each one having the same input structure. This allows for a variation of each variable's coefficient throughout the day but causes each hour to be forecasted independently from the adjacent values.
The autoregressive part of the model significantly reduces the forecasting error by using previous errors to capture effects that are not modeled by the predictors used. This feature is beneficial in forecasting, but in modeling it may cover up the effect of the predictors-both used and missing.
In order to overcome the issues described, the AR model is modified to better suit our purposes. Firstly, the AR components are removed for the analysis of the input variables. Nevertheless, the full model with AR components will be tested at the end, once the variables are properly defined. Secondly, in order to understand how daylight affects load, it is not possible to limit each model to one hour because the variation of daylight within one hour throughout the year is too short. Therefore, the daylight variables need to be considered as global variables. The proposed model is, therefore, a linear regression with exogenous variables that affect differently each hour but also with common variables regarding daylight. It is described in (1): where Y h is a vector with the natural logarithm of the load values at hour h each day, X h is a matrix with the values of all exogenous variables for hour h each day, θ h is the vector of coefficients for each variable (the same variable has a different coefficient for each hour), ε h is a vector containing the residuals, D up and D down are the matrices containing the piecewise variables and ϕ up and ϕ down are the vectors of coefficients for the daylight variables linearization. The model described in (1) will be used to analyze the effect of the daylight variables and design their proper pre-treatment. However, once these variables have been designed, they will be introduced into the initial AR and NARX models.

Load Data
The load data are used both as an input and as an output in the model. For this research, hourly data from the Spanish inland system from 2011 to 2018 [39] have been used. The use of large databases is important when using variables that repeat over long periods of time (temperature, daylight or holidays). Otherwise, the modeled behavior may be partial. In addition, databases may differ in predictability which increases the difficulty of comparing and analyzing accuracy results [3,40]. To avoid this issue, Table 1 includes a measure of approximate entropy (ApEn) as a measure of unpredictability. More details about this metric and its parameters can be found in [41]. Load data are used as an input in two ways: -Long-term trend: Economic growth is the main driver for long-term trends in electricity demand in Spain. The base model uses a quadratic polynomial of time to model these trends, but when larger periods of data are used for training this approach is no longer valid. In order to use more training years, the linear and quadratic terms are substituted by a 52-week moving average of the load. -Recent trend: Load series are highly autocorrelated. Therefore, even if the most relevant predictors are used, a recent value of the series has hidden information from which accuracy of the model may benefit. To this end, the base model includes the most recent known value at the time of the forecast as an input. Nevertheless, this value will not be used in the design stage of our model because it may cover up part of the effect of the variables to be analyzed.
Finally, as in most other models, the output of the model is calculated as the natural logarithm of the load.

Temperature Data
The temperature database that was used contains daily minimum and maximum values for the same 2011-2018 period from 59 weather stations scattered throughout the country [42]. Data series from each station are highly correlated to each other and carry very similar information. The selection of the relevant stations, linearization of the variables using the concept of heating degree days (HDD) and cooling degree days (CDD) and the determination of the lags to use is done by empiric procedure described in [6]. The results for location (Madrid-MAD, Barcelona-BCN, Seville-SEV, Zaragoza-ZAR and Bilbao-BIL) and lags used that were significant at a 0.05 level are described in Table 2: Table 2. Temperature variables significant at a 0.05 level.
The type of day is a very important aspect that determines the load profile. The classification system used is very complex and it is thoroughly defined in [29]. It uses binary variables to categorize days considering the many different circumstances that affect special days in Spain.

Daylight
In order to include variables that represent how the consumers' behavior is affected by available daylight at any given hour, literature suggests [43][44][45] that some hours are affected by sunset time while others are affected by sunrise. Midnight and midday hours are normally considered unaffected. In [46], sunrise and sunset time was used to quantify the effect of daylight savings time in the Spanish systems. The following paragraphs present a similar approach to include these variables in STLF models. Two raw variables are considered to capture the influence of daylight: number of hours to sunrise and number of hours to sunset, as described in (2).
where DL up and DL down are the morning and evening variables, d is the date, h is the hour, t sunrise (d) and t sunset (d) is the sunrise and sunset times for day d in Madrid and L 1 and L 2 are the first and last hours of the morning interval. The model will consider that all hours can be affected by one event, therefore all hours are assigned one variable. The threshold that splits which hours are assigned to sunrise and which are assigned to sunset will be determined empirically. In any case, it should not be a critical parameter because hours adjacent to such threshold should not be affected by either event. Sunrise and sunset times are calculated for the coordinates of Madrid.
In order to understand how the daylight information should be introduced in the model, it is important to understand how it affects the load profile. The clock changes that happen when DST is set on and off are a good opportunity to visualize this. Figure 2 shows how load profile changes on a Sunday before and after the spring time shift. The evening change is much more noticeable than the morning change and midday and midnight hours seem to be unaffected. This means that the relation between load at a given hour and the number of hours to sunrise or sunset is not linear. Therefore, it is necessary to linearize the variables before entering the model.

Linearization of Daylight Variables
In order to linearize the relation between load and our daylight variables, we need to visualize the effect that they have on the load. To this end, a piecewise linearization is a good starting point to understand the relation among variables:

Piecewise Linearization
Piecewise linearization is used to separate each variable in intervals in which the relation between variables could be assumed to be linear and continuous among intervals. The span of sunset and sunrise intervals is separated into n intervals (INT 1 . . . n ) defined by n + 1 thresholds (TH 0 . . . n ). The number of intervals is set empirically to an interval length of 60 min after testing 15-, 30-, 60-and 120-min long intervals. Longer intervals caused larger fitting errors while shorter intervals provoked overfitting.
The linearization is realized by creating two new variables (slope and intercept) for each interval. Both initial variables DL up and DL down are linearized so a total of 2 × n × 2 = 4n variables are created. These variables are described in (3): where M up/dw,i K up/dw,i are the slope and intercepts variables for the ith interval for either the morning or evening function, respectively, and DL up/dw(d,h) is either the morning or evening function. These 2 times 2n variables provide a slope and an intercept for each of the intervals defined. However, this does not imply continuity of the response. In order to force continuity in the domain of the variables, it is necessary to introduce a series of restrictions forcing each adjacent interval to have the same value at the shared threshold. These restrictions are described in (4): where m i and k i are the slope and intercept coefficients for the ith interval. These restrictions reduce the number of variables because they restrict the value of the coefficient. Each restriction subtracts one variable, so the new number of variables is 2 × (2n − (n − 1)) = 2 (n +1). These new variables are described in (5): where Vup/dw,i(d,h) is the ith variable of either the sunrise or sunset group. Each variable takes a value of zero if the hour h does not belong to the corresponding group. The model described in (1) allows us to interpret the effect that the daylight variables have on electricity consumption. By isolating the part of the model that deals with daylight variables it is possible to obtain the corresponding coefficient for each instant. These coefficients are the response of the load to variations in available daylight and they are described in (6):

Sigmoid Linearization
The curves seen in Figure 3 show a stable coefficient for hours far before sunrise or sunset; there is a steep change around each event and again a constant coefficient for hours long after. This behavior can be approximated by the sigmoid function described in (7): where b, L, k and DL 0 are the parameters that determine the size of the step, its position on both axes and the slope in the transition. These parameters have different values for the sunrise and sunset linearizations and can be obtained from fitting the sigmoid response to the piecewise response. The results are shown in (8): where R t (−∞) is the natural logarithm of the average value of the right side of the curve, R t (∞) is the natural logarithm of the average value of the left side and R' t represents the slope.
where Rt(−∞) is the natural logarithm of the average value of the right side of the curve, Rt(∞) is the natural logarithm of the average value of the left side and R't represents the slope. Figure 3 shows the comparison of both linearizations. The curves include the intervals in which each hour may vary throughout the year. Due to the current DST policy, the time to sunrise for any hour varies only two hours from winter to summer while the time to sunset varies four hours throughout the year.

Validation of Linearized Variables
Both linearizations show a similar behavior but the sigmoid linearization requires less variables. In order to assess the quality of each linearization, the modeling error of each model will be compared. As a base reference, a third model without any information regarding daylight will also be tested.  Figure 3 shows the comparison of both linearizations. The curves include the intervals in which each hour may vary throughout the year. Due to the current DST policy, the time to sunrise for any hour varies only two hours from winter to summer while the time to sunset varies four hours throughout the year.

Validation of Linearized Variables
Both linearizations show a similar behavior but the sigmoid linearization requires less variables. In order to assess the quality of each linearization, the modeling error of each model will be compared. As a base reference, a third model without any information regarding daylight will also be tested.
The test consists on an out-of-sample test using data from 2011 to 2017 in which data from every other year is used to model each year. The results are shown in Table 3. Other modeling and forecasting errors from the same database can be found in [6,37] for context. The results show that both linearizations provide a very similar result and that, as expected, their use drastically reduces the modeling error at sunrise and sunset hours. Considering that the sigmoid linearization matches the piecewise accuracy with a fraction of the variables, the sigmoid linearization is chosen.

Other Modeling Improvements
The variations in load as available daylight changes is due to consumers' behavior. Still, it is arguable that their behavior will change in the same way, regardless of the day of the week. In order to allow for this variation for different types of days, the model includes a different set of daylight variables as described in (7) for each type of day. The categories considered after testing other combinations are Mondays, Sundays and holidays, Saturdays, and the rest of the weekdays.

Description of Tests
As it was aforementioned, the aim of the paper is to provide a methodology that improves forecasting accuracy through the use of available daylight information. In order to assess the different possibilities considered, several changes have been made from the original forecasting model in the design and testing stages.
The metric used to compare results is the mean average percentage error (MAPE) as described in (9). The simulation conditions are day-ahead predictions cast at 10 a.m. the day before. For the design phase, the data used is from 2011 to 2017, and the model used corresponds to the AR model stripped from its autoregressive components as described in (1). In order to obtain out-of-sample results, each year is simulated using a model estimated with the rest of the years. In that way, we obtain a 7-year-long period of results. The design phase allows us to determine the best model, which is then tested in the assessment phase.
where MAPE d is the MAPE for day d, F d,h is the forecast for day d at hour h and L d,h is the actual load on day d and hour h.
In the assessment phase, both the AR and neural network (NN) models are used, including all of their designed features. The testing data is a 1-year-long period of 2018, and the model is estimated with the previous data (2011-2017).

Results
This section presents the results not only for the final model but also for the design phase, as a justification of the design decisions and a guide to the application of this methodology to other databases. The benchmark used to compare the results of the model is the current forecast used by the system TSO, as our proposed model is held to the same computational restrictions and deals with the same database. Comparison across techniques applied to different databases may hide differences in the predictability of each data series and may lead to wrong conclusions about the quality of each technique [3]. The computational burden of the proposed model is the same as that of the previous system.

Design Results
Six different input configurations were tested in this stage. The details of each configuration are included in Table 4. The results for the out-of-sample tests are shown in Figure 4. Table 5 also shows the results by grouping the different periods: sunrise, sunset, mid-day and mid-night.  From Figure 4 it can be seen that in the sunrise and sunset intervals both linearizations (traces 3 and 4) are essentially equivalent and that the first improvement comes from using either one of them and then, a second increase in accuracy is obtained by discriminating by type of day.
Considering that the model 5 uses 35 more variables than model 6, the piecewise linearization is discarded for the assessment phase and the type of day discrimination is found to be significant.

Assessment Results
The previous section presented the results that led to the selection of a sigmoid linearization with type-of-day discrimination due to its accuracy and low number of variables used. The errors presented were out-of-sample modeling errors. Nonetheless, the actual assessment of this paper's contribution is to compare the original forecasting model used by the TSO with the proposed one in real conditions. To this end, the AR model is restored to its original configuration and the new variables are added as input. The two configurations are compared using new data from 2018. Table 6 presents the forecasting errors of both models along with the forecasting error of the NN model also used in the TSO's ensemble. This NN model does not include any information regarding daylight other than what was already included in the original model. Table 6 also includes the results of the week after DST is implemented. This particular week is especially sensitive to the modeled phenomenon and, therefore it is a good reference value: The results show that the inclusion of available daylight information causes a reduction of the forecasting error in sunrise and sunset times of 0.12 and 0.18 percentage points respectively. As expected, midday and midnight hours remain unaffected. This accuracy improvement turns the AR model into the best performing model during sunset hours. The overall result indicates that the new AR model is nearly as accurate as the NN model. All these conclusions are even more obvious from the DST week results. The new variables have been introduced in the AR model in order to obtain a significant improvement on sunrise and sunset hours. Table 6 includes results from the original NN model as a comparison for both AR models. However, the question of whether the NN can be improved remains.
Several attempts to include available daylight information into the NN model were carried out. Considering the NN inherent ability to model non-linear behavior, the original variables of the number of hours to sunrise and sunset were included. Also, on a separate attempt, the linearized variables were also tested, yet, none of these tests yielded any improvement in accuracy of the NN model. Figure 5 shows the NN model's result on 2018 compared to the AR models. It shows how the NN model does not show the sunrise and sunset error humps that this paper addresses. This behavior could be explained through the non-linear behavior of the NN model and the high correlation between temperature (already in the model) and daylight. The fact that the AR model became so much more similar to the NN model when the available daylight information was introduced that it leads to thinking that while the NN was able to model this behavior from its original variables, the AR model needed the linearization provided to increase its accuracy. The final forecast emitted by the system is a linear combination of both the AR and NN models. Table 7 shows how the weight in the combination has shifted towards the AR model and the final improvement of the overall system. Such combination is the result of an algorithm implemented to minimize error based on the last 30 days. Further details of this process can be found in [6]. The results of Table 7 include other references as context for the provided results. All the references used refer to the same electric system because comparing accuracy using different databases can be misleading, as some systems are more predictable than others. Nevertheless, the proposed model presents a clear advantage during sunrise and sunset hours over all reported results.
All these results were obtained on an i7 2.2GHz PC with 16Gb of RAM under real-time conditions. A ten-day forecast is emitted every hour with a computing time of 5.3 min on average that includes reading input files, processing the forecast and writing output files.

Conclusions
The possible inputs for STLF models include previous loads, calendar variables, and environmental variables like temperature, humidity or daylight. However, even though the use of electric lighting depends directly on available daylight, its use is scarce in the literature as its own variable. This paper proposes a methodology for including available daylight information into STLF models to achieve the goal of reducing forecasting error during sunrise and sunset hours. It uses as raw variables the number of hours until sunrise and sunset, it describes how sunrise and sunset hours are affected differently and how the non-linear relation between load and the raw variables should be linearized through sigmoid functions without over-complicating the model. The methodology has been tested using out-of-sample data from 2018. The full-year results are compared with other models currently running at the National Transmission System Operator in Spain, proving that the followed approach leads to significant increase in short-term load forecasting accuracy, especially focused on sunrise (1.33% to 1.21%) and sunset (1.56% to 1.38%). The improvement is even more relevant on the weeks in which DST is implemented and daylight changes suddenly: sunrise error falls from 1.90% to 1.37% and sunset error drops from 2.53% to 2.09%. The inclusion of this variables did not cause any improvements in the forecasting model based on neural networks. However, the linear combination of both models used reflected the sunrise and sunset improvement (1.09% vs. 1.12% and 1.30% vs. 1.37%, respectively) and increased the weight of the new AR model (from 37% to 59%).
The methodology described in this paper will be rolled out into the TSO's forecasting system in 2020 as it will help reduce forecasting error especially at peak hours like sunset, in which inaccuracies are particularly costly to the system.

Data Availability Statement:
Restrictions apply to the availability of these data. Data was obtained from REE and AEMET and are available at https://www.esios.ree.es/es and http://www.aemet.es/ es/datos_abiertos.