\section{Introduction}\label{sec:introduction}
\input{sections/introduction}

\section{Methodology}\label{sec:methodology}
\input{sections/methodology}

\section{Case study}\label{sec:case_study}
\input{sections/case_study}

\section{Results}\label{sec:results}
\input{sections/results}

\section{Conclusion and outlook}\label{sec:conclusion}
\input{sections/conclusion}

\section*{Acknowledgement}
The work of H. Cai and P. Heer was supported by the Swiss Federal Office of Energy through the Sustainable Demand-Side Management for the Operation of Buildings (S-DSM) project under the contract number SI/502165-01. The work of T. Zufferey and G. Hug is part of the activities of the Renewable Management and Real-Time Control Platform (ReMaP).

\bibliographystyle{IEEEtran}




\subsection{Original Flexibility Envelope Quantification}\label{sec:original_flex_quant}

Initially, the upper and lower energy bounds are identified by energizing the building to its extremes, which consists of solving two optimization problems as follows \cite{caha2021}:
\begin{subequations}\label{eq:opt_based_quant}
\begin{alignat}{1}
    J_{\mathrm{upper}} &= \min_{u,\delta} \; - w^\intercal u + w_\delta \cdot \lVert \delta \rVert^2,\\
    J_{\mathrm{lower}} &= \min_{u,\delta} \; w^\intercal u + w_\delta \cdot \lVert \delta \rVert^2,
\end{alignat}
\end{subequations}
where $u$ is a vector of active power of all components, $w$ is a vector of weights, $\delta$ is a vector of slack variables introduced to ensure feasibility and $w_\delta\in\mathbb{R}_+$ is the corresponding weight. Optimization constraints include the dynamics of each component, power and energy capacities, state limits and mutual exclusiveness of charging and discharging. The respectively obtained $u$ gives the upper and lower energy bounds. 
The feasible power levels and the corresponding maximum sustained duration can be derived using the procedure shown in Fig.~\ref{fig:envelope_calculation} \cite{caha2021}. When applying this derivation in all the lead times within the horizon, a flexibility envelope spanned by lead time, power levels, and maximum sustained duration is obtained.  An example of flexibility envelope with a 24-hour horizon is shown in Fig.~\ref{fig:envelope}.
In real-time applications, the entire quantification process is periodically repeated to account for the most recent observations.
In the rest of this paper, the prediction and approximation of flexibility envelopes concretely refer to the prediction and approximation of the maximum sustained duration values, which depend on the two other dimensions.

\begin{figure}
    \centering
    \includegraphics[width=\columnwidth]{./figures/flex_concept}
    \caption{Calculation of maximal sustainable duration by (a) first identifying energy bounds and then (b) deriving sustainable duration of chosen power levels.}
    \label{fig:envelope_calculation}
\end{figure}

\subsection{Prediction of Flexibility Envelope}

Each flexibility envelope can be predicted by~\ac{ML} models in order to dispense with the computationally intensive solution of the optimization problem described in (\ref{eq:opt_based_quant}). 
In this work, a simple benchmark model is compared to a linear regressor, a~\ac{KNN} algorithm, a~\ac{SVM}, and an AdaBoost regressor from the scikit-learn Python library. They are presented in detail in the following subsections. Hyperparameter tuning is performed using grid search with 5-fold cross-validation on the training set.
Predictions are performed based on the following input features:
\begin{itemize}
    \item calendar data: time of the day (cyclic encoding)
    \item weather forecast from MeteoSwiss: solar irradiance (in~$\text{W/m}^2$) and temperature (in~$\degree \text{C}$)
    \item residential behavior from historical data: presence at home (binary), non-flexible load (in~$\textrm{kW}$), water draw (in~$\textrm{L/min}$)
    \item comfort requirements: room temperature limits (in~$\degree \text{C}$), domestic hot water temperature limits (in~$\degree \text{C}$)
    \item predicted energy state of flexible appliances: home battery energy (in~$\textrm{kWh}$),~\ac{EV} battery energy (in~$\textrm{kWh}$), \ac{PV} energy production (in~$\textrm{kWh}$), room temperature  (in~$\degree \text{C}$), domestic hot water temperature  (in~$\degree \text{C}$)
\end{itemize}
Note that the input features are a subset of the information required by the base model predictive controller and no additional data are necessary. Furthermore, the linear regressor, KNN, and SVM algorithms require features scaling, which is done by standardization in this work:
\begin{equation}\label{eq:scaling}
    x' = \frac{x - \mu}{\sigma},
\end{equation}
where $x$ and $x'$ are the original and scaled feature values, respectively. In addition, $\mu$ and $\sigma$ are the mean and standard deviation of the original feature values, respectively.

\subsubsection{Linear Regressor}

The linear model outputs a linear combination of the input features:
\begin{equation} \label{eq:lin}
	 \hat{y} = \sum_{i=1}^m \beta_i\cdot x_i
\end{equation}
where $\hat{y}$ is the predicted value, $x_i$ is the corresponding $i^{\textrm{th}}$ input feature, $\beta_i$ is the weight associated with the $i^{\textrm{th}}$ input feature, and $m$ is the total number of features. In this study, the weights are tuned based on the Lasso regression from the Lasso Python class. The training process consists of minimizing a least square loss function with L1 regularization:
\begin{equation}\label{eq:lin_loss}
    J_{\mathrm{lin}} = \min_{\beta} \; \lVert y - X\beta \rVert^2 + \gamma \cdot \lVert \beta \rVert
\end{equation}
where $y$ is a vector of length $n$ containing the target values in the training set, $X$ is an $n \times m$-matrix containing the input features, $\beta$ is a vector of length $m$ containing the weights of the linear regression, and $n$ is the total number of examples in the training set. The regularization parameter $\gamma$ is a hyperparameter whose value is selected from $\{0.2, 0.4, 0.6, 0.8, 1\}$. Both Lasso and Ridge (i.e., L2 regularization) regressions handle collinearity well \cite{dormann2013collinearity}. Nevertheless, the Lasso regression has been chosen over the Ridge regression since it performs better on the training set.

\subsubsection{K-Nearest Neighbors}

The~\ac{KNN} is a non-parametric model whose prediction is a linear combination of the target values of the $k$ nearest neighbors in the training set:
\begin{equation} \label{eq:knn}
	 \hat{y} = \sum_{i=1}^k w_i^{\mathrm{knn}}\cdot y^{\mathrm{knn}}_i,
\end{equation}
where $\hat{y}$ is the predicted value, and $y^{\mathrm{knn}}_i$ and $w_i^{\mathrm{knn}}$ are the target value and the weight associated to the nearest neighbor $i$, respectively. In this work, the model comes from the KNeighborsRegressor Python class, where the number of nearest neighbors $k$ is a hyperparameter whose value is selected from $\{20, 40, 60, 80, 100, 120, 140, 160\}$. Standard Euclidean distance is used as the similarity metric.

\subsubsection{Support Vector Machine}

Used as a regressor, the~\ac{SVM} is also a non-parametric approach that consists of two steps.  In the first step, the original input features are mapped into a higher dimensional space induced by a so-called kernel function. In this new space, the computation occurs in a more efficient manner. In the second step, the algorithm determines the hyper-plane with the narrowest possible margin with a so-called $\varepsilon$-insensitive loss through the input data. Formally, this consists of solving the following convex optimization problem:
\begin{subequations}
\begin{alignat}{3}
&\min_{w,\xi_i,\xi_i^*} \quad && \dfrac{1}{2}\cdot \lVert w \rVert^2 + C \sum_{i=1}^{n} (\xi_i +  \xi_i^*), \\
&\,\,\,\mathrm{s.t.}  && w^\intercal x_i - y_i \geq \varepsilon + \xi_i,\\
&               && y_i-w^\intercal x_i  \geq \varepsilon + \xi_i^*,\\
&               && \xi_i,\xi_i^*,\varepsilon \geq 0, \quad \forall i\in \{1,...,n\},
\end{alignat}
\end{subequations}
where $x_i$ and $y_i$ are a vector of input features (in the higher dimensional space) and the target value corresponding to the $i^{\textrm{th}}$ training sample, respectively. In addition, $w$ is a vector of weights, $C$ is a regularization parameter, $\varepsilon$ defines the margin of insensitive loss, and $n$ is the total number of samples in the training set. Next, $\xi_i$ and $\xi_i^*$ are slack variables introduced to ensure feasibility. The goal is to find a function that stays within an $\varepsilon$-margin from the predicted values $y$, but also exhibits flatness. The first requirement is met through minimization of the $\varepsilon$-insensitive loss. The second requirement is reflected in the L2-regularization of the weights. There is an extensive literature on the SVR algorithm and more information can be found in~\cite{drucker1997support}. In this work, the model comes from the SVR Python class, where $\varepsilon$ is set to $0.1$ and the radial basis function is chosen as Kernel function. Finally, $C$ is a hyperparameter whose value is selected from $\{0.1, 0.2, 0.5, 1, 2, 5\}$.

\subsubsection{AdaBoost Regressor}

The AdaBoost approach is an ensemble method (i.e., multiple weak learners are combined to create a strong learner), where the training of the $(n+1)^\textrm{th}$ weak learner depends on the results of the $n^\textrm{th}$ weak learner. Hence, each weak learner focuses on the error of the previous weak learner. The final prediction is a weighted average of the prediction of all weak learners, where the respective weights are calculated based on the error of each weak learner. More information on the AdaBoost regressor is available in~\cite{freund1997decision}. In this work, the model comes from the AdaBoostRegressor Python class, and the weak learners consist of decision tree regressors. The number of weak learners and the maximal depth of decision trees are hyperparameters whose values are selected from $\{20, 50, 100\}$ and $\{3, 5, \textrm{none}\}$, respectively.

\subsubsection{Benchmark Model}

The \ac{ML}-based predictors are also compared to a simple benchmark model which exploits the time-dependency and the daily pattern of flexibility. Concretely, the benchmark model divides the day into intervals of the same length. For each interval, the average sustained duration is computed for each power level based on the data in the training set. These average values are used as predictions for the sustained duration in the test set:
\begin{equation} \label{eq:benchmark}
	 \hat{y}_j = \frac{1}{|\Omega|}\sum_{i\in\Omega} y_i, \quad \forall j\in\Omega,
\end{equation}
where $\hat{y}_j$ and $y_i$ are the predicted and historical values whose characteristics comply with the condition of $\Omega$, i.e. correspond to a given period of the day. The length of the intervals is a hyperparameter selected from $\{0.25\textrm{h}, 0.5\textrm{h}, 1\textrm{h}, 2\textrm{h}, 4\textrm{h}, 8\textrm{h}\}$ and varies among the various power levels.




\subsection{Approximation of Flexibility Envelope}

The approximation of flexibility envelopes via~\acp{PDF} allows for a reduction of the number of parameters to be transmitted to \acp{DSO}, however at the cost of decreased accuracy. In this work, the flexibility envelopes are approximated based on the 2-dimensional (2D)~\ac{ND}, the 2-dimensional~\ac{SND}, and the 3-dimensional~\ac{GMM}. The 2D models fit a \ac{PDF} to a slice of the flexibility envelope, i.e. at a given lead time within the horizon (see Fig.~\ref{fig:envelope_calculation}(b)), whereas the 3D model directly approximates the entire envelope. Here, the horizon refers to the entire time window (i.e., set of all considered lead times) of the flexibility envelope.

\subsubsection{2D Normal Distribution}

A flexibility envelope at a given lead time tends to follow a normal (or Gaussian) distribution:
\begin{equation}\label{eq:nd}
    \phi(x)  = \frac{1}{\sigma\sqrt{2\pi}}\cdot e^{-\frac{(x-\mu)^2}{2\sigma^2}},
\end{equation}
where $\sigma$ and $\mu$ are the standard deviation and the mean of the~\ac{PDF}.

In this work, the approximation model consists of a weighted sum of two~\acp{ND}, which shows improved performance over a single~\ac{ND}. Formally, the model is defined as:
\begin{equation}\label{eq:2Dsum}
    f(x)  = \sum_{i=1}^2 c_i \cdot \phi_i(x),
\end{equation}
where $\phi_i(x)$ is a~\ac{ND} as defined in~(\ref{eq:nd}), and $c_i$ is the corresponding weight. Parameters are optimized by the~\textit{curve.fit} function from the scipy.optimize Python library, which relies on the Levenberg-Marquardt algorithm~\cite{more1978levenberg}.

\subsubsection{2D Skewed Normal Distribution}

Flexibility envelopes appear to be mostly asymmetric. Indeed, buildings import electricity most of the time such that the flexibility potential is generally skewed to the positive side. Hence, the skewed normal distribution is also considered:
\begin{equation}\label{eq:snd}
    \phi^\textrm{skew}(x)  = \frac{2}{\omega} \cdot \phi\left(\frac{x - \xi}{\omega}\right) \cdot \Phi\left(\alpha \cdot \frac{x - \xi}{\omega}\right),
\end{equation}
where $\phi(x)$ is the~\ac{PDF} of the standard~\ac{ND} as defined in~(\ref{eq:nd}) and $\Phi(x)$ is the~\ac{CDF} of the~\ac{ND}. In addition, $\alpha$, $\xi$, and $\omega$ represent the skewness, the location, and the scale parameters, respectively.

Similar to the previous approximation model, this model considers the superposition of two~\acp{PDF} as defined in~(\ref{eq:2Dsum}), where $\phi_i$ is replaced by $\phi_i^\textrm{skew}$.

\subsubsection{3D Gaussian Mixture Model}

A~\ac{GMM} is the linear combination of a finite number of normal (or Gaussian) distribution components. Its 3D version can be expressed as follows:
\begin{equation}\label{eq:3D-GMM}
    \textrm{GMM}^{\textrm{3D}}(x)  = \sum_{i=1}^K c_i \cdot \phi^{\textrm{3D}}_i(x),
\end{equation}
where $\phi^{\textrm{3D}}_i(x)$ is the $i\textsuperscript{th}$ 3D Gaussian component, and $c_i$ is the respective weight. In contrast to its 2D version, each 3D Gaussian component is defined by 2 means and a $2\times 2$ covariance matrix. In this work, the covariance matrix is diagonal such that only 2 parameters are necessary. In addition, the number of components $K$ can be set such as $18$ to balance between approximation accuracy and computational cost. Equation~(\ref{eq:3D-GMM}) is an extension of the model presented in~(\ref{eq:2Dsum}). Furthermore, a preliminary analysis shows that the performance can be improved if the values at the first and the last lead time are duplicated prior to model fitting. Therefore, this empirical insight is incorporated in the subsequent results.

\subsection{Evaluation Metrics}

For a comprehensive assessment, it is crucial to evaluate a model based on several metrics which can account for different aspects of its performance. Detailed in the following subsections, the error metrics considered in this work are the~\ac{MAE},~\ac{NRMSE}, and $\textrm{R}^2$ score. The number of parameters and the computation time are also considered as key performance indicators in the context of this study.

\subsubsection{Mean Absolute Error}
The~\ac{MAE} accounts for the absolute difference between the estimated and the real sustained duration:
\begin{equation}\label{eq:MAE}
    \text{MAE}  = \frac{\sum_{i=1}^{n} \mid \hat{y}^{i} - y^{i} \mid}{n},
\end{equation}
where $y^{i}$ and $\hat{y}^{i}$ are the real and estimated value of data point $i$, respectively. In addition, $n$ refers to the total number of data points under consideration. The lower the~\ac{MAE}, the better the performance.

\subsubsection{Normalized Root Mean Square Error}

The~\ac{NRMSE} is the normalized version of the~\ac{RMSE} which penalizes larger errors more:
\begin{equation}\label{eq:NRMSE}
    \text{NRMSE}  = \frac{1}{\bar{y}} \cdot \sqrt{\frac{\sum_{i=1}^{n} (\hat{y}^{i} - y^{i})^2}{n}},
\end{equation}
where $\bar{y}$ is the average value over all real data points under consideration. The lower the~\ac{NRMSE}, the better the performance. In contrast to the~\ac{MAE}, the~\ac{NRMSE} enables to fairly compare the outcome at different power levels even if the sustained duration values are not necessarily in the same order of magnitude.

\subsubsection{$\textrm{R}^2$ Score}

The $\textrm{R}^2$ score, also referred to as coefficient of determination, explains how much of the variance in the dataset can be estimated by the model. It is defined as:
\begin{equation}\label{eq:R2}
    \textrm{R}^2  = 1 - \frac{\sum_{i=1}^{n} ( \hat{y}^{i} - y^{i})^2}{\sum_{i=1}^{n} (y^{i} - \bar{y})^2}.
\end{equation}
The coefficient of determination is upper bounded at 1, but has no lower bound. The higher the score, the better the performance. A value of zero corresponds to the performance of a model that always outputs the mean value of the test set.

\subsection{Prediction of Flexibility Envelope}

\begin{table}
\caption{Performance summary of the prediction models}
\label{tab:prediction_summary}
\centering
 \begin{tabular}{|c|c|c|c|c|}
 \hline
 \begin{tabular}{@{}c@{}}Predict.\\model\end{tabular} &
 \begin{tabular}{@{}c@{}}MAE\\(min)\end{tabular} &
 \begin{tabular}{@{}c@{}}NRMSE\\(\%)\end{tabular} &
 $\textrm{R}^2$ score &
 \begin{tabular}{@{}c@{}}Comput.\\time (sec)\end{tabular}\\ [0.5ex] 
 \hline
 Original & - & - & - & $20.75\pm 16.05$  \\
 \hline
 Benchmark & 70 & 46 & 0.67 & 0.15 \\
 Linear & 84 & 61 & 0.49 & 0.19\\
 KNN & 42 & 35 & 0.78 & 0.59\\
 SVM & 42 & 30 & 0.81 & 1.48\\
 AdaBoost & 25 & 21 & 0.89 & 0.52\\
 \hline
\end{tabular}
\end{table}

The benchmark model and the four \ac{ML}-based models have been tuned on the training set and applied to the samples in the test set, which consists of 24 hours of data. Their performance is summarized in Table~\ref{tab:prediction_summary} in terms of~\ac{MAE} (in minute),~\ac{NRMSE} (in percentage), $\textrm{R}^2$ score, and average computation time (in second, per flexibility envelope). AdaBoost exhibits by far the best prediction accuracy according to all error metrics. In contrast, the linear model performs particularly badly, even worse than the benchmark model. In terms of computation time, AdaBoost is the second fastest \ac{ML} model after the linear model, which makes it the most appropriate candidate. In comparison, the original calculation using one CPU process is estimated to take $20.75\pm 16.05$ 
seconds per flexibility envelope most of the time. The reduction in computational efforts is achieved because \ac{ML}-based approaches do not involve an iterative process in searching for the solution. 

Moreover, Figs.~\ref{fig:prediction_MAE},~\ref{fig:prediction_NRMSE}, and~\ref{fig:prediction_R2} provide insights into the error distribution with respect to the power level. First of all, the ranking of performance observed in Table~\ref{tab:prediction_summary} among the various prediction models is respected for a majority of power levels. In terms of absolute error (i.e., MAE), power levels close to zero are naturally more impacted due to the generally higher sustained duration values. In the case of \ac{NRMSE}, however, the accuracy tends to decrease when moving to the sides of the flexibility envelope. A notable exception is the power levels between $-3$~kW and $1$~kW, where the relative error increases substantially. This is principally due to the normalization. The total flexibility in that power range is mainly driven by the PV flexibility which has strong diurnal patterns but a relatively low average over the entire day.
Finally, the $\textrm{R}^2$ score indicates that the variance is generally harder to predict between the power levels of $-2$ kW and $8$ kW. Indeed, such power levels are a composition of small loads and PV output, characterized by high variance. In contrast, the flexibility outside $[-2,8]~\textrm{kW}$ is mainly attributed to the~\ac{EV}, which has more regular patterns.

\begin{figure}[t!]
    \centering
    \includegraphics[width=\columnwidth]{./figures/prediction_MAE}
    \caption{MAE of the flexibility prediction with respect to the power level.}
    \label{fig:prediction_MAE}
\end{figure}

\begin{figure}[t!]
    \centering
    \includegraphics[width=\columnwidth]{./figures/prediction_NRMSE}
    \caption{NRMSE of the flexibility prediction with respect to the power level.}
    \label{fig:prediction_NRMSE}
\end{figure}

\begin{figure}[t!]
    \centering
    \includegraphics[width=\columnwidth]{./figures/prediction_R2}
    \caption{$\textrm{R}^2$ score of the flexibility prediction with respect to the power level.}
    \label{fig:prediction_R2}
\end{figure}

\subsection{Approximation of Flexibility Envelope}

The three approximation models have been applied to all the original flexibility envelopes in the dataset. Their performance is summarized in Table~\ref{tab:approximation_summary} in terms of number of parameters, average~\ac{MAE} (in minute), average~\ac{NRMSE} (in percentage), average $\textrm{R}^2$ score, and average computation time (in second) per flexibility envelope. In terms of accuracy, the 3-dimensional model is clearly outperformed by both 2-dimensional models. Indeed, each slice of envelope (i.e., corresponding to a given lead time) exhibits a distinctive Gaussian shape, but the specific parameters highly vary over the different lead times. Moreover, the skewed normal distribution achieves a slightly higher performance than the standard version. We can observe from Figs.~\ref{fig:approximation_MAE},~\ref{fig:approximation_NRMSE}, and~\ref{fig:approximation_R2} that the 3D model further leads to a larger spread of the approximation errors. It should also be noted that the 2D models, and especially the~\ac{SND}, can almost perfectly account for the variance in the data. In contrast, the 3D-GMM outperforms both 2D approximators in terms of computational time (e.g., $8.7$ times faster than the 2D-ND) and required number of parameters. In fact, the 2D-ND, 2D-SND, and 3D-GMM allow for a reduction of the number of parameters by a factor of $8.83$, $6.63$, and $12.96$, respectively, compared to the original representation.
Hence, the most appropriate algorithm depends on the preference between approximation accuracy, computational burden, and communication load.

\begin{table}[t!]
\caption{Performance summary of the approximation models.}
\label{tab:approximation_summary}
\centering
 \begin{tabular}{|c|c|c|c|c|c|} 
 \hline
 \begin{tabular}{@{}c@{}}Approx.\\model\end{tabular} &
 \begin{tabular}{@{}c@{}}Number of\\parameters\end{tabular} &
 \begin{tabular}{@{}c@{}}MAE\\(min)\end{tabular} &
 \begin{tabular}{@{}c@{}}NRMSE\\(\%)\end{tabular} &
 $\textrm{R}^2$ score &
 \begin{tabular}{@{}c@{}}Comput.\\time (sec)\end{tabular}\\ [0.5ex] 
 \hline
 2D-ND & 132 & 32 & 28 & 0.97 & 3.68\\
 2D-SND & 176 & 23 & 23 & 0.98 & 5.50\\ 
 3D-GMM & 90 & 66 & 72 & 0.83 & 0.63\\ 
 \hline
\end{tabular}
\end{table}

\begin{figure}[t!]
    \centering
    \includegraphics[width=\columnwidth]{./figures/approximation_MAE}
    \caption{Distribution of the MAE for the flexibility envelope approximation.}
    \label{fig:approximation_MAE}
\end{figure}

\begin{figure}[t!]
    \centering
    \includegraphics[width=\columnwidth]{./figures/approximation_NRMSE}
    \caption{Distribution of the NRMSE for the flexibility envelope approximation.}
    \label{fig:approximation_NRMSE}
\end{figure}

\begin{figure}[t!]
    \centering
    \includegraphics[width=\columnwidth]{./figures/approximation_R2}
    \caption{Distribution of the $\textrm{R}^2$ score for the flexibility envelope approximation.}
    \label{fig:approximation_R2}
\end{figure}

\subsection{Overall Comparison}

\begin{figure*}
    \centering
    \includegraphics[width=\textwidth]{./figures/example}
    \caption{Example of the prediction and approximation of the flexibility envelope displayed in Fig.~\ref{fig:envelope}.}
    \label{fig:example}
\end{figure*}


Fig.~\ref{fig:example} shows the outcome of the benchmark prediction model as well as the best performing prediction (i.e., AdaBoost) and approximation (i.e., 2D skewed normal distribution) models for the example flexibility envelope from Fig.~\ref{fig:envelope}. The upper subplots display the predicted and approximated flexibility envelopes, whereas the lower subplots show the difference between the respective upper subplots and the original envelope. While the benchmark prediction model cannot account for the variation among different lead times, both best performing models provide an envelope visually close to the original envelope.

In addition, the best performing prediction and approximation models are integrated into a single pipeline, which is displayed in the right-hand side subplots. Formally, the 2D-SND approximates the envelope predicted by the AdaBoost algorithm. This corresponds to the process desired in this work to address both the computational burden and communication load. 

Table~\ref{tab:pipeline_summary} summarizes the performance of the different models visualized in Fig.~\ref{fig:example}. As expected, the benchmark prediction model exhibits by far the worst performance. In absolute values, the entire envelope can be better predicted than approximated. However, the approximation still leads to lower relative error and higher $\textrm{R}^2$ score. More importantly, although the combination performs slightly worse than both models separately, the respective errors do not directly add up. This is a very encouraging outcome and shows that these various models can be used to considerably speed up the flexibility quantification process and lower the communication load without substantial loss of accuracy.

Finally, a hybrid prediction model has also been considered as an alternative to the overall best performing model. Precisely, instead of predicting all power levels based on the AdaBoost algorithm, model selection is carried out in the training phase such that the best performing prediction model per power level is applied on the test set (see Fig.~\ref{fig:prediction_MAE}). Consequently, AdaBoost is still used in a large majority of the cases but is replaced by the~\ac{SVM} and~\ac{KNN} algorithms for some power levels close to zero. This hybrid approach allows for a MAE of 31 minutes (i.e., $-4$ minutes, compared to the accuracy displayed in Table~\ref{tab:pipeline_summary}), a NRMSE of $26$\% (i.e., $-2$\%), and an $\textrm{R}^2$ score of $0.95$ (i.e., $+0.01$) for the combination.

\begin{table}
\caption{Performance summary of the prediction and approximation models on the flexibility envelope displayed in Fig.~\ref{fig:envelope}}
\label{tab:pipeline_summary}
\centering
 \begin{tabular}{|c|c|c|c|} 
 \hline
 Model & MAE (min) & NRMSE (\%) & $\textrm{R}^2$ score \\ [0.5ex] 
 \hline
 Prediction (Benchmark) & 93 & 55 & 0.79 \\
 Prediction (AdaBoost) & 21 & 25 & 0.95 \\ 
 Approximation (SND) & 25 & 15 & 0.99 \\ 
 Combin. (AdaBoost + SND) & 35 & 28 & 0.94 \\ 
 \hline
\end{tabular}
\end{table}

\subsection{Computational Aspects}

The original flexibility envelope quantification was performed on a  server with an Intel Core i7-8565U CPU and 32~GB of RAM. The optimization problem was solved with CVXPY \cite{diamond2016cvxpy} in Python and MOSEK \cite{mosek} as the solver.
For the machine learning-based models presented in this work, the results are obtained on a different server with a similar specification of Intel Core i7-8565U CPU and 16~GB of RAM, performed in Python using the scikit-learn library \cite{scikit-learn}.

Please note that the time spent on offline training is not accounted for since it does not influence the real-time flexibility quantification process. Moreover, the trained prediction algorithms are valid as long as the system configuration remains unchanged. However, training is again required if a new flexible component is added to the system or if substantial changes in the flexibility quantification model or inhabitant behavior are made.

\section{Introduction}\label{sec:introduction}
\input{sections/introduction}

\section{Methodology}\label{sec:methodology}
\input{sections/methodology}

\section{Case study}\label{sec:case_study}
\input{sections/case_study}

\section{Results}\label{sec:results}
\input{sections/results}

\section{Conclusion and outlook}\label{sec:conclusion}
\input{sections/conclusion}

\section*{Acknowledgement}
The work of H. Cai and P. Heer was supported by the Swiss Federal Office of Energy through the Sustainable Demand-Side Management for the Operation of Buildings (S-DSM) project under the contract number SI/502165-01. The work of T. Zufferey and G. Hug is part of the activities of the Renewable Management and Real-Time Control Platform (ReMaP).

\bibliographystyle{IEEEtran}




\subsection{Original Flexibility Envelope Quantification}\label{sec:original_flex_quant}

Initially, the upper and lower energy bounds are identified by energizing the building to its extremes, which consists of solving two optimization problems as follows \cite{caha2021}:
\begin{subequations}\label{eq:opt_based_quant}
\begin{alignat}{1}
    J_{\mathrm{upper}} &= \min_{u,\delta} \; - w^\intercal u + w_\delta \cdot \lVert \delta \rVert^2,\\
    J_{\mathrm{lower}} &= \min_{u,\delta} \; w^\intercal u + w_\delta \cdot \lVert \delta \rVert^2,
\end{alignat}
\end{subequations}
where $u$ is a vector of active power of all components, $w$ is a vector of weights, $\delta$ is a vector of slack variables introduced to ensure feasibility and $w_\delta\in\mathbb{R}_+$ is the corresponding weight. Optimization constraints include the dynamics of each component, power and energy capacities, state limits and mutual exclusiveness of charging and discharging. The respectively obtained $u$ gives the upper and lower energy bounds. 
The feasible power levels and the corresponding maximum sustained duration can be derived using the procedure shown in Fig.~\ref{fig:envelope_calculation} \cite{caha2021}. When applying this derivation in all the lead times within the horizon, a flexibility envelope spanned by lead time, power levels, and maximum sustained duration is obtained.  An example of flexibility envelope with a 24-hour horizon is shown in Fig.~\ref{fig:envelope}.
In real-time applications, the entire quantification process is periodically repeated to account for the most recent observations.
In the rest of this paper, the prediction and approximation of flexibility envelopes concretely refer to the prediction and approximation of the maximum sustained duration values, which depend on the two other dimensions.

\begin{figure}
    \centering
    \includegraphics[width=\columnwidth]{./figures/flex_concept}
    \caption{Calculation of maximal sustainable duration by (a) first identifying energy bounds and then (b) deriving sustainable duration of chosen power levels.}
    \label{fig:envelope_calculation}
\end{figure}

\subsection{Prediction of Flexibility Envelope}

Each flexibility envelope can be predicted by~\ac{ML} models in order to dispense with the computationally intensive solution of the optimization problem described in (\ref{eq:opt_based_quant}). 
In this work, a simple benchmark model is compared to a linear regressor, a~\ac{KNN} algorithm, a~\ac{SVM}, and an AdaBoost regressor from the scikit-learn Python library. They are presented in detail in the following subsections. Hyperparameter tuning is performed using grid search with 5-fold cross-validation on the training set.
Predictions are performed based on the following input features:
\begin{itemize}
    \item calendar data: time of the day (cyclic encoding)
    \item weather forecast from MeteoSwiss: solar irradiance (in~$\text{W/m}^2$) and temperature (in~$\degree \text{C}$)
    \item residential behavior from historical data: presence at home (binary), non-flexible load (in~$\textrm{kW}$), water draw (in~$\textrm{L/min}$)
    \item comfort requirements: room temperature limits (in~$\degree \text{C}$), domestic hot water temperature limits (in~$\degree \text{C}$)
    \item predicted energy state of flexible appliances: home battery energy (in~$\textrm{kWh}$),~\ac{EV} battery energy (in~$\textrm{kWh}$), \ac{PV} energy production (in~$\textrm{kWh}$), room temperature  (in~$\degree \text{C}$), domestic hot water temperature  (in~$\degree \text{C}$)
\end{itemize}
Note that the input features are a subset of the information required by the base model predictive controller and no additional data are necessary. Furthermore, the linear regressor, KNN, and SVM algorithms require features scaling, which is done by standardization in this work:
\begin{equation}\label{eq:scaling}
    x' = \frac{x - \mu}{\sigma},
\end{equation}
where $x$ and $x'$ are the original and scaled feature values, respectively. In addition, $\mu$ and $\sigma$ are the mean and standard deviation of the original feature values, respectively.

\subsubsection{Linear Regressor}

The linear model outputs a linear combination of the input features:
\begin{equation} \label{eq:lin}
	 \hat{y} = \sum_{i=1}^m \beta_i\cdot x_i
\end{equation}
where $\hat{y}$ is the predicted value, $x_i$ is the corresponding $i^{\textrm{th}}$ input feature, $\beta_i$ is the weight associated with the $i^{\textrm{th}}$ input feature, and $m$ is the total number of features. In this study, the weights are tuned based on the Lasso regression from the Lasso Python class. The training process consists of minimizing a least square loss function with L1 regularization:
\begin{equation}\label{eq:lin_loss}
    J_{\mathrm{lin}} = \min_{\beta} \; \lVert y - X\beta \rVert^2 + \gamma \cdot \lVert \beta \rVert
\end{equation}
where $y$ is a vector of length $n$ containing the target values in the training set, $X$ is an $n \times m$-matrix containing the input features, $\beta$ is a vector of length $m$ containing the weights of the linear regression, and $n$ is the total number of examples in the training set. The regularization parameter $\gamma$ is a hyperparameter whose value is selected from $\{0.2, 0.4, 0.6, 0.8, 1\}$. Both Lasso and Ridge (i.e., L2 regularization) regressions handle collinearity well \cite{dormann2013collinearity}. Nevertheless, the Lasso regression has been chosen over the Ridge regression since it performs better on the training set.

\subsubsection{K-Nearest Neighbors}

The~\ac{KNN} is a non-parametric model whose prediction is a linear combination of the target values of the $k$ nearest neighbors in the training set:
\begin{equation} \label{eq:knn}
	 \hat{y} = \sum_{i=1}^k w_i^{\mathrm{knn}}\cdot y^{\mathrm{knn}}_i,
\end{equation}
where $\hat{y}$ is the predicted value, and $y^{\mathrm{knn}}_i$ and $w_i^{\mathrm{knn}}$ are the target value and the weight associated to the nearest neighbor $i$, respectively. In this work, the model comes from the KNeighborsRegressor Python class, where the number of nearest neighbors $k$ is a hyperparameter whose value is selected from $\{20, 40, 60, 80, 100, 120, 140, 160\}$. Standard Euclidean distance is used as the similarity metric.

\subsubsection{Support Vector Machine}

Used as a regressor, the~\ac{SVM} is also a non-parametric approach that consists of two steps.  In the first step, the original input features are mapped into a higher dimensional space induced by a so-called kernel function. In this new space, the computation occurs in a more efficient manner. In the second step, the algorithm determines the hyper-plane with the narrowest possible margin with a so-called $\varepsilon$-insensitive loss through the input data. Formally, this consists of solving the following convex optimization problem:
\begin{subequations}
\begin{alignat}{3}
&\min_{w,\xi_i,\xi_i^*} \quad && \dfrac{1}{2}\cdot \lVert w \rVert^2 + C \sum_{i=1}^{n} (\xi_i +  \xi_i^*), \\
&\,\,\,\mathrm{s.t.}  && w^\intercal x_i - y_i \geq \varepsilon + \xi_i,\\
&               && y_i-w^\intercal x_i  \geq \varepsilon + \xi_i^*,\\
&               && \xi_i,\xi_i^*,\varepsilon \geq 0, \quad \forall i\in \{1,...,n\},
\end{alignat}
\end{subequations}
where $x_i$ and $y_i$ are a vector of input features (in the higher dimensional space) and the target value corresponding to the $i^{\textrm{th}}$ training sample, respectively. In addition, $w$ is a vector of weights, $C$ is a regularization parameter, $\varepsilon$ defines the margin of insensitive loss, and $n$ is the total number of samples in the training set. Next, $\xi_i$ and $\xi_i^*$ are slack variables introduced to ensure feasibility. The goal is to find a function that stays within an $\varepsilon$-margin from the predicted values $y$, but also exhibits flatness. The first requirement is met through minimization of the $\varepsilon$-insensitive loss. The second requirement is reflected in the L2-regularization of the weights. There is an extensive literature on the SVR algorithm and more information can be found in~\cite{drucker1997support}. In this work, the model comes from the SVR Python class, where $\varepsilon$ is set to $0.1$ and the radial basis function is chosen as Kernel function. Finally, $C$ is a hyperparameter whose value is selected from $\{0.1, 0.2, 0.5, 1, 2, 5\}$.

\subsubsection{AdaBoost Regressor}

The AdaBoost approach is an ensemble method (i.e., multiple weak learners are combined to create a strong learner), where the training of the $(n+1)^\textrm{th}$ weak learner depends on the results of the $n^\textrm{th}$ weak learner. Hence, each weak learner focuses on the error of the previous weak learner. The final prediction is a weighted average of the prediction of all weak learners, where the respective weights are calculated based on the error of each weak learner. More information on the AdaBoost regressor is available in~\cite{freund1997decision}. In this work, the model comes from the AdaBoostRegressor Python class, and the weak learners consist of decision tree regressors. The number of weak learners and the maximal depth of decision trees are hyperparameters whose values are selected from $\{20, 50, 100\}$ and $\{3, 5, \textrm{none}\}$, respectively.

\subsubsection{Benchmark Model}

The \ac{ML}-based predictors are also compared to a simple benchmark model which exploits the time-dependency and the daily pattern of flexibility. Concretely, the benchmark model divides the day into intervals of the same length. For each interval, the average sustained duration is computed for each power level based on the data in the training set. These average values are used as predictions for the sustained duration in the test set:
\begin{equation} \label{eq:benchmark}
	 \hat{y}_j = \frac{1}{|\Omega|}\sum_{i\in\Omega} y_i, \quad \forall j\in\Omega,
\end{equation}
where $\hat{y}_j$ and $y_i$ are the predicted and historical values whose characteristics comply with the condition of $\Omega$, i.e. correspond to a given period of the day. The length of the intervals is a hyperparameter selected from $\{0.25\textrm{h}, 0.5\textrm{h}, 1\textrm{h}, 2\textrm{h}, 4\textrm{h}, 8\textrm{h}\}$ and varies among the various power levels.




\subsection{Approximation of Flexibility Envelope}

The approximation of flexibility envelopes via~\acp{PDF} allows for a reduction of the number of parameters to be transmitted to \acp{DSO}, however at the cost of decreased accuracy. In this work, the flexibility envelopes are approximated based on the 2-dimensional (2D)~\ac{ND}, the 2-dimensional~\ac{SND}, and the 3-dimensional~\ac{GMM}. The 2D models fit a \ac{PDF} to a slice of the flexibility envelope, i.e. at a given lead time within the horizon (see Fig.~\ref{fig:envelope_calculation}(b)), whereas the 3D model directly approximates the entire envelope. Here, the horizon refers to the entire time window (i.e., set of all considered lead times) of the flexibility envelope.

\subsubsection{2D Normal Distribution}

A flexibility envelope at a given lead time tends to follow a normal (or Gaussian) distribution:
\begin{equation}\label{eq:nd}
    \phi(x)  = \frac{1}{\sigma\sqrt{2\pi}}\cdot e^{-\frac{(x-\mu)^2}{2\sigma^2}},
\end{equation}
where $\sigma$ and $\mu$ are the standard deviation and the mean of the~\ac{PDF}.

In this work, the approximation model consists of a weighted sum of two~\acp{ND}, which shows improved performance over a single~\ac{ND}. Formally, the model is defined as:
\begin{equation}\label{eq:2Dsum}
    f(x)  = \sum_{i=1}^2 c_i \cdot \phi_i(x),
\end{equation}
where $\phi_i(x)$ is a~\ac{ND} as defined in~(\ref{eq:nd}), and $c_i$ is the corresponding weight. Parameters are optimized by the~\textit{curve.fit} function from the scipy.optimize Python library, which relies on the Levenberg-Marquardt algorithm~\cite{more1978levenberg}.

\subsubsection{2D Skewed Normal Distribution}

Flexibility envelopes appear to be mostly asymmetric. Indeed, buildings import electricity most of the time such that the flexibility potential is generally skewed to the positive side. Hence, the skewed normal distribution is also considered:
\begin{equation}\label{eq:snd}
    \phi^\textrm{skew}(x)  = \frac{2}{\omega} \cdot \phi\left(\frac{x - \xi}{\omega}\right) \cdot \Phi\left(\alpha \cdot \frac{x - \xi}{\omega}\right),
\end{equation}
where $\phi(x)$ is the~\ac{PDF} of the standard~\ac{ND} as defined in~(\ref{eq:nd}) and $\Phi(x)$ is the~\ac{CDF} of the~\ac{ND}. In addition, $\alpha$, $\xi$, and $\omega$ represent the skewness, the location, and the scale parameters, respectively.

Similar to the previous approximation model, this model considers the superposition of two~\acp{PDF} as defined in~(\ref{eq:2Dsum}), where $\phi_i$ is replaced by $\phi_i^\textrm{skew}$.

\subsubsection{3D Gaussian Mixture Model}

A~\ac{GMM} is the linear combination of a finite number of normal (or Gaussian) distribution components. Its 3D version can be expressed as follows:
\begin{equation}\label{eq:3D-GMM}
    \textrm{GMM}^{\textrm{3D}}(x)  = \sum_{i=1}^K c_i \cdot \phi^{\textrm{3D}}_i(x),
\end{equation}
where $\phi^{\textrm{3D}}_i(x)$ is the $i\textsuperscript{th}$ 3D Gaussian component, and $c_i$ is the respective weight. In contrast to its 2D version, each 3D Gaussian component is defined by 2 means and a $2\times 2$ covariance matrix. In this work, the covariance matrix is diagonal such that only 2 parameters are necessary. In addition, the number of components $K$ can be set such as $18$ to balance between approximation accuracy and computational cost. Equation~(\ref{eq:3D-GMM}) is an extension of the model presented in~(\ref{eq:2Dsum}). Furthermore, a preliminary analysis shows that the performance can be improved if the values at the first and the last lead time are duplicated prior to model fitting. Therefore, this empirical insight is incorporated in the subsequent results.

\subsection{Evaluation Metrics}

For a comprehensive assessment, it is crucial to evaluate a model based on several metrics which can account for different aspects of its performance. Detailed in the following subsections, the error metrics considered in this work are the~\ac{MAE},~\ac{NRMSE}, and $\textrm{R}^2$ score. The number of parameters and the computation time are also considered as key performance indicators in the context of this study.

\subsubsection{Mean Absolute Error}
The~\ac{MAE} accounts for the absolute difference between the estimated and the real sustained duration:
\begin{equation}\label{eq:MAE}
    \text{MAE}  = \frac{\sum_{i=1}^{n} \mid \hat{y}^{i} - y^{i} \mid}{n},
\end{equation}
where $y^{i}$ and $\hat{y}^{i}$ are the real and estimated value of data point $i$, respectively. In addition, $n$ refers to the total number of data points under consideration. The lower the~\ac{MAE}, the better the performance.

\subsubsection{Normalized Root Mean Square Error}

The~\ac{NRMSE} is the normalized version of the~\ac{RMSE} which penalizes larger errors more:
\begin{equation}\label{eq:NRMSE}
    \text{NRMSE}  = \frac{1}{\bar{y}} \cdot \sqrt{\frac{\sum_{i=1}^{n} (\hat{y}^{i} - y^{i})^2}{n}},
\end{equation}
where $\bar{y}$ is the average value over all real data points under consideration. The lower the~\ac{NRMSE}, the better the performance. In contrast to the~\ac{MAE}, the~\ac{NRMSE} enables to fairly compare the outcome at different power levels even if the sustained duration values are not necessarily in the same order of magnitude.

\subsubsection{$\textrm{R}^2$ Score}

The $\textrm{R}^2$ score, also referred to as coefficient of determination, explains how much of the variance in the dataset can be estimated by the model. It is defined as:
\begin{equation}\label{eq:R2}
    \textrm{R}^2  = 1 - \frac{\sum_{i=1}^{n} ( \hat{y}^{i} - y^{i})^2}{\sum_{i=1}^{n} (y^{i} - \bar{y})^2}.
\end{equation}
The coefficient of determination is upper bounded at 1, but has no lower bound. The higher the score, the better the performance. A value of zero corresponds to the performance of a model that always outputs the mean value of the test set.

\subsection{Prediction of Flexibility Envelope}

\begin{table}
\caption{Performance summary of the prediction models}
\label{tab:prediction_summary}
\centering
 \begin{tabular}{|c|c|c|c|c|}
 \hline
 \begin{tabular}{@{}c@{}}Predict.\\model\end{tabular} &
 \begin{tabular}{@{}c@{}}MAE\\(min)\end{tabular} &
 \begin{tabular}{@{}c@{}}NRMSE\\(\%)\end{tabular} &
 $\textrm{R}^2$ score &
 \begin{tabular}{@{}c@{}}Comput.\\time (sec)\end{tabular}\\ [0.5ex] 
 \hline
 Original & - & - & - & $20.75\pm 16.05$  \\
 \hline
 Benchmark & 70 & 46 & 0.67 & 0.15 \\
 Linear & 84 & 61 & 0.49 & 0.19\\
 KNN & 42 & 35 & 0.78 & 0.59\\
 SVM & 42 & 30 & 0.81 & 1.48\\
 AdaBoost & 25 & 21 & 0.89 & 0.52\\
 \hline
\end{tabular}
\end{table}

The benchmark model and the four \ac{ML}-based models have been tuned on the training set and applied to the samples in the test set, which consists of 24 hours of data. Their performance is summarized in Table~\ref{tab:prediction_summary} in terms of~\ac{MAE} (in minute),~\ac{NRMSE} (in percentage), $\textrm{R}^2$ score, and average computation time (in second, per flexibility envelope). AdaBoost exhibits by far the best prediction accuracy according to all error metrics. In contrast, the linear model performs particularly badly, even worse than the benchmark model. In terms of computation time, AdaBoost is the second fastest \ac{ML} model after the linear model, which makes it the most appropriate candidate. In comparison, the original calculation using one CPU process is estimated to take $20.75\pm 16.05$ 
seconds per flexibility envelope most of the time. The reduction in computational efforts is achieved because \ac{ML}-based approaches do not involve an iterative process in searching for the solution. 

Moreover, Figs.~\ref{fig:prediction_MAE},~\ref{fig:prediction_NRMSE}, and~\ref{fig:prediction_R2} provide insights into the error distribution with respect to the power level. First of all, the ranking of performance observed in Table~\ref{tab:prediction_summary} among the various prediction models is respected for a majority of power levels. In terms of absolute error (i.e., MAE), power levels close to zero are naturally more impacted due to the generally higher sustained duration values. In the case of \ac{NRMSE}, however, the accuracy tends to decrease when moving to the sides of the flexibility envelope. A notable exception is the power levels between $-3$~kW and $1$~kW, where the relative error increases substantially. This is principally due to the normalization. The total flexibility in that power range is mainly driven by the PV flexibility which has strong diurnal patterns but a relatively low average over the entire day.
Finally, the $\textrm{R}^2$ score indicates that the variance is generally harder to predict between the power levels of $-2$ kW and $8$ kW. Indeed, such power levels are a composition of small loads and PV output, characterized by high variance. In contrast, the flexibility outside $[-2,8]~\textrm{kW}$ is mainly attributed to the~\ac{EV}, which has more regular patterns.

\begin{figure}[t!]
    \centering
    \includegraphics[width=\columnwidth]{./figures/prediction_MAE}
    \caption{MAE of the flexibility prediction with respect to the power level.}
    \label{fig:prediction_MAE}
\end{figure}

\begin{figure}[t!]
    \centering
    \includegraphics[width=\columnwidth]{./figures/prediction_NRMSE}
    \caption{NRMSE of the flexibility prediction with respect to the power level.}
    \label{fig:prediction_NRMSE}
\end{figure}

\begin{figure}[t!]
    \centering
    \includegraphics[width=\columnwidth]{./figures/prediction_R2}
    \caption{$\textrm{R}^2$ score of the flexibility prediction with respect to the power level.}
    \label{fig:prediction_R2}
\end{figure}

\subsection{Approximation of Flexibility Envelope}

The three approximation models have been applied to all the original flexibility envelopes in the dataset. Their performance is summarized in Table~\ref{tab:approximation_summary} in terms of number of parameters, average~\ac{MAE} (in minute), average~\ac{NRMSE} (in percentage), average $\textrm{R}^2$ score, and average computation time (in second) per flexibility envelope. In terms of accuracy, the 3-dimensional model is clearly outperformed by both 2-dimensional models. Indeed, each slice of envelope (i.e., corresponding to a given lead time) exhibits a distinctive Gaussian shape, but the specific parameters highly vary over the different lead times. Moreover, the skewed normal distribution achieves a slightly higher performance than the standard version. We can observe from Figs.~\ref{fig:approximation_MAE},~\ref{fig:approximation_NRMSE}, and~\ref{fig:approximation_R2} that the 3D model further leads to a larger spread of the approximation errors. It should also be noted that the 2D models, and especially the~\ac{SND}, can almost perfectly account for the variance in the data. In contrast, the 3D-GMM outperforms both 2D approximators in terms of computational time (e.g., $8.7$ times faster than the 2D-ND) and required number of parameters. In fact, the 2D-ND, 2D-SND, and 3D-GMM allow for a reduction of the number of parameters by a factor of $8.83$, $6.63$, and $12.96$, respectively, compared to the original representation.
Hence, the most appropriate algorithm depends on the preference between approximation accuracy, computational burden, and communication load.

\begin{table}[t!]
\caption{Performance summary of the approximation models.}
\label{tab:approximation_summary}
\centering
 \begin{tabular}{|c|c|c|c|c|c|} 
 \hline
 \begin{tabular}{@{}c@{}}Approx.\\model\end{tabular} &
 \begin{tabular}{@{}c@{}}Number of\\parameters\end{tabular} &
 \begin{tabular}{@{}c@{}}MAE\\(min)\end{tabular} &
 \begin{tabular}{@{}c@{}}NRMSE\\(\%)\end{tabular} &
 $\textrm{R}^2$ score &
 \begin{tabular}{@{}c@{}}Comput.\\time (sec)\end{tabular}\\ [0.5ex] 
 \hline
 2D-ND & 132 & 32 & 28 & 0.97 & 3.68\\
 2D-SND & 176 & 23 & 23 & 0.98 & 5.50\\ 
 3D-GMM & 90 & 66 & 72 & 0.83 & 0.63\\ 
 \hline
\end{tabular}
\end{table}

\begin{figure}[t!]
    \centering
    \includegraphics[width=\columnwidth]{./figures/approximation_MAE}
    \caption{Distribution of the MAE for the flexibility envelope approximation.}
    \label{fig:approximation_MAE}
\end{figure}

\begin{figure}[t!]
    \centering
    \includegraphics[width=\columnwidth]{./figures/approximation_NRMSE}
    \caption{Distribution of the NRMSE for the flexibility envelope approximation.}
    \label{fig:approximation_NRMSE}
\end{figure}

\begin{figure}[t!]
    \centering
    \includegraphics[width=\columnwidth]{./figures/approximation_R2}
    \caption{Distribution of the $\textrm{R}^2$ score for the flexibility envelope approximation.}
    \label{fig:approximation_R2}
\end{figure}

\subsection{Overall Comparison}

\begin{figure*}
    \centering
    \includegraphics[width=\textwidth]{./figures/example}
    \caption{Example of the prediction and approximation of the flexibility envelope displayed in Fig.~\ref{fig:envelope}.}
    \label{fig:example}
\end{figure*}


Fig.~\ref{fig:example} shows the outcome of the benchmark prediction model as well as the best performing prediction (i.e., AdaBoost) and approximation (i.e., 2D skewed normal distribution) models for the example flexibility envelope from Fig.~\ref{fig:envelope}. The upper subplots display the predicted and approximated flexibility envelopes, whereas the lower subplots show the difference between the respective upper subplots and the original envelope. While the benchmark prediction model cannot account for the variation among different lead times, both best performing models provide an envelope visually close to the original envelope.

In addition, the best performing prediction and approximation models are integrated into a single pipeline, which is displayed in the right-hand side subplots. Formally, the 2D-SND approximates the envelope predicted by the AdaBoost algorithm. This corresponds to the process desired in this work to address both the computational burden and communication load. 

Table~\ref{tab:pipeline_summary} summarizes the performance of the different models visualized in Fig.~\ref{fig:example}. As expected, the benchmark prediction model exhibits by far the worst performance. In absolute values, the entire envelope can be better predicted than approximated. However, the approximation still leads to lower relative error and higher $\textrm{R}^2$ score. More importantly, although the combination performs slightly worse than both models separately, the respective errors do not directly add up. This is a very encouraging outcome and shows that these various models can be used to considerably speed up the flexibility quantification process and lower the communication load without substantial loss of accuracy.

Finally, a hybrid prediction model has also been considered as an alternative to the overall best performing model. Precisely, instead of predicting all power levels based on the AdaBoost algorithm, model selection is carried out in the training phase such that the best performing prediction model per power level is applied on the test set (see Fig.~\ref{fig:prediction_MAE}). Consequently, AdaBoost is still used in a large majority of the cases but is replaced by the~\ac{SVM} and~\ac{KNN} algorithms for some power levels close to zero. This hybrid approach allows for a MAE of 31 minutes (i.e., $-4$ minutes, compared to the accuracy displayed in Table~\ref{tab:pipeline_summary}), a NRMSE of $26$\% (i.e., $-2$\%), and an $\textrm{R}^2$ score of $0.95$ (i.e., $+0.01$) for the combination.

\begin{table}
\caption{Performance summary of the prediction and approximation models on the flexibility envelope displayed in Fig.~\ref{fig:envelope}}
\label{tab:pipeline_summary}
\centering
 \begin{tabular}{|c|c|c|c|} 
 \hline
 Model & MAE (min) & NRMSE (\%) & $\textrm{R}^2$ score \\ [0.5ex] 
 \hline
 Prediction (Benchmark) & 93 & 55 & 0.79 \\
 Prediction (AdaBoost) & 21 & 25 & 0.95 \\ 
 Approximation (SND) & 25 & 15 & 0.99 \\ 
 Combin. (AdaBoost + SND) & 35 & 28 & 0.94 \\ 
 \hline
\end{tabular}
\end{table}

\subsection{Computational Aspects}

The original flexibility envelope quantification was performed on a  server with an Intel Core i7-8565U CPU and 32~GB of RAM. The optimization problem was solved with CVXPY \cite{diamond2016cvxpy} in Python and MOSEK \cite{mosek} as the solver.
For the machine learning-based models presented in this work, the results are obtained on a different server with a similar specification of Intel Core i7-8565U CPU and 16~GB of RAM, performed in Python using the scikit-learn library \cite{scikit-learn}.

Please note that the time spent on offline training is not accounted for since it does not influence the real-time flexibility quantification process. Moreover, the trained prediction algorithms are valid as long as the system configuration remains unchanged. However, training is again required if a new flexible component is added to the system or if substantial changes in the flexibility quantification model or inhabitant behavior are made.
