Vulnerability Assessment of N-1 Reliable Power Systems to False Data Injection Attacks

1 V ulnerability Assessment of N − 1 Reliable Po wer Systems to F alse Data Injection Attacks Zhigang Chu, Jiazi Zhang, Oliv er K osut, and Lalitha Sankar School of Electrical, Computer and Energy Engineering Arizona State Univ ersity Abstract —This paper studies the vulnerability of large-scale power systems to false data injection (FDI) attacks through their physical consequences. Prior work has shown that an attacker- defender bi-level linear pr ogram (ADBLP) can be used to determine the w orst-case consequences of FDI attacks aiming to maximize the physical power flow on a target line. Understanding the consequences of these attacks requires consideration of power system operations commonly used in practice, specifically real- time contingency analysis (R TCA) and security constrained eco- nomic dispatch (SCED). An ADBLP is formulated with detailed assumptions on attacker’ s knowledge, and a modified Benders’ decomposition algorithm is introduced to solve such an ADBLP . The vulnerability analysis results pr esented f or the synthetic T exas system with 2000 buses show that intelligent FDI attacks can cause post-contingency overflows. Index T erms —F alse data injection attack, cyber -security , vul- nerability of N − 1 reliable power system, bi-level optimization. I . I N T RO D U C T I O N The efficienc y and intelligence of modern electric power systems are increasing rapidly with integration of real-time monitoring, sensing, communication and data processing. This integration is accomplished via a c yber layer consisting of the supervisory control and data acquisition (SCAD A) system in conjunction with the energy management system (EMS). SCAD A monitors the physical system, collects measurements, and sends them to the control center . In the EMS, state esti- mation (SE) estimates the voltage magnitudes and angles from measurements. This estimate along with the subsequent data processing, optimization and communication, specifically real- time contingenc y analysis (R TCA) [1] and security constrained economic dispatch (SCED) [2], allo w for real-time control of the power systems. Howe ver , the integration of the cyber layer also increases the threat of cyber-attacks on power systems that could lead to sev ere physical consequences, as illustrated by the recent cyber -attack in Ukraine (see [3]). Therefore, it is crucial to dev elop techniques to detect and thw art potential attacks, which requires ev aluating system vulnerability to credible at- tacks. Assessing consequences of possible attacks is extremely instructiv e for system operators, and is important for secure power system operations. This paper focuses on f alse data injection (FDI) attacks, wherein a malicious attacker replaces a subset of SCAD A measurements (power flows and injections) with counterfeits. FDI attacks can be designed to target system states [4]–[6], system topology [7], [8], and energy markets [9]. Optimization problems have been proposed to design FDI attacks that aim to maximize line power flow [10], change locational marginal prices [11], or maximize operating cost [12]. Howe ver , the results hav e only been demonstrated for small systems, and do not include the ef fects of N − 1 reliability constraints. N − 1 reliable system operations typically inv olve R TCA to generate security constraints, and SCED to re-dispatch the generators with the security constraints in the most economic sense. In this paper, we focus on the worst-case FDI attacks that maximize power flow on a target line, b ut our goal is to ev aluate vulnerability of significantly larger systems ( i.e. thousands of buses) with A C SE, R TCA, and SCED. The authors of [10] introduce an unobservable FDI attack that re-distributes the loads by changing SCADA measure- ments, to trigger generation re-dispatches that result in phys- ical ov erflo w on a target line. The authors show that the worst-case attack can be found using an attacker-defender bi- lev el linear program (ADBLP), wherein the first lev el models the attacker’ s objectiv e and limitations, and the second level models the system response through DC optimal po wer flo w (OPF). In this paper , we formulate a similar ADBLP with the second lev el modeled as SCED instead of OPF , taking into account the security constraints generated by R TCA. W e also discuss in detail the required knowledge of an attacker to design such an unobservable attack. T echniques to solve ADBLPs with applications to power systems hav e been studied in [13], [14], b ut are limited to scenarios with the same objectiv e for both le vels, and hence, their techniques cannot be applied to either the problem in [10] or its generalization considered here for large power systems. An ADBLP can be reformulated as a mathematical program with equilibrium constraints (MPEC) [15] by replacing the second lev el by its Karush-Kuhn-T ucker (KKT) conditions. Howe ver , MPECs are non-conv ex and hard to solve efficiently in general [16]. The MPEC from the ADBLP can be further reformulated as a mixed-integer linear program (MILP) by rewriting the complementary slackness constraints as mixed- integer constraints. As the system size increases, this MILP be- comes harder to solve due to the increasing number of binary variables. In this paper, we introduce an ef ficient algorithm that utilized Benders’ decomposition to solve ADBLPs. The contributions of this paper are as follows: 1) Knowledge requirement for attacker to design unobserv- able attacks in the presence of R TCA and SCED; 2) Attack design ADBLP modeling SCED as the second lev el problem; 3) Modified Benders’ decomposition algorithm to solve AD- BLPs; 4) Simulations of the designed attacks on the synthetic T exas system to ev aluate the consequences of such attacks. The remainder of this paper is organized as follows. Sec. II describes the power system measurement model and unob- servable attack model. Sec. III discusses attacker’ s requirement to design and launch worst-case unobservable attacks. Sec. IV introduces an ADBLP to find worst-case attacks. Sec.V describes the modified Benders’ decomposition algorithm to solve any ADBLP . Sec. VI presents the simulation results on the T exas system. Sec. VII concludes the paper and considers directions of future work. I I . S Y S T E M A N D A T TAC K M O D E L A. EMS Operation In this paper , we consider an EMS with three core functions operating in the order of SE, R TCA, and SCED. The EMS operating structure is illustrated in Fig. 1. Power system measurement data collected by SCAD A are sent to SE, which estimates the complex voltages after eliminating noise and bad measurements. Giv en the generator set points, SE also estimates the load values in the system. Modern po wer systems typically require N − 1 reliability , i.e., the system must operate with no violations if a contingency occurs (one of the system components, generators or branches, is out of service). R TCA simulates one power flow under each contingency k . W e say a branch has a warning if its power flow is abov e a threshold τ but less than its limit, while a branch has a violation if its power flow exceeds its limit. Note that in base case, the limit is the long-term line rating, while in contingency case it is the short-term rating. Each warning and violation generates one line limit constraint to be modeled in SCED. In contingency cases, these constraints are called security constraints. SCED takes all these constraints and solves an optimization problem to determine the most economic generation dispatch that ensures N − 1 reliable operation. S CA D A S t a t e E s t i m a t or C o n t i n g e n c y A n a l y s i s D C A p p ro x i m a t i o n s S CE D R e a l - ti m e C o n ti n g e n c y A n a l y s i s G e n e r a l i z e d S C ED P h y s i c a l S y s t e m G e n e r a l i z e d S ta te Es ti m a to r G e n e ra t e S e c u r i t y Co n s t ra i n t s B a d D a t a ? B a d D a t a El i m i n a t o r Y e s No M e a s u re m e n t s G e n e ra t i o n D i s p a t c h EM S C o r e F u n c ti o n s Fig. 1: EMS operation with SE, R TCA, and SCED. B. Measurement Model W e model the power system with n b buses, n g generators, and n m measurements. The SCAD A system measurement model is giv en by z = h ( x ) + e (1) where z is the n m × 1 measurement vector; x is the 2 n b × 1 vector of bus voltage magnitudes and angles (states); h ( · ) is the non-linear relationship between measurements and states; e is the n m × 1 vector of measurement noise, whose entries are assumed to be jointly distributed as N (0, R ) where R = diag ( σ 2 1 , σ 2 2 , . . . , σ 2 n m ) . C. Unobservable Attack Model An n m × 1 measurement attack vector a = ¯ z − z is defined to be unobservable to the residual-based bad data detector (BDD) if a = h ( x + c ) − z , where ¯ z is the vector of the false measurements created by the attacker and c is the state attack vector [6]. Gi ven c , an attack subgraph S can be constructed as in [6], such that the non-zero entries of a are all within S . By modifying measurements only in S , the attacker can arbitrarily spoof the states of center buses (load buses corresponding to non-zero entries of c ) without detection. The attack causes the system estimated load to re-distribute between load buses within S , while the total load remain unchanged. I I I . A T TAC K E R A S S U M P T I O N S Among all unobservable FDI attacks, the most dangerous ones are those with serious physical consequences. In this paper , we focus on a class of unobservable attacks where the attacker maliciously changes the SCAD A measurements to maximize the power flo w on a target line, and possibly cause overflo w . The authors of [10] introduce an ADBLP to determine the worst-case unobservable line overflo w attack, wherein the first lev el models the attacker’ s objective and limitations, while the second level models the system response via DCOPF . Assuming the attacker has knowledge of (i) the complete network topology (including line parameters and ratings) and load information, and (ii) the cost, capacity , and operational status of all generators in the system, the authors show that unobserv able attacks found using this optimization successfully result in generation re-dispatches that cause line ov erflo ws on the IEEE R TS 24-b us system. Howe ver , modern po wer systems typically do not use DCOPF to re-dispatch the generation, but rather operates as outlined in Sec. II-A. An attacker who gains kno wledge of EMS operations has an adv antage to accurately predict the system response. In other words, if the attacker is able to perform the same R TCA and SCED as the system does, it can design attacks that maximize the consequences. This is a stronger assumption than that in [10], because in addition to having access to the database of the control center , now the attacker further kno ws the algorithms and assumptions used by the system. While this is a stronger requirement, it is valuable to understand how the system is resilient against such strong adversaries through this worst-case approach. In R TCA, the attacker needs to know the power flow algorithm used to get the same post-contingency flows on all lines, as well as the threshold τ as described in Sec. II, to determine the security constraints to be included in SCED. In SCED, the attacker should kno w how the system models the constraints, as dif ferent system operators may implement SCED differently . W e assume the attacker has full knowledge of R TCA and SCED implementation in the EMS, in particular: 1) Contingency ratings of the branches; 2) Loss handling method; 2 3) Ramp rates and reserve costs of all generators; 4) Reserve policy and requirements; 5) Criteria to determine which base case line limits are to be modeled. This can be the same threshold as τ in post- contingency case, but can also be dif ferent; 6) Branch flow calculation method in both base case and contingency case; 7) Load shedding policy and costs. While it seems unrealistic to gain such kno wledge, it is not entirely impossible, since such complex systems in volv e sophisticated (even nation-state) attackers that can exploit or hav e access to insider knowledge [3], [17]. Again, this is the worst-case assumptions, and therefore, resilience of the system to such w orst-case attacks can serve as an upper bound on risks to the system operations. I V . A D B L P T O F I N D W O R S T - C A S E A T TAC K In this section, we introduce an ADBLP similar to that in [10] to find worst-case line ov erflo w attacks. The first lev el models the attacker’ s objective and limitations, while the sec- ond level models the system response via SCED. W e focus on R TCA that simulates branch contingencies (excluding radial branches), and reports corresponding security constraints to SCED. Contingency k indicates that branch k is out of service. The attacker is assumed to be able to perform R TCA and pick a target line l to maximize its po wer flow when target contingency k t occurs, and possibly create overflo w . W ithout loss of generality , we assume the flow on l is positiv e; if it is not the case, its absolute value can be maximized. In the formulation belo w we assume the attacker aims to maximize post-contingency power flow on the target line, but the base case power flow can also be maximized. Since SCED is DC, the voltage magnitudes are all considered to be 1 p.u., and hence, c is an n b × 1 attack vector on the voltage angles. III. The ADBLP takes the following form: maximize c P l,k t − σ k c k 1 (2a) subject to P l,k t = O TDF l k t ( G B P ∗ G − P D ) (2b) k c k 1 ≤ N 1 (2c) − L S P D ≤ H c ≤ L S P D (2d) { P ∗ G } = arg  min P G ,R G ,P,P k C G ( P G ) + C R R G  (2e) subject to P n g g =1 P Gg = P n b i =1 P Di (2f) ¯ P = P 0 + PTDF ( G B ( P G − P G 0 ) + H c ) (2g) ¯ P k = P k 0 + OTDF k ( G B ( P G − P G 0 ) + H c ) (2h) + LODF k · PTDF k · H c, ∀ k − P max ≤ ¯ P ≤ P max (2i) − P k, max ≤ ¯ P k ≤ P k, max , ∀ k (2j) P G ≥ max { P G 0 − M G T h , P G, min } (2k) P G ≤ max { P G 0 + M G T h , P G, max } (2l) 0 ≤ R G ≤ M G T r (2m) P G + R G ≤ P G, max (2n) P n g g =1 R Gg ≥ P Gg + R Gg , ∀ g (2o) where the variables are: c attack vector , n b × 1 ; ¯ P , ¯ P k vectors of monitored line cyber power flows in base case and under contingency k , respectiv ely; P l,k t physical power flow on target line l under target contingency k t ; P G power output of generators, n g × 1 ; R G spinning reserve of the generators, n g × 1 ; and the parameters are: σ penalty of the l 1 -norm of attack vector c ; G B generators to buses connectivity matrix, n b × n g ; O TDF k outage transfer distrib ution factor matrix under contingency k ; O TDF l k l th row of O TDF k ; N 1 attack vector l 1 -norm limit; L S load shift factor , in percentage; H dependency matrix between power injection mea- surements and states, n b × n b ; P D vector of real loads, n b × 1 ; C G generation cost vector , n g × 1 ; C R reserve cost vector , n g × 1 ; P 0 , P k 0 vectors of pre-SCED monitored line power flows in base case and under contingency k , respec- tiv ely; P G 0 pre-SCED generator outputs, n g × 1 ; PTDF power transfer distribution factor matrix; PTDF k k th row of PTDF; LODF k line outage distribution factors of monitored lines under contingency k ; P max vector of base case line limits; P k, max vector of line limits under contingency k ; P G, min generation lower limits vector , n g × 1 ; P G, max generation upper limits vector , n g × 1 ; M G ramp rates of all generators, n g × 1 ; T h look-ahead time for one period SCED; T r time for spinning reserve requirement. Expression in (2a) captures the attacker’ s objective of max- imizing the po wer flow on line l under tar get contingency k t , and the penalty f actor σ is a small positive number to limit the attack size; constraint (2b) is the calculation of the power flow on line l under target contingency k t ; (2c) models the attack er’ s limited resources. Ideally , l 0 -norm should be used to precisely capture the sparsity of c , but for tractability reasons we use the l 1 -norm as a proxy . Constraint (2d) limits the percentage of load changes at each bus to av oid detection. SCED (2e)-(2o) models the system response to the attack. The objectiv e of the operator (2e) is to minimize the total cost, consisting of generation cost and reserve cost; constraint (2f) is the power balance equation; (2g) is the cyber power flow of the base case monitored lines. Note that this constraint is only modeled for the lines whose pre-SCED power flow is greater than the threshold τ , i.e., | P 0 /P max | ≥ τ . This is under the assumption that the line flows will not change dramatically after the SCED re-dispatch, due to the ramping 3 constraints of the generators. Similarly , (2h) is the cyber power flows on monitored lines under each contingenc y k , where | P k 0 /P k, max | ≥ τ . Here we assume the base case and contingency case monitoring thresholds are the same. In the right hand side of (2h), the first term is the pre-SCED post-contingency flows; the second term is the change of the flows as a result of re-dispatch and false loads; the third term represents the amount of po wer on the monitored lines resulting from the ef fect of false loads on the contingency line k , which is not considered in P k 0 . Constraints (2i) and (2j) are the line limits in base case and contingency case, respectiv ely . The activ e power limits in both base case and contingency cases, P max and P k, max , are approximated from the MV A ratings and reactive flows on the branches by P max = p S 2 max − [max( Q from , Q to )] 2 (3) P k, max = q S 2 k, max − [max( Q k, from , Q k, to )] 2 (4) where S max and S k, max are branch long-term and short-term ratings, respectively; Q from and Q to are the base case reactiv e branch flows at the "from" end and "to" end, respecti vely; Q k, from and Q k, to are those flows in contingency cases. Con- straints (2k) and (2l) are the ramping limits; (2m) is the reserve limit; (2n) is the generation limit. Though the R TCA does not simulate generator contingencies, in SCED it is required that when a generator is out, the reserves of all other generators are sufficient to cover the output of the lost generator . W e assume the SCED does not include a load shedding policy . V . M O D I FI E D B E N D E R S ’ D E C O M P O S I T I O N A L G O R I T H M T O S O LV E A D B L P S ADBLPs with different objectives in the two lev els are in general non-con ve x. The authors of [10] solve their ADBLP by replacing the second level defender’ s problem by its KKT conditions and then con vert the problem into an MILP , but this approach does not apply to large-scale systems due to the numerical difficulty brought on by large number of binary variables. T o the best of our knowledge, there are no existing techniques to solve large-scale ADBLPs efficiently . In this section, we introduce a modified Benders’ decomposition (MBD) algorithm to solve ADBLPs. Benders’ decomposition [18] is an iterati ve approach to solve linear programs in a distributed manner [19]. It is a popular technique to solve optimization problems of large size or with complicating variables. It is also effecti ve in solving complex optimization problems such as stochastic programs and mixed-integer linear programs. In Benders’ decomposition, an optimization prob- lem is decomposed into two sub-problems, wherein v ariables of each sub-problem are treated as constant in the other . The two sub-problems are solved iterativ ely until the solution con verges. Our MBD algorithm modifies the classic Benders’ decomposition algorithm to apply it on any ADBLP . An ADBLP takes the following form (dual v ariable of the defender’ s problem is in parentheses): minimize u c T 1 u + d T 1 v ∗ (5a) subject to A 1 u ≥ b 1 (5b) v ∗ = arg { min v d T 2 v } (5c) subject to A 2 u + A 3 v ≥ b 2 ( β ) (5d) where u and v are the attacker’ s and defender’ s decision variables, respectiv ely . The defender has no control on u , and hence, u in (5d) is treated as a constant in the defender’ s problem. The attacker does not directly control v , but it controls v ∗ by changing u , assuming it has knowledge of the defender’ s objectiv e and constraints. The attack optimization ADBLP (2) fits in the form of (5) where the attack vector c is represented by u and SCED variables P G , R G , P , and P k are represented by v . In the attacker’ s objecti ve function, c T 1 u represents the term − σ k c k 1 , and d T 1 v ∗ represents the term P l,k t in (2a). Equality constraints can be equi v alently written as two inequality constraints. For example, (2f) can be written as 1 T P G ≥ 1 T P D (6a) − 1 T P G ≥ − 1 T P D (6b) which fits the form of (5d). One can similarly map all the constraints in (2) to those in (5). The defender’ s problem (5c)–(5d), which represents the system response (SCED) to a fixed attack vector , has the following dual problem (note that u is treated as constant here since it is the fixed attack vector from the attacker’ s problem): maximize β β T ( b 2 − A 2 u ) (7a) subject to A T 3 β = d 2 (7b) β ≥ 0 . (7c) By weak duality [20], for any feasible primal/dual pair , the dual objectiv e value is always less than the primal one: β T ( b 2 − A 2 u ) ≤ d T 2 v . (8) Since the defender’ s problem is a linear program, it satisfies strong duality . That is, any feasible point ( v , β ) that satisfies β T ( b 2 − A 2 u ) ≥ d T 2 v (9) is an optimal solution to it. Therefore, constraints (5d), (7b), (7c), and (9) guarantee the optimality of the defender’ s prob- lem, and hence, can be used to con vert the ADBLP to a single lev el problem as: minimize u,v ,β c T 1 u + d T 1 v (10a) subject to A 1 u ≥ b 1 (10b) A 2 u + A 3 v ≥ b 2 (10c) A T 3 β = d 2 (10d) β T b 2 − β T A 2 u − d T 2 v ≥ 0 (10e) β ≥ 0 . (10f) The bilinear term β T A 2 u in (10e) is non-conv ex and hard to solve. T o overcome this difficulty , Benders’ decomposition is utilized to decompose this optimization problem into two problems, with u as the variable for the master problem (MP) 4 and v , β as the variables for the slave problem (SP). The MP takes the follo wing form: minimize u,α c T 1 u + α (11a) subject to A 1 u ≥ b 1 (11b) where α is a v ariable introduced to represent d T 1 v ∗ , which will then be updated by adding cuts. The SP is given by: minimize v ,β d T 1 v (12a) subject to β T b 2 − d T 2 v − β T A 2 u ≥ 0 ( δ ) (12b) A 3 v ≥ b 2 − A 2 u ( γ ) (12c) A T 3 β = d 2 ( λ ) (12d) β ≥ 0 . (12e) At the optimal solution of the SP giv en by (12), we hav e d T 1 v ∗ = γ T b 2 + λ T d 2 − γ T A 2 u. (13) An optimality cut can be added to the MP by taking the right hand side of (13): α ≥ γ T b 2 + λ T d 2 − γ T A 2 u. (14) Note that (14) is in the MP , and therefore, u is again a variable. If the SP is infeasible with a gi ven u , slack v ariables s i , i = 1 , 2 , 3 , can be introduced to all of the SP constraint to solve the relaxed SP: minimize v ,β,s i d T 1 v (15a) subject to β T b 2 − d T 2 v − β T A 2 u + s 1 ≥ 0 ( ˆ δ ) (15b) A 3 v + s 2 ≥ b 2 − A 2 u ( ˆ γ ) (15c) A T 3 β + s 3 = d 2 ( ˆ λ ) (15d) β ≥ 0 . (15e) where s i , i = 1 , 2 , 3 are the slack variables introduced to ensure feasibility of the relaxed SP . Then, instead of an optimality cut (14), a feasibility cut is added to the MP: 0 ≥ ˆ γ T b 2 + ˆ λ T d 2 − ˆ γ T A 2 u. (16) The MP and SP can then be solved iterativ ely , with the MP updating u and the SP updating cuts in each iteration. Algorithm 1 Modified Benders’ Decomposition for Bi-level Linear Programs (MBD) 1) Set the iteration number j = 1 and let u (0) = 0 . 2) Solve the SP (12) with u = u ( j − 1) . 3) If the SP is infeasible, solve the relaxed SP (15) and obtain ( ˆ γ ( j ) , ˆ λ ( j ) ) , then add a feasibility cut of form (16) to the MP . Otherwise, solve SP (12) to get ( v ( j ) , β ( j ) , γ ( j ) , λ ( j ) ) , and add an optimality cut of form (14) to the MP . 4) Solve the MP with added cuts and obtain the solution ( u ( j ) , α ( j ) ) . 5) If | d T 1 v ( j ) − α ( j ) α ( j ) | <  , stop. The optimal objective value is obtained as c T 1 u ( j ) + d T 1 v ( j ) . Otherwise, let j = j + 1 and go to step 2). Solving the SP is equi valent to solving the second le vel SCED under attack (2e) − (2o), while the dual variables of the SP provide information on the objectiv e function (2a). Since each cut is formulated linearly on the u domain, adding cuts to the MP does not affect its con ve xity . Thus, MBD is guaranteed to con verge in a finite number of iterations [21]. Howe ver , due to the non-conv exity of the original bi-level optimization problem, global optimal solution cannot be guaranteed [22]. Therefore, the optimal objective v alue obtained by MBD, ˆ P ∗ l,k t , is a lower bound on P ∗ l,k t , the global optimal objective. V I . S I M U L A T I O N R E S U LT S A N D D I S C U S S I O N S In this section, we present physical consequences through simulations of the attacks designed using the ADBLP de- scribed in Sec. IV. W e use the synthetic T exas system with 2000 buses, 3210 branches, and 432 generators [23]. The inputs to the ADBLP described in Sec. IV are obtained from OpenP A [24], a Jav a-based EMS simulation platform that we dev eloped in collaboration with our industry partners IncSys [25] and PowerData [26]. Without attack, the system is operating at steady-state, which means that SCED does not change the generation dispatch between each EMS loop. In the base case power flow solution, the total losses among the system is 2% of the net load. W e assume the SCED handles losses by uniformly increasing all loads by this percentage. R TCA simulates contingencies of all branches whose end bus voltages are both at least 100 kV , except radial branches. Prior to attack, R TCA reports no base case warnings nor violations, and 25 post-contingency warnings. W e e xhausti vely design attacks targeting each of those 25 contingency case warnings and test the attack consequences. In our simulations, the short- term branch limit is assumed to be 115% of the long-term limit, i.e., S k, max = 115% × S max ; the warning threshold τ = 90% ; MBD con ver gence tolerance  = 5 × 10 − 5 ; SCED look ahead time T h = 15 minutes; spinning reserve time T r = 10 minutes. The ADBLP is solved using Matlab with solver CPLEX on a 3.4 GHz PC with 32 GB RAM. A. Approac h for Attack Implementation and System V ulnera- bility Assessment Fig. 2 illustrates the implementation of the attack and the vulnerability assessment approach. For simplicity , we assume that the real loads remain unchanged during the attack period. The physical system behavior and the SCADA measurement collection are simulated by solving an A C po wer flow . The true measurements z 1 from the power flow solution are acquired by the attacker to estimate the states (denoted ˆ x 1 ). It then performs A C power flow-based R TCA to achiev e the security constraints and solves the attack design ADBLP to find the attack vector c . Recall that the second lev el of the ADBLP is a SCED in response to the attack, and by solving it the attacker obtains an estimate on the maximal physical power flow on the target line, which is the optimal objectiv e ˆ P ∗ l,k t . T o implement the designed attack, the attacker then constructs false measurements ¯ z 1 = h ( ˆ x 1 + c ) and injects ¯ z 1 to the system SE instead of the true measurements z 1 . Again, only the measurements in the attack subgraph S are changed. Since the generator outputs are known to the system, the false 5 measurements will cause the SE to estimate a set of false loads. R TCA and SCED are then performed by the system to determine the new optimal generation dispatch P ∗ G in response to the false loads. Once the generators re-dispatch, the attacker again acquires the true measurements z 2 , and estimates the new states ˆ x 2 . It then sends ¯ z 2 = h ( ˆ x 2 + c ) to the system SE to estimate new false loads. The system operator again runs R TCA with the new false loads and observes the cyber power flow ¯ P l,k t . Howe ver , the new dispatch applied on the physical system, will maximize the physical power flow on target line l under target contingency k t , and possibly cause overflo w . The true physical power flow , P l,k t , is obtained by running R TCA with the new dispatch and real loads. B. Results on Maximal Physical P ower Flows Fig. 3 compares physical power flow ˆ P ∗ l,k t predicted by the attacker , the true power flow P l,k t in the physical system, as well as the power flow (cyber) seen by the system operator ¯ P l,k t , as a function of the l 1 -norm constraint N 1 . These power flows are plotted as percentage v alues relativ e to the activ e power limit P l,k, max calculated using (4). The attacker’ s goal is to maximize the power flow on line ‘ln-2025-2055’ when line ‘ln-2054-5236’ is out of service. When the load shift L S = 10% , ˆ P ∗ l,k t and P l,k t increase as N 1 increases. This indicates that the attacks are effecti ve: they successfully cause post-contingency o verflo ws that cannot be seen by the system operators. When L S = 20% , similar results are observed, but ˆ P ∗ l,k t and P l,k t are not monotonically increasing as N 1 increases. This suggests that the MBD algo- rithm provides sub-optimal solutions, because as N 1 increases, the constraints are relaxed, and the optimal solution for a larger N 1 should be at least that of a smaller N 1 . Maximal power flow is higher when a larger load shift is allowed. With L S = 20% , N 1 = 0 . 2 , the power flow is higher than that when L S = 10% , N 1 = 2 , which indicates that in this case load shift is the dominant constraint. The true physical power flow P l,k t is slightly lower than the attacker predicted physical power flo w ˆ P ∗ l,k t . One possible reason for this phenomenon is that the attacker is solving a DC approximation of an A C system, and the reactiv e power flow may change after attack. This could result in a difference in P l,k, max before and after attack. Another possible reason is that the false measurements ¯ z 1 injected by the attacker cause a different set of security constraints than those that the attacker used to solve the attack design ADBLP . The attacker generates the security constraints by running R TCA using the true measurements, but those constraints generated by the system R TCA are based on the false measurements after attack. As a result, the system SCED solution may be different than the attacker predicted re-dispatch. One approach for the attacker to pre vent this situation is to run its own R TCA using the false measurements and include any newly appeared security constraints into the attack design ADBLP , until there are no more new security constraints. Howe ver , this approach has no conv ergence guarantee, and could be too time-consuming to launch the attack in real-time. Note that in order for the attacks to actually cause post- contingency violations requires a particular contingency to P ow e r F l ow S ys t e m S C E D S ys t e m S t at e E s t i m at or S ys t e m R T C A P ow e r F l ow S ys t e m S t at e E s t i m at or A t t ac k e r ’ s R T C A R T C A A t t ack D e s i gn A D B L P C r e at e F al s e M e as u r e m e n t s R T C A C r e at e F al s e M e as u r e m e n t s S ys t e m R e s p on s e t o A t t ac k C yb e r C on s e q u e n c e s P h ys i c al C on s e q u e n c e s A t t ack e r ’ s S t at e E s t i m at or A t t ac k e r ’ s S t at e E s t i m at or Fig. 2: Attack implementation and system vulnerability assessment approach. occur . Thus, the attacker has to create the target contingency itself, or gain insider knowledge about when the contingenc y is likely to occur . Both are plausible for sophisticated attackers. More aggressiv ely , the attacker can aim to create base case ov erflo ws, but the N − 1 reliable constraints may push the system to operate conservati vely . In the synthetic T exas sys- tem, there is no branch whose base case power flow is higher 6 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 90 95 100 105 110 115 Predicted Physical Cyber 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 90 95 100 105 110 115 120 Maximal Power Flow (%) Predicted Physical Cyber Fig. 3: Comparison of attacker predicted, physical, and cyber power flows on line ‘ln-2025-2055’ under contingency ‘ln-2054-5236’, (a) L S = 10% ; (b) L S = 20% . than τ prior to the attack. Thus, to cause base case ov erflow , the attacker has to shift a tremendous amount of load that may easily trigger an alarm at the control center . Moreov er , a large load shift will mov e the system operating condition dramatically with high probability , and thereby create ne w security constraints that are not considered when designing the attack. Thus, the consequences of the attack become unpredictable for the attacker . W e have attempted to design a base case attack targeting branch ‘ln-7058-7095’ that has the highest base case power flo w in percentage, but no overflo w can be found with L S = 90% and N 1 = 20 . With L S = 100% and N 1 = 20 , the attacker’ s predicted power flo w reaches 102 . 29% , but the false measurements create 3197 warnings and 24773 violations at the R TCA solution. C. Results on Attack Resources Fig. 4 illustrates the relationship between maximal power flow and l 0 -norm of the attack vector ( i.e. the number of center buses in the attack) versus the l 1 -norm constraint N 1 for target line ‘ln-2025-2055’ under contingency ‘ln-2054-5236’, with different load shift constraints. As N 1 increases, so does the l 0 - norm of the attack, indicating that l 1 -norm is a valid proxy for l 0 -norm for our problem. If a larger load shift is allowed, the maximal power flow on target line increases, but the resulting l 0 -norm may decrease for the same N 1 . This indicates a trade- off between load shift and attacker’ s resources: as the attacker attempts to av oid detection by minimizing load changes, it will require control over a larger portion of the system to launch a comparable attack. D. Comparison of Physical and Cyber RTCA results Fig. 5 compares the physical and cyber R TCA results after the re-dispatch resulting from an attack on tar get line ‘ln- 2025-2055’ under contingency ‘ln-2054-5236’ with load shift L S = 10% , N 1 = 2 . The cyber post-contingency power flows 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 50 100 150 200 250 300 350 400 450 500 Fig. 4: Comparison of the l 0 -norm of the attack vector for target line ‘ln-2025-2055’ under contingency ‘ln-2054-5236’. on the x-axis represent what the system operator observes, while the y-axis represents the post-contingency power flows in the physical system. There is no point beyond 100% of the x-axis, which indicates that the system operator sees no post- contingency violation after the attack. Therefore, the attack successfully spoofed the operator that the system is in a secure state, while in reality , the tar get line has a 112.2% post-contingency ov erflow . In addition, there are four post- contingency violations that are caused by the same attack, ev en though they are not the attacker’ s targets. 8 0 8 2 8 4 8 6 8 8 9 0 9 2 9 4 9 6 9 8 1 0 0 C y b e r P o w e r F l o w ( % ) 8 0 8 5 9 0 9 5 1 0 0 1 0 5 1 1 0 P h y s i c a l P o w e r F l o w ( % ) V i o l a t i o n V i o l a t i o n Fig. 5: Comparison of the physical and cyber R TCA results after re-dispatch. E. Statistical Results on Attack Consequences As mentioned at the beginning of Sec. VI, we exhausti vely tested attacks targeting the 25 branches with post-contingency warnings. The designed attacks successfully cause overflo ws on 8 out of the 25 target branches. T able I gives the statistical results on attack consequences of these 8 branches. W e deriv ed attacks using l 1 -norm constraints in the range from N 1 = 0 . 2 to N 1 = 2 . The table shows the resulting ranges in maximal power flow and l 0 -norm of the attack vector c across this range. The load shift constraint L S = 10% . The prefix ‘ln’ indicates a transmission line and ‘tx’ indicates a transformer . From the maximal po wer flow range, we can see that some branches are more vulnerable than others, as they have higher 7 ov erflo ws. Thus, the system operators can identify critical lines and critical contingencies for attack protection purposes. For example, they can artificially reduce the line limit to keep the attack from being successful. Measurements around vulnerable branches can be encrypted to prev ent them from being modified. In our ADBLP , the load shift constraint characterizes the detectability of the attack, indicating that load abnormally detectors can help system operators distinguish between natural load changes and possible cyber attacks based on load redistribution. T ABLE I: Statistical Results on Maximal Physical Power Flow and l 0 -norm of the Attack V ector with N 1 ∈ [0 . 2 , 2] T arget Contingency Max PF Range (%) k c k 0 Range ln-6188-7305 ln-7058-7095 101.92–105.08 133–442 ln-6240-6287 ln-6141-6239 102.43–106.76 137–314 ln-7233-7251 tx-6063-6062 105.41–107.90 156–485 ln-1003-1055 ln-3046-3078 102.80–102.94 163–520 ln-2025-2055 ln-2054-5236 107.98–111.00 90–461 ln-2070-5237 ln-2054-5236 101.35–104.35 90–461 ln-1003-1055 ln-1004-3133 102.43–102.56 160–513 ln-7059-7407 ln-7058-7406 100.38–102.24 154–488 V I I . C O N C L U S I O N W e have ev aluated the vulnerability of N − 1 reliable power systems to unobservable FDI attacks via the physical consequences of such attacks. Such N-1 reliable systems are assumed to be operated by an EMS consisting of SE, R TCA, and SCED. The attacker injects intelligently designed f alse measurements to the SE that bypass the bad data detector , and cause the SE to estimate false loads. The SCED re-dispatch resulting from the false loads leads to the power flo w on a target line (picked by the attacker) to be maximized. W e hav e also highlighted the knowledge required by the attacker to design such attacks. In the worst case, the attacker can perform exactly the same R TCA and SCED as the system does, and hence, can approximately predict the system re- sponse to the attacks. Designing these attacks inv olves solving an ADBLP that is non-conv ex and difficult to solve in general. An efficient algorithm based on Benders’ decomposition is introduced to solve the attack design ADBLPs. The designed attacks can successfully cause post- contingency overflo ws on target branches. Moreover , they may create more violations on branches other than the target one. Our vulnerability assessment approach can help system operator identify critical branches and critical contingencies to design protection schemes. Future work will include designing countermeasures to detect, identify , and mitigate such attacks. A C K N O W L E D G M E N T This material is based upon work supported by the National Science Foundation under Grant No. CNS-1449080, and grant S-72 from the Power System Engineering Research Center (PSERC). W e would like to thank the following at ASU: Mr . Andrea Pinceti for creating the base case, Mr . Roozbeh Khodadadeh for help with the test platform, and Prof. Kory Hedman and his team for their support with R TCA and SCED. W e also thank Dr . Robin Podmore (IncSys) and Mr . Christopher Mosier (Powerdata) for the OpenP A software. R E F E R E N C E S [1] A. Mittal, J. Hazra, N. Jain, V . Goyal, D. P . Seetharam, and Y . Sabharwal, “Real time contingency analysis for power grids, ” in ICPP 2011 , 2011. [2] FERC, “Security constrained economic dispatch: definition, practices, issues and recommendations, ” Federal Energy Regulatory Commission, T ech. Rep., 2006. [Online]. A vailable: https://www .ferc.gov/industries/ electric/indus- act/joint- boards/final- cong- rpt.pdf [3] K. Zetter, “Inside the cunning, unprecedented hack of Ukraine’ s power grid, ” March 2016. [Online]. A vailable: http://www .wired.com/2016/03/ inside- cunning- unprecedented- hack- ukraines- power - grid/ [4] Y . Liu, P . Ning, and M. K. Reiter, “False data injection attacks against state estimation in electric power grids, ” in Pr oceedings of the 16th A CM Confer ence on Computer and Communications Security , ser . CCS ’09, Chicago, Illinois, USA, 2009, pp. 21–32. [5] O. Kosut, L. Jia, R. J. Thomas, and L. T ong, “Malicious data attacks on the smart grid, ” IEEE T ransactions on Smart Grid , vol. 2, no. 4, pp. 645–658, 2011. [6] G. Hug and J. A. Giampapa, “V ulnerability assessment of A C state estimation with respect to false data injection cyber-attacks, ” IEEE T ransactions on Smart Grid , vol. 3, no. 3, pp. 1362–1370, 2012. [7] J. Kim and L. T ong, “On topology attack of a smart grid: Undetectable attacks and countermeasures, ” IEEE JSAC , vol. 31, no. 7, pp. 1294– 1305, 2013. [8] J. Zhang and L. Sankar, “Physical system consequences of unobservable state-and-topology cyber -physical attacks, ” IEEE T ransactions on Smart Grid , vol. 7, no. 4, pp. 2016–2025, July 2016. [9] R. Moslemi, A. Mesbahi, and J. M. V elni, “Design of robust profitable false data injection attacks in multi-settlement electricity markets, ” IET Generation, T ransmission Distribution , vol. 12, no. 6, pp. 1263–1270, 2018. [10] J. Liang, L. Sankar, and O. Kosut, “V ulnerability analysis and conse- quences of false data injection attack on power system state estimation, ” IEEE T ransactions on P ower Systems , vol. 31, no. 5, pp. 3864–3872, Sept 2016. [11] L. Jia, J. Kim, R. J. Thomas, and L. T ong, “Impact of data quality on real-time locational mar ginal price, ” IEEE T rans. P ower Systems , v ol. 29, no. 2, pp. 627–636, 2014. [12] Y . Y uan, Z. Li, and K. Ren, “Modeling load redistribution attacks in power systems, ” Smart Grid, IEEE T ransactions on , vol. 2, no. 2, pp. 382–390, June 2011. [13] ——, “Quantitati ve analysis of load redistribution attacks in power systems, ” P arallel and Distributed Systems, IEEE T ransactions on , vol. 23, no. 9, pp. 1731–1738, Sept 2012. [14] D. Alderson, G. Brown, W . Carlyle, and R. W ood, “Solving defender- attacker-defender models for infrastructure defense, ” in 12th INFORMS Computing Society Confer ence , 2011. [15] L. Mathiesen, “Computation of economic equilibria by a sequence of linear complementarity problems, ” Mathematical Pr ogramming Study , vol. 23, pp. 144–162, 1985. [16] Z.-Q. Luo, J.-S. Pang, and D. Ralph, Mathematical pr ograms with equilibrium constraints . Cambridge Univ ersity Press, 1996. [17] “The stuxnet worm: A cyber -missile aimed at iran, ” The Economist, T ech. Rep., 24 September 2010. [Online]. A vailable: http://www . economist.com/blogs/babbage/2010/09/stuxnet_worm [18] J. F . Benders, “Partitioning procedures for solving mixed-v ariables programming problems, ” Numerische Mathematik , no. 4(3), pp. 238– 252, September 1962. [19] A. J. Conejo, R. Minguez, E. Castillo, and R. Garcia-Bertrand, Decom- position T echniques in Mathematical Pr ogr amming . Springer . [20] S. Boyd and L. V andenberghe, Conve x Optimization . Cambridge Univ ersity Press, 2004. [21] A. M. Geoffrion, “Generalized Benders’ decomposition, ” Optimization Theory and Applications , vol. 10, no. 4, 1972. [22] N. V . Sahinidis and I. E. Grossmann, “Con vergence properties of gener- alized Benders’ decomposition, ” Computers and Chemical Engineering , vol. 15, p. 481, 1991. [23] “ACTIVSg2000: 2000-bus synthetic grid on footprint of Texas, ” Sep. 2017. [Online]. A v ailable: https://electricgrids.engr.tamu.edu/ electric- grid- test- cases/activsg2000/ [24] “OpenP A. ” [Online]. A vailable: https://powerdata.com/openpa/ [25] “IncSys. ” [Online]. A vailable: http://www .incsys.com/ [26] “PowerData. ” [Online]. A vailable: https://powerdata.com/ 8