Comparison of Bounds for Optimal PMU Placement for State Estimation in Distribution Grids

TL;DR

This paper investigates optimal sensor placement in large, unbalanced distribution grids to enhance state estimation accuracy, comparing bounds and proposing methods for near-optimal solutions in complex network scenarios.

Contribution

It introduces and compares bounds for optimal PMU placement, leveraging properties like convexity and modularity, to efficiently approximate solutions in large-scale distribution grids.

Findings

Bounds effectively estimate optimal sensor placement performance.

Proposed methods improve state estimation accuracy in large networks.

Validated approach on IEEE benchmark feeders.

Abstract

The lack of measurements in distribution grids poses a severe challenge for their monitoring: since there may not be enough sensors to achieve numerical observability, load forecasts (pseudo-measurements) are typically used, and thus an accurate state estimation is not guaranteed. However, an estimation is required to control distribution grids given the increasing amount of distributed generation. Therefore, we consider the problem of optimal sensor placement to improve the state estimation accuracy in large-scale, 3-phase coupled, unbalanced distribution grids. This is a combinatorial optimization problem whose optimal solution is unpractical to obtain for large networks. We explore the properties of different metrics in the context of optimal experimental design, like convexity and modularity, to propose and compare several tight lower and upper bounds on the performance of the optimal solution. Moreover, we show how to use these bounds to choose near-optimal solutions. We test the method on two IEEE benchmark test feeders, the 123-bus and the 8500-node feeders, to show the effectiveness of the approach.

Equations

$$\begin{array}[]{c}\left[\begin{array}[]{c}I_{\text{src}}\\ I\end{array}\right]=Y\left[\begin{array}[]{c}V_{\text{src}}\\ V\end{array}\right],\;S=\text{diag}(\bar{I})V\end{array}$$

$$(S{)}_{\epsilon }=0,(I{)}_{\epsilon }=0$$

$$\begin{array}[]{c}z_{\text{meas}}\approx C_{\text{meas}}V+\text{diag}(C_{\text{meas}}V)(\omega_{\text{mag}}+j\omega_{\text{ang}})\end{array}$$

$$\begin{array}[]{l}(C_{\text{meas}}V)_{j}=(C_{\text{meas}})_{j,\bullet}V=\\[2.84544pt] \left\{\begin{array}[]{ll}V_{i_{l}}&\mbox{for a voltage measurement}\\[0.0pt] (Y)_{i_{l},\bullet}V&\mbox{for a current measurement}\\[0.0pt] (Y)_{i_{l},m_{l}}(V_{i_{l}}-V_{m_{l}})&\mbox{for a branch-current $i\to m$}\\[0.0pt] &\mbox{measurement}\end{array}\right.\end{array}$$

$$\begin{array}[]{rl}\Sigma_{\text{meas}}&=(\sigma_{\text{mag}}^{2}+\sigma_{\text{ang}}^{2})\text{diag}(\>\mathinner{\!\left\lvert C_{\text{meas}}V\right\rvert}^{2})\\ &\approx(\sigma_{\text{mag}}^{2}+\sigma_{\text{ang}}^{2})\text{diag}(\>\mathinner{\!\left\lvert z_{\text{meas}}\right\rvert}^{2})\end{array}$$

$$\begin{array}[]{l}V_{\text{prior}}=\text{PowerFlow}(S_{\text{psd}})\end{array}$$

$$\begin{array}[]{l}V_{\text{post}}=V_{\text{prior}}+K(z_{\text{meas}}-C_{\text{meas}}V_{\text{prior}})\end{array}$$

$$\begin{array}[]{rl}\Sigma_{\text{post}}=&\Sigma_{\text{prior}}+K(\Sigma_{\text{meas}}+C_{\text{meas}}\Sigma_{\text{prior}}C_{\text{meas}}^{*})K^{*}\\[2.84544pt] &-KC_{\text{meas}}\Sigma_{\text{prior}}-\Sigma_{\text{prior}}C_{\text{meas}}^{*}K^{*}\\[2.84544pt] K=&\arg\min_{K}\text{tr}(\Sigma_{\text{post}})\\[2.84544pt] =&\Sigma_{\text{prior}}C_{\text{meas}}^{*}(C_{\text{meas}}\Sigma_{\text{prior}}C_{\text{meas}}^{*}+\Sigma_{\text{meas}})^{-1}\end{array}$$

$$\begin{array}[]{c}\Sigma_{\text{post}}=F(\Sigma_{F,\text{prior}}^{-1}+(C_{\text{meas}}F)^{*}\Sigma_{\text{meas}}^{-1}C_{\text{meas}}F)^{-1}F^{*}\\[2.84544pt] \Sigma_{F,\text{post}}=(\Sigma_{F,\text{prior}}^{-1}+(C_{\text{meas}}F)^{*}\Sigma_{\text{meas}}^{-1}C_{\text{meas}}F)^{-1}\end{array}$$

$$\begin{array}[]{l}\Sigma_{F,\text{post}}=(\Sigma_{F,\text{prior}}^{-1}+\sum_{i}(C_{\text{meas}}F)_{i,\bullet}^{*}(C_{\text{meas}}F)_{i,\bullet}(\Sigma_{\text{meas}}^{-1})_{i,i})^{-1}\\[2.84544pt] =(\Sigma_{F,\text{prior}}^{-1}+\sum_{i}x_{i}(\tilde{C}_{\text{meas}}F)_{i,\bullet}^{*}(\tilde{C}_{\text{meas}}F)_{i,\bullet}(\Sigma_{\text{meas}}^{-1})_{i,i})^{-1}\\[2.84544pt] =\Sigma_{F,\text{post}}(x)\end{array}$$

$$\underset{x}{\min }m({\Sigma }_{F,\text{post}}(x))\text{ s.t. }h(x)\le 0,x_i\in \{0,1\}\forall i$$

$$\begin{array}[]{l}\textnormal{A-}optimal\mathrel{\mathop{\ordinarycolon}}f_{\textnormal{A}}(x)\mathrel{\mathop{\ordinarycolon}}=\textnormal{tr}(\Sigma_{\textnormal{post}}(x))=\textnormal{tr}(\Sigma_{F,\textnormal{post}}(x))\\[2.84544pt] \textnormal{D-}optimal\mathrel{\mathop{\ordinarycolon}}f_{\textnormal{D}}(x)\mathrel{\mathop{\ordinarycolon}}=\log(\textnormal{det}(\Sigma_{F,\textnormal{post}}(x)))\end{array}$$

$$\begin{array}[]{rl}(\nabla f_{\text{A}}(x))_{i}=&-\text{tr}\Big{(}\Sigma_{F,\text{post}}^{2}(x)(\tilde{C}_{\text{meas}})_{i,\bullet}^{*}(\tilde{C}_{\text{meas}})_{i,\bullet}\Big{)}(\Sigma_{\text{meas}}^{-1})_{i,i}\\[2.84544pt] (\nabla f_{\text{D}}(x))_{i}=&-\text{tr}\Big{(}\Sigma_{F,\text{post}}(x)(\tilde{C}_{\text{meas}})_{i,\bullet}^{*}(\tilde{C}_{\text{meas}})_{i,\bullet}\Big{)}(\Sigma_{\text{meas}}^{-1})_{i,i}\end{array}$$

$${\Lambda }_X:=\sum _{i\in X}({\tilde{C}}_{\text{meas}}{)}_{i,\bullet }^{\ast }({\tilde{C}}_{\text{meas}}{)}_{i,\bullet }({\Sigma }_{\text{meas}}^{-1}{)}_{i,i}⪰0$$

$$\begin{array}[]{ll}\text{submodular:}&\hat{f}(Y\cup\{a\})-\hat{f}(Y)\leq\hat{f}(X\cup\{a\})-\hat{f}(X)\\[2.84544pt] \text{supermodular:}&\hat{f}(Y\cup\{a\})-\hat{f}(Y)\geq\hat{f}(X\cup\{a\})-\hat{f}(X)\end{array}$$

$$x_{\text{opt}}=\arg \underset{x}{\min }f(x)\text{ s.t. }\sum _i{x}_i\le N_{\text{meas}},x_i\in \{0,1\}\forall i$$

$$x_{\text{convex}}=\arg \underset{x}{\min }f(x)\text{ s.t. }\sum _i{x}_i\le N_{\text{meas}},x_i\in [0,1]\forall i$$

$$x_{\text{feas},i}=\left\{\begin{array}[]{l}1\mbox{ if }x_{\text{convex},i}\geq x_{\text{convex},k}\\[1.42271pt] 0\mbox{ otherwise}\end{array}\right.$$

$$k=\arg \underset{i\notin X^{(K-1)}}{\min }\hat{f}(X^{(K-1)}\cup \{i\}),X^{(K)}=X^{(K-1)}\cup \{k\}$$

$$f_{\text{convex}}\le f_{\text{opt}}\le \min (f_{\text{greedy}},f_{\text{feas}})$$

$$\begin{array}[]{c}\hat{f}_{\text{D},\text{greedy}}-\hat{f}_{\text{D}}(\emptyset)\leq(\hat{f}_{\text{D},\text{opt}}-\hat{f}_{\text{D}}(\emptyset))\alpha\\[2.84544pt] \alpha=\left(1-\left(1-\frac{1}{N_{\text{meas}}}\right)^{N_{\text{meas}}}\right)\in(1-e^{-1},1]\end{array}$$

$${\tilde{f}}_{\text{D},\text{greedy}}:=({\hat{f}}_{\text{D},\text{greedy}}-{\hat{f}}_{\text{D}}(\varnothing )){\alpha }^{-1}+{\hat{f}}_{\text{D}}(\varnothing )\le f_{\text{D},\text{opt}}$$

$$x_{\text{opt}}=\arg \underset{x}{\min }f(x)\text{ s.t. }\sum _i{c}_i{x}_i\le b,x_i\in \{0,1\}\forall i$$

$$x_{\text{convex}}=\arg \underset{x}{\min }f(x)s.t.\sum _i{c}_i{x}_i\le b,x_i\in [0,1]\forall i$$

$$k={\arg \max }_{i\in \mathcal{B},i\notin X^{(K-1)}}{x}_{\text{convex},i},X^{(K)}=X^{(K-1)}\cup \{k\}$$

$$f_{\text{convex}}\le f_{\text{opt}}\le \min (f_{\text{greedy2}},f_{\text{feas}})$$

$$\begin{array}[]{c}\tilde{f}_{\text{D},\text{greedy1}}\coloneqq(f_{\text{D},\text{greedy1}}-\hat{f}_{\text{D}}(\emptyset))\beta^{-1}+\hat{f}_{\text{D}}(\emptyset)\leq f_{\text{D},\text{opt}}\\[2.84544pt] \tilde{f}_{\text{D},\text{greedy1a}}\coloneqq(f_{\text{D},\text{greedy1a}}-\hat{f}_{\text{D}}(\emptyset))\beta_{a}^{-1}+\hat{f}_{\text{D}}(\emptyset)\leq f_{\text{D},\text{opt}}\\[2.84544pt] \end{array}$$

$$\begin{array}[]{c}\beta=\left(1-\prod_{i\in X_{\text{greedy1}}}\left(1-\frac{c_{i}}{b}\right)\right)\in(0,1]\\[2.84544pt] \beta_{a}=\left(1-\prod_{i\in X_{\text{greedy1a}}}\left(1-\frac{c_{i}}{b}\right)\right)\in(1-e^{-1},1]\end{array}$$

$$\begin{array}[]{c}\beta\in(1-e^{-\gamma},1],\beta_{a}\in(1-\gamma e^{-\gamma},1]\end{array}$$

$$\underset{\mathcal{A}}{\max }{\tilde{f}}_{\text{D},\text{online}}(\mathcal{A})\le f_{\text{D},\text{opt}}$$

Full text

I Introduction

The operation of a power network requires accurate monitoring of its state: bus voltages, line currents, consumption and generation, to efficiently manage its controllable elements. State Estimation (SE) serves that purpose by estimating a minimum representation of the network state, for example the bus voltage phasors. Taking several measurements and the admittance matrix of the network as parameters, SE typically solves a weighted least-squares problem using an iterative approach like Newton-Raphson [1, 2]. In transmission networks, SE is fundamental, since the volatile and distributed generation, which injects power in different locations of the network, causes fast-changing bidirectional power flows. In contrast, SE has not been so necessary in distribution grids until recently, since these grids used to have a simple radial structure with a single source bus injecting power. But this situation is changing due to the increasing penetration of distributed generation like PV panels, batteries, etc. [3].

The introduction of Phasor Measurement Units (PMUs) improves the monitoring of electrical networks. Some work in the literature is focused on using PMUs to achieve topological observability [4], which can be solved using integer linear programming [5]. Although cheap PMUs [6] are becoming available for mass deployment, their operational and network communication costs [7] may still prevent installing the required minimum number of sensors to achieve topological observability based only on PMUs, especially in distribution grids. If topological observability is not possible, neither is numerical observability ensured, which is required to solve the SE [4]. As a result, conventional Supervisory Control And Data Adquisition (SCADA) measurements need to be combined with PMUs to solve the SE problem [8]. Therefore, authors of recent work [9, 10] formulate the problem of optimal PMU placement in terms of maximising the SE accuracy. This accuracy is typically measured through some metric in the context of optimal design of experiments [11], applied to the covariance matrix of the SE error caused by noisy measurements and unreliable forecasts. As a result, the optimal PMU placement problem is a combinatorial optimization problem with a nonlinear objective function, the SE accuracy metric, and thus cannot be solved using linear programming. Moreover, as the size of the network increases, the optimal solution becomes unpractical to obtain given the number of possible combinations of measurements [9].

In distribution grids, one of the major limitations for SE is the lack of sufficient real-time measurements to achieve observability, even using SCADA, and thus SE algorithms need to rely on pseudo-measurements, such as load forecasts. These pseudo-measurements have typically a large relative noise associated (approximately $50$% [12]), which causes a high uncertainty in the SE estimates [13]. As a result, there is a growing interest in using PMUs in distribution grids [14], and some recent work proposes algorithms to place sensors to satisfy a desired performance, like algorithms using greedy and random combinations of sensors [15, 16] or evolutionary algorithms [17, 18]. However, these approaches have no optimality guarantees. Since optimal solutions are unpractical in large networks, [19] uses the results in [9] to derive lower bounds for the values of the optimal solution.

Our contribution consists in proposing and comparing several bounds for the optimal solution of the PMU placement problem in distribution grids, and extending these bounds to large-scale grids. With these bounds, we will be able to check the gap between any given suboptimal solution and the optimal solution. First, we prove properties like convexity and supermodularity for some metrics in [11]; then, we use these properties to derive a combination of bounds based on convex optimization as in [10], and bounds based on supermodular minimization similar to [9]. These bounds allow us to obtain a close bound for the performance of the optimal solution of the problem under a cardinality constraint, as well as under a budget constraint and sensors with different costs. Additionally, we show how these bounds scale to large networks by using a projected gradient descent algorithm. Moreover, we show through two examples of test feeders, one of them large-scale, how bounds based on the supermodularity property perform better than those based on convexity if the number of sensors to deploy is small, as in distribution networks.

The rest of the paper is structured as follows. Section II contains the nomenclature for the most relevant symbols. Section III-A presents some background about power networks. Section III-B discusses the different types of measurements. Section III-C summarizes the newly proposed methodology for SE [13]. Section IV presents the metrics considered for the problem and their properties. Section V states the optimal sensor placement problem and derives the lower and higher bounds for each optimum value. Section VI shows the effectiveness of the bounds on a test case. Finally, Section VII presents the conclusions. For clarity purposes, proofs can be found in the Appendices A and B. Appendices C and D describe the methods to extend the approach to large grids.

II Nomenclature

[TABLE]

[TABLE]

III State Estimation in Distribution Grids

III-A Distribution Grid Model

A distribution grid consists of buses, where power is injected or consumed, and branches, each connecting two buses. This system can be modeled as a graph $\mathcal{G}=(\mathcal{V},\mathcal{E},\mathcal{W})$ with nodes $\mathcal{V}=\{1,...,N_{\text{bus}}\}$ representing the buses, edges $\mathcal{E}=\{(v_{i},v_{j})\mid v_{i},v_{j}\in\mathcal{V}\}$ representing the branches, and edge weights $\mathcal{W}=\{w_{i,j}\mid(v_{i},v_{j})\in\mathcal{E}\}$ representing the admittance of the branches, which are determined by the length and type of the line cables. In 3-phase networks buses may have up to 3 phases, so that the voltage at bus $i$ is $V_{\text{bus},i}\in\mathbb{C}^{n_{\phi,i}}$, where $n_{\phi,i}\leq 3$ (and the edge weights $w_{i,j}\in\mathbb{C}^{n_{\phi,i}\times n_{\phi,j}}$). The state of the network is then typically represented by the vector bus voltages $V_{\text{bus}}=[V_{\text{src}}^{T},\;V^{T}]^{T}\in\mathbb{C}^{N+3}$, where $V_{\text{src}}\in\mathbb{C}^{3}$ denotes the known voltage of the 3 phases at the source bus, and $V\in\mathbb{C}^{N}$ the voltages in the non-source buses, where $N$ depends on the number of buses and phases per bus. Then, using the Laplacian matrix $Y\in\mathbb{C}^{(N+3)\times(N+3)}$ of the weighted graph $\mathcal{G}$, called admittance matrix [1], the power flow equations to compute the currents $I$ and the power loads $S$ are:

[TABLE]

III-B Measurements

As explained in [13], several different sources of information can be available to solve the SE problem:

Pseudo-measurements, i.e., load estimations $S_{\text{psd}}$ based on predictions and/or known installed load capacity at every bus. Since these pseudo-measurements are estimations rather than actual measurements, we model their uncertainty using a Gaussian noise with a relative large standard deviation (a typical value can be $\sigma_{\text{psd}}\approx 50\%$ [12]). 2. 2.

Virtual measurements, i.e., buses with zero-injections, no loads connected. They can be represented as physical constraints for the voltage states by defining the set of indices of zero-injection buses $\varepsilon=\{i,\cdots,j\}$:

[TABLE] 3. 3.

Real-time PMU measurements, i.e., voltage and current GPS-synchronized measurements of magnitude and phase angle. According to the IEEE standard for PMU [20], these measurements may have a small error. Again, we model this uncertainty using a Gaussian noise with a low standard deviation for the magnitude and the angle, $\sigma_{\text{mag}}\approx 1\%$ and $\sigma_{\text{ang}}\approx 0.01\text{ rad}$ respectively. They can be expressed using a linear approximation [13] with magnitude and angle noise due to the measurements and imperfect synchronization. For a number of $N_{\text{meas}}$ measurements $z_{\text{meas}}\in\mathbb{C}^{N_{\text{meas}}}$ we have

[TABLE]

where $\omega_{\text{mag}},\omega_{\text{ang}}$ are the Gaussian noises with mean [math] and standard deviation $\sigma_{\text{mag}},\sigma_{\text{ang}}$ respectively: $\omega_{\text{mag}}\sim\mathcal{N}(0,\sigma_{\text{mag}}I_{\text{d},N_{\text{meas}}})$, $\omega_{\text{ang}}\sim\mathcal{N}(0,\sigma_{\text{ang}}I_{\text{d},N_{\text{meas}}})$, with $I_{\text{d},n}$ denoting the identity matrix of dimension $n$, and where the matrix $C_{\text{meas}}$ maps state voltages to measurements; it relates the values of the measurements to values of state voltages. Then, for measurement $j$ at phase $l$ of bus $i$ we have

[TABLE]

Since the measurement noises in (3) are small according to the PMUs standard [20], their covariance matrices can be approximated using the measurements:

[TABLE]

III-C State Estimation

Typically, SE consists in finding the voltages that best match the measurements by solving a weighted least-squares problem [1]. As proposed in [13], SE can be decomposed in two parts: First, using the pseudo-measurement estimations or predictions for the loads $S_{\text{psd}}$, we solve the power flow offline to obtain a prior estimate $V_{\text{prior}}$:

[TABLE]

Then, using the real-time PMU measurements $z_{\text{meas}}$, a posterior solution $V_{\text{post}}$ can be derived using a linear filter:

[TABLE]

where the gain matrix $K$ is obtained by minimizing the error covariance $\Sigma_{\text{post}}=\mathbb{E}[(V_{\text{post}}-V)^{*}(V_{\text{post}}-V)]$:

[TABLE]

where $\Sigma_{\text{prior}}$ and $\Sigma_{\text{meas}}$ are the expected error covariance of the prior estimate $V_{\text{prior}}$ and the measurements $z_{\text{meas}}$ respectively.

Remark 1**.**

Extra non-synchronized real-time measurements, like magnitude measurements from SCADA, can be included in the posterior update (6) for a greater improvement of the posterior estimate $V_{\text{post}}$ by using the first-order approximation of the measurement function [13]. If a multi-stage PMU deployment is considered [21], prior installed PMUs can be also considered as extra measurements in the same manner.

Remark 2**.**

As shown in [8], splitting the problem in two steps yields the same first-order approximation as solving the problem in one step. Moreover, the posterior minimum-variance estimator using the linear update (6), is equal to the maximum-likelihood using a weighted least-squares approach [13]. Therefore, we can conclude that for an SE method that assumes Gaussian noises and performs a maximum likelihood estimation, the posterior covariance will be approximately the one in (7), and thus the method developed here for optimal sensor placement can be also extended for other SE techniques satisfying these conditions.

Since $V_{\text{post}}$ in (6) is an unbiased estimator, SE accuracy can be defined as a function of the posterior covariance matrix $\Sigma_{\text{post}}$ in (7), which needs to be minimized to improve the SE accuracy. This motivates a deeper analysis of $\Sigma_{\text{post}}$: if there are zero-injection buses at indices $\varepsilon$ when solving (5), then $\Sigma_{\text{prior}}$ will not be of full rank and thus will not be invertible. For convenience, we consider the restricted subspace $\{V\mid(Y)_{\varepsilon,\bullet}V=0\}$ and the linear transformation proposed in [13] to represent the space of feasible solutions: $V=Fx+V_{0}$ with $x\in\mathbb{C}^{N-\mathinner{\!\left\lvert\varepsilon\right\rvert}}$, where $\mathinner{\!\left\lvert\varepsilon\right\rvert}$ is the cardinality of $\varepsilon$ and $F$ is the null space of $(Y)_{\varepsilon,\bullet}$: $F=\text{ker}((Y)_{\varepsilon,\bullet})\in\mathbb{C}^{N\times N-\mathinner{\!\left\lvert\varepsilon\right\rvert}}$, so that $F^{*}F=I_{\text{d}}$, and $V_{0}$ denotes the voltage under zero loads. Then we have $\Sigma_{\text{prior}}=F\Sigma_{F,\text{prior}}F^{*}$, where $\Sigma_{F,\text{prior}}$ is the covariance of $x$. After some manipulations of (7), the resulting error covariance $\Sigma_{\text{post}}$ ($\Sigma_{F,\text{post}}$ in this subspace) for the posterior estimation $V_{\text{post}}$ can be expressed as

[TABLE]

In the presence of zero-injection buses, $\mathinner{\!\left\lvert\varepsilon\right\rvert}>0$, and then $\Sigma_{\text{post}}\in\mathbb{C}^{N\times N}$ is not full rank: $\text{rank}(\Sigma_{\text{post}})=N-\mathinner{\!\left\lvert\varepsilon\right\rvert}$. It has the same rank as $\Sigma_{F,\text{post}}\in\mathbb{C}^{(N-\mathinner{\!\left\lvert\varepsilon\right\rvert})\times(N-\mathinner{\!\left\lvert\varepsilon\right\rvert})}$, which is full rank since $\Sigma_{F,\text{prior}}$ is full rank. Moreover, the eigenvalues of $\Sigma_{F,\text{post}}$ are all eigenvalues of $\Sigma_{\text{post}}$, because for any eigenvector $v$ of $\Sigma_{F,\text{post}}$, $Fv$ is an eigenvector of $\Sigma_{\text{post}}$ with the same eigenvalue, since $F^{*}F=I_{\text{d}}$. Consequently, we can analyze $\Sigma_{F,\text{post}}$ instead of $\Sigma_{\text{post}}$. These eigenvalues represent the lengths of the axes of the confidence ellipsoid [22]. Moreover, since measurement errors are caused separately by each sensor, they are independent, i.e. $\Sigma_{\text{meas}}$ is diagonal, and we can split $\Sigma_{F,\text{post}}$ by every measurement:

[TABLE]

where $x_{i}\in\{0,1\}$, $x_{i}=1$ if the physical quantity $i$ has a sensor measuring its value, [math] otherwise; and $\tilde{C}_{\text{meas}}$ is the special case of $C_{\text{meas}}F$ with all possible measurements of all types (bus voltage, bus current, and line current) for all nodes and lines in each phase.

Definition 1**.**

In order to improve the accuracy of the SE, the problem of optimal sensor placement consists in minimizing $\Sigma_{\text{post}}(x)$, equivalently $\Sigma_{F,\text{post}}(x)$ in (9), according to a metric $m(\cdot)$ and under a set of constraints $h(\cdot)$ to limit the number of sensors or the total cost:

[TABLE]

For simplicity, we define: $f(x)\equiv m(\Sigma_{F,\text{post}}(x))$.

Remark 3**.**

Problem (10) is equivalent to minimizing the number of sensors or the total cost of the sensors, while enforcing a given accuracy, i.e., satisfying a performance threshold $\lambda$ of the metric: $\min_{\{x\mid f(x)\leq\lambda,\;x_{i}\in\{0,1\}\;\forall i\}}h(x)$, as in [15, 16, 17, 18]. To achieve that, the problem (10) can be solved for different numbers of sensors or budgets until the performance value is below the threshold: $\min_{\{x\mid h(x)\leq 0,\;x_{i}\in\{0,1\}\;\forall i\}}f(x)\leq\lambda$. The representation in (10), as used in [9, 10, 19] is preferred due to the availability of a computationally efficient optimization algorithm to solve the problem, see Appendix C.

IV Metrics for Sensor Placement

There are many possible metrics $f(x)$ available in the context of optimal design of experiments [11]. Concretely, we will focus on the A-optimal and the D-optimal metrics, because of the properties that we will show later. Nonetheless, the results of this paper could be generalized to other metrics with the same properties.

Definition 2**.**

The A-optimal and the D-optimal metrics are defined as follows:

[TABLE]

The A-optimal metric represents the sum of the eigenvalues of $\Sigma_{F,\text{post}}$ and thus the sum of the lengths of axes of the confidence ellipsoid [22]. This is the metric typically used for SE methods, since standard SE maximizes the log-likelihood through a weighted least-squares minimization [1], which is equivalent to the minimum-variance estimator using the trace [13]. The D-optimal metric is the natural logarithm of the product of the eigenvalues of $\Sigma_{F,\text{post}}$ and it is related to the logarithm of the volume of the confidence ellipsoid [22]. These metrics have several properties relevant to the problem of optimal sensor placement:

Convexity: When relaxing $x_{i}$ in (10) to be continuous, $x_{i}\in[0,1]$, the metrics are convex on $x$ [22, Section 7.5] [10], and thus a global optimum can be computed efficiently. Note that despite $\log(x)$ being concave for $x\in\mathbb{R}$, $\log(\text{det}(X))$ is convex for $X\in\mathbb{R}^{n\times n}$ [22]. 2. 2.

Gradient computation: The gradients $\nabla f$ for $f\in\{f_{\text{A}},f_{\text{D}}\}$ can be derived analytically using matrix calculus [23]:

[TABLE]

These expressions will be necessary when developing gradient methods to optimize the placement of sensors in large networks, see Appendix C. 3. 3.

Monotonicity: Moreover, these metrics can be seen as set functions by defining the sets $X=\{i\mid x_{i}=1\}$. Note that there is a bijection between both, i.e., $i\in X\iff x_{i}=1$. Therefore, we will use the notation $X$ and $x$ to denote the set and the vector respectively, and $\hat{f}(X)=f(x)$ to denote the metric functions applied to them. Defining

[TABLE]

the posterior error covariance can be expressed as $\Sigma_{F,\text{post}}(X)=(\Sigma_{F,\text{prior}}^{-1}+\Lambda_{X})^{-1}$.

Proposition 1**.**

$\hat{f}_{\mathrm{A}}(X),\hat{f}_{\mathrm{D}}(X)$* are monotone decreasing.*

Proof.

Consider the sets of sensors $X\subseteq Y$. Since $\Lambda_{Y}-\Lambda_{X}=\Lambda_{Y\setminus X}\succeq 0$, we have $\Lambda_{Y}\succeq\Lambda_{X}$. Then $\Sigma_{F,\text{prior}}^{-1}+\Lambda_{Y}\succeq\Sigma_{F,\text{prior}}^{-1}+\Lambda_{X}\succeq 0$ and thus $(\Sigma_{F,\text{prior}}^{-1}+\Lambda_{X})^{-1}\succeq(\Sigma_{F,\text{prior}}^{-1}+\Lambda_{Y})^{-1}$. Consequently, for any of the functions $\hat{f}\in\{\hat{f}_{\text{A}},\hat{f}_{\text{D}}\}$, $\hat{f}(X)\geq\hat{f}(Y)$. ∎ 4. 4.

Modularity: Submodularity/supermodularity are a sort of concavity/convexity properties when considering set functions. The increment of a set function when adding a new element diminishes (submodular) or increases (supermodular) as the set gets larger.

Definition 3**.**

Given a finite set $\Omega$, a submodular/supermodular function [24] is a set function $\hat{f}\mathrel{\mathop{\ordinarycolon}}2^{\Omega}\to\mathbb{R}$, with $\hat{f}(\emptyset)=0$, where $\emptyset$ is the empty set, so that for any element $a\in\Omega$ and two subsets $X,Y$ so that $X\subseteq Y\subseteq\Omega\setminus\{a\}$ we have

[TABLE]

Proposition 2**.**

$\hat{f}_{\mathrm{D}}(X)-\hat{f}_{\mathrm{D}}(\emptyset)$* is supermodular.*

Proof.

See Appendix A ∎

This means that for $f_{\text{D}}$, there exist bounds between greedy solutions and optimal solutions, as we will see later in (19) and (24). These bounds will allow us to provide a limit on the values of the optimal solution, and thus to check how far the values of other suboptimal solutions can be from the value of the optimal solution.

Remark 4**.**

In [9] was proven that if the rows in $\tilde{C}_{\text{meas}}$ are orthogonal, $\hat{f}_{\mathrm{A}}(X)-\hat{f}_{\mathrm{A}}(\emptyset)$ is supermodular, and lower bounds based on the supermodularity property can also be used for the A-optimal metric. However, the rows of $\tilde{C}_{\text{meas}}$ are in general not orthogonal due to the different types of sensors in (4), the electrical connections represented in $Y$, and the effect of the reduced subspace $F$.

V Optimal Sensor Placement

So far we have looked at the cost function in (10), now we will focus on the constraints. In this section we state the optimal sensor placement problem under two different constraints: a cardinality constraint limiting the number of sensors, which is the most typical approach in the literature [9, 10, 12], and a budget constraint limiting the total cost of the deployed sensors, which may have different costs. For each case, we derive the respective lower and upper bounds on the optimal performance, using the metrics properties derived in Section IV.

V-A Solutions with cardinality constraint

For any of the metrics $f\in\{f_{\text{A}},f_{\text{D}}\}$, the problem of optimal placement of $N_{\text{meas}}$ sensors is

[TABLE]

where $N_{\text{meas}}$ indicated the maximum number of sensors. We denote the value of the optimal solution of (13) $f_{\text{opt}}=f(x_{\text{opt}})$.

Similar to [25, 26] for general estimation problems and in [10] for transmission power networks, we relax the constraints in (13), $x_{i}\in[0,1]$, to get a continuous convex problem using

[TABLE]

We denote the value of the optimal solution of (14) as $f(x_{\text{convex}})=f_{\text{convex}}$. However, this solution will not necessarily be feasible to (13). To create a feasible solution $x_{\text{feas}}$ (with respective value $f(x_{\text{feas}})=f_{\text{feas}}$), we can take the largest $N_{\text{meas}}$ values of $x_{\text{convex}}$, set them to $1$ and the others to [math]:

[TABLE]

where $k$ is such that $\mathinner{\!\left\lvert\{i\mid x_{\text{convex},i}\geq x_{\text{convex},k}\}\right\rvert}=N_{\text{meas}}$. Ties are broken arbitrarily if there is more than one $x_{\text{convex},i}$ with value $x_{\text{convex},k}$, i.e. if $\mathinner{\!\left\lvert\{i\mid x_{\text{convex},i}=x_{\text{convex},k}\}\right\rvert}>1$. This also applies for further possible ties throughout the paper. Another simple way to create a feasible solution would be using a forward greedy sensor selection: at iteration $K$, given the set of selected sensors $X^{(K-1)}$, add a new sensor $k$ such that

[TABLE]

We denote the solution of (16) $x_{\text{greedy}}$ and its value $f(x_{\text{greedy}})=f_{\text{greedy}}$. Then, since $x_{\text{opt}}$ is a feasible suboptimal solution of (14), and $x_{\text{feas}}$ and $x_{\text{greedy}}$ are feasible suboptimal solutions of (13), the following holds for all metrics:

[TABLE]

For the case of $\hat{f}_{\text{D}}$, since $\hat{f}_{\text{D}}(X)-\hat{f}_{\text{D}}(\emptyset)$ is monotone nonincreasing and supermodular ($-\hat{f}_{\text{D}}(X)+\hat{f}_{\text{D}}(\emptyset)$ is monotone nondecreasing submodular), we have an extra lower bound [24] for $f_{\text{opt}}$ of (13):

[TABLE]

so that

[TABLE]

V-B Solutions with budget constraint

Although the problem of sensor placement under a cardinality constraint is simpler and has more desired properties [24], a more realistic and economic representation of the problem would be having a budget constraint on the sensors

[TABLE]

where $b$ represents the budget and $c_{i}$ the cost of installing a sensor at location $i$. This can take into account the extra cost of installing a sensor in a remote area, or locations where specific rights might be required, or different types of sensors, etc. The relaxed convex problem is then

[TABLE]

Similar to (15), a feasible solution $x_{\text{feas}}$ w.r.t. (20) can be built using the convex solution $x_{\text{convex}}$ of (21), by iteratively taking the sensor with highest $x_{\text{convex},i}$: at iteration $K$, and until the set of possible sensors is empty, $\mathcal{B}=\{i\mid c_{i}\leq b-\sum_{j\in X^{(K-1)}}c_{j}\}=\emptyset$, add to the set of selected sensors $X^{(K-1)}$ a new sensor $k\notin X^{(K-1)}$ such that:

[TABLE]

Likewise, we can adapt the forward greedy selection algorithm (16) to take costs into account. Therefore, we can use the the cost-effective forward greedy selection algorithm of [27], which in every iteration adds the sensor with the lowest ratio of objective function improvement divided by sensor cost. In this version we also return intermediate values required for the bounds. For clarity we detail the whole algorithm here:

Using Algorithm 1, we denote the solutions as $X_{\text{greedy1}}=X_{1}$, $X_{\text{greedy2}}=X_{2}$ and $X_{\text{greedy1a}}=X_{1}\cup\{a\}$ with corresponding values $\hat{f}(X_{\text{greedy1}})=f_{\text{greedy1}},\hat{f}(X_{\text{greedy2}})=f_{\text{greedy2}},\hat{f}(X_{\text{greedy1a}})=f_{\text{greedy1a}}$. Since $X_{\text{greedy1}}\subseteq X_{\text{greedy2}}$, $f_{\text{greedy2}}\leq f_{\text{greedy1}}$. Note that in the cardinality constrained case we would have $X_{1}=X_{2}$, and hence $f_{\text{greedy1}}=f_{\text{greedy2}}$; therefore we use $f_{\text{greedy}}$ in that case. Then, as in (17) we know:

[TABLE]

Again, for the case of $\hat{f}_{\text{D}}$, since $\hat{f}_{\text{D}}(X)-\hat{f}_{\text{D}}(\emptyset)$ is nonincreasing supermodular, we have two extra lower bounds for the value $f_{\text{opt}}$ of (20), which are derived from the proofs in [28, 29]:

[TABLE]

with

[TABLE]

Remark 5**.**

If the costs $c_{i}$ are unitary and the budget $b$ equals the number of sensors, i.e $c_{i}=1\;\forall i$ and $b=N_{\text{meas}}$, then $\alpha=\beta$. Consequently, (19) is a particular case of (24).

A drawback of the bound with $\beta$ in (24) is that $\beta$ does not have a lower bound bigger than [math], and thus it may happen that $\beta^{-1}\to\infty$ and the bound becomes trivial. Nonetheless, we can prove the following result:

Proposition 3**.**

Let $\gamma\in(0,1]$ define the percentage of budget used, so that $\sum_{i\in X_{\text{greedy1}}}c_{i}=\gamma b\leq b$, then we have

[TABLE]

Proof.

See Appendix B ∎

Proposition 3 shows that if $\gamma\to 0$, then $\beta\to 0$, but $\beta_{a}\to 1$; so only the bound with $\beta_{a}$ is useful. If $\gamma\to 1$, $\beta\to\beta_{a}$, and the bound with $\beta$ may be better.

Given the supermodularity of $f_{\text{D}}$, see Proposition 2, we can also consider the bound proposed in [27], called online bound, generated by Algorithm 2 to obtain

[TABLE]

Finally we have all bounds:

[TABLE]

Remark 6**.**

Note that since $\mathcal{A}$ may be any set of sensors, $\max_{\mathcal{A}}\tilde{f}_{\mathrm{D},\text{online}}(\mathcal{A})$ is computationally demanding.

VI Test Case

Now we test the bounds on the 123-bus [30] and the 8500-node [31] test feeders. For the budget-constrained case in the 123-bus feeder, we consider two more expensive zones to represent heterogeneous costs: $\tilde{c}_{i}=2$ in the top left blue zone, $\tilde{c}_{i}=1.5$ in the bottom right red zone, $\tilde{c}_{i}=1$ elsewhere, see Fig. 1. Then, sensor costs are normalized $c_{i}=\frac{\tilde{c}_{i}}{\sum_{i}\tilde{c}_{i}}$, so that their average is $1$, and thus the budget is approximately the number of sensors deployed. In the 8500-node feeder, we assign random normal distributed costs: $c_{i}\sim\mathcal{N}(1,0.1)$. For any number of measurements or budget, the sensor locations can be recovered from each solution $x$.

The algorithms are coded in Python and run on an Intel Core i7-6700HQ CPU at 2.60GHz with 16GB of RAM. For the convex optimization problems (14) and (21), we use the function minimize of the scipy.optimize package for the 123-bus test case, and a projected gradient descent algorithm for the 8500-node, see Appendix C, where we use an efficient projection algorithm as in [32] in order to handle the computational complexity of the 8500-node feeder. Additionally, since function evaluations may take a few seconds for the 8500-node feeder, we use matrix algebra results to speed up the computation of the greedy solutions in (16) and Algorithm 1, see Appendix D. Otherwise, it could take several days to evaluate all possible nodes in every step.

Fig. 2 shows the bounds for the 123-bus and the 8500-node test feeders. Both the A,D-optimal metrics are analyzed, under a cardinality and a budget constraint, for different numbers of sensors and budgets respectively. For simplicity, in the online bound $\tilde{f}_{\text{D},\text{online}}(\mathcal{A})$ in (27), we have only considered $\mathcal{A}=\emptyset$ instead of any $\mathcal{A}$. However, we have observed that other options, like using a greedy solution for $\mathcal{A}$, produce similar bounds. The online bound has also been applied to the cardinality constrained problem. The yellow shaded area with squares shows the area between the minimum upper bound and the maximum lower bound, and thus the possible locations of the optimal solution $f_{\{\text{A,D}\},\text{opt}}$.

For the D-optimal metric, in Fig. 2c, 2d, 2g and 2h, it can be observed that the online bound $\tilde{f}_{\text{D},\text{online}}$ (27) outperforms the rest of the lower bounds ($\tilde{f}_{\text{D},\text{greedy}}$ (19), $\tilde{f}_{\text{D},\text{greedy1}}$ (24), $f_{\text{D},\text{convex}}$ (14), (21)) for a small number of sensors under a cardinality and under a budget constraint, especially in the 8500-node feeder. For clarity, $\tilde{f}_{\text{D},\text{greedy1a}}$ (24) has not been plotted, since the results were almost equal to $\tilde{f}_{\text{D},\text{greedy1}}$. For a small number of sensors, it is remarkable how the lower bound based on convex relaxations $f_{\text{D},\text{convex}}$ performs poorly compared to the other lower bounds, as it produces too optimistic results. This is a result of the high SE uncertainty in all nodes of the grid before the PMU placement; as a consequence, the convex optimization allocates partial unfeasible PMUs with $0<x_{i}\ll 1$ in many nodes of the grid, instead of full PMUs with $x_{i}=1$ in a few nodes. This is equivalent to adding more sensors with an error with larger relative standard deviation. This cannot be done, since PMUs have a fixed maximum relative error by design [20], see section III-B.

For the A-optimal metric, in Fig. 2a, 2b, 2e and 2f we do not have the lower bounds based on the supermodularity property; however, the lower bound based on convex relaxations $f_{\text{A},\text{convex}}$ (14), (21) performs better than for the D-optimal metric. It still can be observed that for a small number of sensors this convex bound is far from the greedy and feasible solutions ($f_{\text{A},\text{feas}}$ (15), $f_{\text{A},\text{greedy}}$ (16), $f_{\text{A},\text{greedy2}}$ in Algorithm (1)), but not as much as in the D-optimal metric case, and the bound approaches the values of the solutions quicker.

As expected, for both metrics simple greedy or feasible solutions ($f_{\{\text{A,D}\},\text{greedy}}$, $f_{\{\text{A,D}\},\text{greedy2}}$, $f_{\{\text{A,D}\},\text{feas}}$) perform better than random configurations ($f_{\{\text{A,D}\},\text{rand}}$), even better than the best random configuration out of the 100 samples displayed.

VII Conclusions

We have stated the problem of optimal sensor placement for deploying PMUs to minimize the uncertainty of state estimation in distribution grids, both under a cardinality and a budget constraint. We have analyzed the properties of different metrics, concretely, convexity and supermodularity. Using these properties, we have derived a set of bounds that have enabled us to narrow the gap of the possible value of the optimal solution of this intractable problem. Since the optimal solution is unpractical to obtain in large grids, suboptimal solutions are required instead. Using these bounds, we can compute the maximum gap between any suboptimal solution and the optimal solutions. Moreover, we have observed, that the bounds produced by the supermodularity property are especially relevant when only a small number of sensor can be installed in large networks. In this case, a bound based on convex relaxation produces a too optimistic result.

Future work could include extending these results to take into account network reconfiguration due to different states of switches. Moreover, more exhaustive search algorithms could be developed to obtain solutions closer to the lower bounds. Additionally, it would be interesting to analyze how the choice of the metric affects the effectiveness of the bounds.

Appendix A Proof of Proposition 2

Proof.

The constant term $\hat{f}_{\text{D}}(\emptyset)$ is only necessary to ensure that the function is [math] when $X=\emptyset$. For the rest of the proof this term is not necessary since for any two sets $X,Y$ we have: $(\hat{f}_{\text{D}}(Y)-\hat{f}_{\text{D}}(\emptyset))-(\hat{f}_{\text{D}}(X)-\hat{f}_{\text{D}}(\emptyset))=\hat{f}_{\text{D}}(Y)-\hat{f}_{\text{D}}(X)$.

This proof is very similar to the one used for controlability gramians in [33]. Consider the set of sensors $X\subseteq Y$ and a sensor $a$ so that $\{a\}\notin Y$. Now let us define

[TABLE]

So that we have $\tilde{f}_{\text{D}}(0,X,Y)=\hat{f}_{\text{D}}(X)$, $\tilde{f}_{\text{D}}(1,X,Y)=\hat{f}_{\text{D}}(Y)$, $\tilde{g}_{a}(0,X,Y)=g_{a}(X)$ and $\tilde{g}_{a}(1,X,Y)=g_{a}(Y)$. Now we will prove that $g_{a}(Y)\geq g_{a}(X)$ by computing the gradient of $\tilde{g}_{a}$:

[TABLE]

Since $X\subseteq Y$, we have $\Lambda_{Y\setminus X}\succeq 0$. Since $\Lambda_{\{a\}}\succeq 0$ and $\Lambda_{X\cup\{a\}}=\Lambda_{\{a\}}+\Lambda_{X}$, we have

[TABLE]

and thus

[TABLE]

Consequently, $\frac{\partial\tilde{g}_{a}}{\partial\gamma}(\gamma,X,Y)\geq 0$ and

[TABLE]

so that

[TABLE]

∎

Appendix B Proof of Proposition 3

Proof.

If $\sum_{i\in X_{\text{greedy1}}}c_{i}=\gamma b$, then $\prod_{i\in X_{\text{greedy1}}}\left(1-\frac{c_{i}}{b}\right)$ achieves its maximum at $c_{i}=\frac{\gamma b}{\mathinner{\!\left\lvert X_{\text{greedy1}}\right\rvert}}$ for all $i$, where $|\cdot|$ is the number of elements in a set. Then we have:

[TABLE]

We assume that $c_{i}\leq b$ for all $i$ (sensors with $c_{i}>b$ are discarded since they cannot be installed). Then we know that $b>c_{a}>b-\gamma b$, where $c_{a}$ is the cost of element $a$, and we have:

[TABLE]

∎

Appendix C Efficient Projected Gradient Descent

To solve the convex problem under a cardinality constraint (14) for the 8500-node feeder, we use a projected gradient descent method using the gradient expressions in (12):

[TABLE]

where the suffix $(\cdot)^{(K)}$ denotes the value at iteration $K$. We use $\alpha^{(K)}=\frac{\alpha}{k\mathinner{\!\left\lVert\nabla f(x^{(K)})\right\rVert}}_{2}$ to guarantee convergence of the method [34], where $\alpha$ is a design parameter. For a more efficient projection $\Pi_{\mathcal{X}}(\cdot)$, we use the scaled boxed-simplex projection algorithm proposed in [32], which converges in a finite number of iterations, as opposed to the projection algorithm proposed in [10].

When solving the problem under a budget constraint (21), we can use the change of variables $y_{i}=x_{i}c_{i}$, and the function f_{c}(y)=f\Big{(}\frac{y}{c}\Big{)}, and solve the alternative problem

[TABLE]

using the projected gradient descent method:

[TABLE]

where an efficient implementation of the projection $\Pi_{\mathcal{Y}}(\cdot)$ can be derived by using a modified version of the one in [32] by changing the $1$ to $c_{i}$ for each respective $y_{i}$. The details of the algorithm can be found in [35].

Appendix D Efficient Greedy Computations

The posterior covariance at iteration $K$ can be defined as

[TABLE]

With this expression, we use Woodbury’s matrix identity [23] to obtain a simplified expression that is easier to evaluate:

[TABLE]

Then, we use Sylvester’s determinant identity [23] to get

[TABLE]