Nomenclature

[]

Number of buses

Number of transmission lines

Number of lines connected to bus

Set of buses connected to bus

Set of buses with PMUs installations

Number of PMUs attacked

Binary vector of PMU locations in a network

Binary vector of attacked PMUs

Regression matrix for bus

Noise vector at bus

Noise covariance for the measurement at bus

Voltage phasor at bus

Real and imaginary parts of

Magnitude and angle of

Vector collecting for all buses

Vector collecting for all buses

Vector

Current phasor on line

Magnitude and angle of

Maximum likelihood (ML) estimate of system state

Noiseless measurement at bus

Noisy attacked measurement at bus

ML state estimate using attacked measurements

Expected value of

Biasscaling matrix induced by attack

Attack angle at bus

Vector collecting for all buses

Block diagonal matrix relating and

Vector
I Introduction
Phasor Measurement Units (PMUs) equipped with GPS receivers are installed ubiquitously in smart grids, replacing and augmenting traditional sensors of the Supervisory Control and Data Acquisition (SCADA) systems. The higher sampling rates of PMUs compared to SCADA systems assist the network operator to perform real time Wide Area Monitoring, Protection and Control (WAMPAC)—a unique feature of smarter power grids.
Cyber attacks on PMUs include False Data Injection (FDI) attacks and GPS spoofing attacks. In an FDI attack, the attacker tries to inject false data that are not detectable by bad data detection algorithms into the network; see e.g., the recent works [1, 2]. GPS spoofing is caused by transmitters mimicking the GPS signal with the intention of altering the GPS time estimated by the PMU’s GPS receiver [3, 4]. These attacks maliciously introduce erroneous time stamps, thereby inducing a wrong phase angle in the PMU measurements [5], and are also called time synchronization attacks (TSAs).
The Electric Power Research Institute, in collaboration with the National Electric Sector Cybersecurity Organization Resource, has published a technical report recognizing the vulnerability of PMUs to GPS spoofing under its scenario WAMPAC.12: GPS Time Signal Compromise [6]. This paper considers the effects of GPS spoofing attacks on power grids.
Sensors such as PMUs cannot be placed at every bus location in the power network as they are still costly. Thus they are placed in crucial network locations to ensure the visibility of the network. Consequently, a system operator can generate information about other, nonsensed states using state estimation (SE) routines [7, 8, 9, 10]. The purpose of SE is threefold: protection, operational control, and grid status verification [9]. In a network, measurements from dispersed PMUs are time synchronized at the Phasor Data Concentrator (PDC) and fed to the super PDC for further processing. Thus, systematic attacks on PMUs can have a catastrophic effect on the SE and consequently on the grid.
This paper analyzes the vulnerability of smart grids to GPS spoofing by formulating an optimization problem to identify vulnerable PMUs in the network, and by developing an algorithm for joint state estimation and attack reconstruction. The ensuing section provides a review of the related literature, and details the contributions of this paper.
Ii Literature Review & Paper Contributions
This section categorizes the GPS spoofing attack literature in three groups, namely, the feasibility of spoofing attacks, impacts to the power grid, and countermeasures, and details the contributions of this paper.
The experimental and theoretical feasibility of GPS spoofing attacks is shown in [11, 12, 13, 14, 15]. More specifically, [11] demonstrates the implementation of a laboratorybased GPS spoofer and discusses mechanisms against civilian GPS spoofing. A spoofer implemented in US DOE’s Pacific Northwest National Laboratory is the theme of [12]. The vulnerabilities of PMUs to GPS spoofing are studied in [13] and [14]. The work in [15] manipulates the GPS navigation data to achieve various timing and phase errors in PMU measurements.
The impacts of GPS spoofing attacks on system operations are studied in [16, 15, 13]. Specifically, the effect of GPS spoofing on fault location and voltage stability monitoring algorithm is demonstrated in [15] and [16]. A falsely activated generator trip scenario—as a result of a GPS spoofing attack—is constructed in [13]. In our previous work [17], the statistics of state estimates affected by GPS spoofing are studied.
Countermeasures to spoofing attacks are discussed in [18, 19, 20, 21, 22, 23]. A detection method for TSAs on multiple PMUs using a cross layer detection mechanism is developed in [18]. The work in [19] presents a method to detect hazardous data via signalbased and modelbased methods. Furthermore, a GPS spoofing attack detection scheme based on collaboration among multiple PMUs in a large grid is the theme of [20]. Recommendations for improving spoofing detection in commercial receivers are articulated in [21]. The work in [22] takes advantage of the PMU locations along with the statistics of GPS receivers to detect the spoofing attacks. The work in [23] proposes methods to prevent propagation of cyberattacks in PMU networks.
State estimation and GPS spoofing attack identification in power networks is pursued in [24, 25]. A computationally efficient algorithm to identify an attack on at most one PMU is developed in [24]. The generalized likelihood ratio test (GLRT) is the basis of the method developed in [25], but implementation of GLRT requires the solution of a hard optimization problem with respect to the unknown time delay induced by GPS spoofing. Beyond the two aforementioned works, which specialize in GPS spoofing, recent works that deal with attack identification in power networks include [2, 26, 27, 28]. The majority of these studies require successive measurements across time and detailed dynamical models of the network, including generator parameters. The contributions of this paper are as follows.

A measurement model that explicitly relates the PMU measurements with the network state and the (potential) GPS spoofinginduced phase shifts is developed. This model is leveraged to derive the statistics of the state estimator when GPSspoofed PMU measurements are utilized.

A novel greedy algorithm to identify the most vulnerable PMU locations is developed. The vulnerability is quantified by the state estimation error, which is characterized explicitly in terms of the spoofing attack. Exhaustive search is used to validate the performance of the greedy algorithm. Our analysis is intended to guide the system operator towards increased protection of the network’s most vulnerable PMU locations.

An algorithm to jointly perform state estimation and attack angle reconstruction is developed. The algorithm is based on alternating minimization of a bilinear least squares objective. Lagrangian duality is leveraged so that—despite the nonconvexity—the algorithm features closedform updates. This property renders the algorithm attractive for realtime implementation.

The algorithm is extended to state estimation and attack identification using combined SCADA and GPSspoofed PMU measurements.
Numerical tests indicate that the alternating minimization algorithm can yield smaller state estimation error than the largest normalized residual test (LNRT), which is a classical method for bad data identification [29, Sec. 4.8.4]. The algorithm is also compared to one of the previous GPS spoofing identification approaches [24]. The developed algorithm can identify simultaneous attacks to more than one PMU, in contrast to [24]. Even if 20% of the PMUs are attacked, the algorithm correctly identifies the attacked PMU locations and yields an accurate state estimate, as indicated by tests on standard IEEE transmission networks.
In comparison with [25], this work does not necessitate the solution of a complicated optimization problem with respect to the induced phase shift. The state estimation algorithm is based on a new measurement model that explicitly relates the PMU measurements with the attacked phases, which is different than the additive attack model of [2, 26, 27, 28]. The implication is that the developed algorithm does not require successive measurements across time, and is also able to identify the spoofinginduced attacked phases.
The remainder of the paper is organized as follows. Section III reviews the PMUbased SE. Section IV formulates an optimization problem to identify the most vulnerable PMUs in the network. An algorithm for joint state estimation and attack reconstruction is developed in Section V. Extensions to SE with PMU and SCADA measurements are discussed in VI. Numerical tests on standard IEEE networks are performed in Section VII, and Section VIII concludes the paper.
Iii PMUBased State Estimation
This section outlines the network and measurement model, with and without TSAs, and provides the optimal state estimators using PMU measurements without TSAs.
Iiia Network Model and PMU Measurements
Consider a power network with buses connected via transmission lines. Let be the set of buses connected to bus , and define as the number of lines connected to bus . The system state is the vector of nodal voltages in rectangular coordinates denoted by where and collect the real and imaginary parts and of the complex voltages at buses .
PMUs are installed on select buses of the network; is a binary indicator is equal to 1 if a PMU is installed at bus and 0 otherwise. Vector collects for . The set of buses where PMUs are installed is denoted by A PMU installed at bus measures the bus’s complex voltage as well as the complex currents on all lines that bus is connected to. This collection of measured quantities (in rectangular coordinates) at bus is concatenated in a vector . To make the notation more compact, define as the number of distinct real quantities measured by the PMU at bus .
Let and generically denote the voltage and current phasors at bus and line respectively; and let and denote the corresponding phasor angles. It is convenient for subsequent developments to consider the noiseless version of , which is denoted by :
(1) 
where and are the real and imaginary parts of the complex current injected into line . Note that the current injected into line is different than the current injected into line . To summarize, the noiseless quantities measured at bus comprise the real and imaginary parts of the nodal complex voltage, appended by the real and imaginary parts of the complex currents injected to all lines connected to bus . Using the bus admittance matrix of the network, can be written as a linear function of the system state as . The construction of is provided in [30, 17]. In practice, a PMU at bus measures , which is a noisy version of , i.e.,
where represents an additive Gaussian noise vector that is assumed independent across PMUs and has a known positive definite covariance .
Given that the likelihood of the measurement = is Gaussian, the maximum likelihood (ML) estimate of the system state is given as
(2) 
The optimization in (2) amounts to unconstrained least squares, and can be solved by taking the gradient of the cost function with respect to and setting it to zero, resulting in
(3) 
where it is assumed that the matrix is nonsingular. Invertibility of is tantamount to the state observability, which can be ensured when there is a sufficient number of installed PMUs in the network [30].
Substituting into (3) yields the statistics of the estimator as . That is, the expected value of the estimate is the system state , or in other words the estimator is unbiased. The next section develops a relationship between the measured quantities of bus when there is TSA and their version in the absence of an attack.
IiiB TSAImpacted PMU Measurement Model
As mentioned in [16], TSA affects only the phase of the measurement. Specifically, a TSA on bus introduces a clock offset error [16]. The phase angle error corresponding to the clock offset error for the PMU at bus is denoted by and is given by where Hz. Table I lists typical phase angle errors caused by GPS spoofing attacks reported in the literature.
The noisy attacked PMU measurement at bus is given by
(4) 
Note that there is a potentially different per bus . Combining (4) with (1) and introducing , a linear relationship between and can be derived as:
(5) 
where is a block diagonal matrix^{1}^{1}1Note that depends on , but this dependency is kept implicit in order to keep the notations compact. When
(no attack), we obtain the identity matrix,
. consisting of blocks and each block is the matrix . Therefore, can be thought of as the noiseless measurement that would be available to the PMUenabled bus in the absence of a spoofing attack. The measurement model in (5) connects the state with the attack. It is worth noticing that the measurement is linear in the state, but nonlinear in the attack angle. In Section IV, we use the measurement model (5) to develop a framework which identifies the most vulnerable PMU locations in the network. In Section V, the model in (5) is leveraged to jointly perform the state estimation and attack angle reconstruction from spoofed PMU measurements.Iv Identifying Susceptible PMU Locations
This section formulates an optimization problem to identify the most vulnerable PMUs in the network and develops a computationally attractive solution algorithm. To this end, the statistics of the ML state estimate obtained from the corrupted measurement (5) are characterized next.
Iva Statistics of the Estimates Under Attack
The worstcase scenario is considered, whereby the SE routine does not know (or has not detected) that an attack has occurred. Under this scenario, the SE routine passes the corrupted measurement (5) through the estimator (3). Thus, the state estimate after the attack is given by
(6) 
It should be emphasized that (6) is not the ML estimate that the SE module would derive if the attack magnitude were known. With this observation, we refer to (6) as the attacked ML estimator. In what follows, we derive the statistics of the estimates under attack.
Theorem 1
The attacked ML estimator has the Gaussian distribution
, i.e., it has covariance and expected value(7) 
Proof:
Having derived the distribution of the attacked estimate, the next section details how estimation accuracy metrics typically employed in power systems are affected by the attack.
IvB SE Accuracy Metrics
For any estimator , the mean square error (MSE) matrix is a metric of the estimation accuracy which is introduced in the power systems as early as the seminal work in [31] and has traditionally been utilized in the statistical literature [32].
Definition 1
With denoting the true system state, the mean square error (MSE) matrix is defined by
(8) 
The MSE matrix is formed by the pairwise nodal voltage estimate errors. Its diagonal entries characterize the accuracy of individual nodal voltage estimates. Based on the MSE matrix, it is customary to obtain a scalar metric of the SE accuracy as follows.
Definition 2
The mean square error of an estimator is defined as the trace of the MSE matrix:
(9) 
The MSE sums the squared errors of nodal voltage estimates. The last equality in (9) follows immediately from the linearity of the expectation. The MSE matrix and the MSE are formulated by squaring the difference between the estimate and its true value. This difference yields a bias metric as follows.
Definition 3
The bias of an estimator is defined as
(10) 
Theorem 1 reveals that GPS spoofing introduces a bias in the ML estimator, and is given by where and
(11) 
The matrix analytically characterizes the effect of the attacked PMU phase shift to the estimation bias. The ideal scenario for the bias is to be zero, which occurs in the absence of an attack, yielding and .
The next theorem precisely characterizes the relationship between the MSE and the bias.
Theorem 2
The MSE and the bias of the attacked state estimate satisfy the following relationship:
(12) 
Proof:
The theorem states that as the norm of the bias increases, the MSE increases and vice versa. Recall that depends only on the network topology and the location of the installed PMUs.
Power system protection, control, and status verification are all contingent upon the availability of accurate state estimates. For example, voltage stability assessment (VSA) computes the load change that a system can tolerate before voltage collapse occurs. To perform VSA, the current operating point is obtained from the state estimator [33], and thus corrupted state estimates can have significant impact on the power grid security. It is worth emphasizing that even if a very small number of PMUs are attacked, the estimated voltages at all buses are affected. By leveraging the previously defined SE accuracy metrics, the next section studies how the location and angles of the attacked PMUs affect the quality of the state estimates.
Remark 1
Apart from the MSE and bias, which are general metrics of the SE accuracy, additional indices are used for evaluating the severity of contingencies [34] based on the current network state. Since GPS spoofing disturbs the SE, it also affects the computation of contingency severity indices. The methodology in Sections IVC and IVD can also be applied to these indices.
IvC The Most Vulnerable PMU Location
The objective of this section is to furnish the system operator with an analytical tool to study how the SE accuracy is affected by different attack combinations, and identify the PMUs that can induce the largest bias or MSE if attacked. Theorem 2 reveals that the MSE and the squared norm of the bias are related via an additive constant. Therefore, any of the two metrics can be used to study the effects that the spoofing attacks have on the SE accuracy. An optimization problem to find the attack angle combinations that maximize the norm of the bias is formulated in the sequel.
Specifically, the optimization variables are , which denotes the vector of attacked PMUs, and the attack angles . The optimization problem to maximize the bias, solved by the network operator, is stated next, where be the number of PMUs attacked and an upper bound on the attack angle is considered to potentially account for the attacker’s limited capability to shift the phasor angle:^{2}^{2}2For vectors , notation means , .
(14a)  
subject to  (14b)  
(14c)  
(14d) 
where represents entrywise multiplication.
The two optimization variables in (14) are and . Constraint (14b) imposes bounds on that account for the maximum phase shift induced by the attack. The bounds can be set empirically based on studies such as the ones listed in Table II. Constraint (14c) expresses the binary nature of and represents that an attack can only happen on buses where PMUs are installed. The number of PMUs attacked in the network is captured via constraint (14d). The value of is realistically set to a small number, because it is unlikely that multiple GPS spoofers simultaneously attack at different geographical locations.
Problem (14) is a mixed integer program due to the binary
and thus hard to solve. Furthermore, even if the binary variables are relaxed to intervals
, the resulting problem is still nonconvex due to presence of sinusoids in the objective function.Optimization problem (14) yields the bus locations that if attacked induce the largest bias. The inputs are , , , and . The vector represents the buses with installed PMUs. It is worth emphasizing that optimization problem (14) is solved by the system operator, with the purpose of identifying the attacks that can potentially cause the maximum deviation of the state estimate from its true value. Optimization problem (14) depends on the voltage profile . The voltage profile is typically determined as a result of an optimal power flow routine, which considers the network demand and optimizes economical objectives, subject to reliability constraints (such as line thermal limits). The system operator can thus solve (14) for different voltage profiles , as the nodal injections and loads vary. Section IVD details the solution approach for (14).
IvD Solution Approach
Problem (14) can be optimally solved by enumerating all possible combinations of vector satisfying constraints (14c) and (14d), and solving for each vector the following optimization problem with variable :
(15)  
subject to 
Recall that the bias is a nonconvex function of . We use MATLAB’s nonlinear programming solver fmincon [35] to find a local maximum of (15). As the solution depends on the initial point, it is better to solve (15) for multiple initializations. The solution of (14) is given by and resulting in the largest objective value for (15). The optimal procedure to solve (14) is given in Algorithm 1.
Algorithm 1 uses exhaustive search. It thus has substantial complexity for large networks especially when , because it considers all combinations of PMUs. An alternative approach is to use a greedy algorithm. Consider the case of . The main idea is to use Algorithm 1 with to identify the single most vulnerable PMU and corresponding attack angle, then fix the attack angle on that bus and search for a second PMU bus and attack angle.
Specifically, first Algorithm 1 is run with . Let be the resulting bus index, that is , and let be the corresponding optimal attack angle. The second step is to solve problem (14) with and additional constraints and . This approach effectively reduces the problem of considering all combinations of PMUs into two searches of PMU, which are much simpler. Algorithm 2 summarizes the steps of the greedy approach. The idea can be extended to larger values of . For example, if , Algorithm 2 can be applied to find the first two most vulnerable PMUs (denote them as and with corresponding angles and ). Then, problem (14) with and additional constraints , , , and is solved.
The chief reason why Algorithm 2 has potential for computational effectiveness relative to Algorithm 1 is the complexity of the nonconvex optimization problem (15) that needs to be solved in every iteration of these algorithms. In particular, problem (15) always has one optimization variable for all cases that needs to be solved as Algorithm 2 runs. In contrast, problem (15) always has optimization variables for all cases that needs to be solved as Algorithm 1 runs. Taking advantage of the fact that highly efficient algorithms exist for onedimensional minimization, the reduction in the number of optimization variables plays a significant role.
V State Estimation and Attack Reconstrucion
This section develops an algorithm to jointly estimate the system state and identify the phase shifts for attacked buses.
The network operator has access to the measurement vectors at all buses on which PMUs are installed. Using the model in (5), the problem of jointly estimating and amounts to the following nonlinear least squares problem with variables and , where for clarity, the dependency of on is denoted explicitly:
(16) 
The previous problem is nonlinear due to the trigonometric functions present in the definition of (cf. Section IIIB). To alleviate this nonlinearity, we change the variable by introducing a new optimization variable . In order to be able to uniquely recover , the constraint is added. The resulting problem is a constrained bilinear least squares problem with optimization variables and :
(17a)  
subject to  (17b) 
where the objective function is written as
and is a block diagonal matrix that includes blocks of the term ; the variables and are used interchangeably.
Problem (17) is nonconvex and thus challenging. The nonconvexity arises because the objective is bilinear in the variables and and also because of the quadratic equality constraint (17b). But problem (17) can be efficiently tackled via an alternating minimization (AM) algorithm as explained in the sequel.
In general, the AM algorithm is applicable to minimization with respect to two groups of variables, in this case and . In the first step, minimization with respect to the first group is performed, by assuming the second group is kept fixed. The second step consists of minimization with respect to the second group of variables upon substituting the updated values for the first group of variables. These two steps are repeated until convergence. The convergence criterion is CurrObj PrevObjCurrObj Tolerance , where PrevObj and CurrObj represent the objective function values (17a) before and after the update. The algorithm is initialized by setting for all .
The AM algorithm has the attractive feature that it decreases the objective (17a) in every iteration. Interestingly, the minimizations in this algorithm can be solved in closed form, as shown next.
Va Minimization With Respect to the State
The minimization of (17) with respect to amounts to unconstrained least squares with solution
(18) 
where it is assumed that the matrix is nonsingular.
VB Minimization With Respect to the Attack Angle
The minimization in (17) with respect to takes the following equivalent form:
(19a)  
subject to  (19b) 
where is the th row of and is defined as
(20) 
Problem (19) is nonconvex due to the quadratic equality constraint (19b
). Interestingly, it is possible to solve this problem in closed form. To facilitate the solution, consider the eigenvalue decomposition (EVD) of
given as , where is orthonormal and is a diagonal matrix of nonnegative eigenvalues and . Define further . Notice that and are readily obtained with knowledge of the current iterate and the measurement . The next theorem characterizes the solution to (19).Theorem 3
Proof:
An optimal Lagrange multiplier always exists for (19), as the linear independence constraint qualification holds [36, Sec. 3.1]. The Lagrangian function of (19) is
(23) 
The optimality condition that yields as a function of is given by , which yields
(24) 
Assuming the invertibility of —which will be discussed shortly—(21) is obtained.
To find , substitute (21) into the constraint . This yields the equation where
The function can be written in the following simpler form if is substituted by its EVD :
(25) 
The function is strictly decreasing in in interval . Moreover, the limit of as approaches or is respectively or , therefore, the equation has one real root in the interval .
Eq. (22) is quartic in , that is, it yields four values of as solutions. The solutions of a quartic equation are completely characterized and can be routinely computed; see e.g., [37]. In addition, Theorem 3 ensures that at least one real solution exists, i.e., a situation where (22) has only complex roots never arises. The optimal value of is obtained by substituting each value of first in (21), and then picking the one that yields the smallest objective in (19). Algorithm 3 summarizes the steps to solve (17).
The previous procedure works seamlessly, unless is not invertible, which is unlikely in practice, as explained next. In particular, the following theorem characterizes the invertibility of .
Theorem 4
The matrix in (21) is invertible if or .
Proof:
We will prove the claim by contradiction. Suppose that the matrix is not invertible. Then, the optimal Lagrange multiplier is or .
In practice, it is unlikely that or . The reason is that depends on the noisy (recall that ), and thus the noise must have very particular values in order to yield or . In our numerical tests, we did not encounter any noninvertibility issue.
Remark 2 (Relationship to works on imperfect synchronization)
The modeling of GPS spoofing is related to imperfect PMU synchronization [38], and more broadly to asynchronous sampling in networked sensing systems [39]. But the GPS spoofing attack is different than imperfect PMU synchronization. In the latter, the measurement delay typically translates to less than for Hz frequency; whereas, in GPS spoofing, the attacker can potentially cause a much larger phase shift as documented in Table I. The small phase shift implies that approximations and can be invoked to render linear. These approximations are used in [38]. Problem (16) is more complicated when the attack angles are not assumed small. This section developed tractable algorithms for its solution based on the reformulation (17) that includes the nonconvex constraint (17b).
The work in [39]
includes a general asynchronicity model, but relies upon computing the Fourier transform of the underlying signal to be estimated. As such, it is necessary to acquire a number of signal samples across time. On the other hand, a single set of measurements from all PMUs is sufficient for the present work. The work
[39] also develops an alternating minimization algorithm. It performs full optimization with respect to one set of variables, but takes one step of gradient descent with respect to the other set. On the other hand, the present work performs full minimizations with respect to each set of optimization variables, which becomes possible by exploiting Lagrangian duality and the structure of the particular problem at hand.VC Simplifications Under Diagonal Covariance
Suppose that the covariance
is a diagonal matrix with equal variance for the real and imaginary parts of the voltage and likewise for the real and imaginary parts of every current on the
lines [30, 7, 38]. The covariance is thus assumed to have diagonal entries , , where(27) 
The computation of is very simple in this case, as the following theorem describes.
Theorem 5
The closedform solution for when the covariance has the structure given by (27) is
(28) 
Proof:
Under (27), the matrix is written as where is given by
(29) 
Comments
There are no comments yet.