Dynamic Graph-Based Anomaly Detection in the Electrical Grid

Transmission grids, in order to be observable, have a large number of RTUs and PMUs installed. For the anomaly detection algorithm developed in this paper, processing the large volume of redundant data for anomaly detection is unnecessary and computationally prohibitive. This is because each sensor predominantly captures the relative information of its neighboring sensors as well. Therefore, to create a proper input for anomaly detection, pre-processing techniques can be deployed: Principal Component Analysis (PCA)[46] creates a low-dimensional representation by extracting uncorrelated directions; projection pursuit[47] reduces the input to a low-dimensional projected time series that optimizes the kurtosis coefficient; Independent Component Analysis (ICA)[48] identifies a subset of independent variables. Alternative techniques like cross-correlation analysis[49] also help create a low-dimensional input. A detailed survey regarding dimensionality reduction can be found in[50]. Rather than transforming the redundant input, other algorithms for sensor placement  [51] are also applicable, by suggesting the best several locations of sensors to be installed and providing observability. [30] has shown the selection of a small number of sensor locations with a provably near-optimal probability of detecting an anomaly.

Many graph distance/similarity measures have been proposed in the past that relate to anomaly detection in dynamic graphs. A survey of such measures can be found in [52] [53]. Most of these measures fit in one of the following categories:

1. 1.

   Graph isomorphism and its generalizations: examples include Maximum Common Subgraph (MCS) distance [54], Graph Edit Distance (GED) [54], and variants of GED[55], etc.

2. 2.

   Aggregate statistical measure: preferred for measure for larger graphs, examples include diameter distance [56], clustering based measures [57], and degree distribution[58], etc.

3. 3.

   Iterative methods based on the structural similarity of local neighborhoods: This type of method exchanges node/edge similarities until convergence and computes the similarity between two graphs by coupling the similarity scores of nodes and edges [52].

4. 4.

   Other complex feature-based measures: examples include graph kernel-based similarities [59], modality distance[60], median graph distance[60], etc.

While all of the above graph distance measures have unique advantages, none of them are designed for grid-specific challenges, nor do they capture the implicit physics of the power grid graph. Section III-C will provide a more detailed discussion about grid-specific challenges and develop a novel graph distance measure to meet the needs of the grid anomaly detection.

Line Outage Distribution Factor (LODF) is a sensitivity measure of how much an outage on a line affects real power flow on other lines in the system [61]. This factor can be easily and efficiently calculated by assuming a DC power flow model with lossless lines or a linearized AC power flow model around the operating point and is commonly used to estimate the linear impact of line outage. For an outage on line $k$, LODF $d_{l}^{k}$ gives the ratio between power change $\Delta f_{l}$ on an observed line $l$ and the pre-outage real power $f_{k}$ on the outage line $k$.

Table II shows the symbols used in this paper.

We are given a dynamic graph (grid) $\mathcal{G}_{t}=(\mathcal{V}(t),\mathcal{E}(t))$ at each time tick $t$, where $\mathcal{V}(t)$ denotes the set of nodes (active grid buses), and $\mathcal{E}(t))$ denotes the set of edges (active grid branches). Also we have a fixed set of sensors $\mathcal{S}\subseteq\mathcal{V}$. Each sensor installed on node $i$ can obtain PMU or RTU measurements at each time $t$.

For each sensor on node $i$, we obtain the power flows on all lines adjacent to node $i$, as observed in [30] that using power (rather than current) provides better anomaly detection in practice. For any PMU bus $i$ and edge $e\in\mathcal{N}_{i}$, define the power w.r.t. $i$ along edge $e$ as $s_{ie}(t)=V_{i}(t)\cdot I_{e}(t)^{*}$, where ^∗ is the complex conjugate.

Topology determination: At any time $t$, the given dynamic graph $\mathcal{G}(t)$ can be considered a ‘reference topology’. The observed measurements should be consistent with the reference topology if there are no anomalous measurements and topology. In this work, we use the output from topology estimation (NTP) in the grid energy management system (EMS) as the reference topology. It is a valid reference because the operator believes that this is the best estimate of the grid topology at any time. Notably, this topology is not assumed to be perfect and accurate. Topological errors, along with false measurements, are to be detected by the proposed method.

Consider the simple power grid shown in Figure 2, which evolves over time from $\mathcal{G}_{1}$ to $\mathcal{G}_{2}$ to $\mathcal{G}_{3}$. For simplicity, assume that we have a single sensor, from which we want to detect anomalous events. How do we evaluate whether the current time point ($t=15$) is an anomaly? If the graph had not been changing, we could simply combine all past sensor values to learn a distribution of normal behavior (e.g. fitting a Gaussian distribution as in $\mathcal{N}_{t=15}$), then evaluate the current time point using this Gaussian distribution (e.g. in terms of the number of standard deviations away from the mean).

In the changing graph setting, we still want to learn a model of normal behavior ($\mathcal{N}_{t=15}$), but while taking the graph changes into account. Note that $\mathcal{G}_{2}$ and $\mathcal{G}_{3}$ are only slightly different, while $\mathcal{G}_{1}$ and $\mathcal{G}_{3}$ are very different. Hence, the sensor values coming from $\mathcal{G}_{2}$ (i.e. time $6$ to $10$) should be taken into account more highly when constructing $\mathcal{N}_{t=15}$, as compared to those from $\mathcal{G}_{1}$. Intuitively, the sensor values from $\mathcal{G}_{1}$ are drawn from a very different distribution from the current graph, and thus should not influence our learned model $\mathcal{N}_{t=15}$. In general, the more similar a graph is to the current graph, the more we should take its sensor values into account when learning our current model. This motivates the 3-step process we use:

1. 1.

   Graph Distances: Measure the distance between each past graph and the current graph.

2. 2.

   Temporal Weighting: Weight the past sensor data, where data from graphs that are similar to the current one are given higher weight.

3. 3.

   Anomaly Detection: Learn a distribution of normal behavior ($\mathcal{N}_{t}$) from the weighted sensor values, and measure the anomalousness at the current time based on its deviation from this distribution.

To further clarify the motivation and methodology, we provide an informal definition for the anomaly detection problem and a statistical definition of the anomaly pattern.

In electric grids, the target anomalies correspond to topological errors (unexpected topology changes unknown to the operator) and measurement errors. These are the types of anomalies that our method is proposed for and is likely to detect based on empirical results provided at the end of the manuscript.

In the following sections, we first introduce our domain-aware graph distance measure based on Line Outage Distribution Factors (LODF) [61]. Then, we describe our temporal weighting and anomaly detection framework, which flexibly allows for any given graph distance measure. Finally, we present an alternate distance measure that is locally sensitive, i.e., it accounts for the local neighborhood around a given sensor.

In this section, we describe our proposed graph distance measure to calculate the distance $D(\mathcal{G}_{i},\mathcal{G}_{j})$ between any pair of graphs. For ease of understanding, the rest of Section III uses the example of anomaly detection at $t=15$ in Figure 2 as an extended case study, but our approach can be easily extended to the general case.

Section II-C has provided a short review of existing graph distance measures, but these measures do not apply directly to the grid-specific challenge in this paper. For an anomaly detection algorithm to work well on the power grid applications, the choice of graph distance needs to consider problem-specific challenges along with desirable properties for an anomaly detection algorithm (scalability, sensitivity to change, and ‘importance-of-change awareness’). One grid-specific challenge is that the ideal graph distance should capture ‘grid physics’ rather than only graph structural changes. Specifically, the distance should be sensitive to the redistribution of power flow, not only the addition/deletion of nodes/edges, since the anomaly information is extracted from power flow measurements. Meanwhile, as the ‘importance-of-change awareness’ indicates, grid changes that cause big shifts in power flow (measurements) should result in larger graph distances, than changes that cause minor power impact. Unfortunately, none of the classical graph distances can capture the physics of the power flow and quantify the impact of graph change in terms of power. To handle this, this work proposes a novel design of graph distance by making use of the power sensitivity factor.

Intuitively, the goal is for our graph distance to represent redistribution of line power flow. Critical changes in topology result in large redistributions of power. Thus, the graph distance arising from a topology change should be large if the changed edges can potentially cause large amounts of power redistribution.

Hence, given two graphs $\mathcal{G}_{i}(\mathcal{V}(i),\mathcal{E}(i))$ and $\mathcal{G}_{j}(\mathcal{V}(j),\mathcal{E}(j)$ with different topology, we first define a transition state which takes the union of the two graphs:

Then the topology changes from $\mathcal{G}_{i}$ to $\mathcal{G}_{j}$ can be considered as different line deletions from their base graph $\mathcal{G}_{trans}$. For each single line deletion, e.g. line $p$, we define its contribution $x_{p}$ to graph distance by taking the average of its power impacts on all other lines as measured by LODF:

where $|\mathcal{E}(i)\cup\mathcal{E}(j)|$ denotes the cardinality of set $\mathcal{E}(i)\cup\mathcal{E}(j)$, $d_{l}^{p}$ denotes the LODF coefficient with $p$ as outage line and $l$ as observed line.

Then graph distance $D(\mathcal{G}_{i},\mathcal{G}_{j})$ is given by summing up the contributions of different line deletions from the base graph:

where $\mathcal{E}(i)-\mathcal{E}(j)=\{p|p\in\mathcal{E}(i),p\notin\mathcal{E}(j)\}$ and accordingly, $(\mathcal{E}(i)-\mathcal{E}(j))\cup(\mathcal{E}(j)-\mathcal{E}(i))$ denotes all the edge changes between the two graphs.

This definition uses LODF as a measure of the impact on power flow of the removal of line $p$. Hence, edges with high LODF to many other edges can potentially cause greater changes in power flow, and thus our graph distance measure places greater importance on these edges. Appendix A-A demonstrates the effectiveness of the LODF-based graph distance by comparing it against traditional distance measures for anomaly detection.

In this section, we assume that we are given any distance measurements $D(\mathcal{G}_{i},\mathcal{G}_{j})$ between any pair of graphs $\mathcal{G}_{i}$ and $\mathcal{G}_{j}$, and explain how to use them to assign weights to each previous sensor data. This procedure can take the LODF-based distance defined in the previous subsection as input, but also allows us to flexibly use any given graph distance measure. The proposed Temporal Weighting is given in Algorithm 1.

For the purpose of utilizing previous data from a series of dynamic graphs, Temporal Weighting plays an important role. The resulting weights directly determine how much information to extract from each previous record, thus requiring special care. Intuitively, the weights should satisfy the following principles:

• •

  The larger the distance $d_{t}$, the lower the weight $w_{t}$. This is because high $d_{t}$ indicates that time tick $t$ is drawn from a very different graph from the current one, and thus should not be given high weight when estimating the expected distribution at the current time

• •

  Positivity and Normalization: $\sum_{t}w_{t}=1,w_{t}\geq 0$

To satisfy these conditions, we use a principled optimization approach based on bias-variance trade-off. Intuitively, the problem with using data with high $d_{t}$ is bias: it is drawn from a distribution that is very different from the current one, and that can be considered a biased sample. We treat $d_{t}$ as a measure of the amount of bias. Hence, given weights $w_{1},\cdots,w_{14}$ on previous data (in Figure 2 example), the total bias we incur can be defined as $\sum_{t\in\{1,...,14\}}w_{t}d_{t}$.

We could make the bias low simply by assigning positive weights to only time points from the most recent graph. However, this is still unsatisfactory as it results in a huge amount of variance: since very little data is used to learn $\mathcal{N}_{t=15}$, the resulting estimate has high variance. Multiplying a fixed random variable by a weight $w_{t}$ scales its variance proportionally to $w_{t}^{2}$. Hence, given weights $w=[w_{1},w_{2},\cdots]$, the total amount of variance is proportional to $\frac{1}{2}w^{T}w$, which we define as our variance term.

We thus formulate the following optimization problem as minimizing the sum of bias and variance, thereby balancing the goals of low bias (i.e. using data from similar graphs) and low variance (using sufficient data to form our estimates). We formulate the problem as:

subject to

By writing out its Lagrangian function:

and applying KKT conditions, we can see the optimal primal-dual solution $(w,\lambda^{*},u^{*})$ must satisfy:

Since we have $d_{t}\geq 0$, by further manipulation we have:

Moreover, there is a unique choice of $\lambda^{*}$ such that the resulting weights $w_{t}$ sum up to $1$. This $w_{t}$ against $d_{t}$ relationship is shown in Figure 4. This result is intuitive: as $d_{t}$ increases, the resulting weight we assign $w_{t}$ decreases, and if $d_{t}$ passes a certain threshold, it becomes large enough so that any reduction in variance it could provide is more than offset by its large bias, in which case we assign it a weight of $0$.

Our Temporal Weighting algorithm is in Algorithm 2. During implementation, we adjust the relative importance of bias and variance by normalizing and scaling the graph distances (scaling factor 0.005 works well based on our empirical observation).

Having obtained our weights $w_{t}$, the remaining step is to compute our anomaly score, as shown in Algorithm 3.

We focus on 3 metrics from sensor data as indications of power system anomalies. These metrics were studied in [30] and found to be effective for detecting anomalies in power grid sensor data. In our setting, recall that for each sensor, we can obtain $\Delta s_{i}$ that contains changes of real and reactive power on the adjacent lines, over time. The 3 metrics are:

• •

  Edge anomaly metric: $X_{edge,i}(t)=\max_{l\in{E_{adj}}}\Delta s_{i,l}$ which measures the maximum line flow change among lines connected to the sensor. Let $E_{adj}$ denote the set of lines connected to sensor $i$:

• •

  Average anomaly metric: $X_{ave,i}(t)=\text{mean}\{\Delta s_{i,l}\ |\ l\in{E_{adj}}\},$ which measures the average line flow change on the lines connected to the sensor:

• •

  Diversion anomaly metric: $X_{div,i}(t)=\text{std}\{\Delta s_{i,l}\ |\ l\in{E_{adj}}\},$ which measures the standard deviation of line flow change over all lines connected to the sensor:

Intuitively, for each metric, we want to estimate a model of its normal behavior. To do this, we compute the weighted median and interquartile range (IQR)¹¹1IQR is the difference between 1st and 3rd quartiles of the distribution, and is commonly used as a robust measure of spread. of the detection metric, weighting the time points using our temporal weights $w_{1},\cdots,w_{t}$. (Weighted) median and IQR are preferred choice of distribution parameters over the mean and variance for anomaly detection since they are robust measures of central tendency and statistical dispersion (i.e. they are less likely to be impacted by outliers)[62]. We can then estimate the anomalousness of the current time tick by computing the current value of a metric, then subtracting its weighted median and dividing by its IQR. The exact steps are given in Algorithm 3.

This section focuses on a quantitative analysis of the performance of our method, through statistical error analysis.

Let $T$ denote the width of time window for analysis, then for $\forall$ sensor $s$, the anomalousness of its observation $x_{T+1}$ is evaluated based on its previous data $x_{1},x_{2},...,x_{T}$. The anomaly detection method works by assigning weights $w_{1},w_{2},...,w_{T}$ ($w_{t}\geq 0,\forall t,\sum_{t=1}^{T}w_{t}=1$) to all the previous observations and an alarm is created if $x_{T+1}$ deviates $\sum_{t=1}^{T}w_{t}x_{t}$ by a certain threshold.

Here we investigate the properties of statistical error based on the following definitions and assumptions:

We first demonstrate that the statistical error can be bounded:

Upon obtaining the error bound, here are some intuitive explanations of the upper bound being dependent on $\max_{t}w_{t}$ and $\max_{t}d^{*}_{t}$:

• •

  $\max_{t}w_{t}$: large value for this term indicates that the estimation method depends heavily on a particular prior data point with weight $w_{t}$. This can lead to overfitting and as a result higher error (upper bound) due to high variance.

• •

  $\max_{t}d^{*}_{t}$: large value for this term indicates that a prior data point from a very different distribution has been used for estimation, which can lead to higher error (upper bound) due to high bias.

In this section, we compare DynWatch against baseline anomaly detection approaches, while varying the number of sensors in the grid.

In this subsection, we seek to analyze the scalability of our DynWatch and DynWatch-Local. In reality, PJM, the largest ISO in the U.S., runs ACSE on a 28k bus model, performed every 1 min[10], thus any anomaly detection algorithm that takes significantly less than 1 min may provide valuable information to prevent wrong control decisions in real-time. The following results demonstrate the proposed method’s capability to achieve this.

Here, we generate test cases of different sizes by starting with the case2383 case and duplicating it $3,4,5,\cdots,12$ times. After each duplication, edges are added to connect each node with its counterpart in the last duplication, so that the whole grid is connected. Then for each testcase, we generate 20 dynamic grids and sensor data with 1 randomly chosen sensor and 1200 time ticks, following the same settings as the previous sub-section. Finally, we measure

Figure 9 shows that our method is fast: even on a large case with 60k+ branches, both methods took less than 2s to apply anomaly detection on all 1200 time ticks of the sensor, corresponding to an average of less than 1.7ms per time tick per sensor. The ACTIVSg25k (similar to the largest real-time network in the US) has 32k+ branches, and thus run-time of anomaly detection at each time $t$ will be significantly less than 1 min. The figure also shows that our methods scale close to linearly with the grid size.

In this section, we explore how the proposed Dynwatch algorithm can improve the performance of the standard residual-based ACSE bad-data detection (BDD) method, by testing a type of grid-specific anomaly that SE BDD is known to fail against False Data Injection Attack (FDIA).

False Data Injection Attack (FDIA)[31] is a cyber-attack in which attackers manipulate the value of measurements according to the grid physical model such that the SE outputs incorrect grid estimates while ensuring that its residual does not change by much (ideally remains unchanged).

In this experiment, we construct an attack on a 14 bus network to mislead the operator into thinking that the load reduces by 20%. For any anomalous time tick $t$, this is implemented by simulating power flow with the reduced load and generating measurements based on that.

The time-series measurement data and a comparison of anomaly scores are shown in Figure 1. Results show that the ACSE residual, which is metric for BDD reduced in value when anomalies occur (see the residual decrease in Figure 1 during anomalous operation shown by the dotted red line), implying that the standard SE BDD, along with any other residual-based method, will not be able to detect this coordinated attack.

In addition to standard ACSE BDD, we also implemented the auto-regression (VAR) method to detect grid anomalies. As the VAR algorithm does not consider dynamic graph properties of the power grid, it tends to create alarms on all abrupt measurement changes. This will easily lead to false positives since regular topological changes also result in sudden temporal change. This can be seen in the Figure 1 wherein during regular topology changes (shown in black dotted line) sudden spikes in VAR anomaly score are observed.

In comparison, our proposed method is able to detect all anomalies without False positives (FP). This indicates that the proposed algorithm is more likely to detect anomalies due to complex attack scenarios while being able to reduce the occurrences of false positives.

Here we explore the realistic grid pattern, using a utility-provided real-world dataset (for reasons of confidentiality, we cannot make this dataset public).

Dataset description- 24 hours worth of ACSE output data from a real utility in the Eastern Interconnect of the U.S. The dataset contains all operational data of the grid based on a 0.5-hour interval.

Here we document the statistical findings with the following notable observations:

• •

  Figure 10(a): Topology changes frequently in reality, thus our proposed method, which considers the impact of topology change, provides realistic values in handling time-series sensor data.

• •

  Figure 10(b): Real part of daily loads and generations change smoothly within a range of $(100\pm 8)\%$ times their average values, with the generation following the load demand throughout the day.

• •

  Figure 11: $90\%$ of the time, an individual real part of load changes less than $10\%$ of its maximum daily value within a 30-min interval; less than $17\%$ within a 2-hour interval; and less than $50\%$ within 24 hours.

In the worst possible case of measurement non-synchronization, the measurements and status data can be collected at such a different time that the measurements and reference topology are inconsistent with each other. If this happens, our method should detect it as an anomaly, since they are equivalent to topological or measurement errors.

However, the impact of this event is low. In case “uncorrected” unsynchronized measurements or topology exists, then that is an error that would disrupt the grid operator’s state estimation. In practice, the grid has some mechanisms to prevent this from happening:

1. 1.

   Both status data and analog measurements are sampled very frequently: status data are updated upon change of status in PJM and every 4s in ISO-NE; analog measurements are updated every 10s in PJM. Considering that the frequency of topology change is low, the high acquisition frequency is likely to guarantee all information is up to date.

2. 2.

   Topology estimation (TE) can check the consistency between switch status and analog measurements: e.g., if the switch is observed as ‘open’ but the current measurement on it is non-zero, then this status data is wrong and should be corrected. This processing is within the TE algorithms in the control room.

Therefore, it is likely that we can still obtain a good reference topology not affected by the non-synchronization.

Considering that: 1) the load changes slowly following a daily cycle; 2) conventional generators are re-dispatched every 15 min (less frequent than data acquisition we have sufficient time ticks before each re-dispatch) and each re-dispatch is limited by ramp rate; 3) renewable forecasts are improving and the auxiliary services like batteries are making them more stable and “firm” sources of energy, it’s still reasonable to assume that for any time $t$, its recent previous measurements can provide meaningful information.

What if sudden load change happens? At transmission grid level the load is aggregated at high-voltage nodes because of which abrupt change in load is uncommon. However, there exist rare events like load shedding or an unusual re-dispatch where there might be a sudden redistribution of power in the grid in a very short interval. In case of these rare events, the proposed method may detect the abrupt change as an anomaly. In reality, these events may not be anomalous and will be misclassified as anomalies. During such instances, we will rely on the operator to clear these false positives in their decision-making process where additional trustworthy information is available. This is reasonable since these events are rare and the operators are well aware of the situation of active ongoing load-shedding in a region. In case that the abrupt load change is unexpected to the operator, an indicator for anomalous behavior might be a rather useful one.