Chance constrained day-ahead robust flexibility
needs assessment
for low voltage distribution network

Md Umar Hashmi, Member, IEEE 0000-0002-0193-6703, Arpan Koirala, Graduate Student Member IEEE 0000-0003-4826-7137, Hakan Ergun, Senior Member, IEEE 0000-0001-5171-1986, and Dirk Van Hertem, Senior Member, IEEE 0000-0001-5461-8891
Corresponding author email<EMAIL_ADDRESS>Md U. Hashmi, A. Koirala, H. Ergun,  and D. Van Hertem are with KU Leuven, division Electa & EnergyVille, Genk, Belgium

This work is supported by the H2020 EUniversal project, grant ID: 864334 (<https://euniversal.eu/>) and the energy transition funds project BREGILAB organized by the FPS economy, S.M.E.s, Self-employed and Energy.

For market-based procurement of low voltage (LV) flexibility, DSOs
identify the amount of flexibility needed for resolving probable distribution network (DN) voltage and thermal congestion.
A framework is required to avoid over or under procurement of flexibility in the presence of uncertainty.
To this end, we propose a scenario-based robust chance-constrained (CC) day-ahead flexibility needs assessment (FNA) framework.
The CC level is analogous to the risk DSO is willing to take in flexibility planning.
Multi-period optimal power flow is performed to calculate the amount of flexibility needed to avoid network issues.
Flexibility is defined in terms of nodal power ramp-up and ramp-down and cumulative energy needs over a full day for each node.
Future uncertainties are considered as multiple scenarios generated using multivariate Gaussian distribution and Cholesky decomposition.
These scenarios are utilized to solve the flexibility needs assessment optimal power flow (FNA-OPF) problem.
Zonal clustering of an LV feeder is performed using electrical distance as a measure and spatial partitioning.
The FNA tool calculates ramp-up and ramp-down flexibility's power and energy requirements.
Energy and power needs are often valued differently in many energy markets. We identify the marginal value of flexibility associated with energy and power needs separately.
From numerical results for an LV feeder, it is observed that zonal flexibility needs assessment is more immune to uncertainty than nodal flexibility needs, making it more useful for DSOs to evaluate day-ahead flexibility procurement.
We also propose a Pareto optimal mechanism for selecting CC level to reduce flexibility needs while reducing DN congestion.

Flexibility needs assessment (FNA), scenario generation, chance constrained (CC), optimal power flow (OPF), distribution system operator (DSO), zonal clustering

§ INTRODUCTION

With growing uncertainty in distribution network (DN) operation due to distributed generation (DG) and loads with high simultaneity factor such as electric vehicle charging, it is crucial for DSOs to plan flexible resources to ensure reliable operation of DNs.
ENTSO-e defines flexibility "as the ability of the power system
to cope with variability and uncertainty in demand, generation
and grid capacity, while maintaining a satisfactory level of
reliability at all times" [1].
Flexibility assists power networks to damp fluctuations, which could potentially reduce the reliability of the DN [2, 3].
Authors in [4] propose a computational method to calculate prosumer multi-period flexibility forecast. Using these feasible flexibility spaces, prosumers can participate in the energy market. Reference
[5] shows that small electricity prosumers can utilize their flexibilities for providing grid services while making a profit in day-ahead (DA) market. Authors in [6] schedule flexible resources for DA energy dispatch, which minimizes the energy procurement cost.

[7] presents a case study of Northern Germany, where every MW of additional wind installation will require 0.7 MW of flexibility for avoiding curtailment of renewable generation.
Prior works [8, 9, 10] present market mechanisms for end-user flexible resources for solving low voltage (LV) network issues.
Market mechanism presented in [8] utilizes contractual flexibility by DSOs for H2020 Interflex project.
Authors in [9] propose an efficient design for flexibility markets.

Although flexibility market design is not the focus of the work, interested readers can refer to [11], [12].
In [13] authors observe that locational temporal flexibility mapping is crucial for avoiding future network issues. They propose a mapping tool for Northern Ireland through GIS analysis.

Activation of flexibility in LV DN could aim at (a) mitigating LV DN issues [14] and/or (b) providing services to MV and HV transmission or DNs. The latter would require coordination in the operation of flexible resources, ensuring it does not create new problems [15].
The focus of this work is to utilize flexible resources for solving LV DN issues.
Flexibility needs assessment (FNA) is fairly a new concern for power networks. Recent works such as [16, 17, 18, 19] discuss the need for flexibility provisioning. [16] utilize FNA for prioritization of future investments. Authors in [17, 19] quantifies flexibility provisions based on conditional value at risk.

The goal of this paper is to quantify the flexibility needs of a DN from a DSO's perspective. The DSO utilizes the FNA output provided by the framework proposed in this work for procurement of flexible resources in DA flexibility market implemented in H2020 EUniversal project.
The procurement of such flexible resources in the DA market can be costly.
Under-procuring can lead to network congestion, and over-procuring of flexible resources may not be efficient.
FNA should consider the probable uncertainties which could happen.
Advanced time planning utilizes one or a combination of  (a) probabilistic modelling [20], (b) scenario generation [21], or (c) robust planning [22, 23, 24].

We propose a scenario-based robust DA planning mechanism for quantifying the flexibility needs of a DN.

Chance constrained (CC) levels are utilized to avoid planning for the worst-case scenario.
The DSO also needs to consider how a selection of a CC level would project on probable DN congestion.

Power network clustering is widely utilized for operation and planning of transmission and DNs [25, 26, 27].
In [25, 26] network clustering is used for transmission system operator level resource dispatch and planning.
Authors in [27] use zonal voltage regulation in the context of a LV DN with high penetration of DG.
In this work, we identify zones using spectral partitioning, along with $k$-means clustering for a doubly stochastic measure matrix formulated using a static network admittance matrix.
The proposed framework for zonal clustering partitions the DN in connected segments, a key challenge in other partitioning methods.
These zones are used to aggregate nodal flexibility needs.
Often, exact flexibility needs may not be met by the resources available in the flexibility market. In such a case, zonal FNA can be utilized by DSOs to find alternative resources which could resolve the network congestion.
Advantages of flexibility aggregation are also observed in [28].

In robust optimization, an uncertainty set is defined, over which a robust feasible solution is identified. It is unclear how we can define such an uncertainty set for FNA problem. Thus, we utilize generated scenarios for creating bounds on uncertainty. These scenarios define the uncertainty set for the flexibility needs assessment optimal power flow (FNA-OPF).
The bounded uncertainty and its associated worst-case action selection may lead to too conservative solution

(see Chapter 2 of [29]).
A simple scenario generation model is proposed which takes a point forecast and its associated forecast error as input for generating nodal load and generation scenarios.

The flexibility needs which provides a feasible solution for all scenarios will often lead to an over-procurement of flexible resources. These values of flexibility needs will be robust in true sense. To avoid worst case planning, chance-constrained is used. CC level is analogous to the risk level DSO is willing to consider in flexibility planning.
Chance constraint is associated to the FNA.

The generated scenarios are used as snapshots of Monte Carlo simulations.

FNA-OPF calculates the flexibility needs for solving LV DN issues while minimizing the cost of curtailed load and generation.
Furthermore, spatio-temporal metrics of FNA are proposed.

Numerically, it is observed that variance of the day-ahead FNA can be reduced by zonal aggregation of flexibility needs. This implies aggregation makes it more predictable for the DSO to plan flexibility in a day-ahead setting. The zonal FNA can help DSOs to identify alternative flexible resources which could probably solve DN congestion issues.
We also propose a Pareto optimal mechanism for selecting the level of CC used in FNA-OPF.
The selection of CC level for proposed FNA of a DN should consider DSO's risk for over and under-procurement of flexible resources.
In multi-objective optimization, some of the objectives may be in conflict with others. In such a case, an optimal solution for one objective could lead to unacceptable outcomes for other objectives. In power systems, prior works such as [30, 31, 32] utilize multi-objective optimization. Authors in [32] use epsilon-constrained Pareto front for an electric vehicle scheduling problem of an aggregator, while considering DSO's goal to minimize losses at the same time reducing the charging cost. In our work, we use the Pareto optimal value of chance-constraint to reduce DSO's risk of over-procurement of flexible resources while reducing the probable DN congestion incidents.

The key contributions of the paper are:

$\bullet$ Uncertainty consideration:

We propose a simple scenario generation model using multivariate Gaussian distribution and Cholesky decomposition,  which takes as input PV size at a node, forecasted nodal load profiles, forecast of normalized PV generation,  and the associated forecast errors for load and solar generation. The forecast error creates a boundary of uncertainty. The generated scenarios serve as input for the FNA-OPF.

$\bullet$ Flexibility needs assessment:
The spatial and temporal robust flexibility needs of a DN are identified for avoiding all network incidents for all scenarios generated. For risk-averse FNA of the DN, a chance constrained nodal flexibility level is applied. Further, the nodal FNA is extended to zonal FNA by aggregating flexibility needs in a zone.

$\bullet$ Case studies: Two dedicated numerical case studies are presented. The first case study shows the FNA output of a DN. A Pareto optimal mechanism is provided for tuning CC level. The selection of CC level is projected on to unavoided network incidents using power flow simulations. The FNA of a DN is aggregated for providing zonal FNA.
The second case study quantifies the marginal value of flexibility based on energy or power needs. It is observed that for the DN, power needs are twice as important as the energy needs.

Structure of this paper and the inputs and the outputs of robust chance-constrained FNA under uncertainty.

The structure of the paper and the inputs and outputs of the proposed FNA framework are shown in Fig. <ref>.

In Section <ref>, we detail the methodology to generate scenarios.

In Section <ref>, zonal clustering of a DN based on electrical distance using spatial partitioning is provided.
In Section <ref>, the metrics used in FNA are quantified.
In Section <ref>, the FNA-OPF is detailed.
Section <ref> presents the numerical results.

Section <ref> concludes the paper.

§ SCENARIO GENERATION

The focus of the present paper is to consider probable future uncertainties to provide FNA. This tool is developed as part of the EUniversal project to be implemented by the distribution system operator Mitnetz Storm in Germany. Due to German data limitations, only historical forecast errors and point forecasts are available for scenario generation. This has led us to opt for a simple scenario generation model with data limitation.

The generated scenarios using

multivariate Gaussian distribution and Cholesky decomposition

considers (a) correlation of data points, (b) parameter distribution.
Cholesky factor is utilized to form unique scenarios

emulating different uncertainties.

This structure of scenario generation

can be used for historical data and point forecasts with known forecast errors.

In the present case, large historical data is unavailable, instead a point forecast along with its associated forecast error is known. In such a case, a multivariate Gaussian distribution can be formed by using point time-series forecast as the mean value and the variance of the distribution is proportional to forecast error [33].

§.§ Multivariate distribution and Cholesky decomposition

The multivariate Gaussian distribution is a generalization of the one-dimensional normal distribution to more than one dimension.
For our case, $X_i$ denotes the load profile of consumers connected at time $i$.
The multivariate normal distribution of a $S$-dimensional random vector $\mathbb{X} = (X_1, X_2, .... , X_S)^T$ can be written as
\mathbb{X} \sim \textit{N} (\mu, \Sigma),$

with $S$-dimensional mean vector given as
\mu =  \textbf{E}[\mathbb{X}] =
\mleft[ \begin{array}{c}
\textbf{E}[\mathbb{X}_1] , ..,
\textbf{E}[\mathbb{X}_S]
\end{array} \mright]^T,
and the covariance matrix is given as
\Sigma  =
\mleft[ \begin{array}{ccc}
\texttt{Cov}(\mathbb{X}_1,\mathbb{X}_1)  & .. & \texttt{Cov}(\mathbb{X}_1,\mathbb{X}_S) \\
%  \texttt{Cov}(\mathbb{X}_2,\mathbb{X}_1) & \texttt{Cov}(\mathbb{X}_2,\mathbb{X}_2) & .. & \texttt{Cov}(\mathbb{X}_2,\mathbb{X}_S) \\
: &  \ddots & :\\
\texttt{Cov}(\mathbb{X}_S,\mathbb{X}_1)  & .. & \texttt{Cov}(\mathbb{X}_S,\mathbb{X}_S) \\
\end{array} \mright],
where $\texttt{Cov}(\mathbb{X}_i,\mathbb{X}_j)$ (also denoted as $\Sigma_{i,j}$) denote the covariance between $\mathbb{X}_i$ and $\mathbb{X}_j$ and is given as

\begin{equation}
\Sigma_{i,j} = \textbf{E}[(\mathbb{X}_i - \mu_i)(\mathbb{X}_j - \mu_j)], ~i,j \in\{1,..,S\}.
\end{equation}

Each component $X_i$ has distribution $N(\mu_i, \sigma_i^2)$ with $\Sigma_{ii} = \sigma_i^2$.
The covariance matrix is a square matrix.

To qualify as a covariance matrix, $\Sigma$ must be symmetric ($\Sigma = \Sigma^T$) and positive semidefinite ($x^T\Sigma x \geq 0 $).

In order to ensure that Cholesky decomposition does not fail, a small error is added to the covariance matrix for numerical reasons, e.g. $
\Sigma = \Sigma + \epsilon I_{S},
where $I_{S}$ denotes identity matrix of order $S$, $\epsilon$ denotes small error value.
This ensures that eigenvalues of $\Sigma$ do not decay rapidly, which stabilizes the decomposition. Due to the small magnitude of $\epsilon$, it has inconsequential effects on the samples while ensuring numerical stability [34, 35].

Cholesky decomposition is used to calculate the lower triangular matrix, $L$, for the covariance matrix such that
\Sigma =  L L^T.
$L$ is also called the Cholesky factor.
A lower triangular matrix is convenient because it reduces the calculation of a scenario as $\mu + LZ$ where $Z \sim N(0,I)$ to

\begin{gather*}
X_1 = \mu_1 + l_{11} z_1,\\
X_2 = \mu_2 + l_{21}z_1 + l_{22}z_2, ~\ldots\\
% :\\
X_S = \mu_S + l_{S1}z_1 + l_{S2}z_2 + \ldots +l_{SS}z_S
\end{gather*}

where $z_1, z_2, \ldots, z_S$ are
independent and identically distributed random variables.

The Cholesky factor is utilized to form unique scenarios [36] and is calculated using eigenvalue factorization.
Since $\Sigma$ is symmetric and positive semidefinite, therefore, its eigenvalues $\lambda_1, \ldots, \lambda_S$ are non-negative.
The covariance matrix can be denoted as $\Sigma = V \wedge V^T$, where $V$ is an orthogonal matrix, i.e. $VV^T =I$, $\wedge$ is a diagonal matrix with eigenvalues as diagonal entries. In this case, $L=V\wedge^{0.5}$.

Consider $J$ scenarios are generated for a time-series with $S$ intervals. The scenario matrix is given as

\begin{equation}
S_J =
\mleft[ \begin{array}{ccc}
\mu_1 & ..&\mu_1 \\
\mu_2 & ..&\mu_2 \\
: & ..&:\\
\mu_S & ..&\mu_S
\end{array} \mright]_{S,J}
+ \mleft[ \begin{array}{ccc}
l_{11}& .. & 0 \\
l_{21}& .. & 0 \\
:&  \ddots & :\\
l_{S1}&.. & l_{SS} \\
\end{array} \mright]
\texttt{rnd}(S,J),
\label{eq:scenariosgen}
\end{equation}

where $\texttt{rnd}(S,J) \sim N(0,I)$ denotes uniformly distributed random numbers.
(<ref>) is used to generate $J$ number of unique forecast profiles.

§.§ Temporal and spatially correlated scenarios

The proposed scenario generation tool considers the point forecast as the mean shown in
The covariance matrix is a function of forecast error calculated based on historical data.
The forecast error is assumed to be 1.96 times the standard deviation of the normal distribution. This would ensure that more than 95% of the incidents are covered. The extreme tail events are ignored in the scenario generation.
The scenario generation does not consider the spatial correlation. However, the generated scenarios do not ignore spatial correlation entirely.
The net load of a node consists of a solar generation as DG and load. Since consumption is assumed positive, therefore, the load is strictly non-negative and DG is negative. The load profile is sampled separately depending on load forecast error and its point forecast, both of which are assumed to be known.
On the other hand, DG scenarios are calculated based on its point forecast and associated forecast errors. The DG forecast scenarios are normalized in per kW installed capacity. These scenarios are assumed to be the same for all DGs in the feeder, assuming all DGs are installed in geographical proximity to be not affected by different levels of solar irradiance.
The average spread of 160 Spanish feeders, [37], is a square of size of 158 meters.
Therefore, this assumption would be valid for most DNs.

Fig. <ref> shows input and output of scenario generation.

Scenario generation decoupling load and generation to consider temporal and spatial correlationship.

§ IDENTIFYING ZONES OF A LV DN

Identifying the zones of a LV DN will help the DSO in planning the flexibility needs of a network. Although, the analytical division of a DN is expected to work well, however, due to the sheer number of DN feeders, it becomes crucial to have a standardized framework for dividing DN into zones based on electrical and/or geographical distances.
For example, in the UK, there are more than 1.3 million LV feeders [38].
In this section, we develop a clustering framework to identify the best-suited LV DN zonal partition using electrical distance as a measure.

The proposed zone formation solves following challenges: (i)
The formation of connected zones requires an incidence matrix based measure, therefore, we consider admittance as a measure. (ii) The admittance matrix cannot be used directly, therefore, spectral decomposition of doubly stochastic matrix is used. (iii)  The adequate number of zones is not known a priori. However, to apply unsupervised clustering techniques such as k-means, one should know how many clusters are needed. We use silhouette score as a measure for identifying the best number of clusters.

Although, the different methods mentioned above have been well discussed in the literature, in this work, they have been for the first time combined to perform zonal FNA of a DN.

§.§ Distance measure for clustering

A DN can be represented in standard directed graph notation $G=(\mathscr{N},E, W)$, where $\mathscr{N}$ denotes the set of nodes such that $|\mathscr{N}| = D$ ($D$ denotes the total number of nodes), $E$ denotes edges and $W$ denotes edge weight matrix.
We aim to divide $\mathscr{N}$ into $p$ groups $\{M_1, M_2,...,M_p \}$ such that $M_i \subset \mathscr{N}$ and $M_i \cap M_j = \emptyset, ~ \forall i \neq j$.
The edge weights can be considered as a penalty for cutting that line in the graph while clustering.
It also measures the connection strength of two nodes in a graph.
A good measure would separate the weakly connected nodes and identify nodes that should be clustered together.

A power network can be clustered based on static network parameters such as line impedance/admittance or dynamic parameters such as power flows, line losses, voltage variations [12, 26, 39].
The electrical distance matrix is used for partitioning LV DNs in
[25, 26]. Similar to these prior works, we use the line admittance matrix as the distance measure parameter. The weight $w_{ij} = Y_{ij} = 1/|R_{ij} + jX_{ij}|$, where $Y_{ij}$, $R_{ij}$ and $X_{ij}$ denotes admittance, resistance and impedance between nodes $i$ and $j$ and the edge weight $w: \mathscr{N}\times \mathscr{N} \rightarrow \mathbb{R}^{\geq 0}$ such that (a) $w_{ij} = w_{ji}, ~ \forall i, j$, (b) $w_{ij} = 0, \text{ if } (i,j) \notin E$, (c) $w_{ii} = 0, \forall i$.

§.§ Spectral partitioning

Spectral clustering is used for power network partitioning or creation of zones or network reduction in [27, 26].
Previously described weight matrix

cannot be directly used, as the diagonal elements are all zeros. Authors in [26] use the normalized Laplacian matrix for spectral partitioning.
Tutorial [40] transforms the weight matrix into a doubly stochastic matrix.
A double stochastic matrix is a special type of Markov matrix where not only each row but also each column add to 1.
Spectral properties of doubly stochastic matrices have been detailed in [41].
For this transformed matrix, all eigenvalues are real and smaller than or equal to 1, with one eigenvalue exactly equal to 1.
For identifying $k$ partitions in a graph, the $k$ highest eigenvalues and corresponding orthonormal
eigenvectors are identified.
The eigenvector matrix of the order $N \times k$ is used for DN partitioning, in effect reduces the dimensionality of the problem.
$k$-means clustering is used to partition the spectral data points.

§.§ Goodness of a cluster

Previously, we detailed partitioning a given LV DN into $k$ zones. In this subsection, we deal with identifying the best-suited values of $k$. In a real-world LV DN network partitioning problem, we may not know how many clusters we want. We use performance indices for measuring the goodness of a partition and, based on different values of $k$, identify the best value which fits our needs.
In this work, the goodness of a cluster is measured using the mean silhouette index of the network cluster.
The silhouette coefficient of a node is a confidence indicator of its association in a group $M_x$ [27, 42].
Rousseeuw [43] proposed an interpretation based on value of silhouette coefficient.
The zone selection algorithm is detailed in Appendix <ref>.

§ DAY-AHEAD FLEXIBILITY PLANNING

Flexibility is referred to as resources that can be activated or deactivated, thus in effect increasing or decreasing the consumption based on the grid's needs.
For a DN, nodal voltages and line loadings should remain within operational bounds. Violation of these limits could damage appliances, affect DN operation, or grid elements.
Authors in [11] present two business models for operating electrical flexibility in the context of services provided to transmission and distribution system operators and assisting them to accommodate a greater amount of renewable energy.
Depending on the application of flexibility, flexible resources needs can be measured in terms of response speed (or ramp rate), duration (or energy) and direction (ramp up and ramp down power) [44].
Next, the framework

for calculating the flexibility needs that DSO procures in the energy market is detailed.

§.§ Flexibility: definition

Prior work, [45], used three parameters for defining flexibility: (a) ramp rate, (b) power, and (c) energy.

The flexibility model used in this work utilizes the power and energy parameters. This model resembles the simple storage model used in [46].
The active power flexibility at node $i$ and time $t$ is

given as

\begin{equation}
\Delta P_{i,t}^{\text{flex}} = \Delta P_{i,t}^{\text{flex}+} + \Delta P_{i,t}^{\text{flex}-},
\label{eq:flexdef}
\end{equation}

where $\Delta P_{i,t}^{\text{flex}+}$ and $\Delta P_{i,t}^{\text{flex}-}$ denotes ramp down and ramp up flexibility, respectively.

Ramp down flexibility decreases the nodal load and is analogous to load curtailment.
Ramp up flexibility increases the nodal load and is analogous to generation curtailment.
Ramp up flexibility can be provided by consumer loads such as HVAC, water heaters, pool pumps [47] etc.
Ramp down flexibility can be provided by solar generation, and also by curtailment of flexible loads.
Prosumer energy storage if not fully charged or discharged can provide both ramp up and ramp down flexibility [48].
Note these flexible resources have a limited ramp rate, power and energy to be utilized as flexibility. As an example, the temperature of a water heater can be temporally constrained to be within some bounds of water temperature. If utilized as ramp up flexibility to increase the nodal load, it can only be operated till water temperature hits the upper bound (or the lower bound for ramping down).
Thus, to respect technical device and user comfort constraints, we use power and energy constraints to define flexible resources.
The ramp down and ramp up flexibility energy for node $i$ are given as

\begin{equation}
\begin{split}
E_{i}^{\text{flex}+} = \sum_t \Delta P_{i,t}^{\text{flex}+} \Delta t,~~
E_{i}^{\text{flex}-} = \sum_t \Delta P_{i,t}^{\text{flex}-} \Delta t,
\end{split}
\end{equation}

\begin{gather}
\label{eq:flexrange1}
\Delta P_{i,t}^{\text{flex}+} \in [0, \Delta P_{i,t, \max}^{\text{flex}+}], ~~
% \label{eq:flexrange2}
\Delta P_{i,t}^{\text{flex}-} \in [\Delta P_{i,t, \min}^{\text{flex}-}, 0], ~ \forall t,\\
\label{eq:flexrange3}
E_{i}^{\text{flex}+} = [0, E_{i, \max}^{\text{flex}+}],~~~~
E_{i}^{\text{flex}-} = [E_{i, \min}^{\text{flex}-}, 0].
\end{gather}

Since ramp up and ramp down flexibility levels can compensate each other, therefore, they are dealt with separately in order to correctly assess the DN needs.
Next, we describe the metrics used to quantify flexibility needs.

§.§ Metrics for day-ahead FNA

The limits on ramp up and down power and energy needs are infrastructure constraints that the optimization needs to take into account.
The needs assessment should provide the temporal and locational needs of the DN.
With these inputs and dimensions to the needs assessment problem, the following metrics are utilized:

* Objective function: the optimization minimizes the cost of operation of flexible resources while ensuring network issues are avoided. Assuming ample amount of flexible resource available, there are many solutions that can solve DN issues, some of which could have the same objective function value.
In this work, we utilize a flat penalty factor for penalizing the use of ramp down and ramp up flexibilities. The penalty factors are denoted as $\lambda^{\text{rampDown}}$, $\lambda^{\text{rampUp}}$. These factors are fixed for all nodes across all times.

* Flexibility power needs assessment metrics: Ramp up ($R_{i,t}^{P_{\text{up}}})$ and ramp down ($R_{i,t}^{P_{\text{down}}})$ power (with CC) needs for node or zone $i$ and time $t$: the temporal distribution of ramp up and ramp down power needs will often have a fat-tailed distribution, as shown in Fig. <ref>. Planning for the absolute worst case may lead to drastic over-design of flexibility needs.

* Flexibility energy needs assessment metrics: Ramp up ($R_{i}^{E_{\text{up}}})$ and ramp down ($R_{i}^{E_{\text{down}}})$ energy with chance constrained needs for node or zone $i$ denotes the energy needs of DN.

FNA metrics for ramp up and down power and energy using chance constraint level denoted as $\epsilon_{cc}$

* The nodal and zonal FNA are compared using the average weighted variance of nodal ramp up and down power and energy needs for all $t$.

Tuning CC level ($\epsilon_{cc}$) for FNA: The selection of $\epsilon_{cc}$ is crucial for DSO to appropriately plan for day-ahead flexibilities. Both over and under procurement is not efficient. We project $\epsilon_{cc}$ on mean DN voltage and thermal congestion over all generated scenarios. Next, the Pareto optimal level of CC is utilized for reducing DN flexibility activation and reducing DN congestion events. Pareto optimality is frequently used to select between conflicting goals of a multi-objective optimization problem [49].

§ FLEXIBILITY NEEDS ASSESSMENT

§.§ Notation

A power network consists of nodes ($\mathscr{N}$),
branches ($E$), generators and loads.

Each node $i \in \mathscr{N}$ has two variables, i.e., voltage magnitude ($V_{i,t}$) and phase angle ($\theta_{i,t}$) for time $i$ which depends on power injection and load magnitude. The branch admittance $(i,j) \in E$ together with the nodal voltage differences govern the power flow and losses.
The subset of nodes with loads connected is denoted as $\mathscr{N}_L \subset \mathscr{N}$. The active and reactive power demand in these nodes are denoted as $P^d_{i,t}$ and $Q^d_{i,t}$.
Nodes with generators connected is denoted as $\mathscr{N}_G \subset \mathscr{N}$, have active and reactive power generation denoted as $P^g_{i,t}$ and $Q^g_{i,t}$.

§.§ Optimization formulation

The flexibility needs assessment is performed by solving the optimal power flow problem which minimizes the cost of dispatching flexible resources.
We solve a multi-period optimal power flow problem with penalties associated with generation and load curtailment, refer to (<ref>). Further, multi-period coupled energy curtailment constraint is considered along with instantaneous curtailed power constraint, refer to  (<ref>). The optimization problem is denoted as flexibility needs assessment optimal power flow dispatch or FNA-OPF. FNA-OPF is an adaptation of the resource dispatch tool introduced in [50].

\begin{equation}
\label{eq:porgobjective}
({\text{FNA-OPF}})~\underset{\substack{\Delta P^{\text{flex}}_{i,t}}}{\text{min}}~ \sum_i \Big( \lambda^{\text{rampDown}}  E_i^{\text{flex}+} -
\lambda^{\text{rampUp}} E_i^{\text{flex}-} \Big)
%   \label{objectivefunctionFNA}
\end{equation}

subject to,  (<ref>),  (<ref>)

\begin{equation}
~V_{\min}^i \leq |V_{i,t}| \leq V_{\max}^i, ~ \forall i \in \mathscr{N} , t \in \{1,..,T\},% ~~ \text{voltage limit}
\label{eq:voltage}
\end{equation}

\begin{equation}
(P^g_{i,t}) - (P^d_{i,t} - \Delta P^{\text{flex}}_{i,t}) + jQ^d_{i,t} = \sum_j s_{i,j,t},  ~~\forall~ i,j \in \mathscr{N} ,% ~~ \text{power balance}
\label{eq:powerbalance}
\end{equation}

\begin{equation}
|s_{i,j,t}| \leq s_{i,j}^{\max}, ~~\forall~ i,j \in \mathscr{N}     \label{const:thermal}, % ~~ \text{thermal limit}
\end{equation}

\begin{equation}
P^g_{i,t} \in [P^g_{\min,i},P^g_{\max,i}] , ~~\forall~ i \in \mathscr{N_G}     \label{const:genlimit}, % ~~ \text{thermal limit}
\end{equation}

\begin{equation}
s_{i,j,t}  = \textbf{Y}_{ij}^*V_{i,t}V_{i,t}^* - \textbf{Y}_{ij}^*V_{i,t}V_{j,t}^*, ~ \forall (i,j) \in E \cup E^R, % ~~ \text{Ohm's law}
\label{eq:ohmslaw}
\end{equation}

\begin{equation}
\angle (V_{i,t} V_{j,t}^*) \in [\theta_{i,j}^{\min}, \theta_{i,j}^{\max}], ~~\forall~ i,j \in \mathscr{N}, %~~ \text{phase angle constraint}
\label{const:phase}
\end{equation}

where $s_{i,j,t}$ denote the apparent power flow from node $i$ to $j$ at time $t$.

(<ref>), (<ref>) and (<ref>) denote the voltage constraint for nodes, thermal limit constraint and phase angle constraints for branches, respectively. (<ref>) denotes the nodal balance of active and reactive power in the network.
(<ref>) denotes the generator output power limits.
(<ref>) denotes Ohm's law.

To avoid that both the ramp up and down flexibility described in (<ref>) are non-zero simultaneously they have been defined as $\Delta P_{i,t}^{\text{flex}+}\geq 0, \Delta P_{i,t}^{\text{flex}-}\leq 0$. As such a change in either of them increases the objective function value in (<ref>).
Active power of the generation is not included in the objective function, as
this would influence the amount of flexibility activation.
Especially, optimization with several components which may come in conflict, it is not clear how these components will affect each other.
Furthermore, the active power dispatch is an entirely different problem. The goal of this study is to find out the flexibility needs of a DN provided active power is already dispatched.

§.§ Robust FNA under uncertainty

Robust optimization is a two-stage optimization problem where the inner level minimizes the objective function and the outer level maximizes the inner-optimization within an uncertainty set [22, 23, 24, 29].

For the robust FNA, we utilize scenarios for defining the uncertainty set.
The inner min-problem solves the FNA-OPF for all generated scenarios. The outer level max-problem finds the robust flexibility levels for all time steps at a day-ahead level.
The ramp down and ramp up flexibility needs of a DN is modelled as an empirical cumulative distribution function (ECDF), given as

\begin{equation}
\begin{split}
F_{i,t,+}(z) & = \frac{1}{S} \sum_{j=1}^S \mathbf{I}_{(-\infty,z]}(\Delta P^{\text{flex}+}_{i,t} \text{for scenario } j), \forall i,t, \\
F_{i,t,-}(z) & = \frac{1}{S} \sum_{j=1}^S \mathbf{I}_{(-\infty,z]}(\Delta P^{\text{flex}-}_{i,t} \text{for scenario } j),\forall i,t, \\
\end{split}
\end{equation}

where identity function is given as

\begin{equation}
\mathbf{I}_{(-\infty,z]}(y)=
\begin{cases}
1 ,& \text{if }y\leq z ,\\
0 , & \text{otherwise,}
\end{cases}
\end{equation}

Based on ECDF, the chance constrained robust flexibility needs of the distribution network are identified as below

\begin{equation}
\begin{split}
(\Delta P^{\text{flex}+}_{i,t})^* & = F_{i,t,+}^{-1}(1-\epsilon_{cc}), \\
(\Delta P^{\text{flex}-}_{i,t})^* & = F_{i,t,-}^{-1}(\epsilon_{cc}), \\
\end{split}
\label{ccequation}
\end{equation}

where $\epsilon_{cc} \in [0,1]$ denotes the chance constraint level.

The load profiles and their associated forecast error profiles are utilized by generating 1000 scenarios. These scenarios are used as input for FNA-OPF. Note that to obtain a feasible solution in all the scenarios, we solve FNA-OPF where each load and/or generation can be curtailed. The activated flexibility is analysed as shown in Fig. <ref>. The chance constraint level will denote the amount of risk in the procurement of flexibility. If the $\epsilon_{cc}$ is too low, then the risk for over-procurement is high. However, if the $\epsilon_{cc}$ is too high, then the risk of under-procurement is high, see (<ref>).
Based on a CC level, the nodal ramp up and ramp down power and energy needs of a DN are identified.
These nodal needs' assessment is aggregated into the zones identified in Section <ref>.

§ NUMERICAL CASE STUDIES

§.§ System description

The DN considered is an adaptation of one of the Spanish LV feeders described in [37].
The DN consists of 76 nodes and 75 branches connecting the nodes.
52 loads are connected in 28 DN nodes, implying many nodes have more than one load connected.
The network diagram is shown in Fig. <ref> (in Appendix <ref>).
The network used in this work along with load profiles can be found on GitHub [51].
The numerical simulations are performed using PowerModels.jl in Julia / JuMP
Since the FNA-OPF problem is nonconvex, therefore, IPOPT solver [53] is utilized to solve (<ref>). The feasibility of the solution of FNA-OPF is verified using power flows.

The load and PV distributions are selected so as the network issues are visible for significant periods.
The total load of the DN seen at substation over a day is 2.719 MWh.
Total installed PV is equal to 338 kW peak, with a cumulative PV generation over a day of 1.709 MWh.
Based on the above load profiles, referred to as nominal nodal profiles, 1000 scenarios are generated using Cholesky decomposition and shown in Fig. <ref>.
For the scenario generation, the forecast error for nodal load is assumed to be 30% and PV generation is assumed to be 40%.
Note from Fig. <ref>, that the load perturbations are reduced during aggregation at the substation, however, PV generation due to high simultaneity factor leads to significantly high fluctuations during the day.
For the nominal load profile, the self-sufficiency is 62.85% which means a substantial part of consumer load is met using PV generation.
The voltage set-point at the slack bus is set at 1.01 per unit.

Due to high PV installations and high peak loads during late evening, we observe DN voltage and thermal constraint violations.
For the 1000 scenarios generated, we observe that only 57.03% of optimal power flow calculations are feasible for all time steps using an iteration limit of 3000.
We observe under voltage (voltage $<0.95$)) for 6.165% of the samples in all scenarios, and 3.275% of samples we observe over-voltage (voltage $>1.05$).

Thermal overload is calculated by performing OPF with relaxed network constraints, as OPF output for infeasible cases are not reliable.
Power flow calculations are used for all non-feasible OPF solutions.
We observe thermal over-loads for 0.11% of samples over all samples, i.e., (number of scenarios)$\times$(number of nodes)$\times$(number of time steps in 1 scenario).

Net load seen at transformer with 1000 scenarios. (A) shows the nominal aggregate load profile without PV generation, (B) shows the nominal aggregate load profile with PV generation.

The mean silhouette score for the 76 node network based on the best cluster evaluation described is shown in Fig. <ref>. The mean silhouette score for 12 clusters is 0.8487.
The high silhouette score implies the nodes in a cluster are strongly related to each other and weakly related to other nodes in the other clusters.
Fig. <ref> shows the zones of the DN.
Note that the naming of zones is arbitrarily selected based on k-means centroid.

Silhouette score for 76 node network with 12 zones

§.§ Case study 1: Flexibility needs assessment

FNA-OPF is executed for unbounded power and energy flexibilities.
Thus, we identify the minimum amount of flexibility needed to avoid all network issues in all scenarios.
FNA based on chance constraint level of 5% is shown in Figure <ref>.
The distribution of flexibility needs is fitted using a normal distribution.
Note from Fig. <ref> that the quality of the approximation  depends on the distribution of the flexibility needs for all scenarios.
We utilize an empirical cumulative distribution function for applying chance constraint as normal distribution fit, shown in Fig. <ref>(a). From the numerical results we observed that the error between FNA  calculated by applying chance constraint level on the inverse ECDF and

on an approximate Gaussian fit is less than 1% for most points, however, few instances had an error exceeding 5%.
Since the distribution of FNA of a DN is not known, ECDF is used

for all subsequent results.

Note for ramp up flexibility, the chance constraint level is ($1-$ $\epsilon_{cc}$) level because of the sign reversal.

FNA calculation based on normal distribution fit and CC quartiles. The error between inverse ECDF and quantile of approximate normal distribution fit is 0.33% and 0.02% respectively.

Nodal flexibility power needs with $\epsilon_{cc}=0.05$.

Fig. <ref> shows the nodal ramp up and ramp down flexibility power needs of the DN.
The nodal needs are aggregated into zones of the DN. Zonal flexibility power needs are shown in Fig. <ref>.
The figure shows that the majority of generation and load flexibility is needed in zones 1, 3, 4, 5, 6, and 11. Note that these zones are located at the end of the distribution feeder. The locational disparity of prosumers

can be reduced by implementing a flexibility market.

The participation of flexible resources in such a market can create an additional revenue stream for such prosumers and can alleviate the locational disparity.

From Table <ref> it can be observed that normalized flexibility needs if aggregated in zones will have a reduction in standard deviation (STD) compared to nodal flexibility needs. The mean and STD denotes the parameters for the normal distribution fitted to the histogram of ramp up and down flexibility power needs.
We observe that in this numerical case study that the zonal flexibility needs have an STD which is $43.48\%$ and $57.85\%$ less for ramp up and down nodal flexibility needs.
Thus, zonal flexibility needs assessment is more immune to uncertainty.

The aggregation of uncertain parameters leads to the reduction in variance of uncertainty is a well known concept [54].

Zonal flexibility power needs with $\epsilon_{cc}=0.05$

Comparing aggregate nodal and zonal flexibility variation

2|c|Nodal flex need (kWh)     2|c|Zonal flex need (kWh)

Perf. index    Ramp Up     Ramp Down     Ramp Up     Ramp Down

Mean     - 250.39     328.7     - 250.39     328.7

STD     328.5     256.6     185.8     108.2

STD reduction     -     -     43.48%     57.85%

Table <ref> shows that a significant amount of ramp up and ramp down flexibility is needed for DN, with network issues less than 10% of the time.
In our numerical example, mean ramp up and down needs are around 9.2% and 12.1% respectively.
With growing DG integration, the need for flexibility planning will be crucial for maintaining the reliability of DN.

Fig. <ref> shows the impact of the chance constraint level on the flexibility energy needs. The heavy tail distribution can be avoided by a small value of CC. Choosing $\epsilon_{cc}=$ 5% can reduce ramp up energy needs by 68% and ramp down energy needs 40%, respectively.
The CC level can also be used for identifying the criticality of a  resource in ensuring DN reliability, as shown in the next section.

Flexibility needed with chance constraint level

§.§.§ Projecting CC levels on network issues

The CC levels are analogous to the risk DSOs are willing to take while planning for DN flexibility in a day-ahead setting.
It will be useful to understand how the CC levels project on avoiding probable network issues. We perform traditional OPF calculations with different levels of CC for FNA. The network voltage and thermal congestions are analysed and listed in Table <ref>.
where the mean values of probabilities over 1000 scenarios are provided.
Note for the $\epsilon_{cc}=$ 5%, the mean under-voltage probability reduced from 6.165% for no flexibility to 0.076%, a reduction of 98.77%. The mean over-voltage probability is reduced by 84.2% and thermal over-load probability is reduced by 77.9%. The mean hours of congestion (thermal and/or voltage) for $\epsilon_{cc}=$ 5% is reduced by 93.23%.
As expected, a robust FNA with $\epsilon_{cc}=0$ will eliminate all DN issues.
From Tab. <ref> we observe that CC levels cannot be directly projected on probable DN issues and will be governed by DN layout and the nodal distribution of demand. Next, we discuss the mechanism for selecting the appropriate CC-level for FNA.

Projecting CC levels on probabilities of network congestion

%     1l|c]@l@Under
%     1l|c]@l@Over
%     1l|c]@l@Thermal
load %     c]@l@% of
hours with

flexibility                                                         -       6.165                                                                                3.275                                                                              0.110                                                                                46.36

0                                                          0                                                                                0                                                                                0                                                                                0

With                                                         1          0.066       0.102       0.007       1.108

robust                                                       5                                                          0.076                                                                                0.518                                                                              0.024                                                                                3.138

CC                                                           10      0.203       0.915       0.044       6.196

FNA                                                          20                                                         0.544       1.440       0.082      12.16

planning                                                     30                                                         0.996       1.891       0.117      17.63

40                                                         2.079       2.563       0.153      26.43

§.§.§ Pareto optimal tuning of CC level

The FNA of a DN is governed by network layout and nodal load profiles.
A Pareto front is built for the two conflicting goals of reducing the amount of flexibility needed, and probable network congestion. Fig. <ref> shows that the Pareto optimal value of the CC-level for the considered DN is 5%.

Tuning CC level using Pareto optimality.

§.§ Case study 2: Marginal value of energy and power

Flexible resources are valued differently in different types of markets [55].
The goal of this case study is to quantify the marginal impact of energy and power constraints for flexible resources.
The ramp up and down energy and power constraints are tightened in 10% steps, in order to evaluate:
(a) marginal increase in objective function value of FNA-OPF,
(b) calculate the percentage feasibility of 100 scenarios,
(c) cumulative ramp up and
(d) cumulative ramp down energy needs.
These metrics are shown in Fig. <ref>.

DN flexibility needs with flexibility energy and power bound tightening; (0,0) denotes unconstrained solution of FNA-OPF.

The key observations from this study are:

* The marginal cost for ramp up and down power is more than two times that of the energy. This is observed in Fig. <ref>(a) where the marginal increase in the objective function is more than 487 % compared to the unconstrained objective function value for FNA-OPF. This marginal increase is due to a reduction in flexibility power ranges by 60%. However, for up to 80% reduction in energy constraint over a day, the marginal increase in the objective function is around 186%. Similar trends can be observed for other combinations of simultaneous tightening of power and energy constraints.

* The percentage of feasible FNA-OPF calculations is shown in Fig. <ref>(b). It is expected that the feasibility will reduce with tighter flexibility power and energy constraints.
Note that we limit the levels of power and energy bounds tightening to 0.6 and 0.8 respectively, as the majority of points are infeasible for any further tightening. The solutions with extremely high levels of tightening are unreliable, as in such a case either there is no solution possible or a maximum number of iterations is reached. In both these cases, the solution for FNA-OPF cannot be used.

* Fig. <ref>(b) also indicates the combinatorial nature of the flexibility needs assessment problem. We can observe that there are many combinations of feasible solutions that can solve the same network congestion.

* The ramp up and ramp down energy are shown in Fig. <ref>(c) and (d). It is observed that as power and energy constraints are tightened, an larger amounts of flexibility is needed somewhere else in the network to mitigate solve congestion. As the distance from the congestion location increases, the efficiency of the procured flexibility to solve congestion decreases.

§ CONCLUSION

Power system flexibility will be crucial for reliable operation of distribution network for ensuring a large scale integration of
distributed generation and new loads such as electric vehicle.
An optimization-based framework for calculating the day-ahead flexibility needs of a distribution network is presented in this work.
This flexibility needs assessment tool considers future uncertainties in the form of PV generation and demand scenarios.
The FNA-OPF identifies nodal ramp up and ramp down power and energy needs of a distribution network, which minimizes the cost of dispatching flexible resources.
Based on the distribution of flexibility needs, a chance constraint-based robust needs assessment is performed to avoid under or over procurement of such resources.
In the first case study, this needs assessment is evaluated for a zonal partition of the distribution network. The appropriate number of zones is determined based on electrical distance measure and spatial partitioning.
It is observed that zonal needs assessment is more immune to future uncertainties, and DSOs can utilize this feature along with nodal flex needs assessment for flexibility planning and operation.
The second case study illustrates the impact of the flexibility bounds for power and energy on the flexibility needs and the marginal cost of additional flexibility. We observe that the marginal impact of tighter power flexibility constraints on activation cost is more than twice compared to energy  flexibility constraint.
Finally, a mechanism for DSOs to tune the level of chance constraint for FNA-OPF, considering the flexibility needed and probable network congestion reduction, is demonstrated in this paper.

In future work, we use the need assessment for developing a strategy to select bids in flexibility markets.
The purpose of zonal analysis is to assist DSOs in finding alternative resources in the same zone in case the nodal flexibility needs identified are not available in the flexibility market.

Further, the impact analysis on power and energy flexibility constraints can be utilized to determine the best suitable location and size of flexibility sources in the network.

[1]
(2022) Entso-e 2021 technology factsheets. [Online]. Available:

[2]
H. Nosair and F. Bouffard, “Flexibility envelopes for power system operational
planning,” IEEE Transactions on Sustainable Energy, vol. 6, no. 3,
pp. 800–809, 2015.

[3]
I. F. Abdin and E. Zio, “An integrated framework for operational flexibility
assessment in multi-period power system planning with renewable energy
production,” Applied energy, vol. 222, pp. 898–914, 2018.

[4]
R. Pinto, R. J. Bessa, and M. A. Matos, “Multi-period flexibility forecast for
low voltage prosumers,” Energy, vol. 141, pp. 2251–2263, 2017.

[5]
X. Ayón, J. K. Gruber, B. P. Hayes, J. Usaola, and M. Prodanović, “An
optimal day-ahead load scheduling approach based on the flexibility of
aggregate demands,” Applied Energy, vol. 198, pp. 1–11, 2017.

[6]
E. Hooshmand, A. Rabiee, S. Jalilzadeh, and A. Soroudi, “Optimal flexibility
coordination for energy procurement in distribution networks,” IET
Renewable Power Generation, no. 6, pp. 1191–1203, 2021.

[7]
A. C. Mei$\beta$ner, A. Dreher, K. Knorr, M. Vogt, H. Zarif, L. Jürgens,
and M. Grasenack, “A co-simulation of flexibility market based congestion
management in northern germany,” in 2019 16th International Conference
on the European Energy Market (EEM).1em plus 0.5em minus 0.4emIEEE, 2019, pp. 1–6.

[8]
R. Fonteijn, T. Van Cuijk, P. H. Nguyen, J. Morren, and J. G. Slootweg,
“Flexibility for congestion management: A demonstration of a multi-mechanism
approach,” in 2018 IEEE PES Innovative Smart Grid Technologies
Conference Europe (ISGT-Europe).1em plus 0.5em minus 0.4emIEEE, 2018, pp. 1–6.

[9]
G. Tsaousoglou, J. S. Giraldo, P. Pinson, and N. G. Paterakis, “Mechanism
design for fair and efficient dso flexibility markets,” IEEE
transactions on smart grid, vol. 12, no. 3, pp. 2249–2260, 2021.

[10]
S. S. Torbaghan, N. Blaauwbroek, P. Nguyen, and M. Gibescu, “Local market
framework for exploiting flexibility from the end users,” in 2016 13th
International Conference on the European Energy Market (EEM).1em plus
0.5em minus 0.4emIEEE, 2016, pp. 1–6.

[11]
K. Coninx, G. Deconinck, and T. Holvoet, “Who gets my flex? an evolutionary
game theory analysis of flexibility market dynamics,” Applied energy,
vol. 218, pp. 104–113, 2018.

[12]
M. Zhang, Z. Miao, and L. Fan, “Power grid partitioning: Static and dynamic
approaches,” in 2018 North American Power Symposium (NAPS).1em
plus 0.5em minus 0.4emIEEE, 2018, pp. 1–6.

[13]
O. Agbonaye, P. Keatley, Y. Huang, O. O. Ademulegun, and N. Hewitt, “Mapping
demand flexibility: A spatio-temporal assessment of flexibility needs,
opportunities and response potential,” Applied Energy, vol. 295, p.
117015, 2021.

[14]
P. Arboleya, M. A. Kippke, and S. Kerscher, “Flexibility management in the
low-voltage distribution grid as a tool in the process of decarbonization
through electrification,” Energy Reports, vol. 8, pp. 248–256, 2022.

[15]
A. G. Givisiez, K. Petrou, and L. F. Ochoa, “A review on tso-dso coordination
models and solution techniques,” Electric Power Systems Research,
vol. 189, p. 106659, 2020.

[16]
S. Davidov and J. Curk, “Novel framework for investment prioritisation based
on flexibility needs assessment,” IET Generation, Transmission &
Distribution, vol. 14, no. 25, pp. 6189–6194, 2020.

[17]
A. Laur, J. Nieto-Martin, D. W. Bunn, and A. Vicente-Pastor, “Optimal
procurement of flexibility services within electricity distribution
networks,” European Journal of Operational Research, vol. 285, no. 1,
pp. 34–47, 2020.

[18]
E. Hillberg, A. Zegers, B. Herndler, S. Wong, J. Pompee, J.-Y. Bourmaud,
S. Lehnhoff, G. Migliavacca, K. Uhlen, I. Oleinikova et al.,
“Flexibility needs in the future power system,” 2019.

[19]
X. Yang, C. Xu, H. He, W. Yao, J. Wen, and Y. Zhang, “Flexibility provisions
in active distribution networks with uncertainties,” IEEE Transactions
on Sustainable Energy, vol. 12, no. 1, pp. 553–567, 2020.

[20]
Z. Lu, H. Li, and Y. Qiao, “Probabilistic flexibility evaluation for power
system planning considering its association with renewable power
curtailment,” IEEE Transactions on Power systems, vol. 33, no. 3, pp.
3285–3295, 2018.

[21]
Y. Chen, Y. Wang, D. Kirschen, and B. Zhang, “Model-free renewable scenario
generation using generative adversarial networks,” IEEE Transactions
on Power Systems, vol. 33, no. 3, pp. 3265–3275, 2018.

[22]
A. J. Conejo and X. Wu, “Robust optimization in power systems: a tutorial
overview,” Optimization and Engineering, pp. 1–23, 2021.

[23]
L. Baringo and A. J. Conejo, “Offering strategy via robust optimization,”
IEEE Transactions on Power Systems, vol. 26, no. 3, pp. 1418–1425,

[24]
D. Bertsimas and A. Thiele, “Robust and data-driven optimization: modern
decision making under uncertainty,” in Models, methods, and
applications for innovative decision making.1em plus 0.5em minus
0.4emINFORMS, 2006, pp. 95–122.

[25]
S. Blumsack, P. Hines, M. Patel, C. Barrows, and E. C. Sanchez, “Defining
power network zones from measures of electrical distance,” in 2009
IEEE Power & Energy Society General Meeting.1em plus 0.5em minus
0.4emIEEE, 2009, pp. 1–8.

[26]
R. J. Sánchez-García, M. Fennelly, S. Norris, N. Wright, G. Niblo,
J. Brodzki, and J. W. Bialek, “Hierarchical spectral clustering of power
grids,” IEEE Transactions on Power Systems, vol. 29, no. 5, pp.
2229–2237, 2014.

[27]
J. Ding, Q. Zhang, S. Hu, Q. Wang, and Q. Ye, “Clusters partition and zonal
voltage regulation for distribution networks with high penetration of pvs,”
IET Generation, Transmission & Distribution, vol. 12, no. 22, pp.
6041–6051, 2018.

[28]
S. S. Torbaghan, G. Suryanarayana, H. Höschle, R. D'hulst, F. Geth,
C. Caerts, and D. Van Hertem, “Optimal flexibility dispatch problem using
second-order cone relaxation of ac power flows,” IEEE Transactions on
Power Systems, vol. 35, no. 1, pp. 98–108, 2019.

[29]
A. Ben-Tal and A. Nemirovski, “Robust optimization–methodology and
applications,” Mathematical programming, vol. 92, no. 3, pp.
453–480, 2002.

[30]
E. E. Elattar, A. M. Shaheen, A. M. Elsayed, and R. A. El-Sehiemy, “Optimal
power flow with emerged technologies of voltage source converter stations in
meshed power systems,” IEEE Access, vol. 8, pp. 166 963–166 979,

[31]
Y. Li and Y. Li, “Two-step many-objective optimal power flow based on knee
point-driven evolutionary algorithm,” Processes, vol. 6, no. 12, p.
250, 2018.

[32]
K. Knezović, A. Soroudi, A. Keane, and M. Marinelli, “Robust
multi-objective pq scheduling for electric vehicles in flexible unbalanced
distribution grids,” IET Generation, Transmission & Distribution,
vol. 11, no. 16, pp. 4031–4040, 2017.

[33]
A. Chakraborty, “Generating multivariate correlated samples,”
Computational Statistics, vol. 21, no. 1, pp. 103–119, 2006.

[34]
C. KI Williams, Gaussian processes for machine learning.1em plus
0.5em minus 0.4emTaylor & Francis Group, 2006.

[35]
J. Orduz, “Sampling from a multivariate normal distribution,” blog, March,
2019 [Online]. [Online]. Available: <https://tinyurl.com/bcchm6c9>

[36]
T. Muschinski, G. J. Mayr, T. Simon, N. Umlauf, and A. Zeileis,
“Cholesky-based multivariate gaussian regression,” Econometrics and
Statistics, 2022.

[37]
A. Koirala, L. Suárez-Ramón, B. Mohamed, and P. Arboleya,
“Non-synthetic european low voltage test system,” International
Journal of Electrical Power & Energy Systems, vol. 118, p. 105712, 2020.

[38]
V. Rigoni, L. F. Ochoa, G. Chicco, A. Navarro-Espinosa, and T. Gozel,
“Representative residential lv feeders: A case study for the north west of
england,” IEEE Transactions on Power Systems, vol. 31, no. 1, pp.
348–360, 2015.

[39]
E. Cotilla-Sanchez, P. D. Hines, C. Barrows, S. Blumsack, and M. Patel,
“Multi-attribute partitioning of power networks based on electrical
distance,” IEEE Transactions on Power Systems, vol. 28, no. 4, pp.
4979–4987, 2013.

[40]
J. P. Hespanha, “An efficient matlab algorithm for graph partitioning,”
Santa Barbara, CA, USA: University of California, 2004.

[41]
B. Mourad, “On a spectral property of doubly stochastic matrices and its
application to their inverse eigenvalue problem,” Linear algebra and
its applications, vol. 436, no. 9, pp. 3400–3412, 2012.

[42]
F. Scarlatache, G. Grigoraş, G. Chicco, and G. Cârţină,
“Using k-means clustering method in determination of the optimal placement
of distributed generation sources in electrical distribution systems,” in
2012 13th International Conference on Optimization of Electrical and
Electronic Equipment (OPTIM).1em plus 0.5em minus 0.4emIEEE,
2012, pp. 953–958.

[43]
P. J. Rousseeuw, “Silhouettes: a graphical aid to the interpretation and
validation of cluster analysis,” Journal of computational and applied
mathematics, vol. 20, pp. 53–65, 1987.

[44]
H. Tang and S. Wang, “Energy flexibility quantification of grid-responsive
buildings: Energy flexibility index and assessment of their effectiveness for
applications,” Energy, vol. 221, p. 119756, 2021.

[45]
A. Ulbig and G. Andersson, “Analyzing operational flexibility of electric
power systems,” International Journal of Electrical Power & Energy
Systems, vol. 72, pp. 155–164, 2015.

[46]
M. U. Hashmi, A. Mukhopadhyay, A. Bušić, and J. Elias, “Storage
optimal control under net metering policies,” arXiv preprint
arXiv:2002.01524, 2020.

[47]
Y. Chen, M. U. Hashmi, J. Mathias, A. Bušić, and S. Meyn,
“Distributed control design for balancing the grid using flexible loads,”
in Energy Markets and Responsive Grids.1em plus 0.5em minus
0.4emSpringer, 2018, pp. 383–411.

[48]
A. Bušić, M. U. Hashmi, and S. Meyn, “Distributed control of a fleet
of batteries,” in 2017 American Control Conference (ACC).1em
plus 0.5em minus 0.4emIEEE, 2017, pp. 3406–3411.

[49]
K. Deb and H. Gupta, “Searching for robust pareto-optimal solutions in
multi-objective optimization,” in International conference on
evolutionary multi-criterion optimization.1em plus 0.5em minus
0.4emSpringer, 2005, pp. 150–164.

[50]
M. U. Hashmi, A. Koirala, H. Ergun, and D. Van Hertem, “Flexible and
curtailable resource activation in a distribution network using nodal
sensitivities,” in 2021 International Conference on Smart Energy
Systems and Technologies (SEST).1em plus 0.5em minus 0.4emIEEE, 2021, pp. 1–6.

[51]
“Network and load data Github,” Online,

[52]
C. Coffrin, R. Bent, K. Sundar, Y. Ng, and M. Lubin, “Powermodels. jl: An
open-source framework for exploring power flow formulations,” in 2018
Power Systems Computation Conference (PSCC).1em plus 0.5em minus
0.4emIEEE, 2018, pp. 1–8.

[53]
A. Wächter and L. T. Biegler, “On the implementation of an interior-point
filter line-search algorithm for large-scale nonlinear programming,”
Mathematical programming, vol. 106, no. 1, pp. 25–57, 2006.

[54]
J. Miettinen, H. Holttinen, and B.-M. Hodge, “Simulating wind power forecast
error distributions for spatially aggregated wind power plants,” Wind
Energy, vol. 23, no. 1, pp. 45–62, 2020.

[55]
M. U. Hashmi, W. Labidi, A. Bušić, S.-E. Elayoubi, and T. Chahed,
“Long-term revenue estimation for battery performing arbitrage and ancillary
services,” in IEEE Int. Conference on Communications, Control, and
Computing Technologies for Smart Grids (SmartGridComm).1em plus 0.5em
minus 0.4emIEEE, 2018, pp. 1–7.

§ ZONE SELECTION ALGORITHM

The silhouette coefficient of node $j \in M_x$ is calculated as,
s_j = \frac{b_j - a_j}{\max \{b_j, a_j\}},
where $a_i$ denotes the mean distance between node $j$ and all other nodes in the group $M_x$ and $b_j$ denotes the mean distance between node $j$ and all other nodes not in the group $M_x$.
The mean silhouette score of a group $M_x$ is given as
S_x = \frac{1}{|M_x|} \sum_{i=1}^{|M_x|} s_i,
where $|M_x|$ denotes the number of nodes in the group $M_x$.
The mean silhouette score of a $k$ partitioned network is given as
\text{SC}_T^k = \frac{1}{k} \sum_{x=1}^{k} S_x.
For a given network, vary $k$ to maximize the value of $\text{SC}_T^k$. Note $\text{SC}_T^k$ may be high for a small number of clusters which may not suit the application, therefore, the selection of best $k$ which increases $\text{SC}_T^k$, is a design problem.

The algorithm for selecting the best number of zones that maximizes the mean silhouette score is detailed below.
The output of the algorithm is the ideal number of zones for our flexibility assessment application.

Inputs: Network details, $\text{SC}_T^k=[~]$,
[1]
Calculate admittance matrix, $Y$, for the network,
Make the electrical distance matrix double stochastic,
Set value of $k=2$ clusters.
Calculate $k$ the largest eigenvalues and eigenvectors.
Use $k$-means clustering and calculate $\text{SC}_T^k$ and concatenate,
Increment $k$ till $k\leq N$ and Goto Step 4,
Analyse $\text{SC}_T^k$ vector to select the best suited number of zones. The best suited number of zones depends on (a) how many clusters will meet the required application for which clusters are formed, and (b) silhouette score of the cluster. For example, if the silhouette score is very high for very few numbers of clusters, but it does not meet our purpose, then we may opt for a lower silhouette score. .

§ NETWORK DIAGRAM

The network diagram and the zones are shown in Fig. <ref>.

Network diagram with 12 zones indicated