# Bayesian Error-in-Variables Models for the Identification of Power Networks

Jean-Sébastien Brouillon
Institute of Mechanical Engineering,
École Polytechnique Fédérale de Lausanne,
Lausanne, Switzerland
<EMAIL_ADDRESS>
& Emanuele Fabbiani
Identification and Control of Dynamic Systems Laboratory,
University of Pavia,
Pavia, Italy
<EMAIL_ADDRESS>
& Pulkit Nahata
Institute of Mechanical Engineering,
École Polytechnique Fédérale de Lausanne,
Lausanne, Switzerland
<EMAIL_ADDRESS>
& Florian Dörfler
Automatic Control Laboratory,
Swiss Federal Institute of Technology (ETH),
Zurich, Switzerland
<EMAIL_ADDRESS>
& Giancarlo Ferrari-Trecate
Institute of Mechanical Engineering,
École Polytechnique Fédérale de Lausanne,
Lausanne, Switzerland
<EMAIL_ADDRESS>

(5 May 2021)

###### Abstract

The increasing integration of intermittent renewable generation, especially at
the distribution level, necessitates advanced planning and optimisation
methodologies contingent on the knowledge of the grid, specifically the
admittance matrix capturing the topology and line parameters of an electric
network. However, a reliable estimate of the admittance matrix may either be
missing or quickly become obsolete for temporally varying grids. In this work,
we propose a data-driven identification method utilising voltage and current
measurements collected from micro-PMUs. More precisely, we first present a
maximum likelihood approach and then move towards a Bayesian framework,
leveraging the principles of maximum a posteriori estimation. In contrast with
most existing contributions, our approach not only factors in measurement
noise on both voltage and current data, but is also capable of exploiting
available a priori information such as sparsity patterns and known line
parameters. Simulations conducted on benchmark cases demonstrate that,
compared to other algorithms, our method can achieve significantly greater
accuracy.

_Keywords_ Bayesian inference $\cdot$ Distribution grids $\cdot$ Error-in-
variables $\cdot$ Line admittance estimation $\cdot$ Power systems
identification

##  1 Introduction

A major key to realising green energy systems is the large-scale integration
of renewable energy sources (RESs) in the distribution grid. Nevertheless, RES
proliferation leads to additional risks such as reverse power flows and over-
voltage—especially during periods of peak generation and low consumption [1].
Distribution grid operators are consequently required to put in place
intelligent monitoring and control algorithms in order to maintain the
existing levels of grid safety and reliability [2, 3, 4, 5]. Deploying such
algorithms efficiently requires the knowledge of the topology and the line
parameters of the grid, embedded in its admittance matrix.

An exact estimate of the admittance matrix is hard to obtain for distribution
grids, in particular as topological information and line parameter values
either are unavailable for large chunks of the network, or become obsolete in
the event of a topology change. To circumvent this issue, many recent
contributions work out an up-to-date admittance matrix estimate by utilising
data collected from micro-phasor measurement units ($\mu$PMUs). Although a
more recent development than PMUs, which are commonly deployed on transmission
systems, $\mu$PMUs have already been installed in distribution grids across
America, Asia, and Europe, and their penetration is expected to steadily
increase in the coming years [6].

Due to its increasing relevance, the problem of identifying the topology and
line parameters of a power grid attracted attention in the last few years. In
[7, 8], an approach based on inverter probing is explored. Both works, besides
employing approximate linearized power-flow equations, are restricted to
radial networks. Albeit requiring voltage and current (or power) measurements
at each bus of the grid, identification methods in [9, 10, 11, 12, 13] can be
applied to both radial and meshed structures. In [12, 13], structural
properties of the admittance matrix, such as symmetry and Laplacianity are
used to eliminate redundant admittance matrix parameters. Moreover, [13]
proposes an adaptive Lasso algorithm promoting sparsity.

The tradeoff between voltage stability and the energy of the signal is the
main challenge for the identification of power systems. To compensate for the
lack of signal, online design-of-experiment procedures are presented in [9,
10]; nonetheless, the proposed algorithms require control authority on the
state of the grid, and additional measurements of line power flows. They also
neglect the structural and sparsity properties of the admittance matrix and
stay limited by the small acceptable voltage variations.

All the foregoing works suffer from two limitations. First, they either
completely disregard measurement errors or assume errors solely on certain
measurements. This creates an estimation bias, for $\mu$PMUs introduce an
unavoidable error on all measured electric variables [14, 15]. Second, they do
not capitalize on grid information which may already be available a priori,
for instance sparsity patterns and known network sections and line parameters.
To do away with the first limitation, [16, 17] introduce error-in-variable
(EIV) models taking into consideration all sources of measurement errors. That
notwithstanding, they leave aside all prior information, including structural
properties of the admittance matrix, which can potentially improve grid
identification.

In this paper, we address the limitations of existing works by putting forth a
novel Bayesian grid identification framework, which incorporates EIV models
with an unbiased estimation of the error on both voltage and current data, and
takes advantage of the principles of maximum likelihood estimation (MLE). Our
approach exploits not only the inherent structural properties of the
admittance matrix, but enables to incorporation of grid information known a
priori. We present and describe how to formalize different pieces of
information which may be available to grid operators, partly through data-
driven Bayesian priors. This work has high tutorial value, as it provides a
framework into which new methods can be fitted, and presents not just a formal
description of the problem and the proposed methods, but also practical
solutions to several common issues, from data preprocessing to the numerical
resolution of the optimization problems arising from EIV estimates.

In order to substantiate the efficacy of our method, we conduct simulations on
a large network with realistic voltage and current profiles, and $\mu$PMU
noise levels compatible with the accuracy of actual commercial devices. We
then compare the performance of our proposed methods with other grid
identification procedures proposed in the literature. Our analysis shows not
only that EIV models are needed to obtain reasonable grid estimates, but also
that sparsity needs to be enforced if the topology of the network is unknown.
Moreover, it substantiates the value of injecting prior information into the
estimation algorithm, in case it is available.

The paper proceeds as follows: Sections 2, 3 and 4 define a formal model for
the identification problem, which we solve using likelihood maximization in
Section 5. Sections 6 and 7 introduce prior knowledge using the Bayesian
framework. Section 8 provides numerical methods for the _a posteriori_
optimization and Section 9 presents a realistic simulation and discusses its
outcome. Section 10 concludes and proposed future developments.

###  1.1 Preliminaries and Notation

Let $j=\sqrt{-1}$ denote the imaginary unit. For $x\in\mymathbb{C}^{n}$,
$\overline{x}$ is its complex conjugate and $|x|$ its magnitude, both taken
element-wise. The diagonal matrix of order $n$ associated with $x$ is written
$[x]$. The $\ell_{1}$ and $\ell_{2}$ norms of a vector $x$ are represented by
$\left\lVert x\right\rVert_{1}$ and $\left\lVert x\right\rVert_{2}$,
respectively. Throughout, $\mymathbb{1}_{n}$ and $\mymathbb{0}_{n}$ are
$n$-dimensional vectors of all ones and zeros, whereas $\mathcal{I}_{n}$ and
$\mathcal{O}_{n\times m}$ represent $n$-by-$n$ identity and $m$-by-$n$ zero
matrices, respectively. The unit vector $e_{i},~{}i=1,...,n$ is the $i^{th}$
column of $\mathcal{I}_{n}$. For an $(m,n)$ matrix $A$, $A^{\top}$ denotes its
transpose, $A_{i\cdot}$ its $i^{th}$ row vector, and
$\operatorname*{\textrm{vec}}(A)=[A_{1\cdot}^{\top}\cdots
A_{n\cdot}^{\top}]^{\top}$ the $mn$-dimensional stacked column vector. Given a
square matrix $A$, $\operatorname*{\textrm{ve}}(A)$ is the
$n(n-1)/2$-dimensional vector obtained by removing diagonal and supra-diagonal
elements from $-\operatorname*{\textrm{vec}}(A)$. The Kronecker product
between matrices $A$ and $B$ is $A\otimes B$. Given $n$ elements $x_{n}$,
$[x_{i}]_{i=1}^{n}$ is the stacked column vector $[x_{1}^{\top}\cdots
x_{n}^{\top}]^{\top}$.

Random variables. $\mathcal{X}\sim\mathcal{N}(\mu,\sigma^{2})$ denotes a
Gaussian random variable with expected value
$\operatorname*{E}[\mathcal{X}]=\mu$ and variance
$\operatorname*{Var}[\mathcal{X}]=\operatorname*{Cov}[\mathcal{X},\mathcal{X}]=\sigma^{2}$.

Algebraic graph theory. We denote by
$\mathcal{G}(\mathcal{V},\mathcal{E},\mathcal{W})$ an undirected, weighted,
and connected graph, where $\mathcal{V}$ is the node set of cardinality $n$,
$\mathcal{E}\subseteq(\mathcal{V}\times\mathcal{V})$ the edge set, and
$\mathcal{W}$ collects the edge weights. The adjacency matrix
$W\in\mymathbb{C}^{n\times n}$ has elements $w_{hk}$ corresponding to the
weight of the edges $(h,k)$, if $(h,k)\in\mathcal{E}$ and zero otherwise. The
matrix $L=[W\,\mymathbb{1}_{n}]-W$, $L\in\mymathbb{C}^{n\times n}$, is the
Laplacian matrix associated with $\mathcal{G}$. By definition, a Laplacian
matrix is symmetric and such that $L\mymathbb{1}_{n}=\mymathbb{0}_{n}$.

##  2 Grid model and data collection

###  2.1 Power grid model

An electric distribution network is modeled as an undirected, weighted, and
connected graph $\mathcal{G}(\mathcal{V},\mathcal{E},\mathcal{W})$, where the
nodes in $\mathcal{V}=\\{1,2,\dots,n\\}$ represent buses, either generators or
loads, and edges represent power lines, each connecting two distinct buses and
modeled after the standard lumped $\pi-$model [18]. To each edge
$(h,k)\in\mathcal{E}$ we associate a complex weight equal to the line
admittance $y_{hk}=g_{hk}+jb_{hk}$, where $g_{hk}>0$ is the line conductance
and $b_{hk}\in\mathbb{R}$ the line susceptance.

The network is then completely represented by the $n$-by-$n$ complex
admittance matrix $Y$, with elements $Y_{hk}=-y_{hk}$ for $h\neq k$ and
$Y_{hk}=\sum_{h=1,h\neq k}^{n}y_{hk}+y_{s,h}$, where $y_{s,h}$ is the shunt
element at the $h^{th}$ bus. The admittance matrix $Y$ is typically sparse, as
each bus is connected to few others: notably, this is the case in distribution
grids, which are often characterized by a radial topology. Moreover, for
network where shunt elements are negligible, $Y$ is Laplacian [19].

We consider the network to be either single-phased or phase-balanced, and
operating in sinusoidal regime. To each bus $h\in\mathcal{V}$, we associate a
phasor voltage $v_{h}e^{j\theta_{h}}\in\mathbb{C}$, where $v_{h}>0$ is the
voltage magnitude and $\theta_{h}\in\mathbb{R}$ the voltage angle, and a
phasor current $i_{h}e^{j\phi_{h}}\in\mathbb{C}$, representing the injection
at the bus. We do not assume the presence of a point of common coupling (PCC),
although one or more may be present as long as their fixed $v_{0}$ and
$\theta_{0}$ are known. The current-voltage relation descending directly from
Kirchhoff’s and Ohm’s laws is given by

$i=Yv,$ (1)

where $i\in\mathbb{C}^{n}$ is the vector of nodal current injections, and
$v\in\mathbb{C}^{n}$ the vector of nodal voltages [20].

###  2.2 Data collection

We assume that each bus of the network is equipped with $\mu$PMU, while we do
not require electrical variables to be measured on the lines.

######  Assumption 1.

The network is completely observable, that is current injections and voltages
are measured at each node.

In transmission systems, where the reliability and the economic optimization
of dispatch are primary concerns, synchronized phasor measurements are
provided by PMUs, which sample the magnitude and phase of current and voltage
phasors. Unfortunately, distribution networks are characterized by relatively
small flows of active power and predominantly resistive lines, resulting in
small phase differences between nodes. The accuracy of standard PMUs is not
high enough to reliably appreciate such differences, making the devices
ineffective.

Micro-syncrophasors improve the resolution and accuracy of PMUs by up to two
orders of magnitude, while preserving low costs - see Table 1. As PMUs, these
devices sample the magnitude and phase of current and voltage: state-of-the-
art models achieve frequency up to 120 Hz [21]. While the fast sampling rate
generates a large number of data points in a short time span, the samples are
highly correlated. Indeed, the characteristic frequencies of load variations
are much lower than the sampling rate of $\mu$PMUs. Moreover, due to structure
of the network and the low phase difference between nodes, samples collected
on different buses are highly correlated with each other [13].

Metric | P-class PMU |  $\mu$PMU
---|---|---
Magnitude accuracy [%] | 0.1% | 0.01%
Phase accuracy [rad] | $1.5\cdot 10^{-3}$ | $10^{-4}$
Table 1: PMU and $\mu$PMU accuracy [15, 21]. Magnitude accuracy is a reported
as a percentage of the rated value.

##  3 Problem Statement

Consider a power distribution network as described in Section 2, fulfilling
Assumption 1. The identification problem amounts to reconstruct the admittance
matrix from a sequence of voltage and current measurements corresponding to
different steady states of the system [11, 13].

Let $N$ be the number of samples, and let $v_{t}\in\mymathbb{C}^{n}$ and
$i_{t}\in\mymathbb{C}^{n}$ be the vectors of current injections and voltages
for $t=1,\dots,N$. From (1), one can obtain

$I=YV,$ (2)

where $V=[v_{1},v_{2},\dots,v_{N}]\in\mymathbb{C}^{N\times n}$, and
$I=[i_{1},i_{2},\dots,i_{N}]\in\mymathbb{C}^{N\times n}$.

As described in Section 2.2, the available current and voltage phasors are
corrupted by measurement noise. Therefore, in place of the actual electrical
variables $V$ and $I$, only noisy samples $\tilde{V}$ and $\tilde{I}$ are
available, where

$\displaystyle\tilde{V}$ $\displaystyle=V+\Delta V,$ (3a)
$\displaystyle\tilde{I}$ $\displaystyle=I+\Delta I,$ (3b)

with $\Delta V\in\mymathbb{C}^{N\times n}$ and $\Delta
I\in\mymathbb{C}^{N\times n}$ denoting the complex measurement noise. The
network identification problem then translates into the estimation of $Y$
given $\tilde{V}$ and $\tilde{I}$.

##  4 Noise Model

By design, PMUs and $\mu$PMUs collect current and voltage measurements in
polar coordinates, that is in terms of magnitude and phase [15, Sec. 3.2].
Previous studies have shown, with both theoretical and empirical arguments,
that the measurement noise is approximately Gaussian in polar coordinates,
with zero mean and constant variance [15, Sec. 2.1].

In (1), the admittance matrix $Y$ establishes a linear relationship between
the real and the imaginary parts of $i$ and $v$, but the equation becomes non-
linear if magnitude and phase are considered. To preserve linearity, one needs
to transform the measurements and their associated noise from polar to
Cartesian coordinates. This transformation changes the statistical
distribution of the noise, which becomes non-Gaussian, with a non-zero mean,
and with a time-varying, non-diagonal covariance matrix.

We consider a generic phasor measured by a $\mu$PMU. Without loss of
generality, we will discuss only the case of a voltage phasor; the same
arguments apply to the current. Let $\tilde{v}$ and $\tilde{\theta}$ denote
the measured magnitude and phase, and let $v$ and $\theta$ be the actual
unobservable variables. Then $\tilde{v}=v+\epsilon$ and
$\tilde{\theta}=\theta+\delta$, where
$\epsilon\sim\mathcal{N}(0,\sigma_{\epsilon})$ and
$\delta\sim\mathcal{N}(0,\sigma_{\delta})$ are zero-mean Gaussian variables.
Previous studies suggest that it is safe to introduce the following Assumption
[15, Sec. 2.1].

######  Assumption 2.

The samples taken at two different time steps on the same node, and the
samples taken on two different nodes at the same time are independent.

Using Assumption 2, we can focus on a single sample. Our aim is to write the
measured phasor $\tilde{v}e^{j\tilde{\theta}}$ in Cartesian coordinate
($\tilde{c}+j\tilde{d}$). Highlighting the noise component $\Delta c+j\Delta
d$ gives $\tilde{v}e^{j\tilde{\theta}}=(c+\Delta c)+j(d+\Delta d)$, where

$\displaystyle\Delta c$
$\displaystyle=\tilde{c}-c=v\cos\theta(\cos\delta-1)-v\sin\theta\sin\delta+\epsilon\cos\theta\cos\delta-\epsilon\sin\theta\sin\delta,$
(4a) $\displaystyle\Delta d$
$\displaystyle=\tilde{d}-d=v\sin\theta(\cos\delta-1)+v\cos\theta\sin\delta+\epsilon\sin\theta\cos\delta+\epsilon\cos\theta\sin\delta.$
(4b)

From (4), it can be noted that $\Delta c$ and $\Delta d$ are not distributed
as Gaussian variables, due to interaction terms like $\epsilon\sin\delta$.
However, the first-order Taylor approximation of (4) is a linear combination
of Gaussian variables, suggesting that, for low noise levels, the distribution
of $\Delta c$ and $\Delta d$ is closely approximated by a Gaussian variable.
Similarly to [17], to characterize the noise we compute the expected value of
$\Delta c$ and $\Delta d$

$\displaystyle\operatorname*{E}[\Delta c]$
$\displaystyle=v\cos\theta(e^{-\sigma^{2}_{\delta}/2}-1),$ (5a)
$\displaystyle\operatorname*{E}[\Delta d]$
$\displaystyle=v\sin\theta(e^{-\sigma^{2}_{\delta}/2}-1),$ (5b)

and the associated variance-covariance terms:

$\displaystyle\operatorname*{Var}[\Delta c]$
$\displaystyle=v^{2}e^{-\sigma_{\delta}^{2}}[\cos^{2}\theta(\cosh\sigma_{\delta}^{2}-1)+\sin^{2}\theta\sinh{\sigma_{\delta}^{2}}]+\sigma_{\epsilon}^{2}e^{-\sigma_{\delta}^{2}}[\cos^{2}\theta\cosh\sigma_{\delta}^{2}+\sin^{2}\theta\sinh{\sigma_{\delta}^{2}}],$
(6a) $\displaystyle\operatorname*{Var}[\Delta d]$
$\displaystyle=v^{2}e^{-\sigma_{\delta}^{2}}[\sin^{2}\theta(\cosh\sigma_{\delta}^{2}-1)+\cos^{2}\theta\sinh{\sigma_{\delta}^{2}}]+\sigma_{\epsilon}^{2}e^{-\sigma_{\delta}^{2}}[\sin^{2}\theta\cosh\sigma_{\delta}^{2}+\cos^{2}\theta\sinh{\sigma_{\delta}^{2}}],$
(6b) $\displaystyle\operatorname*{Cov}[\Delta c,\Delta d]$
$\displaystyle=\sin\theta\cos\theta
e^{-2\sigma_{\delta}^{2}}[\sigma_{\epsilon}^{2}+v^{2}(1-e^{\sigma_{\delta}^{2}})].$
(6c)

Unfortunately, the expressions in Eqs. 5 and 6 are of no practical use, as
they rely on the actual unobservable values $v$ and $\theta$. Similar issues
arise in the context of state estimation based on Extended Kalman Filter
(EKF): in particular, research on the filtering of radar signal can be adapted
to our case [22, 23, 24, 25, 15]. Following such developments, we compute the
expectation of the average (5) and the variance (6) conditioned on the
measurements:

$\displaystyle\operatorname*{E}[\Delta c|\tilde{v},\tilde{\theta}]$
$\displaystyle=\tilde{v}\cos\tilde{\theta}(e^{-\sigma^{2}_{\delta}}-e^{-\sigma^{2}_{\delta}/2}),$
(7a) $\displaystyle\operatorname*{E}[\Delta d|\tilde{v},\tilde{\theta}]$
$\displaystyle=\tilde{v}\sin\tilde{\theta}(e^{-\sigma^{2}_{\delta}}-e^{-\sigma^{2}_{\delta}/2}).$
(7b)

The same procedure can be applied to the variances:

$\displaystyle\textrm{Var}[\Delta c|\tilde{v},\tilde{\theta}]$
$\displaystyle=\tilde{v}^{2}e^{-2\sigma_{\delta}^{2}}[\cos^{2}\tilde{\theta}(\cosh
2\sigma_{\delta}^{2}-\cosh\sigma_{\delta}^{2})+\sin^{2}\tilde{\theta}(\sinh
2\sigma_{\delta}^{2}-\sinh\sigma_{\delta}^{2})]+$
$\displaystyle+\sigma_{\epsilon}^{2}e^{-2\sigma_{\delta}^{2}}[\cos^{2}\tilde{\theta}(2\cosh
2\sigma_{\delta}^{2}-\cosh\sigma_{\delta}^{2})+\sin^{2}\tilde{\theta}(2\sinh
2\sigma_{\delta}^{2}-\sinh\sigma_{\delta}^{2})],$ (8a)
$\displaystyle\textrm{Var}[\Delta d|\tilde{v},\tilde{\theta}]$
$\displaystyle=\tilde{v}^{2}e^{-2\sigma_{\delta}^{2}}[\sin^{2}\tilde{\theta}(\cosh
2\sigma_{\delta}^{2}-\cosh\sigma_{\delta}^{2})+\cos^{2}\tilde{\theta}(\sinh
2\sigma_{\delta}^{2}-\sinh\sigma_{\delta}^{2})]+$
$\displaystyle+\sigma_{\epsilon}^{2}e^{-2\sigma_{\delta}^{2}}[\sin^{2}\tilde{\theta}(2\cosh
2\sigma_{\delta}^{2}-\cosh\sigma_{\delta}^{2})+\cos^{2}\tilde{\theta}(2\sinh
2\sigma_{\delta}^{2}-\sinh\sigma_{\delta}^{2})],$ (8b)
$\displaystyle\textrm{Cov}[\Delta c,\Delta d|\tilde{v},\tilde{\theta}]$
$\displaystyle=\sin\tilde{\theta}\cos\tilde{\theta}e^{-4\sigma_{\delta}^{2}}[\sigma_{\epsilon}^{2}+(\tilde{v}^{2}+\sigma_{\epsilon}^{2})(1-e^{\sigma_{\delta}^{2}})].$
(8c)

Eq. 7 suggests that the measurement in Cartesian coordinates are biased, as
the noise has a non-zero average. However, the arguments in Appendix D suggest
that such bias is negligible for realistic noise levels. Moreover, it can
always be computed and substracted from the data. Hence, in the following, the
noise will be considered unbiased.

We finally model the noise on a phasor measurement as

$\begin{bmatrix}\Delta c\\\ \Delta
d\end{bmatrix}\sim\mathcal{N}\left(\mymathbb{0}_{2},\Sigma\right),$ (9)

with the elements of $\Sigma$ defined by (8). The covariance matrix $\Sigma$
is not constant in time, but changes with the actual values of phase and
magnitude: such property will be further discussed in Section 5, while
presenting the estimation methods.

##  5 Frequentist identification

###  5.1 Least squares

From (2) and (3), noisy data are related by the model

$\displaystyle\tilde{I}-\Delta I=(\tilde{V}-\Delta V)Y.$ (10)

For reconstructing $Y$, the works [13, 12] assume $\Delta V=0$ (i.e.
$\tilde{V}=V$) and use Ordinary Least Squares (OLS) as well as its recursive
and regularized variants. For standard OLS, one obtains

$\displaystyle\hat{Y}_{\textrm{OLS}}=\operatorname*{arg\,min}_{\hat{Y}}\|\tilde{I}-\tilde{V}\hat{Y}\|_{F}^{2}.$
(11)

However, if $\Delta V\neq 0$, the OLS introduces a bias. In this case, the
Total Least Squares (TLS) is an unbiased estimator [26]:

$\displaystyle\hat{Y}_{\textrm{TLS}}=\operatorname*{arg\,min}_{\hat{Y}}\min_{I,V}\|[\tilde{V}-V,\tilde{I}-I]\|_{F}^{2}\quad\textrm{s.t.}\;I=V\hat{Y}.$
(12)

Closed-form solutions for both estimators are well-known, and can be written
column-wise, with centered data $\tilde{V}_{c}$ and $\tilde{I_{c}}$, as

$\displaystyle\hat{Y}_{i,\textrm{OLS}}$
$\displaystyle=(\tilde{V}_{c}^{\top}\tilde{V}_{c})^{-1}\tilde{V}_{c}^{\top}\tilde{I}_{i,c},$
(13a) $\displaystyle\hat{Y}_{i,\textrm{TLS}}$
$\displaystyle=(\tilde{V}_{c}^{\top}\tilde{V}_{c}-\sigma_{n+1}^{2}\mathcal{I})^{-1}\tilde{V}_{c}^{\top}\tilde{I}_{i,c},$
(13b)

where $\sigma_{n+1}$ is the smallest singular value of $[V_{c},I_{i,c}]$. The
error covariance, computed in terms of the exact values $V$ [27], is

$\displaystyle\textrm{Var}[\hat{Y}_{i,\textrm{OLS}}]$
$\displaystyle\approx\frac{\sigma_{n+1}^{2}}{N}(V_{c}^{\top}V_{c})^{-1},$
(14a) $\displaystyle\textrm{Var}[\hat{Y}_{i,\textrm{TLS}}]$
$\displaystyle\approx(1+\|\hat{Y}_{i,TLS}\|^{2})\frac{\sigma_{n+1}^{2}}{N}(V_{c}^{\top}V_{c})^{-1}.$
(14b)

It appears that $Y_{i,\textrm{OLS}}$ has a smaller variance than
$Y_{i,\textrm{OLS}}$, but has a bias that grows with $\sigma_{n+1}$. In power
systems, the data covariance $(V_{c}^{\top}V_{c})^{-1}$ is very small compared
to the noise variance (approximated by $\sigma_{n+1}$). Hence, the bias of the
OLS can be very large, making the TLS a more suitable choice. However, the
estimation of a sparse topology can be difficult, because the TLS estimate of
zero elements can be very large due to the large variance. Possible solutions
for this problem are regularization or more complex iterative methods such as
the one presented in [16]. Both methods assume that the samples are
independent and identically distributed (i.i.d). As shown in section 4, µPMU
measurements are not identically distributed.

###  5.2 Maximum likelihood estimator

The high correlation between measurements observed in power systems suggest
that large sample sizes $N$ may be needed in order to obtain a good estimate
of $Y$. Moreover, the estimator should be unbiased and consistent [28, chapter
7, 10]. The Maximum Likelihood Estimator (MLE) presents weak consistency
conditions [29] that are satisfied for linear models. From (2) and (3), noisy
data are related by the model

$\displaystyle\tilde{I}-\Delta I=(\tilde{V}-\Delta V)Y.$ (15)

Considering Gaussian noise $\Delta V$ and $\Delta I$, as described in Section
2.2, the MLE is a weighted Total Least Squares (TLS) estimator and thus shares
most of its properties such as _de_ regularization [26] and unbiasedness [30].

The variables $V$ and $I$ are deterministic, and the noises $\Delta V$ and
$\Delta I$ are independent. The likelihood
$l(\tilde{V},\tilde{I}|V,I,\hat{Y})$ of $(\tilde{V},\tilde{I})$, can be
written as follows.

$\displaystyle l(\tilde{V},\tilde{I}|V,I,\hat{Y})\propto\;$ $\displaystyle
p(V+\Delta V|V,\hat{Y})p(I+\Delta I|I,\hat{Y}),$ (16)
$\displaystyle\textrm{s.t.}\;(\tilde{V}-\Delta V)\hat{Y}=\tilde{I}-\Delta I.$

To work with real variables, we define

$\displaystyle\boldsymbol{v}$
$\displaystyle=\left(\begin{matrix}\Re\left(\operatorname*{\textrm{vec}}(V)\right)\\\
\Im\left(\operatorname*{\textrm{vec}}(V)\right)\end{matrix}\right),\quad\boldsymbol{i}=\left(\begin{matrix}\Re\left(\operatorname*{\textrm{vec}}(I)\right)\\\
\Im\left(\operatorname*{\textrm{vec}}(I)\right)\end{matrix}\right),$ (17)
$\displaystyle\mymathbb{V}$
$\displaystyle=\left(\begin{matrix}\Re(\mathcal{I}_{n}\otimes
V)&-\Im(\mathcal{I}_{n}\otimes V)\\\ \Im(\mathcal{I}_{n}\otimes
V)&\Re(\mathcal{I}_{n}\otimes V)\end{matrix}\right),\quad
y=\left(\begin{matrix}\Re(\operatorname*{\textrm{vec}}(Y))\\\
\Im(\operatorname*{\textrm{vec}}(Y))\end{matrix}\right).$

The same transformations can be applied to $\tilde{V}$, $\tilde{I}$, $\Delta
V$, $\Delta I$, and $\hat{Y}$, resulting in the vectors
$\tilde{\boldsymbol{v}}$, $\tilde{\boldsymbol{i}}$, $\Delta\boldsymbol{v}$,
$\Delta\boldsymbol{i}$, $\hat{y}$, as well as the matrices
$\tilde{\mymathbb}{V}$ and $\Delta\mymathbb{V}$. Note that $\boldsymbol{v}$
and $\mymathbb{V}$ contain the same elements but arranged differently. We will
therefore use them interchangeably when describing optimization problems over
$\Delta\boldsymbol{v}$ or $\Delta\mymathbb{V}$. The matrix $\mymathbb{V}$ is
introduced to represent the product $VY$ with real and vectorized quantities
as shown in (18) below.

Using the vectorized notations (17), we assume that
$\Delta\boldsymbol{v}\sim\mathcal{N}(0,\Sigma_{v})$ and
$\Delta\boldsymbol{i}\sim\mathcal{N}(0,\Sigma_{i})$, as per the approximate
noise model discussed in Section 2.2. The covariance matrices $\Sigma_{v}$ and
$\Sigma_{i}$ are computed from (8) as explained in Appendix A. The likelihood
(16) becomes

$\displaystyle
l(\tilde{\boldsymbol{v}},\tilde{\boldsymbol{i}}|\boldsymbol{v},\boldsymbol{i},\hat{y})=\;$
$\displaystyle
e^{-\Delta\boldsymbol{i}^{\top}\Sigma_{i}^{-1}\Delta\boldsymbol{i}}e^{-\Delta\boldsymbol{v}^{\top}\Sigma_{v}^{-1}\Delta\boldsymbol{v}},$
(18a)
$\displaystyle\textrm{s.t.}\;\tilde{\boldsymbol{i}}-\Delta\boldsymbol{i}=(\tilde{\mymathbb}{V}-\Delta\mymathbb{V})y.$
(18b)

The corresponding log-likelihood is

$\displaystyle\mathcal{L}(\tilde{\boldsymbol{v}},\tilde{\boldsymbol{i}}|\boldsymbol{v},\boldsymbol{i},\hat{y})=$
$\displaystyle-\Delta\boldsymbol{i}^{\top}\Sigma_{i}^{-1}\Delta\boldsymbol{i}-\Delta\boldsymbol{v}^{\top}\Sigma_{v}^{-1}\Delta\boldsymbol{v},$
(19)

subject to (18b). For a fixed, albeit unknown $\boldsymbol{v}$ and
$\boldsymbol{i}$, we use the shorthand notation
$\mathcal{L}(\hat{y},\Delta\boldsymbol{v},\Delta\boldsymbol{i})$. Minimizing
$-\mathcal{L}$ for $\Delta\boldsymbol{i}$, $\Delta\boldsymbol{v}$ and
$\hat{y}$ yields the MLE.

###  5.3 Error covariance analysis

Let $\bar{V}=\frac{1}{N}\mymathbb{1}_{N}^{\top}V$ be the vector of voltage
means, and $\Phi=\mathcal{I}_{n}\otimes(V-\mymathbb{1}_{N}\bar{V})$ be the
centered complex data matrix. The MLE with Gaussian noise as defined in
Section 4 being a weighted TLS estimator, expressions for its error covariance
can be found in the literature [27, 30]. For a time-variant, not identically
distributed case such as (19), [30] proves that the error covariance
$\Sigma_{\textrm{MLE}}$ is equal to the inverse of the Fischer information
matrix $F_{\textrm{MLE}}$ written in [30, (48)]. We adapt this expression for
the complex matrix $\Phi$ using the same transformation to real numbers as for
$V$ in (17), and obtain

$\displaystyle
F_{\textrm{MLE}}=\sum_{q,h,t=1}^{2,n,N}\left(\begin{matrix}\frac{\Re(\Phi_{hN+t})^{\top}\Re(\Phi_{hN+t})}{\Re(z)^{\top}\mathcal{R}_{\Re,qht}\Re(z)}&\frac{\Re(\Phi_{hN+t})^{\top}\Im(\Phi_{hN+t})}{\Re(z)^{\top}\mathcal{R}_{\Re\Im,qht}\Im(z)}\\\
\star&\frac{\Im(\Phi_{hN+t})^{\top}\Im(\Phi_{hN+t})}{\Im(z)^{\top}\mathcal{R}_{\Im,qht}\Im(z)}\end{matrix}\right),$
(20)

where $z=[\operatorname*{\textrm{vec}}(Y)^{\top},-1-j]^{\top}$ and
$\mathcal{R}_{\Re,qht}$, $\mathcal{R}_{\Re\Im,qht}$, and
$\mathcal{R}_{\Im,qht}$ are diagonal and constructed from $\Sigma_{v}$ and
$\Sigma_{i}$ as in [30]. The exact expression is presented in Appendix B.

One should note that $\hat{F}_{\textrm{MLE}}$, computed from noisy data
instead of the exact variables used for $F_{\textrm{MLE}}$, is a good
approximation only if the signal to noise ratio is high enough. Although this
is typically not the case in distribution grids, $F_{\textrm{MLE}}$ can still
be used for theoretical purposes, such as design of experiments [12] or
avoiding unobservable problems (Section 8.1). $F_{\textrm{MLE}}$ also shows
the following properties of the MLE.

######  Lemma 3.

The columns of $\hat{Y}_{\textrm{MLE}}$ are independent. The variance of each
column $(\hat{Y}_{\textrm{MLE}}^{\top})_{h}$ depends only on the same column
$Y_{\cdot h}$ of the exact admittance matrix $Y$.

######  Corollary 4.

With $\textrm{Cov}[\Re(\tilde{V}_{t}),\Im(\tilde{V}_{t})]\approx 0$ for all
$t$ (e.g. with linearized power flow equations), the variances of both real
and imaginary parts of $(\hat{Y}_{\textrm{MLE}}^{\top})_{h}$ are monotone with
their respective values, and their covariance is constant.

###### Proof.

The proofs of both Lemma 3 and Corollary 4 can be found in Appendix B. ∎

######  Remark.

It follows from (8c) that the assumption of Corollary 4 is also often true in
in distribution networks, where voltage angles are very small.

Lemma 3 and Corollary 4 are used in Section 5 and Section 9 to identify and
explain what factors can make the parameter identification very imprecise or
even impossible.

##  6 Bayesian estimation

Line admittances, even if measured, are typically known up to a tolerance.
Some knowledge of $Y$’s structure, such as its sparsity, may also not be
certain or precisely defined. This kind of uncertainty can be modeled via
Bayesian prior distributions.

###  6.1 Maximum _a posteriori_ estimation

Following [31], we describe how to compute Maximum _A Posteriori_ (MAP)
estimates for the error-in-variables model (18). Using Bayes’ rule, the
posterior probability density is

$\displaystyle
p(V,I,\hat{Y}|\tilde{V},\tilde{I})=p(\tilde{V},\tilde{I}|V,I,\hat{Y})\frac{p(V,I)}{p(\tilde{V},\tilde{I})}p(\hat{Y}),$
(21)

where we assume that the line parameters $y$ are independent of the grid state
$(V,I)$ or its measurement $(\tilde{V},\tilde{I})$. The factor
$\frac{p(V,I)}{p(\tilde{V},\tilde{I})}$ can be neglected as it is a quotient
of non-informative priors [31], defined as uniform distributions on the finite
set of feasible voltages and currents. The negative log-posterior minimization
of (21) is then written as

$\displaystyle\min_{\hat{y},\Delta\boldsymbol{v},\Delta\boldsymbol{i}}$
$\displaystyle-\mathcal{L}(\hat{y},\Delta\boldsymbol{v},\Delta\boldsymbol{i})-\log\left(p(\hat{y})\right),$
(22a)
$\displaystyle\tilde{\boldsymbol{i}}-\Delta\boldsymbol{i}=\left(\tilde{\mymathbb{V}}-\Delta\mymathbb{V}\right)\hat{y},$
(22b)

with $\mathcal{L}$ defined by (19). Optimizing (22) provides a maximum _a
posteriori_ (MAP) estimate $\hat{y}_{\textrm{MAP}}$111The optimizer may not be
unique, in this case $\hat{y}_{\textrm{MAP}}$ is one of the elements of the
set of optimizers..

The density function $p(\hat{y})$ can take many forms. If it is Gaussian, then
$-\log(p(\hat{y}))$ corresponds to a weighted ridge regularization [32]. In
this paper, we will focus mainly on the element-wise Laplace distribution
$p(\hat{y}_{h})\propto e^{-\lambda|\hat{y}_{h}|}$, where all elements
$p(\hat{y}_{h})$ of the prior are assumed independent, and therefore
$p(\hat{y})=\prod_{h}p(\hat{y}_{h})$ [33]. One obtains
$-\log(p(\hat{y}))=\lambda\|\hat{y}\|_{1}+\textrm{const}$, where the constant
can be neglected in the optimization problem, obtaining an $\ell_{1}$
regularization term. The $\ell_{1}$ regularization can also be interpreted as
a robustification of the MLE optimization problem
$\min_{\hat{y},\Delta\boldsymbol{v},\Delta\boldsymbol{i}}\;-\mathcal{L}(\hat{y},\Delta\boldsymbol{v},\Delta\boldsymbol{i})$
[34].

Prior distributions are centered on the believed value of the exact admittance
$y_{i}$, which can be different from zero, e.g., in case of an existing line.
More generally, one can also believe that a linear combination of $y$ has a
particular value [35], for example that two lines have the same admittance.
The probability density $p(L\hat{y}-\mu)$ of a linear transformation
$\hat{y}\rightarrow L\hat{y}-\mu$ can describe such a belief. The penalty
function is then

$\displaystyle-\log(p(\hat{y}))=\lambda\|L\hat{y}-\mu\|_{1}.$ (23)

######  Remark.

On one hand, if one has two priors about a line $h$, conditioned on
independent events B and C, then similarly to (21), Bayes’ rule gives that
$p(\hat{y}_{h})$ is proportional to $p(\hat{y}_{h}|B)p(\hat{y}_{h}|C)$. If
$p(\hat{y}_{h}|B)$ and $p(\hat{y}_{h}|C)$ are both Laplace distributions, this
is equivalent to adding a row to $L$ and $\mu$ in (23). The same operation can
be repeated for a larger number of priors on $h$.

On the other hand, if one has no prior about the line $h$, then the non-
informative prior $p(\hat{y}_{h})$ is a uniform distribution on the bounded
support of $\hat{y}_{h}$, and can be factorized away in $p(\hat{y})$, leaving
$p(\hat{y})\propto\prod_{k\neq h}p(\hat{y}_{k})$.

Introducing additional information may reduce the variance of the estimate.
Using the approximation
$\|\hat{y}\|_{1}\approx\hat{y}^{\top}[|y|+\alpha]^{-1}\hat{y}$ similarly to
[36, Appendix A], we can make $\log(p(\hat{y}))$ in (22) smooth. Hence, the
Fisher information is given by the Hessian of the log-likelihood with smoothed
prior.

$\displaystyle F_{\textrm{MAP}}\approx
E\left[\frac{\partial^{2}}{\partial\hat{y}^{2}}(-\mathcal{L}(\hat{y},\Delta
a,\Delta b)-\log(p(\hat{y})))\right].$ (24)

Both the expected value and second derivative operators are distributive. For
$p(\hat{y})\propto e^{-\lambda\|L\hat{y}-\mu\|_{1}}$, $F_{\textrm{MAP}}$ is
therefore approximated as

$\displaystyle F_{\textrm{MAP}}\approx F_{\textrm{MLE}}+\lambda
L^{\top}([|Ly-\mu|]+\alpha\mathcal{I}])^{-1}L.$ (25)

If the data matrix $\mymathbb{V}$ is not full rank so neither is
$F_{\textrm{MLE}}$, (25) can show where prior knowledge is required to make
$F_{\textrm{MAP}}$ full rank.

###  6.2 Prior distribution from known grid parameters

In this section, we will explain how to introduce prior knowledge to the
estimation problem. This knowledge can be directly parameter values, or if
they are all unknown, the Bayesian prior can still include other known grid
properties such as sparsity, constraints on the signs of parameters, or on
ratios between their values.

Introducing the value $\beta_{h}$ known _a priori_ for a single parameter
$y_{h}$ can be done as follows.

$\displaystyle-\log(p(\hat{y}|\beta_{h}))=\lambda\|e_{h}^{\top}(\hat{y}-\beta_{h})\|_{1}.$
(26)

If a set $\mathcal{H}$ of parameters is known, and $\beta$ contains their
values, the log prior distribution becomes

$\displaystyle-\log(p(\hat{y}|\beta))=\lambda\|[L_{h}e_{h}^{\top}]_{h\in\mathcal{H}}(\hat{y}-\beta)\|_{1},$
(27)

where each $L_{h}$ is a diagonal matrix that represents the confidence in the
corresponding value $\beta_{h}$.

Sparsity is the belief that each parameter has a high probability to be zero.
In the Bayesian framework, this translates into a zero-centered distribution,
identical for all parameters, which is a standard $\ell_{1}$ penalty
$\lambda\|\hat{y}\|_{1}$ in the log-space. If some lines are known to exist or
not, the sparsity prior can be removed or strengthened only for them. Let
$\beta_{s}\in\\{0,1,K\\}^{2n^{2}}$ be a sparsity pattern containing a $0$ for
an existing line, an arbitrarily large $K$ for an absent line, and a $1$
otherwise. The corresponding sparsity promoting prior is

$\displaystyle-\log(p(\hat{y}|\beta_{s}))=\lambda\|[\beta_{s}]\hat{y}\|_{1}.$
(28)

Prior knowledge may not directly concern the value of individual parameters.
One could also know the sign of the elements of $y$. The conductance
$\Re(Y_{hk})$ is always non-negative because of energy conservation laws, and
in the absence of series capacitance, the susceptance $\Im(Y_{hk})$ is non-
positive for all lines. The exact Bayesian prior for the sign would be a
uniform distribution with support $[0,y_{\max}]$. Combined with the sparsity,
it yields the exponential distribution. To avoid its discontinuity at zero, we
will use the asymmetric Laplace distribution [37] with a very large parameter
on one side as a smooth approximation. In the log-space, $\log(p(\hat{y}))$
becomes a sum of skewed absolute values written as

$\displaystyle-\log(p(\hat{y}|s))=\lambda\|\hat{y}\|_{1}+K\sum_{h=1}^{2n^{2}}(1-s_{h}\textrm{sgn}(\hat{y}_{h})),$
(29)

where $K$ is an arbitrarily large constant and $s$ is a vector with elements
in $\\{-1,1\\}$ defining the believed signs of $y$.

Constraining the ratios between parameters can be useful in many cases, for
example if two lines are parallel, or if one knows the type of cable used for
a certain line. For a ratio $\hat{\rho}_{hk}=\frac{y_{k}}{y_{h}}$, it gives

$\displaystyle-\log(p(\hat{y}|\hat{\rho}_{hk}))=\lambda\left\|(\rho_{hk}e_{h}-e_{k})^{\top}\hat{y}\right\|_{1}.$
(30)

While the exact value of the resistance $R$ and the reactance $X$ of a line
depends on its length, the type of cable gives their ratio $R$/$X$. A constant
$R$/$X$ ratio for the entire network is common assumption for distribution
grids. Let $\hat{\rho}$ be the estimated $R$/$X$ ratio, a prior enforcing
$\hat{\rho}$ is

$\displaystyle-\log(p(\hat{y}|\hat{\rho}))=\lambda\left\|[\hat{\rho}\mathcal{I}_{n^{2}},-\mathcal{I}_{n^{2}}]\hat{y}\right\|_{1}.$
(31)

We will not provide a direct comparison of all these prior because their
effects are highly dependent of the power grid studied. Moreover, they may not
all be applicable as they rely on specific types of knowledge that may be
absent.

###  6.3 Data-driven prior distributions

The Bayesian prior can model belief from other estimation methods or different
data sets, and is not limited to observations on the system. More
particularly, the MLE can be used to refine the priors from Section 6.2.
Sparsity promotion using (28) can create a bias, as any other $\ell_{1}$
penalties would. To reduce this bias, it can be weighted with
$L=[|\hat{y}_{\textrm{MLE}}|]^{-1}$, similarly to the adaptive Lasso method
[13]. This yieds

$\displaystyle-\log(p(\hat{y}|y_{\textrm{MLE}}))=\lambda\|[|\hat{y}_{\textrm{MLE}}|]^{-1}\hat{y}\|_{1}$
(32)

Note that this method is only asymptotically unbiased (with
$N\rightarrow\infty$). For a finite sample size, it does not cancel the
$\ell_{1}$ penalty’s bias but still improves it [38].

######  Remark.

For constructing the prior, $y_{\textrm{MLE}}$ is considered constant rather
than a random variable, preserving independence towards other potential
priors.

The hyperparameter $\lambda$ in (23) is very important because if it is too
big, then the bias will be unnecessarily large, and if it is too small, some
entries will not be effectively shrunk to zero. However, it can generally not
be tuned using cross validation (which is the standard approach in Lasso)
because $\Sigma_{v}$ and $\Sigma_{i}$ depend on the values of $\boldsymbol{v}$
and $\boldsymbol{i}$. For very sparse systems such as power grids, there may
not exist a $\lambda$ such that the MAP estimate is both sparse enough and
moderately biased. To solve this issue, one can add a prior on
$\|\hat{y}\|_{1}$, centered on its believed value $\hat{\gamma}$. If the sign
of every $y_{i}$ is known and with $s$ from (29), one can construct this prior
as

$\displaystyle-\log(p(\hat{y}|\hat{\gamma}))=\frac{\lambda^{\prime}}{\hat{\gamma}}|s^{\top}\hat{y}-\hat{\gamma}|.$
(33)

If $\lambda^{\prime}>\lambda$, it limits the bias created by a too large
$\lambda$, which makes its tuning much more tolerant to errors and allows a
more aggressive regularization.

The value of $\hat{\gamma}$ can be estimated as
$\|\hat{y}_{\textrm{MLE}}\|_{1}$, but this may be too large because of the
error on zero elements. With a large enough $n$ and a very sparse $y$, this
error comes mainly from zero elements. Let $\hat{y}_{\mymathbb{0}}$ be the
$n_{0}$ estimated parameters supposed to be of sign $s_{\mymathbb{0}}$ but
zero in reality (e.g. the susceptance of an inductive line that does not exist
corresponds to a $-1$ in $s_{\mymathbb{0}}$). Because the MLE is unbiased,
$E[\hat{y}_{\mymathbb{0}}]=\mymathbb{0}_{n_{0}}$ so
$E[s_{\mymathbb{0}}^{\top}\hat{y}_{\mymathbb{0}}]=E[\mymathbb{1}_{n_{0}}^{\top}\hat{y}_{\mymathbb{0}}]=0$.
Hence

$\displaystyle
E[(s_{\mymathbb{0}}-\mymathbb{1}_{\mymathbb{0}})^{\top}(\mymathbb{|}\hat{y}_{\mymathbb{0}}|-\hat{y}_{\mymathbb{0}})+(s_{\mymathbb{0}}+\mymathbb{1}_{\mymathbb{0}})^{\top}(\mymathbb{|}\hat{y}_{\mymathbb{0}}|+\hat{y}_{\mymathbb{0}})]=E[(s_{\mymathbb{0}}+\mymathbb{1}_{\mymathbb{0}})^{\top}(\mymathbb{|}\hat{y}_{\mymathbb{0}}|-\hat{y}_{\mymathbb{0}})+(s_{\mymathbb{0}}-\mymathbb{1}_{\mymathbb{0}})^{\top}(\mymathbb{|}\hat{y}_{\mymathbb{0}}|+\hat{y}_{\mymathbb{0}})],$
(34)

which means that half of the expected error on zero elements has the wrong
sign. Assuming that the error is smaller than the parameter’s value, no non-
zero element should have the wrong sign. This gives

$\displaystyle\hat{\gamma}\approx\|\hat{y}_{\textrm{MLE}}\|_{1}-|\hat{y}_{\textrm{MLE}}^{\top}||s-\textrm{sgn}(\hat{y}_{\textrm{MLE}})|.$
(35)

Even without using (33), this value can be useful to assess if the chosen
$\lambda$ is too large. Note that the same argument applies if
$\hat{y}_{\textrm{MLE}}$ is split into groups of elements (e.g. into columns).
However, if the groups are too small, the variance of each $\hat{\gamma}$ may
be very large, leading to an erroneous prior.

##  7 Structural priors

Due to the structure of power networks, $Y$ has peculiar properties. If the
network does not include phase-shifting transformers and power lines are not
compensated by series capacitors, $Y$ is symmetric. Moreover, for networks
where shunt elements are negligible, $Y$ is Laplacian [19]. Since phase-
shifting transformers are usually employed in transmission systems, and shunt
admittances are negligible for medium-sized grids, with line less than 60 km
long, it is safe to assume that standard distribution networks have a
Laplacian admittance matrix [39].

###  7.1 Formal definition

Under the assumption that $Y$ is Laplacian, entries on and above the main
diagonal of $Y$ can be derived from the elements below the diagonal.
Therefore, in order to avoid redundant variables, one can proceed as in [12]
and use duplication and transformation matrices $D$ and $T$ to remove the
redundant entries from the identification problem and solve for
$\hat{y}_{r}=[\Re(\operatorname*{\textrm{ve}}(\hat{Y}))^{\top},\Im(\operatorname*{\textrm{ve}}(\hat{Y}))^{\top}]^{\top}$
instead. In case some entries of $Y$ are known to be zero, one can derive
variants of $D$ and $T$ and also remove these zero entries from $\hat{y}$ by
following a procedure similar to the one presented in [12, Appendix 2]. In
both cases, $\hat{y}=(\mathcal{I}_{2}\otimes D\cdot T)\hat{y}_{r}$ and the
equation (18b) becomes

$\displaystyle\tilde{\boldsymbol{i}}-\Delta\boldsymbol{i}=(\tilde{\mymathbb}{V}-\Delta\mymathbb{V})(\mathcal{I}_{2}\otimes
D\cdot T)\hat{y},$ (36)

The diagonal entries of $Y$ are often the largest, as they are the sum of all
other entries on the same rows. According to Corollary 4 and its following
remark, the diagonal entries may not only have a large variance, but also
cause one for all other elements. Using the $D$ and $T$ matrices also improves
this last point.

###  7.2 Implications for Bayesian priors

The MAP estimate can be computed using (22) by replacing (22b) by (36). One
can also simply replace $\hat{y}$ by $\hat{y}_{r}$ in all the priors presented
in Section 6. However, the Laplacianity of $Y$ opens additional opportunities.
Similarly to Section 6.2, if the line susceptance is assumed negative,
$\textrm{diag}(Y_{\textrm{MLE}})$ provides the following alternative estimate
for $\hat{\gamma}$.

$\displaystyle\|\hat{y}\|_{1}$ $\displaystyle=\sum_{h=1}^{n}\sum_{k=1,k\neq
h}^{n}\left\|\begin{matrix}\Re(\hat{Y}_{hk})\\\
\Im(\hat{Y}_{hk})\end{matrix}\right\|_{1}+\sum_{h=1}^{n}\left\|\begin{matrix}\Re(\hat{Y}_{hh})\\\
\Im(\hat{Y}_{hh})\end{matrix}\right\|_{1},$
$\displaystyle=2\sum_{h=1}^{n}\left\|[\Re(\hat{Y}_{hh}),\Im(\hat{Y}_{hh})]\right\|_{1}.$
(37)

Note that $\|\hat{y}\|_{1}=4\|\hat{y}_{r}\|_{1}$ because in $\hat{y}$, the
diagonal elements make up half of the norm and the other elements are present
twice. Furthermore, a prior on $\sum_{k=1,k\neq h}^{n}Y_{hk}$ for all
$h\in\mathcal{V}$ can also limit the bias of a sparsity-promoting one. If the
exact value is not available and with a Laplacian $Y$, this prior can be
centered on $Y_{\textrm{hh,MLE}}$. This yields

$\displaystyle-\log(p(\hat{y}_{r}|\textrm{diag}(Y_{\textrm{MLE}})))=$ (38)
$\displaystyle\quad\lambda^{\prime}\sum_{h=1}^{n}\left|\frac{\sum_{k=1,k\neq
h}^{n}\Re(\hat{Y}_{hk})}{\Re(\hat{Y}_{\textrm{MLE},hh})}-1\right|+\left|\frac{\sum_{k=1,k\neq
h}^{n}\Im(\hat{Y}_{hk})}{\Im(\hat{Y}_{\textrm{MLE},hh})}-1\right|.$

More details about the representation of
$p(\hat{y}_{r}|\textrm{diag}(\hat{Y}_{\textrm{MLE}}))$ in the same shape as
(23) are presented in Appendix C.

In the simulation Section 9, we have a large network with a very sparse,
Laplacian admittance matrix. Hence, we will us $D$ and $T$ with the prior
$p(\hat{y}_{r}|s)$, as well as $p(\hat{y}_{r}|\hat{y}_{r,\textrm{MLE}})$ and
$p(\hat{y}_{r}|\textrm{diag}(\hat{Y}_{\textrm{MLE}}))$, which are built using
the non-diagonal and diagonal elements respectively. Both priors are combined
using the remarks in Sections 6.1 and 6.3.

##  8 Numerical methods

Before solving (22) to obtain an estimate, several improvements can be done by
pre-processing available measurements. This section will first explain how
centering and filtering the data, as well as removing hidden nodes can improve
the estimation. Then, we will present and compare different algorithms to
solve the optimization problem.

###  8.1 Data pre-processing

Some networks may have nodes with no load attached. If a node is unloaded,
then the corresponding column of $V$ is a linear combination of the columns
corresponding to neighboring nodes, as it is determined by a simple voltage
divider. In this case, $\Phi$ in (20) does not have full rank, so
$F_{\textrm{MLE}}$ is singular. In other words, the Cramer-Rao bound becomes
infinite, at least for some parameters, which means that they cannot be
reliably estimated reliably by any estimator. In a similar spirit, if all the
nodes are loaded but some nodes have much lower loads than others,
$F_{\textrm{MLE}}$ can be full rank, but some of its eigenvalues can be very
small. When inverting $F_{\textrm{MLE}}$, the small eigenvalues become very
large, which means that the variance of the corresponding estimates will also
be large.

From (25), it is apparent that a prior on lines near an unloaded node may help
with this issue as it potentially adds the missing ranks to
$F_{\textrm{MLE}}$. However, this makes the prior the only source of
information. Another solution is to take out these nodes from the problem, and
identify a reduced matrix $Y_{\textrm{red}}$ such that
$I_{-h}=Y_{\textrm{red}}V_{-h}$. In most applications, such as sensitivity
analysis or control, this reduced matrix is sufficient as it only removes
redundant parameters, while keeping an equivalent model for the remaining
nodes. $Y_{\textrm{red}}$ can be computed using the Kron reduction method
[40].

For a power grid with a rated voltage $V_{\textrm{rated}}$, the matrix
$\tilde{V}^{\top}\tilde{V}$ is almost equal to
$V_{\textrm{rated}}\mymathbb{1}_{n}\mymathbb{1}_{n}^{\top}$. This matrix is
then almost singular with one eigenvalue much larger than all others. This can
be an issue for numerical stability both for minimizing (19) or for solving
(22). If $Y$ is Laplacian, one can use
$\tilde{V}_{c}=\tilde{V}-V_{\textrm{rated}}\mymathbb{1}_{N}\mymathbb{1}_{n}^{\top}$
instead, as the second term is cancelled by
$\mymathbb{1}_{n}^{\top}Y=\mymathbb{0}_{n}^{\top}$. The covariance
$\Sigma_{v}$ still needs to be computed with $V$ and not $V_{c}$ due to the
non-linear transformation (8).

Finally, a low pass (moving average) filter can help reducing the measurement
noise. Measurements from µPMUs are usually very frequent (50 to 120Hz) [41],
and load variation has an average period of a couple minutes. A low pass
filter with a cutoff frequency equal to the one of the load profiles will not
remove a lot of information from the signal. However, as $K$ measurements are
averaged, the noise variance of a filtered measurement is reduced since
$[\sigma_{\varepsilon}^{2},\sigma_{\delta}^{2}]_{\textrm{filtered}}=K^{-1}[\sigma_{\varepsilon}^{2},\sigma_{\delta}^{2}]$.

###  8.2 Optimization algorithms

If $-log(p(\hat{y}))$ is convex, the optimization problem (22) has a convex
cost and bilinear constraints due to the multiplication of $\hat{y}$ and
$\Delta\mymathbb{V}$. Similarly to weighted TLS, no closed-form solution is
available [26].

The most basic algorithm for solving (22) is the alternate block coordinate
descent (BCD), which sets $\Delta\mymathbb{V}$ to constant to solve $\hat{y}$
for the next iteration ($k$) and vice versa, as explained in [42]. With
$c=-\mathcal{L}-\log(p(\hat{y}))$, the update is

$\displaystyle\Delta\mymathbb{V}_{k}$
$\displaystyle=\operatorname*{arg\,min}_{\Delta\mymathbb{V}}\;c\left(\hat{y}_{k-1},\Delta\mymathbb{V},\tilde{\boldsymbol{i}}-\left(\tilde{\mymathbb{V}}\\!-\\!\Delta\mymathbb{V}\right)\hat{y}_{k-1}\right),$
(39a) $\displaystyle\hat{y}_{k}$
$\displaystyle=\operatorname*{arg\,min}_{\hat{y}}\;c\left(\hat{y},\Delta\mymathbb{V}_{k},\tilde{\boldsymbol{i}}-\left(\tilde{\mymathbb{V}}\\!-\\!\Delta\mymathbb{V}_{k}\right)\hat{y}\right).$
(39b)

(39) shows two convex sub-problems that can be solved iteratively. However,
the $\hat{y}$-subproblem (39b) may not admit a closed-form solution, depending
on $p(\hat{y})$. When it does not, (39b) can become a computational
bottleneck.

To improve performance, one can use the approximation

$\|L(\hat{y}-\mu)\|_{1}\approx(\hat{y}-\mu)^{\top}L^{\top}[|L(\hat{y}_{k-1}-\mu)|+\alpha\mymathbb{1}]^{-1}L(\hat{y}-\mu)$

in the expression of $p(\hat{y})$, with a small enough $\alpha$. This
algorithm is called broken adaptive ridge regression (BAR)[43, 44], and
provides a closed-from approximate solution to (39b). If $L$ is diagonal,
another possible alternative is use an operator splitting method such as ADMM
[45, 46]. However, experimental evidence shows that more iterations are needed
for ADMM to converge. A comparison of the three algorithms presented here can
be found in table 2.

Algorithm | iterations to convergence | iterations/second
---|---|---
BCD |  $\sim$10000 | 1.25
BAR |  $\sim$10000 | 30
ADMM |  $\sim$30000 | 28
Table 2: Comparison of the speed of three algorithms on a 9 nodes network with
$N=400$ measurement samples, using a standard $\ell_{1}$ penalty with the same
$\lambda$. In all cases, they are executed on a MacBook Pro with a 2.3GHz
Intel i7 processor running Python 3.8.

######  Remark.

The convergence speed highly depends on the regularization parameter
$\lambda$. To speed it up, a higher penalty can be applied in the first
iterations, and then decrease to the optimal one. This optimal value can be
computed using $\hat{Y}_{\textrm{MLE}}$ as shown in Section 6.3.

If $L$ is diagonal, the proximal gradient method is also usable [46, 47], but
it requires to tune the step size on top of the regularization parameter. This
can take many iterations, especially if the information contained in the data
is very limited. Also, it relies on thresholding providing a closed-form
proximal operator of the $\ell_{1}$ norm [48]. If $L$ is not diagonal, the
proximal operator becomes a piece-wise linear function with a number of pieces
scaling with the number of combinations of signs in $L$. It therefore becomes
quickly prohibitive to compute. In practice, non-diagonal priors are used for
cancelling bias (Section 6.3) or to keep some parameters close to the same
value, which can be needed for a three phased network if some segments have
three times the same line. The proximal Newton method [47] is not usable for
the same reason: it relies on the transformation of the proximal operator by
the Hessian matrix of $c$, which is generally not diagonal.

##  9 Simulation results

We apply the estimation method presented in Section 6, as well as approaches
from other works, to a 56-node network with realistic parameters, admittances,
noise levels and load profiles. In order to obtain the results, the BAR
algorithm is implemented using a hardware-accelerated linear solver from
NVidia’s CUDA tool. The framework is programmed in Python and available on
GitHub [49]. It simulates the voltage and current measurements from the real
network parameters, nominal loads, and given household load profiles. It then
computes the various least squares and regularized MLE estimates from the
simulated measurements.

###  9.1 Setup

To simulate the identification problem, we use the IEEE 123 bus network.
Assuming that the loads on each phase are balanced, its three-phased part is
transformed into a 56 nodes, single-phase feeder using the method and
parameters from [50]. Although our identification method could be applied to a
three-phase network, we will use this simplified version for convenience and
clarity. Additionally, some nodes do not have connected loads. We use the Kron
reduction method described in section 8.1 to remove them and compare the
reduced matrices.

The load profiles for each node are generated with the GENETX generator [51].
It creates random realistic loads for households according to parameters like
penetration of renewable energies (set to 35%) or electric vehicles (40%). The
tool then generates a thousand one-minute-resolution demand profiles, for
households situated in the Netherlands from weeks 12 to 16 of the year. To
create the load profiles for each node (Fig. 1(b)), the demand profiles of
randomly selected households are summed until the nominal power is reached.

(a) Map of the network from [50].

$0$$100$$200$$300$$400$$500$$600$$700$$800$$900$$1{,}000$$1{,}100$$1{,}200$$1{,}300$$1{,}400$$1{,}500$$0$$0.1$$0.2$time
$[\textrm{min}]$Active load $[\textrm{MW}]$ (b) Sample day of generated load
profiles, with node 4’s one highlighted.

Figure 1: Graphic representation of the simulation settings.

Voltage and current values are generated by simulating the network using the
PandaPower library [52] with a measurement frequency of 100Hz. We then add
0.01% of Gaussian noise in polar coordinates, as described in Section 4. Note
that the noise generated by a $\mu$PMU depends on its rating. Assuming that
the $\mu$PMUs are adapted to their nodes, we choose a rating of four times the
nominal power.

Because we only have data on loads every minute and the measurements are taken
at 100Hz, we extrapolate the missing load data using linear interpolation of
consecutive samples. Using all samples collected during an entire month at
100Hz for network identification is not computationally practical. Therefore,
as proposed in Section 8.1, we use a moving average filter with a length of
17000 and under-sample the result by the same amount. This provides 15000
samples of both current and voltages for each node as complex phasors. A
higher amount of data would probably provide a better estimate, but this is
limited by computational complexity. It is also the main limit when the
network size increases.

###  9.2 Results

Thanks to the high accuracy of µPMUs, and to the data preprocessing described
in Section 8.1, the MLE method manages to retrieve an approximate but fair
estimate (Fig. 2(a)). Although the sparsity of the admittance matrix is above
98%, this estimate is dense. Fig. 2(b) shows the MAP estimate with a prior
distribution as described in Section 7.2. It does not use any exact
information, but only the MLE estimate described above as a starting point for
the BAR algorithm. While the following relative Frobenius error (40) is more
than 50% when using the MLE, the MAP yields 10.3%.

$\displaystyle\varepsilon_{F}=\frac{\|\hat{Y}-Y\|_{F}}{\|Y\|_{F}}.$ (40)

We also compare the MAP with the OLS and TLS estimates (used as a baseline),
as well as adaptive Lasso from [13]. We recall that both OLS and Lasso do not
consider errors in the voltage variables. Fig. 3 shows the magnitude of the
estimated admittance for all 55 existing lines, as well as the exact value.
Note that the MLE has been omitted because it is very close to the TLS
estimate, and that this plot does not allow one to compare the density of the
estimated admittance matrix.

(a) MLE estimate.

(b) MAP estimate.

(c) Real reduced admittance matrix.

Figure 2: Heat-maps of real and estimated admittance matrices (elements on the diagonal have been omitted, see Section 7). $0$$5$$10$$15$$20$$25$$30$$35$$40$$45$$50$$55$$0$$0.5$$1$$1.5$$\cdot 10^{4}$exactOLSTLSadaptive LassoMAP Figure 3: Comparison of the magnitude of non-zero elements estimated with different methods. Noise std. | MAP | TLS | Lasso | OLS
---|---|---|---|---
2e-4 | 40.69% | 167.24% | 99.54% | 99.58%
1e-4 | 10.32% | 51.90% | 98.48% | 98.56%
1e-5 | 0.91% | 4.86% | 60.17 | 60.78%
1e-6 | 0.073% | 0.49% | 1.54% | 2.18%
Table 3: Estimation error for various levels of noise, with standard
deviations in phase and magnitude ($\sigma_{\varepsilon}=\sigma_{\delta})$.

Fig. 2(b), corresponding to the noise level in the second row of Table 3
(1e-4), and Fig. 2(c) show that the MAP estimate is not very far from the
actual admittance. However, Table 3 shows that by doubling the noise, the
estimation error grows to 40%. Overall, the error growth is non-linear with
the noise standard deviation, and while it increases fast with more noise, it
also decreases fast with less. We can also observe that the limits of OLS and
Lasso are reached with lower noise.

###  9.3 Prior improvement by measuring selected line admittances

If one needs to further improve the estimate, the remaining solution is to
collect further information on topology or parameters of the electric network
and integrate the results into the Bayesian prior (Section 6.2). To avoid
collecting too much data, the additional information can be focused on
estimates with the largest error. Although the error covariance is practically
not computable (Section 5.3), its properties show that in theory one expects
the estimation error to be concentrated on (i) high admittance elements such
as short lines or switches, (ii) lines connecting a node with a low load, and
(iii) lines with a very low admittance, which may be strongly affected by
regularization and estimated as zero.

The first point follows from Corollary 4, which implies that a high admittance
element also affects the estimation error of the other lines connected to the
same node. Point (ii) is equivalent to the remark in Section 8.1, that if a
node has a low load, the error variances of the admittances of all lines
connected to it may be very high.

These guidelines can help choosing lines to measure to achieve a better
estimate. Table 4 shows the MAP estimation errors for various noise levels and
three different priors. The first is the same as in Section 9.2. The second
incorporates the knowledge of the line from node 4 to node 40 and the lines
connected to nodes 2, 50 and 51. Indeed, these nodes have some of the lowest
loads, hence falling in the category (i). In total, we add information about 7
of the 1540 possible connections, among which 55 actually exist. As a
comparison, the third prior adds the measurements of 7 random other lines to
the first one. Note that inserting prior information has a larger impact if
the estimate is worse.

Noise | no measure | chosen lines | random lines
---|---|---|---
2e-4 | 40.69% | 29.39% | 38.97%
1.7e-4 | 24.88% | 20.41% | 23.83%
1.3e-4 | 14.75% | 12.88% | 14.13%
1e-4 | 10.32% | 9.06% | 9.89%
Table 4: MAP estimation error with exact prior for various levels of noise,
equal in phase and magnitude for current and voltages.

##  10 Conclusion

The penetration of distributed generation and smart devices in the
distribution grid calls for the introduction of advanced control schemes,
which require the exact topology and line parameters. Such information is
often unavailable for distribution networks: as direct measurement is
infeasible, data-driven estimators are needed.

In this work, we proposed to exploit samples collected by micro-PMUs.
Considering a realistic statistical model for the noise affecting both current
and voltage measurements, we built maximum-likelihood and Bayesian estimators.
We argued that the latter can outperform the former, due to their ability to
exploit features of the grid, such as sparsity, as well as available
information on specific lines. Our argument is substantiated by numerical
simulations on benchmark grids: even without any network-specific prior
information, Bayesian methods outperformed state-of-the-art estimators with
realistic noise levels.

Further research may focus on effective methods to selectively collect live
measurements for improving the quality of estimates, as well as the
development of alternative noise models for different sensors and physical
quantities such as power. Formally defining the network’s observability could
lift the limitation to fully observed grids, introducing pseudo-measurements
and Bayesian priors on the grid’s state. Learning a reduced network connecting
only specific nodes could also provide an answer to missing measurements.

##  Appendix A Covariance matrix

In order to solve the Maximum Likelihood problem (18) and all its subsequent
refinements, one needs the covariance matrices $\Sigma_{v}$ and $\Sigma_{i}$.
The construction of the two is identical, thus we will focus on
$\Sigma_{v}=\operatorname*{Cov}[\Delta\boldsymbol{v}]\in{R}^{2nN\times 2nN}$
only.

From Section 4, the $\Sigma_{v}$ is sparse, having non-zero elements only on
three diagonals:

$\displaystyle\Sigma_{v}=\begin{pmatrix}\begin{matrix}\textrm{Var}[\Re(\tilde{V}_{11})]&\lx@intercol\hfil\text{\kern
5.0pt\smash{\raisebox{-4.30554pt}{$$}}}\hfil\lx@intercol\\\ &\ddots&\\\
\lx@intercol\hfil\text{\kern
5.0pt\smash{\raisebox{4.30554pt}{$$}}}\hfil\lx@intercol&\textrm{Var}[\Re(\tilde{V}_{Nn})]\end{matrix}&\vline&\begin{matrix}\textrm{Cov}[\Re(\tilde{V}_{11}),\Im(\tilde{V}_{11})]&\lx@intercol\hfil\text{\kern
5.0pt\smash{\raisebox{-4.30554pt}{$$}}}\hfil\lx@intercol\\\ &\ddots&\\\
\lx@intercol\hfil\text{\kern
5.0pt\smash{\raisebox{4.30554pt}{$$}}}\hfil\lx@intercol&\textrm{Cov}[\Re(\tilde{V}_{Nn}),\Im(\tilde{V}_{Nn})]\end{matrix}\\\
\hline\cr\star&\vline&\begin{matrix}\textrm{Var}[\Im(\tilde{V}_{11})]&\lx@intercol\hfil\text{\kern
5.0pt\smash{\raisebox{-4.30554pt}{$$}}}\hfil\lx@intercol\\\ &\ddots&\\\
\lx@intercol\hfil\text{\kern
5.0pt\smash{\raisebox{4.30554pt}{$$}}}\hfil\lx@intercol&\textrm{Var}[\Im(\tilde{V}_{Nn})]\end{matrix}\end{pmatrix}.$
(41)

The main diagonal of $\Sigma_{v}$ hosts the variance while the $nN$th super-
and sub-diagonals, provide the covariance between the real and the imaginary
part of the measurements. Such particular structure makes it possible to split
$\Sigma_{v}$ into four diagonal blocks and makes it easy to find
$\Sigma_{v}^{-1}$ analytically.

It is also interesting to note that, up to a permutation of the elements in
$\Delta\boldsymbol{v}$, $\operatorname*{Cov}[\Delta\boldsymbol{v}]$ can be
written as a block diagonal matrix where the 2-by-2 blocks are given by (8c).

##  Appendix B Proof in Section 5.3

###  B.1 Proof of Lemma 3

$\Phi=\mathcal{I}_{n}\otimes(V-_{N}\bar{V})$ being a Kronecker product with
identity, $\Phi_{hN+t}^{\top}\Phi_{hN+t}$ is block diagonal for all $h,t$,
with $n$ blocks of size $N$ by $n$. Hence, $F_{\textrm{MLE}}$ is a 2-by-2
block matrix, with each block being a block diagonal matrix with $n$ blocks,
and its inverse has the same shape. The error covariance matrix $\Sigma_{y}$
is then written using the compact notation
$\hat{Y}_{h}^{\Re}+j\hat{Y}_{h}^{\Im}=(\hat{Y}^{\top})_{h,\textrm{MLE}}$ as

$\displaystyle\Sigma_{y}=\begin{pmatrix}\begin{matrix}\textrm{Var}[\hat{Y}_{1}^{\Re}]&\lx@intercol\hfil\text{\kern
5.0pt\smash{\raisebox{-4.30554pt}{$$}}}\hfil\lx@intercol\\\ &\ddots&\\\
\lx@intercol\hfil\text{\kern
5.0pt\smash{\raisebox{4.30554pt}{$$}}}\hfil\lx@intercol&\textrm{Var}[\hat{Y}_{n}^{\Re}]\end{matrix}&\vline&\begin{matrix}\textrm{Cov}[\hat{Y}_{1}^{\Re},\hat{Y}_{1}^{\Im}]&\lx@intercol\hfil\text{\kern
5.0pt\smash{\raisebox{-4.30554pt}{$$}}}\hfil\lx@intercol\\\ &\ddots&\\\
\lx@intercol\hfil\text{\kern
5.0pt\smash{\raisebox{4.30554pt}{$$}}}\hfil\lx@intercol&\textrm{Cov}[\hat{Y}_{n}^{\Re},\hat{Y}_{n}^{\Im}]\end{matrix}\\\
\hline\cr\star&\vline&\begin{matrix}\textrm{Var}[\hat{Y}_{1}^{\Im}]&\lx@intercol\hfil\text{\kern
5.0pt\smash{\raisebox{-4.30554pt}{$$}}}\hfil\lx@intercol\\\ &\ddots&\\\
\lx@intercol\hfil\text{\kern
<EMAIL_ADDRESS>(42)

The Kronecker product constructing $\Phi$ also includes $n(n-1)$ deterministic
zeros to each of their rows, which then have a zero variance and covariance.
Together with assumption 2, this yields

$\displaystyle\mathcal{R}_{\Re,1ht}=\begin{pmatrix}\begin{matrix}[_{n(h-1)}]&\\\
&\textrm{Var}[\Re(\tilde{V}_{t})]\end{matrix}&\text{$$}\\\
\text{$$}&\begin{matrix}[_{n(n-h)}]&\\\
&\textrm{Var}[\Re(\tilde{I}_{ht})]\end{matrix}\end{pmatrix},$ (43)

where $\textrm{Var}[\Re(\tilde{V}_{t})]$ is diagonal from 2. Thus,
$\mathcal{R}_{\Re,1ht}$ is also diagonal. $\mathcal{R}_{\Re,2ht}$ has the same
expression as (43), but with $\Im(\tilde{I}_{ht})$ replacing
$\Re(\tilde{I}_{ht})$. Similar expressions can be derived using
$\textrm{Cov}[\Re(\tilde{V}_{t}),\Im(\tilde{V}_{t})]$ and
$\textrm{Var}[\Im(\tilde{V}_{t})]$ for $\mathcal{R}_{\Re\Im,qht}$ and
$\mathcal{R}_{\Im,qht}$, respectively. With $q\in\\{1,2\\}$ and
$\mathcal{Q}_{qht}$ such that
$\mathcal{Q}_{1ht}=\textrm{Var}[\Re(\tilde{I}_{ht})]$ and
$\mathcal{Q}_{2ht}=\textrm{Var}[\Im(\tilde{I}_{ht})]$, we define

$\displaystyle\mathcal{D}_{\Re,qht}(Y_{\cdot h})$ $\displaystyle=\Re(Y_{\cdot
h})^{\top}\textrm{Var}[\Re(\tilde{V}_{t})]\Re(Y_{\cdot h})+\mathcal{Q}_{qht},$
(44a) $\displaystyle\mathcal{D}_{\Re\Im,qht}(Y_{\cdot h})$
$\displaystyle=\Re(Y_{\cdot
h})^{\top}\textrm{Cov}[\Re(\tilde{V}_{t}),\Im(\tilde{V}_{t})]\Im(Y_{\cdot
h})+\mathcal{Q}_{qht},$ (44b) $\displaystyle\mathcal{D}_{\Im,qht}(Y_{\cdot
h})$ $\displaystyle=\Im(Y_{\cdot
h})^{\top}\textrm{Var}[\Im(\tilde{V}_{t})]\Im(Y_{\cdot h})+\mathcal{Q}_{qht},$
(44c)

Excluding the zeros in (43) yields
$D_{\Re,1ht}(Y_{h})=\Re(z)^{\top}\mathcal{R}_{\Re,1ht}\Re(z)$, and similarly
for $D_{Q,qht}(Y_{h})$ with any $Q\in\\{\Re,\Re\Im,\Im\\}$ and
$q\in\\{1,2\\}$.

Moreover, $\Phi_{hN+t}^{\top}\Phi_{hN+t}$ also has the same sparsity pattern
as the $\mathcal{R}$ matrices (excluding diagonality and the last row and
column). Therefore, for any $k$, the $k^{th}$ $n$-by-$n$ block of
$F_{\textrm{MLE}}$ is only nonzero for the terms of the sum (20) where $h=k$.
This means that

$\displaystyle\textrm{Var}\\!\left[\begin{matrix}\hat{Y}_{h}^{\Re}\\\
\hat{Y}_{h}^{\Im}\end{matrix}\right]\\!=\\!\left(\\!\sum_{q,t=1}^{2,N}\begin{matrix}\frac{\Re(\Phi_{hN+t})^{\top}\Re(\Phi_{hN+t})}{D_{\Re,qht}(Y_{\cdot
h})}&\frac{\Re(\Phi_{hN+t})^{\top}\Im(\Phi_{hN+t})}{D_{\Re\Im,qht}(Y_{\cdot
h})}\\\
\star&\frac{\Im(\Phi_{hN+t})^{\top}\Im(\Phi_{hN+t})}{D_{\Im,qht}(Y_{\cdot
h})}\end{matrix}\right)^{-1}.$ (45)

Equations (42) and (45) respectively show not only that the columns
$(\hat{Y}^{\top})_{h}$ of $\hat{Y}_{\textrm{MLE}}$ are statistically
independent, but also that their variance does not depend on the exact values
of one another, which finishes the proof.

###  B.2 Proof of Corollary 4

Let $Y_{\cdot h}^{1}$ and $Y_{\cdot h}^{2}$ be such that $|\Re(Y_{\cdot
h}^{2})|\geq|\Re(Y_{\cdot h}^{1})|$ and $|\Im(Y_{\cdot
h}^{2})|\geq|\Im(Y_{\cdot h}^{1})|$ element-wise. From (44) and assuming
$\textrm{Cov}[\Re(V_{t}),\Im(V_{t})]=0$, it follows that, for all $q$ and $t$,

$\displaystyle\mathcal{D}_{\Re,qht}(Y_{\cdot h}^{2})$ $\displaystyle\geq
D_{\Re,pit}(Y_{\cdot h}^{1}),$ (46a)
$\displaystyle\mathcal{D}_{\Re\Im,qht}(Y_{\cdot h}^{2})$
$\displaystyle=\mathcal{D}_{\Re\Im,qht}(Y_{\cdot h}^{1}),$ (46b)
$\displaystyle\mathcal{D}_{\Im,qht}(Y_{\cdot h}^{2})$
$\displaystyle\geq\mathcal{D}_{\Im,qht}(Y_{\cdot h}^{1}),$ (46c)

because $\textrm{Var}[\Re(\tilde{V}_{t})]$ and
$\textrm{Var}[\Im(\tilde{V}_{t})]$ are positive diagonal matrices. Using (45),
we then write

$\displaystyle\left(\textrm{Var}\\!\left[\begin{matrix}\hat{Y}_{h}^{2,\Re}\\\
\hat{Y}_{h}^{2,\Im}\end{matrix}\right]\right)^{-1}-\left(\textrm{Var}\\!\left[\begin{matrix}\hat{Y}_{h}^{1,\Re}\\\
\hat{Y}_{h}^{1,\Im}\end{matrix}\right]\right)^{-1}=$ (47)
$\displaystyle\quad\sum_{q,t=1}^{2,N}\left(\begin{matrix}\alpha^{\Re}\Re(\Phi_{hN+t})^{\top}\Re(\Phi_{hN+t})&0\\\
0&\alpha^{\Im}\Im(\Phi_{hN+t})^{\top}\Im(\Phi_{hN+t})\end{matrix}\right),$

with $\alpha^{\Re}=D_{\Re,qht}(Y_{\cdot h}^{2})^{-1}-D_{\Re,qht}(Y_{\cdot
h}^{1})^{-1}$ and $\alpha^{\Im}=D_{\Im,qht}(Y_{\cdot
h}^{2})^{-1}-D_{\Im,pit}(Y_{\cdot h}^{1})^{-1}$. The inequalities in (46) show
that both $\alpha^{\Re}$ and $\alpha^{\Im}$ are non-positive. Hence, the
blocks of (47) are the product of a negative scalar and a quadratic form and
are therefore negative semi-definite. From this definiteness follows that

$\displaystyle\textrm{Var}\\!\left[\begin{matrix}\hat{Y}_{h}^{2,\Re}\\\
\hat{Y}_{h}^{2,\Im}\end{matrix}\right]\succeq\textrm{Var}\\!\left[\begin{matrix}\hat{Y}_{h}^{1,\Re}\\\
\hat{Y}_{h}^{1,\Im}\end{matrix}\right],$ (48)

which finishes the proof.

##  Appendix C Non diagonal Bayesian prior

The goal of the regularization of diagonal elements of $\hat{Y}$ is such that
for all $h$,

$\displaystyle\hat{Y}_{\textrm{MLE},hh}$
$\displaystyle=\frac{1}{2}_{n^{2}}^{\top}\operatorname*{\textrm{vec}}\left(\begin{array}[]{c|c|c}&\hat{Y}_{h,:h-1}^{\top}&\\\
\hline\cr\vspace{-11px}\hfil\\\ \hat{Y}_{h,:h-1}&0&\hat{Y}_{h,h+1:}\\\
\hline\cr\vspace{-11px}\hfil\\\ &\hat{Y}_{h,h+1:}^{\top}&\end{array}\right),$
$\displaystyle=[1,j]\otimes\left(\operatorname*{\textrm{ve}}(e_{h}{}_{n}^{\top}+_{n}e_{h}^{\top})\right)^{\top}\hat{y}_{r},$

where $[\hat{Y}_{h,:h-1},\hat{Y}_{hh},\hat{Y}_{h,h+1:}]=\hat{Y}_{h}$ is the
$h$th row of $\hat{Y}$. The non-diagonal prior is then given by
$\mu_{\textrm{nd}}=_{n}$ and $L_{\textrm{nd}}$ defined as:

$\displaystyle L_{\textrm{nd}}$
$\displaystyle=\frac{\lambda^{\prime}}{\lambda}\\!\left[\left(\begin{matrix}\Re(\hat{Y}_{\textrm{MLE},hh}^{-1})&0\\\
0&\Im(\hat{Y}_{\textrm{MLE},hh}^{-1})\end{matrix}\right)\\!\otimes\operatorname*{\textrm{ve}}(e_{h}{}_{n}^{\top}+_{n}e_{h}^{\top})\right]_{h=1}^{n}.$

##  Appendix D Noise bias in Cartesian coordinates

In Section 4 we stated that the noise bias Eq. 7 is negligible with realistic
levels of accuracy. In order to show this, we adapt the procedure in [25, Sec
III.A]. By using first-order Taylor expansion on Eq. 7, we get:

$\displaystyle\operatorname*{E}[\Delta c|\tilde{v},\tilde{\theta}]$
$\displaystyle\simeq-(\tilde{v}\sigma^{2}_{\delta}/2)\cos\tilde{\theta},$
(49a) $\displaystyle\operatorname*{E}[\Delta d|\tilde{v},\tilde{\theta}]$
$\displaystyle\simeq-(\tilde{v}\sigma^{2}_{\delta}/2)\sin\tilde{\theta}.$
(49b)

Then,

$\lVert[\operatorname*{E}[\Delta
c|\tilde{v},\tilde{\theta}],\operatorname*{E}[\Delta
d|\tilde{v},\tilde{\theta}]]\rVert=\tilde{v}\sigma^{2}_{\delta}/2.$ (50)

Moreover, the minimum eigenvalue $\lambda_{\min}$ of the covariance matrix (8)
is $\min(\sigma^{2}_{\epsilon},\tilde{v}^{2}\sigma^{2}_{\delta})$: therefore,
the minimum standard deviation in the covariance matrix is
$\sigma_{\text{min}}=\sqrt{\lambda_{\text{min}}}=\sqrt{\min(\sigma^{2}_{\epsilon},\tilde{v}^{2}\sigma^{2}_{\delta})}$.

The bias can be considered non-significant if $\lVert[\operatorname*{E}[\Delta
c|\tilde{v},\tilde{\theta}],\operatorname*{E}[\Delta
d|\tilde{v},\tilde{\theta}]]\rVert/\sigma_{\text{min}}$ is small. Adopting the
per-unit system, and using realistic $\mu$PMU accuracy specifications (Table
1), we obtain:

$\frac{\lVert[\operatorname*{E}[\Delta
c|\tilde{v},\tilde{\theta}],\operatorname*{E}[\Delta
d|\tilde{v},\tilde{\theta}]]\rVert}{\sigma_{\text{min}}}\simeq 8.72\cdot
10^{-5}.$ (51)

Thus, the bias is four orders of magnitude smaller than the smallest standard
deviation from the noise covariance matrix and therefore can be safely
neglected.

## References

  * [1] Y. Weng, Y. Liao, and R. Rajagopal.  Distributed energy resources topology identification via graphical modeling.  IEEE Transactions on Power Systems, 32(4):2682–2694, 2017.
  * [2] L. Schenato, G. Barchi, D. Macii, R. Arghandeh, K. Poolla, and A. Von Meier.  Bayesian linear state estimation using smart meters and pmus measurements in distribution grids.  In 2014 IEEE International Conference on Smart Grid Communications (SmartGridComm), pages 572–577, 2014.
  * [3] A. La Bella, S. Negri, R. Scattolini, and E. Tironi.  A two-layer control architecture for islanded ac microgrids with storage devices.  In 2018 IEEE Conference on Control Technology and Applications (CCTA), pages 1421–1426, 2018.
  * [4] Alessio Iovine, Tristan Rigaut, Gilney Damm, Elena De Santis, and Maria Domenica Di Benedetto.  Power management for a dc microgrid integrating renewables and storages.  Control Engineering Practice, 85:59–79, 2019.
  * [5] Alessandra Parisio, Evangelos Rikos, and Luigi Glielmo.  Stochastic model predictive control for economic/environmental operation management of microgrids: An experimental case study.  Journal of Process Control, 43:24–37, 2016.
  * [6] Deepa S Kumar, JS Savier, and SS Biju.  Micro-synchrophasor based special protection scheme for distribution system automation in a smart city.  Protection and Control of Modern Power Systems, 5(1):1–14, 2020\.
  * [7] Guido Cavraro and Vassilis Kekatos.  Graph algorithms for topology identification using power grid probing.  IEEE control systems letters, 2(4):689–694, 2018.
  * [8] Guido Cavraro and Vassilis Kekatos.  Inverter probing for power distribution network topology processing.  IEEE Transactions on Control of Network Systems, 6(3):980–992, 2019\.
  * [9] Xu Du, Alexander Engelmann, Yuning Jiang, Timm Faulwasser, and Boris Houska.  Optimal experiment design for ac power systems admittance estimation.  arXiv preprint arXiv:1912.09017, 2019.
  * [10] M. Angjelichinoski, C. Stefanović, P. Popovski, A. Scaglione, and F. Blaabjerg.  Topology identification for multiple-bus dc microgrids via primary control perturbations.  In 2017 IEEE Second International Conference on DC Microgrids (ICDCM), pages 202–206, June 2017.
  * [11] Ye Yuan, Steven Low, Omid Ardakanian, and Claire Tomlin.  Inverse power flow problem.  arXiv preprint arXiv:1610.06631, 2016.
  * [12] Emanuele Fabbiani, Pulkit Nahata, Giuseppe De Nicolao, and Giancarlo Ferrari-Trecate.  Identification of ac networks via online learning.  arXiv preprint arXiv:2003.06210, 2020.
  * [13] O. Ardakanian, V. W. S. Wong, R. Dobbe, S. H. Low, A. von Meier, C. J. Tomlin, and Y. Yuan.  On identification of distribution grids.  IEEE Transactions on Control of Network Systems, 6(3):950–960, 2019\.
  * [14] Bogdan Pinte, Michael Quinlan, and Karl Reinhard.  Low voltage micro-phasor measurement unit ($\mu$pmu).  In 2015 IEEE Power and Energy Conference at Illinois (PECI), pages 1–4. IEEE, 2015.
  * [15] Styliani Sarri.  Methods and performance assessment of PMU-based real-time state estimation of active distribution networks.  PhD thesis, EPFL, 2016.
  * [16] Jiafan Yu, Yang Weng, and Ram Rajagopal.  Patopa: A data-driven parameter and topology joint estimation framework in distribution grids.  IEEE Transactions on Power Systems, PP, 05 2017.
  * [17] Antoine Wehenkel, Arpan Mukhopadhyay, Jean-Yves Le Boudec, and Mario Paolone.  Parameter estimation of three-phase untransposed short transmission lines from synchrophasor measurements.  IEEE Transactions on Instrumentation and Measurement, 69(9):6143–6154, 2020.
  * [18] Allen J Wood, Bruce F Wollenberg, and Gerald B Sheblé.  Power generation, operation, and control.  John Wiley & Sons, 2013.
  * [19] P. Kundur.  Power System Stability and Control.  CRC Press New York, NY, USA, 2007.
  * [20] Florian Dörfler, John W Simpson-Porco, and Francesco Bullo.  Electrical networks and algebraic graph theory: Models, properties, and applications.  Proceedings of the IEEE, 106(5):977–1005, 2018.
  * [21] Alexandra Von Meier, Emma Stewart, Alex McEachern, Michael Andersen, and Laura Mehrmanesh.  Precision micro-synchrophasors for distribution systems: A summary of applications.  IEEE Transactions on Smart Grid, 8(6):2926–2936, 2017.
  * [22] Simon J Julier and Jeffrey K Uhlmann.  Unscented filtering and nonlinear estimation.  Proceedings of the IEEE, 92(3):401–422, 2004.
  * [23] Mo Longbin, Song Xiaoquan, Zhou Yiyu, Sun Zhong Kang, and Yaakov Bar-Shalom.  Unbiased converted measurements for tracking.  IEEE Transactions on Aerospace and Electronic Systems, 34(3):1023–1027, 1998.
  * [24] Zhansheng Duan, Chongzhao Han, and X Rong Li.  Comments on "unbiased converted measurements for tracking".  IEEE transactions on aerospace and electronic systems, 40(4):1374, 2004.
  * [25] Don Lerro and Yaakov Bar-Shalom.  Tracking with debiased consistent converted measurements versus ekf.  IEEE transactions on aerospace and electronic systems, 29(3):1015–1022, 1993.
  * [26] Ivan Markovsky and Sabine Van Huffel.  Overview of total least-squares methods.  Signal Processing, 87(10):2283–2302, 2007.  Special Section: Total Least Squares and Errors-in-Variables Modeling.
  * [27] Sabine Van Huffel and Vandewalle Joos.  The Total Least Squares Problem: Computational Aspects and Analysis, chapter 8, pages 227–250.  Society for Industial and Applied Mathematics, 1991.
  * [28] Giuseppina Casella and R. Berger.  Statistical Inference.  Duxbury Press, 01 2002.
  * [29] Robert Engle and Daniel McFadden, editors.  Handbook of Econometrics, volume 4.  Elsevier, 1 edition, 1986.
  * [30] John Crassidis and Yang Cheng.  Error-covariance analysis of the total least squares problem.  Journal of Guidance, Control, and Dynamics, 37, 07 2014.
  * [31] Xing Fang, Bofeng Li, Hamza Alkhatib, Wenxian Zeng, and Yibin Yao.  Bayesian inference for the errors-in-variables model.  Studia Geophysica et Geodaetica, 61:35–52, 02 2017.
  * [32] Paul Holland.  Weighted ridge regression: Combining ridge and robust regression methods.  10 1973.
  * [33] Trevor Park and George Casella.  The bayesian lasso.  Journal of the American Statistical Association, 103(482):681–686, 2008.
  * [34] Dimitris Bertsimas and Martin S. Copenhaver.  Characterization of the equivalence of robustification and regularization in linear and matrix regression.  European Journal of Operational Research, 270(3):931–942, 2018\.
  * [35] Minjung Kyung, Jeff Gill, Malay Ghosh, and George Casella.  Penalized regression, standard errors, and bayesian lassos.  Bayesian Analysis, 5:369–412, 06 2010.
  * [36] Michael Osborne, Brett Presnell, and Berwin Turlach.  On the lasso and its dual.  Journal of Computational and Graphical Statistics, 9:319–337, 07 2000.
  * [37] Samuel Kotz, Tomaz J. Kozubowski, and Krzysztof Podgórski.  Asymmetric Laplace Distributions, pages 133–178.  Birkhäuser Boston, Boston, MA, 2001.
  * [38] Hui Zou.  The adaptive lasso and its oracle properties.  Journal of the American Statistical Association, 101(476):1418–1429, 2006.
  * [39] M. Taleb, M. J. Ditto, and T. Bouthiba.  Performance of short transmission lines models.  In 2006 IEEE GCC Conference (GCC), pages 1–7, March 2006.
  * [40] Florian Dörfler and Francesco Bullo.  Kron reduction of graphs with applications to electrical networks.  Computing Research Repository - CORR, 60, 02 2011.
  * [41] Emile Dusabimana and Sung-Guk Yoon.  A survey on the micro-phasor measurement unit in distribution networks.  Electronics, 9(2), 2020.
  * [42] Hao Zhu, Geert Leus, and G.B. Giannakis.  Sparsity-cognizant total least-squares for perturbed compressive sampling.  Signal Processing, IEEE Transactions on, 59:2002 – 2016, 06 2011\.
  * [43] Linlin Dai, Kani Chen, Zhihua Sun, Zhenqiu Liu, and Gang Li.  Broken adaptive ridge regression and its asymptotic properties.  Journal of Multivariate Analysis, 168, 08 2018.
  * [44] Florian Frommlet and Gregory Nuel.  An adaptive ridge procedure for l0 regularization.  PloS one, 11, 05 2015.
  * [45] Stephen Boyd, Neal Parikh, Eric Chu, Borja Peleato, and Jonathan Eckstein.  Distributed optimization and statistical learning via the alternating direction method of multipliers.  Foundations and Trends in Machine Learning, 3:1–122, 01 2011.
  * [46] Amir Beck.  First-Order Methods in Optimization.  SIAM-Society for Industrial and Applied Mathematics, Philadelphia, PA, USA, 2017.
  * [47] Nicholas G. Polson, James G. Scott, and Brandon T. Willard.  Proximal algorithms in statistics and machine learning.  arXiv preprint arXiv:1502.03175, 2015.
  * [48] Masayuki Tanaka and Masatoshi Okutomi.  Unified optimization framework for L2, L1, and/or L0 constrained image reconstruction.  In Computational Imaging II, volume 10222, pages 99 – 108. SPIE, 2017.
  * [49] Jean-Sebastien Brouillon, Emanuele Fabbiani, and Pulkit Nahata.  Identification tool.  https://github.com/donlelef/eiv-grid-id.
  * [50] S. Bolognani, R. Carli, G. Cavraro, and S. Zampieri.  On the need for communication for voltage regulation of power distribution grids.  IEEE Transactions on Control of Network Systems, 6(3):1111–1123, 2019.
  * [51] G. Hoogsteen, A. Molderink, J. L. Hurink, and G. J. M. Smit.  Generation of flexible domestic load profiles to evaluate demand side management approaches.  In 2016 IEEE International Energy Conference (ENERGYCON), pages 1–6, 2016.
  * [52] L. Thurner, A. Scheidler, F. Schafer, J. H. Menke, J. Dollichon, F. Meier, S. Meinecke, and M. Braun.  pandapower - an open source python tool for convenient modeling, analysis and optimization of electric power systems.  IEEE Transactions on Power Systems, 2018.