Optimization of stochastic lossy transport networks and
applications to power grids
arXiv:1712.07411v1 [math.OC] 20 Dec 2017

Alessandro Zocca∗ and Bert Zwart†
December 21, 2017

Abstract
Motivated by developments in renewable energy and smart grids, we formulate a stylized
mathematical model of a transport network with stochastic load fluctuations. Using an
affine control rule, we explore the trade-off between the number of controllable resources
in a lossy transport network and the performance gain they yield in terms of expected
transportation losses. Our results are explicit and reveal the interaction between the level
of flexibility, the intrinsic load uncertainty and the network structure.
Keywords: stochastic network; optimization; lossy transport; total effective resistance; efficiency;
robustness; power grids; renewable energy sources.

1

Introduction

A transport network is an abstract model describing a structure in which some commodity is
transferred from the “source” nodes of the network to the “sink” nodes according to a specified
routing that is determined by some external principle or design, see [10, 46]. Examples of transport
networks are road networks, railways, pipes, and power grids.
In this work we focus in particular on lossy transport networks where a fraction of the
transported good is inevitably lost, having in mind as primary application power systems in
which part of the transported electricity is lost due to heat dissipation in the line conductors.
The main question that we want to address in the present paper is whether these transportation
losses can be reduced in the scenario in which we have no direct control on the routing, but some
of the nodes of the network have controllable loads. This is the case for power systems in which
the line flows are determined by physical laws, but at the same time feature an increasing number
of controllable energy resources, like energy storage devices, smart buildings and appliances, and
electric vehicles.
In the present work we consider a probabilistic model to describe the stochastic fluctuations of
the load in (a subset of) the nodes of the network. This is instrumental to model the stochastic load
fluctuations due to power demand uncertainty and intermittent generation by renewable energy
∗

Department of Computing and Mathematical Sciences, California Institute of Technology, Pasadena, CA
91125, Email: azocca@caltech.edu
†
Centrum Wiskunde and Informatica (CWI), 1098 XG Amsterdam, The Netherlands

1

sources. The current power grids were originally built around conventional power generation
systems and therefore they are not equipped to cope with this massive amount of uncertainty,
especially in power supply. Managing this uncertainty on such a large scale with existing methods
will soon become crucial: in the next decades power grids will have to become more flexible and
robust to reduce the likelihood of contingencies and blackouts, whose social and economic impact
is enormous.
As mentioned earlier, next to the increasing renewable penetration, there is another powerful
trend that is driving this pervasive evolution of power system: the advent of distributed energy
resources. At a high level, all these resources can be seen as “virtual storage/batteries”, in the
sense that over short time-scales they can dynamically reduce their power consumptions and even
inject electricity in the power grid when necessary, see [32, 43]. Even if at the present stage these
resources are not fully incorporated, they have a huge potential: if their penetration increases
and we can actively and optimally control them, then they can make power grids more flexible
and at the same time effectively mitigate the volatile nature of renewable power generation and
allow a higher share of renewable energy sources.
The controllable loads that we consider in this work should be seen as an abstraction of more
concrete examples, such as (i) actual energy storage that is neither full nor empty, (ii) distribution
grids with ample flexible and/or deferrable load [32, 33, 36, 37] or (iii) conventional generators
that can provide balancing services.
The stochastic network model considered in this paper aims to understand the potential that
these controllable resources can have in mitigating the load uncertainty and in particular the
transportation losses due to the stochastic load fluctuations in the network nodes. Specifically, we
consider a network with random sources and sinks modelled by an undirected connected weighted
graph G consisting of n nodes and m edges to which are associated non-negative weights w ∈ Rm
+,
and investigate how much could the average total loss be reduced by operating optimally the
subset B ⊆ V of nodes with controllable loads.
The metric we consider here to quantify the transportation losses due to stochastic fluctuations
is a quadratic function of the load profile that has been introduced in [21], where the authors
showed that is a good approximation for the true total power losses in AC power grids. It is a
scalar quantity that in our setting captures both the correlations in the load fluctuations as well
as the role of the controllable loads in that are available in the network.
Such a metric generalizes the notion of total effective resistance Rtot (G) of the graph G,
also known as Kirchhoff index. This is a key quantity that measures how well connected and
robust a network is and for this reason has been extensively studied and rediscovered in various
contexts, such as complex network analysis [13], theoretical chemistry (for an overview see [47]
and references therein), and probability theory [12].
We then formulate a constrained optimization problem to find the optimal load-sharing factors
for the controllable loads. We prove that the solution is unique and give a closed-form expression
for the optimal control in which appears an interplay between the network structure, the location
of the controllable nodes, and the correlation structure of the load fluctuations.
In particular, our result shows how the correlation structure between load fluctuations affects
the optimal operation of the network resources and this is extremely relevant in power grids with
geographically close wind or solar farms, whose power outputs will be highly correlated. These
insights are derived without making Gaussian assumptions on the distribution of the fluctuations.
We then use this explicit solution to explore the trade-off between the number of controllable
resources available in a network and the performance gain they yield, quantified as transportation
2

loss reduction. The analysis builds on and extends that of the conference paper [48], in which
the scenario with exactly two controllable loads has been considered. In particular, the authors
show therein that the expected total loss due to fluctuations can be reduced by 25% in a line
network by adding one controllable storage device. In the present paper we extended that insight
by showing that for large graphs the total expected loss can be on average reduced by a factor
(1 + 1/k)/2 when k controllable loads are available. This insight holds when each node in the
network is self-sustainable on average.
This suggests that, even if power grids are becoming locally more robust and self-supportive,
transportation losses can be reduced by up to 50 percent, and the number of controllable loads
or balancing services k quantify how close one can get to this reduction. Though our model is
stylized, it provides a simple quantitative estimate on the value of balancing services in a scenario
where each node in the network is self-sustainable on average.
We remark that our stylized model is “static”, in the sense the optimal control we derive does
rebalance the total power mismatch in the network in any scenario, but does not depend on the
realized load fluctuations. Indeed, it depends only on the network structure, on the location of
the controllable loads and on their average covariance structure of the load fluctuations. For this
reason our model is not specific for a precise time-scale and provides insights into the value of
balancing services both in real-time operations as well as long-term planning.
Our work complements a large body of literature on optimal topology design for power
grids [11, 17, 21] and on optimal control of multiple controllable devices and generators where
often the designed controls, often called participation factors are also affine in the stochastic
load fluctuations, see e.g., [1, 5, 6, 18, 19, 22, 30, 34, 39, 40, 41]. Optimal policies for storage
management, especially aiming at the mitigation of the uncertainties in wind generation, have been
explored in [2, 16, 15, 44], where, however, the physical network is not modelled explicitly. Optimal
storage placement can increase network reliability, as shown in [3] by simulation techniques.
Storage can also be used for arbitrage, exploiting temporal price differences [8, 9] and the impact
of storage on energy markets has been studied in [7, 14].
The rest of the paper is organized as follows. We provide a detailed model description in
Section 2. In Section 3 we investigate the optimal load sharing factors in several scenarios. These
results are applied in Section 4 where a scaling law is presented for large networks. In Section 5
we report several numerical experiments. Section 6 concludes.

2

Model description

In this paper we model a lossy transport network as a weighted graph (G, w), where the graph G
is a simple undirected graph G = (V, E) with |V | = n nodes and |E| = m edges and w ∈ Rm
+ is the
collection of edge weights. In the context of power grids, the nodes of G are often referred to as
buses, the edges as transmission lines and the quantity wi,j is the susceptance of the transmission
line connecting buses i and j. Missing edges can be thought as edges with zero weight.
The weighted Laplacian matrix of the graph G is the matrix L ∈ Rn×n defined as
(P
if i = j,
k6=i wi,k
Li,j :=
−wi,j
if i 6= j.
Denote by λ1 , λ2 , . . . , λn its eigenvalues and by v1 , v2 , . . . , vn the corresponding eigenvectors, which
we can take to be pairwise orthogonal. It is well-known that L is a real symmetric positive
3

semi-definite matrix, which means that it has a real non-negative spectrum, i.e., 0 ≤ λ1 ≤ . . . ≤ λn .
All the rows of L sum up to zero and thus the matrix L is singular. Under the assumption
that G is a connected graph, the eigenvalue λ1 = 0 has multiplicity one and the corresponding
eigenvector v1 is the vector with all unit entries, which we denote by 1 ∈ Rn .
Let L+ ∈ Rn×n be the Moore-Penrose pseudoinverse of the weighted Laplacian matrix L.
Using the eigenspace structure of L, the pseudoinverse L+ can be defined as
−1

11T
11T
−
.
L := L +
n
n
+

The matrix L+ is also real, symmetric, and positive semi-definite and its spectrum is 0 < λ−1
n ≤
−1
. . . ≤ λ2 . For further spectral properties of graphs, we refer the reader to [38, 45].
Denote by p ∈ Rn the load profile at the network nodes, where pi is the load at node i for
every i = 1, . . . , n. In the context of power grids the i-th entry of the vector p models the power
generated (if pi > 0) or consumed (if pi < 0) at node i. We say that a load profile p ∈ Rn is
balanced if 1T p = 0.
Given a network with a balanced load profile p ∈ Rn , we define its total loss H = H(p) as
1
H := pT L+ p.
2
The scalar quantity H is a quadratic form of the load profile vector p and, as such, is always
non-negative, thanks to the fact that L+ is a positive semi-definite matrix.
The total loss H will be central in our analysis. H is the most natural choice for an efficiency
metric in a lossy transport network as H quantifies the total losses occurred in the network
when subject to a given load profile. More specifically, for power systems H has been shown
by [21] to be good approximation for total power loss in AC power grids. The quantity H has
also been considered as measure of the network tension, i.e., the aggregate load of the network
in [20, 24], where it is shown to be monotone along a cascade failure, as long as the network
remains connected.
As we will show later, the total loss H is indeed a generalized total effective resistance in
the sense that it quantifies how robust the network G is against a stochastic load profile with a
predefined covariance structure.

2.1

Stochastic loads and load-sharing factors

In this work we are particularly interested in a transport network with a stochastic load profile,
which means that we will take p to be a multivariate random variable.
More precisely, we will take p to be of the form p = µ + ω, where µ ∈ Rn is the nominal
load profile in the network and ω is a n-dimensional random variable modeling the fluctuations.
We henceforth assume that Eω = 0 and denote by Σ ∈ Rn×n the covariance matrix of the load
fluctuations, namely Σi,j = cov(ωi , ωj ) = E[ωi ωj ]. Nodes in which there are no stochastic load
fluctuations can be modelled by setting all the entries equal to zero in the corresponding row
and column of the matrix Σ. We denote by S ⊆ V the subset of nodes with stochastic load
fluctuations and we will henceforth assume that |S| ≥ 1.
We further assume that the load profile is balanced on average, namely 1T µ = 0. This
assumption, however, does not guarantee that every realization of the load profile is balanced:
4

indeed, the total mismatch 1T p in the network is a random variable, which is distributed as the
sum of fluctuations, since
n
X
T
T
T
1 p = 1 (p − µ) = 1 ω =
ωi .
i=1

Let σ 2 := Var
as

P

n
i=1



ωi be the variance of the sum of the load fluctuations, which also rewrites
2

σ =

n
X

Σi,j = 1T Σ1 = tr(Σ11T ).

(1)

i,j=1

Since we assumed that there is at least one node with stochastic load fluctuations, we have σ 2 > 0.
In order to cope with the stochastic fluctuations and, in particular, to keep the network
balanced, we assume that load at each node is controllable: for every i = 1, . . . , n, node i can
deal with (either
Pn generate or store) a controllable fraction αi ∈ R of the realized total mismatch
T
1 (p − µ) = i=1 ωi . In other words, we assume that while using the load-sharing coefficients
α = (α1 , . . . , αn ) ∈ Rn the net load profile p(α) is given by


p(α) = p1 −α1

n
X
i=1

ωi , . . . , pn −αn

n
X





ωi = µ1 +ω1 −α1

i=1

n
X

ωi , . . . , µn +ωn −αn

i=1

n
X


ωi .

i=1

P
For any j = 1, . . . , n, the term αj ni=1 ωi corresponds to the power generated or stored in
the corresponding controllable load when using an affine control responsive to stochastic load
fluctuations. We can rewrite the net power injection p(α) in matrix form as
p(α) = Cα p = Cα (µ + ω) = µ + Cα ω,
where Cα ∈ Rn×n is the matrix defined as

1 − α1 −α1
...
−α1
 −α2 1 − α2 −α2 . . .
−α2

...
..
..

T
.
.
Cα := I − α 1 = 

.
..
..

.
−αn
. . . −αn 1 − αn

(2)





.



Note that the last equality in (2) follows from the fact that 1T µ = 0 since Cα µ = (I − α 1T )µ =
µ − 0 · α = µ. When the load-sharing coefficient α is used, the total mismatch in the network is
equal to
1T p(α) = 1T Cα (µ + ω) = 1T (I − α 1T )ω = (1T − (1T α)1T )ω.
Therefore, since we want the net load profile to be balanced for any realization ω, we need to
impose that
n
X
T
1 α=
αi = 1.
(3)
i=1

We remark that the affine control modelled via the load-sharing vector α is a simplification,
especially when the controllable loads models energy storage. Indeed, it does not incorporate
many details, among which possible ramping constraints or the current state of charge. In our
model such details are omitted on purpose to have a mathematically tractable optimization
5

problem and to better identify the interplay between load uncertainty and storage operations. A
further simplification we make is that we allow to choose α without invoking line limits, as we
also do this for mathematical tractability, see Section 6 for further comments. Note that affine
control is fairly common in power systems operations, see, e.g., [4].

2.2

Expected total loss: definition and properties

For any 1 ≤ k ≤ n, we can model a network where exactly k nodes have controllable loads by
imposing that the remaining n − k nodes have load-sharing coefficients equal to zero, so that for
any such node pi = pi (α). Using the net load profile p(α), the total loss rewrites as
1
H(α) = p(α)T L+ p(α),
2
and therefore {H(α)}α∈Rn is a family of random variables parametrized by the control α. Being
a quadratic form and being L+ a positive semi-definite matrix, it immediately follows that H(α)
is a non-negative random variable for any α ∈ Rn .
The next proposition shows how, leveraging the properties of the matrices L+ and Cα , the
expected total loss EH(α) rewrites as the sum of two contributions, one stochastic and one
deterministic which is not affected by the control α. Furthermore, the expected total loss is
rewritten as a quadratic function the load-sharing vector α. The proof of this proposition is
presented in Appendix A.
Proposition 2.1. Consider a network with nominal load profile µ such that 1T µ = 0 and
zero-mean stochastic fluctuations described by a covariance matrix Σ. Then, the expected total
loss using the control α is given by
1
EH(α) = EHs (α) + µT L+ µ,
2

(4)

where EHs (α) is the expected total loss due to the stochastic fluctuations. Furthermore, EHs (α) ≥
0 for every α ∈ Rn , and the following identity holds:
EHs (α) =

σ2 T +
1
(α L α) − 1T ΣL+ α + tr(ΣL+ ).
2
2

(5)

The first important remark is that the expected total loss EH(α) is a quadratic form in the
vector α ∈ Rn since it can be rewritten as
EH(α) =

σ2 T
α Aα − bT α + c,
2

(6)

where A = L+ is a positive semi-definite matrix, σ 2 > 0, b = L+ Σ1 ∈ Rn , and c = tr(ΣL+ )/2 +
1 T +
µ L µ ∈ R+ . Furthermore, we can already conclude that the nominal load profile µ has no
2
impact on the optimal control α, since it appears only in the constant term c.
A second important observation is that a closer look to (5) reveals that if the covariance
matrix is multiplied by a factor δ > 0, the resulting expected total loss is also δ times the original
one.
Before investigating what is the optimal load-sharing vector for a given network and covariance
structure of the noise, we argue here why the quantity EH(α) can be seen as a generalized notion
of effective resistance. In order to do so, we will first recall some classical definitions.
6

The effective resistance Ri,j between a pair of nodes i and j of the network G is defined as
the electrical resistance measured across nodes i and j when we look at G as electrical network
−1
in which resistors with conductance w1−1 , . . . , wm
are placed at the corresponding network edges.
Equivalently,
+
+
Ri,j := (ei − ej )T L+ (ei − ej ) = L+
i,i + Lj,j − 2Li,j .
The total effective resistance of a graph G is then defined as
n
1X
Rtot (G) :=
Ri,j .
2 i,j=1

The same quantity is also known as Kirchhoff index when the network G is such that all the edge
weights are equal to 1. This is a key quantity that measures how well connected the network
is and for this reason has been extensively studied and rediscovered in various contexts, such
as complex network analysis [13], theoretical chemistry, see [47] and references therein for an
overview, and probability theory [12]. In the context of electrical networks the total effective
resistance Rtot (G) is shown in [17] to be proportional to the average power dissipated in a resistive
network (G, w) when random i.i.d. currents with zero mean and unit variance are injected at the
nodes.
We can look at the network (G, w, α) with controllable loads introduced earlier as a flexible
transport network, where the load-sharing coefficients α1 , . . . , αn can be tuned to respond
optimally to specific stochastic load fluctuations. In this respect, we claim that the quantity
EHs (α) can be seen as a (rescaled) generalized total effective resistance that measures how
“robust” the network (G, w, α) is against stochastic load fluctuations with covariance structure Σ.
To further corroborate this claim, we now show that the expected total loss reduces to the classical
total effective resistance Rtot (G) in the special case where the load-sharing coefficients are all
equal, i.e., αi = 1/n for every i = 1, . . . , n, and the stochastic load fluctuations are i.i.d. random
variables with with zero mean and unit variance, i.e., Σ = I. Using (5), the expected total loss
due to fluctuations rewrites as
1 
n
1
1
Rtot (G)
EHs 1 = 2 (1T L+ 1) − 1T L+ 1 + tr(L+ ) = tr(L+ ) =
,
(7)
n
2n
2
2
2n
where in the last step we use the well-known identity Rtot (G) = n · tr(L+ ) proved by [23] that
relates the total effective resistance of a graph with its spectrum.
Rayleigh’s monotonicity principle [12] states that the pairwise effective resistance Ri,j is
a non-increasing function of the edge weights and, as a consequence, also the total effective
resistance Rtot (G) is. The following proposition, proved in Appendix A, shows that a similar
property also holds for H(α): regardless of the load-sharing vector α, the total power loss does
not increase when edges are added or weights are increased.
Proposition 2.2 (Monotonicity of total power loss). Let G be a weighted connected graph and
let G0 the graph obtained from G by increasing the weight of edge e = (i, j) by β > 0 or by adding
the edge e = (i, j) with weight β > 0. For any load-sharing vector α ∈ Rn and any realization of
the stochastic load fluctuations ω, the following inequality holds
0

HG (α) ≤ HG (α),
0

and, in particular, EHG (α) ≤ EHG (α).
7

3

Optimal load-sharing control

In this section we consider the problem of minimizing the expected total loss EH(α) given a
network G and a stochastic load covariance structure Σ.
Let B ⊆ V be the subset of k nodes with controllable loads. We focus first on the scenario
in which not all the nodes have controllable loads and thus assume 1 ≤ k < n. Besides the
constraint (3), we further add n − k constraints on the optimal load-sharing vector α ∈ Rn
to account for the absence of controllable loads in the nodes in V \ B, obtaining the following
constrained optimization problem in Rn :


minα∈Rn EHs (α)
(8)
s.t. 1T α = 1,


αv = 0, ∀ v ∈ V \ B.
Note we consider only the expected total loss due to stochastic load fluctuations in the objective
function, since in view of (4) it differs only by a constant from the expected total loss.
We henceforth assume that the k nodes with the controllable loads are those with labels
1, . . . , k, i.e., B = {1, . . . , k}. We can make this assumption without loss of generality as it
amounts to a relabelling the network nodes. If this is the case, the rows and columns of matrices
L, L+ , and Σ and the entries of the vector µ are also rearranged accordingly.
e ∈ Rk to
Let PB ∈ {0, 1}n×k be the binary matrix that maps any k-dimensional vector α
n
e = (α
e 1, . . . , α
e k , 0, . . . , 0) ∈ R . Such a matrix can be defined
the n-dimensional vector PB α
component-wise as (PB )i,j := δ{i=j} δ{i≤k} , for i = 1, . . . , n and j = 1, . . . , k, and has the following
structure:
 
I
PB = k ,
O
where Ik is the k × k identity matrix and O ∈ Rn−k×k is a matrix with all entries equal to zero.
In our first main result we present a closed-form expression for the optimal load-sharing
factors of k controllables.
Theorem 3.1 (Optimal load-sharing between k < n controllables). Consider a network with
balanced nominal load profile 1T µ = 0 in which the nodes in B = {1, . . . , k} have controllable
e ∗ 0), with
loads. The solution of the optimization problem (8) is the load-sharing vector α∗ = (α


1 + −1
1 + −1
−1 T + Σ1
∗
e = (LB ) 1 + Ik − (LB ) J (L+
α
∈ Rk ,
(9)
B ) PB L
tB
tB
σ2
+
+ −1
+
T +
T
where L+
B is the k × k principal submatrix of L , i.e., LB := PB L PB , and tB := 1 (LB ) 1 > 0.
If all the nodes with stochastic load fluctuations have controllable loads, i.e., S ⊆ B, then

e∗ =
α

Σ1
,
σ2

(10)

e ∗v = 0 for every v ∈ B \ S.
and, in particular, α
The involved expression (9) for the optimal load-sharing factors reflects the interplay that
exists between the network structure, the location of the controllable loads B, and the correlation
structure of the load fluctuations in determining the losses.
8

In the special case where all the nodes with stochastic load fluctuations have controllable loads
there is a nice interpretation for the optimal load-sharing coefficients: indeed α∗i is proportional
to how much the stochastic fluctuations of node i contribute in relative terms to the variance of
the total mismatch, since
P
n
X
Var(ω
)
+
Cov(ωi , ωj )
1
1
i
Pj6=ni
α∗i = 2 eTi Σ1 = 2
Σi,j =
σ
σ j=1
Var( i=1 ωi )
P
Var(ωi ) + j6=i Cov(ωi , ωj )
.
(11)
=P 
P
n
Var(ω
)
+
Cov(ω
,
ω
)
k
k
j
j6=k
k=1
We further remark that in the case of i.i.d. stochastic fluctuations in all the nodes, the second
term in (9) vanishes, since the vector Σ1 lies in the null space of L+ (being a multiple of 1), and,
therefore, the optimal control is equal to
(L+
)−1 1
e = T B+ −1 .
α
1 (LB ) 1
∗

In particular, it does not depend depend on the variance of the load fluctuations, but only on the
network structure and on the location B of the controllable loads, both encoded in the matrix
L+
B.

3.1

Full controllability

In this subsection we focus on the special case where the load is controllable in every node,
i.e., B = V , which is not covered in Theorem 3.1. Indeed, the proof method does not work in
this scenario due to the non-invertibility of the graph Laplacian L, and for this reason is treated
separately here.
The problem of minimizing the expected loss EH(α) when all the nodes have controllable
loads can be written as an optimization problem on Rn with a single constraint, namely
(
minα∈Rn EH(α)
(12)
s.t. 1T α = 1.
As mentioned earlier, thanks to Proposition 2.1 we can immediately conclude that the optimal
load-sharing vector does not depend on the vector µ, which appears only in the constant term in
the equality (4) for the expected loss.
The next theorem derives an analytical expression for the optimal solution of this optimization
problem.
Theorem 3.2 (Optimal load-sharing between n controllables). Consider a network G with n
nodes and a balanced nominal load profile 1T µ = 0. The solution of the optimization problem (12)
is the load-sharing vector α∗ given by
Σ1
α∗ = 2 .
σ
The highlight of this result is that in the scenario where all nodes have controllable loads, the
optimal control α∗ does not depend on the graph structure, but only on the covariance structure
of the fluctuations.
9

The same interpretation as in the special case S ⊆ B of Theorem 3.1 holds here for the
optimal load-sharing coefficients: α∗i is proportional to how much the stochastic fluctuations of
node i contribute in relative terms to the variance of the total mismatch, see (11).
In particular, when a node i does not have stochastic load fluctuations, then it is optimal not
to use the controllable load in that node not to do anything, since α∗i = 0 in view of the fact that
the i-th row of Σ is identically zero.
It immediately follows from Theorem 3.2 that when the the stochastic load fluctuations are
independent and identically distributed, the optimal load-sharing factors are all equal, namely
α∗ =

1
1.
n

Furthermore, the expected total loss due to stochastic fluctuations when using the optimal
load-sharing coefficient rewrites as


1 T
1
+
+
∗
tr(ΣL ) − 2 1 ΣL Σ1 .
EHs (α ) =
2
σ
In this special case, the fact that EHs (α∗ ) ≥ 0 can equivalently be proved as follows:
p
p
1T ΣL+ Σ1 = tr(ΣL+ ΣJ) ≤ tr((ΣL+ )2 )tr((ΣJ)2 ) ≤ tr(ΣL+ )2 tr(ΣJ)2 = tr(ΣL+ )tr(ΣJ) = σ 2 tr(ΣL+ ).
where both the inequalities leverage in a crucial way that all the matrices Σ, J, L+ are positive
semi-definite. The first inequality follows from the fact that in the space of positive semi-definite
matrices trace is a proper inner-product and thus obeys Cauchy-Schwarz inequality,
p
tr(AB) ≤ tr(A2 ) tr(B 2 )
∀ A, B positive semi-definite matrices.
The second inequality follows from the fact that tr(A2 ) ≤ tr(A)2 for any positive semi-definite
matrix A, obtained by applying Cauchy-Schwarz inequality with A = B.

4

Scaling properties of the expected total loss

In this section we explore the relation between the expected total loss and the number of
controllable loads. Even if the intuition suggests that the expected total loss should be a
decreasing function in the number of controllable loads, but this may not be true in general,
as the total loss depends both on the location of the controllable loads as well as on the the
load-sharing factors. For instance, in Section 5 we present a counterexample of a network in
which by adding the one controllable load and readjusting the load-share factors to be all equal,
the expected total power loss increases.
To get rid of these heterogeneities and obtain a more transparent result for the impact of
the number of controllable loads, we calculate the expected total loss for a fixed number of
controllable loads averaging on all their possible locations and assuming they share equally the
load. Consider an integer 1 ≤ k ≤ n and denote by Bk ⊂ Rn the collection of load-sharing vectors
with exactly k non-zero identical entries (and thus equal to 1/k, in view of (3)), namely
(
)
1X
Bk :=
ei : B ⊆ V, |B| = k ,
k i∈B
10

where ei ∈ Rn is the vector with the i-th entry equal to 1 and zero elsewhere.
Let Hk denote the expected total loss due to stochastic load fluctuations averaging over all
their possible locations of k controllable nodes share equally the load, i.e.,
1 X
Hk :=
EHs (α).
|Bk | α∈B
k


n

This average consists of |Bk | = k terms, one for each possible displacement of k controllable
nodes in a network with n nodes. The following theorem states an explicit expression for Hk that
makes the dependence on the number of controllable node k very explicit, showing in particular
that Hk is, up to a constant, proportional to 1/k.
Theorem 4.1 (Average total loss with k controllable loads). Consider a network G of n nodes
with balanced nominal profile load, i.e., 1T µ = 0. Then,
Hk = C1 +

C2
,
k

where C1 , C2 ∈ R are two constants that do not depend on k given by
1
tr(L+ )
C1 := tr(ΣL+ ) − σ 2
2
2n(n − 1)

and

C2 := σ 2

tr(L+ )
.
2(n − 1)

Both the constants C1 and C2 depend on the graph structure via L+ , on its size n, and on the
covariance matrix Σ. We remark that the constant C2 is always strictly positive, as tr(L+ ) > 0
(L+ being a positive semi-definite matrix) and σ 2 > 0, in view of (1).
We are interested in understanding how the expected total loss scales for large graphs. Assume
we can take a sequence of graphs {Gn }n∈N of growing size, |Vn | = n, and of covariance matrices
{Σn }n∈N with total variance σn2 = 1T Σn 1, so that the limit
(n − 1) tr(Σn L+
n)
n→∞
σn2 tr(L+
n)

γ := lim

exists. Note that the fact that inequality Hn ≥ 0 holds for every n ∈ N (see Proposition 2.1)
guarantees that γ ≥ 0. Under these assumptions, we have that the relative gain of having k
controllable loads with respect to a single one scales as

 

Hk
1
1
1
1
1
1
lim
=
+
=
+
.
2 tr(L+ )
−1
σn
(n−1) tr(Σn L+
n−1
n−1
n→∞ H1
n
n) k
1+γ
1+γ k
1 + limn→∞ n (n−1) tr(Σ L+ )
+
lim
+
n→∞
2
n
σ tr(L )
n

n

n

n

In the scenario where the load fluctuations are independent and identically distributed, regardless
of the graph structure, the asymptotic scaling reads
lim

n→∞

Hk
1
1
= + .
H1
2 2k

as

(13)

(n − 1) tr(Σn L+
(n − 1) σn2 tr(L+
n−1
n)
n)
=
lim
=
lim
= 1.
n→∞
n→∞
n→∞
σn2 tr(L+
n σn2 tr(L+
n
n)
n)
In Subsection 5.5 we show numerically that the way the expected total loss scales on average
with the number of controllable loads as stated in Theorem 4.1 is pretty accurate in general, even
without averaging over all possible locations of the controllable loads.
γ = lim

11

5

Numerical examples

5.1

Impact of covariance structure

The next two figures illustrates how the covariance structure influence the optimal control for
a small network of 12 nodes. Figure 1 covers the case in which all the nodes have stochastic
load fluctuations, i.e., S = V , while in Figure 2 we present a scenario where S ( V . In both
figures the nodes with controllable loads are drawn as squares, whose area is proportional to the
corresponding load-sharing factor.
3

3

3

9

9
8

2

9
8

2

2

4

4

5

5

12

12

5

7

7

8
4

12

7
6

6

6

11

11

(a) i.i.d. load fluctuations
e ∗ = (0.3167, 0.4826, 0.2007)
α

11

(b) Independent but not identically
distributed load fluctuations
e ∗ = (0.3335, 0.4876, 0.1789)
α

(c) Correlated load fluctuations
e ∗ = (0.4109, 0.3841, 0.2050)
α

Figure 1: Example of a stochastic loss transport network in which all nodes have stochastic load
fluctuations but with different covariance structure.
The correlation matrix Σc used for Figure 1 is generated at random in case (c) and is such
that σc2 = 1T Σc 1 = 50.42. The correlation matrix Σb one used in case (b) is obtained using the
diagonal entries of Σc . Note that both in case (a) and (b), the correlation matrices have been
rescaled so that their total variance is still equal to σa2 = σb2 = σc2 . The correlation matrices for
Figure 2 are obtained analogously.
3

3

3
9

9
8

2

9
8

2

4

5

12

7

4

5

12

7

5

12

7

6

6

11

(a) i.i.d. fluctuations
e ∗ = (0.2571, 0.5162, 0.2267)
α

8

2

4

6

11

(b) Independent but not identically
distributed fluctuations
e ∗ = (0.2624, 0.5317, 0.2059)
α

11

(c) Correlated fluctuations
e ∗ = (0.2011, 0.5504, 0.2485)
α

Figure 2: Example of a stochastic loss transport network with different covariance structure for
the load fluctuations. The nodes in black are those with stochastic load fluctuations.

12

5.2

Relative position of controllable loads and stochastic nodes

Figure 3 visualizes the optimal load-sharing factors in the scenario where the subset of stochastic
nodes is fixed, S = {1, 3, 5, 7, 9, 11}, as well as the covariance matrix, but the displacement of the
controllable loads, i.e., the subset B, changes.
12

12
1

1

11

12

1

2

10

10

2

2

3

9

10

3

9

4
8

11

11

3

9

8

4
4

8
5

7
5
6

(a) S ⊂ B = {1, 4, 5, 7, 10, 11}
EHs (α∗ ) = 2.8825

7

5

6

7

6

(b) S 6⊂ B = {1, 3, 5, 7, 9, 11}
EHs (α∗ ) = 3.3552

(c) S ⊂ B c , B = {2, 5, 6, 8, 9, 12}
EHs (α∗ ) = 3.8134

Figure 3: Example of a stochastic loss transport network with the same covariance structure for
the load fluctuations, but different displacement of the controllable loads and the corresponding
expected total loss.
Figure 3(a) presents the scenario in which S ⊆ B and the optimal control for the nodes in
B \ S = {5, 11} is equal to zero, as prescribed by Theorem 3.1. In the other two cases, (b) and
(c), we picked subsets B of controllable loads such that S 6⊆ B and the corresponding value of
the expected losses suggest it is optimal to place the controllable loads in the nodes affected by
stochastic fluctuations.

5.3

Negative load-sharing coefficients

The fact that a load-sharing
corresponding node absorbs
P factor is non-negative means that theP
part of power excess (if ni=1 ωi > 0) and balance out shortages (if ni=1 ωi < 0). In most of the
related work in primary response mechanisms and autonomous generation controls for power
grids, the load-sharing coefficients (often called participation factors) are in fact taken to be
non-negative, i.e., αv ≥ 0 for all v ∈ B. This assumption tacitly implies that all the controllable
generators and storage have “coordinated” actions, in the sense that they either all increase or
all decrease their power output.
In our formulation of the optimization problems (8) and (12) we do not make such an
assumption and load-sharing factors can also be negative, as long as the condition (3) is met.
This is crucial as for certain covariance structures of the load fluctuations (especially when there
are strong negative correlations) it is optimal to have negative load-sharing coefficients in some
nodes: we illustrate this fact for a small network illustrated in Figure 4.

13

2

3

4

1

Figure 4: A small network modeled by a graph with n = 4 nodes and m = 4 edges with unit
weights
Table 1 below lists the optimal load-sharing factors α∗ corresponding to different set of
controllable loads in the network in Figure 4 where the load fluctuations covariance structure is
assumed to be


1
0
0
−0.5
 0
1
0
−0.5 
.
Σ=
 0
0
1
−0.5 
−0.5 −0.5 −0.5
1
The best way for the controllable loads to respond to the negative correlations that the load
fluctuations have in this network is having the controllable load in node 4 taking actions “mirroring”
those of the other three nodes, in the sense that α∗4 < 0 while the load-sharing coefficients of the
other nodes in B are always positive.
B
α∗ (B)
{1, 4}
{6/5, 0, 0, −1/5}
{1, 2, 4}
{2/3, 2/3, 0, −1/3}
{1, 2, 3, 4} {1/2, 1/2, 1/2, −1/2}
Table 1: Sets of controllable loads for the network in Figure 4 with the corresponding optimal
load-sharing factors

5.4

Non-monotonicity of expected power loss when adding controllable loads

An extra controllable load always reduces the expected total power loss if the corresponding
optimal load-sharing vector α∗ is selected, since it corresponds to removing one constraints in the
optimization problem (8). However, if the chosen load-sharing factors of the augmented subset of
controllable loads are not the optimal ones, adding an extra controllable load does not necessarily
reduce the expected total power loss. We illustrate this with an example in which the control is
assumed to be equal-share between the controllable loads in B. Consider the network given in
Figure 5 and assume that the stochastic loads are i.i.d. with unit variance.

14

3
4
2
1
5

Figure 5: A small network modeled by a graph with n = 7 nodes and m = 9 edges with unit
weights
Table 2 below lists the expected total power losses for some subsets B of controllable loads
and compares them with those for some augmented subset B ∪ {4}.
B
{2}
{2, 5}
{2, 6, 7}

EH(B) EH(B ∪ {4})
2.6875
3.0625
2.0625
2.13194
1.79861
1.84375

Table 2: Expected total power losses for the network in Figure 5 assuming equal load-sharing
coefficients for some subsets B of controllable loads and then for the subsets augmented with an
extra node, namely B ∪ {4}.

5.5

Empirical evidence of the scaling law (13)

In the following example we consider the IEEE RTS 96-bus test network and track the expected
total loss while adding one by one controllable loads in random locations and assuming equal
share among the available controllable loads. As illustrated by Figure 6, the theoretical scaling
for the expected total loss scales with the number of controllable loads stated in Theorem 4.1
while averaging on all possible locations is in fact very accurate also for a single instance where
new controllable locations are randomly added.

15

70

70

60

60

50

50

40

40

10

20

30

40

50

60

70

10

(a) I.i.d. load fluctuations

20

30

40

50

60

70

(b) Independent but not identically distributed load
fluctuations

70

60

50

40

10

20

30

40

50

60

70

(c) Correlated load fluctuations

Figure 6: Theoretical scaling (in red) as predicted by Theorem 4.1 vs. empirical total loss (in
black) while adding controllable loads one by one in random locations.

6

Concluding remarks

In this paper we consider a stochastic lossy transport network in which some nodes have
controllable loads and derive a closed-form expression for the optimal control when aiming to
minimize the average total loss. The model is inspired by power systems where distributed energy
resources can be used as virtual storage to mitigate the fluctuations in the power generated by
renewable energy sources. Our analysis (i) uncovers a interplay between the network structure,
the location of the controllable loads and the covariance structure of the load fluctuations and (ii)
give insights in how much the average total loss can be reduced by adding a certain number of
controllable loads to the network.
The stylized mathematical model considered in the present paper gives useful insight in
understanding how the optimal displacement and operations of distributed energy resources are
affected by correlated load fluctuations. For this reason complements the recent efforts in the
electrical engineering community in upgrading the existing models for power grids to account
both for the intrinsic volatility of renewable energy generation and storage capabilities, see,
e.g., [22, 25, 26, 27, 28, 29].
In the present work we do not account for line limits, aiming to get a transparent insight in
the interplay between load uncertainty and optimal storage operations. Nonetheless, we remark
that the quadratic form we consider as a objective function can be rewritten as weighted sum of

16

the line flows squared and for this reason the optimal control should prevents line flows from
being too large and in particular above the corresponding line capacities.
The optimization problems considered in this paper focus on the network efficiency and do
not account for economic factors and neither prevents the excessive usage of controllable loads.
Both these factors can be accounted for by adding an additional term to the objective in the
optimization problems (12) and (8) of the form αT P α + q T α with P is the diagonal matrix with
positive entries P1,1 , . . . , Pn,n on the diagonal and q = (q1 , . . . , qn ) is a vector with non-negative
entries.
Let ξ be a positive real number that weights the relative importance of the penalty/cost term
αT P α + q T α with respect to the average total loss. The generalized optimization problem reads
(
minα∈Rn EH(α) + ξ (αT P α + q T α)
(14)
s.t. 1T α = 1,
The solution of this constrained optimization problem is
 σ2
−1 
−1 
 σ2
∗
+
+
L + ξP
L + ξP
α =
δ1 + I + δJ
(ξq + Σ1),
2
2
−1 −1
with δ := 1
+ ξP
1
, and can be obtained in a similar way of that in Theorem 3.1,
leveraging the fact that the matrix L+ + ξP is positive definite.
Such an extension can be instrumental to explore the trade-off between the best operations
for the network and the corresponding cost or penalties for the excessive usage of the controllable
nodes. This is particularly relevant for instance in the design of primary response mechanisms
and autonomous generation controls for power grids [1, 6, 18, 19, 30, 39, 42]. In these works
only economic criteria are considered to tune the load-sharing coefficients and therefore combining them with a network performance metric such as the expected total loss could be interesting.


T



σ2 +
L
2

Acknowledgment: This research is supported by NWO VICI grant 639.033.413 and NWO
Rubicon grant 680.50.1529.

Appendices
A

Properties of the total power loss

Proof. Proof of Proposition 2.1. Using the expression for Cα in (2), we can rewrite H(α) as
1
1
1
H(α) = (µ + Cα ω)T L+ (µ + Cα ω) = µT L+ µ + ω T CαT L+ Cα ω + µT L+ Cα ω.
2
2
2
Note that one of the terms on the RHS, namely 12 µT L+ µ, is not random and does not depend on
the control α. Define the random variable
1
1
Hs (α) := H(α) − µT L+ µ = ω T CαT L+ Cα ω + µT L+ Cα ω,
2
2
17

(15)

which describes precisely the contribution of the stochastic fluctuations to the transportation
losses. From the fact that L+ is a positive semi-definite matrix it follows that
1 T T +
1
ω Cα L Cα ω = (Cα ω)T L+ (Cα ω) ≥ 0 ∀ ω ∈ Rn .
2
2

(16)

Combining (15), (16), and the fact that E(µT L+ Cα ω) = µL+ Cα Eω = 0 yields that


1 T T +
ω Cα L Cα ω ≥ 0.
EHs (α) = E
2
Applying a classic result for quadratic forms of random vector, see Corollary 3.2b.1 in [31] and
using the fact that Eω = 0, we derive


1 T T +
1
(17)
E
ω Cα L Cα ω = tr(CαT L+ Cα Σ),
2
2
which yields
EHs (α) =

1
tr(CαT L+ Cα Σ),
2

since E(µT L+ Cα ω) = µL+ Cα Eω = 0.
We now derive identity (5). Recall the following well-known properties of the trace of matrix:
(i) The trace is invariant under cyclic permutations, i.e., for any r ∈ N
tr(A1 . . . Ar ) = tr(A2 . . . Ar A1 ) = . . . = tr(Ar A1 . . . Ar−1 ).
(ii) The trace of a matrix and of its transpose coincide, i.e., tr(A) = tr(AT );
(iii) The trace of the outer product of two vectors is their inner product, namely
tr(vwT ) = tr(v ⊗ w) = v T w.
First note that we can rewrite
CαT L+ Cα = (I − α 1T )T L+ (I − α 1T ) = (I − 1 αT )(L+ − L+ α 1T )
= L+ − L+ α 1T − 1 αT L+ + 1 αT L+ α 1T .

(18)

The aforementioned properties of the trace yield
(ii)

(i)

tr(ΣL+ α 1T ) = tr((ΣL+ α 1T )T ) = tr(1αT L+ Σ) = tr(Σ1αT L+ ),

(19)

and
(iii)

(iii)

tr(Σ1αT L+ ) = tr(Σ1(L+ α)T ) = tr((Σ1) ⊗ (L+ α)) = (Σ1)T (L+ α) = 1T ΣL+ α.

(20)

By combining all these equalities and exploiting the linearity of the trace operator, we obtain
(18)

tr(ΣCαT L+ Cα ) = tr(ΣL+ ) − tr(ΣL+ α 1T ) − tr(Σ1 αT L+ ) + tr(Σ1 αT L+ α 1T )
(19)

= tr(ΣL+ ) − 2 · tr(Σ1 αT L+ ) + (αT L+ α) · tr(Σ1 1T )

(20)

= tr(ΣL+ ) − 2 · (1T ΣL+ α) + (αT L+ α) · tr(ΣJ),

where in the second step we also used the fact that αT L+ α is a scalar.
18

Proof. Proof of Proposition 2.2 Assume that e = (i, j) ∈ (V × V ) is the edge with weight
β > 0 that has been added to G or whose edge weight has been increased by β > 0 and let
me = (ei − ej ) ∈ Rn be the corresponding non-weighted incidence vector. In both cases the
Laplacian matrix of the newly obtained graph G0 can be written as
LG0 = LG + β me mTe
and, using the generalized version of the Sherman-Morrison formula in [35], we get
+
L+
G0 = LG −

β −1

1
T +
L+
G me me LG .
+ mTe L+
m
e
G

We can thus rewrite the total loss corresponding to any net power injection vector p(α) as
1
0
HG (α) = p(α)T LG0 p(α)
2
1
1
T +
= p(α)T L+
p(α)T L+
G p(α) −
G me me LG p(α)
−1
2
2(β + mTe L+
m
)
e
G
T +
2
(me LG p(α))
,
= HG (α) −
2(β −1 + mTe L+
G me )
T +
2
and conclude by noticing that mTe L+
G me ≥ 0 and (me LG p(α)) ≥ 0.

B

Proof of Theorem 3.1

In the proof of Theorem 3.1 we will make use of the so-called block matrix inversion formula,
which is stated in the next lemma.
Lemma

 B.1 (Block matrix inversion formula). Consider a matrix with the block structure
A B
. If both A and D − CA−1 B are non-singular matrices, then
C D


A B
C D

−1


=


A−1 − A−1 B(D − CA−1 B)−1 CA−1 −A−1 B(D − CA−1 B)−1
.
−(D − CA−1 B)−1 CA−1
(D − CA−1 B)−1

Proof. Proof of Theorem 3.1. For any vector α ∈ Rn which is such that αi = 0 for every node
e Using this
e ∈ Rk such that α = PB α.
i ∈ V \ B, there exists a unique k-dimensional vector α
T
correspondence and the fact that the nominal load profile is balanced, i.e., 1 µ = 0, we can
rewrite
2
1
(4) σ
e T PBT L+ PB α)
e =
e − 1T ΣL+ PB α
e + tr(ΣL+ ).
EHs (α) = EH(PB α)
(α
2
2
Therefore the n-dimensional optimization problem (8) rewrites as a k-dimensional optimization
problem with a single constraint, namely

2
 min σ (α
e T PBT L+ PB α)
e − 1T ΣL+ PB α
e
k 2
(21)
e
α∈R

T
e = 1.
s.t. 1 α
19

T +
k×k
The matrix L+
is positive definite, as for any vector v ∈ Rk , v 6= 0,
B := PB L PB ∈ R
 
+ v
T +
T T +
> 0,
v LB v = v PB L PB v = (v 0) L
0

where the last inequality follows from the fact that the vector (v 0) is not a multiple of the vector
1 and thus does not lie in the null space of L+ . The optimization problem in (21) has then a
unique solution, since the corresponding Hessian is positive definite.
Let γ ∈ R be the Lagrange multiplier γ associated with the unique equality constraint of the
optimization problem (21). The associated Lagrangian associated is
e γ) =
L(α,

σ2 T T +
e − γ(1T α
e − 1).
e − (PBT L+ Σ1)T α
e (PB L PB )α
α
2

Setting b := PBT L+ Σ1 ∈ Rk , the optimality conditions read
(
e − γ1 = b,
σ 2 L+
Bα
T
e = 1,
1 α
or, equivalently, in matrix form
 2 +
   
e
σ LB −1
α
b
=
.
T
γ
1
1
0

(22)

Being positive definite, L+
B is invertible and its inverse is also positive definite, which means that
+ −1
+
T
tB := 1 (LB ) 1 > 0. In view of the fact that σ 2 t−1
B 6= 0 and LB is invertible, we can use the
block matrix inversion formula given in Lemma B.1 to obtain
1


 1
1
+ −1
+ −1
+ −1
 2 +
−1
(L )
(L ) 1
I − J(LB )
σ LB −1

 σ2 B
tB
tB B
=

.
2
T
σ
1 T + −1
1
0
1 (LB )
tB
tB
The solution of the linear system (22) then reads
1


 1
1
+ −1
+ −1
+ −1
   2 +
−1  
 
(L
J(L
(L
)
I
−
)
)
1
B
B
e
σ LB −1
α
b
 σ2 B
 b
t
t
B
B
=
=

σ2
1 T + −1
γ
1
1
1T
0
1 (LB )
tB
tB
e ∈ Rk is given by
and thus the optimal load-sharing vector α


−1
(L+
1
1
+ −1
∗
B)
e =
α
I − J(LB )
b + (L+
)−1 1
2
σ
t
tB B


1
1 + −1
b
+ −1
+ −1
= (LB ) 1 + (LB )
I − J(LB )
tB
tB
σ2


1 + −1
1 + −1
Σ1
+ −1
+ −1
= (LB ) 1 + (LB ) − (LB ) J(LB )
PBT L+ 2
tB
tB
σ


1
1 + −1
−1
−1 T + Σ1
= (L+
(LB ) J (L+
.
B) 1 + I −
B ) PB L
tB
tB
σ2
20

We now focus on the special case where S ⊆ B and prove identity (10). Rewrite L+ as a block
matrix
 +

LB
L+
C
,
T
(L+
L+
C)
Bc
+
k×k
k×n−k
n−k×n−k
with L+
, L+
and L+
. Note that L+
B ∈ R
C ∈ R
Bc ∈ R
B and LB c are symmetric
matrices, since L+ is. Note that this is consistent with the former definition of L+
B , since
 +

LB
L+
T +
T
T
C
PB L PB = (Ik | O)L+ (Ik | O) = (Ik | O)
(Ik | O) = L+
+
T
B.
(L+
)
L
c
C
B

First note that
−1 T +
(L+
B ) PB L

=

−1
(L+
B ) (Ik | O)



L+
L+
B
C
+
T
(L+
)
L
C
Bc



+
+
+ −1 +
−1
k×n
= (L+
B ) (LB | LC ) = (Ik | (LB ) LC ) ∈ R

From the assumption S ⊆ B it follows that the covariance matrix can be rewritten as


ΣB O
Σ=
,
O O
where ΣB ∈ Rk×k is itself a covariance matrix. Trivially 1T (ΣB | O)1 = tr(ΣB J) = tr(ΣJ) = σ 2 .
Furthermore,




ΣB O
+ −1 T +
+ −1 T + ΣB O
+ −1 +
(LB ) PB L Σ = (LB ) PB L
= (Ik | (LB ) LC )
= (ΣB | O).
O O
O O
e ∗ then rewrites as
The optimal control α


−1
−1
(L+
(L+
∗
−1 T + Σ1
B)
B) J
e =
α
1 + Ik −
(L+
B ) PB L
tB
tB
σ2


(L+ )−1
(L+ )−1 J
1
= B
1 + Ik − B
(ΣB | O) 2
tB
tB
σ
+ −1
+ −1
1
1T (ΣB | O)1
(L )
(L )
= B
1 + (ΣB | O) 2 − B
1
tB
σ
tB
σ2
+ −1
+ −1
(L )
Σ1 (L )
= B
1+ 2 − B
1
tB
σ
tB
Σ1
= 2.
σ

C

Proof of Theorem 3.2

In (6) we already rewrote the objective function of the optimization problem (12) as
EHs (α) =

σ2 T
α Aα − bT α + c.
2

where A = L+ , b = L+ Σ1 ∈ Rn , and c = tr(ΣL+ )/2 ∈ R+ and notice that, for the purpose of
solving the optimization problem (12), we can ignore the constant term c.
21

We can consider the representation with respect to the orthogonal basis of eigenvectors
, v2 , . . . , vn of L+ of a vector α ∈ R, namely
α = β1 1 + β2 v2 + . . . + βn vn ,

(23)

for some β1 , . . . , βn ∈ R. From the constraint 1T α = 1, it immediately follows that β1 = 1/n.
Indeed,
nβ1 = β1 1T 1 = β1 1T 1 + β2 1T v2 + . . . + βn 1T vn = 1T (β1 1 + β2 v2 + . . . + βn vn ) = 1T α = 1.
Let v2 , . . . , vn be an orthonomal basis of eigenvectors of the graph Laplacian L corresponding to
the non-zero eigenvalues λ2 , . . . , λn and define the real coefficients κ2 , . . . , κn as
κi := hΣ1, vi i = 1T Σvi ,

i = 2, . . . , n.

Leveraging the latter definition and the representation (23), we can rewrite
T

+

T

+

α L α = (β1 1+β2 v2 +. . .+βn vn ) L (β1 1+β2 v2 +. . .+βn vn ) =

n
X

β12 1T L+ 1+

βi2 viT L+ vi

i=2

=

n
X
β2
i

i=2

λi

and
T

T

+

T

+

b α = 1 ΣL α = 1 ΣL (β1 1 + β2 v2 + . . . + βn vn ) =

n
X
βi
i=2

λi

T

1 Σvi =

n
X
i=2

βi

κi
.
λi

By combining the latter two identities, the objective function rewrites as
n
n
σ2 T
σ 2 X βi2 X κi
T
EH(α) = α Aα − b α =
−
βi .
2
2 i=2 λi
λi
i=2

(24)

The optimization problem (12) is therefore equivalent to a unconstrained problem in Rn−1 whose
variables are β2 , . . . , βn . Define the function g : Rn−1 → R as
n
n
σ 2 X βi2 X κi
g(β2 , . . . , βn ) :=
−
βi .
2 i=2 λi
λi
i=2

The gradient of the function g is the vector


κi
2 βi
∇g(β2 , . . . , βn ) = σ
−
λi λi i=2,...,n
−1
and its Hessian matrix of g is the diagonal matrix H(g) = σ 2 · diag(λ−1
2 , . . . , λn ). The Hessian
H(g) is constant as it does not depend on β2 , . . . , βn . Furthermore, H(g) is a positive definite
matrix since all its diagonal terms are positive, in view of the fact that λi > 0 for i = 2, . . . , n
and that σ 2 > 0. This means that any critical point satisfying ∇g(β2 , . . . , βn ) = 0 would then
be a minimum for the function g. The optimality condition ∇g(β2 , . . . , βn ) = 0 can be solved
explicitly, yielding
κi
βi∗ = 2 ,
i = 2, . . . , n.
σ

22

,

As a consequence the optimal load-sharing coefficients are unique and they are given by
n
1
1 X
α = 1+ 2
κi vi .
n
σ i=2
∗

Set v1 =

D

and κ1 = hΣ1, v1 i = 1T Σv1 = 1T Σ √1n 1. Note that α∗ rewrites as
!
!
n
n
n
X
X
X
1
κ
1
1
Σ1
1
1
α∗ = 2 √ 1 + 2
κi vi = 2
κi vi =
hΣ1, vi i vi = 2 .
σ n
σ i=2
σ
a i=1
σ
i=1

√1 1
n

Proof of Theorem 4.1

The starting point of the proof are two identities that leverage the properties of the pseudoinverse
L+ of the graph Laplacian. Firstly,
n
X 
 X X 
X

X
X
+
+
+
ei = L
ei = L
ei 1{i∈B}
L
B⊆V : |B|=k

i∈B
+

=L

B⊆V : |B|=k i∈B
n
X
X

ei

i=1

1{i∈B}



i=1 B⊆V : |B|=k

n
X
n
+
ei
=L


−1 
k−1

i=1

B⊆V : |B|=k





n
n−1 +
n − 1 + X 
=
L
ei =
L 1 = 0,
k−1
k
−
1
i=1

(25)

where we use the fact that in a graph with n nodes, each nodes belong to exactly
of k nodes. We further claim that
X  X T  X  n − 2
e i L+
ei =
tr(L+ ),
k
−
1
i∈B
i∈B
B⊆V : |B|=k

with the convention that n−2
= 0. Since
n−1
 X T  X  X
X
e i L+
ei =
eTi L+ ei +
eTi L+ ej ,
i∈B

i∈B

we can rewrite the LHS of (26) as
X  X T  X 
ei L+
ei =
B⊆V : |B|=k

i∈B

i∈B

i∈B

X
B⊆V : |B|=k

n−1
k−1



subsets
(26)

i,j∈B, i6=j

X


eTi L+ ei +

i∈B

X
B⊆V : |B|=k



X


eTi L+ ej .

i,j∈B, i6=j

(27)
The first term on the RHS of (27) can be rewritten as
n
 X
 X
X X
eTi L+ ei =
eTi L+ ei
B⊆V : |B|=k

i∈B

i=1

1{i∈B}



B⊆V : |B|=k

n 
X



 n
n−1 X T +
n−1 T +
=
e L ei =
e L ei
k−1 i
k − 1 i=1 i
i=1

 n


n−1 X +
n−1
=
L =
tr(L+ ),
k − 1 i=1 i,i
k−1
23

while the second term on the RHS of (27) is equal to


X
B⊆V : |B|=k

X

eTi L+ ej



X

=

i,j∈B, i6=j



=

P

X

i6=j

i6=j


n 
X
n−2 T +
=
e L ej
k−2 i
i=1

eTi L+ ej = −tr(L+ ), which immediately follows from

eTi L+ ej +

i6=j
n
X

=

1{i∈B, j∈B}






n−2 X T +
n−2
e L ej = −
tr(L+ ).
k − 2 i6=j i
k−2

In the last step we used the fact that
eTi L+ ej + tr(L+ ) =

X
B⊆V : |B|=k

i6=j



X

eTi L+ ej

ek

T

L+

X

L+
i,i =

i
n
X

k=1

X

eTi L+ ej +

X

eTi L+ ei

i

i6=j


ek = 1T L+ 1 = 0.

k=1

Hence, (27) rewrites as
X
B⊆V : |B|=k

X

ei

T

L

+

i∈B

X



ei =

n − 1





n−2 
n−2
+
−
tr(L ) =
tr(L+ ),
k−1
k−2
k−1

i∈B

which concludes the proof of identity (26).
P
Each load-sharing vector α ∈ Bk can be written as α = k1 i∈B ei , for some B ⊆ V , |B| = k.
In view of Proposition 2.1, the expected total loss due to stochastic fluctuations when using this
load-sharing vector is given by
!T
!
!
X
X
σ2 1 X
1
1
1
EHs (α) =
ei
L+
ei − 1T ΣL+
ei + tr(ΣL+ ).
2 k i∈B
k i∈B
k i∈B
2
Therefore,
1
1
Hk = tr(ΣL+ ) +
2
|Bk |

"
X
B⊆V : |B|=k

2

σ
2

1X
k

i∈B

ei

T

L+

1 X
k

i∈B



ei − 1T ΣL+

1 X
k

ei



#
.

i∈B

Using (25) and (26), we get

n−2
+
+
1
1
n−k
+
2 tr(L ) k−1
+
2 tr(L )

Hk = tr(ΣL ) + σ
=
tr(ΣL
)
+
σ
n
2
2 k2 k
2
2 k · n(n − 1)
+ 
1
tr(L+ ) 1 − nk
1
1 1
+
2 tr(L )
= tr(ΣL+ ) + σ 2
=
tr(ΣL
)
+
σ
−
.
2
2n(n − 1) nk
2
2(n − 1) k n

References
[1] D. Apostolopoulou, P. Sauer, and A. Dominguez-Garcia. Automatic Generation Control and
Its Implementation in Real Time. In 2014 47th Hawaii International Conference on System
Sciences, pages 2444–2452. IEEE, jan 2014.
24

[2] A. Bejan, R. Gibbens, and F. Kelly. Statistical aspects of storage systems modelling in
energy networks. In 46th Annual Conference on Information Sciences and Systems (CISS),
pages 1–6, Princeton, NJ, mar 2012. IEEE.
[3] D. Bhaumik, D. Crommelin, and B. Zwart. A computational method for optimizing storage
placement to maximize power network reliability. In To appear in the Proceeding of the 2016
Winter Simulation Conference, 2015.
[4] D. Bienstock, M. Chertkov, and S. Harnett. Chance-Constrained Optimal Power Flow:
Risk-Aware Network Control under Uncertainty. SIAM Review, 56(3):461–495, 2014.
[5] D. Bienstock, C. Matke, G. Munoz, and S. Yang. Robust linear control of nonconvex battery
operation in transmission systems. Preprint at arXiv:1610.09432, 2016.
[6] M. Chertkov and Y. Dvorkin. Chance Constrained Optimal Power Flow with Primary
Frequency Response. Preprint at arXiv:1703.06724, 2017.
[7] J. Cruise, L. Flatley, and S. Zachary. Impact of storage competition on energy markets.
Preprint at arXiv:1606.05361, 2016.
[8] J. Cruise, R. Gibbens, and S. Zachary. Optimal control of storage for arbitrage, with
applications to energy systems. In 48th Annual Conference on Information Sciences and
Systems (CISS), number i, 2014.
[9] J. Cruise and S. Zachary. The optimal control of storage for arbitrage and buffering, with
energy applications. Preprint at arXiv:1509.05788, 2015.
[10] C. F. Daganzo. Fundamentals of Transportation and Traffic Operations. Emerald Group
Publishing Limited, sep 1997.
[11] D. Deka, H. Nagarajan, and S. Backhaus. Optimal Topology Design for Disturbance
Minimization in Power Grids. Preprint at arXiv:1703.00558, 2017.
[12] P. G. Doyle and J. L. Snell. Random Walks and Electric Networks. American Mathematical
Monthly, 94(January):202, 2000.
[13] W. Ellens, F. Spieksma, P. Van Mieghem, a. Jamakovic, and R. Kooij. Effective graph
resistance. Linear Algebra and its Applications, 435(10):2491–2506, nov 2011.
[14] N. Gast, J.-Y. Le Boudec, A. Proutière, and D.-C. Tomozei. Impact of storage on the
efficiency and prices in real-time electricity markets. Proceedings of the 4th international
conference on Future energy systems - e-Energy ’13, 2:15, 2013.
[15] N. Gast, D.-C. Tomozei, and J.-Y. Le Boudec. Optimal Storage Policies with Wind Forecast
Uncertainties. ACM SIGMETRICS Performance Evaluation Review, 40(3):1–5, 2012.
[16] N. Gast, D.-C. Tomozei, and J.-Y. Le Boudec. Optimal Generation and Storage Scheduling
in the Presence of Renewable Forecast Uncertainties. IEEE Transactions on Smart Grid,
5(3):1328–1339, 2014.

25

[17] A. Ghosh, S. Boyd, and A. Saberi. Minimizing Effective Resistance of a Graph. SIAM
Review, 50(1):37–66, 2008.
[18] S. Guggilam, C. Zhao, E. Dall’Anese, Y. Chen, and S. Dhople. Engineering Inertial and
Primary-frequency Response for Distributed Energy Resources. Preprint at arXiv:1706.03612,
2017.
[19] S. Guggilam, C. Zhao, E. Dall’Anese, Y. Chen, and S. Dhople. Primary frequency response
with aggregated DERs. In 2017 American Control Conference (ACC), pages 3386–3393.
IEEE, may 2017.
[20] L. Guo, C. Liang, and S. Low. Monotonicity Properties and Spectral Characterization
of Power Redistribution in Cascading Failures. Communication, Control, and Computing
(Allerton), 2017 55th Annual Allerton Conference on, 2017.
[21] J. Johnson and M. Chertkov. A majorization-minimization approach to design of power
transmission networks. In 49th IEEE Conference on Decision and Control (CDC), pages
3996–4003. IEEE, 2010.
[22] Y. Kanoria, A. Montanari, D. Tse, and B. Zhang. Distributed storage for intermittent energy
sources: Control design and performance limits. In 49th Annual Allerton Conference on
Communication, Control, and Computing, pages 1310–1317, 2011.
[23] D. J. Klein and M. Randić. Resistance distance. Journal of Mathematical Chemistry,
12(1):81–95, 1993.
[24] C. Lai and S. Low. The redistribution of power flow in cascading failures. In 2013 51st
Annual Allerton Conference on Communication, Control, and Computing (Allerton), pages
1037–1044. IEEE, oct 2013.
[25] W. Lin and E. Bitar. Decentralized control of distributed energy resources in radial distribution systems. In 2016 IEEE International Conference on Smart Grid Communications
(SmartGridComm), pages 296–301. IEEE, nov 2016.
[26] W. Lin and E. Bitar. Performance bounds for robust decentralized control. In 2016 American
Control Conference (ACC), pages 4323–4330. IEEE, 2016.
[27] W. Lin and E. Bitar. Decentralized Stochastic Control of Distributed Energy Resources.
IEEE Transactions on Power Systems, 8950(c):1–1, 2017.
[28] A. Lorca and X. A. Sun. Multistage Robust Unit Commitment with Dynamic Uncertainty
Sets and Energy Storage. IEEE Transactions on Power Systems, 32(3):1678–1688, 2017.
[29] Á. Lorca, X. A. Sun, E. Litvinov, and T. Zheng. Multistage Adaptive Robust Optimization
for the Unit Commitment Problem. Operations Research, 64(1):32–51, feb 2016.
[30] M. Lubin, Y. Dvorkin, and S. Backhaus. A Robust Approach to Chance Constrained
Optimal Power Flow With Renewable Generation. IEEE Transactions on Power Systems,
31(5):3840–3849, sep 2016.

26

[31] A. Mathai and S. Provost. Quadratic forms in random variables: Theory and applications.
Statistics: textbooks and monographs, Vol. 126, 1992.
[32] J. Mathias, A. Bušić, and S. Meyn. Demand Dispatch with Heterogeneous Intelligent Loads.
In 50th Annual Hawaii International Conference on System Sciences (HICSS), Waikoloa,
HI, United States, 2017.
[33] J. Mathias, R. Kaddah, A. Bušić, and S. Meyn. Smart Fridge / Dumb Grid? Demand
Dispatch for the Power Grid of 2020. In 2016 49th Hawaii International Conference on
System Sciences (HICSS), pages 2498–2507. IEEE, jan 2016.
[34] C. Matke, D. Bienstock, G. Muñoz, S. Yang, D. Kleinhans, and S. Sager. Robust optimization
of power network operation: storage devices and the role of forecast errors in renewable
energies. In H. Cherifi, S. Gaito, W. Quattrociocchi, and A. Sala, editors, Complex Networks
& Their Applications V, volume 693 of Studies in Computational Intelligence, pages 809–820.
Springer International Publishing, Cham, 2017.
[35] C. D. Meyer, Jr. Generalized Inversion of Modified Matrices. SIAM Journal on Applied
Mathematics, 24(3):315–323, may 1973.
[36] S. Meyn, P. Barooah, A. Bušić, Y. Chen, and J. Ehren. Ancillary service to the grid using
intelligent deferrable loads. IEEE Transactions on Automatic Control, 60(11):2847–2862,
2015.
[37] S. Meyn, P. Barooah, A. Bušić, and J. Ehren. Ancillary service to the grid from deferrable
loads: The case for intelligent pool pumps in Florida. In 52nd IEEE Conference on Decision
and Control, pages 6946–6953. IEEE, dec 2013.
[38] G. Ranjan, Z.-L. Zhang, and D. Boley. Incremental Computation of Pseudo-Inverse of
Laplacian. In COCOA 2014: 8th International Conference on Combinatorial Optimization
and Applications, volume 5573, pages 729–749. Springer International Publishing, 2014.
[39] L. Roald, G. Andersson, S. Misra, M. Chertkov, and S. Backhaus. Optimal power flow
with wind power control and limited expected risk of overloads. In 2016 Power Systems
Computation Conference (PSCC), volume 60, pages 1–7. IEEE, 2016.
[40] L. Roald, S. Misra, T. Krause, and G. Andersson. Corrective Control to Handle Forecast
Uncertainty: A Chance Constrained Optimal Power Flow. IEEE Transactions on Power
Systems, 32(2):1–1, 2016.
[41] E. Sjodin, D. Gayme, and U. Topcu. Risk-mitigated optimal power flow for wind powered
grids. In 2012 American Control Conference (ACC), pages 4431–4437. IEEE, jun 2012.
[42] K. Sundar, H. Nagarajan, M. Lubin, L. Roald, S. Misra, R. Bent, and D. Bienstock. Unit
commitment with N-1 Security and wind uncertainty. In 2016 Power Systems Computation
Conference (PSCC), pages 1–7. IEEE, jun 2016.
[43] J. A. Taylor and J. L. Mathieu. Uncertainty in Demand ResponseIdentification, Estimation,
and Learning. In The Operations Research Revolution, pages 56–70. INFORMS, sep 2015.

27

[44] P. Van De Ven, N. Hegde, L. Massoulié, and T. Salonidis. Optimal control of end-user energy
storage. IEEE Transactions on Smart Grid, 4(2):789–797, 2013.
[45] P. Van Mieghem. Graph Spectra for Complex Networks. CUP, 2011.
[46] P. Whittle. Networks: Optimisation and Evolution. Cambridge University Press, 2007.
[47] B. Zhou and N. Trinajstić. A note on Kirchhoff index. Chemical Physics Letters, 455(13):120–123, mar 2008.
[48] A. Zocca and B. Zwart. Minimizing heat loss in DC networks using batteries. In 2016 54th
Annual Allerton Conference on Communication, Control, and Computing (Allerton), pages
1306–1313. IEEE, sep 2016.

28

