Curvature analysis of the power grid: congestion, congestion management, and line tripping cascade anticipation

PDF

Download

Open document

Flip pages

Summarize

Download a page range

Download transcript

Copy asset link

Request this asset

Transcript (if available)

Content Curvature Analysis of the Power Grid:
Congestion, Congestion Management, and Line Tripping Cascade
Anticipation
by
Eugenio Grippo
A Dissertation Presented to the
FACULTY OF THE USC GRADUATE SCHOOL
UNIVERSITY OF SOUTHERN CALIFORNIA
In Partial Fulllment of the
Requirements for the Degree
DOCTOR OF PHILOSOPHY
(ELECTRICAL ENGINEERING)
May 2020
Copyright 2020 Eugenio Grippo
Dedication
To my mom, my dad, my brother, my sister, my lovely wife and my beloved grandma
who is no longer with us but I am certainly sure she might be saying \ te felicito nene "...
ii
Acknowledgments
I want to gratefully say thanks to my advisor Prof. Edmond Jonckheere who has always
been there to teach me, help me and push me forward along this usually lonely road of
studying. I also want to extend my gratitude to my very good friend Anantha Karthikeyan
who helped me a lot and gave me courage in many dicult times I went through along
all these years as a PhD student. I remember and say thanks to my hometown professors
Julio Braslavsky, Anibal Zanini and Nicolas Cole who always encouraged me to keep on
studying and helped me at the beginning of my PhD; and I am condent I will see them
again either back home or in an international conference. I do not want to forget my friends
Reza Banirazi and Chi Wang who always helped me when I was in need. Lastly, I would like
to remember a very nice person who helped me a lot at USC and he is no longer with us,
his name is Tim Boston and I will always miss to stop by his oce and stretch his hand and
hear him saying \hello young man". Finally, I want to say thanks to the person that always
helped me to sort many dierent problems I faced in the academia and life itself during all
these years at USC, her name is Diane Demetras and I will always be indebted to her and
her sweet family.
iii
Contents
Dedication......................................................................................................................... ii
Acknowledgments ............................................................................................................. iii
List of Tables .................................................................................................................... vi
List of Figures................................................................................................................... vii
Abstract............................................................................................................................ ix
Introduction...................................................................................................................... 1
1 Essential background................................................................................................... 5
1.1 Graph Theory ..................................................................................................... 5
1.2 Resistive Networks .............................................................................................. 5
1.3 Curvature of Graphs ........................................................................................... 7
1.4 Eective Resistance Curvature of Prototypical Euclidean and Hyperbolic Graphs 8
1.4.1 Euclidean Graphs . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8
1.4.2 Hyperbolic Graphs . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9
1.4.3 A Counterexample: Modied Four-Point Graph . . . . . . . . . . . . 11
1.4.4 Growth-Preferential Attachment Graph . . . . . . . . . . . . . . . . . 11
1.5 Hilbert Space Insight .......................................................................................... 12
1.6 Main Result ........................................................................................................ 16
1.7 Curvature Detection............................................................................................ 20
1.8 Data Network...................................................................................................... 20
2 Curvature and the Power Grid.................................................................................... 23
2.1 Resistive Network Models of Power Flows.......................................................... 24
2.2 Curvature and Congestion .................................................................................. 26
2.2.1 Local Curvature Notion: The Ollivier-Ricci Curvature . . . . . . . . . 27
2.2.2 Wasserstein Distance, Earth Mover's Distance (EMD) and ORC . . . 27
2.2.3 Global Curvature Notion: Eective Resistance Curvature . . . . . . . 28
2.3 New Curvature Centrality Concept..................................................................... 31
2.4 Ollivier-Ricci Curvature Driven Optimal Power Flow (ORC-OPF).................... 32
2.4.1 Negative Curvature Detection via ORC . . . . . . . . . . . . . . . . . 32
2.4.2 FACTS, Storage Devices and the full ORC-OPF . . . . . . . . . . . . 34
2.4.3 Eective Resistance Curvature Driven Optimal Power Flow (ERC-OPF) 38
2.5 Grid Entropy Measure ........................................................................................ 39
2.6 Line Rating Considerations................................................................................. 41
3 Anticipation of Line Tripping Cascade by Grid Curvature and Entropy .................... 44
3.1 Ollivier-Ricci Curvature Approach to Detect Power Grid Areas Under Risk (
Fault / Overheating ).......................................................................................... 44
3.1.1 Negative Curvature Detection via ORC . . . . . . . . . . . . . . . . . 45
iv
3.1.2 ORC Identication of Areas/Lines/Corridors under Risk . . . . . . . 45
3.1.3 First Simulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45
3.1.4 Second Simulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47
3.2 Eective Resistance Curvature Approach to Detect Power Grid Areas Under
Risk ( Fault / Overheating )............................................................................... 48
3.2.1 Negative Curvature Detection via ERC . . . . . . . . . . . . . . . . . 52
3.2.2 ERC Identication of Areas/Lines/Corridors under Risk . . . . . . . 52
3.2.3 First Simulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52
3.2.4 Second Simulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54
4 Local/Nodal Entropy .................................................................................................. 61
5 Conclusions ................................................................................................................. 63
5.1 Research Summary.............................................................................................. 63
5.2 Future Work|Near Term ................................................................................... 65
5.2.1 Future Work|Longer Term . . . . . . . . . . . . . . . . . . . . . . . 68
5.2.2 Potential Biomedical Application . . . . . . . . . . . . . . . . . . . . 69
References......................................................................................................................... 70
v
List of Tables
1 Column indices of nonvanishing entries of rows 1 to 8 of D . . . . . . . . . . 14
2 Ollivier-Ricci Curvature Driven OPF \Main Steps" . . . . . . . . . . . . . . 36
3 Total cost functions values (AC Model with conventional OPF and with ORC-
OPF implementation) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37
4 Total cost functions values (AC Model with conventional OPF and with ERC-
OPF implementation) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39
5 Global curvature values vs grid entropy . . . . . . . . . . . . . . . . . . . . . 40
vi
List of Figures
1 Euclidean tiling (6 sides core cell) . . . . . . . . . . . . . . . . . . . . . . . . 8
2 R
lrp
=R
e
ratio evaluated throughout an Euclidean tiling (6 sides core cell) . . 9
3 Hyperbolic tiling (7 sides core cell) . . . . . . . . . . . . . . . . . . . . . . . 9
4 R
lrp
=R
e
ratio evaluated throughout a hyperbolic tiling (7 sides core cell) . . 10
5 R
lrp
=R
e
ratio evaluated throughout a hyperbolic tiling (8 sides core cell) . . 10
6 A modied four-point grid . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11
7 R
lrp
=R
e
ratio evaluated throughout the modied four-point lattice . . . . . 12
8 R
lrp
=R
e
ratio evaluated throughout a 1000 nodes preferential attachment graph 13
9 Zigzagging needed to construct the null space eigenvectors of D . . . . . . . 17
10 Eigenvalues of the matrices D, T , and D +T = L
G
for a Euclidean tiling
(hexagons) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17
11 Eigenvalues of the matrices D, T , and D +T = L
G
for a hyperbolic tiling
(heptagons) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18
12 Eigenvalues of the matrices D, T , and D +T = L
G
for a hyperbolic tiling
(decagons) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18
13 R
lrp
(1;j) evaluated throughout a Euclidean tiling (6-side core cell) . . . . . . 19
14 R
lrp
(1;j) evaluated throughout a hyperbolic tiling (7-side core cell) . . . . . 19
15 Number of visits per node in a 6-side polygon data network . . . . . . . . . . 21
16 Number of visits per node in a 7-side polygon data network . . . . . . . . . . 22
17 Number of visits per node in a 8-side polygon data network . . . . . . . . . . 22
18 Flat areas of
lrp
=
e
curves for the P -graph of the IEEE300 bus system
indicating negative curvature along the corresponding paths (Note: for clarity
reasons only 6 curves are shown out of 300.) . . . . . . . . . . . . . . . . . . 29
19 Curvature Centrality of critical nodes, (a), for the P -graph of the IEEE300
bus system (Denition 16). The height of each bar represents the number of
times each critical bus appears within all critical lines. . . . . . . . . . . . . 30
20 Curvature values of critical edges . . . . . . . . . . . . . . . . . . . . . . . . 33
21 Curvature ealues of critical edges (No negative curvature values are present.) 37
22 Utilization Factor (
AC
) for the AC model of the IEEE300 bus system with
ORC-OPF implementations and a 700 MW active power constraint . . . . . 42
23 ORC Values Dierence between tripped and no tripped lines (Tripping of line
250 with 700 MW branch limit) . . . . . . . . . . . . . . . . . . . . . . . . . 48
24 ORC values dierence between tripped and no tripped lines (zoom in of lines
150 to 409 of Fig. 23, tripping of line 250 with 700 MW branch limit) . . . . 49
25 ORC values dierence between tripped and no tripped lines (zoom in of lines
1 to 100 of Fig. 23, tripping of line 250 with 700 MW branch limit) . . . . . 49
26 IEEE300 Bus connectivity highlighting aected edges (in red) after the trip
of line 250 (700MW branch limit). The width of each red edge represents the
degree to which each edge has been aected. . . . . . . . . . . . . . . . . . . 50
27 ORC values dierence between tripped and no tripped lines (Tripping of line
322 with 900 MW branch limit.) . . . . . . . . . . . . . . . . . . . . . . . . . 50
28 ORC values dierence between tripped and no tripped lines (zoom in of lines
317 and 320 of Fi. 27 after tripping of line 322 with 900 MW branch limit) . 51
vii
29 Branch utilization factor (Line 322 utilization factor has been highlighted .) . 51
30 IEEE300 Bus connectivity highlighting aected edges (in red) after the trip
of line 322 (900MW branch limit). The width of each red edge represents the
degree to which each edge has been aected. . . . . . . . . . . . . . . . . . . 52
31 Curvature centrality ((a)) dierence between tripped and no tripped lines
(Tripping of line 220 with 700MW branch limit.) . . . . . . . . . . . . . . . . 55
32 Curvature entrality ((a)) dierence between tripped and no tripped lines
(zoom in of nodes 240 to 300 of Fig. 31) . . . . . . . . . . . . . . . . . . . . 56
33 Curvature centrality ((a)) dierence between tripped and no tripped lines
(zoom in of nodes 1 to 150 of Fig. 31) . . . . . . . . . . . . . . . . . . . . . . 56
34 IEEE300 Bus connectivity highlighting aected edges (in red) after the trip
of line 220 (700MW branch limit). The width of each red edge represents the
degree to which each edge has been aected. . . . . . . . . . . . . . . . . . . 57
35 Curvature centrality ((a)) dierence between tripped and no tripped lines
(Tripping of line 266 with 900MW branch limit.) . . . . . . . . . . . . . . . . 57
36 Curvature centrality ((a)) dierence between tripped and no tripped lines
(zoom in of nodes 1 to 150 of Fig. 35) . . . . . . . . . . . . . . . . . . . . . . 58
37 Curvature centrality ((a)) dierence between tripped and no tripped lines
(zoom in of nodes 240 to 300 of Fig. 35) . . . . . . . . . . . . . . . . . . . . 58
38 Branch utilization factor (Line 266 AC utilization factor has been highlighted.) 59
39 IEEE300 Bus connectivity highlighting aected edges (in red) after the trip
of line 266 (900MW branch limit). The width of each red edge represents the
degree to which each edge has been aected. . . . . . . . . . . . . . . . . . . 59
40 Nodal entropy as per Def. 20 (Line 220, connecting buses 119 and 121, has
tripped.) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62
41 Nodal entropy as per Def. 20 (Line 266, connecting buses 169 and 210, has
tripped.) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 63
viii
Abstract
This work is divided in three main parts. The rst one develops the asymptotic relation
between the resistance of the least resistive pathR
lrp
(i;j) and the eective resistanceR
e
(i;j)
between two extreme points i;j of a resistive lattice, as an electrical engineering mean to
check negative curvature of a graph. More specically, the behavior of R
lrp
(i;j)=R
e
(i;j),
understood asymptotically in the sense of a growing size network, is able to discriminate a
Euclidean grid versus a hyperbolic grid. It is shown that the ratio asymptotically saturates
in case of a hyperbolic grid, while it grows without bound in the Euclidean case. In the
case of a Gromov hyperbolic grid quasi-isometric to a tree, the asymptotic saturation of
R
lrp
(i;j)=R
e
(i;j) is easy to prove, but the major contribution of this rst part is to show
that the same result holds under the by far less trivial situation of a hyperbolic grid not
quasi-isometric to a tree. Applications include the role of curvature in congestion control
problems understood in the broad sense as they apply to power grid and data network.
The second part develops a congestion management method for the power grid utilizing
the notion of curvature. It initially uses the curvature concept to detect areas prompt to
congestion (negative curvature areas) and it subsequently applies load balancing techniques
(through FACTS devices) and load (storage devices) deployment to maximize curvature (grid
decongestion) and cost-eectively minimize the generated energy throughout the grid, while
at the same time guaranteeing stability under phase angle and voltage constraints. Two
dierent curvature denitions are compared: Ollivier-Ricci Curvature (ORC) and Eective
Resistance Curvature (ERC), and shown to be locally, globally, resp., relevant. An entropy
concept suitable to power grid is also introduced as a new measure to analyze grid congestion.
The third part utilizes the ORC and the ERC to analyze the new power distribution
of the grid after tripping of main power grid lines. The main idea is to construct a moni-
toring/systematic method that identies the most vulnerable lines and corridors within the
power grid that will be at risk after the tripping of a line. Overhead lines are subject to
thermal limitations; whenever these limits are overpassed, or whenever a fault is detected,
ix
the risk of overheating is imminent and protection switches are activated to protect the rest
of the grid. Under this scenario, the power is suddenly forwarded to dierent grid areas that
might suer a fault/overheating as a consequence of this forwarding; thus, the detection and
identication of these critical areas, lines and corridors is crucial to avoid a cascade eect.
x
Introduction
Resistive networks are ubiquitous in many problems: random walks, power grid (the so-
called DC power
ow equations), data networks, to name but a few. It is expected that the
electrical properties of the network depend strongly on the underlying graph topology. The
most signicant electrical property is the eective resistance R
e
(i;j) between two vertices,
dened as the resistance \seen" between vertexi and vertexj. For example, a random walk
is recurrent if the asymptotic eective resistance is \large", whereas it is transient if the
eective resistance is \small" [1]. Eective resistance has also proved useful in large scale
power
ow in the smart grid [2].
Many real life networks have been found to have properties reminiscent of the Gromov
hyperbolic property, with the diculty that real life networks no matter how huge have nite
size, whereas the Gromov property is an idealization that makes sense only for innite graphs.
A folklore interpretation of a Gromov hyperbolic network is that it is a \fattened" tree.
Certainly, a fattened tree is Gromov hyperbolic, but the converse does not hold true. There
are Gromov hyperbolic graphs that are not fattened trees|case in point, the hyperbolic
tessellation of the Poincar e disk. A simple way to prove that the hyperbolic tessellation of
the Poincar e disk and a fattened tree are not quasi-isometric is that the former has the unit
circle S
1
as Gromov boundary @G, while the latter has a Cantor set structured Gromov
boundary. In either case,j@Gj 2 and such networks are prone to congestion [11]; on the
other hand, networks such thatj@Gj = 1 are not prone to congestion. This exemplies the
subtlety of the problem.
This often neglected Gromov versus fattened tree discrepancy calls for a re-examination
of the \fattened tree" results in the realm of Gromov hyperbolic networks that are not quasi-
isometric to trees. This is precisely one of the purpose of this work. We take a hyperbolic
\tiling" as a model of a 2-dimensional hyperbolic space and prove the asymptotic saturation
of R
lrp
(i;j)=R
e
(i;j).
Thus, the bounded asymptotic evolution of the R
lrp
(i;j)=R
e
(i;j) ratio between two
1
points i;j of the resistive grid associated with a graph might be eciently used to detect
negatively curved behavior in the graph. From there on, congestion control applications go
along two directions: (i) Either leave the network \as is," and adapt the protocol according
to the network curvature; this approach is better suited to wireless networks [3]; or else (ii)
given that such protocol as the Shortest Path Routing wireline protocol is enforced and does
not work well with negatively curved network, modify the link weight so as to make the
curvature non negative; this is the \curvature control" approach [4].
This work focuses on the \curvature control" approach, and it develops a novel method to
overcome the negative curvature/congestion feature of networks. Once the ground work for
negative curvature alleviation is set up, this work further investigates new methods to reduce
transmission costs within the power grid. Indeed, one of the premises of the \smart grid"
is to allow consumers to purchase electricity at the cheapest price, with the drawback that
this creates large transfer of power across the grid with the potential to overload some lines;
thus, an ecient transmission cost reduction procedure will need to contemplate congestion
alleviation. Such a method is proposed in the second part of this work.
After the tripping of a main power grid line occurs, the active power is re-directed through
out the grid; this automatic reaction of the grid is performed by the protection switches to
avoid thermal rating violations. As a consequence, an important amount of active power
could reach vulnerable/critical areas of the grid that could also trip due to thermal limits
violations. This is basically known as a cascade eect. The third part of this work proposes
a curvature approach to detect those vulnerable areas/corridors where overhead conductors
thermal limits could be violated as a consequence of an active power forwarding. The
proposed curvature analysis is also combined with a nodal entropy feature which visually
helps to conrm where the line tripping is ocurring.
This work is organized along the following sections: Section 1|Essential Background|
includes a collection of well suited denitions and results that have been strategically linked
to gradually construct the required background. It includes spectral graph theory together
2
with adjacency matrices and Laplacians.
Section 1.4|Eective Resistance of Prototypical Euclidean and Hyperbolic Grids|builds
on the previous sections to introduce several test cases of lattices with either vanishing or
negative curvature. Section 1.5|Hilbert Space Insight|already highlights the vanishing
versus negative curvature issue in the context of operator perturbation theory. Section 1.6|
Main result|presents a sketch of the proof of the bounded behavior of
R
lrp
(i;j)=R
e
(i;j) in hyperbolic grids. Section 1.8 |Data Networks|presents the rst set
of simulations of the applications that are targeted.
In order to construct a cost reduction procedure within the power grid, the foundational
material for power grid analysis needs to be reviewed (Sec. 2). Therefore, Section 2.1 crafts
an overview of various simplications of the power
ow equations, the \DC power
ow
equations," and derive the P , Q and S-graphs, all of which are resistive network models
of the various power
ows. Negative curvature of such networks then immediately leads to
an overall warning of potential \congestion" or \line overloading" in the underlying grid, as
demonstrated in Sec. 2.2.
The major results come in Section 2.4, where remedies to such congestion situations are
developed in a cost conscientious manner by having the cost of the generation taken into
consideration. Specically, we begin with a curvature-based congestion reduction proce-
dure. Then we develop a curvature-driven algorithm that minimizes a quadratic generation
cost functional. It basically injects loads in very strategically chosen points within the grid
revealed by the curvature procedure, resulting in a drastic reduction in the overall cost
(quadratic function on the generator active power). Next, a combined curvature smooth-
ing/cost minimization optimization is proposed, for both the DC and the AC models, with
the surprising result that the full AC model allows for less con
ict between the curvature
and the cost optimization objectives.
Sec. 3 utilizes the ORC and ERC to detect those areas under risk after a power grid line
is suddenly tripped and the active power is rerouted by the protecton switches of the grid.
3
These vulnerable areas and corridors are ought to be identied in advance to protect them
and avoid the well known cascade eect.
Finally, the last part of this work is summarized in Sec. 5, where conclusions and several
possible future extensions of this work are presented.
4
1 Essential background
1.1 Graph Theory
Let G = (V;E) be a simple graph, where V = (v
1
;v
2
;:::;v
n
) denotes the vertex set and
E = (e
1
;e
2
;:::;e
m
) the edge set. n is called the order of the graph andm is called the size of
the graph. By simple graph, we mean an unweighted, undirected graph containing no loops
nor multiple edges. Each edge e
i
is dened by a pair of verticesfv
j
;v
k
g, which means that
an edge requires exactly two vertices to be dened.
Denition 1: A weighted graph G = (V;E;w) is a graph equipped with a weighting function
w : E!R
+
that associates a weight with every edge in the graph. w might be thought as
the distance between the pair of nodes that dene the edge.
Denition 2: The degree of a vertex is dened as the number of edges connected to the vertex
v and is written as degree(v) =k(v).
Denition 3: The adjacency matrix A is an nn symmetric matrix where a non-diagonal
entry a
ij
= 1 if vertex i and vertex j are connected with an edge and a
ij
= 0 otherwise; the
diagonal entries are set to zero.
Denition 4: The topological graph Laplacian operator L
G
is the symmetric positive semidef-
inite matrix dened as L
G
= diag(degree(v
i
))A.
1.2 Resistive Networks
A resistive network [5] is a weighted graph G = (V;E;w) where w(e) is the resistance of
the edge e.
Denition 5: The eective resistance [6] between two nodes of a weighted graph, R
e
(v
i
;v
j
)
is the electrical resistance seen between the two nodes.
5
Denition 6: The Laplacian of the resistive networks is ann matrix where the o-diagonal
element L
i;j
is dened as1=w(ij) and the diagonal element is dened as L
ii
=
P
j6=i
L
ij
.
In case of a homogeneous resistive networks with edge resistance R, its Laplacian is (1=R)
times its topological graph Laplacian.
Let I = (I
1
;I
2
;:::;I
n
) be the column vector of currents injected at the nodes and
E = (E
1
;E
2
;:::;E
n
) be the column vector of node potentials relative to a common ground.
Recalling Kirchho's laws, it is readily found that
I =
1
R
L
G
E; (1)
where it is assumed that all edges have the same resistance R. It is now straightforward to
conclude that the eective resistance from node i to j within our G graph is given by:
R
e
(i;j) = R(L
ii
1
)
jj
; (2)
where L
ii
denotes the matrix L
G
after removal of column i and row i.
An alternative eective resistance formulation that will be also used through out this
work, is the so called \the two-point resistance" (see [3]). The beauty of this formulation
is that it expresses the resistance between two points in a network as a function of the
eigenvalues
k6=1
and eigenvectors v
k
= (v
k1
;v
k2
;v
k3
;:::;v
kn
) of the graph Laplacian of the
network:
R
e
(i;j) =
n
X
k=2
1

k
jv
ki
v
k;j
j
2
: (3)
Denition 7: The resistance of the least resistive path,R
lrp
(v
i
;v
j
) between vertex nodes (v
i
;v
j
)
in a weighted graph is dened as the smallest sum of the resistances(weights) encountered
along all possible paths.
R
lrp
can be computed eciently using the Bellman-Ford or the Dijkstra algorithm.
6
1.3 Curvature of Graphs
Denition 8 (Curvature of a Surface): To develop the notion of curvature of a surface at a
point p, we draw the tangent plane and the normal to the surface at p. The hight of the
surface above the tangent plane is approximated by a quadratic function, with coecient
matrix with eigenvaluesk
1
,k
2
, the so-called principal curvatures, and their product is named
as the Gauss Curvature K =k
1
k
2
. For a rigorous analysis of curvature of a surface see [7].
Denition 9 (Negatively Curved Surface): A surface X is said to be negatively curved if
K(v)< 0 for every point v2X.
Now that the basic notion of a surface curvature has been dened, we dene the curvature
of a graph.
Denition 10 (Curvature of a Graph): Given a weighted triangulation (T;) of a compact
surface X [8], where the vertex set is given by V = (v
1
;v
2
;v
3
;:::;v
n
), the edge set dened by
E = (e
1
;e
2
;e
3
;:::;e
m
), and the weighting function dened by w :E!R
+
, the curvature K
i
at the vertex v
i
is given by:
K
i
= 2a
i
; (4)
where a
i
is called the cone angle at vertex v
i
, and is given by the sum of all inner angles
having v
i
as a vertex; precisely, the angle at the vertex i of the triangle jik is given by the
cosine law, a
i
= arccos((w(i;j)
2
+w(i;k)
2
w(j;k)
2
)=2w(i;j)w(i;k)).
Denition 11 (Negatively curved graph): A planar graph X that is a triangulation (of some
surface) is said to be negatively curved if K
i
< 0,8i2V (G).
Denition 12: The genus g of as compact surface X is the number of handles that need to be
attached to the sphere S
2
to obtain a homeomorph of X.
Theorem 1 (Gauss-Bonnet theorem for compact surface): Given a compact triangulated sur-
face X with Euler characteristic = 2 2g, we have
7
X
i
K
i
= 2:
1.4 Eective Resistance Curvature of Prototypical Euclidean and
Hyperbolic Graphs
The objective of this section is to evaluate the R
lrp
=R
e
quotient in standard regular
lattice test cases, and compare the curvature outcome (negative if the quotient is bounded
through the lattice) with the degree of the nodes.
1.4.1 Euclidean Graphs
We can clearly see from Fig. 1 and Fig. 2 that, for a 2-dimensional, 6-side regular polygon
lattice, the R
lrp
=R
e
ratio it is not bounded from above, indicating non-negative nature of
the lattice.
Figure 1: Euclidean tiling (6 sides core cell)
8
10
12
14
16
18
20
Rlrp/Reff
R lrp/R eff ratio for different path lenghts
8
10
12
14
Rlrp/Reff
Rlrp/Reff ratio for different path lenghts
0 10 0 200 30 0 400 50 0
0
2
4
6
8
Rlrp/Reff
N umber of nodes in the lattice
0 200 400 600
0
2
4
6
Rlrp/Reff
Node number ~ distance
R lrp/R eff ratio for different path lenghts
Rlrp/Reff ratio for different path lenghts
50 0 600 70 0 800 900
N umber of nodes in the lattice
800 1000 1200 1400
Node number ~ distance
Rlrp /Reff ratio for different paths lengths
Number of nodes in the lattice
Rlrp /
Reff
Figure 2: R
lrp
=R
e
ratio evaluated throughout an Euclidean tiling (6 sides core cell)
1.4.2 Hyperbolic Graphs
From Fig. 3 and Fig. 4, we observe that in a 2-dimensional 7-side regular polygon lattice
the curvature is negative (the R
lrp
=R
e
is bounded from above).

Figure 3: Hyperbolic tiling (7 sides core cell)
Following this line of thoughts, we expect that a tiling with an 8-side building core cell will
also exhibit a bounded evolution of R
lrp
=R
e
. As show in Fig. 5, the simulation experiment
9
Rlrp / Reff ratio for different path lengths
Node number ~ distance
Rlrp /
Reff
Number of nodes in the lattice Number of nodes in the lattice
Figure 4: R
lrp
=R
e
ratio evaluated throughout a hyperbolic tiling (7 sides core cell)
matches the anticipated pattern.
Rlrp / Reff ratio for different path lengths
Rlrp /
Reff
Node number ~ distance
Number of nodes in the lattice
Figure 5: R
lrp
=R
e
ratio evaluated throughout a hyperbolic tiling (8 sides core cell)
10
1.4.3 A Counterexample: Modied Four-Point Graph
If the genus of a compact surface g > 2, the Gauss-Bonnet theorem indicates negative
curvature at some points; less trivially, it can be shown that the surface can be endowed with
a hyperbolic metric and the graph written on the surface would be negatively curved [9].
To illustrate the subtleties involved, consider the graph shown in Fig. 6. It is nonplanar,
with high genus (g> 2), giving the (wrong!) impression of being negatively curved.
Figure 6: A modied four-point grid
If we deploy theR
lrp
=R
e
calculation in this case, we observe from Fig. 7 that the electrical
measure does not show sign of a bounded growing behavior. This fact apparently contradicts
our major claim.
Nevertheless, by inspection of the structure of the graph, we observe that the graph is
pretty meshed; in fact, it is quasi-isometric to the standard square lattice which has vanishing
curvature. The apparent contradiction between g > 2 and 0 curvature can be explained by
the fact that this lattice is not compact.
1.4.4 Growth-Preferential Attachment Graph
This section deploys theR
lrp
=R
e
tool on the so-called \growth-preferential attachment"
graph, a popular model of the build up of the Internet, which follows the experimental
power law distribution. Under proper combination of parameters, this type of networks are
known to have \scaled Gromov-hyperbolic" properties [10]. Fig. 8 clearly shows a bounded
11
Figure 7: R
lrp
=R
e
ratio evaluated throughout the modied four-point lattice
behavior of the R
lrp
=R
e
, indicating that the negative curvature criterion applies to more
general graphs than the test ones.
1.5 Hilbert Space Insight
This paper deals with the eective resistance between two points, extreme in the sense of
(i;j) = arg max
k;`
R
lrp
(k;`), under a spiraling expansion of the network. This is quite dier-
ent from a priori assuming that the network is innite and looking at
lim
R
lrp
(i;j)!1
R
lrp
(i;j)=R
e
(i;j). Nevertheless, the spiraling expansion of the network builds
a Hilbert space operator, which provides some useful asymptotic eigenstructure results for
the nite operator.
The graph Laplacian obtained by assigning a label increasing along a counterclockwise
diverging spiral of vertices, in either the Euclidean or the hyperbolic case, can be decomposed
12
Figure 8: R
lrp
=R
e
ratio evaluated throughout a 1000 nodes preferential attachment graph
as a sum of two operators,
L
G
=T +D;
where T stands for the Toeplitz part, and D is taken as a perturbation departure from the
Toeplitz matrix. The Toeplitz part has symbol () = deg + exp(j) + exp(j). The
(essential) spectrum of the Toeplitz part is well known to be a compact interval [;] :=
([0; 2)) of the real line, and the surprising numerically observed feature is that the spectrum
of L
G
and the one of T coincide as the size of the network increases. This appears to be
the well known perturbation result that the essential spectrum of an operator is invariant
under a compact perturbation. The problem is that, even though D has innitely many
eigenvalues at 0, it is not a compact operator; nevertheless, as we shall prove, some weak
form of the classical spectral perturbation result still holds. This is the main mathematical
contribution of the paper.
The structure of the departure from Toeplitz matrix D is easily seen to be as follows:
13
The diagonal entries are all vanishing, and the upper triangular part is described row-wise
as follows:
1. The rst row consists of two 0's followed by a line of 1's, call it D
1;J
1
, where J
1
=
(i
1;1
;i
1;2
=i
1;1
+ 1;:::;i
1;jJ
1
j)
is the ordered set of consecutive column indices of entries
equal to 1.
2. The second row is zero everywhere, except for those entries indexed byJ
2
= (j
2;1
;:::;j
2;jJ
2
j
),
which equal 1; furthermore, maxJ
1
= minJ
2
.
3. The third row is zero everywhere, except for those entries indexed byJ
3
= (j
3;1
;:::;j
3;jJ
3
j
),
which equal 1; furthermore, maxJ
2
= minJ
3
.
4. etc.
Then the lower triangular part is constructed so as to make D symmetric. Table 1 makes
the J intervals more explicit.
Table 1: Column indices of nonvanishing entries of rows 1 to 8 of D
row# degree = 6 degree = 7 degree = 8
J
1
(3;:::; 7) (3;:::; 8) (3;:::; 9)
J
2
(7;:::; 10) (8;:::; 12) (9;:::; 14)
J
3
(10;:::; 12) (12;:::; 15) (14;:::; 18)
J
4
(12;:::; 14) (15;:::; 18) (18;:::; 22)
J
5
(14;:::; 16) (18;:::; 21) (22;:::; 26)
J
6
(16;:::; 18) (21;:::; 24) (26;:::; 30)
J
7
(18;:::; 19) (24;:::; 27) (30;:::; 34)
J
8
(19;:::; 21) (27;:::; 29) (34;:::; 38)
Denition 1: A Hilbert space operator is compact if it transforms any bounded set into a
compact set.
To prove that the D operator is not compact, we construct a bounded sequencefx
(m)
:
m = 1; 2;:::g,kx
(m)
k = 1, such that the image sequencefDx
(m)
: m = 1; 2;:::g has no
14
converging subsequence. Letfe
k
:k = 1; 2;:::g be the natural basis of the Hilbert space `
2
.
By the Zermelo Axiom of Choice, pick j
m
2J
m
, m = 1; 2;:::. Dene x
(m)
=e
jm
. Then it is
easily seen that Dx
(m)
=e
m
. Note that for any image subsequence such that m
i+1
>m
i
we
have lim sup
i;j!1

Dx
(m
i
)
Dx
(m
j
)

`
2
= 2, so that no such image subsequence is Cauchy.
Since a Hilbert space is complete, no such image subsequence converges.
Theorem 2: The operator D has innitely many eigenvalues at 0.
Before going through the proof, we need some denitions. Ro is a mapping that picks
a column index of D and maps it to the set of row indexes below the diagonal such that
the corresponding entries in D are nonvanishing. Dually, Co picks a row index of D and
maps it to he corresponding column indexes above the diagonal and such that the entries
are nonvanishing.
Proof. (Sketch) We recursively construct a vector x2 `
2
such that Dx = 0 and show that
countably innitely many such vectors exist. Such eigenvectors are indexed by m = 1; 2;:::
in the sense that their rst nonvanishing components are in J
m
. We start with m = 1
and construct the corresponding eigenvector. Pick the rst two nonvanishing components
x
n
1
= x
n
2
= 1 with n
1
;n
2
2 J
1
. Set x = (0;:::; 0;x
n
1
; 0;:::; 0;x
n
2
; 0;:::). This yields
(Dx)
1
= 0. But this yields components indexed by J
Ro(n
1
)
and J
Ro(n
2
)
that are nonzeros,
and that have to be compensated. This is done via components of magnitude 1=2 picked in
J
CoRo(n
1
)
and J
CoRo(n
2
)
. Update the vector x with those new components, but this yields
components further down, in J
RoCoRo(n
1
)
and J
RoCoRo(n
2
)
that are nonzero and need to
be compensated for. The recursive, zigzagging between Ro and Co mapping should now be
obvious. This takes care of the rst null space eigenvector. For the second such eigenvector
start with the rst nonvanishing components in J
2
, etc.
The eigenvalues of the principal top left hand corner submatrix ofD are shown in the red
curve of Figs. 11, 12. Clearly, 0 is a multiple eigenvalue|an innite dimensional eigenvalue
15
for the Hilbert space operator. Since the Hilbert space operator is not compact, it must have
essential eigenvalues away from 0. From the spectrum of the submatrix, we see at least one
area of the red curve tangent to the x-axis at a y value of , meaning an eigenvalue of the
Hilbert space operator accumulating at .
Since the operator D is noncompact, one cannot expect the full essential spectrum of T
and that of T +D to coincide. However, looking at Fig. 12 with the black and blue curves
coinciding over a subset of indices, it appears that a subset of the essential spectrum ofD+T
coincide with an essential subspectrum of D. Precisely,
Theorem 3: Let E
0
(D) be the innite-dimensional eigenspace of D corresponding to the es-
sential spectrumf0g. Then the subsets of the essential spectra ofD andD +T coincide over
E
0
(D).
Proof. Take a sequence x
n
2 E
0
(D) such thatkx
n
k = 1 and x
n
* 0. Clearly Dx
n
! 0.
Take a E
0
(D) element in the eigenspace of T corresponding to in the essential spectrum.
Clearly (TI)x
n
! 0. It follows that (T +D)x
n
= (TI)x
n
+Dx
n
! (T)x
n
! 0.
Hence is in the essential spectrum of T +D.
An early illustration of this theorem is given in Fig. 11, showing that over the eigenspace of
the essential spectrumf0g ofD, the spectra ofT andT +D nearly coincide. The subspectra
of T and T +D do not quite coincide, because E
0
(D) does not completely encompass the
eigenspace of the essential spectrum ofT . Nevertheless, as the degree grows,E
0
(D) expands,
so that over the accumulation point 0 of the spectrum of D, the spectra of T and T +D
coincide much better, as seen from Fig. 12.
1.6 Main Result
As it has been brie
y mentioned in Sec. 1.4, a convenient numbering assignment for this
type of lattices is the one that assigns the node number 1 to the center of the network and
then consecutively labels the nodes in a counterclockwise spiraling manner. This rst entails
16
0 0 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0.....
0 0 0 0 0 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0.....
1 0 0 0 0 0 0 0 1 1 1 1 0 0 0 0 0 0 0 0 0.....
1 0 0 0 0 0 0 0 0 0 0 1 1 1 1 0 0 0 0 0 0 0.....
1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 1 1 0.....
1 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 1 1 0 0 0 0 ....
0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 1 1 0.....
0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1.....
0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0.....
Ro
Co
R1
Figure 9: Zigzagging needed to construct the null space eigenvectors of D
Figure 10: Eigenvalues of the matrices D, T , and D + T = L
G
for a Euclidean tiling
(hexagons)
a polygon of lengthjJ
1
j with R
lrp
(1;J
1
) = 1. This is followed by a polygonal arc of length
jJ
2
j, another polynomial arc of lengthjJ
3
j, etc., all of them having the propertyR
lrp
(1;J
2
) =
R
lrp
(1;J
3
) = ::: = 2, until a second full counterclockwise circle is made with the polygonal
arcs, ending up with the polygonal arc of lengthjJ
Ro(minJ
2
)
j. This is followed by a polygonal
17
Figure 11: Eigenvalues of the matrices D, T , and D + T = L
G
for a hyperbolic tiling
(heptagons)
Figure 12: Eigenvalues of the matrices D, T , and D + T = L
G
for a hyperbolic tiling
(decagons)
18
arc of lengthjJ
Ro(minJ
2
)+1
j,..., ending with a polygonal arc of lengthjJ
Ro(min Ro(minJ
2
+1))
j, etc.
Figure 13: R
lrp
(1;j) evaluated throughout a Euclidean tiling (6-side core cell)
Figure 14: R
lrp
(1;j) evaluated throughout a hyperbolic tiling (7-side core cell)
19
Fig. 13- 14 illustrate the process. They plot R
lrp
(1;j)) against the ascending numbering
of nodesj in the Eucliean and hyperbolic case, respectively. The stair atR
lrp
= 1 has depth
jJ
1
j. Then the stair at R
lrp
(1;J
1
) = 2 has \depth"jJ
1
j +jJ
2
j +::: +jJ
Ro(J
2
)+1
j, etc.
For the 7-side regular lattice of Fig. 14, theR
lrp
(1;j) \stairs" are wider as we depart from
the center; this plays an important role in the calculation of the ratio R
lrp
(1;j)=R
e
(1;j)
between nodes 1 and j, because for a longer range of j nodes, the resistance ratio will only
change if R
e
(1;j) changes, because R
lrp
(1;j) remains constant in that section.
The crucial dierence between deg = 6 and deg = 7 is that in the former case the depth
of the stair improves slowly, while it increases much more dramatically in the latter.
1.7 Curvature Detection
Contrasting Fig. 10 with Figs. 11 and 12, the dierence between 0 curvature and negative
curvature is striking. We now proceed to show another facet of this striking dierence.
Theorem 4: For a regular 2-dimensional 6-side polygon (Euclidean) lattice limited to j ver-
tices labeled from 1, the center, along a counterclockwise diverging spiral, toj, the end vertex,
in a process formalized in the Appendix,
lim
j!1
R
lrp
(1;j)
R
e
(1;j)
=1:
For the same lattice, but built on k-side regular polygons with k 7 making it negatively
curved,
lim
j!1
R
lrp
(1;j)
R
e
(1;j)
<1:
1.8 Data Network
This section analyzes the distribution of packages in data networks of varying topology.
The path taken by the packages in every network simulation is dictated by the Bellman-Ford
20
Visits per node
Node number
Number
of visits
on each
node
Figure 15: Number of visits per node in a 6-side polygon data network
protocol, which basically delivers packages between nodes using minimum length paths. Each
network is ideally supposed to have an innite buer in each node and there are no channel
capacity restrictions. It is also assumed that at each time and at every node a package should
be delivered to another node with uniform probability across all nodes. In other words, every
node acts as a source and as a sink in any time slot. The connectivity of the proposed data
networks are the 2-dimensional 6, 7, and 8-side regular polygon lattices, where each network
has a total of 49 nodes.
The congestion metric in Figs. 15- 17 is the number of packets that has transited through
a node during the run of the experiment. It is observed from Figs. 15- 17 that as we proceed
from 0 curvature to more negative curvature, the congestion at vertex #1 worsens, a fact
consistent with [11].
21
Visits per node
Node number
Figure 16: Number of visits per node in a 7-side polygon data network
Visits per Node
Node number
Figure 17: Number of visits per node in a 8-side polygon data network
22
2 Curvature and the Power Grid
Congestion Management of the Power Grid has been extensively researched within the
last few years; a comprehensive literature review can be found in [33], [34], [35], [36],
[37], [38] and [39]. Within all the algorithms and strategies developed throughout the years
we would mention the followings: Generator Rescheduling (GR), load shedding, Distributed
Generation (DG), Optimal Power Flow (OPF), Flexible Alternating Current Transmission
System (FACTS) devices, implementable via Articial Bee Colony algorithm (ABC), Genetic
Algorithms (GA) and Strength Pareto Evolutionary Algorithm (SPEA), just to name a few.
This section introduces a novel congestion management (CM) technique that basically
grounds its roots in the notion of curvature. The curvature is essentially utilized to detect
areas prompt to congestion (negative curvature areas); once these areas are identied via
the curvature analysis, FACTS devices are deployed to maximize the curvature, but more
practically storage devices are deployed to buses/links to reduce the
ow of active power.
Both deployments are dictated by the curvature analysis. To cost-eectively reduce the
generation of energy throughout the entire grid, an Optimal Power Flow (OPF) algorithm
is implemented to minimize the cost of energy production while maintaining stability under
the usual phase angle and voltage constraints. Therefore, two objectives are simultaneously
reached: grid decongestion and energy generation cost reduction.
Two curvature notions are compared: the local Ollivier-Ricci Curvature (denition in
Sec. 2.2) and the global Eective Resistance Curvature( [42]). Both of them are eectively
used to detect all the critical (negatively curved) edges prompt to transmit most of the active
power within the grid. The Ollivier-Ricci Curvature (ORC) concept traces its origins back
to the earth mover's distance idea developed in the Napoleonic era with the objective of
eciently move earth from one point to the other to level o the landscape. The Earth
Mover's eort is mathematically quantied as the Wasserstein 1-metric distance, and it is
also known as the transportation cost (Gaspard Monge, 1781).
Finally, towards the end of this work, a Markov chain entropy is introduced, compared
23
with the curvature, and applied to the power grid aiming at providing additional global
information on the behavior of the grid.
2.1 Resistive Network Models of Power Flows
Given two buses k and m specied by their voltage magnitude and phase angle pairs
(V
k
;
k
) and (V
m
;
m
), resp., connected by a transmission line with admittanceY
km
=G
km

jB
km
, the power
ow equations are well known as
P
km
= G
km
V
2
k
+B
km
V
k
V
m
sin(
k

m
)
G
km
V
k
V
m
cos(
k

m
);
Q
km
= B
km
V
2
k
B
km
V
k
V
m
cos(
k

m
)
G
km
V
k
V
m
sin(
k

m
);
where P
km
and Q
km
are the active and reactive power, resp.,
owing from bus k to bus m.
Under the standard approximations of a nearly lossless lines (G
km
0) with small phase
angle dierences (
k

m
) and small voltage magnitude dierences (V
k
V
m
), the power

ow equations are simplied to be
P
km
=B
km
V
k
V
m
(
k

m
); Q
km
=V
k
B
km
(V
k
V
m
):
Hence, P
km
can be viewed as the current
owing through a resistor
km
= 1=B
km
V
k
V
m
driven by a voltage potential dierence
k

m
. Active powers injected at some buses are
then modeled as currents injected at the corresponding nodes of the resistive network. Let
us call this resistive network the P -graph.
Similarly, but subject to a discrepancy to be explained soon, Q
km
can be viewed as the
current
owing through a resistor
km
= 1=B
km
V
k
driven by a voltage potential dierence
V
k
V
m
. The discrepancy relative to the active power case is that the resistor is directional,
24

km
6=
mk
. We refer to this directed resistive network as the Q-digraph. Ricci curvature
concepts for such digraphs can be developed in the context of Finsler geometry.
While the P -graph and the Q-digraph allow for a quick snapshot of active and reactive
power load, respectively, the apparent power carried by the line is
p
P
2
km
+Q
2
km
, mandating
some way to combine the two power
ows [2]. Here we combine the two power
ows in the
way their
uctuations, resulting from such renewables as wind farms [47], could overload
the lines. Dene

k
,

m
,

V
k
,

V
m
as the average values and
~

k
(t) :=
k
(t)

k
,
~

m
(t) :=

m
(t)

m
,
~
V
k
(t) := V
k
(t)

V
k
,
~
V
m
(t) := V
m
(t)

V
m
as the
uctuating values. Under
the assumption that P
km
(t) depends more on the phase angles rather than the bus voltages
and that the transmission lines are nearly lossless, the following rst order approximation is
easily obtained [2]:
~
P
km
(t) =B
km

V
k

V
m
cos(

k

m
)

~

k
(t)
~

m
(t)

: (5)
Regarding the reactive power
uctuation, on the other hand, it is assumed that Q
km
(t)
depends more on the bus voltages than the phase angles. Further assuming that the line is
nearly lossless, the rst order approximation
Q
km
(t)

Q
km
=
@Q
km
@V
k

V
k
(t)

V
k

+
@Q
km
@V
m

V
m
(t)

V
m

yields
~
Q
km
(t) =

2B
km

V
k
B
km
V
m
cos(

k

m
)

~
V
k
(t)
B
km

V
k
cos(

k

m
)
~
V
m
(t):
(6)
Here, we apply another approximation, which was veried on the IEEE 300 bus system [2],

V
k

V
m
cos(

k

m
): (7)
25
Then, it follows that
~
Q
km
=B
km

V
k

V
m
cos(

k

m
)

~
V
k
(t)

V
k

~
V
m
(t)

V
m

: (8)
Dene the
uctuating complex power
~
S
km
(t) =
~
P
km
(t)+j
~
Q
km
(t). Combining Eqs. (5)-(8)
yields the approximate
uctuating complex power equations
~
S
km
=B
km

V
k

V
m
cos(

k

m
)

~

k
+j
~
V
k
(t)

V
k

~

m
(t) +j
~
V
m

V
m

| {z }
~
E
km
:
(9)
Clearly the
uctuating complex power can be viewed as a (complex) current
owing through
the resistor

km
= 1=B
km

V
k

V
m
cos(

k

b
) (10)
subject to the complex potential dierence
~
E
km
. We refer to this resistive network asS-graph.
Interestingly, while the resistive network for reactive power is directed, the one derived for
the
uctuating complex power is an undirected graph. This is not surprising in view of the
approximation of Eq.(6) by Eq.(7).
2.2 Curvature and Congestion
This work suggests that one of the main steps towards a succesfull congestion management
method is to be able to eectively detect congestion areas. Having that in mind, this
work proposes a geometric approach to the congestion management problem. Basically, we
propose, as core strategy of the method, the utilization of the curvature notion to detect areas
prompt to congestion (negative curvature areas). The bridge between negative curvature and
congestion in power networks was established in [2] and [42].
As mentioned above, this work will present results utilizing two dierent curvature de-
26
nitions, namely: Ollivier-Ricci Curvature and Eective Resistance Curvature.
2.2.1 Local Curvature Notion: The Ollivier-Ricci Curvature
The Ollivier-Ricci concept has recently been applied to many dierent elds outside
mathematics itself. A clear example is the usage of ORC to dierentiate biological networks
corresponding to cancer cells from normal cells [31], and the detection of changes in brain
structural connectivity in people with ASD (Autism Spectrum Disorders) [32]. It need not
be said that ORC has been long applied in image processing too ([40]); thus, ORC as a rst
step in CM seems natural. The following subsection brie
y denes the ORC for graphs. A
better detailed review of the ORC notion can be found in [26], [27], [28], [29] and [30].
2.2.2 Wasserstein Distance, Earth Mover's Distance (EMD) and ORC
Let H be a discrete metric space equipped with a metric d(:;:), and let c
ij
be the cost
of moving a unit mass from x
i
to x
j
; both x
i
and x
j
belong to H. Denote with p and q
two probability distributions in H. Let
ij
0 be the amount of mass to be transferred
fromx
i
tox
j
. The so-called OMT (Optimal Mass Transportation) is the problem of nding
an optimal transfer of mass fromp toq with minimum cost. This can be formulated as ( [32]):
min

X
i;j
c
i;j

i;j
; (11)
subject to
X
j

i;j
=p
i
; 8i;
X
i

i;j
=q
j
; 8j;

i;i
0; 8i;j;
(12)
27
where i and j are connected via an edge. If the previously formulated problem is solved
with a cost c
i;j
=d(x
i
;x
j
)
r
for any positive integer r, then it is said that the solution of the
optimization problem is the r-Wasserstrein Distance. Moreover, if r = 1, the solution is
called Earth Mover's Distance.
Let now (X;d) be a geodesic metric space equipped with probability measuresfp
x
:x2Xg.
Then the Ollivier-Ricci curvature k(x;y) along the geodesic joining x to y is dened as
W
1
= (1k(x;y))d(x;y); (13)
where W
1
is the EMD distance and d the geodesic distance within the space.
Recall now from Sec. 2.1 that our method will initially calculate a resistive network
model called P -graph; thus, we will have a G = (V;E;w) graph, where V is the set of
nodes/buses, E is the set of lines/edges, and w is the set of resistances/weights. Following
our previous equations, the geodesic distanced(x;y) of the ORC formulation for a graph will
be represented by the minimum number of steps/hops needed to go from x toy. Therefore,
after this recap, a simple and short way to calculate the ORC k(x;y) can be implemented
through a linear programming script, and this is what has been done in the forthcoming
sections.
2.2.3 Global Curvature Notion: Eective Resistance Curvature
We dene the Eective Resistance Curvature of an idealized innite resistive network via
the fraction, already dened in [42],
lim

lrp
(k;m)!1

lrp
(k;m)

e
(k;m)
1: (14)
In the above,
lrp
(k;m) is the resistance of the least resistive path fromk tom in the network,
obtained for example by the Bellman-Ford or the Dijkstra algorithm, and
e
(k;m) is the
eective resistance \seen" at the port km. Precisely, inject a current I at node k and draw
28
0 50 100 150 200 250 300
Number of nodes in the lattice
1
1.5
2
2.5
3
3.5
4
Rlrp/Reff
Rlrp/Reff ratio for different path lenghts
Rlrp/Reff for undirected graphs
Negatively Curved/Congested Lines Indication
Figure 18: Flat areas of
lrp
=
e
curves for theP -graph of the IEEE300 bus system indicating
negative curvature along the corresponding paths (Note: for clarity reasons only 6 curves
are shown out of 300.)
the same current at node m; then,
e
(k;m) := (V
k
V
m
)=I, where V
k
and V
m
represent
voltages induced at nodes k and m, resp.
Denition 13 (Negative Eective Resistance Curvature): An innite network is said to be
negatively curved if (14) is bounded and positively curved otherwise. A nite network is said
to be negatively curved if the fraction (14) is near its lower bound.
This curvature concept is specialized for power
ow problems, although it has some
commonalities with the Gromov [42, 43] and the Ollivier-Ricci [48] concepts. The latter is
a curvature concept, along a path rather than at a vertex, directly related to transport and
hence congestion.
If we construct the P -graph model (see Section II) of the IEEE300 bus network and
compute the various fractions (14) for various buses, we obtain a family of curves as depicted
29
Crtitical Nodes Identification
0 50 100 150 200 250 300
Node # (a)
0
5
10
15
20
25
Most critical nodes/buses
Figure 19: Curvature Centrality of critical nodes, (a), for theP -graph of the IEEE300 bus
system (Denition 16). The height of each bar represents the number of times each critical
bus appears within all critical lines.
in Fig. 18. Recall from [2] and [42] that each curve corresponds to an initial node a and
plots all possibles ratios
lrp
(a;k)=
e
(a;k) versusk6=a. Given a busa, the variousk-buses
are relabeled so that the various ratios
lrp
(a;k)=
e
(a;k) are in increasing order.
Considering
lrp
(a;k)=
e
(a;k) 1, the ratio could reach its lower bound, making the
related curve \
at" with
lrp
(a;k)=
e
(a;k) 1. In this situation, most of the current (in
the resistive model) or power (in the power grid) froma tok will
ow along the least resistive
path, hence overloading the transmission lines along that path. In the Monge-Kantorovich
set up, this means that the transport is along a privileged path; in the Ollivier-Ricci curvature
set up, it means that the curvature is negative.
From Fig. 18, it is clear that the IEEE300 P -graph model has several overloaded lines
corresponding to
at curves. More accurate inspection reveals that there are at least 15 buses
with a
attening behavior along the entire grid, which is further betrayed by the consistent
congestion behavior depicted in Fig. 19 (Sec. 2.3).
30
Conversely, if
lrp
(a;k)=
e
(a;k) is monotone increasing above 1, this implies that there
are many paths of a resistance slightly above
lrp
, and so the current or power will be
distributed along those various paths without overloading some specic ones. In the Monge-
Kantorovich set up, this means that the transport is along many paths; in the Ollivier-Ricci
curvature set up, it means that the curvature is positive.
Going back to Fig. 18, recall that each curve is formed by the points that represent the
values of the ratio
lrp
(a;k)=
e
(a;k) for a xed bus a and varying buses k.
Denition 14 (Critical Buses): A busa is critical if its related
lrp
(a;k)=
e
(a;k) curve tends
to be a
at line for the majority of varying buses k.
Note that various points on a
lrp
(a;k)=
e
(a;k) curve represent ratio values for dierent
paths, named (a; 1); (a; 2); , and so, the curve carries information about many branches,
which form dierent paths. Therefore, the topological information extracted from a
at
curve is directly related to the transmission lines connected to its corresponding bus.
Denition 15 (Critical Lines): A critical line is a transmission line connected to a critical
bus.
The critical lines/buses are responsible for most of the congestion in the grid.
Note that this work will indistinguishably use the critical buses/lines idea for any of the
two curvature notions quoted above; thus, a critical bus/line would be the one associated
with negative curvature connotation regardless of the curvature denition invoked.
2.3 New Curvature Centrality Concept
Dierent graph theoretic centrality measures specialized to the power grid can be found
in the literature [41, 44{46]. The present work, which follows in the footsteps of [2], has
the unique feature that it relates directly to congestion (Def. 14 and Def. 15), rather than
referring to a graph-theoretic feature that can be related to congestion.
31
Denition 16 (Curvature Centrality): The curvature centrality (a) of a critical bus a is the
number ofk-times the ratio
lrp
(a;k)=
e
(a;k) remains \
at" in the (a;k) diagram/curve/ratio.
Therefore, as per Def. 16, every time the ratio
lrp
(a;k)=
e
(a;k) is one (constant) the
curvature centrality stores the bus number a and and keeps a record of the number of times
the bus a produces a constant (equal to one)
lrp
(a;k)=
e
(a;k) ratio for dierent k
0
s. This
record is actually kept for every bus, not only for bus a; thus, an histogram like Fig. 19 is
basically the record of how many times each bus produces a constant
lrp
=
e
ratio along
the entire construction of Fig. 18.
Clearly, the
at areas of the
lrp
=
e
curves are revealing that the
grid/network has serious topological defect that creates high curvature centrality (hence
congestion) as evidenced by Fig. 19. Unfortunately, the topology of the grid can only be
changed at signicant cost, so that we will have to nd less costly alternatives to manipulate
the curvature.
2.4 Ollivier-Ricci Curvature Driven Optimal Power Flow (ORC-
OPF)
2.4.1 Negative Curvature Detection via ORC
We start this section by building up aP -graph of the IEEE 300 bus system; this is done
following Sec. 2.1. Once the resistive network of the power grid is acquired, we calculate
the Ollivier-Ricci curvature of every edge in the model, aiming at detecting those edges
that have negative curvature. Bear in mind that what we actually have is a pure resistive
network model (P -graph) that can be seen as an undirected weighted graph; thus, an ORC
calculation for every edge is straightforward using Sec. 2.2.
Fig. 20 shows the ORC value for every edge in the grid. We actually have 409 edges
within the grid (x-axis), and almost 25 edges show a negative ORC value (y-axis).
32
Ollivier-Ricci Curvature value for every Branch/Edge of the grid/graph
0 50 100 150 200 250 300 350 400 450
Edge/Branch #
-3.5
-3
-2.5
-2
-1.5
-1
-0.5
0
0.5
1
ORC value
Figure 20: Curvature values of critical edges
Observe that this brief section has identied the cause/reason to propose a congestion
management procedure; it has detected the critical edges, which are the lines prompt to be
overloaded within the grid. Therefore, if we are able to smooth over the curvature throughout
the grid, we will be avoiding congestion. This is actually what is done in the next section,
and it is basically a heuristic procedure that makes positively curved those areas of the grid
originally negatively curved.
Note that next sections do not necessarily follow the edge/branch numbering utilized in
the IEEE300 bus system common data. On the other hand, they do follow the bus/node
numbering. Nevertheless, have in mind that the numbering does not change the curvature
analysis outcome.
33
2.4.2 FACTS, Storage Devices and the full ORC-OPF
Recall from Sec. 2.1 that theP -graph is basically a resistive network model composed of

km
= 1=B
km
V
k
V
m
resistances, wherek andm are dierent buses of the grid connected with
an edge, V
k
and V
m
are the voltages at buses k and m resp., and B
km
is the susceptance
value of the line (edge) that joins buses k and m.
Clearly, in order to change the negative curvature value of an edge, we need to change
the resistance value (distance) of such edge, with the aim of obtaining a positive curvature
value for it.
The only possible variables that we have available for changing the resistances are the
voltagesV
k
andV
m
, and the susceptanceB
km
, but if we leave the voltages unchanged for the
grid voltage stability operations, we are only left with the susceptance valueB
km
( 1=X
km
).
This is actually the only choice due to the DC assumptions (mostly the G
km
0 as-
sumption), which sets the format of the
km
to be dependent only on voltages (V
k
and V
m
)
and the susceptance (B
km
).
This seems to be a serious limitation, because just a change in the suceptance might
not be enough to do a nontrivial load balancing of the grid. Here is the approach that will
be taken: for the purpose of the ORC calculation on each edge, we will enforce the DC
assumptions in order to be able to have a relevant P -graph and consequently a curvature
value for each edge. Once the curvature values are computed for all lines and the negative
curvature areas are spotted, we will switch to AC conditions to smooth out the curvature
(which is actually more realistic). This will allow us to have G
km
6= 0, and consequently we
will have a new variable (R
km
) to adjust, which in turn will facilitate the curvature smoothing
process. Actually, while no control can be directly exercised upon the line susceptances, the
apparent susceptances can be modied by, for example, FACTS series compensation to
modify line impedance and static synchronous series compensator (SSSC) that connects an
inductive or capacitive reactance in series with the transmission line.
This curvature smoothing process (making the negative curvature areas positive using
34
FACTS) has been done heuristically and it was done at the same time a collection of loads
were deployed in the surroundings of the critical edges. Clearly, this second stage of the
proposed method allows energy storage by converting the power consumed by the deployed
loads to Gibbs free energy, while at the same time minimizing the overall cost of generating
active and reactive power within the grid.
Once the load balancing is performed and the loads are deployed, a new set of power

ow equations are embedded in a convex optimization algorithm (Algorithm 2.1), with the
objective of minimizing a polynomial cost function of the active and reactive power of each
generator. Clearly, this is also done under the AC assumptions. We illustrate below the
structure of the nonlinear programming algorithm:
Algorithm 2.1: AC Cost Optimization(;V;P
g
;Q
g
;)
min
;V;Pg;Qg
P
gensize
k=1
C
AC
(P
g;k
;Q
g;k
);
subject to
8
>
>
>
>
>
>
>
>
>
>
>
>
>
>
<
>
>
>
>
>
>
>
>
>
>
>
>
>
>
:
F
AC
(;V;P
g
;Q
g
) = 0;

i

i

i
; i = 1;:::;bussize;
V
i
V
i

V
i
; i = 1;:::;bussize;
P
g;k
P
g;k

P
g;k
; k = 1;:::;gensize;
Q
g;k
Q
g;k

Q
g;k
; k = 1;:::;gensize:
return (;V;P
g
;Q
g
)
In the algorithm, P
g;k
and Q
g;k
stand respectively for the active and reactive power
generated by generatork, gensize is the number of generators in the grid,x = [;V;P
g
;Q
g
]
is the optimization state variable where is the phase angle vector carrying the bus phase
angle, V stands for the bus voltages vector, P
g
andQ
g
are respectively the vectors of active
and reactive powers generated by the generators; bussize is the number of buses in the grid;
C
AC
(:) is a degree-2 cost function that weights the cost of generation of each generator k:
35
Table 2: Ollivier-Ricci Curvature Driven OPF \Main Steps"
ORC-OPF
(2.4.2)
identify critical buses
adjust B's and R's (w / FACTS) of critical lines
deploy loads around critical buses
run AC Cost Optimization Algorithm 2.1
C
AC
(P
g;k
;Q
g;k
) =
g;k
(P
g;k
)
2
+
g;k
(P
g;k
)+
+
g;k
(Q
g;k
)
2
+
g;k
(Q
g;k
) +
g;k
:
(15)
F
AC
() = 0 represents the dynamic of the AC power
ow. Finally, (P
g;k
,

P
g;k
), (Q
g;k
,

Q
g;k
),
(
i
;

i
) and (V
i
;

V
i
) are the min and max limits forP
g;k
,Q
g;k
,
i
andV
i
, respectively. Observe
that the cost function is composed of
0
gensize
0
order-two polynomials that could be built
up dierently for each generator; thus, we can weigh (choosing
g;k
,
g;k
,
g;k
,
g;k
and
g;k
)
each generator cost dierently by shaping each polynomial separately.
The combination of the dierent steps made till this point constitutes what we call the
Ollivier-Ricci Curvature Driven OPF (ORC-OPF), and it is summarized in Table 2.
The overall procedure has been implemented in Matlab using the MATPOWER package.
A modied version of the Ulas Yilmaz ORC software has been generated and implemented
based on [40]. The results and a comparison with a standalone AC OPF optimization method
applied to the IEEE 300 bus system are summarized in Table 3. The total amount of storages
deployed were 24, all of them deployed in the surroundings of critical lines. Also, a total of
30 critical lines were adjusted (load balancing) so as to maximize the curvature.
Observe from Table 3 that the optimal cost of generation proposed by the ORC-OPF
spend almost 15% less energy than the optimal generation proposed by the standalone OPF;
thus, these numbers clearly show the eectiveness of the proposed method. Moreover, it can
be claimed now that the structure of the grid is no longer so prompt to congestion as the
original one, because the ORC-OPF method not only has achieved less energy usage, but
36
Table 3: Total cost functions values (AC Model with conventional OPF and with ORC-OPF
implementation)
(dollars/hr) Total Cost Function Value
with conventional OPF 719730.00
with ORC-OPF 613730.00
Ollivier-Ricci Curvature value for every edge after implementation of the ORC-OPF
0 50 100 150 200 250 300 350 400 450
Edge/Branch #
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
ORC value
Figure 21: Curvature ealues of critical edges (No negative curvature values are present.)
it has also smoothed over the negative curvatures areas indicated in Fig. 20. This can be
clearly seen in Fig. 21, where no trace of negative ORC values can be found|and this is
actually the main goal of the method. Notice that Fig. 21 shows more edges than Fig. 20;
the reason is that the edges belonging to the deployed loads added after the ORC-OPF
implementation have been included. We also highlight the fact that in order to calculate
the ORC curvature we are using a graph abstraction, the P -graph. This is basically done
because the ORC needs a graph and a distance to be computed. As a consequence, we utilize
the DC assumptions to eliminate the non-linearities and generate a simple graph/model of
37
the grid. Although we worked with AC assumptions during the load balancing and the OPF
energy optimization, we analyze part of the results (curvature maximization) under the DC
assumptions, because it is the fundamental tool to generate the ORC values that guide the
overall method. We however choose the AC model within the main steps of the procedure,
just to be closer to the real behavior of the grid, where dissipation and non-linearities are
present.
2.4.3 Eective Resistance Curvature Driven Optimal Power Flow (ERC-OPF)
Recall that, under the Eective Resistance Curvature notion, the critical nodes are iden-
tied via the
lrp
(k;m)=
e
(k;m) tool/plots shown in Fig. 18 and such critical nodes are
more clearly shown in Fig. 19 along with their centralities. This type of plots are basically
showing how many transmission lines within the grid are behaving as overloaded paths,
which is captured by the ratio
lrp
(k;m)=
e
(k;m). In order to understand the identication
of these particular paths, recall that the branches of a tree are always congested, and that
their
lrp
(k;m)=
e
(k;m) ratios are always equal to one, because the shortest path resistance
coincides with
e
(eective resistance) in the branches of a tree. It turns out that the Eec-
tive Resistance Curvature notion captures more negative areas/critical buses than the ORC
notion. Therefore, we should expect a better perfomance than the ORC-OPF; and this is ac-
tually the case. Applying exactly the same FACTS/storages allocation procedure done in the
ORC-OPF but utilizing the critical buses depicted in Fig. 19, the overall ERC-OPF method
yields the cost-reduction outcome shown in Table 4. Clearly, if we compare Table 3 and
Table 4, we observe that the ERC-OPF approach utilizes an optimal set of generators that
spend almost 30% less energy than the conventional OPF case, and a 20% less energy that
the ORC-OPF. All these ERC-OPF improvements with respect to the ORC-OPF are easily
explained by the fact that the Eective Resistance Curvature is a global curvature notion,
and that the Ollivier-Ricci Curvature is a local edge curvature calculation; thus, the ERC
captures a proper global snapshot/state of the overall grid in each step of its calculation.
38
The drawback of this ERC-OPC approach is that is more time consuming than the
ORC-OPF due to its higher complexity. Therefore, depending on the application, an a
priori suitable choice of the curvature notion before the implementation of the CM approach
might be pertinent.
Table 4: Total cost functions values (AC Model with conventional OPF and with ERC-OPF
implementation)
(dollars/hr) Cost Function Value
with conventional OPF 719730.00
with ERC-OPF 492280.00
2.5 Grid Entropy Measure
This section develops an entropy-based measure suitable for weighted graphs. In our
case, a weighted graph would be a P -graph resistive network model of the power grid; a
graph that we have utilized to construct the two dierent curvature notions that gave rise
to this CM approach. The idea is to dene a measure that is positively correlated with the
curvature increment of the grid [51]. Recall that the proposed CM method implements a
curvature maximization as a core step towards an ecient cost-eective curvature driven
OPF; thus, to have a topological index that measures this increment would be helpful and
handy for a simple inspection of the new curvature state of the grid. The utilization of graph
metrics to analyze the structure of networks have been extensively used in many dierent
elds such as: computer sciences, biology and mathematical chemistry [53]. Most of this
work founds its origins in the well known work of Shannon in the 1950's. Other authors
further developed the entropy concept and applied it to graphs (Rashevsky, Mowshowitz,
Chen et al, [52]). Here we'll use the R. Kazemi [53] approach:
Denition 17 (Grid Entropy): For an edge weighted graph G = (V, E, w), the entropy of G
is dened by:
39
I(G;w) =
X
uv2E
p
u;v
log(p
u;v
); (16)
where
p
u;v
=w(u;v)=
X
uv2E
w(u;v): (17)
In our case, the weighting function w(uv) will be the resistance value of each edge. In order
to be able to compare the new measure and the curvature we utilize the Global Ollivier-Ricci
curvature value, which is just the average curvature value of the edge curvatures previously
utilized (Fig. 20 and Fig. 21).
Although the IEEE 300 bus system presents negative curvatrue areas, the overall (av-
erage) Ollivier-Ricci curvature value is a positive quantity; thus, after the implementation
of the ORC-OPF we should expect an average curvature increment within the grid. This
is clearly explained by recalling that the CM approach embraces a curvature maximization.
Therefore, we should also expect (as per [51] ) a positive increment in the entropy; this can
be observed in Table 5, where a global Ricci curvature values has been compared against the
entropy of the power grid, before and after the implementations of the ORC-OPF.
Table 5: Global curvature values vs grid entropy
(IEEE300 Bus, AC P -graph ) Global ORC Value Grid Entropy
before ORC-OPF 0.8633 5.2608
after ORC-OPF 0.9676 6.1607
As we mentioned above, the suggested CM approach maximizes the curvature; thus, it
reduces the negative areas of the grid, making a more uniform or homogenous grid. Therefore,
to obtain a higher entropy after the ORC-OPF implementation is natural and physically
expected. If we go back to the second law of thermodynamics, we clearly have a higher
entropy value (in an isolated system) as we move towards an equlibrium; and this is actually
what is happening as we smooth out the curvature of the grid with the CM approach. On
40
the other hand, an initial heterogenous state of the grid (with positive/negative curvature
areas disseminated throughout the grid) cannot have another entropy value than a lower one
compared with a posteriori homogenous state (after grid curvature maximization).
2.6 Line Rating Considerations
As mentioned in previous sections, all of the simulations have taken into account line
rating considerations; this has been accomplished by restricting the active power
owing
through the lines to a maximum of 700MW. This is actually a way to account for the
capacity of the lines, usually determined by the thermal rating [56,58,59].
In particular, the following denitions have been applied:
Denition 18 (DC Utilization Factor): The utilization factor for the branch (k;m) under DC
model assumptions is dened as

DC
=P
k;m
=LC
k;m
;
where LC stands for line capacity, the maximum active power allowed (in MW) through the
branch (k;m).
Denition 19 (AC Utilization Factor): The utilization factor for the branch (k;m) under AC
model assumptions is dened as

AC
=
q
P
2
k;m
+P
2
m;k
=LC
k;m
:
The impact of adding the line limitation constraints

DC;k;m
1;
AC;k;m
1;
in the DC and AC ORC-OPF approach is barely noticeable in the overall nal cost function
41
Branch Utilization
0 50 100 150 200 250 300 350 400 450
Branch #
0
20
40
60
80
100
120
Percentage %
Figure 22: Utilization Factor (
AC
) for the AC model of the IEEE300 bus system with
ORC-OPF implementations and a 700 MW active power constraint
values for a line limit of LC
k:m
= 700MW, although it underlines an important advantage
within the complete method: the ORC-OPF scheme is able to handle realistic line limits.
Fig. 22 shows a line utilization histogram (in percentage) for the IEEE 300 bus system
with a line rating of 700 MW (on the active power of the branches) under AC analysis. As
mentioned earlier, the line rating inclusion within the proposed load balancing encompasses
a promising analysis tool: it would basically consists in the direct inclusion of thermal rating
considerations for bare overhead conductors.
The increase of thermal stress due to variable weather or other conditions [54,55,58,59]
could easily trigger a line overloading that might end up in a blackout (e.g., 1996 Western
North America blackout [56]). As another scenario, when a major line trips, power is rerouted
along other lines that may not have been designed to carry such an amount of power and
hence are likely to be overloaded and trip, leading to a chain reaction eect [57]. Thus, a real
42
time power
ow calculation that includes a DLR (Dynamic Line Rating [58]) could certainly
help to assess power grid functionality.
43
3 Anticipation of Line Tripping Cascade by Grid Cur-
vature and Entropy
Thermal stress in power grids can easily produce a line overloading that can end up in
a blackout (Western North America blackout [56]). A variable weather or other conditions
( [54], [55], [58], [59]) can produce an overheating of bare overhead conductors causing a
line tripping; as a consequence, an important amount of power can be rerouted through lines
that are not designed to carry such an amount of power. Consequently, these lines could
also trip and start a chain reaction eect [57].
The main objective of this section is to try to identify those lines that are prone to be
overloaded after the main line/lines is/are tripped. To accomplish this, the Ollivier-Ricci
Curvature (ORC) concept is utilized and presented as a powerful tool to identify areas under
risk.
The foundations of this section were set in [22], [23], and [24]; recall from those papers
that the ORC is based on the construction of a resistive network model of the grid. The
power of this formulation is that the core structure/inner interconnection of the power grid is
exposed and that the vulnerable areas are precisely identied using the notion of curvature.
3.1 Ollivier-Ricci Curvature Approach to Detect Power Grid Ar-
eas Under Risk ( Fault / Overheating )
This section will basically divide the ORC approach in two parts; the rst part summa-
rizes the plain ORC detection of negative curvature areas [23]. The second part identies
lines/areas under risk after a line tripping scenario utilizing a second ORC calculation over
the new tripped grid. This last calculation is compaired against the initial ORC values and
the under risk lines will be identied as those with a net negative curvature increment.
44
3.1.1 Negative Curvature Detection via ORC
We start this section by building up the S-graph of the IEEE 300 bus system; this is done
following Sec. 2.1. Once the resistive network of the power grid is acquired, we calculate
the Ollivier-Ricci curvature of every edge in the model, aiming at detecting those edges
that have negative curvature. Notice that the resistive network model (S-graph) chosen
can be seen as an undirected weighted graph; thus, an ORC calculation for every edge
is straightforward using Sec. 2.1. Bear in mind that the resistive network models/graphs
mentioned in Sec. 2.1 require the solution of the power
ow equations a priory (bus voltages
and bus voltage angles calculations), and that the model (AC or DC model) adopted to
solve those sets of equations might slightly change the nal resistive network model. Having
said that, this work opted to work with the AC model of the power grid and to utilize the
Optimal Power Flow (OPF) set up to calculate the voltage magnitude and voltage angle of
every single bus. This allows the inclusion of a degree-2 cost function that weights the cost
of generation of each generator, which is actually what is usually done in most of the grids
that consider Congestion Management (CM) within its regular operation. Therefore, this
work clearly proposes a complete power grid model just from the very beginning, aiming at
yielding precision and error reduction.
The AC model of the power grid together with the OPF calculation were implemented
in the MATPOWER environment, and the convex optimization algorithm that solves the
power
ow equations (minimizing a polynomial cost function of the active and reactive power
of each generator) becomes a nonlinear programming as depicted in algorithm 2.1.
3.1.2 ORC Identication of Areas/Lines/Corridors under Risk
3.1.3 First Simulation
After the ORC values of the original power grid are computed (Sec. 3.2.1), a line of the
grid is forced to be tripped to emulate a fault in the grid. Once this is done, a new ORC set
45
of values for every edge is calculated for this new tripped network. Have in mind that the
entire OPF (optimal power
ow) calculation ought to be performed again because the ORC
formula utilizes the voltage and the voltage angle of every bus, and the tripping of just a
single line will change the entire OPF outcome; thus, a new resistive network model will be
obtained after solving a new OPF problem.
Once the second set of ORC values for each edge is obtained (Sec. 3.2.1), the dierence
between the ORC values of the tripped network and those of the non-tripped grid are an-
alyzed. This ORC curvature dierence is depicted in Fig. 23. The performed simulation
tripped line 250 of the grid (connected with buses 152 and 154) and imposed a 700 MW
active power limit in every remaining line. By observing the gure we clearly notice that
some lines are becoming less congested (positive curvature increment) and some of them are
actually incrementing its negative ORC value (Braess' paradox [60]).
So those lines that are prompt to be under risk after a tripping of line 250 are clearly the
ones with highest negative curvature values (recall that the more negative the curvature is,
the more congested the line will be). If we take a second look at Fig. 23 we notice that edges
245, 249, 250 and 251 are the most aected ones after the tripping of line 250. Nevertheless,
lower branches from 1 to 36 , middle branches (194-202 and 216) and upper branches (360
and 362) have also been aected; thus, the tripped line have actually aected most of the
regions of the grid. Fig. 24 and Fig. 25 show the aected edges more closely.
The aected areas can be seen more clearly in Fig. 26, where the IEEE300 Bus system
has been drawn and the aected edges have been highlighted in red with dierent widths to
account for the level of congestion they might acquire if the trip occurs.
Observe from Fig. 23, Fig. 24 and Fig. 25 that lower, middle and upper branches are
aected; this is basically due to the topology of the grid. Bus 216 is connected to lower
branches through bus 198; at the same time, bus 198 is connected to upper branches through
bus 301. Therefore, these buses are clearly forming inner grid corridors that can spread the
excess of active power re-routed due to the tripping of line 250. The identication of this
46
type of corridors is essential if the grid designer is planning to be able to hold an excess of
forwarded active power as a consequence of a cascade eect/line tripping.
3.1.4 Second Simulation
The upcoming simulation will utilize the same AC model of the power grid, the same
OPF formulation/solver and the sameS-graph resisitve network model as of rst simulation
to deploy the ORC analysis; the subtle dierence will be that the active power limit of
the branches will increase to 900 MW and that the tripped line will be line 322. This
line/edge/branch works initially at 90 % of its capacity, which means that it is a branch
that could easily damage adjacent branches if it were to be tripped. The branch capacity is
depicted in Fig. 29, where the AC branch utilization factor developed in Sec. 2.6 is used. The
importance of this branch resides in the fact that one of its end nodes (bus 220) is actually a
PV bus; thus, it dispatches an important amount of active power towards dierent loads (PQ
buses) in the grid. Nevertheless, have in mind that a branch working far from its capacity
limit could also be a dangerous line if tripped due to its potential not convenient topological
location.
Fig. 27 shows the nal result of the ORC analysis for this new scenario. Clearly the
absence of the tripped line (edge 322) shows its eect throughout the entire grid, having
its strongest impact on lines 317 and 320 (Fig. 28) which show the highest negative peaks.
Once again, the key cause of this spread is easily understood by identifying the corridor
that spreads the active power that used to go over the tripped line 322 (buses 220 and 216).
Bus 216 is connected to lower buses through bus 198; thus, the impact is seen all the way
to lower branches. On the other hand, bus 198 is also connected to bus 301; thus, the
forwarding is also seen all the way to higher branches. Therefore, almost every region of the
grid experiences the tripped of line 322; this can be farther appreciated in Fig. 30, where the
complete IEEE300 Bus connectivity is drawn and the aected edges (after the line tripping)
have been highlighted in red with dierent widths to emphasize the level of compromise of
47
Ollivier-Ricci Curvature Difference between tripped and no tripped line grid
0 50 100 150 200 250 300 350 400
Edge/Branch #
-0.1
-0.08
-0.06
-0.04
-0.02
0
0.02
0.04
ORC value
Figure 23: ORC Values Dierence between tripped and no tripped lines (Tripping of line
250 with 700 MW branch limit)
each risky edge.
As can be seen, with this type of risk analysis we can infer a priori what would be the
consequences of loosing a line for any reason (fault, overheating, etc); thus, a proper measure
that could be taken is to re-reinforce those lines that could be aected by the tripping of
any major line that carries power usually within its capacity limits. And this is basically the
message of the proposed approach, because it can help the designer/grid company developer
to determine in advance which path a possible cascade eect will take.
3.2 Eective Resistance Curvature Approach to Detect Power
Grid Areas Under Risk ( Fault / Overheating )
This section follows closely the simulations performed within the previous section, the
main dierence is that a dierent curvature measure (Eective Resistance Curvature, ERC)
is utilized to identify the areas under risk and the corridors that connect them.
48
Ollivier-Ricci Curvature Difference between tripped and no tripped line grid
200 250 300 350 400
Edge/Branch #
-0.08
-0.06
-0.04
-0.02
0
0.02
ORC value
Figure 24: ORC values dierence between tripped and no tripped lines (zoom in of lines 150
to 409 of Fig. 23, tripping of line 250 with 700 MW branch limit)
Ollivier-Ricci Curvature Difference between tripped and no tripped line grid
0 10 20 30 40 50 60 70 80 90
Edge/Branch #
-0.05
-0.04
-0.03
-0.02
-0.01
0
0.01
ORC value
Figure 25: ORC values dierence between tripped and no tripped lines (zoom in of lines 1
to 100 of Fig. 23, tripping of line 250 with 700 MW branch limit)
49
IEEE300 Bus connectivity highlighting affected edges (in red) after the trip of line 250
Figure 26: IEEE300 Bus connectivity highlighting aected edges (in red) after the trip of
line 250 (700MW branch limit). The width of each red edge represents the degree to which
each edge has been aected.
Ollivier-Ricci Curvature Difference between tripped and no tripped lines
0 50 100 150 200 250 300 350 400
Edge/Branch #
-0.025
-0.02
-0.015
-0.01
-0.005
0
0.005
0.01
0.015
0.02
0.025
ORC value
Figure 27: ORC values dierence between tripped and no tripped lines (Tripping of line 322
with 900 MW branch limit.)
50
Ollivier-Ricci Curvature Difference between tripped and no tripped lines
305 310 315 320 325 330 335 340
Edge/Branch #
-0.02
-0.015
-0.01
-0.005
0
0.005
0.01
0.015
0.02
ORC value
X 320
Y -0.02103
Figure 28: ORC values dierence between tripped and no tripped lines (zoom in of lines 317
and 320 of Fi. 27 after tripping of line 322 with 900 MW branch limit)
Branch Utilization
0 50 100 150 200 250 300 350 400
Branch #
0
20
40
60
80
100
120
Percentage %
X 322
Y 90.57
Figure 29: Branch utilization factor (Line 322 utilization factor has been highlighted .)
51
IEEE300 Bus connectivity highlighting affected edges (in red) after the trip of line 322
Figure 30: IEEE300 Bus connectivity highlighting aected edges (in red) after the trip of
line 322 (900MW branch limit). The width of each red edge represents the degree to which
each edge has been aected.
3.2.1 Negative Curvature Detection via ERC
This entire subsection deploys the ERC calculations (Sec. 2.2.3) using theS-graph resis-
tive network model (Sec. 2.1), the AC model of the power grid, and the OPF calculations
performed and implemented in MATPOWER utilizing algorithm 2.1.
3.2.2 ERC Identication of Areas/Lines/Corridors under Risk
3.2.3 First Simulation
Initially, the ERC values of the original power grid are computed (Sec. 3.2.1). As a second
step, a line of the grid is tripped to emulate a fault in the grid and a new ERC set of values
for every node is calculated for this new tripped network. Recall from 3.1.3 that the entire
OPF (optimal power
ow) calculation has to be performed again because the ERC formula
utilizes the voltage and the voltage angle of every bus, and the tripping of just a single line
will change the entire OPF results; thus, a new resistive network model will be obtained
52
after solving a new OPF problem and before a second set of ERC values are computed.
Once the second set of ERC values for each node is obtained (Sec. 3.2.1), the dier-
ence between the ERC values of the tripped network and those of the non-tripped grid are
analyzed. This ERC curvature dierence is depicted in Fig. 31.
The actual simulation imposed a 700 MW active power limit in every line of the grid and
forced the trip of line 220 which is actually connected to buses 119 and 121. Fig. 31 shows
that some nodes are becoming less congested and some of them are actually incrementing
its negative ERC value. In order to properly read Fig. 31 have in mind that a positive (a)
value for node a implies an increment in its negative curvature (denition 16). Therefore,
all negative(a) values are actually saying that those nodes reduce their negative curvature
value (which is actually good). So those nodes that are prompt to be under risk after a
tripping of line 220 are clearly the ones with highest negative curvature values, which means
a higher curvature centrality(a) value. Recall that the more negative the curvature is, the
more congested the node will be.
If we take a second look at Fig. 31 we notice that nodes 268, 269, 280, 281, 283, 284,
285, 286, 287, 289, 290, 291, 296, 297, 298 are the most aected ones (highest (a) values)
after the tripping of line 220, although the impact is also seen at lower buses like 23, 54, 55,
59, 141, 143, 144 and 145. Fig. 32 and Fig. 33 show the aected nodes more closely.
If we translate this into edge numbering, the most negatively aected edges were edges:
11, 21, 22, 23, 24, 25, 26, 27, 29, 30, 31, 32, 33, 34, 35, 36, 115 to 117, 120 to 124, and 242
to 245; thus, the impact is basically seen in lower, middle and upper branches. Nevertheless,
although all these branches seem to be scattered across dierent regions, this does not
mean that they are far apart in distance within the grid; because the bus/edge numbering
adopted does not represent a distance or a measure. Moreover, the bus/edge numbering
could be changed at anytime and the curvature analysis will pick again the same aected
regions/connections as it can be expected. In order to better highlight the latter comment,
Fig. 39 shows the IEEE300 Bus connectivity highlighting aected edges (in red) after the
53
trip of line 220 (with a 700MW branch limit). The width of each red edge represents the
degree to which each edge has been aected.
A deeper analysis/inspection of Fig. 31 and Fig. 39 indicates that all these branch areas
are connected through buses 23 (connected to upper buses 231 and 254), 54, 55 and 59
(connected to middle bus 145). Clearly, these buses are forming inner grid corridors that can
spread the excess of active power rerouted due to the tripping of line 220. The identication
of this type of corridors is essential if the grid designer is planning the grid to be able to
hold an excess of forwarded active power as a consequence of a cascade eect/line tripping;
and gures such as Fig. 31 and Fig. 39 are certainly essential to help the identication of
them.
3.2.4 Second Simulation
This simulation is performed under the same model assumptions (S-graph and AC model
for OPF) as of previous section, although now the active power limit of the branches is
increased to 900 MW and the tripped line is line 266. This line initially works at 64.17 %
of its capacity, which means that it is a branch that could easily damage adjacent branches
if it were to be tripped. This branch capacity is depicted in Fig. 38, where the AC branch
utilization factor developed in Sec. 2.6 is used. This branch (266) is composed of buses 169
and 210, where bus 169 is actually a PV bus; thus, it dispatches an important amount of
active power towards dierent loads (PQ buses) in the grid.
Fig. 35 shows the nal result of the ERC dierence under this new scenario. We observe
that this new histogram shows more serious spread of the active power of the tripped line.
Clearly the absence of the tripped 266 line is also aecting three dierent areas (as in 3.2.3);
the dierence now is that buses 15, 17, 213, 215, 236, 237, 254, 261, 262, 265, 279, 282,
297 and 299 are also included within the aected buses (buses 280, 296 and 298 were not
negatively aected this time).
Therefore, the lower/middle/upper branch area is expanded compaired to 3.2.3; it now
54
Crtitical Nodes Repetition Difference
0 50 100 150 200 250 300
Node #
-40
-30
-20
-10
0
10
20
30
# of repetition difference for each node
Figure 31: Curvature centrality ((a)) dierence between tripped and no tripped lines (Trip-
ping of line 220 with 700MW branch limit.)
includes edges 11, 21-37, 41, 46, 49-58, 64, 67 112-118, 121-126, 242-245, 316-322, 370, 374,
388, 399. The responsible buses for this widening of the risk area are more easily seen in
Fig. 37 and Fig. 36; thus, new corridors (compaired to 3.2.3) are identied and they are
formed by buses 15, 17, 213, 215 and 237.
Fig. 39, as in rst simulation, is also provided as a mean to visually identify the dierent
edges/branches that could be aected as a consequence of the trip of a line (line 266 in this
case). The bus numbering has been omitted to preserve the clarity of the gure. It is worth
to notice that the branch that works close to its capacity limit is not necessarily the most
dangerous line to be tripped; on the other end, a branch that works with a fair line capacity
margin could be topologically located in a region that might be in
icted the worst impact
if the line were to be tripped. This type of conclusion/analysis can only be revealed by the
type of curvature analysis proposed by this work.
As can be seen, with this type of risk analysis we can infer a priori what would be the
consequences of loosing a line for any reason (fault, overheating, etc); thus, a proper measure
55
Crtitical Nodes Repetition Difference
240 250 260 270 280 290 300
Node #
-30
-20
-10
0
10
20
# of repetition difference for each node
X 268
Y 26
Figure 32: Curvature entrality ((a)) dierence between tripped and no tripped lines (zoom
in of nodes 240 to 300 of Fig. 31)
Crtitical Nodes Repetition Difference
0 50 100 150
Node #
-2
-1
0
1
2
3
4
5
6
7
# of repetition difference for each node
X 141
Y 3
Figure 33: Curvature centrality ((a)) dierence between tripped and no tripped lines (zoom
in of nodes 1 to 150 of Fig. 31)
56
IEEE300 Bus connectivity highlighting affected edges (in red) after the trip of line 220
Figure 34: IEEE300 Bus connectivity highlighting aected edges (in red) after the trip of
line 220 (700MW branch limit). The width of each red edge represents the degree to which
each edge has been aected.
Crtitical Nodes Repetition Difference
0 50 100 150 200 250 300
Node #
-50
-40
-30
-20
-10
0
10
20
# of repetition difference for each node
Figure 35: Curvature centrality ((a)) dierence between tripped and no tripped lines (Trip-
ping of line 266 with 900MW branch limit.)
57
Crtitical Nodes Repetition Difference
0 50 100 150
Node #
-6
-4
-2
0
2
4
6
# of repetition difference for each node
Figure 36: Curvature centrality ((a)) dierence between tripped and no tripped lines (zoom
in of nodes 1 to 150 of Fig. 35)
Crtitical Nodes Repetition Difference
255 260 265 270 275 280 285 290 295 300
Node #
-40
-30
-20
-10
0
10
20
# of repetition difference for each node
Figure 37: Curvature centrality ((a)) dierence between tripped and no tripped lines (zoom
in of nodes 240 to 300 of Fig. 35)
58
Branch Utilization
0 50 100 150 200 250 300 350 400
Branch #
0
10
20
30
40
50
60
70
80
90
100
Percentage %
X 266
Y 64.17
Figure 38: Branch utilization factor (Line 266 AC utilization factor has been highlighted.)
IEEE300 Bus connectivity highlighting affected edges (in red) after the trip of line 266
Figure 39: IEEE300 Bus connectivity highlighting aected edges (in red) after the trip of
line 266 (900MW branch limit). The width of each red edge represents the degree to which
each edge has been aected.
59
that could be taken is to re-reinforce those lines that could be aected by the tripping of
any major line that carries power usually within its capacity limits. And this is basically the
message of the proposed approach, because it can help the designer/grid company developer
to determine in advance which path a possible cascading eect will take.
60
4 Local/Nodal Entropy
To further analyze the grid changes after a line tripping, we develop in this section the
notion of nodal entropy. This nodal entropy concept applies to weighted graphs, in particular
in our case to the resistive S-graph network model of the power grid (Sec. 2.1). This local
measure will be positively correlated with the curvature increment of the grid only locally.
The utilization of the entropy measure to analyze networks has been extensively used [53].
Further reading concerning entropy applied to graphs can be found in [52]. Here, we will
follow the entropy approach of R. Kazemi:
Denition 20 (Nodal Entropy): For an edge-weighted graph G = (V;E;w), the nodal entropy
of node k is dened by
I
k
(G;w) =
X
m2N (k)
p
k;m
log(p
k;m
); (18)
where
p
k;m
=
w(k;m)
P
m2N (k)
w(k;m)
: (19)
In our case, the weighting function w(k,m) is the resistance of each edge km andN (k),
the neighborhood of k, is the set of all vertices m such that km2E.
If we compute this local/nodal entropy before and after a line tripping and compare
(substract) the results locally, we will notice that buses near the tripped line will have
undergone a drastic entropy changes as can be seen in Fig. 40 and Fig. 41. Both gures were
obtained within the simulation framework of Sec. 3.2.1 and Sec. 3.2.2, where lines 220 and
266 were tripped; thus, it is evident that buses 119 and 121 (line 220) and buses 169 and
210 (line 266) at the distal ends of the tripped lines undergo an important reduction in their
local/nodal entropy compared to their surrounding/local bus neighbors.
Observe that an excess of active power through line km due to the rerouting of active
power after overheating and tripping ofkm
#
is accompanied by an increment in the dierence
of voltage angles between buses k and m. Therefore, using the S-graph model, we observe
61
that the resistance
km
= 1=(B
km

V
k

V
m
cos(

k

m
)) will also increase, but(km) =w(k;m);
thus, the tripped line will create an uneven distribution off
kn
:n2N (k);n6=m
#
g resulting
in a reduction of the local/nodal entropy of the bus k.
Entropy of nodes
0 50 100 150 200 250 300
Node #
-0.35
-0.3
-0.25
-0.2
-0.15
-0.1
-0.05
0
0.05
Entropy
X 121
Y -0.3396
Figure 40: Nodal entropy as per Def. 20 (Line 220, connecting buses 119 and 121, has
tripped.)
This local/nodal entropy could be further extended to other stochastic probabilities, like
p
km
=
~
P
km
P
m2N(k)
~
P
km
, and it could be used in fault detection. Running this local/nodal entropy
in parallel to a curvature analysis could point out the exact location of an/a attack/line
tripping. This is left for further research.
62
Entropy of nodes
0 50 100 150 200 250 300
Node #
-0.12
-0.1
-0.08
-0.06
-0.04
-0.02
0
0.02
Entropy
X 169
Y -0.1022
Figure 41: Nodal entropy as per Def. 20 (Line 266, connecting buses 169 and 210, has
tripped.)
5 Conclusions
5.1 Research Summary
The present work has proposed, in its early stage, the R
lrp
=R
e
ratio between two nodes
of a nite network as a robust electrical engineering method to detect areas of negatively
curved behavior in a network, areas co-located with areas of congestion. The theory was
developed for regular 2-dimensional spiraling lattices of uniform curvature. However, the
heterogeneous IEEE-300 bus model case-study demonstrated that the R
lrp
=R
e
test is able
to identify areas of congestion. From the formally provable paradigm that congestion betrays
negative curvature, it is hence conjectured that the R
lrp
=R
e
test applies to heterogeneous
networks in detecting areas of vanishing and negative curvature. The formal proof to conrm
the generality of this observation is left for future research.
In its more mature stage, this work has suggested a rst procedure called curvature
63
smoothing to alleviate the negative curvature eects of a power network. This procedure
itself consists of two stages:
1. The rst stage, based on a new curvature-inspired betweenness centrality, identies
critical areas.
2. The second stage alters the resistances/impedances of the lines, which could be achieved
using Flexible AC Transmission System (FACTS) technology.
In the IEEE300 bus system case-study, congestion in the sense of line overloading im-
mediately appears to be an issue. By choosing this non trivial example, and observing that
unrealistic changes of the line parameters are warranted to change the curvature, the anal-
ysis reveals severe restrictions that the topological/combinatorial properties of the power
network impose on what can be achieved in terms of load balancing understood as better
utilization of the transmission network resources.
As second procedure, this work has proposed an optimal cost reduction procedure con-
comitant with the extra load deployment. It consists of two parts:
1. Identication of critical areas.
2. Deployment of additional loads to divert power away from congestion areas and utiliz-
ing the accrued freedom to minimize the cost of generation. The minimization of the
generation cost can be done under the
DC power
ow assumption,
AC power
ow assumption.
Furthermore, the loads can also be designed to store energy for future use.
The second part of the early work was a combination of the curvature smoothing procedure
and the cost reduction procedure. The goal was to keep the power
ow of the lines (P
k;l
)
bounded or at least under control for those problematic lines detected in the body of this
64
work. This resulted in a procedure that drastically reduces the overall cost of the power
ow
based on the curvature knowledge of the grid.
The last part of this work has proposed an Ollivier-Ricci Curvature (Eective Resistance
Curvature) based procedure to track down the eects of a power grid faulty line, taking line
rating into consideration. Once the curvature analysis has identied the stress points, a
second round of curvature computation is implemented after the tripping of a line within
the grid. Those two ORC (ERC) results are stored and compared to do a proper contingency
analysis of the grid after a main line tripping. In order to further help the analysis of the
line tripping consequences a nodal entropy measure was proposed to visually understand the
eects of the line trip in individual/local buses. Clearly, as seen in Sec. 3, more than one
node can be simultaneously negatively aected; thus, it would be up to the grid regula-
tor/grid management to determine how many of those nodes/areas they are nancially able
to protect/reinforce in case a sudden line tripping occurs. Also, the corridor identication
sets up a framework to start protecting those hidden alleys that could carry enough rerouted
active power to aect areas further away from the tripping region/line.
Finally, line rating considerations considered throughout the entire proposed ORC (ERC)
approach of this work have certainly opened the possibility to include thermal rating calcu-
lations, resulting in powerful and complete methods.
5.2 Future Work|Near Term
The future step of this work will be to improve the performance and the eectiveness
of the early Joint Curvature Smoothing and Optimal Cost Reduction method. This will be
achieved by enriching the pool of control actions over the power grid model.
We already proposed to drain power via deployed loads and store it in appropriate energy
reservoirs, but we now propose to re-inject the active power drained by the deployed loads
and available in the reservoirs.
In order to achieve this, the power grid dynamic model should be modied to account for
65
energy injection u. Thus, the new DC curvature-cost optimization equations, taking energy
injection u for granted, will look like:
Algorithm 5.1: DC Cost Optimization(;P
g
)
min
;Pg
P
gensize
k=1
C
DC
(P
g;k
)
subject to
8
>
>
>
>
>
>
<
>
>
>
>
>
>
:
F
DC
(;P
g
; u) = 0

i

i

i
i = 1;:::;bussize
P
g;k
P
g;k

P
g;k
k = 1;:::;gensize
return (;P
g
)
where the linearized power
ow equations F
DC
() = 0 now include a xed energy injectionu
at the critical nodes i
[
's.
A control-enhanced approach will consider the possibility of controlling the re-injected
powers at the i
[
buses; thus, u could become part of the optimization variable set:
Algorithm 5.2: DC Cost Optimization(;P
g
;u)
min
;Pg;u
P
gensize
k=1
C
DC
(P
g;k
) +
P
deployedloads
k=1
B
DC
(u)
subject to
8
>
>
>
>
>
>
>
>
>
>
<
>
>
>
>
>
>
>
>
>
>
:
F
DC
(;P
g
;u) = 0

i

i

i
i = 1;:::;bussize
P
g;k
P
g;k

P
g;k
k = 1;:::;gensize
u
i
u
i
u
i
i = 1;:::;deployedloads
return (;P
g
;u)
whereB
DC
(:) is a degree-2 cost function on the new control action u.
66
Recall that behind this cost optimization procedure a curvature smoothing procedure is
also implemented. Although these two processes are presented as decoupled processes to
simplify the overall exposition of the ideas, a challenging simultaneous implementation of
the methods is warranted.
After these two options are implemented and analyzed, a third approach will include
the energy charging/discharging dynamics in the optimization of the generation cost. This
basically means that another constraint equation should be added in the optimization set
up. Initially the so-called SOC (state of charge for generators) will be proposed, which is the
simplest dynamic model that represents the charging/discharging nature of energy storage
reservoirs:
x[t + 1] =x[t] +u[t],
where u can take either positive or negative values. As it can be seen, this dynamics de-
pends on time (discrete); thus, a nite time horizon optimization should be implemented.
Therefore, a new time dependent optimization set up should be constructed.
Recalling that the curvature-cost optimization achieved better performance with the use
of the AC model, all the previous suggestions will also be implemented in the AC case, aiming
to obtain better results while using the more realistic (AC) model. On the other hand, an
optimization procedure that includes the hitherto unknown dynamical eects revealed by
the data driven approach of [49, 50] could further enhance the curvature/cost reduction
procedure.
Lastly, another step into the future of this work could be the automatic implementation
of a load balancing technique to protect those lines that might be under risk. Clearly the
real implementation of this second approach would have a cost management framework that
will limit the amount of lines to be protected.
Finally, a promising step into the future of this work will be the study of the subtle re-
67
lation between nodal entropy and a tripped line; this work sets up the foundation to further
investigate the possibility to utilize the entropy and the curvature as a fault detection tool.
The global entropy [22] and the nodal entropy are believed to be useful tools to detect an
attack or a fault in the grid with a further practical implementation of the theoretical prob-
ability used in their denition, and a parallel/simultaneous usage of the eective resistance
curvature of the grid.
5.2.1 Future Work|Longer Term
With a second priority and depending on time availability, a further idea of this work
would be to study the possibility to generate the theoretical basis to deploy the curvature
analysis/congestion detection procedure on directed graphs, like the DC reactive power graph.
Although directed graphs arise naturally in many applications, the entire theory developed
in this work can not be applied directly on such a graph. The main barrier that prevents
the application of the method is that the graph Laplacian on directed graphs is no longer
symmetric; thus, the entire R
e
calculation set up no longer holds true and this resistance
calculation is the very rst pillar on which the entire curvature-cost optimization method was
built. On the other hand, theR
lrp
calculation based on the Bellman's principle of optimality
would still hold true on directed graphs. Also a dierent notion of curvature should be
searched, because the Gauss notion of curvature would no longer be appropriate.
To start analyzing this challenging directed graph problem a very novel R
e
formulation
developed in [14] will be investigated. Initially, developed in the eld of study of stochastic
systems evolving on graphs, this new R
e
calculation is based on a covariance matrix which
is a solution of a Lyapunov equation that has the reduced graph Laplacian in his core
formulation. Following the notation in [14], the R
eff
calculation has the following aspect:
R
e
(k;`) =x
(k;k)
+x
(`;`)
2x
(k;`)
;
where
68
X = 2Q
T
Q;

L +

L
T
=I
N1
;

L =QLQ
T
:
HereL is the graph Laplacian,N is the number of nodes of the graph, andQ is a (N 1 N)
matrix.
The main drawback of this formulation is that the R
e
is no longer a distance/metric;
thus, the triangle inequality is no longer enforced. Nevertheless, under certain conditions,
the metric can be restored.
5.2.2 Potential Biomedical Application
A least priority objective of this work will be to apply the curvature analysis in the
biomedical eld. Without diving into details, due to the "out of scope" nature of this ob-
jective, a joint work with the BME department will attempt to conclude if the curvature
detection tool can yield any type of insight at dierentiating cancer networks [15]. Particu-
larly, an attempt at studying the Gene Co-expressions Networks (GCNs) from the curvature
point of view would be performed. In just a few words: \Nodes represent genes and nodes
are connected if the corresponding genes are signicantly co-expressed across appropriately
chosen tissue samples" (taken from Edwin Juarez Rosales's presentation on CAMM Journal
Club, 03/11/2015). It is believed in [15] that cancer networks exhibit higher curvature than
their normal counterparts; thus, another objective of this work would be to shield a dierent
curvature measure to try to assess these type of networks.
69
References
[1] Peter G. Doyle and J. Laurie Snell. Random Walks and Electric Networks.
arXiv:math/0001057 [math.PR] 2000.
[2] R. Banirazi and E. Jonckheere. Geometry of power
ow in negatively curved power
grids. In 49th IEEE Conference on Decision and Control (CDC), pages 6259-6264,
Atlanta, GA, December 15-17 2010.
[3] Chi Wang and E. Jonckheere and R. Banirazi. Interference constrained network per-
formance control based on curvature control. In 2016 American Control Conference
(ACC), pages 6036-6041, Boston, MA, July 6-8, 2016.
[4] M. Lou and E. Jonckheere and Y. Baryshnikov and F. Bonahon and B. Krishnamachari.
Load Balancing by Network Curvature Control. In International Journal of Computers,
Communications and Control (IJCCC), 6 (1), pages 134-149, March, 2011, ISSN 1841-
9836.
[5] T. Jorgensen and J. Pearse. A Hilbert space approach to eective resistance metric. In
Complex Anal. Oper. Theory, pages 975-1030, 2010.
[6] A. ghosh, S. Boyd and A.Saberi. Minimizing eective Resistance of a graph. In SIAM
Review,volume 50, number 1, pages 37-66, Wiley, New York, 2008.
[7] Andrew Pressley. In Elementary Dierential Geometry, Springer International Edition,
2001.
[8] Poonsuk Lohsoonthorn. Hyperbolic Geometry of Networks, USC Thesis, 2008.
[9] Mingji Lou. Trac Pattern in Negatively Curved Networks, USC Thesis, 2008
[10] E. Joncheree, P. Lohsoonthorn and F Bonahon. Scale Gromov Hyperbolic Graphs. In
Journal of Graph Theory, 57: 157-180, 2008, DOI 10.1002/jgt.20275.
70
[11] E. Jonckheree, Mingji Lou, Francis Bonahon and Yuliy Baryshnikov. Euclidean versus
hyperbolic congestion in idealized versus experimental networks. InInternet Mathemat-
ics, 7 (1): 1-27, March 2011.
[12] Joaquin Mur-Amada and Jess Salln-Arasanz. Power
uctuation in a wind farm com-
pared to single turbine. In From turbine to wind farms|Technical requirements and
spin-o products,6, pages 101-132, editor: Gesche Krause, pubisher: Intech, April, 2011,
ISBN 978-953-307-237-1, DOI: 10.5772/15957
[13] E. Grippo and E. Jonckheere. Eective resistance criterion for negative curvature:
application to congestion control. In 2016 IEEE Conference on Control Applications
(CCA), Part of 2016 IEEE Multi-Conference on Systems and Control, pages 129{136,
Buenos Aires, Argentina, Sept. 19-22 2016.
[14] Naomi Ehrich Leonard Gertge Forrest Young, Luca Scardovi. A new notion of eective
resistance for directed graphs - part i: Denition and properties. In IEEE Transactions
on Automatic Control, pages 1727-1736, Jul. 2016.
[15] Yasin Senbabaoglu, Allen Tannenbaum. Graph Curvature for Dierentiating Cancer
Networks. In Scientic Report, doi: 10.1038/srep12323, 2015
[16] CIGRE Working Group B2.43. Guide for Thermal Ratings Calculation of Overhead
Lines. In IEEE Transactions on Power Delivery, December 2014
[17] IEEE Std. IEEE Standard for Calculating the Current-Temperature Relationship of
Bare Overhead Conductors. page 738, 2012
[18] T. Morgan. Eect of Elevated Temperature Operation on the Tensile Strength of Over-
head Conductor. In IEEE Transactions on Power Delivery, Vol. 11, No. 1, pages 345-51,
January 1996.
71
[19] Federal Energy Regulatory Commission and the North American Electric Reliability
Corporation. Arizona-Southern California Outages on September 8, 2011. Causes and
Recommendations. April 2012.
[20] J. Heckenbergerova and J. Hosek. Dynamic thermal rating of power transmission
lines related to wind energy integration. In11th International Conference on Envi-
ronment and Electrical Engineering (EEEIC), Venice, Italy, 18-25 May 2012. (DOI:
10.1109/EEEIC.2012.6221484)
[21] H. Wan and J.D. McCalley. Increasing thermal rating by risk analysis. InIEEE Trans-
actions on Power Systems, volume 14, number 3, pp. 815 - 828, Aug. 1999.
[22] E. Jonckheere, E. Grippo and R. Banirazi. Curvature, Entropy, Congestion Management
and the Power Grid. In3rd IEEE Conference on Control Technology and Applications
(CCTA), pages 535-542 Hong Kong, China, on August 19-21, 2019.
[23] E. Jonckheere and E. Grippo. Ollivier-Ricci Curvature approach to Cost-Eective
Power Grid Congestion Management. In31st Chinese Control and Decision Conference
(CCDC), pages 2118-2123, Nanchang, China, on June 3-5, 2019.
[24] E. Jonckheere, E. Grippo and Reza Banirazi. Congestion Management for Cost-eective
Power Grid Load Balancing using FACTS and Energy Storage Devices allocated via Grid
Curvature Means. InAmerican Control Conference, pages 3909-3915, Philadelphia, PA,
USA on July 10-12, 2019.
[25] Zhongmin Shen. Lectures on Finsler Geometry. World Scientic, May 2001.
[26] M. DoCarmo. Riemannian Geometry. Birkhauser, Boston, 1993.
[27] Y. Ollivier. A visual introduction to Riemannian curvatures and some discrete gener-
alizations. In 50th Seminar of Mathematics Superior vol. 56 of CRM Proceedings and
Lecture Notes, pp. 197-220, 2013.
72
[28] Y. Ollivier. Ricci curvature of Markov chains on metric spaces. In Journal of Functional
Analysis vol. 256, no. 3, pp. 810-864, 2009.
[29] Y. Ollivier. Ricci curvature of metric spaces. In Comptes Rendus Mathematique vol.
345, no. 11, pp. 643-646, 2007.
[30] S. Rachev and L. Ruschendorf. Mass Transportation Problems. In vol. 1 of Probability
and Its Applications, New York: Springer-Verlag, First Ed., 1998.
[31] Romeil Sandhu, Tryphon Georgiou, Ed Reznik, Liangjia Zhu, Ivan Kolesov, Yasin Sen-
babaoglu and Allen Tannenbaum. Graph Curvature for Dierentiating Cancer Net-
works. In Scientic reports. 5. 12323. 10.1038/srep12323, July 14, 2015.
[32] Hamza Farooq, Yongxin Chen, Tryphon T. Georgiou, Allen Tannenbaum and
Christophe Lenglet. Network Curvature as a Hallmark of Brain Structural Connec-
tivity. In bioRxiv 162875; doi: https://doi.org/10.1101/162875 July 12, 2017.
[33] Nurul Idayu Yuso, Abdullah Asuhaimi Mohd Zin, and Azhar Bin Khairuddin. Con-
gestion Management in Power System: A Review. In 2017 3rd International Conference
on Power Generation Systems and Renewable Energy Technologies (PGSRET), Johor
Bahru, Malaysia, April 4-6 2017.
[34] A. Pillay, S. Prabhakar Karthikeyan, and D. P. Kothari. Congestion Management in
Power System: A Review. In International Journal of Electrical Power and Energy
Systems, vol. 70, pp. 83-90, 2015.
[35] K.Vijayakumar. Optimal Location of FACTS Devices for Congestion Management in
Deregulated Power Systems. In International Journal of Computer Applications, Volume
16 No.6, February 2011.
73
[36] Hiren Patel and Ravikumar Paliwal. Congestion management in deregulated power
systems using FACTS devices. In International Journal of Advances in Engineering
and Technology, April 2015.
[37] Vijaya Margaret and K Uma Rao. Demand Response for Residential Loads using Arti-
cial Bee Colony Algorithm to Minimize Energy Cost. In TENCON 2015 - 2015 IEEE
Region 10 Conference, Macao, China, Nov 1-4 2015
[38] Insoon Yang and Asuman E. Ozdaglar. Reducing Electricity Price Volatility via Stochas-
tic Storage Control. In 2016 American Control Conference (ACC), pages 4138-4144,
Boston, MA, USA, July 6-8 2016.
[39] Junjie Qin, Insoon Yang and Ram Rajagopal. Submodularity of Energy Storage Place-
ment in Power Networks. In 2016 IEEE 55th Conference on Decision and Control
(CDC), Las Vegas, NV, USA, Dec. 12- 14 2016.
[40] Y. Rubner, C. Tomasi and L. Guibas. The Earth Mover's Distance as a Metric for
Image Retrieval. In International Journal of Computer Vision, 40(2), page 121, 2000.
[41] Ulrik Brandes and Daniel Fleischer. Centrality measures based on current
ow. In Proc.
22nd Symp. Theoretical Aspects of Computer Sciences, volume 3404/2005 of Lecture
Notes in Computer Science, pages 533{544. Springer Berlin / Heidelberg, 2005.
[42] E. Grippo and E. Jonckheere. Eective resistance criterion for negative curvature:
application to congestion control. In 2016 IEEE Conference on Control Applications
(CCA), Part of 2016 IEEE Multi-Conference on Systems and Control, pages 129{136,
Buenos Aires, Argentina, Sept. 19-22, 2016.
[43] M. Gromov. Asymptotic invariants of innite groups. In M. A. Roller and G. Niblo,
editors, Geometric Group Theory: Proceedings of the Symposium Held in Sussex, 1991.
Cambridge University Press, 1993.
74
[44] Paul Hines and Seth Blumsack. A centrality measure for electrical networks. In HICSS
'08: Proceedings of the Proceedings of the 41st Annual Hawaii International Conference
on System Sciences, page 185, Washington, DC, USA, 2008. IEEE Computer Society.
[45] Zhifang Wang, Anna Scaglione and Robert J. Thomas. Electrical Centrality Measures
for Electric Power Grid Vulnerability Analysis. In 49th IEEE Conference on Decision
and Control , Hilton Atlanta Hotel, Atlanta, GA, USA, December 15-17, 2010.
[46] Ali Pinar, Sandip Roy and Bernie Lesieutre. Power System Extreme Event Detection:
The Vulnerability Frontier. In HICSS '08: Proceedings of the 41st Annual Hawaii
International Conference on System Sciences , page 184, Hawaii, HI, USA, 2008.
[47] Joaquin Mur-Amada and Jess Salln-Arasanz. Power
uctuation in a wind farm com-
pared to single turbine. In Gesche Krause, editor, From turbine to wind farms|
Technical requirements and spin-o products, chapter 6, pages 101{132. Intech, April
2011. ISBN 978-953-307-237-1; DOI: 10.5772/15957.
[48] Y. Ollivier. Ricci curvature on Markov chains on metric spaces. arXiv:0701886v4
[math.PR] 30 Jul 2007, 2007.
[49] L. Shalalfeh, P. Bogdan, and E. Jonckheere. Evidence of long range-dependence in
power grid. In Power and Energy Society General Meeting (PESGM), Boston, MA,
July 2016. Available at http://eudoxus2.usc.edu.
[50] L. Shalalfeh, P. Bogdan, and E. Jonckheere. Kendall's tau of frequency Hurst exponent
as blackout proximity margin. In IEEE International Conference on Smart Grid Com-
munications, Sydney, Australia, November 06-09 2016. 978-1-5090-4075-9/16; available
at http:eudoxus2.usc.edu.
[51] Allen Tannenbaum, Chris Sander, Liangjia Zhu, Romeil Sandhu, Ivan Kolesov, Eduard
Reznik, Yasin Senbabaoglu, Tryphon Georgiou. Ricci curvature and robustness of cancer
networks. http://arxiv.org/abs/1502.04512. 2015.
75
[52] Z. Chen, M. Dehmer, Y. Shi. A note on distance based graph entropies. Entropy 16,
page 54165427, 2014.
[53] Ramin Kazemi. Entropy of weighted graphs with the degree based topological indices
as weights. Communications in Mathematical and in Computer Chemistry, December
2015.
[54] International Council on Large Electric Systems, CIGRE Working Group B2.43. Guide
for Thermal Rating Calculation of Overhead Lines, Technical Brochure 601, CIGRE:
Paris, France, December 2014, ISBN : 978-2-85873-302-6.
[55] IEEE Std. IEEE Standard for Calculating the Current-Temperature Relationship of
Bare Overhead Conductors. page 738, 2012.
[56] T. Morgan. Eect of Elevated Temperature Operation on the Tensile Strength of Over-
head Conductor. In IEEE Transactions on Power Delivery, Vol. 11, No. 1, pages 345-51,
January 1996.
[57] Federal Energy Regulatory Commission and the North American Electric Reliability
Corporation. Arizona-Southern California Outages on September 8, 2011. Causes and
Recommendations. April 2012.
[58] J. Heckenbergerova and J. Hosek. Dynamic thermal rating of power transmission
lines related to wind energy integration. In11th International Conference on Envi-
ronment and Electrical Engineering (EEEIC), Venice, Italy, 18-25 May 2012. (DOI:
10.1109/EEEIC.2012.6221484).
[59] H. Wan and J.D. McCalley. Increasing thermal rating by risk analysis. InIEEE Trans-
actions on Power Systems, volume 14, number 3, pages 815 - 828, Aug. 1999.
76
[60] John Baillieul, Bowen Zhang, and Shuai Wang. The Kirchho-Braess paradox and its
implications for smart microgrids. In54th IEEE CDC, Osaka, Japan, pages 6556-6563,
December 15-18, 2015.
77

Abstract (if available)

Abstract This work is divided in three main parts. The first one develops the asymptotic relation between the resistance of the least resistive path Rₗᵣₚ (i, j) and the effective resistance Reff (i, j) between two extreme points (i, j) of a resistive lattice, as an electrical engineering mean to check negative curvature of a graph. More specifically, the behavior of Rₗᵣₚ (i, j) / Reff (i, j), understood asymptotically in the sense of a growing size network, is able to discriminate a Euclidean grid versus a hyperbolic grid. It is shown that the ratio asymptotically saturates in case of a hyperbolic grid, while it grows without bound in the Euclidean case. In the case of a Gromov hyperbolic grid quasi-isometric to a tree, the asymptotic saturation of Rₗᵣₚ (i, j) / Reff (i, j) is easy to prove, but the major contribution of this first part is to show that the same result holds under the by far less trivial situation of a hyperbolic grid not quasi-isometric to a tree. Applications include the role of curvature in congestion control problems understood in the broad sense as they apply to power grid and data network. ❧ The second part develops a congestion management method for the power grid utilizing the notion of curvature. It initially uses the curvature concept to detect areas prompt to congestion (negative curvature areas) and it subsequently applies load balancing techniques (through FACTS devices) and load (storage devices) deployment to maximize curvature (grid decongestion) and cost-effectively minimize the generated energy throughout the grid, while at the same time guaranteeing stability under phase angle and voltage constraints. Two different curvature definitions are compared: Ollivier-Ricci Curvature (ORC) and Effective Resistance Curvature (ERC), and shown to be locally, globally, resp., relevant. An entropy concept suitable to power grid is also introduced as a new measure to analyze grid congestion. ❧ The third part utilizes the ORC and the ERC to analyze the new power distribution of the grid after tripping of main power grid lines. The main idea is to construct a monitoring/systematic method that identifies the most vulnerable lines and corridors within the power grid that will be at risk after the tripping of a line. Overhead lines are subject to thermal limitations

Linked assets

University of Southern California Dissertations and Theses

Conceptually similar

PDF

Electric vehicle integration into the distribution grid: impact, control and forecast

Asset Metadata

Creator Grippo, Eugenio (author)

Core Title Curvature analysis of the power grid: congestion, congestion management, and line tripping cascade anticipation

School Viterbi School of Engineering

Degree Doctor of Philosophy

Degree Program Electrical Engineering

Publication Date 05/07/2020

Defense Date 05/05/2020

Publisher University of Southern California (original), University of Southern California. Libraries (digital)

Tag cascade effect,complex networks,congestion management,Control,curvature,Gromov,hybrid systems,line tripping,OAI-PMH Harvest,Ollivier-Ricci,power grid,smart grid,vulnerability analysis

Language English

Contributor Electronically uploaded by the author (provenance)

Advisor Jonckheere, Edmond (committee chair), Flashner, Henryk (committee member), Krishnamachari, Bhaskar (committee member)

Creator Email egrippo@usc.edu,leulumlap@gmail.com

Permanent Link (DOI) https://doi.org/10.25549/usctheses-c89-299124

Unique identifier UC11663562

Identifier etd-GrippoEuge-8435.pdf (filename),usctheses-c89-299124 (legacy record id)

Legacy Identifier etd-GrippoEuge-8435.pdf

Dmrecord 299124

Document Type Dissertation

Rights Grippo, Eugenio

Type texts

Source University of Southern California (contributing entity), University of Southern California Dissertations and Theses (collection)

Access Conditions The author retains rights to his/her dissertation, thesis or other graduate work according to U.S. copyright law. Electronic access is being provided by the USC Libraries in agreement with the a...

Repository Name University of Southern California Digital Library

Repository Location USC Digital Library, University of Southern California, University Park Campus MC 2810, 3434 South Grand Avenue, 2nd Floor, Los Angeles, California 90089-2810, USA