Journal of Modern Power Systems and Clean Energy

ISSN 2196-5625 CN 32-1884/TK

网刊加载中。。。

使用Chrome浏览器效果最佳,继续浏览,你可能不会看到最佳的展示效果,

确定继续浏览么?

复制成功,请在其他浏览器进行阅读

Multi-period Power System State Estimation with PMUs Under GPS Spoofing Attacks  PDF

  • Paresh Risbud
  • Nikolaos Gatsis
  • Ahmad Taha
Department of Electrical and Computer Engineering, University of Texas at San Antonio, San Antonio, USA

Updated:2020-07-10

DOI:10.35833/MPCE.2020.000125

  • Full Text
  • Figs & Tabs
  • References
  • Authors
  • About
OUTLINE

Abstract

This paper introduces a dynamic network model together with a phasor measurement unit (PMU) measurement model suitable for power system state estimation under spoofing attacks on the global positioning system (GPS) receivers of PMUs. The spoofing attacks may introduce time-varying phase offsets in the affected PMU measurements. An algorithm is developed to jointly estimate the state of the network, which amounts to the nodal voltages in rectangular coordinates, as well as the time-varying attacks. The algorithm features closed-form updates. The effectiveness of the algorithm is verified on the standard IEEE transmission networks. It is numerically shown that the estimation performance is improved when the dynamic network model is accounted for compared with a previously reported static approach.

I. Introduction

THE smart grid is a massive cyber-physical system (CPS) spanning continents. The cyber part consists of a computerized network including two-way digital communication between devices (e.g., voltage and current sensors, power meters) and the network operation center. The physical part is formed by the power grid itself (generation, transmission, and distribution), which can be as large as a continent. Phasor measurement units (PMUs) equipped with global positioning system (GPS) receivers are installed ubiquitously and replace traditional sensors of the supervisory control and data acquisition (SCADA) systems. PMUs use GPS receivers to accurately time stamp and synchronize measurements. The main advantage of PMUs over SCADA is the higher sampling rate, which enables the operator to perform real-time wide area monitoring, protection and control (WAMPAC).

Sensors such as PMUs that measure nodal voltages, among other quantities, cannot be installed on each bus in the network. Thus, state estimation (SE) routines are performed to gain the visibility of the network. According to [

1], scenario WAMPAC.12: GPS Time Signal Compromise highlights the vulnerability of PMUs to GPS spoofing attacks. In particular, GPS spoofing is a class of cyber attacks in which the attacker mimics the GPS signals in order to alter the GPS time estimated by the GPS receiver of the PMU [2]. This attack is also called time synchronization attack (TSA), since it induces erroneous time stamps, thereby inducing wrong phase angle in the PMU measurements.

Various threat scenarios, security gaps, and mitigation strategies in relation to TSAs have been assessed in [

3] and [4]. Specific selected works representing effects and mitigation of GPS spoofing attacks in power systems include [5]-[20] and these are reviewed next.

A set of studies has investigated the operation at the GPS receiver level and developed methods to detect GPS spoofing and if possible, to provide accurate timing to the PMU. Specifically, [

5] provides a multi-layered approach for reliable GPS-based timing for PMUs against jamming and spoofing. Reference [6] leverages the static nature of PMUs to provide refined signal tracking techniques at the GPS receiver robustifying its operation against jamming and spoofing. Reference [7] develops a multilateration scheme to detect the malicious source of GPS signals by utilizing characteristics of multiple sensors in the grid. Detection techniques against GPS spoofing attacks in power grids based on the visible satellites and the statistics of the receiver clock are implemented in [8]. In an effort to detect attacked PMUs, [9] introduces a trustworthiness measure for PMUs that builds upon antenna- and signal-based techniques at the receiver.

The impacts of GPS spoofing attacks on power grid operation have also been analyzed. In addition to demonstrating the experimental feasibility of spoofing, a false generator trip in a synchrophasor-based automatic control scheme is exhibited in [

10]. A simulation-based approach to studing the effect of GPS spoofing on false alarm and missed generation scenarios in the smart grid is the theme of [11]. An optimization problem to demonstrate the feasibility of GPS spoofing attack is formulated in [12], where the objective is to maximize the PMU clock offset before and after the attack. The impact of GPS spoofing on voltage stability monitoring, fault detection, and event location is detailed in [13]. The effect on synchrophasor assisted load shedding are reported in [14].

More recently, power grid operation such as SE have been upgraded to account for GPS spoofing attacks on PMUs. Specifically, distributed estimation of power system oscillations under GPS spoofing attacks is developed in [

15]. Power system SE resilient to a spoofing attack on a single PMU is furnished in [16], with further computational refinements presented in [17]. In our previous work [18], an alternating minimization (AM) algorithm to jointly perform SE and attack phase angle reconstruction is developed. A sparse error correction for PMU data under GPS spoofing attacks is presented in [19]. A novel GPS spoofing identification and correction algorithm for SE in unbalanced distribution grids is developed in [20]. References [16]-[20] pertain to a static SE setup, whereby a set of measurements from a particular time instant is processed to produce a state estimate corresponding to that time instant. A dynamic SE model under GPS spoofing based on the swing equation is introduced in [21], where the presence of attack is detected using a generalized likelihood ratio test.

The effect of GPS spoofing on PMU measurements bears some resemblance to the problem of imperfect PMU synchronization [

22]. The similarity amounts to the fact that imperfect synchronization induces a time offset in the PMU measurement that in effect translates to a phase offset in the measured voltage or current phasor. The chief difference is that the phase offset is typically much smaller in case of imperfect synchronization, as opposed to GPS spoofing. The latter condition enables the use of small angle approximations [23]-[25] which are not generally valid under spoofing attacks. When dynamic state models with imperfect synchronization are considered, it is assumed that the availability of GPS restores the synchronization [26], which is a situation that is not tenable under GPS spoofing.

Based on a dynamic model for the power system state, this paper develops an SE algorithm that furnishes the phase offsets induced by GPS spoofing across time, in addition to the power system state. There are two approaches to dynamic SE depending on the definition of the state vector [

27], and representative references specifically featuring PMU measurements are mentioned next. The first approach features bus voltages as state variables, and the corresponding simplest dynamic model boils down to a random walk [28]-[30]. In the second approach, dynamic models of generators are considered, wherein rotor angles and rotor speeds are the state variables, following the swing equation [21], [31]-[33]. It is worth noting that when the states are the voltages expressed in rectangular coordinates and PMU measurements are used, SE can be performed efficiently by (linear) weighted least squares or the Kalman filter, whose performance is thoroughly analyzed in [28]. Depending on the variability of power injections, it is also possible to partition the network into three areas, i.e., steady, quasi-steady, and fluctuant, and perform a multi-time scale SE [34].

This paper adopts the random walk model for the network voltage states [

28]. The first contribution is deriving a state and measurement model that explicitly models the GPS spoofing attack angles. The attack may be time varying and start at an unknown time, but these characteristics are considered in the model. The multi-period SE is formulated as a bi-criterion optimization problem, where the measurement and dynamic state equations are contributing to the overall objective function. The states and attack angles are optimization variables and the resulting problem is nonconvex. The second contribution is developing an AM algorithm to compute the solution of the SE problem. The algorithm features closed-form updates and extends the one in [18] to the dynamic case. The effectiveness of the proposed algorithm is demonstrated on standard IEEE transmission networks. Two types of realistic attacks are applied [35], namely, one that shows up suddenly (step attack), and the other that is ramping. The third contribution is numerically demonstrating that the performance of the multi-period SE is superior to prior work that only considers the static nature of voltages, compared to the presented formulations that consider the voltage dynamics.

The remainder of this paper is organized as follows. Section II details the system and attack models. The SE problem is formulated in Section III. The iterative solver is developed in Section IV. Numerical tests and conclusions are presented in Sections V and VI, respectively.

II. System and Attack Models

This section describes the multi-period state and measurement model with and without TSAs on PMU measurements.

A. Attack-free Dynamic Model

Consider a power network with Nb buses connected via Nl transmission lines. Let 𝒩n be the set of buses connected to bus n, and define Ln=|𝒩n| as the number of lines connected to bus n. The nodal voltage at bus n and discrete time period k is written in complex, rectangular, and polar forms as Vn,k=Vn,k,r+jVn,k,i=|Vn,k|ejθn,k, where Vn,k,r and Vn,k,i are the real and imaginary parts of the nodal voltage, respectively; and θn,k is the phase of the nodal voltage. The dynamic equation describing the state evolution over discrete time periods indexed by k follows a random walk model:

vk=vk-1+wkk=1,2,,K (1)

where vk=[vk,rT,vk,iT]TR2Nb×1 is the system state vector at time k, and vectors vk,r and vk,i collect the real and imaginary parts of nodal voltages Vn,k,r and Vn,k,i for n=1,2,,Nb, respectively; wk is the state noise, which follows a Gaussian distribution, that is, wk𝒩 (0,Qk), and Qk is a positive definite covariance matrix at time k; and K is a generic terminal time through which SE is to be performed. The duration of time period k is given by the sampling period Ts. The typical sampling frequency for PMUs ranges between 30 and 120 samples per second, yielding a sampling period Ts of approximately 8-33 ms. The model in (1) is appropriate for the systems where the network dynamics are slow enough compared to the duration of the sampling period [

28].

PMUs are installed on the selected buses of the network, and an is a binary indicator equal to 1 if a PMU is installed at bus n and 0 otherwise. Vector a collects an for n=1,2,,Nb. The set of buses where PMUs are installed is denoted by 𝒩PMU={n{1,2,,Nb}an=1}={n1,n2,,nP}, where P is the number of PMUs installed in the network.

A PMU installed at bus n measures, for all time periods k, the complex nodal voltage as well as the complex currents on all lines that bus n is connected to. This collection of measured quantities (in rectangular coordinates) at bus n for time k is concatenated in a vector zn,kR2(1+Ln). To make the notation more compact, define Mn=2(1+Ln) as the number of distinct real quantities measured by the PMU at bus n for time k.

It is convenient for subsequent developments to consider the noiseless version of zn,k, which is denoted by zn,ktrueRMn:

zn,ktrue=Vn,k,rVn,k,i{Inp,k,r}p𝒩n{Inp,k,i}p𝒩n=|Vn,k|cos(θn,k)|Vn,k|sin(θn,k){|Inp,k|cos(θInp,k)}p𝒩n{|Inp,k|sin(θInp,k)}p𝒩n (2)

where Inp,k,r and Inp,k,i are the real and imaginary parts of the complex current flowing on line (n,p) for time k, respectively; and |Inp,k| and θInp,k are the corresponding current magnitude and phase, respectively.

To summarize, the noiseless quantities measured at bus n𝒩PMU, for discrete time k=1,2,,K, comprise of the real and imaginary parts of the nodal voltage, appended by the real and imaginary parts of the complex currents injected to all lines connected to bus n for time k. Using the bus admittance matrix of the network, zn,ktrue can be written as a linear function of the system state vk as zn,ktrue=Hnvk. The corresponding noisy measurement equation evolving over k discrete time periods is given as:

zn,k=Hnvk+εn,kn𝒩PMU,k=1,2,,K (3)

where εn,k𝒩 (0,Σn,k), and Σn,k is a positive definite measurement noise covariance matrix at time k; and the construction of HnRMn×2Nb is provided in [

36], [37].

B. Attack-aware Dynamic Model

A TSA at node n introduces a generally time-varying delay denoted by δn(t) to all measurements, i.e., voltages and currents, captured by the PMU. To introduce a mathematical model of the attack, let vn(t) be the instantaneous voltage of bus n, and let Vn(t) be the corresponding phasor at continuous time t. The phasor is time-dependent to allow for description of a system with time-varying state. The sampled phasor at period k is thus Vn,k=Vn(kTs). The attacked instantaneous nodal voltage is written as:

vn(t+δn(t))=Re{2Vn(t+δn(t))ej2πf(t+δn(t))} (4)

where f is the system operation frequency; and Re{} is the real part operator. Similar expressions can be written for the line currents. It is worth emphasizing that the same delay δn(t) is introduced across the measurements captured by the PMU at bus n (entries of zn,k). The reason is that GPS spoofing affects the time estimated by the GPS receiver of the PMU, which subsequently affects the time stamp of all measurements of the PMU. We also introduce the sampled version of δn(t), denoted by δn,k=δn(kTs), and the corresponding attack angle Δθn,k=2πfδn,k.

The objective is to formulate a time-varying measurement model that relates the attacked measurements with the network state and the attack in period k. To this end, it is assumed that δn(t)Ts. This assumption is valid as experimentally demonstrated. Other realistic attacks reported in [

10], [12], [13] depict attack angles Δθn,k of 70°, 52° and 60°, respectively. The corresponding delay of the attack is 3.2 ms, 2.4 ms, and 2.8 ms, respectively, which is indeed smaller than Ts. The following theorem characterizes the relationship between the measurement vector, the true state, and the attack at time k.

Theorem 1   The TSA measurement equation at discrete period k is given as:

zn,katk=Γn,kHnvk+εn,k (5)

where zn,katkRMn is the attacked measurement vector; and Γn,kRMn×Mn is a block diagonal matrix consisting of 1+Ln blocks where each block is a 2×2 matrix cos(Δθn,k)-sin(Δθn,k)sin(Δθn,k)cos(Δθn,k).

Proof   The attacked instantaneous voltage at t=kTs is

vn(kTs+δn(kTs))=Re{2Vn(kTs+δn(kTs))ej2πf(kTs+δn(kTs))} (6)

Considering that δn(t)Ts and the network dynamics evolve slowly with respect to the sampling period Ts, the voltage phasor at t+δn(t) can be approximated as Vn(t+δn(t))Vn(t) [

21]. Invoking the latter into (6), it follows that

vn(kTs+δn(kTs))Re{2Vn(kTs)ej2πf(kTs+δn,k)}=Re{2Vn,kej2πfδn,kej2πf(kTs)} (7)

Equation (7) reveals that the measured phasor is

Vn,kej2πfδn,k=|Vn,k|ejθn,kejΔθn,k=|Vn,k|ej(θn,k+Δθn,k) (8)

Extracting the real and imaginary parts of the latter, and repeating for the current measurements, the noisy attacked PMU measurement at bus n for time period k is given by (cf. (2)):

zn,katk=|Vn,k|cos(θn,k+Δθn,k)|Vn,k|sin(θn,k+Δθn,k){|Inp,k|cos(θInp,k+Δθn,k)}p𝒩n{|Inp,k|sin(θInp,k+Δθn,k)}p𝒩n+εn,k (9)

Upon introducing the trigonometric identities cos(a+b)=cos(a)cos(b)-sin(a)sin(b) and sin(a+b)=sin(a)cos(b)+cos(a)sin(b) into (9) and combining with (2) and zn,ktrue=Hnvk, the measurement model in (5) is obtained.

The assumed relationship between δn(t) and Ts enables the approximation Vn(t+δn(t))Vn(t). Thus, when δn(t) is smaller, but not significantly smaller, than Ts, the validity of the assumption and the resulting approximation are eventually determined by the extent in which the network dynamics are slow compared to the duration of the sampling period.

The model in (5) expresses the attacked PMU measurement zn,katk in terms of the spoofing attack Γn,k, the system Hn, and the state vk at time k. This model is leveraged in the next section to formulate the SE problem with spoofed PMU measurements.

III. Multi-period SE

This section presents the joint SE and attack angle reconstruction formulation. Both measurement and state equations are considered in a multi-objective optimization. Let vR2NbK collect vk for k=1,2,,K. Likewise, vector Δθ collects Δθn,k for n𝒩PMU and k=1,2,,K. The optimization is formulated as follows:

(v̂,Δθ̂)=argminv,Δθ(J1(v,Δθ)+J2(v)) (10)

where J1 represents the nonlinear weighted least squares problem of estimating v and Δθ based on the measurement equation (5), as given in (11); and J2 represents the weighted least squares problem of estimating v from the state equation (1), as given in (12).

J1=12n=1Nbk=1Kan||zn,katk-Γn,k(Δθn,k)Hnvk||Σn,k-12 (11)
J2=12k=1K||vk-vk-1||Qk-12 (12)

where the norm notation ||x||P2=xTPx is used. The initial state vector v0R2Nb is considered known, which is an assumption akin to typical multi-period estimation setups that rely on, for example, the Kalman filter.

Problem (10) is nonconvex due to the bilinear term Γn,k(Δθn,k)Hnvk and the sinusoidal dependence of Γn,k on Δθ in J1. Following [

18], a reformulation of (10) that can lead to an efficient solution algorithm is pursued in the sequel. Specifically, a change of variable is introduced as follows:

γn,k=[γn,k,1γn,k,2]T=[cos(Δθn,k)sin(Δθn,k)]T (13)

where γn,kR2×1; n𝒩PMU; and k=1,2,,K. The variable Δθn,k is eliminated, and the matrix Γn,k has blocks of the form γn,k,1-γn,k,2γn,k,2γn,k,1. The variables Γn,k and γn,k are thus used interchangeably.

In order to uniquely map the vector γn,k back to an angle Δθn,k, it is necessary and sufficient to impose the constraint γn,kTγn,k=1. Define γ as the vector including γn,k for all n𝒩PMU, k=1,2,,K. Then, problem (10) becomes equivalent to the following:

(v̂,γ̂)=argminv,γ(J1(v,γ)+J2(v)) (14)

s.t.

γn,kTγn,k=1n𝒩PMU,k=1,2,,K (15)

In order to facilitate the development of an algorithm for the solution of (3), it is supposed that the following condition holds for the measurement model in (3).

Assumption 1 The matrix [Hn1T,Hn2T,,HnPT]T is full column-rank.

This assumption pertains to the observation matrix corresponding to all PMU measurements. It ensures that under normal operation (i.e., no spoofing), the power network is observable and the SE problem has a unique solution [

38]. This condition is readily satisfied with proper placement of PMUs, a topic that has been well researched in the literature.

The transformation of the original unconstrained nonlinear problem (10) into the nonconvex quadratically constrained quadratic program (3) enables the development of an iterative solver with closed-form updates. This is the theme of the next section.

IV. Estimation Algorithm

This section develops an AM algorithm, to jointly solve (15) for v and γ. The algorithm minimizes two sets of variables one after the other. In the first step, the objective is minimized with respect to one set of variables while treating the second set as constant. In the second step, the minimization occurs with respect to the second set of variables upon substituting the updated values of the first set of variables. In this instance, the vectors v and γ constitute the two sets of variables.

The procedure is repeated until convergence. The initialization step includes γn,k=[1,0]T for all n and k. The two minimizations can be performed in closed form, and the related updates are described in the sequel.

A. Minimization with Respect to State

In order to derive the update for v, the objectives J1 and J2 are re-written as explicit functions of the vector v instead of vk. Specifically, J1 is written as:

J1(v,γ)=12n=1Nbk=1Kan||zn,katk-Γn,kHnBkv||Σn,k-12 (16)

where BkR2Nb×2NbK is a matrix such that vk=Bkv. In particular, Bk is constructed as a block diagonal matrix with the 2Nb×2Nb identity matrix in the kth block.

Similarly, J2 can be written as:

J2(v)=12||Ev-v0||Q1-12+k=2K||Akv||Qk-12 (17)

where ER2Nb×2NbK and AkR2Nb×2NbK (k=2,3,,K) are matrices such that v1=Ev and vk-vk-1=Akv. Specifically, matrix E includes the identity matrix in the top 2Nb×2Nb diagonal block and is zero otherwise. Matrix Ak is constructed as a fat matrix with negative identity matrix at the (k-1)th block and identity matrix at the kth block.

The minimization with respect to v is an unconstrained minimization with a convex quadratic objective function. The solution is obtained by solving the first-order optimality condition vJ(v)=0 and is given as:

v̂AM=M-1n=1Nbk=1Kan(Γn,kHnBk)TΣn,k-1zn,katk+ETQ1-1v0 (18)
M=ETQ1-1E+k=2KAkTQk-1Ak+n=1Nbk=1Kan(Γn,kHnBk)TΣn,k-1Γn,kHnBk (19)

The following result asserts that matrix M is indeed invertible.

Theorem 2   Given that Assumption 1 holds and the entries in γ satisfy (15), then matrix M is invertible.

Proof   It follows by the structure of M and the positive definiteness of Qk and Σn,k that M is positive semidefinite. We will show that the third term comprising M is full-rank, and therefore, M is invertible. The third term in M can be written as shown in (20).

k=1KBkTHn1Hn2HnPTΓn1,k000Γn2,k000ΓnP,kTΣn1,k-1000Σn2,k-1000ΣnP,k-1Γn1,k000Γn2,k000ΓnP,kHn1Hn2HnPBk (20)

The central matrix in (20) is full-rank because the convariance matrices Σn,k are positive definite. The matrices Γn,k are full-rank as long as the entries in γ satisfy (15). Finally, the matrix consisting of the Hni is full-rank due to Assumption 1. By successive application of the Sylvester inequality [

39], it follows that the matrix enclosed between BkT and Bk, which has size 2Nb×2Nb, has full rank for k=1,2,,K.

By the construction of Bk, the pre- and post-multiplication by BkT and Bk, respectively, generates a 2NbK×2NbK block diagonal matrix. The kth block in this larger matrix is full-rank, according to the previous explanation. Therefore, the matrix in (19) is full-rank.

The invertibility of M relies upon the γn,k values substituted in (2) satisfying (15). This condition is ensured by the next step.

B. Minimization with Respect to Attack Angles

In order to perform the minimization with respect to γ, only the first of the two objectives in (14) is relevant. Due to the separability of the objective J1(v,γ) (cf. (16)) and constraints per n and k, this minimization can be performed in parallel with respect to the variables pertaining to different PMUs and time periods. The resulting problem is stated for n𝒩PMU and k=1,2,,K as follows:

minγn,k[(zn,katk-Γn,kHnvk)TΣn,k-1(zn,katk-Γn,kHnvk)]s.t.  γn,kTγn,k=1 (21)

Considering that vk is not a variable in (21), the objective in (21) can be rearranged as follows:

minγn,k[(zn,katk-A˜n,kγn,k)TΣn,k-1(zn,katk-A˜n,kγn,k)]s.t.  γn,kTγn,k=1 (22)

where A˜n,k is given by

A˜n,k=hn,1Tvk-hn,2Tvkhn,2Tvkhn,1Tvkhn,Mn-1Tvk-hn,MnTvkhn,MnTvkhn,Mn-1Tvk (23)

where hn,iT is the ith row of Hn (i=1,2,,Mn).

Problem (21) is nonconvex due to the quadratic equality constraints. A problem of this form has been thoroughly analyzed in [

18], where a procedure to obtain a closed form solution based on Lagrangian duality is detailed. Furthermore, when the covariance matrices Σn,k have special structure, the solution of (21) is particularly easy to be obtained. As this is a practically relevant case, this result is stated next.

Theorem 3   Suppose that Σn,k is a diagonal matrix with equal variances for the real and imaginary parts corresponding to each measurement, that is, the diagonal entries of Σn,k satisfy the following for all n and k:

σn,k,12=σn,k,22σn,k,32=σn,k,42σn,k,Mn-12=σn,k,Mn2 (24)

Then, the solution of (21) is given as follows:

γ̂n,k=A˜n,kTΣn,k-1zn,katk||A˜n,kTΣn,k-1zn,katk||2 (25)

Note that the measurement covariance matrix is routinely assumed diagonal, and the particular structure mentioned in Theorem 3 can be found in [

23], [36], [40].

Algorithm 1 describes the AM steps for SE and attack reconstruction, incorporating the results of the present section. The convergence criterion checks whether at least one of two conditions hold: ||Ocurr-Oprev||2ϵ or |Ocurr-Oprev|/|Ocurr|ϵ, where Oprev and Ocurr are the values of the objective function (14) before and after the update, respectively; and ϵ is a predefined tolerance.

Algorithm 1  : SE and attack reconstruction

Result: state estimate v̂AM and attack angle Δθn,k, nNPMU, k=1,2,,K

Input: zn,katk

Initialization: solve (18) for v̂AM by setting γn,k=[1,0]T

repeat

for k=1,2,,K do

  for nNPMU do

   Find the corresponding γ̂n,k via (25)

  end

end

 Update v̂AM using (18)

until convergence or maximum iterations

Algorithm 1 minimizes the objective function after each update of the state and the attack angles. Because the objective function is lower bounded, the sequence of produced objective function values converges to a limit. This point may not be the global minimum, as the problem is nonconvex. Nevertheless, the numerical tests of Section V indicate very favorable SE performance.

C. Rolling Window Implementation

Algorithm 1 can be implemented in a rolling window fashion. Specifically, suppose that the state estimates for periods {1,2,,K} are available, and measurements for periods {K,K+1,,L} are obtained, where KL2K. Then, Algorithm 1 can be applied to provide state estimates for {L-K+1,L-K+2,,L}, where the value vL-K available from the previous window plays the role of the known v0 for the current window. Thus, the windows may be overlapping, and the process continues.

V. Numerical Tests

This section presents the state and attack angle estimation tests using the AM algorithm. The numerical tests are performed on the standard IEEE 14- and 118-bus systems. All network parameters are provided in case files case14.m and case118.m of MATPOWER [

41], from which Hn is computed. The PMU placement vector a for all test cases is obtained using the criterion in [36] based exclusively on the availability of PMU measurements. Table I lists the buses with installed PMUs for each network. Similar to [28], the state noise covariance Qk is diagonal, and its diagonal entries result from standard deviation of 0.001 p.u.. The measurement noise covariance Σn,k is also diagonal resulting from standard deviation of 0.001 p.u. for bus voltage and line current measurements. The measurement noise standard deviation is chosen such that PMU measurements do not violate the IEEE C37.118 standard [42]. In fact, with this choice of standard deviation, the majority of measurements incur a 0.1%-0.2% total variation error [28]. For simplicity, the state and measurement noise covariance matrices are assumed constant across all time periods k.

Table I Optimal PMU Location for IEEE Test Networks
Test case|NPMU|Bus number
IEEE 14-bus 6 2, 4, 6, 7, 10, 14
IEEE 118-bus 94 1-5, 7-19, 21-25, 27-36, 40, 43, 44, 46, 47, 48, 50, 51, 52, 53, 55-60, 64, 65, 66, 67, 68, 70, 71, 73, 75, 76, 77, 80-83, 85-90, 92, 94-104, 106-111, 113-118

The PMU sampling rate is set to 30 samples per second. The simulation considers a time horizon of 35 s. Two realistic attacks, namely, Type I (step) and Type II (ramp) are performed [

35]. The AM algorithm is first tested with attacks on the PMUs located at buses 14 and 7 of the IEEE 14- and 118-bus networks, respectively.

The Type I attack occurs suddenly at a particular time instant and remains constant thereafter. In this test, a Type I attack of 0.5787° appears at 30 s [

12]. It is customary in the GPS community to also express the attack in meters by multiplying the time offset of the attack by the speed of light c=3×108 m/s; the particular attack is thus 8000 m. The value of the attack is chosen so that the measurement exceeds the total variation error of 1% specified by the IEEE C37.118 standard [42]. The Type II attack changes gradually through time. For the attack to be successful, it cannot exceed the distance-equivalent velocity of 400 m/s [10]. Adhering to this value, the attack in this test starts at 10 s from value 0 m and gradually increases to 1000 m at 35 s.

To demonstrate the effectiveness of the algorithm, representative results pertaining to individual buses and the entire network are presented. The true and estimated voltage magnitudes and phase angles at bus 2 of the IEEE 14- and 118-bus networks across the time horizon are given in Fig. 1 and Fig. 2, respectively (bus 2 is not attacked). It is observed that estimated values closely follow the corresponding true values. Figures3 and 4 depict the true and estimated attack angles at buses 14 and 7, respectively, for the Type I and II attacks in the IEEE 14- and 118-bus networks. It is observed that the estimated attack angles closely follow their true values.

Fig. 1 True and estimated voltage magnitudes and phase angles at bus 2 of IEEE 14-bus network across time horizon. (a) Voltage magnitude. (b) Voltage phase angle.

Fig. 2 True and estimated voltage magnitudes and phase angles at bus 2 of IEEE 118-bus network across time horizon. (a) Voltage magnitude. (b) Voltage phase angle.

Fig. 3 True and estimated Type I and Type II attack angles at bus 14 in IEEE 14-bus network across time horizon. (a) Type I. (b) Type II.

Fig. 4 True and estimated Type I and Type II attack angles at bus 7 in IEEE 118-bus network across time horizon. (a) Type I. (b) Type II.

The relative voltage error is defined as ||v̂k-vk||2/||vk||2, where v̂k and vk are the state vector estimated with Algorithm 1 and its true value at period k, respectively. The relative voltage error for the IEEE 14- and 118-bus networks across the time horizon is depicted in Fig. 5. The relative error is of the order of 10-3 and 10-4 for the IEEE 14- and 118-bus networks, respectively. The error in attack angles is defined as ||Δθ̂k-Δθk||2/P, where Δθ̂k and Δθk are the attack angle vector resulting from Algorithm 1 and its true value at period k, respectively. Figure 6 depicts the attack angle errors for Type I and Type II attacks in the IEEE 14-bus network. Likewise, Fig. 7 shows the attack angle errors for both types of attacks in the IEEE 118-bus network.

Fig. 5 Relative voltage errors for IEEE 14- and 118-bus networks across time horizon. (a) IEEE 14-bus network. (b) IEEE 118-bus network.

Fig. 6 Type I and Type II attack angle errors for IEEE 14-bus network across time horizon. (a) Type I. (b) Type II.

Fig. 7 Type I and Type II attack angle errors for IEEE 118-bus network across time horizon. (a) Type I. (b) Type II.

Next, the algorithm is tested with attacks at buses 2, 4, 6, and 14 of the IEEE 14-bus network, and buses 7, 50, 60, and 80 of the IEEE 118-bus network. Only one bus is attacked in each run. The mean relative error is defined as the average of the relative error across the time horizon, that is, 1Kk=1K||v̂k-vk||2/||vk||2. The mean voltage phase angle is the average of the voltage phase angle of a bus across the time horizon, computed from the random walk model of the voltage state. Tables II and III list the mean voltage phase angle for the attacked bus and the resulting mean relative voltage error for the IEEE 14-bus network with Type I attacks of 0.5787° and 5°, respectively. The corresponding values for the IEEE 14-bus network with Type II attack are given in Table IV. Tables V-VII list the mean voltage phase angle for the attacked bus and the resulting mean relative voltage error for the IEEE 118-bus network under Type I and Type II attacks. The values for the mean voltage phase angle and mean relative voltage error in Tables II-VII are computed from single run of the algorithm. However, these values vary only slightly across multiple runs of the algorithm. The results in Tables II-VII reveal that the performance of the algorithm is not sensitive to the location of the attack, the voltage phase angle of the attacked bus, or the attack size.

Table II Mean Relative Voltage Error and Mean Voltage Phase Angle for IEEE 14-bus Network with Type I Attack of 0.5787°
Attacked bus numberMean voltage phase angle (°)Mean relative voltage error
2 -5.0255 0.00044492
4 -10.1426 0.00044767
6 -14.0238 0.00045789
14 -15.8632 0.00042576
Table III Mean Relative Voltage Error and Mean Voltage Phase Angle for IEEE 14-bus Network with Type I Attack of 5°
Attacked bus numberMean voltage phase angle (°)Mean relative voltage error
2 -5.3773 0.00043858
4 -10.1261 0.00044767
6 -14.7917 0.00044817
14 -17.5053 0.00044040
Table IV Mean Relative Voltage Error and Mean Voltage Phase Angle for IEEE 14-bus Network with Type II Attack
Attacked bus numberMean voltage phase angle (°)Mean relative voltage error
2 -5.3774 0.00044050
4 -8.9737 0.00040892
6 -14.5895 0.00044480
14 -16.1468 0.00043225

The quality of the state and attack angle estimation is compared to that of [

18], which performs estimation in a single period and does not account for the dynamic network model. Figures8,9,10 depict results from Algorithm 1 (single-period SE) and [18] (multi-period SE) plotted every 5th sample for the setup corresponding to the IEEE 14-bus network and attack on bus 14. The relative voltage error and attack angle error have been previously defined, while the state error norm depicted in Fig. 9 is defined as ||v̂k-vk||2. In the Type I attack error plots, the error in single-period SE is higher than that in multi-period SE before the attack occurs (at 30 s). Due to the step attack, the error in multi-period SE increases, but remains overall less than that in single-period SE. In Type II attack error plots, the error in single-period SE is greater than that in multi-period SE.

Fig. 8 Relative voltage error and attack angle error for Type I attack in IEEE 14-bus network across time horizon. (a) Relative voltage error. (b) Attack angle error.

Fig. 9 State error norms for Type I and Type II attacks in IEEE 14-bus network across time horizon. (a) Type I. (b) Type II.

Fig. 10 Relative voltage error and attack angle error for Type II attack in IEEE 14-bus network across time horizon. (a) Relative voltage error. (b) Attack angle error.

Table V Mean Relative Voltage Error and Mean Voltage Phase Angle for IEEE 118-bus Network with Type I Attack of 0.5787°
Attacked bus numberMean voltage phase angle (°)Mean relative voltage error
7 14.1334 0.00017028
50 20.8352 0.00016969
60 24.2264 0.00017338
80 29.0855 0.00016888
Table VI Mean Relative Voltage Error and Mean Voltage Phase Angle for IEEE 118-bus Network with Type I Attack of 5°
Attacked bus numberMean voltage phase angle (°)Mean relative voltage error
7 12.7141 0.00016947
50 18.9251 0.00017090
60 23.4543 0.00017410
80 27.3213 0.00017176
Table VII Mean Relative Voltage Error and Mean Voltage Phase Angle for IEEE 118-bus Network with Type II Attack
Attacked bus numberMean voltage phase angle (°)Mean relative voltage error
7 13.0204 0.00016892
50 20.8746 0.00017052
60 22.0586 0.00017233
80 27.8337 0.00017249

Finally, a comparison between Algorithm 1 and the Kalman filter is performed for the IEEE 14-bus network with Type I and Type II attacks at bus 14. The Kalman filter is run according to [

28]. Both the Algorithm 1 and the Kalman filter take the sequence of measurements zn,katk as input and produce the sequence of state estimates v̂k. Figure 11 depicts the resulting state error norm for the Kalman filter and Algorithm 1. It is observed that at the onset of the attack and thereafter, the Kalman filter yields larger state error norm than Algorithm 1 for both of the attack types. This is to be expected, as the Kalman filter is not designed to mitigate spoofing attacks.

Fig. 11 State error norms for Type I and Type II attacks in IEEE 14-bus network across time horizon using Algorithm 1 and Kalman filter. (a) Type I. (b) Type II.

VI. Conclusion

This paper puts forth a dynamic model which relates the measurement, state vector, and GPS spoofing attacks. The resulting nonconvex multi-period SE problem is solved by an AM algorithm that jointly estimates the state and reconstructs the attack. Two realistic attack scenarios are considered to validate the aforementioned algorithm on standard IEEE transmission networks. The numerical tests indicate that the estimation quality under GPS spoofing attacks is improved by considering the dynamic model of the network, as opposed to static estimation approaches.

The developed dynamic model and resulting algorithm are applicable to networks whose dynamics are slow enough compared to the sampling period, while relying on the condition that the delay induced by spoofing attack is smaller than the sampling period. It is an interesting direction to develop models for dynamic SE under spoofing attacks in more general setups. Furthermore, power networks may have available readings from SCADA and PMU systems. Developing algorithms for mitigating spoofing attacks in networks where both types of measurements are combined is an additional research direction.

References

1

T. Popovic, C. Blask, M. Carpenter et al., “Electric sector failure scenarios and impact analyses–version 3.0,” Electric Power Research Institute, Washington, Tech. Rep., Dec. 2015. [百度学术

2

A. Jafarnia-Jahromi, A. Broumandan, J. Nielsen et al., “GPS vulnerability to spoofing threats and a review of anti spoofing techniques,” International Journal of Navigation and Observation, vol. 2012, pp. 1-16, May 2012. [百度学术

3

D. Schmidt, K. Radke, S. Camtepe et al., “A survey and analysis of the GNSS spoofing threat and countermeasures,” ACM Computing Surveys, vol. 48, no. 4, pp. 1-31, May 2016. [百度学术

4

B. Moussa, M. Debbabi, and C. Assi, “Security assessment of time synchronization mechanisms for the smart grid,” IEEE Communications Surveys Tutorials, vol. 18, no. 3, pp. 1952-1973, Feb. 2016. [百度学术

5

L. Heng, J. J. Makela, A. D. Dominguez-Garcia et al., “Reliable GPS-based timing for power systems: a multi-layered multi-receiver architecture,” in Proceedings of Power and Energy Conference at Illinois, Champaign, USA, Feb. 2014, pp. 1-7. [百度学术

6

D. Chou, L. Heng, and G. Gao, “Robust GPS-based timing for phasor measurement units: a position-information-aided approach,” in Proceedings of 27th International Technical Meeting of the Satellite Division of the Institute of Navigation, Tampa, USA, Sept. 2014, pp. 1261-1269. [百度学术

7

D.-Y. Yu, A. Ranganathan, T. Locher et al., “Short paper: detection of GPS spoofing attacks in power grids,” in Proceedings of ACM Conference on Security and Privacy in Wireless & Mobile Networks, Oxford, United Kingdom, Jul. 2014, pp. 99-104. [百度学术

8

F. Zhu, A. Youssef, and W. Hamouda, “Detection techniques for data level spoofing in GPS-based phasor measurement units,” in Proceedings of International Conference on Selected Topics in Mobile Wireless Networking, Cairo, Egypt, Apr. 2016, pp. 1-8. [百度学术

9

Y. Fan, Z. Zhang, M. Trinkle et al., “A cross-layer defense mechanism against GPS spoofing attacks on PMUs in smart grids,” IEEE Transactions on Smart Grid, vol. 6, no. 6, pp. 2659-2668, Nov. 2015. [百度学术

10

D. P. Shepard, T. E. Humphreys, and A. A. Fansler, “Evaluation of the vulnerability of phasor measurement units to GPS spoofing attacks,” International Journal of Critical Infrastructure Protection, vol. 5, no. 3, pp. 146-153, Dec. 2012. [百度学术

11

I. Akkaya, E. A. Lee, and P. Derler, “Model-based evaluation of GPS spoofing attacks on power grid sensors,” in Proceedings of Workshop Modeling and Simulation of Cyber-Physical Energy Systems, Berkeley, USA, May 2013, pp. 1-6. [百度学术

12

X. Jiang, J. Zhang, B. J. Harding et al., “Spoofing GPS receiver clock offset of phasor measurement units,” IEEE Transactions on Power Systems, vol. 28, no. 3, pp. 3253-3262, Aug. 2013. [百度学术

13

Z. Zhang, S. Gong, A. D. Dimitrovski et al., “Time synchronization attack in smart grid: impact and analysis,” IEEE Transactions on Smart Grid, vol. 4, no. 1, pp. 87-98, Mar. 2013. [百度学术

14

J. G. Sreenath, S. Mangalwedekar, A. Meghwani et al., “Impact of GPS spoofing on synchrophasor assisted load shedding,” in Proceedings of IEEE PES General Meeting, Portland, USA, Aug. 2018, pp. 1-5. [百度学术

15

Y. Wang and J. P. Hespanha, “Distributed estimation of power system oscillation modes under attacks on GPS clocks,” IEEE Transactions on Instrumentation and Measurement, vol. 67, no. 7, pp. 1626-1637, Jul. 2018. [百度学术

16

X. Fan, L. Du, and D. Duan, “Synchrophasor data correction under GPS spoofing attack: a state estimation based approach,” IEEE Transactions on Smart Grid, vol. 9, no. 5, pp. 4538-4546, Sept. 2018. [百度学术

17

X. Fan, S. Pal, D. Duan et al., “Closed-form solution for synchrophasor data correction under GPS spoofing attack,” in Proceedings of IEEE PES General Meeting, Portland, USA, Aug. 2018, pp. 1-5. [百度学术

18

P. Risbud, N. Gatsis, and A. Taha, “Vulnerability analysis of smart grids to GPS spoofing,” IEEE Transactions on Smart Grid, vol. 10, no. 4, pp. 3535-3548, Jul. 2019. [百度学术

19

S. D. Silva, T. Hagan, J. Kim et al., “Sparse error correction for PMU data under GPS spoofing attacks,” in Proceedings of IEEE Global Conference on Signal and Information Processing, Anaheim, USA, Nov. 2018, pp. 902-906. [百度学术

20

Y. Zhang, J. Wang, and J. Liu, “Attack identification and correction for PMU GPS spoofing in unbalanced distribution systems,” IEEE Transactions on Smart Grid, vol. 11, no. 1, pp. 762-773, Jan. 2020. [百度学术

21

P. Pradhan, K. Nagananda, P. Venkitasubramaniam et al., “GPS spoofing attack characterization and detection in smart grids,” in Proceedings of IEEE Conference on Communications and Network Security, Philadelphia, USA, Oct. 2016, pp. 391-395. [百度学术

22

P. Castello, C. Muscas, P. Pegoraro et al., “Trustworthiness of PMU data in the presence of synchronization issues,” in Proceedings of IEEE International Instrumentation and Measurement Technology Conference, Houston, USA, May 2018, pp. 1-5. [百度学术

23

P. Yang, Z. Tan, A. Wiesel et al., “Power system state estimation using PMUs with imperfect synchronization,” IEEE Transactions on Power Systems, vol. 28, no. 4, pp. 4162-4172, Nov. 2013. [百度学术

24

J. Du, S. Ma, Y. Wu et al., “Distributed Bayesian hybrid power state estimation with PMU synchronization errors,” in Proceedings of IEEE Global Communications Conference, Austin, USA, Dec. 2014, pp. 3174-3179. [百度学术

25

J. A. Bazerque, U. Ribeiro, and J. Costa, “Synchronization of phasor measurement units and its error propagation to state estimators,” in Proceedings of IEEE PES Innovative Smart Grid Technologies Latin America, Montevideo, Uruguay, Oct. 2015, pp. 508-513. [百度学术

26

M. Todescato, R. Carli, L. Schenato et al., “PMUs clock desynchronization compensation for smart grid state estimation,” in Proceedings of IEEE Conference on Decision and Control, Melbourne, Australia, Dec. 2017, pp. 793-798. [百度学术

27

L. Hu, Z. Wang, I. Rahman et al., “A constrained optimization approach to dynamic state estimation for power systems including PMU and missing measurements,” IEEE Transactions on Control Systems Technology, vol. 24, no. 2, pp. 703-710, Mar. 2016. [百度学术

28

S. Sarri, L. Zanni, M. Popovic et al., “Performance assessment of linear state estimators using synchrophasor measurements,” IEEE Transactions on Instrumentation and Measurement, vol. 65, no. 3, pp. 535-548, Mar. 2016. [百度学术

29

M. B. D. C. Filho and J. C. S. de Souza, “Forecasting-aided state estimation–part I: panorama,” IEEE Transactions on Power Systems, vol. 24, no. 4, pp. 1667-1677, Nov. 2009. [百度学术

30

A. Monticelli, “Electric power system state estimation,” Proceedings of IEEE, vol. 88, no. 2, pp. 262-282, Feb. 2000. [百度学术

31

J. Zhang, G. Welch, G. Bishop et al., “A two-stage Kalman filter approach for robust and real-time power system state estimation,” IEEE Transactions on Sustainable Energy, vol. 5, no. 2, pp. 629-636, Apr. 2014. [百度学术

32

N. Zhou, D. Meng, Z. Huang et al., “Dynamic state estimation of a synchronous machine using PMU data: a comparative study,” IEEE Transactions on Smart Grid, vol. 6, no. 1, pp. 450-460, Jan. 2015. [百度学术

33

J. Zhao, M. Netto, and L. Mili, “A robust iterated extended Kalman filter for power system dynamic state estimation,” IEEE Transactions on Power Systems, vol. 32, no. 4, pp. 3205-3216, Jul. 2017. [百度学术

34

Y. Guo, B. Zhang, W. Wu et al., “Multi-time interval power system state estimation incorporating phasor measurements,” in Proceedings of IEEE PES General Meeting, Denver, USA, Jul. 2015, pp. 1-5. [百度学术

35

A. Khalajmehrabadi, N. Gatsis, D. Akopian et al., “Realtime rejection and mitigation of time synchronization attacks on the global positioning system,” IEEE Transactions on Industrial Electronics, vol. 65, no. 8, pp. 6425-6435, Aug. 2018. [百度学术

36

V. Kekatos, G. Giannakis, and B. Wollenberg, “Optimal placement of phasor measurement units via convex relaxation,” IEEE Transactions on Power Systems, vol. 27, no. 3, pp. 1521-1530, Aug. 2012. [百度学术

37

P. Risbud, N. Gatsis, and A. Taha, “Assessing power system state estimation accuracy with GPS-spoofed PMU measurements,” in Proceedings of 7th IEEE Conference on Innovative Smart Grid Technologies, Minneapoils, USA, Sept. 2016, pp. 1-5. [百度学术

38

A. Gomez-Exposito, A. J. Conejo, and C. Canizares, Electric Energy Systems: Analysis and Operation. Boca Raton: CRC Press, 2018. [百度学术

39

R. Horn and C. R. Johnson, Matrix Analysis, 2nd ed. Cambridge: Cambridge University Press, 2013. [百度学术

40

A. Gomez-Exposito, P. Rousseaux, C. Gomez-Quiles et al., “On the use of PMUs in power system state estimation,” in Proceedings of 17th Power System Computation Conference, Stockholm, Sweden, Aug. 2011. [百度学术

41

R. Zimmerman, C. Murillo-Saanchez, and R. Thomas, “MATPOWER: steady-state operations, planning, and analysis tools for power systems research and education,” IEEE Transactions on Power Systems, vol. 26, no. 1, pp. 12-19, Feb. 2011. [百度学术

42

IEEE Standard for Synchrophasor Measurements for Power Systems, IEEE Std C37.118.1-2011 (Revision of IEEE Std C37.118-2005), 2011. [百度学术