Dynamic Distribution State Estimation Using Synchrophasor Data

TL;DR

This paper presents a fast, robust, and real-time algorithm for distribution state estimation using streaming synchrophasor data, leveraging time-varying optimization and prediction-correction methods.

Contribution

It introduces a novel online algorithm for dynamic distribution state estimation that is computationally efficient and capable of tracking states under measurement outliers.

Findings

Algorithm provably tracks state variables accurately.

Computational efficiency achieved without matrix-inverse calculations.

Robustness to measurement outliers demonstrated.

Abstract

The increasing deployment of distribution-level phasor measurement units (PMUs) calls for dynamic distribution state estimation (DDSE) approaches that tap into high-rate measurements to maintain a comprehensive view of the distribution-system state in real time. Accordingly, this paper explores the development of a fast algorithmic framework by casting the DDSE task within the time-varying optimization realm. The time-varying formulation involves a time-varying robustified least-squares approach, and it naturally models optimal trajectories for the estimated states under streaming of measurements. The formulation is based on a linear surrogate of the AC power-flow equations, and it includes an element of robustness with respect to measurement outliers. The paper then leverages a first-order prediction-correction method to achieve simple online updates that can provably track the state variables from heterogeneous measurements. This online algorithm is computationally efficient as it relies on the Hessian of the cost function without computing matrix-inverse. Convergence and bounds on the estimation errors of proposed algorithm can be analytically established.

Equations

$$\mathsf{diag}({\mathbf{H}}^{\text{T}}({\mathbf{i}}^{\Delta }{)}^{\ast })\mathbf{v}+{\mathbf{s}}^Y=\mathsf{diag}(\mathbf{v}){\mathbf{i}}^{\ast },$$

$${\mathbf{s}}^{\Delta }=\mathsf{diag}\left(\mathbf{H}\mathbf{v}\right)({\mathbf{i}}^{\Delta }{)}^{\ast },\mathbf{i}={\mathbf{Y}}_{L0}{\mathbf{v}}_0+{\mathbf{Y}}_{LL}\mathbf{v},$$

$$\mathbf{z}:=[\Re \{\mathbf{v}{\}}^{\text{T}},\Im \{\mathbf{v}{\}}^{\text{T}}{]}^{\text{T}}.$$

$$\tilde{\mathbf{z}}$$

$$\left[\begin{array}{c}{\mathbf{y}}_v^{(k)}\\ {{\mathbf{y}}_u^Y}^{(k)}\\ {{\mathbf{y}}_u^{\Delta }}^{(k)}\end{array}\right]=$$

$$+\left[\begin{array}{c}{\mathbf{m}}_v\\ 0\\ 0\end{array}\right]+\left[\begin{array}{c}{\mathbf{n}}_v^{(k)}\\ {{\mathbf{n}}_u^Y}^{(k)}\\ {{\mathbf{n}}_u^{\Delta }}^{(k)}\end{array}\right]$$

$${\mathbf{y}}^{(k)}={\mathbf{G}}^{(k)}{\mathbf{u}}^{(k)}+\bar{\mathbf{m}}+{\mathbf{n}}^{(k)}$$

$${\mathbf{G}}^{(k)}:=\left[\begin{array}{cc}{\mathbf{J}}_v{{\mathbf{M}}^{\mathbf{Y}}}^{(k)}& {\mathbf{J}}_v{{\mathbf{M}}^{\mathbf{\Delta }}}^{(k)}\\ {\mathbf{J}}^Y& 0\\ 0& {\mathbf{J}}^{\Delta }\end{array}\right].$$

$${\ell }^{(k)}(\mathbf{u}):=$$

$$+\mathsf{H}\left({{\mathbf{y}}_u^{\Delta }}^{(k)}-{\mathbf{J}}^{\Delta }{\mathbf{u}}^{\Delta }\right)$$

$$\displaystyle\mathsf{H}_{i}(\epsilon_{i}):=\left\{\begin{array}[]{ll}-\delta\epsilon_{i}-\delta^{2}/2,&\textrm{if~{}}\epsilon_{i}<-\delta\\ |\epsilon_{i}|^{2}/2,&\textrm{if~{}}|\epsilon_{i}|\leq\delta\\ \delta\epsilon_{i}-\delta^{2}/2,&\textrm{if~{}}\epsilon_{i}>\delta\end{array}\right.$$

$${\mathsf{P}}^{(k)}(\mathbf{u}):\underset{\mathbf{u}\in {\mathbb{R}}^{N_{\varphi }}}{\min }{f}^{(k)}(\mathbf{u}):={\ell }^{(k)}(\mathbf{u})+r^{(k)}(\mathbf{u}),k\in \mathbb{N}$$

$${\nabla }_x{f}^{(k+1)}({\mathbf{x}}^{(k+1|k)})=(1-\gamma ){\nabla }_x{f}^{(k)}({\mathbf{x}}^{(k)}),$$

$$(1-\gamma ){\nabla }_x{f}^{(k)}({\mathbf{x}}^{(k)})={\nabla }_x{f}^{(k+1)}({\mathbf{x}}^{(k+1|k)})\approx {\nabla }_x{f}^{(k)}({\mathbf{x}}^{(k)})$$

$$+{\nabla }_{xx}{f}^{(k)}({\mathbf{x}}^{(k)})({\mathbf{x}}^{(k+1|k)}-{\mathbf{x}}^{(k)})+h{\nabla }_{tx}{f}^{(k)}({\mathbf{x}}^{(k)}).$$

$$\begin{gathered} {\mathbf{x}}^{(k+1|k)}={\mathbf{x}}^{(k)}-{\nabla }_{xx}{f}^{(k)}({\mathbf{x}}^{(k)}{)}^{-1}\times \\ \left(\gamma {\nabla }_x{f}^{(k)}({\mathbf{x}}^{(k)})+h{\nabla }_{tx}{f}^{(k)}({\mathbf{x}}^{(k)})\right), \end{gathered}$$

$${\mathbf{x}}^{(k+1|k)}={\mathrm{argmin}}_{\mathbf{x}}$$

$${\hat{f}}^{(k)}(\mathbf{x}):=$$

$$-$$

$${\hat{\mathbf{x}}}_{p+1}=$$

$$\displaystyle+h\nabla_{tx}f^{(k)}(\bm{x}_{k})\big{]},~{}~{}~{}p=0,\ldots,P-1$$

$${\hat{\mathbf{x}}}_{c+1}={\hat{\mathbf{x}}}_c-\beta {\nabla }_x{f}^{(k+1)}({\hat{\mathbf{x}}}_c)({\hat{\mathbf{x}}}_c-{\mathbf{x}}^{(k)})$$

$${\bar{\mathbf{u}}}_{p+1}$$

$$+h{\nabla }_{t\mathbf{u}}{f}^{(k-1)}({\hat{\mathbf{u}}}^{(k-1)})+\gamma \left({\nabla }_{\mathbf{u}}{f}^{(k-1)}({\hat{\mathbf{u}}}^{(k-1)})\right)]$$

$${\hat{\mathbf{z}}}^{(k|k-1)}=\left[\begin{array}{cc}{{\mathbf{M}}^{\mathbf{Y}}}^{(k-1)}& {{\mathbf{M}}^{\mathbf{\Delta }}}^{(k-1)}\end{array}\right]{\hat{\mathbf{u}}}^{(k|k-1)}+\mathbf{m}.$$

$${\bar{\mathbf{u}}}_{c+1}$$

$${\hat{\mathbf{z}}}^{(k)}=\left[\begin{array}{cc}{{\mathbf{M}}^{\mathbf{Y}}}^{(k)}& {{\mathbf{M}}^{\mathbf{\Delta }}}^{(k)}\end{array}\right]{\hat{\mathbf{u}}}^{(k)}+\mathbf{m}.$$

$$\Vert {\nabla }_{t\mathbf{u}}{f}^{(k)}(\mathbf{u})\Vert \le C_0,\text{for all }k\in \mathbb{N}.$$

$$\nu \le \Vert {\nabla }_{\mathbf{u}\mathbf{u}}{f}^{(k)}(\mathbf{u})\Vert \le L,\text{ for all }k\in \mathbb{N}.$$

$$\begin{gathered} \Vert {\nabla }_{\mathbf{u}\mathbf{u}}{f}^{(k)}(\mathbf{u})\Vert =\Vert {\nabla }_{\mathbf{u}\mathbf{u}}{\ell }^{(k)}(\mathbf{u})+{\nabla }_{\mathbf{u}\mathbf{u}}{r}^{(k)}(\mathbf{u})\Vert \le \\ \Vert {{\mathbf{G}}_v^{(k)}}^{\text{T}}{\mathbf{G}}_v^{(k)}\Vert +\Vert {{\mathbf{J}}^Y}^{\text{T}}{\mathbf{J}}^Y\Vert +\Vert {{\mathbf{J}}^{\Delta }}^{\text{T}}{\mathbf{J}}^{\Delta }\Vert +L_r. \end{gathered}$$

$$\alpha <2/L,\beta <2/L,$$

$${\varrho }_{\mathrm{P}}=\max \{|1-\alpha \nu |,|1-\alpha L|\},{\varrho }_{\mathrm{C}}=\max \{|1-\beta \nu |,|1-\beta L|\}.$$

$$\tau_{0}:=\varrho_{\mathrm{C}}^{C}\left[\varrho_{\mathrm{P}}^{P}+(\varrho_{\mathrm{P}}^{P}+1)\Big{(}1-\gamma+\gamma\frac{2L}{\nu}\Big{)}\right]<1.$$

$$\underset{k\to \infty }{\mathrm{limsup}}\Vert {\mathbf{u}}^{(k)}-{\mathbf{u}}^{(\ast ,k)}\Vert =O({\varrho }_{\mathrm{C}}^{C}h).$$

Full text

I Introduction

Recently, power distribution networks have witnessed an increasing connection of renewable energy sources, electric vehicles, energy storage systems, among other distributed energy resources. These transformations have propelled the development and deployment of advanced sensing, communications, and control technologies. In particular, the synchrophasor technology based on distribution-level phasor measurement units (PMUs) [1] has equipped distribution systems operators with synchronized, low latency and high-resolution measurements that can be collected on a fast time scale. It is timely to design efficient and effective distribution situational awareness modules that can leverage the availability of PMU data therein.

Distribution dynamic state estimation (DDSE) is a fundamental tool that can enable general distribution operations and control tasks. This is similar to the operational paradigm for wide-area transmission grids [2]. Nonetheless, distribution networks have unique characteristics compared to transmission systems, such as unbalanced loads and higher resistance-to-reactance ratios [3]. Specifically, the unbalanced nature of distribution systems necessitates the multi-phase modeling at increased dimension and coupling. Moreover, nodes with zero injections commonly exist while historic load data are often included to improve redundancy. The high variation of accuracy from both types of data can lead to numerical conditioning issues, along with the higher resistance-to-reactance ratios. Thus, the dynamic SE approaches traditionally developed for transmission grids may not be directly applicable for distribution systems with PMU data.

Traditionally, distribution SE research has been limited to the static setting, constrained by meter availability and data rates. Several efforts have been focused on addressing the numerical conditioning issue using e.g., branch current formulation [4], conversion to current measurements [5, 6], zero-injection information [7, 8]. More recent work has considered the incorporation of synchrophasor data for static SE, but using the purely linear SE assuming sufficient observability from PMU data only [9, 10, 11]. This is not yet a reality for distribution networks, and thus legacy meters should still be included for distribution SE. Motivated by the availability of fast PMU measurements, dynamic SE methods has been developed relying on the recursive Kalman filter updates for distribution systems [12, 13]. Note that dynamic SE in transmission systems also uses variations of Kalman filtering recursions as the “workhorse” algorithm; see e.g., [14, 15]. Please see a recent review on distribution SE in [3].

This paper aims to develop an efficient and effective DDSE solution technique that can address the unique characteristics of distribution networks and fast sampling rates of PMU data. We leverage the linearized multi-phase AC power flow model recently developed in [16]. The model can account for unbalanced operation, as well as for wye and delta connections. Based on this model, we formulate a new DDSE problem by advocating a time-varying optimization formalism. The time-varying formulation involves a time-varying robustified least-squares approach, and it naturally models optimal trajectories for the estimated states under streaming of measurements. Based on the proposed time-varying optimization model for the DDSE task, this paper proposes an algorithmic framework to track the state of the distribution system by leveraging running prediction-correction methodologies [17, 18]. Prediction-correction methods involves two phases, sequentially implemented at each time step: i) a prediction phase where, based on the measurements collected up the current time instant, the algorithm attempts to predict the optimal solution of the next time period by exploring into intrinsic temporal correlations of the cost function; and, ii) once a new datum/measurement becomes available, the correction phase refines the predicted solution. To facilitate the development of computationally affordable algorithms, the paper considers first-order prediction-correction (FOPC) methods that rely on the Hessian of the cost function, instead of requiring the computation of its inverse [18]. FOPC methods are very attractive for DDSE problems where measurements are collected at high frequency by distribution-level PMUs or even distributed energy resources to enhance real-time situational awareness. In particular, FOPC have merits in the following two cases:

$\bullet$ The DDSE task has a given computational budget to estimate the state before a new measurement is collected and processed. In this case, in par with the computational demand, FOPC outperforms traditional iterative algorithms.

$\bullet$ FOPC can perform a prediction of the state while waiting for the measurement to be transmitted from the PMUs to the state estimator; once the measurement is received, the correction step can be performed. The prediction stage is shown to enable substantial improvements in terms of tracking performance.

Relative traditional approaches based on Kalman filtering, it is worth pointing out that: (i) FOPC provides an appreciable flexibility to include a variety of performance metrics in the cost function of the problem; for example, the paper will demonstrate how the DDSE can be easily robustified by modifying the cost function. (ii) KF requires covariance matrices; in a DDSE setting where heterogeneous measurements are collected at different rates, it is practically challenging to obtain accurate estimates of noise covariance matrices. And, (iii) KF typically needs an Hessian inverse computations, which can be too computationally burdensome; on the other hand, FOPC relies on first-order updates. Overall, FOPC is naturally data-driven, while KF is grounded on models for the dynamics and the noises. Hence, using the preferred time-varying optimization algorithms, our proposed FOPC-based DDSE methods are more computationally efficient and flexible to incorporate various types of measurements.

II Modeling and Problem Statement

We consider a generic multi-phase unbalanced distribution system with multiphase nodes collected in the set ${\cal N}\cup\{0\}$, ${\cal N}:=\{1,\ldots,N\}$, and distribution line segments collected in the set of edges ${\cal E}:=\{(m,n)\}$. Node [math] denotes the three-phase slack bus, i.e., the point of connection of the distribution grid with the rest of the electrical system. At each multiphase node, the loads can be either wye- or delta-connected [19], with the number of each type being $N^{Y}$, $N^{\Delta}$ respectively.

We briefly introduce the AC power-flow model for multiphase distribution systems (a comprehensive description can be found in [19, 16]). To this end, let ${\bf v}$ be a vector collecting the line-to-ground voltages in all phases of the nodes in ${\cal N}$; similarly, vector ${\bf i}$ collects all the phase net current injections, ${\bf i}^{\Delta}$ the phase-to-phase currents in all the delta connections, and vectors ${\bf s}^{Y}$ and ${\bf s}^{\Delta}$ collect the net complex powers injected for wye- and delta-connected loads, respectively. Let $N_{\phi}$ denote the total number of single-phase connections; for example, if all the nodes are three-phase, it follows that $N_{\phi}=3N$. With these definitions in place, the AC power-flow equations can be compactly written as:

[TABLE]

where $[{\bf Y}_{L0},{\bf Y}_{LL}]\in\mathbb{C}^{N_{\phi}\times(3+N_{\phi})}$ is the submatrix of the admittance matrix ${\bf Y}$ by eliminating the slack-bus rows, while ${\bf H}$ is a $N_{\phi}\times N_{\phi}$ block-diagonal matrix mapping ${\bf i}^{\Delta}$ to line currents; see [16] for a detailed description. For real-valued notations, define the $2N_{\phi}\times 1$ rectangular-form voltage vector

[TABLE]

Similarly, define the vectors ${\bf u}^{Y}:=[\Re\{{\bf s}^{Y}\}^{\textsf{T}},\Im\{{\bf s}^{Y}\}^{\textsf{T}}]^{\textsf{T}}$ and ${\bf u}^{\Delta}:=[\Re\{{\bf s}^{\Delta}\}^{\textsf{T}},\Im\{{\bf s}^{\Delta}\}^{\textsf{T}}]^{\textsf{T}}$ for power variables.

The proposed approach leverages a linearized AC power flow model to facilitate the development of computationally affordable algorithms that can be implemented in real time. For any given complex powers ${\bf s}^{Y},{\bf s}^{\Delta}$ and its voltage solution ${\bf v}$, this paper leverages a fixed-point approximation of (1) as detailed in [16], to obtain a linearized model that exactly touches upon $({\bf s}^{Y},{\bf s}^{\Delta},{\bf v})$ and $(\mathbf{0},\mathbf{0},{\bf w})$, where ${\bf w}$ is the zero-load voltage solution. Notice that, relative to alternative linearization techniques (e.g., [20, 21, 22, 23]), the approach in [16] accounts for both wye and delta connections.

Accordingly, consider the following linearized model of (1):

[TABLE]

where the model parameters ${\bf M}^{Y}\in\mathbb{R}^{2N_{\phi}\times 2N^{Y}_{\phi}}$, ${\bf M}^{\Delta}\in\mathbb{R}^{2N_{\phi}\times 2N^{\Delta}_{\phi}}$, and $\mathbf{m}\in\mathbb{R}^{2N_{\phi}\times 1}$, with ${\bf M}:=[{\bf M}^{Y},{\bf M}^{\Delta}]$ and ${\bf u}:=[({\bf u}^{Y})^{\textsf{T}},({\bf u}^{\Delta})^{\textsf{T}}]^{\textsf{T}}$. The vector $\mathbf{m}:=[\Re\{\mathbf{w}\}^{\textsf{T}},\Im\{\mathbf{w}\}^{\textsf{T}}]^{\textsf{T}}$ ensures that $(\mathbf{0},\mathbf{0},{\bf w})$ always satisfies the model (3) for ${\bf u}=\mathbf{0}$. In addition, matrices ${\bf M}^{\Delta}$ and ${\bf M}^{Y}$ are computed based on the network parameters and a given voltage profile. It is worth re-iterating that a linear model is leveraged in this paper to obtain the time-varying convex optimization model described in the next section and, based on that, synthesize dynamic algorithmic solutions with provable tracking properties.

II-A Problem statement

We introduce the measurement model for our $\mu$PMU assisted DDSE problem. Assume that the temporal axis is discretized as $t_{k}=hk$, where $k=0,1,\ldots,$ and the sampling period is $h=t_{k+1}-t_{k}$. Assume that $\mu$PMUs are located at a subset of nodes ${\cal M}_{\bf v}\subset{\cal N}$; let $M_{\bf v}$ denote the number of line-to-ground voltage measurements collected from the multi-phase nodes in ${\cal M}_{\bf v}$. These $\mu$PMUs can obtain accurate measurements of the voltages in rectangular coordinates, and they can produce measurements in real time. Furthermore, net injected powers from wye and delta connections are measured at the multi-phase nodes in ${\cal M}_{\bf u}^{Y}$ and ${\cal M}_{\bf u}^{\Delta}$, respectively, with ${\cal M}_{\bf u}:={\cal M}_{\bf u}^{Y}\bigcup{\cal M}_{\bf u}^{\Delta}\subset{\cal N}$. Accordingly, let $M_{\bf u}$ denote the number of nodes in ${\cal M}_{\bf u}$.

In this setting, at each time instant $t_{k}$ a set of new measurements are collected and processed for DSSE; the interval $h$ can be small – even on the order of seconds – if the fast-acting measurement capabilities of $\mu$PMUs and distributed energy resources are leveraged. In particular, a (subset of) the following quantities are measured at every time $t_{k}$, $k\in\mathbb{N}$:

• •

${\bf y}_{v}^{(k)}\in\mathbb{R}^{2M_{\bf v}\times 1}$: measurements of the line-to-ground voltages at all the phases of the nodes ${\cal M}_{\bf v}$. The measurement model for the $\mu$PMUs is ${\bf y}_{v}^{(k)}={\bf z}_{{\cal M}_{\bf v}}^{(k)}+{\bf n}_{v}^{(k)}$ with the measurement noise ${\bf n}_{v}^{(k)}$.

• •

${\bf y}_{u}^{(k)}\in\mathbb{R}^{2M_{\bf u}\times 1}$: measurements of the net active and reactive powers from wye and/or delta connections at nodes ${\cal M}_{\bf u}$. The measurement model is ${\bf y}_{u}^{(k)}={\bf u}_{{\cal M}_{\bf u}}^{(k)}+{\bf n}_{u}^{(k)}$ with the measurement noise ${\bf n}_{u}^{(k)}$.

Using (3), the measurement equation per time $t_{k}$ is:

[TABLE]

where: ${\bf J}_{v}$ is a suitable permutation matrix which selects rows of ${\bf z}$ to form ${\bf z}_{{\cal M}_{\bf v}}$; and, similarly for ${\bf J}^{Y}$ and ${\bf J}^{\Delta}$, selecting the measured loads. Recall that ${\bf M^{Y}}^{(k)}$ and ${\bf M^{\Delta}}^{(k)}$ are time-variant is the last iterate ${\bf z}^{(k)}$ is utilized as the voltage profile for the fixed-point linearization. In our algorithmic development later on, ${\bf z}^{(k)}$ is approximated by its estimated value $\widehat{{\bf z}}^{(k)}$ at each time step.

Remark 1** (Heterogeneous measurements).**

Although we model the DDSE problem for only voltage phasor data from PMUs and power data from pseudo-measurements or smart meters, it can be generalized to encompass a variety of measurements available in distribution systems, as reviewed in [3]. First, the zero-injection constraint is always satisfied by model (4) as only load buses are included by the power variables. Second, as detailed in [16, 24], the linearized model (3) is generalizable to voltage magnitude, line currents, and power flows. Accordingly, real-time measurements of these variables can be included in (4). To fit various sensing frequency of different measurements, one can set $h$ as the fastest sampling time (typically from PMUs) and maintain the values of other slower measurements until new datum arrives.

Next, consider rewriting the measurement model (4) in the following compact form:

[TABLE]

where ${\bf y}^{(k)}:=[({\bf y}_{v}^{(k)})^{\textsf{T}},({{\bf y}_{u}^{Y}}^{(k)})^{\textsf{T}},({{\bf y}^{\Delta}_{u}}^{(k)})^{\textsf{T}}]^{\textsf{T}}$, $\bar{\mathbf{m}}:=[\mathbf{m}_{v}^{\textsf{T}},\mathbf{0}^{\textsf{T}},\mathbf{0}^{\textsf{T}}]^{\textsf{T}}$, ${\bf n}^{(k)}:=[({\bf n}_{v}^{(k)})^{\textsf{T}},({{\bf n}^{Y}_{u}}^{(k)})^{\textsf{T}},({{\bf n}^{\Delta}_{u}}^{(k)})^{\textsf{T}}]^{\textsf{T}}$, and

[TABLE]

The data $\mathbf{y}_{u}$ is usually collected or generated at much lower quality compared to the high-resolution $\mu$PMU data. In fact, $\mathbf{y}_{u}$ is either collected from meters or generated from historic load information. Thus, to account for different granularities and precisions, we define the instantaneous error mismatch loss function $\ell^{(k)}({\bf z})$ at time $t_{k}$ as:

[TABLE]

with $\mathbf{G}_{v}^{(k)}:=[{\bf J}_{v}{\bf M^{Y}}^{(k)},{\bf J}_{v}{\bf M^{\Delta}}^{(k)}]$, and the Huber loss function \mathsf{H}({\mbox{\boldmath\epsilon}}):=\sum_{i}\mathsf{H}_{i}(\epsilon_{i}) written as

[TABLE]

where $\delta>0$ is a positive parameter determined by the load data quality. The Huber loss function is utilized to reject possible outliers, or down-weight data with substantial measurement errors [25].

Remark 2** (Weighted error objective).**

To accommodate varying data quality, non-uniform weights can be assigned to different types of measurements. With the high accuracy of voltage data from PMUs, one can use a large positive weight for the voltage error term in (4). Note that this weighted error objective does not affect the problem structure and thus is not included specifically by the algorithmic developments.

Using the error mismatch $\ell^{(k)}$, we formulate the following state estimation problem at time $t_{k}$:

[TABLE]

where $r^{(k)}({\bf u})$ is a (possibly time-varying) regularization function that renders the overall cost function $f^{(k)}({\bf u})$ globally strongly convex. Notice further that the problem (12) is unconstrained; however, possible prior information on the minimum and maximum values of the vector ${\bf u}$ can be naturally incorporated in the proposed approach.

Problem (12) models a time-varying state estimation task under streaming of measurements, and it implicitly defines an optimal trajectory $\{{\bf u}^{(k,\star)}\}_{k\in\mathbb{N}}$ for the estimation task. Accordingly, the objective here is to track $\{{\bf u}^{(k,\star)}\}_{k\in\mathbb{N}}$ by processing the incoming measurements in real time. One way to obtain $\{{\bf u}^{(k,\star)}\}_{k\in\mathbb{N}}$ is to solve the problem $\mathsf{P}^{(k)}({\bf u})$ to convergence (i.e., a batch solution) at each time step $t_{k}$. However, in a real-time setting with an asynchronous streaming of measurements, a batch solution of (12) might not be achievable within an interval $h$ due to underlying communication and computational complexity requirements. Thus, the objective of the paper is to develop an algorithmic solution to generate a sequence $\{{\bf u}^{(k)}\}$ of approximate optimizers for the time-varying problem $\{\mathsf{P}^{(k)}({\bf u})\}$, which eventually converges to the optimization trajectory $\{{\bf u}^{(k,\star)}\}$. Accordingly, the next section presents a running prediction-correction method to solve (12) in real time. Before doing so, a remark on the strong convexity of $\{\mathsf{P}^{(k)}({\bf u})\}$ is in order.

Strong convexity can ensure that the optimizer trajectory for the sequence of problems $\{\mathsf{P}^{(k)}({\bf u})\}$ is unique. In addition, it allows us to establish convergence for the proposed online algorithms. Regarding the regularization function, possible choices are exemplified next:

i) $r^{(k)}({\bf u})=\frac{a}{2}\|{\bf u}-{\bf u}_{\textrm{pr}}^{(k)}\|_{2}^{2}$, $a>0$ and where ${\bf u}_{\textrm{pr}}^{(k)}$ is a priori guess on the load profile;

ii) $r^{(k)}({\bf u})=\frac{a}{2}\|{\bf u}\|_{2}^{2}$, and it is time invariant;

iii) if the Hessian of the cost function is available, or, a subspace tracking method is in place, the regularization function can be set to $r^{(k)}({\bf u})=\frac{a}{2}{\bf u}^{\textsf{T}}{\bf Q}^{(k)}{\bf u}$, where the positive eigenvalues of the matrix ${\bf Q}^{(k)}$ are in the null space of the Hessian of $\ell^{(k)}({\bf u})$. The options i) and ii) would involve a deviation from optimal solutions that one would have obtained by minimizing $\ell^{(k)}({\bf u})$, with the magnitude of a possible deviation dependent on the parameter $a$; see [26]. The option iii) would not perturb the optimal solution, but the overall solution would incur a higher computational complexity.

III Dynamic state estimation

We consider First-Order Prediction-Correction (FOPC) method [18] to solve the DDSE problem at hand. Inspired by Kalman filtering approaches, prediction-correction approaches allow one to solve a broad class of time-varying convex optimization objectives in a dynamic setting by involving two stages: a prediction phase and a correction phase. In the prediction phase, the algorithm attempts to approach the optimal solution of the next time period (without new observations) by tapping into intrinsic temporal correlations of the cost function; on the other hand, in the correction phase the predicted vector is corrected using the latest measurement. This mechanism improves the response time to external dynamics and shows a good convergence result when the objective function changes smoothly over time, which is the case in the DDSE setting.

III-A FOPC algorithm

Consider a continuously time-varying unconstrained optimization with objective $f(\bm{x};t)$ to model external dynamics, along with the sampled counterpart $\{f^{(k)}:=f(\bm{x};t_{k})\}$ for $t_{k}=kh$, $k\in\mathbb{N}$. The goal is to produce a trajectory $\{\bm{x}^{(k)}\}$ such that $\bm{x}^{(k)}\approx\bm{x}^{*}(t_{k})$, where $\bm{x}^{*}(t_{k})$ denotes the optimal solution at time $t_{k}$. In order to predict the solution at time $t_{k+1}$, a strategy is to find $\bm{x}^{(k+1|k)}$ that satisfies the condition

[TABLE]

with $\gamma\in[0,1]$. Varying $\gamma$, this condition imposes optimality ($\gamma=1$), or the fact that the estimate $\bm{x}^{(k+1|k)}$ is no worse than $\bm{x}^{(k)}$ in terms of suboptimality ($\gamma=0$) even when the function changes. The choice $\gamma=1$ combines moving towards the optimizer while moving with the varying objective function; the choice $\gamma=0$ represents a rigid motion with the objective.

Notwithstanding the choice of $\gamma$, condition (13) cannot be computed at time $t_{k}$ without information about $f^{({k+1})}$. Instead, consider the following Taylor approximation:

[TABLE]

By solving this equation, we have a recursion of the form:

[TABLE]

which (as anticipated) combines a Newton’s step with a rigid motion with the objective function.

Concerning over the cost of computing the inverse of the Hessian in a possibly small time interval $h$, the FOPC method further involves a first-order update to solve (14). Specifically, the prediction solution of (14) is sought by constructing an equivalent quadratic optimization as follows

[TABLE]

where

[TABLE]

Thus, one can replace the update (15) with the gradient descent solution for (17); i.e., each iteration $p$ is given by:

[TABLE]

where integer $P$ is the number of prediction steps that one can afford within an interval $h$, and $\alpha>0$ the stepsize to be designed later on. Hence, the predicted solution is set to $x^{(k+1|k)}=\widehat{\bm{x}}_{P}$.

Initializing at $x^{(k+1|k)}$, the correction phase further involves $C$ first-order gradient steps; that is, for each $c=0,\ldots,C-1$

[TABLE]

with the stepsize $\beta>0$. The final corrected estimate is set as $\bm{x}^{(k+1)}=\widehat{\bm{x}}_{C}$. The number of steps $P$ and $C$ for the prediction and correction stages, respectively, are selected based on the duration of the interval $h$.

The complete running FOPC algorithm for the DDSE problem (12) is tabulated as Algorithm 1. Selection strategies for $P,C$ as well as the stepsizes $\alpha$ and $\beta$ will be elaborated in the next subsection. Notice that, once the prediction $\widehat{{\bf u}}^{(k|k-1)}$ and the corrected estimate $\widehat{{\bf u}}^{(k)}$ are obtained, the voltage vectors can be readily calculated using the linearized model (3). Alternatively, given the estimate of the power injections, the voltages can be calculated by solving the AC power flow equations (see, e.g., the fixed-point power flow method in [16]).

III-B Online tracking results

This subsection describes in which sense the sequence of approximate optimizers $\{\widehat{{\bf u}}^{(k)}\}_{k\in\mathbb{N}}$ generated by Algorithm 1 tracks the sampled solution trajectory ${\bf u}^{(*,k)}$. To this end, we will adapt some of the results of [18] to the DDSE problem.

The following assumptions on the cost function of (12) and its time variations are presupposed.

AS 1**.**

The cost function $f^{(k)}({\bf u})$ is $\nu$-strongly convex for all $k\in\mathbb{N}$. The regularized function $r^{(k)}({\bf u})$ is $L_{r}$-strongly smooth for all $k\in\mathbb{N}$.

AS 2**.**

The time variation of the gradient of the cost $f^{(k)}({\bf u})$ is upper bounded as

[TABLE]

Assumption 1 guarantees that the solution trajectory is unique. This is why one may need a regularizer $r^{(k)}({\bf u})$. As explained in Section II, this assumption can be verified by utilizing a suitable regularization function. Assumption 2 makes sure that the functional changes are bounded and thus it is possible to bound the errors arising from a time-varying problem. Note that, in the sequel, the exact value $C_{0}$ will not be needed for determining the parameters such as stepsizes, but we will only need the bounded condition.

First of all, we show that the Hessian of the cost function is lower and upper bounded uniformly in time, as follows.

Proposition 1**.**

Under Assumption 1, the Hessian $\nabla_{{\bf u}{\bf u}}f^{(k)}({\bf u})$ of the cost function (12) is lower and upper bounded uniformly in time as

[TABLE]

Proof:

The lower bound follows directly from Assumption 1. The upper bound follows by direct computation. The Hessian of $f^{(k)}({\bf u})$ is bounded as

[TABLE]

By properly defining $L$ as the upper bound of the right-hand term, the Hessian is upper bounded. ∎

We are now ready for the online tracking result.

Theorem 1** (Adapted from Theorem 3 of [18]).**

Consider the sequence $\{\widehat{{\bf u}}^{(k)}\}_{k\in\mathbb{N}}$ generated by Algorithm 1, and let Assumptions 1-2 hold true. Let ${\bf u}^{(*,k)}$ be the optimizer of (12) at time $t_{k}$. Choose stepsizes $\alpha$ and $\beta$ as

[TABLE]

and define the following non-negative quantities

[TABLE]

Further, select the number of correction steps $C$ in a way that

[TABLE]

Then, the sequence $\{\|{\bf u}^{(k)}-{\bf u}^{(*,k)}\|\}_{k\in\mathbb{N}}$ converges linearly with rate $\tau_{0}$ to an asymptotic error bound, and

[TABLE]

Theorem 1 asserts that the sequence $\{\widehat{{\bf u}}^{(k)}\}_{k\in\mathbb{N}}$ generated by Algorithm 1 converges to (and tracks) the sampled solution trajectory ${\bf u}^{(*,k)}$ up to an asymptotic bound. This bound is linearly related to the sampling period $h$ and exponentially decreasing with $C$. Theorem 1 requires standard bounds on the stepsizes (27) and a condition on the number of prediction and correction steps (29).

If one chooses the parameter $\gamma=0$ (i.e., rigid motion only), then the second condition boils down to

[TABLE]

Since both $\varrho_{\mathrm{P}}$ and $\varrho_{\mathrm{C}}$ are less than $1$ by construction, the condition (31) is not very restrictive. In fact, fixing a level of prediction $\bar{P}$, then the number of correction steps one has to perform is

[TABLE]

For reasonable values such as $\varrho_{\mathrm{C}}=\varrho_{\mathrm{P}}=0.8$ and $\bar{P}=4$, it implies that $C\geq 3$ is sufficient. When $\gamma>0$, Condition (29) gets more restrictive, while the tracking error accuracy may benefit from the Newton’s step in (15).

The choice of $\gamma$ is important to trade-off convergence region, requirements for prediction and correction steps, and conditioning on the measurement matrix. On one hand, if $\gamma=0$, then the convergence region is bigger, Condition (29) is less restrictive, which is good when one can afford only a small number of prediction and correction steps (in a fast sampling scenario); however $\gamma=0$ could be more affected by a badly conditioned measurement matrix. On the other hand, if $\gamma=1$, you get a Newton step in the prediction that helps in case the measurement matrix is badly conditioned, the convergence region gets smaller, and the number of required prediction and correction steps gets higher.

IV Numerical Simulations

Numerical tests have been performed using the IEEE 37-bus and 123-bus test feeders (see e.g., [27] for a description of the feeders) on a standard laptop with Intel*®* CoreTM i7-7500 CPU @2.70Hz. The 37-bus test feeder consists of 32 nodes (phases) that are connected to non-zero loads, all delta-connected. The 123-bus test feeder is a popular case with 72 non-zero load nodes and various single-, two-, and three-phase lines, with a mix of delta and wye connections. The load profiles were generated from a real dataset that the National Renewable Energy Laboratory produced from real consumption data received from a utility company in California; the data includes 55 load consumption trajectories over the course of 24-hour, at a time-resolution of 6 seconds. The trajectories for the active power are described in more details in [28]; 5 representative trajectories are shown in Fig. 1. A constant power factor of 0.95 has been postulated to create the trajectories for the reactive power. The load profiles were randomly chosen for each load node.

As for the measurement settings, the DDSE algorithm used both the load power data and PMU voltage data. We assume the load active/reactive power injections are observed everywhere in the system. Nonetheless, since they are typically collected by smart meters, the load profiles in Fig. 1 were downsampled to generate the measurement data at a slower time resolution of 10 minutes. Specifically, the power measurements were kept as the average value within every 10-minute window. Accordingly, the parameter for the Huber’s loss function in (11) was chosen to be $\delta=8{e}{-4}$, coinciding with the deviation level between the down-sampled power measurement and the actual value. Since PMUs are not typically installed everywhere, only selected buses are assumed to be equipped with high-quality voltage measurements at the same resolution as the load data. Addictive Gaussian noises with a small standard deviation of $1{e}{-5}$ were added to the actual voltage profiles to reflect the high sensing capability of PMUs. Accordingly, we set the weights of voltage measurements to be $1{e}{3}$ to match the small noise variation therein.

IV-A 37-Bus Test Feeder

This test feeder was used to compare the FOPC-based DDSE algorithm under different computational settings, with the performance in terms of both the tracking error with the instantaneous optimum and estimation error with the actual voltage state. For simplicity of implementation, the regularization term in (12) was chosen to be $r(\mathbf{u})=\frac{1}{2}\|\mathbf{u}\|_{2}$. This sufficiently small term guarantees convergence and yet does not degrade significantly the performance of the estimation problem (as we will show). To enhance the effectiveness of the prediction step under potentially ill-conditioning issue (see Sec. III), we set the parameter $\gamma=0.9$ to be close to 1.

Based on these parameter settings, the norm of the Hessian matrix as in Prop. 1 can be bounded within the interval $[3.5e4,~{}4e4]$. Thus, the stepsize parameters $\alpha$ and $\beta$ were chosen as $1e{-4}$ according to Theorem 1. Accordingly, the two parameters therein, $\varrho_{\mathrm{P}}$ and $\varrho_{\mathrm{C}}$, are around $0.65$, both less than 1. To satisfy the convergence condition in (29), it is sufficient to have $C=5$ steps of correction even if no prediction is performed ($P=0$), as $\tau_{0}=0.8$ in this case.

IV-A1 Fixed number of correction steps

We first show the advantage to add prediction phase before new measurements are processed. As mentioned in Sec. I, FOPC can perform a prediction of the state while waiting for the measurement to be transmitted from the PMUs; once the measurement is received, a fixed number of correction steps can be performed. Thus, this test uses a fixed $C=5$ steps of correction and compares the results from $P=0,5,10$ steps of prediction.

Fig. 2 plots the relative tracking error with the instantaneous optimal state ${\bf u}^{(*,k)}$. Clearly, a larger number of $P$ does help the tracking of the optimal solution, with more noticeable change in mismatch error from $P=0$ to $P=5$ steps of prediction. It also makes the relative error trajectory more quickly to reach the steady-state level of below $0.1$. Hence, the prediction phase has been shown to improve the tracking error performance for the time-varying DDSE problem.

Furthermore, we compare the relative estimation errors for both the system-wide power states and corresponding voltage outputs as compared to the ground-truth values, as plotted in Fig. 3. As for the power states, the relative estimation error increases from the tracking error level. This could be due to the approximation error of the linearized model adopted by the objective function and the measurement error from low-resolution power data. Meanwhile, the effectiveness of prediction phase is still evident in improving the estimation error and convergence rate. More interestingly, the voltage estimation error is very minimal at the level of below 1e-3 for all scenarios, corroborating that the PMU voltage data is instrumental for recovering the feeder voltage profile even under highly uncertain power measurements.

IV-A2 Fixed computational time

This test compares the FOPC performance under a total computational time constraint. The computational time in the prediction phase is mainly spent on a one-shot computation of Hessian $\nabla_{{\bf u}{\bf u}}f^{(k)}(\widehat{{\bf u}}^{(k)})$ in (20), while that of the correction phase grows linearly with $C$ as it needs to compute the gradient in every step. We choose two sets of $(P,C)$ values: $(8,3)$ and $(0,6)$, both taking roughly a total of $0.3$ms per iteration. For the case of $(P,C)=(8,3)$, we verified again the convergence condition in (29) with $\tau_{0}\approx 0.8$.

Fig. 4 and Fig. 5 plot the relative tracking and estimation error trajectories for the scenarios, respectively, as in the last test. Interestingly, even with smaller $C$ to compensate for the one-shot Hessian computation, the improvement of including the prediction phase can be demonstrated. Both the tracking and estimation errors for the power states are lower for the case of $(P,C)=(8,3)$. The voltage error is again very small thanks to the high-quality PMU data.

IV-A3 Varying number of PMUs

We further compare the performance when different number of PMUs are installed in the system. Fig. 6 plots the relative error of estimating both the ground-truth power and voltage variables. Clearly, more high-quality voltage data can significantly improve the estimation error performance. For the 37-bus case, it seems that 3 PMUs are sufficient for estimating the system-wide voltage, as there is no noticeable improvement with 5 PMUs. This is not necessarily the case when estimating the power state, due to the high uncertainty of power measurements.

IV-B 123-Bus Test Feeder

Last, we tested the FOPC method on the 123-bus case to demonstrate its scalability. The parameter settings follow from those in the 37-bus tests, with the estimation error comparisons for fixed $C$ given in Fig. 7. Similar results have been observed, corroborating the improvement of including the prediction phase and high-quality voltage data from PMUs.

V Conclusions

This paper presented a distribution state estimation algorithm that can dynamically incorporate fast and accurate PMU voltage data. The first-order prediction-correction (FOPC) algorithm is proposed to solve the time-varying optimization problem of DDSE using a linearized power flow model. Compared to existing recursive updates, the FOPC iterations are computationally simple and require no specific modeling of system transition, suitable for the time-critical DDSE problem where the load dynamics is difficult to model. Numerical tests have shown that the data-driven prediction phase of FOPC is effective in reducing the mismatch error in tracking the power state variable. With the availability of high-quality voltage data, the voltage estimation error is significantly small. Future work includes exploring more diverse types of distribution system measurements and large-scale system validations using real data.