Phase retrieval is a prevalent problem in digital signal processing and experimental physics that consists of estimating a complex signal from magnitude measurements. This paper expands the classical phase retrieval framework to electric power systems with unknown network models and limited access to observations of voltage magnitudes, active power injections, and reactive power injections. The proposed method recovers the phase angles and the power-phase angle submatrices of the AC power flow Jacobian matrix. This is made possible by deriving topology and parameter-free expressions for the structural symmetries of the power flow Jacobian that do not depend on the phase angles. These physical laws provide structural constraints for the proposed phase retrieval method. The paper then presents sufficient conditions for guaranteed recovery of the voltage phase angles, which also depend solely on voltage magnitudes, active power injections, and reactive power injections. The method offers two significant benefits: both estimating the voltage phase angles and recovering the power flow Jacobian matrix—a basis for approximating the power flow equations. Simulations on widely studied open-source test networks validate the findings.

§ INTRODUCTION

Flows of active and reactive power in AC electric power systems are directly related to the nodal voltage phase angles. Thus, there is great interest in knowing the phase angles for many network planning, control, and protection applications. To this end, broad efforts have taken place to deploy phasor measurement units (PMUs), which use communication networks and protocols to synchronize measurements of phase angles.

Nonetheless, due to their complexity and cost, PMU deployment remains limited in some regions and grid types, such as in transmission system boundaries, transmission systems in rural environments, and distribution systems. Another method to determine the voltage phase angles is through state estimation. However, conventional state estimation techniques require detailed and complete knowledge of the network topology and model parameters. This is rarely possible in practice <cit.>, particularly for combined transmission and distribution systems. In addition, network models are often inaccurate and difficult to maintain.

Therefore, we propose a method that enables the recovery of the voltage phase angles in electric power systems where the network model is unknown or unavailable. We extend previous works on recovering sensitivity matrices that relate the bus complex injections to bus voltage phasors to show that the structure of the AC power flow Jacobian enables the recovery of the power-to-phase angle sensitivities—and consequently, the phase angles—when neither observations of the phase angles nor the network topology are available.

§.§ Problem description

Each bus $i$ in an electricity network is modeled with a complex power injection and voltage phasor, denoted as

\begin{equation}
\label{eq:pmu_signal_def}
v_i \angle{\theta_i} \in \C, \quad
p_i + j q_i \in \C, \quad i=1,\dots,n. \tag{Phasors}%\tag{Synchrophasors}
\end{equation}

We often lack full information about these complex variables, but receive measurements of the form

\begin{equation}
\label{eq:ami_signal_definition}
\quad  v_i \in \R, \quad p_i + jq_i \in \C, \quad i =1,\dots,n, \tag{Magnitudes}
\end{equation}

where $v_i, p_i$,  and $q_i$ are the voltage magnitude, active power injection, and reactive power injection at bus $i$, respectively. The measurements (<ref>) may emerge from substation sensors in the context of transmission systems, or from sensors in “advanced metering infrastructure" (AMI) in the context of distribution systems. Critically, these measurements lack information about the voltage phase angles $\theta_1,\dots,\theta_n$.

Motivated by the contrast between the measurements in (<ref>) and (<ref>), we extend the classical phase retrieval problem to electric power systems with unknown network models, through the lens of a recovered AC power flow Jacobian matrix  <cit.>. Phase retrieval is a rich digital signal processing theory <cit.> that consists of estimating a complex-valued signal $\vx \in \C^n$ given observations of real-valued magnitude measurements $\vb \in \R^m$. In the case of a linear system of equations described by a matrix $\mA \in \C^{m \times n}$, classical phase retrieval is concerned with the problem

\begin{equation}
\label{eq:high-level-phir}
\operatorname{find} \quad \vx \in \C^{n} \quad \text{subject to:} \quad | \mA \vx | = \vb,
\end{equation}

where $\vb \in \R^m$ are the element-wise magnitudes of the measurements $\mA \vx$.

Importantly, we will show that the structure of our particular problem allows us to circumvent full knowledge of the analogous $\mA$—i.e., the power flow Jacobian matrix—by exploiting the symmetries in its structure that arise due to electrical circuit physics. In Section <ref> and in Section <ref>, this structure is reviewed. We derive a novel representation of the symmetries within this structure in Section <ref>, which then leads to the recovery of the  voltage phasors and the full power flow Jacobian in Section <ref>. This representation depends on neither the phase angles nor the network model. Moreover, we then use this representation to develop sufficient conditions for guaranteed voltage phase retrieval and power flow Jacobian invertibility in Section <ref>.

§.§ Contributions

In summary, this work proposes a novel application of phase retrieval—electric power systems—and synthesizes theory from both fields for this application. Specifically, considering only observations of the voltage magnitudes, active power injections, and reactive power injections as in (<ref>), we show that it is possible to recover the voltage phasors as in (<ref>).  This is achieved by deriving and exploiting a novel expression for the physical law describing symmetries in the power flow Jacobian matrix. The proposed method  contributes solutions to several classes of power system estimation problems from only the measurements in (<ref>). Specifically, the contributions of this paper are:

* Estimating the voltage phase angles $\vtheta$ as in (<ref>), or their perturbations $\Delta \vtheta$, from only the measurements in (<ref>), and developing sufficient conditions that guarantee voltage phase retrieval.

* Estimating the power-voltage phase angle sensitivity matrices $\DPDTH, \DQDTH \in \mathbb{R}^{n \times n}$, which are  two of the four blocks of the power flow Jacobian $\mJ(\vx) \in \R^{2n \times 2n}$, which is defined for a given operating point of the network $\vx \triangleq[ \vtheta^T,\vv^T ]^T \in \mathbb{R}^{2n}$.

* Estimating the current injection phase angles from observations of the current magnitudes described by the complex power measurements in (<ref>).

[01]\(j\)The imaginary unit ($j^2 = -1$)
[02]\(A,a\) Scalar
[05]\( \setA \)Set
[06]\(\mathds{1}\)Vector of all ones
[07]\(\mId\)The identity matrix
[08]\(\frac{\partial \vf}{\partial \vx}\)$m \times n$ Jacobian matrix of $\vf: \R^n \mapsto \R^m$
[09]\(\frac{\partial f_i}{\partial x_k}\)$(i,k)$-th entry of the Jacobian matrix $\frac{\partial \vf}{\partial \vx}$
[10]\(\norm{\va}_2\)Euclidean/$\ell_2$ norm of vector $\va$, $\norm{\va}_2 = \left(\sum_{i=1}^n a_i^2 \right)^{\frac{1}{2}}$
[11]\(\norm{\mA}_2\)Spectral norm/2-norm of matrix $\mA$, $\norm{\mA}_2 = \max_{\vx \neq \boldsymbol{0}} \frac{\norm{\mA \vx}_2}{\norm{\vx}_2}$
[12]\(\lambda(\mA)\)Eigenvalues of matrix $\mA$
[13]\(\sigma(\mA)\)Singular values of matrix $\mA$
[14]\(\diag(\vx)\)Diagonal matrix, with entries of vector $\vx$ along the diagonal
[15]\(\trace{\mA}\)Sum of the diagonal entries of matrix $\mA$
[16]\(\rank(\mA)\)Rank of matrix $\mA$
[17]$\lvert \cdot \rvert$Element-wise absolute value
[19]\((\cdot)^*\)Complex conjugate transpose
[20]\((\cdot)^{\dagger}\)Moore-Penrose pseudoinverse
[21]\((\cdot) \circ (\cdot)\)Element-wise multiplication
[22]\(\conj{(\cdot)}\)Complex conjugate of a complex argument
[23]\(\Re{\cdot}\)Real component of a complex argument
[24]\(\Im{\cdot}\)Imaginary component of a complex argument

§ PRELIMINARIES

The benefits of having voltage phase information as in (<ref>) are well-known <cit.>. Knowledge of voltage phase information allows for improved power system model validation, control, and  protection <cit.>. It also enables improvements in grid reliability, planning, and investment decision-making. Some of these benefits include increasing the accuracy of state estimators in distribution <cit.> and transmission <cit.> systems. However, the cost and complexity of PMUs have resulted in limited deployment. Approximately 3000 PMUs have been deployed in the United States as of 2022 <cit.>, up from 2000 PMUs as of 2015 and 500 PMUs as of 2009 <cit.>.

§.§ Phaseless electricity network modeling techniques

Motivated by the heterogeneous deployment of PMUs, a growing literature has developed on phaseless—or “non-PMU"—measurement-based modeling techniques for electricity networks. This literature focuses on power systems where the network model is unknown, and where measurements are limited to $\vv,\vp,\vq$ as in (<ref>).

This literature has developed methods for estimating the unknown topology and electrical parameters—both for single-phase<cit.> and multiphase unbalanced <cit.> networks. The phaseless estimation of the bus impedance matrix <cit.> has recently seen experimental validation in hardware.

Additionally, this literature has developed methods for estimating complex power injections using voltage magnitude measurements. These methods largely center on the voltage-power submatrices of the inverse power flow Jacobian <cit.>, and have also seen recent experimental validation in hardware <cit.>.

Our work contributes to this burgeoning literature. It provides the first set of methods and sufficient conditions for the recovering the voltage phase angles as in (<ref>) given only measurements of the form (<ref>)—even when the network model is unknown.

§.§ Review of the phase retrieval problem

Optimization problems of the form (<ref>) are known as phase retrieval problems. Classically, the goal of phase retrieval is to recover the phase of a complex signal $\vx \in \C^n$ given a design matrix $\mA \in \C^{m \times n}$ and real-valued magnitude measurements  $|\mA \vx | = \vb$, where $\vb \in \R^{m}$, and $|\cdot|$ denotes the elementwise magnitudes of a vector. In this paper, we make the noteworthy distinction that the design matrix (the power flow Jacobian) is not assumed to be known, in contrast with the classical approach to phase retrieval.

One of the earliest investigations of this problem was in experimental physics, where <cit.> developed an iterative method for recovering the phase of a wave function from the magnitude of its Fourier Transform coefficients, provided by intensity measurements. The authors of <cit.> note wide applications for their method in X-ray crystallography, which has been discussed more recently in <cit.>. Since then, the phase retrieval problem and its various formulations have seen applications in electron microscopy <cit.>; optical and medical imaging <cit.>; experimental physics <cit.>; and digital signal processing. We refer the interested reader to <cit.> and <cit.> for multidisciplinary surveys of the formulations and applications of phase retrieval.

§.§.§ Classical phase retrieval

The classical phase retrieval problem, where $\mA$ is assumed to be known, (in our case, it will be recovered) can be understood as a quadratic optimization problem with complex decision variables:

\begin{equation}
\label{eq:basic_ph_retr_prob}
\minimize_{\vx \in \C^n, \vy \in \C^m} \quad \norm{\mA \vx - \vy}_2^2
\quad \text{subject to:}
\quad |\vy| = \vb.
\end{equation}

The introduction of the decision variables $\vy \in \C^m$ in (<ref>) is sometimes referred to as “lifting" the magnitude measurements $\vb \in \R^m$ to $\C^m$, which allows us to recover the complex vector $\mA \vx \in \C^m$. It is noted in <cit.> that the problem (<ref>) can be recast as solving for the complex signal $\vx \in \C^n$ and for the phases of the magnitude measurements $\vb$. The phases are modeled as rectangular complex numbers $\vu \in \C^m$ that are constrained to lie on the unit circle, i.e., $|u_i| = 1 \ \ \forall i=1,\dots,n$. With this framework, the phase retrieval program (<ref>) can be rewritten as

\begin{equation}
\label{eq:basic-ph-retr-prob-fixed-phase}
\minimize_{\vx \in \C^n, \vu \in \C^m} \quad \norm{\mA \vx - \diag(\vb) \vu}^2_2 \quad \text{subject to:} \quad |\vu| = \mathds{1}.
\end{equation}

The authors of <cit.> note that, for any fixed guess of the phase $\hat{\vu} \in \C^m$ of the magnitude measurements $\vb$, an estimate $\hat{\vx} \in \C^n$ of the complex signal is unique and is given as

\begin{equation}
\hat{\vx} = \mA^{\dagger} \diag(\vb) \hat{\vu} = (\mA^* \mA)^{-1} \mA^* \diag(\vb) \hat{\vu},
\end{equation}

where $(\cdot)^*$ is the complex conjugate transpose and $(\cdot)^{\dagger}$ is the Moore-Penrose pseudoinverse.

With this in hand, we can substitute $\vx$ with $\mA^{\dagger} \diag(\vb) \vu$, and solve the quadratic program

\begin{equation}
\label{eq:generalized-high-level-phir}
\minimize_{\vu \in \C^m : |\vu| = \mathds{1}} \vu^* \mM \vu, \ \ \text{subject to:} \ \   \mM = \diag(\vb)(\mId - \mA \mA^{\dagger})\diag(\vb),
\end{equation}

where $\mM \succeq 0$. However, the modulus constraints in (<ref>) are computationally challenging. They are equivalent to the quadratic equality constraints $\Re{u_i}^2 + \Im{u_i}^2 =1$ for all $i$, which produces a non-convex feasible region.

§.§.§ The PhaseCut algorithm

This paper is inspired by a well-known approach to the phase retrieval problem known as PhaseCut. The algorithm, developed in <cit.>, casts the phase retrieval problem in terms of a semidefinite programming (SDP) relaxation of the maximum cut problem, which finds the maximum number of edges in a network that form a partition of the network into two disjoint subsets. The purpose of PhaseCut essentially boils down to the fact that the intractable modulus constraints in (<ref>) can be avoided via a relaxed SDP. This is achieved by modeling the phase $\vu \in \C^m$ of the signal $|\mA \vx |$ as a rank-one matrix $\mX \triangleq \vu \vu^*$. With this matrix decision variable, (<ref>) can be written as

\begin{equation}
\label{eq:semidefinite_relaxed_phase_retrieval}
\minimize_{\mX \succeq 0 } \trace{ \mX \mM }, \quad \text{subject to:} \quad \rank(\mX) = 1, \quad \diag(\mX) = \mathds{1},
\end{equation}

which is relaxed by dropping the $\rank(\mX)=1$ constraint. Although we do not explicitly use this SDP relaxation in this work, its precursor quadratic formulation, (<ref>) is central to the forthcoming results.

§.§ Review of Newton-Raphson power flow

§.§.§ Power flow equations

We consider an arbitrary electric power network with $n$ “PQ" buses—i.e., buses where the complex power injections are specified. The well-known power flow equations in polar form for buses $i =1,\dots,n$ are

\begin{align}
\label{eq:active_power_balance}
p_i &= v_i \sum_{k=1}^{n}  v_k (G_{ik} \cos\theta_{ik} + B_{ik}\sin\theta_{ik}) = p^g_i - p^d_i, \\
\label{eq:reactive_power_balance}
q_i &= v_i \sum_{k=1}^{n}  v_k (G_{ik} \sin\theta_{ik} - B_{ik}\cos \theta_{ik}) = q^g_i - q^d_i,
\end{align}

where $\theta_{ik}\triangleq\theta_i -\theta_k$ is the phase angle differences between buses $i,k$, $p^g_i,q^g_i$ are the generated active and reactive powers at bus $i$, $p^d_i,q^d_i$ are the demanded active and reactive powers  at bus $i$, $v_i,v_k$ are the voltage magnitudes at buses $i,k$, and $G_{ik},B_{ik}\in \mathbb{R}$ are the real and imaginary components of the $(i,k)$-th entry $Y_{ik} \triangleq G_{ik} + j B_{ik} \in \C$ of the nodal admittance matrix $\mY = \mG + j \mB \in \C^{n \times n}$.

By convention, the state of the system is defined as $\vx \triangleq [\vtheta^T,\vv^T]^T \in \R^{2n}$, where $\vv \in \R^n$ are the voltage magnitudes and $\vtheta \in (-\pi,\pi]^n$ are the voltage phase angles. The complex voltage at each bus takes the form $\overline{v}_i \triangleq v_i e^{j \theta_i} = v_i(\cos(\theta_i) + j \sin(\theta_i)) \in \C$. The bus active and reactive power injection vector $\vg(\cdot): \R^{2n}\mapsto \R^{2n}$, i.e., the vectorization of the power flow equations (<ref>), is defined as $\vg(\vx) \triangleq [\vp(\vx)^T,\vq(\vx)^T]^T \in \R^{2n}$. The solution $\vx^*$ of the non-linear system of equations (<ref>) satisfies $\vg(\vx^*) - \left[(\vp^g - \vp^d )^T,(\vq^g - \vq^d)^T\right]^T = \boldsymbol{0}.$

§.§.§ Newton-Raphson power flow in polar coordinates

The Newton-Raphson power flow algorithm is a well-studied iterative approach to solving the non-linear power flow equations (<ref>). This method consists of repeatedly solving the linear system of equations

\begin{equation}
\renewcommand*{\arraystretch}{1.5}
\label{eq:nr_pf_definition}
\begin{bmatrix}
\Delta \vp\\
\Delta \vq
\end{bmatrix} =
\begin{bmatrix}
\DPDTH(\vx) & \DPDV(\vx)\\
\DQDTH(\vx) & \DQDV(\vx)
\end{bmatrix}
\begin{bmatrix}
\Delta \vtheta\\
\Delta \vv
\end{bmatrix} = \mJ(\vx) \Delta \vx,
\end{equation}

where—around a specified operating point—$\Delta \vp,\Delta \vq\in \R^n$ are small perturbations in the active and reactive power injections respectively, and $\Delta \vtheta,\Delta \vv \in \R^n$ are small perturbations in the voltage phase angles and magnitudes, respectively. The matrix $\mJ(\vx) \in \R^{2n \times 2n}$ in (<ref>) is the power flow Jacobian, which contains the partial derivatives of the power flow equations (<ref>) with respect to the voltage magnitudes $\vv$ and phase angles $\vtheta$.

In this work, we consider an estimate of the power flow Jacobian $\mJ(\vx_t)$ computed at a time point $t$. We use $\mJ(\vx_t)$ to provide an approximation <cit.> of the AC power flow equations (<ref>) around the operating point measured at a time point $t$ as

\begin{equation}
\label{eq:first-order-linearization}
\vg(\vx_{t+1}) \approx \vg(\vx_{t}) + \mJ(\vx_t)(\vx_{t+1} - \vx_{t}).
\end{equation}

Here, the power flow Jacobian $\mJ(\vx_t)$ denotes the sensitivity of the injections to changes in the state from the operating point measured at time $t$. We consider the notion of observed or predicted perturbations in the injections and state around the operating point measured at time point $t$. Denote these perturbations as $\Delta \vg(\vx_t) \triangleq \vg(\vx_{t+1}) - \vg(\vx_t) = [\vp_{t+1}^T,\vq_{t+1}^T]^T - [\vp_{t}^T,\vq_{t}^T]^T$ and $\Delta \vx_t \triangleq \vx_{t+1} - \vx_t = [\vtheta_{t+1}^T,\vv_{t+1}^T]^T - [\vtheta_{t}^T,\vv_{t}^T ]^T$.

In this paper, the network model is unknown, and PMU measurements as in (<ref>) are entirely unavailable. Thus, we can neither measure nor solve for the phase angles, and the quantities $\Delta \vtheta_t$, $\DPDTH(\vx_t)$, and $\DQDTH(\vx_t)$ are unknown. We will assume that the power-voltage magnitude sensitivities $\DPDV(\vx_t)$ and $\DQDV(\vx_t)$ are known in the sense that they can be estimated from prior measurements, which is a well-studied problem (e.g., see <cit.>).

In this paper, we ask the following question: can we recover the phase angles $\vtheta_t$ as in (<ref>) and the power flow Jacobian submatrices $\DPDTH(\vx_t), \DQDTH(\vx_t)$ as in (<ref>), using only the measurements $\vv_t,\vp_t,\vq_t$ as in (<ref>)?  This question is illustrated in Fig. <ref>.

Illustration of voltage phase retrieval showing model-free power flow Jacobian reconstruction technique

§.§ The structure of the power flow Jacobian

§.§.§ Power-voltage magnitude sensitivity matrices

Forming the entries of the block submatrices of the power flow Jacobian requires that we differentiate the power flow mismatch equations (<ref>) with respect to the voltage magnitudes and phase angles at each bus. First, consider the derivatives of the net active power and reactive power injections at bus $i$ with respect to the voltage magnitudes at bus $k$:

\begin{align}
\label{eq:dpdv_jac_block}
\frac{\partial p_i}{\partial v_k} =
\begin{dcases}
v_i (G_{ik}\cos\theta_{ik} + B_{ik} \sin\theta_{ik} ) & i \neq k,\\
2 v_i G_{ii}  +  \sum_{k: k\neq i} v_k(G_{ik} \cos\theta_{ik} + B_{ik}\sin\theta_{ik})  & i = k.
\end{dcases}
\\
\label{eq:dqdv_jac_block}
\frac{\partial q_i}{\partial v_k} =
\begin{dcases}
v_i (G_{ik}\sin\theta_{ik} - B_{ik} \cos\theta_{ik} ) & i \neq k,\\
-2v_iB_{ii} + \sum_{k : k\neq i} v_k(G_{ik} \sin\theta_{ik} - B_{ik}\cos\theta_{ik}) & i = k.
\end{dcases}
\end{align}

The diagonal terms can be alternatively expressed as $\frac{\partial p_i}{\partial v_i} = \frac{p_i(\vx)}{v_i}+ G_{ii} v_i $ and $\frac{\partial q_i}{\partial v_i} = \frac{q_i(\vx)}{v_i} - B_{ii}v_i$.

§.§.§ Power-voltage phase angle sensitivity matrices

Similarly to the power-voltage magnitude sensitivities in Section <ref>, by differentiating the net active and reactive power injections with respect to the voltage phase angles, we obtain

\begin{align}
\label{eq:dpdth_jac_block}
\frac{\partial p_i}{\partial \theta_k} =
\begin{dcases}
v_i v_k (G_{ik} \sin\theta_{ik} - B_{ik}\cos\theta_{ik}) & i \neq k,\\
-v_i \sum_{k: k\neq i} v_k (G_{ik}\sin\theta_{ik} - B_{ik} \cos \theta_{ik}) & i = k,
\end{dcases}\\
\label{eq:dqdth_jac_block}
\frac{\partial q_i}{\partial \theta_k} =
\begin{dcases}
-v_i v_k(G_{ik}\cos \theta_{ik} + B_{ik}\sin\theta_{ik}) & i \neq k,\\
v_i \sum_{k: k\neq i} v_k (G_{ik} \cos\theta_{ik} + B_{ik} \sin \theta_{ik}) & i = k.
\end{dcases}
\end{align}

The diagonal terms can alternatively be expressed as $\frac{\partial p_i}{\partial \theta_i}= -q_i(\vx) - B_{ii}v_i^2 $ and $\frac{\partial q_i}{\partial \theta_i} = p_i(\vx) - G_{ii}v_i^2 $.

§ ANALYTICAL RESULTS

§.§ Model-free recovery of the power flow Jacobian

There are well-known mathematical symmetries that relate the power flow Jacobian block submatrices. These symmetries arise due to the physical laws of the power flow equations <cit.>. The canonical expressions  for these symmetries <cit.> depend on the network topology and electrical parameters. In this section, we propose a topology and parameter-free extension of these known symmetries that:

* Has a valuable engineering implication, supporting system identification without voltage phase angles.

* Has a valuable theoretical implication, by providing structural constraints for applying the phase retrieval framework—informed by the physics of Newton-Raphson power flow.

Specifically, we show that we can represent $\DPDTH$ and $\DQDTH$ as a function of the power-to-voltage magnitude submatrices of the power flow Jacobian $\DPDV$ and $\DQDV$, the voltage magnitudes $\boldsymbol{v}$, and the active and reactive power injections $\boldsymbol{p},\boldsymbol{q}$. Representing $\DPDTH$ and $\DQDTH$ in this way is valuable because:

* It removes dependence on direct measurements of the voltage phase angles, as in (<ref>).

* It requires only measurements as in (<ref>), and estimates for the power-voltage magnitude sensitivities; which is a well-studied problem <cit.>.

This result is achieved through the following Lemma:

Consider an electric power network with PQ buses $\{1,\dots,n\}$. Let $\vx = [\vtheta^T,\vv^T]^T$ be a guess for the network state, producing power injections $\vp(\vx),\vq(\vx)$. For all buses $i,k \in \{1,\dots,n\}$, the entries $\frac{\partial p_i}{\partial \theta_k}$ and $\frac{\partial q_i}{\partial \theta_k}$ of the power-voltage phase angle sensitivity matrices, $\DPDTH(\vx),\DQDTH(\vx) \in \mathbb{R}^{n \times n}$, as defined in (<ref>), can be written as

\begin{align}
\label{eq:thm1_dpdth}
\frac{\partial p_i}{\partial \theta_k} = \begin{dcases}
v_k \frac{\partial q_i}{\partial v_k} & i \neq k\\
v_i \frac{\partial q_i}{\partial v_i} - 2 q_i  & i = k,
\end{dcases}
\\
\label{eq:thm1_dqdth}
\frac{\partial q_i}{\partial \theta_k} = \begin{dcases}
-v_k \frac{\partial p_i}{\partial v_k} & i \neq k\\
-v_i \frac{\partial p_i}{\partial v_i} + 2p_i & i = k.
\end{dcases}
\end{align}

The relationship for the off-diagonal terms is well-known <cit.>, and is found by substituting (<ref>) and (<ref>) into (<ref>) and (<ref>), respectively:

\begin{align}
\frac{\partial p_i}{\partial \theta_k} =
\left(
v_i(G_{ik} \sin \theta_{ik} - B_{ik} \cos \theta_{ik}
\right)
}_{\triangleq \frac{\partial q_i}{\partial v_k}}, \quad i \neq k,\\
\frac{\partial q_i}{\partial \theta_k} =
\left(
v_i(G_{ik} \cos \theta_{ik} + B_{ik} \sin \theta_{ik})
\right)
}_{\triangleq \frac{\partial p_i}{\partial v_k}}, \quad  i\neq k.
\end{align}

The topology-free expressions for the diagonal terms $\frac{\partial p_i}{\partial \theta_i}$ and $\frac{\partial q_i}{\partial \theta_i}$ in (<ref>) and (<ref>) are found by manipulating the power flow equations (<ref>) and the definition of the power-voltage magnitude sensitivities (<ref>). To show this, first consider the expression (<ref>). We want to show that this holds for the case when $i=k$. Substituting the definition of the diagonal reactive power-voltage magnitude sensitivities $\frac{\partial q_i}{\partial v_i}$, (<ref>), and the reactive power balance equation, (<ref>), suppose that we can write $\frac{\partial p_i}{\partial \theta_i}$ for $i=1,\dots,n$ as

\begin{equation}
\begin{split}
%\frac{\partial p_i}{\partial \theta_i} =
v_i\frac{\partial q_i}{\partial v_i}  - 2 q_i =
-2v_i^2 B_{ii}  + v_i \sum_{k : k\neq i} v_k(G_{ik} \sin\theta_{ik} - B_{ik}\cos\theta_{ik})
- 2v_i \sum_{k=1}^{n}  v_k (G_{ik} \sin\theta_{ik} - B_{ik}\cos \theta_{ik}).
\end{split}
\end{equation}

Next, pull out the $i$th term from the argument of the sum over $k=1$ to $n$,

\begin{equation}
% v_i(G_{ii}  \underbrace{\sin\theta_{ii}}_{=0} - B_{ii}\underbrace{\cos\theta_{ii}}_{=1}),
%no underbrace
v_i(G_{ii}  \sin\theta_{ii} - B_{ii}\cos\theta_{ii}) = -v_iB_{ii},
\end{equation}

which allows us to express $\frac{\partial p_i}{\partial \theta_i}=   v_i \frac{\partial q_i}{\partial v_i} - 2q_i$ as

\begin{equation}
\label{eq:dpdth_proof_penultimate_step}
\begin{split}
-2v_i^2 B_{ii}  &+ v_i
\sum_{k: k \neq i} v_k(G_{ik} \sin\theta_{ik} - B_{ik}\cos\theta_{ik}) %\dots
- 2v_i \sum_{k: k \neq i}  v_k (G_{ik} \sin\theta_{ik} - B_{ik}\cos \theta_{ik}) +2 v_i^2 B_{ii}.
\end{split}
\end{equation}

Simplifying (<ref>), we obtain (<ref>) for $i=k$,

\begin{align}
%\frac{\partial p_i}{\partial \theta_i} =
- v_i \sum_{k: k \neq i}  v_k (G_{ik} \sin\theta_{ik} - B_{ik}\cos \theta_{ik}) \triangleq \frac{\partial p_i}{\partial \theta_i},
\end{align}

as desired, which verifies that the expression (<ref>) is true for $i=k$.

Analogously, we can see that the expression (<ref>) for the reactive power-voltage phase angle sensitivity matrix is true by substituting in the expression (<ref>) for $\frac{\partial p_i}{\partial v_i}$, and the active power balance equation (<ref>). Then, applying the same logic as before, suppose that we can write $ \frac{\partial q_i}{\partial \theta_i}$ as:

\begin{align}
\frac{\partial q_i}{\partial \theta_i}  = \ -v_i \frac{\partial p_i}{\partial v_i} &+ 2p_i,\\
\begin{split}
= -2 v_i^2 G_{ii}  &- v_i \sum_{k : k\neq i} v_k(G_{ik} \cos\theta_{ik} + B_{ik}\sin\theta_{ik})
+2v_i \sum_{k=1}^{n}  v_k (G_{ik} \cos\theta_{ik} + B_{ik}\sin\theta_{ik}),
\end{split}\\
\begin{split}
= -2 v_i^2 G_{ii}  &- v_i \sum_{k : k\neq i} v_k(G_{ik} \cos\theta_{ik} + B_{ik}\sin\theta_{ik})
+2v_i \sum_{k: k \neq i}  v_k (G_{ik} \cos\theta_{ik} + B_{ik}\sin\theta_{ik}) + 2v_i^2G_{ii}.
\end{split}
\end{align}

Simplifying (<ref>) yields (<ref>) for $i=k$ as

\begin{equation}
v_i \sum_{k: k \neq i}  v_k (G_{ik} \cos\theta_{ik} + B_{ik}\sin\theta_{ik}) \triangleq \frac{\partial q_i}{\partial \theta_i},
\end{equation}

which verifies the desired expression (<ref>) for $i=k$.

Notice that by Lemma <ref> the power-voltage phase angle sensitivity matrices can be expressed as functions $\DPDTH,\DQDTH : \R^n \times \R^n \mapsto \R^{n \times n}$ of the form

(,) = () - 2 (),

(,) = -()  + 2 (),

which are implicitly parameterized by $\DQDV$ and $\DPDV$.

Furthermore, the matrices (<ref>) are functions that depend solely on (<ref>), and depend on neither the network topology nor the voltage phase angles. The matrices $\DPDV,\DQDV$ can also be efficiently estimated from (<ref>) using well-studied methods <cit.>.

Computationally, we can verify expressions (<ref>) and (<ref>)  by considering a simple 2-bus network, shown in the “flat start" configuration at the left of Fig. <ref>. We solve the system (<ref>) for the state $\vx_k = [\theta_2^{(k)}, v_2^{(k)}]^T$ iteratively for $k=1,\dots$ using the usual expressions (<ref>) and the phaseless expressions (<ref>) for $\frac{\partial p_2}{\partial \theta_2}$ and $\frac{\partial q_2}{\partial \theta_2}$. The equivalence of the mismatches yielded by this procedure is illustrated on the right of Fig. <ref>. Further numerical evidence for the exactness of the expressions (<ref>) in a large-scale test network is available in Appendix <ref>.

Equivalence of the Newton-Raphson iterations (<ref>) using the $\theta$-free expressions (<ref>) and standard expressions (<ref>) for $\DPDTH,\DQDTH$ (right) for a simple two bus test case (left).

§.§ Phase retrieval via power flow Jacobian recovery

We propose a specialized variant of the quadratic phase retrieval program (<ref>) for electric power networks. The program leverages the structural symmetries of the power flow Jacobian described in Section <ref>.

§.§.§ Formulating the objective

Let $\vvbar_t \in \C^n$ be observed voltage phasors as in (<ref>), where $\overline{v}_{t}^{(i)} = v_{t}^{(i)} e^{j \theta_t^{(i)}} \triangleq v_t^{(i)} \angle \theta_t^{(i)} \in \C $ denotes the voltage phasor at bus $i \in \{1,\dots,n\}$. Suppose that we observe measurements $\vp_{t},\vq_{t} \in \R^n$ of active and reactive power injections, respectively; and of voltage magnitudes  $\vv_{t} \triangleq |\vvbar_{t}| \in \R^n$, as in (<ref>) and Fig. <ref>. We want to recover the phase angle perturbations $\Delta \vtheta_t \in (-\pi,\pi]^n$ using the observed perturbations $\Delta \vp_t$, $\Delta \vq_t$, and $\Delta \vv_t$. These are used in the linear approximation (<ref>), taken around the network operating point.

Notice that we can equivalently represent a voltage phase angle $\theta_i$ as a complex number $u_i \in \C$. This is done by requiring that $u_i$ lies on the unit circle, and that its argument is the voltage phase angle $\theta_i$, i.e.,

\begin{equation}
|u_i|=1, \quad \textcolor{black}{ \theta_i = \operatorname{atan2}\left( \Im{u_i},\Re{u_i}\right) \quad i=1,\dots,n,}
\end{equation}

where $\operatorname{atan2}(\Im{u_i},\Re{u_i})$ is the argument, or phase angle[
The phase $\theta = \operatorname{atan2}(y,x)$ of $x + jy \in \C$ is defined as $\theta= \arctan(\frac{y}{x})\mathds{1}[x \neq 0] + (1-2\mathds{1}[y<0])(\pi \mathds{1}[x <0] + \frac{1}{2}\pi \mathds{1}[x =0])$, and is undefined if $x=y=0$. The operator $\mathds{1}[\cdot] =1$ if $[\cdot]$ is true, and $\mathds{1}[\cdot] =0$ otherwise.
of $u_i \in \C$.
Drawing inspiration from the classical phase retrieval algorithm (<ref>), at each iteration $t$ we want to find the voltage phase angles $\vtheta_t$ that satisfy Euler's formula for a complex number on the unit circle $\vu_t$, i.e., $\vvbar_t = \diag(\vv_t) \vu_t(\vtheta_t),$ where

\begin{equation}
u_t^{(i)}\left(\theta_t^{(i)}\right) =  \cos\theta_t^{(i)} +j \sin\theta_t^{(i)}, \quad \left|u_t^{(i)}\right|=1, \ i=1,\dots,n.
\end{equation}

We will use this representation to directly apply the logic of the classical phase retrieval problem

\begin{equation}
\operatorname{find} \quad  \Delta \vtheta_t \in \R^n
\quad \text{subject to:} \quad \big| \diag(\Delta \vv_t) \vu_t(\Delta \vtheta_t) \big| = \Delta \vv_t,
%\quad \textsf{st:} \quad |\mJ(\vx_t)^{-1}\vg(\vx_t)| = \Delta \vv_t,%|\mJ(\vx_k) [\Delta \vtheta_k,\Delta \vv_k]^T| = \vg(\vx_k).
\end{equation}

which consists of solving for a given iteration/sample $t$:

\begin{equation}
\label{eq:generic-nr-phase-recovery}
\minimize_{\Delta \vtheta_t} \quad \norm{\Delta \vg(\vx_t) - \mJ(\vx_t)\Delta \vx_t}_2^2, \quad \text{subject to:} \quad \big| \diag(\Delta \vv_t) \vu_t(\Delta \vtheta_t) \big| = \Delta \vv_t.
\end{equation}

where $\Delta \vg(\vx_t) \triangleq \big[ \Delta \vp_t^T,\Delta \vq_t^T\big]^T$ denotes the observed perturbations $\big[ \Delta \vp_t^T,\Delta \vq_t^T\big]^T = \big[\vp_{t+1}^T, \vq_{t+1}^T\big]^T - \big[\vp_{t}^T, \vq_{t}^T\big]^T$ and $\Delta \vx_t \triangleq \big[ \Delta \vtheta_t^T,\Delta \vv_t^T\big]^T$ denotes the grid state change $ \big[ \Delta \vtheta_t^T,\Delta \vv_t^T\big]^T = \big[\vtheta_{t+1}^T,\vv_{t+1}^T\big]^T - \big[\vtheta_{t}^T,\vv_{t}^T \big]^T$.

§.§.§ Phase retrieval by exploiting the structure of the power flow Jacobian

Importantly, the formulation (<ref>) assumes knowledge of the power-phase angle sensitivities $\DPDTH(\vx_t),\DQDTH(\vx_t)$, which is impossible without knowledge of the network model or measurements of the form (<ref>). To resolve this, we can appeal to
Lemma <ref>. This allows us to reconstruct the full power flow Jacobian without knowledge of the voltage phase angles by considering the matrix decision variables $\DPDTH(\vv_t,\vq_t),\DQDTH(\vv_t,\vp_t) \in \R^{n \times n}$ in the phase retrieval program as in (<ref>)

The structural symmetries of the power flow Jacobian (<ref>) can be understood as variants of the structural constraints and the real-valued phase constraints described in <cit.> and <cit.>, respectively.  The phase retrieval program for samples $t=1,\dots,$ can then be written as

\begin{equation}
\label{eq:nr-phret-program-with-symmetry-constraints}
\begin{split}
\renewcommand*{\arraystretch}{1.5}
\minimize_{\Delta \vtheta_t,
\DPDTH,\DQDTH}   \quad
\norm{
\begin{bmatrix}
\Delta \vp_t\\
\Delta \vq_t
\end{bmatrix}
\begin{bmatrix}
\DPDTH(\vv_t,\vq_t) & \DPDV\\
\DQDTH(\vv_t,\vp_t) & \DQDV
\end{bmatrix}
\begin{bmatrix}
\Delta \vtheta_t\\
\Delta \vv_t
\end{bmatrix}
}_2^2 \quad
\text{subject to:} \quad \eqref{eq:augmented-phase-sensitivities}, %\quad | \diag(\Delta \vv_t) \vu_t(\Delta \vtheta_t)| = \Delta \vv_t.
\end{split}
\end{equation}

where, as $\Delta \vtheta_t \in (-\pi,\pi]^n$ by construction, we have dropped the modulus constraints. Modern solvers can solve the program (<ref>)  efficiently. Note that if $\DPDTH(\vv_t,\vq_t)$ is known and invertible, there is a unique solution for the phase angles. Dropping the arguments for simplicity, we have

\begin{equation}
\label{eq:delta-theta-recovery-expression}
\Delta \vtheta_t = \left(\DPDTH\right)^{-1}\left( \Delta \vp_t - \DPDV \Delta \vv_t \right),
%-\left(\DPDTH\right)^{-1} \DPDV \Delta \vv_t + \left(\DPDTH\right)^{-1} \Delta \vp_t,
\end{equation}

because $ \Delta \vp_t = \DPDTH \Delta \vtheta_t + \DPDV \Delta \vv_t$, so $\DPDTH \Delta \vtheta_t = - \DPDV \Delta \vv_t + \Delta \vp_t$, and we have $\Delta \vtheta_t = -\DPDTH^{-1} \DPDV \Delta \vv_t + \DPDTH^{-1} \Delta \vp_t$ as desired. Similarly, if $\DQDTH(\vv_t,\vp_t)$ is invertible, by applying analogous steps for $\Delta \vq_t = \DQDTH \Delta \vtheta_t + \DQDV \Delta \vv_t$, we have

\begin{equation}
\label{eq:delta-theta-reactive-recovery-expression}
\Delta \vtheta_t = \left( \DQDTH\right)^{-1}\left(  \Delta \vq_t -\DQDV \Delta \vv_t \right).
\end{equation}

The expressions in (<ref>) and the program (<ref>) enable the phase angle submatrices of the power flow Jacobian—and thus, the phase angle behavior itself—to be recovered using just the data in (<ref>). Furthermore, we will demonstrate the effectiveness of the  phase retrieval program (<ref>) using estimated power-voltage magnitude sensitivities $\DPDV$ and $\DQDV$ in the numerical results in Section <ref>.

§.§.§ Recovering the current phase angles

For any guess for the voltage phase $\hat{\vtheta} \in (-\pi,\pi]^n$, let $\vvbar \triangleq \vv \angle \hat{\vtheta} \in \C^n$ denote the estimated complex voltage phasors. The complex current injection $\vell \in \C^n$ is related to $\vvbar$ and the complex power injection  $ \vs \triangleq \vp + j \vq \in \C^n$ through Ohm's law as $\vs = \vvbar \circ \conj{\vell}$, where $(\cdot)\circ(\cdot)$ denotes the elementwise multiplication of two vectors, and where $\conj{(\cdot)}$ denotes the elementwise complex conjugate.

Once the voltage phase angles are recovered, all other quantities in the above relationship are given by the measurements in (<ref>). Therefore, we can compute an estimate for the complex current injections as
$\vell = \conj{\diag\left(\vvbar\right)^{-1} \vs}$. While we do not focus our estimation results on the current injections, we verify that this formulation is physically valid in Appendix <ref>.

§.§ Recovery guarantees based on Gershgorin discs

The Gershgorin Circle Theorem provides regions of the complex plane where the eigenvalues of a matrix are guaranteed to lie. It has been applied in voltage stability assessment for siting and control of distributed energy resources  <cit.>, in an impedance-based stability requirement <cit.>, in small signal stability of swing equations <cit.>, and many other electric power systems applications.

The motivations of using the Gershgorin Circle Theorem  to provide phase retrieval guarantees are twofold. First, the Theorem's  structure is intuitively similar to the power flow Jacobian symmetry expression (<ref>). Second, the power flow Jacobian submatrices are empirically diagonally dominant, as we show numerically in Fig. <ref> of the subsequent section. This property is favorable for applying the Gershgorin Circle Theorem.

Below, in Section <ref>, we apply Gershgorin to yield  sufficient conditions for the guaranteed recovery of the voltage phase angles. In Section <ref>, we apply a block matrix variant of Gershgorin to yield a sufficient condition for invertibility of the power flow Jacobian. The sufficient conditions of both Theorems provide these guarantees without requiring any knowledge of the voltage phase angles $\vtheta$.

§.§.§ Sufficient condition for guaranteed voltage phase retrieval

Let $\setB \subset \{1,\dots,n\}$ be a subset of PQ buses under study. Let $\vx = [\vtheta^T, \vv^T]^T \in \R^{2|\setB|}$ be a guess for the network state, let $\vp(\vx),\vq(\vx) \in \R^{|\setB|}$ be the observed power injections, and let $\DPDV(\vx),\DQDV(\vx) \in \R^{|\setB| \times |\setB|}$ be the power-voltage magnitude sensitivity matrices. For any perturbations in (<ref>), if for all $i \in \setB$

\begin{align}
\label{eq:thm1_row_offdiag}
|q_i| &> \frac{1}{2} v_i \left( \sum_{k \in \setB \setminus\{i\}} \left| \frac{\partial q_k}{\partial v_i}\right| - \left| \frac{\partial q_i}{\partial v_i}\right|\right)\\
\label{eq:thm1_column_offdiag} \textsf{or} \quad
|q_i| &> \frac{1}{2}\left(\sum_{k \in \setB \setminus \{i\}} v_k \left| \frac{\partial q_i}{\partial v_k}\right| - v_i \left| \frac{\partial q_i}{\partial v_i} \right| \right),
\end{align}

then there exists unique voltage phase angle perturbations $\Delta \vtheta$ such that (<ref>) holds. Analogously, if for all $i \in \setB$

\begin{align}
|p_i| &> \frac{1}{2}v_i \left( \sum_{k \in \setB \setminus\{i\}} \left| \frac{\partial p_k}{\partial v_i}\right| - \left| \frac{\partial p_i}{\partial v_i} \right|\right)\\
\textsf{or} \quad |p_i| &> \frac{1}{2} \left(\sum_{k \in \setB \setminus\{i\}} v_k \left| \frac{\partial p_i}{\partial v_k}\right| - v_i \left| \frac{\partial p_i}{\partial v_i}\right| \right),
\end{align}

then there exists unique voltage phase angle perturbations $\Delta \vtheta$ such that (<ref>) holds.

The proof of Theorem <ref> can be found in Appendix <ref>.

Theorem <ref> is valuable because it depends solely on measurements of the form (<ref>), and the matrices $\DPDV,\DQDV$, which can be efficiently estimated from (<ref>) using well-studied methods. Theorem <ref> gives explicit bounds on the magnitude of the active or reactive power injection at each bus in a chosen set of buses $\setB$, such that the phase angles at these buses can be guaranteed to be uniquely recovered.

Furthermore, the bounds are each functions of only the voltage magnitudes and one of either the active or reactive power injections, respectively.  If the conditions in Theorem <ref> hold, then $\Delta \vtheta$ can be uniquely recovered via either the expression in (<ref>) or (<ref>), which are equivalent by construction.

§.§.§ Sufficient condition for power flow Jacobian invertibility and guaranteed voltage phase retrieval

We can establish conditions for the invertibility of the full power flow Jacobian without knowledge of the phase angles by appealing to an analogue of the Gershgorin Circle Theorem for block matrices <cit.>. To yield this result, we introduce the following assumption:

The matrices $\DPDTH(\vv,\vq)$ and $\DQDV$ are full rank.

Assumption <ref> is not restrictive, as it only requires that the diagonal blocks of the Jacobian are invertible, without certification of the invertibility of the entire Jacobian. Moreover, we show in Table <ref> of Appendix <ref> that the assumption holds for all the test cases we study, including very large-scale, realistic test cases. With Assumption <ref> in hand, we can use the matrices (<ref>) to yield the following condition, which guarantees the non-singularity of the full Jacobian given the measurements in (<ref>).

Let $\setB \subset \{1,\dots,n\}$ be a subset of PQ buses under study with a guess for the network state $\vx$ and injections $\vp(\vx),\vq(\vx)$ as in Theorem <ref>. Consider the voltage phase angle sensitivity matrix functions $\DPDTH(\vv,\vq)$ and $\DQDTH(\vv,\vp)$ as defined in (<ref>). Suppose that Assumption <ref> holds. Then, if

\begin{align}
\label{eq:thm2_row_offdiag}
\norm{\left(\DPDTH(\vv,\vq)\right)^{-1} \DPDV}_2 < 1 \quad \textsf{and} \quad \norm{\left(\DQDV\right)^{-1} \DQDTH(\vv,\vp)}_2 <1,\\
\label{eq:thm2_column_offdiag}
\textsf{or,} \quad  \norm{\left(\DPDTH(\vv,\vq)\right)^{-1} \DQDTH(\vv,\vp)}_2 <1  \quad \textsf{and} \quad \norm{\left(\DQDV\right)^{-1} \DPDV}_2 < 1,
\end{align}

then the power flow Jacobian $\mJ(\vx) \in \R^{2|\setB| \times 2|\setB|}$ is non-singular, i.e., $\rank\left(\mJ(\vx)\right) = 2|\setB|$. Thus, the network is not in a state of voltage collapse due to Jacobian singularity, (e.g., saddle node bifurcation), and there exists unique voltage phase angle perturbations $\Delta \vtheta$ such that $[\Delta\vtheta^T, \Delta \vv^T]^T = \mJ(\vx)^{-1}[\Delta \vp^T,\Delta \vq^T]^T$ for any perturbations in (<ref>).

The proof of Theorem <ref> can be found in Appendix <ref>.

Theorem <ref> is a stronger condition than that of Theorem <ref>. If Theorem <ref> holds, then $\mJ(\vx)$ is guaranteed to be invertible. This is also a sufficient condition for the network to not be in a state of voltage collapse due to power flow Jacobian singularity, see <cit.> for more details. However, note that Theorem <ref> requires all the measurements in (<ref>) for each condition. In contrast, each condition of Theorem <ref> only requires the voltage magnitudes and one of either the reactive power or active power injections.

§ NUMERICAL RESULTS

In this section, we numerically validate the proposed algorithms developed in this paper. A Julia package based on  <cit.> has been developed for this work, and is publicly available.

We first present numerical analyses of the recovery guarantees in Section <ref>.
Then, we present numerical experiments for the phase retrieval program (<ref>). In Section <ref>, we  perform phase retrieval for a single time point at a nominal operating point specified by the network's AC power flow solution, with noise added. For this reason, the subscript $t$ is often dropped. In Section <ref>, actual time-series measurement data are used to apply the algorithm sequentially, including nonidealities such as asynchronous measurements, delays, and noise.

§.§ Validation of recovery guarantees

In addition to the proofs of Theorem <ref> and Theorem <ref> that are available in  Appendix <ref> and Appendix <ref>, respectively, we report numerical validation of these results for several Matpower <cit.> test networks. This includes the widely used  <cit.> test network, and the large-scale, realistic Texas A&M synthetic grids <cit.>. The quantities used for validating the theorem were computed with  <cit.>.

In Table <ref>, we show the percentage of the PQ buses that are in the largest subset of PQ buses $\setB$ that satisfy Theorem <ref>. For the buses that do not satisfy Theorem <ref>, we show the maximum length $r_{\sf worst}>0$ of the radii segment of the discs that crosses into the left-hand side of the complex plane.

In Table <ref>, we show the percentage of the PQ buses that are in the largest subset of PQ buses $\setB$ that satisfy Theorem <ref> (see the end of the proof of Theorem <ref> in Appendix <ref> for more details). For each test network, regardless of whether Theorem <ref> is satisfied for the entire set of PQ buses, we show the maximum spectral norm of the four matrix arguments in (<ref>), (<ref>).

Case     # PQ Buses     % Satisfying Thm. <ref>     $r_{\sf worst}$

9     100.0%     —

13     100.0%     —

24     95.83%     $1.4 \times 10^{-14}$

40     100.0%     —

64     100.0%     —

77     94.81%     8.32

162     96.91%     0.088

444     94.37%     3.014

1608     84.83%     28.31

Analysis of Theorem <ref> for PQ buses of various test cases

Case     # PQ Buses     % Satisfying Thm. <ref>     $\sigma_{\sf max}$

9     100.0%     0.876

13     100.0%     0.401

24     95.83%     1.437

40     100.0%     0.444

64     100.0%     0.473

77     100%     0.954

162     100%     0.698

444     99.77%     1.090

1608     99.69%     1.180

Analysis of Theorem <ref> for PQ buses of various test cases

§.§ Recovery of voltage phase angles and power-voltage phase angle sensitivities

§.§.§ Simulated measurements

To test the phase retrieval program (<ref>) on several standard Matpower <cit.> test cases, we simulate noisy measurements. We used  <cit.> to solve the AC power flow equations for each test case, yielding a nominal solution for the grid state $\vx^* = [\vtheta^T,\vv^T]^T$, with corresponding power flow Jacobian $\mJ(\vx^*)$. To provide a simulation of the measurements in (<ref>), we computed observed injection perturbations as

\begin{equation}
\label{eq:noisy-mismatches}
\Delta \vg^{\sf obs}(\vx^*) \triangleq \mJ(\vx^*) \vx^* + \begin{bmatrix} \vxi^p \\ \vxi^p \end{bmatrix},
\end{equation}

where the entries of the noise vectors $\vxi^p,\vxi^q  \in \R^n$ are distributed as $\xi^p_i,\xi_i^q \sim \setN(0,\sigma_{\sf meas}^2)$ for all $i \in \{1,\dots,n\}$. We then computed the state perturbations as $\Delta \vx^{\sf obs} \triangleq \mJ(\vx^*)^{-1}\Delta \vg^{\sf obs}$.

§.§.§ Recovery of voltage phase angles

We solved the phase retrieval program (<ref>) for various Matpower <cit.> test networks: the IEEE RTS 24-bus test case, the IEEE 30-bus test case, and the Grid Modernization Laboratory Consortium Reliability Test System () case <cit.>. In Figure <ref>, we show the estimated PQ bus voltage phase angles at the AC power flow solution of  for multiple noise levels, as in (<ref>). The ground truth bus voltage phase angles are shown in black. The results verify that the phase retrieval program (<ref>) is effective and robust to noise. Additional experiments for other test cases are available in Appendix <ref>.

Performance of (<ref>) by measurement noise level. The shaded ribbon indicates a $\pm 2 e$ error interval where $e$ is the absolute bus voltage phase error.

§.§.§ Recovery of the power-voltage phase angle sensitivity matrices

The structural symmetries of the power flow Jacobian (<ref>)  that were derived in Section <ref> allow us to recover the power-voltage phase angle sensitivity matrices, i.e., the two unknown blocks of the power flow Jacobian. This is achieved by simultaneously recovering the phase angle perturbations and the Jacobian submatrix decision variables $\DPDTH,\DQDTH$ in the phase retrieval program (<ref>), subject to the structural constraints (<ref>). In Fig. <ref>, we show the optimal values of the power-voltage phase angle sensitivity matrices obtained by solving the phase retrieval program (<ref>), compared with the deterministic matrices defined in (<ref>).

Recovery of the power-phase angle submatrices $\DPDTH,\DQDTH$ of the power flow Jacobian for the  network via the phase retrieval program (<ref>) with $\sigma_{\sf meas} = 0.1$. Bus indices correspond to PQ bus ordering of the basic network format of  <cit.>. Frobenius norm relative error is displayed on the second row.

§.§.§ Robustness to errors in the estimates of the power-voltage magnitude sensitivity matrices

We introduce errors to the power-voltage magnitude sensitivity matrices to simulate the application of the phase retrieval program (<ref>) using sensitivity matrices $\DPDV,\DQDV$ that were pre-estimated from past measurements of (<ref>). These errors are simulated by perturbing the deterministic, model-derived matrices with an additive Gaussian random matrix as

\begin{equation}
\label{eq:simulated-jacobian-estimation-error}
\DPDV^{\sf obs} = \DPDV(\vx^*) + \mXi^{\sf pv}, \quad
\textsf{and} \quad \DQDV^{\sf obs} = \DQDV(\vx^*) + \mXi^{\sf qv},
\end{equation}

where $\Xi^{\sf pv}_{i,k},\Xi^{\sf qv}_{i,k} \sim \setN(0,\sigma^2_{\sf jac})$ for all $i,k \in \{1,\dots,n\}$. Matrices $\DPDV^{\sf obs}$ and $\DQDV^{\sf obs}$ are used as inputs in the phase retrieval program (<ref>) and the power-phase angle sensitivity functions (<ref>) in Section <ref>.

§.§ Impact of measurement nonidealities

In practice, several nonidealities may be present in measurements of the form (<ref>), such as:

* The measurements may not be instantaneous. They may be averaged measurements over a particular sampling period, such as one hour.

* The measurements may not be time-synchronized. They may possess delays and asynchronous sampling.

Below, we assess the impact of these nonidealities on the proposed phase retrieval scheme. In Section <ref>, we subsample measurements provided by <cit.>, and also propose a model to analyze the impact of random delays and asynchronous measurements on the method. Simulations are then provided in Section <ref>.

Ground truth (blue) and estimated (orange dashed) voltage phase angles at 15 min. granularity, juxtaposed with ground truth 5 min. granularity voltage phase angles (black dots). Shown are PQ bus 2 “AGRICOLA" (left) and PQ bus 22 “BARKLA" (right) of the  network model.

§.§.§ Modeling asynchronous measurements and non-ideal sampling rates

In practice, measurements of the form (<ref>) regularly have delays and may be asynchronous. A simple model of this phenomenon can be obtained via a Bernoulli random vector $\vd \in \{0,1\}^n$, where the entries are i.i.d. as $d_i \sim \textsf{Bern}(p) $ for nodes $i=1,\dots,n$. In this model, the measurement data stream of node $i$ is delayed  with probability $p$, and is not delayed with probability $1-p$.

Results of a survey of 31 distribution utilities in 2021 found that $>$60% of respondents have access to AMI at sampling rates at least as frequent as 15 minutes <cit.>, while a smaller percentage had 5 minute sampling rates. We model this nonideality by subsampling 5-minute ground truth measurements to 15-minute granularity.

§.§.§ Simulations

In Fig. <ref>, the impacts of time delays, data rates, and measurement error are depicted. The ground-truth 5-minute time-series of active and reactive power load measurements are those provided by the  network model <cit.>. These ground truth measurements were used to sequentially solve the AC power flow equations, and obtain the voltage magnitudes and phase angles at 5-min intervals over the course of a year.  We then recover the voltage phase angles by applying the program (<ref>) sequentially using measurements of the form (<ref>). The measurements were corrupted in the following ways:

* The 5-minute ground truth $\vv_t,\vp_t,\vq_t$ measurements are subsampled to 15 min. granularity measurements.

* Node $i$ is delayed by 15 minutes relative to the true time-series with probability $p=0.1$, modeled by the Bernoulli random variable $\vd$. If $d_i=1$, then $\texttt{circshift}(\cdot)$ is applied to each measurement stream of bus $i$.

* Additive white Gaussian noise is added as in Section <ref>, with $\sigma_{\sf meas} = 0.1$.

The 5-min granularity “ground-truth" phase angles at the AC power flow solutions provided by  are shown as black dots. The title depicts the relative percentage error between the recovered phase angle time series (interpolated back to 5-min granularity) and the 5-min ground truth.

§.§ Baseline comparison with known topology

We compare the method with a standard state estimation approach where the topology is known. We achieve this comparison by using exact closed form expressions <cit.> for the complex-power to voltage magnitude and complex power-voltage phase angle sensitivities that depend on knowledge of the topology $\mY$. These matrices $\DSDV,\DSDTH \in \C^{n \times n}$ are defined as

\begin{align}
\label{eq:dsdv-def} \DSDV &\triangleq \diag(\vvbar)\left( \diag(\conj{\mY} \conj{\vvbar}) + \conj{\mY} \diag(\conj{\vvbar})\right) \diag(\vv)^{-1},\\
\label{eq:dsdth-def} \DSDTH &\triangleq j \diag(\vvbar)\left( \diag(\conj{\mY} \conj{\vvbar}) - \conj{\mY} \diag(\conj{\vvbar})\right),
\end{align}

where $\conj{(\cdot)}$ denotes the complex conjugate. Using (<ref>) and (<ref>), closed form expressions for the power flow Jacobian submatrices are given as
$\DPDV= \Re{\DSDV}$, $\DQDV =\Im{\DSDV}$,  $\DPDTH=\Re{\DSDTH}$, and $\DQDTH = \Im{\DSDTH}$.

The above matrices can then be used to build a deterministic (model-based) Jacobian $\mJ^{\sf model} \in \R^{2n \times 2n}$ to directly solve a standard state estimation problem of the form

\begin{equation}
\minimize_{\Delta \vtheta} \quad \norm{\mJ^{\sf model} [\Delta \vtheta^T, \Delta \vv^T]^T - \Delta \vg^{\sf obs} }_2^2,
\end{equation}

where $\DPDTH,\DQDTH$ are known and are not decision variables, and $\Delta \vv \triangleq \Delta \vv^{\sf obs}$. This stands in contrast with the program (<ref>), where $\DPDTH,\DQDTH$ are unknown and are included as decision variables.

In Fig. <ref>, we benchmark the performance of the phase retrieval program (<ref>) that does not require a model against this standard state estimation approach with a known model. The quantities $\sigma_{\sf meas}$ and $\sigma_{\sf jac}$ refer to the standard deviation of the additive Gaussian noise applied to $\Delta \vg^{\sf obs}$ and $\DPDV^{\sf obs},\DQDV^{\sf obs}$, respectively, as denoted in (<ref>) and (<ref>), respectively. Note that in Fig. <ref>, we are applying noise to both the measurements and $\DPDV,\DQDV$.

Impact of () model uncertainty on recovered phase angle relative error vs. measurement noise level. Shaded regions indicate $\pm$1 standard deviation of the relative errors computed over 20 bootstraps.

§ DISCUSSION

The proposed extension of the phase retrieval framework to electric power systems has several potential consequences, key limitations, and directions for future work. Below, we discuss these items in Section <ref>, Section <ref>, and Section <ref>, respectively.

§.§ Commentary on results

The proposed methods are a useful development from both a practical and theoretical perspective. These methods have important applications for numerous data-driven modeling scenarios in electric power systems where we only have access to measurements of the form (<ref>). The most obvious is in distribution systems, where it is typical to only have access to data like (<ref>). Applications in transmission systems are also appropriate, as large-scale PMU deployment remains an ongoing infrastructure upgrade challenge.

The results in Lemma <ref> and its application in the program (<ref>) are valuable results for data-driven power flow analysis. This is because they show that the state evolution matrices of the unmeasured voltage phase angle states can be recovered without knowledge of the voltage phase angles. The method can also be shown to be provably correct via the conditions based on spectral theory developed in Theorem <ref> and Theorem <ref>.

§.§ Limitations

The largest limitation of this research is that it primarily assumes a balanced power system that can be represented as a single-phase equivalent system.[Appendix <ref> develops a preliminary extension to multiphase unbalanced networks.] An additional limitation is the fact that the phase retrieval program (<ref>) relies on the structure of the power flow Jacobian, which is a time-varying matrix. While the matrix can be recursively estimated, e.g., as in <cit.>, or via Kalman filtering; this is a key computational limitation.  A potential solution to this problem is to reformulate it in terms of the bus admittance matrix. Furthermore, the algorithm is only suitable for the PQ buses of the system or a system that only has PQ buses.

§.§ Future work

There are several key directions for future work to improve the practical application of the proposed methods.

§.§.§ Multiphase extensions

One of the most clear research needs is generalizing the proposed algorithms to multiphase unbalanced networks. This can facilitate practical implementation beyond single-phase equivalent models in distribution systems and in other unbalanced systems. We provide preliminary results that indicate that the proposed algorithms can be generalized to multiphase unbalanced networks in Appendix <ref>.

§.§.§ Additional phase retrieval formulations

Furthermore, recent advances in the phase retrieval literature can be synthesized with the problem setting described in this paper. For example, <cit.> and <cit.> develop convex formulations for the phase retrieval problem. These formulations may provide numerical and theoretical advantages in comparison to the approach provided by PhaseCut <cit.>, which our work is largely based on. However, <cit.> requires an appropriate guess for the phase vector, and the conditions for the effectiveness of the algorithms are based on the quality and construction of this initialization vector, which may be difficult to obtain in this application setting. An important direction for future work is to determine how other aspects of the phase retrieval framework <cit.> can be applied to electric power systems and related to the AC power flow equations.

§.§.§ Partial/estimated topology knowledge

The limitation imposed by the time-varying nature of the power flow Jacobian could potentially be resolved by developing an extension of this work based on the nodal admittance matrix. This would require estimation of the network topology and parameters, e.g., as in <cit.>. The nodal admittance matrix is, in principle, not time-varying. This could possibly allow measurements of the form (<ref>) to be recovered without updating the power flow Jacobian.

§.§.§ Uncertain reactive power measurements

If reactive power measurements $\vq$ are unavailable, remotely readable reactive power control parameters—which are mandated to be available in the IEEE 1547-2018 standard <cit.>—can be used to infer reactive power measurements when only the $\vv,\vp$ variables of (<ref>) are available. Alternatively, distribution operators can use historical power factors and other heuristics, which are widely used in practice <cit.>, to infer the reactive power injections. There is an opportunity for future work to develop methods to exploit knowledge of known/estimated reactive power control parameters to obtain this knowledge.

§ CONCLUSION

This work proposed theory and algorithms to recover voltage phasors, current phasors, and power-voltage phase angle sensitivity matrices in electric power systems where the network model is unknown. The method requires only the voltage magnitudes, active power injections, and reactive power injections. The methods are directly inspired by the classical phase retrieval framework from the signal processing literature.

To achieve these outcomes, this research developed a novel representation of the structure of the power flow Jacobian matrix (Lemma <ref>) that depends on neither the topology of the network nor the voltage phase angles. The proposed representation  of the structure of the power flow Jacobian was then used in Section <ref> to provide a novel approach to recovering the voltage phase angles and a basis for an approximation of the power flow equations as a function of the phase angles. The method relies neither on partial observability of the phase angles nor on topological information. The proposed representation of the power flow Jacobian was also used to construct sufficient conditions (Theorem <ref>, Theorem <ref>) that certify when the method will work well. Consistent with the method itself, the conditions only require the measurements in (<ref>).

Foreseeably, the proposed theory and algorithms could have significant practical engineering implications. The results of this research can enable engineers to “fill in the information gap" in electric power networks that lack deployment of PMU measurements, as in (<ref>). It could also potentially save significant cost, by allowing fewer PMU measurements to be placed. As PMU deployment continues, these results can assist network control paradigms by retrieving comparable information to that provided by PMUs in areas that lack access to them.

§ ACKNOWLEDGEMENT

Samuel Talkington sincerely thanks Daniel Turizo for his feedback on the manuscript and valuable discussions.

§ PROOF OF THEOREM 1

For any matrix $\mA \in \C^{n \times n}$, by the Gershgorin Circle Theorem the eigenvalues of $\mA$ are guaranteed to lie in the union of the $i=1,\dots,n$ Gershgorin discs $\setG_i(\mA)$ of the matrix, i.e.,

\begin{equation}
\label{eq:apdx:thm1:gershgorin-thm}
\lambda_i(\mA) \in \bigcup_{i=1}^n \setG_i(\mA), \quad  i=1,\dots,n,
\end{equation}

\begin{equation}
\setG_i(\mA) \triangleq \{w \in \C : | w - A_{ii}| \leq \sum_{k: k\neq i} |A_{ik}|\} \subseteq \C.
\end{equation}

It suffices to show that $\DPDTH$ is invertible to show that (<ref>) holds. By (<ref>), the eigenvalues of $\DPDTH$ satisfy

\begin{equation}
\lambda_i\left(\DPDTH\right) \in \bigcup_{i \in \setB} \left\{w\in \C: \left|w - \frac{\partial p_i}{\partial \theta_i} \right| \leq \sum_{k \in \setB \setminus \{i\}} \left|\frac{\partial p_k}{\partial \theta_i}\right| \right\} \quad \forall i \in \setB.
\end{equation}

Therefore, if

\begin{equation}
\label{eq:thm1-step1}
\forall i \in \setB, \quad \left|\frac{\partial p_i}{\partial \theta_i} \right| > \sum_{k \in \setB \setminus \{i\}} \left| \frac{\partial p_k}{\partial \theta_i}\right|,
\end{equation}

then 0 fails to be an eigenvalue of $\DPDTH$, in which case, $\DPDTH$ is invertible and (<ref>) holds.

By Lemma <ref>, the symmetries in the entries of the power flow Jacobian (<ref>) hold, and we know that $\frac{\partial p_i}{\partial \theta_i} = v_i\frac{\partial q_i}{\partial v_i} - 2q_i$ for all $i \in \setB$, and $\frac{\partial p_i}{\partial \theta_k} = v_k \frac{\partial q_i}{\partial v_k}$ for all $k \in \setB\setminus\{i\}$.

Applying the power flow Jacobian structure (<ref>) to (<ref>), we obtain

\begin{equation}
\label{eq:thm1-step2}
\forall i \in \setB \quad |v_i| \sum_{k \in \setB \setminus \{i\}} \left| \frac{\partial q_k}{\partial v_i}\right| < \left|v_i\frac{\partial q_i}{\partial v_i} -2q_i \right|.
\end{equation}

Moreover, by applying the triangle inequality to (<ref>), we have

\begin{equation}
\label{eq:thm1-penultimate-row-offdiag}
\forall i \in \setB \quad |v_i| \sum_{k \in \setB \setminus \{i\}} \left| \frac{\partial q_k}{\partial v_i}\right|  < \left|v_i\frac{\partial q_i}{\partial v_i} \right| + 2 \left|q_i \right|.
\end{equation}

By construction $v_i>0$, so $|v_i| = v_i$. Rearranging (<ref>), if

\begin{equation}
\label{eq:thm1-apdx-row-condition}
\forall i \in \setB \quad |q_i| > \frac{1}{2} v_i \left( \sum_{k \in \setB \setminus \{i\}}
\left| \frac{\partial q_k}{\partial v_i}\right| - \left| \frac{\partial q_i}{\partial v_i}\right|
\right),
\end{equation}

then 0 fails to be an eigenvalue of $\DPDTH$. This means that $\DPDTH$ is invertible, and therefore, it is guaranteed that we can uniquely recover the phase angle deviations $\Delta \vtheta$ as (<ref>).

This completes the proof for the row off-diagonal condition in (<ref>). The column off-diagonal condition (<ref>) can be shown through the same application of the power flow Jacobian symmetry and the triangle inequality, and rearranging appropriately. We have that

\begin{align}
\forall i \in \setB \quad  \left| \frac{\partial p_i}{\partial \theta_i}\right| &= \left| v_i \frac{\partial q_i}{\partial v_i}  + 2 q_i\right|\\
&< v_i \left| \frac{\partial q_i}{\partial v_i} \right| + 2 |q_i|\\
&> \sum_{k \in \setB \setminus \{i\}}\left| v_k \frac{\partial q_i}{\partial v_k} \right| =  \sum_{k \in \setB \setminus \{i\}} v_k\left|  \frac{\partial q_i}{\partial v_k} \right|.
\end{align}

Rewriting (<ref>), we obtain the desired condition that if

\begin{equation}
\label{eq:thm1-apdx-column-condition}
\forall i \in \setB \quad |q_i| > \frac{1}{2}\left(\sum_{k \in \setB \setminus \{i\}} v_k \left| \frac{\partial q_i}{\partial v_k}\right| - v_i \left| \frac{\partial q_i}{\partial v_i} \right| \right),
\end{equation}

then 0 fails to be an eigenvalue of $\DPDTH$, and we can uniquely recover the phase angle perturbations $\Delta \vtheta$ as (<ref>).

Note that if either the row off-diagonal bound (<ref>) or the column off-diagonal bound (<ref>) holds for all $i \in \setB$, then 0 is guaranteed to not be an eigenvalue of $\DPDTH$, which means that $\DPDTH$ is guaranteed to be invertible. This means that the recovery of $\Delta \vtheta$ via (<ref>) is guaranteed to exist and be unique.

Identical logic can be applied to yield the bounds for the active power injection (<ref>). By Lemma <ref> and the triangle inequality, if for all $i \in \setB$

\begin{align}
\left|\frac{\partial q_i}{\partial \theta_i} \right| &= \left| -v_i \frac{\partial p_i}{\partial v_i} + 2 p_i\right|,\\
&\leq v_i \left| \frac{\partial p_i}{\partial v_i}  \right| + 2 |p_i|,\\
&> \sum_{k \in \setB \setminus\{i\}} \left| \frac{\partial q_k}{\partial \theta_i}\right|,\\
&= |v_i| \sum_{k \in \setB \setminus\{i\}} \left| \frac{\partial p_k}{\partial v_i}\right| = v_i \sum_{k \in \setB \setminus\{i\}} \left| \frac{\partial p_k}{\partial \theta_i}\right|,
\end{align}

or the corresponding column off-diagonal condition holds, then 0 is not an eigenvalue of $\DQDTH$, and there exists unique phase angle perturbations such that (<ref>) holds.

§ PROOF OF THEOREM 2

Case     Assumption <ref> holds?

All test networks satisfy Assumption <ref>.

The Gershgorin Circle Theorem can be generalized to block matrices, as shown in <cit.> and <cit.>. Consider a block matrix $\mM$ of the form

\begin{equation}
\mM = [\mA_{ij}] \quad \textsf{with blocks} \quad \mA_{ij} \in \C^{m \times m}, \quad i,j=1,\dots,n,
\end{equation}

where $\rank(\mA_{ii}) = m$ for all $i=1,\dots,n$. According to <cit.>, if $\lambda \in \C$ is an eigenvalue of $\mM$, then there exists  at least one  $i \in \{1,\dots,n\}$ such that

\begin{equation}
\label{eq:echeverria-block-gershgorin}
\sum_{j \neq i}^n \norm{(\mA_{ii} - \lambda \mId)^{-1} \mA_{ij}}_2 \geq 1,
\end{equation}

where $\norm{\mA}_2$ is the matrix 2-norm, or spectral norm, which is defined as $\norm{\mA}_2 \triangleq \max_{\vx \neq \boldsymbol{0}} \frac{\norm{\mA \vx}}{\norm{\vx}_2} = \sigma_{\sf max}(\mA)$, where $\sigma_{\sf max}(\mA)$ is the largest singular value of $\mA$.

Furthermore, every eigenvalue $\lambda$ of $\mM$ lies in the unions of these discs.

\begin{equation}
\lambda(\mM) \in \bigcup_{i =1}^n \left\{ w \in \C : \sum_{j\neq i}^n \norm{(\mA_{ii} - w \mId)^{-1} \mA_{ij}}_2 \geq 1\right\}.
\end{equation}

If all the blocks $\mA_{ij}$ have dimension $1 \times 1$, then this statement reduces to the classical Gershgorin Circle Theorem. Suppose that $\mM$ takes the form of a $2n \times 2n$ block matrix that satisfies Assumption <ref>, which we write as

\begin{equation}
\renewcommand*{\arraystretch}{1.5}
\mM \triangleq \begin{bmatrix}
\mA_{11} & \mA_{12}\\
\mA_{21} & \mA_{22}
\end{bmatrix},
\end{equation}

where $\rank(\mA_{11}) = \rank(\mA_{22}) = n$.
Then, applying (<ref>), if 0 is an eigenvalue of $\mM$, then

\begin{equation}
\norm{\mA_{11}^{-1} \mA_{12}}_2 \geq 1  \quad \textsf{or}  \quad \norm{\mA_{22}^{-1} \mA_{21}}_2 \geq 1.
\end{equation}

The contraposition is that if

\begin{equation}
\label{eq:apdx:thm2:blockgershgorinsufficient}
\norm{\mA_{11}^{-1} \mA_{12}}_2 < 1  \quad \textsf{and}  \quad \norm{\mA_{22}^{-1} \mA_{21}}_2 <1,
\end{equation}

then 0 is guaranteed to not be an eigenvalue of $\mM$. That is, (<ref>) is sufficient to certify the invertibility of the block matrix $\mM$.

We can complete the proof by applying the above logic to the power flow Jacobian. Let $\mJ \in \R^{2|\setB|\times 2|\setB|}$ be the power flow Jacobian for PQ buses $\setB$ under study. Recall that we defined $\vx = [\vtheta^T, \vv^T]^T \in \R^{2|\setB|}$ as a guess for the network state. We want to establish when the entire block matrix

\begin{equation}%\arraycolsep=1.4pt\def\arraystretch{1.5}
\renewcommand*{\arraystretch}{1.5}
\mJ(\vx) \triangleq \begin{bmatrix}
\DPDTH(\vx) & \DPDV(\vx)\\
\DQDTH(\vx) & \DQDV(\vx)
\end{bmatrix}
\end{equation}

is invertible—with $\vtheta$ unknown. This is equivalent to showing that $0$ is not an eigenvalue of $\mJ$. By Assumption <ref>, $\DPDTH$ and $\DQDV$ are invertible[Assumption <ref> empirically holds for all test cases studied, including large synthetic real-world networks, as shown in Table <ref>.]. Then by <cit.>, if $\lambda \in \C$ is an eigenvalue of $\mJ$, then either

\begin{equation}
\norm{\left(\DPDTH - \lambda \mId\right)^{-1} \DPDV}_2 \geq 1
\end{equation}

\begin{equation}
\norm{\left(\DQDV - \lambda \mId\right)^{-1} \DQDTH}_2 \geq 1.
\end{equation}

Let $\mJ_{\#}: \R^{|\setB|} \times \R^{|\setB|} \times \R^{|\setB|} \mapsto \R^{2 |\setB| \times 2 |\setB|}$ be a matrix-valued function of the measurements in (<ref>). Consider the power-voltage phase angle sensitivity matrix functions $\DPDTH(\cdot,\cdot),\DQDTH(\cdot,\cdot): \R^{|\setB|} \times \R^{|\setB|} \mapsto \R^{|\setB| \times |\setB|}$ as defined in (<ref>). By parameterizing $\mJ_{\#}$ by $\vtheta$, we can write

\begin{equation}
\renewcommand*{\arraystretch}{1.5}
\mJ_{\#}(\vv,\vp,\vq | \vtheta)  \triangleq \begin{bmatrix}
\DPDTH(\vv,\vq) & \DPDV(\vv|\vtheta)\\
\DQDTH (\vv,\vp) & \DQDV(\vv|\vtheta)
\end{bmatrix},
\end{equation}

which is a function of (<ref>). Then, by Lemma <ref>, we have $\lambda \in \lambda(\mJ) \iff \lambda \in \lambda( \mJ_{\#})$.

Therefore, it suffices to show that 0 is not an eigenvalue of $ \mJ_{\#}$ to complete the proof.  Applying <cit.> again, if $\lambda$ is an eigenvalue of $\mJ_{\#}$, then either

\begin{equation}
\norm{\left(\DPDTH(\vv,\vq) - \lambda \mId\right)^{-1} \DPDV}_2 \geq 1
\end{equation}

\begin{equation}
\norm{\left(\DQDV - \lambda \mId\right)^{-1} \DQDTH(\vv,\vp)}_2 \geq 1.
\end{equation}

Finally, applying Lemma <ref> and forming the contrapositive, if

\begin{equation}
\label{eq:apdx2:proofthm2:final}
\norm{\left(\DPDTH(\vv,\vq)\right)^{-1} \DPDV}_2 < 1, \quad \textsf{and} \quad  \norm{\left(\DQDV\right)^{-1} \DQDTH(\vv,\vq)}_2 <1,
\end{equation}

then 0 fails to be an eigenvalue of $\mJ$. Note that the inequalities (<ref>) depend solely on (<ref>). This completes the proof for (<ref>). Analogous arguments can be applied to $\mJ_{\#}^T$ to yield the second set of inequalities (<ref>), because $\sigma(\mJ) = \sigma(\mJ^T)$ and $\lambda(\mJ) = \lambda(\mJ^T)$.

As in Theorem <ref>, conditions (<ref>) and (<ref>) allow us to certify Theorem <ref> for a subset of PQ buses $\setB$ under study. This can be achieved by exploiting the properties of singular values of principle submatrices. Let $\mS_{\bullet} \in \R^{(|\setB| -b) \times (|\setB| -b)}$ be a principal submatrix of one of the Jacobian submatrix products,

\begin{equation}
\label{eq:jac-submatrix-products}
\mS \in \left\{ \DPDTH(\vv,\vq)^{-1} \DPDV, \ \ \DQDV^{-1} \DQDTH(\vv,\vp), \ \ \DPDTH(\vv,\vq)^{-1}\DQDTH(\vv,\vp), \ \ \DQDV^{-1} \DPDV \right\},
\end{equation}

with the corresponding rows and columns of $b < |\setB|$ buses removed. Note that $\mS$ is one of the arguments of $\norm{\cdot}_2$ for the spectral norm  conditions in (<ref>) or (<ref>). Then, by the submatrix property of singular values <cit.>, $\forall i=1,\dots,|\setB| -b,$ the singular values of $\mS_{\bullet} $ satisfy $\sigma_{i + 2b}(\mS) \leq \sigma_i(\mS_{\bullet}) \leq \sigma_i(\mS)$ for all Jacobian submatrix products $\mS$ as defined in (<ref>). Thus, the conditions (<ref>) and (<ref>) remain sufficient for the $|\setB|-b$ buses.

§ ADDITIONAL NUMERICAL RESULTS

§.§ Extension to three-phase unbalanced networks

Consider a three-phase unbalanced distribution network with $n$ buses. The compound admittance matrix of the network, $\mY \in \C^{3n \times 3n}$, is a block matrix $\mY= \left[\mM_{i,j}\right]_{i,j=1}^n$. Block $\mM_{i,j} \in \C^{3 \times 3}$ describes the mutual admittance between the three phases of branch $i \to j$ with entries $\mM_{i,j} = -\big[y^{\phi}_{i,j}\big]_{\phi \in \setF_{ij}}$. Set $\setF_{ij}$ denotes the phases of the branch $i \to j$.  The diagonal blocks $\mM_{i,i}$ are computed as $\mM_{i,i} = \mY_i -\sum_{j:j\neq i} \mM_{i,j}$, where $\mY_i \in \C^{3 \times 3}$ is the equivalent admittance of the loads connected at bus $i$, and the admittance-to-ground. We refer the reader to <cit.> for more comprehensive information on modeling three-phase networks with the bus admittance matrix.

As a proof-of-concept to verify that the proposed symmetry in Lemma <ref> and the phase retrieval method (<ref>) is applicable in the unbalanced case, we compute the power-voltage phase angle sensitivities using the compound admittance matrix $\mY$. According to <cit.>, we can express the complex-power to voltage magnitude and complex power-voltage phase angle sensitivities as complex matrices $\DSDV,\DSDTH \in \C^{n \times n}$ of the form

\begin{align}
\DSDV &= \diag(\vvbar)\left( \diag(\conj{\mY} \conj{\vvbar}) + \conj{\mY} \diag(\conj{\vvbar})\right) \diag(\vv)^{-1},\\
\DSDTH &= j \diag(\vvbar)\left( \diag(\conj{\mY} \conj{\vvbar}) - \conj{\mY} \diag(\conj{\vvbar})\right),
\end{align}

where $\vvbar \in \C^n$ are the complex voltage phasors, $\vv = |\vvbar|$ are the voltage magnitudes, and $\conj{(\cdot)}$ denotes the complex conjugate. We then have $\DPDV,\DQDV = \Re{\DSDV},\Im{\DSDV}$, and $\DPDTH,\DQDTH=\Re{\DSDTH},\Im{\DSDTH}$.

These closed-form expressions allow us to verify that the structural symmetries derived in Lemma <ref> hold in the unbalanced case. In Table <ref>, we show the relative error between the deterministic $\DPDTH,\DQDTH$ sensitivities and the model-free power-voltage phase angle sensitivity functions $\DPDTH(\vv,\vq),\DQDTH(\vv,\vp)$ as defined in  (<ref>). The results are shown for a three-phase unbalanced distribution test network with three buses.

Data used to compute the values in Table <ref> were provided by  <cit.>. The case used was the three-phase unbalanced . The  polar AC power flow model was used to obtain a solution. The voltage magnitudes, active powers, and reactive powers were then computed at the AC power flow solution, and used to generate $\DPDTH(\vv,\vq)$ and $\DQDTH(\vv,\vp)$ as defined in  (<ref>).

Quantity     Value

$\norm{\DPDTH -\DPDTH(\vv,\vq)}_F \Big/ \norm{\DPDTH}_F$      $2.510 \times 10^{-8}$

$\norm{\DQDTH - \DQDTH(\vv,\vp)}_F \Big/ \norm{\DQDTH}_F$     $1.725 \times 10^{-7}$

Verification that the structure expressions (<ref>) of Lemma <ref> hold for multiphase unbalanced networks.

§.§ Verifying Lemma <ref> in deterministic Newton-Raphson power flow

2*Iteration      4c| Relative error $\frac{\|(\cdot) - \hat{(\cdot)}\|}{\|(\cdot)\|} \times 10^{16}$
$\vtheta$      $\vell $     $\DPDTH$     $\DQDTH$

t=1     $4.922 $     $1.513 $     $2.407 $     5.096

t=2     $5.324 $     1.907      2.525     3.381

t=3     $7.444 $     1.979     3.020     5.118

t=4     $4.377 $     1.514      2.066     3.001

t=5     $4.681 $     1.497     2.904      6.780

t=6     $3.530 $     1.182     2.973     5.583

Relative error ($\times 10^{16}$) between the phaseless and deterministic Newton-Raphson power flow model with no noise for the  test case <cit.>.

The phaseless sensitivity functions $\DPDTH(\vv,\vq),\DQDTH(\vv,\vp)$ as defined in (<ref>) are equivalent to the deterministic matrices $\DPDTH,\DQDTH$ in (<ref>), (<ref>). This can be shown by using $\DPDTH(\vv,\vq),\DQDTH(\vv,\vp)$ to iteratively solve the Newton-Raphson power flow model (<ref>) in the noiseless case without considering $\vtheta$ as a grid state. To this end, we compare the quantities $\DPDTH,\DQDTH,$ $\vtheta,$ and $\vell$ at every Newton-Raphson iteration $t=1,\dots,$ when solving the  test case with both the phaseless and deterministic $\DPDTH,\DQDTH$ in Table <ref>. The results  indicate that—with no observations of the phase angles—comparative results for the Newton-Raphson power flow model are achieved within machine precision at each iteration. This demonstrates that the results in Lemma <ref> are physically valid.

Performance of (<ref>) by measurement noise level. The shaded ribbon indicates a $\pm 2 e$ error interval where $e$ is the absolute bus voltage phase error.

§.§ Recovery of voltage phase angles and power-voltage phase angle sensitivities

Additional numerical simulations are shown in this appendix for the same phase retrieval technique (<ref>) used for the simulations in Section <ref>. In particular, we report results for the IEEE 30-bus test system and the IEEE 24-bus reliability test system (RTS).

In Fig. <ref>, we show the estimation of the bus voltage phase angles at the AC power flow solution using (<ref>) for additional test cases with multiple noise levels. The experiments are the same as was done in Section  <ref>, which yielded the prior results in Fig. <ref> for the  network model.   Fig. <ref> shows the power-voltage phase angle sensitivity recovery results for additional test cases, in the same manner as was done in Section  <ref>, which yielded the prior results in Fig. <ref> for the  network model.

Recovery of the power-voltage phase angle submatrices of the AC power flow Jacobian for additional test cases via the
phase retrieval program (<ref>) with $\sigma_{\sf meas} = 0.1$. Bus indices correspond to PQ bus ordering of the basic network format of  <cit.>. Frobenius norm relative error is displayed on the second row.