1

Real-Time Minimization of Average Error in the
Presence of Uncertainty and Convexification of
Feasible Sets
arXiv:1612.07287v1 [math.OC] 21 Dec 2016

Andrey Bernstein, Member, IEEE, Niek J. Bouman, Member, IEEE, Jean-Yves Le Boudec, Fellow, IEEE

Abstract—We consider a two-level discrete-time control framework with real-time constraints where a central controller issues
setpoints to be implemented by local controllers. The local
controllers implement the setpoints with some approximation
and advertize a prediction of their constraints to the central
controller. The local controllers might not be able to implement
the setpoint exactly, due to prediction errors or because the
central controller convexifies the problem for tractability. In this
paper, we propose to compensate for these mismatches at the level
of the local controller by using a variant of the error diffusion
algorithm. We give conditions under which the minimal (convex)
invariant set for the accumulated-error dynamics is bounded,
and give a computational method to construct this set. This can
be used to compute a bound on the accumulated error and hence
establish convergence of the average error to zero. We illustrate
the approach in the context of real-time control of electrical grids.
Index Terms—Hierarchical control, real-time control, error
diffusion, convex dynamics, robust set-invariance, power grids

I. I NTRODUCTION
We consider a two-level discrete-time control framework
with real-time constraints, consisting of several local controllers and a central controller. Such a framework was recently proposed in the context of real-time control of electrical
grids [1]. The task of a local controller is (i) to implement
setpoints issued by the central controller, and (ii) to advertise
a prediction of the constraints on the feasible setpoints to the
central controller. In turn, the central controller uses these
advertisements to compute next feasible setpoints for the local
controllers. These setpoints correspond to the solution to an
optimization problem posed by the central controller. Due
to real-time constraints, the central controller is restricted
to operate on convex feasible sets and continuous variables.
Hence, the advertisement of the local controllers is in the form
of convex sets.
Formally, the interaction between the local controllers and
the central controller is assumed to be as in Algorithm 1.
Example 1. Consider a central grid controller, whose task is
to control the grid and the resources connected to it in realtime. Consider a single-phase (or balanced) system, where the
setpoints are pairs (P, Q) ∈ R2 that represent the requested
active (P ) and reactive (Q) power consumption/production.
As an example of the resources, consider a photovoltaic (PV)
The authors are with EPFL, Lausanne, Switzerland.
Manuscript received XXXXXXXX, 2015; revised XXXXXXX, 2015.

Algorithm 1 Interaction between local controllers and central
controller
1: Set n = 0.
2: loop
3:
At time step n, every local controller:
(a) Receives a setpoint request xn sent by the central
controller.
(b) Implements an approximation yn of xn . The implemented setpoint yn is constrained to lie in some set
Sn ∈ S, where S is the collection of all possible
feasible sets of the local controller.
(c) Performs a prediction of its feasible set Sn+1 ∈ S
that will be valid at step n+1, and advertises to the
central controller the convex hull An+1 = chSn+1 .
4:
Upon receiving the advertisements from all the local
controllers, the central controller chooses the setpoints
xn+1 ∈ An+1 for every local controller and sends them
over a communication network.
5:
n := n + 1.
6: end loop

plant and a heater system with a finite number of heating states.
In the case of the PV, the local controller predicts the set of
feasible (implementable) power setpoints Sn+1 in step 3(c) of
Algorithm 1, but the actual set at the time of the implementation may differ from the prediction due to high volatility
of solar radiation. In the case of the heater system, the set
Sn+1 is a finite set. Thus, the corresponding local controller
sends a convex hull of Sn+1 in step 3(c) of Algorithm 1, and
consequently may receive a non-implementable setpoint from
the central controller in the next step.
This example illustrates the source of the potential difference between yn and xn in step 3(b) of Algorithm 1. More
generally, the local controller might not be able to exactly
implement the requested setpoint, because of two reasons:
1) Due to convexification of the feasible set. The actual
set of implementable setpoints might be non-convex. For
example, suppose that the local controller is controlling a
collection of “on-off” devices, which would correspond
to a discrete set of implementable setpoints.
2) Because of uncertainty. The actual set of implementable
setpoints might differ from its prediction, for example,
because of external disturbances.
In this paper, we propose to use the metric of total accu-

2

mulated error between requested and implemented setpoints.
Our goal is to analyze the greedy algorithm for the choice
of yn in step 3(b) of Algorithm 1, namely the algorithm
that performs online minimization of the accumulated error.
(This algorithm is also called the error diffusion algorithm
and is well-known in the context of image processing and
digital printing.) We show conditions on the collection S under
which the accumulated error is bounded for all n and propose
computational methods to find tight bounds. As a consequence,
the average error converges to zero. In Section I-A, we discuss
how our approach compares to the existing literature.
In general, the boundedness of accumulated error (and
thus convergence of the average error to zero) is a relevant
metric in any application where the integral of the control
variable is an important quantity. In this paper, we illustrate
our approach in the context of real-time control of electrical grids [1], where the setpoints are active and reactive
power injections/absorptions, and the integral thereof is the
consumed/produced energy. In this application, “real-time”
means a period time in the order of 100 ms. Thus, the
proposed framework with a fast and efficient central controller
is a natural solution, as discussed in detail in [1]. Other
applications in which boundedness of accumulated error is
relevant include signal processing, digital printing, scheduling,
and assignment problems [2], [3].
A. Related Work
The problem of the local controller can be viewed as
controlling the quantity en defined recursively by
en+1 = en + xn − yn ,
where xn is the requested setpoint by the central controller
and yn is an implemented setpoint by the local controller;
see Algorithm 1. Note that in this context, yn is an input (or
control) variable, while xn can be viewed as an exogenous
disturbance variable.
The task of the local controller is to find a policy that
achieves en ∈ Q for all n, for some bounded set Q. At first
glance, this is a standard goal in control theory, and involves
a robust control invariant set [4]. Moreover, when the set Sn
is non-convex or discrete, this problem is closely related to
control of Mixed Logic Dynamical (MLD) Systems (e.g., [5],
[6], [7]). The classical approach to design a robust controller
in this setting is based on formulating a corresponding MPC
and solving MILP or MIQP problems. However, in real-time
applications, solving MIQP/MILP at every time step may be
not feasible. Moreover, our problem has the following features
which are not present in the classical setting:
1) The feasible set Sn is not known in advance and depends
on the history of the process up to time step n.
2) The feasible set Sn can be uncertain in the sense that it
is not known to the central controller at the time of the
decision making, and thus prediction errors might arise.
3) We allow to control the set of possible “disturbances” An
– it is the advertisement sent by the local controller.
Thus, the methods from the classical control theory do not
apply directly.

A different approach, which we pursue in the present paper,
is to directly analyze a specific algorithm for the choice of
yn , namely a greedy algorithm which, at each time step, minimizes the next-step accumulated error. In particular, it chooses
yn = VorSn (xn + en ), where VorS (z) is the closest point
to z in S. The classical version of this algorithm is known
as error diffusion (or Floyd-Steinberg dithering) in the field
of image processing and digital printing, where the variables
are one-dimensional [8], [9], [10], [11]. The extension to the
general d-dimensional case was considered over the recent
years in several papers. In [12], the problem of a single (fixed)
feasible set S is considered in the special case where S are
the corner points of a polytope A = chS, and an algorithm
to construct a minimal invariant set for the corresponding
dynamical system is proposed. However, the boundedness
of this set is not guaranteed in general. [2] show how to
construct bounded invariant sets for that problem, and extend
the results to a finite collection of polytopes A1 , . . . , AK .
Namely, they show that there exists a bounded set Q that is
simultaneously invariant for the dynamical systems defined
with respect to A1 , . . . , AK . [13] extends the results of [2] to
the case where the polytope may change from step to step,
and argues that there exists a bounded invariant set Q for
this changing dynamical system. Moreover, the conditions are
extended to an infinite collection of polytopes, provided that
the set of face-normals of this collection is finite. Other papers
in this line of research consider specific applications [14], [15],
[3] and/or other special cases [16]. Finally, the optimality of
the error diffusion algorithm was recently analyzed in [3].
B. Our Contribution
Our paper extends the state-of-the-art on the general errordiffusion algorithm, with the following main contributions that
are relevant to our control application:
• We consider general non-convex feasible sets S rather
than corner points of a polytope.
• We consider an uncertain case, in which the feasible set
is not known at the time of advertisement and hence
predicted. Specifically, in this case, it may happen that
A=
6 chS.
• We propose a computational method for constructing the
minimal invariant set for the accumulated error dynamics
in the case of finite collection of feasible sets. We also
show some important special cases in which the minimal
invariant set can be computed explicitly. As a result, we
obtain tight bounds on the accumulated error.
II. N OTATION
Throughout the paper, k · k denotes the `2 norm. We use
N and N>0 to refer to the natural numbers including and
excluding zero, respectively. For arbitrary n ∈ N>0 , we write
[n] for the set {1, . . . , n}.
Let d ∈ N>0 . For arbitrary sets U, V ∈ Rd , U +V represents
the Minkowski sum of U and V, which is defined as U +
V := {u + v | u ∈ U, v ∈ V}. Likewise, U − V represents
the Minkowski difference, defined as U − V := {u − v) | u ∈
U, v ∈ V}. We let chU denote the convex hull of the set U,

3

and by ∂chU we denote the boundary of the convex hull of
U.
For any compact set V, we define the diameter of V as

A. Perfect Prediction

diamV := max{kv − wk : v, w ∈ V}.

Recall from Algorithm 1 that since An = chSn , the setpoint
request xn lies in chSn , while the local controller implements
a setpoint yn according to (2). Therefore, the dynamics for the
accumulated error variable en (1) is given by

Let S ⊂ Rd be an arbitrary non-empty closed set. Any
mapping VorS (x) that satisfies

en+1 = en + xn − VorSn (en + xn ).

VorS (x) = c,

where c ∈ arg min kρ − xk.
ρ∈S

is called a closest-point (or projection) operator onto S. The
Voronoi cell associated with the set S and a point c ∈ S is
defined as
VS (c) = {x ∈ Rd : kx − ck ≤ kx − c0 k, ∀c0 ∈ S}.
For any set V ⊆ Rd , we denote the intersection of V with
VS (c) by
V ∩ VS (c) := V{c}
whenever the set S is clear from the context.
Finally, throughout the paper, a “set” (or “subset”) denotes
a “closed set” (or “closed subset”) unless specified otherwise
explicitly.
III. P ROBLEM D EFINITION
Fix the dimension d ∈ N. Let S be a collection of subsets of
Rd . This collection represents the set of all possible feasible
sets of the local controller. Recall that the interaction between
the local controllers and the central controller is given in
Algorithm 1.
In this paper, we focus on two cases for the prediction step
3(c):
(i) Perfect prediction, namely An+1 = chSn+1 , and
(ii) Persistent prediction, namely An+1 = chSn .
The performance metric considered is the accumulated error
defined recursively by
en+1 = en + xn − yn .

(1)

We analyze the greedy algorithm for the choice of yn ,
namely the algorithm that chooses yn ∈ Sn so that ken+1 k is
minimized. That is,
yn = arg min ken + xn − yk = VorSn (en + xn )

(2)

y∈Sn

which is the closest point to en + xn in Sn . This algorithm is
also known as error diffusion.
As mentioned in the introduction, our goal in this paper is
to:
1) Find conditions on the collection S under which the
accumulated error en is bounded for all n.
2) Propose computational methods to find tight bounds.

Similarly to [2], [13], for any non-empty set S ⊂ Rd and any
x ∈ chS we define the map
GS,x : Rd → Rd

e 7→ e + x − VorS (e + x).

The dynamics for the accumulated error (3) can be then
expressed as
en+1 = GSn ,xn (en ).
(4)
Definition 1 (Invariance). We say that a set Q ⊆ Rd is Ginvariant with respect to a set S ⊂ Rd if
∀ x ∈ chS,

GS,x (Q) ⊆ Q.

We say that Q is G-invariant with respect to a collection S if
it is G-invariant with respect to every S ∈ S.
Remark 1. Our definition of invariance is a special case of
robust positively invariant sets in robust set-invariance theory
[4]. Indeed, for any dynamical system sn+1 = f (sn , wn ),
a robust positively invariant set is any set Ω that satisfies
f (Ω, w) ⊆ Ω for all w ∈ W, where W is a set of all possible
disturbances. In our case, we can view the pair (Sn , xn ) as a
disturbance wn that lies in the set
W := {(S, x) : S ∈ S, x ∈ chS}.
The following observation makes the connection between
invariant sets and boundedness of the accumulated error.
Observation 1. Let S be a collection of subsets of Rd .
Consider the dynamics (4), where Sn ∈ S for all n. Let Q
be a G-invariant set with respect to the collection S. Then, if
e0 ∈ Q, it holds that en ∈ Q for all n ≥ 1.
In particular, it follows from Observation 1 that if Q is a
bounded invariant set that contains the origin, and e0 ∈ Q, then
the accumulated error is bounded for all n by maxv∈Q kvk.
Our next goal is to find minimal G-invariant sets in order
to obtain tight bounds for the accumulated error. In particular,
we provide conditions for the existence of bounded minimal
G-invariant sets, and a method to compute these sets. To that
end, we first define the following set-operators similarly to
[12].
Definition 2 (Set operators induced by a collection of sets).
Fix d ∈ N. For any set S ⊆ Rd , define:
[
gS (Q) :=
(chS + Q){c} − c.
c∈S

IV. M AIN R ESULTS
In this section, we present our main results for the two
cases in the prediction step 3(c) of Algorithm 1. The proofs
are deferred to Section VII.

(3)

For a collection S, let
gS (Q) :=

[
S∈S

gS (Q).

4

We have the following alternative characterization of invariance in terms of the operator gS .
Proposition 1.
(i) For any Q ⊆ Rd ,
Q ⊆ gS (Q).

(ii) Q ⊆ Rd is G-invariant with respect to S if and only if
gS (Q) ⊆ Q.
The next two theorems show how to find minimal invariant
sets and provide conditions on their boundedness.
Theorem 1 (Minimal Invariant Set). Let Q ⊆ Rd . The iterates
gnS (Q) := gS (gn−1
(Q))
S

for n ∈ N, n ≥ 1,

g0S (Q) := Q,

0

are monotonic, in the sense that gnS (Q) ⊆ gnS (Q) for all n ≤
n0 , and the limit set
[
n
g∞
gnS (Q)
S (Q) := lim gS (Q) =
n→∞

n≥0

is the minimal G-invariant set with respect to S that contains
Q.
Theorem 2 (Bounded Invariant Set). If the collection S is
such that chS := {chS, S ∈ S} is a collection of polytopes
such that:
(i) The sizes of the polytopes are uniformly bounded;
(ii) The set N of outgoing normals to the faces of the
polytopes is finite; and
(iii) The bounded Voronoi cells of S are uniformly bounded;
then the minimal G-invariant set with respect to S that
contains the origin is bounded.
In addition, we have the following variants of Theorems 1
and 2 on the existence of minimal convex invariant sets.
Definition 3 (Convex set operator induced by a collection of
sets). Fix d ∈ N. For any set S ⊆ Rd , define:
!
[
GS (Q) := ch(gS (Q)) = ch
(chS + Q){c} − c .
c∈S

For a collection S, let
GS (Q) :=

[
S∈S

GS (Q).
d

Proposition 2. A convex set Q ⊆ R is invariant with respect
to S if and only if
Q = GS (Q)
Theorem 3 (Minimal Convex Invariant Set). Let Q ⊆ Rd .
The following statements hold:
(i) The iterates
GnS (Q) := GS (Gn−1
(Q))
S
G0S (Q) := chQ,

for n ∈ N, n ≥ 1,
0

are monotonic, in the sense that GnS (Q) ⊆ GnS (Q) for
all n ≤ n0 , and the limit set
[
n
G∞
GnS (Q)
S (Q) := lim GS (Q) =
n→∞

n≥0

is the minimal convex invariant set with respect to S that
contains Q.
(ii) Under the conditions of Theorem 2, G∞
S ({0}) is a
bounded set.
We note that the choice between the convex iteration of
Theorem 3 and the original iteration of Theorem 1 reflects the
trade-off between a) performing convex hull at every iteration
to keep the number of vertices describing the iterate at a
minimum, and thereby lowering the cost of computing the
rest of each iteration (intersection, union, and Minkowskisum operations), and b) not computing convex hull in each
step (hence saving this computational cost), at the expense of
having a (possibly) non-convex iterate and a potential increase
in the number of vertices needed to represent it, which in turn
could lead to an increased computational cost of the iteration
as a whole.
1) Computational Method for Computing an Invariant Set:
The iteration of Theorem 3 can be turned into a computational
method (an algorithm that does not necessarily terminate) by
augmenting the iteration with the stopping rule that corresponds to the invariance property (Proposition 10): ending the
iteration when the vertex-representation of Gn+1
(Q) equals
S
that of GnS (Q).
We have implemented the method for the special case of
point sets in C++ with the help of the CGAL library [17],
the source code is available online [18]. To prevent loss of
precision during the iterations, and to be able to perform
exact equality tests, we use exact rational arithmetic, instead
of floating-point arithmetic. Note that this choice restricts all
vertices to have coordinates in Q.
A problem with exact rational arithmetic is that a vertex
coordinate might approach a mixed number (a sum of an
integer and a proper fraction) whose fractional part has small
numerator and denominator (by “small” we mean just a few
digits), like 1/2, 5 13 , 300 56 , etc., but never reach it in finite
time. To mitigate this problem, we apply, in every iteration,
the following conditional rounding function to each coordinate
of every vertex of GnS (Q): for a given  ∈ R and a finite set
of proper fractions with small numerators and denominators
X ⊂ Q, we define
R:Q→Q
(
bqc + t
q 7→
q

if |t − (q mod 1)| ≤ ,
otherwise,

where t := arg minx∈X |t − (q mod 1)|.
Recall that it immediately follows from the stopping rule
that if the method converges, it means that it has found an
invariant set. Hence, we may in principle perturb the set in
an arbitrary way after each iteration in an attempt to aid
convergence. The “rounding trick” outlined above works well
in practice.
Note, however, that by perturbing the set GnS (Q) (through
R) during the iteration, we cannot guarantee anymore that
the method finds the minimal invariant set. Nonetheless, if the
method converges and coordinate rounding occurs only just
before convergence, in other words, if the last “rounding-free”
iterate is δ-close (measured by a suitable metric for sets) to

5

1.0

1.0

0.5
0.5
-1

-0.5

1

2

3

4

5
-1.0 -0.5
-0.5

-1.0
Fig. 1. A single-point-set example. The minimal invariant error set is shown
in light gray.

the invariant set found by the method, then that invariant set
is a δ-close approximation to the minimal invariant set, by the
monotonicity the iterates.
An additional benefit of applying R is that the method is
likely to find an invariant set whose vertex-coordinates have
small representation.
2) Numerical Examples in R2 : In a first example, we let
S = {S}, with S = X × Y, where X = {−1, 1, 3, 5} and
Y := {−1, 1}, i.e., eight points on a rectangular grid. Figure 1
shows these vertices and the minimal invariant error set (in
gray), which was found by our computational method (based
on Theorem 3) after one iteration.
In our second example, we consider S0 = {S1 , S2 , S3 }, with

0.5 1.0

-1.0
Fig. 2. The set of setpoints S1 (all points). The gray and lightgray points
indicate the setpoints that are removed from S1 to construct S2 and S3 .

1

-3

-2

-1

1

2

3

-1
S''=8S1 <
-2

S''=8S2 <

-3

S''=8S3 <

-4

S1 := {(−1, −1), (0, −1), (1, −1), (1, 0), (1, 1), (0, 1), (−1, 1), (−1, 0)},
S2 := S1 \ {(0, −1)},

S3 := S2 \ {(−1, −1)}.
In words: the set S1 is a collection of points that are placed
equidistantly on a rectangle; see Figure 2, and the set S2 and
S3 respectively are created by removing one resp. two points
from S1 . In particular, it holds that S3 ⊂ S2 ⊂ S1 . This could
correspond to a setting in practice where the local controller
can implement points from S1 most of the time, but once
in a while one particular setpoint, (and sometimes even an
additional particular setpoint) becomes temporarily infeasible.
The minimal invariant error set, shown in Figure 3, has 7
vertices and was found after 177 iterations, with  = 10−8
as rounding parameter. From this figure, we see that the
minimal invariant error set corresponding to S0 (the “joint”
error set) is significantly larger than the minimal invariant
error sets corresponding to singletons S00 = {Si } for all
i ∈ {1, 2, 3}. Hence, if having a small invariant error-set is
of central importance in an application, then one could, in
this particular example, decide to only use the setpoints in S3
at the cost of providing a smaller feasible set (on average) to
the central controller.
B. Persistent Prediction
We now analyze the case of using a persistent predictor
in the local controller (step 3(c) of Algorithm 1), where the
advertised predicted feasible set is given by An+1 = chSn .
We then have the following error dynamics
en+1 = en + xn − VorSn (en + xn )
where xn ∈ An = chSn−1 and Sn , Sn−1 ∈ S. Observe
that this dynamics involves a pair of sets rather then a single
set, and hence the previous results cannot be applied directly.

Fig. 3.
The minimal invariant error set for S0
=
{S1 , S2 , S3 };
a
convex
polygon
with
corner-vertices
{(−3, − 72 ), (−1, − 92 ), (1, − 92 ), (1, 12 ), ( 12 , 1), (− 12 , 1), (−3, 12 }.
In
the same figure, the invariant error sets are drawn corresponding to the cases
00
where S = {Si }, for all i ∈ {1, 2, 3}. Those error sets are much smaller
than the “joint” invariant error set.

Fortunately, we can reformulate this dynamics in terms of the
modified request zn = en + xn similarly to the approach in
[2], [12]. For this new state variable, we have that
zn+1 = zn + xn+1 − VorSn (zn ),

(5)

where xn+1 ∈ An+1 = chSn . Thus, the following operator
can be defined.
Definition 4. For any finite non-empty set S ⊂ Rd and any
x ∈ chS we define the map
FS,x : Rd → Rd

z 7→ z + x − VorS (z)

The dynamics (5) can be then expressed as
zn+1 = FSn ,xn+1 (zn )

(6)

and F -invariance with respect to S is given by Definition 1
by replacing G with F .
The following proposition makes the connection between
F -invariant sets and boundedness of the accumulated error.
Proposition 3. Let S be a collection of subsets of Rd . Consider
the dynamics (6), where Sn ∈ S for all n. Let D be an F invariant set with respect to the collection S, and assume that
z0 ∈ D.
(i) Then it holds that
ken k ≤ max max kvk
S∈S v∈D−chS

6

for all n ≥ 1.
(ii) If in addition S ⊆ D for all S ∈ S, we have that ken k ≤
diamD for all n ≥ 1.
Proof: Part (i) of the proposition follows trivially by the
invariance of D, and the fact that en = zn − xn for zn ∈ D
and xn ∈ chS for some S ∈ S.
For part (ii), observe that
en+1 = xn + en − VorSn (xn + en ) = zn − VorSn (zn )

where both zn ∈ D and VorSn (zn ) ∈ Sn ⊆ D. Hence,
ken+1 k ≤ diamD.
We next state our main result that provides conditions for the
existence of bounded minimal F -invariant sets, and a method
to compute these sets. Hence we obtain tight bounds for the
accumulated error in the case of persistent prediction. To that
end, we first define the following set-operators similarly to
[12].
Definition 5 (Set operators). Fix d ∈ N. For any set S ⊆ Rd ,
define:
[
pS (D) := chS +
D{c} − c, PS (D) := ch(pS (D)).
c∈S

Also, for a collection S, let
[
pS (D) :=
pS (D),

PS (D) := ch(pS (D)).

S∈S

Also, define the iterates of the above
pnS (D), PnS (D), pnS (D) and PnS (D), as before.

operators,

Proposition 4.
(i) For any set D ⊆ Rd , D ⊆ pS (D), and D is F -invariant
with respect to S if and only if
pS (D) ⊆ D.

(ii) For any convex set D ⊆ Rd , D ⊆ PS (D), and convex
D is F -invariant with respect to S if and only if
PS (D) ⊆ D.

Theorem 4. Let S be a collection of sets S ⊂ Rd , and let D
be any given set. The following statements hold:
(i) The iterates pnS (D), PnS (D) are monotonically nondecreasing, and the sets
n
p∞
S (D) := lim pS (D) and
n→∞

n
P∞
S (D) := lim PS (D)
n→∞

are the minimal F -invariant set and the minimal convex
F -invariant set, respectively, containing the set D.
(ii) Under the conditions of Theorem 2, both p∞
S ({0}) and
P∞
S ({0}) are bounded sets.
V. A PPLICATION E XAMPLES AND N UMERICAL
I LLUSTRATION
In this section, we give concrete examples of the design
of local controllers for some important types of resources in
the context of real-time control of electrical grid. We also
perform numerical simulation of the COMMELEC system [1]
that includes one central controller (grid agent) and several
local controllers (resource agents), and show how boundedness of the accumulated error helps to achieve better overall
performance of the closed-loop COMMELEC system.

A. Discrete Resource in One Dimension
In this section, we show how our general results apply to
local controllers that can only implement setpoints from a discrete and one-dimensional set. As an application, we consider
a heating system consisting of a finite number of heaters that
each can either be switched on or off (see Section V-B below).
In particular, this system can only produces/consumes real
power, and thus its set of feasible setpoints is one-dimensional.
We focus here on deterministic systems, hence we are in the
case of perfect prediction.
Consider a finite collection S = {S k }k∈[K] of finite nonempty subsets S k ⊆ R, where K ∈ N. For any S ∈ S, whose
elements we label as s1 < s2 < . . . < s|S| , we let
(
0
if |S| = 1
∆S :=
maxi∈[|S|−1] si+1 − si if |S| > 1.
denote the maximum stepsize of S. We also let
∆S := max ∆S
S∈S

(7)

denote the maximum step size of the collection S.
Theorem 5. Let S be a finite collection of finite non-empty
subsets of R. Then the set [−∆S /2, ∆S /2] is the minimal
G-invariant set with respect to S (as per Definition 1) that
contains the origin. Hence, if e0 ∈ [−∆S /2, ∆S /2], the
accumulated error under the greedy algorithm (2) is bounded
by |en | ≤ ∆S /2 for all n ≥ 1.
Proof: We use the iteration of Theorem 1 to show that
g∞
({0})
= [−∆S /2, ∆S /2]. For any S ∈ S, it is easy to see
S
that
[
gS ({0}) =
chS{c} − c = [−∆S /2, ∆S /2].
c∈S

Indeed, consider the two points in S that attain the maximum
step size, namely si and si+1 such that si+1 − si = ∆S . Then
[

chS{si } − si
chS{si+1 } − si+1 = [−∆S /2, ∆S /2]
and the rest of the terms chS{c} − c are contained in
[−∆S /2, ∆S /2]. Therefore, the first iteration yields
[
gS ({0}) =
[−∆S /2, ∆S /2] = [−∆S /2, ∆S /2] := Q.
S∈S

For the second iteration, for any S ∈ S, consider
[
g2S ({0}) = gS (Q) =
(chS + Q){c} − c.

(8)

c∈S

Denote chS = [s1 , s|S| ]. Observe that for any si ∈ S, i 6=
1, |S|,
(chS + Q){si } = VS (si )

= [si − (si − si−1 )/2, si + (si+1 − si )/2]
⊆ [si − ∆S /2, si + ∆S /2]

by the definition of the Voronoi cell of si and of the maximal
step size of the collection S. For c = s1 ,
(chS + Q){s1 } = [s1 − ∆S /2, s1 + (s2 − s1 )/2]
⊆ [s1 − ∆S /2, s1 + ∆S /2],

7

and similarly
(chS + Q){s|S| } = [s|S| − (s|S| − s|S|−1 )/2, s|S| + ∆S /2]
⊆ [s|S| − ∆S /2, s|S| + ∆S /2].

Hence, for all c ∈ S, (chS + Q){c} − c ⊆ [−∆S /2, ∆S /2].
Substituting in (8) yields
g2S ({0}) = gS (Q) ⊆ Q
and consequently gS (Q) = Q by Proposition 1. Therefore,
[
g2S ({0}) = gS (Q) =
gS (Q) = Q.
S∈S

Thus the iteration has converged and, by Theorem 1,
g∞
S ({0}) = Q := [−∆S /2, ∆S /2] is the minimal invariant
set.
B. Example: Resource Agent for Heating a Building
In this section, we present a concrete local-controller example: we will design a local controller for managing the
temperature in a building with several rooms. The reason for
showing this example is twofold. First, we wish to give a
concrete example of a local controller that controls a load
that can only implement power setpoints from a discrete set.
Second, the local-controller design shows a concrete usage
example of the COMMELEC framework [1], and might serve
as a basis for an actual implementation.
In this section and next section, we will use the terminology
of the COMMELEC framework and denote the local controller
as “resource agent”, while the advertised feasible set of power
setpoints An as “P Q profile”, where P and Q denote active
and reactive power, respectively.
The heating system’s objective is to keep the rooms’ temperatures within a certain range. For rooms whose temperature
lies in that range, there is some freedom in the choice of the
control actions related to those rooms. The resource agent’s job
is to monitor the building and spot such degrees of freedom,
and expose them to the grid agent, which can then exploit
those for performing Demand Response.
Our example is inspired by [19], which also considers the
problem of controlling the temperature in the rooms of a
building using multiple heaters. We address two issues that
were not addressed in [19]:
1) We show that by rounding requested setpoints into implementable setpoints using the error diffusion algorithm we
obtain a resource agent with bounded accumulated-error.
2) We prevent the heaters from switching on and off with
the same frequency as COMMELEC’s control frequency,
which is crucial in an actual implementation.
1) Simple Case: a Single Heater: For simplicity, we first
analyze a scenario with only one heater. The main aspects
of our proposed design (as mentioned above) are in fact
independent of the number of heaters, and we think that those
aspects are more easily understood in this simple case. We
will generalize our example to an arbitrary number of heaters
in Section V-B2.

a) Model and Intended Behavior: We model the heater
as a purely resistive load (it does not consume reactive power)
that can be either active (“on”) or inactive (“off”). It consumes
Pheat > 0 Watts while being active, and zero Watts while being
inactive.
From the perspective of the resource agent, the heater has
a state that consists of two binary variables: sn ∈ {0, 1},
which corresponds to whether the heater is on (sn = 1) or
off (sn = 0), and `n ∈ {0, 1} indicates whether the heater
is “locked”, in which case we cannot switch on or switch off
the heater. Formally, if `n = 1 then sn+1 = sn necessarily
holds. Hence, `n exposes a physical constraint of the heater,
namely that it cannot (or should not) be switched on and off
with arbitrarily high frequency. In the typical case where the
minimum switching period of the heater is (much) larger than
the COMMELEC’s control period (≈ 100 ms), the heater will
“lock” immediately after a switch, i.e., assuming `n = 0,
setting sn+1 such that sn+1 6= sn will induce `n0 = 1 for
every n0 ∈ [n + 1, n + N ], after which `n+N +1 = 0. Here,
N ∈ N represents the minimum number of timesteps for which
the heater cannot change its state from on to off or vice versa.
Suppose that the heater is placed in a room, and that the
temperature of this room is a scalar quantity. (We do not aim
here to model heat convection through the room or anything
like that.) The temperature in the room, denoted as Tn , should
remain within predefined “comfort” bounds,
Tn ∈ [Tmin , Tmax ] n = 1, 2, . . .
where Tmin , Tmax ∈ R. If Tn is outside this interval (and only
if `n = 0), then the resource agent should take the trivial
action, i.e., ensure that the heater is active if Tn < Tmin , and
inactive if Tn > Tmax . The more interesting case is if Tn lies
in [Tmin , Tmax ] (again, provided that `n = 0), as this gives
rise to some flexibility in the heating system: the degree of
freedom here is whether to switch the heater on or off, which
obviously directly corresponds to the total power consumed
by the heating system. The goal is to delegate this choice
to the grid agent, which we can accomplish by defining an
appropriate COMMELEC advertisement.
b) Defining the Advertisement and the Rounding Behavior: Let the discrete set of implementable real-power setpoints
at time n be defined as


{0}
if (`n = 0 ∧ Tn > Tmax )∨





(`n = 1 ∧ sn = 0),

Sn := {−Pheat , 0} if `n = 0 ∧ Tmin ≤ Tn ≤ Tmax ,



{−Pheat }
if (`n = 0 ∧ Tn < Tmin )∨




(`n = 1 ∧ sn = 1),
where ∧ and ∨ stand for “and” and “or”, respectively. Note
that Sn only contains non-positive numbers, by the convention
in COMMELEC that consuming real power corresponds to
negative values for P . We define the P Q profile as a perfect
prediction of chSn , namely An = chSn , and assume that
the resource agent uses the error diffusion algorithm (2) to
implement setpoints yn . We then have the following immediate
corollary of Theorem 5.

8

P

Corollary 1. The accumulated error of the single-heater
resource agent as defined in this section is bounded by 12 Pheat .
2) General Case: an Arbitrary Number of Heaters: Here,
we extend the single-heater case to a setting with r heaters, for
r ∈ N, r ≥ 1 arbitrary. As we will see, also this multi-heater
case can be analyzed using Theorem 5.
Like in the single-heater case, we assume that each heater is
purely resistive. We furthermore assume that heater i consumes
Piheat Watts of power when active (and zero power when
inactive), for every i ∈ [r]. Also similarly to the singleheater case, we assume that each heater is placed in a separate
(i)
room, whose (scalar) temperature is denoted as Tn . Not
surprisingly, our objective shall now be to keep the temperature
in each room within the predefined comfort bounds, i.e.,
Tn(i) ∈ [Tmin , Tmax ] ∀n ∈ N>0 , ∀i ∈ [r].
In the one-heater case, the only degree of freedom is the
choice to switch that heater on or off. In case of multiple
heaters, there is potentially some freedom in choosing which
subset of the heaters to activate, and note that there will
typically1 be an exponential number of those subsets (exponential in the number of heaters). Each subset corresponds
to a certain total power consumption, i.e., a power setpoint.
As in the single-heater case, the P Q profile will be defined
as the convex hull of the collection of these setpoints. When
the grid agent requests some setpoint from the P Q profile,
the resource agent has to select an appropriate subset whose
corresponding setpoint is closest (in the Euclidean sense) to
the requested setpoint. Note that there can be several subsets of
heaters that correspond to the same setpoint. A simple method
to resolve this ambiguity would be, for example, to choose the
subset consisting of the coldest rooms, however, as this topic
is beyond the scope of this work, we leave the choice of such
a selection method to the resource-agent designer.
When going from the single-heater setting to a multipleheaters scenario, we merely need to re-define Sn , which we
will name Sen here to avoid confusion with the single-heater
case. The definition of the P Q profile, and rule for computing
yn given in Section V-B1b also apply to the multi-heater case,
provided that all occurrences of Sn in those definitions are
replaced by Sen .
(i)
(i)
For every i ∈ [r], let sn and `n represent the state
variables sn and `n (as defined in the single-heater case)
(i)
for the i-th heater. Let Ln := {i ∈ [r] : `n = 1}
denote the set of rooms whose heater is locked at timestep
(i)
n. Furthermore, let Cn := {i ∈ [r] : Tn < Tmin } and
(i)
Wn := {i ∈ [r] : Tmin ≤ Tn ≤ Tmax }. Informally speaking,
Cn contains the rooms that are “too cold”, and Wn the rooms
whose temperatures are within the comfort bounds.
If A ⊆ [r], we write A for the complement with respect to
[r], i.e. A := [r] \ A.
Let
X

Sen := an −
Piheat : I ⊆ Ln ∩ Wn
(9)
i∈I

1 Provided

that not too many heaters are locked.

Pmax
x
T (x)

d

te

S ra

ϕ
Q
0

Fig. 4. The parameterized collection of sets {T (x) : x ∈ [0, Pmax ]}, and
its relation to the rated-power constraint of the PV converter.

represent the set of implementable (active) power setpoints,
with
X
X
heat
Pjheat .
s(i)
−
an := −
n Pi
i∈Ln

j∈Ln ∩Cn

As in the one-heater example, we use Theorem 5 to bound
the accumulated error of the resource agent. To this end, let e
S
denote the collection of all possible sets Sen (9). It is easy to see
that this is a finite collection. Further, the maximum stepsize
of this collection (7) is given by ∆eS = maxi∈[r] Piheat . This
gives us the following corollary.
Corollary 2. The accumulated error of the multiple-heaters
resource agent as defined in this section is bounded by
( 12 maxi∈[r] Piheat ).
C. Example: Resource Agent for a Photovoltaic (PV) System
Here, we explain how we can apply Theorem 4 to devise
a local controller (or a resource agent in the terminology of
[1]) for a PV system with bounded accumulated-error.
Let Srated and φmax denote the rated power of the converter
and angle corresponding to the minimum power factor, respectively. We suppose that these quantities are given (they
correspond to physical properties of the PV system), and that
Srated ≥ 0 and 0 ≤ φmax < π. Note that for any power setpoint
(P, Q), the rated power imposes the constraint P 2 + Q2 ≤
2
Srated
; the angle φmax imposes that arctan(Q/P ) ≤ φmax .
Let us now choose Pmax , ϕ ∈ R such that 0 ≤ Pmax ≤
max
Srated , 0 ≤ ϕ ≤ φmax , and Pcos
ϕ = Srated , and let
|Q|
≤ tan ϕ}
P
be a triangle-shaped set in the P Q plane; see Figure 4 for
an illustration. Note that for any combination of Pmax and ϕ,
the triangle T (x) for any x ∈ [0, Pmax ] is fully contained in
the disk that corresponds to the rated-power constraint, and,
moreover, the two upper corner points of T (Pmax ) lie on the
boundary of that disk.
Let pmax
be the maximum real power available at timestep
n
n ∈ N (typically determined by the solar irradiance). Using
pmax
n , we define the set of implementable points at timestep n
as
Sn := T (min(pmax
n , Pmax )).
T (x) := {(P, Q) ∈ R2 : 0 ≤ P ≤ x,

9

Theorem 6. Let Pmax , ϕ, and Sn (for every n ∈ N) be defined
as above for a given PV system. Consider a resource agent for
this PV system that: (i) uses a persistent predictor to advertise
An+1 := Sn , and (ii) implements setpoints according to the
greedy (error diffusion) algorithm. Then, D = T (Pmax ) is the
minimal F -invariant set with respect to S = {T (x)}0≤x≤Pmax .
Therefore, if e0 ∈ D, the accumulated error for this resource
agent is bounded by


Pmax
, 2 Pmax tan ϕ
ken k ≤ diamD = max
cos ϕ
for all n ∈ N.
The proof of Theorem 6 relies on the following general
result.
d
Lemma 1. Let
T S be a collection of non-empty subsets of R .
Assume that S∈S S =
6 ∅. Let
[
chS.
D :=
S∈S

If
T D is F -invariant with respect to S then, for any s0 ∈
S∈S S, D is the minimal F -invariant set with respect to S
that contains s0 .
T
Proof: Let s0 ∈ S∈S S. We use the iteration of Theorem 4 to prove that D is the minimal F -invariant set that
contains s0 . For the first iteration, for every S ∈ S, we have
by Definition 5 that
[
pS ({s0 }) := chS +
{s0 }{c} − c
c∈S

=

chS + s0 − s0 = chS,

where the second equality follows by the fact that s0 ∈ S.
Therefore,
[
[
pS ({s0 }) :=
pS ({s0 }) =
chS := D.
S∈S

S∈S

Now for the second iteration,
p2S ({s0 }) = pS (D) = D
by the invariance of D. Thus the iteration has converged, and
by Theorem 4, D is the minimal F -invariant set with respect
to S that contains s0 .
T
Proof of Theorem 6: Observe that S∈S S = {0} and
[
[
[
chS =
S=
T (x) = T (Pmax ).
S∈S

S∈S

0≤x≤Pmax

Hence, it remains to show that T (Pmax ) is F -invariant with
respect to S. The minimality property then follows immediately from Lemma 1. By Proposition 4, it is enough to show
that
pS (T (Pmax )) ⊆ T (Pmax ).
In other words, we want to show that for any S ∈ S,
[
S+
T (Pmax ){c} − c ⊆ T (Pmax ).
c∈S

First, note that it is straightforward to characterize the different
types of Voronoi cells of S (see also Figure 5): for interior

points, the Voronoi cell is the point itself; for non-corner
points on the boundary, the Voronoi cell is the outwardpointing ray that emanates from that point and is normal to
the facet; for corner points on the boundary, the Voronoi cell
is a cone, namely the union of all rays that emanate from that
corner point, whose directions vary (continuously) between the
normals of the adjacent facts.
Now, it is not hard to
S see that for any x ∈ [0, Pmax ], we
can construct G(x) := c∈T (x) T (Pmax ){c} − c as shown in
Figure 5. From this construction it then immediately follows
that T (x) + G(x) = T (Pmax ), which proves F -invariance of
T (Pmax ) with respect to S. Hence, from Proposition 3 we then
have that ken k ≤ diamD for all n ∈ N. Because T (Pmax )
is an isosceles triangle, its diameter is either Pmax / cos ϕ (the
length of one of its legs) for ϕ ≤ π6 or 2 Pmax tan ϕ otherwise
(the length of its base).
D. Simulation
As in [20], we take a case study that makes reference to the
low voltage microgrid benchmark defined by the CIGRÉ Task
Force C6.04.02 [21]. For the full description of the case study
and the corresponding agents design, the reader is referred to
[20].
There are two modifications compared to the original case
study: (i) The PV agents are updated with the algorithm
described in Section V-C, and (ii) instead of using an uncontrollable load in the case study, we use a resistive heaters
system, and the corresponding agent is implemented according
to the methods described in Section V-B.
We simulate a rather extreme scenario involving a highly
variable solar irradiance profile. That is, we let the irradiance
vary according to a square wave with a period of 300 ms. This
will cause the PV agent’s P Q profile to be highly variable.
We let the cost function of the PV agent be the same as in
[20]; this cost function encourages to maximize active-power
output. The cost function of the heater is set to a quadratic
function, whose minimum lies at half the heater power, namely
at −7.5 kW. With respect to the locking behavior of the heater,
we let it lock for one second after a switch.
The results are shown in Figures 6–11. For comparison,
we run the same scenario with resource agents for which
the accumulated error might grow unboundedly. I.e., those
RAs do not apply the error-diffusion technique described in
this paper, instead, they just project the request to the closest
implementable setpoint, like in [20].
Figure 6 and 8 illustrate how the use of the error-diffusion
algorithm in the heater agent affects the convergence properties
of the closed-loop system. Recall that as in [20], the central
controller (i.e., the grid agent) is using a gradient-descent
algorithm to compute power setpoints. It can be seen from
Figure 6 that without error-diffusion, the grid agent “gets
stuck” on a setpoint that is not close to the optimal value,
and the implemented setpoint is also suboptimal. On the other
hand, with error-diffusion, both requested and implemented
setpoints are optimal on average, as can also be seen in
Figure 8. Similar behavior can be observed in Figure 9 and
11 for the PV agent. Moreover, Figure 9 shows how the errordiffusion algorithm helps to improve utilization of renewables.

10

P
P

Pmax
P

VorS (w)
VorS (v)

x

v

w

v2

S = T (x)
u
VorS (u) = {u}

x

v3

G(x)
Q
(0, 0)

T (x)

ϕ

ϕ
Q

Q
v1 = p1 = (0, 0)

(0, 0)

Fig. 5. Left – The three Voronoi cell-types for the set S, 1) a singleton for any interior point (u), 2) a cone for any corner point (v), and 3) a ray (normal
to the corresponding facet) for any non-corner point on the boundary (w). Right – Construction of the set G(x).

Without Error Diffusion

−10

−15
0

0.5

1

Without Error Diffusion

400

Req.
Impl.

−5

1.5

2
Time [s]

2.5

3

3.5

4

Real Power [kW]

Real Power [kW]

0

300
200
100
0
−100

With Error Diffusion

0

1

2

3

4

0

0.5

1

1.5

2
Time [s]

2.5

3

3.5

4

Fig. 6. Sequence of requested vs. implemented setpoints belonging to the
heater agent. In the top figure (no error diffusion) the grid agent “gets stuck”
on a setpoint that is not close to the optimal value. In the bottom figure (with
error diffusion), the heater switches on and off with an appropriate duty cycle
(the switching frequency is limited by the locking-duration parameter).

In particular, there is less curtailment of the PV, hence more
energy is produced from this renewable energy source. Finally,
the simulation results corroborate the theoretical analysis by
showing delivery of the requested power on average, i.e.,
energy, which can be valuable for an application like a virtual
power plant (Figure 7, 8, 10 and 11).
VI. D ISCUSSION
In this section, we comment and emphasize several important issues related to our results.
We first discuss the following conjecture related to Theorem
1 and 2.
Conjecture: Consider a case where the collection S
contains a single point set S, which is shown in Figure 12.
One vertex, p ∈ S, lies inside the convex hull of the four
other points in S, whose coordinates are (±10, ±10). When
p is moved towards the boundary of chS, it is easy to see
that the corresponding (bounded) Voronoi cell (shown in dark
gray) can become arbitrarily large. (Note that if p is placed
on the boundary, its Voronoi cell becomes an unbounded set).
Interestingly, we observe that the minimal invariant set (light
gray) also becomes larger as p is moved towards ∂chS, in
particular, it exactly covers the Voronoi cell (when the latter is

Real Power [kW]

−15

6

5
Time [s]

6

7

8

9

10

7

8

9

10

5
0
−5
−10
−15
−20

0

1

2

3

4

Fig. 7. Plots of the accumulated error of the heater agent. In the absence of
error diffusion (top figure), the accumulated error grows linearly with time.
With error diffusion (bottom figure), the accumulated error is bounded from
above and below.

Without Error Diffusion
0
Real Power [kW]

−10

5
Time [s]

With Error Diffusion

10

Req.
Impl.

−5

−5

−10
Req.
−15

0

1

2

3

Impl.
4

5
Time [s]

Objective
6

7

8

9

10

8

9

10

With Error Diffusion
0
Real Power [kW]

Real Power [kW]

0

−5

−10
Req.
−15

0

1

2

3

Impl.
4

5
Time [s]

Objective
6

7

Fig. 8. Time-averaged sequences of requested and implemented setpoints
belonging to the heater agent. In the absence of error diffusion (top figure), the
implemented setpoint does not converge to the objective, where the objective
is the minimizer of the cost function. With error diffusion (bottom figure),
the sequence of implemented setpoint converges towards the objective.

11

60

Without Error Diffusion

Real Power [kW]

10
Req.

8

Impl.

50

Max. Power

6

40

4
2

30

0
0

0.5

1

1.5

2
Time [s]

2.5

3

3.5

4

With Error Feedback

5

Real Power [kW]

10
8

-10

6
Req.

Impl.

15

10

Max. Power

-5

p

5
5

10

-5

-10 -5
-5

-10

-10

4
2

20

10

p

10

5

10

-10 -5

p
5 10

-10

0
0

0.5

1

1.5

2
Time [s]

2.5

3

3.5

4

Fig. 9. Sequence of requested vs. implemented setpoints belonging to the
PV agent. The bottom figure shows that error diffusion helps to maximize
utilization of the PV.
Without Error Diffusion

Real Power [kW]

0
−20
−40
−60
−80
−100

0

1

2

3

Real Power [kW]

4

5
Time [s]

6

5
Time [s]

6

7

8

9

10

7

8

9

10

With Error Diffusion

5
0
−5
−10
−15

0

1

2

3

4

Fig. 10. Plot of the accumulated error of the PV agent. Like in the heateragent case, the accumulated error grows unboundedly in the absence of error
diffusion (top figure), whereas the accumulated error is bounded with error
diffusion.
Without Error Diffusion

Real Power [kW]

6

Req.

5

shifted by −p) in the cases corresponding to the three positions
of p shown in Figure 12. This leads us to conjecture that an
invariant region must cover all bounded Voronoi cells, when
shifted by their corresponding vertex.
This conjecture is also corroborated by condition (3) in
Theorem 2. Indeed, in order to proof boundedness of the
minimal invariant set, we require that all bounded Voronoi
cells are uniformly bounded.
Finally, we would like to point out that our result for
the case of the PV system presented in Section V-C can be
generalized to any local controller with collection S having
the following properties:
•

•

•

Impl.

4
3

Every S ∈ S is convex.
The sets are monotonic in the sense that for any S, S 0 ∈
S, either S ⊆ S 0 or S 0 ⊆ S. Thus, there exists Smax ∈ S
such that S ⊆ Smax for all S ∈ S.
For every S ∈ S:
[
S+
(Smax ){c} − c ⊆ Smax ,
c∈S

2

or in other words, for every x ∈ S and y ∈ Smax

1
0

0

1

2

3

4

5
Time [s]

6

7

8

9

10

With Error Diffusion

6
Real Power [kW]

Fig. 12. A case where the bounded Voronoi cell (in dark gray), corresponding
to a point p in the interior of the convex hull of all points (chS), can grow
arbitrarily large as p is moved towards the boundary of chS. The corner points
of S are at coordinates (±10, ±10). The three plots correspond to three
positions of p: px = 0 and py ∈ {0, 5, 9}, where px and py represent the xand y-coordinate of p. In each plot, the minimal invariant set corresponding
to S is shown in light gray, and exactly covers the bounded Voronoi cell that
corresponds to p (when the latter is translated by −p) in these examples; we
conjecture that the minimal invariant set must cover all bounded Voronoi cells
(when translated appropriately).

Req.

5

x + y − VorS (y) ∈ Smax .
Under these conditions, it is easy to show that the set Smax
is the minimal F -invariant set with respect to S. The proof is
identical to that of Theorem 6.

Impl.

4
3
2
1
0

0

1

2

3

4

5
Time [s]

6

7

8

9

10

Fig. 11. Time-averaged sequence of requested vs. implemented setpoints
belonging to the PV agent. Also this plot shows that error diffusion helps to
increase the utilization of the PV.

VII. P ROOFS
The proofs in this section extend the previous results to our
more general case. In particular, the proof of minimality is
an extension of the proofs given in [12], while the proof of
boundedness is an extension of proofs in [2], [13].

12

A. Preliminaries
We next list well-known results that are useful in our proofs.
Proposition 5 (Convex hull is non-decreasing). For all sets
A and B with A ⊆ B,
chX ⊆ chY
Corollary 3. For all sets A and B,
ch(A ∪ B) ⊇ chA ∪ chB
Proposition 6 (The Minkowski sum is distributive over
unions). Let G be a group and let A, B and C be arbitrary
subsets of G. Then,
A + (B ∪ C) = (A + B) ∪ (A + C).
Proposition 7 (Minkowski sums and convex hulls commute).
For arbitrary sets A and B,
ch(A + B) = chA + chB.
B. Proofs for the Case of Perfect Prediction
In this section, we present the proofs for the results of
Section IV-A. For clarity of exposition, we first present proofs
for the case of single set S, and then extend to a collection S.
1) Single-Set Case: Recall Definition 2 of the operator gS .
Lemma 2. For any A ⊆ Rd , we have that
A ⊆ gS (A).
Proof: Let v ∈ A. There exists c ∈ S such that c =
Vor(c + v). In particular, c is one of the vertices of chS.
Therefore, v = (c + v) − c = x − Vor(x) for x = c + v ∈
(chS + A){c} and Vor(x) = c, implying that v ∈ gS (A).

(ii) (Additivity) gS (A ∪ B) = gS (A) ∪ gS (B).
Proof: Property (i) is straightforward. To prove (ii), first
note that by the monotonicity property (i), we have that
gS (A) ⊆ gS (A ∪ B) and gS (B) ⊆ gS (A ∪ B), hence
gS (A) ∪ gS (B) ⊆ gS (A ∪ B). For the other direction,
let v ∈ gS (A ∪ B). Then, v = x − Vor(x) for some
x ∈ chS + A ∪ B = (chS + A) ∪ (chS + B), where the
last equality follows by Proposition 6. Therefore, v ∈ gS (A)
or v ∈ gS (B). This implies that gS (A) ∪ gS (B) ⊇ gS (A ∪ B)
and completes the proof of the Lemma.
Definition 6 (Iterates of gS (A)). In the setting of Definition
2, we let
g0S (A) := A
and
gnS (A) := gS (gn−1
(A)) for n ∈ N, n ≥ 1.
S
Proposition 9 (Monotonicity of the iterates). Let A ⊆ Rd .
Then the iterates gnS (A) are monotonic, in the sense that
0
gnS (A) ⊆ gnS (A) for all n ≤ n0 . Thus, there exists
[
n
g∞
gnS (A).
S (A) := lim gS (A) =
n→∞

Proof: This proposition is a direct consequence of Lemma
3 (i) and Lemma 2.
Theorem 7. Let A ⊆ Rd . The set g∞
S (A) is the minimal Ginvariant set containing the set A.
that

Proof: We first prove the invariance of g∞
S (A). We have

gS (g∞
S (A))

Proposition 8. A ⊆ Rd is G-invariant with respect to S if
and only if
A = gS (A)
Proof: Note that by Lemma 2, we only need to consider
the condition A ⊇ gS (A) for the purpose of the proof.
(⇒) Assume that A ⊇ gS (A). Let x ∈ chS + A. Also, by
the definition of Voronoi cells, x ∈ VS (Vor(x)). Thus,
x ∈ (chS + A){Vor(x)}

and
x − Vor(x) ∈ (chS + A){Vor(x)} − Vor(x) ⊆ gS (A) ⊆ A,
where the last set inclusion follows by the hypothesis. Hence,
by Definition 1, A is an invariant set.
(⇐) Assume that A is an invariant set. Also, assume by the
way of contradiction that A 6⊇ gS (A). Namely, there exists
v ∈ gS (A) such that v ∈
/ A. But every v ∈ gS (A) can be
written as v = x − Vor(x) for some x ∈ chS + A. In other
words, there exists x ∈ chS + A such that x − Vor(x) ∈
/ A, a
contradiction to Definition 1.
Lemma 3 (Properties of gS ).
(i) (Monotonicity) If A ⊆ B then gS (A) ⊆ gS (B).

n≥0

= gS 


[

n≥0

=

[
n≥0

=

[
n≥1

=

A∪

gnS (A)

gS (gnS (A))
gnS (A)
[
n≥1

gnS (A) = g∞
S (A),

where the second equality follows by Lemma 3 (ii). Thus, by
Proposition 8, g∞
S (A) is invariant. It also contains A by its
definition.
To prove minimality, let Q be any invariant set containing
A. Then
Q ⊇ gS (Q) ⊇ gS (A),
where the first inclusion follows by Proposition 8 and the second inclusion follows by the monotonicity property (Lemma
3 (i)). Applying these rules recursively, we obtain that
Q ⊇ gnS (A),

n ≥ 0.

This implies that
Q⊇

N
[
n=0

gnS (A),

N ≥0

13

invariance follows by

and hence
Q ⊇ g∞
S (A).



This completes the proof of the Theorem.
2) Convex Case: We next prove a variant of Theorem 7 for
the case of convex minimal G-invariant sets.

gS (G∞
S (A))

= gS 

n≥0

=

[
n≥0

⊆

Definition 7 (Iterates of GS (A)). In the setting of Definition
3, we let

=

G0S (A) := A

[
n≥0

[
n≥1

= A∪

and


[

GnS (A)

gS (GnS (A))
GS (GnS (A))
GnS (A)
[
n≥1

GnS (A) = G∞
S (A).

GnS (A) := GS (Gn−1
(A)) for n ∈ N, n ≥ 1.
S

To prove minimality, let Q be any invariant convex set
containing A. Then

Proposition 10. A convex set A ⊆ Rd is invariant if and only
if

Q = GS (Q) ⊇ GS (A),

A = GS (A)

where the first inclusion follows by Proposition 10 and the second inclusion follows by the monotonicity property (Lemma
4). Applying these rules recursively, we obtain that

Proof: Observe that by Lemma 2, we have for any A that
A ⊆ gS (A) ⊆ GS (A).
Thus, we next focus on the condition A ⊇ GS (A).
(⇒) Assume that A ⊇ GS (A). Then clearly A ⊇ gS (A),
and by Proposition 8, A is an invariant set.
(⇐) Assume that A is a convex invariant set. Also, assume
by the way of contradiction that A 6⊇ GS (A). Namely, there
exists v ∈ GS (A) such
/ A. But every v ∈ P
GS (A) can
P that v ∈
be written as v = i βi vi with vi ∈ gS (A) and i βi = 1,
β ≥P
0. By the invariance of A, we thus have that vi ∈ A, but
v = i βi vi ∈
/ A, a contradiction to convexity of A.
Lemma 4 (Monotonicity of GS ). If A ⊆ B then GS (A) ⊆
GS (B).
Proof: The result follows by Lemma 3 (i) and the
monotonicity of the convex hull (Lemma 5).
Proposition 11 (Monotonicity of the iterates). Let A ⊆ Rd
be a convex set. Then the iterates GnS (A) are monotonic, in
0
the sense that GnS (A) ⊆ GnS (A) for all n ≤ n0 . Thus, there
exists
[
n
G∞
GnS (A).
S (A) := lim GS (A) =
n→∞

n≥0

Proof: By Lemma 2, we have that A ⊆ GS (A). The result
then follows by applying the monotonicity property of GS
(Lemma 4) for this inclusion recursively, and using Definition
6.
Theorem 8. Let A ⊆ Rd be a convex set. The set G∞
S (A) is
the minimal convex invariant set containing the set A.
Proof: The proof is similar to that of Theorem 7. The

Q ⊇ GnS (A),

n ≥ 0.

This implies that
Q⊇

N
[
n=0

GnS (A),

N ≥0

and hence
Q ⊇ G∞
S (A).
This completes the proof of the Theorem.
3) Multiple-Set Case: Recall Definition 2 of the operator
gS . We next give the proofs of Proposition 1 and Theorem
1 presented in Section IV-A. We also prove the convex
counterpart of Theorem 1, namely Theorem 3 (i). The proof
of boundedness of the minimal invariant sets (namely, proofs
of Theorem 2 and Theorem 3 (ii)) is given in Section VII-B4.
Proof of Proposition 1: By Lemma 2, A ⊆ gS (A) for
every S ∈ S. Hence, A ⊆ gS (A) and part (i) of the proposition
follows.
We next prove the two directions of part (ii).
(⇒) Assume A ⊇ gS (A). Then, also A ⊇ gS (A) for all
S ∈ S. Thus, by Theorem 7, A is invariant for all S ∈ S,
hence jointly invariant.
(⇐) Assume A is invariant for
S all S ∈ S. Then A ⊇ gS (A)
all S ∈ S, therefore also A ⊇ S∈S gS (A) = gS (A).
To prove Theorem 1, we will need the following auxiliary
results.
Lemma 5 (Properties of gS ).
(i) (Monotonicity) If A ⊆ B then gS (A) ⊆ gS (B).
(ii) (Additivity) gS (A ∪ B) = gS (A) ∪ gS (B).
Proof: This result follows trivially by Lemma 3 and by
using the definition of gS as the union of gS over S ∈ S.

14

Proposition 12 (Monotonicity of the iterates). Let A ⊆ Rd .
Then the iterates gnS (A) are monotonic, in the sense that
0
gnS (A) ⊆ gnS (A) for all n ≤ n0 . Thus, there exists
g∞
S (A)

:= lim

n→∞

gnS (A)

=

[

By the definition of g, we have that
[
gS (Qr ) =
(chS + Qr ) ∩ VS (c) − c
c∈S

"

gnS (A).

#
[

=

n≥0

c∈Sc

(chS + Qr ) ∩ VS (c) − c

[

"
Proof: This proposition is a direct consequence of Lemma
5 (i) and Proposition 1 (i).
Proof of Theorem 1: The proof is exactly the same as
that of Theorem 7 by using the properties of g∞
S established
in Lemma 5 and Proposition 12.
Finally, the proof of Theorem 3 (i) is exactly the same as
that of Theorem 8 given the results for the operator gS (A)
proved above.
4) Proof of Boundedness of the Minimal G-Invariant Set:
In this section, we prove Theorem 2 and Theorem 3 (ii). To
that end, we use the results from [2], [13], where it was proven
that there exists a convex bounded invariant set that contains
the origin in the special case where a set S is the collection
of corner points (vertexes) of chS.
For any S ∈ S, consider a partition of S imposed by chS
into two disjoint sets:
1) Sc : the “corner points” of S, namely Sc ⊆ S such that
Sc ⊆ ∂ chS.
2) Snc : the “inner points” of S, namely Snc ⊆ S such that
Snc 6⊆ ∂ chS.
Note that in this definition, we include in the corner points
also any points in S that lie on the edges of chS. It is easy
to see that the Voronoi cells associated with Snc are bounded,
while the Voronoi cells associated with Sc are unbounded.
Moreover, chS = chSc .
We have the following basic result about the Voronoi cells
of S and Sc .
Lemma 6. We have that VS (c) ⊆ VSc (c) for all c ∈ Sc .
Proof: Let c ∈ Sc and x ∈ VS (c). By the definition of
the Voronoi cell
0

kx − ck ≤ kx − c k,

∀c ∈ S.

∀c0 ∈ Sc ,

implying that x ∈ VSc (c).
In [2], [13], it was shown that for the collection Sc =
{Sc , S ∈ S}, under the hypothesis of Theorem 2, for any r ≥ 0
large enough, one can construct a bounded convex invariant set
Qr which contains a ball with radius r centered in the origin.
That is
∀Sc ∈ Sc ,

gSc (Qr ) = Qr .

We next extend this proof to the case where Snc 6= ∅.

c∈Snc

(chS + Qr ) ∩ VS (c) − c

"
⊆

#
[
c∈Sc

(chS + Qr ) ∩ VSc (c) − c

[

"

#
[
c∈Snc

(chS + Qr ) ∩ VS (c) − c
#

"
= gSc (Qr )

[

[
c∈Snc

(chS + Qr ) ∩ VS (c) − c ,

where the set inclusion follows by Lemma 6. Now observe
that there exists r0 ≥ 0 such that for all r ≥ r0 ,
(chS + Qr ) ∩ VS (c) = VS (c),

c ∈ Snc

as the Voronoi cells for c ∈ Snc are bounded. Thus, for r ≥ r0 ,
[
(chS + Qr ) ∩ VS (c) − c
c∈Snc

is a bounded set. Moreover, since the Voronoi cells for c ∈ Sc
are unbounded,
[
(chS + B(0, r)) ∩ VSc (c) − c = Rd .
lim
r→∞

c∈Sc

Hence, there exists r1 ≥ r0 , such that for all r ≥ r1
[
(chS + Qr ) ∩ VS (c) − c
c∈Snc

⊆
⊆

[
c∈Sc

[
c∈Sc

(chS + B(0, r)) ∩ VSc (c) − c
(chS + Qr ) ∩ VSc (c) − c

= gSc (Qr ).

0

However, since Sc ⊆ S, it is also true that
kx − ck ≤ kx − c0 k,

#
[

(10)

(11)

By the uniform boundedness of VS (c) for c ∈ Snc and
S ∈ S, it follows that there exists r1∗ < ∞ such that (11)
holds for all r ≥ r1∗ and all S ∈ S. Thus, for all S ∈ S,
gS (Qr ) ⊆ gSc (Qr ) = Qr , where the last equality follows by
the invariance of Qr for Sc (10). Therefore, for r ≥ r1∗ ,
[
gS (Qr ) =
gS (Qr ) ⊆ Qr ,
S∈S

implying that Qr is a G-invariant set with respect to S by
Proposition 1. Since Qr is a bounded convex set, the results
of both Theorem 2 and Theorem 3 (ii) follow.
C. Proofs for the Case of Persistent Prediction
In this section, we present the proofs for the results of
Section IV-B.

15

1) Single-Set Case: We note that the proofs for minimality
of F -invariant sets for the case of a single set S (namely,
involving the operator pS in Definition 5) are given in [12]
for the case where a set S is the collection of corner points
(vertexes) of chS. This proof extends exactly in the same way
as that of G-invariance provided in Section VII-B, thus we
present the next results without a proof.
Lemma 7 (Properties of pS ; Extension of Lem. 3.11 in [12]).
(i) (Monotonicity) If A ⊆ B then pS (A) ⊆ pS (B).
(ii) (Additivity) pS (A ∪ B) = pS (A) ∪ pS (B).
Proposition 13 (Extension of Lemma 3.2 and Proposition 3.12
in [12]). Any D ⊆ Rd satisfies D ⊆ pS (D). Moreover, D is
F -invariant with respect to S if and only if pS (D) ⊆ D.
Proposition 14 (Monotonicity of the iterates; Extension of
Lemma 3.13 in [12]). Let D ⊆ Rd . Then the iterates pnS (D)
0
are monotonic, in the sense that pnS (D) ⊆ pnS (D) for all
n ≤ n0 . Thus, there exists
[
n
p∞
pnS (D).
S (D) := lim pS (D) =
n→∞

0

the sense that PnS (D) ⊆ PnS (D) for all n ≤ n0 . Thus, there
exists
[
n
P∞
PnS (D).
S (D) := lim PS (D) =
n→∞

Proof: By Proposition 13, we have that D ⊆ PS (D). The
result then follows by applying the monotonicity property of
PS (Lemma 8) for this inclusion recursively.
Theorem 10. Let D ⊆ Rd be a convex set. The set P∞
S (D)
is the minimal convex F -invariant set containing the set D.
Proof: The proof is similar to that of Theorem 8. The
invariance follows by


[

pS (P∞
PnS (D)
S (D)) = pS
n≥0

=
⊆
=

Proof: Observe that by Proposition 13, we have for any
D that
D ⊆ pS (D) ⊆ PS (D).
Thus, we next focus on the condition D ⊇ PS (D).
(⇒) Assume that D ⊇ PS (D). Then clearly D ⊇ pS (D),
and by Proposition 13, D is an F -invariant set.
(⇐) Assume that D is a convex F -invariant set. Also,
assume by the way of contradiction that D 6⊇ PS (D).
Namely, there exists v ∈ PS (D) such
/ D.PBut every
P that v ∈
v ∈ PS (D) can be written as v = i βi vi with i βi = 1,
β ≥ 0, and vi ∈ pS (D) ⊆ D by the F -invariance
of D.
P
Therefore, we have that vi ∈ D, but v =
β
v
∈
/
D, a
i i i
contradiction to convexity of D.
Lemma 8 (Monotonicity of PS ). If A ⊆ B then PS (A) ⊆
PS (B).
Proof: The result follows by Lemma 7 (i) and the
monotonicity of the convex hull (Lemma 5).
Proposition 16 (Monotonicity of the iterates). Let D ⊆ Rd
be a convex set. Then the iterates PnS (D) are monotonic, in

[
n≥0

[
n≥1

Theorem 9 (Extension of Corollary 3.15 in [12]). Let D ⊆
Rd . The set p∞
S (D) is the minimal F -invariant set containing
the set D.

Proposition 15. A convex set D ⊆ Rd is F -invariant if and
only if
D = PS (D)

[
n≥0

n≥0

2) Convex Case: The case of the convex iteration (namely,
the one involving the operator PS in Definition 5) is not
treated in [12]. Hence, we next present the proofs for the
following results, which are similar to the ones given in
Section VII-B2.

n≥0

pS (PnS (D))
PS (PnS (D))
PnS (D)

= D∪

[
n≥1

PnS (D) = P∞
S (D).

To prove minimality, let Q be any invariant convex set
containing D. Then
Q ⊇ PS (Q) ⊇ PS (D),
where the first inclusion follows by Proposition 15 and the second inclusion follows by the monotonicity property (Lemma
8). Applying these rules recursively, we obtain that
Q ⊇ PnS (D),

n ≥ 0.

This implies that
Q⊇

N
[
n=0

PnS (D),

N ≥0

and hence
Q ⊇ P∞
S (D).
This completes the proof of the Theorem.
3) Multiple-Set Case: Recall Definition 5 of the operators
pS and PS . Also, recall the proofs for the case of G-invariance
given in Section VII-B3. Observe that, given the results of
Sections VII-C1 and VII-C2 for a single-set case, the proof
of Proposition 4 is exactly the same as that of Proposition 1.
Also, the proof of Theorem 4 (i) is exactly the same as these
of Theorem 1 and Theorem 3 (i).
4) Proof of Boundedness of the Minimal F -Invariant Set:
The proof of Theorem 4 (ii) follows exactly that of the
boundedness of the minimal G-invariant set provided in full
detail in Section VII-B4. It is omitted here to due to space
constraints.

16

R EFERENCES
[1] 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.
[2] R. L. Adler, B. Kitchens, M. Martens, C. Pugh, M. Shub, and C. Tresser,
“Convex dynamics and applications,” Ergodic Theory and Dynamical
Systems, vol. 25, pp. 321–352, 4 2005.
[3] D. Coppersmith, T. Nowicki, G. Paleologo, C. Tresser, and
C. Wu, “The optimality of the online greedy algorithm in
carpool and chairman assignment problems,” ACM Trans. Algorithms,
vol. 7, no. 3, pp. 37:1–37:22, Jul. 2011. [Online]. Available:
http://doi.acm.org/10.1145/1978782.1978792
[4] F. Blanchini, “Set invariance in control,” Automatica, vol. 35,
no. 11, pp. 1747 – 1767, 1999. [Online]. Available:
http://www.sciencedirect.com/science/article/pii/S0005109899001132
[5] M. S. Branicky, V. S. Borkar, and S. K. Mitter, “A unified framework for
hybrid control: model and optimal control theory,” IEEE Transactions
on Automatic Control, vol. 43, no. 1, pp. 31–45, Jan 1998.
[6] A. Bemporad and M. Morari, “Control of systems integrating logic,
dynamics, and constraints,” Automatica, vol. 35, no. 3, pp. 407–427,
Mar. 1999.
[7] J. Lygeros, C. Tomlin, and S. Sastry, “Controllers for
reachability specifications for hybrid systems,” Automatica,
vol. 35, no. 3, pp. 349 – 370, 1999. [Online]. Available:
http://www.sciencedirect.com/science/article/pii/S0005109898001939
[8] R. W. Floyd and L. Steinberg, “An adaptive algorithm for spatial grey
scale,” in Proc. Dig. SID International Symp., Los Angeles, California,
1975, p. 36–37.
[9] R. Gray, “Oversampled sigma-delta modulation,” Communications,
IEEE Transactions on, vol. 35, no. 5, pp. 481–489, May 1987.
[10] D. Anastassiou, “Error diffusion coding for a/d conversion,” Circuits
and Systems, IEEE Transactions on, vol. 36, no. 9, pp. 1175–1186, Sep
1989.
[11] R. Adams and R. Schreier, “Stability theory for delta-sigma modulators,”
in Delta-Sigma Data Converters: Theory, Design, and Simulation, S. R.
Norsworthy, R. Schreier, and G. C. Temes, Eds. New Jersey: Wiley,
1996.
[12] T. Nowicki and C. Tresser, “Convex dynamics: properties of invariant
sets,” Nonlinearity, vol. 17, no. 5, p. 1645, 2004. [Online]. Available:
http://stacks.iop.org/0951-7715/17/i=5/a=005
[13] C. Tresser, “Bounding the errors for convex dynamics on one or
more polytopes,” Chaos, vol. 17, no. 3, 2007. [Online]. Available:
http://scitation.aip.org/content/aip/journal/chaos/17/3/10.1063/1.2747053
[14] R. Adler, B. Kitchens, M. Martens, A. Nogueira, C. Tresser, and C. W.
Wu, “Error bounds for error diffusion and related digital halftoning
algorithms,” in Circuits and Systems, 2001. ISCAS 2001. The 2001 IEEE
International Symposium on, vol. 2, May 2001, pp. 513–516 vol. 2.
[15] R. L. Adler, B. P. Kitchens, M. Martens, C. P. Tresser, and C. W. Wu,
“The mathematics of halftoning,” IBM J. Res. Dev., vol. 47, no. 1, pp. 5–
15, Jan. 2003. [Online]. Available: http://dx.doi.org/10.1147/rd.471.0005
[16] R. Adler, T. Nowicki, G. Swirszcz, C. Tresser, and S. Winograd,
“Error diffusion on acute simplices: Geometry of invariant sets,” Chaotic
Modeling and Simulation (CMSIM), vol. 1, pp. 3–26, 2015.
[17] “CGAL, Computational Geometry Algorithms Library,” version 4.9,
http://www.cgal.org.
[18] N. J. Bouman, “Convex dynamics research software,” 2016,
https://github.com/niekbouman/convex dynamics.
[19] G. Costanzo, A. Bernstein, L. Reyes Chamorro, H. Bindner, J.-Y.
Le Boudec, and M. Paolone, “Electric Space Heating Scheduling for
Real-time Explicit Power Control in Active Distribution Networks,” in
IEEE PES ISGT Europe 2014, 2014.
[20] L. Reyes-Chamorro, A. Bernstein, J.-Y. Le Boudec, and M. Paolone, “A
composable method for real-time control of active distribution networks
with explicit power setpoints. Part II: Implementation and validation,”
Electric Power Systems Research, vol. 125, pp. 265–280, 2015.
[21] S. Papathanassiou, N. Hatziargyriou, and K. Strunz, “A benchmark low
voltage microgrid network,” in Proceedings of the CIGRÉ Symposium
“Power Systems with Dispersed Generation: technologies, impacts on
development, operation and performances”, Apr. 2005, Athens, Greece.

