# Identifying Secure Operating Ranges for
DER Control using Bilevel Optimization

Kshitij Girigoudar, and Line A. Roald The authors are with the Department of
Electrical and Computer Engineering, University of Wisconsin-Madison, USA
(e-mail<EMAIL_ADDRESS>roald@wisc.edu).

###### Abstract

Active distribution grids are accommodating an increasing number of
controllable electric loads and distributed energy resources (DERs). A
majority of these DERs are managed by entities other than the distribution
utility, such as individual customers or third-party aggregators, who control
the loads and DERs without consideration of any distribution grid constraints.
This makes it challenging for a distribution system operator (DSO) to allow
third-party aggregators and transmission operators to fully exploit the
flexibility offered by these resources while also ensuring that distribution
grid constraints such as voltage magnitude limits are not violated. In this
paper, we develop a bilevel optimization-based framework to determine the
aggregate power flexibility that can be obtained from an unbalanced
distribution grid while ensuring that there is no disaggregation solution that
leads to grid constraint violations. The results are a set of constraints and
operating rules that are easy to communicate, and which provide the entities
that procure flexibility from DERs (e.g. transmission operators or third-party
aggregators) with the ability to freely implement their own disaggregation
strategy without intervention from the DSO. The proposed approach is tested on
two unbalanced distribution feeders and our simulation results indicate that
it is possible to determine a wide range of aggregate power flexibility, as
long as a simple set of rules for DER control activation are followed.

###### Index Terms:

Aggregate power flexibility, bilevel optimization, distribution grids, optimal
power flow, solar PV, strong duality.

## Nomenclature

Sets & Indices

$\mathcal{N}$

Set of nodes excluding slack node

${\Phi}$

Set of phases, ${\Phi}=\\{a,b,c\\}$

$i$

Node index, $i\in\mathcal{N}$

$\phi$

Phase index, $\phi\in\Phi$

Parameters

$n$

Number of single-phase nodes, $n=3|\mathcal{N}|$

${v}_{\text{d}i}^{\phi},v_{\text{q}i}^{\phi}$

Rectangular components of current voltage phasor with magnitude
$|v_{i}^{\phi}|$

$\underline{v},\overline{v}$

Lower and upper voltage magnitude limits

${V}_{\text{d}i},{V}_{\text{q}i}\in\mathbb{R}^{3}$

Rectangular components of current three-phase voltage phasor $V_{i}$

${p}_{i}^{\phi},q_{i}^{\phi}$

Active and reactive components of current power injection

${P}_{i},Q_{i}\in\mathbb{R}^{3}$

Active and reactive components of current three-phase power injection

${p}_{\text{G},i}^{\phi},q_{\text{G},i}^{\phi}$

Active and reactive components of current power generation

${p}_{\text{L},i}^{\phi},q_{\text{L},i}^{\phi}$

Active and reactive components of current power load demand

$\underline{p}_{\text{G},i}^{\phi},\overline{p}_{\text{G},i}^{\phi}$

Lower and upper limits for active power generation

$\underline{p}_{\text{L},i}^{\phi},\overline{p}_{\text{L},i}^{\phi}$

Lower and upper limits for active power load demand

$\Delta\underline{p},\Delta\overline{p}$

Lower and upper limits for maximum available active power flexibility

${pf}_{\text{L},i}^{\phi}$

Load constant power factor setting

$|{s}_{\text{G},i}^{\phi}|$

Solar PV inverter apparent power capacity

${\gamma}_{\text{G},i}^{\phi}$

Solar PV inverter constant power ratio setting

${P}_{G},P_{L}\in\mathbb{R}^{n}$

Active power generation and load vectors at current operating point

Variables

$\Delta\bm{p}^{-},\Delta\bm{p}^{+}$

Lower and upper limit variables for aggregate active power flexibility

$\Delta\bm{p}_{\mathbf{G},i}^{\phi},\Delta\bm{p}_{\mathbf{L},i}^{\phi}$

Active power generation and load flexibility variables

$\Delta\bm{P}_{\mathbf{G}},\\!\Delta\bm{P}_{\mathbf{L}}\\!\in\\!\mathbb{R}^{n}$

Vector of active power generation and load flexibility variables

$\bm{q}_{\mathbf{G},i}^{\phi},\bm{q}_{\mathbf{L},i}^{\phi}$

Reactive power generation and load variables

$\bm{Q}_{\mathbf{G}},\bm{Q}_{\mathbf{L}}\in\mathbb{R}^{n}$

Vectors of reactive power generation and load variables

$\bm{P},\bm{Q}\in\mathbb{R}^{n}$

Vectors of active and reactive power injection variables

$\bm{v}_{\mathbf{d}i}^{\phi},\bm{v}_{\mathbf{q}i}^{\phi}$

Rectangular components of voltage variable with magnitude
$|\bm{v}_{i}^{\phi}|$

$\bm{V}_{\mathbf{d}i},\bm{V}_{\mathbf{q}i}\in\mathbb{R}^{3}$

Rectangular components of three-phase voltage variable

$\bm{V}_{\mathbf{d}},\bm{V}_{\mathbf{q}}\in\mathbb{R}^{n}$

Vectors of voltage variables in rectangular form

$\bm{\gamma}_{\mathbf{G},i}^{\phi}$

Solar PV inverter power ratio variable

$\overline{\bm{q}}_{\mathbf{G},i}^{\phi}$

Solar PV inverter maximum available reactive power variable

$\bm{X}_{\textbf{u}},\bm{X}_{\textbf{l}}$

Vectors of upper- and lower-level variables of bilevel problem

##  I Introduction

Today’s distribution grids are experiencing an increasing amount of smart
electric loads and distributed energy resources (DERs), such as electric
vehicles, smart home appliances, rooftop solar photovoltaic (PV) systems and
energy storage. With growing penetration of dispatchable DERs, distribution
grids have evolved from being considered as passive loads to take on a more
active role in system operation [1]. Importantly, dispatchable DERs owned by
individual customers can provide flexibility in the form of energy services
such as operating reserves to transmission systems. However, each customer-
owned DER can only contribute a comparatively small amount of flexibility. To
achieve a high enough level of flexibility to act as a reserve, aggregation
and coordinated control of many DERs is required. Such aggregation and control
services are typically provided by DER aggregators. In this paper, we refer to
them as _third-party aggregators_ to highlight that they are not part of a
distribution utility.

Typically, a distribution feeder comprises of DERs managed by various entities
whose control objectives can be different from each other. The flexibility
offered by DERs relies on assumptions about who controls them. DERs managed
directly by their owners, who choose not to enroll their DERs with a third-
party aggregator, may act as uncontrollable loads and may increase variability
rather than offering any flexibility. On the one hand, DERs can be controlled
individually by the distribution system operator (DSO), who can use them to
manage distribution grid constraints. Furthermore, third-party aggregators can
control multiple DERs to achieve an aggregated response across a large
population of DERs, typically without awareness of their exact geographical
location and without any consideration of the distribution grid constraints.
To comply with the recently announced FERC Order 2222 [2] in the USA,
independent system operators (ISOs) will be required to allow third-party DER
aggregators to participate in electricity markets. This may lead to a
significant increase in the number of aggregators and the overall potential
for DER flexibility. The focus of this paper is to propose a method that can
be used to determine limits on the combined DER flexibility that can be
offered from within a distribution feeder without violating any distribution
grid constraints.

There are numerous works in literature determining the flexibility of
individual DERs [3, 4, 5] or aggregate power flexibility of distribution grids
[6, 7, 8, 9, 10, 11]. Most of these methods assume that network constraints
such as voltage magnitude limits are non-binding, and hence do not consider
them. Approaches that do consider network constraints include geometric
methods, where the aggregate power flexibility is determined using polytopic
projection or calculating the Minkowski sum of the polytope sets, as discussed
for balanced distribution grids in [12, 13, 14] and for unbalanced three-phase
network models in [15, 16, 17, 18]. However, when assessing the total
available range of flexibility, these works assume that the DSO directly
controls all DERs, i.e. the DSO acts as an aggregator and can leverage DER
resources to both provide transmission grid services and manage distribution
grid constraints. Other works model the interactions between aggregators and
DSOs [19, 20, 21, 22], or between end consumers and aggregators [22, 23],
including methods that provide DSOs with the opportunity to change control
signals from the aggregator that would lead to violations[24].

Existing methods rely on strong assumptions on the frequency of communication
between the DSO, ISO and aggregators or either require or assume an allocation
of flexibility among DERs (i.e. which DERs are going to provide the desired
flexibility). In a practical setting, where multiple aggregators and DER
owners may be competing for the same limited grid capacity, it is problematic
that the methods employed to determine the safe amount of DER flexibility may
create limits that favor certain DERs over others. This is particularly when
the DSO does not have any information about the cost of contracting
flexibility from individual DERs. To address this issue, we propose a
different approach that focuses on assessing a range of aggregate DER
flexibility, while making more practical assumptions regarding communication
between the DSO, aggregator and ISO. We propose a hierarchical decision making
process where an upper-level authority (i.e. DSO) takes decisions based on the
unknown, and thus assumed to be worst-case, responses by a lower-level entity
(i.e. aggregators). This framework naturally lends itself to modelling as a
bilevel optimization problem. A key feature of the proposed approach is that
it does not prescribe or limit which DERs need to be providing flexibility, as
long as the cumulative amount of actuation stays within the allowable range of
aggregate DER flexibility. Thus, the method places minimal limits on both how
much flexibility should be procured from each of the aggregators active within
the feeder (allowing ISOs to procure flexibility from the aggregator that has
the most competitive offer), and does not limit how aggregators choose to
actuate the DERs that they control to achieve a certain set-point after they
have been contracted to provide flexibility.

The main contributions of the paper are as follows: 1) We propose a bilevel
formulation to determine aggregate DER flexibility range that can be offered
by distribution grids which can guarantee that _any_ disaggregation solution
(including the worst-case condition) is feasible, thus enabling more safe
disaggregation without further coordination with the DSO. 2) We present a
computationally tractable, strong duality based reformulation of the bilevel
problem which can be solved efficiently using an iterative approach. 3) The
proposed method is tested on two unbalanced distribution feeders, where we
investigate how the range of the feeder’s aggregate power flexibility can be
maximized by enforcing a set of simple rules for DER actuation.

The rest of the paper is organized as follows: Section II describes the
problem setup. Section III presents the bilevel problem formulation. Section
IV details the analytical reformulation of the problem and an iterative
solution approach to solve the problem efficiently. We present the numerical
results on two test cases in Section V. Section VI concludes.

##  II Problem Setup

We next introduce the problem setup by providing a motivating example, and
giving an overview of the bilevel optimization framework.

###  II-A Motivating Example

To comply with the FERC Order 2222 [2], ISOs in the USA will be required to
allow DER aggregators to participate in electricity markets. While the
clearing of electricity markets by ISOs currently only considers transmission
grid constraints, increased levels of generation and reserve provision from
DERs located at the distribution grid level may require that ISOs include
distribution grid constraints in their optimal dispatch algorithms. Here, we
present a method that enables DSOs to identify limits on how much DER
flexibility the ISO can procure from a given feeder without causing constraint
violations. We demonstrate the intended use of this method with a small,
illustrative example.

Consider the grid depicted in Fig. 1 where generators and a distribution grid
are connected to the transmission system. On the transmission side of the
grid, two generators G1 and G2 operated by separate companies bid into the
electricity market. The amount of power provided by each of them is determined
by the ISO by solving a transmission system dispatch problem. This problem
considers the transmission line constraints defined by the minimum and maximum
power transfer capacity $p^{-},p^{+}$, which bounds the total amount of power
that the ISO can get from the two generators combined. This example is a
simple illustration of the current state of the art in ISO market clearing.

Now, consider the distribution side of the grid where two DER aggregators A1
and A2 control DERs at multiple nodes, as illustrated by the red and green
circles. These aggregators can now bid the flexibility offered by the DERs
into the electricity market. However, including the operational limits of each
individual DER (analogous to individual generator constraints) along with a
detailed model of the distribution grid (analogous to the transmission grid
model) into the optimization problem would significantly increase the
computational complexity [25, 26, 27]. As an alternative, we propose a method
to enable the DSO to identify aggregate active power flexibility limits
$\Delta p^{-},\Delta p^{+}$ around the current operating point. These limits
represent the amount of active power flexibility that can be provided from the
combined set of DERs (controlled both by A1 and A2) without causing internal
distribution grid constraints to be violated. By integrating these limits in
their dispatch algorithms, the ISO can determine a secure amount of power
flexibility to procure from the aggregators.

Figure 1: Example of transmission-distribution interaction.

###  II-B Bilevel Modeling Framework

To enable the DSO to identify a secure aggregate power flexibility range, we
propose a bilevel optimization approach.

####  II-B1 Bilevel problem overview

The upper-level (or leader) problem represents the decision making of the DSO,
and maximizes the aggregate change in DER power injections relative to the
current operating point while guaranteeing that the worst-case voltage
magnitudes (obtained by solving a set of lower level problems) stay within
acceptable limits. The decision variables $\bm{X}_{\mathbf{u}}$ of this upper-
level problem are the control actions determined by the DSO. Importantly, this
set of decision variables include $\Delta\bm{p}^{+}$ and $\Delta\bm{p}^{-}$,
the upper and lower limits on the aggregate active power flexibility, that can
safely be provided from the feeder. Each lower-level (or follower) problem
solves a three-phase optimal power flow (OPF) problem for a given node
$i\in\mathcal{N}$ in the network with objective to find the worst-case
(largest or smallest) voltage magnitude achievable at that node. Here,
$\mathcal{N}$ represent the set of three-phase nodes in the feeder (excluding
the substation, which acts as a slack node) and $n=3\cdot|\mathcal{N}|$ is the
total number of single-phase nodes in the network. The lower-level problem
formulation is similar to the model proposed in [28], which we extend to
include the unbalanced three-phase power flow model derived in [29]. The
decision variables $\bm{X}_{\mathbf{l}}$ of the lower-level problem are
assumed to be independently determined by the aggregator(s) and may therefore
be equal to the values that would lead to worst case conditions in the grid.
Note that there are many lower level problems. Specifically, given a network
with $n$ single-phase nodes, we need $2n$ lower-level problems, i.e. $n$
problems identifying smallest achievable voltage magnitude and $n$ problems
finding largest achievable voltage magnitude for each single-phase node.

####  II-B2 Input data requirements

To solve the bilevel problem, we require information regarding the
distribution grid parameters as well as information about the current
operating state. Furthermore, the location and feasible range of set-points
for individual DERs is assumed to be known. Note that this DER information
does not have to be real-time information, but could simply be device
characteristics provided by the aggregators when new DERs enroll in their
programs.

####  II-B3 Rules on DER activation

In order to obtain practical ranges of flexibility which are not overly
conservative, we enforce a simple rule on DER control actuation. This rule
states that the active power injection of all DERs in the feeder have to
either be increased or decreased (i.e., we cannot have a scenario where some
DERs reduce their power injections and others increase them). This rule is
needed to ensure that limiting the aggregate active power response is
sufficient to limit voltage magnitude violations in the feeder. Without it,
DERs in one part of the feeder can arbitrarily increase their injections as
long as the increases are offset by an equally large decrease in injections
elsewhere. If this is allowed, even a zero change in the overall active power
injection at the substation can cause voltage violations. To integrate this
rule in our optimization problem, we consider two cases:

  1. (a)

_Positive (+) case:_ The positive case assumes that the active power
injections at any node $i\in\mathcal{N}$ connected to phase
$\phi\in\\{a,b,c\\}$ is controlled to a value that is greater than the current
active power injection. For loads, this implies that demand is reduced while
for DERs such as solar PV inverters, it implies that active power generation
is increased. Thus, we include the following constraints,

$\displaystyle\Delta\bm{p}_{\mathbf{L},i}^{\phi}\leq
0,~{}~{}\Delta\bm{p}_{\mathbf{G},i}^{\phi}\geq
0,~{}\forall_{\phi\in\\{a,b,c\\},i\in\mathcal{N}}.$ (1)

Here $\Delta\bm{p}_{\mathbf{L},i}^{\phi},\Delta\bm{p}_{\mathbf{G},i}^{\phi}$
represent the deviation of the active power of the individual loads and solar
PV inverters at node $i\in\mathcal{N}$ connected to phase
$\phi\in\\{a,b,c\\}$, respectively, from the current operating point.

  2. (b)

_Negative (-) case:_ This negative case assumes that the active power
injections at every node $i\in\mathcal{N}$ connected to phase $\phi$ is
controlled to a value that is lower than the current operating point. This
gives rise to the following constraints for flexibility activation,

$\displaystyle\Delta\bm{p}_{\mathbf{L},i}^{\phi}\geq
0,~{}~{}\Delta\bm{p}_{\mathbf{G},i}^{\phi}\leq 0,~{}$
$\displaystyle\forall_{\phi\in\\{a,b,c\\},i\in\mathcal{N}}.$ (2)

Enforcing this rule increases the number of lower-level problems we need to
solve, as the worst-case conditions determined for the positive and negative
cases will be different. Given a network with $n$ single-phase nodes, we now
need $4n$ lower-level problems, i.e. $2n$ each for the positive and negative
case.

####  II-B4 DER reactive power control

We further assume that all DERs are equipped with smart inverters which can
provide reactive power support as per IEEE Standard 1547-2018 [30] and that
the DSO is allowed to determine the reactive power setpoints for these
inverters. By making this assumption, it is possible to significantly increase
the amount of aggregate power flexibility that DERs are allowed to provide, as
demonstrated by our case study results in Section V. In this work, we consider
three reactive power modes (i.e. constant power factor mode, constant reactive
power mode and voltage-reactive power mode) in which the inverters can
operate. The detailed modelling of the inverter modes is discussed later in
Section III-D.

##  III Bilevel Problem Formulation

In this section, we describe our model of the distribution grid and the
flexibility offered by DERs such as solar PV inverters and controllable loads,
before presenting the full bilevel problem formulation.

###  III-A Notation

We consider a network where the distribution substation with index $i=0$ is
chosen as the slack node. $\mathcal{N}$ denotes the set of non-slack nodes.
Without loss of generality, we assume all nodes have three phases with the set
of phases defined by $\Phi=\\{a,b,c\\}$, and total number of single-phase
nodes is $n=3\cdot|\mathcal{N}|$. Apart from the substation, the only other
generators in the network are single-phase, residential solar PV systems. To
simplify notation, we assume that there is one solar PV inverter and one load
at each single-phase node. If there is no solar PV or load at some node, the
corresponding entries are set to zero. All the optimization variables are
denoted using bold symbols and all vectors as well as matrices are represented
using capital letters. Given a vector $X\in\mathbb{C}^{n}$, the matrix
${dg}(X)\in\mathbb{C}^{(n\times n)}$ is a diagonal matrix with elements of $X$
on its diagonal. We use $\odot$ to denote the element-wise product of two
vectors. Given a complex phasor $x$, the complex conjugate is represented by
$x^{*}$.

###  III-B Current Operating Point

At any three-phase node $i\in\mathcal{N}$, we assume that we have access to
the current voltage phasor ${V}_{i}=V_{\text{d}i}+jV_{\text{q}i}$, with
$V_{\text{d}i}=[v_{\text{d}i}^{a}~{}~{}v_{\text{d}i}^{b}~{}~{}v_{\text{d}i}^{c}]^{\top}$
and
$V_{\text{q}i}=[v_{\text{q}i}^{a}~{}~{}v_{\text{q}i}^{b}~{}~{}v_{\text{q}i}^{c}]^{\top}$
representing the real and imaginary components, respectively. Similarly, we
can define the three-phase active and reactive power injections by
${P}_{i}=[p_{i}^{a}~{}~{}p_{i}^{b}~{}~{}p_{i}^{c}]^{\top}$ and
${Q}_{i}=[q_{i}^{a}~{}~{}q_{i}^{b}~{}~{}q_{i}^{c}]^{\top}$, respectively. The
power injections are expressed as the difference between the generation and
load demand using

$\displaystyle{p}_{i}^{\phi}={p}_{\text{G},i}^{\phi}-{p}_{\text{L},i}^{\phi},~{}$
$\displaystyle\forall_{\phi\in\Phi,i\in\mathcal{N}},$ (3a)
$\displaystyle{q}_{i}^{\phi}={q}_{\text{G},i}^{\phi}-{q}_{\text{L},i}^{\phi},~{}$
$\displaystyle\forall_{\phi\in\Phi,i\in\mathcal{N}},$ (3b)

where ${p}_{\text{G},i}^{\phi},{q}_{\text{G},i}^{\phi}$ are the active and
reactive power generation by solar PV inverters;
${p}_{\text{L},i}^{\phi},{q}_{\text{L},i}^{\phi}$ denote the respective active
and reactive load demand. The power balance is maintained by the substation,
which supplies the difference between the load and power generation as well as
other losses in the grid. We denote the current active and reactive power
injection at the substation by
${P}_{\text{G},0},{Q}_{\text{G},0}\in\mathbb{R}^{3}$. Note that the assumption
of perfect knowledge of the current system state can be relaxed to include a
range of possible operating points.

###  III-C System Modeling in Lower-Level Problem

We next describe the model used in the lower-level problem to represent the
worst-case operating condition after DER control actuation.

####  III-C1 Voltage representation

The voltage variables in the bilevel problem are expressed in rectangular
form. For any node $i\in\\{0,\mathcal{N}\\}$, the real and imaginary
components of voltage are given by
$\bm{V}_{\mathbf{d}i}=[\bm{v}_{\mathbf{d}i}^{a}~{}~{}\bm{v}_{\mathbf{d}i}^{b}~{}~{}\bm{v}_{\mathbf{d}i}^{c}]^{\top}$
and
$\bm{V}_{\mathbf{q}i}=[\bm{v}_{\mathbf{q}i}^{a}~{}~{}\bm{v}_{\mathbf{q}i}^{b}~{}~{}\bm{v}_{\mathbf{q}i}^{c}]^{\top}$,
respectively. The distribution substation is assumed to be perfectly balanced,
i.e.

$\displaystyle\bm{V}_{\mathbf{d}0}+j\cdot\bm{V}_{\mathbf{q}0}$
$\displaystyle=\begin{bmatrix}e^{j0^{\circ}}&e^{-j120^{\circ}}&e^{j120^{\circ}}\end{bmatrix}^{\top}.$
(4)

For every other node $i\in\mathcal{N}$ with phase $\phi$, the voltage
magnitude variable is given by
$|\bm{v}_{i}^{\phi}|=\sqrt{(\bm{v}_{\mathbf{d}i}^{\phi})^{2}+(\bm{v}_{\mathbf{q}i}^{\phi})^{2}}$.
To avoid introducing a non-linear constraint for the voltage magnitude, we
introduce a linear approximation for this relationship. Given the current
voltage magnitude
$|{v}_{i}^{\phi}|=\sqrt{{(v_{\text{d}i}^{\phi})}^{2}+{(v_{\text{q}i}^{\phi})}^{2}}$,
we use a First-order Taylor approximation to get the following linear
equation:

$\displaystyle{(v_{\text{d}i}^{\phi})}^{2}\\!+\\!{(v_{\text{q}i}^{\phi})}^{2}\\!+\\!2v_{\text{d}i}^{\phi}\bm{v}_{\mathbf{d}i}^{\phi}\\!+\\!2v_{\text{q}i}^{\phi}\bm{v}_{\mathbf{q}i}^{\phi}={{|{v}_{i}^{\phi}|}^{2}\\!+\\!2|{v}_{i}^{\phi}||\bm{v}_{i}^{\phi}|}.$
(5)

Note that we do not enforce voltage magnitude constraints in the lower level
problem, as these are accounted for in the upper level problem, as further
discussed in Section III-E.

####  III-C2 Load modeling

The total demand from load at node $i$ in phase $\phi$ is constrained by the
lower and upper limits
$\underline{p}_{\text{L},i}^{\phi},\overline{p}_{\text{L},i}^{\phi}$ for the
active power demand from this node, i.e.

$\displaystyle\underline{p}_{\text{L},i}^{\phi}\leq{p}_{\text{L},i}^{\phi}+\Delta\bm{p}_{\mathbf{L},i}^{\phi}\leq\overline{p}_{\text{L},i}^{\phi},~{}\forall_{\phi\in\Phi,i\in\mathcal{N}}.$
(6)

Assuming that loads operate with a constant power factor
$pf_{\text{L},i}^{\phi}$, the reactive power demand is given by

$\displaystyle\bm{q}_{\mathbf{L},i}^{\phi}={\frac{\sqrt{1-{(pf_{\text{L},i}^{\phi})}^{2}}}{{pf_{\text{L},i}^{\phi}}}}\cdot\left({p}_{\text{L},i}^{\phi}+\Delta\bm{p}_{\mathbf{L},i}^{\phi}\right),~{}\forall_{\phi\in\Phi,i\in\mathcal{N}}.$
(7)

####  III-C3 Solar PV active power modeling

The constraint on solar PV generation at node $i$ in phase $\phi$ is given by

$\displaystyle\underline{p}_{\text{G},i}^{\phi}\leq{p}_{\text{G},i}^{\phi}+\Delta\bm{p}_{\mathbf{G},i}^{\phi}\leq\overline{p}_{\text{G},i}^{\phi},~{}\forall_{\phi\in\Phi,i\in\mathcal{N}},$
(8)

where $\underline{p}_{\text{G},i}^{\phi},\overline{p}_{\text{G},i}^{\phi}$
denote the respective lower and upper limits for the active power generation
of the inverter. The reactive power modeling of solar PV inverters is
discussed later.

####  III-C4 Aggregate active power flexibility

The total active power flexibility $\Delta\bm{p}$ provided by the distribution
grid is given by the sum of the flexibility from individual loads and solar PV
inverters, and is bounded by lower and upper limits of the aggregate power
flexibility provided from the feeder $\Delta\bm{p}^{-},\Delta\bm{p}^{+}$. This
gives rise to the following constraint,

$\displaystyle\Delta\underline{p}\\!\leq\\!\Delta\bm{p}^{-}\\!\leq\\!\underbrace{\sum_{i\in\mathcal{N}}\sum_{\phi\in\Phi}\left(\Delta\bm{p}_{\mathbf{G},i}^{\phi}-\Delta\bm{p}_{\mathbf{L},i}^{\phi}\right)}_{\Delta\bm{p}}\\!\leq\\!\Delta\bm{p}^{+}\\!\leq\\!\Delta\overline{p}.$
(9)

Here,
$\Delta\overline{p}={\sum_{i\in\mathcal{N}}\sum_{\phi\in\Phi}(\overline{p}_{\text{G},i}^{\phi}-{p}_{\text{G},i}^{\phi}+\overline{p}_{\text{L},i}^{\phi}-{p}_{\text{L},i}^{\phi})}$
is the maximum available aggregate flexibility that can be offered by the
network, and
$\Delta\underline{p}={\sum_{i\in\mathcal{N}}\sum_{\phi\in\Phi}(\underline{p}_{\text{G},i}^{\phi}-{p}_{\text{G},i}^{\phi}+\underline{p}_{\text{L},i}^{\phi}-{p}_{\text{L},i}^{\phi})}$
is the minimum available aggregate flexibility. In each lower level problem,
the constraint (9) is combined with the condition that the active power
injections either increase (1) or decrease (2). Note that setting
$\Delta\bm{p}=0$ forces
$\Delta\bm{p}_{\mathbf{L},i}^{\phi}=\Delta\bm{p}_{\mathbf{G},i}^{\phi}=0$ and
recovers the current operating point.

####  III-C5 Power Flow

We use a fixed-point power flow interpretation in this paper to linearize the
power flow equations because it exhibits better global approximation accuracy
when compared to other linear models such as first-order Taylor approximation
[31]. With this approximation, the vector of voltage variables
$\bm{V}_{\mathbf{d}}=[\bm{V}_{\mathbf{d}i\in\mathcal{N}}^{\top}]^{\top},\bm{V}_{\mathbf{q}}=[\bm{V}_{\mathbf{q}i\in\mathcal{N}}^{\top}]^{\top}$
is expressed as

$\displaystyle\bm{V}_{\mathbf{d}}=\real\\{{Z}_{1}\Big{\\}}+\real\\{Z_{2}\Big{\\}}\cdot\bm{P}+\imaginary\\{Z_{2}\Big{\\}}\cdot\bm{Q},$
(10a)
$\displaystyle\bm{V}_{\mathbf{q}}=\imaginary\\{{Z}_{1}\Big{\\}}+\imaginary\\{Z_{2}\Big{\\}}\cdot\bm{P}-\real\\{Z_{2}\Big{\\}}\cdot\bm{Q}.$
(10b)

Here, $Z_{1}\in\mathbb{R}^{n}$ and $Z_{2}\in\mathbb{R}^{n\times n}$ are fixed
and calculated using the current operating point as described in [29].
$\bm{P}\in\mathbb{R}^{n}$ includes the active power flexibility
$\Delta\bm{P}_{\mathbf{L}}=[\Delta\bm{P}_{\mathbf{L}i\in\mathcal{N}}^{\top}]^{\top},\Delta\bm{P}_{\mathbf{G}}=[\Delta\bm{P}_{\mathbf{G}i\in\mathcal{N}}^{\top}]^{\top}$
from loads and PV inverters, respectively, and $\bm{Q}\in\mathbb{R}^{n}$
includes the reactive power injections of PV inverters
$\bm{Q}_{\mathbf{G}}=[\bm{Q}_{\mathbf{G}i\in\mathcal{N}}^{\top}]^{\top}$ and
load reactive power demand
$\bm{Q}_{\mathbf{L}}=[\bm{Q}_{\mathbf{L}i\in\mathcal{N}}^{\top}]^{\top}$, as
defined by

$\displaystyle\bm{P}={P}_{\text{G}}+\Delta\bm{P}_{\mathbf{G}}-({P}_{\text{L}}+\Delta\bm{P}_{\mathbf{L}}),$
(11a) $\displaystyle\bm{Q}=\bm{Q}_{\mathbf{G}}-\bm{Q}_{\mathbf{L}},$ (11b)

where
$P_{\text{G}}=[P_{\text{G}i\in\mathcal{N}}^{\top}]^{\top},P_{\text{L}}=[P_{\text{L}i\in\mathcal{N}}^{\top}]^{\top}$
are the respective active power generation and load demand vectors at the
current operating point.

###  III-D Reactive Power Control from Solar PV Inverters

Following the IEEE Standard 1546-2018 [30], we consider three modes in which
the PV inverters can operate. We next describe some constraints that are
common to all inverters, before discussing the three reactive power control
modes. Note that the set-points provided by the DSO are different for each
reactive power control mode.

####  III-D1 Inverter Constraints

For all three modes, the solar PV inverter at node $i$ connected to phase
$\phi$ should not exceed the inverter apparent power capacity
$|{s_{\text{G},i}^{\phi}}|$. In its original form, the apparent power
constraint is quadratic. To avoid introducing non-linear constraints, we outer
approximate the circular feasible region using linear constraints [32]. The
resulting constraints are given by

$\displaystyle 0\leq
p_{\text{G},i}^{\phi}+\Delta\bm{p}_{\mathbf{G},i}^{\phi}\leq|{s_{\text{G},i}^{\phi}}|,$
$\displaystyle\forall_{\phi\in\Phi,i\in\mathcal{N}},$ (12a)
$\displaystyle-|{s_{\text{G},i}^{\phi}}|\leq\bm{q}_{\mathbf{G},i}^{\phi}\leq|{s_{\text{G},i}^{\phi}}|,$
$\displaystyle\forall_{\phi\in\Phi,i\in\mathcal{N}},$ (12b) $\displaystyle
p_{\text{G},i}^{\phi}+\Delta\bm{p}_{\mathbf{G},i}^{\phi}+\bm{q}_{\mathbf{G},i}^{\phi}\leq\sqrt{2}\cdot|{s_{\text{G},i}^{\phi}}|,$
$\displaystyle\forall_{\phi\in\Phi,i\in\mathcal{N}},$ (12c) $\displaystyle
p_{\text{G},i}^{\phi}+\Delta\bm{p}_{\mathbf{G},i}^{\phi}-\bm{q}_{\mathbf{G},i}^{\phi}\leq\sqrt{2}\cdot|{s_{\text{G},i}^{\phi}}|,$
$\displaystyle\forall_{\phi\in\Phi,i\in\mathcal{N}}.$ (12d)

####  III-D2 Constant power factor mode

This mode assumes that the DSO provides a target power factor setting
$pf_{\text{G},i}^{\phi}$ to the DERs. The reactive power injection
$\bm{q}_{\mathbf{G},i}^{\phi}$ is given by

$\displaystyle\bm{q}_{\mathbf{G},i}^{\phi}=\bm{\gamma}_{\mathbf{G},i}^{\phi}\cdot\left({p}_{\text{G},i}^{\phi}+\Delta\bm{p}_{\mathbf{G},i}^{\phi}\right),\hskip
47.00003pt\forall_{\phi\in\Phi,i\in\mathcal{N}},$ (13a)
$\displaystyle\\!{\frac{-\sqrt{1\\!-{(pf_{\text{G},i}^{\phi})}^{2}}}{{pf_{\text{G},i}^{\phi}}}}\\!\leq\\!\bm{\gamma}_{\mathbf{G},i}^{\phi}\\!\leq\\!{\frac{\sqrt{1\\!-{(pf_{\text{G},i}^{\phi})}^{2}}}{{pf_{\text{G},i}^{\phi}}}},~{}\forall_{\phi\in\Phi,i\in\mathcal{N}},$
(13b)

where $\bm{\gamma}_{\mathbf{G},i}^{\phi}$ is the power ratio of the inverter
and is a decision variable in the upper-level optimization problem. Note that
(13a) consist of bilinear terms involving the upper-level variables
$\bm{\gamma}_{\mathbf{G},i}^{\phi}$ and lower-level variables
$\Delta\bm{p}_{\mathbf{G},i}^{\phi}$. For a bilevel problem, the upper-level
variables are considered as parameters in the lower-level problem, and as a
result, (13) is linear for the lower-level problem.

####  III-D3 Constant reactive power mode

In this mode, the reactive power injection $\bm{q}_{\mathbf{G},i}^{\phi}$ is
specified by the DSO and hence, considered as an upper-level variable. In our
work, we assume that PV inverters are operating within some power factor range
[30] defined by the power ratio $\gamma_{\text{G},i}^{\phi}$ and we enforce
limits on the reactive power using

$\displaystyle\bm{q}_{\mathbf{G},i}^{\phi}\geq-\gamma_{\text{G},i}^{\phi}\\!\cdot\\!\left({p}_{\text{G},i}^{\phi}+\Delta\bm{p}_{\mathbf{G},i}^{\phi}\right),$
$\displaystyle\forall_{\phi\in\Phi,i\in\mathcal{N}},$ (14a)
$\displaystyle\bm{q}_{\mathbf{G},i}^{\phi}\leq\gamma_{\text{G},i}^{\phi}\\!\cdot\\!\left({p}_{\text{G},i}^{\phi}+\Delta\bm{p}_{\mathbf{G},i}^{\phi}\right),$
$\displaystyle\forall_{\phi\in\Phi,i\in\mathcal{N}}.$ (14b)

####  III-D4 Voltage-reactive power mode

In this mode, the reactive power is a function of the voltage magnitude, as
illustrated in Fig. 2 for a PV inverter at node $i$ connected to phase $\phi$.
The maximum available reactive power $\overline{\bm{q}}_{\mathbf{G},i}^{\phi}$
is specified by the DSO and hence, included in our bilevel task as an upper-
level decision variable. When operating outside the grey shaded region in Fig.
2, where the voltage magnitude is violating the specified limits
$\underline{v},\overline{v}$, the inverter either injects or absorbs the
maximum available reactive power $\overline{\bm{q}}_{\mathbf{G},i}^{\phi}$. In
the grey shaded region, the reactive power $\bm{q}_{\mathbf{G},i}^{\phi}$
depends on the voltage magnitude $|\bm{v}_{i}^{\phi}|$ and can be expressed as

$\displaystyle\bm{q}_{\mathbf{G},i}^{\phi}=\overline{\bm{q}}_{\mathbf{G},i}^{\phi}-2\overline{\bm{q}}_{\mathbf{G},i}^{\phi}\cdot\left(\frac{|\bm{v}_{i}^{\phi}|-\underline{v}}{\overline{v}-\underline{v}}\right),$
$\displaystyle\forall_{\phi\in\Phi,i\in\mathcal{N}},$ (15a) $\displaystyle
0\leq\overline{\bm{q}}_{\mathbf{G},i}^{\phi}\leq|s_{\text{G},i}^{\phi}|,$
$\displaystyle\forall_{\phi\in\Phi,i\in\mathcal{N}}.$ (15b)

Here, the actual upper bound of the maximum available reactive power
$\overline{\bm{q}}_{\mathbf{G},i}^{\phi}$ is given by
$\sqrt{|s_{\text{G},i}^{\phi}|^{2}-({p}_{\text{G},i}^{\phi}+\Delta\bm{p}_{\mathbf{G},i}^{\phi})^{2}}$.
Since this expression is nonlinear, we relax the upper bound as defined in
(15b).

We assume that the inverter is always operating in the gray shaded region
shown in Fig. 2. This assumption is reasonable since we enforce constraints on
the voltage magnitude $|\bm{v}_{i}^{\phi}|$ to be within the limits
$\underline{v},\overline{v}$ in the upper-level problem. Similar to the
constant power factor mode, the reactive power constraints in (15a) consist of
bilinear terms involving the upper-level variables
$\overline{\bm{q}}_{\mathbf{G},i}^{\phi}$ and lower-level variables
$|\bm{v}_{i}^{\phi}|$.

Figure 2: Reactive power-voltage magnitude characteristic.

###  III-E Bilevel Optimization Problem

The objective of the bilevel problem is to find the maximum range of aggregate
power flexibility while ensuring that the grid is secure even in the worst-
case conditions. In this section, we discuss the procedure to formulate the
bilevel problem when inverters are operating in constant power factor mode and
DSO controls the power ratio $\bm{\gamma}_{\mathbf{G},i}^{\phi}$ of each
inverter. The same procedure can be followed for bilevel problem formulations
with inverters operating in other reactive power modes. The upper-level
variable vector $\bm{X}_{\mathbf{u}}$ and lower-level variable vector
$\bm{X}_{\mathbf{l}}$ are denoted by

$\displaystyle\bm{X}_{\mathbf{u}}:=\\{\Delta\bm{p}^{+},\Delta\bm{p}^{-},\bm{\gamma}_{\mathbf{G},i}^{\phi}\forall_{\phi\in\Phi,i\in\mathcal{N}}\\},$
$\displaystyle\bm{X}_{\mathbf{l}}:=\left\\{\\{\bm{v}_{\mathbf{d}i}^{\phi},\bm{v}_{\mathbf{q}i}^{\phi},|{\bm{v}_{i}^{\phi}}|,\Delta\bm{p}_{\mathbf{G},i}^{\phi},\Delta\bm{p}_{\mathbf{L},i}^{\phi},\bm{q}_{\mathbf{L},i}^{\phi},\bm{q}_{\mathbf{G},i}^{\phi}\\}_{\forall_{\phi\in\Phi,i\in\mathcal{N}}}\right\\}.$

Since we need to represent both the positive and negative active power
activation cases and must also determine worst-case voltage magnitude (minimum
and maximum) achievable, we need to consider 4 lower level problems at every
single-phase node. Therefore, we use four sets of lower-level variables
corresponding to four scenarios summarized in Table I.

Each of the lower-level problems include system, inverter and active power
flexibility constraints. The set of system constraints is given by

$\displaystyle\mathbb{S}({\bm{X}}_{\mathbf{l}}):=\begin{cases}\text{Voltage
constraints}~{}\eqref{eq:Vlim_lin},\\\ \text{Load reactive power
constraints}~{}\eqref{eq:load_q},\\\ \text{Power
flow}~{}\eqref{eq:FP_rec}.\end{cases}$

The set of inverter constraints (for the constant power factor mode) is given
by

$\displaystyle\mathbb{I}(\bm{\gamma}_{\mathbf{G},i}^{\phi},{\bm{X}}_{\mathbf{l}}):=\text{Inverter
constraints}~{}\eqref{eq:inv_lim_lin},\eqref{eq:inv_cons_pf},$

while the active power flexibility constraints for either positive or negative
case are represented by

$\displaystyle\mathbb{F}(\Delta\bm{p}^{+},\Delta\bm{p}^{-},{\bm{X}}_{\mathbf{l}}):=\text{Positive
case}~{}\eqref{eq:P_pos},\eqref{eq:P_ss_constr}~{}\text{or}$
$\displaystyle\mathbb{F}(\Delta\bm{p}^{+},\Delta\bm{p}^{-},{\bm{X}}_{\mathbf{l}}):=\text{Negative
case}~{}\eqref{eq:P_neg},\eqref{eq:P_ss_constr}.$

TABLE I: Notation of lower-level problem variables corresponding to different scenarios Scenario | Active power | Worst case | Lower-level
---|---|---|---
(#) | activation case | voltage magnitude | variable set
1 | Positive | Minimum | ${\bm{X}}_{\mathbf{l},1}$
2 | Positive | Maximum | ${\bm{X}}_{\mathbf{l},2}$
3 | Negative | Minimum | ${\bm{X}}_{\mathbf{l},3}$
4 | Negative | Maximum | ${\bm{X}}_{\mathbf{l},4}$

The bilevel problem to determine the aggregate power flexibility limits
$\Delta\bm{p}^{+},\Delta\bm{p}^{-}$ can then be formulated as

$\displaystyle\max_{\bm{X}_{\mathbf{u}}}~{}$
$\displaystyle\Delta{\bm{p}^{+}}-\Delta{\bm{p}^{-}}$ ($\text{P}^{\pm}$) s.t.
$\displaystyle\underline{v}\leq|\bm{{v}}_{i}^{\phi}|_{1},|\bm{{v}}_{i}^{\phi}|_{2}\leq\overline{v},\hskip
45.00006pt\forall_{\phi\in\Phi,i\in\mathcal{N}},$
$\displaystyle\underline{v}\leq|\bm{{v}}_{i}^{\phi}|_{3},|\bm{{v}}_{i}^{\phi}|_{4}\leq\overline{v},\hskip
45.00006pt\forall_{\phi\in\Phi,i\in\mathcal{N}},$
$\displaystyle\text{where}~{}\forall{s\in\\{1,2,3,4\\}},$
$\displaystyle|\bm{{v}}_{i}^{\phi}|_{s}=\max_{{\bm{X}}_{\mathbf{l},s}}~{}(-1)^{s}\cdot|\bm{{{v}}}_{i}^{\phi}|,\hskip
45.00006pt\forall_{\phi\in\Phi,i\in\mathcal{N}},$ $\displaystyle\hskip
18.99995pt\text{s.t.
}\,\mathbb{S}({\bm{X}}_{\mathbf{l},s}),~{}\mathbb{I}(\bm{\gamma}_{\mathbf{G},i}^{\phi},{\bm{X}}_{\mathbf{l},s}),$
$\displaystyle\hskip
36.0001pt\mathbb{F}(\Delta\bm{p}^{+},\Delta\bm{p}^{-},{\bm{X}}_{\mathbf{l},s}),$

Similar problems can be formulated for other types of inverter control by
replacing
$\mathbb{I}(\bm{\gamma}_{\mathbf{G},i}^{\phi},{\bm{X}}_{\mathbf{l}})$ with
either (12), (14) for constant reactive power mode or (12), (15) for voltage-
reactive power mode.

##  IV Problem Reformulation and Solution Method

We next use the strong-duality theorem for the lower-level problem to reduce
the bilevel problem into a single-level problem. Given that the lower-level
problem has a finite optimal solution, the dual feasible set is non empty and
strong duality holds for every primal and dual feasible pairs.

###  IV-A Single-level Reformulation

To derive the strong duality based reformulation of the bilevel problem
$\text{P}^{\pm}$, we first rewrite the bilevel problem in a simplified and
more abstract form with one follower to get

$\displaystyle\max_{\Delta{\bm{p}},\bm{\gamma}}~{}$
$\displaystyle\Delta{\bm{p}}$ (16a) s.t.
$\displaystyle\bm{\hat{x}}_{\mathbf{l}}\leq b_{1},$ (16b)
$\displaystyle\text{where}~{}\bm{\hat{x}}_{\mathbf{l}}=\max_{\bm{x}_{\mathbf{l}}}~{}\bm{x}_{\mathbf{l}}$
(16c) $\displaystyle\hskip 59.00002pt\text{s.t. }a_{2}\bm{x}_{\mathbf{l}}\leq
b_{2},$ (16d) $\displaystyle\hskip
75.0001pta_{3}\bm{x}_{\mathbf{l}}+d_{3}\bm{\gamma}\bm{x}_{\mathbf{l}}\leq
b_{3},$ (16e) $\displaystyle\hskip
75.0001pta_{4}\bm{x}_{\mathbf{l}}+d_{4}\Delta{\bm{p}}\leq b_{4},$ (16f)

where $\bm{x}_{\mathbf{l}}$ is the follower problem variable and
$\bm{\gamma},\Delta{\bm{p}}$ represent the upper-level variables. Therefore,
the follower problem (16c)-(16f) is linear in $\bm{x}_{\mathbf{l}}$ even
though we have bilinear terms in (16e).

By taking dual of the lower-level problem (16c)-(16f) and applying the strong
duality theorem [33], we reduce the bilevel problem (16) to a single-level
optimization problem given by

$\displaystyle\max_{\Delta{\bm{p}},\bm{\gamma},\bm{x}_{\mathbf{l}},\bm{\lambda}}~{}$
$\displaystyle\Delta{\bm{p}}$ (17a) s.t.
$\displaystyle\bm{{x}}_{\mathbf{l}}\leq b_{1},$ (17b) $\displaystyle
a_{2}\bm{x}_{\mathbf{l}}\leq b_{2},$ (17c) $\displaystyle
a_{3}\bm{x}_{\mathbf{l}}+d_{3}\bm{\gamma}\bm{x}_{\mathbf{l}}\leq b_{3},$ (17d)
$\displaystyle a_{4}\bm{x}_{\mathbf{l}}+d_{4}\Delta{\bm{p}}\leq b_{4},$ (17e)
$\displaystyle\bm{x}_{\mathbf{l}}\geq{\begin{bmatrix}b_{2}&b_{3}&(b_{4}-d_{4}\Delta\bm{p})\end{bmatrix}}\bm{\lambda},$
(17f)
$\displaystyle{\begin{bmatrix}a_{2}&(a_{3}+d_{3}\bm{\gamma})&a_{4}\end{bmatrix}}\bm{\lambda}=1,$
(17g) $\displaystyle\bm{\lambda}\geq 0,$ (17h)

where (17c)-(17e) and (17g)-(17h) are the primal and dual constraints for the
lower level problem, respectively, and (17f) enforces strong duality.

A similar approach can be used to derive the single-level reduction of
$\text{P}^{\pm}$ which comprises of multiple follower problems. We can also
choose other inverter reactive power control modes and formulate the bilevel
problem and the single-level reformulated problem in the same way as described
above for the constant power factor mode. Note that the presence of bilinear
terms $\bm{\gamma}\bm{x}_{\mathbf{l}}$ in the primal follower problem
constraint (17d), $\Delta\bm{p}\bm{\lambda}$ in the strong duality constraint
(17f) and $\bm{\gamma}\bm{\lambda}$ in the dual equality constraint (17g) make
the problem challenging to solve. We will next discuss an iterative solution
approach to handle the bilinear terms and solve (17) efficiently.

###  IV-B Solution Method

As the number of follower problems in (16) increases, the number of bilinear
(and non-convex) terms seen in the single-level reformulated problem (17) also
increase. Typically, the voltage magnitude violations in a distribution grid
occur only at a few nodes, which in turn determine the aggregate power
flexibility of the system. Thus, instead of solving the bilevel problem with
all $4n$ follower problems included, we use an iterative approach where we
identify the nodes where voltage violations are most likely to occur and only
include the lower-level problems corresponding to these nodes in the bilevel
optimization problem (17). The iterative approach is illustrated in Fig. 3 and
comprises of the following steps:

Figure 3: Iterative solution method to solve bilevel problem. Figure 5: Worst-
case upper (orange) and lower (blue) aggregate power flexibility limits for
every single-phase node in the IEEE-13 node feeder. From left to right:
Results of different inverter control modes.

####  IV-B1 Worst case

We first identify the aggregate flexibility range by assuming that the
inverter set-points can take on worst-case values (instead of being chosen by
the DSO). In this case, the inverter set-points become lower level variables.
Solving the resulting problem is straightforward for the constant reactive
power mode with linear inverter constraints (14). However, the inverter
constraints (13a), (15a) for constant power factor and voltage-reactive power
modes include bilinear terms for this case. To avoid this problem of non-
convexity, we fix the inverter set-points. If the lower-level objective is to
maximize voltage magnitude, the worst-case inverter setpoints are determined
as:

  * •

_Constant power factor mode_ : We set $\bm{\gamma}_{\mathbf{G},i}^{\phi}$ for
every inverter to its upper limit defined in (13b). This corresponds to the
situation where all inverters are boosting the voltage profile in the network
by injecting the full available reactive power, leading to higher voltage
magnitudes.

  * •

_Voltage-reactive power mode_ : We set
$\overline{\bm{q}}_{\mathbf{G},i}^{\phi}$ to its lower limit defined in (15b).
This represents the condition where inverters provide no reactive power
support when voltage magnitude is close to or violating the upper limit
$\overline{v}$.

If the lower-level objective is to minimize the voltage magnitude, we set the
corresponding variables to the other extreme.

With these assumptions, $\Delta\bm{p}^{+},\Delta\bm{p}^{-}$ are the only
variables that are shared by all the lower-level problems, i.e. there are no
other variables or constraints linking the $4n$ follower problems. This allows
us to solve each follower problem separately. Specifically, we solve $4n$
worst case problems with one leader and one follower each to obtain $2n$
aggregate power flexibility upper limits for the positive case (with
$\Delta\bm{p}^{-}=0$) and $2n$ aggregate power flexibility lower limits for
the negative case (with $\Delta\bm{p}^{+}=0$). The difference between the
minimum upper limit and the maximum lower limit is the worst case aggregate
power flexibility range for the network corresponding to the worst-case
inverter control. Note that each of the $4n$ small worst case problems has
only one bilinear term in the strong duality inequality constraint. It is
faster to solve the $4n$ smaller problems than one large problem with $4n$
followers, especially if solved in parallel.

####  IV-B2 Ideal case

Next, we solve the ideal case problem P± as defined in Section III-E. This
problem cannot be broken down into several bilevel problems because we have
multiple upper-level variables shared by the lower-level problems. To reduce
the number of bilinear terms, we initially only include the two lower-level
problems corresponding to the two problems that determined the worst case
aggregate power flexibility range in step 1). This represents a relaxation of
the original optimization problem and by solving it, we obtain an upper bound
on the new aggregate power flexibility range and proposed inverter setpoints.
Note that this ideal case problem comprises of multiple constraints with
bilinear terms. Apart from the dual and strong duality constraints, bilinear
terms are also present in the inverter constraints (13a) and (15a) for the
constant power factor and voltage-reactive power mode, respectively, which
make the problem non-convex. Commercial solvers such as Gurobi can handle such
constraints efficiently using spatial branching [34] if the number of bilinear
terms are small. We use this approach for solving the ideal case problem.
Alternatively, it is possible to employ piecewise McCormick envelopes to relax
the bilinear terms [35], but we defer this to future work.

####  IV-B3 Feasibility check

The solution obtained by solving the ideal case problem in step 2) may allow
for voltage violations at some single-phase nodes in the network since not all
follower problems were considered inside the problem. Therefore, we perform a
feasibility check where we fix the aggregate power flexibility range and
inverter setpoints to the values determined in step 2) and solve all $4n$
lower-level problems defined in $\text{P}^{\pm}$. If any of the lower-level
problems have maximum or minimum voltage magnitudes outside the limits, we go
back to step 2) and add this follower problem to the ideal case problem.

We iterate between step 2) and 3) in this way until the feasibility check
determines that all voltage magnitudes obtained from the $4n$ lower level
problems are within limits.

##  V Case Study

In this section, we evaluate our proposed approach using the IEEE 13-node
feeder [36] and one of the taxonomic distribution feeders from Pacific
Northwest National Laboratory (PNNL) [37]. The optimization problem is
implemented in Julia and the optimization problem is solved using Gurobi [34].
All simulations were run on a Windows 10 PC with 3.00 GHz Intel Xeon processor
and 16 GB RAM. For our analysis, we solve three optimization problems based on
the three inverter reactive power control modes. For the constant power factor
mode, we use the inverter power factor $pf_{\text{G},i}^{\phi}=0.9$ to define
the upper and lower bounds in (13b) and for the constant reactive power mode,
we use the power ratio $\gamma_{\text{G},i}^{\phi}=0.48$ to define upper and
lower bounds in (14). The load power factor $pf_{\text{L},i}^{\phi}$ at node
$i$ connected to phase $\phi$ is set to $0.95$.

###  V-A IEEE 13-node Feeder

Fig. 4 shows the IEEE 13-node feeder with 75 single-phase solar PV
installations at seven nodes, illustrated by house blocks where each block
corresponds to aggregation of 5 houses. The maximum apparent power rating of
each single-phase solar PV inverter is 60 kVA. We have chosen a current
operating point where the PV penetration level, calculated as the ratio of
total PV generation (in kW) to the total rated load (in kW), is 45%. The
voltage limits are set to $\underline{v}=0.9$ p.u. and $\overline{v}=1.1$ p.u.
The taps of the voltage regulator connecting nodes 650 and 630 are set to high
values and as a result, the system is more prone to overvoltage conditions
compared to undervoltage scenarios. So, we will only discuss results which
focus on identifying aggregate power flexibility limits
$\Delta{\bm{p}^{-}},\Delta{\bm{p}^{+}}$ that mitigate overvoltages in the
feeder (i.e. the lower-level problem objective is to maximize voltage
magnitude). The maximum available aggregate power flexibility in the system is
$\Delta\overline{p}=-\Delta\underline{p}=1.64$ MW.

Figure 4: Modified IEEE-13 node feeder. Figure 6: Optimal setpoints (blue
dots) for different inverter reactive power control modes for the IEEE-13 node
feeder. The red dashed lines represent the maximum and minimum limits of the
control setpoints.

####  V-A1 Worst case aggregate power flexibility limits

The first step of our iterative solution approach is to solve the worst case
problem for each single-phase node separately and identify the location where
overvoltages are most likely going to occur. Fig. 5 shows the aggregate power
flexibility limits obtained by solving the worst case problem for the three
inverter reactive power control modes. For the constant power factor mode, the
lower limit on the aggregate power flexibility $\Delta\bm{p}^{-}$ (blue
triangles) is equal to $-1.64$ MW for most of the single-phase nodes, but much
closer to zero for some. In contrast, the upper limit $\Delta\bm{p}^{+}$
(orange triangles) is much lower than $1.64$ MW for most single-phase nodes.
For the constant reactive power mode, the aggregate power flexibility range
(i.e. the distance between the blue and orange triangles) is zero for multiple
single-phase nodes. For the voltage-reactive power mode, the aggregate power
flexibility range at almost all the single-phase nodes is equal to the maximum
available aggregate power flexibility. From the above results, we see that the
choice of reactive power control mode has a very significant impact on the
worst-case aggregate power flexibility range. However, for all three inverter
reactive power control modes, the smallest aggregate power flexibility range
was obtained when maximizing the voltage magnitude at phase $b$ of node 675
(single-phase node #46 in Fig. 5). Furthermore, the average computation time
to solve the worst case problem with one follower is less than 0.1 seconds.

####  V-A2 Actual aggregate power flexibility limits

We next solve for the true aggregate power flexibility limits by following the
procedure described in Section IV-B. The ideal case problem starts with one
lower-level problem corresponding to node 675 at phase $b$, and additional
problems are iteratively added following the feasibility check. Table II
summarizes the number of iterations and the aggregate power flexibility limits
obtained after termination of the iterative process.

First, we observe that the algorithm only requires a few iterations,
regardless of which converter control mode is used. While only one iteration
(i.e. only a single follower problem) is required for the problem with
constant reactive power mode, the constant power factor and voltage-reactive
power modes require that more follower problems are added to the ideal case
problem. With any of the inverter control modes, the computation time to solve
the ideal case problem in the last iteration is less than a second. This
indicates that the iterative approach can solve the problem efficiently.

Next, we observe that the aggregate power flexibility limits are same for all
inverter reactive power control modes and they are equal to the maximum
available aggregate power flexibility in the system. This indicates that if
the DSO is able to control the reactive power settings of the inverters, then
the amount of flexibility offered by the DERs connected to the grid can be
leveraged to the maximum technical limit.

TABLE II: Aggregate power flexibility and number of iterations for the IEEE-13
node feeder with different inverter control modes

Inverter mode $\Delta{\bm{p}^{-}}$ $\Delta{\bm{p}^{+}}$ No. of (MW) (MW)
Iterations Constant power factor -1.64 1.64 3 Constant reactive power -1.64
1.64 1 Voltage-reactive power -1.64 1.64 3

Fig. 6 illustrates the optimal inverter setpoints for the different inverter
control modes. For the constant power factor and constant reactive power
modes, the reactive-power injections for most of the solar PV systems are
negative (reflected as either a negative power ratio or a negative reactive
power injection). This is as expected since the absorption of the reactive
power will lead to lowering of the voltage magnitudes thereby avoiding
overvoltage conditions. For the voltage-reactive power mode, the maximum
reactive power limit is close to the apparent power limit of the inverters
(red dashed line) to ensure that the reactive power absorbed by the solar PV
systems is high, which again helps keep the voltage magnitudes within limits.

####  V-A3 Linear approximation accuracy

Figure 7: Worst-case voltage magnitudes obtained by solving nonlinear (blue
circles) and linear (orange crosses) versions of the lower-level problem for
IEEE 13-node feeder. The red dashed lines are the upper voltage magnitude
limits. From left to right: Results of different inverter control modes.

To obtain a linear lower-level optimization problem, we linearized multiple
constraints. To analyze the accuracy of the inverter setpoints obtained by
solving the linearized problem, we formulate the nonlinear counterpart of the
lower-level problem by making the following changes. (i) We replace the
linearized inverter constraints (12) with the original quadratic constraints.
(ii) We replace the first order Taylor approximation of the relationship
between the voltage variables defined in (5) with the original quadratic
constraints. (iii) We use the nonlinear power flow equation in rectangular
form [38] instead of the fixed-point power flow equations (10).

We then fix the values of the upper level variables (i.e. the aggregate power
flexibility limits and inverter setpoints) to the solution identified above
and solve the nonlinear subproblem for all single-phase nodes to check if the
worst-case voltage magnitudes are within the limits. We also compare the
resulting voltage magnitudes with the voltage magnitudes that are obtained
from the bilevel problem with the linearized formulation.

Figure 9: R2-12-47-2 feeder results obtained by solving ideal case problem. We
show the setpoints determined for different inverter modes. The red dashed
lines are maximum and minimum limits on the upper-level variables.

Fig. 7 shows the worst-case voltage magnitudes obtained by solving the
nonlinear and linear versions of the problem, for all different inverter
control modes. Across the three different control modes, we see that the
approximation accuracy is high, i.e. the nonlinear voltage magnitudes (blue
circles) and linear voltage magnitudes (orange crosses) are close to each
other. The maximum error is $\approx$ 0.006 p.u., and occurs for the constant
reactive power mode. The reason for the higher inaccuracy in this case may be
that the solution is further away from the initial operating point which was
used to linearize the power flow equations in (10). Importantly, for all three
inverter reactive power control modes, the set-points determined by the
linearized problem keep the voltage magnitudes within the limits even we solve
the nonlinear problem.

###  V-B PNNL Taxonomic Feeder: R2-12-47-2

Figure 8: Modified R2-12-47-2 taxonomic feeder [39] visualized using [40].

To investigate scalability of our solution method, we next run simulations on
the R2-12-47-2 feeder shown in Fig. 8 which comprises of 820 single-phase
nodes [39]. Solar PV inverters are connected to 556 single-phase nodes in the
feeder, and each inverter is rated at 5 KVA to achieve a PV penetration level
of 20% of the total rated load. The voltage limits for this test case are
tighter and set to $\underline{v}=0.95$ p.u. and $\overline{v}=1.05$ p.u. The
maximum available aggregate power flexibility before considering voltage
constraints is $\pm 4.86$ MW. In this case, we solve the bilevel problem to
determine the aggregate power flexibility and inverter setpoints that mitigate
both overvoltages and undervoltages in the feeder.

####  V-B1 Worst case aggregate power flexibility limits

We first solve the worst-case problem for every single-phase node in the
network. Different from the IEEE-13 node feeder results, it was observed that
all three inverter reactive power control modes provide an aggregate power
flexibility range of zero at multiple single-phase nodes in the network.
Furthermore, the average computation time to solve the worst-case problem at
each single-phase nodes is about 2 seconds. Since there are multiple single-
phase nodes with aggregate power flexibility range of zero, we randomly pick
one of the them as the first set of lower-level problems to be included in the
ideal case problem.

####  V-B2 Actual aggregate power flexibility limits

Next, we compare the results obtained after termination of the iterative
process where the ideal case problem is solved repeatedly until the
feasibility check is passed. Similar to the IEEE 13-node results, the
aggregate power flexibility limits determined for all inverter modes are equal
to the maximum allowable aggregate power flexibility $\pm 4.86$ MW in the
system. It is interesting to note that the feasibility check passed after the
first iteration for all three inverter modes. In addition, the maximum
computation time to solve the ideal case problem with one follower was
approximately one minute.

The inverter setpoints for the different reactive power modes are shown in
Fig. 9. Compared to the IEEE-13 node feeder results, most of the power ratio
setpoints for the constant power factor mode are no longer close to the lower
limit. Instead, we see that multiple inverters are injecting reactive power
into the grid. This might be because the lower PV penetration level reduces
the probability of overvoltages and increases the risk of undervoltages. The
inverter setpoints obtained for the reactive power mode are much closer to
zero than the lower limit which is similar to the reactive power setpoints
obtained for the IEEE-13 node feeder in Fig. 6(b). For the voltage-reactive
power mode, the maximum reactive power setpoints are closer to the lower limit
as opposed to the IEEE-13 node results shown in Fig. 6(c) where the setpoints
are near the upper limit. This again indicates that a high reactive power
absorption is not required for this test case.

##  VI Conclusion

This paper explores the coordination between transmission systems,
distribution grids and DER aggregators. The main goal is to identify the
aggregate power flexibility range that can be provided from distribution-
connected resources to the transmission system without causing any constraint
violations within the distribution grid itself. We formulate a bilevel
optimization task where the upper-level problem represents the decision making
of the DSO, which maximizes the aggregate power flexibility and determines DER
inverter reactive power setpoints that ensure no voltage violations occur,
while the lower-level problem determines the worst-case disaggregation
strategy (i.e. assuming that whoever controls the DERs have no information or
interest in enforcing internal distribution grid constraints). Our case
studies demonstrate that it is important to allow the DSO to choose the DER
inverter reactive power setpoints in order to utilize the maximum available
flexibility in the system. These setpoints will ensure that the grid is secure
even during worst-case conditions. The results for the IEEE 13-bus feeder
showed that we are able to obtain high-quality solutions for the inverter
setpoints by solving the single-level, strong duality based reformulated
problem. Our results on the larger taxonomic feeder demonstrate that the
proposed approach can be applied to large, realistic distribution feeders. In
both cases, we are able to enable DERs to use their full range of flexibility,
even though our formulation assumes worst-case behavior of aggregators.

For future work, we will focus on integrating DERs such as electric vehicles
and battery storage systems with temporally-coupled constraints in our
optimization framework. We also plan to investigate the impact of limited
information about the current operating point on the aggregate flexibility
range obtained by the bilevel problem.

## References

  * [1] C. D’Adamo, S. Jupe, and C. Abbey, “Global survey on planning and operation of active distribution networks-update of CIGRE c6. 11 working group activities,” in _CIRED 2009-20th International Conference and Exhibition on Electricity Distribution-Part 1_. IET, 2009, pp. 1–4.
  * [2] FERC, “FERC order no. 2222: A new day for distributed energy resources.” [Online]. Available: https://www.ferc.gov/media/ferc-order-no-2222-fact-sheet
  * [3] Y. Yi and G. Verbic, “Operating envelopes under probabilistic electricity demand and solar generation forecasts,” _arXiv e-prints_ , pp. arXiv–2207, 2022.
  * [4] N. Nazir and M. Almassalkhi, “Convex inner approximation of the feeder hosting capacity limits on dispatchable demand,” in _2019 IEEE 58th Conference on Decision and Control (CDC)_. IEEE, 2019, pp. 4858–4864.
  * [5] ——, “Grid-aware aggregation and realtime disaggregation of distributed energy resources in radial networks,” _IEEE Transactions on Power Systems_ , vol. 37, no. 3, pp. 1706–1717, 2021.
  * [6] J. Silva, J. Sumaili, R. J. Bessa, L. Seca, M. A. Matos, V. Miranda, M. Caujolle, B. Goncer, and M. Sebastian-Viana, “Estimating the active and reactive power flexibility area at the TSO-DSO interface,” _IEEE Transactions on Power Systems_ , vol. 33, no. 5, pp. 4741–4750, 2018.
  * [7] X. Xu, Q. Lyu, M. Qadrdan, and J. Wu, “Quantification of flexibility of a district heating system for the power grid,” _IEEE Transactions on Sustainable Energy_ , vol. 11, no. 4, pp. 2617–2630, 2020.
  * [8] E. Polymeneas and S. Meliopoulos, “Aggregate modeling of distribution systems for multi-period OPF,” in _2016 Power Systems Computation Conference (PSCC)_. IEEE, 2016, pp. 1–8.
  * [9] L. Ageeva, M. Majidi, and D. Pozo, “Analysis of feasibility region of active distribution networks,” in _2019 International Youth Conference on Radio Electronics, Electrical and Power Engineering (REEPE)_. IEEE, 2019, pp. 1–5.
  * [10] M. Heleno, R. Soares, J. Sumaili, R. J. Bessa, L. Seca, and M. A. Matos, “Estimation of the flexibility range in the transmission-distribution boundary,” in _2015 IEEE Eindhoven PowerTech_. IEEE, 2015, pp. 1–6.
  * [11] Z. Tan, H. Zhong, Q. Xia, C. Kang, X. S. Wang, and H. Tang, “Estimating the robust PQ capability of a technical virtual power plant under uncertainties,” _IEEE Transactions on Power Systems_ , vol. 35, no. 6, pp. 4285–4296, 2020.
  * [12] F. L. Müller, J. Szabó, O. Sundström, and J. Lygeros, “Aggregation and disaggregation of energetic flexibility from distributed energy resources,” _IEEE Transactions on Smart Grid_ , vol. 10, no. 2, pp. 1205–1214, 2017.
  * [13] L. Zhao, W. Zhang, H. Hao, and K. Kalsi, “A geometric approach to aggregate flexibility modeling of thermostatically controlled loads,” _IEEE Transactions on Power Systems_ , vol. 32, no. 6, pp. 4721–4731, 2017.
  * [14] H. Hao, B. M. Sanandaji, K. Poolla, and T. L. Vincent, “Aggregate flexibility of thermostatically controlled loads,” _IEEE Transactions on Power Systems_ , vol. 30, no. 1, pp. 189–198, 2014.
  * [15] X. Chen, E. Dall’Anese, C. Zhao, and N. Li, “Aggregate power flexibility in unbalanced distribution systems,” _IEEE Transactions on Smart Grid_ , vol. 11, no. 1, pp. 258–269, 2019.
  * [16] S. Wang and W. Wu, “Aggregate flexibility of virtual power plants with temporal coupling constraints,” _IEEE Transactions on Smart Grid_ , vol. 12, no. 6, pp. 5043–5051, 2021.
  * [17] B. Cui, A. Zamzam, and A. Bernstein, “Network-cognizant time-coupled aggregate flexibility of distribution systems under uncertainties,” in _2021 American Control Conference (ACC)_. IEEE, 2021, pp. 4178–4183.
  * [18] X. Chen and N. Li, “Leveraging two-stage adaptive robust optimization for power flexibility aggregation,” _IEEE Transactions on Smart Grid_ , vol. 12, no. 5, pp. 3954–3965, 2021.
  * [19] C. Zhang, Q. Wang, J. Wang, P. Pinson, J. M. Morales, and J. Østergaard, “Real-time procurement strategies of a proactive distribution company with aggregator-based demand response,” _IEEE Transactions on Smart Grid_ , vol. 9, no. 2, pp. 766–776, 2016.
  * [20] C. Zhang, Q. Wang, J. Wang, P. Pinson, and J. Østergaard, “Real-time trading strategies of proactive disco with heterogeneous dg owners,” _IEEE Transactions on Smart Grid_ , vol. 9, no. 3, pp. 1688–1697, 2016.
  * [21] P. Sheikhahmadi, S. Bahramara, A. Mazza, G. Chicco, and J. P. Catalão, “Bi-level optimization model for the coordination between transmission and distribution systems interacting with local energy markets,” _International Journal of Electrical Power & Energy Systems_, vol. 124, p. 106392, 2021.
  * [22] S. Riaz, H. Marzooghi, G. Verbič, A. C. Chapman, and D. J. Hill, “Generic demand model considering the impact of prosumers for future grid scenario analysis,” _IEEE Transactions on Smart Grid_ , vol. 10, no. 1, pp. 819–829, 2017.
  * [23] M. Zugno, J. M. Morales, P. Pinson, and H. Madsen, “A bilevel model for electricity retailers’ participation in a demand response market environment,” _Energy Economics_ , vol. 36, pp. 182–197, 2013.
  * [24] S. Ross and J. Mathieu, “Strategies for network-safe load control with a third-party aggregator and a distribution operator,” _IEEE Transactions on Power Systems_ , vol. 36, no. 4, pp. 3329–3339, 2021.
  * [25] Z. Li, Q. Guo, H. Sun, and J. Wang, “Coordinated economic dispatch of coupled transmission and distribution systems using heterogeneous decomposition,” _IEEE Transactions on Power Systems_ , vol. 31, no. 6, pp. 4817–4830, 2016\.
  * [26] C. Lin, W. Wu, X. Chen, and W. Zheng, “Decentralized dynamic economic dispatch for integrated transmission and active distribution networks using multi-parametric programming,” _IEEE Transactions on Smart Grid_ , vol. 9, no. 5, pp. 4983–4993, 2017.
  * [27] C. Lin, W. Wu, B. Zhang, B. Wang, W. Zheng, and Z. Li, “Decentralized reactive power optimization method for transmission and distribution networks accommodating large-scale DG integration,” _IEEE Transactions on Sustainable Energy_ , vol. 8, no. 1, pp. 363–373, 2016.
  * [28] D. Molzahn and L. A. Roald, “Grid-aware versus grid-agnostic distribution system control: A method for certifying engineering constraint satisfaction,” in _Proceedings of the 52nd Hawaii International Conference on System Sciences_ , 2019.
  * [29] K. Girigoudar and L. A. Roald, “Linearized three-phase optimal power flow models for distribution grids with voltage unbalance,” in _2021 60th IEEE Conference on Decision and Control (CDC)_. IEEE, 2021, pp. 4214–4221.
  * [30] D. G. Photovoltaics and E. Storage, “IEEE standard for interconnection and interoperability of distributed energy resources with associated electric power systems interfaces,” _IEEE Std_ , pp. 1547–2018, 2018.
  * [31] A. Bernstein, C. Wang, E. Dall’Anese, J.-Y. Le Boudec, and C. Zhao, “Load flow in multiphase distribution networks: Existence, uniqueness, non-singularity and linear models,” _IEEE Transactions on Power Systems_ , vol. 33, no. 6, pp. 5832–5843, 2018.
  * [32] X. Chen, W. Wu, and B. Zhang, “Robust restoration method for active distribution networks,” _IEEE Transactions on Power Systems_ , vol. 31, no. 5, pp. 4005–4015, 2015.
  * [33] B. Zeng and Y. An, “Solving bilevel mixed integer program by reformulations and decomposition,” _Optimization online_ , pp. 1–34, 2014.
  * [34] Gurobi Optimization, LLC, “Gurobi Optimizer Reference Manual,” 2022. [Online]. Available: https://www.gurobi.com
  * [35] P. M. Castro, “Tightening piecewise mccormick relaxations for bilinear problems,” _Computers & Chemical Engineering_, vol. 72, pp. 300–311, 2015\.
  * [36] K. Schneider, B. Mather, B. Pal, C.-W. Ten, G. Shirek, H. Zhu, J. Fuller, J. Pereira, L. Ochoa, L. De Araujo _et al._ , “Analytic considerations and design basis for the IEEE distribution test feeders,” _IEEE Transactions on Power Systems_ , vol. 33, no. 3, pp. 3181–3188, 2017.
  * [37] K. P. Schneider, Y. Chen, D. P. Chassin, R. G. Pratt, D. W. Engel, and S. E. Thompson, “Modern grid initiative distribution taxonomy final report,” Pacific Northwest National Lab.(PNNL), Richland, WA (United States), Tech. Rep., 2008.
  * [38] J. Arrillage, C. P. Arnold, and B. Harker, “Computer modelling of electrical power systems,” 1983.
  * [39] PNNL, “GridLAB-D taxonomy feeder.” [Online]. Available: https://github.com/gridlab-d/Taxonomy_Feeders
  * [40] NRECA, “Open modeling framework.” [Online]. Available: https://github.com/dpinney/omf


David A. Cox, John B. Little, and Henry K. Schenck.
Toric varieties, volume 124 of Graduate Studies in
American Mathematical Society, Providence, RI, 2011.

[10]
William Fulton.
Introduction to toric varieties, volume 131 of Annals of
Mathematics Studies.
Princeton University Press, Princeton, NJ, 1993.
The William H. Roever Lectures in Geometry.

[11]
Catherine Greenhill and Brendan D. McKay.
Asymptotic enumeration of sparse multigraphs with given degrees.
SIAM J. Discrete Math., 27(4):2064–2089, 2013.

[12]
R. C. King, C. Tollu, and F. Toumazet.
The hive model and the polynomial nature of stretched
Littlewood-Richardson coefficients.
Sém. Lothar. Combin., 54A:Art. B54Ad, 19, 2005/07.

[13]
Jesús A De Loera and Bernd Sturmfels.
Algebraic unimodular counting.
Mathematical programming, 96:183–203, 2003.

[14]
Jesús A. De Loera, Raymond Hemmecke, Jeremiah Tauzer, and Ruriko Yoshida.
Effective lattice point counting in rational convex polytopes.
Journal of Symbolic Computation, 38(4):1273 – 1302, 2004.
Symbolic Computation in Algebra and Geometry.

[15]
Brendan D. McKay, Ian M. Wanless, and Nicholas C. Wormald.
Asymptotic enumeration of graphs with a given upper bound on the
maximum degree.
Combin. Probab. Comput., 11(4):373–392, 2002.

[16]
Eckhard Meinrenken and Reyer Sjamaar.
Singular reduction and quantization.
Topology, 38(4):699–762, 1999.

[17]
Marni Mishna, Mercedes Rosas, and Sheila Sundaram.
Vector partition functions and Kronecker coefficients.
Journal of Physics A: Mathematical and Theoretical,
54(20):205204, May 2021.

[18]
Marni Mishna and Stefan Trandafir.
Estimating and computing Kronecker coefficients: a vector partition
function approach.
<https://arxiv.org/abs/2210.12128>, 2022.

[19]
Tadao Oda.
Convex bodies and algebraic geometry—toric varieties and
applications. I.
World Sci. Publishing, Singapore, 1988.

[20]
Etienne Rassart.
A polynomiality property for Littlewood-Richardson coefficients.
J. Combin. Theory Ser. A, 107(2):161–179, 2004.

[21]
Alexander Schrijver.
Theory of linear and integer programming.
Wiley-Interscience Series in Discrete Mathematics. John Wiley &
Sons, Ltd., Chichester, 1986.
A Wiley-Interscience Publication.

[22]
Bernd Sturmfels.
On vector partition functions.
Journal of Combinatorial Theory, Series A, 72(2):302 – 309,

[23]
Sven Verdoolaege, Matthias Koeppe, and Kevin M Woods.
An implementation of the projection theorem of Barvinok and
In DTAI seminar, Date: 2008/03/10-2008/03/10, Location: Leuven,
Belgium, 2008.

[24]
Sven Verdoolaege and Kevin Woods.
Counting with rational generating functions.
J. Symbolic Comput., 43(2):75–91, 2008.