Detection of false data injection attacks using unscented Kalman filter

Abstract

It has recently been shown that state estimation (SE), which is the most important real-time function in modern energy management systems (EMSs), is vulnerable to false data injection attacks, due to the undetectability of those attacks using standard bad data detection techniques, which are typically based on normalized measurement residuals. Therefore, it is of the utmost importance to develop novel and efficient methods that are capable of detecting such malicious attacks. In this paper, we propose using the unscented Kalman filter (UKF) in conjunction with a weighted least square (WLS) based SE algorithm in real-time, to detect discrepancies between SV estimates and, as a consequence, to identify false data attacks. After an attack is detected and an appropriate alarm is raised, an operator can take actions to prevent or minimize the potential consequences. The proposed algorithm was successfully tested on benchmark IEEE 14-bus and 300-bus test systems, making it suitable for implementation in commercial EMS software.

1 Introduction

Due to the ever-increasing reliance on modern cyber infrastructures in power systems, cyber security has recently been considered to be among the most important issues in modern power systems. Supervisory control and data acquisition (SCADA) systems are vulnerable to attacks that are directed not only at data communication infrastructures but also to those directed at control centers and even remote terminal units (RTUs). The SCADA communication network is very diverse and consists of fiber optics, microwave and satellite connections, while the exchanged data is often unencrypted - leaving substantial space for potential attacks. Even though phasor measurement unit (PMU) based measurements are regarded as generally more secure than SCADA measurements, they are also susceptible to malicious attacks, as explained in [1]. For example, it was stated in [2] that among 245 reported incidents - the energy sector led all others with 79 (32%) incidents. A large number of those threats targeted SCADA devices with the intention of gaining unauthorized access. Earlier, an experimental attack was conducted by the researchers at the US Department of Energy’s Idaho Lab, which resulted in the self-destruction of one generating unit, emphasizing the impact of cyber threats [3]. State estimation (SE), which is among the most important real-time applications of commercial energy management system (EMS) software, is thus indirectly vulnerable to cyber-attacks, which may result in severe system instability, suboptimal operation, financial losses, and even loss of human life.

To identify and remove bad measurements—so they cannot compromise the accuracy of estimated results, bad data detection (BDD) algorithms have been devised as an integral part of SE algorithms [4]. In the early BDD algorithms, after the existence of bad measurements was confirmed, only the single measurement with the largest normalized residual was removed in one iteration of the SE algorithm. Multiple loops were needed until all bad measurements were eliminated and the detection test was passed. Later, a new approach was devised, in which multiple bad measurements were identified and removed at the same time, thus substantially reducing the SE execution time [4]. Both of these algorithms, which are presented in [4], rely on normalized residual techniques and show good results only in cases of independent and non-interacting bad measurements.

For correlated bad measurements, novel techniques have been developed, including hypothesis testing identification (HTI), which again relies on normalized residuals [5]. This and similar techniques are still among the most commonly used in commercial EMS software. However, the concept of “false data injection (FDI)” attacks was recently introduced in [6], proving that measurements can be manipulated in a way that does not trigger SE BDD modules. Given the knowledge of the system configuration and element parameters, an attacker can create undetectable malicious attacks. The development of techniques for identification and mitigation of such attacks is essential for the secure operation of large power grids.

The substantial amount of research carried out on FDI attacks can be classified into three categories.

• 1) Vulnerability of SE. The weaknesses of BDD algorithms have been explored from the perspective of an attacker as well as ways to construct malicious attacks with minimal resources and maximum impacts on power grids [6,7,8,9,10,11,12,13,14,15,16,17,18]. Although network topology and electric parameters are regarded as completely or partially known by the attacker, in some research [6, 10], they were not required, whereas others relied only on the measurement matrix [13,14,15,16,17,18].

• 2) Consequences of an FDI attack. FDI attacks on SE have been analyzed from the perspective of EMS applications such as contingency analysis, optimal power flow and automatic generation control [19,20,21,22].

• 3) Countermeasures development. These studies have focused on the detection of stealthy attacks and protection of the power system. New enhanced BDD algorithms have been proposed, optimal PMU placement strategies created, and discussions undertaken regarding ways to improve the protection of communication systems [12, 23,24,25,26,27,28,29,30,31,32,33,34].

In [30], Kullback-Leibler distance was proposed for tracking the dynamics of measurement variations to detect FDI attacks. The statistical behavior of the state estimation process, through cumulative sum based approach (CUSUM), was used for FDI attacks detection in [31]. Method based on the short-term state forecasting which considers temporal correlations between nodal states was proposed in [32]. Detection scheme using two physical system parameters and their behavior was proposed in [33]. In [34] to detect FDI attacks, a cosine similarity matching technique was used and tested in power systems for which estimated (expected) measurement values were obtained using a Kalman filter, by comparing them with actual measurements. The underlying SE observation model was regarded as linear, as was the state transition matrix (approximated by the identity matrix), enabling a Kalman filter to be used for SE.

However, in most commercial SE algorithms, the utilized measurement observation model is highly nonlinear, especially regarding the line current and transformer tap position measurements, implying the possible use of Kalman filter extensions intended for nonlinear systems - extended or unscented Kalman filters (EKF or UKF, respectively). Kalman filtering techniques were first proposed for use in dynamic state estimation in [35], and significant efforts have been subsequently undertaken to improve their performance [36,37,38,39,40,41]. In an EKF, which is probably the most widely used dynamic SE algorithm, the state distribution is approximated by Gaussian random variables (GRVs), which are propagated through a “first order” linearization of the nonlinear system, and EKF can therefore be seen as “first order” approximation of the optimal filter [36]. Enhancements to the EKF were proposed in [37, 38] to improve its performance and robustness by incorporating the high nonlinearity of measurement functions. Furthermore, an iterated EKF based on the generalized maximum likelihood approach was proposed in [39] to estimate dynamic states during disturbances. In [40], a UKF was introduced as a tool for dynamic state estimation, through a combination an unscented transformation and Kalman filter. In the UKF technique, the minimal set of selected sample points, which capture the true mean and covariance of the GRVs, is transferred through the complete nonlinear system, obtaining the posterior mean and covariance from the “third order” approximation of the optimal filter [41].

In this paper, we investigated the detection of FDI attacks by using a UKF to predict and update SVs starting from the previously known state and compared them with the results acquired from a typically used WLS-based SE algorithm. It should be noted that the UKF used in the proposed algorithm can be replaced by other nonlinear filters without any structural changes. We show that the SVs under attack significantly deviated between the UKF and WLS-based SE. The measurements influencing the suspicious SVs are those that may be under malicious attack.

The main contributions of the proposed method are as follows: ① a derivation of the time-variant transition function (necessary for the UKF prediction step) by a combination of power flow equations with load/generation very short-term forecasts and generator schedules; ② false data detection using the normalized SV residuals obtained from WLS-based SE and UKF estimates, as well as the UKF state covariance matrix; ③ an analysis of the most critical scenario, where the attacker may gain access to the complete network model and set of measurements; ④ the efficient synergy of WLS-based SE and UKF algorithms for the real-time (online) detection of FDI attacks.

The paper is organized as follows. In Section 2, we formulate the problem of detecting FDI attacks, in which the advisory agent has an unlimited number of resources at its disposal. In Section 3, the proposed UKF based detection algorithm is explained, and the BDD process using UKF and SE results is covered in Section 4. The results of numerical simulations for benchmark power systems (with 14 and 300 buses), empowered with the data needed for forecast analysis, are provided in Section 5, and concluding remarks are presented in Section 6. The algorithmic details of the UKF are provided in Appendix A.

2 Problem formulation

We consider a scenario in which an attacker has gained access to an unlimited number of measurements by breaking into a SCADA system (Case 2 in Fig. 1), primary domain controller (PDC) (Case 5) and control center, or by directly tampering with the PMUs (Case 4), RTUs (Case 3) and communications network (Case 1).

Fig. 1

FDI attack detection by comparing WLS-based SE and UKF results

With the idea of analyzing the most critical scenario, which may not be the most typical and frequent in practice, we suppose that the attacker has sufficient knowledge of the power system (topology and physics) as described by a nonlinear SE measurements-SVs model

$$\varvec{z} = \varvec{h}(\varvec{x}) + \varvec{e}$$

(1)

or by a linear SE model (assumed only for analysis in this section)

$$\varvec{z} = \varvec{Hx} + \varvec{e}$$

(2)

where z and x are the M-dimensional measurement ($z_{j} \in {\varvec{R}}^{M \times 1}$) and N-dimensional state vectors ($x_{i} \in {\varvec{R}}^{N \times 1}$), M and N are the total number of measurements and SVs; $\varvec{h}(\varvec{x})$ is the vector of nonlinear functions of the system state vector x – observation model function; e is the measurement error vector (assumed to be zero mean multivariate Gaussian noise with covariance R); H is the Jacobian matrix and $\varvec{H} = \frac{{{\partial }\varvec{h}(\varvec{x})}}{{\partial \varvec{x}}}$.

2.1 Analysis of compromised measurement cases

Typically, the goal of an attacker is to change either the SV(s) ($x_{i}$) to the desired value by $\delta x_{i}$ or to change some particular measurement (probably power flow) to an arbitrary scalar value. To change specific measurements and remain undetected by the conventional bad data detection algorithms, the adversary would need to modify the corresponding SV (the one affecting the desired power flow measurement, and a minimum number of other measurements). Therefore, without loss of generality, we suppose that SV will be changed as a goal of the attack.

When we assume a realistic scenario where an attacker cannot obtain access to every measurement in the system (some can be regarded as secure), the measurements can be divided into two groups: ①fully protected measurements ($1,\;2,\; \cdots \;,\;m$), denoted by ‘*’; ② unprotected measurements ($1,\;2,\; \cdots \;,\;M-m$).

The measurement vector z can then be written as:

$$\varvec{z} = \left[ {\begin{array}{*{20}c} {z_{1}^{*} } \\ \vdots \\ {z_{m}^{*} } \\ {z_{1} } \\ \vdots \\ {z_{{M-m}} } \\ \end{array} } \right]$$

(3)

Likewise, the state vector can be divided into the protected SVs that correspond to at least one of the secure measurements ($x_{k}^{*}$) and the unprotected ones ($x_{N - k}$ is the total k SV which is fully protected – subvector $x_{k}^{*}$ in (4)), where all the corresponding measurements may be compromised.

If the attacker plans to change the i^th SV ($x_{i}$) by $\delta x_{i}$ while remaining undetected, then all the measurements corresponding to the rows with non-zero elements (total $\ell$ measurements) in the ith column of matrix H need to be changed ($z_{\ell } + \delta z_{\ell }$). The attacker can change only the unprotected SVs ($x_{N - k}$) without detection because all the corresponding measurements are those that can be tampered with $z_{{M - m}}^{{}}$ (the measurement error vector e in (2) is ignored).

$$\left[ {\begin{array}{*{20}c} {z_{1}^{*} } \\ \vdots \\ {z_{m}^{*} } \\ {z_{1} } \\ \vdots \\ z_{\ell } + \delta z_{\ell } \\ \vdots \\ \\ {z_{M - m} } \\ \end{array} } \right] = \left[ {\begin{array}{*{20}c} \cdots & {H_{1i} } & \cdots \\ {} & \vdots & {} \\ \cdots & {H_{ii} } & \cdots \\ {} & \vdots & {} \\ \cdots & {H_{Mi} } & \cdots \\ \end{array} } \right]\left[ {\begin{array}{*{20}c} {x_{1}^{*} } \\ \vdots \\ x_{k}^{*} \\ x_{1}^{{}} \\ \vdots \\ x_{i} + \delta x_{i} \\ \vdots \\ {x_{N - k} } \\ \end{array} } \right]$$

(4)

The Jacobian matrix used by an adversary may be represented as a sum of the actual Jacobian H and some increment δH that is used to encapsulate all the uncertainties the attacker may have. The measurement-SVs model with multiple SVs under attack ($\varvec{x}_{N - k}$) can now be written as:

$$\begin{aligned} \left[ {\begin{array}{*{20}c} {\varvec{z}_{m}^{*} } \\ {\varvec{z}_{M - m} } \\ \end{array} } \right] & = \left[ {\begin{array}{*{20}c} {\varvec{H}_{m,k}^{{}} } & {\varvec{H}_{m,N - k}^{{}} } & {\mathbf{0}} \\ {\varvec{H}_{M - m,k}^{{}} } & {\varvec{H}_{M - m,N - k}^{{}} } & {\varvec{H}_{M - m,N - k}^{{}} + \delta \varvec{H}_{M - m,N - k} } \\ \end{array} } \right]\left[ {\begin{array}{*{20}c} {\varvec{x}_{k}^{*} } \\ {\varvec{x}_{N - k}^{{}} } \\ {\delta \varvec{x}_{N - k} } \\ \end{array} } \right] \\ & = \left[ {\begin{array}{*{20}c} {\varvec{H}_{m,k}^{{}} } & {\varvec{H}_{m,N - k}^{{}} } \\ {\varvec{H}_{M - m,k}^{{}} } & {\varvec{H}_{M - m,N - k}^{{}} } \\ \end{array} } \right]\left[ {\begin{array}{*{20}c} {\varvec{x}_{k}^{*} } \\ {\varvec{x}_{N - k}^{{}} } \\ \end{array} } \right] + \left[ {\begin{array}{*{20}c} {\mathbf{0}} \\ {(\varvec{H}_{M - m,N - k}^{{}} + \delta \varvec{H}_{M - m,N - k} )\delta \varvec{x}_{N - k} } \\ \end{array} } \right] \\ & = \left[ {\begin{array}{*{20}c} {\varvec{H}_{m,k}^{{}} } & {\varvec{H}_{m,N - k}^{{}} } \\ {\varvec{H}_{M - m,k}^{{}} } & {\varvec{H}_{M - m,N - k}^{{}} } \\ \end{array} } \right]\left[ {\begin{array}{*{20}c} {\varvec{x}_{k}^{*} } \\ {\varvec{x}_{N - k}^{{}} } \\ \end{array} } \right] + \left[ {\begin{array}{*{20}c} {\mathbf{0}} & {\mathbf{0}} \\ {\varvec{H}_{M - m,N - k}^{{}} } & {\mathbf{1}} \\ \end{array} } \right]\left[ {\begin{array}{*{20}c} {\delta \varvec{x}_{N - k} } \\ {\delta \varvec{H}_{M - m,N - k} \delta \varvec{x}_{N - k} } \\ \end{array} } \right] \\ & = \left[ {\begin{array}{*{20}c} {\varvec{H}_{m,k}^{{}} } & {\varvec{H}_{m,N - k}^{{}} } & {\mathbf{0}} & {\mathbf{0}} \\ {\varvec{H}_{M - m,k}^{{}} } & {\varvec{H}_{M - m,N - k}^{{}} } & {\varvec{H}_{M - m,N - k}^{{}} } & {\mathbf{1}} \\ \end{array} } \right]\left[ {\begin{array}{*{20}c} {\varvec{x}_{k}^{*} } \\ {\varvec{x}_{N - k}^{{}} } \\ {\delta \varvec{x}_{N - k} } \\ {\delta \varvec{z}_{M - m} } \\ \end{array} } \right] \\ \end{aligned}$$

(5)

where $\delta \, \varvec{z}_{{M{ - }m}} = \delta \varvec{H}_{{M{ - }m,\,N{ - }k}} \delta \varvec{x}_{{N{ - }k}}$ and ${\mathbf{1}}$ is an identity matrix.

Note that in (5), the state vector has been expanded with the additional sets of variables ($\delta \varvec{x}_{N - k}$ and $\delta \varvec{z}_{M - m}$) that, represent malicious SV changes and measurement changes (resulting as a consequence of the attacker’s incomplete knowledge of Jacobian matrix H), respectively.

Unfortunately, the resulting set of (5) is generally undetermined and cannot be solved even for the simplest scenario with the maximum measurement redundancy ($M > 3N$) and only one SV under attack ($x_{i}$) – irrespective of the size of $\varvec{x}_{k}^{*}$. Therefore, a WLS-based SE cannot by itself be enhanced to detect FDI attacks but only in synergy with some other method. In consequence, a new algorithm that combines a WLS-based SE and forecast based UKF is proposed in Section 3.

2.2 Nonlinear problem formulation

Based on the analysis provided in Section 2.1, a critical case wherein all the measurements are susceptible to attacks ($m\text{ = 0}$ and $k\text{ = 0}$ in (4)) is considered. Let $\delta \varvec{z}$ be the non-zero false data attack vector from (4) that is created by the adversary, who has a complete knowledge of the nonlinear observation model function $\varvec{h}(\varvec{x})$. The modified set of measurements can then be written as:

$$\varvec{z^{\prime}} = \varvec{z} + \delta \varvec{z} + \varvec{e} = \varvec{h}(\varvec{x} + \delta \varvec{x}) + \varvec{e} = \varvec{h}(\varvec{x^{\prime}}) + \varvec{e}$$

(6)

The typically used BDD techniques are based on a normalized residual approach, where the $\ell^{2}$-norm of the measurement residual is compared against a threshold τ.

$$\left\| {\varvec{z^{\prime}} - \varvec{h}(\varvec{x^{\prime}})} \right\| = \left\| {\varvec{h}(\varvec{x} + \delta \varvec{x}) + \varvec{e} - \varvec{h}(\varvec{x} + \delta \varvec{x})} \right\| = \left\| \varvec{e} \right\| < \tau$$

(7)

Therefore, malicious attacks created in this manner cannot be detected using standard BDD techniques [6].

3 Proposed algorithm

A flow-chart of the proposed algorithm for the detection of FDI attacks is shown in Fig. 2, wherein the different steps are described in sequence.

Fig. 2

Proposed algorithm for detection of FDI attack

3.1 UKF-based states prediction and update

The prediction of the SVs from the (k – 1)^th time instant to the those at the k^th time instant

$$\varvec{x}_{k|k - 1} = \varvec{f}(\varvec{x}_{k - 1} ) + \varvec{w}_{k}$$

(8)

is performed by the UKF prediction step (described briefly in Appendix A), where the linearized transition model can be derived from increments to the power flow equations:

$$\left[ {\begin{array}{*{20}c} {\Delta \varvec{P}} \\ {\Delta \varvec{Q}} \\ \end{array} } \right] = \left[ {\begin{array}{*{20}c} {\Delta \varvec{P}_{g} - \Delta \varvec{P}_{p} } \\ {\Delta \varvec{Q}_{g} - \Delta \varvec{Q}_{p} } \\ \end{array} } \right] = \left[ {\begin{array}{*{20}c} \varvec{E} & \varvec{T} \\ \varvec{K} & \varvec{L} \\ \end{array} } \right]\left[ {\begin{array}{*{20}c} {\Delta\varvec{\theta}} \\ {\Delta \varvec{V}} \\ \end{array} } \right] = \varvec{J}\Delta \varvec{x}$$

(9)

where E, T, K, L are corresponding submatrices of Jacobian matrix J for the power flow equations, with elements $E_{ij} = \left. {\frac{{\partial P_{i} }}{{\partial \theta_{j} }}} \right|_{{{\theta}_{0} ,{V}_{0} }}$, $T_{ij} = \left. {\frac{{\partial P_{i} }}{{\partial V_{j} }}} \right|_{{{\theta}_{0} ,{V}_{0} }}$, $K_{ij} = \left. {\frac{{\partial Q_{i} }}{{\partial \theta_{j} }}} \right|_{{{\theta}_{0} ,{V}_{0} }}$ and $L_{ij} = \left. {\frac{{\partial Q_{i} }}{{\partial V_{j} }}} \right|_{{{\theta}_{0} ,{V}_{0} }}$. Note that the only equation for active and reactive power in the Slack-bus is excluded (i, j ≠ Slack), because $\theta_{Slack} = 0$.

For the assumed constant power factor ($\cos \varphi$), we have:

$$\left[ {\begin{array}{*{20}c} {\Delta \varvec{P}} \\ {\Delta \varvec{Q}} \\ \end{array} } \right] = \left[ {\begin{array}{*{20}c} {\Delta \varvec{P}} \\ {\text{tan }\varphi \Delta \varvec{P}} \\ \end{array} } \right]$$

(10)

where

$$\Delta P_{i} = s_{i} P_{k - 1,i}$$

(11)

and $P_{k - 1,i}$ represents the per unit injection change (referenced to the base case values (1 p.u.) at time instant $t_{k - 1}$) taken from the active power generation/load curve, whereas $s_{i}$ represents the correlation coefficient between $\Delta P_{i}$ and $P_{k - 1,i}$, which is also obtained from the generation/load curve. These curves are acquired as a result of the application of short-term forecasts and generator schedules.

Finally, we obtain the following time-variant transition model for UKF

$$\varvec{x}_{k} = \varvec{x}_{k - 1} + \Delta \varvec{x} = \varvec{x}_{k - 1} + \varvec{J}_{k - 1}^{{\text{inv}}} \left[ {\begin{array}{*{20}c} {\Delta \varvec{P}} \\ {\text{tan }\varphi \Delta \varvec{P}} \\ \end{array} } \right]$$

(12)

or

$$\varvec{x}_{k} = \varvec{F}_{k - 1} \varvec{x}_{k - 1}$$

(13)

where $\varvec{F}_{k - 1} = {\mathbf{1}} + (\varvec{J}_{k - 1}^{{i\text{nv}}} )^{\prime}$is the linearized UKF transition function. The elements of the $(\varvec{J}_{k - 1}^{{i\text{nv}}} )^{\prime}$ matrix are defined as:

$$(\varvec{J}_{k - 1,ij}^{{{\text{inv}} }} ) ^{\prime }= J_{k - 1,ij}^{\text{inv}} \frac{{s_{i} P_{k - 1,i} }}{{x_{k - 1,i} }}$$

(14)

The derived transition function is generally a full matrix, because changes in load and generation buses influence multiple SVs. Being time-varying, it captures sudden load/generation changes in power systems to a much higher degree than constant transition functions.

The update step of the UKF is then performed using a highly non-linear observation model $\varvec{h}(\varvec{x})$ acquired with the WLS-based SE algorithm, as explained in Appendix A.

3.2 PMU measurement buffering

The state estimation calculation is usually triggered either periodically at periods of 1-5 minutes or after a topological or significant analogue measurement change. Depending on the SCADA pooling time and measurement dynamics, the expected SE triggering period is usually 30 seconds to 1 minute – which is much longer than the refresh rates of PMU-based measurements. Due to these differences, measurement buffering was proposed in [42] For $N^{{\text{PMU}}}$ measurements that arrive in the time interval between two WLS-based SE executions ($t_{k - 1} \le t_{i} \le t_{k - 1} + \Delta t^{{\text{SE}}}$) the mean and variance respectively, are

$$\mu_{z} = \frac{1}{{N^{{\text{PMU}}} }}\sum\limits_{i = 1}^{{N^{{\text{PMU}}} }} {z(t_{i} )}$$

(15)

$$\text{var} (\mu_{z} ) = \frac{1}{{\text{(}N^{{\text{PMU}}} \text{)}^{2} }}\sum\limits_{i = 1}^{{N^{{\text{PMU}}} }} {\text{var} \{ z(t_{i} )\} } = \frac{{\sigma^{2} }}{{N^{{\text{PMU}}} }}$$

(16)

3.3 WLS-based SE

Because the SE model is nonlinear, the function $\varvec{h}(\varvec{x})$ is linearized and the resulting linearized SE is formulated as a WLS problem. The WLS estimator minimizes the objective function by satisfying the first-order optimality conditions.

$$\frac{{\partial \varvec{j}(\varvec{x})}}{{\partial \varvec{x}}} = 0$$

(17)

where $\varvec{j}(\varvec{x}) = \Delta \varvec{z}^{\text{T}} \varvec{R}^{ - 1} \Delta \varvec{z}$ is the objective function; $\Delta \varvec{z} = \varvec{z} - \varvec{h}(\varvec{x})$ is the vector of measurement residuals; $\varvec{R}$ is the covariance matrix of measurement error vector.

The minimization problem can be solved using an iterative scheme, changing the point of linearization (for simplicity, the iteration index is removed).

$$(\varvec{H}^{\text{T}} \varvec{R}^{ - 1} \varvec{H})\Delta \hat{\varvec{x}} = \varvec{H}^{\text{T}} \varvec{R}^{ - 1} \Delta \varvec{z}$$

(18)

$$\varvec{x} = \varvec{x} + \Delta \hat{\varvec{x}}$$

(19)

$$\hat{\varvec{z}} = \varvec{h}(\varvec{x} + \Delta \hat{\varvec{x}})$$

(20)

where the Jacobian matrix H is assumed to be constant inside one SE cycle.

4 FDI attack detection using UKF

Knowing the initial network condition at time instant k-1 (vector of SVs $\varvec{x}_{k - 1}$) and assuming that the forecasted data are available for all the power consumers and generators for the following time instant k, a UKF can be used to determine the network condition on the basis of the transition matrix and available telemetered measurements. A transition matrix ($\varvec{F}_{k - 1}$) can be created using forecasted data, as part of load flow calculation for future time instant k using (10)-(14).

In a commercial EMS environment, an SE calculation is usually triggered in four ways: ① periodically, on user defined time intervals that are primarily in the range of 1–5 minutes; ② after significant analogue measurement(s) change(s); ③ after every topological or change in transformer tap position; ④ on the request of a user. Therefore, it can be expected that the period between time instants k – 1 and k is no longer than 5 minutes. For such short time periods, it can be assumed that the trust factor for the forecasted injections (generations minus loads) is high, resulting in small values of the process noise vector $\varvec{w}_{k}$ in (A1).

When an FDI attack occurs, the maliciously changed group of measurements would swing the estimates of the SV(s) provided by the WLS-based SE algorithm in the desired direction. On the other hand, the estimates acquired from the UKF would only partially swing, due to the binding effect of the transition matrix ($\varvec{F}_{k - 1}$) and the small process noise vector ($\varvec{w}_{k}$).

For the purpose of detecting FDI attacks and identifying false data, a normalized SV residuals based approach is proposed. A normalized SV residual ($r_{i}$) is calculated as the absolute value of a difference between a SV value estimated by the WLS-based SE ($x_{i}^{{\text{SE}}}$) and that calculated by the UKF ($x_{i}^{{\text{UKF}}}$) divided by the standard deviation acquired from the UKF covariance matrix $\varvec{C}$ in (A4).

$$r_{i} = \frac{{\left| {x_{i}^{{\text{SE}}} - x_{i}^{{\text{UKF}}} } \right|}}{{\sqrt {C_{i,i} } }}\quad \;i = 1,\;2, \cdots ,N$$

(21)

It should be noted that the traditional BDD algorithms, which are typically based on the normalized measurement residuals, are used to filter out bad measurements in both the WLS-based SE and UKF in the proposed algorithm. These filtered estimates are later used to calculate the residuals (21).

To determine whether the regarded SVs are under attack, the normalized SV residuals ($r_{i}$ in (21)) are compared to a predefined threshold. When an FDI attack is confirmed, the proposed algorithm generates an alarm in the control room to warn a system operator.

Because the proposed method relies on forecast results, false positive detections may occur. Note that false positives are the cases where the proposed algorithm detects an attack although there is none, and false negatives are the cases when the algorithm fails to detect an existing attack. Detection depends on the specified threshold, which should provide the minimum number of false positives and false negatives. Threshold selection is very important because higher threshold values would tolerate larger forecast errors, while at the same time limiting the minimum attack intensity that could be detected. Lower threshold values would, on the other hand, increase the probability of false positive detections. Being system oriented and forecast quality dependent, the following steps are proposed for threshold definition:

• 1) Generation of daily-based load curves obtained from the application of very short-term load forecasts for multiple day types (week day, weekend, and holiday) and seasons (spring, summer, autumn, and winter).

• 2) Iterative adjustment (decrease/increase) of the threshold value, starting from the initial one.

• 3) Application of FDI attacks of different intensities.

• 4) Attack detection using the proposed algorithm.

Only sudden changes in generation/load, that are not captured in a very short-term forecast as well as unplanned system events such as equipment outages (for example, overhead line outages), could produce differences between UKF and WLS-based SE estimates that are not the consequence of a malicious attack and would therefore lead to false positive detections. However, such detections could easily be disregarded after the system operator becomes aware of them. Orchestrating FDI attacks during such unplanned events to evade detections is regarded as highly unlikely. An attacker would need to have access to all the measurements in the region over a long period and then attack one of the SVs affected by the outage during a short time frame (inside a 5 minute’ interval, before the next execution of a very short-term forecast).

5 Application

To evaluate the proposed algorithm, simulations were conducted on the modified 14-bus test system [43] shown on Fig. 3, and on an IEEE 300-bus test system [44]. The SVs calculated using the WLS-based SE and UKF were compared in cases with and without a malicious attack on the smaller 14-bus test system. On the same system, the sensitivity of the algorithm to different transition process noise values was tested. The ability of the proposed method to detect attacks of different intensities and its sensitivity to forecast accuracy was demonstrated on the 300-bus test case.

Fig. 3

Single-line diagram of modified 14-bus test system

5.1 Modified 14-bus test system

A single-line diagram of the modified 14-bus test system is shown in Fig. 3. Compared to the original model in [43], generation was removed from buses 3 and 6 and added to buses 5, 11 and 13. The types of all the generating units and utilized daily load curves are specified in Fig. 3. The initial voltage phasors and nodal active and reactive power injections were also slightly changed from their original values to compensate for the addition/removal of generating units shown in Table 1.

Table 1 Active and reactive power injections into system buses

High measurement redundancy was assumed (with active and reactive power flow measurements on all branches, shunt power injection measurements, and voltages at every bus). PMU-based voltage measurements were placed on the generator (generation, PV and slack) buses. Forecasted generation curves for solar and wind units, and schedules for thermal and hydro units were used in Fig. 4 in conjunction with forecasted load curves in Fig. 5 to estimate near-term network condition.

Fig. 4

Forecasted generation curves

Fig. 5

Forecasted daily load curves

In Fig. 6, the WLS-based SE and UKF estimated voltage magnitudes at bus 12 (load, PQ bus) were plotted for the analyzed 6-hour period when no attack was introduced. Estimation calculations were conducted every 30 seconds, which is an expected average SE sampling period for utilities with modern EMSs. The maximum normalized residuals for all the voltage magnitudes are shown on the figure as well. It can be easily concluded that as the UKF estimates followed the WLS-based SE estimates, the normalized residuals were below the predefined threshold, and no FDI attack was detected.

Fig. 6

Obtained results for modified 14-bus test system without FDI attack

An FDI attack aimed at the voltage at bus 12 was started after 1 hour, as shown in Fig. 7. The intensity of the voltage magnitude change targeted by the attacker was 0.1 p.u. It can be observed that although the WLS-based SE calculated voltage immediately reached the targeted voltage, the voltage calculated by the UKF only reached the targeted value after more than one hour, thus enabling the FDI attack to be detected, which can be seen by analysing the maximum normalized residuals over time.

Fig. 7

Obtained results for modified 14-bus test system with an FDI attack starting after one hour

The influence of the process noise vector, $\varvec{w}_{k}$ in (A1), was also analysed, and was determined to be critical to the FDI detection process. In Fig. 8, the UKF estimated voltages for three different process noise values were plotted; the smallest process noise was 10 times smaller than the medium process noise and 100 times smaller than the largest process noise.

Fig. 8

UKF estimated voltages at bus 12 using three different process noise values

From the presented results, it was concluded that for large process noise, the UKF results closely followed the WLS-based SE estimates, thus reducing the time available for FDI attack detection, whereas for lower process noise, the differences between the WLS-based SE and UKF results were much more significant and remained so for longer periods of time. This conclusion emphasizes the necessity for high-quality generation and load forecasts, and for short SE execution periods.

5.2 IEEE 300-bus test system

Furthermore, FDI attacks of different intensities were aimed at bus 4 (load, PQ type), which affected all the necessary measurements to change the SE results without being detected by the typically used BDD algorithms. The voltages estimated using the WLS-based SE method in the first execution after the attack commenced were compared with those acquired from the UKF, as shown in Table 2. It was concluded that for attacks of greater strength (targeted malicious voltage magnitude changes in p.u.), the differences between the WLS-based SE and UKF estimates were significant, as was the normalized residual, which overtopped the defined threshold by almost 20 times. The maximum normalized residuals decreased with decreasing attack intensity - the 0.01 p.u. attack was the last to be successfully detected (last row in Table 2). Depending on the accuracy of measurement equipment in the power system, voltage magnitudes may be estimated using the WLS-based SE algorithm alone, with errors exceeding 1% (even in the absence of a malicious attack). The consequences of such errors are minor, which leads to the conclusion that the proposed method is instantly capable of detecting any FDI attack that may create significant consequences to the network operation.

Table 2 Detection of FDI attacks of different strengths

In addition, the sensitivity of the proposed method was tested against different very short-term forecast accuracies. The expected mean absolute percentage error (MAPE) for very short-term forecast (5-minute in advance) in the modern EMS systems is significantly lower than 1%. Therefore, the maximum normalized SV residuals in the absence of an FDI attack were obtained for forecast accuracies ranging from 0% to 1% MAPE, as shown in Table 3. The minimum strengths of the FDI attacks that were successfully detected without triggering any false positive detections (see in Table 3) expectedly increased with the increase of forecast MAPE. However, even for the higher forecast errors the proposed method proved capable of detecting FDI attacks of significant strengths. The range of the FDI attack strengths that could be detected would further increase with the less strict selection of the detection threshold (allowing false positive detections).

Table 3 Detection of FDI attacks for different forecast accuracies

One potential drawback of the proposed approach occurs in cases when the malicious attack is comprised of multiple low strength attacks aimed at the same SV, driving its value towards the targeted value in small increments over a substantial period of time. For attacks constructed in that manner, the differences between WLS-based SE and UKF calculated values could be below the threshold in the time instant of the incremental attack. In addition, the UKF estimation would reach the WLS-based SE after a limited number of iterations but before the following incremental attack - allowing the attack to be continuously undetected, as shown in Fig. 9.

Fig. 9

Obtained results for IEEE 300 bus test system under an incremental FDI attack over a 24-hour time period

On the other hand, constructing this type of FDI attack is much more difficult because a prolonged period of time is needed for the attack to attain the desired strength necessary to produce any significant consequences. Furthermore, multiple meaningful changes to all the affected measurements are needed. Such practical complications make this type of an attack an unlikely scenario. Nevertheless, our subsequent research will focus on these slow, incremental FDI attacks.

6 Conclusion

In this paper, the synergy between a traditional WLS-based SE and UKF was used for detection of FDI attacks in power systems. An unlimited number of compromised measurements was supposed, and a more realistic nonlinear observation model was used. To implement the UKF prediction step, a transition function was derived using a power flow model in combination with load/generation short-term forecasts and generator schedules.

The proposed method was tested on benchmark IEEE 14-bus and 300-bus test systems, and all FDI attacks of significant strength were successfully identified for reasonable forecast quality. The only possible failure cases which have been identified are as follows: ① detection of lower strength FDI attacks in systems with bad forecast quality; ② false positive detections in cases of sudden changes in generation/load, that are not captured within a very short-term forecast and unplanned system events (such as equipment outages); ③ detection of attacks consisting of multiple small intensity attacks executed over a long time period. While the third case is highly unlikely and impractical, the second one could easily be disregarded by the system operator.

By relying on a more realistic and accurate nonlinear measurement observation model, and a fast, easily implemented and robust UKF, the proposed solution can be successfully deployed as part of an SE package in commercial EMS software.

Appendix A

The transition of a system state vector from the (k-1)^th to k^th time instant can be described as:

$$\varvec{x}_{k} = \varvec{f}(\varvec{x}_{k - 1} ) + \varvec{w}_{k}$$

(A1)

where $\varvec{x}_{k}$ is the N-dimensional state estimation vector at k^th time instant; $\varvec{f}( \cdot )$ is the N-dimensional nonlinear set of state transition functions; $\varvec{w}_{k}$ is the N-dimensional process noise vector, assumed to be zero mean multivariate Gaussian noise with covariance $\varvec{Q}_{k}$, at k^th time instant; N is the number of unknown system SVs (the same as in (1)).

In this paper, a linearized form of the transition model function $\varvec{f}( \cdot )$ is proposed (F in (13)).

The observation model function (in our case, the SE model) is represented by the following equation (which is the same as (1) for the k^th time instant).

$$\varvec{z}_{k} = \varvec{h}(\varvec{x}_{k} ) + \varvec{e}_{k}$$

(A2)

The UKF uses the unscented transformation sampling technique to perform a nonlinear transformation. A minimal set of sampling points (known as sigma points) is selected and transformed using the nonlinear function, whereas the new mean and covariance are formed out of those transformed points. The basic underlying idea is that it is easier to approximate a Gaussian distribution, than it is to approximate an arbitrary nonlinear function.

The state estimation vector at the (k - 1)th time instant ($\varvec{x}_{k - 1}$) and its covariance matrix ($\varvec{C}_{k - 1}$) are augmented with the mean (0) and covariance matrix $\varvec{Q}_{k}$ of the process noise vector $\varvec{w}_{k}$, respectively.

$$\varvec{x}_{k - 1}^{a} = \left[ {\varvec{x}_{k - 1}^{\text{T}}\, {\mathbf{0}}} \right]^{\text{T}}$$

(A3)

$$\varvec{C}_{k - 1}^{a} = \left[ {\begin{array}{*{20}c} {\varvec{C}_{k - 1} } & {\mathbf{0}} \\ {\mathbf{0}} & {\varvec{Q}_{k} } \\ \end{array} } \right]$$

(A4)

$$n^{a} = 2N$$

(A5)

where $n^{a}$ represents the dimension of the augmented state vector ($\varvec{x}_{k - 1}^{a}$), and $\varvec{C}_{k - 1}^{a}$ represents the augmented state covariance matrix.

The $2n^{a} + 1$ sigma points ($\varvec{\chi}_{k - 1|k - 1}^{{}}$) and weight are, respectively, defined as:

$$\varvec{\chi}_{k - 1|k - 1}^{i} = \left\{ \begin{array} {ll}\varvec{x}_{k - 1}^{a} &\quad i = 0 \hfill \\ \varvec{x}_{k - 1}^{a} + \varvec{B}_{i} &\quad i = 1, 2,\ldots , n^{a} \hfill \\ \varvec{x}_{k - 1}^{a} - \varvec{B}_{{i - n^{a} }} &\quad i = n^{a} + 1, n^{a} +2,\ldots , 2n^{a} \hfill \\ \end{array} \right.$$

(A6)

$$W_{i}^{s} = \left\{ \begin{array}{ll} \frac{\kappa }{{n^{a} + \kappa }}&\quad i = 0 \hfill \\ \frac{1}{{2\left( {n^{a} + \kappa } \right)}}&\quad i = 1, 2,\ldots , 2n^{a} \hfill \\ \end{array} \right.$$

(A7)

$$W_{i}^{c} = \left\{ \begin{array}{ll} \frac{\kappa }{{n^{a} + \kappa }} + (1 - \alpha^{2} + \beta )&\quad i = 0 \hfill \\ \frac{1}{{2(n^{a} + \kappa )}}&\quad i = 1,2, \ldots , 2n^{a} \hfill \\ \end{array} \right.$$

(A8)

where $\varvec{B}_{i}$ is the ith column of the matrix and $\varvec{B} = \sqrt {\left( {n^{a} + \kappa } \right)\varvec{C}_{k - 1}^{a} }$; $\alpha$ is the constant that determines the spread of the sigma points around the $\varvec{x}_{k - 1}^{a}$ (usually set to a small number, e.g., 0.001); $\beta$ is the constant used to incorporate the prior knowledge of the distribution of state vector ($\beta = 2$ for Gaussian distributions); $\kappa$ is the scaling parameter, usually calculated as $\kappa = n^{a} (\alpha^{2} - 1)$; $W_{i}^{s} \;(W_{i}^{c} )$ is the weight associated with the ith sigma point and used to calculate the predicted state vector (covariance matrix).

The predicted sigma points (${\varvec{\upchi}}_{k|k - 1}^{i}$) are acquired by propagating the previously calculated sigma points using the state transition function $\varvec{f}( \cdot )$, to obtain the predicted state vector with mean ($\varvec{x}_{k|k - 1}$) and covariance ($\varvec{C}_{k|k - 1}$), respectively

$$\varvec{\chi}_{k|k - 1}^{i} = \varvec{f}(\varvec{\chi}_{k - 1|k - 1}^{i} )$$

(A9)

$$\varvec{x}_{k|k - 1} = \mathop \sum \limits_{i = 0}^{{2n^{a} }} W_{i}^{s}\varvec{\chi}_{k|k - 1}^{i}$$

(A10)

$$\varvec{C}_{k|k - 1} = \mathop \sum \limits_{i = 0}^{{2n^{a} }} W_{i}^{c} \left( {\varvec{\chi}_{k|k - 1}^{i} - \varvec{x}_{k|k - 1} } \right)\left( {\varvec{\chi}_{k|k - 1}^{i} - \varvec{x}_{k|k - 1} } \right)^{\text{T}}$$

(A11)

Afterwards, the obtained predicted state vector and state covariance matrix should be augmented with the mean (0) and covariance matrix ($\varvec{R}_{k}$) of the observation noise vector, respectively.

$$\varvec{x}_{k|k - 1}^{a} = \left[ {\varvec{x}_{k|k - 1}^{\text{T}}\, {\mathbf{0}}} \right]^{\text{T}}$$

(A12)

$$\varvec{C}_{k|k - 1}^{a} = \left[ {\begin{array}{*{20}c} {\varvec{C}_{k|k - 1} } & {\mathbf{0}} \\ {\mathbf{0}} & {\varvec{R}_{k} } \\ \end{array} } \right]$$

(A13)

$$n^{a} = 2N + N = 3N$$

(A14)

After a new set of sigma points ($\varvec{\chi}_{k|k - 1}^{i}$) are defined using (A6)–(A8), they are propagated through the observation model (h) to determine the predicted measurement vector and the predicted measurement covariance matrix, respectively.

$$\hat{\varvec{z}}_{k} = \mathop \sum \limits_{i = 0}^{{2n^{a} }} W_{i}^{s} \bar{\varvec{z}}_{k}^{i}$$

(A15)

$$\varvec{C}_{{\varvec{z}_{k} \text{,}\varvec{z}_{k} }} = \mathop \sum \limits_{i = 0}^{{2n^{a} }} W_{i}^{c} \left( {\bar{\varvec{z}}_{k}^{i} - \hat{\varvec{z}}_{k} } \right)\left( {\bar{\varvec{z}}_{k}^{i} - \hat{\varvec{z}}_{k} } \right)^{\text{T}}$$

(A16)

where

$$\bar{\varvec{z}}_{k}^{i} = \varvec{h}\text{(}\varvec{\chi}_{k|k - 1}^{i} \text{)}$$

(A17)

The state-measurement cross-covariance matrix is determined as:

$$\varvec{C}_{{\varvec{x}_{k} \text{,}\varvec{z}_{k} }} = \mathop \sum \limits_{i = 0}^{{2n^{a} }} W_{i}^{c} \left( {{\varvec{\upchi}}_{k|k - 1}^{i} - \varvec{x}_{k|k - 1} } \right)\left( {\bar{\varvec{z}}_{k}^{i} - \hat{\varvec{z}}_{k} } \right)^{\text{T}}$$

(A18)

and is used to calculate the Kalman gain

$$\varvec{K}_{k} = \varvec{C}_{{\varvec{x}_{k} \text{,}\varvec{z}_{k} }} \varvec{C}_{{\varvec{x}_{k} \text{,}\varvec{z}_{k} }}^{ - 1}$$

(A19)

The state vector is updated by adding the innovation (${\mathbf{z}}_{k} - {\hat{\mathbf{z}}}_{k}$) weighted by the Kalman gain to the predicted state

$$\varvec{x}_{k|k} = \varvec{x}_{k|k - 1} + \varvec{K}_{k} (\varvec{z}_{k} - \hat{\varvec{z}}_{k} )$$

(A20)

The updated covariance matrix is determined using

$$\varvec{C}_{k|k} = \varvec{C}_{k|k - 1} - \varvec{K}_{k} \varvec{C}_{{z_{k} z_{k} }} \varvec{K}_{k}^{\text{T}}$$

(A21)