Frugal day-ahead forecasting of multiple local electricity loads by aggregating adaptive models

This paper focuses on day-ahead electricity load forecasting for substations of the distribution network in France; therefore, the corresponding problem lies between the instability of a single consumption and the stability of a countrywide total demand. Moreover, this problem requires to forecast the loads of over one thousand substations; consequently, it belongs to the field of multiple time series forecasting. To that end, the paper applies an adaptive methodology that provided excellent results at a national scale; the idea is to combine generalized additive models with state-space representations. However, extending this methodology to the prediction of over a thousand time series raises a computational issue. It is solved by developing a frugal variant that reduces the number of estimated parameters: forecasting models are estimated only for a few time series and transfer learning is achieved by relying on aggregation of experts. This approach yields a reduction of computational needs and their associated emissions. Several variants are built, corresponding to different levels of parameter transfer, to find the best trade-off between accuracy and frugality. The selected method achieves competitive results compared to individual models. Finally, the paper highlights the interpretability of the models, which is important for operational applications.


Introduction
Electricity consumption forecasting is essential for numerous activities: managing the electricity network, investment and production planning, trading on the electricity markets, and reducing power wastage.To suit these activities, forecasts are performed at different horizons, from short-term (hours, days) to long-term (years), and on different scales, from individual to national.New variability in the electricity load has recently emerged due to several factors: new decentralized production units (like solar panels), new actors with the opening of the French electricity market, new uses (like electric mobility), and the COVID-19 pandemic.They bring essential changes in the electricity consumption, and we argue that forecasting models must be updated to take them into account.
Numerous forecasting approaches have been proposed in the past before this new variability, and we refer to the Global Energy Forecasting Competitions (GEFCom) for an overview [Hong et al., 2014[Hong et al., , 2016[Hong et al., , 2019]].Classical time series methods such as auto-regressive integrated moving-average (ARIMA, Huang and Shih [2003], Chodakowska et al. [2021]) and exponential smoothing [Jalil et al., 2013] have been used to forecast at the very-short term (hours-ahead).Machine learning approaches provide better forecasts using explanatory variables, especially for larger horizons.Indeed, judging that electricity demand is a human activity, it can be predicted using explanatory variables, the most important ones being weather and calendar data.In particular, Generalized Additive Models (GAMs, Wood [2017]) have been widely applied to forecast the electricity consumption [Pierrot and Goude, 2011, Goude et al., 2014, Fasiolo et al., 2021].
More recently, a state-space approach has been investigated to adapt GAMs over time [ de Vilmarest, 2022a]; this online model allows GAMs to adapt to new variabilities, such as the COVID-19 pandemic, to improve forecasts [Obst et al., 2021, de Vilmarest andGoude, 2022].Finally, the variety of models motivates predicting the demand with a combination of forecasters, yielding a final prediction better than any individual model; that is the aggregation of experts [Cesa-Bianchi and Lugosi, 2006], which has also shown good results on load forecasting [Gaillard and Goude, 2015].Aggregation of experts is an online method which combines predictions given by experts in a weighed sum where the weights evolve over time.
In this paper, we study day-ahead electricity load forecasting on the distribution network in France.The electricity consumption is measured on about 2,200 substations located at the frontier between the high voltage grid and the distribution network.More precisely, we have access to the data of 1,344 of them.We thus consider forecasting at a local scale.The high number of substations motivates the need for a methodology to forecast a large number of time series efficiently; we look for the right trade-off between good accuracy and frugality in the number of parameters and computational complexity.Frugal machine learning can be motivated by economic reasons, for instance to avoid costly data collection or transmission, and by environmental goals like the reduction of computers energy consumption.An example where frugality is key is the implementation of machine learning algorithms in portable devices such as smartwatches [Evchenko et al., 2021]: transmission of data is difficult and intensive use of machine learning algorithms is limited by computational and battery saving constraints.In our case, the limitation comes from the learning process which becomes infeasible in a reasonable time with such a large number of time series to predict.
We present this challenge from a transfer learning point of view that we represent with diagrams in Fig. 1.We train individual models on a few time series, relying on GAMs and their adaptive variant based on state-space models.Then, we transfer these models learned on a small fraction of the time series to all of them, instead of training an individual model for each time series.We use aggregation of experts as a transfer tool.Precisely, we apply our base models on each time series even if the models were not trained on them, and then aggregating these weak forecasters yields a procedure able to scale to a large number of time series.We build several aggregations involving different kinds of models, where transfer learning occurs at multiple levels.We discuss the computational cost for each aggregation type and focus on their frugal nature.Indeed, optimization of individual models to a large number of time series is very costly, especially the performing models we want to use, while aggregation of experts is not.We show that our methods compete with individual models that we detail, even during the first French lockdown due to COVID-19, while achieving low computational cost.Transfer learning through aggregation of experts has also been studied in a hierarchical setting, using regional load data as sources and national load data as target [Gaucher et al., 2021].
Following this introduction, Section 2 establishes the methods used in our work.In Section 3 we introduce the data and the models applied.We analyze and present in Section 4 the results of the forecasting tasks.Finally, we conclude the paper and propose future work.

Methodology
In this section we first characterize our problem in the transfer learning framework.The transfer is obtained by aggregation of experts, and we describe its principle.Finally, we introduce the different models we consider and the resulting aggregation methods.

Transfer learning context
We start by expressing the transfer learning context using the definitions and vocabulary used for transfer learning methods [Zhuang et al., 2020].First, we define the notion of domain and task.A domain D = {X , P (X)} is composed of a feature space X and a marginal distribution P (X) of an instance set X = {x | x i ∈ X , i = 1, . . ., D}.A task T = {Y, f } is composed of a label space Y and a Figure 1: Transfer learning through aggregation of experts in two steps : (a) we train individual models for E time series for a small E compared to the number of time series, then (b) to forecast any time series k, we apply models M 1 , M 2 , . . ., M E and we combine all the resulting forecasts using aggregation of experts.decision function f .Given m S ∈ N + source domains and tasks, and target domain and task, transfer learning is the use of the knowledge within the source domains and tasks to improve the performance of the target decision function.
Each substation represents one target whose decision function varies from the other substations' decision functions.Thus, we have m T = 1, 344 transfer learning tasks, each corresponding to one forecasting task.It results in a multi-task transfer learning situation.On the other hand, judging that we will use several sources to do the transfer, we have m S > 1, resulting in a multi-source transfer learning method.
In our case, we are in a homogeneous and inductive transfer learning situation.Indeed, a homogeneous transfer learning scenario occurs when the source and target feature spaces are equal, as well as the source and target label spaces.The explanatory variables (calendar and meteorological variables and past electricity load) are equally available for all substations and range thus in the same feature space.The mutual label space is R + corresponding to the electricity measurement range.Concerning inductive transfer learning, it happens when labeled data are available in the target domain to induce the target decision functions, which is our case.

Transfer learning using aggregation of experts
A natural idea would be to determine a subsample of m S data sets, build one individual model per data set, and select the best among the m S models for each substation we want to forecast.The information on the best model is unavailable; one way to estimate it is through the aggregation of experts.
In the aggregation of experts context [Cesa-Bianchi and Lugosi, 2006], we want to predict a bounded sequence of observations y 1 , . . ., y n ∈ [0, B] (B is unknown) using E forecasting models called experts.For each time step t ∈ {1, . . ., n}, they provide E forecasts (ŷ e t ) E e=1 of the observation y t .The aggregation ŷt = E e=1 pe,t ŷe t is then computed where the weights (p e,t ) E e=1 are updated online according to past performances of each expert.Forecast error is measured with a convex loss function t (y t , •).The goal is to minimise the so-called regret R T = 1 T T t=1 t (y t , ŷt ) − 1 T T t=1 t (y t , ŷ t ) where ŷ t is given by an oracle model which can use unavailable information to build a forecast difficult to beat.The regret is the difference between the error suffered by the aggregation and the error of the oracle.The latter can be the best-fixed convex combination of all the experts or the best-fixed expert (constant over time).We use the ML-Poly algorithm proposed by Gaillard et al. [2014] and implemented in the R package opera [Gaillard and Goude, 2016].This algorithm tracks the best expert or the best convex combination of experts by giving more weight to an expert that will generate a low regret.This makes this algorithm particularly interesting as no parameter tuning is needed.

Adaptive experts
In this paragraph, we detail the forecasting methods we aggregate.We choose GAMs adapted by Kalman filtering as proposed by Obst et al. [2021].We consider in this section the forecasting of one time series.

Generalized additive models
We consider generalized additive models (GAM, Wood [2017]) where the response variable y t is expressed as the following sum: where β 0 is an intercept, ε t is the model error at time t, (x t,d ) D d=1 are the D explanatory variables available at time t and (f d ) D d=1 are linear or nonlinear smooth functions called GAM effects.An effect f d is expressed as a projection where (B d,k ) m d k=1 is a spline basis of dimension m d and (β d,k ) m d k=1 are the corresponding coefficients estimated by a ridge regression, where we minimize the following criterion: The penalty term controls the second derivatives f d to force the effects to be smooth.We denote C 1 as the computational cost of GAM estimation.

Adaptation by Kalman filter
To adapt a GAM, we use a multiplicative correction of the following GAM effects vector f , where f j is a normalized version of f j obtained by subtracting the mean on the train set and dividing by the standard deviation.
The adaptation is obtained assuming that a state-space property is satisfied.Precisely, we estimate a vector θ t called state under the assumption that where (ε t ) and (η t ) are Gaussian white noises of respective variance / covariance σ 2 and Q.
Starting from a Gaussian prior and assuming the variances σ 2 and Q are known, the Kalman filter [Kalman, 1960] achieves the estimation of the state θ t .This is a Bayesian method where at each step, the state posterior distribution is obtained as a Gaussian distribution θ t | (x s , y s ) s<t ∼ N ( θt , P t ).
The Kalman filter depends on the initial parameters θ1 , P 1 (the prior) and the variances σ 2 and Q.We choose these hyper-parameters by maximizing the likelihood on the training set.We use an iterative greedy procedure [ de Vilmarest, 2022a, chapter 5].We refer to this choice of hyper-parameters as the dynamic setting.Its computational cost is noted C 2 and can be large; it is detailed in the experimental study.Thanks to transfer learning, we avoid the costly estimation of these hyper-parameters for each time series.

Models considered
As previously mentioned, the estimation of the GAM and of the Kalman variances is computationally costly when applied to 1,344 time series.We neglect the computational time due to Kalman updates, and we propose to reduce the number of GAMs and Kalman variances that are estimated.We define hybrid experts M i,j,k , where the different parts of the forecaster are trained on different time series.GAM is trained on data set i, Kalman variances optimized (in the dynamic setting) on data set j, and resulting adapted model applied on data set k, where i, j, and k range in {1, . . ., m T }.We denote M i,∅,k a non-adapted GAM and M i,0,k an adapted GAM in the static setting where no variances are optimized.
There are three basic models without transfer learning where the models are optimized on the same substation's data we forecast: M k,k,k , M k,0,k , and M k,∅,k .It is the scenario of traditional machine learning.We also consider three models involving transfer learning.The most simple transfer is GAM transfer: a source data set is used to train the GAM, and the resulting GAM is applied to a different target data set.It corresponds to the model M i,∅,k , i = k.An immediate improvement of the previous model is the adaptation of the transferred GAM with a Kalman filter optimized on the same data set used to train the GAM.It results in models M i,i,k , i = k.The adaptation step helps the GAM transfer and improves the basic GAM.Finally, we consider the Kalman filter transfer case: the data set we want to forecast is used to train a GAM adapted by a Kalman filter optimized on another data set.These are models M k,j,k , j = k.

Aggregation models
As previously said, the estimation of Kalman variances in the dynamic setting is computationally costly (e.g.training of M k,k,k model), and we want to avoid individual training for all time series.Moreover, although the computational cost of GAM effects estimation is smaller, we possibly want to train only a few of them to reduce the cost and the number of model parameters.To do so, we build three aggregation methods based on the three previous models where transfer learning occurs, provided in Section 2.3.3.We obtain AGG GAM TL from the aggregation of n 1 M i,∅,k models, AGG GAM-Kalman TL from n 2 M i,i,k models, and AGG Kalman TL from n 3 M k,j,k models.The data sets used to train GAMs and Kalman filters in each aggregation method are randomly chosen among the 1,344 data sets.
Table 1 details the three aggregation methods, and the two individual adapted GAMs.It specifies if a transfer is involved, the type of the model used, and the computational costs corresponding to GAM and Kalman variances estimation.More precisely, the computational costs correspond to GAM and Kalman variances estimations.Judging that GAM, Kalman filter, and aggregation applications are computationally cheap, we overlook them.For each aggregation method, parameter n i is the unique hyper-parameter of the method and corresponds to its number of sources m S .It must be as small as possible compared to m T = 1344 to force the frugality of the method.
Table 1: Description of the individual adapted models and aggregation methods considered in the case study.C 1 and C 2 correspond to the computational costs of GAM and Kalman variances estimation, respectively.m T is the number of targets in the transfer learning context, and n i is the number of experts in the three aggregation methods.To give an order of magnitude, in our application, parameters n i are inferior to 10 while m T = 1, 344.

Data and model presentation
We first introduce French local electricity consumption and explanatory variables.We then detail the models in the practical point of view: definition of GAMs formula, training of Kalman filter and construction of the aggregation methods.To guarantee data confidentiality, substations' identities are not provided, and the electricity loads represented in the different figures are normalized by the average load.

Electricity load data
The data are provided by Enedis, the operator in charge of the electric power distribution in France.
The data are composed of 1,344 time series, each of which is the electricity consumption of one substation represented in blue in Fig. 2. They cover metropolitan France and reflect thus its local electricity consumption.Forecasting each substation's consumption involves m T = 1, 344 forecasting tasks.The data are available from June 1 st , 2014, to December 31 st , 2021, with a 30-minute temporal resolution.We forecast electricity consumption for the next day with information up to the previous day; that's why we fix an operational constraint on data availability with a lag of 48 hours.Although most of the 1,344 time series present classical temporal and meteorological patterns, there are some counter-intuitive and contrary variations.We refer to Fig. 3 for comparing one substation with basic behavior and two substations with unusual behaviors.

Explanatory variables
For explanatory variables, we choose meteorological and calendar variables, a common practice in electricity load forecasting.Weather variables are obtained from the French weather forecaster Météo-France.They are composed of temperature, cloudiness, and wind measurements from 27 weather stations.These stations are unequally spread over metropolitan France and are represented in red in Fig. 2. Weather data are available in the same range as electricity consumption with a 3 hours temporal resolution.We transform these data into 30 minutes of temporal resolution thanks to linear interpolation.Temperature is highly correlated to electricity load with a different impact for cold and hot regions, as shown in Fig. 4. The two patterns are due to the use of electrical heating and air-conditioning systems in France, which are highly electric consuming.On the other hand, calendar variables gather indicators of holidays by region, bank holidays, and working days, as well as the instant of the day, time of year, and days of the week.

Time Segmentation
The data have an availability range of 7.5 years.We train and validate the models on the data up to December 31 st , 2019, and we test afterward.The test set includes COVID-19 observations and three lockdown periods in France.We propose to consider the second and third lockdowns as normal periods because the electricity consumption strongly varies only during the first lockdown, see Fig.

Generalized Additive Model formula
To build a GAM, we must determine its formula: which explanatory variables and spline bases to use.To do so, 2018 is used as a test period, and the different tests have been achieved with a forwardbackward heuristic.We assume that a single GAM formula could be used for all target forecasting tasks.Indeed, although substations' behaviors are different, they can be explained by the same explanatory variables.Finally, we apply a GAM by instant of the day; that is, we deal with 48 GAMs for each forecasting task.We get the unique following GAM formula: where at each time step t, • y t is the electricity load for the considered instant, • DayT ype t is a categorical variable indicating the type of the day.There are five categories: Monday, Tuesday to Thursday, Friday, Saturday, and Sunday, • BankHoliday t is a binary variable indicating whether the day t is a bank holiday or a vacation day, • T oY t is the time of year whose values grow linearly from 0 on the 1st of January midnight to 1 on the 31st of December 23h30, • W orkingDay t is a binary variable indicating whether the day t is a working day, i.e., not a weekend day or a bank holiday, • T rend t is the number of the current observation, • T emp t is the temperature of the closest weather station, • T emp95 t and T emp99 t are exponentially smoothed T emp t variable of factor α = 0.95 and 0.99.E.g. for α = 0.95 at a given time step t, T emp95 t = αT emp95 t−1 + (1 − α)T emp t , • T empM in t and T empM ax t are the minimal and maximal value of T emp t at the current day, • Load2D t and Load1W t are loads of 2 days before and the load of the week before, • ε t is Gaussian noise with 0 mean and constant variance.
This model is a short-term model in terms of electricity past consumption availability.It is thus called ST GAM.We also consider a model in a mid-term context which is the same model without the Load2D t and Load1W t effects.It is called MT GAM.
We also examine the case where ε t is an auto-correlated error term.In that case, we use an ARIMA (autoregressive integrated moving average) model by selecting the best model with AIC criterion in the family of ARIMA(p,d,q).Correcting residuals of mid and short-term GAM with an ARIMA model achieves the same performance; therefore, we choose the more straightforward mid-term formula.It can be explained by the redundancy of correcting auto-correlated residuals with an ARIMA model and the linear lag terms in the short-term GAM.We call the following model MT GAM + ARIMA.
Thin plate spline basis with low dimensions represent all the effects except f 1 and f 2 .Indeed, time of year has a cyclic impact (see Fig. 4 (a)); therefore, we use a cyclic cubic splines basis of dimension 20.GAM effects estimation is done using the Generalized Cross Validation criterion and takes a few tens of seconds in practice.Thus, it is computationally reasonable to train individual GAM for many forecasting tasks.Finally, we train the GAMs using R and the package mgcv [Wood, 2015].

Kalman filtering
We examine in this section Kalman filtering adaptation of the previous short-term GAM.In the static setting, we can apply one individual adapted GAM to each forecasting task as no Kalman variances estimation is necessary.It is called GAM + Kalman Static.In this case, we are in a traditional machine learning context where the model is optimized on the same data set we forecast.
Execution time is essential in the dynamic setting.This is an unavoidable operational constraint, and we choose to calculate Kalman variances on only a subsample of the 1,344 substations.To select this subsample, we first use the short-term GAM predictions for 2018 to sort the 1,344 substations by forecast performance.The 44 substations corresponding to the 44 worst performances are put aside as we assume they won't provide interesting experts.We then randomly draw 6 substations in groups of 100 substations.We finally get a subsample of 78 substations representative of the forecasting difficulty.The parallel computing of the 78 corresponding sets of Kalman variances on a 36-core virtual machine on 2014-2018 data takes about 1.6 days.The computation of 1,344 individual Kalman variances would thus last about 28 days in a similar setting.
We calculate two sets of Kalman variances on the same substations subsample: one on 2014-2018 data and the other on 2014-2019 data.The first set is used to set the aggregation hyper-parameters n i with 2019 as the validation set, and the second set is used to forecast 2020-2021 as the test set.We thus have at most 78 GAMs adapted by their individual Kalman variances.In that case, we call the corresponding model GAM + Kalman Dynamic, and we show that the aggregation methods achieve equivalent performances.
Kalman filters estimation and application are achieved using R and the package viking [de Vilmarest, 2022b].

Aggregation
We refer to the Aggregation models section and Table 1 for a remainder on the three aggregation methods.We focus here on their application.
For each aggregation method, the unique hyper-parameter is n i , the number of experts in the aggregation, which is also the number of sources in the transfer learning task.As said previously, we want m S to be the smallest possible compared to m T = 1, 344.We proceed with a grid search on 2019 data as the validation set.We perform 10 forecasts of 182 representative substations for different values of n i and compute the 10 corresponding medians.The 182 representative substations are obtained in the same way we obtain the 78 substations used to optimize Kalman variances, but here we pick 14 substations by groups of 100.We represent results in Fig. 6.We can see that aggregation is gradually robust under the randomness of substations selection to compute GAM effects and Kalman variances.Moreover, there is an elbow phenomenon after n 1 = n 2 = 9 for AGG GAM TL and AGG GAM-Kalman TL, and after n 3 = 6 for AGG Kalman TL.These values of n i are very little compared to m T = 1, 344.

Experiments
In this section, we first provide model analysis thanks to visualization plots and we then compare performances on substations data sets thanks to a comparative metric.

Model dynamics
We first analyze the aggregation methods at the expert scale and then at the GAM effect scale.Concerning aggregations of experts, we want to know which expert is important in the mixture and when.To do so, we study the distribution of the weights of each expert.Indeed, the higher the weight, the more important the corresponding expert is in the forecast.One expert can be important in the forecast for specific observations and not for others.To give an example, we represent in Fig. 7 boxplots of the weights of the 1,344 AGG GAM-Kalman TL aggregations for three interesting instants of the day where the consumption behaviors are very different.Experts 3, 4, 5, and 6 are essential in the same way according to the three instants of the day, with medians near the uniform weight of 1/9.On the other hand, the other experts show variations according to the instant of the day.We can see that, for each instant, each expert is useless at least once (weight close to 0) and of great importance at least once (weight close to 0.5 or bigger).
We can focus on each forecast to see which experts are used and when in an aggregation.We display in Fig. 8 the evolution of weights from AGG GAM-Kalman TL aggregation.These figures represent the evolution of weights according to midday observations of 2 substations.We observe the same previous observations: some experts are called very little, while others are important.The new information these figures provide is the evolution of experts' importance during aggregation.For  example, in Fig. 8 (a), Exp4's weight increases and becomes of paramount importance, with the weight being superior to 1/2 for the last observations.In Fig. 8 (b), Exp7 vanishes rapidly while the impact of Exp1 increases.Experts 5 and 9 contribute very poorly to the first observations and are quickly absent in the aggregation.
We also analyze the aggregation methods at the effect scale.To that end, we visualize the evolution of state coefficients from Kalman filtering; precisely, we plot the D curves corresponding to the D GAM effects.We provide an example in Fig. 9 (a) with the evolution of state coefficients of an adapted GAM M k,k,k during the first lockdown.That example is representative of what we observe for other state coefficients' evolution.Let's focus first on the coefficients associated with Bias and f 2 (T oY ), which is the effect modeling the time of year during working days.Their coefficients are more significant compared to others and evolve rapidly.Concerning the bias coefficients, we think they evolve to balance the gap between the source load used to train the GAM and the target load, which can be of different scales.Then some state coefficients evolve significantly, like BankHoliday or DayT ype, while others seem time-invariant.
If we select AGG Kalman TL, we can express the aggregation as a GAM adapted by hybrid state coefficients ( where at time t we denote θ t,d = E e=1 pe,t θt,d,e the hybrid state coefficients, E the number of experts, D the number of GAM effects, pe,t the weight associated with expert e, and θt,d,e the adaptation coefficient associated with the GAM effect f d of expert e.This new expression allows for gathering information on both adaptation and aggregation and improving the interpretation of the whole model.We represent in Fig. 9 (b) the evolution of these hybrid state coefficients during the first lockdown for the same target forecasting task as in Fig. 9  n 3 M k,j,k models coefficients.We see that both types of coefficients differ in amplitude and dynamic.Bias coefficients are still important compared to the others and evolve dynamically, but it's not the case for f 2 coefficients.Moreover, most of the hybrid state coefficients evolve contrary to the previous ones.Finally, we can make the same observations of this example for the other substations.
If we consider now any aggregation method, we can express the whole mixture as a GAM using the linearity of GAM effects, GAMs, the adaptation procedure, and the aggregation of experts.The new coefficients of the spline basis are composed of the basic coefficients of the spline basis, the adaptation coefficients, and the aggregation weights.

Numerical results
To compute numerical performances of the models, we use the normalized mean absolute error (NMAE) defined as It is the mean absolute error between the ground truth y and the forecast ŷ normalized by the absolute mean of y.This metric allows us to compare forecasts of different scales, which is necessary in our case study.We choose the NMAE rather than the mean absolute percentage error because our time series are sometimes very close to 0. More precisely, we use the percentage of NMAE by multiplying it by 100.Table 2 provides the numerical performances of the 1,344 forecasting tasks.For each model, the first, second, and third quartiles of the 1,344 performances are provided during the three validation periods.Kalman variances and GAMs have been trained on 2014-2019 for each test period.We display with the same setting in Table 3 the numerical performances of the forecasting tasks corresponding to the substations subsample chosen to optimize the set of Kalman variances.
We start by comparing the performances of the individual and aggregation models displayed in Table 2.Each time we cite improvements in the performances, we refer to the three respective validation periods.We can see that the information on the past electricity load added in the GAM formula is of great importance, with an improvement between ST GAM and MT GAM median performances of 1.84%, 2.83%, and 3.13%.There is also an improvement when MT GAM residuals are modeled with an ARIMA model: there is an improvement between MT GAM and MT GAM + ARIMA median performances of 0.77%, 1.46%, and 1.70%.GAM + Kalman Static provides slightly weaker performances than MT GAM + ARIMA, and the latter is thus the individual benchmark to beat.
Concerning the aggregation methods, AGG GAM TL shows the worst performance.However, this model is still interesting as its performances approach the ST GAM performances while its computational cost is much lower.On the other hand, the two other aggregation methods show excellent and close performances.AGG Kalman TL is better than AGG GAM-Kalman TL outside of the first lockdown (improvement in the median of 0.03% and 0.07%), while it is the opposite for the first lockdown (decrease of 0.16%).The gap between these two last models is not significant.As AGG GAM-Kalman TL is less complex than AGG Kalman TL, we consider it the best compromise between forecasting performance and computational cost.Finally, compared to MT GAM + ARIMA, its median performances improve by 0.62%, 1.54%, and 0.84%.We represent its 1,344 NMAE scores Fig. 10.
As explained in Section 2.2, aggregation aim at beating two oracles: the best-fixed expert and the best-fixed convex combination of experts.We call them Expert Oracle and Convex Oracle, respectively, and we compute their for AGG GAM-Kalman TL.The aggregation competes with the best-fixed expert and the best-fixed combination of experts.We observe the same for the two other aggregation Concerning the test periods, we can see that the 2020 forecasts worsen during the first French lockdown with a downgrade of 2.18% the two MT GAM + ARIMA median performances.AGG GAM-Kalman TL and AGG Kalman TL lower downgrades: 1.26% and 1.45% in the median, respectively; that is, the online nature of Kalman filtering and aggregation of experts allow adaptation to the extreme change in electricity consumption.On the other hand, ST GAM and MT GAM performances worsen in 2021 as these models are not updated while the consumption behavior evolves.The downgrades are significant: 1.37% and 2.36% in the median, respectively.The performance gaps between 2020 out of the first lockdown and 2021 are tighter for adaptative models, especially for the two best aggregation methods: 0.22% and 0.18% in the median, respectively.Once again, it shows that the aggregation methods adapt well to the data distribution evolution.
Finally, we provide in Table 3 performances of 69 forecasting tasks.Precisely, they correspond to the 78 substations we used to optimize the set of Kalman variances out of the 9 substations data sets selected to train the experts involved in the three aggregation methods.That is we avoid the lucky case where Kalman variances are optimized on a data set and used to forecast this data set.We can thus compare the aggregation methods with GAM + Kalman Dynamic in the most common case.As previously said, GAM + Kalman Dynamic achieves great performance forecasting the national electricity load and is thus a model we want to compete with.We can see that it is the case for AGG GAM-Kalman TL and AGG Kalman TL; therefore, these 2 methods are accurate comparing to the state-of-the-art model.

Conclusion
In this paper, we propose a frugal method to forecast multiple electricity loads.We avoid the estimation of an individual model per time series and we obtain a scalable methodology with limited resources based on aggregation of experts.The chosen experts are GAMs adapted by Kalman filtering, which are state-of-the-art models in electricity consumption forecasting.The aggregation methods allow the transfer of GAMs and Kalman filters separately or simultaneously.We show that they provide good forecasts compared to classical individual models, especially during the first French lockdown.Moreover, they compete with individual state-of-the-art GAM adapted by Kalman filtering.On the other hand, they don't need human intervention nor expertise and are thus simple to use.The method is frugal in terms of parameter estimation, and the computational cost has been discussed extensively.
Finally, it benefits from the interpretability of GAM and the aggregation of experts.We see various ways to extend our work.First, we think we can improve the random selection of the time series used to train GAMs and Kalman variances.We may characterize clusters of substations according to some characteristics (geographic, weather, type of consumers) and select one representative per cluster.The second way of improvement is the inclusion of new explanatory variables at the same local scale.Geo-tracking and communication data reflect human behavior and are therefore helpful for electric consumption forecasting [Gaucher et al., 2021].Finally, the computational complexity of the method introduced essentially depends on the estimation of a fixed number of models; therefore, we could apply the method to a larger data set, for instance the electrical consumption at a finer granularity.

Figure 2 :
Figure 2: Location of the 1,344 substations (in blue) and 27 weather stations (in red).

Figure 3 :
Figure3: Electricity loads in March for three substations: (a) averaged daily load and (b) averaged weekly load.Substation A has the expected daily and weekly patterns: the load is ordinarily high during daylight, especially when people are at home, low during the night, and similar during weekdays, with a drop during weekends.On the other hand, substation B has opposed daily behavior at night, and load is increasing during the week and maximal during weekends.Substation C shows a high load in the mid-afternoon, and load is still important at night.Concerning its weekly load, it drops Friday in addition to the weekend.

Figure 4 :
Figure 4: Dependence of 2015 electricity consumption at 6 p. m. to (a) time of year and (b) temperature.Blue values for a substation in the north of France and red values for one in the south of France.

Figure 5 :
Figure 5: Weekly moving average of the electricity load of one substation in 2019, 2020, and 2021.

Figure 6 :
Figure 6: Grid search of the number of experts: (a) n 1 for AGG GAM TL, (b) n 2 for AGG GAM-Kalman TL, and (c) n 3 for AGG Kalman TL.

Figure 7 :
Figure 7: Boxplot representation of the weights of the 1,344 AGG GAM-Kalman TL aggregations for three instants of the day.2021 target forecasting tasks.
Figure 8: Weights associated with the experts of AGG GAM-Kalman TL aggregation at midday for (a) substation D and (b) substation E. Experts are denoted Expi, and the forecasting period is 2020 out of the first lockdown.

Figure 9 :
Figure 9: Evolution of (a) state coefficients of a Kalman TL model ( θt,d − θ0,d ) D d=1 and (b) hybrid state coefficients of a AGG Kalman TL aggregation ( θ t,d − θ 0,d ) D d=1 , during the first lockdown and for the same target forecasting task.

Figure 10 :
Figure 10: Sorted performances of AGG Kalman TL for the 1,344 forecasting tasks in 2021.