Inverse Power Flow Problem
This paper formulates an inverse power flow problem which is to infer a nodal admittance matrix (hence the network structure of a power system) from voltage and current phasors measured at a number of buses. We show that the admittance matrix can be …
Authors: Ye Yuan, Steven Low, Omid Ardakanian
1 In v erse Po wer Flo w Problem Y e Y uan, Senior Member , IEEE, Ste ven H. Low , F ellow , IEEE, Omid Ardakanian, Member , IEEE, Claire J. T omlin, F ellow , IEEE Abstract —This paper f ormulates an in verse power flo w prob- lem which is to infer a nodal admittance matrix (hence the network structur e of a power system) from voltage and current phasors measured at a number of b uses. W e show that the admittance matrix can be uniquely identified from a sequence of measurements corresponding to different steady states when every node in the system is equipped with a measurement device, and a Kron-reduced admittance matrix can be determined even if some nodes in the system are not monitored (hidden nodes). Furthermore, we propose effective algorithms based on graph theory to unco ver the actual admittance matrix of radial systems with hidden nodes. W e provide theoretical guarantees for the reco vered admittance matrix and demonstrate that the actual admittance matrix can be fully recover ed even from the Kron- reduced admittance matrix under some mild assumptions. Sim- ulations on standard test systems confirm that these algorithms are capable of providing accurate estimates of the admittance matrix from noisy sensor data. Index T erms —Inv erse Power Flow Problem, System Identifi- cation, Kron Reduction, Phasor Measurement Units. I . I N T RO D U C T I O N T HE power industry has witnessed profound changes in recent years. These changes are mostly driven by the widespread adoption of distributed energy resources (DER), activ e participation of customers in emerging energy markets, and rapid deployment of measurement, communication, and control infrastructure resulting in an unprecedented level of visibility and controllability , especially for distribution grids. Despite the increased amount of uncertainty , these changes offer opportunities for system operators to improve power system stability and efficienc y by lev eraging advanced op- timization and control techniques. Most of these techniques require the knowledge of the network topology in real time. The in verse power flo w (IPF) problem we define in this paper concerns the estimation of the nodal admittance matrix from synchronized measurements of voltage and current pha- sors (i.e., magnitudes and phase angles) which can be obtained from phasor measurement units (PMUs) or con ventional super- visory control and data acquisition (SCADA) technology . The IPF problem underlies sev eral crucial smart grid applications, affecting real-time system operation and long-term planning, the most important of which are: Y e Y uan is with School of Artificial Intelligence and Automation, Huazhong University of Science and T echnology , Wuhan, China. Ste ven H. Low is with Department of Computer and Mathematical Sciences and Electrical Engineering, California Institute of T echnology , Pasadena, United States. Omid Ardakanian is with Department of Computing Science, Uni- versity of Alberta, Edmonton, Canada. Claire J. T omlin is with Department of Electrical Engineering and Computer Sciences, Univ ersity of California, Berkeley , United States. Correspondence to slow@caltech.edu . This v ersion corrects the proof to Proposition 3 and adds details about constructing the admittance matrix of the overall network from its components. 1) State Estimation [ 1 ] combines the knowledge of the admittance matrix with a set of known state-variables to determine the unknown ones, e.g., voltage magnitude and phase angle of unobserved buses, thereby building a real-time model of the network for control. 2) Optimization & Contr ol [ 2 ] techniques determine a sequence of operations that can transition the power system from one steady state to another steady state that meets certain stability and efficiency targets. 3) Event Detection [ 3 ] concerns identifying faults, line outages, and other critical ev ents, such as transformer tap changes, capacitor and switching operations from changes in the real-time network model. 4) Cybersecurity [ 4 ] is the problem of identifying the potential vulnerabilities of a power system and designing strategies to protect it from the potential cyber attacks using telemetry data along with information about its topology . In this paper , we lay out a theoretical framework for the IPF problem. Using the bus injection model (BIM), we propose efficient algorithms to identify the admittance matrix. In par- ticular , we show that when the system has no hidden nodes, the admittance matrix can be uniquely identified from a sequence of complex voltage and current measurements corresponding to different steady states. Should there be some hidden nodes in the network, we show that a reduced admittance matrix (from Kron reduction [ 5 ]) can be determined; we develop a method based on graph decomposition, maximal clique search- ing and composition for identifying the admittance matrix of the original system for radial networks. Po wer flo w simulations are performed on the IEEE 14-bus system to illustrate the theoretical results and e valuate their sensiti vity to measurement noise introduced by transducers. The paper is outlined as follo ws: after surveying related work in Section II , we formulate the IPF problem and propose a solution for the case that the system is fully observable in Section III . When the system has hidden nodes, we propose ef- ficient algorithms to solve the IPF problem for radial networks with theoretical guarantee in Section IV . W e ev aluate the identification accurac y of the proposed algorithms in Section V and conclude the paper by presenting directions for future work in Section VI . I I . R E L A T E D W O R K The av ailability of high-precision, high-sample-rate mea- surements of transmission and distribution networks in recent years has given impetus to research on topology and admit- tance identification. The IPF problem that identifies both topology and admit- tance matrix has been studied extensiv ely in transmission 2 networks [ 6 ], [ 7 ], [ 8 ], [ 9 ] as well as distribution grids [ 10 ], [ 11 ] using single-phase a.c. and d.c. power flow models. For example, the topology identification problem is cast as a sparse subspace learning problem in [ 6 ] and an efficient algorithm is proposed to estimate the admittance matrix of the underlying power system from the measured power injection of different buses. In [ 12 ], the topology of an urban (mesh) distribution network is inferred from voltage magnitude, and real/reactiv e power measurements carried out by smart meters at the end- nodes. A graphical model is built to describe the probabilistic relationships between dif ferent voltage measurements using lasso. In [ 10 ] a graphical learning based approach is dev el- oped to estimate the radial grid topology from nodal voltage measurements. The learning algorithm is based on condi- tional independence tests for continuous variables o ver chordal graphs. An efficient algorithm for topology identification of a power system is also proposed in [ 8 ] drawing on ideas from compressiv e sensing and graph theory . The authors assume that power and phase angle measurements are av ailable for all nodes. Different algorithms hav e been developed in the literature to make this inference, using v arious techniques such as weighted least square, maximum-likelihood/maximum-a-posteriori esti- mation, minimum spanning tree, sparse recov ery , Lasso/Group Lasso, blind identification, quickest change detection theory , as well as graphical model learning. There are three limitations in the current literature that we propose ov ercome. First, most of the literature focuses on topology identification or change detection, but there is not much work on joint topology and parameter identification, with notable exceptions of [ 6 ], [ 13 ], [ 14 ]. Second, most papers require measurements at every node in the network, with the exceptions of [ 15 ], [ 14 ], [ 16 ], [ 17 ], [ 18 ], [ 19 ]. In particular , [ 14 ] learns the topology with parameters from a stochastic perspecti ve, the true topology can only be found in probability , even when the number of samples is lar ge; [ 17 ], [ 18 ] assume that perturbed data are av ailable (therefore a special inv erter is assumed) to identify the network, which could be strong in practice; [ 19 ] proposed a method based on recursive grouping to estimate the topology and branch impedance for networks that may ha ve hidden nodes, howe ver , without guarantee. I I I . I P F W I T H O U T H I D D E N N O D E S In this section we study the IPF problem when voltage and current phasor measurements are available at every bus in the system. W e formulate the identification problem as a constrained least squares problem and then con vert it to an equiv alent unconstrained least squares problem. W e note that Y has a certain structure that can be exploited when solving the IPF problem— (a) Y is a symmetric but not Hermitian complex matrix (i.e., Y ∈ S N ) and (b) Y encodes the topology of a connected graph (or a connected tree for radial networks). A. Pr oblem formulation Let C denote the set of complex numbers, R the set of real numbers, and N the set of integers. For A ∈ C n × n , Re ( A ) and Im ( A ) denote matrices with the real and imaginary parts of A , respecti vely . Let S n ⊆ C n × n be the set of all n × n complex symmetric (not necessarily Hermitian) matrices. The transpose of a matrix A is denoted A T and its Hermitian (complex conjugate) transpose is denoted A H . A [ i, j ] denotes the element of A located at ith ro w and jth column. W e define I as the identity matrix with an appropriate dimension and define 1 as an all- 1 column vector with an appropriate dimension. A po wer system can be modeled by an undirected connected graph G = ( N , E ) where N := { 1 , 2 , . . . , N } represents the set of buses, and E ⊆ N × N represents the set of lines, each connecting two distinct buses. A bus j ∈ N can be a load bus, a generator bus, or a swing bus. Let V j be the complex voltage at bus j and s j be the net complex po wer injection (generation minus load) at that bus. W e use s j to denote both the complex number p j + i q j ( i ≜ √ − 1 ) and the real pair ( p j , q j ) depending on the context. For each line ( i, j ) ∈ E , we denote its series admittance by y ij . The bus admittance matrix of this system is denoted Y , which is an N × N complex- valued matrix whose off-diagonal elements are Y ij = − y ij and diagonal elements are Y ii = − P j = i Y ij , assuming that there is no shunt element (this assumption can be relaxed). Hence, the current injection vector can be expressed as I = Y V . The formulated IPF problem is: giv en voltage and current measurements of different steady-states, i.e., V i ( k ) and I i ( k ) for k = 1 , . . . , K and i ∈ M = { 1 , . . . , N } , how to recover the true admittance matrix Y . Specially , when M is a subset of N , what part of the true admittance matrix Y can be recov ered under what condition. B. Identification algorithm In this section, we consider the case when M = { 1 , . . . , N } and propose a solution to the IPF . W e will relax this full measurement condition, i.e., when M ⊂ { 1 , . . . , N } in the next section. For k ∈ { 1 , . . . , K } , the Kirchhof f ’ s laws for a gi ven time index yields I ( k ) = Y V ( k ) . Rewriting this formula in vector form for all time indices yields the following equation for a bus i : I i (1) I i (2) . . . I i ( K ) | {z } I K i = V 1 (1) V 2 (1) . . . V N (1) V 1 (2) V 2 (2) . . . V N (2) . . . . . . . . . . . . V 1 ( K ) V 2 ( K ) . . . V N ( K ) | {z } V K Y i 1 Y i 2 . . . Y iN | {z } Y i . (1) The admittance matrix Y can be obtained from solving the optimization problem below: ˆ Y K,l 2 ≜ arg min Y V K Y − I K F (2) s.t.: Y ∈ S N , Y ii = − X j = i Y ij , ∀ i. in which I K is a K × N matrix, i.e., I K = I K 1 I K 2 . . . I K N . Define vec ( Y ) = Y 11 Y 21 . . . Y N 1 Y 12 Y 22 . . . Y N N T , 3 Fig. 1. An illustrativ e example of a three-node power system with two lines connecting bus 1 to buses 2 and 3 respectively . and apply the vec operator to the objecti ve function, we obtain: min vec ( Y ) ∈ C N 2 × 1 I ⊗ V K vec ( Y ) − vec ( I K ) 2 (3) s.t.: Y ∈ S N , Y ii = − X j = i Y ij , ∀ i, where ⊗ is the Kronecker product. This holds because vec ( AB C ) = ( C T ⊗ A ) v ec ( B ) . Let svec : S N → C ( N 2 − N ) / 2 × 1 be a mapping from a symmetric complex matrix to a complex vector defined as: svec ( Y ) = Y 21 Y 31 . . . Y N 1 Y 32 Y 42 . . . Y N N − 1 T . It can be readily seen that svec is a bijection for any matrix Y that satisfies a) Y ∈ S N and b) 1 T Y = 0 . Based on this definition, we have vec ( Y ) = Γ svec ( Y ) , where Γ ∈ R N 2 × ( N 2 − N ) / 2 maps svec ( Y ) to the v ectorized admittance matrix as illustrated below . Example 1. F or the network depicted in F igur e 1 , the Γ matrix has the following form Y 11 Y 21 Y 31 Y 12 Y 22 Y 32 Y 13 Y 23 Y 33 | {z } vec ( Y ) = − 1 − 1 0 1 0 0 0 1 0 1 0 0 − 1 0 − 1 0 0 1 0 1 0 0 0 1 0 − 1 − 1 | {z } Γ Y 21 Y 31 Y 32 | {z } svec ( Y ) . Based on the definition of Γ in the abo ve equation, the constrained optimization problem can be con verted to an unconstrained one: min svec ( Y ) ∈ C ( N 2 − N ) / 2 × 1 I ⊗ V K Γ | {z } F svec ( Y ) − vec ( I K ) 2 , (4) in which I denotes an identity matrix and F ≜ ( I ⊗ V K )Γ . W e define ˜ F = Re ( F ) − Im ( F ) Im ( F ) Re ( F ) , and ˜ b = Re ( vec ( I K )) Im ( vec ( I K )) . The optimization problem ( 4 ) can be written as an uncon- strained quadratic program in the real domain: min ˜ f ( Y ) ∈ R ( N 2 − N ) × 1 ˜ F ˜ f ( Y ) − ˜ b 2 , (5) in which ˜ f ( Y ) ≜ [ svec ( Y R ) T svec ( Y I ) T ] T , Y R = Re ( Y ) , and Y I = Im ( Y ) . This least square problem yields a solution provided that ˜ M has full column rank: ˜ f ( Y ) = ˜ F T ˜ F − 1 ˜ F T ˜ b. (6) W e compute the solution of the original optimization prob- lem ( 4 ) from the solution of the optimization problem ( 5 ) by taking the in verse map of ˜ f . A sufficient condition to guarantee the exactness of the solution is that V K has full column rank. When V K does not ha ve full column rank, we can characterize the part of the admittance matrix that is identifiable (see [ 3 ]). Proposition 1 (Exactness) . If V K has full column rank, the optimization pr oblem ( 5 ) has a unique solution given by ( 6 ) . Pr oof: Since Γ ∈ R N 2 × ( N 2 − N ) / 2 and has full column rank (this can be checked easily), there exists a matrix Γ † such that Γ † Γ = I . For the Kronecker product I ⊗ V K ∈ C K N × N 2 , I ⊗ V K has full column rank when V K has full column rank; therefore, ˜ F and F ha ve full column rank gi ven the fact that rank ( ˜ F ) = 2 rank ( F ) . Finally , we prove by contradiction that if ˜ F has full column rank, the solution to the optimization problem ( 5 ) is unique. Suppose there exists ˜ f ( Y 1 ) and ˜ f ( Y 2 ) ( ˜ f ( Y 1 ) = ˜ f ( Y 2 ) ) such that ˜ F ˜ f ( Y 1 ) = ˜ b and ˜ F ˜ f ( Y 1 ) = ˜ b , then ˜ F ˜ f ( Y 1 ) − ˜ f ( Y 2 ) = 0 , which contradicts the full column rank assumption. Remark 1. The approac h can be easily extended to the case of nonzer o shunt elements where 1 T Y = 0 by changing the definitions of svec ( Y ) , Γ and ˜ f ( · ) . Specifically if Y includes shunt elements then svec ( Y ) will include diagonal elements ( Y 11 , . . . , Y N N ) . W e can add the element-wise positivity constraint to this problem if the conductance and susceptance of each line are positiv e 1 . min ˜ f ( Y ) ≥ 0 ˜ F ˜ f ( Y ) − ˜ b 2 . (7) The above problem is kno wn as nonne gativ e least squares, which is a con ve x optimization problem and a global mini- mizer can be solved using different methods, such as the acti ve set method [ 20 ]. I V . I P F W I T H H I D D E N N O D E S In the previous section we solve the IPF problem when voltage and current measurements are av ailable at all buses. In this section we consider the case when voltage and current measurements are av ailable only at a subset of the buses. In a distrib ution system, for e xample, measurements are typi- cally av ailable at the substation and customer meters but not throughout the grid. In Section IV -A , we show that, in the presence of hidden nodes, the algorithm presented in Section III can identify a Kron-reduced admittance matrix ¯ Y , defined in ( 11 ) below , for the nodes where measurements are av ailable (Corollary 1 ). In Section IV -B we show that when the network is a tree 1 The conductance of a line is always positive, the susceptance can be negati ve or positive depending on whether the line is inductive or capacitiv e. 4 then it is indeed possible to uniquely identify the original admittance matrix Y from its Kron reduction under reasonable assumptions. A. Kr on-r educed admittance matrix ¯ Y W e call a bus/node a measur ed bus/node if measurements of its voltage and current injection are av ailable for identification. W e call it a hidden bus/node otherwise. Let M and H represent the set of measured and hidden nodes respecti vely . Assumption 1. W e make the following assumptions: 1 .1 The underlying graph G is connected. 1 .2 Hidden nodes have zer o injection I i ( k ) = 0 for all i ∈ H . Note that the second assumption, i.e., the injected current is zero at e very hidden node, is reasonable because if generators or loads are connected to a node, then the current they inject or draw is typically measured. Let H be the number of hidden nodes. W ithout loss of generality we label the buses so that the first N − H b uses are measured and the last H buses are hidden, i.e., M := { 1 , . . . , M ≜ N − H } and H := { N − H + 1 , . . . , N } . W e partition the bus admittance matrix Y into four sub-matrices: Y = Y 11 Y 12 Y 21 Y 22 = G 11 G 12 G 21 G 22 + i B 11 B 12 B 21 B 22 = G + i B . (8) Here Y 11 ∈ S N − H describes the connectivity among the measured nodes, Y 12 = Y T 21 ∈ C ( N − H ) × H the connectivity between the measured and the hidden nodes, and Y 22 ∈ S H the connecti vity among the hidden nodes. For i ∈ M , ( I i ( k ) , V i ( k ) , k = 1 , . . . , K ) denote the current and voltage measurements at bus i at time k . T o simplify notation, we index the entries of Y 22 , not by i, j = 1 , . . . , H , but by i, j = N − H + 1 , . . . , N . W e index the entries of Y 12 by i = 1 , . . . , N − H and j = N − H + 1 , . . . , N and similarly for Y 21 = Y T 12 , as well as submatrices G ij , B ij in ( 8 ). For i ∈ H , I i ( k ) = 0 , ∀ k , and V i ( k ) is the voltage at bus i but is not available for identification. T o simplify notation, W e abuse V 1 ( k ) to denote both the voltage at bus 1 at time k and the vector of all voltages at measured buses at time k , depending on the context; similarly for V 2 ( k ) and I 1 ( k ) . Then I 1 ( k ) 0 = Y 11 Y 12 Y 21 Y 22 V 1 ( k ) V 2 ( k ) , ∀ k (9) If Y 22 is in vertible then eliminating V 2 from ( 9 ) yields a relation I 1 ( k ) = ¯ Y V 1 ( k ) , ∀ k (10) between currents and voltages at measured nodes through the Kr on-reduced admittance matrix ¯ Y ∈ S m defined as: ¯ Y ≜ Y 11 − Y 12 Y − 1 22 Y T 12 (11) for the set of measured nodes. In the rest of this subsection we first justify the in vertibility of Y 22 and hence the definition of ¯ Y . Proposition 1 then implies that ¯ Y can be identified from voltage and current measurements. Moreov er ¯ Y is the best we can identify for general networks because multiple admittance matrices Y may reduce to the same ¯ Y . Assumption 2. The admittance matrix Y ≜ G + i B defined in ( 8 ) satisfies: 2 .1 Series impedances of the lines ar e r esistive and inductive: G [ i, j ] ≤ 0 and B [ i, j ] ≥ 0 for any i = j ; 2 .2 Diagonal dominance: G 22 [ i, i ] ≥ − P j = i G 22 [ i, j ] and − B 22 [ i, i ] ≥ P j = i B 22 [ i, j ] hold for any i . Lemma 1. Under Assumption 1 and 2 , G 22 ⪰ 0 , − B 22 ⪰ 0 and G 22 − B 22 ≻ 0 , when H < N . Pr oof: From the Gershgorin Theorem and Assump- tion 2 .2, all eigen values of the submatrix G 22 lie in the right- half plane including the origin and all eigenv alues of B 22 lie in the left-half plane including the origin. T ogether with the fact that G 22 and B 22 are symmetric, we hav e G 22 ⪰ 0 and − B 22 ⪰ 0 . This implies that G 22 − B 22 ⪰ 0 . W e now show that indeed G 22 − B 22 ≻ 0 . Suppose for the sake of contradiction that there exists a nonzero x ∈ R H such that x T ( G 22 − B 22 ) x = 0 . Denote by A 22 := G 22 − B 22 and A 21 := G 21 − B 21 so that A 22 [ i, i ] = X i,j ∈H : j = i ( − A 22 [ i, j ]) + X i ∈H ,j ∈M ( − A 21 [ i, j ]) Then x T ( G 22 − B 22 ) x = X i,j ∈H : i = j A 22 [ i, j ] x i x j − A 22 [ i, j ] x 2 i + X i ∈H ,j ∈M ( − A 21 [ i, j ]) x 2 i = X i,j ∈H : i 0 and − A 21 [ i, j ] = − G 21 [ i, j ] + B 21 [ i, j ] > 0 for all i ∼ j . Therefore, for x T ( G 22 − B 22 ) x = 0 in ( 12 ), we must ha ve: 1. x i = x j if i ∼ j is a connection between hidden nodes in H ; 2. x i = 0 for any hidden node i ∈ H connected to at least one observed node j ∈ M . Since the network is connected, for e very hidden node i ∈ H , there is a path that connects the hidden node to an observed node k ∈ M . For all nodes j on this path from i to k , the abo ve properties implies that x j = 0 . Since this is true for all hidden nodes, we hav e x = 0 , a contradiction. Hence G 22 − B 22 ≻ 0 . Recall that the network is connected and has N nodes of which H are hidden. Proposition 2. Under Assumptions 1 and 2 , if H < N then Y 22 is in vertible. 5 Pr oof: W e prov e that 0 is not an eigen value of Y 22 . Suppose for the sake of contradiction that there exists a nonzero vector ( v + i w ) such that ( G 22 + i B 22 )( v + i w ) = 0 , i.e., G 22 v − B 22 w = 0 and G 22 w + B 22 v = 0 , This implies ( G 22 + B 22 ) v + ( G 22 − B 22 ) w = 0 , ( G 22 + B 22 ) w − ( G 22 − B 22 ) v = 0 . (13) Since G 22 − B 22 ≻ 0 by Lemma 1 , we can eliminate v from ( 13 ) to get ( G 22 − B 22 ) + ( G 22 + B 22 )( G 22 − B 22 ) − 1 ( G 22 + B 22 ) w = 0 . Multiplying w T on the left we have w T ( G 22 − B 22 ) w + ˆ w T ( G 22 − B 22 ) − 1 ˆ w = 0 , where ˆ w = ( G 22 + B 22 ) w . This contradicts G 22 − B 22 ≻ 0 unless w = 0 . But w = 0 implies that ( G 22 − B 22 ) v = 0 from ( 13 ), meaning v = 0 . This is a contradiction and hence Y 22 is in vertible. Proposition 2 guarantees that Y 22 is in vertible and hence the Kron-reduced admittance ¯ Y in ( 11 ) is well defined under Assumption 2 . Moreov er , because of ( 10 ), the algorithm in Section III can identify the Kron-reduced admittance matrix ¯ Y from voltage and current measurements (Proposition 1 ). Corollary 1. Suppose Assumptions 1 and 2 and the condition in Pr oposition 1 hold. If H < N , then the Kr on-r educed admittance matrix ¯ Y can be identified fr om voltage V K 1 and curr ent I K 1 measur ements at the measur ed nodes. An admittance matrix Y ∈ S N specifies a unique weighted undirected graph G ( Y ) = ( N ( Y ) , E ( Y )) with N := { 1 , . . . , N } and E ⊆ N × N such that there is an edge ( i, j ) if and only if Y [ i, j ] = 0 . Its Kron reduction ¯ Y specifies a unique weighted graph ¯ G := G ( ¯ Y ) = ( M , ¯ E ) that can be obtained from G through Algorithm 1 . Algorithm 1 Graph Condensation Algorithm 1: Input: a graph G = ( N , E ) with N nodes and a set of measured nodes M = { 1 , 2 , . . . , M } and admittance matrix Y 2: f or v = M + 1 : N do 3: Remov e hidden node v from N = N − { v } and all edges from E that are incident on v ; 4: For all node pairs w and l that are neighbors of v , add an edge between w and l to E ; 5: Update the admittance matrix Y = Y / Y [ i, i ] using eq. ( 15 ). 6: end for 7: r eturn ¯ G = G and ¯ Y = Y . Each iterative step in the algorithm remov es a hidden node i ∈ H and connects all its neighbors to each other . This step can be represented algebraically as an update on the admittance matrix to compute its Schur complement of Y [ i, i ] . Specifically we can partition possibly after permutation an admittance matrix Y of an appropriate dimension into the form: Y = Y ( i, i ) Y ( i, i ] Y ( i, i ] T Y [ i, i ] , where Y [ i, i ] ∈ C is the i th diagonal elements of Y , and Y ( i, i ) , Y ( i, i ] are shown in Eq. ( 14 ). Then each iterativ e step of Algorithm 1 updates the admittance matrix by ( Y [ i, i ] is always inv ertible due to Proposition 2 ): Y / Y [ i, i ] = Y ( i, i ) − Y ( i, i ] Y − 1 [ i, i ] Y ( i, i ] T . (15) This step is repeated until all the hidden nodes are removed from the original graph, producing the Kron-reduced graph ¯ G = {M , ¯ E } [ 5 ] and its admittance matrix ¯ Y (sho wn in Fig. 2 ). Giv en an admittance matrix Y , each partition ( Y 11 , Y 12 , Y 22 ) in ( 8 ) defines uniquely a Kron-reduced matrix ¯ Y := ¯ Y ( Y 11 , Y 12 , Y 22 ) gi ven by ( 11 ). This mapping is clearly not injectiv e in general, i.e., giv en an M × M symmetric matrix ¯ Y ∈ S M (possibly with ¯ Y 1 = 0 ) there are generally multiple N × N symmetric matrices Y that can be partitioned into (non- unique) ( Y 11 , Y 12 , Y 22 ) whose Kron reductions are the giv en ¯ Y , as long as N ≥ M . Example 2. Consider a (Kr on-reduced) admittance matrix for a two-node network ( θ = 0 ): ¯ Y = θ − θ − θ θ . The following 3 × 3 admittance matrice Y with the given partition has ¯ Y as its Kron reduction: Y := Y 11 Y 12 Y T 12 Y 22 = θ + θ ′ − θ − θ ′ − θ θ 0 − θ ′ 0 θ ′ , for arbitrary nonzer o θ ′ . The underlying network is shown in F ig. 1 with Node 3 as the hidden node, so that ¯ Y corresponds to the Kron-r educed admittance matrix for Nodes 1 and 2 in F ig. 1 . In this case the hidden node has degr ee 1. Another 3 × 3 admittance matrices ˜ Y that also has ¯ Y as its Kr on r eduction is: ˜ Y := ˜ Y 11 ˜ Y 12 ˜ Y T 12 ˜ Y 22 = θ 1 0 − θ 1 0 θ 2 − θ 2 − θ 1 − θ 2 θ 1 + θ 2 , as long as ( θ 1 , θ 2 ) satisfies θ 1 θ 2 θ 1 + θ 2 = θ . The network underlying ˜ Y is isomorphic to that in F ig. 1 with Node 1 being the hidden node, so that ¯ Y is the Kr on-reduced admittance matrix for Nodes 2 and 3 in F ig. 1 . In this case the hidden node has de gr ee 2. Example 2 sho ws that in general only the Kron-reduced admittance matrix ¯ Y is identifiable from measurements at the measured nodes. For arbitrary networks it is impossible to identify the original admittance matrix Y whose Kron reduction yields ¯ Y . W e next show the surprising result that, when the underlying network is a tree and e very hidden nodes 6 Y ( i, i ) = Y [1 , 1] . . . Y [1 , i − 1] Y [1 , i + 1] . . . Y [1 , N ] . . . . . . . . . . . . . . . . . . Y [ i − 1 , 1] . . . Y [ i − 1 , i − 1] Y [ i − 1 , i + 1] . . . Y [ i − 1 , N ] Y [ i + 1 , 1] . . . Y [ i + 1 , i − 1] Y [ i + 1 , i + 1] . . . Y [ i + 1 , N ] . . . . . . . . . . . . . . . . . . Y [ N , 1] . . . Y [ N , i − 1] Y [ N , i + 1] . . . Y [ N , N ] , Y ( i, i ] = Y [1 , i ] . . . Y [ i − 1 , i ] Y [ i + 1 , i ] · · · Y [ N , i ] . (14) Y ¯ Y Kron Reducti on with resp ect to M AAACFXicbVDLSgMxFM34rPU16tJNsBVcSJmpgi4LbgQRqtgHtEPJpLdtaCYzJBmlDP0JN/6KGxeKuBXc+Tdm2llo64HA4Zx7yT3HjzhT2nG+rYXFpeWV1dxafn1jc2vb3tmtqzCWFGo05KFs+kQBZwJqmmkOzUgCCXwODX94kfqNe5CKheJOjyLwAtIXrMco0Ubq2MdXMhT4FroxTQX8wPQAS1ARUI11iIvtgOgBJTy5Hhc7dsEpORPgeeJmpIAyVDv2V7sb0jgAoSknSrVcJ9JeQqRmlMM4344VRIQOSR9ahgoSgPKSSaoxPjRKF/dCaZ7QeKL+3khIoNQo8M1keqOa9VLxP68V6965lzARxRoEnX7Ui3kaN60Id5k06fnIEEIlM7diOiCSUG2KzJsS3NnI86ReLrknpfJNuVA5zerIoX10gI6Qi85QBV2iKqohih7RM3pFb9aT9WK9Wx/T0QUr29lDf2B9/gDi5Z6W Graph Condensation Algorithm 1 AAACB3icbVDLSsNAFJ3UV62vqEtBBovgqiR1YZeVLnRZwT6gDWUymTRDJzNhZiKU0J0bf8WNC0Xc+gvu/BunaRbaeuDC4Zx7uZzjJ4wq7TjfVmltfWNzq7xd2dnd2z+wD4+6SqQSkw4WTMi+jxRhlJOOppqRfiIJin1Gev6kNfd7D0QqKvi9nibEi9GY05BipI00sk9vJEoi2BI8IFzlIrxmYyGpjmLojuyqU3NywFXiFqQKCrRH9tcwEDiNCdeYIaUGrpNoL0NSU8zIrDJMFUkQnqAxGRjKUUyUl+U5ZvDcKAEMhTTDNczV3xcZipWaxr7ZjJGO1LI3F//zBqkOG15GeZJqwvHiUZgyqAWclwIDKgnWbGoIwiY6xRBHSCKsTXUVU4K7HHmVdOs197JWv6tXm42ijjI4AWfgArjgCjTBLWiDDsDgETyDV/BmPVkv1rv1sVgtWcXNMfgD6/MHPzGY3g== Y/ Y [1 , 1] ... Fig. 2. T wo equiv alent schemes to compute the Kron reduced ¯ Y . has a degree ≥ 3 , then the original admittance matrix Y can indeed be discovered ev en in the presence of hidden nodes. B. Radial networks: exact identification Consider a radial network and suppose we hav e identified a Kron-reduced admittance matrix ¯ Y from partial voltage and current measurements. In this section we dev elop a novel algorithm to compute the original admittance matrix Y from ¯ Y under the following additional assumptions. Assumption 3. The admittance matrix Y satisfies: 3 .1 The underlying graph G ( Y ) is a tr ee. 3 .2 Every hidden node has a degr ee ≥ 3 . 3 .3 Ther e is no shunt element in Y , i.e., Y 1 = 0 . Remark 2. Assumption 3 .2 is a necessary condition for identification as shown in Example 2 where the hidden node has a degr ee 1 or 2. Remark 3. Assumption 3 .3 can be further r elaxed as demon- strated in the full version [ 21 ]. W e start with some definitions. Consider an undirected graph G = ( N , E ) where N := { 1 , . . . , N } is the set of nodes and E ⊆ N × N is the set of edges. A complete graph is one in which all nodes are adjacent. A subgraph of G is a graph G ′ = ( N ′ , E ′ ) with N ′ ⊆ N and E ′ ⊆ E . A clique of G is a complete subgraph of G . A maximal clique of G is a clique that is not a subgraph of another clique of G . W e say G is a tr ee if there is exactly one path between e very two nodes. A for est is a disjoint union of trees. For our purposes, an admittance matrix Y is a complex symmetric matrix (we usually assume Y satisfies Assumption 3 .3). W e sometimes refer to G ( Y ) as the underlying graph of Y and write G := ( N , E ) when Y is clear from the context. Consider two N × N admittance matrices Y 1 and Y 2 . W e define two functions of ( Y 1 , Y 2 ) and their underlying graphs. First Y 3 := Y 1 + Y 2 is also an admittance matrix and its underlying graph G ( Y 3 ) = ( N ( Y 3 ) , E ( Y 3 )) is the graph with the same set of nodes and edges in both graphs, i.e., E ( Y 3 ) := E ( Y 1 ) ∪ E ( Y 2 ) . When the matrices are clear from the context, we also denote the function Y 3 = Y 1 + Y 2 by G 3 = G 1 ⊕ G 2 . Note that if Y 1 and Y 2 satisfy Assumption 3 .3, so does Y 3 . Second define the N × N matrix Y 4 := Y 1 \ Y 2 as a function of ( Y 1 , Y 2 ) by: Y 4 [ i, j ] := Y 1 [ i, j ] if i ∼ j and ( i, j ) ∈ E 2 − P j Y 4 [ i, j ] if i = j 0 otherwise The underlying graph G ( Y 4 ) is a subgraph of G ( Y 1 ) where edges in G ( Y 2 ) have been removed. When the matrices are clear from the conte xt, we also denote the function Y 4 = Y 1 \ Y 2 by G 4 = G 1 / G 2 . Note that Y 4 satisfies Assumption 3 .3 by definition. Fix an (unknown) admittance matrix Y and assume its un- derlying graph G := G ( Y ) is a tree. Suppose its Kron-reduced admittance matrix ¯ Y and its underlying graph ¯ G := G ( ¯ Y ) are giv en. For example ¯ Y is obtained according to Corollary 1 from partial v oltage and current measurements at measured nodes in M . Next, we will propose a recursiv e algorithm to recover Y from ¯ Y . Specifically , W e can decompose ¯ G to two graphs G 1 and G 2 ( ¯ Y 1 and ¯ Y 2 correspondingly) with distinct properties in Section IV -B1 . Secondly , we further introduce a partition of Y in Section IV -B2 and a corresponding parameterization of Y in Section IV -B3 . Thirdly , we can compute these parameters from known quantity in ¯ Y in Section IV -B4 . Finally , the overall recursiv e algorithm to recov er Y is proposed in IV -B5 . 1) Decomposition of ¯ G : Let ¯ E 1 denote the subset of all edges of ¯ G that are between measured nodes in the original graph G , and ¯ E 2 denote the subset of all edges of ¯ G that hav e been added by Step 4 of the graph condensation Algorithm 1 when hidden nodes are remov ed from G . By definition of ¯ E 1 , ¯ E 2 , we have ¯ G = ( M , ¯ E 1 ∪ ¯ E 2 ) . Lemma 2. Under Assumption 1 and 3 , ¯ E 1 ∩ ¯ E 2 = ∅ . Pr oof: If there exists an edge ( i, j ) ∈ ¯ E 1 ∩ ¯ E 2 , then ( i, j ) must be an edge in the original graph G and nodes i and j must also be connected through a path consisting of only hidden nodes. This creates a loop and contradicts that G is a tree. Hence ¯ E 1 ∩ ¯ E 2 = ∅ . Since ¯ G = ( M , ¯ E 1 ∪ ¯ E 2 ) , Lemma 2 motiv ates decomposing ¯ G into two subgraphs, G 1 := ( M , ¯ E 1 ) and G 2 := ( M , ¯ E 2 ) , both defined on M of measured nodes but with disjoint edge sets. While the graph ¯ G := G ( ¯ Y ) is defined by the Kron- reduced admittance matrix ¯ Y , at this point the graphs G 1 , G 2 are only defined in terms of the graph ¯ G (in fact in terms of G ) 7 and are not associated with any admittance matrices. Define the matrices: ¯ Y 1 := Y 11 − diag { 1 T Y 11 } (16a) ¯ Y 2 := diag { 1 T Y 11 } − Y 12 Y − 1 22 Y T 12 . (16b) The key observation, stated in the next result, is that G 1 , G 2 hav e simple structures, that the matrices defined in ( 16 ) are indeed admittance matrices, and that G 1 , G 2 are the underlying graphs of these admittance matrices. Even though we do not know the submatrices Y 11 , Y 12 , Y 22 of Y , the simple structures of G 1 , G 2 allow us to compute ¯ Y 1 , ¯ Y 2 as we will see. Theorem 1 (Separability) . Suppose the admittance matrix Y satisfies Assumptions 1 , 2 and 3 . Then 1. G 1 is a forest. 2. G 2 = ⊕ i C i for some C i ar e edge-disjoint maximal cliques each with mor e than 2 nodes. 2 3. G 1 = G ( ¯ Y 1 ) and G 2 = G ( ¯ Y 2 ) so that ¯ G = G 1 ⊕ G 2 . Pr oof: For the first assertion, G 1 is a forest since it is a subgraph of the tree G . For the second assertion G 2 is a collection of maximal cliques C i due to Step 4 of the graph condensation Algorithm 1 . T o show that the maximal clique (in each) C i is of size at least 3, suppose C i consists of m i (measured) nodes and, in the original graph G , these m i mea- sured nodes “surround” h i hidden nodes, i.e., the neighbors of each of these hidden nodes are either hidden nodes or nodes in C i . Let d j denote the degrees of hidden nodes j = 1 , . . . , h i . These m i + h i nodes form a (connected) subtree of G with exactly m i + h i − 1 edges. Since m i of these edges are between measured and hidden nodes and h i − 1 edges are between hidden nodes, we must have P h i j =1 d j = m i + 2( h i − 1) and hence m i = 2 + P h i i =1 ( d i − 2) . Since h i ≥ 1 and d i ≥ 3 (Assumption 3 .2), we have m i ≥ 3 . T o show that C i and C j are edge-disjoint, suppose for the sake of contradiction that there is an edge ( k , l ) in both C i and C j . By the definition of G 2 , ( k , l ) is not an edge in the original graph G . Since nodes k , l are both in C i , there is a path from k to l in G that consists of only hidden nodes connected to measured nodes in the maximal clique C i . Since nodes k , l are both in C j , there is disjoint path from k to l in G that consists of a set of hidden nodes connected to nodes in C j . These two paths form a loop in G , a contradiction. Hence C i and C j do not share any edge in G 2 . For the third assertion, note that the matrix ¯ Y 1 defined in ( 16a ) is symmetric and hence an admittance matrix. The diagonal entry Y 11 [ i, i ] of Y 11 is the negativ e sum of the off-diagonal entries of the i th row/column of the original admittance matrix Y (plus any shunt element at b us i ), so that the i th entry of 1 T Y 11 is equal to the i th row sum of Y 12 (plus any shunt element at bus i ). Hence ¯ Y 1 satisfies Assumption 3 .3 if Y does. Moreover , by the definition of G 1 , the edges in E 1 correspond exactly to the off-diagonal entries of Y 11 that are 2 Strictly speaking, each C i is a subgraph of G 2 with M as its node set. It consists of a single maximal clique and the remaining isolated nodes in M . W e will abuse notation and use C i to both refer to this subgraph of G 2 and to the maximal clique in C i . nonzero. This implies that the graph G ( ¯ Y 1 ) that underlies the admittance matrix in ( 16a ) is indeed G 1 . The matrix ¯ Y 2 defined in ( 16b ) is also symmetric and hence is an admittance matrix. If Y satisfies Assumption 3 .3, then 1 T Y 11 = − 1 T Y T 12 , 1 T Y 12 = − 1 T Y 22 . This implies diag { 1 T Y 11 } = diag { 1 T Y 12 Y − 1 22 Y T 12 } , i.e., ¯ Y 2 defined in ( 16b ) satisfies Assumption 3 .3 when Y does. Next we show that G 2 = G ( ¯ Y 2 ) . From ( 16 ) ¯ Y = Y 11 − Y 12 Y − 1 22 Y T 12 = ¯ Y 1 + ¯ Y 2 , and hence ¯ G := G ( ¯ Y ) = G 1 ⊕ G 2 with G ( ¯ Y 1 ) = G 1 . Therefore we have G ( ¯ Y 2 ) = G 2 . This concludes the proof. Remark 4. F r om the thir d assertion, we have shown that, once G 1 and G 2 ar e obtained fr om G , ¯ Y 1 and ¯ Y 2 defined in ( 16 ) can be obtained. There are many algorithms for solving the clique problem, such as the Bron-K erbosch algorithm, which we adopt in Algorithm 2 . Algorithm 2 Graph Decoupling Algorithm 1: Input: a condensed graph ¯ G 2: Set G ′ = ¯ G , i = 1 . 3: while G ′ has a clique with more than two nodes do 4: Use Bron-Kerbosch Algorithm to find a clique ( ≥ 3 nodes, together with other isolating nodes) C i in G ′ , 5: Let G ′ = G ′ / C i , i = i + 1 , 6: end while 7: r eturn G 2 = ⊕ i C i , G 1 = ¯ G / G 2 and the corresponding ¯ Y 1 and ¯ Y 2 . 2) P artition of Y : Next we propose an algorithm to obtain Y 11 , Y 22 and Y 12 , and therefore the original admittance matrix Y . The decomposition of ¯ G into G 1 and G 2 guaranteed by Theorem 1 allows us to partition the set M into a subset of internal measured nodes that are not connected to any hidden nodes and a disjoint subset of boundary measured nodes that connect to some hidden nodes. W e can similarly partition H into a subset of internal hidden nodes that are not connected to any measured nodes and the disjoint subset of boundary hidden nodes that connect to some measured nodes. The decomposition of ¯ G into G 1 and G 2 identifies only the types of measured nodes, but not those of hidden nodes. W e can hence arrange the original admittance matrix Y into the follo wing structure (only the upper triangular submatrix is shown as Y is symmetric): Y = Y 11 Y 12 Y 22 =: Y 11 , 11 Y 11 , 12 0 0 Y 11 , 22 Y 12 , 21 0 Y 22 , 11 Y 22 , 12 Y 22 , 22 . (17a) Here, for Y 11 , the submatrix Y 11 , 11 corresponds to connecti vity among the internal measured nodes, Y 11 , 22 corresponds to 8 connectivity among the boundary measured nodes, and Y 11 , 12 corresponds to connectivity between the internal and boundary measured nodes. Similarly , for Y 22 , the submatrix Y 22 , 11 corresponds to connectivity among the boundary hidden nodes, Y 22 , 22 to that among the internal hidden nodes, and Y 22 , 12 to that between the internal and boundary hidden nodes. The submatrix Y 12 , 21 corresponds to connectivity between the set of boundary measured nodes and the set of boundary hidden nodes. Denote the in verse Y − 1 22 by: X 22 := Y − 1 22 =: X 22 , 11 X 22 , 12 X T 22 , 12 X 22 , 22 . W e hav e ¯ Y = Y 11 − Y 12 Y − 1 22 Y T 12 = Y 11 , 11 Y 11 , 12 Y T 11 , 12 Y 11 , 22 − 0 0 0 Y 12 , 21 X 22 , 11 Y T 12 , 21 , (17b) where X 22 , 11 = ( Y 22 , 11 − Y 22 , 12 Y − 1 22 , 22 Y T 22 , 12 ) − 1 from ( 17a ) and the W oodbury formula. Specifically , gi ven the definition of Schur complement det( Y 22 , 11 − Y 22 , 12 Y − 1 22 , 22 Y T 22 , 12 ) det Y 22 , 22 = det Y 22 , and from the in vertibility of Y 22 (shown in Proposition 2 ), the right-hand side of the above equation is nonzero. Therefore det( Y 22 , 11 − Y 22 , 12 Y − 1 22 , 22 Y T 22 , 12 ) cannot be zero, and as a result, the in vertibility of X 22 , 11 can be guaranteed. Since we can compute ¯ Y from partial voltage and current measurements, we can identify submatrices Y 11 , 11 and Y 11 , 12 for internal measured nodes from ¯ Y according to ( 17b ). The edges in E 1 correspond to the off-diagonal entries of [ Y 11 , 11 Y 11 , 12 ] as well as Y T 11 , 12 , and they form a forest (Theo- rem 1 ). The edges in E 2 correspond to the off-diagonal entries of Y 11 , 22 − Y 12 , 21 X 22 , 11 Y T 12 , 21 , and they form a collection of cliques. Recall that both G 1 and G 2 hav e M as their node set; see the example in Fig. 3 . In the rest of this subsection we focus on identifying the remaining submatrices Y 11 , 22 , Y 12 , 21 as well as Y 22 (or specifically , Y 22 , 11 , Y 22 , 12 , Y 22 , 22 ) of Y . For this purpose we assume without loss of generality that all measured nodes are boundary measured nodes, i.e., the ro ws and columns corresponding to submatrices Y 11 , 11 and Y 11 , 12 as well as their contributions to the diagonal entries of Y 11 , 22 hav e been remov ed from Y . Then Y = Y 11 Y 12 Y 22 =: Y 11 , 22 Y 12 , 21 0 Y 22 , 11 Y 22 , 12 Y 22 , 22 . (18) Our goal is to identify Y in ( 18 ) giv en its Kron-reduction: ¯ Y = Y 11 , 22 − Y 12 , 21 X 22 , 11 Y T 12 , 21 . Theorem 1 .2 implies that the underlying Kron-reduce graph G ( ¯ Y ) is a disjoint collection of maximal cliques C i among boundary measured nodes. By hidden nodes in a maximal clique C i of the Kron-reduced graph ¯ G , we mean the nonempty set of hidden nodes in the original graph G that are connected either to the measured nodes in C i or other hidden nodes in C i . A measured node can be in multiple cliques C i though C i are edge-disjoint (Theorem 1 .2). Lemma 3. Suppose the admittance matrix Y satisfies Assump- tions 1 , 2 and 3 . A measur ed node can connect to only one hidden node in any cliques C i of which it is a member . Pr oof: If a measured node connects to more than one hidden node in a maximal cliques C i , there exists a loop since there is a path between any two hidden nodes in C i , hence a contradiction. W e further assume, without loss of generality , that G ( ¯ Y ) consists of a single clique; otherwise, we can repeatedly apply Algorithm 3 below to each clique separately to determine the corresponding submatrices and then combine them to obtain Y 22 and Y 12 . The general case has been discussed and elaborated in the Appendix. Remark 5. W ith this further assumption, Lemma 3 guarantees that Y 12 has exactly one nonzer o element in each r ow . 3) P arameterization of Y : Recall that there are M (bound- ary) measured nodes, indexed by 1 , . . . , M , so that Y 11 , 22 in ( 18 ) is M × M . Suppose there are H b boundary hidden nodes, index ed by M + 1 , . . . , M + H b , and H i := H − H b internal hidden nodes, indexed by M + H b + 1 , . . . , M + H . Then Y 22 , 11 in ( 18 ) is H b × H b and Y 22 , 22 is H i × H i . Suppose each measured node i ∈ { 1 , . . . , M } is connected to the hidden node h ( i ) ∈ { M + 1 , . . . , M + H b } by a line with series admittance y ih ( i ) . From Remark 5 we know there is a unique h ( i ) for each i , but voltage and current measurements only identify the identity of each measured node i , but not the hidden node h ( i ) it is connected to (nor the v alues of H , H b , H i ). The next result suggests a method to identify all measured nodes that are connected to the same boundary hidden node. Proposition 3. Suppose the admittance matrix Y satisfies Assumptions 1 , 2 and 3 . T wo measured nodes i and j are connected to the same hidden node if and only if the off- diagonal entries of r ows i and j of ¯ Y 2 ar e pr oportional, i.e., ther e exists γ ( i, j ) = 0 such that ¯ Y 2 [ i, k ] ¯ Y 2 [ j, k ] = γ ( i, j ) , k = i, j, k = 1 , . . . , M . Pr oof: As discussed above, each row of Y 12 , 21 has exactly a single nonzero entry . For each measured node i , let h ( i ) denote the hidden node to which i is adjacent and y ih ( i ) denotes the series admittance of line ( i, h ( i )) . Then Y 12 , 21 =: − y 1 h (1) e T h (1) . . . − y M h ( M ) e T h ( M ) where the unit vector e i ∈ { 0 , 1 } H b is the column vector with a single 1 in the i th entry and 0 elsewhere. Denote the ( i, j ) th entry of X 22 , 11 by β ij . Then the i th ro w of matrix Y 12 , 21 X 22 , 11 Y T 12 , 21 in ( 16b ) is y ih ( i ) β h ( i ) h (1) y 1 h (1) · · · β h ( i ) h ( M ) y M h ( M ) . The necessity of the proposition can be prov en as follows: if i and j are connected to the same hidden node, i.e., h ( j ) = 9 Fig. 3. Illustration of admittance matrices and their underlying graphs: ( Y , G ( Y )) , ( ¯ Y , ¯ G ) , ( ¯ Y 1 , G 1 ) , and ( ¯ Y 2 , G 2 ) . h ( i ) , then b T h ( j ) = b T h ( i ) . Therefore, for k = i, j , we have from ( 16b ) ¯ Y 2 [ i, k ] = γ ( i, j ) ¯ Y 2 [ j, k ] for γ ( i, j ) := y ih ( i ) /y j h ( j ) . W e now present the sufficiency proof. Fix i = j , and suppose there exists γ ′ ( i, j ) = 0 such that ¯ Y 2 [ i, k ] = γ ′ ( i, j ) ¯ Y 2 [ j, k ] for all k = i, j . Then γ β h ( i ) h ( k ) = β h ( j ) h ( k ) , k = i, j (19) where γ := y ih ( i ) / γ ′ ( i, j ) y j h ( j ) . W e hav e to pro ve that h ( i ) = h ( j ) . Suppose h ( i ) = h ( j ) . Assume without loss of generality that h ( i ) = 1 and h ( j ) = 2 (corresponding to nodes M + 1 and M + 2 respectiv ely). By Assumptions 3 .3, Y is symmetric and hence ¯ Y , X 22 := Y − 1 22 , and X 22 , 11 are symmetric. Then ( 19 ) means that the matrix X 22 , 11 is of the form X 22 , 11 = β 11 β 12 β 13 · · · β 1 H b β 12 β 22 γ β 13 · · · γ β 1 H b β 13 γ β 13 β 33 · · · β 3 H b . . . . . . . . . . . . . . . β 1 H b γ β 1 H b β 3 H b · · · β H b H b (20a) =: A 11 A 12 A T 12 A 22 . (20b) In particular rank ( A 12 ) ≤ 1 . Suppose more than one measured nodes are connected to each of the hidden nodes h ( i ) and h ( j ) . Specifically let h ( i ) = h ( i ′ ) and h ( j ) = h ( j ′ ) for distinct i ′ , j ′ . Then letting k = i ′ and then k = j ′ in ( 19 ) implies that β h ( i ) h ( i ) = β h ( i ) h ( i ′ ) = 1 γ β h ( j ) h ( i ′ ) = 1 γ β h ( j ) h ( i ) β h ( j ) h ( j ) = β h ( j ) h ( j ′ ) = γ β h ( i ) h ( j ′ ) = γ β h ( i ) h ( j ) . This implies that β 21 = β 12 = γ β 11 and β 22 = γ β 12 in ( 20a ), i.e., the first two rows of X 22 , 11 are linearly dependent. This contradicts the in vertibility of X 22 , 11 and hence h ( i ) = h ( j ) if more than one measured nodes are connected to each of h ( i ) and h ( j ) . Suppose then at least one of h ( i ) , h ( j ) is connected to a single measured node (i.e., i or j ) so that γ β h ( i ) h ( k ) = β h ( j ) h ( k ) in ( 19 ) may not hold for k = i or k = j or both. Let X − 1 22 , 11 =: B 11 B 12 B T 12 B 22 (21) where B 11 ∈ C 2 × 2 and B 22 ∈ C ( H b − 2) × ( H b − 2) . W e no w prov e rank ( B 12 ) ≤ 1 through direct calculation. W e then deriv e a contradiction first for the case when rank ( B 12 ) = 0 and then when rank ( B 12 ) = 1 . This implies that h ( i ) = h ( j ) , i.e., i and j are connected to the same (boundary) hidden node. Proof of rank ( B 12 ) ≤ 1 . From ( 20 ) and ( 21 ), the entries of B 12 are giv en by: for k ≥ 3 , [ B 12 ] 1 k = ( − 1) k +1 det ( X 22 , 11 ) β 12 b T β 22 γ b T γ c − k B − k (22a) [ B 12 ] 2 k = ( − 1) k det ( X 22 , 11 ) β 11 b T β 12 γ b T c − k B − k (22b) where b T is the first row of A 12 , c − k is the H b − 3 dimensional (column) vector which is the first column of A T 12 with its k th entry removed, and B − k is the ( H b − 3) × ( H b − 2) matrix which is A 22 with its k th row remov ed. For ease of reference we state the following calculation as a lemma without proof. Lemma 4. Given any scalars α, β , γ ∈ C , any column vectors b, c , and any matrix B of matching sizes, we have the following 10 determinant α b T β γ b T c B = ( γ α − β ) b T B . Applying Lemma 4 to ( 22 ) we hav e [ B 12 ] 1 k = ( − 1) k +1 a 1 d k , k ≥ 3 [ B 12 ] 2 k = ( − 1) k a 2 d k , k ≥ 3 where a 1 := γ β 12 − β 22 det ( X 22 , 11 ) , a 2 := γ β 11 − β 12 det ( X 22 , 11 ) , d k := b T B − k i.e., B 12 is of the form B 12 = a 1 d 3 − a 1 d 4 · · · ( − 1) H b − 1 a 1 d H b − 2 − a 2 d 3 a 2 d 4 · · · ( − 1) H b a 2 d H b − 2 (23) This shows that rank ( B 12 ) ≤ 1 . Case 1: rank ( B 12 ) = 0 . The ke y observation is that, X − 1 22 , 11 = Y 22 / Y 22 , 22 is the admittance matrix of the Kron- reduced network consisting of only boundary hidden nodes where two boundary hidden nodes h ( k ) and h ( k ′ ) are adjacent if and only if, in the original netw ork, h ( k ) and h ( k ′ ) are either adjacent or are connected by a path with only internal hidden nodes (ev en though Y 22 may not hav e zero row sums because each diagonal entry of Y 22 may include line admittances that can be interpreted as a shunt admittance in the Kron-reduced network). Note that h ( k ) and h ( k ′ ) cannot both be adjacent and connected by a path with only internal hidden nodes because, otherwise, there is a loop in the original network. Then B 12 describes the connectivity of nodes h ( i ) , h ( j ) with other boundary hidden nodes in the Kron-reduced network. When rank ( B 12 ) = 0 (either a 1 = a 2 = 0 or d k = 0 for all k ≥ 3 in ( 23 )), then h ( i ) and h ( j ) are connected at most to each other in the Kron-reduced network represented by Y 22 / Y 22 , 22 . If either h ( i ) or h ( j ) is connected to a single measured node, this violates Assumptions 3 .2 that every hidden node has degree at least 3 in the original (not Kron- reduced) network. Case 2: rank ( B 12 ) = 1 . Then either a 1 > or a 2 > 0 , or both. Suppose one of h ( i ) = 1 and h ( j ) = 2 , say , h ( i ) is connected to more than one measured node (and h ( j ) is connected to a single measured node). The argument above shows that γ β h ( i ) h ( i ) = β h ( j ) h ( i ) = β h ( i ) h ( j ) and hence a 2 = 0 and [ B 12 ] 2 k = 0 for all k ≥ 3 . This means that h ( j ) = 2 is connected only to a single measured node and at most the boundary hidden node h ( i ) in the Kron-reduced network represented by Y 22 / Y 22 , 22 , violating Assumptions 3 .2 that e very hidden node has degree at least 3 in the original network. Suppose then both h ( i ) and h ( j ) are connected to a single measured node each. Then a 1 > 0 and a 2 > 0 because otherwise, either [ B 12 ] 1 k = 0 or [ B 12 ] 2 k = 0 for all k ≥ 3 , violating Assumptions 3 .2. But ( 23 ) implies that [ B 12 ] 1 k = 0 if and only if [ B 12 ] 2 k = 0 , i.e., h ( i ) and h ( j ) are connected to exactly the same set of boundary hidden nodes in the Kron- reduced network represented by Y 22 / Y 22 , 22 . If h ( i ) and h ( j ) are connected to each other and no common boundary hidden node or if they are connected to only 1 common hidden node (and not to each other), then Assumptions 3 .2 is violated. If they are connected to 2 or more common boundary hidden nodes, then there must be a loop in the original (not Kron- reduced) network, violating Assumptions 3 .1. This completes the proof of suf ficiency . Note that if there are only M = 3 measured nodes then Assumptions 3 .1 and 3 .2 imply that all of them must be connected to the same boundary hidden node. Giv en the Kron-reduced admittance matrix ¯ Y 2 , Proposition 3 allows us to group together the (boundary) measured nodes that are connected to the same (boundary) hidden node. This also identifies the number of boundary hidden nodes, e ven though we do not know (yet) the number or identity of internal hidden nodes nor the connectivity among the nodes. W e can re-arrange the submatrix matrix Y 12 , 21 into a form easier for identification. Specifically let measured nodes 1 , . . . , k 1 be connected to hidden node M + 1 , measured nodes k 1 + 1 , . . . , k 2 to hidden node M + 2 , . . . , measured nodes k H b − 1 + 1 , . . . , k H b := M to hidden node M + H b . Note that Proposition 3 yields the values for H b and ( k 1 , k 2 , . . . , k H b = M ) even though it provides no information about the value of H , the total number of hidden nodes. T o simplify notation, denote the series admittance y ih ( i ) of line ( i, h ( i )) by y i . Then Y 12 = Y 12 , 21 0 where Y 12 , 21 is M × H b and can be arranged into the following simple form: Y 12 , 21 = − y 1 0 · · · 0 . . . . . . . . . . . . − y k 1 0 · · · 0 0 − y k 1 +1 · · · 0 . . . . . . . . . . . . 0 − y k 2 · · · 0 . . . . . . . . . . . . 0 0 · · · − y k H b =: − ˆ y 1 0 · · · 0 0 − ˆ y 2 · · · 0 . . . . . . . . . . . . 0 0 · · · − ˆ y H b , where for i = 1 , . . . , H b , ˆ y i is a ( k i − k i − 1 ) -dimensional column vector corresponding to k i − k i − 1 measured nodes that are connected to the hidden node M + i . Since Y has zero row sum by Assumption 3 .3, the diagonal matrix diag { 1 T Y 11 } = diag { 1 T Y 11 , 22 } = diag ( y i = y ih ( i ) , i = 1 , . . . , M ) . W e ha ve Y 12 Y − 1 22 Y T 12 = diag ( ˆ y j ) X 22 , 11 diag ( ˆ y T j ) = β 11 ˆ y 1 ˆ y T 1 β 12 ˆ y 1 ˆ y T 2 · · · β 1 H b ˆ y 1 ˆ y T H b β 21 ˆ y 2 ˆ y T 1 β 22 ˆ y 2 ˆ y T 2 · · · β 2 H b ˆ y 2 ˆ y T H b . . . . . . . . . . . . β H b 1 ˆ y H b ˆ y T 1 β H b 2 ˆ y H b ˆ y T 2 · · · β H b H b ˆ y H b ˆ y T H b . Then the admittance matrix corresponding to the graph G 2 in 11 Theorem 1 is: ¯ Y 2 = diag ( ˆ y 1 ) 0 · · · 0 diag ( ˆ y 2 ) · · · 0 . . . . . . diag ( ˆ y M ) − β 11 ˆ y 1 ˆ y T 1 β 12 ˆ y 1 ˆ y T 2 · · · β 1 H b ˆ y 1 ˆ y T H b β 22 ˆ y 2 ˆ y T 2 · · · β 2 H b ˆ y 2 ˆ y T H b . . . . . . β H b H b ˆ y H b ˆ y T H b . (24) Recall that we hav e already identified the Kron-reduced admit- tance matrix ¯ Y 2 , i.e., we know every entry of ¯ Y 2 on the left- hand side of ( 24 ). W e now explain how to identify ( β ij , i, j = 1 , . . . , H b ) and ( y i = y ih ( i ) , i = 1 , . . . , M ) on the right-hand side of ( 24 ). In particular , ( y i = y ih ( i ) , i = 1 , . . . , M ) yields Y 12 of the original admittance matrix Y . 4) Computation of parameters in Y : Let ¯ Y 2 ,k 1 be the diagonal submatrix consisting of the first k 1 rows and columns of ¯ Y 2 corresponding to the first k 1 measured nodes connected to the first hidden node M + 1 : ¯ Y 2 ,k 1 := diag ( ˆ y 1 ) − β 11 ˆ y 1 ˆ y T 1 = y 1 · · · 0 . . . . . . y k 1 − β 11 y 1 . . . y k 1 y 1 · · · y k 1 . (25) W e first explain how to identify ( β 11 , ˆ y 1 ) on the right-hand side of ( 25 ) from the kno wledge of ¯ Y 2 ,k 1 on the left-hand side of ( 25 ). The identification of other β ii , ˆ y i corresponding to k i − k i − 1 measured nodes connected to the hidden node M + i from the diagonal blocks ¯ Y 2 ,k i := diag ( y k i − 1 − 1 , . . . , y k i ) − β ii ˆ y i ˆ y T i can be done similarly . Case 1: k 1 ≥ 2 . In this case, hidden node M +1 is connected to two or more measured nodes indexed by i = 1 , . . . , k 1 . Consider the first two measured nodes and the corresponding 2 × 2 principal submatrix of Y 2 ,k 1 : for i, j = 1 , 2 ¯ Y 2 ,k 1 [ i, j ] = y i − β 11 y 2 i if i = j − β 11 y i y j if i = j (26) This leads to the following equations in ( β 11 , y 1 , y 2 ) : y 1 − β 11 y 2 1 = ¯ Y 2 ,k 1 [1 , 1] =: a 1 − β 11 y 1 y 2 = ¯ Y 2 ,k 1 [1 , 2] =: a 2 y 2 − β 11 y 2 2 = ¯ Y 2 ,k 1 [2 , 2] =: a 3 yielding: y 1 = a 1 a 3 − a 2 2 a 2 + a 3 , y 2 = a 1 a 3 − a 2 2 a 1 + a 2 , β 11 = − a 2 ( a 1 + a 2 )( a 2 + a 3 ) ( a 1 a 3 − a 2 2 ) 2 . (27) T o identify other ( y j , j > 2) , note that − β 11 y 1 y j = Y 2 ,k 1 [1 , j ] , j = 3 , . . . , k 1 yielding y j = − Y 2 ,k 1 [1 , j ] β 11 y 1 , where β 11 and y 1 are given by ( 27 ). Once ˆ y 1 , . . . , ˆ y k j are found, we can calculate from off-diagonal entries of ¯ Y 2 all β ij from ( 24 ). Case 2: Once we have recov ered the coefficients for hidden boundary nodes with at least two connections to measured nodes in Case 1, next, we can treat these recovered hidden nodes as measured nodes and repeat the above procedure until no hidden node is left. A key step is to construct a new Kron reduced matrix once parts of the admittance matrix have been found. Let the original Y ha ve the following partition as in ( 18 ): Y = Y 11 , 22 Y 12 , 21 0 Y 22 , 11 Y 22 , 12 Y 22 , 22 . The Kron reduced admittance matrix can be decomposed to ¯ Y 1 and ¯ Y 2 . Specifically , ¯ Y 2 has the following form: ¯ Y 2 = diag { 1 T Y 11 , 22 } − Y 12 , 21 X 22 , 11 Y T 12 , 21 = diag { 1 T Y 11 , 22 } − Y 12 , 21 0 Y 22 , 11 Y 22 , 12 Y T 22 , 12 Y 22 , 22 − 1 Y T 12 , 21 0 . Based on the results in Case 1, one can recov er diag { 1 T Y 11 , 22 } , Y 12 , 21 and X 22 , 11 . Since ¯ Y 1 is known from Algorithm 2 , diag { 1 T Y 11 , 22 } allows us to compute Y 11 from the equality ( 16a ) and the partition in ( 18 ): Y 11 = ¯ Y 1 + diag { 1 T Y 11 , 22 } . (28) Hence the entire rows and columns of Y corresponding to the boundary measured nodes are known after ( 28 ). W e can then focus on the submatrices Y 22 , 11 , Y 22 , 12 , Y 22 , 22 corresponding to only the boundary and internal hidden nodes, i.e., we can reduce the unkno wn admittance matrix Y to the new smaller admittance matrix below , which amounts to restricting attention to the subgraph without the boundary measured nodes. Y = Y 22 , 11 Y 22 , 12 Y 22 , 22 . The Kron reduced admittance matrix of this new (unknown) admittance matrix Y can then be obtained from the knowledge of X 22 , 11 : ¯ Y := Y 22 , 11 − Y 22 , 12 Y − 1 22 , 22 Y T 22 , 12 = X − 1 22 , 11 . Moreov er , we have identified the set of boundary hidden nodes. Applying Theorem 1 , Algorithm 2 and Proposition 3 to this new ¯ Y allo ws us to identify a set of internal hidden nodes to which this set of boundary hidden nodes are connected. Moreov er , we can treat the set of boundary hidden nodes as boundary measured nodes and the newly identified internal hidden nodes as boundary hidden nodes. Therefore, ev en though we do not kno w the number or the identity of the r emaining internal hidden nodes, we can partition the new (unknown) admittance matrix Y into the form at the beginning of Case 2 and therefore repeat the computation on this new (smaller) admittance matrix recursively , strictly reducing the 12 number of internal hidden nodes in each iteration until the set of internal hidden nodes becomes null. Case 3: For any hidden node that connects to one or zero measured node, these hidden nodes will ev entually hav e more than one connection to measured nodes once the other hidden nodes hav e been recov ered and therefore can be recovered. It is easy to show that there will nev er exist a scenario that all the hidden nodes have at most 1 connection to measured nodes for a tree graph. T o see this, note that for any clique, H ≥ M as ev ery hidden node connects to a different measured node. On one hand, the sum of all hidden nodes’ degrees is greater than 3 H under Assumption 3 . On the other hand, it is at most 2( H − 1) + M , which is twice the sum of all edges between hidden nodes and the number of connections between hidden nodes and measured nodes. Ho wev er , 2( H − 1) + M < 3 H , a contradiction. Case 4: If all hidden nodes are hidden boundary nodes, i.e., Y 22 = Y 22 , 11 , then Y 22 = X − 1 22 , 11 and hence the entire admittance matrix Y can be identified. If there are hidden nodes that are not hidden boundary nodes, we can treat hidden boundary nodes as measured nodes now and repeat the above procedure based on Case 2. 5) Over all r ecursive algorithm: The overall identification procedure is summarized in Algorithm 3 3 . Algorithm 3 Recover Y from ¯ Y 1: Input: ¯ Y 1 and ¯ Y 2 2: f or each pair of nodes ( j, k ) do 3: Compute γ [ j, k ] from ¯ Y 2 . 4: end for 5: Solv e for diag { 1 T Y 11 , 22 } , Y 12 , 21 and X 22 , 11 from ( 25 ), set ˆ Y = ¯ Y 1 + diag { 1 T Y 11 , 22 } Y 12 , 21 Y T 12 , 21 X − 1 22 , 11 := ˆ Y 11 ˆ Y 12 ˆ Y T 12 ˆ Y 22 and set ¯ Y 2 = X − 1 22 , 11 6: if the graph corresponding to ¯ Y 2 , i.e., G ( ¯ Y 2 ) is not radial then 7: for each pair of nodes ( j, k ) do 8: Compute γ [ j, k ] from ¯ Y 2 . 9: end for 10: Solve for diag { 1 T Y 11 , 22 } , Y 12 , 21 and X 22 , 11 from ( 25 ) and set ˆ Y = ˆ Y 11 ˆ Y 12 0 ˆ Y T 12 diag { 1 T Y 11 , 22 } Y 12 , 21 0 Y T 12 , 21 X − 1 22 , 11 . 11: Set ˆ Y 11 = ˆ Y 11 ˆ Y 12 ˆ Y T 12 diag { 1 T Y 11 , 22 } , ˆ Y 12 = 0 Y 12 , 21 . 12: Set ¯ Y = X − 1 22 , 11 and apply Algorithm 2 to obtain ¯ Y 1 and ¯ Y 2 . 13: end if 14: r eturn Y = ˆ Y 3 For notational simplicity , we assume without loss of generality that all measured nodes are boundary measured nodes. Y et, this assumption can easily be relaxed. 1 2 3 4 5 6 7 8 9 10 11 12 13 14 bus number 1 2 3 4 5 6 7 8 9 10 11 12 13 14 bus number 0 1 2 3 4 5 6 7 8 # 10 -5 Fig. 4. The identification error when there is no hidden state. The color of a cell located at row i and column j represents the value of | Y [ i, j ] − ˆ Y [ i, j ] | . It can be seen that the identification error using measurements of 15 time slots ( K = 15 ) is quite small compared to the absolute value of the elements of Y . V . S I M U L AT I O N S A. Simulated data In this section we implement the proposed algorithms in MA TLAB and ev aluate their identification accuracy by per- forming simulations in MA TPO WER [ 22 ]. The optimization problems are solved using the CVX toolbox [ 23 ]. W e run power flow analysis on the IEEE 14 -bus test system, repre- senting a portion of a po wer system in the Midwestern U.S., which has 14 buses, 11 aggregated loads, and 5 generators, 3 of which are synchronous compensators used for reactive power support [ 24 ]. T o validate our identification algorithms in three-phase distribution systems, we have also performed simulations on IEEE test distribution networks, which are not presented here (see our previous work [ 25 ] for the identifica- tion results in distribution systems). W e assume that PMUs are installed at selected buses and that they can precisely measure the voltage and current magnitudes and phase angles that we obtain from power flow calculations, unless stated otherwise. For each scenario, we run 100 steady state simulations, each pertaining to a time slot, to determine the voltage and current magnitude and phase angle of e very bus, while varying the real and reacti ve po wer demand of the loads across the time slots. Specifically , for a giv en time slot, the real and reactiv e power consumption of a constant PQ load are computed by multiplying a scaling factor drawn from a uniform distrib ution ov er the interval [0 . 8 , 1 . 2] by the real and reactiv e power consumption data provided in [ 24 ]. W e obtain the admittance matrix of this system using a built- in function of the power flow simulator . It turns out that the absolute v alues of nonzero complex elements of the admittance matrix are between 1 . 86 and 40 . 06 , reminding the readers that a complex number’ s absolute value is its distance from zero in the complex plane. W e first consider the scenario that ev ery bus is equipped with a PMU. Assuming that the self admittance of bus 7 , i.e., 13 1 2 3 4 5 6 7 8 9 10 11 12 13 14 bus number 1 2 3 4 5 6 7 8 9 10 11 12 13 14 bus number 0 0.005 0.01 0.015 Fig. 5. The identification error when white Gaussian noise is added to both complex voltage and current measurements. Similar to the previous case, the errors are sufficiently small. the transformer bus, is known, Fig. 4 shows the identification error , defined as | Y − ˆ Y | , and the vertical color bar indicates the mapping of data v alues into colors. Hence, the color of a cell located at row i and column j represents the value of | Y [ i, j ] − ˆ Y [ i, j ] | . It can be seen that the identification error using measurements of 15 time slots ( K = 15 ) is quite small compared to the absolute v alue of the elements of Y , and this error does not vary much if we use more observations. W e next analyze the sensiti vity of the proposed algorithm to the measurement error which is typically introduced by the transducers. T o this end, white Gaussian noise with a signal-to-noise ratio of 125 is added to both complex voltage and current measurements. The signal-to-noise ratio is chosen such that the measurement accuracy lies within the reported range for existing PMU technology . Fig. 5 sho ws the absolute identification error for this case. Similar to the previous case, the errors are sufficiently small. In general, we observe that the identification error increases as we decrease the signal-to- noise ratio and it becomes really lar ge when the signal-to-noise ratio drops below 100 ; at this point we say that the Y matrix cannot be identified from data. B. Experimental data Next, we test the proposed IPF algorithm on the IEEE 5- bus benchmark sho wn in the left panel of Fig. 6 . The network contains 4 dynamic loads, 2 generators and 7 transmission lines. The system admittance matrix is obtained from [ 26 ] with the absolute values of nonzero complex elements of the admittance matrix are between 3 . 80 and 31 . 32 . Simulation is performed on O P A L - RT real time simulator [ 27 ] in our lab shown in the right panel of Fig. 6 . The generator is represented by a sixth order model with turbine gov ernors and excitation systems. Assume that the sampling frequency of PMU is 20 Hz and the dynamic loads change ov er time. W e collect PMU data (nodal voltages and injected currents) of all buses for 40 seconds. Note that the frequency fluctuates with changes of the load demands, which is consistent with field measurements; Fig. 6. Left: The network topology of the 5-bus benchmark systems. Right: The OP AL-R T real-time simulation platform that are used for data generation. PMU has 0 . 49% a verage total v ector error in simulation. Fig. 7 shows the identification error, defined as | Y − ˆ Y | (element- wise absolute value), and the vertical color bar indicates the mapping of data values into colors. One can see that the IPF is capable of identifying the topology correctly and estimating system parameters with small error from noisy data with the maximum parameter estimation error 0 . 2126 . 1 2 3 4 5 Bus Number 1 2 3 4 5 Bus Number 0 0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16 0.18 0.2 Fig. 7. The identified error of admittance matrix along buses for the 5-bus benchmark example we tried. C. Hidden node case Giv en a graph shown in Fig. 8 (left), if sensors are deployed at nodes M = { 1 , 2 , 6 , 7 , 8 , 9 , 12 } , we can use Kron reduction to obtain the graph shown in Fig. 8 (right). In this example, the hidden nodes are H = { 3 , 4 , 5 , 10 , 11 } , from which nodes { 3 , 5 , 10 , 11 } hav e degree less than 3 , and therefore, cannot 14 Fig. 8. An example of a 12 -bus network. Left: The original netw ork topology . Green circles denote the observed nodes while yellow ones denote the hidden nodes. Right: The topology of the condensed graph obtained from Kron reduction. The red lines denote the edges in G 1 and the grey ones denote the ones in G 2 . S t e p 1 : D e c om pos i t i on Cl i que T re e Com pos i t i on S t e p 3 : S t e p 2 : E s t i m a t e # of hi dde n nod e s + + Fig. 9. A step-by-step illustration of the proposed algorithm. Step 1 is described in Algorithm 2 that one separates the graph to a number of cliques and a forest. Step 2 is described in Algorithm 3 that one recov ers the admittance submatrix from every clique. Step 3 combines all the recovered graphs (together with admittance submatrices) to recover the whole admittance matrix. be identified from data. W e now illustrate ho w the proposed algorithm can identify the actual admittance matrix including Node 4 . The first step is to decompose the graph corresponding to the estimated ¯ Y matrix to two graphs using Algorithm 2 , one is a collection of cliques and the other one is a tree and some isolated nodes. The second step is to apply Algorithm 3 to identify the original admittance matrices for all cliques. The final step is to take the union of the cliques and the graph obtained in the second step. These steps are sho wn in Fig. 9 . Indeed, e ven the original graph does not satisfy the assump- tion on the degree of hidden nodes. W e can recov er a graph with smallest number of nodes that is a) consistent with ¯ Y and b) satisfies Assumption 3 . Consider another graph with the following topology in Fig. 10 , the admittance matrix has the following form 4 : 4 For notational simplicity , we used Y ij to denote the Y [ i, j ] element in this section. 1 6 8 5 4 7 2 3 Fig. 10. An example of a 8 -bus network, where 5 nodes are measured, and there exists three hidden nodes. Y = Y 11 0 0 0 0 − Y 16 0 0 Y 22 0 0 0 0 − Y 27 0 Y 33 0 0 0 − Y 37 0 Y 44 0 0 0 − Y 48 Y 55 0 0 − Y 58 Y 66 − Y 67 − Y 68 Y 77 0 Y 88 . Through direct computation, we have ¯ Y = Y 11 − X 11 Y 2 16 − X 12 Y 16 Y 27 − X 12 Y 16 Y 37 − X 13 Y 16 Y 48 − X 13 Y 16 Y 58 Y 22 − X 22 Y 2 27 − X 22 Y 27 Y 37 − X 23 Y 27 Y 48 − X 23 Y 27 Y 58 Y 33 − X 22 Y 2 37 − X 22 Y 37 Y 48 − X 23 Y 37 Y 58 Y 44 − X 33 Y 2 48 − X 33 Y 48 Y 58 Y 55 − X 33 Y 2 58 . For the ‘zero shunt element’ case, we follow Algorithm 3 . For starters, we compute γ γ = ∗ 0 0 0 0 ∗ 1 0 0 ∗ 0 0 ∗ 1 ∗ . It means that nodes 2 and 3 are connected to the same hidden node, and nodes 4 and 5 are connected to the same hidden node. Thus, it is straightforward to compute Y 22 , Y 33 , Y 44 , Y 55 , X 22 , X 23 , X 33 from ¯ Y and furthermore Y 27 = − Y 22 , Y 37 = − Y 33 , Y 48 = − Y 44 , Y 58 = − Y 55 . Next, we obtain a “new” ¯ Y by including two new hid- den boundary nodes as measured nodes in eq. ( 29 ). where X 12 Y 16 X 13 Y 16 can be computed from the corresponding row of ¯ Y and ¯ Y new , 11 = Y 11 − X 11 Y 2 16 + X 12 Y 16 X 13 Y 16 X 22 , X 23 X 23 , X 33 − 1 X 12 Y 16 X 13 Y 16 . It is noted that nodes 2 , 3 , 4 , 5 hav e no connections to non hidden boundary nodes, then they can be remov ed. ¯ Y new = " Y 11 − Y − 1 66 Y 2 16 Y − 1 66 Y 16 Y 67 Y 68 . Y 77 0 0 Y 88 − Y − 1 66 Y 67 Y 68 Y 67 Y 68 # . The last equality due to matrix inv erse lemma. W e can repeat 15 ¯ Y new = Y new , 11 0 X 12 Y 16 X 13 Y 16 X 22 , X 23 X 23 , X 33 − 1 diag { Y 22 , Y 33 , Y 44 , Y 55 } − Y 27 0 − Y 37 0 0 − Y 48 0 − Y 58 X 22 , X 23 X 23 , X 33 − 1 , (29) the procedure, by noting that the size of ¯ Y new is three, and furthermore solve for Y 11 , Y 77 , Y 88 , Y 66 , Y 67 , Y 68 , Y 16 , there- fore recover the original admittance matrix Y . V I . C O N C L U S I O N S This paper presents a framew ork for the in verse po wer flow problem which identifies the admittance matrix of a power system from synchronized voltage and current measurements pertaining to a subset of its buses. The algorithms proposed in this work can identify the graph topology together with its associated admittance matrix with guarantee for radial networks; and it can further jointly address state estimation and topology identification problems with theoretical guarantees, if certain conditions are met. These findings are supported by high-fidelity power flow simulations performed on standard test systems. In future work, we plan to extend our framework to three phase power flow models, which takes the mutual coupling between phases into account, develop efficient algorithms for identifying the admittance matrix of radial distrib ution systems with fe w measurement nodes, and analyze the sensiti vity of the identification results to non-stationary measurement errors. V I I . A C K N O W L E D G E M E N T W e thank Dr . W ei Zhou (Huazhong Univ ersity of Science and T echnology) for useful discussion. Y e Y uan was supported by the National Natural Science Foundation of China under Grant 92167201. A P P E N D I X Section IV describes in detail how to identify each max- imal clique in isolation and claims that admittance matrices obtained from these maximal cliques can be put together to construct the admittance matrix Y of the entire network. Here, we elaborate on how this is done. Theorem 1 implies that the graph G ( ¯ Y ) underlying the Kron-reduced admittance matrix ¯ Y contains a collection of edge-disjoint maximal cliques {C j ( ¯ Y ) } . The boundary ob- served nodes in each C j ( ¯ Y ) are connected to each other through a tree T j ( Y ) of hidden nodes in the original network G ( Y ) . As mentioned previously , e very boundary observed node in C j ( ¯ Y ) is connected to e xactly one hidden node in T j ( Y ) for , otherwise, there exists a loop in G ( Y ) , contradicting that G ( Y ) is a tree. Theorem 1 does not explicitly describe how these maximal cliques are connected in the graph G 2 . T wo maximal cliques can be connected in exactly one of three ways, corresponding to ho w the boundary observed nodes in the original graph G ( Y ) are connected through internal observed or hidden nodes. For each of these three cases, we now explain how to deriv e the (Kron-reduced) admittance matrix of ev ery maximal clique and how to put them together to construct the overall admittance matrix Y after each maximal clique has been identified in isolation. Consider two maximal cliques C j ( ¯ Y ) and C k ( ¯ Y ) in the graph G 2 in Theorem 1 . The three ways they can be connected are illustrated in Figure 11 . A single graph G ( Y ) may contain maximal cliques spanning cases 1, 2, and 3. Case 1: Figure 11 (a). C j ( ¯ Y ) and C k ( ¯ Y ) are disconnected when they are separated in the original graph G ( Y ) by internal observed nodes ( Y in equation (18) assumes there are no internal observed nodes). The Kron-reduced matrix ¯ Y takes the form ( × denotes nonzero entries): ¯ Y = × × × × × × × × × × × × × (30) where the submatrix of ¯ Y , shaded in red and denoted by W 11 , corresponds to the maximal clique C j ( ¯ Y ) and the submatrix of ¯ Y , shaded in blue and denoted by W 22 , corresponds to C k ( ¯ Y ) . The first rows and columns correspond to the internal measured nodes and their connection to the two blocks of boundary measured nodes. The admittance matrix ¯ Y j of the maximal clique C j ( ¯ Y ) (i.e., if C j ( ¯ Y ) were the entire graph) is: ¯ Y j := W 11 − diag 1 T W 11 , (31a) i.e., the diagonal entries of W 11 is modified so that ¯ Y j has zero row/column sums. Here 1 denotes a vector of all 1s of appropriate size. Similarly the admittance matrix ¯ Y k of the C k ( ¯ Y ) is: ¯ Y k := W 22 − diag 1 T W 22 . (31b) For each maximal clique C i ( ¯ Y ) in isolation, Section IV explains how to identify the admittance matrix Y i from its Kron reduction ¯ Y i . The admittance matrix Y i represents the topology and line admittances of the tree T i ( ¯ Y ) of hidden nodes and the observed nodes in C i ( ¯ Y ) . For Case 1, suppose the admittance matrices Y j and Y k identified from their Kron- 16 (a) Disconnected maximal cliques (b) Connected by a line in Y 11 , 22 (c) Connected by a boundary observed node Fig. 11. Examples of maximal cliques in ¯ G 2 . Shaded nodes are measured nodes and unshaded nodes are hidden nodes. (These examples do not satisfy Assumptions 3 .2). reductions ¯ Y j and ¯ Y k respectiv ely are: Y j =: Y j 11 , 22 Y j 12 , 21 0 Y j 22 , 11 Y j 22 , 12 Y j 22 , 22 (32a) Y k =: Y k 11 , 22 Y k 12 , 21 0 Y k 22 , 11 Y k 22 , 12 Y k 22 , 22 . (32b) Then the admittance matrix Y of the overall network is shown in ( 33 ), where the top-left submatrix Y 11 of Y has the same structure as ¯ Y in ( 30 ). The submatrices Y 11 , 11 and Y 11 , 12 correspond to the subgraph of internal measured nodes and their connecti vity to boundary measured nodes. They can be identified directly from the gi ven Kron-reduced admittance matrix ¯ Y of the overall network. The submatrices ˆ Y j 11 , 22 and ˆ Y k 11 , 22 are obtain from Y j 11 , 22 , Y k 11 , 22 and Y 11 , 12 by modifying the diagonal entries of Y j 11 , 22 and Y k 11 , 22 respectiv ely so that Y has zero row and column sums, i.e., " ˆ Y j 11 , 22 0 0 ˆ Y k 11 , 22 # := Y j 11 , 22 0 0 Y k 11 , 22 − diag 1 T Y 11 , 12 . Case 2: Figure 11 (b). C j ( ¯ Y ) and C k ( ¯ Y ) are connected by a line in Y 11 , 22 between two boundary observed nodes. The Kron-reduced matrix ¯ Y is of the form: ¯ Y = × × × × × × × × × × =: W 11 W 12 W T 12 W 22 (34) where the submatrix W 11 shaded in red corresponds to the maximal clique C j ( ¯ Y ) and the submatrix W 22 shaded in blue corresponds to C k ( ¯ Y ) . The admittance matrices ¯ Y j and ¯ Y k of the maximal cliques C j ( ¯ Y ) and C k ( ¯ Y ) respecti vely are giv en by ( 31 ). Suppose the admittance matrices Y j and Y k identified from their Kron-reduced matrices ¯ Y j and ¯ Y k respectiv ely are giv en by ( 32 ). Then the admittance matrix Y of the ov erall network is Y = ˆ Y j 11 , 22 W 12 Y j 12 , 21 0 0 0 W T 12 ˆ Y k 11 , 22 0 Y k 12 , 21 0 0 Y j 22 , 11 0 Y j 22 , 12 0 0 Y k 22 , 11 0 Y k 22 , 12 Y j 22 , 22 0 0 Y k 22 , 22 =: Y 11 Y 12 Y T 12 Y 22 , where the submatrix Y 11 of Y has the same structure as ¯ Y in ( 34 ). The submatrices ˆ Y j 11 , 22 and ˆ Y k 11 , 22 are obtain from Y j 11 , 22 , Y k 11 , 22 and W 12 by modifying the diagonal entries of Y j 11 , 22 and Y k 11 , 22 respectiv ely so that Y has zero row and column sums, i.e., ˆ Y j 11 , 22 := Y j 11 , 22 − diag ( W 12 1 ) ˆ Y k 11 , 22 := Y k 11 , 22 − diag 1 T W 12 . Case 3: Figure 11 (c). C j ( ¯ Y ) and C k ( ¯ Y ) are connected by a shared boundary observed node. The Kron-reduced matrix ¯ Y is of the form: ¯ Y = × × × × × × × × × × × × (35) where the submatrix W 11 shaded in red corresponds to the maximal clique C j ( ¯ Y ) and the submatrix W 22 shaded in blue corresponds to C k ( ¯ Y ) . The admittance matrices ¯ Y j and ¯ Y k of the maximal cliques C j ( ¯ Y ) and C k ( ¯ Y ) respecti vely are giv en by ( 31 ). The two maximal cliques C j ( ¯ Y ) and C k ( ¯ Y ) are connected by a single shared node. Suppose without loss of generality that the last row/column of ¯ Y j and the first row/column of 17 Y = Y 11 , 11 Y 11 , 12 0 0 0 0 ˆ Y j 11 , 22 0 Y j 12 , 21 0 0 0 0 ˆ Y k 11 , 22 0 Y k 12 , 21 0 0 Y j 22 , 11 0 Y j 22 , 12 0 0 Y k 22 , 11 0 Y k 22 , 12 Y j 22 , 22 0 0 Y k 22 , 22 =: Y 11 Y 12 Y T 12 Y 22 . (33) ¯ Y k correspond to the shared node. Suppose the admittance matrices Y j and Y k identified from ¯ Y j and ¯ Y k respectiv ely are Y j =: ˆ Y j 11 , 22 ˆ r j T 11 , 22 Y j 12 , 21 0 ˆ r j 11 , 22 d j 11 , 22 r j 12 , 21 0 Y j 22 , 11 Y j 22 , 12 Y j 22 , 22 Y k =: d k 11 , 22 ˆ r k 11 , 22 r k 12 , 21 0 ˆ r kT 11 , 22 ˆ Y k 11 , 22 Y k 12 , 21 0 Y k 22 , 11 Y k 22 , 12 Y k 22 , 22 , where ( ˆ r j 11 , 22 , d j 11 , 22 , r j 12 , 21 ) and ( ˆ r k 11 , 22 , d k 11 , 22 , r k 12 , 21 ) are the rows corresponding to the shared node. For the example in Figure 11 (c), ˆ Y j 11 , 22 is 2 × 2 and ˆ r j 11 , 22 is 1 × 2 ; both ˆ Y k 11 , 22 and ˆ r k 11 , 22 are 1 × 1 . T o construct the admittance matrix Y of the o verall netw ork, the rows in Y j and Y k corresponding to this shared node is combined into a single row: h ˆ r j 11 , 22 d j 11 , 22 + d k 11 , 22 ˆ r k 11 , 22 r j 12 , 21 r k 12 , 21 i The admittance matrix Y is then Y = ˆ Y j 11 , 22 ˆ r j T 11 , 22 0 Y j 12 , 21 0 0 0 ˆ r j 11 , 22 d j 11 , 22 + d k 11 , 22 ˆ r k 11 , 22 r j 12 , 21 r k 12 , 21 0 0 0 ˆ r kT 11 , 22 ˆ Y k 11 , 22 0 Y k 12 , 21 0 0 Y j 22 , 11 0 Y j 22 , 12 0 0 Y k 22 , 11 0 Y k 22 , 12 Y j 22 , 22 0 0 Y k 22 , 22 =: Y 11 Y 12 Y T 12 Y 22 , where submatrix Y 11 of Y has the same structure as ¯ Y in ( 35 ). R E F E R E N C E S [1] K. Dehghanpour, Z. W ang, J. W ang, Y . Y uan, and F . Bu, “ A survey on state estimation techniques and challenges in smart distribution systems, ” IEEE T ransactions on Smart Grid , vol. 10, no. 2, pp. 2312–2322, 2018. [2] C. Le Floch, S. Bansal, C. J. T omlin, S. J. Moura, and M. N. Zeilinger, “Plug-and-play model predictive control for load shaping and voltage control in smart grids, ” IEEE T ransactions on Smart Grid , vol. 10, no. 3, pp. 2334–2344, 2017. [3] O. Ardakanian, V . W . W ong, R. Dobbe, S. H. Low , A. von Meier , C. J. T omlin, and Y . Y uan, “On identification of distribution grids, ” IEEE T ransactions on Control of Network Systems , vol. 6, no. 3, pp. 950– 960, 2019. [4] P . Zhuang, R. Deng, and H. Liang, “False data injection attacks against state estimation in multiphase and unbalanced smart distribution systems, ” IEEE T ransactions on Smart Grid , vol. 10, no. 6, pp. 6000– 6013, 2019. [5] F . Dorfler and F . Bullo, “Kron Reduction of Graphs W ith Applications to Electrical Networks, ” IEEE T ransactions on Circuits and Systems I: Re gular P apers , vol. 60, no. 1, pp. 150–163, Jan 2013. [6] X. Li, H. V . Poor, and A. Scaglione, “Blind topology identification for power systems, ” in Proc. IEEE SmartGridComm , 2013. [7] V . Kekatos, G. B. Giannakis, and R. Baldick, “Grid topology iden- tification using electricity prices, ” in IEEE PES General Meeting — Confer ence Exposition , July 2014, pp. 1–5. [8] M. Babakmehr, M. G. Simoes, M. B. W akin, and F . Harirchi, “Com- pressiv e Sensing-Based T opology Identification for Smart Grids, ” IEEE T ransactions on Industrial Informatics , vol. 12, no. 2, pp. 532–543, April 2016. [9] Y . Y uan, X. T ang, W . Zhou, W . Pan, X. Li, H.-T . Zhang, H. Ding, and J. Goncalves, “Data dri ven discovery of cyber physical systems, ” Nature communications , vol. 10, no. 1, pp. 1–9, 2019. [10] D. Deka, S. Backhaus, and M. Chertko v , “Estimating distribution grid topologies: A graphical learning based approach, ” in 2016 P ower Systems Computation Confer ence (PSCC) , 2016, pp. 1–7. [11] Y . Liao, Y . W eng, M. Wu, and R. Rajagopal, “Distrib ution grid topology reconstruction: An information theoretic approach, ” in North American P ower Symposium , Oct 2015, pp. 1–6. [12] Y . Liao, Y . W eng, and R. Rajagopal, “Urban Distribution Grid T opology Reconstruction via Lasso, ” in IEEE PES General Meeting , July 2016, pp. 1–6. [13] J. Y u, Y . W eng, and R. Rajagopal, “Patopa: A data-driven parameter and topology joint estimation frame work in distribution grids, ” IEEE T ransactions on P ower Systems , vol. 33, no. 4, pp. 4335–4347, 2017. [14] S. Park, D. Deka, S. Backhaus, and M. Chertkov , “Learning with end- users in distribution grids: T opology and parameter estimation, ” IEEE T ransactions on Control of Network Systems , vol. 7, no. 3, pp. 1428– 1440, 2020. [15] H. Zhu and G. B. Giannakis, “Sparse ov ercomplete representations for efficient identification of power line outages, ” IEEE T ransactions on P ower Systems , vol. 27, no. 4, pp. 2215–2224, Nov ember 2012. [16] D. Deka, S. Backhaus, and M. Chertkov , “Structure learning in power distribution networks, ” IEEE T ransactions on Contr ol of Network Sys- tems , vol. 5, no. 3, pp. 1061–1074, September 2018. [17] G. Cavraro and V . Kekatos, “Graph algorithms for topology identifi- cation using power grid probing, ” IEEE contr ol systems letters , vol. 2, no. 4, pp. 689–694, 2018. [18] G. Cavraro, V . Kekatos, and S. V eeramachaneni, “V oltage analytics for power distrib ution network topology verification, ” IEEE T ransactions on 18 Smart Grid , vol. 10, no. 1, pp. 1058–1067, 2017. [19] H. Li, Y . W eng, Y . Liao, B. Keel, and K. E. Brown, “Distribution grid impedance & topology estimation with limited or no micro-pmus, ” International Journal of Electrical P ower & Ener gy Systems , vol. 129, p. 106794, 2021. [20] S. Wright and J. Nocedal, “Numerical optimization, ” Springer Science , vol. 35, pp. 67–68, 1999. [21] Y . Y uan, O. Ardakanian, S. Lo w , and C. T omlin, “On the in verse power flow problem, ” arXiv preprint , 2016. [22] R. D. Zimmerman, C. E. Murillo-Sanchez, and R. J. Thomas, “MA T - PO WER: Steady-State Operations, Planning, and Analysis T ools for Power Systems Research and Education, ” IEEE T ransactions on P ower Systems , vol. 26, no. 1, pp. 12–19, Feb 2011. [23] M. Grant and S. Boyd, “CVX: Matlab software for disciplined conv ex programming, version 2.1, ” http://cvxr .com/cvx , Mar. 2014. [24] IEEE, “14 Bus T est Case, ” https://www .ee.washington.edu/research/ pstca/pf14/pg tca14bus.htm . [25] O. Ardakanian, Y . Y uan, R. Dobbe, A. von Meier, S. Low , and C. T omlin, “Event detection and localization in distribution grids with phasor measurement units, ” in IEEE P ower Energy Society General Meeting , July 2017, pp. 1–5. [26] “IEEE 5 bus test system data, ” http://shodhganga.inflibnet.ac.in/ bitstream/10603/26549/14/14 appendix.pdf , [Online; accessed 20-Jan- 2019]. [27] “OP AL-R T real-time digital simulator , ” https://www .opal- rt.com/ simulator- platform- op5600/ , [Online; accessed 20-Jan-2019]. Y e Y uan (M’13, SM’20) recei ved the B.Eng. degree from the Department of Automation, Shanghai Jiao T ong University , Shanghai, China, in 2008, and the M.Phil. and Ph.D. degrees from the Department of Engineering, University of Cambridge, Cambridge, U.K., in 2009 and 2012, respecti vely . He has been a Full Professor at the Huazhong Uni versity of Science and T echnology , Wuhan, China since 2016. Prior to that, he was a Postdoctoral Researcher at UC Berkeley , a Junior Research Fellow at Darwin Col- lege, Uni versity of Cambridge. His research interests include system identification and optimization with applications to cyber- physical systems. Steven Low (F’08) is the F . J. Gilloon Professor of the Department of Computing & Mathematical Sciences and the Department of Electrical Engi- neering at Caltech and an Honorary Professor of the University of Melbourne. Before that, he was with A T&T Bell Laboratories, Murray Hill, NJ, and the Univ ersity of Melbourne, Australia. He has held honorary/chaired professorship in Australia, China and T aiwan. He was a co-recipient of IEEE best paper awards, an awardee of the IEEE INFOCOM Achiev ement A ward and the ACM SIGMETRICS T est of Time A ward, and is a Fello w of IEEE, ACM, and CSEE. He was well- known for work on Internet congestion control and semidefinite relaxation of optimal power flo w problems in smart grid. His research on networks has been accelerating more than 1TB of Internet traffic every second since 2014. His research on smart grid is providing large-scale cost effecti ve electric vehicle charging to workplaces. He recei ved his B.S. from Cornell and PhD from Berkeley , both in EE. Omid Ardakanian (M’15) is an Assistant Professor at the University of Alberta, Canada. He received his B.Sc. from Sharif University of T echnology in 2009, and M.Math. and Ph.D. from the University of W a- terloo in 2011 and 2015. From 2015 to 2017, he was an NSERC Postdoctoral Fellow at UC Berkeley and the University of British Columbia. He recei ved best paper awards at A CM e-Energy , ACM BuildSys, and IEEE PES General Meeting. His research focuses on the design and implementation of intelligent networked systems. Claire J. T omlin (F’10) is the Charles A. Desoer Professor of Engineering in the Department of Elec- trical Engineering and Computer Sciences (EECS), Univ ersity of California Berkeley (UC Berkeley). She was an Assistant, Associate, and Full Professor in Aeronautics and Astronautics at Stanford Uni ver- sity from 1998 to 2007, and in 2005, she joined UC Berkeley . Claire works in the area of control theory and hybrid systems, with applications to air traffic management, UA V systems, energy , robotics, and systems biology . She is a MacArthur Foundation Fellow (2006), an IEEE Fellow (2010), and in 2017, she was awarded the IEEE Transportation T echnologies A ward. In 2019, Claire was elected to the National Academy of Engineering and the American Academy of Arts and Sciences.
Original Paper
Loading high-quality paper...
Comments & Academic Discussion
Loading comments...
Leave a Comment