Table Of ContentJOURNALOFLATEXCLASSFILES,VOL.13,NO.9,SEPTEMBER2014 1
From State Estimation to Network Reconstruction
Farnaz Basiri, Jose Casadiego, Marc Timme, Dirk Witthaut
Abstract—We develop methods to efficiently reconstruct the is lacking. Manual switches still exist in many distribution
topologyandlineparametersofapowergridfromthemeasure- grids, transmission lines can get lost in contingency cases and
ment of nodal variables. We propose two compressed sensing
targeted attacks can cause entire regions of a grid to collapse.
algorithms that minimize the amount of necessary measurement
We present an example where the connectivity in two regions
resources by exploiting network sparsity, symmetry of connec-
tions and potential prior knowledge about the connectivity. The of the grid is a priori unknown, but can be reconstructed from
algorithmsarereciprocaltoestablishedstateestimationmethods, asinglesnapshotofthenodalvariables.Infuturepowergrids,
where nodal variables are estimated from few measurements the network’s structure will become even more variable: The
7
given the network structure. Hence, they enable an advanced
1 transmission of high-voltage direct current (HVDC) lines in
grid monitoring where both state and structure of a grid are
0 hybridpowergridscanbeactivelycontrolledandtheeffective
subject to uncertainties or missing information.
2 line parameters can be regulated with flexible alternating cur-
IndexTerms—Stateestimation,Powersystemcontrol,SCADA
n rent transmission system (FACTS) devices. Then it becomes
systems, Compressed sensing
a essential to include lacking information or uncertainties of
J
the grid topologies in any grid monitoring tool. The methods
8 I. INTRODUCTION presented in this article shall contribute to the development of
2
such advanced monitoring algorithms [3]–[5].
The secure and reliable operation of complex power grids
] requiresapreciseknowledgeofthegridtopologyandthestate
Y II. THEDCAPPROXIMATION
of connected generation and transmission elements. A central
S The DC approximation describes the flow of real power
tool is ‘state estimation’, i.e the estimation of the current
. flow in AC power grids in a linearized way. It is derived
s state of the nodes from few, possibly noisy measurements
c fromthefullnonlinearloadflowequationsunderthefollowing
and the grid data. Various methods and algorithms have been
[ additional simplifying assumptions [2], [6]–[8]:
developed in this field, see e.g. [1] for a review. A prime
1 example is the estimation of the voltage phase angles from • All bus voltage magnitudes |Vs| are close to 1 p.u. and
v power flow measurements (see e.g. [2] for basic textbook approximated to be exactly |Vs|=1 p.u.
4
examples). State estimation algorithms typically assume that • Voltageangledifferencesacrossbranchesaresmall,such
8
0 the network structure, i.e. the location and parameters of that sin(ϕs−ϕr)≈ϕs−ϕr and cos(ϕs−ϕr)≈1.
9 transmission lines and transformers, are known. However, • Branchescanbeconsideredlossless.Inparticular,branch
0 parts of the information might be lacking or uncertain, for resistances and charging capacitances are negligible.
1. example for contingencies or attacks, such that the question The approximate real power flow from node s to node r over
0 arises: Is it possible to reconstruct the network structure from a transmission line with reactance Xsr is then given by
7 local measurements only? |V ||V |
1 P = s r sin(ϕ −ϕ )≈X−1(ϕ −ϕ ) (1)
: In this article we develop two algorithms that enable the sr Xsr s r sr s r
v faithful reconstruction of the network structure or parts of it
where ϕ is the phase angle at bus s. The real power balance
Xi from the measurement of nodal variables. Special attention at the nosde s reads
is paid to the efficiency of these methods: We show how
r N N
a to reduce the number of measurements using methods from (cid:88) (cid:88)
P = P = X−1(ϕ −ϕ ). (2)
compressed sensing exploiting the structural properties of s sr sr s r
r=1 r=1
power grids. We here focus on the DC approximation to
For notational convenience we summarize these equations
formulate the fundamental ideas of network reconstruction
in matrix form. We define the vector of all voltage angles and
most clearly.
power injections,
The presented algorithms can find several applications in
the monitoring and operation of (smart) power grids. While ϕ=(ϕ ,...,ϕ )T ∈RN,P =(P ,...,P )T ∈RN,
1 N 1 N
it is unlikely that the entire network structure is unknown,
where N is the number of buses or nodes in the grid and the
situations regularly occur where parts of the information
superscript ‘T’ denotes the transpose of a vector or matrix.
The DC approximation then yields the condition
F. Basiri is with the Forschungszentrum Ju¨lich, Institute for Energy and
ClimateResearch-SystemsAnalysisandTechnologyEvaluation(IEK-STE),
Bϕ=P, (3)
52428Ju¨lich.
J.CasadiegoandM.TimmearewiththeMaxPlanckInstituteforDynamics wherethenodalsusceptancematrixB ∈RN×N haselements
andSelf-Organization(MPIDS),37077Go¨ttingen,Germany.
D.WitthautiswiththeForschungszentrumJu¨lich,InstituteforEnergyand (cid:88)N
C52li4m28ateJu¨Rliecshe,arGcherm-SaynsyteamndsAthnealIynssitsitauntedfToerchTnhoeloorgeyticEavlaPluhaytsioicns,(IUEnKiv-SeTrsEit)y, Bnk = j=1Xn−j1 if k =n; (4)
ofCologne,50937Ko¨ln,Germany. −X−1 if k (cid:54)=n.
nk
JOURNALOFLATEXCLASSFILES,VOL.13,NO.9,SEPTEMBER2014 2
ThematrixBisaLaplacianmatrix,whichhasonezeroeigen- IV. NETWORKRECONSTRUCTION
value with eigenvector (1,1,...,1)T [9]. This eigenvector
Let us come back to the defining equation for the DC
represents a global shift of all voltage angles without physical
approximation
significance. To avoid this complication, one often fixes the
Bϕ=P. (9)
voltage angle at a reference bus and excludes this bus from
thecalculation,i.e.thecorrespondingrowandcolumnwillbe In state estimation, we assume that the nodal phase angles
missing in B. We finally note that mathematically equivalent ϕ are difficult to measure and shall be estimated from the
models of flow are used to describe hydraulic networks [10] knowledge of the two other quantities in this equation, B
or vascular networks of plants [11]. and P. But we can also reverse the problem and ask what
happens if we don’t know the nodal susceptance matrix B.
Can we efficiently reconstruct its entries from measurements
III. STATEESTIMATION
of the two remaining quantities P and ϕ? This problem may
State estimation is a mathematical procedure by which the appear hypothetical at a first glance, as the network structure
state of an electric power system is extracted from a set and parameters are generally known. But there are situations
of measurements [1]. In general, any measurement can be where at least parts of the network structure are unknown:
expressed as a function of the system state. Let Switchescanbeopenorclosed,transmissionlinescanundergo
failures, or plans may be inaccurate. Even more, in times
z =h (x) (5)
i i of war or terrorism we may be extremely unsure about the
denote a measured quantity, where x is the system state and physical integrity of parts of the grid. Can we reconstruct the
h is a function specific to the measured quantity z . Assume structure of the grid from measurements only?
i i
that M measurements are taken. Then, all measurements can In the following we develop a mathematical theory of
be written in compact form network reconstruction. To introduce the method, we start
fromthehypotheticalsituationthatB iscompletelyunknown,
z =h(x) (6)
whereas P and ϕ can be measured perfectly. We then study
morerealisticcases,whereonly partsofthegridtopologyare
where
unknown.
x∈RN×1 is the system state, We first notice that a single measurement of P and ϕ will
z ∈RM×1 is a vector of measured quantities and obviously not be enough to reconstruct the entire grid topol-
ogy. However, we may repeat our measurement at different
his a vector function, i.e. a mappingRN →RM
times t ,t ,...,t such that we obtain a large number of
1 2 M
Typically more measurements than the number of state vari- conditions on the entries of B of the form
ables to be determined are taken, i.e. M > N, such that (cid:88)
B ϕ (t )=P (t ) (10)
rk k m r m
equation (6) is overdetermined. It is then solved in a least-
k
squaresfashion,i.e.onecalculatesxsuchthat(cid:107)z−h(x)(cid:107)2 is
2 which hold for all rows r = 1,...,N and all measurement
minimized. For underdetermined equations (6) typically many
timestepsm=1,...,M.Furthermoreweknowthattherow-
solutions exist.
sums of the nodal susceptance matrix (4) vanish,
It is often helpful to work with a simplified DC approx-
(cid:88)
imation model for the measurement equations in analyzing B =0. (11)
rk
the inherent limitations of various methods related solely to
k
the measurement configuration [12]. The state of the grid is
In total we thus have N ×(M +1) conditions which we can
described by the voltage phase angles x = ϕ and shall be
use to reconstruct the unknown entries of B.
estimated from measurements. As explicated above, the real
powerflowmeasuredfrombusstor canbeapproximatedby
A. Row-wise reconstruction
P =X−1(ϕ −ϕ )+e (7)
sr sr s r sr In a first approach we aim to reconstruct the nodal admit-
whereϕ isthephaseangleatbussande isameasurement tance matrix row-by-row (cf. [13]–[16]). To simplify notation
s sr
we collect all measurement conditions (10) and the condition
error. Similarly, a power injection measurement
(11) for a given row r and rewrite them in matrix form
N
(cid:88)
Ps = Psr+e(cid:48)s (8) ϕ1(t1) ··· ϕN(t1) B Pr(t1)
. . r1 .
r=1 .. .. .. = .. (12)
at a given bus s can be expressed as a sum of flows along ϕ1(tM) ··· ϕN(tM)B. Pr(tM)
incident branches to that bus. The measurement vector z 1 ··· 1 rN 0
(cid:124) (cid:123)(cid:122) (cid:125)
thus consists of a subset of the power flows Psr and power (cid:124) (cid:123)(cid:122) (cid:125) =:Br (cid:124) (cid:123)(cid:122) (cid:125)
injections P . In the DC approximation the vector function =:Φ =:Pr
s
h are all linear, which strongly simplifies the problem such Here,B ∈RN×1 denotesthetransposeoftherthrowofthe
r
thatitisespeciallysuitabletointroducenewconceptssuchas nodaladmittancematrixB,whichwewanttoreconstruct.The
network reconstruction. matrix Φ ∈ R(M+1)×N summarizes all measurement results
JOURNALOFLATEXCLASSFILES,VOL.13,NO.9,SEPTEMBER2014 3
for the voltage phase angles at all time steps, while the vector where
P ∈R(M+1)×1 summarizesthevaluesofthepowerinjection ϕc(t1)
r
onlyfornoder.Thelastrowofthematrixequationrepresents Pred =P −(cid:88)N B K ... (17)
the condition (11) for the row sum. r r c=1 rc rcϕc(tM)
Now we have to distinguish, whether equation (12) is
1
(over)determined or underdetermined. If the number of lin-
early independent equations is larger or equal to N we and Φrred is the submatrix of Φr obtained by deleting all
generally have enough information to directly compute Br. columns c for which Krc = 1 and Brred is the submatrix of
In particular there is exactly one solution if rank(Φ¯|P¯r) = Br obtained by deleting all rows r(cid:48) for which Kr,r(cid:48) =1. The
rank(Φ¯) = N, which can be calculated using Gaussian dimension of the reduced linear system of equations (16) is
elimination.Nevertheless,matrixinversioncanbenumerically smallerthanintheoriginalproblem(12),suchthatthecorrect
ill-conditioned, especially for systems with large N [17]. solution can typically be found from fewer measurements.
If the data includes some measurement noise, the system
of equations (12) can be overdetermined when measuring C. Reconstruction from power flow measurements
M +1 > N times such that it is solvable in a least squares
The reconstruction scheme introduced above requires the
fashion, i.e. we have to determine the vector B ∈ RN×1
r knowledge of the voltage phase angle at all nodes of the grid.
minimizing the 2-norm
This data can in principle be obtained using phasor measure-
M+1 ment units (PMUs), but these are typically very expensive.
(cid:88)
(cid:107)ΦBr−Pr(cid:107)22 = (ΦBr−Pr)2k. (13) The measurement of other quantities such as real power flows
k=1 is typically much simpler.
Toefficientlyreconstructthegridtopology,wewanttorelyon Fortunately,weusuallydonothavetoreconstructtheentire
asfewmeasurementsaspossible.IfM+1<N,wegenerally grid in most cases. In a typical application we know a lot
have an underdetermined set of equations which admits many aboutthegridandhavetoreconstructonlypartsofthematrix
possible solutions. Is it possible to obtain the correct solution B. Then the measurement of power injections and real power
also in this case? Do we have more information about B flows of a certain subset of transmission lines can be enough
r
which we can exploit? Indeed, we know that a power grid to perform the reconstruction of the remaining data using
is typically very sparse – a single substation is connected to the methods described above. Technically, this amounts to
onlyfewothersubstations.Thuswechoosetheonesolutionto combining methods of classical state estimation and network
equation(12)whichminimizesthenumberofnon-zeroentries reconstruction.
of B . Unfortunately, the direct minimization is computa- So assume that we can measure the power injection P
r n
tionally hard in general. In 2006 Candes, Romberg, Tao and for all nodes n ∈ {1,...,N} at the time steps t ,...,t .
1 M
Donohoshowedthatanefficientreconstructionisnevertheless In addition we have knowledge about the transmission line
possibleusingaconvexsurrogateforsparsity:the1-norm[18], parameters,inparticularX ,andmeasuretherealpowerflow
sr
[19]. Under weak conditions, the correct sparse solution can P (t )=X−1(cid:0)ϕ (t )−ϕ (t )(cid:1) (18)
be calculated efficiently by minimizing the 1-norm sr m sr r m s m
N for a subset L of all transmission lines at all time steps
(cid:88)
(cid:107)Br(cid:107)1 = |Brk| (14) t1,...,tM. If this set of equations is fully determined or even
k=1 overdeterminedwereconstructthestatevectorϕ(tm)asinthe
subject to the constraint (12). This problem can be mapped to classical state estimation problem described in section III. As
a linear program which we explicate in the following section. before the resulting estimates for the nodal phase angles for
all time steps t ,...,t are then summarized in the matrix
1 M
B. Partial reconstruction Φ. In addition we have partial knowledge about the nodal
susceptance matrix, in particular we already know the entries
Inapplicationswewillrarelyencounterthesituationthatthe
structure of a power grid is completely unknown. In a typical B =B =−X−1 (19)
rs sr sr
applicationweknowthevalueofB ofalltransmissionlines
rs
for all (s,r) ∈ L. This information can then be used to
(s,r) in one part of the grid and have to reconstruct only the
reconstruct the remaining entries of B as described in the
remaining entries of the matrix B – for instance we might
previous section.
want to monitor the position of manual switches at remote
places of the grid.
To keep track of our knowledge about the network, we V. IMPLEMENTATION
define the matrix K ∈RN×N with entries A. Least-squares solution
(cid:26)
1 if the value of B is known for line (r,s) If we have many measurements or prior knowledge avail-
K := rs
rs 0 otherwise. able, then the system of equations (12) can be be overdeter-
(15) mined. In this case equation must be solved in a least squares
We again proceed row by row. For the rth row we reduce fashion, i.e. we have to solve
equation (12) to
ΦredBred =Pred, (16) min(cid:107)ΦBr−Pr(cid:107)22, (20)
r r r Br
JOURNALOFLATEXCLASSFILES,VOL.13,NO.9,SEPTEMBER2014 4
which leads to [15]
Start
B =(ΦTΦ)−1ΦT P (21)
r r
Measurement of Φ and 𝑃
providedthatΦisfullranksuchthattheinverseexists.Least-
squares solutions are already implemented in many numerical
Add constraint ∑ 𝐵 =0
solvers, for instance in the MATLAB function mldivide. 𝑐 𝑟𝑐
B. Minimizing the 1-norm no
𝑟𝑘(Φ)<𝑁 ?
The key to an efficient reconstruction of B from an un-
derdetermined system of equations is the minimization of the Linear System
yes underdetermined
1-norm [18], [19]. This problem can be mapped to a linear
underdetermined
program which can be solved efficiently, i.e. the computation Calculate null space 𝑊of Φ
time scales at most polynomially with the input length.
So we consider the linear system of equations (12) and 𝑟 = 1
assume that it is solvable and underdetermined. The solutions
span an affine subspace of of RN with dimension D = Calculate Specific solution 𝐵𝑟(sp)
N −rank(Φ). All solutions can be written as
Solve linear program
Br =B(rsp)+Wy, (22) min ‖𝑊𝑦+𝐵(sp)‖
𝑦 𝑟
1
where the columns of the matrix W ∈ RN×D form a basis
for the kernel (nullspace) of Φ and y ∈ RD is a vector of Set 𝐵(𝑟,:)=[𝑊𝑦+𝐵(sp)]𝑇 𝐵=Φ\𝑃𝑟
parameters. B(sp) is a specific solution to equation (12); in 𝑟
r
the implementation it is obtained using the MATLAB function 𝑟 = 𝑟+1
mldivide. We now search for the vector y such that the 1-
yes
norm (14) assumes its minimum. This optimization problem 𝑟≤𝑁
is rewritten as a linear program [20]
no
B(sp)+Wy (cid:22)s
min1Ts such that r , (23) Return 𝐵
s B(sp)+Wy (cid:23)−s
r
where s ∈ RN is an auxiliary variable, 1 ∈ RN is a vector
Fig.1. Flowchartforrow-wisereconstruction.MATLABcodeforthecentral
of ones and (cid:22) and (cid:23) denotes entry-wise comparison. In the partislistedinFig.2.
standard form of many commercial solvers the optimization
problem reads
If the system is underdetermined we proceed row-by-row and
minfTx such that Ax(cid:22)b. (24)
reconstruct B as described in section V-B using the function
x r
null to compute the matrix W and the function linprog
using the further auxiliary variables
to solve the linear program (24).
(cid:18) (cid:19) (cid:18) (cid:19)
y W −1l
x= , A= ,
s −W −1l D. Iterative reconstruction
b=(cid:18) B(rsp) (cid:19), f =(cid:18)0(cid:19), In practice we typically have some prior knowledge about
−B(sp) 1 the entries of the matrix B as outlined in section IV-B. An
r
algorithm to perform the reconstruction exploiting the prior
where 1l ∈ RN×N is the identity matrix and 0 ∈ RD is a
knowledge is illustrated by the flow-chart in Fig. 3
vector of zeros.
Mostinterestingly,wecanalsousethisalgorithmtogreatly
improve the convergence of the reconstruction algorithm. For
C. Row-Wise Reconstruction
intermediate values of M we typically face the situation that
Wefirstdescribetheimplementationoftherow-wiserecon- some rows are successfully reconstructed while others are
struction of B from the measurement data without any prior not. If we have not reconstructed B successfully, but gained
knowledge about the entries of B, whose basic idea has been further knowledge, we may use it in the following. Having
outlined in section IV-A. The algorithm to solve this problem successfullyreconstructedtherthrowofB,wealsoknowthe
is illustrated by the flow-chart in Fig. 1. We implement this entriesB =B forallotherrowsr(cid:48) duetothesymmetry
r(cid:48),r r,r(cid:48)
algorithm in Matlab, the central part of the program code is of the matrix. We thus propose to perform the reconstruction
listed in Fig. 2. iteratively. In each step of the reconstruction process we use
First data is collected and the matrix Φ and the vectors the initial knowledge about the entries of the matrix B and
P are formed. Then the problem is solved depending on the also the knowledge gained in previous steps.
r
rank of the matrix Φ. If the system is (over)determined, we The algorithm shown in Fig. 3 starts with the input of the
directly solve it for B using the Matlab-function mldivide. measurement results and the prior knowledge of the system
JOURNALOFLATEXCLASSFILES,VOL.13,NO.9,SEPTEMBER2014 5
% Case 1: Overdetermined system: Sta rt
if rank(Phi) >= N
Input measurement results Φ,𝑃
Breconstruct = Phi \ Pt; 𝑟
and prior knowledge 𝐵, 𝐾new
% Case 2: Underdetermined system: Set 𝐾=𝐾new,𝑟=0
else
Proceed to next row 𝑟=𝑟+1
% Calculate nullspace of Phi
W = null(Phi); Set Φ𝑟red , 𝑃𝑟red, 𝑁𝑟red=𝑁−Σ𝑐𝐾𝑟𝑐
D = size(W,2); yes
𝑁red=0 ? Linear System
𝑟 underdetermined
no
% treat all rows separately: underdetermined
no
for r=1:N 𝑟(Φred)<𝑁red
yes
% one specific solution Solve linear program 𝐵red=Φred\𝑃red
min ‖𝐵red‖ s.t. Φred 𝐵red=𝑃red 𝑟 𝑟 𝑟
Br_sp = Phi \ Pt(:,r); 𝑟 1 𝑟 𝑟
no
Successful?
% Choose sparse solution
f = [zeros(D,1); ones(N,1)];
A = [W, -eye(N); - W, -eye(N)]; Set 𝐾𝑟n𝑐ew=1 ∀𝑐=1,…𝑁
Add results 𝐵red to matrix 𝐵
b = [-Br_sp; +Br_sp]; 𝑟
y = linprog(f,A,b);
no
y = y(1:D); 𝑟=𝑁?
yes
yes
% Write result to matrix B 𝐾new=𝐾 ?
Breconstruct(r,:) = (W*y+Br_sp)’; no Output
no 𝐾new=1 ∀ 𝑟,𝑐n o? yes Warning
𝑟𝑐
end
end Return 𝐵
Fig.2. MATLABcodeforrow-wisereconstruction Fig.3. Flowchartfortheiterativereconstructionalgorithm.
substations.Henceweseekforasolutionwhichismostsparse,
which is encoded in the matrices B and K. It solves the
i.e. which minimizes
reconstruction problem iteratively, processing row-by-row in
each iteration. For each row, we first calculate the reduced (cid:107)Br(cid:107)0 =number of non-zero entries ofBr. (25)
matrices Φred and Pred as defined in section IV-B. Then we
r r However, we cannot minimize (cid:107)B (cid:107) directly as this is a
attempt to reconstruct the vector Bred as described above: If r 0
r computationally hard problem. The ingenious contribution of
the problem is (over)determined we solve the linear system
Candes, Romberg, Tao and Donoho [18], [19] was the proof
of equations using the Matlab-function mldivide, if it is
that the most sparse solution can be obtained by minimiz-
underdetermined we minimize the 1-norm (cid:107)Bred(cid:107) subject to
r 1 ing (cid:107)B (cid:107) with high probability if a sufficient amount of
the constraint (16). If the reconstruction was successful, we r 1
measurement data is available. Otherwise we find a solution
add the results to the matrix B and update the knowledge
which minimizes the 1-norm, but has a large number of non-
matrix K. To facilitate the bookkeeping we define two ma-
zero entries (cid:107)B (cid:107) . In algorithm in Fig. 3 we thus adopt the
trices K and Knew: We set K = Knew at the beginning r 0
definition that the reconstruction of the rth row is assumed
of each iteration and only modify Knew during the step. The
to be successful if (a) the linear system is overdetermined
iteration stops when the reconstruction has been successfully
or (b) the reconstructed row B is sufficiently sparse, i.e.
completed, i.e. Knew = 1 for all r,c = 1...,N, or when r
rc (cid:107)B (cid:107) ≤ d for some upper limit d . A typical value
no further progress has been made, i.e. Knew has not been r 0 max max
used in the following is d =15.
modified during the last step. max
One unsolved problem remains: If we do not know the
VI. APPLICATIONSANDPERFORMANCE
matrix B a priori, how can we know if the reconstruction of
A. Reconstruction from time series
the rth row has been successful? In the underdetermined case
we cannot decide whether the reconstructed values B are We demonstrate the applicability of network reconstruction
rc
definitelycorrect–butwecandecideiftheyarereasonablein from time series measurements for a test grid taken from [21]
termsoftheconnectivityofthegrid.Powergridsaregenerally illustrated in Figure 4 (a). This data set includes hourly data
verysparse:Asinglesubstationisconnectedtoonlyfewother for demand and generation P (t) for one year and all N =
k
JOURNALOFLATEXCLASSFILES,VOL.13,NO.9,SEPTEMBER2014 6
(a) 300
(a)
stn 200
e
m
e
ru 100
s
a
e
m
fo re 0100 101 102 103
(b) P(t)km[MW] 120000000 −−2100000 bmun m 60 (b)
φ(t)kmdegree] −1100 −1100 uminim 40
[ 0 20 40 60 0 20 40 60 20
time t [hour] time t [hour]
(c) m m
400
0
100 101 102 103
s 300 grid size N
or
err
of 200 Fig. 5. Performance of (a) the row-wise reconstruction algorithm shown
er inFig.1and(b)theiterativealgorithmshowninFig.3.Wehavetestedthe
b
m algorithmsforseveraltestgridsfrom[22],[23]forsyntheticpowerinjection
nu 100 data,choosingeithertheϕk (bluecircles)orthePk (redcrosses)randomly
(see text for details). Shown is the minimum number of measurements M
0 needed to obtain a faithful reconstruction as a function of the grid size N.
20 25 30 35 40 45 50 55 60 Resultshavebeenaveragedover10randomrealizations.
number of measurements M
Fig. 4. Example of the complete grid reconstruction for a test grid with exploitthesparsityofthematrixB asanadditionalstructural
N = 117 nodes. (a) Test grid from [21]. It is assumed that the network
information.
topology,i.e.thelocationandparametersofthetransmissionlines(grey)are
completely unknown. (b) The reconstruction is based on the time series of
thenodalpowerinjectionsPk(t)andvoltageanglesϕk(t).Twotimeseries B. Minimum Measurement Requirements and Scaling
areshown.(c)PerformanceoftheiterativealgorithmshowninFig.3.Perfect
reconstructionisobtainedwithM ≥55measurements.Weusethetestgrid The algorithms presented in figures 1 and 3 allow for a
andtheloadandgenerationtimeseriesofthescenario2013from[21].For
reconstructionofthenetworktopologyalsointheunderdeter-
thesakeofsimplicitywehaveremovedtheimport/exportnodesandrescaled
thegenerationtoexactlymatchtheload. mined case. But how many measurements M are needed
min
for a faithful reconstruction and how does this number scale
with the grid size N?
117 nodes. The nodal voltage angles ϕk(t) are then obtained To systematically study the efficiency of the reconstruction
by solving the DC approximation (3) for all time steps. The algorithms, we consider various test grids of different size
resulting time series of the power injections and angles are taken from [22], [23] and use randomized synthetic data for
shown in Figure 4 (b) for two nodes as an example. the time series. We consider two different types of time series
The algorithm shown in Figure 3 can now reconstruct data. First we draw the voltage phase angles ϕ (t) uniformly
k
the entire network topology - i.e. all entries of the nodal at random from the interval [−π/8,+π/8] (type I). This
susceptance matrix B - from the time series data. To evaluate scenario is not realistic but close to the original mathematical
theperformanceofthealgorithmwevarytheamountofinput work[18],[19],whererandomGaussiansamplingvectorsare
data.ForeachvalueofM werunthealgorithmandcountthe considered. Second, we choose the power injections P (t)
k
number of successfully reconstructed entries and the number at random from a normal distribution with mean zero and
ofreconstructionerrors.Toaccountforsmallnumericalerrors, standard deviation as in the original test grid (type II). For
we say that an entry (r,k) of the nodal susceptance matrix is both types and both reconstruction algorithms (row-wise and
successfully reconstructed if iterative), we run the reconstruction algorithm as a function
|Breconstructed−Btrue|<(cid:15), (26) of the number of measurements M taken into account and
r,k r,k determine the minimum number of measurements M for
min
where the numerical tolerance is chosen as (cid:15) = 10−3 × a faithful reconstruction, i.e. zero reconstruction errors. Each
max |B |. numerical experiment is repeated 10 times.
ij ij
Figure 4 (c) shows that M = 55 measurements are First, we observe that the iterative algorithm allows for a
min
sufficient to faithfully reconstruct the entire nodal suscep- faithfulreconstructionfromunderdeterminedinputdataforall
tance matrix B. Most importantly, a faithful reconstruction is grids under consideration and both data types. The minimum
possible in the strongly underdetermined case: The minimum number of measurements M is much smaller than the grid
min
number of measurements M = 55 is less than half of the sizeN inallcases.Therow-wisealgorithmperformslesswell
min
matrix dimension N = 117. This is possible because we can as it does not exploit the symmetry of the matrix B.
JOURNALOFLATEXCLASSFILES,VOL.13,NO.9,SEPTEMBER2014 7
(a) pfloowwer (b) 1 prior knowledge grid topology is encoded in the matrix K, which is illustrated
[MW] in figure 6 (b). White entries indicate where we do not know
37
? 29.6 mber 10 the entries of the matrix B.
22.2 bus nu 20 staTteheestmimisastiinogniannfdornmeatwtioonrkisrecfoounnstdrubcytiocno.mFbirisnti,ntgheclvaoslstiacgael
14.8
phase angles for all nodes of the network are reconstructed as
?
7.4 30 describedinsectionIII,theresultsbeingshowninfigure6(c).
1 10 20 30
0 bus number Second, the iterative reconstruction algorithm shown in figure
(c) phase (d) reconstructed network 3 is applied using the measured power injections, the prior
angles 1 knowledge B and K as well as the estimated voltage phase
φ [°]
k
6 anglesasinputdata.Thealgorithmthenfaithfullyreconstructs
4.8 mber 10 all the missing information about the grid topology. The
3.6 bus nu 20 rfiegcuornest6ru(cdt)edexvaacltuleysmoaftcthhethlienetrusuesvceaplutaens.ceMsoXsts−ri1mpshoortwanntliyn,
2.4
the full reconstruction is possible already for M = 1 in this
1.2 30
1 10 20 30 case, i.e. with a single snapshot of the power injections P,
0 bus number
making use of the prior knowledge.
Fig.6. Reconstructionofmissinginformationfromflowmeasurements(a)
The IEEE 30-bus test grid. Real power flows are measured on the colored VII. CONCLUSIONANDOUTLOOK
lines. In two regions the grid topology is assumed to be unknown. (b) The
prior knowledge matrix K. Grey shading indicates that the respective line Wehaveintroducedtwoalgorithmstoreconstructthestruc-
parameterBjk isknown.Thewhiteregioniscompletelyunknown.(c)The ture of a power grid from nodal measurements only. One
voltagephaseanglesϕkarereconstructedfromthepowerflowmeasurements algorithmisstraightforwardandrow-based,reconstructingthe
usingstandardstateestimationmethods.(d)Thecompletelygridtopologyis
perfectlyreconstructedusingthealgorithmshowninFig.3. lines of the grid node by node. The second is iterative and
takesintoaccountbothpriorknowledgeaboutthepresenceof
absence of lines as well as the knowlegde generated by that
Second, efficient reconstruction is possible with both algo- algorithm during previous steps. We have demonstrated how
rithms for the type I input data. The theory of compressed these algorithmscan be usedto reconstruct theentire network
sensing shows that a reconstruction of sparse vectors is pos- structure from time series or missing information on the grid
sible with high probability from a number of measurements topology from a single snapshot.
scaling only logarithmically with the problem dimension N The presented algorithms exploit several structural prop-
if some conditions are satisfied [18], [19]. This logarithmic erties of power grids to reduce the number of necessary
scaling is well confirmed by the results of our numerical measurements. Power grids are typically very sparse, i.e.
experiments. For the iterative algorithm the scaling seems to each substation is connected to only few other substations.
be even slower. Hence, methods from compressed sensing can be used which
However, the proofs for the efficiency of compressed sens- allowforafaithfulreconstructionalsointheunderdetermined
ing depend on two features of the sampling vectors: isotropy case. The iterative reconstruction algorithm makes use of the
and incoherence. Roughly speaking, isotropy means that all symmetry of the nodal susceptance matrix. Loosely speaking
directions in the vector space are sampled equally well. the algorithm solves the simple parts of the reconstruction
Incoherenceguaranteesthattheredundancyintheinformation problem first and than uses the gained information for the
gained by two measurements is small. Both features are no remaining parts. We have shown that this trick leads to a vast
longer guaranteed for type II input data. In fact, we observe reduction of the measurement resources.
that the row-wise reconstruction algorithm does not perform The present work presents a step towards hybrid state
well in this case. The minimum number of measurements estimation/network reconstruction algorithms. In future smart
M is only slightly smaller than the problem dimension N grids one can face the situation where both nodal variables
min
such that the iterative algorithm performs much better. and the grid topology are subject to uncertainties or lacking
information. A pure state estimation starting from perfect
knowledgeofthenetworkstructureisthennolongersufficient.
C. Reconstruction from flow measurements
Advanced grid monitoring algorithms must be able of coping
The iterative algorithm allows to take into account prior with missing information in both state and structure.
knowledge to efficiently reconstruct the missing information
about the grid topology. Such a problem can arise in practice ACKNOWLEDGEMENTS
whenalargerpartofthegridissubjecttodamagesorattacks. We thank D. Gross and R. Kueng for inspiring discus-
Anexampleofsuchasituationisdepictedinfigure6(a).Itis sions.WegratefullyacknowledgesupportfromtheHelmholtz
assumed that we have no reliable information about the grid Association (via the joint initiative “Energy System 2050
in two areas indicated by the question marks. However, we - A Contribution of the Research Field Energy” and the
do have information about the connections of all remaining grant no. VH-NG-1025 to D.W.) and the Federal Ministry of
nodes and we assume that we can measure the real power Education and Research (BMBF grant nos. 03SF0472B and
flowalongthecoloredsolidlines.Thepriorknowledgeofthe 03SF0472E to D.W. and M.T.).
JOURNALOFLATEXCLASSFILES,VOL.13,NO.9,SEPTEMBER2014 8
REFERENCES JoseCasadiegoreceivedhisLicentiatedegreeinPhysicsfromtheUniversity
of Carabobo (Venezuela) in 2009. He completed his PhD studies in the
[1] A. Meliopoulos, B. Fardanesh, and S. Zelingher, “Power system state
International Max Planck Research School for Physics of Biological and
estimation:modelingerroreffectsandimpactonsystemoperation,”in
Complexsystems(Go¨ttingen)in2016workingonmodel-freeapproachesfor
Proceedings of the 34th Annual Hawaii International Conference on
revealinginteractionsincomplexnetworksfromcollectivedynamics.Henow
SystemSciences(HICSS-34). IEEE,2001,p.2012.
works as a Postdoctoral Researcher at the Network Dynamics Group of the
[2] A. J. Wood, B. F. Wollenberg, and G. B. Sheble´, Power Generation,
MaxPlanckInstituteforDynamicsandSelf-Organization.
OperationandControl. NewYork:JohnWiley&Sons,2014.
[3] Y.Yuan,G.-B.Stan,S.Warnick,andJ.Goncalves,“Robustdynamical
networkstructurereconstruction,”Automatica,vol.47,no.6,pp.1230
–1235,2011.
[4] J. Yu, G. G. Karady, and L. Gu, “Applications of embedded hvdc in
powersystemtransmission,”in2012PowerEngineeringandAutomation
Conference,2012,pp.1–6.
[5] S. Bolognani, R. Carli, and M. Todescato, “State estimation in power
distributionnetworkswithpoorlysynchronizedmeasurements,”in53rd
IEEEConferenceonDecisionandControl,2014,pp.2579–2584.
[6] J.J.GraingerandW.D.StevensonJr.,PowerSystemAnalysis. New
York:McGraw-Hill,1994.
[7] K.Purchala,L.Meeus,D.V.Dommelen,andR.Belmans,“Usefulness
of dc power flow for active power flow analysis,” in IEEE Power
Engineering Society General Meeting, 2005, 2005, pp. 454–459 Vol.
1.
[8] D. Van Hertem, J. Verboomen, K. Purchala, R. Belmans, and W. L.
Kling, “Usefulness of DC power flow for active power flow analysis
withflowcontrollingdevices,”in8thIEEInternationalConferenceon
ACDCPowerTransmission(ACDC2006). IEE,2006,pp.58–62.
Marc Timme studied physics and applied mathematics at the Universities
[9] M.Newman,Networks:AnIntroduction. NewYork:OxfordUniversity
of Wu¨rzburg (Germany), Stony Brook (New York, USA) and Go¨ttingen
Press,2010.
(Germany).HeholdsanM.A.inPhysicsfromStonyBrookandaDoctorate
[10] N.HwangandR.Houghtalen,FundamentalsofHydraulicEngineering
in Theoretical Physics (Go¨ttingen). After postdoctoral and visiting stays at
Systems. UpperSaddleRiver,NJ:PrenticeHall,1996.
Cornell University (New York, USA) and the National Research Center of
[11] E.Katifori,G.J.Szo¨llo¨si,andM.O.Magnasco,“DamageandFluctu-
Italy (Sesto Fiorentino) he is heading the Max Planck Research Group on
ationsInduceLoopsinOptimalTransportNetworks,”Phys.Rev.Lett.,
Network Dynamics at the Max Planck Institute for Dynamics and Self-
vol.104,p.48704,2010.
Organization.HeisAdjunctProfessorattheUniversityofGo¨ttingen.
[12] A.AburandA.G.Exposito,Powersystemstateestimation:theoryand
implementation. BocaRaton:CRCPress,2004.
[13] M. Timme, “Revealing Network Connectivity from Response Dynam-
ics,”Phys.Rev.Lett.,vol.98,no.22,p.224101,2007.
[14] S. G. Shandilya and M. Timme, “Inferring network topology from
complexdynamics,”NewJ.Phys.,vol.13,no.1,p.013004,2011.
[15] M. Timme and J. Casadiego, “Revealing networks from dynamics: an
introduction,”J.Phys.AMath.Theor.,vol.47,no.34,p.343001,2014.
[16] X. Han, Z. Shen, W.-X. Wang, and Z. Di, “Robust Reconstruction of
ComplexNetworksfromSparseData,”Phys.Rev.Lett.,vol.114,no.2,
p.028701,2015.
[17] J. Stoer and R. Bulirsch, Introduction to Numerical Analysis. New
York:Springer,1993.
[18] D.L.Donoho,“Compressedsensing,”Inf.Theory,IEEETrans.,vol.52,
no.4,pp.1289–1306,2006.
[19] E.J.Candes,J.K.Romberg,andT.Tao,“Stablesignalrecoveryfrom
incomplete and inaccurate measurements,” Comm. Pure Appl. Math.,
vol.59,no.8,pp.1207–1223,2006.
[20] S. Boyd and L. Vandenberghe, Convex optimization. Cambridge:
CambridgeUniversityPress,2004.
[21] H.Barrios,A.Roehder,H.Natemeyer,andA.Schnettler,“Abenchmark Dirk Witthaut received his Diploma (MSc) and PhD in Physics from the
case for network expansion methods,” in Proc. IEEE PowerTech 2015 TechnicalUniversityofKaiserslautern,Germany,in2004and2007,respec-
Eindhoven,2015,pp.1–6. tively. He has been a Postodoctoral Researcher at the Niels Bohr Institute
[22] R. D. Zimmerman, C. E. Murillo-Sanchez, and R. J. Thomas, “MAT- inCopenhagen,DenmarkandattheMaxPlanckInstituteforDynamicsand
POWER: Steady-State Operations, Planning and Analysis Tools for Self-OrganizationinGo¨ttingen,GermanyandaGuestLecturerattheKigali
Power Systems Research and Education,” IEEE Trans. Power Syst., Institue for Science and Technology in Rwanda. Since 2014 he is heading
vol.26,p.12,2011. a Research Group at Forschungszentrum Ju¨lich, Germany and he is Junior
[23] S. Fliscounakis, P. Panciatici, F. Capitanescu, and L. Wehenkel, “Con- ProfessorattheUniversityofCologne,Germany.
tingencyrankingwithrespecttooverloadsinverylargepowersystems
takingintoaccountuncertainty,preventiveandcorrectiveactions,”IEEE
Trans.PowerSyst.,vol.28,p.4909,2013.
FarnazBasiriFarnazBasirireceivedtheB.Sc.andM.Sc.degreeinelectrical
powerengineeringfromtheRWTHAachenUniversity,Aachen,Germanyin
2013 and 2015, respectively. Since 2015 she joined the Reaserch Group at
Reaserch Center Ju¨lich, Germany as a Research Associate and is currently
pursuingthePh.D.degree.
