Aggregation and Disaggregation of Energetic Flexibility from Distributed Energy Resources

TL;DR

This paper presents a scalable method using zonotopic sets to quantify, aggregate, and distribute flexibility from distributed energy resources, enabling efficient control and market participation.

Contribution

It introduces a generic, scalable approach for describing and pricing flexibility of distributed resources using zonotopic sets, including algorithms for fair control distribution.

Findings

Efficient aggregation of flexibility for large numbers of resources.

Algorithm for fair and computationally efficient control distribution.

Method for computing aggregate regulation power bid-curves.

Abstract

A variety of energy resources has been identified as being flexible in their electric energy consumption or generation. This energetic flexibility can be used for various purposes such as minimizing energy procurement costs or providing ancillary services to power grids. To fully leverage the flexibility available from distributed small-scale resources, their flexibility must be quantified and aggregated. This paper introduces a generic and scalable approach for flexible energy systems to quantitatively describe and price their flexibility based on zonotopic sets. The description proposed allows aggregators to efficiently pool the flexibility of large numbers of systems and to make control and market decisions on the aggregate level. In addition, an algorithm is presented that distributes aggregate-level control decisions among the individual systems of the pool in an economically fair and computationally efficient way. Finally, it is shown how the zonotopic description of flexibility enables an efficient computation of aggregate regulation power bid-curves.

Tables

Table 1. TABLE I : PEV parameter ranges.

Variable	Description	Value range
$H$	Planning horizon	24 h
$\overline{p}$	Maximum rated power	3 kW
$\underset{¯}{p}$	Minimum rated power	-3 kW
$\overline{e}-\underset{¯}{e}$	Battery capacity	[20,40] kWh
$\left(e_0-\underset{¯}{e}\right)/\left(\overline{e}-\underset{¯}{e}\right)$	Initial state-of-charge	[0.2,0.8]

Table 2. TABLE II : Solving time of LP-based disaggregation ( V-A ).

$J$ =	$N$ = 24	48	96	168	288	672
${10}^3$	2.9 s	15.9 s	34.0 s	88.9 s	871.7 s	M
$2\cdot {10}^3$	5.9 s	45.7 s	145.4 s	670.5 s	1603.7 s	M
$5\cdot {10}^3$	26.3 s	155.9 s	931.5 s	2647.4 s	M	M
${10}^4$	99.8 s	449.9 s	3811.1 s	7319.2 s	M	M

Table 3. TABLE III : Solving time of subgradient-based disaggregation.

$J$ =	$N$=24	48	96	168	288	672
${10}^3$	0.13 s	0.13 s	0.17 s	0.24 s	0.32 s	0.70 s
$2\cdot {10}^3$	0.12 s	0.15 s	0.22 s	0.29 s	0.42 s	0.88 s
$5\cdot {10}^3$	0.15 s	0.20 s	0.27 s	0.47 s	0.76 s	1.71 s
${10}^4$	0.19 s	0.29 s	0.48 s	0.82 s	1.30 s	2.97 s

Table 4. TABLE IV : Baseline costs for different types of feasible sets.

	$\overline{r}=0$ kW	20 kW	40 kW	60 kW
Polytopes	-77.9 EUR	-59.5 EUR	-41.0 EUR	-22.2 EUR
Zonotopes	-77.7 EUR	-58.6 EUR	-39.1 EUR	-18.5 EUR
A-A Boxes	-43.1 EUR	-25.8 EUR	-8.6 EUR	8.7 EUR

Equations

$${\underline{p}}_k$$

$${\underline{e}}_k-e_0\le t_s\sum _{i=1}^k{p}_i\le {\bar{e}}_k-e_0,k=1,\dots ,N,$$

$${\underline{r}}_k\le (p_k-p_{k-1})/t_s\le {\bar{r}}_k,k=2,\dots ,N.$$

$${\underline{x}}_{k+1}\le A_k{x}_k+B_k{u}_k+C_k{p}_k\le {\bar{x}}_{k+1},k=1,\dots ,N-1,$$

$$Z=\{x\in {\mathbb{R}}^N:x=c+G\beta ,-\bar{\beta }\le \beta \le \bar{\beta }\},$$

$$g^{(i)}$$

$$g^{(N+j)}$$

$$\Lambda :=\frac{1}{q}\sum _{i=1}^q\frac{{\Delta }_{Z,i}}{{\Delta }_{P,i}}\in [0,1],$$

$${\mathcal{P}}^{(\text{agg})}:=$$

$$=$$

$$\displaystyle\begin{split}Z^{(\text{agg})}:=\,&Z(G,c^{(\text{agg})},\bar{\beta}^{(\text{agg})})\\ =\,&Z\bigg{(}G,\sum\limits_{j\in\mathcal{J}}{c^{(j)}},\sum\limits_{j\in\mathcal{J}}{\bar{\beta}^{(j)}}\bigg{)}.\end{split}$$

$$T^{(j)}(p^{(j)}):=f^{(j)}(p^{(j)})-t_s{v}^{(j{)}^{\top }}{p}^{(j)}.$$

$$T^{(j)}({\beta }^{(j)})$$

$$=\sum _{i=1}^g{T}_i^{(j)}({\beta }_i^{(j)})+T_{\text{fix}}^{(j)},$$

$$T^{(\text{agg})}({\beta }^{(\text{agg})})=\sum _{i=1}^g{T}_i^{(\text{agg})}({\beta }_i^{(\text{agg})})+T_{\text{fix}}^{(\text{agg})},$$

$$T_i^{(\text{agg})}({\bar{\beta }}_i^{(\text{agg})})=\sum _{j\in \mathcal{J}}{T}_i^{(j)}({\bar{\beta }}_i^{(j)})+T_{\text{fix}}^{(\text{agg})}.$$

$$\hat{R}(\bar{r}):=\underset{\{p^{(j)},r^{(j)}{\}}_{j\in \mathcal{J}}}{\min }\sum _{j\in \mathcal{J}}\left\{T^{(j)}(p^{(j)})+{\hat{v}}^{\top }{p}^{(j)}\right\},$$

$${\bar{r}}_{\max }=$$

$$G^{PE}{\beta }^{(j)}\ge 0,j\in \mathcal{J},$$

$$G^{PE}\sum _{j\in \mathcal{J}}{\beta }^{(j)}\ge \bar{r}1,$$

$$-{\bar{\beta }}^{(j)}\le {\beta }^{(j)}\le {\bar{\beta }}^{(j)},j\in \mathcal{J}.$$

$$\underset{{\beta }^{(\text{agg})}}{\min }$$

$$-\sum _{j\in \mathcal{J}}{\beta }_{\text{rem}}^{(j)}\le {\beta }^{(\text{agg})}\le \sum _{j\in \mathcal{J}}{\beta }_{\text{rem}}^{(j)}.$$

Full text

I Introduction

Increasing shares of distributed renewable energy resources lead to a growing need for energetic flexibility in power grids. On the one hand, renewable energy sources, such as wind and solar power, introduce large amounts of variability and uncertainty to power grids because of their inherent intermittency and limited controllability. On the other hand, substituting well-controllable traditional generating units, such as nuclear, gas, and pumped-hydro storage power plants, by renewables reduces the rotational inertia of the power system and the supply of highly reliable regulation power [1]. Consequently, future power grids will require additional sources of flexibility to guarantee safe and reliable operation.

The flexibility of distributed, small-scale energy systems, on both the demand and the supply side, can serve as an additional resource to the large-scale supply flexibility that has been traditionally used to balance the grid. Sources of small-scale flexibility include heating, ventilation, and air-conditioning systems, plug-in electric vehicles and stationary batteries, as well as micro-generation units. However, quantifying and making use of the flexibility of large numbers of distributed systems is challenging, in particular because of the computational complexity. To take full advantage of the flexibility available and to reduce the complexity of planning, trading and control, the aggregate flexibility of an entire group of systems must be computed and represented in a concise and compact form. The importance of aggregation as a key enabler for incorporating large numbers of flexible systems has been highlighted in numerous works, e.g. , [2, 3, 4, 5, 6]. The tasks of collecting, aggregating and controlling the flexibility of a group of systems are performed by an entity called aggregator [7, 8]. The aggregator concludes contracts with individual flexible energy resources, defining the way of communication and control between the systems and the aggregator, the specifics of the flexibility offered by the systems, as well as the details on how the systems are reimbursed for the flexibility provided [7, 2, 3, 9]. Thus, the aggregator acts as an intermediary between flexible resources and wholesale markets, molding individual flexibilities into tradable products.

Various approaches have been proposed to characterize the energetic flexibility both of individual systems and of entire populations. We consider flexibility from the perspective of the power grid and define it as the ability of a system to adjust the volume and/or timing of its electric energy intake from or of its energy output to the power grid. Consequently, the set of all power trajectories that can be tracked by the system over a given time horizon is a natural description of its flexibility and is referred to as the feasible set [4, 10]. Among the approaches for characterizing the feasible set of individual systems are power nodes [11, 12], generalized batteries [13, 5, 14], and resource polytopes [4, 10, 6].

Regarding the description of the flexibility of an entire population of systems, two main types of approaches can be identified. Top-down approaches attempt to capture the aggregate flexibility directly based, for example, on the probabilistic properties of the underlying systems [15, 16, 17]. In contrast, here we consider a bottom-up approach in which the flexibility of each individual system is described first and aggregated subsequently. The particular description of individual feasible sets determines the computational complexity of aggregating them. In general, computing the exact aggregate feasible set is computationally intractable [4, 6]. Consequently, attempts to compute the exact aggregate feasible set are restricted to special cases, such as symmetric constraints [13] or binary resource polytopes [4]. Other approaches use set descriptions that provide inner- and/or outer-approximations of the exact aggregate feasible set. Outer-approximations are derived based on polytopes [8, 18] or semi-definite constraints [19]. While some outer-approximations can be computed efficiently, the fact that they can contain infeasible power trajectories is a major drawback. In contrast, inner-approximations provide sufficient conditions for feasibility. Approaches for computing inner-approximations are based on zonotopics sets [10] or on projection to polytopic sets [6]. Both inner- and outer-approximations for battery models are provided in [14, 13].

Compared to flexibility aggregation, disaggregation, i.e. , distributing an aggregate control signal among the systems of a population, has attracted limited attention. Hierarchical control schemes have been proposed in [20, 4]. Other approaches use a priority-stack [13] or randomized [16] controllers.

The main contribution of this work is to introduce zonotopes, a subclass of polytopes, as approximations of the original feasible set of a flexible system. Zonotopes can be aggregated efficiently with regard to memory requirements and computational complexity and can capture time-variant, asymmetric power, energy and ramp constraints. In addition, an efficient disaggregation algorithm based on zonotopes is presented that distributes a given aggregate-level signal among the systems in an economically fair way. Moreover, an aggregator can use the algorithm to compute regulation power bids.

Section II introduces the idea of using zonotopes to approximate the flexibility of individual systems. Optimal zonotopic approximations are computed in Sec. III. Sec. IV discusses the aggregation of flexibility, followed by Sec. V where a disaggregation algorithm based on flexibility costs is derived. In Sec. VI, zonotopic descriptions of flexibility are used to perform feasibility checks and to compute bid curves for regulation power. A conclusion and outlook are provided in Sec. VII.

II Description of flexibility

II-A Polytopic feasible sets

Consider a single energy resource over a finite discrete-time horizon comprising $N$ time steps each of duration $t_{s}$. By definition, the feasible set ${\mathcal{P}\subseteq\mathbb{R}^{N}}$ comprises all the power trajectories ${p\in\mathbb{R}^{N}}$ the system is able to follow. This set is defined by the dynamics and constraints the system must respect. The constraints most common to flexible energy resources are the following.

1) Power constraints: In general, the constant power $p_{k}$ the system draws from (${p_{k}>0}$) or feeds into (${p_{k}<0}$) the power grid during time step $k$, is bounded, i.e. ,

[TABLE]

2) Energy constraints: Most flexible energy resources owe their flexibility to an energy buffer with limited capacity. Thus, the cumulative electric energy $e_{k}$ the system can consume or generate up to and including time slot $k$ is bounded, i.e. ,

[TABLE]

where $e_{0}$ is the initial energy level of the system.

3) Ramp-rate constraints: Often, there are limits on the rate at which power can vary, i.e. ,

[TABLE]

4) State constraints: Various energy resources have internal states that are driven by the electric power consumption or generation of the system. Examples are the state-of-charge of a battery or the temperature of a heating system. Assuming linear state dynamics, the state constraints can be written as

[TABLE]

where $u_{k}$ summarizes all inputs other than power.

The set of all power trajectories ${p:=[p_{1},\dots,p_{N}]^{\top}}$ that satisfy constraints (1)-(4) is given by the convex polytope, referred to as the resource polytope in [4], ${P:=\{p\in\mathbb{R}^{N}:Ap\leq b\}}$, where $A,\,b$ summarize the constraint matrices and the limits of (1)-(4), respectively.

Many energy systems are abstracted as idealized energy buffers subjected to power and energy constraints (1) and (2) only, see [18, 6], among others. Such systems are referred to as PE-systems, and their feasible sets are called PE-polytopes.

II-B Zonotopic feasible sets

The description of flexibility by polytopes poses serious computational challenges when making decisions for an entire group of systems, cf. Sec. IV. To alleviate this difficulty we propose to inner-approximate individual polytopic feasible sets by a subclass of polytopes, known as zonotopes, that can be aggregated efficiently. Zonotopes are centrally symmetric objects. This allows a zonotope $Z$ to be expressed in terms of its center $c\in\mathbb{R}^{N}$ and generators ${g^{(i)}\in\mathbb{R}^{N}}$, ${\|g^{(i)}\|_{2}=1}$, ${i=1,\dots,g,}$ as

[TABLE]

where the generators are summarized in the generator matrix ${G:=[g^{(1)},\dots,g^{(g)}]\in\mathbb{R}^{N\times g}}$. The scaling factor ${\beta\in\mathbb{R}^{g}}$ is symmetrically bounded by $\bar{\beta}$. The shorthand $Z(G,c,\bar{\beta})$ is used to denote such a zonotope. We propose a particular generator matrix such that the associated family of zonotopes is well suited to approximate PE-polytopes. For arbitrary dimension ${N\in\mathbb{N}}$, a total of ${2N-1}$ generators are constructed as

[TABLE]

with ${i=1,\dots,N}$ and ${j=1,\dots,N-1}$. The corresponding generator matrix is ${G^{PE}\in\mathbb{R}^{N\times(2N-1)}}$. Figure 1 provides a two-dimensional example of such a zonotope. It is shown in Appendix A that zonotopes with generator matrix $G^{PE}$ can reconstruct every possible facet of a PE-polytope. However, the approximation quality is influenced by the fact that zonotopes are centrally symmetric by construction, and, in general, will feature additional facets that do not exist in a PE-polytope.

III Computing zonotopic approximations

III-A Computing optimal zonotopes

In this section, it is shown how to compute zonotopes that inner-approximate a given polytope and that are optimal with regard to specific objectives. A key ingredient to all the optimization problems formulated below is a set of constraints that enforce that a zonotope is included in a polytope: ${Z(G,c,\bar{\beta})\subseteq P(A,b)\ \Leftrightarrow\ Ac+|AG|\bar{\beta}\leq b}$, cf. [10] for more details. Given a resource polytope $P(A,b)$ and a generator matrix $G$, different objectives can be used to compute zonotopes that inner-approximate the polytope [10]. The choice of objective depends on the purpose for which the aggregate feasible set will be used. Here we consider the objective illustrated in Fig. 1. The distances between parallel zonotope facets are given as ${\Delta_{Z}=2|F^{\top}G|\bar{\beta}}$, where ${F\in\mathbb{R}^{N\times q}}$ is a matrix that contains the normal vectors of all possible zonotope facets as columns, i.e. , ${F:=[f^{(1)},\dots,f^{(q)}]^{\top}\in\mathbb{R}^{N\times q}}$, with ${||f^{(i)}||_{2}=1}$, ${i=1,\dots,q}$. Zonotopes with generator matrix $G^{PE}$ can have at most ${q=N^{2}+N}$ pairs of parallel facets, cf. Appendix A. Similarly, the extension of the polytope in the direction of the $i^{\text{th}}$ zonotope facet normal is denoted by $\Delta_{P,i}$, ${i=1,\dots,q}$, and computed in a preprocessing step by solving $q$ Linear Programs (LP). Let ${\Delta_{P}:=[\Delta_{P,1},\dots,\Delta_{P,q}]^{\top}}$. The zonotope approximates the polytope perfectly if and only if ${\Delta_{Z}=\Delta_{P}}$. Thus, a convenient and informative approximation quality measure is

[TABLE]

where ${\Lambda=0}$ implies that the zonotope is a point, and ${\Lambda=1}$ is reached iff the polytope is approximated perfectly by the zonotope. The inner-approximating zonotope maximizing (8) can be found by solving the following LP: {maxi!} {c,¯β}w^⊤¯β \addConstraintAc+—AG—¯β≤b \addConstraint¯β≥0,

with ${w:=(2/q)(1/\Delta_{P})^{\top}|F^{\top}G|}$. The objective (8) expresses (8) as a linear function of $\bar{\beta}$. Constraints (8)-(8) guarantee that the zonotope is included in the polytope.

Problem (8) may result in degenerate zonotopes with ${\Delta_{Z,i}=0}$ for some ${i\in\{1,\dots,q\}}$ even though ${\Delta_{P,i}>0}$. This is undesirable because such zonotopes lack a significant portion of the polytope flexibility. To prevent such degeneracy, (8) can be augmented by constraints enforcing the zonotope to contain a full-dimensional object, e.g. , a box. Let $C_{\max}$ denote the maximum-volume axis-aligned cube in $P$. Further, let $B_{\max}$ be the box with maximum cumulative edge length from the set of axis-aligned boxes ${\{B:C_{\max}\subseteq B\subseteq P\}}$. We augment (8) by convex constraints enforcing that ${B_{\max}\subseteq Z(G,c,\bar{\beta})}$.

III-B Approximation results

How well a zonotope can approximate a given polytope depends on a number of factors, including the choice of generator matrix, the shape of the polytope (in particular its symmetry properties), and, in case the approximation is based on optimization, the objective function chosen. Moreover, the number of possible facets of the zonotope and the polytope grows with the number of dimensions $N$ and also affects the approximation quality. Here, the approximation quality of zonotopes is assessed for resource polytopes of plug-in electric vehicles (PEVs) over a planning horizon of 24 h with 2 h sampling time. It is assumed that the PEVs are subject to power and energy constraints only. Trip constraints can easily be translated into power constraints by setting both the upper and lower power bounds to zero for all time slots that fall into trip periods. Table I provides the parameter ranges from which the parameters of 100 PEVs are sampled uniformly. For every system, an optimal inner-approximating axis-aligned box $B^{\ast}$ and zonotope $Z^{\ast}$ are computed via (8) using ${G=I_{N}}$ and ${G=G^{PE}}$, respectively. Averaged over the population, the quality of approximation (8) for the box and the zonotope are ${\Lambda_{B}=0.31}$ and ${\Lambda_{Z}=0.63}$, respectively.

IV Aggregation of flexibility

IV-A Aggregation setup

To take full advantage of the flexibility available from entire groups of flexible energy resources, their individual flexibilities are pooled by aggregators. Aggregators act as intermediaries between flexible systems and wholesale markets. Consider an aggregator that manages the electric energy consumption/production of a population of systems indexed by ${j\in\mathcal{J}:=\{1,\dots,J\}}$. Every system provides the aggregator with a description of its feasible set ${\mathcal{P}^{(j)}\subseteq\mathbb{R}^{N}}$ for a given planning horizon comprising $N$ time steps. The feasible set of a system is either known directly from the physical properties and constraints of the system or has to be estimated from measurement data [21]. The aggregator decides how to optimally use the flexibility available from all the systems in the population. The aggregator might trade electricity on wholesale energy markets or offer dedicated services to flexibility markets, e.g. , regulation power. Finally, the aggregator is required to assign to every system a reference schedule ${p^{(j)}\in\mathcal{P}^{(j)}}$, ${j\in\mathcal{J}}$, that the system is committed to follow. This assignment is done either ahead of time, e.g. , in the case of a day-ahead reference schedule, or during time of delivery, e.g. , upon the activation of regulating power. Given the electricity consumption/production of an individual system during the period of delivery, money flows either from the aggregator to the system or vice versa according to the terms agreed upon by the system and the aggregator.

IV-B Aggregation of feasible sets

A power trajectory ${p^{(\text{agg})}\in\mathbb{R}^{N}}$ that can be followed collectively by the population of $J$ systems is called aggregate feasible. The set of aggregate feasible power trajectories is denoted by $\mathcal{P}^{(\text{agg})}$ and referred to as the aggregate feasible set. It is given by the Minkowski sum of all individual feasible sets [13, 19]:

[TABLE]

The next sections look into computing $\mathcal{P}^{(\text{agg})}$ for the special cases where the individual feasible sets are given either by general polytopes or by zonotopes.

1) Aggregation of polytopic flexibility: When the feasible sets $\mathcal{P}^{(j)}$ are given by polytopes the aggregate feasible set (9) is also a polytope [22]. Different methods to compute the Minkowski sum of polytopes exist, for example by computing the convex hull of the pairwise sum of all vertices of the polytopes or by projecting an augmented polytope to a lower-dimensional subspace [23, 22, 6]. However, no efficient algorithm exists to compute the Minkowski sum of general polytopes in hyper-plane representation [23, 22]. Consequently, attempts to compute the Minkowski sum of resource polytopes are restricted to special types of polytopes (symmetric, binary hyperplane normal vectors) [13, 20], are limited to low dimensions $N$, or compute approximations of the true sum [6, 14].

2) Aggregation of zonotopic flexibility: The Minkowski sum of zonotopes can be computed efficiently. The aggregate feasible set $Z^{(\text{agg})}$ of the individual feasible sets $Z(G,c^{(j)},\bar{\beta}^{(j)})$ is given explicitly by

[TABLE]

The structure of (10) is particularly convenient for aggregators. First, $Z^{(\text{agg})}$ is a zonotope because zonotopes are closed under Minkowski addition. Second, if systems are added to or removed from the population or if individual feasible sets change, the aggregate feasible set can be adjusted efficiently by updating the sums $c^{(\text{agg})}$ and $\bar{\beta}^{(\text{agg})}$ accordingly.

V Disaggregation

V-A Problem formulation

In the aggregation scheme introduced above, the aggregator is responsible for trading a sufficient amount of electric energy such that every system can follow a feasible power trajectory ${p^{(j)}\in\mathcal{P}^{(j)}}$, ${j\in\mathcal{J}}$. The contract between the aggregator and each individual system $j$ defines the price $v_{k}^{(j)}$ at which the system buys/sells electric energy from/to the aggregator for every time step ${k=1,\dots,N}$. The aggregator assigns a reference trajectory $p^{(j)}$ to every system, for which the system pays the amount ${t_{s}v^{(j)^{\top}}p^{(j)}}$, ${v^{(j)}:=[v^{(j)}_{1},\dots,v^{(j)}_{N}]^{\top}}$, to the aggregator. While the systems guarantee that every ${p^{(j)}\in\mathcal{P}^{(j)}}$ can be implemented, they might prefer certain trajectories over others. For example, the owner of a building might prefer those power trajectories that keep the indoor air-temperature within certain limits. Similarly, the operator of a gas power plant might favor constant power trajectories over ones that vary quickly because of ramping costs. To express particular preferences and additional operating costs, every system defines a convex flexibility cost function ${f^{(j)}:\mathcal{P}^{(j)}\rightarrow\mathbb{R}}$. If we assume that the aggregator reimburses the individual systems for their flexibility costs, the total costs $T^{(j)}(p^{(j)})$ of having system $j$ implementing the power trajectory $p^{(j)}$ are given by the difference of the flexibility costs and the energy costs:

[TABLE]

If $T^{(j)}$ is negative, money flows from the system owner to the aggregator. If, in contrast, the costs are positive, the system owner receives money from the aggregator.

The disaggregation problem the aggregator needs to solve is how to best distribute a given trajectory ${p^{(\text{agg})}\in\mathcal{P}^{(\text{agg})}}$ into exactly $J$ trajectories ${p^{(1)},\dots,p^{(J)}}$ such that they are feasible and sum up to $p^{(\text{agg})}$. This disaggregation is not unique in general, and the aggregator may want to find the one that is optimal with regard to a particular objective. A profit-oriented aggregator will minimize the total costs by solving {mini} {p^(j)}_j∈J∑_j∈J T^(j)(p^(j))

\addConstraint∑_j∈J p^(j)= p^(agg) \addConstraintp^(j)∈P^(j), j∈J.

V-B Disaggregation using zonotopes

We consider feasible sets $\mathcal{P}^{(j)}$ given by the zonotopes ${\{p^{(j)}\in\mathbb{R}^{N}:p^{(j)}=c^{(j)}+G\beta,\,-\bar{\beta}^{(j)}\leq\beta\leq\bar{\beta}^{(j)}\}}$, ${j\in\mathcal{J}}$. The particular structure of the Minkowski sum of zonotopes (10) allows an efficient disaggregation algorithm whose complexity scales particularly well with the number of systems in the population. To exploit this advantage, we restrict attention to flexibility costs that can be expressed as a function of $\beta$. In particular, we consider cost functions ${f^{(j)}(\beta^{(j)}):\mathbb{R}^{g}\rightarrow\mathbb{R}}$ that can be written as ${f^{(j)}(\beta^{(j)})=f^{(j)}_{1}(\beta^{(j)}_{1})+\dots+f^{(j)}_{g}(\beta^{(j)}_{g})}$, where each ${f^{(j)}_{i}:\mathbb{R}\rightarrow\mathbb{R}}$ is piece-wise linear convex, and $\beta_{i}^{(j)}$ refers to the $i^{\text{th}}$-element of the vector ${\beta^{(j)}\in\mathbb{R}^{g}}$. Substituting $f^{(j)}(p^{(j)})$ with $f^{(j)}(\beta^{(j)})$ and $p^{(j)}$ with ${c^{(j)}+G\beta^{(j)}}$ allows to rewrite the total costs (11) in terms of $\beta$ as

[TABLE]

where $g^{(i)}$ is the $i^{\text{th}}$ generator, and ${T^{(j)}_{\text{fix}}:=-t_{s}v^{(j)^{\top}}c^{(j)}}$ is a fixed cost independent of $\beta^{(j)}$. The components $T_{i}^{(j)}(\beta_{i}^{(j)})$ are piece-wise linear convex. For the purpose of illustration, Fig. 2 shows two example cost functions $T_{i}^{(1)}$ and $T_{i}^{(2)}$. The breakpoint interval lengths and the slopes of the line segments are denoted by $l^{(j)}_{i,t}$ and $q^{(j)}_{i,t}$, respectively, for ${t=1,\dots,m^{(j)}_{i}}$ with $m^{(j)}_{i}$ being the number of line segments.

We define the aggregate cost $T^{(\text{agg})}(\beta^{(\text{agg})})$ as the total cost of the cheapest disaggregation of $\beta^{(\text{agg})}$ into ${\beta^{(1)},\dots,\beta^{(J)}}$: {mini}[3] {β^(j)}_j∈J ∑_j∈JT^(j)(β^(j))

T^(agg)(β^(agg)):= \addConstraint∑_j∈Jβ^(j)= β^(agg) \addConstraint-¯β^(j)≤β^(j)≤¯β^(j), j∈J.

The structure of the individual costs (12) implies that $T^{(\text{agg})}$ can be written as

[TABLE]

with ${T_{\text{fix}}^{(\text{agg})}:=T_{\text{fix}}^{(1)}+\dots+T_{\text{fix}}^{(J)}}$. Each component $T^{(\text{agg})}_{i}$ is piece-wise linear convex and can be constructed by concatenating the line segments of ${T_{i}^{(1)},\dots,T_{i}^{(J)}}$ sorted by ascending slopes $q_{i,t}^{(j)}$, see Fig. 2. The offset of $T_{i}^{(\text{agg})}$ is given by

[TABLE]

The variable $s_{i,t}$ stores the system index that belongs to line segment $t$ of $T_{i}^{(\text{agg})}$: ${s_{i,1}=2}$ for the example in Fig. 2. Up to its offset, each cost component $T_{i}^{(\text{agg})}$, ${i=1,\dots,g}$, is fully defined by the ordered list of triples ${L_{i}:=\{(l_{i,1}^{(\text{agg})},q_{i,1}^{(\text{agg})},s_{i,1}),\dots,(l_{i,M_{i}}^{(\text{agg})},q_{i,M_{i}}^{(\text{agg})},s_{i,M_{i}})\}}$, with ${M_{i}:=m_{i}^{(1)}+\dots+m_{i}^{(J)}}$. These lists are key prerequisites for the disaggregation algorithm discussed below.

V-C Subgradient-based disaggregation algorithm

Given an aggregate trajectory ${p^{(\text{agg})}\in Z^{(\text{agg})}}$, the feasible sets $Z(G,c^{(j)},\bar{\beta}^{(j)})$, ${j\in\mathcal{J}}$, and the aggregate cost function $T^{(\text{agg})}(\beta^{(\text{agg})})$, the disaggregation problem is to solve {mini} β T^(agg)(β)

\addConstraintp^(agg)= c^(agg)+ G β \addConstraint-¯β^(agg)≤β≤¯β^(agg).

The key difference between the above problem and its original formulation (V-A) is that in (V-C) the number of decision variables ($g$), the number of equality constraints ($N$), and the number of inequality constraints ($2g$) are independent of the population size $J$.

Problem (V-C) is convex with piece-wise linear objective and can be written as an LP using an epigraph reformulation, with the objective substituted by a set of linear inequalities. This is inconvenient for the application discussed here because the number of inequalities can be large. Instead, problem (V-C) is solved via a projected subgradient method [24]. The disaggregation algorithm proposed here comprises four steps:

Initialization: Find $\beta(0)$ in the feasible set of (V-C). 2. 2.

Subgradient steps: In iteration $k$, compute a subgradient of $T^{(\text{agg})}$ at $\beta(k)$ and execute a subgradient step. 3. 3.

Projection: Project the new ${\beta(k+1)}$ back onto the feasible set of (V-C). Repeat 2) and 3) until converged, and denote the best $\beta$ found so far by $\beta^{(\text{agg})\ast}$. 4. 4.

Assignment: Distribute $\beta^{(\text{agg})\ast}$ into $\beta^{(1)\ast},\dots,\beta^{(J)\ast}$.

The next sections elaborate on the above steps.

1) Initialization: The subgradient method requires a feasible starting point ${\beta(0)\in\mathbb{R}^{g}}$. It can be obtained by projecting an arbitrary ${\beta\in\mathbb{R}^{g}}$ onto the feasible set of (V-C).

2) Subgradient step: The subgradient method iteratively executes a subgradient step that updates the decision variable as ${\beta(k+1)=\beta(k)-\alpha(k)\nabla(k)}$, where the subgradient and the step size of iteration $k$ are denoted by $\nabla(k)$ and $\alpha(k)$, respectively. A subgradient of $T^{(\text{agg})}$ at $\beta(k)$ is computed element-wise according to line 6 in Fig. 3. Several convergence results for the projected sugradient method exist in the literature [24], including different step size rules that influence the convergence properties. Here, ${\alpha(k)=a/k}$, ${a>0}$, is used, for which the subgradient method converges to the optimal value, i.e. , ${\lim_{k\rightarrow\infty}T(\beta(k))=T^{\ast}}$ [24]. Setting ${a=(\bar{\beta}^{(\text{agg})}_{1}+\dots+\bar{\beta}^{(\text{agg})}_{g})/g}$ proved to be effective.

3) Projection: The $\beta$ updated by the subgradient step in general does not satisfy all constraints of (V-C). To make it feasible, it is projected onto to the feasible set using the Euclidean projection $Q(\cdot)$, see line 9 in Fig. 3.

The subgradient method is not a descent method. Thus, at every iteration, the best objective value found so far is stored. The subgradient algorithm is stopped once the termination criterion in line 15 in Fig. 3 is met. The optimal $\beta^{(\text{agg})\ast}$ and the corresponding objective value $V(k^{\ast})$ are returned.

4) Distribution to individual values: Once the subgradient method has terminated, the best choice of the decision variable $\beta^{(\text{agg})\ast}$ has to be distributed to individual ${\beta^{(1)\ast},\dots,\beta^{(J)\ast}}$. This is done by the algorithm presented in Fig. 4. Finally, the corresponding individual minimum-cost power trajectories are given by ${p^{(j)\ast}=c^{(j)}+G\beta^{(j)\ast},\ j\in\mathcal{J}}$.

V-D Performance of subgradient-based disaggregation

To assess the performance of the subgradient-based disaggregation, its solving time is compared with the solving time of the original LP formulation (V-A) for different populations of PEVs sampled uniformly from the parameter ranges provided in Tab. I. The flexibility costs $f^{(j)}$ are assumed linear with random costs. The optimization models and the subgradient algorithm are implemented in Matlab®. Optimization problems are solved by IBM ILOG CPLEX®111IBM, ILOG, and CPLEX are trademarks of International Business Machines Corp., registered in many jurisdictions worldwide. Other product and service names might be trademarks of IBM or other companies. v.12.6 run on a desktop PC featuring an Intel®  Core i5-2400 CPU @ 3.1 GHz and 16 GB of RAM.

The original disaggregation problem (V-A) is convex and can, in theory, be solved efficiently. However, the problem grows with the population size $J$ and the number of time steps $N$. In the case of PE-systems, problem (V-A) comprises $NJ$ decision variables, $N$ equalities, and ${J(4N-2)}$ inequalities. Therefore, even in the simplest case of linear costs, the solving time and memory requirements of (V-A) can be significant. Table II summarizes the solving time for different problem sizes. For some combinations of $N$ and $J$, the problem size exceeds the maximum memory that Matlab can allocate (abbreviated by M). We note that memory is exhausted for problem sizes of practical interest; for example, participation in the weekly secondary reserve market with 15 min intervals would require the solution of a problem with ${N=672}$.

The solving time of the subgradient-based disaggregation depends on the termination parameters $h$ and $\epsilon$. The values ${h=5}$ and ${\epsilon=10^{-3}}$ have proved to yield good results in practice: the projected subgradient method terminates with a relative optimality gap of less than ${2.5\cdot 10^{-3}}$ for all experiments executed. The maximum solving time for different choices of $J$ and $N$ are provided in Tab. III. The subgradient-based disaggregation clearly outperforms the LP-based disaggregation (V-A). Because the number of decision variables and constraints are independent of $J$, the problem can be solved for large populations and long planning horizons.

VI Use-cases for zonotopic flexibility

VI-A Feasibility

Checking the feasibility of an aggregate trajectory $p^{(\text{agg})}$ for polytopic feasible sets requires solving the disaggregation problem (V-A), which can be computationally expensive. In contrast, checking the feasibility in the case of zonotopes is easy because the aggregate feasible set can be computed explicitly via (10). It holds that ${p^{(\text{agg})}\in Z^{(\text{agg})}}$ if and only if ${|F^{\top}(p^{(\text{agg})}-c^{(\text{agg})})|\leq|F^{\top}G|\bar{\beta}^{(\text{agg})}}$, where ${F\in\mathbb{R}^{N\times q}}$ is the matrix with the normal vectors of all possible zonotope facets as columns, cf. Appendix (A) and [25]. That is, checking whether a given trajectory is aggregate feasible amounts to validating a set of $q$ inequalities.

VI-B The cost of offering regulation power

Aggregators can use flexibility to offer ancillary services. Here, a symmetric constant regulation power service is considered. Offering $\bar{r}$ units of regulation power requires that in every time step of the planning horizon, a deviation of $\bar{r}$ from the baseline power $p^{(\text{agg})}$ is possible. In the power space, this translates into reserving an axis-aligned cube with edge length $2\bar{r}$ and center $p^{(\text{agg})}$ in the aggregate feasible set. A lower bound $\bar{r}_{\max}$ on the maximum regulation power available from the population of systems is {maxi!} ¯r,{p^(j), r^(j)}_j∈J¯r

\addConstraintB(p^(j),r^(j))⊆P^(j), j∈J \addConstraintr^(j)≥0, j∈J \addConstraint∑_j∈Jr^(j)≥¯r1

with ${B(p^{(j)},r^{(j)}):=\{p\in\mathbb{R}^{N}:p^{(j)}-r^{(j)}\leq p\leq p^{(j)}+r^{(j)}\}}$. For a fixed amount of regulation power ${\bar{r}\in[0,\bar{r}_{\max}]}$ and an expected energy wholesale price $\hat{v}$, the aggregator can find the cheapest baseline by solving

[TABLE]

subject to constraints (VI-B)-(VI-B). The expected cost of offering $\bar{r}$ units of regulation power is ${\hat{R}(\bar{r})-\hat{R}(0)}$. This quantity can serve as a lower bound on the regulation power bid price.

However, problem (14) with polytopic feasible sets has the same complexity as the disaggregation problem (V-A) and thus is inconvenient for large populations and long planning horizons. In contrast, for zonotopic feasible sets, an approximative solution to (14) can be found efficiently by using the subgradient-based disaggregation algorithm of Sec. V. In a first step, $\bar{r}_{\max}$ is computed by solving (VI-B) reformulated for zonotopes with generator matrices $G^{PE}$ as

[TABLE]

Denote by ${\beta^{(j)}_{\max},\,j\in\mathcal{J}},$ the optimizers of (15). The constraints guarantee that every feasible set $Z(G^{PE},c^{(j)},\bar{\beta}^{(j)})$ includes a box with semi-edge lengths given by $G\beta^{(j)}_{\max}$. To keep a certain amount ${\bar{r}=\eta\bar{r}_{\max}}$, ${\eta\in[0,1]}$, of regulation power available, every system reserves a portion ${\eta|\beta^{(j)}_{\max}|}$. The remaining ${\beta^{(j)}_{\text{rem}}:=\bar{\beta}^{(j)}-\beta^{(j)}_{\max}}$ defines a zonotope within which the baseline $p^{(j)}$ can be chosen freely. The cheapest aggregate baseline is ${p^{\ast}=c^{(\text{agg})}+G^{PE}\beta^{(\text{agg})\ast}}$ with ${\beta^{(\text{agg})\ast}}$ minimizing

[TABLE]

This problem is solved efficiently by the disaggregation algorithm of Sec. V.

Table IV provides the expected baseline costs $\hat{R}$ for different amounts of regulation power $\bar{r}$ reserved. The costs are computed for a population of 100 PEVs with the parameters sampled from Tab. I. Historic EPEX day-ahead prices are used for $\hat{v}$. Considered are the cases where i) problem (14) is solved using the original polytopic feasible sets, problem (16) is solved for zonotopic feasible sets, and iii) problem (16) is solved for feasible sets approximated by axis-aligned boxes. While the zonotopic approximation is suboptimal compared with the solution found on the original polytopic sets, it still yields significantly lower baseline costs than feasible sets approximated by boxes do. These results indicate that zonotopes are a reasonable compromise between the complex, but accurate polytopic feasible sets and the simple, but inaccurate approximations by boxes.

VII Conclusion

The description, aggregation, and disaggregation of the energetic flexibility of energy resources are important and challenging tasks. This paper introduced a generic and scalable approach for flexible energy systems to quantitatively describe and price flexibility based on zonotopic sets. The description can be used to aggregate flexibility efficiently, disaggregate control decisions, and solve aggregate-level problems such as computing the minimum baseline costs under regulation capacity constraints. It has been shown that zonotopic flexibility yields a reasonable compromise between accuracy of approximation and computational complexity. Future work will include the incorporation of additional features, such as ramp rates, in the zonotope representation. Further, the applicability of distributed optimization methods for solving various disaggregation problems will be investigated.

Appendix A Facets of PE-Zonotopes

In this section we characterize the facets of zonotopes defined by (6)-(7). It follows that zonotopes can reconstruct every possible facet of a PE-polytope.

Theorem 1**.**

Every facet of a full dimensional zonotope ${Z\subseteq\mathbb{R}^{N}}$ defined by (6)-(7) has normal vector ${\sum\{e_{i}:j\leq i\leq k\}}$ for some ${1\leq j\leq k\leq N}$ ($e_{i}$ is the unit vector of coordinate $i$).

Proof. Let $F$ be a facet of $Z$, $H$ the affine hyperplane generated by $F$, and $c$ a normal vector of $F$ so that ${F=\{z\in Z:zc=M\}}$ with ${M=\max\{zc:z\in Z\}}$. Let ${1\leq l\leq N}$ be a coordinate which maximizes $|c_{l}|$. We know that ${c_{l}\neq 0}$. Let $j$ be the minimum, and $k$ the maximum value in ${\{1,\ldots,N\}}$ for which ${c_{i}=c_{l}}$ for all ${j\leq i\leq k}$. $Z$ is centrally symmetric, so without loss of generality we can assume that $c_{l}>0$.

The following claim follows from the definition (6)-(7) of $Z$.

Claim 1**.**

Let ${z=c+\sum_{i=1}^{g}\beta_{i}g^{(i)}\in Z}$. If ${z\in F}$ then

${\beta_{N+j-1}=\overline{\beta}_{N+j}}$* if ${j>1}$,* 2. 2.

${\beta_{N+k}=-\overline{\beta}_{N+k}}$* if ${k<N}$.*

For a set ${X\subseteq\mathbb{R}^{N}}$ and ${I\subseteq\{1,\ldots,N\}}$ we denote by $X|_{I}$ the projection of $X$ onto $\mathbb{R}^{I}$. Let ${I_{1}=\{j,\ldots,k\}}$ and ${I_{2}=\{1,\ldots,N\}-I_{1}}$. ${I_{2}=\emptyset}$ is possible. For ${x,y\in F}$ we define ${\gamma(x,y)=x|_{I_{1}}\times y|_{I_{2}}}$, the combination of $x$ and $y$ along the partition ${I=I_{1}\dot{\cup}I_{2}}$. Let ${x=c+\sum\beta^{x}_{i}g^{i}}$, and ${y=c+\sum\beta^{y}_{i}g^{i}}$. By Claim 1 one can write $\gamma(x,y)$ as ${c+\sum\beta_{i}g^{(i)}\in F}$ where

${\beta_{i}=\beta^{x}_{i}}$ if ${i\in I_{1}}$ and ${\beta_{i}=\beta^{y}_{i}}$ if ${i\in I_{2}}$, 2. 2.

${\beta_{N+i}=\beta^{x}_{N+i}}$ if ${i\geq j}$ and ${i<k}$, 3. 3.

${\beta_{N+i}=\beta^{y}_{N+i}}$ if ${i<j-1}$ or ${i>k}$, 4. 4.

${\beta_{N+j-1}=\overline{\beta}_{N+j}}$ if ${j>1}$, 5. 5.

${\beta_{N+k}=-\overline{\beta}_{N+k}}$ if ${k<N}$.

It follows that ${\gamma(x,y)\in Z}$. By symmetry, ${\gamma(y,x)\in Z}$. Clearly, ${\gamma(x,y)+\gamma(y,x)=x+y}$ and thus $\gamma(x,y),\,\gamma(y,x)\in F$.

Let ${v_{m}\in F|_{I_{m}}\subseteq\mathbb{R}^{I_{m}}}$ for ${m=1,2}$. By definition, there exists a vector ${z_{m}=c+\sum\beta_{i}^{m}g^{(i)}\in F}$ so that ${z_{m}|_{I_{m}}=v_{m}}$ for ${m=1,2}$. It is easy to see that ${v_{1}\times v_{2}=\gamma(z_{1},z_{2})\in F}$. It follows that ${F=F|_{I_{1}}\times F|_{I_{2}}}$, and also that $F|_{I_{m}}$ is the set of points maximizing $c|_{I_{m}}$ in $Z|_{I_{m}}$.

${\dim(F|_{I_{1}})\leq|I_{1}|-1}$ because the vectors $g_{N+i}$ for ${j\leq i\leq k}$ are generators, and clearly ${\dim(F|_{I_{2}})\leq|I_{2}|}$. Now ${N-1=\dim(F)=\dim(F|_{I_{1}})+\dim(F|_{I_{2}})\leq|I_{1}|+|I_{2}|-1=N-1}$, and thus equality most hold throughout. $F|_{I_{2}}$, being the set of points maximizing $c|_{I_{2}}$ in $Z_{I_{2}}$, can only be full-dimensional if ${c|_{I_{2}}=0}$. After normalizing, we get that $c$ has the form ${\sum\{e_{i}:j\leq i\leq k\}}$.∎

Theorem 2**.**

Let ${Z\subseteq\mathbb{R}^{N}}$ be a zonotope defined by (6)-(7). If ${\overline{\beta}>0}$ then for all ${j\leq i\leq k}$ the objective vector ${c=\sum\{e_{i}:j\leq i\leq k\}}$ defines a facet.

Proof. Let ${F=\{z\in Z:zc=M\}}$ with ${M=\max\{zc:z\in Z\}}$, and $H$ the affine hyperplane generated by $F$. Clearly, the following vectors are generators of $H$:

$g^{(i)}$ for ${i<j}$ and ${i>k}$, 2. 2.

$g^{(N+i)}$ for ${i\geq j}$ and ${i<k}$.

The rank of this set of vectors is ${N-1}$, so $H$ is a hyperplane and $F$ is a facet.

It follows from Theorems 1 and 2 that the maximum number of facets of a zonotope is ${N^{2}+N}$ and this bound is reached if ${\overline{\beta}>0}$.

Acknowledgment

The authors gratefully acknowledge the fruitful discussions with Stefan Wörner, Ulrich Schimpel, and Marco Laumanns.