From State Estimation to Network Reconstruction

TL;DR

This paper introduces compressed sensing algorithms for power grid topology and parameter reconstruction, enabling efficient monitoring of grid structure and state even with limited measurements and uncertain information.

Contribution

It presents novel compressed sensing methods that leverage network sparsity and prior knowledge to reconstruct power grid topology and parameters from minimal measurements.

Findings

Algorithms effectively reconstruct grid topology and parameters.

Methods outperform traditional state estimation in resource efficiency.

Applicable to uncertain or incomplete grid information.

Abstract

We develop methods to efficiently reconstruct the topology and line parameters of a power grid from the measurement of nodal variables. We propose two compressed sensing algorithms that minimize the amount of necessary measurement resources by exploiting network sparsity, symmetry of connections and potential prior knowledge about the connectivity. The algorithms are reciprocal to established state estimation methods, where nodal variables are estimated from few measurements given the network structure. Hence, they enable an advanced grid monitoring where both state and structure of a grid are subject to uncertainties or missing information.

Equations

$$P_{sr}=\frac{|V_s||V_r|}{X_{sr}}\sin ({\phi }_s-{\phi }_r)\approx X_{sr}^{-1}({\phi }_s-{\phi }_r)$$

$$P_s=\sum _{r=1}^N{P}_{sr}=\sum _{r=1}^N{X}_{sr}^{-1}({\phi }_s-{\phi }_r).$$

$$\mathbf{\phi }$$

$$\mathbf{B}\mathbf{\phi }=\mathbf{P},$$

$$B_{nk}=\left\{\begin{array}[]{lll}\displaystyle\sum\nolimits_{j=1}^{N}X_{nj}^{-1}&\mbox{if}&k=n;\\[5.69054pt] -X_{nk}^{-1}&\mbox{if}&k\neq n.\end{array}\right.$$

$$z_i=h_i(\mathbf{x})$$

$$\mathbf{z}=\mathbf{h}(\mathbf{x})$$

$$\mathbf{x}\in {\mathbb{R}}^{N\times 1}\text{mbox}isthesystemstate,$$

$$\mathbf{z}\in {\mathbb{R}}^{M\times 1}\text{mbox}isavectorofmeasuredquantitiesand$$

$$\mathbf{h}\text{mbox}isavectorfunction,i.e.amapping{\mathbb{R}}^N\to {\mathbb{R}}^M$$

$$P_{sr}=X_{sr}^{-1}({\phi }_s-{\phi }_r)+e_{sr}$$

$$P_s=\sum _{r=1}^N{P}_{sr}+e_s^{\prime }$$

$$\mathbf{B}\mathbf{\phi }=\mathbf{P}.$$

$$\sum _k{B}_{rk}{\phi }_k(t_m)=P_r(t_m)$$

$$\sum _k{B}_{rk}=0.$$

$$\underset{=:\mathbf{\Phi }}{\underbrace{\left(\begin{array}{ccc}{\phi }_1(t_1)& \cdots & {\phi }_N(t_1)\\ ⋮& & ⋮\\ {\phi }_1(t_M)& \cdots & {\phi }_N(t_M)\\ 1& \cdots & 1\end{array}\right)}}\underset{=:{\mathbf{B}}_r}{\underbrace{\left(\begin{array}{c}{B}_{r1}\\ ⋮\\ B_{rN}\end{array}\right)}}=\underset{=:{\mathcal{P}}_r}{\underbrace{\left(\begin{array}{c}{P}_r(t_1)\\ ⋮\\ P_r(t_M)\\ 0\end{array}\right)}}$$

$$\Vert \mathbf{\Phi }{\mathbf{B}}_r-{\mathcal{P}}_r{\Vert }_2^2=\sum _{k=1}^{M+1}(\mathbf{\Phi }{\mathbf{B}}_r-{\mathcal{P}}_r{)}_k^2.$$

$$\Vert {\mathbf{B}}_r{\Vert }_1=\sum _{k=1}^N|{\mathbf{B}}_{rk}|$$

$$K_{rs}:=\left\{\begin{array}[]{l l}1&\;\mbox{if the value of $B_{rs}$ is known for line $(r,s)$}\\ 0&\;\mbox{otherwise}.\end{array}\right.$$

$${\mathbf{\Phi }}_r^{\mathrm{red}}{\mathbf{B}}_r^{\mathrm{red}}={\mathcal{P}}_r^{\mathrm{red}},$$

$${\mathcal{P}}_r^{\mathrm{red}}={\mathcal{P}}_r-\sum _{c=1}^N{B}_{rc}{K}_{rc}\left(\begin{array}{c}{\phi }_c(t_1)\\ ⋮\\ {\phi }_c(t_M)\\ 1\end{array}\right)$$

$$P_{sr}(t_{m})=X_{sr}^{-1}\big{(}\varphi_{r}(t_{m})-\varphi_{s}(t_{m})\big{)}$$

$$B_{rs}=B_{sr}=-X_{sr}^{-1}$$

$$\underset{{\mathbf{B}}_r}{\min }\Vert \mathbf{\Phi }{\mathbf{B}}_r-{\mathcal{P}}_r{\Vert }_2^2,$$

$${\mathbf{B}}_r=({\mathbf{\Phi }}^T\mathbf{\Phi }{)}^{-1}{\mathbf{\Phi }}^T{\mathcal{P}}_r$$

$${\mathbf{B}}_r={\mathbf{B}}_r^{(\mathrm{sp})}+\mathbf{W}\mathbf{y},$$

$$\min_{\boldsymbol{s}}\boldsymbol{1}^{T}\boldsymbol{s}\quad\mbox{such that}\quad\begin{array}[]{l}\boldsymbol{B}_{r}^{\rm(sp)}+\boldsymbol{W}\boldsymbol{y}\preceq\boldsymbol{s}\\ \boldsymbol{B}_{r}^{\rm(sp)}+\boldsymbol{W}\boldsymbol{y}\succeq-\boldsymbol{s}\\ \end{array},$$

$$\underset{\mathbf{x}}{\min }{\mathbf{f}}^T\mathbf{x}\text{mbox}suchthat\mathbf{A}\mathbf{x}⪯\mathbf{b}.$$

$$\mathbf{x}$$

$$\mathbf{b}$$

$$\Vert {\mathbf{B}}_r{\Vert }_0=\text{mbox}numberofnon-zeroentriesof{\mathbf{B}}_r.$$

$$|B_{r,k}^{\mathrm{reconstructed}}-B_{r,k}^{\mathrm{true}}|<ϵ,$$

Full text

I Introduction

The secure and reliable operation of complex power grids requires a precise knowledge of the grid topology and the state of connected generation and transmission elements. A central tool is ‘state estimation’, i.e the estimation of the current state of the nodes from few, possibly noisy measurements and the grid data. Various methods and algorithms have been developed in this field, see e.g. [1] for a review. A prime example is the estimation of the voltage phase angles from power flow measurements (see e.g. [2] for basic textbook examples). State estimation algorithms typically assume that the network structure, i.e. the location and parameters of transmission lines and transformers, are known. However, parts of the information might be lacking or uncertain, for example for contingencies or attacks, such that the question arises: Is it possible to reconstruct the network structure from local measurements only?

In this article we develop two algorithms that enable the faithful reconstruction of the network structure or parts of it from the measurement of nodal variables. Special attention is paid to the efficiency of these methods: We show how to reduce the number of measurements using methods from compressed sensing exploiting the structural properties of power grids. We here focus on the DC approximation to formulate the fundamental ideas of network reconstruction most clearly.

The presented algorithms can find several applications in the monitoring and operation of (smart) power grids. While it is unlikely that the entire network structure is unknown, situations regularly occur where parts of the information is lacking. Manual switches still exist in many distribution grids, transmission lines can get lost in contingency cases and targeted attacks can cause entire regions of a grid to collapse. We present an example where the connectivity in two regions of the grid is a priori unknown, but can be reconstructed from a single snapshot of the nodal variables. In future power grids, the network’s structure will become even more variable: The transmission of high-voltage direct current (HVDC) lines in hybrid power grids can be actively controlled and the effective line parameters can be regulated with flexible alternating current transmission system (FACTS) devices. Then it becomes essential to include lacking information or uncertainties of the grid topologies in any grid monitoring tool. The methods presented in this article shall contribute to the development of such advanced monitoring algorithms [3, 4, 5].

II The DC approximation

The DC approximation describes the flow of real power flow in AC power grids in a linearized way. It is derived from the full nonlinear load flow equations under the following additional simplifying assumptions [6, 2, 7, 8]:

• •

All bus voltage magnitudes $|V_{s}|$ are close to 1 p.u. and approximated to be exactly $|V_{s}|=1$ p.u.

• •

Voltage angle differences across branches are small, such that $\sin(\varphi_{s}-\varphi_{r})\approx\varphi_{s}-\varphi_{r}$ and $\cos(\varphi_{s}-\varphi_{r})\approx 1$.

• •

Branches can be considered lossless. In particular, branch resistances and charging capacitances are negligible.

The approximate real power flow from node $s$ to node $r$ over a transmission line with reactance $X_{sr}$ is then given by

[TABLE]

where $\varphi_{s}$ is the phase angle at bus $s$. The real power balance at the node $s$ reads

[TABLE]

For notational convenience we summarize these equations in matrix form. We define the vector of all voltage angles and power injections,

[TABLE]

where $N$ is the number of buses or nodes in the grid and the superscript ‘$T$’ denotes the transpose of a vector or matrix. The DC approximation then yields the condition

[TABLE]

where the nodal susceptance matrix $\boldsymbol{B}\in\mathbb{R}^{N\times N}$ has elements

[TABLE]

The matrix $\boldsymbol{B}$ is a Laplacian matrix, which has one zero eigenvalue with eigenvector $(1,1,\ldots,1)^{T}$ [9]. This eigenvector represents a global shift of all voltage angles without physical significance. To avoid this complication, one often fixes the voltage angle at a reference bus and excludes this bus from the calculation, i.e. the corresponding row and column will be missing in $\boldsymbol{B}$. We finally note that mathematically equivalent models of flow are used to describe hydraulic networks [10] or vascular networks of plants [11].

III State Estimation

State estimation is a mathematical procedure by which the state of an electric power system is extracted from a set of measurements [1]. In general, any measurement can be expressed as a function of the system state. Let

[TABLE]

denote a measured quantity, where $\boldsymbol{x}$ is the system state and $h_{i}$ is a function specific to the measured quantity $z_{i}$. Assume that $M$ measurements are taken. Then, all measurements can be written in compact form

[TABLE]

where

[TABLE]

Typically more measurements than the number of state variables to be determined are taken, i.e. $M>N$, such that equation (6) is overdetermined. It is then solved in a least-squares fashion, i.e. one calculates $\boldsymbol{x}$ such that $\|\boldsymbol{z}-\boldsymbol{h}(\boldsymbol{x})\|_{2}^{2}$ is minimized. For underdetermined equations (6) typically many solutions exist.

It is often helpful to work with a simplified DC approximation model for the measurement equations in analyzing the inherent limitations of various methods related solely to the measurement configuration [12]. The state of the grid is described by the voltage phase angles $\boldsymbol{x}=\boldsymbol{\varphi}$ and shall be estimated from measurements. As explicated above, the real power flow measured from bus $s$ to $r$ can be approximated by

[TABLE]

where $\varphi_{s}$ is the phase angle at bus $s$ and $e_{sr}$ is a measurement error. Similarly, a power injection measurement

[TABLE]

at a given bus $s$ can be expressed as a sum of flows along incident branches to that bus. The measurement vector $\boldsymbol{z}$ thus consists of a subset of the power flows $P_{sr}$ and power injections $P_{s}$. In the DC approximation the vector function $\boldsymbol{h}$ are all linear, which strongly simplifies the problem such that it is especially suitable to introduce new concepts such as network reconstruction.

IV Network Reconstruction

Let us come back to the defining equation for the DC approximation

[TABLE]

In state estimation, we assume that the nodal phase angles $\boldsymbol{\varphi}$ are difficult to measure and shall be estimated from the knowledge of the two other quantities in this equation, $\boldsymbol{B}$ and $\boldsymbol{P}$. But we can also reverse the problem and ask what happens if we don’t know the nodal susceptance matrix $\boldsymbol{B}$. Can we efficiently reconstruct its entries from measurements of the two remaining quantities $\boldsymbol{P}$ and $\boldsymbol{\varphi}$? This problem may appear hypothetical at a first glance, as the network structure and parameters are generally known. But there are situations where at least parts of the network structure are unknown: Switches can be open or closed, transmission lines can undergo failures, or plans may be inaccurate. Even more, in times of war or terrorism we may be extremely unsure about the physical integrity of parts of the grid. Can we reconstruct the structure of the grid from measurements only?

In the following we develop a mathematical theory of network reconstruction. To introduce the method, we start from the hypothetical situation that $\boldsymbol{B}$ is completely unknown, whereas $\boldsymbol{P}$ and $\boldsymbol{\varphi}$ can be measured perfectly. We then study more realistic cases, where only parts of the grid topology are unknown.

We first notice that a single measurement of $\boldsymbol{P}$ and $\boldsymbol{\varphi}$ will obviously not be enough to reconstruct the entire grid topology. However, we may repeat our measurement at different times $t_{1},t_{2},\ldots,t_{M}$ such that we obtain a large number of conditions on the entries of $\boldsymbol{B}$ of the form

[TABLE]

which hold for all rows $r=1,\ldots,N$ and all measurement time steps $m=1,\ldots,M$. Furthermore we know that the row-sums of the nodal susceptance matrix (4) vanish,

[TABLE]

In total we thus have $N\times(M+1)$ conditions which we can use to reconstruct the unknown entries of $\boldsymbol{B}$.

IV-A Row-wise reconstruction

In a first approach we aim to reconstruct the nodal admittance matrix row-by-row (cf. [13, 14, 15, 16]). To simplify notation we collect all measurement conditions (10) and the condition (11) for a given row $r$ and rewrite them in matrix form

[TABLE]

Here, $\boldsymbol{B}_{r}\in\mathbb{R}^{N\times 1}$ denotes the transpose of the $r$th row of the nodal admittance matrix $\boldsymbol{B}$, which we want to reconstruct. The matrix $\boldsymbol{\Phi}\in\mathbb{R}^{(M+1)\times N}$ summarizes all measurement results for the voltage phase angles at all time steps, while the vector $\boldsymbol{\mathcal{P}}_{r}\in\mathbb{R}^{(M+1)\times 1}$ summarizes the values of the power injection only for node $r$. The last row of the matrix equation represents the condition (11) for the row sum.

Now we have to distinguish, whether equation (12) is (over)determined or underdetermined. If the number of linearly independent equations is larger or equal to $N$ we generally have enough information to directly compute $\boldsymbol{B}_{r}$. In particular there is exactly one solution if ${\rm rank}(\bar{\boldsymbol{\Phi}}|\bar{\boldsymbol{\mathcal{P}}_{r}})={\rm rank}(\bar{\boldsymbol{\Phi}})=N$, which can be calculated using Gaussian elimination. Nevertheless, matrix inversion can be numerically ill-conditioned, especially for systems with large $N$ [17]. If the data includes some measurement noise, the system of equations (12) can be overdetermined when measuring $M+1>N$ times such that it is solvable in a least squares fashion, i.e. we have to determine the vector $\boldsymbol{B}_{r}\in\mathbb{R}^{N\times 1}$ minimizing the 2-norm

[TABLE]

To efficiently reconstruct the grid topology, we want to rely on as few measurements as possible. If $M+1<N$, we generally have an underdetermined set of equations which admits many possible solutions. Is it possible to obtain the correct solution also in this case? Do we have more information about $\boldsymbol{B}_{r}$ which we can exploit? Indeed, we know that a power grid is typically very sparse – a single substation is connected to only few other substations. Thus we choose the one solution to equation (12) which minimizes the number of non-zero entries of $\boldsymbol{B}_{r}$. Unfortunately, the direct minimization is computationally hard in general. In 2006 Candes, Romberg, Tao and Donoho showed that an efficient reconstruction is nevertheless possible using a convex surrogate for sparsity: the 1-norm [18, 19]. Under weak conditions, the correct sparse solution can be calculated efficiently by minimizing the 1-norm

[TABLE]

subject to the constraint (12). This problem can be mapped to a linear program which we explicate in the following section.

IV-B Partial reconstruction

In applications we will rarely encounter the situation that the structure of a power grid is completely unknown. In a typical application we know the value of $B_{rs}$ of all transmission lines $(s,r)$ in one part of the grid and have to reconstruct only the remaining entries of the matrix $\boldsymbol{B}$ – for instance we might want to monitor the position of manual switches at remote places of the grid.

To keep track of our knowledge about the network, we define the matrix $\boldsymbol{K}\in\mathbb{R}^{N\times N}$ with entries

[TABLE]

We again proceed row by row. For the $r$th row we reduce equation (12) to

[TABLE]

where

[TABLE]

and $\boldsymbol{\Phi}_{r}^{\rm red}$ is the submatrix of $\boldsymbol{\Phi}_{r}$ obtained by deleting all columns $c$ for which $K_{rc}=1$ and $\boldsymbol{B}_{r}^{\rm red}$ is the submatrix of $\boldsymbol{B}_{r}$ obtained by deleting all rows $r^{\prime}$ for which $K_{r,r^{\prime}}=1$. The dimension of the reduced linear system of equations (16) is smaller than in the original problem (12), such that the correct solution can typically be found from fewer measurements.

IV-C Reconstruction from power flow measurements

The reconstruction scheme introduced above requires the knowledge of the voltage phase angle at all nodes of the grid. This data can in principle be obtained using phasor measurement units (PMUs), but these are typically very expensive. The measurement of other quantities such as real power flows is typically much simpler.

Fortunately, we usually do not have to reconstruct the entire grid in most cases. In a typical application we know a lot about the grid and have to reconstruct only parts of the matrix $\boldsymbol{B}$. Then the measurement of power injections and real power flows of a certain subset of transmission lines can be enough to perform the reconstruction of the remaining data using the methods described above. Technically, this amounts to combining methods of classical state estimation and network reconstruction.

So assume that we can measure the power injection $P_{n}$ for all nodes $n\in\{1,\ldots,N\}$ at the time steps $t_{1},\ldots,t_{M}$. In addition we have knowledge about the transmission line parameters, in particular $X_{sr}$, and measure the real power flow

[TABLE]

for a subset $\mathcal{L}$ of all transmission lines at all time steps $t_{1},\ldots,t_{M}$. If this set of equations is fully determined or even overdetermined we reconstruct the state vector $\boldsymbol{\varphi}(t_{m})$ as in the classical state estimation problem described in section III. As before the resulting estimates for the nodal phase angles for all time steps $t_{1},\ldots,t_{M}$ are then summarized in the matrix $\boldsymbol{\Phi}$. In addition we have partial knowledge about the nodal susceptance matrix, in particular we already know the entries

[TABLE]

for all $(s,r)\in\mathcal{L}$. This information can then be used to reconstruct the remaining entries of $\boldsymbol{B}$ as described in the previous section.

V Implementation

V-A Least-squares solution

If we have many measurements or prior knowledge available, then the system of equations (12) can be be overdetermined. In this case equation must be solved in a least squares fashion, i.e. we have to solve

[TABLE]

which leads to [15]

[TABLE]

provided that $\boldsymbol{\Phi}$ is full rank such that the inverse exists. Least-squares solutions are already implemented in many numerical solvers, for instance in the MATLAB function mldivide.

V-B Minimizing the 1-norm

The key to an efficient reconstruction of $\boldsymbol{B}$ from an underdetermined system of equations is the minimization of the 1-norm [18, 19]. This problem can be mapped to a linear program which can be solved efficiently, i.e. the computation time scales at most polynomially with the input length.

So we consider the linear system of equations (12) and assume that it is solvable and underdetermined. The solutions span an affine subspace of of $\mathbb{R}^{N}$ with dimension $D=N-{\rm rank}({\boldsymbol{\Phi}})$. All solutions can be written as

[TABLE]

where the columns of the matrix $\boldsymbol{W}\in\mathbb{R}^{N\times D}$ form a basis for the kernel (nullspace) of ${\boldsymbol{\Phi}}$ and $\boldsymbol{y}\in\mathbb{R}^{D}$ is a vector of parameters. $\boldsymbol{B}_{r}^{\rm(sp)}$ is a specific solution to equation (12); in the implementation it is obtained using the MATLAB function mldivide. We now search for the vector $\boldsymbol{y}$ such that the 1-norm (14) assumes its minimum. This optimization problem is rewritten as a linear program [20]

[TABLE]

where $\boldsymbol{s}\in\mathbb{R}^{N}$ is an auxiliary variable, $\boldsymbol{1}\in\mathbb{R}^{N}$ is a vector of ones and $\preceq$ and $\succeq$ denotes entry-wise comparison. In the standard form of many commercial solvers the optimization problem reads

[TABLE]

using the further auxiliary variables

[TABLE]

where \mbox{\mbox{1}!\mbox{l};}\in\mathbb{R}^{N\times N} is the identity matrix and $\boldsymbol{0}\in\mathbb{R}^{D}$ is a vector of zeros.

V-C Row-Wise Reconstruction

We first describe the implementation of the row-wise reconstruction of $\boldsymbol{B}$ from the measurement data without any prior knowledge about the entries of $\boldsymbol{B}$, whose basic idea has been outlined in section IV-A. The algorithm to solve this problem is illustrated by the flow-chart in Fig. 1. We implement this algorithm in Matlab, the central part of the program code is listed in Fig. 2.

First data is collected and the matrix $\boldsymbol{\Phi}$ and the vectors $\boldsymbol{\mathcal{P}}_{r}$ are formed. Then the problem is solved depending on the rank of the matrix ${\boldsymbol{\Phi}}$. If the system is (over)determined, we directly solve it for $\boldsymbol{B}$ using the Matlab-function mldivide. If the system is underdetermined we proceed row-by-row and reconstruct $\boldsymbol{B}_{r}$ as described in section V-B using the function null to compute the matrix $\boldsymbol{W}$ and the function linprog to solve the linear program (24).

V-D Iterative reconstruction

In practice we typically have some prior knowledge about the entries of the matrix $\boldsymbol{B}$ as outlined in section IV-B. An algorithm to perform the reconstruction exploiting the prior knowledge is illustrated by the flow-chart in Fig. 3

Most interestingly, we can also use this algorithm to greatly improve the convergence of the reconstruction algorithm. For intermediate values of $M$ we typically face the situation that some rows are successfully reconstructed while others are not. If we have not reconstructed $\boldsymbol{B}$ successfully, but gained further knowledge, we may use it in the following. Having successfully reconstructed the $r$th row of $\boldsymbol{B}$, we also know the entries $B_{r^{\prime},r}=B_{r,r^{\prime}}$ for all other rows $r^{\prime}$ due to the symmetry of the matrix. We thus propose to perform the reconstruction iteratively. In each step of the reconstruction process we use the initial knowledge about the entries of the matrix $\boldsymbol{B}$ and also the knowledge gained in previous steps.

The algorithm shown in Fig. 3 starts with the input of the measurement results and the prior knowledge of the system which is encoded in the matrices $\boldsymbol{B}$ and $\boldsymbol{K}$. It solves the reconstruction problem iteratively, processing row-by-row in each iteration. For each row, we first calculate the reduced matrices $\boldsymbol{\Phi}_{r}^{\rm red}$ and $\boldsymbol{\mathcal{P}}_{r}^{\rm red}$ as defined in section IV-B. Then we attempt to reconstruct the vector $\boldsymbol{B}_{r}^{\rm red}$ as described above: If the problem is (over)determined we solve the linear system of equations using the Matlab-function mldivide, if it is underdetermined we minimize the 1-norm $\|\boldsymbol{B}_{r}^{\rm red}\|_{1}$ subject to the constraint (16). If the reconstruction was successful, we add the results to the matrix $\boldsymbol{B}$ and update the knowledge matrix $\boldsymbol{K}$. To facilitate the bookkeeping we define two matrices $\boldsymbol{K}$ and $\boldsymbol{K}^{\rm new}$: We set $\boldsymbol{K}=\boldsymbol{K}^{\rm new}$ at the beginning of each iteration and only modify $\boldsymbol{K}^{\rm new}$ during the step. The iteration stops when the reconstruction has been successfully completed, i.e. $K^{\rm new}_{rc}=1$ for all $r,c=1\ldots,N$, or when no further progress has been made, i.e. $K^{\rm new}$ has not been modified during the last step.

One unsolved problem remains: If we do not know the matrix $\boldsymbol{B}$ a priori, how can we know if the reconstruction of the $r$th row has been successful? In the underdetermined case we cannot decide whether the reconstructed values $B_{rc}$ are definitely correct – but we can decide if they are reasonable in terms of the connectivity of the grid. Power grids are generally very sparse: A single substation is connected to only few other substations. Hence we seek for a solution which is most sparse, i.e. which minimizes

[TABLE]

However, we cannot minimize $\|\boldsymbol{B}_{r}\|_{0}$ directly as this is a computationally hard problem. The ingenious contribution of Candes, Romberg, Tao and Donoho [18, 19] was the proof that the most sparse solution can be obtained by minimizing $\|\boldsymbol{B}_{r}\|_{1}$ with high probability if a sufficient amount of measurement data is available. Otherwise we find a solution which minimizes the 1-norm, but has a large number of non-zero entries $\|\boldsymbol{B}_{r}\|_{0}$. In algorithm in Fig. 3 we thus adopt the definition that the reconstruction of the $r$th row is assumed to be successful if (a) the linear system is overdetermined or (b) the reconstructed row $\boldsymbol{B}_{r}$ is sufficiently sparse, i.e. $\|\boldsymbol{B}_{r}\|_{0}\leq d_{\rm max}$ for some upper limit $d_{\rm max}$. A typical value used in the following is $d_{\rm max}=15$.

VI Applications and Performance

VI-A Reconstruction from time series

We demonstrate the applicability of network reconstruction from time series measurements for a test grid taken from [21] illustrated in Figure 4 (a). This data set includes hourly data for demand and generation $P_{k}(t)$ for one year and all $N=117$ nodes. The nodal voltage angles $\varphi_{k}(t)$ are then obtained by solving the DC approximation (3) for all time steps. The resulting time series of the power injections and angles are shown in Figure 4 (b) for two nodes as an example.

The algorithm shown in Figure 3 can now reconstruct the entire network topology - i.e. all entries of the nodal susceptance matrix $\boldsymbol{B}$ - from the time series data. To evaluate the performance of the algorithm we vary the amount of input data. For each value of $M$ we run the algorithm and count the number of successfully reconstructed entries and the number of reconstruction errors. To account for small numerical errors, we say that an entry $(r,k)$ of the nodal susceptance matrix is successfully reconstructed if

[TABLE]

where the numerical tolerance is chosen as $\epsilon=10^{-3}\times\max_{ij}|B_{ij}|$.

Figure 4 (c) shows that $M_{\rm min}=55$ measurements are sufficient to faithfully reconstruct the entire nodal susceptance matrix $\boldsymbol{B}$. Most importantly, a faithful reconstruction is possible in the strongly underdetermined case: The minimum number of measurements $M_{\rm min}=55$ is less than half of the matrix dimension $N=117$. This is possible because we can exploit the sparsity of the matrix $\boldsymbol{B}$ as an additional structural information.

VI-B Minimum Measurement Requirements and Scaling

The algorithms presented in figures 1 and 3 allow for a reconstruction of the network topology also in the underdetermined case. But how many measurements $M_{\rm min}$ are needed for a faithful reconstruction and how does this number scale with the grid size $N$?

To systematically study the efficiency of the reconstruction algorithms, we consider various test grids of different size taken from [22, 23] and use randomized synthetic data for the time series. We consider two different types of time series data. First we draw the voltage phase angles $\varphi_{k}(t)$ uniformly at random from the interval $[-\pi/8,+\pi/8]$ (type I). This scenario is not realistic but close to the original mathematical work [18, 19], where random Gaussian sampling vectors are considered. Second, we choose the power injections $P_{k}(t)$ at random from a normal distribution with mean zero and standard deviation as in the original test grid (type II). For both types and both reconstruction algorithms (row-wise and iterative), we run the reconstruction algorithm as a function of the number of measurements $M$ taken into account and determine the minimum number of measurements $M_{\rm min}$ for a faithful reconstruction, i.e. zero reconstruction errors. Each numerical experiment is repeated 10 times.

First, we observe that the iterative algorithm allows for a faithful reconstruction from underdetermined input data for all grids under consideration and both data types. The minimum number of measurements $M_{\rm min}$ is much smaller than the grid size $N$ in all cases. The row-wise algorithm performs less well as it does not exploit the symmetry of the matrix $\boldsymbol{B}$.

Second, efficient reconstruction is possible with both algorithms for the type I input data. The theory of compressed sensing shows that a reconstruction of sparse vectors is possible with high probability from a number of measurements scaling only logarithmically with the problem dimension $N$ if some conditions are satisfied [18, 19]. This logarithmic scaling is well confirmed by the results of our numerical experiments. For the iterative algorithm the scaling seems to be even slower.

However, the proofs for the efficiency of compressed sensing depend on two features of the sampling vectors: isotropy and incoherence. Roughly speaking, isotropy means that all directions in the vector space are sampled equally well. Incoherence guarantees that the redundancy in the information gained by two measurements is small. Both features are no longer guaranteed for type II input data. In fact, we observe that the row-wise reconstruction algorithm does not perform well in this case. The minimum number of measurements $M_{\rm min}$ is only slightly smaller than the problem dimension $N$ such that the iterative algorithm performs much better.

VI-C Reconstruction from flow measurements

The iterative algorithm allows to take into account prior knowledge to efficiently reconstruct the missing information about the grid topology. Such a problem can arise in practice when a larger part of the grid is subject to damages or attacks. An example of such a situation is depicted in figure 6 (a). It is assumed that we have no reliable information about the grid in two areas indicated by the question marks. However, we do have information about the connections of all remaining nodes and we assume that we can measure the real power flow along the colored solid lines. The prior knowledge of the grid topology is encoded in the matrix $\boldsymbol{K}$, which is illustrated in figure 6 (b). White entries indicate where we do not know the entries of the matrix $\boldsymbol{B}$.

The missing information is found by combining classical state estimation and network reconstruction. First, the voltage phase angles for all nodes of the network are reconstructed as described in section III, the results being shown in figure 6 (c). Second, the iterative reconstruction algorithm shown in figure 3 is applied using the measured power injections, the prior knowledge $\boldsymbol{B}$ and $\boldsymbol{K}$ as well as the estimated voltage phase angles as input data. The algorithm then faithfully reconstructs all the missing information about the grid topology. The reconstructed values of the line susceptances $X_{sr}^{-1}$ shown in figure 6 (d) exactly match the true values. Most importantly, the full reconstruction is possible already for $M=1$ in this case, i.e. with a single snapshot of the power injections $\boldsymbol{P}$, making use of the prior knowledge.

VII Conclusion and Outlook

We have introduced two algorithms to reconstruct the structure of a power grid from nodal measurements only. One algorithm is straightforward and row-based, reconstructing the lines of the grid node by node. The second is iterative and takes into account both prior knowledge about the presence of absence of lines as well as the knowlegde generated by that algorithm during previous steps. We have demonstrated how these algorithms can be used to reconstruct the entire network structure from time series or missing information on the grid topology from a single snapshot.

The presented algorithms exploit several structural properties of power grids to reduce the number of necessary measurements. Power grids are typically very sparse, i.e. each substation is connected to only few other substations. Hence, methods from compressed sensing can be used which allow for a faithful reconstruction also in the underdetermined case. The iterative reconstruction algorithm makes use of the symmetry of the nodal susceptance matrix. Loosely speaking the algorithm solves the simple parts of the reconstruction problem first and than uses the gained information for the remaining parts. We have shown that this trick leads to a vast reduction of the measurement resources.

The present work presents a step towards hybrid state estimation/network reconstruction algorithms. In future smart grids one can face the situation where both nodal variables and the grid topology are subject to uncertainties or lacking information. A pure state estimation starting from perfect knowledge of the network structure is then no longer sufficient. Advanced grid monitoring algorithms must be able of coping with missing information in both state and structure.

Acknowledgements

We thank D. Gross and R. Kueng for inspiring discussions. We gratefully acknowledge support from the Helmholtz Association (via the joint initiative “Energy System 2050 - A Contribution of the Research Field Energy” and the grant no. VH-NG-1025 to D.W.) and the Federal Ministry of Education and Research (BMBF grant nos. 03SF0472B and  03SF0472E to D.W. and M.T.).