# Probabilistic Load Forecasting Based on Adaptive Online Learning

Verónica Álvarez, Santiago Mazuelas, and José A. Lozano Manuscript received
XXXX, YYYY, and revised AAAA, BBBB. This research was supported in part by the
Spanish Ministry of Economy and Competitiveness MINECO under Ramon y Cajal
Grant RYC-2016-19383, Project PID2019-105058GA-I00, and Project
TIN2017-82626-R, the BCAM’s Severo Ochoa Excellence Accreditation
SEV-2017-0718, the Basque Government through the ELKARTEK and BERC 2018-2021
programmes, Basque Government under the grant “Artificial Intelligence in BCAM
number EXP. 2019/00432”, and the Iberdrola Foundation under the 2019 Research
Grants. V. Álvarez is with the BCAM-Basque Center for Applied Mathematics,
Bilbao 48009, Spain (e-mail: valvarez@bcamath.org). S. Mazuelas is with the
BCAM-Basque Center for Applied Mathematics, and IKERBASQUE-Basque Foundation
for Science, Bilbao 48009, Spain (e-mail: smazuelas@bcamath.org). J.A. Lozano
is with the BCAM-Basque Center for Applied Mathematics, Bilbao 48009, Spain
(e-mail: jlozano@bcamath.org).

###### Abstract

Load forecasting is crucial for multiple energy management tasks such as
scheduling generation capacity, planning supply and demand, and minimizing
energy trade costs. Such relevance has increased even more in recent years due
to the integration of renewable energies, electric cars, and microgrids.
Conventional load forecasting techniques obtain single-value load forecasts by
exploiting consumption patterns of past load demand. However, such techniques
cannot assess intrinsic uncertainties in load demand, and cannot capture
dynamic changes in consumption patterns. To address these problems, this paper
presents a method for probabilistic load forecasting based on the adaptive
online learning of hidden Markov models. We propose learning and forecasting
techniques with theoretical guarantees, and experimentally assess their
performance in multiple scenarios. In particular, we develop adaptive online
learning techniques that update model parameters recursively, and sequential
prediction techniques that obtain probabilistic forecasts using the most
recent parameters. The performance of the method is evaluated using multiple
datasets corresponding with regions that have different sizes and display
assorted time-varying consumption patterns. The results show that the proposed
method can significantly improve the performance of existing techniques for a
wide range of scenarios.

###### Index Terms:

Load forecasting, probabilistic load forecasting, online learning, hidden
Markov model.

##  I Introduction

Load forecasting is crucial for multiple energy management tasks such as
scheduling generation capacity, planning supply and demand, and minimizing
energy trade costs [1, 2, 3, 4]. The importance of load forecasting is
significantly expanding in recent years due to the increasing development of
power systems and smart grids [4]. In addition, accurate load forecasting has
a beneficial impact in environment and economy by reducing energy waste and
purchase [5].

Forecasting methods are enabled by exploiting consumption patterns related to
multiple factors such as past loads, hours of day, days of week, holidays, and
temperatures [6, 7, 8, 9, 4]. Accurate forecasting is hindered by intrinsic
uncertainties in load demand and dynamic changes in consumption patterns [10,
11, 12, 13]. These problems are becoming more relevant in recent years due to
the integration of renewable energies, electric cars, and microgrids [14, 15,
16, 17, 18]. Uncertainties in load demand cannot be assessed by methods that
obtain single-value forecasts, and dynamic changes in consumption patterns
cannot be captured by methods based on offline learning of static models.
Methods that obtain probabilistic forecasts are essential for optimal
stochastic decision making (e.g., unit commitment [19]) [20]. Methods based on
online learning of adaptive models are necessary to harness dynamic changes in
consumption patterns [21]. Figure 1 illustrates how changes in consumption
patterns affect the performance of methods based on offline and online
learning. The top figure shows the two-peak consumption pattern that both
methods exploit to train their models, and the middle subfigure shows how both
methods accurately forecast until a flatter pattern emerges in day $t_{0}+3d$.
Then, methods based on offline learning cannot adapt to the new pattern, while
methods based on online learning correctly adapt to such change.

Fig. 1: The offline learning method cannot capture the change from two-peak to
flatter pattern at day $t_{0}+3d$, while the online learning method can
harness such change.

Most conventional techniques for load forecasting obtain single-value
forecasts based on offline learning. Such techniques can be classified in
three main groups: techniques based on statistical methods (e.g., linear
regression (LR) [9], autoregressive moving average (ARMA) [22, 23],
autoregressive integrated moving average (ARIMA) [24], and seasonal
autoregressive integrated moving average (SARIMA) [25, 26]), techniques based
on artificial intelligence (e.g., deep learning [27], neural networks [28],
and support vector machine (SVM) [29, 9]), and techniques based on weighted
combinations of several forecasts [30, 31, 32]. Existing techniques that
obtain probabilistic forecasts are based on offline learning, while those
based on online learning obtain single-value forecasts. Current probabilistic
methods are based on Gaussian process (GP) [33] and quantile regression (QR)
[34, 35]. Online learning methods adapt to dynamic changes in consumption
patterns by adjusting offline learning algorithms. In particular, such methods
retrain regularly the static models of conventional techniques such as ARMA
[36], update the weights in combined forecasts [37], or update the functions
in the additive model [38].

In this paper, we present techniques for adaptive probabilistic load
forecasting (APLF) that can harness changes in consumption patterns and assess
consumption uncertainties. Specifically, the main contributions of the paper
are as follows:

  * •

We model the data using hidden Markov models and develop online learning
techniques for APLF that update HMM parameters recursively.

  * •

We develop sequential prediction techniques for APLF that obtain probabilistic
forecasts using the most recent HMM parameters.

  * •

We describe in detail the efficient implementation of the steps for online
learning and probabilistic prediction in APLF method.

  * •

We quantify the performance improvement provided by the methods presented in
comparison with existing load forecasting techniques under multiple scenarios.

The rest of this paper is organized as follows. Section II describes the
problem of load forecasting and introduces the performance metrics. In Section
III, we present the theoretical results for APLF using HMMs. Section IV
compares procedures of offline learning and online learning, and describes in
detail the implementation of APLF. The performance of APLF and existing
techniques is compared in Section V under multiple scenarios. Finally, Section
VI draws the conclusions.

### Notations

$N\left(x;\mu,\sigma\right)$ denotes the Gaussian density function of the
variable $x$ with mean $\mu$ and standard deviation $\sigma$;
$p\left(x|y\right)$ denotes the probability of variable $x$ given variable
$y$; $p(x,y)$ denotes the joint probability of variables $x$ and $y$;
$\mathds{1}\\{\cdot\\}$ denotes the indicator function; bold lowercase letters
represent vectors; bold capital letters represent matrices; $\mathbf{I}_{K}$
denotes the $K\times K$ identity matrix; $\mathbf{0}_{K}$ denotes the zero
vector of length $K$; $[\cdot]$ denotes vectors; and $[\cdot]^{\text{T}}$ and
$\mathbb{E}\\{\cdot\\}$ denote the transpose and expectation of its argument.

##  II Problem Formulation

Load forecasting methods estimate future loads given past loads and factors
that affect future loads such as hours of day, days of week, and weather
forecasts. Forecasting techniques determine a prediction function that assigns
instance vectors $\mathbf{x}$ (predictors) to target vectors $\mathbf{y}$
(responses).

Instance vectors $\mathbf{x}$ are composed by past loads and observations
related to future loads (e.g., weather forecasts), and target vectors
$\mathbf{y}$ are composed by future loads. We denote load by $s$ and load
forecast by $\hat{s}$, with $s_{t}$ and $\hat{s}_{t}$ being the load and the
load forecast at time $t$. In addition, for each time $t$, we denote by
$\mathbf{r}_{t}$ the observations vector at time $t$ that can include data
such as weather forecasts $w_{t}$. Then, for a prediction horizon $L$ (e.g.,
$24$ hours, $30$ minutes) and prediction times $t+1,\leavevmode\nobreak\
t\leavevmode\nobreak\ +\leavevmode\nobreak\ 2,\leavevmode\nobreak\
...,\leavevmode\nobreak\ t+L$, the instance vector composed by past loads and
observations is given by $\mathrm{\mathbf{x}}\leavevmode\nobreak\
=\leavevmode\nobreak\
\left[s_{t-m},...,s_{t},\mathbf{r}_{t+1}^{\text{T}},...,\mathbf{r}_{t+L}^{\text{T}}\right]^{\text{T}}$,
the target vector composed by $L$ future loads is given by
$\mathbf{y}\leavevmode\nobreak\ =\leavevmode\nobreak\
[s_{t+1},s_{t+2},...,s_{t+L}]^{\text{T}}$, and the vector of load forecasts is
given by $\hat{\mathbf{y}}\leavevmode\nobreak\ =\leavevmode\nobreak\
[\hat{s}_{t+1},\hat{s}_{t+2},...,\hat{s}_{t+L}]^{\text{T}}$. Furthermore, each
time $t$ is categorized by a calendar variable $c(t)\leavevmode\nobreak\
\in\leavevmode\nobreak\ \left\\{1,2,...,C\right\\}$ that describes time
factors affecting load demand such as hour of day, day of week, month of year,
and holiday. The calendar variable is used to model separately loads
corresponding with each calendar type $c(t)\leavevmode\nobreak\
\in\leavevmode\nobreak\ \left\\{1,2,...,C\right\\}$ as described in Section
III. Conventional techniques such as LR [9] and SVMs [29] use instance vectors
composed by past loads, observations, and calendar variables. The proposed
APLF method uses instance vectors composed by one past load and observations.

Forecasting methods determine prediction functions using training samples
formed by pairs of vectors $\mathbf{x}$ and $\mathbf{y}$. Offline learning
algorithms determine a static prediction function $f$ using training samples
obtained up to time $t_{0}$,
$(\mathbf{x}_{1},\mathbf{y}_{1}),(\mathbf{x}_{2},\mathbf{y}_{2}),...,(\mathbf{x}_{t_{0}},\mathbf{y}_{t_{0}})$,
while online learning algorithms determine prediction functions $f_{t}$ using
all available training samples at $t\geq t_{0}$,
$(\mathbf{x}_{1},\mathbf{y}_{1}),(\mathbf{x}_{2},\mathbf{y}_{2}),...,(\mathbf{x}_{t},\mathbf{y}_{t})$.
Therefore, the static prediction function $f$ cannot adapt to changes in
consumption patterns that occur after time $t_{0}$, while prediction functions
$f_{t}$ can adapt to patterns’ changes using the latest information (see
Figure 1).

Performance of forecasting algorithms is evaluated in terms of accuracy using
the absolute value of prediction errors:

$e=\left|s-\hat{s}\right|$ (1)

while probabilistic performance can be evaluated using metrics such as pinball
losses [39] and calibration [40]. Overall prediction errors are commonly
quantified using root mean square error (RMSE) given by

$\text{RMSE}=\sqrt{\mathbb{E}\left\\{\left|s-\hat{s}\right|^{2}\right\\}}$

and mean average percentage error (MAPE) given by

$\text{MAPE}=100\cdot\mathbb{E}\left\\{\frac{\left|s-\hat{s}\right|}{s}\right\\}.$

The pinball loss of the $q$-th quantile forecast $\hat{s}^{\left(q\right)}$ is
given by

$L\big{(}s,\hat{s}^{\left(q\right)}\big{)}=\left\\{\begin{matrix}q\left(s-\hat{s}^{\left(q\right)}\right)&\text{if}&s\geq\hat{s}^{\left(q\right)}\\\
\left(1-q\right)\left(\hat{s}^{\left(q\right)}-s\right)&\text{if}&s<\hat{s}^{\left(q\right)}\end{matrix}\right.$

and the overall pinball loss is commonly quantified by the average over all
quantiles. The calibration of the $q$-th quantile forecast
$\hat{s}^{\left(q\right)}$ is given by

$C(q)=\mathbb{E}\left\\{\mathds{1}_{s\leq\hat{s}^{\left(q\right)}}\right\\}$

and quantifies the probability with which the load is less than the quantile
forecast $\hat{s}^{(q)}$. Finally, the expected calibration error (ECE) is
given by

$\text{ECE}=\mathbb{E}\left\\{|q-C(q)|\right\\}$

and quantifies the overall calibration error of probabilistic forecasts.

##  III Models and Theoretical Results

This section first describes the HMM that models loads and observations, we
then develop the techniques for online learning and probabilistic forecasting.

Fig. 2: Hidden Markov model for sequences $\\{s_{t}\\}_{t\geq 1}$ and
$\\{\mathbf{r}_{t}\\}_{t\geq 1}$ characterized by conditional distributions
$p\left(s_{t}|s_{t-1}\right)$ and
$p\left(\mathrm{\mathbf{r}}_{t}|s_{t}\right)$.

We model the relationship between the loads $\left\\{s_{t}\right\\}_{t\geq 1}$
and observations $\left\\{\mathbf{r}_{t}\right\\}_{t\geq 1}$ using HMMs also
known as state-space models [41, 42]. Such models allow to predict hidden
states from past states and observations, and are determined by the
conditional distribution $p\left(s_{t}|s_{t-1}\right)$ that represents the
relationship between two following loads, and the conditional distribution
$p\left(\mathrm{\mathbf{r}}_{t}|s_{t}\right)$ that represents the relationship
between each load and observation vector (see Figure 2). We model the sequence
of loads and observations as a non-homogeneous HMM so that both conditional
distributions change in time. Such dynamic modelling allows to adapt to
changes in consumption patterns.

The conditional distributions $p\left(s_{t}|s_{t-1}\right)$ and
$p\left(\mathbf{r}_{t}|s_{t}\right)$ are modeled using Gaussian distributions
with mean $\mathbf{u}^{\text{T}}\boldsymbol{\eta}$ and standard deviation
$\sigma$, where $\mathbf{u}$ is a known feature vector and parameters
$\boldsymbol{\eta}$, $\sigma$ are different for each calendar type
$c(t)\in\\{1,2,...,C\\}$. For each $c=c(t)$, we denote by
$\mathbf{u}_{s}^{\text{T}}\boldsymbol{\eta}_{s,c}$ and $\sigma_{s,c}$ the mean
and standard deviation that determine the conditional distribution of load at
time $t$ given load at time $t-1$ that is

$\displaystyle p\left(s_{t}|s_{t-1}\right)$
$\displaystyle=N\left(s_{t};\mathbf{u}_{s}^{\text{T}}\boldsymbol{\eta}_{s,c},\sigma_{s,c}\right)$
(2)

with $\boldsymbol{\eta}_{s,c}\in\mathbb{R}^{2}$, $\sigma_{s,c}\in\mathbb{R}$,
and $\mathbf{u}_{s}^{\text{T}}=\left[1,s_{t-1}\right]^{\text{T}}$. In
addition, for each $c=c(t)$, we denote by
$\mathbf{u}_{r}^{\text{T}}\boldsymbol{\eta}_{r,c}$ and $\sigma_{r,c}$ the mean
and standard deviation that determine the conditional distribution of load at
time $t$ given observations at time $t$. Hence, assuming there is no prior
knowledge available for the loads, we have that

$\displaystyle p\left(\mathbf{r}_{t}|s_{t}\right)$ $\displaystyle\propto
p\left(s_{t}|\mathbf{r}_{t}\right)=N\left(s_{t};\mathbf{u}_{r}^{\text{T}}\boldsymbol{\eta}_{r,c},\sigma_{r,c}\right)$
(3)

with $\boldsymbol{\eta}_{r,c}\in\mathbb{R}^{R}$, $\sigma_{r,c}\in\mathbb{R}$,
and $\mathbf{u}_{r}=u_{r}\left(\mathbf{r}_{t}\right)\in\mathbb{R}^{R}$. The
proposed method can consider general functions $u_{r}(\cdot)$ and observations
$\mathbf{r}_{t}$. In cases where the observations vector is high dimensional,
APLF method can use a function $u_{r}(\cdot)$ that reduces the dimensionality
of observations, that is, if $\mathbf{r}\in\mathbb{R}^{N}$,
$u_{r}(\mathbf{r})\in\mathbb{R}^{R}$ can be the result of applying a
dimensionality reduction method such as principal component analysis (PCA). In
the experimental results of Section V, we use a simple function $u_{r}(\cdot)$
that returns binary vectors and encodes weather shifts.

Using the above models, at each time $t$ the HMM describing the sequences of
loads and observations is characterized by parameters

$\Theta=\\{\boldsymbol{\eta}_{s,c},\sigma_{s,c},\boldsymbol{\eta}_{r,c},\sigma_{r,c}:c=1,2,...,C\\}$
(4)

where $\boldsymbol{\eta}_{s,c},\sigma_{s,c}$ characterize the conditional
distribution $p(s_{t}|s_{t-1})$ and $\boldsymbol{\eta}_{r,c},\sigma_{r,c}$
characterize the conditional distribution $p(\mathbf{r}_{t}|s_{t})$ for times
$t$ with calendar type $c=c(t)$.

The parameters for each calendar type and conditional distribution can be
obtained by maximizing the weighted log-likelihood of all the loads obtained
at times with the same calendar type. Specifically, if
$s_{t_{1}},s_{t_{2}},...,s_{t_{n}}$ are loads obtained at times with calendar
type $c\in\\{1,2,...,C\\}$, i.e., $c=c(t_{1})=c(t_{2})=...=c(t_{n})$, and
$\mathbf{u}_{t_{1}},\mathbf{u}_{t_{2}},...,\mathbf{u}_{t_{n}}$ are the
corresponding feature vectors for parameters
$\boldsymbol{\eta}_{s,c},\sigma_{s,c}$ as given in (2) or for parameters
$\boldsymbol{\eta}_{r,c},\sigma_{r,c}$ as given in (3), the exponentially
weighted log-likelihood of loads up to time $t_{i}$, for $i=1,2,...,n$, is
given by

$L_{i}\left(\boldsymbol{\eta},\sigma\right)=\underset{j=1}{\overset{i}{\sum}}\lambda^{i-j}\log
N(s_{t_{j}};\mathbf{u}_{t_{j}}^{\text{T}}\boldsymbol{\eta},\sigma)$ (5)

where weights $\lambda^{i-j}$, $j=1,2,...,i$, allow to increase the influence
of the most recent data using a parameter $\lambda\in\left(0,1\right)$ that is
commonly known as forgetting factor. The maximization of (5) is a convex
optimization problem since $L_{i}(\boldsymbol{\eta},\sigma)$ is a concave
function because Gaussian distributions are log-concave and $\lambda^{i-j}>0$
for any $i,j$. In addition, the maximum of (5) is unique as long as its
Hessian is negative definite which happens for any $i\geq i_{0}$ such that

$\mathbf{H}_{i_{0}}=\underset{j=1}{\overset{i_{0}}{\sum}}\lambda^{{i_{0}}-j}\mathbf{u}_{t_{j}}\mathbf{u}_{t_{j}}^{\text{T}}$
(6)

is a non-singular matrix.

The next Theorem shows that the maximization of the weighted log-likelihood in
(5) can be solved recursively using parameters given by

$\displaystyle\boldsymbol{\eta}_{i}$
$\displaystyle=\boldsymbol{\eta}_{i-1}+\frac{\mathbf{P}_{i-1}\mathbf{u}_{t_{i}}}{\lambda+\mathbf{u}_{t_{i}}^{\text{T}}\mathbf{P}_{i-1}\mathbf{u}_{t_{i}}}\left(s_{t_{i}}-\mathbf{u}_{t_{i}}^{\text{T}}\boldsymbol{\eta}_{i-1}\right)$
(7) $\displaystyle\sigma_{i}$
$\displaystyle=\sqrt{\sigma_{i-1}^{2}-\frac{1}{\gamma_{i}}\left(\sigma_{i-1}^{2}-\frac{\lambda\left(s_{t_{i}}-\mathbf{u}_{t_{i}}^{\text{T}}\boldsymbol{\eta}_{i-1}\right)^{2}}{\lambda+\mathbf{u}_{t_{i}}^{\text{T}}\mathbf{P}_{i-1}\mathbf{u}_{t_{i}}}\right)}$
(8)

with

$\displaystyle\mathbf{P}_{i}$
$\displaystyle=\frac{1}{\lambda}\left(\mathbf{P}_{i-1}-\frac{\mathbf{P}_{i-1}\mathbf{u}_{t_{i}}\mathbf{u}_{t_{i}}^{\text{T}}\mathbf{P}_{i-1}}{\lambda+\mathbf{u}_{t_{i}}^{\text{T}}\mathbf{P}_{i-1}\mathbf{u}_{t_{i}}}\right)$
(9) $\displaystyle\gamma_{i}$ $\displaystyle=1+\lambda\gamma_{i-1}.$ (10)

######  Theorem 1.

Let $i_{0}$ be an integer such that the matrix $\mathbf{H}_{i_{0}}$ given by
(6) is non-singular, and $\boldsymbol{\eta}_{i}^{*}\in\mathbb{R}^{K}$,
$\sigma_{i}^{*}\in\mathbb{R}$ be parameters that maximize the weighted log-
likelihood given by (5) with
$N(s_{t_{j}};\mathbf{u}_{t_{j}}^{\text{T}}\boldsymbol{\eta}_{i}^{*},\sigma_{i}^{*})\leq
M$ for any $j\leq i\leq n$.

If parameters $\boldsymbol{\eta}_{i}\in\mathbb{R}^{K}$,
$\sigma_{i}\in\mathbb{R}$ are given by recursions in (7)-(10) for $i>0$ with
$\boldsymbol{\eta}_{0}=\mathbf{0}_{K}$, any $\sigma_{0}$,
$\mathbf{P}_{0}=\mathbf{I}_{K}$, and $\gamma_{0}=0$. Then, we have that

$L_{i}\left(\boldsymbol{\eta}_{i}^{*},\sigma_{i}^{*}\right)-L_{i}\left(\boldsymbol{\eta}_{i},\sigma_{i}\right)\leq\pi
M^{2}\left\|\boldsymbol{\eta}_{i}^{*}\right\|^{2}\lambda^{i}$ (11)

for any $i\geq i_{0}$.

In addition, if parameters $\boldsymbol{\eta}_{i}\in\mathbb{R}^{K}$,
$\sigma_{i}\in\mathbb{R}$ are given by recursions in (7)-(10) for $i>i_{0}$
with

$\displaystyle\boldsymbol{\eta}_{i_{0}}$
$\displaystyle=\mathbf{P}_{i_{0}}\Bigg{(}\underset{j=1}{\overset{i_{0}}{\sum}}\lambda^{i_{0}-j}s_{t_{j}}\mathbf{u}_{t_{j}}\Bigg{)}$
(12) $\displaystyle\sigma_{i_{0}}$
$\displaystyle=\sqrt{\frac{1}{\gamma_{i_{0}}}\Bigg{(}\underset{j=1}{\overset{i_{0}}{\sum}}\lambda^{{i_{0}}-j}s_{t_{j}}^{2}-\underset{j=1}{\overset{i_{0}}{\sum}}\lambda^{i_{0}-j}s_{t_{j}}\mathbf{u}_{t_{j}}^{\text{T}}\boldsymbol{\eta}_{i_{0}}\Bigg{)}}$
(13) $\displaystyle\mathbf{P}_{i_{0}}$
$\displaystyle=(\mathbf{H}_{i_{0}})^{-1}$ (14) $\displaystyle\gamma_{i_{0}}$
$\displaystyle=\underset{j=1}{\overset{i_{0}}{\sum}}\lambda^{i_{0}-j}.$ (15)

Then, we have that $\boldsymbol{\eta}_{i}=\boldsymbol{\eta}_{i}^{*}$ and
$\sigma_{i}=\sigma_{i}^{*}$ for any $i\geq i_{0}$.

###### Proof.

See appendix A. ∎

The first part of the above result shows that parameters given by recursion
(7)-(10) for $i>0$ with $\boldsymbol{\eta}_{0}=\mathbf{0}_{K}$, any
$\sigma_{0}$, $\mathbf{P}_{0}=\mathbf{I}_{K}$, and $\gamma_{0}=0$,
approximately maximize the weighted log-likelihood. In addition, the log-
likelihood difference with respect to the maximum given by (11) decreases
exponentially fast with the number of iterations $i$ since $\lambda\in(0,1)$.
The second part of the above result shows that parameters given by recursion
(7)-(10) for $i\geq i_{0}$ with $\boldsymbol{\eta}_{i_{0}},\sigma_{i_{0}}$
given by (12)-(15), maximize the weighted log-likelihood for any $i\geq
i_{0}$. Note that the parameters are updated in recursions (7)-(10) by adding
a correction to the previous parameters $\boldsymbol{\eta}_{i-1}$ and
$\sigma_{i-1}$. Such correction is proportional to the fitting error of the
previous parameter
$s_{t_{i}}-\mathbf{u}_{t_{i}}^{\text{T}}\boldsymbol{\eta}_{i-1}$ so that
parameters are updated depending how well they fit the most recent data.

Recursion (7) for parameters describing means is similar to that used for the
recursive minimization of weighted least squares [43]. The main technical
novelty of the above result lies in recursion (8) for parameters describing
standard deviations, and inequality (11) describing the near-optimality of
parameters initialized with $\boldsymbol{\eta}_{0}=\mathbf{0}_{K}$, any
$\sigma_{0}$, $\mathbf{P}_{0}=\mathbf{I}_{K}$, and $\gamma_{0}=0$. Techniques
for least squares only allow to recursively obtain parameters describing
means, while Theorem 1 allows to recursively obtain parameters describing both
means and standard deviations. Such generalization is of practical relevance
since it allows to obtain probabilistic models determined by time-changing
means and standard deviations. In addition, existing techniques address the
possible singularity of matrix (6) during the initial steps by adding a
regularization term in (5), but such approach cannot be used to obtain
standard deviations. The bound in (11) guarantees that parameters given by
recursions (7) and (8) are close to be optimal when initialized with
$\boldsymbol{\eta}_{0}=\mathbf{0}_{K}$, any $\sigma_{0}$,
$\mathbf{P}_{0}=\mathbf{I}_{K}$, and $\gamma_{0}=0$ and are optimal when
initialized as given by (12)-(15).

Fig. 3: APLF method obtains load forecasts together with reliable uncertainty
assessments of load demand.

(a) Offline learning algorithms use $t_{0}$ training samples to train the
model that is later used to obtain all forecasts.

(b) APLF method updates the model every time a new sample is obtained. Then,
the model used for forecasting adapts to the most recent data.

Fig. 4: Block diagrams for offline learning algorithms and APLF method.

The above Theorem provides the procedure for obtaining parameters $\Theta$
described in (4). Specifically, Theorem 1 allows to update parameters
$\boldsymbol{\eta}_{s,c}$, $\sigma_{s,c}$ and $\boldsymbol{\eta}_{r,c}$,
$\sigma_{r,c}$ every time new loads and observations are obtained
corresponding with calendar type $c$. Such parameters are updated using their
previous values and the states variables given by (9) and (10). In the
following, we denote $\mathbf{P}_{s,c}$ and $\gamma_{s,c}$ (resp.
$\mathbf{P}_{r,c}$ and $\gamma_{r,c}$) the state variables required to update
parameters $\boldsymbol{\eta}_{s,c}$ and $\sigma_{s,c}$ (resp.
$\boldsymbol{\eta}_{r,c}$ and $\sigma_{r,c}$) for $c=1,2,...,C$. We denote by
$\Gamma$ the list composed by those state variables, that is

$\Gamma=\\{\mathbf{P}_{s,c},\gamma_{s,c},\mathbf{P}_{r,c},\gamma_{r,c}:c=1,2,...,C\\}.$
(16)

In addition, we denote $\lambda_{s}$ (resp. $\lambda_{r}$) the forgetting
factor required to update parameters $\boldsymbol{\eta}_{s,c}$ and
$\sigma_{s,c}$ (resp. $\boldsymbol{\eta}_{r,c}$ and $\sigma_{r,c}$), for
$c=1,2,...,C$.

The previous result describes how to update HMM parameters using the most
recent data, the next result shows how to obtain probabilistic forecasts using
the HMM characterized by parameters $\Theta$.

######  Theorem 2.

If $\left\\{s_{t},\mathbf{r}_{t}\right\\}_{t\geq 1}$ is an HMM characterized
by parameters $\Theta$ as described in (4). Then, for $i=1,2,...,L$

$p\left(s_{t+i}|s_{t},\mathbf{r}_{t+1},...,\mathbf{r}_{t+i}\right)=N\left(s_{t+i};\hat{s}_{t+i},\hat{e}_{t+i}\right)$
(17)

where means $\hat{s}_{t+1},\hat{s}_{t+2},...,\hat{s}_{t+L}$ and standard
deviations $\hat{e}_{t+1},\hat{e}_{t+2},...,\hat{e}_{t+L}$ can be computed by
the following recursions

$\displaystyle c$
$\displaystyle=c\left(t+i\right),\hat{\mathbf{u}}_{s}=[1,\hat{s}_{t+i-1}]^{\text{T}},\mathbf{u}_{r}=u_{r}\left({\mathbf{r}}_{t+i}\right)$
$\displaystyle\hat{s}_{t+i}$
$\displaystyle=\frac{\hat{\mathbf{u}}_{s}^{\text{T}}\boldsymbol{\eta}_{s,c}\sigma_{r,c}^{2}+\mathbf{u}_{r}^{\text{T}}\boldsymbol{\eta}_{r,c}\left(\sigma_{s,c}^{2}+\left(\mathbf{v}^{\text{T}}\boldsymbol{\eta}_{s,c}\right)^{2}\hat{e}_{t+i-1}^{2}\right)}{\sigma_{r,c}^{2}+\sigma_{s,c}^{2}+\left(\mathbf{v}^{\text{T}}\boldsymbol{\eta}_{s,c}\right)^{2}\hat{e}_{t+i-1}^{2}}$
(18) $\displaystyle\hat{e}_{t+i}$
$\displaystyle=\sqrt{\frac{\sigma_{r,c}^{2}\left(\sigma_{s,c}^{2}+\left(\mathbf{v}^{\text{T}}\boldsymbol{\eta}_{s,c}\right)^{2}\hat{e}_{t+i-1}^{2}\right)}{\sigma_{r,c}^{2}+\sigma_{s,c}^{2}+\left(\mathbf{v}^{\text{T}}\boldsymbol{\eta}_{s,c}\right)^{2}\hat{e}_{t+i-1}^{2}}}$
(19)

for $\mathbf{v}=\left[0,1\right]^{\text{T}}$, $\hat{s}_{t}=s_{t}$,
$\hat{e}_{t}=0$, and $i=1,2,...,L$.

###### Proof.

See appendix B. ∎

The above Theorem allows to recursively obtain load forecasts $\hat{s}_{t+i}$
together with estimates of their accuracy $\hat{e}_{t+i}$, and a complete
probabilistic description of such forecasts
$N\left(s_{t+i};\hat{s}_{t+i},\hat{e}_{t+i}\right)$, for $i=1,2,...,L$ that
allows to quantify the probability of forecast intervals (see Figure 3). In
order to obtain such $L$ forecasts, Theorem 2 uses instance vectors
$\mathbf{x}$ composed by one past load and observations together with the most
recent parameters. Specifically, for each $i=1,2,...,L$, recursions (18) and
(19) allow to obtain the probabilistic forecast at time $t+i$,
$N\left(s_{t+i};\hat{s}_{t+i},\hat{e}_{t+i}\right)$, using 1) the
probabilistic forecast
$N\left(s_{t+i-1};\hat{s}_{t+i-1},\hat{e}_{t+i-1}\right)$ at previous time
$t+i-1$; 2) observations vector $\mathbf{r}_{t+i}$ at time $t+i$; and 3)
parameters $\boldsymbol{\eta}_{s,c},\sigma_{s,c}$ and
$\boldsymbol{\eta}_{r,c},\sigma_{r,c}$ corresponding with calendar type at
time $t+i$, $c=c(t+i)$.

The results in this section provide theoretical guarantees for the learning
and prediction steps of APLF method. The next section describes APLF in
comparison with approaches based on offline learning, and details the
implementation steps for learning and prediction using APLF.

##  IV Implementation

Offline learning methods obtain one model, while online learning methods
obtain a sequence of models. In particular, APLF method learns a model every
time a new sample is obtained. Figure 4 describes the block diagram for load
forecasting based on offline learning and based on APLF method.

Load forecasting methods based on offline learning train a model using a set
of samples. Later, such methods predict the load for each new instance vector
using the learned model. At learning, $t_{0}$ training samples are used to
obtain the model following different approaches. For instance, ARMA algorithm
calculates parameters for the autoregressive, moving average, and error terms
[22], and techniques based on artificial intelligence calculate parameters
determining a regression function [9, 29, 44]. At prediction, the learned
model and the instance vector $\mathbf{x}_{t}$ are used to obtain load
forecasts $\hat{\mathbf{y}}_{t}$ at time $t$, for $t>{t_{0}}$. The model that
uses the prediction step is the same for all $t$, and the actual loads
$\mathbf{y}_{t}$ are only used to quantify the prediction error (see Figure
4(a)).

Load forecasting methods based on online learning train models regularly using
the most recent samples. Later, such methods predict for each new instance
vector using the latest learned model. At learning, training samples and
possibly the previous model are used to obtain a new model following different
approaches. For instance, the method proposed in [36] recalculates parameters
of the ARMA algorithm, the method proposed in [37] recalculates the weights in
combined forecasts, and the adaptive forgetting factor (AFF) method [38]
updates the functions in the additive model. At prediction, the latest learned
model and the instance vector $\mathbf{x}_{t}$ are used to obtain load
forecasts $\hat{\mathbf{y}}_{t}$ at time $t$.

APLF is a forecasting method based on online learning that updates model
parameters using recursions in Theorem 1 and obtain probabilistic forecasts as
given by Theorem 2. At learning, APLF obtains the new model using the instance
vector $\mathbf{x}_{t}$, the actual loads $\mathbf{y}_{t}$, and the previous
model. At prediction, APLF uses the latest model to obtain load forecasts
$\hat{\mathbf{y}}_{t}=\left[\hat{s}_{t+1},\hat{s}_{t+2},...,\hat{s}_{t+L}\right]^{\text{T}}$
and estimated errors
$\hat{\mathbf{e}}_{t}=\left[\hat{e}_{t+1},\hat{e}_{t+2},...,\hat{e}_{t+L}\right]^{\text{T}}$
that determine probabilistic forecasts as given in (17). The model used for
forecasting adapts to the most recent data, and the actual loads
$\mathbf{y}_{t}$ are not only used to quantify the error, but also to update
the model (see Figure 4(b)).

Algorithm 1 and Algorithm 2 detail the efficient implementation of the
learning and prediction steps of APLF. The corresponding source code in Python
and Matlab languages is publicly available on the web
https://github.com/MachineLearningBCAM/Load-forecasting-IEEE-TPWRS-2020. The
running times of Algorithm 1 and Algorithm 2 are amenable for real-time
implementation with very low latency since APLF has memory complexity
$O(CR^{2})$, the learning step has computational complexity $O(LR^{3})$, and
the prediction step has computational complexity $O(LR)$. Note that the values
of $R$ are small in practice, for instance we use $R=3$, $L=24$, and $C=48$ in
the numerical results of Section V. Algorithm 1 follows recursions given in
Theorem 1 for parameters $\Theta$ and state variables $\Gamma$ described in
(4) and (16), respectively. Specifically, such algorithm updates parameters
$\boldsymbol{\eta}_{s,c},\sigma_{s,c}$ and
$\boldsymbol{\eta}_{r,c},\sigma_{r,c}$ as well as state variables
$\mathbf{P}_{s,c},\gamma_{s,c}$ and $\mathbf{P}_{r,c},\gamma_{r,c}$ using
instance and actual loads vectors with calendar type $c$.

Algorithm 1 Learning step for APLF

Input: $\Theta$ model parameters

$\Gamma$ state variables

$\lambda_{s}$, $\lambda_{r}$ forgetting factors

$\mathbf{x}_{t}=\left[s_{t},\mathbf{r}_{t+1}^{\text{T}},\mathbf{r}_{t+2}^{\text{T}},...,\mathbf{r}_{t+L}^{\text{T}}\right]^{\text{T}}$
new instance vector

$\mathbf{y}_{t}=\left[s_{t+1},s_{t+2},...,s_{t+L}\right]^{\text{T}}$ new loads
vector

$t$ time

Output: $\Theta$ updated model parameters

$\Gamma$ updated state variables

for $i=1,2,...,L$ do

$c\leftarrow c\left(t+i\right)$

$\mathbf{u}_{s}\leftarrow\left[1,s_{t+i-1}\right]^{\text{T}}$

$\mathbf{u}_{r}\leftarrow u_{r}\left(\mathbf{r}_{t+i}\right)$

for $j=s,r$ do $\displaystyle\mathbf{P}_{j,c}$
$\displaystyle\leftarrow\frac{1}{\lambda_{j}}\left(\mathbf{P}_{j,c}-\frac{\mathbf{P}_{j,c}\mathbf{u}_{j}\mathbf{u}_{j}^{\text{T}}\mathbf{P}_{j,c}}{\lambda_{j}+\mathbf{u}_{j}^{\text{T}}\mathbf{P}_{j,c}\mathbf{u}_{j}}\right)$
$\displaystyle\gamma_{j,c}$ $\displaystyle\leftarrow
1+\lambda_{j}\gamma_{j,c}$ $\displaystyle\sigma_{j,c}$
$\displaystyle\leftarrow\sqrt{\sigma_{j,c}^{2}-\frac{1}{\gamma_{j,c}}\left(\sigma_{j,c}^{2}-\frac{\lambda_{j}\left(s_{t+i}-\mathbf{u}_{j}^{\text{T}}\boldsymbol{\eta}_{j,c}\right)^{2}}{\lambda_{j}+\mathbf{u}_{j}^{\text{T}}\mathbf{P}_{j,c}\mathbf{u}_{j}}\right)}$
$\displaystyle\boldsymbol{\eta}_{j,c}$
$\displaystyle\leftarrow\boldsymbol{\eta}_{j,c}+\frac{\mathbf{P}_{j,c}\mathbf{u}_{j}}{\lambda_{j}+\mathbf{u}_{j}^{\text{T}}\mathbf{P}_{j,c}\mathbf{u}_{j}}\left(s_{t+i}-\mathbf{u}_{j}^{\text{T}}\boldsymbol{\eta}_{j,c}\right)$

Algorithm 2 follows initializations and recursions given in Theorem 2 using
the parameters $\Theta$ described in (4) and the new instance vector.
Specifically, such algorithm obtains $L$ load forecasts and $L$ estimates of
their accuracy using the latest parameters
$\boldsymbol{\eta}_{s,c},\sigma_{s,c}$ and
$\boldsymbol{\eta}_{r,c},\sigma_{r,c}$ and the corresponding instance vector.
Note that the proposed method can predict at any time of the day and use
general prediction horizons $L$. In addition, these prediction times and
horizons can change from one day to another just by modifying the
corresponding inputs in Algorithm 2.

Algorithm 2 Prediction step for APLF

Input: $\Theta$ model parameters

$\mathbf{x}_{t}=\left[s_{t},\mathbf{r}_{t+1}^{\text{T}},\mathbf{r}_{t+2}^{\text{T}},...,\mathbf{r}_{t+L}^{\text{T}}\right]^{\text{T}}$
new instance vector

$t$ time

Output:
$\hat{\mathbf{y}}_{t}=\left[\hat{s}_{t+1},\hat{s}_{t+2},...,\hat{s}_{t+L}\right]^{\text{T}}$
load forecasts

$\hat{\mathbf{e}}_{t}=\left[\hat{e}_{t+1},\hat{e}_{t+2},...,\hat{e}_{t+L}\right]^{\text{T}}$
estimated errors

$N(s_{t+i};\hat{s}_{t+i},\hat{e}_{t+i})$, $i=1,2,...,L$

probabilistic forecasts

$\hat{s}_{t}\leftarrow s_{t}$

$\hat{e}_{t}=0$

for $i=1,2,...,L$ do

$c\leftarrow c\left(t+i\right)$

$\hat{\mathbf{u}}_{s}\leftarrow\left[1,\hat{s}_{t+i-1}\right]^{\text{T}}$

$\mathbf{u}_{r}\leftarrow u_{r}\left(\mathbf{r}_{t+i}\right)$

Obtain load forecast $\hat{s}_{t+i}$ using equation (18)

Obtain prediction error $\hat{e}_{t+i}$ using equation (19)

##  V Numerical results

Fig. 5: Consumption patterns can significantly change in two consecutive years. Table I: RMSE and MAPE of prediction errors for APLF and $11$ state-of-the-art techniques on $7$ datasets. Method | Large-size region | Medium-size region | Small-size region
---|---|---|---
Belgium | New Engld. | GEFCom12 | GEFCom2014 | Dayton | 400 builds. | 100 builds.
[GW] |  [$\%$] | [GW] |  [$\%$] | [MW] |  [$\%$] | [MW] |  [$\%$] | [GW] |  [$\%$] | [kW] |  [$\%$] | [kW] |  [$\%$]
LR | $1.47$ | $11.8$ | $1.73$ | $8.0$ | $5.48$ | $20.6$ | $0.30$ | $15.0$ | $0.46$ | $15.5$ | $0.06$ | $9.1$ | $0.07$ | $14.3$
SARIMA | $0.81$ | $5.5$ | $1.22$ | $5.4$ | $3.42$ | $11.6$ | $0.25$ | $12.0$ | $0.22$ | $7.7$ | $0.05$ | $9.9$ | $0.07$ | $16.5$
QR | $1.05$ | $9.2$ | $1.17$ | $5.6$ | $5.46$ | $25.5$ | $0.29$ | $14.7$ | $0.20$ | $7.1$ | $0.08$ | $17.6$ | $0.12$ | $26.8$
GP | 0.52 | 4.0 | $0.89$ | $4.2$ | $2.50$ | 8.5 | 0.24 | 10.6 | $0.19$ | $6.3$ | 0.04 | 6.8 | 0.05 | 11.2
SVM | $0.69$ | $4.7$ | $1.11$ | $5.5$ | $3.28$ | $12.2$ | 0.24 | $12.3$ | 0.17 | 5.7 | 0.04 | $8.0$ | 0.05 | $11.4$
DRN | $1.74$ | $13.0$ | 0.52 | 2.3 | 2.17 | $7.6$ | $0.27$ | $19.5$ | $0.31$ | $11.3$ | 0.04 | $7.0$ | $0.07$ | $14.7$
AR | $0.66$ | $5.1$ | $1.28$ | $5.6$ | $3.94$ | $16.2$ | $0.30$ | $18.5$ | $0.38$ | $13.6$ | 0.04 | $8.9$ | $0.07$ | $16.9$
ARNFS | $1.08$ | $9.4$ | $1.95$ | $10.9$ | $4.41$ | $17.4$ | $0.33$ | $18.4$ | $0.31$ | $14.1$ | 0.04 | $8.1$ | 0.05 | $11.4$
ARRFFS | $1.18$ | $10.3$ | $2.05$ | $11.2$ | $4.54$ | $17.2$ | $0.34$ | $17.7$ | $0.28$ | $10.3$ | $0.08$ | $17.5$ | $0.10$ | $23.4$
SFDA | $1.14$ | $8.9$ | $1.41$ | $10.2$ | $5.04$ | $16.8$ | $0.35$ | $14.5$ | $0.39$ | $21.8$ | 0.06 | 13.0 | 0.08 | 18.6
AFF | $0.95$ | $6.7$ | $1.23$ | $5.8$ | $2.91$ | $10.7$ | $0.25$ | $15.2$ | $0.26$ | $9.6$ | 0.05 | 10.2 | 0.07 | 14.6
APLF | 0.33 | 2.3 | 0.86 | 3.9 | 2.15 | 8.1 | 0.20 | 9.6 | 0.16 | 5.5 | 0.03 | 6.3 | 0.05 | 11.0

This section first describes the datasets used for the experimentation, and
then compares the performance of APLF method with respect to that of existing
techniques. The first set of numerical results quantifies the prediction
errors, while the second set of numerical results evaluates the performance of
probabilistic load forecasts and analyzes the relationship between training
size and prediction error.

(a) New England.

(b) Dayton.

(c) $100$ buildings.

Fig. 6: Load demand presents different consumption patterns and variability in
the three regions with different sizes.

(a) New England.

(b) Dayton.

(c) $100$ buildings.

Fig. 7: CDFs of prediction errors provide detailed performance comparison
between APLF method and conventional techniques in three regions with
different sizes.

Seven publicly available datasets are selected for numerical experimentation.
The datasets correspond with regions that have different sizes and display
different consumption patterns that change over time. Such changes can be
observed in Figure 5 that shows load demand per hour of day and day of year
during 2004 and 2005 in Dayton (US). This figure shows that consumption
patterns change significantly not only for different seasons but also between
consecutive weeks and between consecutive years. Therefore, methods based on
static models often obtain inferior accuracies since they cannot adapt to
dynamic changes in consumption patterns.

We group the seven datasets by size of the region: large, medium, and small.
Two datasets belong to large-size regions: load demand in Belgium from
2017-2019 made available by Elia group, and load demand in New England from
2003-2014 made available by ISO-NE organization. Three datasets belong to
medium-size regions: Global Energy Forecasting Competition 2012 (GEFCom2012)
dataset from 2004-2007 [45], Global Energy Forecasting Competition 2014
(GEFCom2014) dataset from 2005-2011 [39], and load demand in Dayton from
2004-2016 made available by PJM interconnection. Finally, two datasets belong
to small-size regions that correspond with load demand for $400$ and $100$
buildings in New South Wales from 2013 and are made available by the
Australian Government as part of the project Smart Grid Smart Cities.

In the numerical results, training for offline learning algorithms is done
using different number of samples depending on the datasets. Training sizes
for offline learning algorithms are chosen taking into account the length of
datasets: two years of data are used for training in New England and
GEFCom2012 datasets; one year of data are used for training in Belgium,
GEFCom2014, and Dayton datasets; and $3/4$ of a year of data are used for
training in $400$ and $100$ buildings datasets. Prediction for all algorithms
is done using the rest of the data. Specifically, at $11$ a.m. of each day,
the forecasting methods obtain future loads for a prediction horizon of $L=24$
hours hence, every vector of load forecasts is formed by forecasts obtained
from $1$ to $L$ hours ahead.

APLF results are obtained using the following implementation details. The
instance vector composed by past loads and observations is given by
$\mathbf{x}=\left[s_{t},\mathbf{r}_{t+1}^{\text{T}},...,\mathbf{r}_{t+L}^{\text{T}}\right]^{\text{T}}$.
The observations vector $\mathbf{r}_{t}$ contains the temperature $w_{t}$ at
time $t$ and the mean of past temperatures $\bar{w}_{c\left(t\right)}$ at
calendar type $c(t)$, i.e.,
$\mathbf{r}_{t}=\left[w_{t},\bar{w}_{c\left(t\right)}\right]^{\text{T}}$. The
observations vector $\mathbf{r}_{t}$ is represented by the feature vector
$u_{r}\left(\mathbf{r}_{t}\right)$ for a function $u_{r}(\cdot)$ that encodes
temperature shifts. Specifically, such function determines the vector
$u_{r}\left(\mathbf{r}_{t}\right)=\left[1,\alpha_{1},\alpha_{2}\right]^{\text{T}}$,
where $\alpha_{1}$ (resp. $\alpha_{2}$) takes value $1$ if the temperature is
above (resp. below) certain thresholds and takes value $0$ otherwise, that is

$\displaystyle\alpha_{1},\alpha_{2}=\left\\{\begin{matrix}1,0&\text{if}&w_{t}-\bar{w}_{c\left(t\right)}&>&W_{1}&\text{and}\\\
&&w_{t}>W_{2}&\text{or}&w_{t}<W_{3}&\\\
0,1&\text{if}&w_{t}-\bar{w}_{c\left(t\right)}&<&-W_{1}&\text{and}\\\
&&w_{t}>W_{2}&\text{or}&w_{t}<W_{3}&\\\ 0,0&&\text{otherwise}&&&&&\\\
\end{matrix}\right.$

where we take threshold values $W_{1}=20^{\circ}F$, $W_{2}=80^{\circ}F$, and
$W_{3}=20^{\circ}F$ for all datasets. The calendar information $c(t)$
specifies the type of hour: $c(t)$ from $1$ to $24$ indicates the hour of day
of weekdays, and $c(t)$ from $25$ to $48$ indicates the hour of day of
weekends and holidays, i.e., $c\left(t\right)\in\left\\{1,...,C\right\\}$ with
$C=48$. Then, APLF method obtains parameters
$\boldsymbol{\eta}_{s,c},\sigma_{s,c}$ and
$\boldsymbol{\eta}_{r,c},\sigma_{r,c}$ for $C=48$ calendar types as given by
Algorithm 1 initialized with $\boldsymbol{\eta}_{0}=\mathbf{0}_{K}$, any
$\sigma_{0}$, $\mathbf{P}_{0}=\mathbf{I}_{K}$, and $\gamma_{0}=0$. Such
parameters are updated by taking forgetting factors as
$\lambda_{s}\leavevmode\nobreak\ =\leavevmode\nobreak\ 0.2$ and
$\lambda_{r}=0.7$ for any calendar type and for all datasets.111Possible
numerical instabilities of state matrices $\mathbf{P}$ in Algorithm 1 are
addressed by their reinitialization in case their trace becomes larger than
$10$ similarly to methods based on recursive least squares [46]. Forgetting
factors ($\lambda_{s},\lambda_{r}$) and threshold values ($W_{1},W_{2},W_{3}$)
are the hyper-parameters of APLF method. Hyper-parameters’ values can be
selected by methods such as cross-validation over a grid of possible values.
For simplicity, we select values for hyper-parameters by inspection over one
dataset and then we use the same values in all datasets. The numerical results
corroborate the robustness of APLF to the choice of hyper-parameters since we
use the same hyper-parameters’ values in all numerical results and with
significantly different datasets.

APLF method is compared with $11$ state-of-the-art techniques based on
statistical methods, artificial intelligence, and weighted combination of
several forecasts. Three techniques are based on statistical methods: LR [9],
SARIMA [25], and QR [34]; six techniques are based on artificial intelligence:
GP [33], SVM [29], deep residual network (DRN) [27], and three versions of AR-
NARX method [9] based on linear regression (AR) [9], fixed size least squares
SVM using the Nyström method (ARNFS) [44], and fixed size least squares SVM
using Random Features (ARRFFS) [47]; finally, two techniques are based on
weighted combination of several forecasts: secondary forecasting based on
deviation analysis (SFDA) [31] and AFF [38].

Table II: Pinball loss and ECE for APLF and $2$ state-of-the-art techniques on $7$ datasets. Region | QR | GP | APLF
---|---|---|---
Pinball loss | ECE | Pinball loss | ECE | Pinball loss | ECE
Belgium | 0.34 [GW] | 0.08 | 0.14 [GW] | 0.19 |  0.11 [GW] | 0.07
New Engld. | 0.70 [GW] | 0.07 | 0.24 [GW] | 0.09 |  0.22 [GW] | 0.07
GEFCom12 | 1.03 [MW] | 0.06 | 0.78 [MW] | 0.14 |  0.77 [MW] | 0.12
GEFCom14 | 0.06 [MW] | 0.60 |  0.05 [MW] | 0.15 | 0.06 [MW] | 0.19
Dayton | 0.09 [GW] | 0.06 |  0.04 [GW] | 0.12 |  0.04 [GW] | 0.05
400 builds. | 0.02 [kW] | 0.10 |  0.01 [kW] | 0.05 |  0.01 [kW] | 0.07
100 builds. | 0.03 [kW] | 0.08 |  0.01 [kW] | 0.03 |  0.01 [kW] | 0.08
Fig. 8: CDFs of pinball losses compare the probabilistic performance of APLF
method with state-of-the-art probabilistic techniques.

(a) New England.

(b) Dayton.

(c) $100$ buildings.

Fig. 9: Calibration plots describe the reliability of probabilistic forecasts
of APLF method in comparison with state-of-the-art probabilistic techniques in
three regions with different sizes.

In the first set of numerical results we quantify the prediction error of APLF
in comparison with $11$ conventional techniques for the $7$ datasets. Overall
prediction errors RMSE and MAPE assessing prediction errors are given in Table
I. Such table shows that existing techniques such as DRN and QR can achieve
high accuracies in certain large-size regions using sizeable training datasets
(e.g., New England dataset), however such techniques become inaccurate in
other datasets such as those corresponding with small-size regions and smaller
training datasets (e.g., $100$ buildings dataset). In addition, Table I also
shows that the online learning method AFF achieves higher accuracies than
multiple offline learning algorithms such as LR, ARRFFS, and SFDA. Figures 6
and 7 provide more detailed comparisons using $5$ representative existing
techniques (AR, SARIMA, QR, SVM, and AFF), and $3$ datasets that correspond
with regions of assorted sizes. Figure 6 shows two days of load demand and
load forecasts in the three regions and Figure 7 shows the empirical
cumulative distribution functions (CDFs) of the absolute value of prediction
errors. Table I, and Figures 6 and 7 show that the proposed APLF method
achieves high accuracies in comparison with conventional techniques in every
dataset studied. In particular, Figure 7 shows that high errors occur with low
probability for APLF method. For instance, in New England dataset, the error
of APLF method is less than $0.8$ GW with probability $0.8$, while the $4$
other methods reach errors of around $1.3$ GW with such probability.

Fig. 10: The length of the training dataset does not affect the accuracy of
online learning algorithms but significantly affects the accuracy of offline
learning algorithms.

In the second set of numerical results we quantify the probabilistic
performance of APLF in comparison with QR [34] and GP [33] and we study the
relationship between training size and prediction error. Pinball loss and ECE
assessing probabilistic forecasts are given in Table II. Such table shows that
APLF achieves high performance in terms of both pinball loss and ECE, while GP
sometimes achieves poor results in terms of ECE and QR achieves poor results
in terms of pinball loss. Figures 8 and 9 provide more detailled
quantification of the probabilistic performance of different methods. Figure 8
shows the empirical CDFs of pinball losses of APLF method, QR, GP , and the
benchmark for the GEFCom2014 dataset [39]. These CDFs show that the
probability of high pinball losses is significantly lower for APLF method. In
particular, the CDFs in Figure 8 show that APLF and GP have a similar median
pinball loss of around $0.04$ MW. However, APLF has pinball losses less than
$0.08$ MW with probability $0.8$, while GP reaches pinball losses of $0.16$ MW
with probability $0.8$. Figure 9 shows the correspondence between the
calibration $C(q)$ of probabilistic forecasts and the quantile $q$ for the
datasets used in Figures 6 and 7. These calibration plots show that GP and QR
tend to obtain forecast quantiles higher than the true quantiles, while APLF
obtains more unbiased probabilistic forecasts. In particular, the true load is
higher than the 50 quantile forecast load with probability very near 50 $\%$
for APLF. Figure 9 shows that APLF obtains improved calibrations especially in
the lower quantiles.

Finally, Figure 10 shows the RMSE obtained by APLF method and the $5$
conventional techniques shown in Figures 6 and 7 for different sizes of
training sets using GEFCom2012 dataset. These results are obtained computing
RMSEs $20$ times for each size of training set. The samples used for training
in these numerical results are different at each experiment and testing sets
always contain two years of data. As can be observed from Figure 10, the
accuracy of online learning algorithms does not significantly change with the
length of the training dataset, while offline learning algorithms require
large training datasets to achieve accurate results.

APLF method achieves remarkable results both in terms of single-value and
probabilistic forecasts, and adapts to different consumption patterns in every
region studied even where variability in load demand is more significant.
Numerical results confirm that APLF better captures dynamic changes in
consumption patterns than existing methods.

##  VI Conclusion

The paper proposes techniques for adaptive probabilistic load forecasting
(APLF) that can adapt to changes in consumption patterns and assess
consumption uncertainties. We developed online learning techniques that update
model parameters using a novel recursive algorithm, and we developed
prediction techniques that obtain probabilistic forecasts using the most
recent parameters. In addition, we described the theoretical guarantees and
efficient implementation of the online learning and probabilistic prediction
steps for APLF. The paper also compared the accuracy of the proposed APLF with
convectional techniques in multiple datasets. These datasets represent
challenging scenarios with different sizes and different consumption patterns
that change over time. The experimental results showed the performance
improvement of APLF method in terms of prediction errors and probabilistic
forecasts. The results showed that APLF can improve performance in a wide
range of scenarios, and consider general load-related observations and
prediction horizons.

##  Appendix A Proof Of Theorem 1

###### Proof.

We first prove that for any $i>0$ optimal parameters
$\boldsymbol{\eta}_{i}^{*}$, $\sigma_{i}^{*}$ satisfy

$\displaystyle\underset{j=1}{\overset{i}{\sum}}\lambda^{i-j}\mathbf{u}_{t_{j}}\mathbf{u}_{t_{j}}^{\text{T}}\boldsymbol{\eta}_{i}^{*}$
$\displaystyle=\mathbf{q}_{i}$ (20)
$\displaystyle\gamma_{i}{\sigma_{i}^{*}}^{2}$
$\displaystyle=\underset{j=1}{\overset{i}{\sum}}\lambda^{{i}-j}s_{t_{j}}^{2}-\mathbf{q}_{i}^{\text{T}}\boldsymbol{\eta}_{i}^{*}$
(21)

while parameters $\boldsymbol{\eta}_{i}$, $\sigma_{i}$ and matrix
$\mathbf{P}_{i}$ given by recursions (7)-(10) with
$\boldsymbol{\eta}_{0}=\mathbf{0}_{K}$, any $\sigma_{0}$,
$\mathbf{P}_{0}=\mathbf{I}_{K}$, and $\gamma_{0}=0$ satisfy

$\displaystyle\mathbf{P}_{i}^{-1}\boldsymbol{\eta}_{i}$
$\displaystyle=\mathbf{q}_{i}$ (22) $\displaystyle\gamma_{i}{\sigma_{i}}^{2}$
$\displaystyle=\underset{j=1}{\overset{i}{\sum}}\lambda^{{i}-j}s_{t_{j}}^{2}-\mathbf{q}_{i}^{\text{T}}\boldsymbol{\eta}_{i}$
(23) $\displaystyle\mathbf{P}_{i}^{-1}$
$\displaystyle=\lambda^{i}\mathbf{I}_{K}+\underset{j=1}{\overset{i}{\sum}}\lambda^{i-j}\mathbf{u}_{t_{j}}\mathbf{u}_{t_{j}}^{\text{T}}$
(24)

where
$\mathbf{q}_{i}=\underset{j=1}{\overset{{i}}{\sum}}\lambda^{{i}-j}s_{t_{j}}\mathbf{u}_{t_{j}}$.
Then, in the second step of the proof we obtain bound in (11) using equations
(20)-(24). Finally, we prove that parameters $\boldsymbol{\eta}_{i}$, and
$\sigma_{i}$ given by recursions (7)-(10) with $\boldsymbol{\eta}_{i_{0}}$,
$\sigma_{i_{0}}$ given by (12)-(15) satisfy
$\boldsymbol{\eta}_{i}=\boldsymbol{\eta}_{i}^{*}$ and
$\sigma_{i}=\sigma_{i}^{*}$, for $i\geq i_{0}$.

Parameters $\boldsymbol{\eta}_{i}^{*}$ and $\sigma_{i}^{*}$ satisfy equations
(20) and (21), respectively, because they maximize the log-likelihood in (5).
The differentiable function $L_{i}\left(\boldsymbol{\eta},\sigma\right)$ is
concave since Gaussian functions are log-concave. Then,
$L_{i}(\boldsymbol{\eta},\sigma)$ has a maximum achieved by parameters that
result in zero derivatives. Since

$\displaystyle L_{i}\left(\boldsymbol{\eta},\sigma\right)$
$\displaystyle=-\underset{j=1}{\overset{{i}}{\sum}}\lambda^{{i}-j}\frac{\big{(}s_{t_{j}}-\mathbf{u}_{t_{j}}^{\text{T}}\boldsymbol{\eta}\big{)}^{2}}{2\sigma^{2}}+\lambda^{{i}-j}\log\big{(}\sigma\sqrt{2\pi}\big{)}$

we have that

$\displaystyle\frac{\partial
L_{i}(\boldsymbol{\eta},\sigma)}{\partial\boldsymbol{\eta}}$
$\displaystyle=\underset{j=1}{\overset{{i}}{\sum}}\lambda^{{i}-j}\frac{\mathbf{u}_{t_{j}}(s_{t_{j}}-\mathbf{u}_{t_{j}}^{\text{T}}\boldsymbol{\eta})}{{\sigma}^{2}}$

that becomes zero for $\boldsymbol{\eta}_{i}^{*}$ given by (20), and

$\displaystyle\frac{\partial
L_{i}(\boldsymbol{\eta},\sigma)}{\partial\sigma}=\underset{j=1}{\overset{i}{\sum}}\lambda^{{i}-j}\frac{\big{(}s_{t_{j}}-\mathbf{u}_{t_{j}}^{\text{T}}\boldsymbol{\eta}\big{)}^{2}}{{\sigma}^{3}}-\lambda^{{i}-j}\frac{1}{\sigma}$

that becomes zero for $\sigma_{i}^{*}$ given by (21) since $\gamma_{i}$ given
by (10) equals $\gamma_{i}=\underset{j=1}{\overset{i}{\sum}}\lambda^{{i}-j}$.

By induction, we prove that $\boldsymbol{\eta}_{i}$ and $\sigma_{i}$ given by
recursions (7) and (8) satisfy equations (22) and (23) for
$\boldsymbol{\eta}_{0}=\mathbf{0}_{K}$, any $\sigma_{0}$,
$\mathbf{P}_{0}=\mathbf{I}_{K}$, and $\gamma_{0}=0$. Firstly, we prove it for
$i=1$. From (9), we have that

$\displaystyle\mathbf{P}_{1}$
$\displaystyle=\frac{1}{\lambda}\left(\mathbf{I}_{K}-\frac{\mathbf{u}_{t_{1}}\mathbf{u}_{t_{1}}^{\text{T}}}{\lambda+\mathbf{u}_{t_{1}}^{\text{T}}\mathbf{u}_{t_{1}}}\right)=\left(\lambda\mathbf{I}_{K}+\mathbf{u}_{t_{1}}\mathbf{u}_{t_{1}}^{\text{T}}\right)^{-1}$
(25)

applying the matrix inversion Lemma and using that
$\mathbf{P}_{0}=\mathbf{I}_{K}$. Hence, from (7), (8), and (10), we get

$\displaystyle\mathbf{P}_{1}^{-1}\boldsymbol{\eta}_{1}$
$\displaystyle=\left(\lambda\mathbf{I}_{K}+\mathbf{u}_{t_{1}}\mathbf{u}_{t_{1}}^{\text{T}}\right)\frac{s_{t_{1}}\mathbf{u}_{t_{1}}}{\lambda+\mathbf{u}_{t_{1}}^{\text{T}}\mathbf{u}_{t_{1}}}=s_{t_{1}}\mathbf{u}_{t_{1}}$
$\displaystyle\sigma_{1}^{2}$ $\displaystyle=\frac{\lambda
s_{t_{1}}^{2}}{\lambda+\mathbf{u}_{t_{1}}^{\text{T}}\mathbf{u}_{t_{1}}}=\frac{\lambda
s_{t_{1}}^{2}+s_{t_{1}}^{2}\mathbf{u}_{t_{1}}^{\text{T}}\mathbf{u}_{t_{1}}-s_{t_{1}}^{2}\mathbf{u}_{t_{1}}^{\text{T}}\mathbf{u}_{t_{1}}}{\lambda+\mathbf{u}_{t_{1}}^{\text{T}}\mathbf{u}_{t_{1}}}$
$\displaystyle=s_{t_{1}}^{2}-s_{t_{1}}\mathbf{u}_{t_{1}}^{\text{T}}\frac{s_{t_{1}}\mathbf{u}_{t_{1}}}{\lambda+\mathbf{u}_{t_{1}}^{\text{T}}\mathbf{u}_{t_{1}}}=s_{t_{1}}^{2}-s_{t_{1}}\mathbf{u}_{t_{1}}^{\text{T}}\boldsymbol{\eta}_{1}$

since $\boldsymbol{\eta}_{0}=\mathbf{0}_{K}$, $\mathbf{P}_{0}=\mathbf{I}_{K}$,
and $\gamma_{0}=0$.

If (22) and (23) hold for $i-1$, then for $i$ we have that

$\displaystyle\mathbf{P}_{i}$
$\displaystyle=\left(\lambda{\mathbf{P}_{i-1}}^{-1}+\mathbf{u}_{t_{i}}\mathbf{u}_{t_{i}}^{\text{T}}\right)^{-1}$
(26)

applying the matrix inversion Lemma to equation (9). Therefore, using the
recursion of $\boldsymbol{\eta}_{i}$ in (7), we have that

$\displaystyle\mathbf{P}_{i}^{-1}\boldsymbol{\eta}_{i}=$
$\displaystyle\lambda{\mathbf{P}_{i-1}}^{-1}\boldsymbol{\eta}_{i-1}+\frac{\lambda\mathbf{u}_{t_{i}}}{\lambda+\mathbf{u}_{t_{i}}^{\text{T}}\mathbf{P}_{i-1}\mathbf{u}_{t_{i}}}\left(s_{t_{i}}-\mathbf{u}_{t_{i}}^{\text{T}}\boldsymbol{\eta}_{i-1}\right)$
$\displaystyle+\mathbf{u}_{t_{i}}\mathbf{u}_{t_{i}}^{\text{T}}\boldsymbol{\eta}_{i-1}+\frac{\mathbf{u}_{t_{i}}\mathbf{u}_{t_{i}}^{\text{T}}\mathbf{P}_{i-1}\mathbf{u}_{t_{i}}}{\lambda+\mathbf{u}_{t_{i}}^{\text{T}}\mathbf{P}_{i-1}\mathbf{u}_{t_{i}}}\left(s_{t_{i}}-\mathbf{u}_{t_{i}}^{\text{T}}\boldsymbol{\eta}_{i-1}\right)$
$\displaystyle=$
$\displaystyle\lambda\mathbf{q}_{i-1}+\frac{\lambda\mathbf{u}_{t_{i}}}{\lambda+\mathbf{u}_{t_{i}}^{\text{T}}\mathbf{P}_{i-1}\mathbf{u}_{t_{i}}}\left(s_{t_{i}}-\mathbf{u}_{t_{i}}^{\text{T}}\boldsymbol{\eta}_{i-1}\right)$
$\displaystyle+\mathbf{u}_{t_{i}}\mathbf{u}_{t_{i}}^{\text{T}}\boldsymbol{\eta}_{i-1}+\frac{\mathbf{u}_{t_{i}}\mathbf{u}_{t_{i}}^{\text{T}}\mathbf{P}_{i-1}\mathbf{u}_{t_{i}}}{\lambda+\mathbf{u}_{t_{i}}^{\text{T}}\mathbf{P}_{i-1}\mathbf{u}_{t_{i}}}\left(s_{t_{i}}-\mathbf{u}_{t_{i}}^{\text{T}}\boldsymbol{\eta}_{i-1}\right)$
(27) $\displaystyle=$
$\displaystyle\lambda\mathbf{q}_{i-1}+\mathbf{u}_{t_{i}}\left(s_{t_{i}}-\mathbf{u}_{t_{i}}^{\text{T}}\boldsymbol{\eta}_{i-1}\right)+\mathbf{u}_{t_{i}}\mathbf{u}_{t_{i}}^{\text{T}}\boldsymbol{\eta}_{i-1}$
$\displaystyle=$
$\displaystyle\underset{j=1}{\overset{i-1}{\sum}}\lambda^{i-j}s_{t_{j}}\mathbf{u}_{t_{j}}+s_{t_{i}}\mathbf{u}_{t_{i}}=\mathbf{q}_{i}$

where the equality (27) is obtained by using the induction hypothesis. Using
the recursion of $\sigma_{i}$ in (8), we have that

$\displaystyle\gamma_{i}{\sigma_{i}}^{2}=$
$\displaystyle\left(\gamma_{i}-1\right){\sigma_{i-1}}^{2}+\left(s_{t_{i}}-\mathbf{u}_{t_{i}}^{\text{T}}\boldsymbol{\eta}_{i-1}\right)$
$\displaystyle\cdot\left(s_{t_{i}}-\frac{s_{t_{i}}\mathbf{u}_{t_{i}}^{\text{T}}\mathbf{P}_{i-1}\mathbf{u}_{t_{i}}}{\lambda+\mathbf{u}_{t_{i}}^{\text{T}}\mathbf{P}_{i-1}\mathbf{u}_{t_{i}}}-\frac{\lambda\mathbf{u}_{t_{i}}^{\text{T}}\boldsymbol{\eta}_{i-1}}{\lambda+\mathbf{u}_{t_{i}}^{\text{T}}\mathbf{P}_{i-1}\mathbf{u}_{t_{i}}}\right)$
$\displaystyle=$
$\displaystyle\underset{j=1}{\overset{i-1}{\sum}}\lambda^{{i}-j}s_{t_{j}}^{2}-\lambda\mathbf{q}_{i-1}^{\text{T}}\boldsymbol{\eta}_{i-1}+\left(s_{t_{i}}-\mathbf{u}_{t_{i}}^{\text{T}}\boldsymbol{\eta}_{i-1}\right)$
(28)
$\displaystyle\cdot\left(s_{t_{i}}-\frac{s_{t_{i}}\mathbf{u}_{t_{i}}^{\text{T}}\mathbf{P}_{i-1}\mathbf{u}_{t_{i}}}{\lambda+\mathbf{u}_{t_{i}}^{\text{T}}\mathbf{P}_{i-1}\mathbf{u}_{t_{i}}}-\frac{\lambda\mathbf{q}_{i-1}^{\text{T}}\mathbf{P}_{i-1}\mathbf{u}_{t_{i}}}{\lambda+\mathbf{u}_{t_{i}}^{\text{T}}\mathbf{P}_{i-1}\mathbf{u}_{t_{i}}}\right)$
(29) $\displaystyle=$
$\displaystyle\underset{j=1}{\overset{i-1}{\sum}}\lambda^{{i}-j}s_{t_{j}}^{2}+s_{t_{i}}^{2}-\left(\underset{j=1}{\overset{i-1}{\sum}}\lambda^{i-j}s_{t_{j}}\mathbf{u}_{t_{j}}^{\text{T}}+s_{t_{i}}\mathbf{u}_{t_{i}}^{\text{T}}\right)$
$\displaystyle\cdot\left(\boldsymbol{\eta}_{i-1}+\frac{\mathbf{P}_{i-1}\mathbf{u}_{t_{i}}}{\lambda+\mathbf{u}_{t_{i}}^{\text{T}}\mathbf{P}_{i-1}\mathbf{u}_{t_{i}}}\left(s_{t_{i}}-\mathbf{u}_{t_{i}}^{\text{T}}\boldsymbol{\eta}_{i-1}\right)\right)$
(30)

where the equality (29) is obtained by using the induction hypothesis. Then,
we obtain (23) from (30) by using the recursion for $\boldsymbol{\eta}_{i}$ in
(7).

In addition, by induction we prove that $\mathbf{P}_{i}$ given by recursion
(9) satisfies the equation (24). The case $i=1$ is proved in the equation
(25). If (24) holds for $i-1$, then for $i$ by using (26) and the induction
hypothesis we have that

$\displaystyle\mathbf{P}_{i}$
$\displaystyle=\bigg{(}\lambda\Big{(}\lambda^{i-1}\mathbf{I}_{K}+\underset{j=1}{\overset{i-1}{\sum}}\lambda^{i-1-j}\mathbf{u}_{t_{j}}\mathbf{u}_{t_{j}}^{\text{T}}\Big{)}+\mathbf{u}_{t_{i}}\mathbf{u}_{t_{i}}^{\text{T}}\bigg{)}^{-1}$

that directly leads to (24).

To obtain the bound in (11), we first use the definition of $L_{i}$ and
equations (21) and (23) to obtain

$\displaystyle L_{i}$
$\displaystyle\left(\boldsymbol{\eta}_{i}^{*},\sigma_{i}^{*}\right)-L_{i}\left(\boldsymbol{\eta}_{i},\sigma_{i}\right)=\frac{\gamma_{i}}{2}\log\left(\frac{\sigma_{i}^{2}}{{\sigma_{i}^{*}}^{2}}\right)$
$\displaystyle=\frac{\gamma_{i}}{2}\log\left(1+\frac{\sigma_{i}^{2}-{\sigma_{i}^{*}}^{2}}{{\sigma_{i}^{*}}^{2}}\right)\leq\frac{\gamma_{i}}{2}\frac{\Big{|}\sigma_{i}^{2}-{\sigma_{i}^{*}}^{2}\Big{|}}{{\sigma_{i}^{*}}^{2}}$
$\displaystyle=\frac{1}{2}\left|\frac{\mathbf{q}_{i}^{\text{T}}\left(\boldsymbol{\eta}_{i}^{*}-\boldsymbol{\eta}_{i}\right)}{{\sigma_{i}^{*}}^{2}}\right|.$
(31)

We then use following inequalities

$\displaystyle\left|\mathbf{q}_{i}^{\text{T}}\left(\boldsymbol{\eta}_{i}^{*}-\boldsymbol{\eta}_{i}\right)\right|$
$\displaystyle\leq\lambda^{i}\left\|{\boldsymbol{\eta}_{i}^{*}}\right\|\left\|\boldsymbol{\eta}_{i}\right\|$
(32) $\displaystyle\frac{1}{{\sigma_{i}^{*}}^{2}}$ $\displaystyle\leq 2\pi
M^{2}$ (33)

where the inequality (32) is obtained using equations (20) and (22) because

$\displaystyle\big{|}$
$\displaystyle\mathbf{q}_{i}^{\text{T}}\left(\boldsymbol{\eta}_{i}^{*}-\boldsymbol{\eta}_{i}\right)\big{|}=\Bigg{|}\mathbf{q}_{i}^{\text{T}}\Bigg{(}\bigg{(}\underset{j=1}{\overset{i}{\sum}}\lambda^{i-j}\mathbf{u}_{t_{j}}\mathbf{u}_{t_{j}}^{\text{T}}\bigg{)}^{-1}-\mathbf{P}_{i}\Bigg{)}\mathbf{q}_{i}\Bigg{|}$
$\displaystyle=\Bigg{|}\mathbf{q}_{i}^{\text{T}}\bigg{(}\underset{j=1}{\overset{i}{\sum}}\lambda^{i-j}\mathbf{u}_{t_{j}}\mathbf{u}_{t_{j}}^{\text{T}}\bigg{)}^{-1}\bigg{(}\mathbf{P}_{i}^{-1}-\underset{j=1}{\overset{i}{\sum}}\lambda^{i-j}\mathbf{u}_{t_{j}}\mathbf{u}_{t_{j}}^{\text{T}}\bigg{)}\mathbf{P}_{i}\mathbf{q}_{i}\Bigg{|}$
$\displaystyle=\Bigg{|}\lambda^{i}\mathbf{q}_{i}^{\text{T}}\bigg{(}\underset{j=1}{\overset{i}{\sum}}\lambda^{i-j}\mathbf{u}_{t_{j}}\mathbf{u}_{t_{j}}^{\text{T}}\bigg{)}^{-1}\mathbf{P}_{i}\mathbf{q}_{i}\Bigg{|}=\lambda^{i}\left|{\boldsymbol{\eta}_{i}^{*}}^{\text{T}}\boldsymbol{\eta}_{i}\right|$

and the inequality (33) is obtained due to the fact that

$\displaystyle\log\left(\frac{1}{\sigma_{i}^{*}}\right)-\log\sqrt{2\pi}\leq\log
M\Rightarrow\left|\frac{1}{\sigma_{i}^{*}}\right|\leq\sqrt{2\pi}M$

because

$L_{i}\left(\boldsymbol{\eta}_{i}^{*},\sigma_{i}^{*}\right)=\gamma_{i}(-\log\sigma_{i}^{*}-\log\sqrt{2\pi})$

and

$L_{i}\left(\boldsymbol{\eta}_{i}^{*},\sigma_{i}^{*}\right)\leq\gamma_{i}\log
M$

since
$N(s_{t_{j}};\mathbf{u}_{t_{j}}^{\text{T}}\boldsymbol{\eta}_{i}^{*},\sigma_{i}^{*})\leq
M$ for any $j\leq i\leq n$.

Substituting inequalities (32) and (33) in (31), we have that

$\displaystyle
L_{i}\left(\boldsymbol{\eta}_{i}^{*},\sigma_{i}^{*}\right)-L_{i}\left(\boldsymbol{\eta}_{i},\sigma_{i}\right)\leq$
$\displaystyle\pi
M^{2}\lambda^{i}\left\|{\boldsymbol{\eta}_{i}^{*}}\right\|\left\|\boldsymbol{\eta}_{i}\right\|$

that leads to bound in (11) using the definition of $\boldsymbol{\eta}_{i}$
given by (22) and the following inequalities

$\displaystyle\Bigg{\|}\mathbf{P}_{i}\underset{j=1}{\overset{i}{\sum}}\lambda^{i-j}\mathbf{u}_{t_{j}}\mathbf{u}_{t_{j}}^{\text{T}}\boldsymbol{\eta}_{i}^{*}\Bigg{\|}$
$\displaystyle\leq\Bigg{\|}\mathbf{P}_{i}\underset{j=1}{\overset{i}{\sum}}\lambda^{i-j}\mathbf{u}_{t_{j}}\mathbf{u}_{t_{j}}^{\text{T}}\Bigg{\|}\left\|\boldsymbol{\eta}_{i}^{*}\right\|\leq\left\|\boldsymbol{\eta}_{i}^{*}\right\|$

where the last inequality is obtained because for any $i$ such that the matrix
(6) is not singular, we have that

$\displaystyle\bigg{(}\underset{j=1}{\overset{i}{\sum}}\lambda^{i-j}\mathbf{u}_{t_{j}}\mathbf{u}_{t_{j}}^{\text{T}}\bigg{)}^{-1}\bigg{(}\underset{j=1}{\overset{i}{\sum}}\lambda^{i-j}\mathbf{u}_{t_{j}}\mathbf{u}_{t_{j}}^{\text{T}}+\lambda^{i}\mathbf{I}_{K}\bigg{)}\succeq\mathbf{I}_{K}$

which implies

$\displaystyle\bigg{(}\underset{j=1}{\overset{i}{\sum}}\lambda^{i-j}\mathbf{u}_{t_{j}}\mathbf{u}_{t_{j}}^{\text{T}}+\lambda^{i}\mathbf{I}_{K}\bigg{)}^{-1}\underset{j=1}{\overset{i}{\sum}}\lambda^{i-j}\mathbf{u}_{t_{j}}\mathbf{u}_{t_{j}}^{\text{T}}$
$\displaystyle=\mathbf{P}_{i}\underset{j=1}{\overset{i}{\sum}}\lambda^{i-j}\mathbf{u}_{t_{j}}\mathbf{u}_{t_{j}}^{\text{T}}\preceq\mathbf{I}_{K}.$

Now, we proof by induction that for any $i\geq i_{0}$ parameters
$\boldsymbol{\eta}_{i}$ and $\sigma_{i}$ given by recursions (7)-(10) with
$\boldsymbol{\eta}_{i_{0}}$, $\sigma_{i_{0}}$, $\mathbf{P}_{i_{0}}$, and
$\gamma_{i_{0}}$ given by (12)-(15) satisfy
$\boldsymbol{\eta}_{i}=\boldsymbol{\eta}_{i}^{*}$,
$\sigma_{i}=\sigma_{i}^{*}$, and
$\mathbf{P}_{i}=(\underset{j=1}{\overset{i}{\sum}}\lambda^{i-j}\mathbf{u}_{t_{j}}\mathbf{u}_{t_{j}}^{\text{T}})^{-1}$.
Firstly, for $i=i_{0}$ the assertions are obtained directly from (12)-(15)
since $\mathbf{H}_{i_{0}}$ is non-singular and $\boldsymbol{\eta}_{i}^{*}$ and
$\sigma_{i}^{*}$ satisfy 20 and 21, respectively.

If $\boldsymbol{\eta}_{i-1}=\boldsymbol{\eta}_{i-1}^{*}$ and
$\sigma_{i-1}=\sigma_{i-1}^{*}$ hold, then for $i$ we have that

$\displaystyle\mathbf{P}_{i}=(\lambda\mathbf{P}_{i-1}^{-1}+\mathbf{u}_{t_{i}}\mathbf{u}_{t_{i}}^{\text{T}})^{-1}=\bigg{(}\underset{j=1}{\overset{i}{\sum}}\lambda^{i-j}\mathbf{u}_{t_{j}}\mathbf{u}_{t_{j}}^{\text{T}}\bigg{)}^{-1}$
(34)

applying the matrix inversion Lemma to equation (9). From (7), we get

$\displaystyle\boldsymbol{\eta}_{i}=$
$\displaystyle\mathbf{P}_{i-1}\mathbf{q}_{i-1}+\frac{\mathbf{P}_{i-1}\mathbf{u}_{t_{i}}}{\lambda+\mathbf{u}_{t_{i}}^{\text{T}}\mathbf{P}_{i-1}\mathbf{u}_{t_{i}}}\big{(}s_{t_{i}}-\mathbf{u}_{t_{i}}^{\text{T}}\mathbf{P}_{i-1}\mathbf{q}_{i-1}\big{)}$
$\displaystyle=$
$\displaystyle\frac{\mathbf{P}_{i-1}}{\lambda+\mathbf{u}_{t_{i}}^{\text{T}}\mathbf{P}_{i-1}\mathbf{u}_{t_{i}}}\mathbf{q}_{i}=\mathbf{P}_{i}\mathbf{q}_{i}$

by replacing the induction hypothesis and using (34) together with the matrix
inversion Lemma. Then, the result for $\sigma_{i}=\sigma_{i}^{*}$ can be
obtained analogously to the steps in (28)-(30). ∎

##  Appendix B Proof of Theorem 2

The proof uses the following Lemma.

######  Lemma 3.

Let $N\left(x;a,b\right)$, $N\left(y;\alpha x,\beta\right)$ be two Gaussian
density functions, then

$\displaystyle N\left(x;a,b\right)N\left(y;\alpha x,\beta\right)=$
$\displaystyle N\left(x;\frac{a\beta^{2}+\alpha
yb^{2}}{\beta^{2}+\alpha^{2}b^{2}},\sqrt{\frac{b^{2}\beta^{2}}{\beta^{2}+\alpha^{2}b^{2}}}\right)N\left(y;a\alpha,\sqrt{\beta^{2}+\alpha^{2}b^{2}}\right).$

###### Proof.

$\displaystyle N\left(x;a,b\right)$ $\displaystyle N\left(y;\alpha
x,\beta\right)$ $\displaystyle=$ $\displaystyle\frac{1}{2\pi
b\beta}\exp\left\\{\frac{-\left(x-a\right)^{2}}{2b^{2}}+\frac{-\left(y-\alpha
x\right)^{2}}{2\beta^{2}}\right\\}$ $\displaystyle=$
$\displaystyle\frac{1}{2\pi
b\beta}\exp\left\\{-\frac{x^{2}-2x\frac{a\beta^{2}+\alpha
yb^{2}}{\beta^{2}+\alpha^{2}b^{2}}+\frac{a^{2}\beta^{2}+y^{2}b^{2}}{\beta^{2}+\alpha^{2}b^{2}}}{2\frac{b^{2}\beta^{2}}{\beta^{2}+\alpha^{2}b^{2}}}\right\\}.$

Then, the result is obtained since the above expression equals

$\displaystyle\frac{1}{2\pi
b\beta}\exp\left\\{-\frac{\left(x-\frac{a\beta^{2}+\alpha
yb^{2}}{\beta^{2}+\alpha^{2}b^{2}}\right)^{2}}{2\frac{b^{2}\beta^{2}}{\beta^{2}+\alpha^{2}b^{2}}}-\frac{\left(y-\alpha
a\right)^{2}}{2\left(\beta^{2}+\alpha^{2}b^{2}\right)}\right\\}$

by completing the squares. ∎

Proof of Theorem 2.

In the following, $s_{t:t+i}$ and $\mathbf{r}_{t+1:t+i}$ denote the sequences
$\left\\{s_{t},s_{t+1},...,s_{t+i}\right\\}$ and
$\left\\{\mathbf{r}_{t+1},\mathbf{r}_{t+2},...,\mathbf{r}_{t+i}\right\\}$
respectively, for any $i$.

We proceed by induction, for $i=1$ we have that

$\displaystyle p$
$\displaystyle\left(s_{t+1}|s_{t},\mathbf{r}_{t+1}\right)\propto
p\left(s_{t+1},s_{t},\mathbf{r}_{t+1}\right)$
$\displaystyle=p(\mathbf{r}_{t+1}|s_{t+1},s_{t})p(s_{t+1}|s_{t})p(s_{t})$
$\displaystyle\propto
p\left(\mathbf{r}_{t+1}|s_{t+1}\right)p\left(s_{t+1}|s_{t}\right)$ (35)
$\displaystyle\propto
N\left(s_{t+1};\mathbf{u}_{r}^{\text{T}}\boldsymbol{\eta}_{r,c},\sigma_{r,c}\right)N\left(s_{t+1};\mathbf{u}_{s}^{\text{T}}\boldsymbol{\eta}_{s,c},\sigma_{s,c}\right)$
(36)

where proportionalty relationships are due to the fact that $s_{t}$ and
$\mathbf{r}_{t+1}$ are known, (35) is obtained because the conditional
distribution of $\mathbf{r}_{t+1}$ depends only on $s_{t+1}$ since
$\\{s_{t},\mathbf{r}_{t}\\}_{t\geq 1}$ form a HMM, and (36) is obtained
because we model conditional distributions as Gaussian given by (2) and (3).

Using the previous Lemma, (36) leads to (17) with $\hat{s}_{t+1}$ and
$\hat{e}_{t+1}$ given by (18) and (19), respectively, since
$\hat{s}_{t}=s_{t}$ and $\hat{e}_{t}=0$.

If the statements hold for $i-1$, then for $i$ we have that

$\displaystyle p$
$\displaystyle\left(s_{t+i}|s_{t},\mathbf{r}_{t+1:t+i}\right)\propto
p\left(s_{t+i},s_{t},\mathbf{r}_{t+1:t+i}\right)$ $\displaystyle=$
$\displaystyle\displaystyle\int
p\left(s_{t},s_{t+i-1:t+i},\mathbf{r}_{t+1:t+i}\right)ds_{t+i-1}$ (37)
$\displaystyle=$ $\displaystyle\displaystyle\int
p\left(s_{t},s_{t+i-1:t+i},\mathbf{r}_{t+1:t+i-1}\right)p\left(\mathbf{r}_{t+i}|s_{t+i}\right)ds_{t+i-1}$
(38) $\displaystyle=$ $\displaystyle\hskip
2.84544ptp\left(\mathbf{r}_{t+i}|s_{t+i}\right)$
$\displaystyle\cdot\displaystyle\int
p\left(s_{t},s_{t+i-1},\mathbf{r}_{t+1:t+i-1}\right)p\left(s_{t+i}|s_{t+i-1}\right)ds_{t+i-1}$
(39) $\displaystyle\propto$ $\displaystyle\hskip
2.84544ptp\left(\mathbf{r}_{t+i}|s_{t+i}\right)$
$\displaystyle\cdot\displaystyle\int
p\left(s_{t+i-1}|s_{t},\mathbf{r}_{t+1:t+i-1}\right)p\left(s_{t+i}|s_{t+i-1}\right)ds_{t+i-1}$
$\displaystyle\propto$ $\displaystyle\hskip
2.84544ptN\left(s_{t+i};\mathbf{u}_{r}^{\text{T}}\boldsymbol{\eta}_{r,c},\sigma_{r,c}\right)\cdot$
$\displaystyle\displaystyle\int
N\left(s_{t+i-1};\hat{s}_{t+i-1},\hat{e}_{t+i-1}\right)N\left(s_{t+i};\mathbf{u}_{s}^{\text{T}}\boldsymbol{\eta}_{s,c},\sigma_{s,c}\right)ds_{t+i-1}$
(40)

where proportionalty relationships are due to the fact that $s_{t}$ and
$\mathbf{r}_{t+1:t+i}$ are known, (37) is obtained by marginalizing, (38) and
(39) are obtained by using the properties of HMMs, and (40) is obtained by
replacing the induction hypothesis and modelling conditional distributions as
Gaussian given by (2) and (3). Then, the result is obtained by applying the
previous Lemma to (40) twice, and substituting
$\mathbf{u}_{s}=[1,s_{t+i-1}]^{\text{T}}$. $\hfill\square$

## References

  * [1] H. Seifi and M. S. Sepasian, _Electric power system planning: issues, algorithms and solutions_. Springer Science & Business Media, 2011.
  * [2] H. L. Willis, _Power distribution planning reference book_. CRC press, 2004.
  * [3] D. K. Ranaweera, G. G. Karady, and R. G. Farmer, “Economic impact analysis of load forecasting,” _IEEE Transactions on Power Systems_ , vol. 12, no. 3, pp. 1388–1392, 1997.
  * [4] R. Weron, _Modeling and forecasting electricity loads and prices: A statistical approach_. John Wiley & Sons, 2007, vol. 403.
  * [5] D. K. Ranaweera, G. G. Karady, and R. G. Farmer, “Economic impact analysis of load forecasting,” _IEEE Transactions on Power Systems_ , vol. 12, no. 3, pp. 1388–1392, 1997.
  * [6] M. U. Fahad and N. Arbab, “Factor affecting short term load forecasting,” _Journal of Clean Energy Technologies_ , vol. 2, no. 4, pp. 305–309, 2014.
  * [7] S. K. Sheikh and M. Unde, “Short term load forecasting using ANN technique,” _International Journal of Engineering Sciences & Emerging Technologies_, vol. 1, no. 2, pp. 97–107, 2012.
  * [8] J. W. Taylor and R. Buizza, “Neural network load forecasting with weather ensemble predictions,” _IEEE Transactions on Power Systems_ , vol. 17, no. 3, pp. 626–632, 2002.
  * [9] M. Espinoza, J. A. Suykens, R. Belmans, and B. De Moor, “Electric load forecasting,” _IEEE Control Systems Magazine_ , vol. 27, no. 5, pp. 43–57, 2007.
  * [10] D. Jang, J. Eom, M. J. Park, and J. J. Rho, “Variability of electricity load patterns and its effect on demand response: A critical peak pricing experiment on Korean commercial and industrial customers,” _Energy Policy_ , vol. 88, pp. 11–26, 2016.
  * [11] G. Huebner, D. Shipworth, I. Hamilton, Z. Chalabi, and T. Oreszczyn, “Understanding electricity consumption: A comparative contribution of building factors, socio-demographics, appliances, behaviours and attitudes,” _Applied Energy_ , vol. 177, pp. 692–702, 2016.
  * [12] T. Kronenberg, “The impact of demographic change on energy use and greenhouse gas emissions in Germany,” _Ecological Economics_ , vol. 68, no. 10, pp. 2637–2645, 2009.
  * [13] G. Ditzler, M. Roveri, C. Alippi, and R. Polikar, “Learning in nonstationary environments: A survey,” _IEEE Computational Intelligence Magazine_ , vol. 10, no. 4, pp. 12–25, 2015.
  * [14] N. Amjady, F. Keynia, and H. Zareipour, “Short-term load forecast of microgrids by a new bilevel prediction strategy,” _IEEE Transactions on Smart Grid_ , vol. 1, no. 3, pp. 286–294, 2010.
  * [15] F. Javed, N. Arshad, F. Wallin, I. Vassileva, and E. Dahlquist, “Forecasting for demand response in smart grids: An analysis on use of anthropologic and structural data and short term multiple loads forecasting,” _Applied Energy_ , vol. 96, pp. 150–160, 2012.
  * [16] W. Kong, Z. Y. Dong, D. J. Hill, F. Luo, and Y. Xu, “Short-term residential load forecasting based on resident behaviour learning,” _IEEE Transactions on Power Systems_ , vol. 33, no. 1, pp. 1087–1088, 2017.
  * [17] E. Gelenbe and Y. Caseau, “The impact of information technology on energy consumption and carbon emissions,” _Ubiquity_ , vol. 2015, no. June, pp. 1–15, 2015.
  * [18] R. Baran and L. F. L. Legey, “The introduction of electric vehicles in Brazil: Impacts on oil and electricity consumption,” _Technological Forecasting and Social Change_ , vol. 80, no. 5, pp. 907–917, 2013.
  * [19] T. Shiina and J. R. Birge, “Stochastic unit commitment problem,” _International Transactions in Operational Research_ , vol. 11, no. 1, pp. 19–32, 2004.
  * [20] T. Hong and S. Fan, “Probabilistic electric load forecasting: A tutorial review,” _International Journal of Forecasting_ , vol. 32, no. 3, pp. 914–938, 2016.
  * [21] S. Shalev-Shwartz _et al._ , “Online learning and online convex optimization,” _Foundations and Trends in Machine Learning_ , vol. 4, no. 2, pp. 107–194, 2012.
  * [22] S. S. Pappas, L. Ekonomou, D. C. Karamousantas, G. Chatzarakis, S. Katsikas, and P. Liatsis, “Electricity demand loads modeling using auto regressive moving average (ARMA) models,” _Energy_ , vol. 33, no. 9, pp. 1353–1360, 2008.
  * [23] S.-J. Huang and K.-R. Shih, “Short-term load forecasting via ARMA model identification including non-Gaussian process considerations,” _IEEE Transactions on Power Systems_ , vol. 18, no. 2, pp. 673–679, 2003.
  * [24] N. Amjady, “Short-term hourly load forecasting using time-series modeling with peak load estimation capability,” _IEEE Transactions on Power Systems_ , vol. 16, no. 3, pp. 498–505, 2001.
  * [25] N. Mohamed, M. H. Ahmad, and Z. Ismail, “Short term load forecasting using double seasonal ARIMA model,” in _Proceedings of the Regional Conference on Statistical Sciences_ , vol. 10, 2010, pp. 57–73.
  * [26] Y. Chakhchoukh, P. Panciatici, and L. Mili, “Electric load forecasting based on statistical robust methods,” _IEEE Transactions on Power Systems_ , vol. 26, no. 3, pp. 982–991, 2010.
  * [27] K. Chen, K. Chen, Q. Wang, Z. He, J. Hu, and J. He, “Short-term load forecasting with deep residual networks,” _IEEE Transactions on Smart Grid_ , vol. 10, no. 4, pp. 3943–3952, 2018.
  * [28] H. S. Hippert, C. E. Pedreira, and R. C. Souza, “Neural networks for short-term load forecasting: A review and evaluation,” _IEEE Transactions on Power Systems_ , vol. 16, no. 1, pp. 44–55, 2001.
  * [29] B.-J. Chen and M.-W. Chang, “Load forecasting using support vector machines: A study on EUNITE competition 2001,” _IEEE Transactions on Power Systems_ , vol. 19, no. 4, pp. 1821–1830, 2004.
  * [30] S. Fan and L. Chen, “Short-term load forecasting based on an adaptive hybrid method,” _IEEE Transactions on Power Systems_ , vol. 21, no. 1, pp. 392–401, 2006.
  * [31] Y. Wang, Q. Xia, and C. Kang, “Secondary forecasting based on deviation analysis for short-term load forecasting,” _IEEE Transactions on Power Systems_ , vol. 26, no. 2, pp. 500–507, 2010.
  * [32] T. Senjyu, P. Mandal, K. Uezato, and T. Funabashi, “Next day load curve forecasting using hybrid correction method,” _IEEE Transactions on Power Systems_ , vol. 20, no. 1, pp. 102–109, 2005.
  * [33] J. Lourenço and P. Santos, “Short term load forecasting using Gaussian process models,” _Proceedings of Instituto de Engenharia de Sistemas e Computadores de Coimbra_ , 2010.
  * [34] B. Liu, J. Nowotarski, T. Hong, and R. Weron, “Probabilistic load forecasting via quantile regression averaging on sister forecasts,” _IEEE Transactions on Smart Grid_ , vol. 8, no. 2, pp. 730–737, 2015.
  * [35] Y. Wang, Q. Chen, N. Zhang, and Y. Wang, “Conditional residual modeling for probabilistic load forecasting,” _IEEE Transactions on Power Systems_ , vol. 33, no. 6, pp. 7327–7330, 2018.
  * [36] L. D. Paarmann and M. D. Najar, “Adaptive online load forecasting via time series modeling,” _Electric Power Systems Research_ , vol. 32, no. 3, pp. 219–225, 1995.
  * [37] A. Laouafi, M. Mordjaoui, S. Haddad, T. E. Boukelia, and A. Ganouche, “Online electricity demand forecasting based on an effective forecast combination methodology,” _Electric Power Systems Research_ , vol. 148, pp. 35–47, 2017.
  * [38] A. Ba, M. Sinn, Y. Goude, and P. Pompey, “Adaptive learning of smoothing functions: Application to electricity load forecasting,” in _Advances in Neural Information Processing Systems_ , 2012, pp. 2510–2518.
  * [39] T. Hong, P. Pinson, S. Fan, H. Zareipour, A. Troccoli, and R. J. Hyndman, “Probabilistic energy forecasting: Global energy forecasting competition 2014 and beyond,” _International Journal of Forecasting_ , vol. 32, no. 3, pp. 896–913, 2016.
  * [40] T. Gneiting, F. Balabdaoui, and A. E. Raftery, “Probabilistic forecasts, calibration and sharpness,” _Journal of the Royal Statistical Society: Series B (Statistical Methodology)_ , vol. 69, no. 2, pp. 243–268, 2007.
  * [41] S. Mazuelas, Y. Shen, and M. Z. Win, “Belief condensation filtering,” _IEEE Transactions on Signal Processing_ , vol. 61, no. 18, pp. 4403–4415, 2013.
  * [42] O. Cappé, E. Moulines, and T. Rydén, _Inference in hidden Markov models_. Springer Science & Business Media, 2006.
  * [43] W. Liu, J. C. Principe, and S. Haykin, _Kernel adaptive filtering: a comprehensive introduction_. John Wiley & Sons, 2011, vol. 57.
  * [44] M. Espinoza, J. A. Suykens, and B. De Moor, “Load forecasting using fixed-size least squares support vector machines,” in _Computational Intelligence and Bioinspired Systems_ , vol. 3512. Springer, 2005, pp. 1018–1026.
  * [45] T. Hong, P. Pinson, and S. Fan, “Global energy forecasting competition 2012,” _International Journal of Forecasting_ , vol. 30, no. 2, pp. 357–363, 2014.
  * [46] A. P. Liavas and P. A. Regalia, “On the numerical stability and accuracy of the conventional recursive least squares algorithm,” _IEEE Transactions on Signal Processing_ , vol. 47, no. 1, pp. 88–96, 1999.
  * [47] A. Rudi and L. Rosasco, “Generalization properties of learning with random features,” in _Advances in Neural Information Processing Systems_ , 2017, pp. 3215–3225.

  *[LR]: linear regression
  *[SVM]: support vector machine
  *[GP]: Gaussian process
  *[QR]: quantile regression
  *[APLF]: adaptive probabilistic load forecasting
  *[HMM]: hidden Markov model
  *[HMMs]: hidden Markov model
  *[SVMs]: support vector machine
  *[AFF]: adaptive forgetting factor
  *[SFDA]: secondary forecasting based on deviation analysis


EPJ Web of Conferences Lattice2017  11institutetext: CSSM, Department of
Physics, The University of Adelaide, Adelaide, SA, Australia 5005
22institutetext: RIKEN Advanced Institute for Computational Science, Kobe,
Hyogo 650-0047, Japan 33institutetext: UniversitÃt Hamburg, Regionales
Rechenzentrum, 20146 Hamburg, Germany

# An update on the BQCD Hybrid Monte Carlo program

Taylor Ryan Haar 11 Yoshifumi Nakamura 22 Hinnerk StÃben Speaker, 33
<EMAIL_ADDRESS>
###### Abstract

We present an update of BQCD, our Hybrid Monte Carlo program for simulating
lattice QCD. BQCD is one of the main production codes of the QCDSF
collaboration and is used by CSSM and in some Japanese finite temperature and
finite density projects. Since the first publication of the code at Lattice
2010 the program has been extended in various ways. New features of the code
include: dynamical QED, action modification in order to compute matrix
elements by using Feynman-Hellman theory, more trace measurements (like
${\mathrm{T}r}(D^{-n})$ for $\kappa$, $c_{\mathrm{SW}}$ and chemical potential
reweighting), a more flexible integration scheme, polynomial filtering, term-
splitting for RHMC, and a portable implementation of performance critical
parts employing SIMD.

##  1 Introduction

BQCD is a Hybrid Monte Carlo program for simulating lattice QCD with dynamical
Wilson fermions. It was first published at Lattice 2010 Nakamura:2010qh  and
has been used by several groups: the QCDSF-UKQCD collaboration Stuben:2000wm ;
Gockeler:2007rm ; Cundy:2009yy ; Horsley:2015vla ; Chambers:2014qaa ;
Chambers:2015bka , CSSM Chambers:2014qaa ; Chambers:2015bka ; Haar:2016bwe ;
Hollitt:2017tkf , Japanese finite density Jin:2013wta ; Jin:2015taa  and
finite temperature Jin:2017jjp  projects, and the RQCD collaboration
Bali:2017pdv .

Here we report on extensions and optimizations that were made meanwhile and
give an update on compute performance. The code and a manual can be downloaded
from bqcd-download . New features of the program are:

  * •

_Actions:_ hopping term with chemical potential, clover $O(a)$ improved Wilson
action plus a CPT breaking term, QCD+QED, QCD+Axion. See section 2.

  * •

_Algorithms:_ polynomial filtering, a generalized multiscale integration
scheme, truncated RHMC, Zolotarev approximation. See section 4.

  * •

_Compute performance optimizations_ : explicit vectorization with SIMD
intrinsics, improvement of MPI communication. See section 5.

##  2 Actions

###  2.1 Gauge actions

Implemented are the Wilson gauge action and an improved gauge action, see
Nakamura:2010qh .

###  2.2 Fermion actions

At the time of Nakamura:2010qh  the program could simulate the standard Wilson
fermion action $S^{\mathrm{Wilson}}_{F}$, $O(a)$ clover improved Wilson
fermions, and the SLiNC fermion action Cundy:2009yy . The new version can also
simulate:

  * •

the hopping term with chemical potential $\mu$ Jin:2013wta ; Jin:2015taa

$\begin{split}S_{F}=\sum_{x}\Bigg{\\{}\bar{\psi}(x)\psi(x)&-\kappa\sum_{i}^{3}\left[\bar{\psi}(x)U_{i}^{\dagger}(x-\hat{i})(1+\gamma_{i})\psi(x-\hat{i})+\bar{\psi}(x)U_{i}(x)(1-\gamma_{i})\psi(x+\hat{i})\right]\\\
&-\kappa\left[\bar{\psi}(x)U_{4}^{\dagger}(x-\hat{4})(1+\gamma_{4})e^{-\mu}\psi(x-\hat{4})+\bar{\psi}(x)U_{4}(x)(1-\gamma_{4})e^{\mu}\psi(x+\hat{4})\right]\Bigg{\\}}\end{split}$
(1)

  * •

the clover $O(a)$ improved Wilson action plus a CPT breaking term with
coefficient $\lambda$ and a $4\times 4$ matrix $H$ Chambers:2014qaa ;
Chambers:2015bka

$S_{F}=S^{\mathrm{Wilson}}_{F}-\frac{i}{2}\kappa\,c_{\mathrm{S}W}\sum_{x}\big{[}\bar{\psi}(x)\sigma_{\mu\nu}F_{\mu\nu}(x)\psi(x)+\kappa\,\lambda\bar{\psi}(x)H\psi(x)\big{]}$
(2)

###  2.3 QCD+QED

The program can simulate QCD+QED Horsley:2015vla  using the action

$S=S_{G}+S_{A}+\sum_{q}S_{F}^{q}\;.$ (3)

$S_{G}$ is an SU(3) gauge action, $S_{A}$ is the non-compact U(1) gauge action

$S_{A}=\frac{\beta_{\mathrm{QED}}}{2}\sum_{x,\mu<\nu}\left[A_{\mu}(x)+A_{\nu}(x+\hat{\mu})-A_{\mu}(x+\hat{\nu})-A_{\nu}(x)\right]^{2}\;,$
(4)

and the fermion action for flavour $q$ is

$\begin{split}S_{F}^{q}&=\sum_{x}{\Bigg{\\{}}\kappa_{q}\sum_{\mu}\left[\overline{q}(x)(\gamma_{\mu}-1)e^{-iQ_{q}A_{\mu}(x)}\tilde{U}_{\mu}(x)q(x+\hat{\mu})\right.\\\
&-\left.\overline{q}(x)(\gamma_{\mu}+1)e^{iQ_{q}A_{\mu}(x-\hat{\mu})}\tilde{U}^{\dagger}_{\mu}(x-\hat{\mu})q(x-\hat{\mu})\right]+\;\overline{q}(x)q(x)-\frac{1}{2}\kappa_{q}c_{\mathrm{S}W}\sum_{\mu,\nu}\overline{q}(x)\sigma_{\mu\nu}F_{\mu\nu}(x)q(x)\Bigg{\\}}\;,\end{split}$
(5)

where $Q_{u}=+2/3$, $Q_{d}=Q_{s}=-1/3$ and $\tilde{U}_{\mu}$ is a singly
iterated stout link.

###  2.4 QCD+Axion

The program can simulate QCD+Axion thiscontrib235  using the action

$S=S_{G}+S_{a}+\sum_{q}S_{F}^{q}\;.$ (6)

$S_{a}$ is scalar action for the axion field $\phi_{a}$

$S_{a}=\kappa_{a}\sum_{x}\sum_{\mu}\Big{(}\phi_{a}(x)-\phi_{a}(x+\mu)\Big{)}\phi_{a}(x)\,,\\\
$ (7)

and the fermion action for flavour $q$ in the case of Wilson fermions is

$\begin{split}S_{F}^{q}&=\sum_{x}\bar{q}(x)\Big{[}1+(\kappa_{q}\lambda_{q}+f_{\mathrm{inv}}\phi_{a})\gamma_{5}\Big{]}q(x)\\\
&-\kappa_{q}\sum_{x,\mu}\Big{[}\bar{q}(x)(1-\gamma_{\mu})U_{\mu}(x)q(x+a\hat{\mu})+\bar{q}(x-a\hat{\mu})(1+\gamma_{\mu})U_{\mu}^{\dagger}(x-a\hat{\mu})q(x)\Big{]}\,,\end{split}$
(8)

where

$\kappa_{q}={1\over 2am_{q}+8}\,,\quad\lambda_{q}=i2am_{q}{\theta\over
N_{f}}\,,\quad f_{\mathrm{inv}}=i2\kappa_{q}m_{q}{\sqrt{\kappa_{a}}\over
f_{a}N_{f}}\,.$ (9)

##  3 Measurements

The following quantities can be measured with BQCD: plaquettes (quadratic and
rectangular), topological charge (cooling method), Polyakov loop, Wilson flow,
traces of the fermion matrix ($\mathrm{Tr}\,(M^{-1})$,
$\mathrm{Tr}\,(\gamma_{5}M^{-1})$, $\mathrm{Tr}\,(M^{\dagger}M)^{-1}$), quark
determinant with chemical potential, smallest and largest eigenvalue of the
Dirac matrix, meson and baryon propagators.

##  4 Algorithms

In addition to nested integrators for multiple time scales, a generalized
integration scheme Haar:2016bwe  has been implemented. This is where separate
integration schemes for each action term are superimposed onto a single time
step evolution. This allows the integration step-sizes for each action term to
be completely independent of the others.

Rational Hybrid Monte Carlo (RHMC) is implemented with rational approximations
from the Remez algorithm. In the case of approximating
$(W^{\dagger}W)^{-1/2}$, an alternative rational function is available, namely
the Zolotarev optimal rational approximation (see e.g. Chiu:2002eh  for an
explanation).

Alongside Hasenbusch filtering, there are two new filtering methods, one of
which applies exclusively to RHMC:

  * •

_Polynomial filtering_ applies to both RHMC and standard HMC fermions, and is
the application of a polynomial filter $P(W^{\dagger}W)$ to split the fermion
action into several terms:

$S_{PF}=\overline{\psi}_{1}P(W^{\dagger}W)\psi_{1}+\overline{\psi}_{2}P(W^{\dagger}W)^{-1}W^{\dagger}W\psi_{2}.$
(10)

  * •

_Term-splitting for RHMC_ splits the sum in the rational approximation
$R(W^{\dagger}W)$ for RHMC into several terms, giving action

$S_{tRHMC}=\overline{\psi}_{1}R_{1,t}(W^{\dagger}W)\psi_{1}+\overline{\psi}_{2}R_{t+1,N}(W^{\dagger}W)\psi_{2}$
(11)

where

$R_{i,j}(K)=c_{n}^{\delta_{i1}}\sum_{k=i}^{j}\frac{W^{\dagger}W+a_{k}}{W^{\dagger}W+b_{k}},$
(12)

$a_{k},b_{k}$ are ordered decreasing.

BQCD has a wide range of iterative solvers: _cg_ , BiCGstab, GMRES, GCRODR,
multishift _cg_ , block multishift _cg_.

##  5 Optimization of compute performance

###  5.1 SIMD vectorization and MPI

In addition to parallelization with MPI and OpenMP a third level of parallel
implementation was introduced for solvers: SIMD vectorization with SIMD
intrinsic functions. The SIMD implementation is generic, i.e. it works for any
size of SIMD vectors. In order to achieve this, the data layout of arrays had
to be changed. All arrays (for gauge, spin-colour and clover fields) now have
SIMD vectors as the smallest structure. In Fortran notation the gauge field is
defined in the following way

> old: | `complex(8) :: u(3, 3, volume/2)`
> ---|---
> new: | `real(8) :: u(SIMDsize, re:im, 3, 3, volume/2 / SIMDsize)`

and the new layout of the spin-color field is

> old: | `complex(8) :: a(4, 3, volume/2)`
> ---|---
> new: | `real(8) :: a(SIMDsize, re:im, 2, 3, volume/2 / SIMDsize, 2)`

where the $4$ spin components of the spin-colour field are split into $2+2$
components which optimizes MPI communication in $t$-direction. The clover
arrays, for which a packed format is used, were changed accordingly.

At the single core level the SIMD code is about 2 times faster than the
corresponding Fortran code. With this speed-up computations are increasingly
dominated by communication and improvement of MPI communication becomes
important. Hence, the following MPI optimizations where made:

  * •

The overhead introduced by the reduction to two-component spinors was
minimized. Previously the projection was done for the whole local volume, now
it is only done for boundary sites, and there is no projection in the
$t$-direction needed any more.

  * •

All MPI ’buffers’ are consecutive in memory and aligned to SIMD vector
boundaries.

  * •

Communication can overlap with computation. This is implemented with MPI plus
OpenMP, where the _master thread_ communicates while the other threads
compute.

In tables 1 and 2 performance figures are listed for machines and lattices
that are currently used in production. The optimized code runs between 1.3 and
1.7 times faster.

Table 1: Double precision performance of the _cg_ -solver of BQCD on a Cray XC40 (24 cores per node). |  $48^{3}\times 96$ lattice |  $64^{3}\times 96$ lattice
---|---|---
| Fortran | SIMD | Fortran | SIMD
| per core | overall | per core | overall | per core | overall | per core | overall
#cores | Mflop/s | Tflop/s | Mflop/s | Tflop/s | Mflop/s | Tflop/s | Mflop/s | Tflop/s
1536 |  930 |  1.4 | 1557 |  2.4 | 785 |  1.2 | 1028 |  1.6
3072 |  949 |  2.9 | 1516 |  4.7 | 795 |  2.4 | 1231 |  3.8
6144 | 1222 |  7.5 | 1558 |  9.6 | 876 |  5.4 | 1419 |  8.7
9216 | 1253 | 11.5 | 1775 | 16.4 | – | – | – | –
12288 | 1274 | 15.7 | 1678 | 20.6 | 927 | 11.4 | 1572 | 19.3
Table 2: Double precision performance of the _cg_ -solver of BQCD on an IBM BlueGene/Q (8192 cores, a _midplane_ , is the smallest partition that has a fully wired torus network and with our SIMD implementation the largest possible partition for the $48^{3}\times 96$ lattice, where the lattice volume per core is $48\times 3^{3}$). |  $48^{3}\times 96$ lattice |  $64^{3}\times 96$ lattice
---|---|---
| Fortran | SIMD | Fortran | SIMD
| per core | overall | per core | overall | per core | overall | per core | overall
#cores | Mflop/s | Tflop/s | Mflop/s | Tflop/s | Mflop/s | Tflop/s | Mflop/s | Tflop/s
8192 | 539 | 4.4 | 912 | 7.5 | 525 |  4.3 | 752 |  6.2
16384 | – | – | – | – | 596 |  9.8 | 783 | 12.8
32768 | – | – | – | – | 503 | 16.5 | 771 | 25.3

###  5.2 QUDA

BQCD can run on GPUs by employing the QUDA library Clark:2009wm . QUDA has a
BQCD interface to its _cg_ and multishift _cg_ solvers.

## Acknowledgements

We would like to thank Gerrit Schierholz, Roger Horsley, Waseem Kamleh, Paul
Rakow and James Zanotti for support, stimulating discussions and bug reports.
The computations were performed on a Cray XC40 of the North-German
Supercomputing Alliance (HLRN) and the IBM BlueGene/Q at JÃlich Supercomputer
Centre (JSC).

## References

  * (1) Y. Nakamura, H. StÃben, PoS LATTICE2010, 040 (2010), 1011.0199
  * (2) H. StÃben (UKQCD, QCDSF), Nucl. Phys. Proc. Suppl. 94, 273 (2001), hep-lat/0011045
  * (3) M. GÃckeler et al. (QCDSF), PoS LAT2007, 041 (2007), 0712.3525
  * (4) N. Cundy et al., Phys. Rev. D79, 094507 (2009), 0901.3302
  * (5) R. Horsley et al., JHEP 04, 093 (2016), 1509.00799
  * (6) A.J. Chambers et al. (QCDSF/UKQCD, CSSM), Phys. Rev. D90, 014510 (2014), 1405.3019
  * (7) A.J. Chambers et al., Phys. Rev. D92, 114517 (2015), 1508.06856
  * (8) T. Haar, W. Kamleh, J. Zanotti, Y. Nakamura, Comput. Phys. Commun. 215, 113 (2017), 1609.02652
  * (9) S. Hollitt, P. Jackson, R. Young, J. Zanotti, PoS INPC2016, 272 (2017)
  * (10) X.Y. Jin, Y. Kuramashi, Y. Nakamura, S. Takeda, A. Ukawa, Phys. Rev. D88, 094508 (2013), 1307.7205
  * (11) X.Y. Jin, Y. Kuramashi, Y. Nakamura, S. Takeda, A. Ukawa, Phys. Rev. D92, 114511 (2015), 1504.00113
  * (12) X.Y. Jin, Y. Kuramashi, Y. Nakamura, S. Takeda, A. Ukawa, Phys. Rev. D96, 034523 (2017), 1706.01178
  * (13) G.S. Bali, S. Collins, A. Cox, A. SchÃfer, Phys. Rev. D96, 074501 (2017), 1706.01247
  * (14) https://www.rrz.uni-hamburg.de/bqcd
  * (15) G. Schierholz, Y. Nakamura, _Dynamical QCD+Axion simulation: First results_ , in _Proceedings,35th International Symposium on Lattice Field Theory (Lattice2017): Granada, Spain_, to appear in EPJ Web Conf.
  * (16) T.W. Chiu, T.H. Hsieh, C.H. Huang, T.R. Huang, Phys. Rev. D66, 114502 (2002), hep-lat/0206007
  * (17) M.A. Clark, R. Babich, K. Barros, R.C. Brower, C. Rebbi, Comput. Phys. Commun. 181, 1517 (2010), 0911.3191