Distributed Estimation and Detection of Anomalies in Control Systems

Abstract

Methods and Systems for detecting anomalies in a control area of a control system. Estimating for the control area, a first state from a historical state over a first time period using a model of dynamics, and defining a transition of the first state as a function of control inputs, the first state includes a generator state for each generator, the control inputs include a network state for each bus, a mechanical input to each generator or power consumptions at the buses. Updating estimated first state, by connecting measurements of rotor frequency of each generator and measurements of the network states on the buses with the generator state of each generator, to obtain a second state over a second time period later than the first time period, and detecting anomalies based on a statistic deviation of the second state from its corresponding prediction derived from the first state.

Description

FIELD

The present disclosure relates generally to electric power systems, and more particularly to electrical anomaly detection in electric power systems.

BACKGROUND

Anomaly detection has been applied to control systems to increase cyber-physical security and optimize operations of a control system. The principles typically used in anomaly detection can include identifying normal behavior and a threshold selection procedure for identifying anomalous behavior. Usually, the challenge is to develop a method that is able to detect the abnormalities specific to an industry's needs.

An anomaly such as a cyber-physical attack may be undetectable from tampered measurements if there is a set of normal operating conditions consistent with the tampered measurements. Conventional cyber-physical attack detection approaches are centralized implemented.

In cyber-physical attack security applications one of the major problems is distinguishing between normal circumstance and “anomalous” or “abnormal” circumstances. For example, malfunction mechanical devices can be viewed as abnormal modifications to normal programs. The detection of anomalous activities is a difficult problem in which the detection of anomalous activities is disadvantaged by not having appropriate data and/or because of the variety of different activities that need to be monitored. Additionally, protective measures based on established conventional practices are vulnerable to activities designed specifically to undermine these assumptions.

Therefore, there is a need for developing more advanced distributed methods for real-time estimation and detection of anomalies in control systems.

SUMMARY

The present disclosure relates to methods and systems of dynamic state estimation of Electric Power Systems (EPS) for detecting anomalies in control areas of a control system having multiple control areas. The disclosure includes methods and systems relating to the detection of anomalies in a control area positioned within a multiple of neighbouring control areas in the control system that takes into account both process and measurement anomalies.

According to embodiments of the present disclosure, the methods and systems for dynamic state estimation are based on partitioning a control area from the multiple neighboring control areas in the control system, and detecting the anomalies for that individual control area.

Specifically, the present disclosure is based on the realization that dynamics of the control area can be decoupled if the voltage phase angles on the buses are treated as control inputs to a model of dynamics of the control area defining a transition of the state of the control area as a function of control inputs. At least one aspect regarding the decoupling dynamics of the control area from the control system, is that it allows for detecting anomalies within an individual control area rather than for the entire control system.

This is significant, because convention control system state variables include states of generators and buses, i.e., a rotor angle and a rotor frequency of generators and voltage phase angles on buses for the entire control system and not for a single control area. The states of the convention control system are globally coupled because the phase angle on a bus is related to the states of all buses within the control system. Cyber-attacks target entire computer information systems and infrastructures such as control systems, by manipulating with measurements of the control system as part of cyberterrorism. Resulting in the cyber-attacks being undetectable from tampered measurements of convention control systems, if there is a set of normal operating conditions consistent with those measurements when detection is based on the entire central system.

We realized that by avoiding centralized processing, and focusing on a single control area we make it more difficult for cyber-attackers to stage a global attack on the entire control system.

Our realization of decoupling dynamics of the control area from the entire control systems is contrary to the technical engineering principles of the structured network of the conventional control system, due to the fact that system networks of electrical components used to supply, transfer and use electric power are a dynamically coupled system. In other words, the conventional control system has a state on a bus of the generator that depends on states of all buses in the entire conventional control system, so that to detect a cyber-attack on one control area, the state of the entire conventional control system is required. The reason for this requirement is that the model of dynamics for the entire conventional control system's state transition is coupled together, and the measurement model of the conventional control system is also coupled together, which requires that any test for detecting a cyber-attack must be for the entire conventional control system and not for a single control area, as according the present disclosure. Another realization that the methods and systems of the present disclosure are based on is that the state of generators can be estimated according to a measurement model of the first control area, i.e. control area, connecting measurements of the rotor frequency of each generator and measurements of the voltage phase angles on the buses of the control area with the rotor angle and the rotor frequency of each generator of the control area. For example, we are able to estimate a first state of the control area from a historical state of the control area over a first state time period using a model of dynamics of the control area. In order to determine a transition of the first state of control area we used a function of control inputs, where the state of the control area includes a rotor angle and a rotor frequency for each generator in the control area. Based on the control inputs including one or combination of phase angles on the buses of the control area, phase angles on some neighbouring buses of neighbouring control areas, a mechanical input to each generator in the control area, or power consumptions at the buses in the control area. Essentially, embodiments of the present disclosure are based on tracking the dynamic relationship between local states and local operational measurements to detect the anomalies or cyber-attack in a distributed fashion with limited neighbouring communications. Wherein, no new measurements or devices are needed or required to be added to the EPS system, which provides an unencumbered implementation for the methods and systems of the present disclosure into the control systems industry.

We further realized that a cyber-attack can be detected in a distributed fashion based on a statistic deviation of the state of the local area, control area, from a normal state of the local area even considering uncertainties of the model of dynamics and the model of measurements. Wherein we use a Kalman filter based approach to deal with the level of ambiguity under which we assess how likely a cyber-attack has occurred introduced by the statistic deviation or noises. The Kalman filter can be distributed applied for each control area through iterative liner solvers and a neighborhood-approximation to an estimated state covariance matrix. In comparison, the estimated state covariance used in a centralized Kalman filter is, in general full, due to the global coupling of the bus voltages with the generator rotor angles for an interconnected power system.

We also formulate the dynamic state estimation problem using a structure-preserving model that preserves the sparse coupling of the dynamics. Both the generator state, i.e., rotor angle and frequency of each generator, and the network state, i.e., voltage phase angle at each bus, are considered in the estimation. Wherein we decompose the problem in such a way that each control area solves a problem of reduced dimension, and a decoupled state-space formulation is also used to avoid communication requirements between all generators.

For example, we updated the first state of the control area according to a measurement model of the control area, by connecting measurements in the control area of the rotor frequency of each generator and measurements of the phase angles on the buses, with the rotor angle and the rotor frequency of each generator in the control area, to obtain a second state of the control area over a second state time period that is later than the first state time period. Wherein we were able to detect the cyber-attack based on a statistic deviation of the second state of the control area from its corresponding prediction derived from the first state of the control area.

The local attack detection statistic is designed jointly with the dynamic state estimation in order to limit communication requirements. The communication requirements consist of buses' sharing their weighted measurement residuals within a user-specified neighborhood. The size of the neighborhood can be adjusted to allow for a tradeoff between accuracy and communication overhead.

According to another embodiment of the disclosure, a method for detecting anomalies in a control area of a control system with multiple control areas. The control area includes a set of generators in communication with a set of buses and some neighboring buses of neighboring controls area of the control system. The method including accessing a memory in communication with at least one processor, wherein the memory includes stored historical states of the control area. Estimating, by the processor, a first state of the control area over a first state time period from the historical states of the control area. Using a model of dynamics of the control area, and defining a transition of the first state of the control area as a function of control inputs. Wherein the first state of the control area includes a generator state for each generator in the control area. Further, wherein the control inputs include one or combination of: a network state for each bus in the control area; a network state for some neighboring buses of neighboring controls areas; a mechanical input to each generator in the control area or power consumptions at the buses in the control area. Updating, by the processor, the estimated first state of the control area according to a measurement model of the control area, by: connecting measurements of the rotor frequency of each generator; a weighted combination of measurements of the network states on the buses in the control area and some neighboring buses of neighboring controls areas; with the generator state of each generator in the control area, to obtain a second state of the control area over a second state time period later than the first state time period. Detecting, by the processor, the anomalies based on a statistic deviation of the second state of the control area from its corresponding prediction derived from the first state of the control area.

According to another embodiment of the disclosure, a system for detecting anomalies in a control area of a power system with multiple control areas. The control area includes a set of generators in communication with a set of buses and some neighboring buses of neighboring controls area of the control system. The system including a computer readable memory to store historical states of the control area, current states of the control area, a model of dynamics of the control area and a measurement model of the control area. A set of sensors arranged to sense measurements in the control area for the set of generators, the set of buses and to measure one or combination of rotor frequencies for each generator, voltage phase angles on the buses, a mechanical input to each generator, or power consumptions at the buses. A processor in communication with the computer readable memory configured is to: estimate a first state of the control area from a historical state of the control area over a first state time period using the stored model of dynamics of the control area, and defining a transition of the first state of the control area as a function of control inputs. Wherein the first state of the control area includes a rotor angle and a rotor frequency for each generator in the control area. Further, wherein the control inputs include one or combination of phase angles on the buses of the control area, a mechanical input to each generator in the control area, and power consumptions at the buses in the control area. Update the estimated first state of the control area according to the stored measurement model of the control area, by connecting measurements of the rotor frequency of each generator in the control area and a weighted combination of measurements of the voltage phase angles on the buses of the control area and some neighboring buses of the neighboring control areas, with the rotor angle and the rotor frequency of each generator in the control area, to obtain a second state of the control area over a second state time period later than the first state time period. Detect the anomalies based on a statistic deviation of the second state from its corresponding prediction derived from the first state of the control area.

According to another embodiment of the disclosure, a detector for detecting anomalies in a control area of an electric power system (EPS) with multiple control areas. The control area includes a set of generators in communication with a set of buses and some neighboring buses of neighboring controls area of the control system. The detector includes acquiring a plurality of measurements from sensors configure for sensing the set of generators and the set of buses over a first state time period. Acquiring a respective second plurality of measurements from sensors configure for sensing the set of generators and the set of buses over a second state time period. At least one processor having a computer readable memory configured to receive the plurality of measurements sensed by the sensors over the first state time period. Estimate a first state of the control area from a historical state of the control area over the first state time period using a model of dynamics of the control area, and defining a transition of the first state of the control area as a function of control inputs. Wherein the first state of the control area includes a rotor angle and a rotor frequency for each generator in the control area, wherein the control inputs include one or combination of phase angles on the buses of the control area, a mechanical input to each generator in the control area, and power consumptions at the buses in the control area. Receive the plurality of measurements sensed by the sensors over the second state time period. Update the estimated first state of the control area according to a measurement model of the control area, by connecting sensed measurements of the rotor frequency of each generator in the control area and a weighted combination of sensed measurements of the phase angles on the buses of the control area and some neighboring buses of the neighboring control areas, with the rotor angle and the rotor frequency of each generator in the control area, to obtain a second state of the control area over the second state time period later than the first state time period. Detect the anomalies based on a statistic deviation of the second state of the control area from its corresponding prediction derived from the first state of the control area.

Further features and advantages of the present disclosure will become more readily apparent from the following detailed description when taken in conjunction with the accompanying Drawing.

BRIEF DESCRIPTION OF THE DRAWINGS

The present disclosure is further described in the detailed description which follows, in reference to the noted plurality of drawings by way of non-limiting examples of exemplary embodiments of the present disclosure, in which like reference numerals represent similar parts throughout the several views of the drawings. The drawings shown are not necessarily to scale, with emphasis instead generally being placed upon illustrating the principles of the presently disclosed embodiments.

FIG. 1A is a schematic of a power network in which anomalies and/or cyber-physical attacks are detected, according to one embodiment of the present disclosure;

FIG. 1B is a schematic for partitioning the power network into multiple control areas that anomalies and/or cyber-physical attacks are detected independently for each control area with limited communications with neighbors, according to one embodiment of the present disclosure;

FIG. 1C is a schematic for a measurement configuration of the power network, according to some embodiments of the present disclosure;

FIG. 2A is a schematic block diagram of an example of a control system that may be used with one or more embodiments described herein, according to some embodiments of the present disclosure;

FIG. 2B is a block diagram of a method for distributed estimation and detection of anomalies in an electric power system (EPS), according to one embodiment of the present disclosure;

FIG. 2C is a schematic of an exemplar EPS controlled system, according to some embodiments of the present disclosure;

FIG. 3A is a schematic block diagram for distributed estimation and detection of cyber-physical attacks for predetermined period of covariance matrix update, according to some embodiments of the present disclosure;

FIG. 3B is a schematic block diagram for determining local area attack statistic and declare possible attack, according to some embodiments of the present disclosure;

FIG. 4A is a schematic histogram of the local attack statistic for Area I, d_I before attack at area I, according to some embodiments of the present disclosure;

FIG. 4B is a schematic histogram of the local attack statistic for Area I, d_I during attack at area I, according to some embodiments of the present disclosure;

FIG. 4C is a schematic histogram of the local attack statistic for Area I, d_I after attack at area I, according to some embodiments of the present disclosure;

FIG. 4D is a schematic histogram of the local attack statistic for Area IV, d_IV before attack at area I, according to some embodiments of the present disclosure;

FIG. 4E is a schematic histogram of the local attack statistic for Area IV, d_IV during attack at area I, according to some embodiments of the present disclosure; and

FIG. 4F is a schematic histogram of the local attack statistic for Area IV, d_IV after attack at area I, according to some embodiments of the present disclosure.

While the above-identified drawings set forth presently disclosed embodiments, other embodiments are also contemplated, as noted in the discussion. This disclosure presents illustrative embodiments by way of representation and not limitation. Numerous other modifications and embodiments can be devised by those skilled in the art which fall within the scope and spirit of the principles of the presently disclosed embodiments.

DETAILED DESCRIPTION

The following description provides exemplary embodiments only, and is not intended to limit the scope, applicability, or configuration of the disclosure. Rather, the following description of the exemplary embodiments will provide those skilled in the art with an enabling description for implementing one or more exemplary embodiments. Contemplated are various changes that may be made in the function and arrangement of elements without departing from the spirit and scope of the subject matter disclosed as set forth in the appended claims.

Specific details are given in the following description to provide a thorough understanding of the embodiments. However, understood by one of ordinary skill in the art can be that the embodiments may be practiced without these specific details. For example, systems, processes, and other elements in the subject matter disclosed may be shown as components in block diagram form in order not to obscure the embodiments in unnecessary detail. In other instances, well-known processes, structures, and techniques may be shown without unnecessary detail in order to avoid obscuring the embodiments. Further, like reference numbers and designations in the various drawings indicated like elements.

Also, individual embodiments may be described as a process which is depicted as a flowchart, a flow diagram, a data flow diagram, a structure diagram, or a block diagram. Although a flowchart may describe the operations as a sequential process, many of the operations can be performed in parallel or concurrently. In addition, the order of the operations may be re-arranged. A process may be terminated when its operations are completed, but may have additional steps not discussed or included in a figure. Furthermore, not all operations in any particularly described process may occur in all embodiments. A process may correspond to a method, a function, a procedure, a subroutine, a subprogram, etc. When a process corresponds to a function, the function's termination can correspond to a return of the function to the calling function or the main function.

Furthermore, embodiments of the subject matter disclosed may be implemented, at least in part, either manually or automatically. Manual or automatic implementations may be executed, or at least assisted, through the use of machines, hardware, software, firmware, middleware, microcode, hardware description languages, or any combination thereof. When implemented in software, firmware, middleware or microcode, the program code or code segments to perform the necessary tasks may be stored in a machine readable medium. A processor(s) may perform the necessary tasks.

Definition of Terms

We use v_k to denote the k-th entry of a vector v, and the (i,j)-th entry of a matrix M is given by M_ij. The i-th row of a matrix M is denoted M_i. The transpose of a vector or matrix X is denoted by X^T. The inverse of a matrix M is denoted by M⁻¹. The entry-wise matrix multiplication is denoted ⊗. A matrix M is positive definite is denoted M0, and AB means A−B0.

The present disclosure is directed to dynamic state estimation in electric power systems (EPS), for detecting anomalies in control areas of a control system having multiple control areas. In particular, the detection of anomalies in a single control area positioned within a multiple of neighbouring control areas in the control system that takes into account both process and measurement anomalies.

Specifically, the present disclosure is based on realizing that dynamics of the control area can be decoupled if the phase angles on the buses are treated as control inputs to a model of dynamics of the control area defining a transition of the state of the control area as a function of control inputs. At least one aspect regarding the decoupling dynamics of the control area from the control system, is that it allows for detecting anomalies within an individual control area rather than for the entire control system. Thus, by avoiding centralized processing, and focusing on a single control area we make it more difficult for cyber-attackers to stage a global attack on the entire control system.

The methods and systems of the present disclosure are also based on that the state of generators can be estimated according to a measurement model of the first control area, i.e. control area, connecting measurements of the rotor frequency of each generator and measurements of the voltage phase angles on the buses of the control area with the rotor angle and the rotor frequency of each generator of the control area. For example, we estimate a first state of the control area from a historical state of the control area over a first state time period using a model of dynamics of the control area.

In order to determine a transition of the first state of control area we use a function of control inputs, where the state of the control area includes a rotor angle and a rotor frequency for each generator in the control area. Based on the control inputs including one or combination of phase angles on the buses of the control area and neighbouring border buses, a mechanical input to each generator in the control area or power consumptions at the buses in the control area. In other words, we are essentially tracking the dynamic relationship between local states and local operational measurements to detect the anomalies or cyber-attack in a distributed fashion.

We then update the estimated first state of the control area according to a measurement model of the control area, by connecting measurements of the rotor frequency of each generator in the control area and measurements of the weighted phase angles on the buses of the control area, with the rotor angle and the rotor frequency of each generator in the control area, to obtain a second state of the control area over a second state time period later than the first state time period. Lastly, we detect the anomalies based on a statistic deviation of the second state of the control area from its corresponding prediction derived from the first state of the control area.

Formulation for State Space Model and Attack Detection of a Power System

FIG. 1 is a schematic of a power network in which anomalies and/or cyber-physical attacks are detected, according to one embodiment of the present disclosure. The system includes a set of generators 110, 140 and connected to a set of loads 130, 139 through transmission lines 150, 160. A generator 110, or load 130 is connected to the system through a bus 120, and a line 150 is connected to the system through two terminal buses.

The neighbourhood distance is used in the present disclosure to determine the communication requirements between generators or portions of power system. The distance of neighbourhood, that is the distance between a pair of neighboring generators or loads in term of number of hops, is measured by the number of connected branches residing in the path between the neighboring buses that the generators or loads connect to. For example, generator 110 is connected to generator 140 through 3 transmission lines, 150, 160 and 170, therefore generator 110 and generator 140 is 3-hop neighbors.

The power system is partitioned to be a set of control areas in the present disclosure. Two adjacent control areas may have common transmission lines, but there are no overlapping buses between two control areas.

FIG. 1B is a schematic for partitioning the exemplar power network into multiple control areas that cyber-physical attack detected independently for each area with limited communications with neighbors. In FIG. 1B, the original power system as shown in FIG. 1A is partitioned 4 control areas, including area I, 155, area II, 165, area III, 175 and area IV, 180.

The anomalies and/or cyber-physical attacks are detected based on dynamic state estimation of the power system. The state estimation may use a set of measurements or input information acquired by a set of sensors from the system.

FIG. 1C is a schematic for measurement configuration of the power system, according to embodiments of the present disclosure. The measurements may include the rotor frequencies of generators, 159, and voltage phase angles of buses 169. The generator mechanical inputs of generators, 179 and load power consumption at buses, 189 can also be available for estimating system states and detecting anomalies.

FIG. 2A is a schematic block diagram of an example of a control system 200 that may be used with one or more embodiments described herein. The control system 200 may comprise of at least one processor 201, and a memory 215 connected to the processor 201, as well as a set of sensors 230. Further, the control system 200 includes communication from neighbouring control areas 225 regarding neighbouring areas' states and attributes. Contemplated may include one or more network interfaces (e.g., wired, wireless, PLC, etc.) used with the system. Also contemplated may bet at least one power supply (i.e. battery, plug-in, etc.) may be used with the control system.

Still referring to FIG. 2A, the memory 215 can be implemented within the processor 201 and/or external to the processor 201. In some embodiments, the memory 215 stores historical states of generators and buses of EPS including a first control area, received neighbouring control area states and attributes and determined control inputs of the first control area.

Similarly, the memory 215 can store the state dynamic model 220, measurement model 240, anomaly detection model 250, and/or instructions to the processor 201 of how to use the state dynamics model 220, measurement model 240 and anomaly detection model 250. For example, in various embodiments, the processor 201 determines estimations of the updated generator rotor angles and frequencies corresponding to the updated current state using the measurement model that relates the measurements to state variables.

Still referring to FIG. 2A, in some embodiments, the processor 201 can include functionality blocks including a state dynamics model block 220, a measurement model block 240, and an anomaly detection model block 250. The data receiving, estimation execution, and results transmitting can be implemented in the processor 201.

Further, the processor 201 may comprise hardware elements or hardware logic adapted to execute software programs and manipulate data structures. The memory 215 may include an operating system, portions of which can typically resident in memory 215 and executed by the processor 201, functionally organizes the control system by, inter alia, invoking operations in support of software processes and/or services executing on the device.

As may be apparent to those skilled in the art that other processor and memory types, including various computer-readable media, may be used to store and execute program instructions pertaining to the techniques described herein. Also, while the description illustrates various processes, contemplated may be that various processes can be embodied as modules configured to operate in accordance with the techniques herein (e.g., according to the functionality of a similar process). Further, while the processes have been shown separately, those skilled in the art will appreciate that processes may be routines or modules within other processes.

Still referring to FIG. 2A, the set of sensors 230 can include phasor measurement units which provide measurements for voltage phase angle for a bus. The sensors 230 can also include remote measurement units which provide the rotor frequency for a bus.

The system 200 can also include a screen or display 270. In some embodiments, the result of the anomaly detection can be rendered on the display 270 or submitted to different applications that can be internal or external to the system. For example, the results can be sent to a real-time security control and prevention application which can be added to the processor 201.

FIG. 2B is a block diagram of a method 202 for distributed estimation and detection of anomalies in an electric power system (EPS), according to one embodiment of the present disclosure. The method including partitioning the EPS into multiple neighbouring control areas 207, wherein the multiple neighbouring control areas include a first control area. Further, estimating current states of generators in the first control area using the states of the generators via a state dynamics model 221 of the generators of the first control area, wherein the states of the generators is from historical data of the EPS 209. Further, receiving, over a communication channel, a current state and attribute of at least some of neighbouring buses located in the neighboring control areas 226. Wherein the neighbouring buses are the buses of neighbouring control areas that connects to at least one boarder bus of the first control area via at least one transmission line within pre-determined neighbourhood distance. Determining measurements of states from generators and buses in the first control area 233. Then, updating the current state of the generators in the first control area using a measurement model 236 on a basis of the estimated current states of the generators in the first control area, the received states and attributes from the neighbouring control areas, and the measurements of states from at least some generators and buses in the first control area. Finally, determining anomaly based on statistic of deviation of measurement from the states of first control area using anomaly detection model 261, so as to update the current state and anomaly status 272. Wherein at least some steps of the method are performed by one or more processor and the updated current state and detection results may be displayed on a display, i.e. monitor or some other related electronic display device.

FIG. 2C shows a schematic of an exemplar EPS controlled according to one embodiment of the present disclosure. The power generation facilities 110 are coupled to substations 120. Associated with substations 120 is a regional control module 276. The regional control module manages power production, distribution, and consumption within its region. Also associated with each region are industrial/commercial/residential loads 130 representative of industrial plant or large commercial enterprises or residential loads. According to some embodiments of the disclosure, each regional control module 276 using one or more applications is operable to manage the power production and distribution within its region, and monitoring the healthy and security status of the region of grid under its control.

In some implementations, power producing entity, such as the power generation plants 110, interfaces with the regional grid via a local control module 272. The local control module 272 can standardize control command responses with each of the plurality of power providers. By offering to the regional control module 276 a standardized response from each of the plurality of power producing entities, the regional control module can actively manage the power grid in a scalable manner. The regional control module 276 is further aware of the electricity producing capacity within the region and the limitations to the power grid. The regional control module 276 understands topology with respect to the power providers and power consumers and its own ability to distribute the power.

Each regional control module 276 is communicatively coupled to a control system 277 via, e.g., a wide area network 274. The wide area network can be the Internet or other means to communicate data among remote locations. Additionally, or alternatively, the data can be exchanged between the control system 277 and the regional control modules 276 via a local area network or Intranet. To that end, the control system 277 includes a transceiver 280 for exchanging data between the control system and regional control modules 276 via the network 274. Also, control system 277 includes one or several processors 211A and 211B to balance amounts of electricity passing through an electrical grid.

The control system 277 is operable to manage the interaction of several regional control modules 276 and the power providers under their control. As previously described, each regional control module 276 using applicable applications can dynamically manage the power consumers and power providers within its control. As demand within a certain region managed by a regional control module 276 increases or decreases, the regional control module 276 needs to act to compensate for power production within a particular region. To that end, the regional control module 276 makes a decision about supplying or requesting the electricity from the grid. The control system 277 receives, transmits or retransmits such request to balance amount of electricity going in or off the grid.

Different embodiments of the disclosure monitor one or combination of the voltage phase angles on the buses of the EPS and/or the rotor frequencies of generators in the EPS. For example, one embodiment determines anomalies from the updated states of the EPS, and then the control system 277 can issue a command to the regional control module 276 to control their production or loads or network topology accordingly. For example, the regional control module 276 can issue a command to a local control module 272 of a generation facility 110 to block its automatic generation control to avoid executing incorrect load flowing function, when an anomaly 15 found.

Structure-Preserving Model

Consider a power system that includes a set of {1, 2, . . . , n} buses and a set of transmission lines connecting between the buses. Let the subset of buses with generators be denoted G with |G|=n_g. Let the subset of buses with loads be denoted D with |D|=n_d, and n_g+n_d=n.

Assumed that the transient voltage behind its internal reactance of a generator at bus i, i∈G is constant, then the dynamics of the generator can be described as:

$\begin{array}{cc}{\stackrel{.}{\delta }}_i\left(t\right)={\omega }_i\left(t\right),M_i{\stackrel{.}{\omega }}_i\left(t\right)=P_i^G\left(t\right)-\frac{E_iV_i}{Z_i}\sin \left({\delta }_i\left(t\right)-{\theta }_i\left(t\right)\right)-D_i{\omega }_i\left(t\right),& \left(1\right)\end{array}$

where E_i and δ_i are the internal voltage magnitude and the rotor angle of generator at bus i, M_i, D_i and Z_i are the inertia, the damping coefficient and the transient reactance of generator at bus i, V_i and θ_i are the voltage magnitude and the phase angle of bus i, P_i^G is the equivalent power injections of mechanical power input of generator at bus i.

The power flow equation at generator bus i, i∈G is as follows:

$\begin{array}{cc}{P}_i^D=\frac{E_iV_i}{Z_i}\sin \left({\theta }_i-{\delta }_i\right)+\sum _{j=1}^nV_iV_jY_{\mathrm{ij}}^B\sin \left({\theta }_i-{\theta }_j\right)+\sum _{j=1}^nV_iV_jY_{\mathrm{ij}}^G\cos \left({\theta }_i-{\theta }_j\right),Q_i^D=\frac{1}{Z_i}V_i^2-\frac{E_iV_i}{Z_i}\cos \left({\theta }_i-{\delta }_i\right)+\sum _{j=1}^nV_iV_jY_{\mathrm{ij}}^G\sin \left({\theta }_i-{\theta }_j\right)-\sum _{j=1}^nV_iV_iY_{\mathrm{ij}}^B\cos \left({\theta }_i-{\theta }_j\right),& \left(2\right)\end{array}$

where, P_i^D and Q_i^D are the equivalent active and reactive power injections for power demands of generator bus i, Y_ij^G and Y_ij^B are the elements of bus conductance matrix and bus susceptance matrix at the row corresponding to bus i and the column corresponding to bus j.

Analogously, for any load bus i, i∈D (i.e. i∉G), the power flow equation at the bus is

$\begin{array}{cc}{P}_i^D=\sum _{j=1}^nV_iV_jY_{\mathrm{ij}}^B\sin \left({\theta }_i-{\theta }_j\right)+\sum _{j=1}^nV_iV_jY_{\mathrm{ij}}^G\cos \left({\theta }_i-{\theta }_j\right),Q_i^D=\sum _{j=1}^nV_iV_jY_{\mathrm{ij}}^G\sin \left({\theta }_i-{\theta }_j\right)-\sum _{j=1}^nV_iV_jY_{\mathrm{ij}}^B\cos \left({\theta }_i-{\theta }_j\right),& \left(3\right)\end{array}$

It can be assumed for a practical power system that all angular differences for each generator and line are small, the resistance at each transmission lines is ignorable, and the voltage at each generator or generator is close to its nominal rated value, and therefore we can have: sin(δ_i−θ_i)≈δ_i−θ_i, cos(δ_i−θ_i)≈1, sin(θ_i−θ_j)≈θ_i−θ_i, cos(θ_i−θ_j)≈1, Y_ij^G=0 and E_i=V_i=1 per unit. With these assumptions, the equations (1)-(3) can be linearized around the network steady state condition to yield following the dynamic linearized swing equation and the algebraic DC power flow equation:

$\begin{array}{cc}{\stackrel{.}{\delta }}_i\left(t\right)={\omega }_i\left(t\right),i\in GM_i{\stackrel{.}{\omega }}_i\left(t\right)=-\frac{1}{Z_i}{\delta }_i\left(t\right)-D_i{\omega }_i\left(t\right)+\frac{1}{Z_i}{\theta }_i\left(t\right)+P_i^G\left(t\right),i\in G0=\frac{1}{Z_i}{\delta }_i\left(t\right)-\frac{1}{Z_i}{\theta }_i\left(t\right)-\sum _{j=1}^nY_{\mathrm{ij}}^B\left({\theta }_i\left(t\right)-{\theta }_j\left(t\right)\right)+P_i^D\left(t\right),i\in G0=-\sum _{j=1}^nY_{\mathrm{ij}}^B\left({\theta }_i\left(t\right)-{\theta }_j\left(t\right)\right)+P_i^D\left(t\right),i\notin G,& \left(4\right)\end{array}$

Equation (4) is usually called as the linearized structure-preserving model. It can be also expressed in compact form as:

$\begin{array}{cc}\left[\begin{array}{ccc}I& 0& 0\\ 0& M^{\mathrm{INER}}& 0\\ 0& 0& 0\end{array}\right]\stackrel{.}{X}\left(t\right)=-\left[\begin{array}{ccc}0& -I& 0\\ L^{\mathrm{GG}}& D^{\mathrm{DAMP}}& L^{\mathrm{GD}}\\ L^{\mathrm{DG}}& 0& L^{\mathrm{DD}}\end{array}\right]X\left(t\right)+{\left[0^TP^G\left(t\right)P^D\left(t\right)\right]}^T,& \left(5\right)\end{array}$

where M^INER and D^DAMP are diagonal matrices, M_ii^INER=M_i, D_ii^DAMP=D_i; P^G (t)=└P₁^G (t), . . . , P_ng^G(t)┘, and P^D(t)=[P₁^D(t), . . . , P_n^D(t)]. The state X(t)=[δ(t)^T ω(t)^T θ(t)^T] consists of the rotor angle δ and rotor frequency ω at every generator and the voltage phase angle θ at every bus. The matrix

$L=\left[\begin{array}{cc}{L}^{\mathrm{GG}}& L^{\mathrm{GD}}\\ L^{\mathrm{DG}}& L^{\mathrm{DD}}\end{array}\right]$

is a Laplacian matrix with a sparsity structure related to the underlying network. We note that L^GG is a diagonal matrix,

$L_{\mathrm{ii}}^{\mathrm{GG}}=\frac{1}{Z_i},L^{\mathrm{DD}}$

is related to the bus susceptance matrix,

$L_{\mathrm{ii}}^{\mathrm{DD}}=\{\begin{array}{cc}\sum _{j=1,j\ne i}^nY_{\mathrm{ij}}^B+\frac{1}{Z_i},& i\in G\\ \sum _{j=1,j\ne i}^nY_{\mathrm{ij}}^B,& i\notin G\end{array},L_{\mathrm{ij}}^{\mathrm{DD}}=-Y_{\mathrm{ij}}^B\mathrm{and}L^{\mathrm{GD}}={\left(L^{\mathrm{DG}}\right)}^T,L_{\mathrm{ij}}^{\mathrm{GD}}=\{\begin{array}{cc}-\frac{1}{Z_i}& j=i,i\in G\\ 0& \mathrm{otherwise}\end{array}.$

The control input U(t)=[0^T P^G(t) P^D(t)]^T is given by the mechanical power input at the generators, P^G, and the electrical power demand at each bus, P^D, including those connected to a generator.

For a practical application, there exist both process and measurement noise. Additionally, the measurements are taken to be a discrete-time process rather than a continuous time process, which is more realistic for digitally sampled systems. The process noise v(t) and measurement noise n[k] are assumed to be Gaussian noise processes. With the introduction of noise, we cannot apply a deterministic filter for attack detection. For the noisy setting, we appeal to Kalman filtering techniques. However, since the dynamic system matrix in (5) is a singular matrix, the dynamic model in (5) is not directly applicable for Kalman filtering. One solution is to eliminate θ from the state dynamics through the equation

θ(t)=(L^DD)⁻¹(P^D(t)−L^DGδ(t)).  (6)

Although L^DD is sparse, its inverse is not, which leads to a coupling of the dynamics for ω(t) amongst all generators that is not present in the original formulation (5). Such a global coupling makes developing a distributed solution with reasonable communication requirements infeasible.

Dynamic Model

To solve the problem described above, we treat the bus voltage phase angle as a control input rather than eliminating it from the dynamics. The state given by the generator variables can be described as:

x(t)=└x₁(t)x₂(t) . . . x_ng(t)┘^T,  (7)

x _i(t)=[δ_i(t)ω_i(t)]^T.  (8)

The dynamics for the generator can be described as follows:

$\begin{array}{cc}{\stackrel{.}{\delta }}_i\left(t\right)={\omega }_i\left(t\right)+v_{{\delta }_i}\left(t\right),& \left(9\right)\\ {\stackrel{.}{\omega }}_i\left(t\right)=\frac{-1}{M_iZ_i}{\delta }_i\left(t\right)-\frac{D_i}{M_i}{\omega }_i\left(t\right)+\frac{1}{M_i}\left(P_i^G\left(t\right)+\frac{1}{Z_i}{\theta }_i\left(t\right)\right)+v_{{\omega }_i}\left(t\right),& \left(10\right)\end{array}$

where v_δi(t) and v_ωi(t) are the process noises on δ_i and ω_i, respectively. Therefore, these equations can be collected in matrix-form:

{dot over (x)}(t)=Āx(t)+u(t)+v(t),  (11)

where Ā is block diagonal, and its non-zero entries are set as:

$\stackrel{\_}{A}\left({\stackrel{.}{\delta }}_i\left(t\right),{\omega }_i\left(t\right)\right)=1.0,\stackrel{\_}{A}\left({\stackrel{.}{\omega }}_i\left(t\right),{\delta }_i\left(t\right)\right)=\frac{-1}{M_iZ_i},\stackrel{\_}{A}\left({\stackrel{.}{\omega }}_i\left(t\right),{\omega }_i\left(t\right)\right)=\frac{D_i}{M_i}.$

The control input is

$\begin{array}{cc}u\left(t\right)=\lfloor u_{{\delta }_1}\left(t\right)u_{{\omega }_1}\left(t\right)\dots u_{{\delta }_{n_g}}\left(t\right)u_{{\omega }_{n_g}}\left(t\right)\rfloor .& \left(12\right)\\ u_{{\delta }_i}\left(t\right)=0,u_{{\omega }_i}\left(t\right)=\frac{1}{M_i}\left(P_i^G\left(t\right)+\frac{1}{Z_i}{\theta }_i\left(t\right)\right).& \left(13\right)\end{array}$

The process noise is assumed to be v(t)˜N(0, Q) and uncorrelated in time. After digital sampling with sampling period T, the continuous-time dynamic system is converted to a discrete-time system:

$\begin{array}{cc}x\left[k+1\right]=\mathrm{Ax}\left[k\right]+\mathrm{Bu}\left[k\right]+v\left[k\right],\mathrm{where}& \left(14\right)\\ A=e^{\stackrel{\_}{\mathrm{TA}}},B={\int }_0^{\stackrel{\_}{T}}e^{\tau \stackrel{\_}{A}}d\tau ,& \left(15\right)\end{array}$

and v[k]˜N(0, Q). The sampled process noise covariance matrix Q is related to the un-sampled covariance matrix Q via

$\begin{array}{cc}Q={\int }_0^{\stackrel{\_}{T}}e^{\stackrel{\_}{A}\tau }\stackrel{\_}{Q}e^{{\stackrel{\_}{A}}^T\tau }d\tau .& \left(16\right)\end{array}$

To evaluate the integrals in (15)-(16), we use a second-order approximation:

$A\approx I+\stackrel{\_}{\mathrm{TA}}+\frac{{\stackrel{\_}{T}}^2}{2!}{\stackrel{\_}{A}}^2,B\approx \stackrel{\_}{T}I+\frac{{\stackrel{\_}{T}}^2}{2!}\stackrel{\_}{A},Q\approx \stackrel{\_}{\mathrm{TQ}}+\frac{{\stackrel{\_}{T}}^2}{2!}\left(\stackrel{\_}{\mathrm{AQ}}+{\stackrel{\_}{\mathrm{QA}}}^T\right).$

The voltage phase angles are usually measured by a Phasor Measurement Unit (PMU). Since the sampling rate of the PMUs is high enough to track the system dynamics, we assume the second-order approximation is accurate.

Measurement Model

The measurements used in this present disclosure include the measurements of the generator rotor frequencies ω, and the measurements of the bus voltage phase angles θ. Similar to the dynamics, using θ directly in the measurement equations yields the following measurement model:

{circumflex over (ω)}_i[k]=ω_i[k]+n_ωi[k],  (17)

{circumflex over (θ)}_i^k[k]=[−(L^DD)⁻¹L^DG]_iδ[k](L^DD)_i⁻¹P^D[k]+n_θi[k],  (18)

where, P^D[k] is a known input. This formulation is not desirable for distributed processing since (L^DD)⁻¹ couples the measurement of the phase angle {circumflex over (θ)} at any given bus to the rotor angle δ at all generators and the electrical power demand P^D at all buses. We handle this by instead considering measurement of

{circumflex over (θ)}_i*[k]≙L_i^DD{circumflex over (θ)}[k]−P_i^D[k]  (19)

The relation to the state is:

$\begin{array}{cc}{\theta }_i^{*}\left[k\right]=\frac{1}{Z_i}{\delta }_i\left[k\right]+n_{{\theta }_i^{*}}\left[k\right]& \left(20\right)\end{array}$

where n_θi*[k] is the measurement noise.

Let y[k]=[y₁[k] . . . y_n[k]]^T, where y_i[k] is the measurement set local to bus i, and defined as:

$\begin{array}{cc}{y}_i\left[k\right]=\{\begin{array}{c}\left[{\hat{\omega }}_i\left[k\right]{\hat{\theta }}_i^{*}\left[k\right]\right]i\in G\\ {\hat{\theta }}_i^{*}\left[k\right]i\notin G\end{array},& \left(21\right)\end{array}$

We can have the following decoupled measurement model:

y[k]=Hx[k]+n[k].

The measurement matrix, H can be decoupled (i.e. block diagonal), since the measurements y_i involve only local variables ω_i and δ_i as shown in equations (17) and (20). The non-zero elements for the measurement matrix H are:

$H\left({\hat{\omega }}_i,{\omega }_i\right)=1.0,i\in G\mathrm{and}H\left({\hat{\theta }}_i^{*},{\delta }_i\right)=\frac{1}{Z_i},i\in G.$

This new formulation maintains the original sparse, localized coupling inherent to power systems rather than a global coupling. The quantity θ_i* in (19) is a linear combination of voltage phase angles at neighboring buses and thus can be calculated in a distributed fashion with limited communication, since L_i^DD only have non-zero element for the entry corresponding to the buses that have direct connection with bus i. Furthermore, only the local electric power demand P_i^D is needed at each bus rather than the global P^D.

However, measurement of θ* introduces the following complication. The measurement noise is Gaussian with n_ω˜N(0,R_ω) and n_θ˜N(0, R_θ), where R_ω and R_θ are assumed to be diagonal. However, R_θ*=L^DDR_θ(L^DD)^T is not diagonal. R_θ* introduces coupling and communication requirements that depend on neighboring rather than global information. The sparsity pattern of L^DD is the same as the adjacency matrix of the underlying network. The matrix L_ij^DD is nonzero if and only if bus i and bus j are neighbors, and [L^DDR_θ(L^DD)^T]_ij is nonzero if and only bus i and bus j are at most 2-hop neighbors. In summary, we have a new formulation that transfers all of the coupling to the measurement covariance matrices and process-measurement covariance.

Kalman Filtering with Correlated Process and Measurement Noise

Due to the presence of θ as a control input in (13) and in the measurements in (19), the process noise v[k] is correlated with the measurement noise n[k+1]. Let

E{v[k]n[j] ^T } Cδ _j,k+1,  (22)

where δ is the Kronecker delta, that is δ is 1 if j equals to (k+1), and 0 otherwise. Before discretization the control input to δ_i is zero, and the control input to ω_i depends only on θ_i. After discretization, the matrix B in (15) is block-diagonal introducing a coupling between u_δi and θ_i. Therefore, we specify C to have non-zero entries only at:

C(u_δi,{circumflex over (θ)}_j*) if j=i or j∈N_i

C(u_ωi,{circumflex over (θ)}_j*) if j=i or j∈N_i.

N_i is the set of 1-hop neighboring buses of bus i.

The Kalman gain and state covariance estimation update formulas with process and measurement noise correlated according to c are given below:

{circumflex over (x)} ⁻ [k]=A{circumflex over (x)} ⁺ [k−1]+Bu[k−1],  (23)

P ⁻ [k]=AP ⁺ [k−1]A^T+Q.  (24)

{circumflex over (x)} ⁺ [k]={circumflex over (x)} ⁻ [k]+K[k](y[k]−H{circumflex over (x)}⁻[k]),  (25)

S[k]=HP ⁻ [k]H ^T +HC+C ^T H ^T +R,  (26)

K[k]=(P⁻[k]H^T+C)S[k]⁻¹,  (27)

P ⁺ [k]=P ⁻ [k]−K[k](HP⁻[k]+C^T).  (28)

Equations (23) and (24) represent the dynamic update step, in which the estimated state, {circumflex over (x)}⁺[k−1], and estimated covariance, P⁺[k−1], are updated according to the system dynamics.

Equations (25)-(28) represent the measurement update step, in which the predicted state estimate, {circumflex over (x)}⁻[k], and the predicted covariance estimate, P⁻[k], are updated with the measurements to produce the current estimate, {circumflex over (x)}⁺[k] and P⁺[k].

Cyber-Physical Attack Model and Detection

The vectors f₁ and f₂ are introduced into state space model to represent additive state and measurement anomaly, or attack vectors, respectively:

x[k+1]=Ax[k]+Bu[k]+v[k]+f₁[k],  (29)

y[k]=Hx[k]+n[k]+f ₂ [k].  (30)

In terms of hypothesis testing, the aim of the present disclosure is to distinguish between

H ₀(No Attack, or Anomaly):∀k,f₁[k]=0,f₂[k]=0

H ₁(Attack, or Anomaly):There exist {k*} for which f₁└k*┘≠0,f₂└k*┘≠0.  (31)

Bad data from a faulty sensor can be viewed as one particular attack, or anomaly in this framework.

Our criterion for detecting attacks or anomalies is a statistic based on the output of the dynamic state estimation. It is assumed that the attack, or anomaly vectors f₁ and f₂ in Equations (29)-(30) are unknown. A sliding window attack statistic is proposed based on the measurement residual. The global attack statistic at time k is defined as follows:

$\begin{array}{cc}d\left[k\right]=\sum _{j=k-w+1}^k{\gamma }^T\left[j\right]S^{-1}\left[k\right]\gamma \left[j\right],& \left(32\right)\end{array}$

where the sliding window is of length w and m is the number of measurements. The measurement residual, i.e. Kalman innovation, is a zero mean Gaussian random variable, and the statistic d[k]˜χ² (win) is a chi-squared random variable with wm degrees of freedom. This is due to the fact that the covariance matrix of global innovation γ[k] is equal to S[k].

This fact can be proven as follows: the global innovations covariance matrix can be defined as:

$\begin{array}{cc}E\left[\gamma \left[k\right]{\gamma }^T\left[k\right]\right]=E\lfloor \left(y\left[k\right]-H{\hat{x}}^{-}\left[k\right]\right){\left(y\left[k\right]-H{\hat{x}}^{-}\left[k\right]\right)}^T\rfloor & \left(33a\right)\\ =E\left\{{\left[H\left(x\left[k\right]-{\hat{x}}^{-}\left[k\right]\right)+n\left[k\right]\right]\left[H\left(x\left[k\right]-{\hat{x}}^{-}\left[k\right]\right)+n\left[k\right]\right]}^T\right\}& \left(33b\right)\\ ={\mathrm{HP}}^{-}\left[k\right]H^T+R+E\left[H\left(x\left[k\right]-{\hat{x}}^{-}\left[k\right]\right)n^T\left[k\right]\right]+E\lfloor n\left[k\right]{\left(x\left[k\right]-{\hat{x}}^{-}\left[k\right]\right)}^TH^T\rfloor .& \left(33c\right)\end{array}$

From equation (23) we have that

$\begin{array}{cc}x\left[k\right]-{\hat{x}}^{-}\left[k\right]=\left(\mathrm{Ax}\left[k-1\right]+\mathrm{Bu}\left[k-1\right]+v\left[k-1\right]\right)-\left(A{\hat{x}}^{+}\left[k-1\right]+\mathrm{Bu}\left[k-1\right]\right)& \left(34a\right)\\ =A\left(x\left[k-1\right]-{\hat{x}}^{+}\left[k-1\right]\right)+v\left[k-1\right],& \left(34b\right)\\ E\left[\left(x\left[k\right]-{\hat{x}}^{-}\left[k\right]\right)n^T\left[k\right]\right]=E\left[A\left(x\left[k-1\right]-{\hat{x}}^{+}\left[k-1\right]\right)n^T\left[k\right]\right]+E\left[v\left[k-1\right]n^T\left[k\right]\right]=C,& \left(35\right)\end{array}$

where the last line follows from the fact measurement noise at time k is uncorrelated with the measurements and state at the previous time step. Plugging into (33c), the desired result can be obtained as:

E[γ[k]γ ^T [k]]=HP ⁻ [k]H ^T +R+HC+C ^T H ^T =S[k]  (36)

Distributed Implantation of Estimation and Online Attack Detection

According to embodiments of the present disclosure, the attack detection is implemented jointly with distributed dynamic state estimation. The distributed and dynamic nature of the invented method facilitates detecting attacks in an online fashion as new measurements become available, making this particularly attractive for monitoring the health of critical cyber-physical systems, such as power grids. In addition to eliminating the need for communication with a centralized control center, each control area is solved as a problem of reduced dimension with respect to the original global problem. In this aim, local states and local measurements are used for each control area.

As shown in Equation (27), the Kalman filtering needs an inversion for the innovations covariance matrix S[k]. In order to achieve distributed implementation of dynamic state estimation, it is solved using an iterative (distributed) inversion method instead of a direct (centralized) inversion method. To avoid using full communication scheme that is all generators communicate with each other, a proximity-based limited communication scheme, i.e. l-hop communication is used. That is only generators at most l-hops away communicate with each other. The parameter l is a pre-determined neighborhood distance.

Distributed Dynamic State Estimation

The present disclosure uses a criterion for detecting attacks that is a statistic based on the output of dynamic state estimation. Dynamic state estimation is distributed solved for each control area using local state and measurement models and limited communications with adjacent control areas.

Assumed the network buses are partitioned into a set of N control areas. The state local to control area I is the generator rotor angle and frequency at the buses in the control area I, x_I=[δ_I ω_I]^T. There is no overlap between states of neighboring areas.

The measurements local to control area I, y_I are the frequencies of all generators contained in the control area and the weighted phase angles θ_I* at buses contained in the control area. Based on the characteristics of matrix L^DD, the θ_I* for the border buses of control area I depends on measurements of phase angles at neighboring buses at neighboring control areas, so control area I needs exchanging the measurements of neighboring buses with their neighboring control areas.

One remarkable feature of the Kalman filter is that the estimated covariance matrices, such as P⁻[k], P⁺[k] and S[k] do not depend on the measurements. Therefore, they can be computed in advance offline. It is stressed here again that the system dynamics matrix A and measurement matrix H are both decoupled (i.e. no mixing is introduced between states in different areas). Assuming buses are labeled consecutively across control areas, therefore the state-space model can be written as:

x _I [k+1]=A_Ix_I[k]+B_Iu_I[k]+v_I[k],

y _I [k]=H _I x _I [k]+n _I [k]

I={1,2, . . . N}.  (37)

Since the power network is an interconnected system, it can expect that a need to exchange information between control areas. Indeed, this need is reflected in the fact that the Kalman gain K[k] is not a block-diagonal matrix. The innovation, or measurement residual, for control area I is defined as:

γ_I[k]y_I[k]−H_I{circumflex over (x)}_I⁻[k].  (38)

The measurement update to the local estimate is then

{circumflex over (x)} _I ⁺ [k]={circumflex over (x)} _I ⁺ [k]+K _I [k]γ[k],  (39)

where K_I[k] are the rows of K[k] corresponding to control area I. Since K[k] is not a block-diagonal matrix, control area I needs to communicate the entries of γ[k] with adjacent control areas to update its local estimate {circumflex over (x)}_I⁺±[k]. Because the formula for K[k] in (27), K[k]=(P⁻[k]H^T+C)S⁻¹[k] contains a matrix inversion, S[k]⁻¹, it is difficult to calculate in a distributed way if a direct inversion method is used. Instead, the present disclosure uses an iterative inversion method to achieve distributed implementation for dynamic state estimation.

Consider the linear system

S[k]β[k]=γ[k],  (40)

Then, the measurement update in (39) is given as

{circumflex over (x)} _I ⁺ [k]={circumflex over (x)} _I ⁺ [k]+(P⁻[k]H^T+C)_Iβ[k]  (41)

The key to dealing with the inverse in a distributed way is to iteratively solve (40) for β[k] without explicitly calculating the inverse and use the result in (41). The present disclosure uses the damped Jacobi method to achieve a fully distributed solution to (40).

The matrix S[k] can be decomposed into the difference of a diagonal matrix, S^D[k] and a matrix containing the remaining off-diagonal entries, S^E[k]:

S[k]=S ^D [k]−S ^E [k].  (42)

Iteratively solving for β[k] using the damped Jacobi method amounts to finding the fixed point of

β^t+1[k]=β^t[k]+α(S^D[k])⁻¹(γ[k]−S[k]β^t[k]),  (43)

where α is the damping parameter. Since S^D[k] is diagonal, its inverse (S^D[k])⁻¹ is diagonal, and each block can be computed locally. The sparsity of S[k] determines the entries from β^t[k] that need to be communicated with neighboring areas.

The damped Jacobi method in (43) can be proven to converge if the damping factor is selected according to

$\alpha <{\min }_i\frac{2S_{\mathrm{ii}}\left[k\right]}{\sum _{j\ne i}S_{\mathrm{ij}}\left[k\right]}.$

Since S[k] is the covariance matrix for the innovations, it is symmetric and positive semi-definite. Furthermore, by the standard assumptions for Kalman filtering, S[k] is an invertible matrix, thus S[k] positive definite, that is S[k]0. The damped Jacobi method can be converged if and only

$0\prec S\left[k\right]\prec \frac{2}{\alpha }S^D\left[k\right].$

Therefore, a sufficient condition is to choose α so that

$\left(\frac{2}{\alpha }S^D\left[k\right]-S\left[k\right]\right)$

is diagonally dominant. Diagonal dominance requires that

$\begin{array}{cc}{\left[\frac{2}{\alpha }S^D\left[k\right]-S\left[k\right]\right]}_{\mathrm{ii}}=\left(\frac{2}{\alpha }-1\right)S_{\mathrm{ii}}\left[k\right]>\sum _{j\ne i}{\left[\frac{2}{\alpha }S^D\left[k\right]-S\left[k\right]\right]}_{\mathrm{ij}}=\sum _{j\ne i}S_{\mathrm{ij}}\left[k\right],& \left(44\right)\end{array}$

Then we can have:

$\begin{array}{cc}\alpha <\frac{2S_{\mathrm{ii}}\left[k\right]}{\sum _{j\ne i}S_{\mathrm{ij}}\left[k\right]}\forall i.& \left(45\right)\end{array}$

The result utilizes the fact that S_ii[k]>0 since it is the value of a variance.

Since the damped Jacobi method is iterative, an inner-loop of iterations can be introduced for each outer-loop k of the Kalman filter. The number of inner-loop iterations T_inner is a tunable parameter that can be chosen to achieve a specified tolerance on the difference between consecutive inner iterations. After completing pre-determined number of iterations, T_inner, we obtain β^Tinner, which is multiplied by (P⁻[k]H^T+C) in order to calculate (41). Since P⁻[k] is in general a full matrix, in order to avoid communication between all generators, we propose the following masking approximation:

{tilde over (P)} ⁻ [k] N _l ⊗P ⁻ [k].  (46)

where N_l is a mask matrix for “l-hop” neighbors, and defined as:

$\begin{array}{cc}{\left[N_l\right]}_{\mathrm{ij}}\stackrel{\mathrm{def}}{=}\{\begin{array}{cc}1,& i,j\mathrm{are}l\text{-}\mathrm{hop}\mathrm{neighbors}\\ 0,& \mathrm{Otherwise}\end{array}& \left(47\right)\end{array}$

Thus, {tilde over (P)}_ij⁻[k] is nonzero if and only if the buses corresponding to x_i and x_j are at most l-hops away (e.g., direct neighbors are l-hop neighbors.) By tuning l, there is a tradeoff between accuracy of estimation and communication requirements.

The sparsity of S[k] in (26) determines the entries from β^t[k] that need to be communicated with neighboring areas. We analyze the communication requirements of our distributed method in terms of the sparsity patterns of relevant matrices and the l-hop neighbor approximation in (47).

Using (43), β_I is iteratively solved, then neighbors need to communicate their entries of the vector β according to the sparsity pattern of S[k]. For l≥2, using a l-hop neighbor mask in (47), the sparsity pattern of S[k] has non-zero entries only at pairs of measurements corresponding to buses that are at most l-hops away. Therefore, at most l-hop neighbor communication is needed for each control area.

This can be proved as follows: The matrix S[k] consists of the sum of four terms. Consider the sparsity pattern of the first item, H{tilde over (P)}⁻[k]H^T:

$\begin{array}{cc}{\left(H{\tilde{P}}^{-}\left[k\right]H^T\right)}_{\mathrm{ij}}=\sum _m^{}H_{{\theta }_i^{*},m}\sum _l^{}{\tilde{P}}_{ml}^{-}\left[k\right]H_{{\omega }_j,l}& \mathrm{ }\left(49\right)\\ =\sum _m^{}H_{{\theta }_i^{*},m}{\tilde{P}}_{m,{\omega }_j}\left[k\right]=H_{{\theta }_i^{*},{\delta }_i}{\tilde{P}}_{{\delta }_i,{\omega }_j}\left[k\right]& \left(50\right)\end{array}$

For concreteness, take row i of H to correspond to the measurement θ_i* and take row j to correspond to the measurement {circumflex over (ω)}_j.

$\begin{array}{cc}{\left(H{\tilde{P}}^{-}\left[k\right]H^T\right)}_{\mathrm{ij}}=\sum _m^{}H_{\mathrm{im}}\sum _l^{}{\tilde{P}}_{ml}^{-}\left[k\right]H_{\mathrm{jl}}.& \left(48\right)\end{array}$

is nonzero if and only if {tilde over (P)}_δi,ωj is not zero. From the definition of l-hop mask in (47), this is equivalent to having bus i and bus j be at most l-hop neighbors. A similar argument follows for the other measurement types. In conclusion, the calculation of (H{tilde over (P)}⁻[k]H^T) requires l-hop neighbor communication. In contrast, without the masking approximation, the estimate of the state covariance matrix P⁻ is full, and each generator would need to communicate with every other generator. The second item is the measurement noise covariance matrix:

$\begin{array}{cc}R=\left[\begin{array}{cc}{R}_{\omega }& 0\\ 0& L^{\mathrm{DD}}{R_{\theta }\left(L^{\mathrm{DD}}\right)}^T\end{array}\right],& \left(51\right)\end{array}$

where R_ω and R_θ are assumed to be diagonal. The entry L_ij^DD is nonzero if and only bus i and bus j are 1-hop neighbors, and [L^DDR_θ(L^DD)^T]_ij is nonzero if and only bus i and bus j are at most 2-hop neighbors. Therefore, R requires at most 2-hop neighbor communication. The third and fourth items are related to analyze the matrix HC. The entry [HC]_ij is nonzero if and only if the states that measurement i depends upon have overlap with the control inputs correlated to measurement j. Measurement i can either be {circumflex over (ω)}_i or {circumflex over (θ)}_i* which depends on ω_i or θ_i, respectively. From (51), a measurement of ω_j is not correlated to the control inputs, so [HC]_ij is zero for columns j corresponding to measurements of ω. Measurements of θ_i* are correlated with the control input at bus i and its 1-hop neighboring buses. Therefore, HC only requires neighbor-to-neighbor communication. The same argument holds for C^TH^T. In summary, at most l-hop neighbor communication is needed due to calculation of (H{tilde over (P)}⁻[k]H^T).

For Equation (41), there are additional communication requirements for calculating ({tilde over (P)}⁻[k]H^T+C) of the method, after calculating β. The sparsity of P⁻H^T is such that the columns corresponding to measurements of {circumflex over (θ)}* at a load bus are zero. In order to calculate the entry of vector └P⁻¹[k]H^Tβ┘ corresponding to measurement {circumflex over (ω)}_i the entries of β corresponding to the measurements of ω and θ* at all other generators are needed. If an l-hop neighbor mask is applied, then only the entries of β corresponding to measurements at generators at most l-hops away are needed. In summary, to approximately calculate local entries of the vector [P⁻¹[k]H^T+C]β, areas must communicate local entries of β with at most their l-hop neighbors.

Distributed Attack Detection

During the dynamic state estimation, the quantity β_I[j]=[S⁻¹[k]γ[j]]_I is calculated locally in each control area I. Therefore no additional communication is required to calculate the local attack statistic

$\begin{array}{cc}{d}_I\left[k\right]=\sum _{j=k-w+1}^k{\gamma }_I^T\left[j\right]{\beta }_I\left[j\right].& \left(52\right)\end{array}$

If one had access to the global detection statistic, a classic chi-squared detection test could be used. For real-time attack detection in large networks, it is not feasible to collect

$\sum _I^{}d_I\left[k\right]$

over all areas. Instead, the present disclosure proposes that each area base its attack detection on its local attack statistic information. If S⁻¹[k] were block diagonal, then d_I[k] would be distributed as a chi-squared random variable with wm_I degrees of freedom, where m_I is the number of measurements in area I. However, in general S⁻¹[k] is a full matrix, and d_I[k] does not have an easily characterized distribution. Instead, this invent compares the real-time local attack statistic with an analytical mean and variance of d_I[k] under hypothesis H₀ (No Attack) and sliding window w=1 to determine the possible cyber-attack. Using w=1, the mean and variance of the local attack statistic without an attack can be quantified as follows:

$\begin{array}{cc}E\left[d_I\left[k\right]\right]=\sum _{i\in I}^{}\sum _{l=1}^nS_{\mathrm{il}}^{-1}\left[k\right]S_{\mathrm{il}}\left[k\right],& \left(53a\right)\\ \mathrm{var}\left[d_I\left[k\right]\right]=E\left[d_I^2\left[k\right]\right]-({E\left[d_I\left[k\right]\right)}^2,& \left(53b\right)\\ \mathrm{where}& \\ E\left[d_I^2\left[k\right]\right]=\sum _{i,j\in I}^{}\sum _{l,m=1}^nS_{\mathrm{il}}^{-1}\left[k\right]S_{\mathrm{jm}}^{-1}\left[k\right]\left(S_{\mathrm{il}}\left[k\right]S_{\mathrm{jm}}\left[k\right]+S_{\mathrm{ij}}\left[k\right]S_{lm}\left[k\right]+S_{\mathrm{im}}\left[k\right]S_{\mathrm{lj}}\left[k\right]\right).& \left(53c\right)\end{array}$

Given the analytical value for the variance, a threshold τ_I is set such that if |d_I[k]|>τ_I(var[d_I[k]])| an attack is declared. For example, the threshold can be a multiple of Var(d_I). A nice feature of the invented method is that different false alarm probabilities can be set per area based on the areas' noise characteristics, and extra information about the location is available since the local partial sums of the global attack variable is monitored in the present disclosure.

Procedure for Distributed Estimation and Attack Detection

FIG. 3A is a schematic block diagram for distributed estimation and detection of cyber-physical attacks for predetermined period of covariance matrix update; and FIG. 3B is a schematic block diagram for determining local area attack statistic and declare possible attack;

In FIG. 3A, the distributed estimation and attack detection of a power system is implemented through the following steps:

Step-310: determine the system dynamics matrix A, measurement relationship matrix H, state and measurement correlation matrix C and state covariance matrix Q and measurement covariance matrix R.

Step-320: Set the total number of time step for covariance update of attack detection K, the threshold number of multiple of no attack variance for each area I, τ_I, pre-determined neighborhood distance l, and initialize covariance matrix P⁺[0]. The threshold τ_I is set per area based on the areas' noise characteristics, and extra information about the location is available.

Step-330: Calculate offline the covariance matrices, including {P⁻[k]}_k=0^K, {P⁺[k]}_k=0^K, and {S[k]}_k=0^K, and communicated to each control area according to:

P ⁻ [k]=AP ⁺ [k−1]A^T+Q

S[k]=HP ⁻ [k]H ^T +HC+C ^T H ^T +R

P ⁺ [k]=P ⁻ [k]−(P⁻[k]H^T+C)S⁻¹[k](HP⁻[k]+C^T)

Step-340: execute on-line attack detection based on a local cyber attack static for each control area for each detection step. Each time step k determines local attack statistic and declare an attack for each area following the steps described in FIG. 3B:

Step-345: Calculate predicted state estimate for each control area I at time step k:

x _I ⁻ [k]=A _I {circumflex over (x)} _I ⁺ [k−1]+B_Iu_I[k−1]

This step does not require any communication since A is block diagonal and [u]_I only depends on local information.

Step-350: Read measurements for local control area and communicate with neighboring areas to get measurements for neighboring buses within 1-hop neighborhood distance in adjacent control areas.

Step-360: Iterate T_inner steps to determine covariance weight vector β_I for each control area and communicate with its l-hop neighbors of covariance weight vectors during the iterations according to:

β_I^t+1=β_I^t+α[S^D[k]]_I⁻¹(γ_I−S_I,I∪NIl[k]β_I∪NIl^t),

wherein N_I^l is the set of neighboring buses of buses in the control area I within l-hop neighborhood at adjacent control areas, l is a pre-determined neighborhood distance. S_I,I∪NIl[k] is the entries of S[k] at rows corresponding to control area I, and columns corresponding to control area I and its l-hop neighbors, N_I^l. β_I∪NIl^t is the entries of β^t corresponding to control area I and its l-hop neighbors, N_I^l.

Step-370: Using a neighbor-limited approximation to P⁻[k], {tilde over (P)}⁻[k] and communication β_I^Tinner with its l-hop neighbors, each area calculates

{circumflex over (x)} _I ⁺ [k]={circumflex over (x)} _I ⁻ [k]+({tilde over (P)}⁻[k]H^T+C)_I,I∪NIlβ_I∪NIl^Tinner.

wherein {tilde over (P)}_ij⁻[k] is nonzero and equals to P_ij⁻[k], if and only if the buses corresponding to x_i and x_j are at most l-hops away. ({tilde over (P)}⁻[k]H^T+C)_I,I∪NIl is the entries of ({tilde over (P)}⁻[k]H^T+C) at rows corresponding to control area I, and columns corresponding to control area I and its 1-hop neighbors, N_I^l.

Step-380: Use β_I to update the local attack detection statistic d_I[k]:

$d_I\left[k\right]=\sum _{j=k-w+1}^k{\gamma }_I^T\left[j\right]{\beta }_I\left[j\right],$

wherein w is the length of sliding window and set as 1.

Step-390: check if |d_I[k]|>τ_IVar(d_I[k]). Var(d_I[k]) is a standard variance of the local attack statistic without an attack, and defined according to:

$\mathrm{var}\left[d_I\left[k\right]\right]=E\left[d_I^2\left[k\right]\right]-({E\left[d_I\left[k\right]\right)}^2E\left[d_I\left[k\right]\right]=\sum _{i\in I}^{}\sum _{l=1}^nS_{\mathrm{il}}^{-1}\left[k\right]S_{\mathrm{il}}\left[k\right]E\left[d_I^2\left[k\right]\right]=\sum _{i,j\in I}^{}\sum _{l,m=1}^nS_{\mathrm{il}}^{-1}\left[k\right]S_{\mathrm{jm}}^{-1}\left[k\right]\left(S_{\mathrm{il}}\left[k\right]S_{\mathrm{jm}}\left[k\right]+S_{\mathrm{ij}}\left[k\right]S_{lm}\left[k\right]+S_{\mathrm{im}}\left[k\right]S_{\mathrm{lj}}\left[k\right]\right)$

wherein S⁻¹[k] is the inverse of S[k]. If yes, an attack is declared.

Example

FIG. 4A-4F gives an example for the local attack detection statistic behavior in the system shown in FIG. 1. In each of the six figures of FIG. 4, the histogram of the attack statistic without an attack present and with an attack present is overlaid. Each event in the histogram corresponds to different values for the correlated process and measurement noise. The figures give the attack statistic at three different time steps (before, during, and after the attack) and in two different control areas. The attack is a corruption of the signal reading the power demand at bus 120 in Control Area I in FIG. 1A. The power demand at bus 120 is taken to be multiple times its regular value in the state estimation. The threshold of multiple for no-attack variance τ_I is 1.

FIG. 4 A-4C show the attack statistics in Control Area I, where the attack takes place. The analytic values for the mean and variance are given using the covariance matrix from our simulations. The vertical line is the analytic mean, and the dashed lines are at plus and minus 1 analytic standard deviation. As shown in the plots, before the attack the histogram for d_I matches exactly with or without an attack. At the time step where the attack occurs, the variance of the statistic d_I is increased with respect to the case when no attack is present.

FIG. 4D-4F show the histogram for the attack statistic d_IV in neighboring control area IV. It can be seen that during the attack, the variance of d_IV is only slightly increased. This points to another potential advantage of using the local attack statistic d_I in (52) rather than the global attack statistic d in (32). Since the variance of the local attack statistic where the attack is taking place is increased with respect to the local statistics in other areas, this suggests that using the local attack statistic helps not only in detecting the presence of an attack but also in identifying where the attack is taking place.

Using the method described above, it can be determined that a possible attack is at control area I, but control area IV.

According to embodiments of the present disclosure a method for detecting anomalies or cyber attacks in a control system or an electric power system. The electric power system includes a set of generators and loads connected to buses, and each pair of buses are connected with branches. Wherein a voltage phase angle of each bus and a rotor frequency of each generator are measured at a pre-determined sampling frequency. Wherein the electric power system is partitioned into a set of control areas in which no common buses exist between the control areas. Wherein a mechanical power for each generator and an active power demand for each load are given for each sampling interval corresponding to the sampling frequency. The method comprising: determining a local state model and a local measurement model for each control area based on measurements at the control area and neighboring buses of adjacent control areas within a pre-determined neighborhood distance. Determining dynamic state estimations for each control area for each sampling interval using iterative inversion solution for covariance matrix and communications among adjacent control areas within the pre-determined neighborhood distance. Determining a local anomalies or cyber-attack statistic for each control area according to results of determined dynamic state estimations. Estimating the possibility for a possible anomalies or cyber attack in the control area by comparing the determined local cyber-attack statistic with a variance of measurement residuals for the control area without anomalies or cyber attack, and declaring possible anomalies, or cyber attacks, when the local statistic is greater than pre-determined multiples of non-anomalies or the non-cyber attack variance of measurement residuals.

According to aspects of the method, wherein the local state for control area I is the generator rotor angle and frequency at the generators in the control area I, x_I=[δ_I ω_I]^T. Wherein x_I is the vector of local states, δ_I and ω_I are the vectors of rotor angle and frequency of all generators in the area; wherein the local measurements for control area I, y_I=[{circumflex over (ω)}_I θ_I*]^T, {circumflex over (ω)}_I are the rotor frequencies of all generators contained in the control area, and {circumflex over (δ)}_I* are the weighted bus phase angles for all buses contained in the control area including the terminal buses of generators. Wherein the weighted bus phase angle at bus i is determined as a linear combination of voltage phase angles at neighboring buses and an local electric power demand P_i^D:

{circumflex over (θ)}_i*L_i^DD{circumflex over (θ)}−P_i^D,

L^DD is a matrix related bus power injection to bus phase angle for all buses in the system according to the susceptances of transmission lines and transient reactance of generators in the system; L_i^DD is the row corresponding to bus i of matrix L^DD; {circumflex over (θ)} is the vector of voltage phase angle for all buses; {circumflex over (θ)}_i* is determined using measurements at the control area and measurements with their neighbors in other control areas within 1-hop neighborhood distance.

According to other aspects of the method, further comprising: determining the system dynamics matrix A, input control matrix B, measurement relationship matrix H, state and measurement correlation matrix C, state covariance matrix Q, and measurement covariance matrix R. Specifying the number of time steps for a period of covariance update K, the multiple threshold of standard variance of no-attack statistic for each control area τ_I, I={1, . . . ,N}, N is the total number of control areas. Calculating offline the matrices of predicted state covariance {P⁻[k]}_k=0^K, estimated state covariance {P⁺[k]}_k=0^K, and innovation covariance {S[k]}_k=0^K for each time step k in the set of time steps of covariance update period, and communicating the matrices to each control area. Executing on-line cyber-physical attack detection based on a local cyber attack statistic for each control area I, I={1, . . . , N} at each time step k, k∈{1, 2, . . . , K}.

According to other aspects of the method, wherein the system dynamics matrix A, input control matrix B, state covariance matrix Q and measurement covariance matrix R are determined according to:

$A=e^{\stackrel{\_}{\mathrm{TA}}},B={\int }_0^{\stackrel{\_}{T}}e^{\tau \stackrel{\_}{A}}d\tau ,Q={\int }_0^{\stackrel{\_}{T}}a^{\stackrel{\_}{A}\tau }\stackrel{\_}{Q}e^{{\stackrel{\_}{A}}^T\tau }d\tau ,R=\left[\begin{array}{cc}{R}_{\omega }& 0\\ 0& L^{\mathrm{DD}}{R_{\theta }\left(L^{\mathrm{DD}}\right)}^T\end{array}\right],$

wherein Q is sampled process noise covariance matrix, Q is pre-determined un-sampled covariance matrix of states, T is the length of each sampling interval of phase measuring units, Ā is an un-sampled matrix of system dynamics related the state derivatives {dot over (x)}(t) with states x(t), A is a sampled matrix of system dynamics related the states at next time step with states at current time steps, Ā are set with non-zero values at

$\stackrel{\_}{A}\left({\stackrel{.}{\delta }}_i\left(t\right),{\omega }_i\left(t\right)\right)=1.0,\stackrel{\_}{A}\left({\stackrel{.}{\omega }}_i\left(t\right),{\delta }_i\left(t\right)\right)=\frac{-1}{M_iZ_i},\stackrel{\_}{A}\left({\stackrel{.}{\omega }}_i\left(t\right),{\omega }_i\left(t\right)\right)=-\frac{D_i}{M_i}{\omega }_i\left(t\right)·A^T$

is the transpose of matrix A, B is an sampled matrix related the states with system inputs at the current time steps, R_ω and R_θ are the pre-determined covariance matrices for measurements of rotor frequency and bus phase angle at all generators and all buses; wherein the state x(t) includes the generator variables for all generators, x(t)=└x₁(t) x₂(t) . . . x_ng(t)┘^T, and the state for generator i is defined as x_i(t)=[δ_i(t) ω_i(t)]^T, and the control input u(t) includes inputs for all generators, u(t)=└u_δ1(t) u_ω1(t) . . . u_δng(t) u_ωng(t)┘, and the control input corresponding to the rotor angle and frequency of the generator at bus i are defined as:

$u_{{\delta }_i}\left(t\right)=0,u_{{\omega }_i}\left(t\right)=\frac{1}{M_i}\left(P_i^G+\frac{1}{Z_i}{\theta }_i\right),$

where n_g is the number of generators in the system, M_i, Z_i and P_i^G are the inertia, the transient reactance, and the mechanical power input of the generator at bus i; wherein the entries of measurement correlation matrix C are set with non-zero values at C(u_i,δ, {circumflex over (θ)}_j*) and C(u_i,ω, {circumflex over (θ)}_j*) if j=i or j∈N_i, N_i is the set of 1-hop neighboring buses of bus i; wherein entries of measurement relationship matrix H are set with non-zero values at H({circumflex over (ω)}_i,ω_i)=1.0 and

$H\left({\hat{\theta }}_i^{*},{\delta }_i\right)=\frac{1}{Z_i},i\in G;$

According to other aspects of the method, the matrices of predicted state covariance, estimated state covariance, and innovation covariance for time step k are determined according to:

P ⁻ [k]=AP ⁺ [k−1]A^T+Q,

S[k]=HP ⁻ [k]H ^T +HC+C ^T H ^T +R,

P ⁺ [k]=P ⁻ [k]−(P⁻[k]H^T+C)S⁻¹[k](HP⁻[k]+C^T),

k∈{1, . . . ,K};

According to other aspects of the method, determining local attack statistic and declare an attack for each area at time step k, further comprising: calculating predicted state estimate for each control area I at time step k. Reading measurements for local control area and communicate with neighboring areas to get measurements for neighboring buses within 1-hop neighbors. Iterating T_inner steps to determine covariance weight vector β_I for each control area and communicate with its l-hop neighbors of covariance weight vectors during the iterations. Using a neighbor-limited approximation to P⁻[k], {tilde over (P)}⁻[k] and communication β_I^Tinner, with its l-hop neighbors to determine estimated state for each control area I at time step k. Using β_I^Tinner to update the local attack detection statistic d_I[k]. Checking if |d_I[k]|>τ_I(Var(d_I[k])). If yes, an attack is declared.

According to other aspects of the method, wherein the predicted state estimate for each control area I at time step k are determined according to:

x _I ⁻ [k]=A _I {circumflex over (x)} _I ⁺ [k−1]+B_Iu_I[k−1]

x_I⁻[k] and {circumflex over (x)}_I⁻[k−1] are the vectors of predicted local state at time step k, and estimated local state at time step (k−1); A_I and B_I are the sub-matrices of system dynamics and measurement matrices corresponding to the buses of control area I, u_I[k−1] is the vector of control inputs of control area I at time step (k−1);

According to other aspects of the method, wherein the weight vector of innovation covariance for control area I, β_I is solved iteratively according to:

β_I^t+1=β_I^t+α[S^D[k]]_I⁻¹(γ_I−S_I,I∪NIl[k]β_I∪NIl^t)

t={0,1, . . . ,T_inner}

wherein T_inner is the total number of iterations, S^D[k] is the diagonals of matrix S[k], α is the damping parameter and determined as a value less than the minimum of

$\frac{2S_{\mathrm{ii}}\left[k\right]}{\sum _{j\ne i}^{}S_{\mathrm{ij}}\left[k\right]}$

for all row i of matrix S[k]. N_I^l is the set of neighboring buses of buses in the control area I within l-hop neighborhood distance. S_I,I∪NIl[k] is the entries of S[k] at rows corresponding to control area I, and columns corresponding to control area I and its l-hop neighbors, N_I^l. β_I∪NIl^l is the entries of β^t corresponding to control area I and its l-hop neighbors, N_I^l.

According to other aspects of the method, the estimated state for control area I is determined according to:

{circumflex over (x)} _I ⁺ [k]={circumflex over (x)} _I ⁻ [k]+({tilde over (P)}⁻[k]H^T+)_I,I∪NIlβ_I∪NIl^Tinner.

wherein {tilde over (P)}⁻[k] is a neighbor-limited approximation to P⁻[k], the entry of {tilde over (P)}⁻[k], {tilde over (P)}_ij⁻[k] is nonzero and equals to P_ij⁻[k], if and only if the buses corresponding to x_i and x_j are at most 1-hops away. ({tilde over (P)}⁻[k]H^T+C)_I,I∪NIl, is the entries of ({tilde over (P)}⁻[k]H^T+C) at rows corresponding to control area I, and columns corresponding to control area I and its 1-hop neighbors, N_I^l. β^Tinner is determined the weight vector of innovation covariance.

According to other aspects of the method, wherein the local attack detection statistic d_I[k] is determined according to:

$d_I\left[k\right]=\sum _{j=k-w+1}^k{\gamma }_I^T\left[j\right]{\beta }_I\left[j\right];$

wherein γ_I[k] is the innovation vector of measurement residual determined as:

γ_I[k]y_I[k]−H_I{circumflex over (x)}_I⁻[k]

wherein y_I[k] is the vector of local measurements for control area I at time step k, {circumflex over (x)}_I⁻[k] is the vector of predicted states for control area I at time step k, w is the slide window size, and set as 1;

According to other aspects of the method, wherein the standard variance of the local attack statistic without an attack, Var(d_I[k]) is determined according to:

$\mathrm{var}\left[d_I\left[k\right]\right]=E\left[d_I^2\left[k\right]\right]-({E\left[d_I\left[k\right]\right)}^2E\left[d_I\left[k\right]\right]=\sum _{i\in I}^{}\sum _{l=1}^nS_{\mathrm{il}}^{-1}\left[k\right]S_{\mathrm{il}}\left[k\right]E\left[d_I^2\left[k\right]\right]=\sum _{i,j\in I}^{}\sum _{l,m=1}^nS_{\mathrm{il}}^{-1}\left[k\right]S_{\mathrm{jm}}^{-1}\left[k\right]\left(S_{\mathrm{il}}\left[k\right]S_{\mathrm{jm}}\left[k\right]+S_{\mathrm{ij}}\left[k\right]S_{lm}\left[k\right]+S_{\mathrm{im}}\left[k\right]S_{\mathrm{lj}}\left[k\right]\right)$

Wherein S⁻¹[k] is the inverse of innovation covariance matrix S[k] at time step k, n is the total number of buses in the system.

According to a system for detecting cyber-physical attacks in an electric power system, wherein the electric power system includes a set of generators and loads connected to buses, and each pair of buses are connected with branches. Wherein a voltage phase angle of each bus and a rotor frequency of each generator are measured at a sampling frequency. Wherein the electric power system is partitioned into a set of control areas in which no common buses exist between the control areas. Wherein mechanical power for each generator and active power demand for each load are given for each sampling interval corresponding to the pre-determined sampling frequency. The system comprising: determining a local state model and a local measurement model for each control area based on measurements at the control area and neighboring buses of adjacent control areas within 1-hop neighborhood distance. Determining dynamic state estimations for each control area for each sampling interval using iterative inversion solution for state covariance matrix and communications among adjacent control areas within the pre-determined neighborhood distance. Determining a local cyber-attack statistic according to results of determined dynamic state estimations. Declaring a possible cyber attack in the control area by comparing the determined cyber-attack statistic with a variance of measurement residuals for the control area without an attack.

Claims

1. A method for detecting anomalies in a control area of a control system with multiple control areas, wherein the control area includes a set of generators in communication with a set of buses and some neighboring buses of neighboring controls area of the control system, comprising: accessing a memory in communication with at least one processor, wherein the memory includes stored historical states of the control area;estimating, by the processor, a first state of the control area over a first state time period from the historical states of the control area, using a model of dynamics of the control area, and defining a transition of the first state of the control area as a function of control inputs, wherein the first state of the control area includes a generator state for each generator in the control area, wherein the control inputs include one or combination of: a network state for each bus in the control area; a network state for some neighboring buses of neighboring controls areas; a mechanical input to each generator in the control area; or power consumptions at the buses in the control area;updating, by the processor, the estimated first state of the control area according to a measurement model of the control area, by: connecting measurements of the rotor frequency of each generator; a weighted combination of measurements of the network states on the buses in the control area and some neighboring buses of neighboring controls areas; with the generator state of each generator in the control area, to obtain a second state of the control area over a second state time period later than the first state time period; anddetecting, by the processor, the anomalies based on a statistic deviation of the second state of the control area from its corresponding prediction derived from the first state of the control area.

2. The method of claim 1, wherein the generator state includes a rotor angle and a frequency of each generator in the set of generators in the control area.

3. The method of claim 1, wherein the network state includes a voltage phase angle at each bus in the set of buses in the control area.

4. The method of claim 1, wherein the first state time period and the second state time period each cover approximately the same time interval.

5. The method of claim 1, wherein the control system is an electric power system (EPS).

6. The method of claim 1, wherein the anomalies include cyber-attacks and power consumption anomalies.

7. The method of claim 1, wherein the second state of the control area are real-time observations.

8. The method of claim 1, further comprising: determining whether, the statistic deviation of the second state from its corresponding prediction derived from the first state of the control area is associated with a cyber attack based, at least in part, on historical cyber attack data stored in a memory in communication with the processor.

9. The method of claim 1, further comprising: determining whether the statistic deviation of the second state from its corresponding prediction derived from the first state of the control area is associated with a cyber attack based, at least in part, on historical anomalies associated with historical states of the control area indicating there is not a cyber attack.

10. The method of claim 9, further comprising: receiving an input from a surface of a user interface in communication with the processor, by a user, of the historical anomalies associated with historical states of the control area that are not indicative of a cyber attack.

11. The method of claim 1, further comprising: determining, by the processor, whether the statistic deviation of the second state from its corresponding prediction derived from the first state of the control area is associated with a cyber attack based, at least in part, on exceeding a cyber attack threshold;providing a notification responsive to the statistic deviation of the second state from its corresponding prediction derived from the first state of the control area exceeding the cyber attack threshold.

12. The method of claim 1, wherein the model of dynamics of the control area includes using an iterative inversion solution for a covariance matrix and communications among adjacent control areas within the pre-determined neighborhood distance.

13. The method of claim 1, wherein the measurement model of the control area is based on measurements at the control area and neighboring buses of adjacent control areas within a pre-determined neighborhood distance.

14. The method of claim 1, wherein the detecting of the anomalies includes: determining a function of a covariance matrix of the deviation of the second state from its corresponding prediction derived from the first state of the control area;comparing the function of the covariance matrix with a threshold determined as a multiple of a variance of the function of the covariance matrix allowed by uncertainties of the model of dynamics and the model of measurements; anddetecting a cyber-attack when the function of the covariance matrix is above the threshold.

15. A system for detecting anomalies in a control area of a power system with multiple control areas, wherein the control area includes a set of generators in communication with a set of buses and some neighboring buses of neighboring controls area of the control system, comprising: a computer readable memory to store historical states of the control area, current states of the control area, a model of dynamics of the control area and a measurement model of the control area;a set of sensors arranged to sense measurements in the control area for the set of generators, the set of buses and to measure one or combination of rotor frequencies for each generator, voltage phase angles on the buses, a mechanical input to each generator, or power consumptions at the buses;a processor in communication with the computer readable memory configured to: estimate a first state of the control area from a historical state of the control area over a first state time period using the stored model of dynamics of the control area, and defining a transition of the first state of the control area as a function of control inputs, wherein the first state of the control area includes a rotor angle and a rotor frequency for each generator in the control area, wherein the control inputs include one or combination of phase angles on the buses of the control area, a mechanical input to each generator in the control area, and power consumptions at the buses in the control area;update the estimated first state of the control area according to the stored measurement model of the control area, by connecting measurements of the rotor frequency of each generator in the control area and a weighted combination of measurements of the voltage phase angles on the buses of the control area and some neighboring buses of the neighboring control areas, with the rotor angle and the rotor frequency of each generator in the control area, to obtain a second state of the control area over a second state time period later than the first state time period; anddetect the anomalies based on a statistic deviation of the second state from its corresponding prediction derived from the first state of the control area.

16. The system of claim 15, wherein the first state time period and the second state time period each cover approximately the same time interval and the second state of the control area are real-time observations.

17. The system of claim 15, wherein the detected anomalies assist in managing the management of the EPS.

18. The method of claim 15, wherein the detecting of the anomalies includes: determining a function of a covariance matrix of the deviation of the second state from its corresponding prediction derived from the first state of the control area;comparing the function of the covariance matrix with a threshold determined as a multiple of a variance of the function of the covariance matrix allowed by uncertainties of the model of dynamics and the model of measurements; anddetecting a cyber-attack when the function of the covariance matrix is above the threshold.

19. A detector for detecting anomalies in a control area of an electric power system (EPS) with multiple control areas, wherein the control area includes a set of generators in communication with a set of buses and some neighboring buses of neighboring controls area of the control system, comprising: acquire a plurality of measurements from sensors configure for sensing the set of generators and the set of buses over a first state time period;acquire a respective second plurality of measurements from sensors configure for sensing the set of generators and the set of buses over a second state time period;at least one processor having a computer readable memory configured to: receive the plurality of measurements sensed by the sensors over the first state time period;estimate a first state of the control area from a historical state of the control area over the first state time period using a model of dynamics of the control area, and defining a transition of the first state of the control area as a function of control inputs, wherein the first state of the control area includes a rotor angle and a rotor frequency for each generator in the control area, wherein the control inputs include one or combination of phase angles on the buses of the control area, a mechanical input to each generator in the control area, and power consumptions at the buses in the control area;receive the plurality of measurements sensed by the sensors over the second state time period;update the estimated first state of the control area according to a measurement model of the control area, by connecting sensed measurements of the rotor frequency of each generator in the control area and a weighted combination of sensed measurements of the phase angles on the buses of the control area and some neighboring buses of the neighboring control areas, with the rotor angle and the rotor frequency of each generator in the control area, to obtain a second state of the control area over the second state time period later than the first state time period; anddetect the anomalies based on a statistic deviation of the second state of the control area from its corresponding prediction derived from the first state of the control area.

20. The detector of claim 19, wherein the first state time period and the second state time period each cover approximately the same time interval and the second state of the control area are real-time observations.