Uncertainty Error Modeling for Non-Linear State Estimation With Unsynchronized SCADA and $\mu$PMU Measurements

Classical power system SE follows a Weighted Least Squares (WLS) framework [8]. A power system with $n$ buses and $d$ measurements can be modeled as a set of non-linear algebraic equations in the following measurement model:

where $\mathbf{z}\in\mathbb{R}^{1\times d}$ is the measurement vector, $\mathbf{x}\in\mathbb{R}^{1\times N}$ is the state variables vector, $h:\mathbb{R}^{1\times N}\rightarrow\mathbb{R}^{1\times d}$ is a continuous non-linear differentiable function, and $\mathbf{e}\in\mathbb{R}^{1\times d}$ is the measurement error vector. Each measurement error $e_{i}$ is assumed to follow a zero mean Gaussian distribution. $d$ is the number of measurements while $N=2n-1$ is the number of unknown state variables, i.e., the complex voltages at each bus.

In the traditional WLS approach, the best state vector estimate in (1) is determined by minimizing the weighted norm of the residual, represented with the cost function $J(\mathbf{x})$:

where $\mathbf{W}$ is the inverse covariance matrix of the measurements, otherwise known as the weight matrix. The solution to (2) is found through the iterative Newton-Raphson method. The linearized form of (1) then becomes:

where $\mathbf{H}=\frac{\partial h}{\partial\mathbf{x}}$ is the Jacobian matrix of $h$ at the current state estimate $\mathbf{x}^{*}$, $\Delta\mathbf{z}=\mathbf{z}-h(\mathbf{x}^{*})=\mathbf{z}-\mathbf{z}^{*}$ is the measurement vector correction and $\Delta\mathbf{x}=\mathbf{x}-\mathbf{x}^{*}$ is the state vector correction. The WLS solution can further be understood geometrically [10] as the projection of $\Delta\mathbf{z}$ onto the Jacobian space $\mathfrak{R}$($H$) by a linear projection matrix $\mathbf{K}$, i.e. $\Delta\hat{\mathbf{z}}=\mathbf{K}\Delta\mathbf{z}$:

The SE can further be interpreted by visualizing the geometrical position of the measurement error in relation to the Jacobian range space $\mathfrak{R}$($H$). Through decomposing the measurement vector space into a direct sum of $\mathfrak{R}$($H$) and $\mathfrak{R}($H$)^{\perp}$, it is possible to also decompose the measurement error vector $\mathbf{e}$ into two components:

The component $\mathbf{e_{D}}$ is the detectable error, which is equivalent to the residual in the classical WLS model. The component $\mathbf{e_{U}}$ is the undetectable component of the residual. $\mathbf{e_{D}}$ lies in the space orthogonal to the Jacobian range space, while $\mathbf{e_{U}}$ is hidden in the Jacobian space [8]. To formalize the impact of the undetectable component $\mathbf{e_{U}}$, the measurement Innovation Index ($II$) is introduced [11]:

A low $II$ indicates that a large component of the error is not reflected in the residual. The composed measurement error ($CME$) can then be formulated in terms of the residual and $II$, after which the normalized $CME^{N}$ is obtained:

where $\sigma_{i}$ is the $i$-th measurement standard deviation.

The SE yields $CME$ values based on measurements taken throughout the power system. Real measurements in this work are defined as those obtained from SCADA and $\mu$PMUs. In addition, per the methodology in [12], synthetic measurements (SM) are artificially created in low redundancy areas considering measurement $II$ and $n$-tuple of critical measurements. This serves to maintain a high global redundancy level ($GRL$), i.e., the number of measurements divided by the number of state variables, by increasing the $\chi^{2}$ degrees of freedom, improving the reliability of the $\chi^{2}$ hypothesis test [8].

$\chi^{2}$ hypothesis testing is used for bad data detection during gross error (GE) analytics. The $CME$-based objective function value (8) is compared to a $\chi^{2}$ threshold. The value of this threshold depends on a chosen probability $p$ and the degrees of freedom $d$, taken to be the number of measurements fed to the SE:

If the value of $J_{CME}$ is greater than or equal to the $\chi^{2}$ threshold, then a GE is detected. For this work, a two-step SE approach is used, illustrated in Fig. 1. The classical first step presented in [13] uses an empirical approach wherein all measurements are weighted equally proportional to the measurement magnitude, after which the $J_{CME}(\hat{\mathbf{x}})$ is compared to $\chi^{2}_{d,p}$. Only if $J_{CME}(\hat{\mathbf{x}})\geq\chi^{2}_{d,p}$ does the procedure advance to the second SE step, during which the SE is repeated—this time with meter precision values for each measurement type per state-of-the-art methodologies [3]. This serves to strike a balance between GE detectability and deciphering the measurement(s) in error. Next, following [14], bad data is identified through analysis of the $CME^{N}$. This work seeks to improve upon the first-step of the classical two-step SE by replacing the empirical assumption with a continually updating model of the unsynchronized measurement uncertainty error.

Due to the asynchronicity between $\mu$PMU and SCADA measurements, as well as the co-asynchronicity among SCADA, load dynamics cause SCADA measurements to become obsolete, thus providing an inaccurate snapshot of the system state. This uncertainty error must be included into the error vector $\mathbf{e}$ from (1) to better provide accurate state variables, error detection, and GE analytics. The drift in the system state due to load dynamics can be modelled as Ornstein–Uhlenbeck (OU) stochastic processes, which have been used to successfully reflect time-varying load profiles [15]. OU load modelling relies on the assumption that, although load demand profiles do fluctuate over time, they rarely deviate sharply from their previous value. Further, the OU is a mean-reverting process, making it an appropriate model for ascribing the uncertainty error to outdated measurement variance only.

A stochastic differential equation is used to represent the $i$-th load apparent power, $S_{i}(t)=P_{i}(t)+jQ_{i}(t)$, where $P_{i}(t)$ and $Q_{i}(t)$ are the $i$-th load active and reactive power respectively:

where $\theta_{i}$ $>$ 0 is the decay rate controlling to what degree the load varies from its initial value, $\sigma_{i}$ $>$ 0 controls the scale of the noise, and $W_{i}(t)$ denotes a complex Wiener process whose sample paths follow a continuous Gaussian process [16]. Statistical analysis of $\mu$PMU data can be used to estimate these parameters [15]. The load variations can then be expressed as a discrete time process for every $\Delta$t $\mu$PMU sample [6]:

where $\tau$ $<$ t is the time when unsynchronized SCADA is received, $t=k\Delta t$ and $\tau=q\Delta t$ for $k$ and $q$ time indices, and $\zeta_{i}$ is a complex random variable of zero mean and variance $\frac{\sigma_{i}^{2}}{2\theta_{i}}(1-e^{-2\theta_{i}(t-\tau)})$. Dynamic load uncertainty can substantially increase the $J_{CME}(\hat{\mathbf{x}})$ in the absence of actual measurement error, especially when measurements are unsynchronized, as illustrated in Fig. 2. The variance of the discrete-time Gaussian variable $S_{i}[k]$ can then be found recursively by:

where $\gamma=e^{-2\theta_{i}\Delta t}$. The updated variances of the unsynchronized SCADA measurements then constitute the time-varying inverse covariance matrix of the measurements, $\mathbf{W}(t)$. The improvement in the $J_{CME}(\hat{\mathbf{x}})$ due to the updating $\mathbf{W}(t)$ is observed in Fig. 3, as the $J_{CME}(\hat{\mathbf{x}})$ should not exceed the $\chi^{2}$ threshold due to load dynamics only.

The classical WLS model in (1) does not consider the possibility of errors attributed to dynamic load measurement uncertainty. This model can be augmented by considering, for any $i$-th measurement, a noiseless true measurement $z_{i}^{true}[k]$ at time $k$, which is related to a previous value $z_{i}^{true}[q]$ plus random load variation $\Delta z_{i}^{true}[k]$, given by (12), whereas the noisy and out-of-date unsynchronized measurement, $z_{i}^{meas}[q]$, is represented in (13):

where $e_{i}[q]$ is the measurement noise of the same Gaussian distribution as in (1). At a given instance $k$ when the first step of the non-linear SE is performed, the precise value of $z_{i}^{true}[k]$ is unobtainable; the best available measurement in time is $z_{i}^{meas}[q]$. However, it is possible to approximate $z_{i}^{meas}[k]$ by assuming that $\Delta z_{i}^{true}[k]$ is attributed only to load variations:

The model in (14) is used for approximating the unsynchronized SCADA measurement variances, which are used to update the time-varying $\mathbf{W}(t)$ weight matrix.

The proposed work also seeks to formalize its assertion of preserving measurement fidelity. The work in [6] uses an OU process load modeling approach with variance update equations, however it relies on linearizing SCADA and AMI measurements. Unlike the proposed work, this linearization process introduces yet an additional error component exclusive to the LSE approach, which this work calls the linearization truncation error $\epsilon^{l}$. To formalize this additive error, the Taylor series of the LSE measurement model can be developed as:

By setting (15) equal to the standard SE measurement model $z_{i}=h_{i}(x)+e_{i}$, one obtains:

By defining the linearized form of $z_{i}$ as $z_{i}^{l}=h_{i,0}+\frac{\partial h_{i}}{\partial x}\Delta x$, one finally obtains:

where $e_{i}$ is the measurement error and $\epsilon_{i}^{l}$ is the LSE additive truncation error, which contributes to SE inaccuracies. The proposed model does not require measurement function linearization, and thus preserves measurement fidelity by avoiding the error component $\epsilon_{i}^{l}$. Further, in the referenced LSE approach, additional error is introduced in the calculation of linearized measurement covariances, which must be approximated using the truncated Taylor series of their functions [6]. The proposed work thus further preserves measurement fidelity by circumventing measurement covariance approximation.