arXiv:1803.06921v1 [] 14 Mar 2018

Approximating Flexibility in Distributed Energy
Resources: A Geometric Approach
Soumya Kundu, Karanjit Kalsi

Scott Backhaus

Optimization and Control Group
Pacific Northwest National Laboratory
Richland, WA 99352 USA
Email: {soumya.kundu, karanjit.kalsi}@pnnl.gov

Information Systems and Modeling Group
Los Alamos National Laboratory
Los Alamos, NM 87545 USA
Email: backhaus@lanl.gov

Abstract—With increasing availability of communication and
control infrastructure at the distribution systems, it is expected
that the distributed energy resources (DERs) will take an
active part in future power systems operations. One of the
main challenges associated with integration of DERs in grid
planning and control is in estimating the available flexibility
in a collection of (heterogeneous) DERs, each of which may
have local constraints that vary over time. In this work, we
present a geometric approach for approximating the flexibility of
a DER in modulating its active and reactive power consumption.
The proposed method is agnostic about the type and model of
the DERs, thereby facilitating a plug-and-play approach, and
allows scalable aggregation of the flexibility of a collection of
(heterogeneous) DERs at the distributed system level. Simulation
results are presented to demonstrate the performance of the
proposed method.
Index Terms—Demand response, load aggregation, Minkowski
sum, polynomial optimization.

I. I NTRODUCTION
Traditionally the bulk of the responsibility of maintaining the (real-time) balance between load and generation has
rest upon the conventional generators, e.g. spinning reserves,
which provide various forms of grid ancillary services such as
inertia support, regulation and ramping. As the power systems
transition towards a greener grid with larger penetration of
renewable generation, a large fraction of which is expected
to be distributed, there is a general consensus that various
forms of distributed energy resources (DERs), including flexible and responsive electrical loads, need to be coordinated
and controlled in real-time to provide grid support. This
requires an appropriate understanding of the loads behavior,
including their physical models, at the details useful for realtime coordination. It is important that such models capture
the necessary information, such as the associated dynamics
(thermal dynamics of a residential air-conditioner), measure
of the available flexibility (reactive power capacity of a gridconnected inverter), the end-user constraints (desired stateof-charge of an electric vehicle battery), etc., while being
tractable for real-time operation.

Modeling of an ensemble of flexible loads for ancillary
services (in particular, frequency regulation and ramping) have
been explored in the literature in recent years [1]–[8]. The
methods proposed in these articles are applicable to ensembles
of similar loads (albeit with heterogeneous parameters), such
as a collection of residential air-conditioners, or a collection
of plug-in electric vehicles. Such types of aggregate flexibility
models are suitable for a transmission system operator which
views the net demand flexibility available at the distributionlevel as a lumped-model. Aggregation at this level (tens of
thousands of loads) do not explicitly take into account the
operational constraints (line-flow and voltage limits) at the
mid/low-voltage distribution systems. As the fraction of flexible loads increase, however, chances of violation of operational
constraints due to control of flexible demand will increase.
In order to efficiently coordinate tens of thousands of
flexible loads in distribution systems, while also satisfying
line-flow and node voltage constraints, hierarchical modeling
and control frameworks [9], [10] become attractive. Load
aggregators, referred to as aggregate device controller (or
ADC), help in keeping the size of control problem tractable
by aggregating the neighboring DERs locally. This aggregation
can be done at the level of a couple of service transformers
(tens of residential customers). It is the responsibility of the
ADC to capture the aggregated flexibility, in terms of active
and reactive power, of the local DERs which are likely to
be of different types (e.g. a collection of air-conditioners,
electric water-heaters, batteries, solar photovoltaic inverers
wind inverters) and ratings.
The rest of the paper is organized as follows. In Section II
we describe the problem of aggregating flexibility of heterogeneous DERs. Section III presents some key concepts that form
the basis of our work. In Section IV we discuss in details
the geometric programming approach to the aggregation of
flexibility, while exploring the metrics of quality of approximation in Section V. Numerical results illustrating the concept
are presented in Section VI, before concluding the article in
Section VII.
II. P ROBLEM D ESCRIPTION

This work was supported by the United States Department of Energy under
the Grid Modernization Lab Consortium initiative.

Consider a hierarchical distribution system operation framework in which a distribution system operator (DSO) com-

municates with aggregate device controllers (ADCs) at every
control period and solves an optimal dispatch problem to
schedule the DERs. It is assumed that the ADCs are so
placed that the operational constraints (line-flow and voltage
limits) are trivially satisfied at the ADC-level. The primary
responsibilities of the ADC are two-fold: 1) before the start
of each control period (5-15 min), the ADC estimates the net
flexibility (in active and reactive power) available among the
DERs under its control, and exchanges that information to the
DSO; 2) during the control period, the ADC coordinates the
DERs in real-time to track the active and reactive power setpoints dispatched by the DSO. In this paper we restrict our
discussion only on the first part, i.e. the problem of estimating
the net flexibility in a group of (dissimilar) DERs.
For the purpose of this article, any energy (consuming or
generating) resource which offers certain flexibility in active
(p) and/or reactive (q) power, possibly over a reasonably short
time window (such as a 5-15 min long control period), is
considered as a DER. We use the notation F to represent
the flexibility domain as a collection of (p, q)-points that are
physically admissible by the DER (possibly via some local
device-level control). Note that the flexibility domain (F )
could be a continuous or a discrete domain. For example, a
solar photovoltaic (PV) inverter that can modulate its active
and reactive power over a continuous range will have a
continuous flexibility domain, while the flexibility domain
for switching loads, such as an air-conditioner or electric
water-heater, will be discrete. Moreover this flexibility is timevarying, and depends on exogenous parameters as well as enduser preferences.
Fig. 1 shows examples of flexibility domains for certain
types of DERs, with positive (negative) values of p and q
denoting consumption (generation). Discrete flexibility domain
of a switching load (e.g. air-conditioner) that operates in
two discrete operational states (‘on’ and ‘off’) is shown in
Fig. 1(b), while the rest of the plots represent continuous
flexibility domains. Batteries (Fig. 1(a) offer full four-quadrant
flexibility, while PV (Fig. 1(c)) and wind (Fig. 1(d)) inverters
offer flexibility only on the left half-plane (active power
generation). It must be noted that the flexibility offered by
the DERs is dynamic, and change based on end-usage and
exogenous influence. For example, if there is a cloudy sky,
the PV inverter output might only be restricted to a small
fraction of its rated generation. Similarly, the air-conditioner
may be forced to operate mostly in ‘on’-state, if the outside
air-temperature is high. Finally, the flexibility domains should
be able to capture DER uncertainties, which could be modeled
in the robust sense via a conservative estimate of the available
flexibility.
In this article, we will focus our discussion on the challenges
of aggregating the flexibility of a heterogeneous mix of DERs
at the ADC, while the issues regarding the dynamic evolution
of flexibility and the associated uncertainties will be addressed
in future work. Let us assume that there are N DERs with an
ADC. The flexibility domain of the i-th DER is denoted by
Fi , such that its active and reactive power consumption (with

negative value signifying net generation)
(pi , qi ) ∈ Fi ,

i ∈ {1, 2, . . . , N } .

The goal of the flexibility aggregation task is to find the net
flexibility domain F in the form of a Minkowski sum of the
individual DER flexibility domains (Fi ), such that,


PN
N
p = i=1 pi 

]
PN
Fi = (p, q)
.
F :=
q = i=1 qi


i=1
(pi , qi ) ∈ Fi ∀i

Calculating the exact Minkowski sum of the individual flexibility domains is computationally complex, especially as the
number of DERs increases. Thus, from a practical point-ofview, a desirable approach is to construct approximations of
the aggregate flexibility domain at the ADC in a scalable
way. In the rest of the paper, we will discuss methods of
computing the approximated aggregate flexibility at the ADC
using geometric optimization procedures.
III. BACKGROUND
A. Homothetic Transformation
In a recent work, [8], authors used a geometric optimizationbased approach to compute the inner and outer polytopic
approximation of the aggregated flexibility of an ensemble
of thermostatically-controlled loads, in the two-dimensional
space of control variable and active power consumption. We
propose to apply a similar approach to approximate the aggregated flexibility (Fi ) in active and reactive power consumption.
Before we explain the approach, let us first introduce the
concept of homothets.
Definition 1: [8] A homothet of a compact convex domain
F ⊂ Rn is defined as the family of domains which can be
expressed as H[α, β; F ] := α F + β, where α is a positive
scalar and β ∈ Rn . Henceforth, the scalar α will be referred
to as the ‘scaling factor’ and β as the ‘translational vector’.
Thus a homothetic transformation is a transformation applied on the state-space (Rn ) using uniform scaling and
translation. The following result is useful in computing the
inner and outer approximations of a set of homothets of a
given domain (extension of the following result for more than
two homothets follows trivially),
Theorem 1: [11] Consider a compact convex domain F ⊂
Rn , and two of its homothets H[αi , βi F ] , αi > 0 , βi ∈ Rn
∀i ∈ {1, 2} . The U
Minkowski sum of the homothets is given
2
by H[α, β; F ] = i H[αi , βi ; F ] , with the scaling factor as
α = α1 + α2 and the translational vector as β = β1 + β2 .
B. Sum-of-Squares Programming
Definition 2: Any n-variate polynomial that can be expressed as sum of squared polynomials, is called a sum-ofsquares (SOS) polynomial. We denote the ring of all SOS
polynomials in x ∈ Rn by Σ[x] .
Positivstellensatz theorem allows one to translate a set of
semi-algebraic constraints into SOS feasibility conditions. In
particular, Putinar’s Positivstellensatz theorem says:

q

q

q

q

p

(a) batteries

(b) switching loads

p

(c) PV inverters

p

(d) wind inverters

Figure 1: Representation of continuous and discrete flexibility domains for an individual DER of certain types.
Theorem 2: [12], [13] Consider a compact domain K =
{x ∈ Rn | gi (x) ≥ 0 ∀i = 1, . . . , m }, where gi (x) ∀i are polynomials and gi (·) ≥ 0 ∀i define compact domains. Then a
polynomial f (x) is positive on K if and only
P if there exist
SOS polynomials σi (x) such that f (x) − i σi (x)gi (x) is
SOS.
For every SOS polynomial f (x) of degree 2d (d is a positive
integer), there exists a positive semi-definite matrix Ξ such that
the Gramm matrix representation f (x) = z(x)T Ξ z(x) holds,
where z(x) is a vector of monomials in x of degree less than
or equal to d [14]. Thus each SOS problem can be cast into
an equivalent semidefinite programming (SDP) problem and
solved via SOSTOOLS [15], in conjunction with SDP solvers
such as SeDuMi [16]. Next we demonstrate the application of
SOS programming to compute the approximations of flexibility domains using prototypes.
IV. AGGREGATE F LEXIBILITY: G EOMETRIC A PPROACH
In this section we describe how one can use Theorem 1 to
design an algorithm to compute the outer and inner approximation of the aggregated flexibility without explicitly computing
their Minkowski sum. The idea is to first define a ‘prototype’
domain F 0 ⊂ R2 , either as a polytope (e.g. a box constraint)
or a generic convex domain, and obtain the outer and inner
approximations of the individual flexibility domains (Fi ⊂ R2 )
of the i-th ADC as the homothets of the prototype. For each
i ∈ {1, 2, . . . , N } , we seek to find positive scalars αi , αi ,
and vectors βi , βi ∈ R2 , such that




(1)
H αi , βi ; F 0 ⊆ Fi ⊆ H αi , βi ; F 0 .
Computation of the outer and inner approximations would be
discussed in details in the following sub-section. Note that,
Proposition 1: If the flexibility domains of a set of DERs
are outer and inner approximated by homothets of a compact
convex prototype domain F 0 ⊂ R2 as in (1), then their
aggregate flexibility domain F ⊂ R2 is outer and inner
approximated as follows,




H α, β; F 0 ⊆ F ⊆ H α, β; F 0
where α =

N
X
i=1

αi , β =

N
X
i=1

βi , α =

N
X

αi , β =

i=1

Proof The proof follows trivially from Theorem 1.

N
X
i=1

βi .

Let us now discuss how to construct the outer and inner
approximations of flexibility domains using prototypes. In this
paper, we consider the two types of flexibility domains: either
1) the flexibility region is a compact domain with a boundary
defined by piece-wise polynomials, or 2) the flexibility domain
is a collection of admissible points in the (p, q)-space. Note
that, we allow non-convex representations of the flexibility
domains. In this paper, we consider flexibility domains that
can be represented in the generic form of:

F=

K
[

k=1

Fk

(3)


k
where F k = (p, q) g1k (p, q) ≥ 0, . . . , gm
k (p, q) ≥ 0 , (4)
where gjk (p, q) ∀k ∈ {1, 2, . . . , mk } ∀i are polynomials. This
representation can be used to represent the continuous and
discrete flexibility domains of the types depicted in Fig. 1, as
explained in the following examples.
Example 1: (BATTERIES) The flexibility domain of a battery, with a maximum charge/discharge rate of pmax and the
apparent power rating of s > pmax , is given by
o
n
p
F = (p, q) p ∈ [−pmax , pmax ], |q| ≤ s2 − p2 ,
which could be expressed in the form of (3), with g1 (p, q) =
s2 −p2 −q 2 , and g2 (p, q) = (pmax )2 −p2 .
Example 2: (PV I NVERTERS) The flexibility domain of a
battery unit, with a maximum active power generation of pmax
and the apparent power rating of s > pmax , is given by
n
o
p
F = (p, q) p ∈ [−pmax , 0], |q| ≤ s2 −p2 ,
which could be expressed in the form of (3), with g1 (p, q) =
s2 −p2 −q 2 , and g2 (p, q) = −pmax p−p2 .
Example 3: (W IND I NVERTERS) The flexibility domain of
a wind inverter, with a maximum active power
√ generation of
pmax and the apparent power √
ratings s1 > αpmax (due to
rotor current limits) and s2 > αpmax (due to stator current

limits) , for some α > 0 is given by [17]–[19]
[
[
F = F1 F2 F3

F 1 = (p, q) p ∈ [−p0 , 0] , q ∈ [−q 0 , q 0 ]


q
2
max
0
2
2
F = (p, q) p ∈ [−p
, −p ), 0 ≤ q ≤ s2 − αp


q
F 2 = (p, q) p ∈ [−pmax , −p0 ), − s21 − αp2 ≤ q ≤ 0
where p0 and q 0 are much smaller than the rated capacities.
The flexibility can be expressed in the form of (3), with
g11 (p, q) = −p2 −p p0 , g21 (p, q) = (q 0 )2 −q 2 ,

g12 (p, q) = −p2 −p p0 +pmax −p0 pmax ,

g22 (p, q) = q, g32 (p, q) = s22 −αp2 −q 2 ,

g13 (p, q) = g12 (p, q), g23 (p, q) = −q, g33 (p, q) = s21 −αp2 −q 2 .
Example 4: (A IR -C ONDITIONERS) The flexibility domain
of a residential air-conditioner with an active power consumption rating of pmax (equal to the power consumed in ‘on’ state)
is represented by
[
F = F1 F2
F 1 = {(0, 0)} , F 2 = {(pmax , γpmax )}

where γ > 0 is related to the power factor. One possible way
of representing the flexibility in the form of (3) is by choosing
2
g11 (p, q) = −p2 , g21 (p, q) = −q 2 , g12 (p, q) = − (p − pmax ) and
2
max 2
g2 (p, q) = − (q − γp
) .
In the scope of this work, we will focus our attention to
prototype domains that are convex polygons, expressed as



 
a1p a1q
b1 

 


p
 ..

 .. 
0
F = (p, q) A
≤ b , A = .
,
b
=

 .  . (5)
q




anp anq
bn

Example 5: (U NIT S QUARE) A prototype domain of the
shape of a unit-square centered around the origin, F 0 =
{(p, q) | |p| ≤ 1, |q| ≤ 1 }, can be expressed in the form of (5)

T

T
1 −1 0 0
with A =
, b= 1 1 1 1 .
0 0 1 −1
A. Constructing Outer Approximation

An outer approximation of the compact flexibility domain in
(3), using a homothet of the prototype in (5) amounts to finding

T
positive scalar α and 2-dimensional vector β = βp βq
such that the following set inclusion condition is satisfied,

 
 0  



βp
p
p
0
0 0
0
=α 0 +
H α, β; F = (p, q)
, (p , q ) ∈ F
q
q
βq

 
 
β
p
(6)
= (p, q) A
≤ α b+A p ⊇ F .
q
βq
This condition translates into a set of semi-algebraic conditions
as follows:
∀i ∈ {1, . . . , n} : aip (βp −p) + aiq (βq −q)+α bi ≥ 0 on F .

After relaxing the inequality ‘≥’ to strictly inequality ‘>’ (note
that for convex domains, for any α satisfying the inequality
condition, we can always find an infinitesimally larger α
that satisfies the strict inequality), we can use Theorem 2 to
construct the following SOS optimization problem
minimize α > 0,

(7)

s.t.,
βp ∈ R, βq ∈ R,

k
σij
∈ Σ[p, q] ∀i ∈ {1, . . . , n}, j ∈ {1, . . . , mk }, k ∈ {1, . . . , K}



 βp −p

Pmk


a1p a1q
+α b1 − j=1 σ1j gjk ∈ Σ[p, q]


β
−q

q

..
∀k :
.







k
P

β
−p
m
p

 anp anq
+α bn − j=1 σnj gjk ∈ Σ[p, q]
βq −q

Note that the SOS conditions in the above SOS problem are
k
) and hence can
affine in the decision variables (α, βp , βq , σij
be solved directly. Note that the value of α represents how
large the homothet is, and therefore by solving for the minimal
α we ensure that the outer approximation is the tightest.
Remark 1: The problem formulated in (7) can be used
to compute the outer approximation for both the continuous
flexibility domains as well as the discrete flexibility domains
of the type shown in Fig.1(b) (and described in Example 4).
B. Constructing Inner Approximation

Solving for the inner approximation is not as straightforvector
ward. Given a scaling factor α and a translational

β , checking whether the homothet H α, β; F 0 is included
completely inside the flexibility domain F can be done by
solving the following feasibility problem
given α > 0, βp ∈ R, βq ∈ R,
(8a)
 

 
β
p
(8b)
check: ∀(k, i) gik ≥ 0 on (p, q) A
≤ α b+A p
q
βq
which can be cast into an SOS feasibility problem (for some
given α > 0 , βp ∈ R and βq ∈ R):

 i ∈ {1, . . . , mk },
k
j ∈ {1, . . . , n},
(9)
find σij ∈ Σ[p, q] ∀

k ∈ {1, . . . , K}

s.t.












 βp −p

k
a
a
+α bj ) ∈ Σ[p, q]
σ
(
jp
jq
j=1 1j
βq −q
..
∀k :
.







P
−p
β

p
n

k
k

a
a
+α bj ) ∈ Σ[p, q]
g
−
σ
(
jp
jq
 mk
j=1 mk j
βq −q
g1k −

Pn

Notice that the SOS constraints in (9) are bilinear in the
k
decision polynomial variables σij
and the homothet parameters
α and β. Therefore this cannot be solved directly to find out
the best inner approximation, and can only be used to check
whether a given homothet is an inner approximation or not.

In this paper, we propose a heuristic to find the homothet
parameters that define the best inner approximation. Note that
given a β the best inner approximation is the largest homothet
that is contained inside the flexibility domain, which can be
solved by running a bisection search for the maximum α for
which (9) is feasible. The value of this largest α can be further
improved if we can identify a direction in which to translate
the homothet, by detecting the edges of the polytope which
touch the flexibility domain boundary (i.e. binding edges). The
best inner approximation is computed in a two-stage iterative
process as described below.
1) Step 1: Find α given β . Given a β we can run a bisection
search for the largest positive scalar α that solves (9).
2) Step 2: Update β . Once we have found the largest
homothet H α, β; F 0 of the prototype polygon, that is contained wholly inside the flexibility domain, we can identify
the edges of the homothet which lie on the boundary of the
flexibility domain. The edges represent the binding constraints
that limit the expansion of the inner homothet. Thus we
choose a direction b̂ that is a vector sum of all the directions
orthogonal to the binding edges, and choose a new b ← b+εb̂ ,
for some sufficiently small scalar ε . Checking which edges are
binding can be done by solving appropriate sum-of-squares
problems, which we have omitted for brevity.
These two steps are repeated until a convergence
in the

value of α is achieved. Let us denote by H α∗ , β ∗ ; F 0 the
largest homothet contained inside the flexibility domain,
such


that there exists no α > α∗ and β for which H α, β; F 0  is
an inner approximation. Let us assume that H α0 , β 0 ; F 0 is
an inner approximation for some translational vector β 0 and
scaling factor α0 . Then the success (convergence) of the above
two-stage algorithm depends on the following condition:
Assumption 1: α(ǫ) increases monotonically with ǫ ∈ [0, 1] ,
where α(ǫ) is the largest scaling factor corresponding
to the

translational vector β(ǫ) = β 0 + ǫ β ∗ − β 0 .
Note that the inner approximations do not exist for DERs
that have discrete flexibility domains. Hence the above method
is only applicable to continuous flexibility domains.
V. Q UALITY

OF

A PPROXIMATION

In order for the approximation of the flexibility domains to
be used in network-level resource optimization, it is important
that the approximations are close to the actual flexibility
domains. In this section, we propose metric(s) to quantify
the closeness of the approximation. Such metrics could be
useful to choose from a set of prototypes to approximate
the flexibility domain of a given DER. Furthermore, such
metrics could be passed up to the network optimizer, along
with the approximating homothets, to facilitate better resource
allocation. In this paper, we discuss two such possible metrics.
First one is related to the worst case error in a dispatched point
(from network operator) and a feasible solution:
Definition 3: The distance metric π d (F , F ) for two sets
F ⊃ F is defined as the maximal distance between any point
in F \F and its nearest feasible point in F .

Proposition 2: The distance metric for
 the outer
 and inner
0
approximations
of
the
type
H
=
H
α,
β;
F
and H =


H α, β; F 0 , where F 0 is a polygon, is given by
π d (H, H) = max (α − α) vi + β − β
i

2

where vi ∀i is the vector representing the vertices of F 0 .
Proof Since homothetic transformation include only scaling
and translation, the edges in both the inner and outer homothets are parallel to each other. Therefore the maximal distance
between the two polygons occur at the respective vertices.
The second metric is related to the likelihood of a dispatched
point lying inside the actual flexibility domain:
Definition 4: The area metric π a (F , F ) for the outer approximated and actual flexibility domains, F ⊃ F , is defined
as the likelihood that an optimally dispatched point in F also
lies inside F .
To compute this metric one needs to assign a probability
of each point in the approximated domains being dispatched
by the network operator. In absence of any prior knowledge
of the probability, one can start by assuming that each point
in the approximated domain is equally likely to be dispatched.
Under such an assumption, we can compute the area metric
as follows:
Proposition 3: The area metric for the outer
 and inner
0
H
=
H
α,
β;
F
and H =
approximations
of
the
type


2
H α, β; F 0 , is given by π a (H, H) = (α/α) .


Proof Note that the area under a homothet H α, β; F 0 is
equal to α2 times the area under the prototype F 0 . The rest
follows directly.
Note that the larger the value of the area metric, the better
are the domain approximations, while for the distance metric
the smaller value implies better approximation. These metrics
are applicable directly to continuous flexibility domains. For
discrete domains, alternative methods need to be devised to
compute the metrics (not covered in this article).
VI. N UMERICAL R ESULTS
In this section, we present some numerical results to demonstrate the aggregation of flexibility domains of different types
of DERs. Fig. 2 shows the plots of outer and inner approximations of a collection of five DERs (consisting of PV inverters,
wind inverters and batteries), using unit square as a prototype,
including the approximation of their Minkowski sum. The
values of the distance metric and the area metric for all these
approximations are also computed and shown on the plots.
Fig. 3 shows the approximations using a regular unit hexagon
as the prototype, for the same set of devices (only a couple
of them shown on the figure). In this particular example, the
approximations using a regular hexagon as the prototype turns
out to be better than the approximations using unit square
as the prototype, according to both the area metric and the
distance metric. For example, the value of the area metric of
the approximation of the aggregated flexibility turns out to be

d

=0.43

4

reactive power consumption, q [kVAR]

6
4
2
0
-2
-4
-6

-4

-2

0

2

4

6

d

=0.18

15
actual
outer
inner

2
1
0
-1
-2
-3

-4

-2

0

2

a

=0.14
actual
outer
inner

10

5

0

-5

-10

-15
-20

4

=10.39,

-15

-10

-5

0

5

active power consumption, p [kW]

active power consumption, p [kW]

active power consumption, p [kW]

(a) battery (DER 1)

(b) solar PV inverter (DER 2)

(c) wind inverter (DER 3)

d

10

=3.62,

a

d

=0.50

4
actual
outer
inner

8
6
4
2
0
-2
-4
-6
-8
-10
-10

a

3

-4
-6

8

=3.06,

-5

0

5

d

=0.12

40
actual
outer
inner

3
2
1
0
-1
-2
-3
-4
-8

10

=3.68,

a

reactive power consumption, q [kVAR]

-6

reactive power consumption, q [kVAR]

reactive power consumption, q [kVAR]

a

actual
outer
inner

-8
-8

reactive power consumption, q [kVAR]

=3.00,

reactive power consumption, q [kVAR]

d

8

-6

-4

-2

0

a

=0.26
actual
outer
inner

30
20
10
0
-10
-20
-30
-40
-50

2

=23.73,

10

-40

-30

-20

-10

0

10

active power consumption, p [kW]

active power consumption, p [kW]

active power consumption, p [kW]

(d) battery (DER 4)

(e) wind inverter (DER 5)

(f) aggregated flexibility

20

30

Figure 2: Outer and inner approximations of a collection of PV inverters, wind inverters and batteries, using unit square as a
prototype.
a

d

=0.26

5
actual
outer
inner

3

reactive power consumption, q [kVAR]

reactive power consumption, q [kVAR]

=3.09,

2
1
0
-1
-2
-3
-4
-4

-3

-2

-1

0

1

2

3

=3.39,

a

d

=0.15

40
actual
outer
inner

4
3

reactive power consumption, q [kVAR]

d

4

2
1
0
-1
-2
-3
-4
-5
-6

-4

-2

0

2

4

=21.41,

a

=0.30
actual
outer
inner

30
20
10
0
-10
-20
-30
-40
-40

-30

-20

-10

0

10

active power consumption, p [kW]

active power consumption, p [kW]

active power consumption, p [kW]

(a) solar PV inverter (DER 2)

(b) wind inverter (DER 5)

(c) aggregated flexibility

20

30

Figure 3: Outer and inner approximations of the same set of DERs, using regular unit hexagon as a prototype.
higher for the regular hexagon prototype (π a = 0.30) than
for the square prototype (π a = 0.26), while the value of the
distance metric is lower with regular hexagon as the prototype
(π d ≈ 21.4) than with a square prototype (π d ≈ 23.7).
Finally Fig. 4 shows the results from approximating the
aggregated flexibility of a collection of discrete loads (residential air-conditioners) via a triangle as a prototype. While
the approximation of an individual DER results in a very
sparse approximated flexibility domain, the approximation of
the aggregated flexibility domain tends to be much denser.
While a specific quantification of the accuracy (or, quality) of
the approximation for discrete DERs is not provided in this
article, Fig. 4 shows qualitatively how the larger population
size of an ensemble of switching loads can lead to better
approximation of the aggregated flexibility.

VII. C ONCLUSIONS
We envision a distribution system optimization framework
in which thousands of these DERs are scheduled to operate the
system efficiently. Heterogeneity of the DERs, manifested in
the different forms of their flexibility domains, pose a serious
challenge to the scalability of such an optimization problem.
Furthermore, any such optimization framework needs to be
able to accommodate a plug-and-play approach whereby any
new DER can be easily integrated into the system. In this
work, we presented a geometric approach to homogenizing the
flexibility domains of different types of DERs via homothets
of some prototype domain. We further looked into the choice
of the prototype domain, and the choice of the parameters of
the homothets, such that the outer (and inner, when applicable)
approximations of the flexibility domains are good. Simulation
results are provided to illustrate the approach. Future efforts

reactive power consumption, q [kVAR]

2
actual
outer

1.8
1.6
1.4
1.2
1
0.8
0.6
0.4
0.2
0
0

1

2

3

4

5

6

7

8

active power consumption, p [kW]

(a) individual air-conditioner

(b) collection of air-conditioners

Figure 4: Outer approximations of a collection of ten residential air-conditioners with discrete flexibility domains, using a
triangle as a prototype.

will concentrate on extending this framework to model the
dynamic evolution of the flexibility, as well as modeling the
uncertainties in the flexibility of DERs.
ACKNOWLEDGMENT
The authors would like to thank Dr. Andrey Bernstein of
National Renewable Energy Laboratory, USA, for his helpful
remarks and discussions on this work.
R EFERENCES
[1] D. S. Callaway, “Tapping the energy storage potential in electric loads to
deliver load following and regulation, with application to wind energy,”
Energy Conversion and Management, vol. 50, no. 5, pp. 1389–1400,
2009.
[2] S. Kundu and I. A. Hiskens, “Hysteresis-based charging control of plugin electric vehicles,” in Proceedings of the IEEE 51st Annual Conference
on Decision and Control, 2012, pp. 5598–5604.
[3] C. Perfumo, E. Kofman, J. H. Braslavsky, and J. K. Ward, “Load
management: Model-based control of aggregate power for populations of
thermostatcally controlled loads,” Energy Conversion and Management,
vol. 55, pp. 36–48, Mar. 2012.
[4] J. L. Mathieu, S. Koch, and D. S. Callaway, “State estimation and
control of electric loads to manage real-time energy imbalance,” IEEE
Transactions on Power Systems, vol. 28, no. 1, pp. 430–440, Feb 2013.
[5] W. Zhang, J. Lian, C. Y. Chang, and K. Kalsi, “Aggregated modeling
and control of air conditioning loads for demand response,” IEEE
Transactions on Power Systems, vol. 28, no. 4, pp. 4655–4664, Nov
2013.
[6] 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, 2015.
[7] B. M. Sanandaji, T. L. Vincent, and K. Poolla, “Ramping Rate Flexibility
of Residential HVAC Loads,” IEEE Transactions on Sustainable Energy,
vol. 7, no. 2, pp. 865–874, April 2016.

[8] L. Zhao, W. Zhang, H. Hao, and K. Kalsi, “A geometric approach
to aggregate flexibility modeling of thermostatically controlled loads,”
arXiv preprint arXiv:1608.04422, 2016.
[9] D. S. Callaway and I. A. Hiskens, “Achieving controllability of electric
loads,” Proceedings of the IEEE, vol. 99, no. 1, pp. 184–199, 2011.
[10] A. Bernstein, L. Reyes-Chamorro, J.-Y. Le Boudec, and M. Paolone,
“A composable method for real-time control of active distribution
networks with explicit power setpoints. part i: Framework,” Electric
Power Systems Research, vol. 125, pp. 254–264, 2015.
[11] R. Schneider, Convex bodies: the Brunn–Minkowski theory. Cambridge
University Press, 2013, no. 151.
[12] M. Putinar, “Positive polynomials on compact semi-algebraic sets,”
Indiana University Mathematics Journal, vol. 42, no. 3, pp. 969–984,
1993.
[13] J.-B. Lasserre, Moments, Positive Polynomials and Their Applications.
World Scientific, 2009, vol. 1.
[14] P. A. Parrilo, “Structured semidefinite programs and semialgebraic
geometry methods in robustness and optimization,” Ph.D. dissertation,
Caltech, Pasadena, CA, 2000.
[15] A. Papachristodoulou, J. Anderson, G. Valmorbida, S. Prajna,
P. Seiler, and P. A. Parrilo, “SOSTOOLS: Sum of squares
optimization toolbox for MATLAB,” 2013, available from
http://www.eng.ox.ac.uk/control/sostools.
[16] J. F. Sturm, “Using SeDuMi 1.02, a MATLAB toolbox for optimization over symmetric cones,” Optimization Methods and Software, vol. 11-12, pp. 625–653, Dec. 1999, software available at
http://fewcal.kub.nl/sturm/software/sedumi.html.
[17] T. Lund, P. Sørensen, and J. Eek, “Reactive power capability of a wind
turbine with doubly fed induction generator,” Wind energy, vol. 10, no. 4,
pp. 379–394, 2007.
[18] J. Tian, C. Su, and Z. Chen, “Reactive power capability of the wind
turbine with doubly fed induction generator,” in Industrial Electronics
Society, IECON 2013 - 39th Annual Conference of the IEEE, Nov 2013,
pp. 5312–5317.
[19] J. A. Martin and I. A. Hiskens, “Reactive power limitation due to windfarm collector networks,” in PowerTech, 2015 IEEE Eindhoven, June
2015, pp. 1–6.

