Bad Data Analysis

One of the essential state estimation routines is the bad data analysis, which follows after obtaining the weighted least-squares (WLS) estimator. Its main task is to detect and identify measurement errors, and eliminate them if possible. This is usually done by processing the measurement residuals [5, Ch. 5], and typically, the largest normalized residual test is used to identify bad data. The largest normalized residual test is performed after we obtained the solution of the state estimation in the repetitive process of identifying and eliminating bad data measurements one after another [24].

Additionally, the Chi-squared test, which can precede the largest normalized residual test, serves to detect the presence of bad data and quickly determine if the largest normalized residual test should be performed [5, Sec. 5.4].

To initiate the process, let us construct the PowerSystem type and formulate the AC model:

system = powerSystem()

addBus!(system; label = 1, type = 3, active = 0.5)
addBus!(system; label = 2, type = 1, reactive = 0.3)
addBus!(system; label = 3, type = 1, active = 0.5)

@branch(resistance = 0.02, susceptance = 0.04)
addBranch!(system; label = 1, from = 1, to = 2, reactance = 0.6)
addBranch!(system; label = 2, from = 1, to = 3, reactance = 0.7)
addBranch!(system; label = 3, from = 2, to = 3, reactance = 0.2)

addGenerator!(system; label = 1, bus = 1, active = 3.2, reactive = 0.2)

acModel!(system)

Following that, we introduce the Measurement type, which represents a set of measurement devices $\mathcal M$:

monitoring = measurement(system)

addWattmeter!(monitoring; label = "Watmeter 1", bus = 3, active = -0.5)
addWattmeter!(monitoring; label = "Watmeter 2", from = 1, active = 0.2)
addWattmeter!(monitoring; label = "Watmeter 3", bus = 3, active = 3.1)

addVarmeter!(monitoring; label = "Varmeter 1", bus = 2, reactive = -0.3)
addVarmeter!(monitoring; label = "Varmeter 3", from = 1, reactive = 0.2)

addPmu!(monitoring; label = "PMU 1", bus = 1, magnitude = 1.0, angle = 0.0)
addPmu!(monitoring; label = "PMU 2", bus = 3, magnitude = 0.9, angle = -0.2)

Let the WLS estimator $\hat {\mathbf x}$ be obtained by solving the AC state estimation:

analysis = gaussNewton(monitoring)
stateEstimation!(analysis)

Chi-Squared Test

Next, we perform the chi-squared test to check for the presence of outliers:

chi = chiTest(analysis; confidence = 0.96)

At this stage, JuliaGrid uses the objective value obtained from the AC state estimation:

\[ f(\hat{\mathbf x}) = \sum_{i=1}^k\cfrac{[z_i - h_i(\hat{\mathbf x})]^2}{v_i},\]

This value is stored in the ChiTest type as:

julia> chi.objective69608.89418770038

Next, retrieve the value from the Chi-squared distribution corresponding to the detection confidence and $(k - s)$ degrees of freedom, where $k$ is the number of measurement functions and $s$ is the number of state variables. This provides the value of $\chi^2_{p}(k - s)$:

julia> chi.treshold10.025519286562865

Then, the bad data test can be defined as:

\[ f(\hat{\mathbf x}) \geq \chi^2_{p}(k - s).\]

If the inequality is satisfied, bad data is suspected in the measurement set:

julia> chi.detecttrue

Largest Normalized Residual Test

As observed from the Chi-squared test, bad data is present in the measurement set. We then perform the largest normalized residual test to identify the outlier and remove it from the measurements:

outlier = residualTest!(analysis; threshold = 4.0)

In this step, we employ the largest normalized residual test, guided by the analysis outlined in [5, Sec. 5.7]. To be more precise, we compute all measurement residuals based on the obtained estimate of state variables:

\[ r_i = z_i - h_i(\hat {\mathbf x}), \;\;\; i \in \mathcal M.\]

The normalized residuals for all measurements are computed as follows:

\[ \bar{r}_i = \cfrac{|r_i|}{\sqrt{S_{ii}\Sigma_{ii}}} = \cfrac{|r_i|}{\sqrt{C_{ii}}}, \;\;\; i \in \mathcal M.\]

In this equation, we denote the diagonal entries of the residual covariance matrix $\mathbf C \in \mathbb{R}^{k \times k}$ as $C_{ii} = S_{ii}\Sigma_{ii}$, where $S_{ii}$ is the diagonal entry of the residual sensitivity matrix $\mathbf S$ representing the sensitivity of the measurement residuals to the measurement errors.

The subsequent step involves selecting the largest normalized residual, and the $j$-th measurement is then suspected as bad data and potentially removed from the measurement set $\mathcal{M}$:

\[ \bar{r}_{j} = \text{max} \{\bar{r}_{i}, \forall i \in \mathcal{M} \}.\]

Users can access this information using the variable:

julia> outlier.maxNormalizedResidual260.80802083468797

If the largest normalized residual, denoted as $\bar{r}_{j}$, satisfies the inequality:

\[ \bar{r}_{j} \ge \epsilon,\]

the corresponding measurement is identified as bad data and subsequently removed. In this example, the bad data identification threshold is set to $\epsilon = 4$. Users can verify the satisfaction of this inequality by inspecting:

julia> outlier.detecttrue

This indicates that the measurement labeled as:

julia> outlier.label"Watmeter 3"

is removed from the measurement set and marked as out-of-service.

Subsequently, we can solve the system again, but this time without the removed measurement:

analysis = gaussNewton(monitoring)
stateEstimation!(analysis)

Following that, we check for outliers once more:

outlier = residualTest!(analysis; threshold = 4.0)

To examine the value:

julia> outlier.maxNormalizedResidual1.0990106756466569

As this value is now less than the threshold $\epsilon = 4$, the measurement is not removed, or there are no outliers. This can also be verified by observing the bad data flag:

julia> outlier.detectfalse

Residual Covariance Matrix

In AC state estimation, the residual covariance matrix $\mathbf C$ is given by:

\[ \mathbf C = \mathbf S \bm \Sigma = \bm \Sigma - \mathbf J (\hat {\mathbf x}) [\mathbf J (\hat {\mathbf x})^T \bm \Sigma^{-1} \mathbf J (\hat {\mathbf x})]^{-1} \mathbf J (\hat {\mathbf x})^T,\]

while for DC state estimation and state estimation using only PMUs, it is computed as:

\[ \mathbf C = \mathbf S \bm \Sigma = \bm \Sigma - \mathbf H [\mathbf H^T \bm \Sigma^{-1} \mathbf H]^{-1} \mathbf H^T.\]

It is important to note that only the diagonal entries of $\mathbf C$ are required. The main computational challenge lies in computing the inverse on the right-hand side of the equations. The JuliaGrid package employs a computationally efficient sparse inverse method, obtaining only the necessary elements.