Power System State Estimation via Feasible Point Pursuit: Algorithms and Cramer-Rao Bound

TL;DR

This paper introduces feasible point pursuit algorithms for power system state estimation, providing improved solutions for nonconvex problems and establishing a Cramér-Rao bound as a performance benchmark.

Contribution

The paper proposes novel FPP-based solvers for power flow and PSSE, outperforming existing methods and converging to stationary points, with a new CRLB derivation for PSSE.

Findings

FPP-based solvers outperform traditional methods in accuracy.

The algorithms converge to stationary points of the WLS problem.

Numerical tests show superior performance on IEEE systems.

Abstract

Accurately monitoring the system's operating point is central to the reliable and economic operation of an electric power grid. Power system state estimation (PSSE) aims to obtain complete voltage magnitude and angle information at each bus given a number of system variables at selected buses and lines. Power flow analysis is a special case of PSSE, and amounts to solving a set of noise-free power flow equations. Physical laws dictate quadratic relationships between available quantities and unknown voltages, rendering general instances of power flow and PSSE nonconvex and NP-hard. Past approaches are largely based on gradient-type iterative procedures or semidefinite relaxation (SDR). Due to nonconvexity, the solution obtained via gradient-type schemes depends on initialization, while SDR methods do not perform as desired in challenging scenarios. This paper puts forth novel \emph{feasible point pursuit} (FPP)-based solvers for power flow and PSSE, which iteratively seek feasible solutions for a nonconvex quadratically constrained quadratic programming (QCQP) reformulation of the weighted least-squares (WLS) problem. Relative to the prior art, the developed solvers offer superior performance at the cost of higher complexity. Furthermore, they converge to a stationary point of the WLS problem. As a baseline for comparing different estimators, the Cram{\' e}r-Rao lower bound (CRLB) is derived for the fundamental PSSE problem in this paper. Judicious numerical tests on several IEEE benchmark systems showcase markedly improved performance of our FPP-based solvers for both power flow and PSSE tasks over popular WLS-based Gauss-Newton iterations and SDR approaches.

Tables

Table 1. TABLE I: Empirical success rate on IEEE test systems with θ = 0.1 ​ π 𝜃 0.1 𝜋 \theta=0.1\pi .

Test case	$5$-bus	$9$-bus	$14$-bus	$24$-bus	$30$-bus	$39$-bus
FPP-based	$100\%$	$100\%$	$100\%$	$100\%$	$100\%$	$100\%$
SDR-based	$0$	$29\%$	$40\%$	$2\%$	$94\%$	$97\%$
GN-based	$100\%$	$100\%$	$100\%$	$100\%$	$87\%$	$64\%$

Table 2. TABLE II: Empirical success rate on IEEE test systems with θ = 0.3 ​ π 𝜃 0.3 𝜋 \theta=0.3\pi .

Test case	$5$-bus	$9$-bus	$14$-bus	$24$-bus	$30$-bus	$39$-bus
FPP-based	$100\%$	$100\%$	$100\%$	$100\%$	$100\%$	$100\%$
SDR-based	$4\%$	$33\%$	$15\%$	$10\%$	$0$	$0$
GN-based	$100\%$	$55\%$	$33\%$	$21\%$	$0$	$5\%$

Equations

$$\mathbf{i}=\mathbf{Y}\mathbf{v}.$$

$$\displaystyle Y_{mn}:=\left\{\begin{array}[]{ll}-y_{mn},{}&(m,n)\in\mathcal{E}\\ {y}_{nn}^{g}+\sum_{k\in\mathcal{N}_{n}}y_{nk},{}&m=n\\ 0,{}&{\text{otherwise}}\end{array}\right.$$

$$I_{mn}=y_{mn}^s{V}_m+y_{mn}(V_m-V_n)$$

$$\mathbf{p}+j\mathbf{q}=\mathrm{diag}(\mathbf{v})\overline{\mathbf{i}}=\mathrm{diag}(\mathbf{v})\overline{\mathbf{Y}}\overline{\mathbf{v}}$$

$$P_{mn}^f+j{Q}_{mn}^f$$

$$=\left(y_{mn}^s+{\overline{y}}_{mn}^s\right)V_m{\overline{V}}_m-{\overline{y}}_{mn}{V}_m{\overline{V}}_n$$

$$z_{\ell }={\mathbf{v}}^{\mathcal{H}}{\mathbf{H}}_{\ell }\mathbf{v},1\le \ell \le L$$

$${\mathbf{Y}}_n$$

$${\mathbf{Y}}_{mn}^f$$

$${\mathbf{Y}}_{mn}^t$$

$${\mathbf{H}}_{V,n}:={\mathbf{e}}_n{\mathbf{e}}_n^{\mathcal{T}}⪰0,\forall n\in \mathcal{N}$$

$${\mathbf{H}}_{P,n}:=\frac{1}{2}\left({\mathbf{Y}}_n+{\mathbf{Y}}_n^{\mathcal{H}}\right),{\mathbf{H}}_{Q,n}:=\frac{j}{2}\left({\mathbf{Y}}_n-{\mathbf{Y}}_n^{\mathcal{H}}\right)$$

$$\displaystyle\bm{H}_{P,mn}^{f}:=\frac{1}{2}\left(\bm{Y}_{mn}^{f}+\big{(}\bm{Y}_{mn}^{f}\big{)}^{\mathcal{H}}\right)$$

$$\displaystyle\bm{H}_{Q,mn}^{f}:=\frac{j}{2}\left(\bm{Y}_{mn}^{f}-\big{(}\bm{Y}_{mn}^{f}\big{)}^{\mathcal{H}}\right)$$

$$\displaystyle\bm{H}_{P,mn}^{t}:=\frac{1}{2}\left(\bm{Y}_{mn}^{t}+\big{(}\bm{Y}_{mn}^{t}\big{)}^{\mathcal{H}}\right)$$

$$\displaystyle\bm{H}_{Q,mn}^{t}:=\frac{j}{2}\left(\bm{Y}_{mn}^{t}-\big{(}\bm{Y}_{mn}^{t}\big{)}^{\mathcal{H}}\right).$$

$$\mathrm{find}$$

$$\mathrm{subject}\mathrm{to}$$

$$z_{\ell }={\mathbf{v}}^{\mathcal{H}}{\mathbf{H}}_{\ell }\mathbf{v}+{\eta }_{\ell },1\le \ell \le L$$

$$\hat{\mathbf{v}}:=\arg \underset{\mathbf{v}\in {\mathbb{C}}^N}{\min }\sum _{\ell =1}^L{w}_{\ell }{\left(z_{\ell }-{\mathbf{v}}^{\mathcal{H}}{\mathbf{H}}_{\ell }\mathbf{v}\right)}^2$$

$$\underset{\mathbf{v}\in {\mathbb{C}}^N,\{s_{\ell }{\}}_{l=1}^L}{\mathrm{minimize}}$$

$$\mathrm{subject}\mathrm{to}$$

$$\underset{\mathbf{v}\in {\mathbb{C}}^N,\{s_{\ell }{\}}_{l=1}^L}{\mathrm{minimize}}$$

$$\mathrm{subject}\mathrm{to}$$

$$\underset{\mathbf{v}\in {\mathbb{C}}^N,\mathbf{s}\in {\mathbb{R}}^L}{\mathrm{minimize}}$$

$$\mathrm{subject}\mathrm{to}$$

$${\mathbf{v}}^{\mathcal{H}}\left(-{\mathbf{H}}_{\ell }\right)\mathbf{v}\le -z_{\ell }+s_{\ell },1\le \ell \le L.$$

$${\mathbf{v}}^{\mathcal{H}}{\mathbf{H}}_{\ell }^{(+)}\mathbf{v}+{\mathbf{v}}^{\mathcal{H}}{\mathbf{H}}_{\ell }^{(-)}\mathbf{v}\le z_{\ell }+s_{\ell }$$

$$(\mathbf{v}-\mathbf{y}{)}^{\mathcal{H}}{\mathbf{H}}_{\ell }^{(-)}(\mathbf{v}-\mathbf{y})\le 0.$$

$${\mathbf{v}}^{\mathcal{H}}{\mathbf{H}}_{\ell }^{(-)}\mathbf{v}\le 2\Re \left\{{\mathbf{y}}^{\mathcal{H}}{\mathbf{H}}_{\ell }^{(-)}\mathbf{v}\right\}-{\mathbf{y}}^{\mathcal{H}}{\mathbf{H}}_{\ell }^{(-)}\mathbf{y}.$$

$${\mathbf{v}}^{\mathcal{H}}{\mathbf{H}}_{\ell }^{(+)}\mathbf{v}+2\Re \left\{{\mathbf{y}}^{\mathcal{H}}{\mathbf{H}}_{\ell }^{(-)}\mathbf{v}\right\}\le z_{\ell }+{\mathbf{y}}^{\mathcal{H}}{\mathbf{H}}_{\ell }^{(-)}\mathbf{y}+s_{\ell }.$$

$$\bm{v}^{\mathcal{H}}\big{(}\!-\!\bm{H}_{\ell}^{(-)}\big{)}\bm{v}\!-2\Re\!\left\{\!\bm{y}^{\mathcal{H}}\bm{H}_{\ell}^{(+)}\bm{v}\!\right\}\!\leq\!-z_{\ell}-\bm{y}^{\mathcal{H}}\bm{H}_{\ell}^{(+)}\bm{y}+s_{\ell}.$$

$$\left\{{\mathbf{v}}_{k+1},{\mathbf{s}}_{k+1}\right\}:=\arg \underset{\mathbf{v},\mathbf{s}}{\min }\sum _{l=1}^L{w}_{\ell }{s}_{\ell }^2$$

$$\mathrm{subject}\mathrm{to}$$

$${\mathbf{v}}^{\mathcal{H}}{\mathbf{H}}_{\ell }^{(+)}\mathbf{v}+2\Re \left\{{\mathbf{y}}^{\mathcal{H}}{\mathbf{H}}_{\ell }^{(-)}\mathbf{v}\right\}\le z_{\ell }+{\mathbf{y}}^{\mathcal{H}}{\mathbf{H}}_{\ell }^{(-)}\mathbf{y}+s_{\ell }$$

$${\mathbf{v}}^{\mathcal{H}}{\mathbf{H}}_{\ell }^{(-)}\mathbf{v}+2\Re \left\{{\mathbf{y}}^{\mathcal{H}}{\mathbf{H}}_{\ell }^{(+)}\mathbf{v}\right\}\ge z_{\ell }+{\mathbf{y}}^{\mathcal{H}}{\mathbf{H}}_{\ell }^{(+)}\mathbf{y}-s_{\ell }$$

$$\forall \ell =1,2,\dots ,L$$

$${\overline{\mathbf{H}}}_{\ell }:=\left[\begin{array}{cc}\Re ({\mathbf{H}}_{\ell })& -\Im ({\mathbf{H}}_{\ell })\\ \Im ({\mathbf{H}}_{\ell })& \Re ({\mathbf{H}}_{\ell })\end{array}\right]\in {\mathbb{R}}^{2N\times 2N},1\le \ell \le L.$$

Full text

I Introduction

Recognized as the greatest engineering achievement of the $21$st century [1], the electric power grid is a complex cyber-physical system comprising multiple subsystems, each with a transmission infrastructure to deliver electricity from power generators to distribution networks to customers. Accurately monitoring the operational condition of a power grid is crucial to various system control and optimization tasks, which include unit commitment, optimal power flow (OPF), and economic dispatch [2, 3]. To enable such an accurate monitoring, a set of system variables are specified (and enforced) or measured at selected buses and lines for determining or estimating the system’s operating point, namely complex voltages at all buses of the grid. These two tasks correspond to the so-termed power flow analysis and power system state estimation (PSSE), respectively. Both are central to monitoring, control, and future planning of electricity networks.

In power engineering, power flow analysis is a numerical analysis of the normal steady-state flow of electric power over the grid, that is crucial for planning future power system expansions (e.g., designing components such as generators, lines, transformers, and capacitors), as well as in determining the best operation of the existing systems [3]. The goal of power flow analysis is to obtain complete voltage magnitude and angle information at each bus for specified or enforced load and generator active power and voltage conditions [3]. Once this information is available, other system variables including active and reactive power flows as well as generator reactive power outputs can be analytically obtained.

Power flow analysis amounts to solving a set of quadratic equations given by the nonlinear AC power flow model obeying Ohm’s and Kirchhoff’s laws. Solving power flow equations for both transmission and distribution systems is known to be NP-hard [4]. Due to the nonlinear nature, several numerical solvers have been developed to obtain a solution that is within an acceptable tolerance [3, 5]. Past solvers include the Gauss-Seidel and Newton-Raphson iterative algorithms [3], and the semidefinite relaxation (SDR) [5, 6]. The Gauss-Seidel method is reported as the earliest devised power flow solver [3]. On the other hand, the Newton-Raphson algorithm iteratively seeks improved approximations to the zeros of real-valued functions, featuring quadratic convergence whenever the initial point lands within a small neighborhood of the zeros [7]. As convergence of both algorithms relies heavily on the initial point, they may diverge if the initialization is not reliable [6]. With a carefully designed objective function and sufficiently small angle differences across lines, the SDR approaches have been shown capable of recovering the true power flow solution provided that the set of available specifications includes all voltage magnitudes, and the active power flows over a spanning tree of the network [5, 6].

The task of PSSE can be described as estimating the voltage magnitudes and angles at all buses across the network from a subset of supervisory control and data acquisition (SCADA) measurements including active and reactive power injections and flows (at both the sending and receiving ends), as well as squared voltage magnitudes [2]. Since its appearance in the $1970$s [8], PSSE has become a prerequisite for supervisory control, system planning, and economic dispatch [8]. Nonlinear SCADA measurements however, render the PSSE problem nonconvex and NP-hard in general [9].

PSSE solvers so far are largely based on Gauss-Newton iterations and SDR heuristics. The “workhorse” Gauss-Newton method for nonconvex optimization has two limitations [10, Sec. 1.5], i.e., sensitivity to the initial guess, and lack of convergence guarantees. SDR-based approaches on the other hand solve first for a matrix variable that can be computationally expensive [9, 11, 12, 6, 5, 13]. SDR’s performance degrades when the data-size is relatively small, or when the data do not include all voltage magnitudes [6, 5]. For PSSE of large-scale networks, distributed Gauss-Newton and SDR implementations have been reported in [9], [14], [15].

Solving power flow equations and the PSSE can be shown equivalent to solving nonconvex QCQPs, which is generally NP-hard [16]. Many heuristics have been put forward. A feasible point pursuit (FPP) algorithm developed in [17] was shown to enjoy improved performance over the SDR-based methods. The FPP heuristic has been employed for solving the OPF problem [18], where the resulting solver was empirically shown more effective for multi-phase transmission networks than popular SDR- and moment relaxation-based ones [18].

Building on our precursors [17, 19] and inspired by the inherent nonconvex challenge, the goal of this work is to develop power flow and PSSE solvers capable of attaining or approximating the global optimum at manageable computational complexity. Starting with the WLS formulation, the power flow and PSSE tasks are equivalently reformulated as a nonconvex QCQP, which can be readily tackled by FPP. We further show that our FPP-based solvers converge to a stationary point of the WLS problem. As a baseline for comparing different SE approaches, the Cramér-Rao lower bound (CRLB) is derived for the basic PSSE problem under additive white Gaussian noise (AWGN). This is achieved by means of Wirtinger’s calculus for functional analysis over complex domains. Finally, numerical experiments using several IEEE benchmark systems corroborate the superior performance of our proposed solvers over existing methods for both power flow and PSSE tasks.

Regarding notation, matrices (vectors) are denoted by upper- (lower-) case boldface letters, and $\overline{(\cdot)}$, $(\cdot)^{\mathcal{T}}$, and $(\cdot)^{\mathcal{H}}$ stand for complex conjugate, transpose, and conjugate-transpose, respectively. Calligraphic letters are reserved for sets, e.g., $\mathcal{N}$. Symbol ${\Re}\{\cdot\}$ (${\Im}\{\cdot\}$) takes the real (imaginary) part of a complex-valued object, and ${\rm diag}(\bm{x})$ is a diagonal matrix holding in order entries of $\bm{x}$ on its diagonal.

II System Modeling and Problem Statement

An electric transmission network having $N$ nodes (buses) and $E$ edges (lines) can be represented by a graph $\mathcal{G}:=\{\mathcal{N},\mathcal{E}\}$, whose nodes $\mathcal{N}:=\{1,2,\ldots,N\}$ correspond to buses, and whose edges $\mathcal{E}:=\{(m,n)\}\subseteq\mathcal{N}\times\mathcal{N}$ correspond to transmission lines. For every bus $n\in\mathcal{N}$, let $V_{n}:=|V_{n}|e^{j\theta_{n}}$ be the nodal complex voltage, whose magnitude and phase are given by $|V_{n}|$ and $\theta_{n}$, respectively; likewise for the complex current injection $I_{n}:=|I_{n}|e^{j\phi_{n}}$. Let also $S_{n}:=P_{n}+jQ_{n}$ be the corresponding complex power injection, in which $P_{n}$ and $Q_{n}$ are the active and reactive power injection, respectively. For every line $(m,n)\in\mathcal{E}$, let $I_{mn}$ denote the complex current flowing from bus $m$ to $n$, and $S_{mn}^{f}:=P_{mn}^{f}+jQ_{mn}^{f}$ the complex power flow from bus $m$ to $n$ seen at the sending end, where $P_{mn}^{f}$ and $Q_{mn}^{f}$ are the active and reactive power flow, respectively; and likewise for the receiving-end (active and reactive) power flow $P_{mn}^{t}$ and $Q_{mn}^{t}$.

The AC power flow model dictates that system variables $\{P_{n}\}$, $\{Q_{n}\}$, $\{P_{mn}^{f}\}$, $\{Q_{mn}^{f}\}$, $\{P_{mn}^{t}\}$, $\{Q_{mn}^{t}\}$, and $\{|V_{n}|^{2}\}$ are quadratic functions of the state vector $\bm{v}$. Clearly, this holds true for the squared voltage magnitude understood as $|V_{n}|^{2}=V_{n}\overline{V}_{n}$. To specify the relationship between power quantities and $\bm{v}$, introduce $\bm{Y}\in\mathbb{C}^{N\times N}$ to represent the bus admittance matrix, which is in general symmetric. Ohm’s law in conjunction with Kirchhoff’s law reads as

[TABLE]

It is worth mentioning that $\bm{Y}$ is sparse, thus enabling efficient computations in large-size power networks, and its $(m,n)$-th entry is given by

[TABLE]

where $y_{mn}$ denotes the admittance of line $(m,n)\in\mathcal{E}$, ${y}_{nn}^{g}$ the admittance to the ground at bus $n\in\mathcal{N}$, and $\mathcal{N}_{n}$ the set of neighboring buses directly connected to bus $n$. For $m\neq n$, let ${y}_{mn}^{s}$ be the shunt admittance at bus $m$ associated with line $(m,n)$. Recall from Ohm’s and Kirchhoff’s laws that the current flowing from bus $m$ to $n$ can be expressed as

[TABLE]

whereby the reverse-direction current $I_{nm}$ can be given symmetrically. Due to ${y}^{s}_{mn}\neq 0$ in general, it holds $I_{mn}\neq-I_{nm}$.

The AC model also asserts $P_{n}+jQ_{n}=V_{n}\overline{I}_{n},\,\forall n\in\mathcal{N}$. Appealing again to (1) leads to the next matrix-vector form

[TABLE]

where both active and reactive power injections are quadratically related to $\bm{v}$. Likewise, the sending-end active and reactive power flow over line $(m,n)\in\mathcal{E}$ can be written as

[TABLE]

where the second equality is obtained by substituting $I_{mn}$ in (6) into the first. Hence, $P_{mn}^{f}$ and $Q_{mn}^{f}$ can also be expressible as quadratic functions of $\bm{v}$. By symmetry, this quadratic relationship also holds for $P_{mn}^{t}$ and $Q_{mn}^{t}$.

To perform either power flow analysis or PSSE, a total of $L$ system variables are specified or measured by the system operator. The nonlinear AC networks have available the next seven types of quantities: $|V_{n}|^{2}$, $P_{n}$, $Q_{n}$, $P_{mn}^{f}$, $Q_{mn}^{f}$, $P_{mn}^{t}$, and $Q_{mn}^{t}$. If $\mathcal{N}_{V}$, $\mathcal{N}_{P}$, $\mathcal{N}_{Q}$, $\mathcal{E}_{P}^{f}$ ($\mathcal{E}_{Q}^{f}$), and $\mathcal{E}_{P}^{t}$ ($\mathcal{E}_{Q}^{t}$) denote the selected sets of buses/lines over which actual quantities of the corresponding type are available, the elaborated quadratic relationships prompt us to define the $L\times 1$ data vector \bm{z}:=\big{[}\{|V_{n}|^{2}\}_{n\in\mathcal{N}_{V}},\;\{P_{n}\}_{n\in\mathcal{N}_{P}},\;\{Q_{n}\}_{n\in\mathcal{N}_{Q}},\;\{P_{mn}^{f}\}_{(m,\,n)\in\mathcal{E}_{P}^{f}}, \{Q_{mn}^{f}\}_{(m,n)\in\mathcal{E}_{Q}^{f}},\,\{P_{mn}^{t}\}_{(m,n)\in\mathcal{E}_{P}^{t}},\,\{Q_{mn}^{t}\}_{(m,n)\in\mathcal{E}_{Q}^{t}}\big{]}^{\mathcal{T}}\in\mathbb{R}^{L}, whose entries can be succinctly given by

[TABLE]

where $\{\bm{H}_{\ell}\}_{\ell=1}^{L}$ are some coefficient matrices to be specified. For this purpose, let $\left\{\bm{e}_{n}\in\mathbb{R}^{N}\right\}_{n=1}^{N}$ be the canonical basis of $\mathbb{R}^{N}$, and introduce also the admittance-dependent matrices

[TABLE]

For $|V_{n}|^{2}=V_{n}\overline{V}_{n}=\bm{v}^{\mathcal{H}}\bm{e}_{n}\bm{e}_{n}^{\mathcal{T}}\bm{v}$, it is clear that the corresponding $\bm{H}_{n}$ in (9) is

[TABLE]

which are rank-$1$. By taking separately the real and imaginary parts of (7) and (8), we obtain the $\{\bm{H}_{\ell}\}$ associated with the active and reactive power injections for all buses $n\in\mathcal{N}$

[TABLE]

and with sending-end and receiving-end active and reactive power flow at all lines $(m,n)\in\mathcal{E}$

[TABLE]

It is worth stressing that all $\{\bm{H}_{\ell}\}$ in (11) and (12) are sparse, low-rank, and Hermitian, but they are non-definite in general. The power flow and PSSE problems are formulated in order next.

II-A Power flow analysis

Power flow analysis deals with specified power quantities, which are enforced for optimally operating an electric power grid. Specifically, given $L$ perfectly known specifications $\{z_{\ell}\}_{\ell=1}^{L}$ and valid network parameters $\{\bm{H}_{\ell}\}_{\ell=1}^{L}$ as in (9), the goal of power flow analysis is to decide the state vector $\bm{v}\in\mathbb{C}^{N}$ that satisfies all specifications, namely,

[TABLE]

Recall that each bus in a power system is classified as a PQ, PV, or slack (reference) bus based on the constraints imposed per bus. PQ buses, which often correspond to loads, specify and enforce only active and reactive power injection $P_{n}$ and $Q_{n}$ on bus $n$. On the other hand, the PV buses, which are typically associated with generators, enforce active power injection $P_{n}$ and voltage magnitude $|V_{n}|$. For the slack bus, its voltage phase is fixed at $\theta_{n}=0$, by convention. With $\theta_{n}=0$, the power flow problem in (13) is equivalent to solving for $2N-1$ real-valued unknowns from $L$ quadratic equations. The classical power flow problem considers the particular case where the $L=2N-1$ specifications are enforced only at the PV, PQ, and slack buses as opposed to a combination of buses and lines [5].

II-B Power system state estimation

PSSE on the other hand deals with noisy observations acquired by the SCADA system adhering to

[TABLE]

where $\eta_{\ell}$ accounts for the zero-mean distributed measurement error with known variance $\sigma_{\ell}^{2}$, henceforth assumed independent across meters. The goal of PSSE is, given SCADA measurements $\{z_{\ell}\in\mathbb{R}\}_{\ell=1}^{L}$ and also parameters $\{\bm{H}_{\ell}\}_{\ell=1}^{L}$, estimate the state vector $\bm{v}\in\mathbb{C}^{N}$.

Adopting the WLS criterion, the SE task can be cast as that of solving the following nonlinear LS problem

[TABLE]

where entries of the weight vector $\bm{w}:=[w_{1}~{}\cdots~{}w_{L}]^{\mathcal{T}}$ are often taken as $w_{\ell}:=1/\sigma_{\ell}^{2}$ for known $\sigma_{\ell}^{2}$ values. The WLS estimate $\hat{\bm{v}}$ coincides with the maximum likelihood one when the error vector $\boldsymbol{\eta}:=[\eta_{1}~{}\cdots~{}\eta_{L}]^{\mathcal{T}}$ obeys the multivariate Gaussian distribution $\mathcal{N}(\bm{0},{\rm diag}(\boldsymbol{\sigma}^{2}))$ with $\bm{\sigma}^{2}:=[\sigma_{1}^{2}~{}\cdots~{}\sigma_{L}^{2}]^{\mathcal{T}}$. Unfortunately, due to the quadratic terms $\{\bm{v}^{\mathcal{H}}\bm{H}_{\ell}\bm{v}\}$ inside the squares, the WLS SE problem is nonconvex. Minimizing nonconvex objectives, which typically exhibit many stationary points, is NP-hard in general [16]. Hence, solving the problem in (15) is indeed challenging.

PSSE approaches so far can be grouped as convex and nonconvex ones. The latter includes the “workhorse” Gauss-Newton method, which is also typically employed in practice: Upon linearizing the error function in the LS cost around a given estimate, the minimizer of the norm of the resulting linearized approximation is used to initialize the next iteration [10, Sec. 1.5]. Minimizing nonconvex functions, Gauss-Newton iterations can be problematic due to: i) its sensitivity to the initial point; and, ii) lack of convergence guarantee to even a stationary point [10]. Convex approaches via SDR [9, 6] express all data $\{z_{\ell}\}$ as linear functions of the outer-product $\bm{V}:=\bm{v}\bm{v}^{\mathcal{H}}\in\mathbb{C}^{N\times N}$. Problem (15) is then convexified by dropping the nonconvex constraint ${\rm rank}(\bm{V})=1$. SDR-based methods seldom yield solutions of rank-$1$ in the noisy case. Further eigen-decomposition or randomization procedures are required to recover the estimator $\hat{\bm{v}}$ from the SDR solution $\hat{\bm{V}}$. Performance of SDR solutions degrades when the data size is small, or when the set of measurements does not include the voltage magnitude at all buses, as will be demonstrated by our numerical results in Sec. V.

III Feasible Point Pursuit based Solvers

In this section, the FPP-based power flow and PSSE solvers will be developed based on procedures distinct from existing iterative optimization and SDR-based SE approaches. To this end, some basics of FPP are first reviewed. For nonconvex QCQPs, FPP iteratively solves a series of convexified QCQPs obtained with successive convex inner-restrictions of the original nonconvex feasibility set, and with additive slacks to approximate the feasible solutions of the original nonconvex QCQP [17]. Specifically, starting with an initial guess, FPP first decomposes the quadratic terms in all nonconvex constraints into their convex and nonconvex parts by means of eigen-decomposition, which can be efficiently carried out offline; then it linearizes the nonconvex parts around the current iterate to obtain a restricted convex QCQP. Due to restriction of the feasibility set, the convexified QCQP may be infeasible. To sustain feasibility, a slack variable is introduced for each relaxed constraint, with a convex penalty on the slack variables added to the cost function, which can enforce sparing use of slacks to produce solutions of minimal constraint violation. The minimizer of the regularized convex QCQP subproblem is taken as the next iterate, which will be used as the linearization point of the nonconvex components at the next iteration. This successive convex approximation and feasibility-restoring procedure is repeated until a certain stopping criterion is met. Further details of FPP can be found in [17], [18].

Note that the power flow problem (13) consists of quadratic equality constraints, which are not in the standard QCQP form. To apply FPP, equalities are relaxed to inequalities, while penalizing the slack variables $\bm{s}:=\{s_{\ell}\geq 0\}_{\ell=1}^{L}$, yielding

[TABLE]

where other choices of the convex penalty function $f(\cdot)$ include the (weighted) $\ell_{1}$ or $\ell_{\infty}$ norm. Problem (16) is equivalent to the original power flow formulation (13) when the latter is feasible. To see this, assume that the set of power flow equations in (13b) admits (possibly more than one) feasible solutions. Clearly at the optimum of (16), the objective reduces to zero, the slack variables $\{s_{\ell}\}_{\ell=1}^{L}$ take zero values, and all equalities in (16b) are achieved, thus yielding a feasible solution to the set of power flow equations in (13).

Similarly, our PSSE formulation in (15) minimizes a quartic polynomial of $\bm{v}$. To use FPP, problem (15) is reformulated as

[TABLE]

where the slack variables $\bm{s}:=\{s_{\ell}\geq 0\}_{\ell=1}^{L}$ in this case relate to the deviations between noisy measurements $\{z_{\ell}\}_{\ell=1}^{L}$ and the actual quantities $\{\bm{v}^{\mathcal{H}}\bm{H}_{\ell}\bm{v}\}_{\ell=1}^{L}$. Problem (17) can be similarly shown equivalent to (15). Other convex penalty functions $f(\cdot)$ in (17a) can also be selected. In particular, if the error vector follows the multivariate Laplace distribution, i.e., $\bm{\eta}\sim{\rm Laplace}(\bm{0},\,\bm{b})$ with $\bm{b}:=[b_{1}~{}\cdots~{}b_{L}]^{\mathcal{T}}$ collecting all scaling parameters, minimizing the $\ell_{1}$-based function $f(\bm{s})=\sum_{\ell=1}^{L}w_{\ell}s_{\ell}$ with $w_{\ell}=1/b_{\ell}$ in (17) produces the maximum likelihood estimate [2, 6].

Apparently, the reformulated power flow and PSSE problems are of the same form [cf. (16) and (17)], except for a minor difference in the cost functions. Setting unit weights $w_{\ell}=1$ in (17) reduces problem (17) to (16). Without loss of generality, we will hereafter focus on the PSSE formulation (17), and develop the novel FPP solver. The power flow problem can be readily handled with all weights being $w_{\ell}=1$.

In this direction, let us first convert problem (17) into a standard QCQP. Note that constraints (17b) can be replaced by two sets of inequalities to arrive at

[TABLE]

Problem (18) is nonconvex even for (semi)definite coefficient matrices $\{\bm{H}_{\ell}\}_{\ell=1}^{L}$. Next we demonstrate how to take advantage of FFP to solve the problem at hand in detail.

As discussed in Sec. II, there are two types of $\{\bm{H}_{\ell}\}$ matrices, one corresponding to the squared voltage magnitude, and the other to power quantities. Type-I $\{\bm{H}_{\ell}\}$ are positive semidefinite [cf. (10)], while Type-II are non-definite [cf. (11) and (12)]. For ease of exposition, let us introduce the FPP constraint convexification procedure using one nonconvex quadratic constraint in (18). Along the lines of FPP, consider the term $\bm{v}^{\mathcal{H}}\bm{H}_{\ell}\bm{v}$ in (18b) for some $\bm{H}_{\ell}$ in (11), which can be decomposed into its convex and nonconvex components as

[TABLE]

where $\bm{H}_{\ell}^{(+)}$ and $\bm{H}_{\ell}^{(-)}$ represent the positive semidefinite (convex) and negative semidefinite (nonconvex) parts of $\bm{H}_{\ell}$ in (19), respectively. For the nonconvex source $\bm{v}^{\mathcal{H}}\bm{H}_{\ell}^{(-)}\bm{v}$ in (19), an inner linear restriction will be derived next.

The following inequality holds for any $\bm{y}\in\mathbb{C}^{N}$ due to the negative semidefiniteness of $\bm{H}_{\ell}^{(-)}$

[TABLE]

Upon expanding the left-hand-side and rearranging terms, one arrives at

[TABLE]

Key to the FPP algorithm is replacing the nonconvexity stemming from $\bm{H}_{\ell}^{(-)}$ in (19) or (18b) by its inner linear approximation at some given point $\bm{y}$ to yield

[TABLE]

The strategy in selecting the linearization point $\bm{y}$ will be discussed shortly. In the same fashion, the nonconvex quadratic constraints in (18c) can be replaced by

[TABLE]

Heed that the flexibility introduced by the slacks $\{s_{\ell}\}_{\ell=1}^{L}$ always restores the feasibility of the relaxed constraints, which contributes to improved performance of FPP over other convexification approaches [17]. In the presence of noise, the minimum values required for $\{s_{\ell}\geq 0\}_{\ell}^{L}$ to satisfy (21) and (22) depend on the measurement error contained in $\{z_{\ell}\}_{\ell=1}^{L}$.

The FPP method replaces all nonconvex constraints in (18b) by their convex restriction (21), and those in (18c) by (22) to derive a convexified QCQP regularized with slack variables to ensure feasibility. Minimizing some convex penalty function of the slacks $\{s_{\ell}\}_{\ell=1}^{L}$ not only minimizes the fitting error between $\{z_{\ell}\}$ and $\{\bm{v}^{\mathcal{H}}\bm{H}_{\ell}\bm{v}\}$, but also enforces sparing use of slacks and promotes solutions of minimal constraint violation.

In a nutshell, the developed FPP-based PSSE solver can be understood as follows. Starting with an initial point $\bm{v}_{0}$ (typically the flat voltage profile point, i.e., all-ones vector), our FPP-based solver successively tackles a sequence of convexified QCQPs with the linearization point being the current iterate $\bm{v}_{k}$, which is the $\bm{v}$-minimizer obtained by solving a convexified QCQP at the previous iteration. Hence, assuming available the $\bm{v}$-minimizer $\bm{v}_{k}$ at the $(k+1)$-st iteration, our FPP-based solver boils down to solving the following convexified QCQP subproblem

[TABLE]

where $\bm{y}:=\bm{v}_{k}$ is the $\bm{v}$-minimizer of (23) at the $k$-th iteration. The QCQP in (23) is convex, which can be solved in polynomial time using off-the-shelf solvers [20]

The FPP-based PSSE solver is summarized in Alg. 1. The following three properties of our FPP-based solver are worth highlighting.

Remark 1* *(Power flow analysis).

Cast as a special instance of PSSE, the power flow problem in (13) can be solved by our developed FPP-based PSSE solver with unit weights $w_{\ell}=1$.

Remark 2* *(Bad data removal).

Besides the $\ell_{2}$-norm in (23a), other convex penalty functions can be used to fit different (noisy) data models. In particular, adopting the weighted $\ell_{1}$-norm (i.e., replacing $s_{\ell}^{2}$ with $|s_{\ell}|$) yields the weighted least-absolute-value estimator known for bad data cleansing [21].

Remark 3* *(Synchrophasors).

Synchrophasors, if available, can be easily incorporated into the developed PSSE formulation (23). To see this, letting $\bm{\zeta}_{n}=\bm{\Phi}_{n}\bm{v}+\bm{\epsilon}_{n}$ collect the noisy PMU data at bus $n$, hybrid estimation exploiting both nonlinear SCADA measurements and linear PMU ones can be achieved [22] with an additional data-fitting term for the PMU data in (23a), namely, $\sum_{n\in\mathcal{P}}\|\bm{\zeta}_{n}-\bm{\Phi}_{n}\bm{v}\|_{2}^{2}$, where $\mathcal{P}$ denotes the subset of the PMU-instrumented buses.

On the theoretical side, the following result establishes convergence of our developed FPP-based solvers to a stationary point of the WLS formulation.

Proposition 1** (Global convergence of FPP-based solvers).**

Let $\{\bm{v}_{k}\}_{k=0}^{\infty}$ be any sequence generated by the FPP-based solver in Alg. 1. Then, all limit points of $\{\bm{v}_{k}\}_{k=0}^{\infty}$ are stationary points of the WLS problem in (15).

Proof.

As elaborated in Sec. III, solving problem (18) is equivalent to solving problem (15). The nonconvex QCQP of complex-valued vector $\bm{v}\in\mathbb{C}^{N}$ in (18) can be equivalently posed as a QCQP of the expanded real-valued vector ${\bm{u}}:=[\Re(\bm{v})^{T}~{}\Im(\bm{v})^{T}]^{T}\in\mathbb{R}^{2N}$, where the associated quadratic matrices $\{\overline{{\bm{H}}}_{\ell}\}$ are given as

[TABLE]

Accordingly, each constraint in (18) can be re-expressed as the difference between two convex functions. To see this, consider e.g. constraint (18b), which can be rewritten as

[TABLE]

where $\overline{\bm{H}}_{\ell}^{(+)}$ and $\overline{\bm{H}}_{\ell}^{(-)}$ are the positive and negative semidefinite parts of $\overline{\bm{H}}_{\ell}$, hence rendering terms $\bm{u}^{\mathcal{T}}\overline{\bm{H}}_{\ell}^{(+)}\bm{u}-s_{\ell}$ and $\bm{u}^{\mathcal{T}}(-\overline{\bm{H}}_{\ell}^{(-)})\bm{u}$ both convex. Alg. 1 is tantamount to an application of the convex-concave procedure [23, 24] to the reformulated QCQP in the real domain. Hence, the sequence generated by Alg. 1 converges to a stationary point of (15) by appealing to the results in [25, Thm. 10]. ∎

IV Cramér-Rao Bound for PSSE

According to standard results from estimation theory [26], the variance of any unbiased estimator is lower bounded by the Cramér-Rao lower bound (CRLB). Appreciating its key role as a performance benchmark across different estimators, this section establishes the CRLB for the fundamental PSSE problem. The CRLB analysis of PSSE however, entails finding derivatives (gradient and Hessian) of a real-valued function with respect to multiple complex-valued variables. To address this challenge, we call for advanced complex analysis tools based on the so-termed Wirtinger derivative and Wirtinger’s calculus, which are detailed in the Appendix. The following result provides a closed-form CRLB for any unbiased PSSE solver under the AWGN model in (14), which can be directly used to assess the performance of other PSSE solvers.

Proposition 2**.**

Consider estimating the unknown state vector $\bm{v}\in\mathbb{C}^{N}$ from noisy data $\{z_{\ell}\}_{\ell=1}^{L}$ obeying the model in (14), where the noise $\eta_{\ell}$ is assumed Gaussian distributed with mean zero and variance $\sigma_{\ell}^{2}$, and is also independent across meters. Then the covariance matrix of any unbiased estimator $\hat{\bm{v}}$ obeys

[TABLE]

where the Fisher information matrix is given by

[TABLE]

Furthermore, $\bm{F}$ has at least rank-$1$ deficiency even when all possible SCADA measurements are available.

The proof of Prop. 2 is deferred to the Appendix. Even though the Fisher information matrix (FIM) in (26) is rank deficient, the pseudo-inverse of FIM qualifies itself as a valid yet generally looser lower bound on the mean-square error (MSE) of any unbiased estimator [27]. This lower bound is often attainable in practice, and is predictive of optimal estimator performance [27], as will be demonstrated by our numerical tests in Sec. V. The derived CRLB in (25) will be employed to benchmark and compare performance of different PSSE solvers next.

V Simulated Tests

In the section, we compare the proposed FPP-based solvers in Alg. 1 with existing alternatives including the WLS via Gauss-Newton iterations (GN-based), and the SDR-based solver (SDR-based) [9] for both power flow and PSSE tasks on several IEEE benchmark systems [28]. Throughout, all reported numerical results were obtained by averaging over $100$ independent Monte Carlo realizations. The three PSSE solvers from noisy measurements are compared in terms of the mean-square error $\sum_{i=1}^{100}\left\|\hat{\bm{v}}_{i}-\bm{v}\right\|_{2}^{2}/100$, where $\hat{\bm{v}}_{i}$ is the returned estimate at the $i$-th realization, and $\bm{v}$ the actual voltage profile. In the absence of noise, performance of the power flow solvers is assessed through the empirical success rate over $100$ trials. A success is declared for a trial if the returned power flow solution $\hat{\bm{v}}$ incurs a relative violation on the given set of $L$ power flow equations, given by $\sum_{\ell=1}^{L}(z_{\ell}-\hat{\bm{v}}^{\mathcal{H}}\bm{H}_{\ell}\hat{\bm{v}})^{2}/\sum_{\ell=1}^{L}z_{\ell}^{2}$ less than $10^{-3}$. (The reason why $\|\bm{v}-\hat{\bm{v}}\|_{2}^{2}$ is not used is due to existence of possibly multiple solutions $\bm{v}$ satisfying the set of power flow equations.)

Different system quantities and voltage profiles were generated via the MATLAB-based toolbox MATPOWER [29]. The Gauss-Newton method was implemented using the SE function ‘doSE.m’ in MATPOWER, which was modified to terminate either upon convergence, or, when the condition number of the approximate linearization exceeds $10^{5}$ flagging explosion of the iterates [9]. The SDR- and FPP-based solvers were realized via the optimization modeling package YALMIP [30], as well as the interior-point solver SeDuMi [20]. Furthermore, the flat-voltage profile point was used as the initial guess for the Gauss-Newton and FPP approaches. In order to fix the phase ambiguity, the phase generated at the reference bus is set to [math] in all tests. The FPP solver stops either when a maximum number $100$ of iterations are reached, or when the objective value improvement between two consecutive iterations becomes smaller $10^{-5}$. All experiments were conducted on an Intel CPU @ $3.4$ GHz ($32$ GB RAM) computer.

To evaluate the performance of the FPP-based solver for power flow analysis, the first experiment simulates noiseless data corresponding to the classical power flow problem. That is, a total of $L=2N-1$ system variables were specified at the PV, PQ, and slack buses to solve for $2N-1$ real-valued unknowns in $\bm{v}\in\mathbb{C}^{N}$ with the reference bus’s phase fixed at [math]. The actual voltage magnitude of each bus was uniformly distributed over $[0.9,\,1.1]$, and its angle over $[-\theta,\,\theta]$ with $\theta=0.1\pi$ and $0.3\pi$. Empirical success rate results on several IEEE benchmark systems were reported in Tables I and II for $\theta=0.1\pi$ and $0.3\pi$, respectively. Apparently, our developed FPP-based power flow solver solves exactly the classical power flow problem in all simulated tests, while the SDR-based one fails with high probability. The Gauss-Newton method performs well when the initial point lies close to the actual solution due to small $\theta$ in Table I, while it diverges frequently for large $\theta$ values in Table II.

The second experiment compares the MSE performance of various approaches relative to the analytical Cramér-Rao bound in (25) on the IEEE $14$-bus test system [28]. The actual voltage magnitude and angle of each bus were generated uniformly over $[0.9,\,1.1]$, and $[-0.4\pi,\,0.4\pi]$, respectively. Initially, all voltage magnitudes as well as all sending-end and receiving-end active power flow were taken, which corresponds to the base case $3$ in the $x$-axis of Fig. 1. To demonstrate the SE performance evolution of various approaches with respect to the increasing number of measurements, additional types of measurements were included in a deterministic manner detailed next. All seven types of SCADA measurements in (10)-(12) were ordered as $\{|V_{k}|^{2},P_{mn}^{f},P_{mn}^{t},Q_{mn}^{f},Q_{mn}^{t},P_{n},Q_{n}\}$. Each $x$-axis value in Fig. 1 implies that the number of ordered types of measurements was used in the experiment to obtain the mean-square errors. For instance, $5$ on the $x$-axis corresponds to the case where the first $5$ types of measurements (i.e., all $|V_{k}|^{2},P_{mn}^{f},P_{mn}^{t},Q_{mn}^{f},Q_{mn}^{t}$) were used; and likewise for all other $x$-axis values. Measurement noise was randomly and independently generated from Gaussian distribution having zero-mean and standard deviation $0.1$. The SDR estimator was recovered from the SDR solution by picking the minimum-cost vector over the eigenvector and $5,000$ zero-mean Gaussian randomizations with covariance matrix being the SDR solution. The MSE as well as the CRLB versus the types of measurements available are shown in Fig. 1, corroborating the near-optimal performance relative to the CRLB and robustness of our developed FPP-based PSSE solver.

The last experiment on the IEEE $30$-bus benchmark system simulates a high signal-to-noise ratio and complete-data scenario, where all voltage magnitude as well as all active power flow at both sending- and receiving-ends were measured to be advantageous to the SDR-based method [6]. Independent zero-mean Gaussian noise was assumed to have standard deviations $0.05$ for power measurements and $0.02$ for voltage measurements. The actual voltage magnitude and angle of each bus were generated uniformly at random over $[0.9,\,1.1]$, and $[-0.4\pi,\,0.4\pi]$, respectively. Fig. 2 depicts the average magnitude and angle estimation errors of three PSSE schemes across buses. The curves in Fig. 2 demonstrate the merits of the FPP-based PSSE solver in this scenario.

VI Conclusions

Motivated by the inherent nonconvexity of the power flow and PSSE tasks and leveraging recent advances in handling nonconvex QCQPs, this work first reformulated power flow and PSSE as a nonconvex QCQP. The resulting nonconvex QCQP was subsequently solved by the FPP algorithm. The novel FPP-based solvers were shown to converge to a stationary point of the WLS formulation. To fairly compare different PSSE solvers from noisy data, the CRLB for PSSE assuming an AWGN model was derived based on Wirtinger’s calculus for functions over complex domains. Extensive numerical tests showed markedly improved performance of our FPP-based solver for both power flow and PSSE tasks at the price of increased runtime over competing Gauss-Newton- and SDR-based alternatives on a variety of IEEE test systems.

Pertinent future research directions include developing distributed implementations for large-scale power networks by exploiting the natural low-rank and sparsity structure present in the coefficient matrices $\{\bm{H}_{\ell}\}$. Another possibility consists of leveraging state-of-the-art approaches for tackling random quadratic systems of equations to solve the power flow and PSSE problems [31]. Generalizing feasible point pursuit algorithms to other nonconvex power grid control tasks such as stochastic energy management [32], and distribution system-level power flow and PSSE [13] constitute meaningful directions for future research as well.

VII Appendix

Proof of Proposition 2: For the AWGN model in (14) with $\bm{\eta}\sim\mathcal{N}(\bm{0},\operatorname{diag}(\bm{\sigma^{2}}))$, the data likelihood can be written as

[TABLE]

and the negative log-likelihood $f(\bm{v})=-\ln p(\bm{z};\bm{v})$ is

[TABLE]

The Fisher information matrix is defined as the Hessian of the objective function $f(\bm{v})\in\mathbb{R}$ with respect to the variable vector $\bm{v}\in\mathbb{C}^{N}$. So the task of deriving the Cramér-Rao bound amounts to finding the Hessian of a real-valued function with respect to a complex-valued vector. Recall from Wirtinger’s calculus that $f(\bm{v})$ can be equivalently rewritten as $f(\bm{v},\overline{\bm{v}})$ [33]. Upon introducing the conjugate coordinates $[\bm{v}^{\mathcal{T}}~{}\overline{\bm{v}}^{\mathcal{T}}]^{\mathcal{T}}\in\mathbb{C}^{2N}$, the so-called Wirtinger derivative is [33]

[TABLE]

Our definitions here follow the convention in multivariate calculus that derivatives are denoted by row vectors, and gradients by column vectors. For brevity, let $\phi_{\ell}(\bm{v},\overline{\bm{v}}):=z_{\ell}-\overline{\bm{v}}^{\mathcal{T}}\bm{H}_{\ell}\bm{v}$. Accordingly, the derivatives of $f$ in (27) can be obtained as

[TABLE]

where the partial derivatives of $\phi_{\ell}$ can be found as

[TABLE]

In the conjugate coordinate system, the complex Hessian is defined as

[TABLE]

whose blocks are given by

[TABLE]

After substituting (29) and (30) into the last equations, and with some tedious algebraic manipulations, the first block of ${\boldsymbol{\mathcal{H}}}$ can be obtained as

[TABLE]

The other blocks can be derived in a similar fashion. Upon omitting algebraic details, the remaining three blocks can be obtained as follows

[TABLE]

Evaluating the Hessian ${\boldsymbol{\mathcal{H}}}$ in (31) [and its blocks in (32)-(35)] at the true value of $\bm{v}$, and taking the expectation with respect to the noise vector $\bm{\eta}$, it is easy to verify that $\mathbb{E}\left[\phi_{\ell}(\bm{v},\overline{\bm{v}})\right]=0$. Hence, the $\phi_{\ell}$-related terms disappear, so the FIM $\bm{F}:=\mathbb{E}[{\boldsymbol{\mathcal{H}}}]\in\mathbb{C}^{2N\times 2N}$ can be expressed as [34]

[TABLE]

where $\bm{G}:=[\bm{g}_{1}~{}\cdots~{}\bm{g}_{L}]\in\mathbb{C}^{2N\times L}$ is introduced to show the rank-deficiency of $\bm{F}$, whose $\ell$-th column is given as

[TABLE]

To demonstrate the rank-$1$ deficiency of $\bm{F}$, it suffices to find a nonzero vector $\bm{d}\in\mathbb{C}^{2N}$ such that $\bm{F}\bm{d}=\bm{0}$. To this end, consider the vector $\bm{d}:=\left[\bm{v}^{\mathcal{T}}~{}-\overline{\bm{v}}^{\mathcal{T}}\right]^{\mathcal{T}}\neq\bm{0}$. It is straightforward to check that for all $\ell=1,\,2,\,\ldots,\,L$

[TABLE]

therefore giving rise to $\bm{F}\bm{d}=\sum_{\ell=1}^{L}\bm{g}_{\ell}\left(\bm{g}_{\ell}^{\mathcal{H}}\bm{d}\right)=\bm{0}$. That is, for any nonzero $\bm{v}$, there always exists a nonzero vector $\bm{d}=\left[\bm{v}^{\mathcal{T}}~{}-\overline{\bm{v}}^{\mathcal{T}}\right]^{\mathcal{T}}$ lying in the null space of $\bm{F}$, hence verifying the rank-$1$ deficiency of $\bm{F}$. This concludes the proof.