\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.


\section{Introduction}

\paragraph{Anatomical versus Diagnostic} Anatomical (organ-level) and diagnostic (pathology-level) segmentation are fundamentally different, with respect to their clinical objectives. For instance, within the clinical workflow of prostate MRI, segmentation of the prostate anatomy \{whole-gland (WG), transitional zone (TZ), peripheral zone (PZ)\} is used to estimate its volume, boundaries and geometry --enabling the calculation of prostate-specific antigen (PSA) density, guiding treatment planning and future interventions \cite{NEJM,PSAdPlanning,RadioTherap}. Meanwhile, segmentation of clinically significant prostate cancer (csPCa) is primarily used to detect the number of malignant lesions present in the prostate gland (if any), and characterize each instance for diagnosis, risk stratification and/or targeted biopsies (similar to PI-RADS guidelines \cite{PIRADSv2}). In other words, while the former anatomical task leans towards \textit{semantic} single-object quantification, the latter diagnostic task can be framed as \textit{panoptic} multi-object detection \cite{Panoptic,PanopticPCa}. Despite these differences, in medical image computing, these objectives are often treated and trained as near-identical segmentation tasks using U-Nets, and evaluated in a similar manner using the Dice Similarity Coefficient (DSC) \cite{DSC1,DSC2,DSC3,Jaeger}. Recent studies confirm that while DSC can be appropriate for single-object quantification, its inability to measure multi-object detection makes it unsuitable for lesion localization \cite{CommonPitfalls,DSC_LIM,DSCB,Panoptic}. We hypothesize that such considerations should not only be made during evaluation, but should also be incorporated at train-time by using task-specific loss function(s).

\vspace{-2mm}

\paragraph{Segmentation Loss Functions} We can distinguish between loss functions that minimize mismatch in underlying distributions and those that minimize mismatch in segmentation regions/boundaries \cite{LossOdyssey}. Distribution-based loss functions, such as the widely adopted cross-entropy (CE) loss, minimize dissimilarity between the target distribution (ground-truth) and its approximation (model predictions), over the course of training. Balanced cross-entropy (BCE) loss is an extension of CE, introducing class weights $\alpha$ to address imbalanced class frequencies in the dataset. Focal loss (FL) represents a further extension of BCE, introducing an additional hyperparameter $\gamma$ to differentiate between easy and difficult examples \cite{FocalLoss}. When $\gamma=0$, FL reduces down to CE or BCE depending on the value of $\alpha$. On the other hand, region or boundary-based loss functions are derived from integrals over segmentation regions (e.g. soft Dice loss \cite{VNet}) or distance metrics in the space of contours (e.g. boundary loss \cite{BoundaryLoss}). Unlike distribution-based losses, which can benefit from class re-weighting, the sum of Dice and boundary loss (DB) is implicitly impervious to class imbalance, and in turn, recommended for highly imbalanced datasets \cite{BoundaryLoss}.

\vspace{-1mm}

\section{Experiments and Analysis}

\vspace{-1.5mm}

\paragraph{Materials} We used 200 prostate bpMRI (T2W, high b-value DWI, ADC) exams from the publicly-available ProstateX dataset \cite{ProstateX}, paired with voxel-level delineations of WG, TZ, PZ and csPCa \cite{ProstateXAnnot}. All images were resampled to $0.5\times0.5\times3.0$ mm$^3$ resolution, center-cropped to $160\times160\times20$ voxels and intensity-normalized (T2W, DWI: \textit{z}-score; ADC: linear) \cite{E2E}, prior to usage.

\vspace{-1mm}

\paragraph{Bayesian Segmentation Model} We used a probabilistic adaptation \cite{ProbUNet} of the deep attentive 3D U-Net, developed and validated specifically for prostate bpMRI in our previous work \cite{E2E}. Monte-Carlo dropout nodes were added to capture both \textit{epistemic} and \textit{aleatoric} uncertainty during inference (as recommended by \citet{MCProb}). Cosine annealing learning rate \cite{CosineAnnLRWR} (decaying from $10^{\textnormal{-}4}$ to $10^{\textnormal{-}7}$) and \textit{AMSGrad} optimizer \cite{AMSBound} were used to train the model. Data augmentations comprised of additive Gaussian noise (standard deviation $0$-$0.5$), horizontal flip, rotation ($\pm$$7.5$\textdegree), translation ($0$-$15$\% horizontal and/or vertical shifts) and scaling ($0$-$15$\%) centered along the axial plane. Fig. 1 illustrates the train-time schematic of the model, and its source code is made publicly-available at \textcolor{magenta}{\url{https://github.com/DIAGNijmegen/prostateMR_3D-CAD-csPCa}}.

\vspace{-1mm}

\paragraph{Experiments} We investigated four different loss policies (CE, BCE, FL, DB) as the segmentation loss function ($\textbf{\textit{L}}_\textbf{S}$) in our model for anatomical (\{WG, TZ, PZ\} in T2W MRI) and diagnostic (csPCa in bpMRI) segmentation. For the former task, $\alpha$ is set as [$0.05$, $0.30$, $0.65$] as per the empirical class distribution of the dataset, and for the latter task, $\alpha$ is set as [$0.75$, $0.25$] as per the findings of previous studies \cite{UCLA,E2E}. In both cases, $\gamma$ is set to its default value of $2.00$ as used in \cite{FocalLoss,UCLA,E2E}. All metrics were computed in 3D over 3 independent runs $\times$ 5-fold cross-validation $\times$ mean of 100 executions of probabilistic inference per image. Identical data splits were maintained for all configurations across both tasks.

\begin{figure}[b!]
\centering
\includegraphics[width=0.975\textwidth]{figures/fig1.pdf}
\caption{Train-time schematic for the Bayesian segmentation model. $\textbf{\textit{L}}_\textbf{KL}$ denotes the Kullback--Leibler divergence loss between prior distribution $\textbf{\textit{P}}$ and posterior distribution $\textbf{\textit{Q}}$. $\textbf{\textit{L}}_\textbf{S}$ denotes the segmentation loss between prediction $\textbf{\textit{p}}$ and ground-truth $\textbf{\textit{Y}}$. For each execution of the model, one sample \textbf{z} $\in \textbf{\textit{Q}}$ (train-time) or \textbf{z} $\in \textbf{\textit{P}}$ (test-time) is drawn to predict one segmentation mask $\textbf{\textit{p}}$ \cite{ProbUNet}.}
\label{fig1}
\end{figure}

\newpage

\section{Results and Discussion}

\begin{figure}[h!]
  \centering
  \begin{tabular}{ccc}
  \centering
  \includegraphics[width=0.29\textwidth]{figures/nll_z.pdf} & \includegraphics[width=0.29\textwidth]{figures/sdsc.pdf} & \includegraphics[width=0.29\textwidth]{figures/hdsc.pdf} \\
  {\small \hspace{0.04\textwidth} \textbf{(a)}} & {\small \hspace{0.04\textwidth} \textbf{(c)}} & {\small \hspace{0.04\textwidth} \textbf{(d)}} \\
  \end{tabular}

  \vspace{1mm}

  \begin{tabular}{ccc}
  \centering
  \includegraphics[width=0.29\textwidth]{figures/nll_l.pdf} & \includegraphics[width=0.29\textwidth]{figures/auroc.pdf} & \includegraphics[width=0.29\textwidth]{figures/froc.pdf} \\
  {\small \hspace{0.04\textwidth} \textbf{(b)}} & {\small \hspace{0.04\textwidth} \textbf{(e)}} & {\small \hspace{0.04\textwidth} \textbf{(f)}} \\
  \end{tabular}
  \centering
  \includegraphics[width=0.9\textwidth]{figures/legend.pdf}
  \caption{Anatomical \textbf{(top-row)} and diagnostic \textbf{(bottom-row)} segmentation of \{WG, TZ, PZ\} in T2W MRI and csPCa in bpMRI, respectively. For each $\textbf{\textit{L}}_\textbf{S}$ configuration, we evaluated $\textbf{(a, b)}$ model calibration with negative log-likelihood (NLL), $\textbf{(c, d)}$ segmentation quality with mean soft/hard Dice Similarity Coefficient (DSC) across all 600 observations (3 runs $\times$ 5 folds $\times$ 40 validation samples), $\textbf{(e)}$ patient-level diagnosis with Area Under Receiver Operating Characteristic (AUROC), $\textbf{(f)}$ lesion-level localization with Free-Response Receiver Operating Characteristic (FROC). Transparent areas in $\textbf{(f)}$ FROC indicate $95$\% confidence intervals. Arrows ($\uparrow$ or $\downarrow$) across $\textbf{(a-f)}$ indicate whether the corresponding metric should be maximized/minimized for ideal performance. }

  \vspace{-4mm}
\label{fig2}
\end{figure}

\paragraph{Calibration} We find that distribution-based loss functions (CE, BCE, FL) produce implicitly well calibrated predictions, in comparison to DB \cite{TMI_1} (as shown in Fig. 2 (a, b) and \textit{Appendix A}). Among them, FL exhibits notably better calibration. Our findings concur with that of \citet{FLeffects}, where the authors demonstrate that FL is formulated as such, that it minimizes dissimilarity between ground-truth and predicted distributions, while \textit{increasing entropy} of the latter --thereby limiting NLL overfitting and inducing both calibration and regularization. Such a property can also lead to more useful uncertainty estimates \cite{FLeffects,UnCa}.

\vspace{-3mm}

\paragraph{Anatomical Segmentation} From Fig. 2 (c), it is clear that while CE, BCE and DB achieve similar segmentation quality with respect to raw \textit{softmax} predictions (\{WG: $81$-$89$; TZ: $82$-$85$; PZ: $71$-$74$\} mean DSC), FL performs worse (\{WG: $70$; TZ: $64$; PZ: $50$\} mean DSC). When considering binarized \textit{argmax}(\textit{softmax}) predictions, as seen in Fig. 2 (d), overall FL performance is only marginally lower than that of CE, BCE and DB (\{WG: $88$-$90$; TZ: $84$-$86$; PZ: $72$-$76$\} mean DSC). We attribute this decline in FL performance to its high predictive uncertainty near class boundaries (as seen in \textit{Appendix A}), which limits its ability to produce precise contour definitions.

\vspace{-3mm}

\paragraph{Diagnostic Segmentation} From Fig. 2 (e, f), we observe that while all loss functions are similar in patient-level diagnosis ($76$-$78$ mean AUROC), results indicate that distribution-based loss functions, especially FL, perform substantially better than region/boundary-based losses (DB) at lesion detection (with an average increase of $15.7\%.$ in maximum detection sensitivity and $0.25$ in partial area under FROC between $0.1$-$2.5$ false positives per patient). We attribute this improvement to the implicit train-time calibration induced by FL, which can facilitate better risk stratification and higher detection sensitivities across difficult diagnostic tasks (e.g. via retaining low-confidence lesion predictions, as opposed to muting them out --refer to \textit{Appendix A}). We believe that the growing success of csPCa detection models trained using FL derivatives (such as the \textit{FocalNet} \cite{UCLA}, among others \cite{E2E,ISBI_Cao,ProCDet}) can, in part, be attributed to this phenomenon rather than architectural enhancements or class weighting alone (note BCE performs near-identical to CE, if not worse, across both tasks). On the other hand, while DB loss is a natural fit for maximizing DSC, it is ill-posed to optimize panoptic segmentation objectives \cite{Panoptic,CommonPitfalls}. Subsequently, its miscalibration (refer to Fig. 2 (a, b)) or polarized predictions (refer to \textit{Appendix A}), translate to a relatively flat FROC curve with lower maximum detection sensitivity (refer to Fig. 2 (f)) --presumably detection of low-confidence, difficult or small lesions are skipped in favour of maximizing confidence and overlap of clear lesions \cite{CommonPitfalls}.

In conclusion, we recommend distribution-based loss functions (in particular, FL) for diagnostic or panoptic segmentation tasks. For anatomical or semantic segmentation tasks, we observe that although distribution-based losses improve calibration, (with the exception of FL) they are equivalent to region/boundary-based losses in terms of standard/hard DSC. Hence, CE, BCE, DB or presumably their composites (as applied by the \textit{nnU-Net} \cite{nnUNet} and investigated by J. Ma et al. \cite{LossOdyssey}), can be recommended. Further analyses are required, using larger datasets and multiple medical imaging modalities, to draw out definitive conclusions.

\vspace{-2mm}

\section*{Broader Impact}

\vspace{-2mm}

Prostate cancer is one of the most prevalent cancers in men worldwide \cite{PCaStat2019}. In the absence of experienced radiologists, its multifocality, morphological heterogeneity and strong resemblance to numerous non-malignant conditions in MR imaging, can lead to low inter-reader agreement ($<50\%$) and sub-optimal interpretation \cite{mpMRI,LimitedPIRADS2,LimitedPIRADS3}. The development of automated, reliable detection algorithms has therefore become an important research focus in medical image computing ---offering the potential to support radiologists with consistent quantitative analysis, improve diagnostic accuracy, and in turn, minimize unnecessary biopsies in patients \cite{UnnBiop,RadPRvsDL}.

To the best of our knowledge, this study has no foreseeable negative societal effects.

\vspace{-2mm}

\begin{ack}

\vspace{-2mm}

This research is supported, in parts, by the European Union H2020: ProCAncer-I project (EU grant 952159), Health{\textasciitilde}Holland (LSHM20103), ContextVision AB (Linköping, Sweden) and Siemens Healthineers (CID: C00225450).
\end{ack}

{
\fontsize{8.75pt}{8.75pt}\selectfont
\bibliographystyle{unsrtnat}