Initialization-free distributed coordination for economic
dispatch under varying loads and generator commitment ?
Ashish Cherukuri a

arXiv:1409.4382v1 [math.OC] 15 Sep 2014

a

Jorge Cortés a

Department of Mechanical and Aerospace Engineering, University of California, San Diego, CA, 92093, USA

Abstract
This paper considers the economic dispatch problem for a network of power generating units communicating over a strongly
connected, weight-balanced digraph. The collective aim is to meet a power demand while respecting individual generator
constraints and minimizing the total generation cost. We design a distributed coordination algorithm consisting of two
interconnected dynamical systems. One block uses dynamic average consensus to estimate the evolving mismatch in load
satisfaction given the generation levels of the units. The other block adjusts the generation levels based on the optimization
objective and the estimate of the load mismatch. Our convergence analysis shows that the resulting strategy provably converges
to the solution of the dispatch problem starting from any initial power allocation, and therefore does not require any specific
procedure for initialization. We also characterize the algorithm robustness properties against the addition and deletion of units
(capturing scenarios with intermittent power generation) and its ability to track time-varying loads. Our technical approach
employs a novel refinement of the LaSalle Invariance Principle for differential inclusions, that we also establish and is of
independent interest. Several simulations illustrate our results.
Key words: distributed optimization; economic dispatch; power networks; dynamic average consensus; invariance principles

1

Introduction

The advent of renewable energy sources and their integration into electricity grids is making power generation
and distribution an increasingly decentralized problem.
The large-scale and highly dynamic nature of the resulting grid optimization problems makes traditional centralized, top-down approaches impractical because they
rely on the assumption of a fixed, limited number of generation units. To solve these problems efficiently, there is
a need to design distributed algorithms that can handle
dynamic loads, provide plug-and-play capabilities, are
robust against transmission and generation failures, and
adequately preserve the privacy of the entities involved.
These considerations motivate us to consider here the design of distributed algorithmic solutions to the economic
dispatch (ED) problem, where a group of power generators aims to meet a power demand while minimizing
the total generation cost (the summation of individual
costs) and respecting the individual generators’ capacity
constraints. We are interested in the synthesis of strate? A preliminary version of this work appears at the 2014
Allerton Conference on Communication, Control, and Computing, Monticello, Illinois.
Email addresses: acheruku@ucsd.edu (Ashish
Cherukuri), cortes@ucsd.edu (Jorge Cortés).

Preprint submitted to Automatica

gies that solve the ED problem starting from any initial
power allocation, can handle time-varying loads, and are
robust against intermittent power generation caused by
unit addition and deletion.
Literature review: The ED problem has been traditionally solved in a centralized manner, see e.g. [Chowdhury and Rahman, 1990] and references therein. Since
distributed algorithmic solutions to grid optimization
problems are envisioned as part of the future smart
grid [Farhangi, 2010], this has motivated the appearance
of a number of distributed strategies for the ED problem in the literature. While there exists a broad variety
in the assumptions made, a majority of the works rely
on the specific form of the solutions of the optimization
problem and propose consensus-based algorithms. Many
works consider convex, quadratic cost functions for the
power generators and perform consensus over their incremental costs under undirected [Zhang and Chow,
2012, Kar and Hug, 2012] or directed [DominguezGarcia et al., 2012, Binetti et al., 2014a] communication topologies. Some works consider general convex
cost functions, like we do here, but either do not consider capacity constraints on the generators [Mudumbai
et al., 2012], assume the initial power allocation to meet
the total load [Cherukuri and Cortés, 2013, Pantoja
et al., 2014], or require feedback on the power mismatch

5 November 2021

gradient” dynamics, that solves the ED problem starting from any initial power allocation. This strategy has
two components: one component seeks to optimize the
network generation cost while keeping constant the total power generated; the other component is a feedback
correction term driven by the error between the desired
total load and the network generation. This latter term
is responsible for ensuring that the algorithm trajectories asymptotically satisfy the load satisfaction constraint irrespective of the initial power allocation. These
observations set the basis for our second contribution,
which is the synthesis of a distributed coordination algorithm, termed “dynamic average consensus + Laplaciannonsmooth-gradient” dynamics, with the same convergence guarantees. Our design consists of two coupled
dynamical systems: a dynamic average consensus algorithm to estimate the mismatch between generation and
desired load in a distributed fashion and a distributed
Laplacian-nonsmooth-gradient dynamics that employs
these estimates to dynamically allocate the unit generation levels. The convergence analysis of both the centralized and distributed algorithms relies on a combination
of tools from algebraic graph theory, nonsmooth analysis, set-valued dynamical systems, and dynamic average consensus, and most notably on a refined version
of the LaSalle Invariance Principle for differential inclusions, which constitutes our third contribution. Roughly
speaking, the application of the known LaSalle Invariance Principle would only establish asymptotic convergence towards the network satisfaction of the total load.
Instead, the use of the refined version allows us, for each
algorithm, to establish global convergence of the trajectories to the solutions of the ED problem. Our final contribution is the formal characterization of the robustness properties of the distributed algorithm. Building on
the observation that the mismatch dynamics between
network generation and total load is exponentially convergent and input-to-state stable, we establish the algorithm ability to track time-varying loads and its robustness in scenarios with intermittent power generation.

from the shift in frequency due to primary droop control [Zhang et al., 2014]. Along with load and capacity
constraints, Binetti et al. [2014a], Loia and Vaccaro
[2013] consider transmission losses, and Binetti et al.
[2014b] additionally take into account valve-point loading effects and prohibited operating zones. However,
these constraints make the problem nonconvex and
prevent these works from obtaining theoretical guarantees on the algorithm convergence properties. In [Du
et al., 2012], the authors propose best-response dynamics for a potential-game formulation of the nonconvex
ED problem, but the implementation requires all-to-all
communication among the generators. Xiao and Boyd
[2006], Johansson and Johansson [2009] introduce distributed methods to solve resource allocation problem
very similar to the ED problem, but without taking
into account individual agent constraints. Instead,
these are incorporated in the formulation of Simonetto et al. [2012], but the proposed algorithm arrives
at suboptimal solutions of the optimization problem.
Our algorithm design and analysis rely on dynamic average consensus and differential inclusions. In dynamic
average consensus, see e.g. [Freeman et al., 2006, Kia
et al., 2014] and references therein, each agent has access to a time-varying input signal and interacts with
its neighbors in order to track the average of the input
signals across the network. We build on our previous
work [Cherukuri and Cortés, 2013], which requires a
proper algorithm initialization, and employ tools from
dynamic average consensus to synthesize a coordination strategy that converges from any initial condition.
Regarding analysis, our technical approach builds on
Lyapunov stability tools for differential inclusions and
nonsmooth systems, see e.g. [Bacciotti and Ceragioli,
1999, Cortés, 2008] and references therein. Of particular
importance is the work [Arsie and Ebenbauer, 2010] for
differential equations, that provides a way to further
refine the description of omega-limit sets of trajectories
by employing more than one LaSalle-type function.
Statement of contributions: We start with the formal definition of the ED problem for a network of power generators communicating over a strongly connected, weightbalanced digraph. The optimization problem is convex
as the individual cost functions are smooth and convex,
the load satisfaction is a linear constraint, and the individual generators’ capacities prescribe convex inequality
constraints. Our formulation is a simplification of the ED
problem in its full generality, which in practice may have
additional constraints (e.g., transmission losses, line capacity constraints, valve-point loading effects, ramp rate
limits, prohibited operating zones) that make it nonconvex. However, our developments show that obtaining a
provably correct algorithmic solution for the formulation here of the ED problem given our performance requirements (distributed, convergent irrespective of initial condition, able to handle time-varying loads, and
robust to intermittent power generation) is challenging.
Our first contribution is the design of a centralized algorithm, termed “load mismatch + Laplacian-nonsmooth-

2

Preliminaries

This section introduces basic concepts and preliminaries.
We begin with some notational conventions. Let R, R≥0 ,
R>0 , Z≥1 denote the real, nonnegative real, positive real,
and positive integer numbers, resp. For r ∈ R we denote
n
Hr = {x ∈ Rn | 1>
n x = r}. The 2- and ∞-norms on R
and their respective induced norms on Rn×n are denoted
with k·k and k·k∞ , resp. We let B(x, δ) = {y ∈ Rn | ky−
xk < δ}. For D ⊂ Rn , D denotes its closure. For x ∈ Rn ,
xi ∈ R denotes its i-th component. Given vectors x, y ∈
Rn , x ≤ y if and only if xi ≤ yi for all i ∈ {1, . . . , n}.
We denote 1n = (1, . . . , 1) ∈ Rn .A set-valued map f :
Rn ⇒ Rm associates to each point in Rn a set in Rm . For
a symmetric matrix A ∈ Rn×n , λmin (A) and λmax (A)
denote the minimum and maximum eigenvalues of A.
Finally, we let [u]+ = max{0, u} for u ∈ R.
2

2.1

Lx ky −y 0 k, for all y, y 0 ∈ B(x, ). A function f : Rn → R
is regular at x ∈ Rn if, for all v ∈ Rn , the right and generalized directional derivatives of f at x in the direction
of v coincide, see [Cortés, 2008] for definitions of these
notions. A function that is continuously differentiable at
x is regular at x. Also, a convex function is regular. A
set-valued map H : Rn ⇒ Rn is upper semicontinuous at
x ∈ Rn if, for all  ∈ (0, ∞), there exists δ ∈ (0, ∞) such
that H(y) ⊂ H(x) + B(0, ) for all y ∈ B(x, δ). Also, H
is locally bounded at x ∈ Rn if there exist , δ ∈ (0, ∞)
such that kzk ≤  for all z ∈ H(y) and y ∈ B(x, δ).
Given a locally Lipschitz function f : Rn → R, let Ωf
be the set (of measure zero) of points where f is not
differentiable. The generalized gradient ∂f : Rn ⇒ Rn is

Graph theory

We present basic notions from algebraic graph theory following [Bullo et al., 2009]. A directed graph (or digraph)
is a pair G = (V, E), with V = {1, . . . , n} the vertex set
and E ⊆ V × V the edge set. A path is a sequence of vertices connected by edges. A digraph is strongly connected
if there is a path between any pair of vertices. The sets
of out- and in-neighbors of v are, resp., N out (v) = {w ∈
V | (v, w) ∈ E} and N in (v) = {w ∈ V | (w, v) ∈ E}. A
weighted digraph G = (V, E, A) is composed of a digraph
(V, E) and an adjacency matrix A ∈ Rn×n
≥0 with aij > 0
iff (i, j) ∈ E. ThePweighted out- and in-degree
Pn of i are,
n
in
resp., dout (i) =
j=1 aij and d (i) =
j=1 aji . The
Laplacian matrix is L = Dout −A, where Dout is the diagonal matrix with (Dout )ii = dout (i), for all i ∈ {1, . . . , n}.
Note that L1n = 0. If G is strongly connected, then 0 is
a simple eigenvalue of L. G is undirected if L = L> . G
is weight-balanced if dout (v) = din (v), for all v ∈ V iff
>
1>
n L = 0 iff L + L ≥ 0. Note that any undirected graph
is weight-balanced. If G is weight-balanced and strongly
connected, then 0 is a simple eigenvalue of L + L> . In
such case, one has for x ∈ Rn ,
x> (L + L> )x ≥ λ2 (L + L> ) x −

2
1 >
(1n x)1n ,
n

∂f (x) = co{ lim ∇f (xi ) | xi → x, xi ∈
/ S ∪ Ωf },
i→∞

where co denotes convex hull and S ⊂ Rn is any set
of measure zero. The map ∂f is locally bounded, upper semicontinuous, and takes non-empty, compact, and
convex values. A critical point x of f satisfies 0 ∈ ∂f (x).
Given a set-valued map H : Rn ⇒ Rn , a differential
inclusion on Rn is
ẋ ∈ H(x).
(3)
A solution of (3) on [0, T ] ⊂ R is an absolutely continuous map x : [0, T ] → Rn that satisfies (3) for almost all
t ∈ [0, T ]. If H is locally bounded, upper semicontinuous,
and takes non-empty, compact, and convex values, then
existence of solutions is guaranteed. The set of equilibria
of (3) is Eq(H) = {x ∈ Rn | 0 ∈ H(x)}. A set S ⊂ Rn is
weakly (resp., strongly) positively invariant under (3) if,
for each x ∈ S, at least a solution (resp., all solutions)
starting from x is (resp., are) entirely contained in S.
For dynamics with uniqueness of solution, both notions
coincide and are referred as positively invariant. Given
a locally Lipschitz function f : Rn → R, the set-valued
Lie derivative LH f : Rn ⇒ R of f with respect to (3) is

(1)

with λ2 (L+L>) the smallest non-zero eigenvalue of L+L>.
2.2

Dynamic average consensus

Here, we introduce notions on dynamic average consensus following [Kia et al., 2014]. Consider n ∈ Z≥1 agents
communicating over a strongly connected, weightbalanced digraph G whose Laplacian is denoted as L.
Each agent is associated with a state xi ∈ R and an input signal t 7→ ui (t) ⊂ R that is measurable and locally
essentially bounded. The aim is to provide a distributed
dynamics such that thePstate of each agent xi (t) tracks
n
the average signal n1 i=1 ui (t) asymptotically. This
can be achieved via the dynamics Xdac : R2n → R2n ,

LH f (x) = {a ∈ R | ∃v ∈ H(x) s.t. ζ > v = a for all
ζ ∈ ∂f (x)}.

ẋ = −αx − βLx − v + νu,
v̇ = αβLx,

For a trajectory t 7→ ϕ(t), ϕ(0) ∈ Rn of (3), the evolution
of f along it satisfies

where α, β, ν > 0 are design parameters and v ∈ Rn
is an auxiliary state. If the initial condition satisfies 1>
n v(0) = 0 and the time-derivatives of the input signals are bounded, then one can show, cf. [Kia
et al., 2014, Corollary
4.1], that the error signal
Pn
t 7→ xi (t) − n1 i=1 ui (t) is ultimately bounded for
each i ∈ {1, . . . , n}. Moreover, this error vanishes if the
input signal converges to a constant value.
2.3

d
f (ϕ(t)) ∈ LH f (ϕ(t))
dt
for almost all t ≥ 0. The omega-limit set of the trajectory, denoted Ω(ϕ), is the set of all points y ∈ Rn for
which there exists a sequence {tk }∞
k=1 with tk → ∞ and
limk→∞ ϕ(tk ) = y. If the trajectory is bounded, then the
omega-limit set is nonempty, compact, connected, and
weakly invariant. These tools allow us to characterize
the asymptotic behavior of solutions of differential inclusions. In Appendix A we develop a novel refinement
of the LaSalle Invariance Principle for differential inclusions, see e.g., [Cortés, 2008], which is suitable for the

Nonsmooth analysis and differential inclusions

We review here some notions from nonsmooth analysis and differential inclusions following [Cortés, 2008]. A
function f : Rn → Rm is locally Lipschitz at x ∈ Rn
if there exist Lx ,  ∈ (0, ∞) such that kf (y) − f (y 0 )k ≤
3

analysis of the coordination algorithms.
3

Our design strategy relies on the following reformulation of the ED problem without inequality constraints.
Consider the modified ED problem

Problem statement

This section presents the network model and the economic dispatch problem we set out to solve in a distributed and robust fashion. Consider n ∈ Z≥1 power
generators communicating over a strongly connected
and weight-balanced digraph G = (V, E, A). Each generator corresponds to a vertex in the digraph and an
edge (i, j) represents the ability of generator j to send
information to generator i. The cost of power generation for unit i is measured by fi : R → R≥0 , assumed to
be convex and continuously differentiable. Representing the power generated by unit i by Pi ∈ R, the total
cost incurred by the network with the power allocation
P = (P1 , . . . , Pn ) ∈ Rn is measured by f : Rn → R≥0 as
f (P ) =

n
X

minimize

f  (P ),

(5a)

subject to

1>
n P = Pl ,

(5b)

where the objective function is
f  (P ) =

n
1 X
fi (Pi ) + ( ([Pi − PiM ]+ + [Pim − Pi ]+ )).
 i=1
i=1

n
X

This corresponds to each generator i ∈ {1, . . . , n} having
the modified local cost
1
fi (Pi ) = fi (Pi ) + ([Pi − PiM ]+ + [Pim − Pi ]+ ).


fi (Pi ).

Note that fi is convex, locally Lipschitz, and continuously differentiable on R except at Pi = Pim and
Pi = PiM . Moreover, the total cost f  is convex, locally Lipschitz, and regular. According to our previous
work [Cherukuri and Cortés, 2013, Proposition 5.2], the
solutions to the original (4) and the modified (5) ED
problems coincide for  ∈ R>0 such that

i=1

Note that f is convex and continuously differentiable.
The generators aim to minimize the total cost P
f (P ) while
n
meeting the total power load Pl ∈ R>0 , i.e., i=1 Pi =
Pl . Each generator has an upper and a lower limit on the
power it can produce, Pim ≤ Pi ≤ PiM for i ∈ {1, . . . , n}.
Formally, the economic dispatch (ED) problem is

<
minimize

f (P ),

(4a)

subject to 1>
n P = Pl ,
P

m

≤P ≤P

(4b)
M

.

1
2 maxP ∈FED k∇f (P )k∞

.

(6)

Throughout the paper, we assume the parameter  satisfies this condition. A useful fact is that P ∗ ∈ Rn is a
solution of (5) if and only if there exists µ ∈ R such that

(4c)

The constraint (4b) is the load condition and (4c) are the
box constraints. The set of allocations satisfying the box
constraints is FB = {P ∈ Rn | P m ≤ P ≤ P M }. Further, we denote the feasibility set of (4) as FED = FB ∩
HPl = {P ∈ Rn | P m ≤ P ≤ P M and 1>
n P = Pl } and
∗
∗
the set of solutions as FED
. Since FED is compact, FED
is compact. Note that P M ∈ FED implies FED = {P M }.
Similarly P m ∈ FED implies FED = {P m }. Therefore,
we assume P M and P m are not feasible.
Our objective is to design a distributed coordination algorithm that allows the team of generators to solve the
ED problem (4) starting from any initial condition, can
handle time-varying loads, and is robust to intermittent
power generation.
Remark 3.1 (Additional practical constraints): We do
not consider here, for simplicity, other constraints on
the ED problem such as transmission losses, transmission line capacities, valve-point loading effects, ramp
rate limits, and prohibited operating zones. As our
forthcoming treatment will show, the design and analysis of algorithmic solutions to the ED problem without
these additional constraints is already quite challenging
given our performance requirements. Nevertheless, Remark 5.4 later comments on how to adapt our algorithm
to deal with more general scenarios.
•

µ1n ∈ ∂f  (P ∗ )
4

and

∗
1>
n P = Pl .

(7)

Robust centralized algorithmic solution

This section presents a robust strategy to make the
network power allocation converge to the solution set
of the ED problem starting from any initial condition.
Even though this algorithm is centralized, its design
provides enough insight to tackle later the design of
a distributed algorithmic solution. Consider the “load
mismatch + Laplacian-nonsmooth-gradient” (abbreviated lm+L∂) dynamics, represented by the set-valued
map Xlm+L∂ : Rn ⇒ Rn ,
Ṗ ∈ −L∂f  (P ) +

1
(Pl − 1>
n P )1n ,
n

(8)

where L is the Laplacian associated to the strongly connected and weight-balanced communication digraph G.
For each generator, the first term seeks to minimize the
total cost while leaving unchanged the total generated
power. The second term is a feedback element that seeks
to drive the units towards the satisfaction of the load.
The first term is computable using information from its

4

f  (P (t)) → ∞ (as f  is radially unbounded). Therefore,
there exist a sequence of times {tk }∞
k=1 with tk → ∞
such that for all k ∈ Z≥1 ,

neighbors but the second term requires them to know
the aggregated state of the whole network, which makes
it not directly implementable in a distributed manner.
The next result shows that the trajectories of (8) converge to the set of solutions of the ED problem.
Theorem 4.1 (Convergence of the trajectories of
Xlm+L∂ to the solutions of ED problem): The trajectories
of (8) starting from any point in Rn converge to the set
of solutions of (4).

1>
n P (tk ) − Pl <

This implies that there exists a sequence {ζk }∞
k=1 with
ζk ∈ ∂f  (P (tk )) such that, for all k ∈ Z≥1 ,
1
>
− ζk> Lζk + (Pl − 1>
n P (tk ))(1n ζk ) > 0
n
 L + L> 
1
⇒ −ζk>
ζk +
1 > ζk > 0
(10)
2
nk n
>
2
1
λ2 (L + L )
1
ζk − (1>
⇒−
1> ζk > 0,
+
n ζk )1n
2
n
nk n

PROOF. Our proof strategy proceeds by applying the
refined LaSalle Invariance Principle for differential inclusions established in Appendix A, cf. Proposition A.1.
Consider the following function V1 : Rn → R≥0 ,
V1 (P ) =

1
and max LXlm+L∂ f  (P (tk )) > 0. (9)
k

1
2
(Pl − 1>
n P) .
2

where we have used (9) in the first implication and (1)
in the second. Next, we consider two cases depending on
>
whether (a) 1>
n ζk is bounded or (b) 1n ζk → ∞. In
case (a), taking the limit k → ∞ in the last inequality
of (10), we get

The set-valued Lie derivative of V1 along Xlm+L∂ is
2
LXlm+L∂ V1 (P ) = {−(Pl − 1>
n P ) } = {−2V1 (P )}.

Thus, starting at any P (0) ∈ Rn , the trajectory of Xlm+L∂
satisfies V1 (P (t)) = V1 (P (0))e−2t and its omega-limit
set (provided the trajectory is bounded, a fact that we
assume is true for now and establish later) is contained
in HPl . In the notation of Proposition A.1, HPl plays the
role of the closed submanifold S of Rn . We next show
that the hypotheses of this result hold. In the notation of
the Lemma A.1, the function f  , the map (P, ζ) 7→ −Lζ,
and the set-valued map P ⇒ −L∂f  (P ) play the role of
W , g, and F , respectively (our choice of F is because
the dynamics Xlm+L∂ takes the form Ṗ ∈ −L∂f  (P )
on S = HPl ). Notice that ζ 7→ −Lζ is a continuous map
and, since G is strongly connected and weight-balanced,
we have ζ > (−Lζ) = − 12 ζ > (L + L> )ζ ≤ 0 for any ζ ∈
∂f  (P ). Therefore, Lemma A.1(i) is satisfied. Moreover,
if ζ > (−Lζ) = 0 for some ζ ∈ ∂f  (P ), then ζ ∈ span{1n }.
Since for P ∈ HPl , we have

lim

k→∞

ζk −

1 >
(1 ζk )1n = 0.
n n

(11)

Since, kP (t)k → ∞ and 1>
n P (t) → Pl , there exist i, j ∈
{1, . . . , n} such that Pi (tk ) → ∞ and Pj (tk ) → −∞.
∗
Let P ∗ ∈ FED
and µ1n ∈ ∂f  (P ∗ ) for some µ ∈ R.
Then, without loss of generality, we assume that Pi∗ ≤
Pi (tk ) ≤ Pi (tk+1 ) and Pj∗ ≥ Pj (tk ) ≥ Pj (tk+1 ) for all k.
This fact along with the expression of ∂fi : R ⇒ R,

1

{∇fi (Pi ) −  }


1

[∇fi (Pi ) −  , ∇fi (Pi )]

∂fi (Pi ) = {∇fi (Pi )}


[∇fi (Pi ), ∇fi (Pi ) + 1 ]


{∇f (P ) + 1 }
i
i


if Pi < Pim ,
if Pi = Pim ,
if Pim < Pi < PiM ,
if Pi = PiM ,
if Pi > PiM .

gives us the following property for all k ∈ Z≥1 ,

LXlm+L∂ f  (P ) = {−ζ > Lζ | ζ ∈ ∂f  (P )},

min ∂fi (Pi (tk )) ≥ µ, max ∂fj (Pj (tk )) ≤ µ,
min ∂fi (Pi (tk+1 )) ≥ max ∂fi (Pi (tk )),
max ∂fj (Pj (tk+1 )) ≤ min ∂fj (Pj (tk )).

we deduce 0 ∈ LXlm+L∂ f  (P ), i.e., Lemma A.1(ii) holds.
The application of Lemma A.1 then yields that Proposition A.1(ii) holds too. In addition, from the above analysis, note that if 0 ∈ LXlm+L∂ f  (P ) for some P ∈ HPl , then
there exists µ ∈ R such that µ1n ∈ ∂f  (P ) and, from (7),
P is a solution of (4). Therefore, {P ∈ HPl | 0 ∈
LXlm+L∂ f  (P )} is the set of solutions of the ED problem
and belongs to a level set of f  , which establishes that
Proposition A.1(i) also holds.
To be able to apply Proposition A.1 and conclude the
proof, it remains to show that the trajectories of Xlm+L∂
are bounded. We reason by contradiction, i.e., assume
there exists a trajectory t 7→ P (t), P (0) ∈ Rn of Xlm+L∂
such that kP (t)k → ∞. From the analysis above, we
know that along this trajectory 1>
n P (t) → Pl and

(12a)
(12b)
(12c)

Note that the limit (11) yields limk→∞ |(ζk )i − (ζk )j | =
0. On the other hand, from (12b)-(12c), we obtain
|(ζk )i − (ζk )j | ≤ |(ζk+1 )i − (ζk+1 )j | for all k. Therefore,
we obtain (ζk )i = (ζk )j for all k and from (12a), we
get µ = (ζk )i = (ζk )j for all k. From (12b)-(12c), this
further implies that µ ∈ ∂fi (x) for all x ∈ [Pi∗ , ∞) and
that µ ∈ ∂fj (x) for all x ∈ (−∞, Pj∗ ]. Using this fact,
one can construct an unbounded set of solutions to the
ED problem in the following manner. First, fix all the
components of P ∗ except i and j. Now pick any x ∈ R≥0

5

and consider Pi∗ + x and Pj∗ − x. From what we have reasoned so far, all such points that we obtain by varying
x are solutions to the ED problem as they satisfy (7).
∗
This contradicts the fact that FED
is bounded.
In case (b), assume without loss of generality that
>
1>
n ζk → ∞ (the argument for 1n ζk → −∞ follows similarly). As reasoned above, there exists j ∈ {1, . . . , n}
such that Pj (tk ) → −∞ and there exists µ ∈ R such
that (ζk )j ≤ µ for all k ∈ Z≥1 . Using this fact, we upper
bound the left hand side of the inequality (10) by
>

that this dynamics is distributed over the communication graph G. For each i ∈ {1, . . . , n}, zi plays the
role of an estimator associated to i which aims to track
the average signal t 7→ n1 (Pl − 1>
n P (t)). This observation justifies substituting the feedback term in (8) by
z ∈ Rn , giving rise to the “dynamic average consensus +
Laplacian-nonsmooth-gradient” dynamics, abbreviated
dac+L∂ for convenience, mathematically represented by
the set-valued map Xdac+L∂ : R3n ⇒ R3n ,
Ṗ ∈ −L∂f  (P ) + ν1 z,
ż = −αz − βLz − v + ν2 (Pl er − P ),
v̇ = αβLz,

2

1
1 >
λ2 (L + L )
ζk − (1>
ζk )1n +
(1 ζk )
2
n n
nk n
2
1
λ2 (L + L> ) 
1 >
(ζk )j − (1>
≤−
(1 ζk )
+
n ζk )
2
n
nk n

> 
2
1
λ2 (L + L )
1 >
µ − (1>
≤−
ζk ) +
(1 ζk ), (13)
2
n n
nk n

−

where ν1 > 0 is a design parameter. Unlike (8), this dynamics is distributed, as each agent only needs to interact with its neighbors to implement it.
5.1

where the last inequality is valid for all but a finite
number of k. Hence, as 1>
n ζk → ∞, there is k̄ ∈ Z≥1
such that the expression in (13) is negative for k ≥ k̄,
contradicting (10). Thus, we conclude the trajectories
are bounded. 2

Convergence analysis

Here we characterize the asymptotic convergence properties of the dac+L∂ dynamics. We start by establishing
an important fact on the omega-limit set of any trajectory of (15) with initial condition in Rn × Rn × H0 .
Lemma 5.1 (Characterizing the omega-limit set of
the trajectories of the dac+L∂ dynamics): The omegalimit set of any trajectory of (15) with initial condition
(P0 , z0 , v0 ) ∈ Rn ×Rn ×H0 is contained in HPl ×H0 ×H0 .

From the proof above, it is interesting to note that the
feedback term (8) drives the mismatch between generation and load to zero at an exponential rate, no matter
what the initial power allocation. This is a good indication of its robustness properties: time-varying loads or
scenarios with generators going down and coming back
online can be handled as long as the rate of these changes
is lower than the exponential rate of convergence associated to the load satisfaction. We provide a formal characterization of these properties for the distributed implementation of this strategy in the next section.
5

(15a)
(15b)
(15c)

PROOF. From (15c), note that 1>
n v̇ = 0. Since v0 ∈
H0 , this implies that 1>
n v(t) = 0 for all t ≥ 0. Now,
define ζ(t) = 1>
n P (t) − Pl and note that
>
ζ̇(t) = 1>
n Ṗ (t) = ν1 1n z(t),

Robust distributed algorithmic solution

where we have used (15a), and

This section presents a distributed strategy to solve
the ED problem starting from any initial power allocation. We build on the centralized design presented in
Section 4. We also formally characterize the robustness
properties against addition and deletion of generators
and time-varying loads.
Given the discussion on the centralized nature of the
dynamics (8), the core idea of our design is to employ a
dynamic average consensus algorithm that allows each
unit in the network to estimate the mismatch in load
satisfaction. To this end, we assume the total load Pl is
only known to one generator r ∈ {1, . . . , n} (its specific
identity is arbitrary). Following Section 2.2, consider the
dynamics,

ζ̈(t) = ν1 1>
n ż(t)
= ν1 1>
n (−αz(t) − βLz(t) − v(t) + ν2 (Pl ek − P (t))
= −ν1 α(1>
n z(t)) − ν1 ν2 ζ(t) = −αζ̇(t) − ν1 ν2 ζ(t),

where we have used (15b). We write this system as a
first-order one by defining x1 = ζ and x2 = ζ̇ to get
"

ẋ1
ẋ2

#

"
=

1

#" #
x1

−ν1 ν2 −α

x2

0

.

(16)

Evaluating the Lie derivative of the positive definite,
radially unbounded function V2 (x1 , x2 ) = ν1 ν2 x21 + x22
along the above dynamics and applying the LaSalle Invariance Principle [Khalil, 2002], we deduce that x1 (t) →
0 and x2 (t) → 0 as t → ∞, that is, 1>
n P (t) → Pl and
1>
z(t)
→
0.
Since
the
system
(16)
is
linear,
the convern
gence is exponential. 2

ż = −αz − βLz − v + ν2 (Pl er − P ),
v̇ = αβLz,
where er ∈ Rn is the unit vector along the r-th direction and α, β, ν2 > 0 are design parameters. Note

6

tory t 7→ (P (t), ξ1 (t), ξ2 (t)) starting in D(Rn × Rn × H0 )
belongs to HPl × H0 × H0 .
Our next step is to show that the hypotheses of Proposition A.1 are satisfied where HPl × H0 × H0 and V3 play
the role of the closed submanifold S of R3n and the function W , respectively. To do so, we resort to Lemma A.1.
Define the continuous function g : R3n × R3n → R3n by

The next result builds on this fact and Proposition A.1
to establish that the trajectory of power allocations under (15) converges to the solution set of the ED problem.
Theorem 5.2 (Convergence of the dac+L∂ dynamics to
the solutions of ED problem): For α, β, ν1 , ν2 > 0 with
ν 2 λmax (L> L)
ν1
+ 2
< λ2 (L + L> ),
>
βν2 λ2 (L + L )
2α

(17)

g(P, ξ1 , ξ2 , ζ̂) = (−Lζ̂1 + ν1 ξ1 , −βLξ1 − ξ2 ,

the trajectories of (15) starting from any point in Rn ×
∗
∗
Rn × H0 converge to the set Faug
= {(P, z, v) ∈ FED
×
n
{0} × R | v = ν2 (Pl er − P )}.

ν1 ν2 ξ1 − αξ2 − ν2 Lζ̂1 ),

and note that the dynamics (19) can be expressed as
Xdac+L∂ (P, ξ1 , ξ2 ) = {g(P, ξ1 , ξ2 , ζ̂) | ζ̂ ∈ ∂V3 (P, ξ1 , ξ2 )}.
For (P, ξ1 , ξ2 ) ∈ HPl × H0 × H0 and ζ̂ ∈ ∂V3 (P, ξ1 , ξ2 ),

PROOF. Our proof strategy is based on the refined
LaSalle Invariance Principle for differential inclusions
established in Appendix A, cf. Proposition A.1. Before
justifying that all its hypotheses are satisfied, we reformulate the expression for the dynamics to help simplify
the analysis. Consider first the change of coordinates,
(P, z, v) 7→ (P, z, v̄), with v̄ = v − ν2 (Pl er − P ). The
set-valued map Xdac+L∂ then takes the form

ζ̂ > g(P, ξ1 , ξ2 , ζ̂) = −ζ > Lζ + ν1 ζ > ξ1 − βν1 ν2 ξ1> Lξ1
− αkξ2 k2 − ν2 ξ2> Lζ,

where we have used that ζ = ζ̂1 ∈ ∂f  (P ), ζ̂2 = ν1 ν2 ξ1 ,
and ζ̂3 = ξ2 . Since the digraph G is strongly connected
and weight-balanced, we apply (1) and the fact that
1>
n ξ1 = 0 to bound the above expression as

Xdac+L∂ (P, z, v̄) = {(−Lζ + ν1 z, −(α + βL)z − v̄,
(αβL + ν1 ν2 )z − ν2 Lζ) ∈ R3n | ζ ∈ ∂f  (P )}.

1
1
− λ2 (L+L> )kηk2 +ν1 η > ξ1 − βν1 ν2 λ2 (L+L> )kξ1 k2
2
2
− αkξ2 k2 − ν2 ξ2> Lη = γ > M γ,

The change of coordinates shifts the equilibrium of the
consensus dynamics to the origin. Under the additional
change of coordinates (P, z, v̄) 7→ (P, ξ1 , ξ2 ), with
" #
ξ1
ξ2

"
=

I 0

#" #
z

αI I

v̄

,

(20)

>
> > >
where η = ζ − n1 (1>
n ζ)1n , γ = [η , ξ1 , ξ2 ], and

(18)



1
− 12 ν2 L>
− 21 λ2 (L + L> )I
2 ν1 I


1
.
M =
− 12 βν1 ν2 λ2 (L + L> )I
0


2 ν1 I
1
− 2 ν2 L
0
−αI

the set-valued map Xdac+L∂ takes the form
Xdac+L∂ (P, ξ1 , ξ2 ) = {(−Lζ + ν1 ξ1 , −βLξ1 − ξ2 , (19)
ν1 ν2 ξ1 − αξ2 − ν2 Lζ) ∈ R3n | ζ ∈ ∂f  (P )}.

Reasoning with the Schur complement [Boyd and Vandenberghe, 2009], M ∈ R3n×3n is negative definite if

This extra change of coordinates makes it easier to identify the candidate Lyapunov function V3 : R3n → R≥0 ,

1
− λ2 (L + L> )I−
2
"
#−1 "
#
1
i − 1 βν ν λ (L + L> )I 0
h
ν
I
1
2
2
1
2
2
1
1
>
2 ν1 I − 2 ν2 L
0
−αI
− 12 ν2 L

1
V3 (P, ξ1 , ξ2 ) = f  (P ) + (ν1 ν2 kξ1 k2 + kξ2 k2 ).
2
For convenience, denote the overall change of coordinates by D : R3n → R3n ,

1
ν1
ν2
= − λ2 (L + L> )I +
I + 2 L> L,
>
2
2βν2 λ2 (L + L )
4α

(P, ξ1 , ξ2 ) = D(P, z, v) = (P, z, v + αz − ν2 (Pl er − P )).

is negative definite. This latter fact is implied by (17).
As a consequence, ζ̂ > g(P, ξ1 , ξ2 , ζ̂) ≤ 0 and so,
Lemma A.1(i) holds. Moreover, ζ̂ > g(P, ξ1 , ξ2 , ζ̂) = 0 if
and only if η = ξ1 = ξ2 = 0, which means ζ ∈ span{1n }.
Using this fact along with the definition of the set-valued
Lie derivative and the characterization of optimizers (7), we deduce that ζ̂ > g(P, ξ1 , ξ2 , ζ̂) = 0 if and only
if (a) 0 ∈ LXdac+L∂ V3 (P, ξ1 , ξ2 ) and (b) P is a solution of

Our analysis now focuses on proving that, in the new
coordinates, the trajectories of (15) converge to the set
∗
∗
F aug = D(Faug
) = FED
× {0} × {0}.

Note that D(HPl ×H0 ×H0 ) = HPl ×H0 ×H0 and therefore, from Lemma 5.1, the omega-limit set of a trajec-

7

Next, we consider two cases, depending on whether the
sequence {P (tk )} is (a) bounded or (b) unbounded.
In case (a), the sequence {(ξ1 (tk ), ξ2 (tk ))} must be
unbounded. Since M is negative definite, we have
γk> M γk ≤ λmax (M )k(ξ1 (tk ), ξ2 (tk ))k2 . Thus, (22) implies that

the ED problem. Fact (a) implies that Lemma A.1(ii)
holds and hence, Proposition A.1(ii) holds too. Fact
(b) implies that over the set HPl × H0 × H0 , we have
0 ∈ LXdac+L∂ V3 (P, ξ1 , ξ2 ) if and only if (P, ξ1 , ξ2 ) ∈ F aug .
Since, F aug belongs to a level set of V3 , we conclude
that Proposition A.1(i) holds too.
To be able to apply Proposition A.1 and conclude the
proof, it remains to show that the trajectories starting from D(Rn × Rn × H0 ) are bounded. We reason
by contradiction, i.e., assume there exists a trajectory t 7→ (P (t), ξ1 (t), ξ2 (t)), with initial condition
(P (0), ξ1 (0), ξ2 (0)) ∈ D(Rn × Rn × H0 ) of Xdac+L∂
such that k(P (t), ξ1 (t), ξ2 (t)k → ∞. Since V3 is radially
unbounded, this implies V3 (P (t), ξ1 (t), ξ2 (t)) → ∞. Additionally, from Lemma 5.1, we know that 1>
n P (t) → Pl
and 1>
n ξ1 (t) → 0. Thus, there exists a sequence of times
{tk }∞
k=1 with tk → ∞ such that for all k ∈ Z≥1 ,
1>
n ξ1 (tk ) < 1/k,
max LXdac+L∂ V3 (P (tk ), ξ1 (tk ), ξ2 (tk )) > 0.

ν1 >
1 ζk
nk n
βν1 ν2
+
λ2 (L + L> ) > 0.
2nk 2

λmax (M )k(ξ1 (tk ), ξ2 (tk ))k2 +

Now, from the expression of ∂f  , since {P (tk )} is
bounded, the sequence {ζk } must be bounded. Combining these facts with λmax (M ) < 0, one can find k̄ ∈ Z≥1
such that the above inequality is violated for all k ≥ k̄,
which is a contradiction. For case (b), we use the bound
γk> M γk ≤ λmax (M )kηk k2 to deduce from (22) that

(21a)
(21b)

λmax (M )kηk k2 +

βν1 ν2
ν1 >
1n ζk +
λ2 (L + L> ) > 0.
nk
2nk 2

One can then use a similar argument as laid out in the
proof of Theorem 4.1, considering the two cases of 1>
n ζk
being bounded or unbounded, arriving in both cases at
similar contradictions. This concludes the proof. 2

Note that (21b) implies that there exists a sequence

{ζk }∞
k=1 with ζk ∈ ∂f (P (tk )) such that
− ζk> Lζk + ν1 ζk> ξ1 (tk ) − βν1 ν2 ξ1 (tk )> Lξ1 (tk )

− αkξ1 (tk )k2 − ν2 ξ2 (tk )> Lζk > 0,

for all k ∈ Z≥1 , where we have used the fact that an
element of LXdac+L∂ V3 (P, ξ1 , ξ2 ) has the form given in (20).
Letting ηk = ζk − n1 (1>
n ζk )1n , we use (1) to deduce from
the above inequality that

Note that as a consequence of the above result, the
dac+L∂ dynamics does not require any specific preprocessing for the initialization of the power allocations.
Each generator can select any generation level, independent of the other units, and the algorithm guarantees
1
1
>
convergence to the solutions of the ED problem.
− λ2 (L + L> )kηk k2 + ν1 ηk> ξ1 (tk ) + ν1 (1>
ζ
)(1
ξ
(t
))
k
1
k
n
n
2
n
Remark 5.3 (Distributed selection of algorithm design
1
1
2
− βν1 ν2 λ2 (L + L> )kξ1 (tk ) − (1>
parameters): The convergence of the dac+L∂ dynamn ξ1 (tk ))1n k
2
n
ics relies on a selection of the parameters α, β, ν1 and
2
>
− αkξ1 (tk )k − ν2 ξ2 (tk ) Lηk > 0.
ν2 ∈ R>0 that satisfy (17). Checking this inequality requires knowledge of the spectrum of matrices related
Further, using the expression
to the Laplacian matrix, and hence the entire network
structure. Here, we provide an alternative condition that
1 >
1
2
2
2
implies (17) and can be checked by the units in a disξ
(t
))1
k
=
kξ
(t
)k
−
(1
ξ
(t
))
,
kξ1 (tk ) − (1>
1 k
n
1 k
1 k
n n
n n
tributed way. Let nmax be an upper bound on the number
of units, dout
max be an upper bound on the out-degree of
the inequality can be rewritten as
all units, and amin be a lower bound on the edge weights,
γk> M γk +

1
>
ν1 (1>
n ζk )(1n ξ1 (tk ))
n
βν1 ν2
2
+
λ2 (L + L> )(1>
n ξ1 (tk )) > 0,
2n

n ≤ nmax , max dout (i) ≤ dout
max , min aij ≥ amin . (23)
i∈V

A straightforward generalization of [Mohar, 1991, Theorem 4.2] for weighted graphs gives rise to the following
lower bound on λ2 (L + L> ),

where γk> = [ηk> , ξ1 (tk )> , ξ2 (tk )> ]. Using now the
bound (21a), we arrive at the inequality,
ν1 >
βν1 ν2
γk> M γk +
1n ζk +
λ2 (L + L> ) > 0.
nk
2nk 2

(i,j)∈E

4amin
≤ λ2 (L + L> ).
n2max

(22)

(24)

On the other hand, using properties of matrix norms [Bern-

8

5.2

stein, 2005, Chapter 9], one can deduce
√
λmax (L> L) = kLk2 ≤ ( nkLk∞ )2
√
2
out 2
≤ (2 ndout
max ) ≤ 4nmax (dmax ) .

In this section, we study the robustness properties of
the dac+L∂ dynamics in the presence of time-varying
load signals and intermittent power unit generation. Our
analysis relies on the exponential stability of the mismatch dynamics between total generation and load established in Lemma 5.1, which implies that (16) is inputto-state stable (ISS) [Khalil, 2002, Lemma 4.6], and consequently robust against arbitrary bounded perturbations. The following result provides an explicit, exponentially decaying, bound for the evolution of any trajectory
of (16). While the rate of decay can also be determined
by computing the eigenvalues of matrix defining the dynamics, here we employ a Lyapunov argument to obtain
also the value of the gain associated to the rate.
Lemma 5.5 (Convergence rate of the mismatch dynamics (16)): Let R ∈ R2×2 be defined by

(25)

Using (24)-(25), the left-hand side of (17) can be upper
bounded by
ν1
ν22 λmax (L> L)
+
βν2 λ2 (L + L> )
2α
2
ν1 n2max
2ν 2 nmax (dout
max )
≤
+ 2
.
4amin βν2
α
Further, the right-hand side of (17) can be lower
bounded using (24). Putting the two together, we obtain
the new condition
2
4amin
2ν 2 nmax (dout
ν1 n2max
max )
< 2 ,
+ 2
4amin βν2
α
nmax

Robustness analysis

"
#
α2 + ν1 ν2 + (ν1 ν2 )2
α
1
.
R=
2αν1 ν2
α
1 + ν1 ν2

(26)

which implies (17). The network can ensure that this
condition is met in various ways. For instance, if the
bounds nmax , dout
max , and amin are not available, the network can implement distributed algorithms for max- and
min-consensus [Ren and Beard, 2008] to compute them
in finite time. Once known, any generator can select α,
β, ν1 and ν2 satisfying (26) and broadcast its choice. Alternatively, the computation of the design parameters
can be implemented concurrently with the determination of the bounds via consensus by specifying a specific
formula to select them that is guaranteed to satisfy (26).
Note that the units necessarily need to agree on the parameters, otherwise if each unit selects a different set of
parameters, the dynamic average consensus would not
track the average input signal.
•
Remark 5.4 (Distributed loads and transmission
losses): Here we expand on our observations in Remark 3.1 regarding the inclusion of additional constraints on the ED problem. Our algorithmic solution
can be easily modified to deal with the alternative scenarios studied in [Zhang et al., 2011, Kar and Hug, 2012,
Binetti et al., 2014a, Loia and Vaccaro, 2013], where
each generator has the knowledge of the load at the corresponding bus that it is connected to and the total load
is the aggregate of these individual loads. Mathematically, denoting the load demanded at generator
Pn bus i by
PiL ∈ R, the total load is given by Pl = i=1 PiL . For
this case, replacing the vector Pl er by P L in the dac+L∂
dynamics (15b) gives an algorithm that solves the ED
problem for the load Pl . Our solution strategy can also
handle transmission losses as modeled in [Binetti et al.,
2014a], where it is assumed that each generator i can
estimate the power loss in the transmission lines adjacent to it. With those values available, the generator
could add them to the quantity PiL , which would make
the network find a power allocation that takes care of
the transmission losses.
•

Then R  0 and any trajectory t 7→ x(t) of the dynam−c2 t
ics
kx(0)k, where c1 =
p (16) satisfies kx(t)k ≤ c1 e
λmax (R)/λmin (R) and c2 = 1/2λmax (R).
PROOF. Let A ∈ R2×2 be the system matrix of (16).
Then, one can see that A> R + RA = −I, i.e., V4 (x) =
x> Rx is a Lyapunov function for (16). Note that
λmin (R)kxk2 ≤ V4 (x) ≤ λmax (R)kxk2 .

(27)

From the Lyapunov equation, we have LAx V4 (x) =
−kxk2 ≤ − λmax1 (R) V4 (x), which implies V4 (x(t)) ≤

e−1/λmax (R) V4 (x(0)) along any trajectory t 7→ x(t)
of (16). Again using (27), we get
kx(t)k2 ≤

λmax (R) −1/λmax (R)
e
kx(0)k2 ,
λmin (R)

which concludes the claim.

2

In the above result, it is interesting to note that the convergence rate is independent of the specific communication digraph (as long as it is weight-balanced). We use
next the exponentially decaying bound obtained above
to illustrate the extent to which the network can collectively track a dynamic load (which corresponds to a
time-varying perturbation in the mismatch dynamics)
and is robust to intermittent power generation (which
corresponds to perturbations in the state of the mismatch dynamics).
5.2.1 Tracking dynamic loads
Here we consider a time-varying total load given by a
twice continuously differentiable trajectory R≥0 3 t 7→
9

The trajectory invariance routine ensures that the
dynamics (16) is the appropriate description for the evolution of the load satisfaction mismatch. This, together
with the ISS property established in Lemma 5.5, implies
that the mismatch effect in power generation caused by
addition/deletion events vanishes exponentially fast. In
particular, if the number of addition/deletion events is
finite, then the set of generators converge to the solution
of the ED problem. We formalize this next.
Proposition 5.7 (Convergence of dac+L∂ dynamics
under intermittent power generation): Let nmax be the
maximum number of generators that can contribute to
the power generation at any time. Let Σnmax be the set of
digraphs that are strongly connected and weight-balanced
and whose vertex set is included in {1, . . . , nmax }. Let
σ : [0, ∞) → Σnmax be a piecewise constant, rightcontinuous switching signal described by the set of
switching times {t1 , t2 , . . . } ⊂ R≥0 , with tk ≤ tk+1 ,
each corresponding to either an addition or a deletion
σ
event. Denote by Xdac+L∂
the switching dac+L∂ dynamics corresponding to σ, defined by (15) with L replaced
by L(σ(t)) for all t ≥ 0, and assume agents execute the
trajectory invariance routine when they leave or
join the network. Then,
(i) at any time t ∈ {0} ∪ {t1 , t2 , . . . }, if the variables (P (t), z(t)) for the generators in σ(t) satisfy
>
1>
n P (t) − Pl ≤ M1 and 1n z(t) ≤ M2 for some
M1 , M2 > 0, then the magnitude of the mismatch
between generation and load becomes less than or
equal to ρ > 0 in time

Pl (t) and show how the total generation of the network
under the dac+L∂ dynamics tracks it. We assume the
signal is known to an arbitrary unit r ∈ {1, . . . , n}. In
this case, the dynamics (16) takes the following form
"

ẋ1

"

#
=

ẋ2

1

#" #
x1

−ν1 ν2 −α

x2

0

"
+

0
−αṖl − P̈l

#
.

Using Lemma 5.5, one can compute the following bound
on any trajectory of the above system
kx(t)k ≤ c1 e−c2 t kx(0)k +

c1
sup αṖl (s) + P̈l (s) .
c2 s∈[0,t]

In particular, for a signal with bounded Ṗl and P̈l , the
mismatch between generation and load, i.e., x1 (t) is
bounded. Also, the mismatch has an ultimate bound
as t → ∞. The following result summarizes this notion
formally. The proof is straightforward application of
Lemma 5.5 following the exposition of input-to-state
stability in [Khalil, 2002].
Proposition 5.6 (Power mismatch is ultimately
bounded for dynamic load under dac+L∂ dynamics):
Let R≥0 3 t 7→ Pl (t) be twice continuously differentiable
and such that
sup Ṗl (t) ≤ d1 ,
t≥0

sup P̈l (t) ≤ d2 ,
t≥0

for some d1 , d2 > 0. Then, the mismatch 1>
n P (t) − Pl (t)
between load and generation is bounded along the trajectories of (15) and has ultimate bound cc12 (αd1 + d2 ), with
c1 , c2 given in Lemma 5.5. Moreover, if Ṗl (t) → 0 and
P̈l (t) → 0 as t → ∞, then 1>
n P (t) → Pl (t) as t → ∞.

tρ =

1  c1 (M1 + ν1 M2 ) 
ln
,
c2
ρ

provided no event occurs in the interval (t, t + tρ );
(ii) if the number of events is finite, say N , then the
σ
trajectories of Xdac+L∂
converge to the set of solutions of the ED problem for the group of generators
in σ(tN ) provided (17) is met for σ(tN ).
Note that the generators can ensure that the condition (17), required for the convergence of the dac+L∂ dynamics, holds at all times even under addition and deletion events, if they rely on verifying that (26) holds and
the bounds (23) are valid for all the topologies in Σnmax .

5.2.2 Robustness to intermittent power generation
Here, we characterize the algorithm robustness against
unit addition and deletion to capture scenarios with
intermittent power generation. Addition and deletion
events are modeled via a time-varying communication
digraph, which we assume remains strongly connected
and weight-balanced at all times. When a unit stops generating power (deletion event), the corresponding vertex and its adjacent edges are removed. When a unit
starts providing power (addition event), the corresponding node is added to the digraph along with a set of
edges. Given the intricacies of the convergence analysis
for the dac+L∂ dynamics, cf. Theorem 5.2, it is important to make sure that the state v remains in the set H0 ,
irrespectively of the discontinuities caused by the events.
The following routine makes sure that this is the case.
trajectory invariance: When a unit i joins the
network at time t, it starts with vi (t) = 0. When a
unit i leaves the network at time t, it passes a token
with value vi (t) to one of its in-neighbors j ∈ N in (i),
who resets its value to vj (t) + vi (t).

6

Simulations in a IEEE 118 bus system

This section illustrates the convergence of the dac+L∂
dynamics to the solutions of the ED problem (4) starting from any initial power allocation and its robustness properties. We consider the IEEE 118 bus system [IEEE 118 bus], that consists of 54 generators.
The cost function of each generator i is quadratic,
fi (Pi ) = ai + bi Pi + ci Pi2 , with coefficients belonging to
the ranges ai ∈ [6.78, 74.33], bi ∈ [8.3391, 37.6968], and
ci ∈ [0.0024, 0.0697]. The communication topology is the
digraph G described in Table 1. We choose the design parameters as ν1 = 1, ν2 = 1.3, α = 10, β = 40,  = 0.0086,
which satisfy the conditions (6) and (17) for G. The to10

G

digraph over 54 vertices consisting of a directed cycle through vertices 1, . . . , 54 and bi-directional edges
{(i, id54 (i + 5)), (i, id54 (i + 10)), (i, id54 (i + 15)), (i, id54 (i + 20))} for each i ∈ {1, . . . , 54}, where
id54 (x) = x if x ∈ {1, . . . , 54} and x − 54 otherwise. All edge weights are 0.1.

Ĝ

obtained from G by replacing the directed cycle with an undirected one keeping the edge weights same

Ĝ\{4,11,25,45}

obtained from Ĝ by removing the vertices {4, 11, 25, 45} and the edges adjacent to them

Ĝ\{4,25,27}
obtained from Ĝ by removing the vertices {4, 25, 27} and the edges adjacent to them
Table 1
Definition of the digraphs G, Ĝ, Ĝ\{4,11,25,45} , and Ĝ\{4,25,27} .
4700

tal load is 4600 for the first 150 seconds and 4200 for the
next 150 seconds, and is known to unit 3. Figure 1(a)-(c)
depicts the evolution of the power allocation, total cost,
and the mismatch between the total generation and
load under the dac+L∂ dynamics starting at the initial
condition (P (0), z(0), v(0)) = (0.5 ∗ (P m + P M ), 0, 0).
Note that the generators initially converge to a power
allocation that meets the load 4600 and minimizes the
total cost of generation. Later, with the decrease in
desired load to 4200, the network decreases the total
generation while minimizing the total cost.
Next, we consider a time-varying total load given by a
constant plus a sinusoid, Pl (t) = 4300 + 100 sin(0.05t).
With the same communication topology, design parameters, and initial condition as above, Figure 1 (d)-(f) illustrates the behavior of the network under the dac+L∂
dynamics. As established in Proposition 5.6, the total
generation tracks the time-varying load signal and the
mismatch between these values has an ultimate bound.
Additionally, to illustrate how that the mismatch vanishes if the load becomes constant, we show in Figure 2
a load signal that consists of short bursts of sinusoidal
variation that decay exponentially. The difference between generation and load becomes smaller and smaller
as the load tends towards a constant signal.
Our final scenario considers addition and deletion of generators. The initial communication topology is the undirected graph Ĝ described in Table 1. The design parameters and the initial condition are the same as above. The
total load is 4200 and is same at all times. For the first 100
seconds, the power allocations converge to a neighborhood of a solution of the ED problem for the set of generators in Ĝ. At time t = 100s, the units {4, 11, 25, 45} stop
generating power and leave the network. We select these
generators because of their substantial impact in the total power generation. After this event, the resulting communication graph is Ĝ\{4,11,25,45} , cf. Table 1. The generators implement the trajectory invariance routine,
after which the dac+L∂ dynamics drives the mismatch to
zero and minimizes the total cost. At t2 = 200s, another
event occurs, the units {11, 45} get added to the network
while the generator 27 leaves. The resulting communication topology is Ĝ\{4,25,27} , cf. Table 1. After executing the trajectory invariance routine, the dynamics
converges eventually to the optimizers of the ED problem for the set of generators in Ĝ\{4,25,27} , as shown in
Figure 1(g)-(i). This example illustrates the robustness

total generation
total load

4500
4300
4100
3900
0

150

300

450

600

750

900

Fig. 2. Evolution of the total power generation for the IEEE
118 bus example under the dac+L∂ dynamics for the communication digraph G, design parameters ν1 = 1, ν2 = 1.3,
α = 10, β = 40 and  = 0.0086, and time-varying total load.
The example depicts the input-to-state stability of the mismatch dynamics.

of the dac+L∂ dynamics against intermittent generation
by the units, as formally established in Proposition 5.7.
In addition to the presented examples, we also successfully simulated scenarios of the kind described in Remark 5.4, where the total load is not known to a single
generator and is instead the aggregate of the local loads
connected to each of the generator buses, but we do not
report them here for space reasons.
7

Conclusions

We have designed a novel provably-correct distributed
strategy that allows a group of generators to solve the
economic dispatch problem starting from any initial
power allocation. Our algorithm design combines elements from average consensus to dynamically estimate
the mismatch between generation and desired load and
ideas from distributed optimization to dynamically allocate the unit generation levels. Our analysis has shown
that the mismatch dynamics between total generation
and load is input-to-state stable and, as a consequence,
the coordination algorithm is robust to initialization errors, time-varying load signals, and intermittent power
generation. Our technical approach relies on tools from
algebraic graph theory, dynamic average consensus,
set-valued dynamical systems, and nonsmooth analysis,
including a novel refinement of the LaSalle Invariance
Principle for differential inclusions that we have stated
and proved. Future work will explore the study of the
preservation of the generator box constraints under the
proposed coordination strategy, the extension to scenarios that involve additional constraints, such as transmission losses, transmission line capacity constraints,
ramp rate limits, prohibited operating zones, and valvepoint loading effects, and the study of the stability and

11

4

500

8

400

x 10

4800

7.5

4600

7

4400

6.5

4200

total generation
total load

300
200
100
0
0

50

100

150

200

250

300

6
0

50

(a) Power allocation

100

150

200

250

300

(b) Total cost

50

100

150

200

250

300

(c) Total mismatch

4

500

8

400

x 10

4600

7.5

4400

300

7

200

4200

100

6.5

0
0

6
0

50

100

150

200

250

300

50

(d) Power allocation

100

150

200

250

300

4000
0

total generation
total load
50

(e) Total cost

100

150

200

250

300

(f) Total mismatch

4

500

8

400

7.5

x 10

7
300

6.5
6

200

5.5

100
0
0

4000
0

5
50

100

150

200

(g) Power allocation

250

300

4.5
0

50

100

150

200

(h) Total cost

250

300

4800
4600
4400
4200
4000
3800
3600
3400
3200
0

total generation
total load

50

100

150

200

250

300

(i) Total mismatch

Fig. 1. Evolution of the power allocation, the total cost, and the total mismatch between generation and load under the dac+L∂
dynamics for the IEEE 118 bus example in different scenarios. In the first case (a)-(c), the communication topology is G, the
load is initially 4600 and later 4200, and parameters are ν1 = 1, ν2 = 1.3, α = 10, β = 40, and  = 0.0086. In the second scenario
(d)-(f), the digraph and the parameters remain the same but the load is time-varying given by Pl (t) = 4300 + 100 sin(0.05t).
In the last case (g)-(i), the parameters remain the same, the communication graph is initially the graph Ĝ. At t = 100s,
units {4, 11, 25, 45} leave the network, resulting in the communication topology Ĝ\{4,11,25,45} , and the remaining agents run
the trajectory invariance routine. Later, at t = 200s, units {11, 45} join the network while unit 27 leaves it, resulting
in the communication topology Ĝ\{4,25,27} . After implementing the trajectory invariance routine, the dac+L∂ dynamics
eventually converges to an optimizer of the ED problem for the network Ĝ\{4,25,27} .

convergence properties of algorithm designs that combine our approach here with traditional primary and
secondary generator controllers.

A distributed auction-based algorithm for the nonconvex economic dispatch problem. IEEE Transaction on Industrial Informatics, 10(2):1124–1132, 2014b.
S. Boyd and L. Vandenberghe. Convex Optimization. Cambridge
University Press, 2009. ISBN 0521833787.
F. Bullo, J. Cortés, and S. Martı́nez. Distributed Control of
Robotic Networks. Applied Mathematics Series. Princeton University Press, 2009. ISBN 978-0-691-14195-4. Electronically
available at http://coordinationbook.info.
A. Cherukuri and J. Cortés. Distributed generator coordination
for initialization and anytime optimization in economic dispatch. IEEE Transactions on Control of Network Systems,
2013. Conditionally accepted. Available at http://carmenere.
ucsd.edu/jorge.
B. H. Chowdhury and S. Rahman. A review of recent advances
in economic dispatch. IEEE Transactions on Power Systems,
5(4):1248–1259, November 1990.
J. Cortés. Discontinuous dynamical systems - a tutorial on solutions, nonsmooth analysis, and stability. IEEE Control Systems Magazine, 28(3):36–73, 2008.
A. D. Dominguez-Garcia, S. T. Cady, and C. N. Hadjicostis.
Decentralized optimal dispatch of distributed energy resources.
In IEEE Conf. on Decision and Control, pages 3688–3693,

References
A. Arsie and C. Ebenbauer. Locating omega-limit sets using
height functions. Journal of Differential Equations, 248(10):
2458–2469, 2010.
J. P. Aubin and A. Cellina. Differential Inclusions, volume 264 of
Grundlehren der mathematischen Wissenschaften. Springer,
New York, 1984.
A. Bacciotti and F. Ceragioli. Stability and stabilization of discontinuous systems and nonsmooth Lyapunov functions. ESAIM:
Control, Optimisation & Calculus of Variations, 4:361–376,
1999.
D. S. Bernstein. Matrix Mathematics. Princeton University Press,
Princeton, NJ, 2005.
G. Binetti, A. Davoudi, F. L. Lewis, D. Naso, and B. Turchiano.
Distributed consensus-based economic dispatch with transmission losses. IEEE Transactions on Power Systems, 29(4):1711–
1720, 2014a.
G. Binetti, A. Davoudi, D. Naso, B. Turchiano, and F. L. Lewis.

12

and Ebenbauer, 2010] for differential equations. Our
motivation for developing this refinement comes from
the need to provide the necessary tools to tackle the
convergence analysis of the coordination algorithms
presented in Sections 4 and 5. Nevertheless, the results
stated here are of independent interest.
Proposition A.1 (Refined LaSalle Invariance Principle for differential inclusions): Let F : Rn ⇒ Rn be upper
semicontinuous, taking nonempty, convex, and compact
values at every point x ∈ Rn . Consider the differential
inclusion ẋ ∈ F (x) and let t 7→ ϕ(t) be a bounded solution whose omega-limit set Ω(ϕ) is contained in S ⊂ Rn ,
a closed embedded submanifold of Rn . Let O be an open
neighborhood of S where a locally Lipschitz, regular function W : O → R is defined. Assume the following holds,
(i) the set E = {x ∈ S | 0 ∈ LF W (x)} belongs to a
level set of W ,
(ii) for any compact set M ⊂ S with M ∩ E = ∅, there
exists a compact neighborhood Mc of M in Rn and
δ < 0 such that supx∈Mc max LF W (x) ≤ δ.
Then, Ω(ϕ) ⊂ E.

Hawaii, USA, December 2012.
L. Du, S. Grijalva, and R. G. Harley. Potential-game theoretical
formulation of optimal power flow problems. In IEEE Power
and Energy Society General Meeting, San Diego, CA, July
2012. Electronic Proceedings.
H. Farhangi. The path of the smart grid. IEEE Power and Energy
Magazine, 8(1):18–28, 2010.
R. A. Freeman, P. Yang, and K. M. Lynch. Stability and convergence properties of dynamic average consensus estimators.
In IEEE Conf. on Decision and Control, pages 398–403, San
Diego, CA, 2006.
IEEE 118 bus. http://motor.ece.iit.edu/data/JEAS IEEE118.doc.
B. Johansson and M. Johansson. Distributed non-smooth resource
allocation over a network. In IEEE Conf. on Decision and
Control, pages 1678–1683, Shanghai, China, December 2009.
S. Kar and G. Hug. Distributed robust economic dispatch in
power systems: A consensus + innovations approach. In IEEE
Power and Energy Society General Meeting, San Diego, CA,
July 2012. Electronic proceedings.
H. K. Khalil. Nonlinear Systems. Prentice Hall, 3 edition, 2002.
ISBN 0130673897.
S. S. Kia, J. Cortés, and S. Martı́nez. Dynamic average consensus under limited control authority and privacy requirements.
International Journal on Robust and Nonlinear Control, 2014.
To appear.
V. Loia and A. Vaccaro. Decentralized economic dispatch in smart
grids by self-organizing dynamic agents. IEEE Transactions
on Systems, Man & Cybernetics: Systems, 2013. To appear.
B. Mohar. Eigenvalues, diameter, and mean distance in graphs.
Graphs and Combinatorics, 7(1):53–64, 1991.
R. Mudumbai, S. Dasgupta, and B. B. Cho. Distributed control
for optimal economic dispatch of a network of heterogeneous
power generators. IEEE Transactions on Power Systems, 27
(4):1750–1760, 2012.
A. Pantoja, N. Quijano, and K. M. Passino. Dispatch of distributed generators under local-information constraints. In
American Control Conference, pages 2682–2687, Portland, OR,
June 2014.
W. Ren and R. W. Beard. Distributed Consensus in MultiVehicle Cooperative Control. Communications and Control
Engineering. Springer, 2008. ISBN 978-1-84800-014-8.
A. Simonetto, T. Keviczky, and M. Johansson. A regularized
saddle-point algorithm for networked optimization with resource allocation constraints. In IEEE Conf. on Decision and
Control, pages 7476–7481, Hawaii, USA, December 2012.
L. Xiao and S. Boyd. Optimal scaling of a gradient method
for distributed resource allocation. Journal of Optimization
Theory & Applications, 129(3):469–488, 2006.
W. Zhang, W. Liu, X. Wang, L. Liu, and F. Ferrese. Online
optimal generation control based on constrained distributed
gradient algorithm. IEEE Transactions on Power Systems,
2014. To appear.
Z. Zhang and M. Chow. Convergence analysis of the incremental
cost consensus algorithm under different communication network topologies. IEEE Transactions on Power Systems, 27(4):
1761–1768, 2012.
Z. Zhang, X. Ying, and M. Chow. Decentralizing the economic
dispatch problem using a two-level incremental cost consensus
algorithm in a smart grid environment. In North American
Power Symposium, Boston, MA, August 2011. Electronic Proceedings.

A

Before proceeding with the proof of the result, we establish an auxiliary result.
Lemma A.2 Under the hypotheses of Proposition A.1,
the sets Ω(ϕ) and E have nonempty intersection.

PROOF. By contradiction, assume Ω(ϕ) ∩ E = ∅.
Then, using the hypothesis (ii) in Proposition A.1, there
exists δ < 0 such that supx∈Ω(ϕ) max LF W (x) ≤ δ. Let
x ∈ Ω(ϕ). Since this set is weakly positively invariant,
there exists a trajectory t 7→ ϕ̃(t) of the differential
inclusion with ϕ̃(0) = x such that ϕ̃(t) ∈ Ω(ϕ). Since
d
dt W (ϕ̃(t)) ∈ LF W (ϕ̃(t)) for almost all t ≥ 0, we get
W (ϕ̃(t)) − W (x) ≤ δt. This is in contradiction with the
fact that t 7→ ϕ̃(t) belongs to the compact set Ω(ϕ),
where W is lower bounded. 2
We are now ready to prove Proposition A.1.
PROOF. [Proof of Proposition A.1] We consider two
cases, depending on whether the set Ω(ϕ) (a) is or (b)
is not contained in a level set of W . In case (a), given
any x ∈ Ω(ϕ), there exists a trajectory of F starting at
x that remains in Ω(ϕ) (because of the weak positive
invariance of the omega-limit set). If x 6∈ E, then by
the hypotheses (ii), there exists a compact neighborhood
Mx of x in Rn and δ < 0 such that supy∈Mx LF W (y) ≤
δ. Since Ω(ϕ) ⊂ S, the trajectory of F starting at x
remains in the set Mx ∩ S for a finite time, say t1 . Over
the time interval [0, t1 ], we have W (t) − W (0) ≤ δt.
This, however, is in contradiction with the fact that the
trajectory belongs to Ω(ϕ) which is contained in a level
set of W . Therefore, x ∈ E, and since this point is generic,
we conclude Ω(ϕ) ⊂ E.
Next, we consider case (b) and reason by contradiction, i.e., assume that Ω(ϕ) is not contained in E (see

Refined LaSalle Invariance Principle for differential inclusions

In this section we provide a refinement of the LaSalle
Invariance Principle for differential inclusions, see
e.g., [Cortés, 2008], by extending the results of [Arsie

13

B

W
Ω(ϕ)

E

bm

bPM

B

E

P

S

Ω(ϕ)

P
U

P

S

Fig. A.1. Illustration (adapted from [Arsie and Ebenbauer,
2010, Figure 1]) depicting various elements involved in the
case (b) of the proof of Proposition A.1.

x∈U

set-valued

sup max LF W (x) > δ.

x∈Mc

Note that this implies that supx∈Mc max LF W (x) ≥ 0.
Now, for each k ∈ Z≥1 , consider the compact neighborhood Mk = M + B(0, k1 ) of M. From the above, we deduce the existence of a sequence {xk }∞
k=1 with xk ∈ Mk
such that

bP
M = sup W (x).

lim max LF W (xk ) = ` ≥ 0.

x∈UP

k→∞

Note that the neighborhoods U and U P can be chosen
such that the set Ω(ϕ) \ (U ∪ U P ) is nonempty, compact,
and its intersection with E is empty. Along with this,
one can select  in such a way that bm > bP
M and from
assumption (ii) we get
sup

max LF W (x) ≤ δ < 0,

(A.1)

Since the whole sequence belongs to the compact set
M1 , there exists a subsequence, which we denote with
the same indices for simplicity, such that
lim xk = x̃ ∈ M.

k→∞
x∈B \(U ∪U P )

Lie

PROOF. We reason by contradiction, i.e., assume that
for all compact neighborhoods Mc of M in Rn and all
δ < 0, we have

Figure A.1). Given  > 0, let B ⊂ O be a compact
neighborhood of Ω(ϕ) in Rn such that d(B , Ω(ϕ)) ≤ .
Let U be an open neighborhood of E in Rn and define
U = U ∩ B . Note that U is nonempty because Ω(ϕ) ∩ E
is nonempty by Lemma A.2. Since Ω(ϕ) is not contained
in a level set of W but E is by hypotheses (i), we can
choose P ∈ Ω(ϕ) \ E such that W (P ) 6= W (E). Without
loss of generality, assume W (P ) < W (E) (the reasoning
is analogous for the other case). Select an open neighborhood U P of P in Rn and define UP = U P ∩B . Define
the following quantities
bm = inf W (x),

of

Here we present an auxiliary result for the convergence
analysis of the algorithms of Sections 4 and 5.
Lemma A.1 (Continuity property of set-valued Lie
derivatives): Let W : Rn → R be a locally Lipschitz and
regular function. Let g : Rn × Rn → Rn be a continuous
function and define the set-valued map F : Rn ⇒ Rn by
F (x) = {g(x, ζ) | ζ ∈ ∂W (x)}. Assume that
(i) S is an embedded submanifold of Rn such that
ζ > g(x, ζ) ≤ 0 for all x ∈ S and all ζ ∈ ∂W (x),
(ii) for any x ∈ S, if ζ > g(x, ζ) = 0 for some ζ ∈ ∂W (x),
then x ∈ E = {z ∈ S | 0 ∈ LF W (z)}.
Then, for any compact set M ⊂ S with M ∩ E = ∅, there
exists a compact neighborhood Mc of M in Rn and δ < 0
such that supx∈Mc max LF W (x) ≤ δ.

U

U
UP

Continuity properties
derivatives

(A.1)

(A.2)

From (A.1), there exists a sequence ζk ∈ ∂W (xk ) such
that

(in the case W (P ) > W (E), we would reason with the
quantities bM = supx∈U W (x) and bP
m = inf x∈UP W (x)).
Since Ω(ϕ) is the omega-limit set of ϕ and B is a compact neighborhood of Ω(ϕ), there exists t1 > 0 such
that ϕ(t1 ) ∈ UP and ϕ(t) ∈ B for all t ≥ t1 . Moreover,
since Ω(ϕ) ∩ E is nonempty, there must also exist t2 > t1
such that ϕ(t2 ) ∈ U . From continuity of the trajectory
we deduce that there exist times t∗1 , t∗2 ∈ (t1 , t2 ), t∗1 < t∗2
such that ϕ(t∗1 ) and ϕ(t∗2 ) lie on the boundary of the
compact set B \ (U ∪ UP ), with ϕ(t∗1 ) belonging to the
closure of UP and ϕ(t∗2 ) to the closure of U . However,
∗
this is not possible as W (ϕ(t∗2 )) ≥ bm > bP
M ≥ W (ϕ(t1 ))
∗ ∗
and, in the interval [t1 , t2 ], the trajectory belongs to
B \ (U ∪ UP ), where the function W can only decrease
due to (A.1), which is a contradiction. 2

lim ζk> g(xk , ζk ) ≥ 0.

k→∞

(A.3)

Since ∂W is upper semicontinuous with compact values, the set ∂W (M1 ) is compact, cf. [Aubin and Cellina, 1984, Proposition 3, p. 42]. This implies that the
sequence {ζk } belongs to the compact set ∂W (M1 ) and
so, there exists a subsequence, denoted again by the same
indices for simplicity, such that ζk → ζ̃. Since ∂W is upper semicontinuous and takes closed values, we deduce
from [Aubin and Cellina, 1984, Proposition 2, p. 41] that
ζ̃ ∈ ∂W (x̃). From (A.2) and (A.3), since g is continuous,
we obtain ζ̃ > g(x̃, ζ̃) ≥ 0. By assumption (i), this implies ζ̃ > g(x̃, ζ̃) = 0. Assumption (ii) then implies x̃ ∈ E,
which along with (A.2) contradicts M ∩ E = ∅. 2
14

