Computation of Ultra-Short-Term Prediction Intervals of the Power Prosumption in Active Distribution Networks

Plouton Grammatikos,
Fabrizio Sossan,
Jean-Yves Le Boudec,
and Mario Paolone
F. Sossan is at the Haute École spécialisée de Suisse occidentale Sion, P. Grammatikos, M. Paolone (Distributed Electrical Systems Laboratory) and J.Y. Le Boudec (Laboratory for Communications and Applications) are at the École Polytechnique Fédérale de Lausanne, Switzerland. Email<EMAIL_ADDRESS>{plouton.grammatikos, mario.paolone<EMAIL_ADDRESS>
Microgrids and, in general, active distribution networks require ultra-short-term prediction, i.e., for sub-second time scales, for specific control decisions. Conventional forecasting methodologies are not effective at such time scales. To address this issue, we propose a non-parametric method for computing ultra short-term prediction intervals (PIs) of the power prosumption of generic electrical-distribution networks. The method groups historical observations into clusters according to the values of influential variables. It is applied either to the original or to the differentiated power-prosumption time series. The clusters are considered statistically representative pools of future realizations of power prosumption (or its derivative). They are used to determine empirical PDFs and, by extracting the quantiles, to deliver PIs for respective arbitrary confidence levels. The models are validated a posteriori by carrying out a performance analysis that uses experimentally observed power-prosumption for different building types, thus allowing the identification of the dominant model.

prosumption, forecast, prediction intervals, electrical load, microgrids.

§ INTRODUCTION

After being the mainstream framework for the integration and coordination of distributed generation, the concepts of an active distribution network (ADN) and a microgrid recently came to prominence to tackle the challenges caused by the large-scale integration of variable renewable generation.
ADNs comprise low-voltage (LV) or medium-voltage (MV) electrical grids with systems in place to control a combination of distributed energy resources (DERs), such as generators, loads, and storage devices <cit.>.

Due to the low level of aggregation, the ADN requirements for electrical-power prosumption[In power systems, the prosumption indicates the aggregated power provided or consumed by users that have the capability to generate electricity by means of user-owned distributed generation locally.] forecasting are different than for conventional large interconnected grids.
One example relates to the possible violation of the ampacity rating of transformers, power converters, and lines due to sudden changes in the prosumption associated with the highly stochastic nature of prosumers[A node that can both absorb or inject power due to prosumption.]. A significant change in the prosumers' renewable-power generation, or a spike in load, can create power-flow variations that can exceed the rating of transformers. Whereas, a spike in current could cause the line relays to trip. This example is particularly relevant in the presence of photovoltaics (PV) and electrical-vehicle (EV) charging stations (CS) in the grid (e.g., <cit.>). The former can exhibit power variations of even 60% of their capacity in under a second <cit.>, whereas the latter can cause significant load changes of hundreds of kW within a few seconds <cit.>. The prediction of the prosumption can be integrated into various real-time (RT) (e.g., <cit.>) and model-predictive control (MPC) (e.g., <cit.>) frameworks to ensure the safe operation of the ADNs and the optimal usage of their resources.

A further example is based on the capability of microgrids to operate autonomously. When connected to the external grid, microgrids can provide ancillary services to the upper grid layer <cit.>. Whereas, in case of contingencies, they can operate islanded to the main grid and can enhance the resiliency of the supply to the local load.
The islanding maneuver (i.e., the operation sequence for bringing a microgrid from connected to off-grid) requires a prediction of the prosumption in the range of the fundamental frequency period (20 ms) in order to correctly set the gains of the slack resource droop control (e.g., <cit.>).

A third example, which is of high importance both in distribution and transmission networks, is that of voltage sags <cit.>. Voltage sags are defined as a sudden reduction of the voltage between 90% and 10% of the nominal value and can last from 10 ms up to 1 minute. They are caused primarily by power-system faults, such as short circuits, or by the start-up of large motors and can cause system outages if not treated in time. For distribution systems, where voltage sags typically last between 90-2000 ms <cit.>, in order to mitigate the voltage drop, a control framework equipped with a forecasting tool acting in the sub-second range could have a timely reaction to the voltage sag by injecting an optimally computed active/reactive power into the grid.

As ultra-short-term power-prosumption forecasts are actionable for fundamental decisions in the context of ADN/microgrid operation and their RT control, we note that well-established forecasting methodologies (e.g., developed using several techniques, such as regression-based model, artificial neural network <cit.>) are not suited to this purpose because, besides referring to point predictions, they were developed considering a high level of aggregation and forecasting horizon from 15 minutes and up. Moreover, in certain applications, such as robust optimization (e.g., <cit.>), worst-case analysis is required; therefore, point predictions are inadequate.

As stated in <cit.> and further supported by the experimental measurements of this paper, when decreasing the aggregation level and measurements sampling time, the power-prosumption volatility and noise level become prominent because consumer behaviors tend not to cancel out.

As the current state of the art appears to be inadequate to deliver ultra-short-term power-prosumption prediction intervals (PIs) of ADN prosumers, we propose an adaptive non-parametric method based on pattern recognition. The algorithm is designed to be computationally efficient, thus allowing for the delivery of high-time resolution probabilistic PIs in RT and at a high sampling rate with low computational overhead. The model is initially trained using a time series of the aggregated-power prosumption without requiring any knowledge of the nature and number of loads/generators present in the network. With an efficient updating and aging procedure, it is then continuously updated as new measurements become available.

The paper is organized as follows. In sec:statement, the problem of estimating the power prosumption PI is stated, along with a review of methods already developed in the context of power system applications. In Section <ref>, we describe the proposed PI models whose performance is analyzed in sec:results by using experimental data for different building types. Finally,  in sec:conclusions, we summarize the findings.

§ PROBLEM STATEMENT

As stated in the previous section, most of the existing literature on power-prosumption forecasting is concerned with point predictions, specifically the problem of estimating the expected realization of the power prosumption for a given look-ahead time. Whereas, we target the computation of PIs; in other words, we predict, with a given confidence level, the interval where the future power-prosumption realization is expected to lie.
Denoting the PI at the target confidence level $\alpha$ as the couple $({P}^{\downarrow\alpha},{P}^{\uparrow\alpha})$ composed by the lower and upper bound of the interval, we address the problem that consists in finding the one-step-ahead PI as a function of a sequence of $n$ historical power-prosumption measurements until the time instant $i$, specifically:

\begin{align}
\left( {P}^{\downarrow\alpha}_{i+1|i}, {P}^{\uparrow\alpha}_{i+1|i} \right) = f(P_{i}, \dots, P_{i+1-n})
\end{align}

where $i$ is the current time interval, and $f$ is a PI estimation model.

When using parametric point predictors (such as autoregressive integrated moving average (ARIMA) models), we can determine the PIs by estimating the variance of the model residuals and computing the quantiles for the prescribed confidence level. This procedure can be performed under the hypothesis of Gaussian iid (independent and identically distributed) model residuals.

In cases where this hypothesis does not hold, non-parametric methods could be considered. For example, in order to determine the PIs of the power output of a wind farm, the authors of <cit.> apply quantile regression to characterize the historical residuals of a state-of-the-art point-prediction model.

The same concept, but developed using fuzzy inference instead of quantile regression, is described in <cit.>. As far as forecasting the electrical-power prosumption is concerned, non-parametric methods have been proposed lately in <cit.> and <cit.>.
In the former work, an artificial neural network (ANN) with an empirically chosen number of layers was trained using historical data to provide a 30-minute-ahead PI for a given confidence level according to the values of selected data features.
The latter work is concerned with predicting the minimum and maximum bounds of the power consumption  by applying empirical-mode decomposition and support-vector regression to an interval-valued time signal obtained from a one-hour historical sample of power consumption. Both methods target a prediction horizon that is too long for the requirements of RT ADN/microgrid operation discussed in the introduction.

§ COMPUTATION OF PREDICTION INTERVALS

As it will be exhaustively described in the rest of this section, the estimation methodology for PIs consists in grouping historical power-prosumption measurements into clusters according to the value of selected influential variables. At the time of delivering a PI, the values of the influential variables are determined, thus allowing for the selection of the appropriate cluster that is finally used as the empirical PDF (probability distribution function) of the next realization. The algorithm is designed to deliver a PI in rolling RT with a minimum report rate of 20ms (longer values are also analyzed).

The algorithm operation sequence is sketched in fig:stages. The first phase, called batch training, consists of the off-line training of the estimation model that uses historical data. In the second phase, the one-step-ahead PI is delivered. And finally, in the third phase, the new progressively available measurements are used for the on-line model training. This methodology is applied in two flavors, specifically on the original and differentiated time-series, as discussed in the following two sections.

block = [rectangle, draw, fill=blue!0,
text width=15em, text centered, rounded corners, minimum height=3em, text width=8em]
block1 = [rectangle, draw, fill=blue!5,
text width=15em, text centered, minimum height=2em]
input = [ellipse, draw, fill=red!5, text width=5em, text centered, minimum height=3em, node distance=3.5cm]
output = [ellipse, draw, fill=green!5, text width=5em, text centered, minimum height=3em, node distance=3.5cm]
line = [draw, -latex']

[node distance = 1.25cm, auto]

[block] (init) $i = $ initial time slot index;
[block1, below of=init, node distance=1cm] (batch) Off-line batch training;

[below of=init, node distance=1.5cm] (loop) ;
[block1, below of=batch, node distance=1.0cm] (estimation) On-line PI Estimation;

[block, below of=estimation, node distance=1.0cm] (passing) Rolling index $i=i+1$;
[block1, below of=passing, node distance=1.0cm] (online) On-line training;
[input, left of=batch] (AA) Historical data;
[output, left of=estimation] (AB) A PI is produced;
[input, left of=online] (AC) New realization;

[line] (init) – (batch);
[line] (batch) – (estimation);
[line] (estimation) – (passing);
[line] (passing) – (online);
[line] (online) – +(2.5,0) – (+2.5,0 |- loop) – (loop);

[line] (AA) – (batch);
[line] (estimation) – (AB);
[line] (AC) – (online);

Operation sequence of the PI estimation models. The batch training is performed off-line, whereas the PIs computation and on-line training are performed in rolling RT. Ellipses denote the input and output of each phase.

§.§ PI Model A

§.§.§ Off-Line Batch Training

We consider $n$ historical power-prosumption measurements $P_{i+1-n}, \dots, P_{i}$ and the respective time stamp $t_{i+1-n}, \dots, t_{i}$. Power-prosumption measurements are discretized in time and amplitude with user-defined discretization steps $\Delta T$ and $\Delta P$, respectively.
We define $c$ as the operator that determines a non-negative integer $l$, said label, according to the power-prosumption value and timestamp in the following way:

\begin{align}
l_j = c(P_j, t_j), \;\; l_j \in \{0,1,\dots,L-1\}\label{eq:classification}.
\end{align}

where $j$ is a generic time slot, and $L$ is the total number of labels.
The process used to operate the classification $c$ is described in detail in sec:classification.
The label in (<ref>) is utilized to group the historical power-prosumption measurements $\mathcal{P}_{i,n}$ into $L$ clusters, denoted as $\mathcal{G}^0_i, \dots, \mathcal{G}^{L-1}_i$. Each cluster contains all the historical power-prosumption measurements for which the previous observation was of the respective given label. For example, the cluster $\mathcal{G}^0_i$ contains the measurements until the time slot $i$ for which the respective previous realization was with label $0$, $\mathcal{G}^1_i$ those for which the respective previous realization was with label 1, and so on. Formally, the clusters are defined as:

\begin{align}
\begin{aligned}
\mathbb{G}^l_i = \{ P_{j+1}:c(P_{j}, t_{j}) = l, \; j=i-n,\dots,i\},\\
\end{aligned}\label{eq:clusters}
\end{align}

Let $\mathcal{G}^0_i, \dots, \mathcal{G}^{L-1}_i$ be the normalized histogram of each cluster computed as

\begin{align}
\mathcal{G}^l_i(x) = \frac{1}{|\mathbb{G}^l_i|}\sum_{p \in \mathbb{G}^l_i} \delta(x-p), \;\; l=0,\dots, L-1\label{eq:histograms}
\end{align}

where $|\cdot|$ denotes the set cardinality (, the number of elements it contains) and $\delta$ is the Dirac measure:

\begin{align}
\delta(x) = \begin{cases}
1 & x=0 \\
0 & \text{otherwise}.
\end{cases}
\end{align}

As the power prosumption is bounded, say between $P_\text{min}$ and $P_\text{max}$, histograms in (<ref>) are defined over a finite domain. Specifically, the domain is as:

\begin{align}
x \in \mathbb{X} = \{P_\text{min}, P_\text{min}+\Delta P, P_\text{min}+2\Delta P, \dots,  P_\text{max}\}. \label{eq:xdomain}
\end{align}

The value of $\Delta P$ is chosen as a trade-off between accuracy and computational efficiency. Indeed, the smaller the $\Delta P$ is, the more accurate the prediction of the PI is. This aspect is made clear below. However, choosing a very small step will require more memory to store all the measurements, thus resulting in a slower computation of the PI.

§.§.§ On-Line PI Estimation

At time $i$, the objective of the PI estimator is to determine the PI for the time slot $i+1$ at a given arbitrary target confidence level, said $\alpha$. The underlying idea is to assume the clusters (<ref>) as a statistically representative pool of possible realizations of the one-step-ahead power-prosumption realization. Therefore, the normalized histograms (<ref>) are assumed to be discrete PDFs and used to extract the symmetric quantiles corresponding to the $\alpha$ confidence level. Let

\begin{align}
l_i = c(P_i, t_i) \label{eq:currentlabel}
\end{align}

be the label calculated with the information at the current time instant. The PI lower and upper bounds are determined as

\begin{align}
{P}^{\downarrow\alpha}_{i+1|i} = (1-\alpha)/2 \text{~~quantile~of~~} \mathcal{G}_i^{l_i},\label{eq:upperquantile}\\
{P}^{\uparrow\alpha}_{i+1|i} = (1+\alpha)/2 \text{~~quantile~of~~}  \mathcal{G}_i^{l_i} \label{eq:lowerquantile}.
\end{align}

For the sake of clarity, the lower and upper quantiles in the expressions above are approximated by, respectively,  (also see fig:quantilecomp):

\begin{align}
\inf_{x \in \mathbb{X}} \left\{x : F_{}(x) \geq (1-\alpha)/2 \right\}\label{eq:lowerquantilecomp}
\end{align}

\begin{align}
\sup_{x \in \mathbb{X}} \left\{x : F_{}(x) \leq (1-\alpha)/2 \right\},\label{eq:upperquantilecomp}
\end{align}

where $F$ denotes the discrete CDF (cumulative distribution function) of $\mathcal{G}_i^l$ calculated by computing its cumulative sum.

Exemplification of the procedure (<ref>)-(<ref>) to approximate the quantiles for (<ref>)-(<ref>). In this case, the target confidence level $\alpha$ is 80%.

There are two advantages to this approach. First, histograms carry the complete information over the empirical PDFs, thus allowing for computing PIs at arbitrary confidence levels by training only one model. Second, it overcomes the problem of quantile crossing that arises, for example, in <cit.>, from treating the bounds at a given confidence level as two different time series.

§.§.§ On-Line Training

as time passes, new measurements become available and can be included to improve future PI estimates. Once the outcome $P_{i+1}$ is known, the normalized histogram associated with the label $l_i$ is updated with the new information, and the other histograms stay the same. Formally, the training procedure is as follows:

\begin{align}
\mathcal{G}_{i+1}^{l}(x) =
\begin{cases}
\hfil \phi \mathcal{G}_i^{l}(x) + (1-\phi) \delta(x-P_{i+1}) & l = l_i \\
\hfil \mathcal{G}_i^{l}(x) & l \neq l_i ,
\end{cases}
%\dfrac{a \mathcal{H}_i^{l_i}(x) + \delta(x-P_{i+1})}{a+1}
\end{align}

\begin{align}
\label{eq:forgetting-factor}
\phi=\frac{T_\phi/T}{T_\phi/T+1}
\end{align}
is called the forgetting factor, $T$ is the measurement period and $T_\phi$ is called the forgetting time constant. The forgetting factor controls how much past measurements influence the computation of PIs. Specifically, each new measurement has the same weight in the computation as all the measurements in the past $T_\phi$ seconds. The adoption of such a forgetting factor is important in order to track changes in the composition of prosumers' load/generation patterns.

§.§ PI Model B

§.§.§ Off-Line Batch Training

we apply the same principles described for Model A but on the once differentiated power-prosumption training data-set. The differentiated time series is denoted as:

\begin{align}
\label{eq:power-differentials}
B_{j} = P_j - P_{j-1},\;\; j=i-n+1,\dots,i.
\end{align}

The observations clusters are now calculated as follows:

\begin{align}
\begin{aligned}
\mathbb{H}^l_i = \{ B_{j+1}:c(P_{j}, t_{j}) = l, \; j=i-n+1,\dots,i\},\\
\end{aligned}
\end{align}

and are used to determine the normalized histograms of the differentiated power prosumption time series, which are denoted as $\mathcal{H}^0_i,\dots,\mathcal{H}^{L-1}_i$.

§.§.§ On-Line PI Estimation

The one-step-ahead PI bounds are computed as:

\begin{align}
{P}^{\downarrow\alpha}_{i+1|} = P_i + (1-\alpha)/2 \text{ quantile of } \mathcal{H}_i^{l_i},\\
{P}^{\uparrow\alpha}_{i+1|i} = P_i + (1+\alpha)/2 \text{ quantile of } \mathcal{H}_i^{l_i},
\end{align}

i.e., the current power-prosumption plus two back-off terms representing the expected power prosumption variation with respect to the current realization $P_i$.

§.§.§ On-Line Training

Once the prosumption $P_{i+1}$ is known and the power difference $B_{i+1}$ is computed, the normalized histogram corresponding to the current label $l_i$ is updated by adding the new differenced value, whereas the others stay the same,

\begin{align}
\mathcal{H}_{i+1}^{l}(x) =
\begin{cases}
\hfil \phi \mathcal{H}_i^{l}(x) + (1-\phi) \delta(x-B_{i+1}) & l = l_i \\
\hfil \mathcal{H}_i^{l}(x) & l \neq l_i.
\end{cases}
\end{align}

§.§ Classification According to Influential Variables

The assignment performed through the function $c(\cdot)$ in (<ref>) is realized by first clustering the historical measurements into groups based on similarities with respect to chosen influential variables and, then, assigning a cluster label to each new measurement. Needless to say, influential (or explanatory) variables are quantities that have an influence on the power prosumption. In general, they can be discovered using numerical methods (like analysis of variance, correlation analysis, or other procedures <cit.>) or identified by exploiting any empirical knowledge on the observed process (in our case the structure of prosumption of a given node).

Visualization of k-means clustering for a node with an EV charging station with 8 clusters and using the power level and the time of day as influential variables. Positive power is consumption.

In this work, influential variables are chosen by adopting the latter approach. The chosen variables are (i) power magnitudes and (ii) time of day (in seconds). The first variable accounts that, in a limited capacity feeder, prosumption variations depend on the same power-prosumption magnitude. For example, when the consumption is large, a load disconnection is more likely than a load insertion because many loads are already active, and vice-versa. The second variable is supposed to capture the different power-prosumption patterns that might occur during the day.

The historical measurements are grouped, by using a k-means approach, into $L$ clusters, according to the values of the influential variables. An example of the clustering method is shown on fig:k-means-clustering for a node with a level-3 EV charging station. The k-means algorithm was trained with two weeks of historical measurements by using both the power level and the time of day as influential variables and by grouping the measurements into $L=8$ clusters, each of which is labeled with a different color in the figure. The value of $L$ is chosen by the user and is fixed a priori. Its influence on the performance of the algorithm is studied in Section <ref>.

During the off-line and on-line training phases of the algorithm, the label of measurement $(P_i, t_i)$ is computed as follows:

\begin{equation}
l_i = \argmin {l=0,...,L-1} d(e(P_i, t_i), c_l)
\label{eq:label-computation}
\end{equation}

where $e(P_i, t_i)$ is the point in the $m-d$ space whose coordinates are the  $m$ influential variables chosen in the classification scheme,  $c_l$ is the center of cluster $l$ and $d(\cdot, \cdot)$ is the Euclidean distance between two points.

The different classification schemes are introduced with the objective of performing an a-posteriori validation of the selection process of influential variables and the number of clusters. Indeed, as stated earlier in this section, influential variables are assigned by exploiting the empirical knowledge of the process: By comparing the performance of different classification schemes (in sec:results), it is possible to infer whether the progressively more complex classification schemes are meaningful or not.

§.§ Implementation Aspects and Complexity

The main design requirement of the proposed algorithm is to deliver PIs in RT in order to, for example, activate enough capacity in inertia-less microgrids or to assist in the decision process of setting the droop controller of slack generators in the islanding maneuver <cit.>.
Given the large PIs reporting rate, computational complexity is a central aspect and hence is addressed in this section. The batch-training phase does not have RT requirements because it is performed off-line. This consists in labeling each observation of the training data-set by applying the discussed classification algorithm. As the number of labels $L$ is fixed by design, the overall complexity of the classification procedure, given by eq:label-computation, for one observation is constant time, or $O(1)$. Iterating it over a set of $N$ training data is an operation with linear time complexity, or $O(N)$.

The computation of the PIs and the on-line training are performed in rolling RT.
The former operation requires computing the label in (<ref>), which is $O(1)$, and the PI bounds by (<ref>)-(<ref>), which involve a minimum and maximum search over the set $\mathbb{X}$, a problem with log time complexity with respect to the set cardinality, $O(\log_2|\mathbb{X}|)$, e.g., using a binary search that can exploit the monotonicity of the discrete CDFs. However, as the cardinality of $\mathbb{X}$ is fixed by design in (<ref>), the complexity of the problem can be regarded as constant time.
Therefore, delivering PIs and performing on-line training are procedures whose complexities do not scale with the size of the problem.

The training data and progressively incoming measurements are encoded in $L$ normalized histograms. Each of them is stored using $2 \times |\mathbb{X}|$ doubles, specifically the height and value of each bin. For example, assuming a discretization of 1024 levels (10 bit), preserving the information for a 1 year at 20 ms of resolution with the proposed method requires 128 kB per label (considering a double representation of 64-bit), while storing the individual values would require approximately 14 Gb.

§ PERFORMANCE EVALUATION

§.§ A real case application: university buildings

To test the performance of the proposed models, we consider three sequences of power-prosumption measurements that were recorded from different points inside the MV distribution network of the EPFL campus.

The first one consists of an office building with a maximum consumption of 80 kW and that is equipped with a 30 kVA roof PV; the second one includes an office load with a maximum consumption of 30 kW and a 150 kW level 3 EV charging station; and the third one is a heat pump with a maximum consumption of 1.5 MW.

The measurements are with a resolution of 20 ms and are provided by a PMU-based metering infrastructure that has been deployed on the university campus (see <cit.>). We consider 45 days of historical power-prosumption measurements that span the period of September-October 2022. In each case, we consider two weeks of training data; they are used to construct the clusters (as explained in sec:classification) and to perform the off-line training of the algorithm. Then, the proposed PI estimation models are operated for a month (with on-line training), and the estimated PIs are validated against the latter data-set, at 20 ms resolution. Each month consists of approximately 130 million data points. The evaluation is performed in a simulated environment coded in C++. The simulations are executed in a Windows Server with 128GB RAM and an Intel Xeon Gold 6130 CPU at 2.10GHz.

§.§ Performance Metrics

We introduce the following metrics to allow for a quantitative comparison between the performance of models and classification schemes. The first is the PI normalized averaged width (PINAW), which is as follows:

\begin{align}
\text{PINAW} = \dfrac{1}{N} \sum_{j=1}^{N} ({P}^{\uparrow\alpha}_j - {P}^{\downarrow\alpha}_j)/P_\text{nom}.
\end{align}

The second metric is the PI coverage probability (PICP),  the percentage of power-prosumption realization that falls inside the predicted PI. It is as follows:

\begin{align}
\text{PICP} = \dfrac{1}{N} \sum_{j=1}^{N} b^\alpha_j
\end{align}

\begin{align}
b^\alpha_j = \begin{cases}
1,  & {P}^{\downarrow\alpha}_j \le P_j \le {P}^{\uparrow\alpha}_j \\
0, & \text{otherwise}.
\end{cases}
\end{align}

Because there is a trade-off between the width of the PI and the accuracy of the model, it is imperative that we define a third metric to quantify it. The metric we chose for this purpose is a modification of the coverage width-based criterion (CWC) proposed in <cit.> and is defined as follows:

\begin{align}
\label{eq:cwc}
\text{CWC} = \text{PINAW} \max(1, e^{-\mu \frac{\text{PICP}-\alpha}{1-\alpha}})
\end{align}

where $\mu$ is a user-defined parameter that quantifies the trade-off between PICP and PINAW. For our experiments, we chose $\mu=\frac{log(10)}{10}$. This means that a deviation of one order of magnitude in the error rate penalizes the width of the interval by $10$ times. The same result can be achieved with a linear penalty with slope $0.9$, as can be seen in fig:cwc_plot. The $x$-axis is the percentage deviation of the error rate from the target error rate, whereas the $y$-axis is the ratio $\frac{\text{CWC}}{\text{PINAW}}$. We observe that the exponential function penalizes less severely the error rates within one order of magnitude in comparison to the linear function. However, with our choice, higher error rates are more penalized.

[xmin=-2, xmax=17.5, ymin=0, ymax=17,
axis lines=middle,
x label style=at=(current axis.right of origin),anchor=north,
y label style=at=(current axis.above origin),anchor=south,
[red,domain=0:12]  exp(ln(10)/10*x) node[above]$\mu=\frac{log(10)}{10}$;
[red,domain=-3:0] 1;
[blue,domain=0:14.5] 1+0.9*x node[right]slope $0.9$;
[dashed] (0, 10) – (10, 10);
[dashed] (10, 0) – (10, 10);

Visual representation of the CWC. $\text{error}=1-\text{PICP}$ is the error rate of the algorithm, and $\text{target}=1-\alpha$ is the target error rate. The red line indicates eq:cwc with the chosen value of $\mu$, and the blue line indicates a linear penalty.

§.§ Clustering Using Power Levels

To evaluate the algorithm, we consider four target confidence levels $\alpha$, namely 99, 99.9, 99.99, and 99.999%. The parameters to be analyzed are (i) the model A or B, (ii) the number of clusters $L$ and (iii) the forgetting factor $\phi$ or, equivalently, the forgetting constant $T_\phi$. Each combination of model, $L$ and $T_\phi$ is called a configuration. To choose the best configuration for each building, we run the algorithm for both models A and B, with all possible combinations of values of $L$ in the set $\{1, 8, 64, 256, 512, 1024\}$ and $T_\phi$ in the set $\{1, 60, 3600, 21600, 86400, 604800\}$ seconds[The values correspond to 1 second, 1 minute, 1 hour, 6 hours, 1 day and 1 week respectively.]. For all the experiments, the value of the quantization step $\Delta P$ was chosen such that the size of the domain $\mathbb{X}$ (see eq:xdomain) is equal to 2000 points. It was experimentally observed that the values of PIs computed using larger domains did not change within three significant digits. Also, the average time needed to compute the PI and do one cycle of on-line training was around $50\mu s$; thus, the method meets the RT requirements.

For each combination of models, of a number of clusters, and time constants, and for each confidence level $\alpha$, we compute the three metrics PINAW, PICP, and CWC. Tables <ref>, <ref>, and <ref> show the top $5$ configurations, i.e., those that achieve the smallest CWC, for each target confidence level, for the three buildings, respectively.

Concerning the office building and the one hosting the EV charging station, model B outperforms model A for every confidence level. Also, for confidence levels up to $99.9\%$, the algorithm is able to achieve the target confidence. For higher confidence levels, however, the algorithm misses the target confidence by at most one order of magnitude. Given our choice of the parameter $\mu$, this means that the computed CI is penalized by a maximum factor of $10$ in the computation of CWC.

Regarding the optimal configuration, fewer clusters, together with larger forgetting time constants, are better choices. More specifically, the optimal number of clusters $L$ does not need to be larger than 8, whereas a value of $T_\phi$ between one day and one week is the optimal choice. By taking a closer look at Tables <ref> and <ref>, we observe a trade-off between the values of $L$ and $T_\phi$ in the performance of the algorithm. If we focus on table:elg_cwc_99, for example, we notice that we could choose, without affecting the value of CWC by more than $1\%$, either a combination of $8$ clusters and a forgetting time constant of one week, or $256$ clusters and one day, or $1024$ clusters and $6$ hours. In fact, this table showcases that the choice of $L$ and $T_\phi$ has a minor effect on the computation of PIs for the office building. The performance of the algorithm is mainly influenced by the differentiation of the measurements performed by model B. For the charging station, however, more clusters perform up to $30\%$ worse compared to the optimal case of only $8$ clusters.

The results are quite different for the heat pump in table:cct_cwc. First of all, the algorithm can only approach the target confidence levels within one order of magnitude. Unlike the other two buildings, the optimal results are achieved by a combination of few clusters (ideally only one) and a small forgetting time constant (less than an hour). This implies that the consumption of the heat pump changes more rapidly than that of the office building and the charging station, and that the measurements older than one hour do not influence the computation of PIs. We also observe that model A outperforms model B for large confidence levels. Even though the value of PINAW computed by model A is up to $5$ times larger than the one computed by model B, model B fails to achieve the target confidence level by up to two orders of magnitude, which results in an exponential increase in the metric CWC.

Overall, we conclude that the optimal choice of model, $L$ and $T_\phi$ depends on the characteristics of the prosumer. Nodes with low volatility, such as an office building, benefit from a long memory of up to one week, whereas nodes characterized by rapid power changes require the use of short memory. The clustering of power measurements seems to benefit mainly nodes with clearly distinct power-levels, such as those at a charging station. But in any case, the number of clusters does not need to be more than $8$. Finally, the differentiation of power (i.e., model B) improves the computation of PIs, provided that the power differentials have low volatility, as is the case for the office building and the charging station.

§.§ Clustering Using Power Levels and Time of Day

In this section, we consider whether adding the time of day (TOD) as a feature of the algorithm, in addition to the power level (P), would improve the performance of the algorithm. The idea is that the clustering based on TOD will capture the patterns that the prosumption exhibits during the day and that this might affect the computation of the PIs. We re-run the experiments of the previous section, with the addition of the TOD. The best configuration is again the one that achieves the smallest value of CWC for each target confidence level. The comparison between the two cases, specifically (i) using only power level (P) and (ii) using power level and the time of day (P+TOD), is shown on table:time-of-day.

We observe that the best configuration for the two clustering strategies is the same in most cases. For the office building, the introduction of the TOD does not affect the value of CWC by more than $1\%$, which confirms once again that the clustering method does not affect the performance of the algorithm. For the charging station, the value of CWC differs by $5-10\%$ between the two clustering methods. However, the introduction of the new feature could either improve or worsen the performance, depending on the confidence level. Therefore, we cannot conclude whether it consistently performs better. As far as the heat pump is concerned, the usage of the TOD in the clustering has no effect on the performance, because the algorithm performs better when all measurements are put in one cluster.

The results indicate that this feature does not play a critical role in the computation of PIs. Perhaps other features, apart from the power level and the time of day, could influence the performance of the algorithm. The effect of additional features could be studied in future research.

§.§ Effect of the Measurement Period

[Office building]

[Charging station]

[Heat pump]

Performance evaluation as a function of the measurement period

A crucial objective of this work is to find out how the proposed algorithm scales as the measurement period increases. To test the algorithm on different measurement periods, we integrate the available data-set. In particular, given the measurements $P_j^{20}, j=1..N$ at $20ms$ resolution, the measurements at resolution $T~(ms)$ are recomputed as:

\begin{equation}
P_i^{T} = \frac{20}{T} \sum_{j=(i-1)\frac{T}{20}+1}^{i\frac{T}{20}} P_j^{20}, i=1..\frac{20}{T}N
\end{equation}

where $T$ is assumed to be a multiple of $20ms$. The forgetting factor $\phi$ is scaled according to $T$, as in eq:forgetting-factor. For this and the following section, we do clustering based only on the power level.

We perform again simulations with varying cluster numbers $L$ and forgetting time constants $T_\phi$ for different confidence levels. In fig:resolution-performance, we plot the values of PINAW and PICP as a function of the measurement period from $20ms$ up to $300s$ for the three buildings considered. Each point on the graphs corresponds to a different configuration (the one that achieves the smallest CWC), as shown in table:config-resolution.

The graph for the office building showcases that the algorithm performs well even for larger measurement periods, provided that the target confidence level is less than $99.99\%$. Indeed, for low confidence levels, the error rate is kept below the target, and the average width of the PI is less than $20\%$ of the nominal value.

For larger confidence levels, however, the algorithm fails because either the target level is not achieved or the width of the PI is so high that it becomes useless for grid control. It is worth noting that at a resolution of five minutes, the result of the algorithm is identical for all target confidence levels above $99.99\%$. This implies that, when the measurements are sparse, we might need a larger history to meet high confidence levels.

Similar observations can be made for the charging station. The main difference is that the maximum measurement period with an acceptable performance depends on the target confidence level. With target confidence of $99\%$, the algorithm computes a low PI width for measurement periods up to $10sec$; whereas, with $99.9-99.999\%$, it has acceptable performance but only when the resolution is sub-second. For higher confidence, the algorithm cannot predict accurate PIs for periods larger than $20ms$.

Looking at the results for the heat pump on fig:resolution-performance, we notice a break in the trend of increasing PINAW. To understand why this happens, we look at table:resolution-cct. The best configuration for confidence levels $99.99-99.999\%$ changes from using model A at $20ms$ resolution to model B at higher resolutions. Model A has been shown in sec:clustering-results to generate, in general, larger PIs than model B. This increase in the width results in a lower error rate, which in turn might reduce the value of CWC, which is the sole metric used to compare configurations to one another. The value of CWC ultimately depends on the choice of the parameter $\mu$, which might affect the trends in the graphs of fig:resolution-performance.

§.§ Confidence Level Uncertainty

[Office building]

[Charging station]

[Heat pump]

Box plot of model accuracy versus target confidence level at $20ms$ resolution

The metrics of sec:metrics evaluate, over the full one-month period, the average performance of the algorithm. However, it would be interesting to see how the algorithm performs over time. For this purpose, we split the one-month evaluation period into six-hour windows and compute the error rate achieved by the algorithm in each window. Hence, for each confidence level, we compute a histogram of $120$ estimations of the error rate. We then depict the statistical measures of the histograms using the box plots of fig:box_plots. The blue box indicates the 25th and 75th percentiles, the red line is the median, the red crosses are outliers, and the small square indicates the target error rate. We computed the box plots only for the best configuration for each building and for confidence level at $20ms$ resolution, which are shown in table:config-resolution.

If the algorithm is consistent in predicting accurate PIs, then the median of the box plot should be close to the target error-rate. We see that this is indeed the case for the office building. For the largest confidence level, our algorithm performs even better than expected, albeit with many outliers. For the charging station, the performance is sometimes worse than expected, but the target error-rate is nevertheless contained within the 25th and the 75th percentile of the box plot. This is not always the case for the heat pump, hence the algorithm cannot consistently estimate accurate PIs if there is high volatility in the measurements.

The results showcase that our algorithm computes accurate
PIs, irrespective of the building, provided that the target confidence level is less than $99.9\%$. For higher confidence levels, the performance is still acceptable, but the consistency of the predictions depends on the dynamics of the node.

§ CONCLUSION

Motivated by the requirements of ADNs real-time control, we have presented a non-parametric method for computing ultra-short-term PIs (prediction intervals) of the power prosumption in generic ADNs (buildings).
The method consists in grouping historical measurements into clusters, according to the value of selected influential variables. The clusters are considered statistically representative pools of future power-prosumption realizations and are used to extract PIs at arbitrary confidence levels by calculating the quantiles from the respective PDF.

The proposed method has been applied to the original and once-differentiated power-prosumption time series, and different influential variables have been considered. The performance of the method was tested for different types of prosumers by using experimental measurements from an MV distribution network.
The performance analysis enabled us to make an a-posteriori selection of the parameters of the algorithm. The algorithm was shown to compute relatively narrow PIs for the studied prosumers, for time resolutions from $20ms$ up to a few minutes in some cases, provided that the target confidence level is below $99.9\%$.

A final statement concerns the computational complexity, which becomes a relevant concern especially when considering densely sampled time series and the high reporting rate for the predictions.
We have shown that the proposed algorithm performs the PI computation and on-line training in constant time hence is scalable.