A Bi-Level Model for Detecting and Correcting Parameter Cyber-Attacks in Power System State Estimation

Abstract

Power system state estimation is an important component of the status and healthiness of the underlying electric power grid real-time monitoring. However, such a component is prone to cyber-physical attacks. The majority of research in cyber-physical power systems security focuses on detecting measurements False-Data Injection attacks. While this is important, measurement model parameters are also a most important part of the state estimation process. Measurement model parameters though, also known as static-data, are not monitored in real-life applications. Measurement model solutions ultimately provide estimated states. A state-of-the-art model presents a two-step process towards simultaneous false-data injection security: detection and correction. Detection steps are ${\chi }^2$ statistical hypothesis test based, while correction steps consider the augmented state vector approach. In addition, the correction step uses an iterative solution of a relaxed non-linear model with no guarantee of optimal solution. This paper presents a linear programming method to detect and correct cyber-attacks in the measurement model parameters. The presented bi-level model integrates the detection and correction steps. Temporal and spatio characteristics of the power grid are used to provide an online detection and correction tool for attacks pertaining the parameters of the measurement model. The presented model is implemented on the IEEE 118 bus system. Comparative test results with the state-of-the-art model highlight improved accuracy. An easy-to-implement model, built on the classical weighted least squares solution, without hard-to-derive parameters, highlights potential aspects towards real-life applications.

1. Introduction

The Power System State Estimator (PSSE) is a major tool for real-time grid monitoring. The end-goal of PSSE is to estimate the system states, typically buses complex voltages, given a set of measurements. Several protection schemes and grid functionality rely on the output of the State Estimation (SE) process. The main inputs of PSSE are measurements set and model parameters. The first is a collection of different system measurement types. For example, circuit breaker status, real and reactive power flows, real and reactive power injection, and voltage magnitudes. The measurement model parameters represent the components of the underlying physical system. With any perturbation in the measurement set and/or model parameters, the PSSE will result in a wrong estimate of the system states. Much research has addressed measurement cyber-attacks. These are usually modeled as False Data Injection (FDI). Measurement model parameters cyber-attacks, on the other hand, have limited research in the field of power systems. In fact, these parameters are considered static and without error during the SE process. Hence, no monitoring scheme is presented in real-life applications. These parameters are prone to cyber-attacks. The cyber-attack in this context could be in the form of an external entity who is able to access the database and alter some of those parameters, or an internal entity who is able to gain super user privileges to change the database [1,2,3,4,5]. The former is a class of cyber-attack called Remote to User attack (R2U) while the latter is known as User to Root attack (U2R).

In the literature, research on detecting FDI pertaining to SE measurements is much explored [3,6,7,8,9]. The work in [10,11,12,13] investigated FDI attack in measurements only. Moreover, a DC model state estimation is considered. The work in [14] considered attack into states in addition to measurement FDI. However, the DC model assumes that states are linearly related to measurements. In addition, voltage magnitudes are assumed to be 1 pu. Such assumption is not accurate in some studies where accurate system model is needed. The AC state estimation, on the other hand, provides an accurate model compared to DC state estimation, since the relationship between states and measurements is non-linear. The work in [15] proposed a convexification framework for the AC state estimation based on semi-definite programming (SDP) for solving cyber attack pertaining measurements sensors. In addition to modeling solutions, Machine Learning (ML) based solutions are also presented [16,17,18]. The problem of detecting cyber-attacks in the measurement model parameters, on the other hand, has been much less considered [19]. Further, the presented solutions considered that cyber-attacks on measurements have been already corrected [20]. However, if a simultaneous attack happened, i.e., on measurements and parameters, how can a measurement correction be made? Existing work towards parameter cyber-attacks [21,22,23,24] considers a two-step approach: detection and correction. In the detection step, the measurements’ residual is analyzed and a pattern is extracted. An attack to a line parameter would result in the normalized residual of the measurements associated with that line to have a higher value compared to the other measurements, assuming no FDI attack [25]. In the correction step, the line’s parameters are corrected in an iterative process using WLS in conjunction with Taylor series expansion. After correction, a SE routine is executed again to check if the normalized residual test does not detect errors. Otherwise, the correction routine is repeated until SE does not flag. In [26], errors on system parameters are addressed while estimating system states. Hence, an augmented objective function is built on the minimization of measurement and parameter residuals. While it is effective to have such a state estimator in a single level model, and eliminating post-processing detection algorithms, the work in [26] assumed errors in parameters are varied in a small range, not considering the possibility of R2U and U2R attacks that enables an adversary to alter those parameters in any range. In addition, the final estimate is sensitive to Gaussian noise in the measurements set and extended redundancy due to the increase size of the state vector.

The aforementioned solutions come with the cost that a non-linear system is linearized using Taylor series expansion and solved in an iterative process to estimate the system parameters. In addition, a simultaneous parameter attack would result in estimating all suspicious parameters under attack in a sequential order. Thus, the correction of one attack depends on the other. Hence, the choice of what attack to correct first might influence the result while there is no guarantee of convergence to the correct physical solution.

In this work, a simultaneous cyber-attack detection and correction bi-level model is presented, towards the solution of previously mentioned state-of-the-art limitations. The bi-level model combines the two steps in a single optimization framework. The presented framework takes advantage of the temporal and spatio characteristic of the grid. In addition, the formulated optimization problem eliminates the effect of the presence of measurements Gaussian noise on parameter correction. Hence, the contribution of this paper towards the state-of-the-art are two-fold:

• An explicit mathematical bi-level model for detecting and correcting cyber-attack pertaining state estimator static data.
• Using the temporal and spatio characteristics of the grid to eliminate non-linearity in parameter correction and providing a sliding-window for an online monitoring scheme of the measurement model parameters.

The remainder of this paper is organized as follows. Section 2 presents background theory on the SE and measurement and parameter attack modelling. Bi-level model and framework is presented in Section 3. Section 4 presents a case study and concluding remarks are provided in Section 5.

2. Background

2.1. State Estimation

AC State estimation aims solving a non-linear algebraic differentiable set of equations that have the following form [27]:

$$\mathbf{z}=\mathbf{h}\left(\mathbf{x}\right)+\mathbf{e}.$$

(1)

where $\mathbf{z}\in {\mathbb{R}}^m$ is the measurement vector, $\mathbf{x}\in {\mathbb{R}}^N$ is the state variables vector (typically voltage magnitudes V and voltage angles $\theta$), $\mathbf{h}\left(\mathbf{x}\right)$:${\mathbb{R}}^m\to {\mathbb{R}}^N,(m>N)$ is a non-linear differentiable function that relates the states to the measurements, $\mathbf{e}$ is the measurement error vector assumed with zero mean, standard deviation $\sigma$ and having Gaussian probability distribution, and $N=2n-1$ is the number of unknown state variables and n is the number of buses in the system. Hence, in classical Weight Least Square State Estimation (WLS SE), the approach consists of solving the following minimization problem:

$$\begin{array}{cc}\hfill {\displaystyle \underset{\mathbf{x}}{min}}& \multicolumn{1}{c}{J\left(\mathbf{x}\right)={[\mathbf{z}-\mathbf{h}\left(\mathbf{x}\right)]}^T\mathbf{W}[\mathbf{z}-\mathbf{h}\left(\mathbf{x}\right)].}\end{array}$$

(2)

where $\mathbf{W}$ is a diagonal weight matrix composed by the inverse of the squared values of measurement standard deviations ($\sigma$): $\mathbf{W}$ = $diag({[{\sigma }_1^{-2},\dots ,{\sigma }_m^{-2}]}^T)$. $J\left(x\right)$ index is a norm in the measurements vector space.

The measurement model in (1) relies on two data sets: measurements set and grid graph, i.e, connectivity and system parameters. If corrupted data is used, then the obtained solution will mislead the operators who monitor the grid. Corrupted data could be attributed to measurement(s) and/or system parameters (database). Given the non-linear relationship, it would be a difficult task to distinguish the source of bad data when there is a simultaneous attack [25]. Hence, in this work, the way is paved for the model to be able to clearly distinguish the error source in the measurement model seen in (1), i.e., is the FDI on the measurement set, system model parameters, or both, and how to correct this?

The database in this context is the model representation of different components that compose the physical power grid. For instance, a typical model of a long transmission system line is represented by the $\pi$-model. Hence, in SE, this model contributes to the bus admittance matrix, i.e., $Y_{bus}$ through its parameters such as line conductance $g_{km}$, line susceptance $b_{km}$ and shunt admittance $b_{km}^{sh}$. Depending on the system under study, a combination of those parameters might be considered. For instance, in short and medium transmission lines, $b_{km}^{sh}$ has a negligible effect on the voltage. Hence, it could be excluded from the model. For long transmission lines, however, $b_{km}^{sh}$ is important for estimating the voltage. The challenging scenario is when all parameters are included. Therefore, with any perturbation in these parameters, the state estimator might lead to a solution that does not depict the true underlying physical system. The task would be even more challenging when both measurements and parameters have contributed to estimate an untrue states, i.e., V and $\theta$, how one could identify the source of erroneous with confidence?

The classical WLS model in (2) minimizes the residual. The work in [28] proved, however, that the error in (1) has a unique decomposition; detectable and undetectable components. The error can be written as follows

$${\parallel e\parallel }^2={\parallel e_D\parallel }^2+{\parallel e_U\parallel }^2.$$

(3)

where $e_D$ is the detectable error while $e_U$ is the undetectable error. Hence, the Innovation concept, i.e., $II$, is used to quantify the undetectable part as follows:

$${II}_i=\frac{\parallel e_D^i\parallel }{\parallel e_U^i\parallel }=\frac{\sqrt{1-P_{ii}}}{\sqrt{P_{ii}}}.$$

(4)

where $P_{ii}$ is the ith entry in the projection matrix P. The P matrix is obtained based on the Jacobian matrix $H=\frac{\partial h}{\partial x}$ and measurements weight W calculated as follows:

$$P=H{\left(H^{T}WH\right)}^{-1}{H}^{T}W.$$

(5)

Hence, the error in (3) is then composed by using the Innvoation Index in (4) to obtain the Composed Measurement Error $CME$ in its normalized form for each measurement i as follows:

$$CM{E}_i^N=\frac{r_i}{{\sigma }_i}\left(\sqrt{1+\frac{1}{{I{I}_i}^2}}\right).$$

(6)

where $r_i$ is the ith measurement mismatch which is the detectable part of the error, and ${\sigma }_i$ is the standard deviation of the ith measurement. Therefore, the minimization problem in (2) should minimize the composed error in (6) instead of the residual [21].

2.2. Bi-Level Optimization

Bi-level optimization is a mathematical programming framework where a constraint in an optimization problem is another optimization problem. The main optimization problem is generally called upper (leader) model while the constraint which is another optimization problem is called lower (follower) model. This type of optimization framework arises in situation where hierarchical decision-making is involved. In other words, a decision from one task affects the decisions of the other task and vice versa. This framework has two types of variables, the upper-level variables and the lower-level variables [29].

3. Framework

The SE process is run every 60–90 s to monitor the status of the grid [27]. After every run, an estimate of system states (typically complex bus voltages) and measurements are obtained. Processing these outputs would yield valuable temporal information considering the next run. Hence, this paper addresses the following question: knowing prior states and database, can one retrieve the current database? To address this question, a model is constructed based on the non-linear algebraic equations used in AC SE.

3.1. Preliminaries

Consider a transmission line connecting bus k and bus m, and represented in a $\pi$-model. With the line admittance $y_{km}$, the conjugate of the complex power flow through that line can be written as [27]:

$$S_{km}^{*}=P_k-j{Q}_k=E_k^{*}{I}_{km}.$$

(7)

where $E_k$ is the complex voltage at bus k, $I_{km}$ is the complex current flowing from bus k to bus m, and the * indicates the conjugate of the complex quantity. Using $I_{km}=(E_k-E_m)y_{km}$, we can write the complex power as:

$$S_{km}^{*}=V_k^2{y}_{km}{V}_k{e}^{-j{\theta }_k}(V_k{e}^{j{\theta }_k}-V_m{e}^{j{\theta }_m})+j{b}_{km}^{sh}.$$

(8)

where $y_{km}$ is the admittance between bus k and bus m, $V_k$ and $V_m$ are the magnitudes of the complex voltages at bus k and m, respectively, ${\theta }_k$ and ${\theta }_m$ are the angles of the complex voltages at bus k and m, respectively, and $b_{km}^{sh}$ is the shunt admittance of the line connecting bus k and bus m. Expanding the right hand side of (7) and decomposing the expression into real and imaginary parts, one can obtain the following:

$$P_{km}=(V_k^2-V_k{V}_{m}cos{\theta }_{km})g_{km}-(V_k{V}_{m}sin{\theta }_{km})b_{km}.$$

(9)

$$Q_{km}=(-V_k{V}_{m}sin{\theta }_{km})g_{km}+(-V_k^2)b_{km}^{sh}+(V_k{V}_{m}cos{\theta }_{km}-V_k^2)b_{km}.$$

(10)

where $g_{km}$ is the real part of the line admittance connecting bus k and bus m, i.e., $\Re \left\{y_{km}\right\}$, and $b_{km}$ is the imaginary part of the line admittance connecting bus k and bus m, i.e., $\Im \left\{y_{km}\right\}$. Equations (9) and (10) represent the real and reactive power flows in the line connecting bus k and bus m, respectively. With the real power flow from bus m to bus k, i.e., $P_{mk}$ (by changing bus index), one can express the real power loss of the same line as:

$$\begin{array}{cc}\hfill P_{km}^{loss}& =P_{km}+P_{mk}\hfill \\ & =g_{km}(V_k^2+V_m^2-2V_k{V}_{m}cos{\theta }_{km})\hfill \\ & =g_{km}{|E_k-E_m|}^2.\hfill \end{array}$$

(11)

where $E_k$ and $E_m$ are complex voltages at bus k and bus m, respectively. Similar procedure to (11), an expression of the reactive power loss in the line can expressed as:

$$\begin{array}{cc}\hfill Q_{km}^{loss}& =Q_{km}+Q_{mk}\hfill \\ & =-(V_k^2+V_m^2)b_{km}^{sh}-(V_k^2+V_m^2-2V_k{V}_{m}cos{\theta }_{km})b_{km}\hfill \\ & =-(V_k^2+V_m^2)b_{km}^{sh}-({|E_k-E_m|}^2)b_{km}.\hfill \end{array}$$

(12)

Equations (9)–(12) are the basic equations that govern a line connecting bus k and bus m from the SE perspective. In AC SE, Equations (9) and (10) are used when $P_{km}$ and/or $Q_{km}$ are present in the measurements set. The relationship between $g_{km}$ and $b_{km}$ can be derived from the actual impedance of the line as:

$$\begin{array}{cc}\hfill y_{km}& =g_{km}+j{b}_{km}=\frac{1}{z_{km}}=\frac{1}{R_{km}+j{X}_{km}}.\hfill \end{array}$$

(13)

$$\begin{array}{cc}\hfill g_{km}& =\Re \{\frac{1}{z_{km}}\}=\frac{R_{km}}{R_{km}^2+X_{km}^2}.\hfill \end{array}$$

(14)

$$\begin{array}{cc}\hfill b_{km}& =\Im \{\frac{1}{z_{km}}\}=-\frac{X_{km}}{R_{km}^2+X_{km}^2}.\hfill \end{array}$$

(15)

If (14) is divided by (15), the following expression is obtained:

$$\begin{array}{c}\hfill \frac{b_{km}}{g_{km}}=-\frac{X_{km}}{R_{km}}.\end{array}$$

(16)

Therefore, (16) correlates line conductance to its susceptance. Further, the term $\left(\frac{X_{km}}{R_{km}}\right)$ is known as $X/R$ ratio and commonly used in short circuit studies. In transmission systems, this ratio is higher compared to distribution systems. In addition, this ratio is a characteristic of the line that indicates the tangent angle between line resistance and line inductance. Hence, having this factor will eliminate the non-linearity in retrieving the original measurement model parameters that will be presented in Section 3.3.

3.2. Cyber-Attack Model

With the mathematical concepts presented in Section 3.1, a FDI in line parameters can be modeled. Consider a line connecting buses k and m has a FDI in its parameter model. Then, this FDI can be modeled as follows:

$$\begin{array}{c}\hfill g_{km}^{pert}=g_{km}+\mathsf{\Delta }{g}_{km}.\end{array}$$

(17)

$$\begin{array}{c}\hfill b_{km}^{pert}=b_{km}+\mathsf{\Delta }{b}_{km}.\end{array}$$

(18)

$$\begin{array}{c}\hfill b_{km}^{sh,pert}=b_{km}^{sh}+\mathsf{\Delta }{b}_{km}^{sh}.\end{array}$$

(19)

where $g_{km}$, $b_{km}$, and $b_{km}^{sh}$ are the true line parameters, $\mathsf{\Delta }{g}_{km}$, $\mathsf{\Delta }{b}_{km}$, and $\mathsf{\Delta }{b}_{km}^{sh}$ are the deviation (due to attack) in line parameters, and $g_{km}^{pert}$, $b_{km}^{pert}$, and $b_{km}^{sh,pert}$ are the perturbed quantities. By substituting (17)–(19) into (9) and (10) one can derive:

$$P_{km}^{pert}=(V_k^2-V_k{V}_{m}cos{\theta }_{km})g_{km}^{pert}-(V_k{V}_{m}sin{\theta }_{km})b_{km}^{pert}.$$

(20)

$$Q_{km}^{pert}=(-V_k{V}_{m}sin{\theta }_{km})g_{km}^{pert}+(-V_k^2)b_{km}^{sh,pert}+(V_k{V}_{m}cos{\theta }_{km}-V_k^2)b_{km}^{pert}.$$

(21)

where $P_{km}^{pert}$ and $Q_{km}^{pert}$ are the attacked (deviated) real and reactive power, considering values obtained in (9) and (10), respectively, due to a FDI in line parameter(s). Note that the voltages at buses k and m are the same as the ones estimated to obtain $P_{km}$ and $Q_{km}$ in (9) and (10). Hence, with this notion, the system operators can make use of data already available from SE to further secure the state estimator routine over time. In addition, it can be viewed as a filtering stage prior to run SE routine to validate system database after initialization. Hence, any flag from SE after validating system database would be identified to measurement set considering a previously defined confidence level.

3.3. Bi-Level Optimization Model

Having established the necessary mathematical concepts in Section 3.1 and Section 3.2, an optimization framework for estimating measurement model parameters (i.e., $g_{km}$, $b_{km}$ and $b_{km}^{sh}$) for any line connecting bus k and bus m can be formulated. The framework hypothesis that a free of attack SE output sample exists. Let us label this sample with $t^{-}$. Hence, at time $t^{-}$, system states are estimated (i.e., $E_k^{t^{-}}$ and $E_m^{t^{-}}$). If {$P_{km}$ or $P_{mk}$} and {$Q_{km}$ or $Q_{mk}$} are part of the measurement set, then estimated measurements $h_{km}^P$ and $h_{km}^Q$ are already available. If not, an estimated measurement out of $\{P_{km},P_{mk}\}$ and an estimated measurement out of $\{Q_{km},Q_{mk}\}$ are generated after SE is converged. This step can be augmented to the existing SE routine without a major modification. Therefore, the bi-level model can be derived as:

$$\begin{array}{cc}\hfill \underset{x_u}{min}& \sum _{i=1}^m\left(1+\frac{1}{{I{I}_i(x_u,x_l)}^2}\right)W_{ii}{r}_i^2(x_u,x_l).\hfill \end{array}$$

(22)

$$\begin{array}{cc}\hfill \mathrm{s}.\mathrm{t}.r_i& =z_i-h_i(x_u,x_l),\forall i=1,2,\dots m\hfill \end{array}$$

(23)

$$\begin{array}{cc}\hfill g_{km}& =g_{km}^{pert}-\mathsf{\Delta }{g}_{km},\forall km\in \mathcal{L}\hfill \end{array}$$

(24)

$$\begin{array}{cc}\hfill b_{km}& =b_{km}^{pert}-\mathsf{\Delta }{b}_{km},\forall km\in \mathcal{L}\hfill \end{array}$$

(25)

$$\begin{array}{cc}\hfill b_{km}^{sh}& =b_{km}^{sh,pert}-\mathsf{\Delta }{b}_{km}^{sh},\forall km\in \mathcal{L}\hfill \end{array}$$

(26)

$$\begin{array}{cc}\hfill x_l& \in \Psi \left(x_u\right)\hfill \end{array}$$

(27)

where ${\mathbf{x}}_{\mathbf{u}}$ is the decision variable vector for the upper-level optimization problem, i.e., voltage magnitude V and voltage angle $\theta$ for all buses, and ${\mathbf{x}}_{\mathbf{l}}$ is the decision variable vector for the lower-level optimization problems, i.e., deviations in system database $\mathsf{\Delta }{g}_{km},\mathsf{\Delta }{b}_{km},$ and $\mathsf{\Delta }{b}_{km}^{sh}$ for all lines. The variable $\mathcal{L}$ is the set of lines in the system under study, and $\Psi \left(x_u\right)$ is a parameterized range constraint for the lower-level decision vector $x_l$. Such constraint is obtained through the lower-level (follower) optimization problem defined as follows:

$$\begin{array}{cc}\hfill \underset{x_l}{min}& \mathsf{\Delta }{g}_{km}+\mathsf{\Delta }{b}_{km}+\mathsf{\Delta }{b}_{km}^{sh}\hfill \end{array}$$

(28)

$$\begin{array}{cc}\hfill \mathrm{s}.\mathrm{t}.\mathsf{\Delta }{g}_{km}& =g_{km}^{pert}-g_{km}\hfill \end{array}$$

(29)

$$\begin{array}{cc}\hfill \mathsf{\Delta }{g}_{km}& =g_{km}^{pert}-g_{km}\hfill \end{array}$$

(30)

$$\begin{array}{cc}\hfill \mathsf{\Delta }{b}_{km}& =b_{km}^{pert}-b_{km}\hfill \end{array}$$

(31)

$$\begin{array}{cc}\hfill \mathsf{\Delta }{b}_{km}^{sh}& =b_{km}^{sh,pert}-b_{km}^{sh}\hfill \end{array}$$

(32)

$$\begin{array}{cc}\hfill b_{km}& =g_{km}{\left(\frac{X}{R}\right)}^{ratio,km}\hfill \end{array}$$

(33)

$$\begin{array}{cc}\hfill P_{km}& =\left(f_{P_{km}}^g\right)g_{km}+\left(f_{P_{km}}^b\right)b_{km}\hfill \end{array}$$

(34)

$$\begin{array}{cc}\hfill Q_{km}& =\left(f_{Q_{km}}^g\right)g_{km}+\left(f_{Q_{km}}^b\right)b_{km}+\left(f_{Q_{km}}^{b^{sh}}\right)b_{km}^{sh}\hfill \end{array}$$

(35)

$$\begin{array}{cc}\hfill P_{mk}& =\left(f_{P_{mk}}^g\right)g_{km}+\left(f_{P_{mk}}^b\right)b_{km}\hfill \end{array}$$

(36)

$$\begin{array}{cc}\hfill Q_{mk}& =\left(f_{Q_{mk}}^g\right)g_{km}+\left(f_{Q_{mk}}^b\right)b_{km}+\left(f_{Q_{mk}}^{b^{sh}}\right)b_{km}^{sh}\hfill \end{array}$$

(37)

$$\begin{array}{cc}\hfill P_{km}^{pert,loss}& =P_{km}+P_{mk}+(f_{P_{km}}^g+f_{P_{mk}}^g)\mathsf{\Delta }{g}_{km}+(f_{P_{km}}^b+f_{P_{mk}}^b)\mathsf{\Delta }{b}_{km}\hfill \end{array}$$

(38)

$$\begin{array}{cc}\hfill P_{km}^{pert,loss}& =\left({|E_k^{t^{-}}-E_m^{t^{-}}|}^2\right)g_{km}^{pert}\hfill \end{array}$$

(39)

$$\begin{array}{cc}\hfill Q_{km}^{pert,loss}& =Q_{km}+Q_{mk}+(f_{Q_{km}}^b+f_{Q_{mk}}^b)\mathsf{\Delta }{b}_{km}+(f_{Q_{km}}^{b_{km}^{sh}}+f_{Q_{mk}}^{b_{km}^{sh}})\mathsf{\Delta }{b}_{km}^{sh}\hfill \end{array}$$

(40)

$$\begin{array}{cc}\hfill Q_{km}^{pert,loss}& =\left(f_{Q_{km}}^g\right)g_{km}^{pert}+\left(f_{Q_{km}}^b\right)b_{km}^{pert}+\left(f_{Q_{km}}^{b^{sh}}\right)\mathsf{\Delta }{b}_{km}^{sh,pert}\hfill \end{array}$$

(41)

$$\begin{array}{cc}\hfill g_{km},b_{km}^{sh}& \ge 0\hfill \end{array}$$

(42)

$$\begin{array}{cc}\hfill b_{km}& \le 0\hfill \end{array}$$

(43)

In the upper-level model, the weighted norm of the error at time $t^{-}$ is minimized [28]. After $\mathsf{\Delta }t$ seconds, the inner-level model, the parameters delta $g_{km}^{pert}$, $b_{km}^{pert}$, and $b_{km}^{sh,pert}$, which are the current status of the database at time $t=\mathsf{\Delta }t$ + $t^{-}$, which the system operator would like to check, are optimized. The variables $g_{km}$, $b_{km}$, and $b_{km}^{sh}$ are the unknown true states of the database that we seek to obtain. The ${\left(X/R\right)}^{\left(ratio\right)}$ is the known ratio of line inductance to line resistance. The function $f_{mea{s}_{km}}^{param}$ is a function evaluation of the coefficient associated with the given parameter $param$ from bus k to bus m for the specified measurement type $meas$ as (9) and (10). The inner model is evaluated using the states V and $\theta$ of the two buses connecting line $km$ at previous time $t^{-}$. $P_{km}^{pert,loss}$ and $Q_{km}^{pert,loss}$ are losses in the line evaluated given the states at time $t^{-}$ and the current status of system database at time t. In (34)–(37), only one estimated measurement of each type at time $t^{-}$ is required. The other two can be free to be obtained by the chosen optimization solver.

From the previous bi-level model, line parameters can be obtained independently from each other. This allows the system operators to take advantage of parallel computation. In addition, the inner optimization problem is linear in its decision variables. The constraint (33) ensures the optimal solution of parameter values are unique and correspond to the correct physical solution. Hence, any off-the-shelf solver can be used to seek solution. The flowchart of the presented framework is shown in Figure 1.

As illustrated in Figure 1, from the prospective of SE, the process starts by uploading data of measurements and system model parameters. The SE routine is executed by system operator every often to monitor the grid. The framework presented in Section 3.3 is initialized with a true sample that is free from measurement and parameter errors. This sample is labeled as $t^{-}$. Then, for a sample t, SE routine is performed. On such, the bi-level model is executed. If SE detects an error, then the presented inner (lower) level model in Section 3.3 is performed to check if the error source is due to measurement model parameters. To do so, current status of the measurement model parameters are sent to the presented model to be executed. After execution, if errors in line parameters are above certain threshold, defined considering a level of confidence, then the corresponding line is updated to the solution obtained by the model in Section 3.3 and SE is executed again. Otherwise, no parameter error is detected [19]. If error is detected after updating measurement model parameters, then the source of this error is due to errors in measurement. In such case, [30] is run. After correcting errors from data at sample t, the base data in the presented model can be updated if the sample t is trusted by system operator. Considering Figure 1, the contribution of this work towards the WLS SE state-of-the-art process is highlighted with the boxes colored in green.

4. Case Study

The presented bi-level model was validated using the IEEE 118-bus system. By using the MATLAB package MATPOWER [31], 21,600 samples (i.e., one day’s worth) of measurements were generated with Gaussian noise based on a common daily load profile that contains temporal information of a power system’s changing state. The measurement set includes real and reactive power flows, power injections, and all voltage magnitudes, resulting in 712 measurements with Global Redundancy Level (GRL = $m/N$) of 3.029, which relates the number of measurements (m) to the number of states (N) to be estimated. Measurement’ standard deviations are considered as 1% of their absolute values. For optimization, Gurobi solver [32] is used for solving the bi-level model. All simulations are conducted on a personal Apple Mac computer: macOS High Sierra 32 GB RAM 1876 MHz DDR3, 4 GHz Intel Core i7.

Towards validation, five independent 100 Monte Carlo simulations were conducted for a selected sample. In each simulation, a line is selected randomly to have cyber-attacks, modeled as FDI in model parameters, i.e, $g_{km}$, $b_{km}$, and $b_{km}^{sh}$. The size of the cyber-attacks is drawn from a uniform distribution between $\pm 5$% and $\pm 40$% of their actual values. The optimization framework presented in Section 3.3 and Figure 1 is conducted after each attack. Case study results are presented in Figure 2. Figure 2 shows that the absolute error after correcting line parameters is less than an order of 3.

To further evaluate the accuracy and performance of the presented bi-level model, around 20% of the samples (out of 21,600) are selected randomly to be compromised with parameter cyber-attacks. Each of those samples, a random line is selected to have a FDI parameter attack. The attack is in the same range as those performed for the aforementioned simulations. The confusion matrix for the SE output using ${\chi }^2$ test as a detection method is illustrated in

Table 1. SE performance result for parameter attack detection ${\chi }^2$ test [19].

		Prediction Outcome
Actual value		Normal	Anomaly
	Normal	17336	0	Normal
		sample	sample
	Anomaly	3252	1012	Anomaly
		sample	sample
		Normal	Anomaly

. The ${\chi }^2$ threshold is calculated based on two parameters: number of measurements and confidence level. In this test, the number of measurements is 712 and the confidence level is chosen to be 95% [19]. As seen, a substantial amount of samples were not detected by ${\chi }^2$ test. Meanwhile, the presented bi-level model was executed after each SE run. All anomaly samples were not only detected, but also corrected in a single optimization run. Observed errors in correction were similar to the results shown in Figure 2. The execution time of the proposed model was monitored for all anomaly samples. On average, for 170 lines, the total execution time was 0.3964 s with a standard deviation of 0.0533 s. It is worth mentioning that these reported statistics are without using parallel computation. Hence, a lower execution time could be achieved with parallelism.

The $CM{E}^N$ methodology presented in [30] for parameter attack processing, which is the composed measurement error $CME$ in its normalized form, is also explored in the comparative test case scenarios. An anomaly sample is selected and the resulted $CM{E}^N$ of the measurements were listed in a descending order based on their absolute values for a threshold value of 3. In this sample, the underlying true attack is on line connecting bus 94 and bus 95. The result is shown in

Table 2. Parameter Cyber-attacks Identification [30].

Measurement	From Bus	To Bus	${\mathit{CME}}^{\mathit{N}}$
Real Power Flow	96	95	10.093
Reactive Power Flow	95	96	9.5299
Reactive Power Flow	94	95	7.9748
Reactive Power Flow	94	96	7.7034
Real Power Flow	94	95	6.3127
Real Power Flow	94	96	5.8595
Real Power Injection	95	95	4.0285

. Based on the strategy presented in [30], the attack is characterized as a parameter attack. However, not a specific line is determined as the one that is compromised. Instead, a region where the attack might be at could be inferred. Hence, the superiority of the proposed framework is that it can identify and correct the attack in a single process. In addition, it can be used as a pre-processing step prior executing SE routine.

For stealthy attack, a line is selected and its parameters, i.e., g, b and $b_{sh}$ are attacked gradually from 0 to 20% of their values. The performance index J as well as the $CME$ in its normalized form ($CM{E}^N$) are recorded. The results are shown in Figure 3 and Figure 4. In Figure 3, the performance index $J\left(x\right)$ (colored in blue) increased with the increase size of the attack in the line’s parameters under attack. In this case, even though the performance index $J\left(x\right)$ increased, the ${\chi }^2$ test still did not detect the error. For identification, the $CM{E}^N$ is obtained for every attack and the absolute error is calculated and presented in Figure 4. As shown, due to the increase size of the attack in a single line, the error is spread into multiple estimation of measurements. After each attack scenario, the bi-level model is performed. The error due to correction of parameters is calculated and shown in Figure 5.

The same scenario of the previous stealth attack is simulated for multiple lines in this case. The results are shown in Figure 6 and Figure 7. As seen from the figures, a similar trend has occurred. However, the errors in measurement estimation are increased. The bi-level model is performed and lines are corrected. The observed error in correction for the stealthy attacks is presented in Figure 8.

5. Conclusions

This paper presents a bi-level model for correcting parameter FDI cyber-attacks on the SE process. The presented model combines the two processed that are usually performed by SE for detection and correction into a single process for parameter attack processing. The presented model can be used as a post-state estimation cyber-attack processing or prior to validate the database of measurement model parameters and measurements. Meanwhile, the framework can be used as an online tool due to the capability of performing parallel computations. In addition, most the information needed in this framework is already available among the data set used by SE. Comparative test results on the IEEE 118-bus system show that the presented model is able to correct parameters with high accuracy, while further processing measurement cyber-attacks. The existing state estimator software can be adjusted to incorporate the presented framework without major modifications, enabling the current work to be utilized by utilities. The model can be solved by solvers that do not require sophisticated features.