From Static to Dynamic Anomaly Detection with Application to Power System Cyber Security
FFrom Static to Dynamic Anomaly Detectionwith Application to Power System Cyber Security
KAIKAI PAN, PETER PALENSKY, AND PEYMAN MOHAJERIN ESFAHANI
Abstract.
Developing advanced diagnosis tools to detect cyber attacks is the key to security of powersystems. It has been shown that multivariate data injection attacks can bypass bad data detection schemestypically built on static behavior of the systems, which misleads operators to disruptive decisions. In thisarticle, we depart from the existing static viewpoint to develop a diagnosis filter that captures the dynamicssignatures of such a multivariate intrusion. To this end, we introduce a dynamic residual generator approachformulated as robust optimization programs in order to detect a class of disruptive multivariate attacks thatpotentially remain stealthy in view of a static bad data detector. We investigate two possible desired features:(i) a non-zero transient and (ii) a non-zero steady-state behavior of the residual generator in the presence ofan attack. In case (i), the problem is reformulated as a finite, but possibly non-convex, optimization program.We further develop a linear programming relaxation that improves the scalability, and as such practicality,of the diagnosis filter design. In case (ii), it turns out that the resulting robust program admits an exactconvex reformulation, yielding a Nash equilibrium between the attacker and the residual generator. Thisassertion has an interesting implication: the proposed approach is not conservative in the sense that theadditional knowledge of the worst-case attack does not improve the diagnosis performance. To illustrate ourtheoretical results, we implement the proposed diagnosis filter to detect multivariate attacks on the systemmeasurements deployed to generate the so-called Automatic Generation Control signals in a three-area IEEE39-bus system. Introduction
The digital transformation of our power system does not only lead to better observability, flexibility andefficiency, but also introduces a phenomenon that is new to power system controls: cyber security threats.NIST [7] defines five functions for protecting Information and Communication Technology (ICT): (i) Identify,(ii) Protect, (iii) Detect, (iv) Respond, (v) Recover. It would be naive to think an ICT system can be perfectlyprotected in order to address the issues raised by (iii)-(v). This paper focuses on (iii) Detection for supervisorycontrol and data acquisition (SCADA) systems, which are in charge of transmitting measurement and controlsignals between power system substations and control centers. Such SCADA systems are notorious for beingbased on legacy ICT, and are a popular target for adversaries [13, 6] nowadays. The consequences of asuccessful attack on SCADA systems can be catastrophic to an economy and society in general [24, 17]. Inthis light, it is of utmost importance to detect these attacks and respond accordingly. Notably, if the maliciousattacks can be detected sufficiently fast, the corrupted signals can be disconnected or corrected by resilientcontrols, preventing further severe damage [34].
Literature on anomaly detection.
Traditionally, SCADA systems deploy bad data detection (BDD) tofilter out possible erroneous measurements due to sensor failures or anomalies [33]. The BDD process capturesonly a snapshot of the steady states of system trajectories, and thus only exploits possible static impact of
Date : September 24, 2019.The authors are with the Delft University of Technology, The Netherlands (email: { K.Pan, P.Palensky,P.MohajerinEsfahani } @tudelft.nl). a r X i v : . [ m a t h . O C ] S e p ROM STATIC TO DYNAMIC ANOMALY DETECTION 2 intrusions. Although this method can perform successfully in detecting basic attacks, it may fail in thepresence of the so-called stealthy multivariate attacks that carefully launch synthesized false data injectionsgiven full knowledge of the system model [15].It was first explored in [20] that such an attack can perturb the state estimation function without triggeringalarms in BDD. Since then vulnerability and impact analysis of stealthy attacks on power systems have beena prominent subject in the literature. A typical notion to quantify the vulnerability to stealthy attacks isdirectly concerned with the level of efforts required to alter specific measurements [12, 27]. Without advanceddiagnosis tools, tampering measurements remains undetected, causing state deviations, equipment damages oreven cascading failures [18]. Techniques proposed to deal with stealthy attacks include statistical methods suchas sequential detection using Cumulative Sum (CUSUM)-type algorithms [16], and measurements consistencyassessment under certain observability assumptions [35]. A detection method that leverages online informationis described in [3], which is applicable by ensuring the availability and accuracy of load forecasts and generationschedules. In [19], a mechanism is introduced to formulate the detection scheme as a matrix separationproblem, but it only recovers intrusions among corrupted measurements over a particular period of time.These techniques are essentially static detection methods that may be confined by certain prior assumptionson the distribution of measurement errors. Despite an extensive and ongoing literature focusing on the staticpart of BDD mechanism, the following question remains largely unexplored:
Would it be possible to detect stealthy multivariate attacks in a real-time operation by exploiting the attackimpact on the dynamics of system trajectories during the transient?
The importance of an appropriate answer to this question has been reinforced thanks to recent advancesin sensing technology in the modern power systems. Our main objective in this article is to address thisquestion.
Related work.
Detection methods concerning system dynamics have primarily emerged under the topic of fault detection and isolation filters . A subclass of these schemes is the observer-based approach applied initiallyto linear models [21]; see also [9] for a comprehensive summary of the large body of literature. The authorsin [25] further extend the modeling framework to general linear differential-algebraic equations (DAEs),enhancing the applicability of such methods particularly for power system applications due to the commongoverning physical laws in this setting. Recently, a variant of observer-based methods is also investigatedin [1] so as to deal with unknown natural exogenous inputs.An inherent shortcoming of many observer-based approaches is that the degree of the resulting diagnosisfilter is effectively the same as the system dynamics, which may yield an unnecessarily complex filter inlarge-scale power systems. To our best of knowledge, there are relatively much fewer studies in the literatureon the design of the reduced-order observers where the conditions for a minimum order existence need to besatisfied [9, 10]. The closest approach in the literature is [23] where a scalable optimization-based filter designis developed for high-dimensional nonlinear control systems. However, the proposed method opts for mainlydealing with a single fault scenario, and may not be as effective in case of smart multivariate adversarialinputs.An effective approach toward security and modeling the interaction between attackers and detectors buildson the rich framework of game theory. Recently, the authors in [32] propose a two player mixed strategy gameto address a dynamic resource-planning problem between an attacker targeting the communication equipmentand a defender protecting the control network. Similar frameworks have also been deployed to model thedynamics of information flow between an advanced persistent threat and a detector [30, 29].
ROM STATIC TO DYNAMIC ANOMALY DETECTION 3
Our contributions.
The main objective of this article is to develop a diagnosis filter to detect multivariatedata injection attacks in a real-time operation. For this purpose, considering a class of disruptive multivariateattack scenarios (Definition 2.5), we first characterize the attack impact on power system dynamics througha set of differential equations. Having transferred the dynamics into the discrete-time domain, we furtherrestrict the diagnosis filter to a family of dynamic residual generators that entirely decouples the contributionsof the attacks from the system states and natural disturbances. In order to identify an admissible multivariateattack scenario, we propose an optimization-based framework to robustify the diagnosis filter with respectto such attacks, i.e., aiming to design a filter whose residual (output) is sensitive to any plausible disruptivemultivariate attacks. The main contributions of this article are as follows:(i) Unlike the existing literature, we go beyond a static viewpoint of anomaly detection to capture theattack impact on the dynamics of system trajectories. To this end, we characterize the diagnosis filterdesign approach as a robust optimization program. It is guaranteed that while the filter residual isdecoupled from system states and disturbances, it still remains sensitive to all admissible disruptivemultivariate attacks even if the attacker has full knowledge about the diagnosis filter architecture(Definition 4.1 and the program (18)).(ii) To detect attacks during the transient behavior, we reformulate the resulting robust program as afinite, possibly non-convex, optimization program (Theorem 4.3). To improve the scalability of theproposed solution, we further propose a linear programming relaxation which is highly tractable forlarge scale systems (Corollary 4.4). It is guaranteed that if the optimal value of the relaxed programis positive, the resulting diagnosis filter is able to detect any admissible disruptive attack scenarios,which may remain stealthy through the lens of a static detector.(iii) We further explore the steady-state behavior of the diagnosis filter in the presence of a plausible attackscenario (Lemma 4.6). In this case, we develop an exact convex reformulation of the resulting robustprogram. As a byproduct, we show that the proposed solution is indeed a Nash equilibrium (saddlepoint) between the attacker and the residual generator (Theorem 4.7). An interesting implication ofsuch a Nash equilibrium is that the information of the attack signal may not necessarily improve theperformance of the diagnosis filter. In other words, if the proposed convex optimization fails to havea desirable feasible solution, it then implies that there exists a disruptive stealthy attack where theexact knowledge of the attack signal still does not help design a successful residual generator.In addition to the above theoretical results, we validate the performance and effectiveness of the proposeddiagnosis filter on a multi-area IEEE 39-bus system. Numerical results illustrate that the diagnosis filtersuccessfully generates a residual “alert” in the presence of multivariate attacks that are stealthy in a staticviewpoint, even in a noisy environment with imprecise measurements.Section 2 introduces the problem of power system cyber security, and the challenges posed by multivariateattacks are highlighted. Section 3 discusses a model instance of power system dynamics under attacks onmeasurements. Our diagnosis filter design is proposed in Section 4 where an optimization framework isintroduced, and numerical simulations are reported in Section 5.
Notation.
The symbols R , N represent the set of real numbers and integers, respectively. Given a matrix A ∈ R m × n , A (cid:62) denotes its transpose, and the space Im( A ) represents its range space. Throughout the paper,the matrix I is the identity matrix with an appropriate dimension. Given a column vector a ∈ R m , diag( a )denotes an m × m diagonal matrix with the elements of vector a sitting on the main diagonal and the restof the elements being zero. We also denote by diag[ A , A , . . . , A k ] a block matrix whose main diagonalelements are the matrices A , A , . . . , A k . Given a vector a ∈ R m , the associated (cid:96) ∞ − norm is denoted by (cid:107) a (cid:107) ∞ = max i ≤ m | a i | . ROM STATIC TO DYNAMIC ANOMALY DETECTION 4
Electrical Power GridElectrical Power Grid d [ k ] u [ k ] Y [ k ]Diagnosis FilterDiagnosis Filter r [ k ]Feedback ControllerFeedback Controller Communication NetworkCommunication Network f [ k ] Figure 1.
Schematic block diagram of the system model.2.
Problem Statement
Static detection and system modeling
For a power grid, measurements are collected by remote sensors and transmitted through a SCADAnetwork. The typical BDD is conducted to detect the erroneous measurements at each time instance. We cansee this as a static process: it only concerns the system states X [ k ] ∈ R n X and measurements Y [ k ] ∈ R n Y attime step k ∈ N , which can be described by Y [ k ] = CX [ k ] + D f f [ k ] , (1)where C ∈ R n Y × n X is the measurement matrix, and f [ · ] ∈ R n f represents the data injection attacks onmeasurements. Note that the matrix D f characterizes which measurement is vulnerable to attacks. Itis customary to define a residual signal for a static detector, r S [ k ] := Y [ k ] − ˆ Y [ k ], where ˆ Y [ · ] denotesthe estimated measurements. In the traditional weighted least squares estimation, the estimate of state is( C (cid:62) C ) − C (cid:62) Y [ k ], assuming that C has full column rank with high measurement redundancy. Then themeasurements estimate is C ( C (cid:62) C ) − C (cid:62) Y [ k ], and the residual signal can be further expressed as r S [ k ] = (cid:0) I − C ( C (cid:62) C ) − C (cid:62) (cid:1) Y [ k ] . (2)Such an anomaly detector has shown a good effectiveness in detecting erroneous data and basic attacks [8].However, in the face of coordinated attacks on multiple measurements, this static detector can fail. In thisarticle, motivated by this shortcoming, we take a dynamic design perspective where we shift the emphasis onan attack as a static process to its effects on power system dynamics. In particular, we opt for differentiatingthe attack impact on the systems trajectories from natural disturbances such as load deviations.To model its impact on the dynamics, let us consider a more general modeling framework in Figure 1.The electrical grid is operated by a digital controller that receives measurements as inputs and sends controlsignals to the actuators through communication networks. These transmitted data are applied in discrete-time samples. On the power grid side, the input d [ k ] ∈ R n d represents natural disturbances. On the controllerside, a control signal u [ k ] ∈ R n u is computed given the measurements Y [ k ]. Note that with the closed-loop ROM STATIC TO DYNAMIC ANOMALY DETECTION 5 control, the corruptions f [ k ] on the measurements would affect the system dynamics. The dynamics of theclosed-loop system is (cid:40) X [ k + 1] = A x X [ k ] + B d d [ k ] + B u u [ k ] ,Y [ k ] = CX [ k ] + D f f [ k ] , (3)where A x , B d and B u are constant matrices. Let us highlight the difference between the dynamical system (3)and the respective static counterpart (1). In fact, the time independence of the first equation in (3) describesthe dynamics of the system, while the algebraic equation (1) represents the relation on each time instance anddescribes a static relation between the states and outputs. The aim of this study is to exploit such dynamicsinformation in (3) in order to design a diagnosis filter to detect stealthy multivariate attacks. To illustrate theattack impact on the system dynamics, we can simply consider the feedback controller as a linear operatorsuch that u [ k ] = GY [ k ] where G ∈ R n u × n Y is a matrix gain. By defining the closed-loop system matrices A cl := A + B u GC and B f := B u GD f , we can reformulate (3) into (cid:40) X [ k + 1] = A cl X [ k ] + B d d [ k ] + B f f [ k ] ,Y [ k ] = CX [ k ] + D f f [ k ] . (4) Remark 2.1 (Dynamic feedback controller) . The restriction to only a static feedback controller u [ k ] = GY [ k ] to transfer from (3) to (4) is without loss of generality. Namely, the proposed framework is rich enough tosubsume a dynamic controller architecture as well. Indeed, when the controller has certain dynamics, itsuffices to augment the system dynamics (3) with the controller states and outputs. We refer to Appendix 2.1,for such a detailed analysis. Remark 2.2 (Attacks impact on the dynamics of system trajectories) . In light of (4) , matrices B f , D f capture the attack impact on the power system dynamics, mapping attacks f [ · ] to the system states andmeasurements respectively. In the following, we show that the state-space description (4) is a particular case of DAE model. Byintroducing a time-shift operator q : qX [ k ] → X [ k + 1], one can fit (4) into H ( q ) x [ k ] + L ( q ) y [ k ] + F ( q ) f [ k ] = 0 , (5)where x := [ X (cid:62) d (cid:62) ] (cid:62) represents the unknown signals of system states and disturbances; y := Y containsall the available data for the operator. Let n x and n y be the dimensions of x [ · ], y [ · ]. We denote n r as thenumber of rows in (5). Then H, L, F are polynomial matrices in terms of the time-shift operator q with n r rows and n x , n y , n f columns separately, by defining, H ( q ) := (cid:34) − qI + A cl B d C (cid:35) , L ( q ) := (cid:34) − I (cid:35) , F ( q ) := (cid:34) B f D f (cid:35) . Challenge: multivariate attacks
We start this subsection with an existing result characterizing the set of stealthy multivariate attacks thatcan bypass the static detector.
Lemma 2.3 (Stealthy attack values [20, Theorem 1]) . Consider the measurement equation (1) and the staticdetector with the respective residual function (2) . Then, an attack f [ · ] remains stealthy, i.e., it does notcause any additional residue to (2) , if it takes values from the set F := (cid:8) f [ k ] ∈ R n f : D f f [ k ] ∈ Im ( C ) , k ∈ N (cid:9) , (6) ROM STATIC TO DYNAMIC ANOMALY DETECTION 6
One can observe that a stealthy attack D f f [ · ] described in (6) has the knowledge of the system model (1)through the range space of C . That is, it represents a tampered value D f f [ k ] = C ∆ X where ∆ X ∈ R n X canbe any injected bias influencing certain sensor measurements. Such multivariate attacks would also challengethe detector design as they may neutralize the diagnosis filter outputs. Assumption 2.4 (Stationary attacks) . Throughout this article, we consider attacks f [ · ] that are time-invariant, i.e., f [ k ] = 0 for all k ≤ k min ; f [ k ] = f ∈ F for all k > k min . Namely, the attack occurs as aconstant bias injection f on measurements during the system operations at a specific unknown time instance k min , and it remains unchanged since then. Advanced attacks also pursue a maximized impact on the system dynamics. Thus, an adversary wouldtry to inject “ smart ” false data, possibly with large magnitudes, in such a way that it causes the maximumdamage. The next definition opts to formalize this class of attacks.
Definition 2.5 (Disruptive stealthy attack) . Consider a set of vectors F b := [ f , f , . . . , f d ] representing afinite basis for the set of stealthy attacks (6) , i.e., the set F defined in (6) can equivalently be represented by F = (cid:40) F (cid:62) b α = d (cid:88) i =1 α i f i (cid:12)(cid:12)(cid:12) α = [ α , α , · · · , α d ] (cid:62) ∈ R d (cid:41) . We call a signal f ∈ F disruptive stealthy attack if its corresponding coefficients α is a polytopic set, i.e., itbelongs to A := (cid:8) α ∈ R d | Aα ≥ b (cid:9) , (7) where A ∈ R n b × d and b ∈ R n b are given matrices. We emphasize that the subsequent analysis and the proposeddiagnosis filter design only rely on the convexity of the set A . Namely, the choice (7) may be adjusted accordingto the application at hand, as long as the convexity of the set is respected. Cyber Security of Power Systems: AGC modeling
In this section, we first go through a modeling instance of power system dynamics in the form of (4):Automatic Generation Control (AGC) closed-loop system under attacks. This model will be used to validateour diagnosis filter. Figure 2 depicts the diagram of a three-area IEEE 39-bus system. AGC is a feedbackcontroller that tunes the setpoints of participated generators (e.g., G11 of Area 1) to maintain the frequencyas its nominal value and the tie-line (e.g., L1-2 between Area 1 and 2) power as the scheduled one.In the work of AGC, a linearized model is commonly used for the load-generation dynamics [28]. For athree-area system, the frequency dynamics in Area i can be written as∆ ˙ ω i = 12 H i (∆ P m i − ∆ P tie i − ∆ P l i − D i ∆ ω i ) , (8a)where H i is the equivalent inertia constant; D i is the damping coefficient and ∆ P l i denotes load deviations.Here ∆ P tie i , ∆ P m i represent the total tie-line power exchanges from Area i and the total generated powerin Area i , i.e., ∆ P tie i = (cid:80) j ∈E i ∆ P tie i,j where E i denotes the set of areas that connect to Area i , and∆ P m i = (cid:80) G i g =1 ∆ P m i,g where G i denotes the number of participated generators in Area i , and we have∆ ˙ P m i,g = − T ch i,g (∆ P m i,g + 1 S i,g ∆ ω i − φ i,g ∆ P agc i ) , (8b)∆ ˙ P tie i,j = T ij (∆ ω i − ∆ ω j ) , (8c)where T ch i,g is the governor-turbine’s time constant; S i,g denote the droop coefficient; T ij is the synchronizingparameter between Area i and j . Note that ∆ P agc i is the signal from AGC for the participated generators to ROM STATIC TO DYNAMIC ANOMALY DETECTION 7
Area 3Area 1 Area 2
G31G11 G22G21
G23
L1-3-1 L1-3-2 L1-2 L2-3
G12 G32
Figure 2.
Three-area 39-bus system: the measurements of the tie-lines (in red) L1-3, L1-2,L2-3 are attacked.track the load changes, and φ i,g is the participating factor, i.e., (cid:80) G i g =1 φ i,g = 1. After receiving the frequencyand tie-line power measurements, the area control error (ACE) is computed for an integral action, ACE i = B i ∆ ω i + (cid:88) j ∈E i ∆ P tie i,j , (8d)∆ ˙ P agc i = − K I i ACE i , (8e)where B i is the frequency bias and K I i represents the integral gain. Based on the equations (8), the linearizedmodel of Area i can be presented as the state equation˙ X i ( t ) = A ii X ( t ) + B i,d d i ( t ) + (cid:88) j ∈E i A ij X j ( t ) , (9)where X i is the state vector; d i := ∆ P l i denotes load deviations. Recall Remark 2.1 that (9) is an augmentedmodel for the closed-loop AGC system that X i consists of not only the electrical grid states (e.g., frequency,generator output and tie-line power) but also the controller state ∆ P agc i , i.e., X i := (cid:104) { ∆ P tie i,j } j ∈E i ∆ ω i { ∆ P m i,g } G i ∆ P agc i (cid:105) (cid:62) . Besides in (9), A ii is the system matrix of Area i ; A ij is a matrix whose only non-zero element is − T ij in row1 or 2 and column 3; B i,d is the matrix for load deviations.In addition to (9), we assume a measurement model with high redundancy that the measurements of eachtie-line power (∆ P tie i,j ) and the total tie-lines’ power (∆ P tie i ), the frequency (∆ ω i ), each generator output(∆ P m i,g ) and the total generated power (∆ P m i ), and the AGC controller output (∆ P agc i ) are all available.Besides, vulnerabilities within SCADA networks may allow cyber intrusions. Thus the output equation is Y i ( t ) = C i X ( t ) + D i,f f i ( t ) , (10)where Y i is the system output and C i is the output tall-matrix with full column rank. Here f i denotesmultivariate attacks and the matrix D i,f quantifies which output is attacked. In the aforementioned section, ROM STATIC TO DYNAMIC ANOMALY DETECTION 8 due to the feedback loop, attacks on the measurements would also affect the frequency dynamics. Hence thestate equation (9) during attacks becomes˙ X i ( t ) = A ii X ( t ) + B i,d d i ( t ) + B i,f f i ( t ) + (cid:88) j ∈E i A ij X j ( t ) , where B i,f is the matrix that relates attacks to system states.Using the state equations of each area, the continuous-time model of the three-area system can be obtained,˙ X ( t ) = ˜ A cl X ( t ) + ˜ B d d ( t ) + ˜ B f f ( t ) , (11)where X is the vector consisting of groups of dynamic states in each area; d is the vector for all areas’ loaddeviations; f denotes all the attack signals in the three-area, namely, X = (cid:104) X (cid:62) X (cid:62) X (cid:62) (cid:105) (cid:62) , d = (cid:104) ∆ P l ∆ P l ∆ P l (cid:105) (cid:62) , f = (cid:104) f (cid:62) l f (cid:62) f (cid:62) (cid:105) (cid:62) . In (11), ˜ A cl is the closed-loop system matrix; ˜ B d , ˜ B f are constant matrices that relate load deviations andattacks to system states. For the three-area system, these matrices are˜ A cl = A A A A A A A A A , ˜ B d = diag (cid:104) B ,d , B ,d , B ,d (cid:105) , ˜ B f = diag (cid:104) B ,f , B ,f , B ,f (cid:105) . We can also obtain the output equation of the system, Y ( t ) = CX ( t ) + D f f ( t ) , (12)where Y is the system output vector containing all the three areas’ outputs; C is the output matrix; D f quantifies all the vulnerable signals. Similarly, these matrices are Y = (cid:104) Y (cid:62) Y (cid:62) Y (cid:62) (cid:105) (cid:62) , C = diag (cid:104) C , C , C (cid:105) , D f = diag (cid:104) D ,f , D ,f , D ,f (cid:105) . To obtain the sampled discrete-time model as (4), (11) and (12) must be discretized. We deploy a zero-orderhold (ZOH) discretization for a given sampling period T s [26], A cl = e ˜ A cl T s , B d = (cid:90) T s e ˜ A cl ( T s − t ) ˜ B d d t. (13)Note that the attack matrix ˜ B f has the same matrix transformation as ˜ B d , resulting B f . The above approxi-mation is exact for a ZOH and (13) corresponds to the analytical solution of the discretization. Therefore, theabove model can be described in the form of (4) which again can be fitted into the DAE (5). In Appendix 2.2,we provide the detailed description of the involved parameters of the three-area 39-bus system as well as theattack scenarios on the AGC measurements.4. Robust Dynamic Detection
Preliminaries for diagnosis filter construction
An ideal detection aims to implement a non-zero mapping from the attack to the diagnostic signal whiledecoupled from system states and disturbances, given the available data y [ · ] in the control center. In thepower system dynamics described via a set of DAE, we restrict the diagnosis filter to a type of dynamicresidual generator in the form of linear transfer functions, i.e., r D [ k ] := R ( q ) y [ k ] where r D is the residual The inputs signals d ( · ) and f ( · ) in (11) are assumed to be piecewise constant within the sampling periods. ROM STATIC TO DYNAMIC ANOMALY DETECTION 9 signal of the diagnosis filter and R ( q ) is a transfer operator. Note that y [ · ] is associated with the polynomialmatrix L ( q ) in (5). We propose a formulation of transform operator R ( q ) as R ( q ) := a ( q ) − N ( q ) L ( q ) , where N ( q ) is a polynomial vector with the dimension of n r and a predefined order d N . To make R ( q )physically realizable, stable dynamics a ( q ) with sufficient order need to be added as the denominator whereall the roots are strictly contained in the unit circle. Note that, unlike the observer-based methods, here d N can be much less than the dimension of system dynamics. Then N ( q ) and a ( q ) are the two variables for adiagnosis filter design. By multiplying a ( q ) − N ( q ) in the left of (5), we have r D [ k ] = a ( q ) − N ( q ) L ( q ) y [ k ] = − a ( q ) − N ( q ) H ( q ) x [ k ] (cid:124) (cid:123)(cid:122) (cid:125) (I) − a ( q ) − N ( q ) F ( q ) f [ k ] (cid:124) (cid:123)(cid:122) (cid:125) (II) , (14)where term (I) in (14) is due to x [ · ] of system states and natural disturbances. Term (II) is the desiredcontribution from the attacks f [ · ]. In view of this diagnosis filter description, we introduce a class of residualgenerator which is sensitive to disruptive stealthy attacks as defined in Definition 2.5. Definition 4.1 (Robust residual generator) . Consider a linear residual generator represented via a poly-nomial vector N ( q ) . This residual generator is robust with respect to disruptive stealthy attacks introducedin Definition 2.5 if (cid:40) ( I ) N ( q ) H ( q ) = 0 , ( II ) N ( q ) F ( q ) F b α (cid:54) = 0 , ∀ α ∈ A , (15) where the basis matrix F b and the set A are the same as the ones in Definition 2.5. In the next step, we show that the polynomial equations (15) in Definition 4.1 can be characterized as afeasibility problem of a finite robust program.
Lemma 4.2 (Linear program characterization) . Consider the polynomial matrices H ( q ) = (cid:80) i =0 H i q i , N ( q ) := (cid:80) d N i =0 N i q i and F ( q ) = F , where H i ∈ R n r × n x , N i ∈ R n r , and F ∈ R n r × n f are constant ma-trices. Then, the family of robust residual generators in (15) is characterized by (cid:40) ( I ) ¯ N ¯ H = 0 , ( II ) (cid:13)(cid:13) ¯ N V ( α ) (cid:13)(cid:13) ∞ > , ∀ α ∈ A , (16) where (cid:107) · (cid:107) ∞ denotes the infinite vector norm, and ¯ N := (cid:104) N N · · · N d N (cid:105) , ¯ H := H H · · · H H ...... . . . . . . · · · H H , V ( α ) := F F b α · · · F F b α ...... . . . · · · F F b α . Proof.
The proof follows a similar line of arguments as [23, Lemma 4.2]. The key step is to observe that N ( q ) H ( q ) = ¯ N ¯ H [ I, qI, · · · , q d N +1 I ] (cid:62) , and N ( q ) F F b α = ¯ N V ( α )[ I, qI, · · · , q d N I ] (cid:62) . The rest of the prooffollows rather straightforwardly, and we omit the details for brevity. (cid:3) Robust diagnosis filter: transient behavior
In light of (16), we can define a symmetric set for the design variable ¯ N of the dynamic residual generator, N := { ¯ N ∈ R ( d N +1) n r | ¯ N ¯ H = 0 , (cid:107) ¯ N (cid:107) ∞ ≤ η } . (17) ROM STATIC TO DYNAMIC ANOMALY DETECTION 10
The second constraint in the set is added to avoid possible unbounded solutions. To design a robust residualgenerator, we aim to find an ¯ N ∈ N that for all α ∈ A , (16) can be satisfied. To this end, a naturalreformulation of the residual synthesis is to consider an objective function as the second quantity in (16)influenced by the parameters N and the attacker action α , i.e., J ( ¯ N , α ) := (cid:107) ¯ N V ( α ) (cid:107) ∞ . A successful scenariofrom an attacker viewpoint is to minimize this objective function given a residual generator. Therefore, wetake a rather conservative viewpoint where the attacker may have complete knowledge of the system modeland even the residual generator parameters, and exploits it so as to synthesize a stealthy attack. We thenreformulate the diagnosis filter design as the robust optimization program, γ (cid:63) := max ¯ N ∈N min α ∈A (cid:110) J ( ¯ N , α ) := (cid:107) ¯ N V ( α ) (cid:107) ∞ (cid:111) . (18)The optimal value γ (cid:63) of the robust reformulation (18) is indeed an indication whether the attack stillremains stealthy in the dynamic setting, i.e., if γ (cid:63) > N (cid:63) yields a diagnosis filter inthe form of (14) which detects all the admissible attacks introduced in Definition 2.5. However, if γ (cid:63) = 0, thenit implies that for any possible detectors (static or dynamic) there exists a stationary disruptive attack thatremains stealthy. In the next step, we show that the robust program (18) can be equivalently reformulatedas a finite (non-convex) optimization problem. Theorem 4.3 (Finite reformulation of (18)) . The robust optimization (18) can be equivalently described viathe finite optimization program γ (cid:63) = max ¯ N, β, λ b (cid:62) λ s.t. d N (cid:80) i =0 ( β i − β i +1 ) N i F F b = λ (cid:62) A, (cid:62) β = 1 , β ≥ , ¯ N ∈ N , λ ≥ , (19) where β = [ β , β , · · · , β d N +1 ] (cid:62) is an R d N +2 -valued auxiliary variable.Proof. See Appendix 1.1. (cid:3)
The exact reformulation program (19) for (18) is unfortunately non-convex due to the bilinearity betweenthe variables β and N i in the first constraint. In the following corollary, we suggest a convex relaxationof the program by restricting the feasible set of the variable β to a 2 d N + 2 finite possibilities where β =[0 , · · · , , · · · , (cid:62) in which the only non-zero element of the vector is the i -th element. Corollary 4.4 (Linear program relaxation) . Given i ∈ { , . . . , d N + 2 } , consider the linear program γ (cid:63)i := max ¯ N, λ b (cid:62) λ s.t. ( − i N (cid:98) i/ (cid:99) F F b = λ (cid:62) A, ¯ N ∈ N , λ ≥ , (LP i ) where (cid:98) · (cid:99) is the ceiling function that maps the argument to the least integer. Then, the solution to theprogram (LP i ) is a feasible solution to the exact robust design reformulation (19) , and max { i ≤ d N +2 } γ (cid:63)i ≤ γ (cid:63) .In particular, if for any i ∈ { , . . . , d N +2 } we have γ (cid:63)i > , then the solution to LP i offers a robust residualgenerator detecting all admissible disruptive attacks introduced by Definition 2.5. Corollary 4.4 suggests that the maximum optimal value of { γ (cid:63) , γ (cid:63) , · · · , γ (cid:63) d N +2 } and its corresponding ¯ N (cid:63) provide a suboptimal solution to the original robust design (18). ROM STATIC TO DYNAMIC ANOMALY DETECTION 11
We note that the focus of this article is on stationary (time-invariant) attacks. It is also important tohighlight that the robust design perspective (18) allows the attacker to know the system model and filterparameters. In such a setting, the detection procedure could be much more difficult if the attacker wouldbe able to dynamically adapt the attacks over the time, i.e., the attack signal is time-varying. In fact, in amultivariate attack scenario, one can construct a disruptive time-varying attack bypassing any linear residualgenerators. The next remark alludes more to this situation.
Remark 4.5 (Time-varying stealthy attacks) . Consider a multivariate attack f = [ f f · · · f n f ] (cid:62) whereeach element is a time-varying signal f i = f i [ k ] . Then, the residual (14) can be rewritten as a ( q ) r D [ k ] = − n f (cid:88) i =1 (cid:16) N ( q ) F i f i [ · ] (cid:17) [ k ] , (20) where F = [ F F · · · F n f ] represents the attack dynamics matrix. One can inspect that when the time-varyingrelation (cid:80) n f i =1 (cid:0) N ( q ) F i f i [ · ] (cid:1) [ k ] = 0 holds for every k , for instance when f n f [ k ] = − (cid:0) N ( q ) F n f (cid:1) − n f − (cid:88) i =1 (cid:16) N ( q ) F i f i [ · ] (cid:17) [ k ] , then the residual outcome (20) stays zero for all k , and as such, the attack remains undetected. The proposed robust design in (18) does not necessarily enforce a non-zero steady-state residual of thediagnosis filter under multivariate attacks. Namely, the design perspective of (18) focuses on detection ofattacks during the transient behavior without any requirements on long-term behavior of the residual. Indeed,the residual signal r D may return to zero value after a successful reaction to the attack occurrence. A morestringent perspective is to require a non-zero steady-state behavior under any admissible attack scenario in α ∈ A . This extension is addressed in the next subsection. Robust diagnosis filter: steady-state behavior
In order to design a diagnosis filter with non-zero steady-state residual “alert” when a multivariate at-tack occurs, the robust optimization (18) can be modified by a more conservative (smaller) objective func-tion J ( ¯ N , α ) := | ¯ N ¯ F α | where ¯ F := (cid:104) F F b F F b · · · F F b (cid:105) (cid:62) . (21)A similar treatment as the preceding subsection can establish a framework for computational purposes. Thenext lemma follows similar objective as in Lemma 4.2 with a more demanding requirement of the non-zerolong-term residual behavior. Lemma 4.6 (Non-zero steady-state residual characterization) . For the polynomial matrices H ( q ) , N ( q ) and F ( q ) as defined in Lemma 4.2, the family of dynamic residual generators with non-zero steady-state residualunder multivariate attacks can be characterized by the algebraic relations (cid:40) ( I ) ¯ N ¯ H = 0 , ( II ) | ¯ N ¯ F α | > , ∀ α ∈ A , (22) where ¯ F is defined in (21) , and the matrices ¯ N , ¯ H are as defined in Lemma 4.2.Proof. Recall that N ( q ) H ( q ) = ¯ N ¯ H [ I, qI, · · · , q d N +1 I ] (cid:62) . Thus if ¯ N ¯ H = 0, the diagnosis filter be-comes r D [ k ] = − a ( q ) − N ( q ) f [ k ]. Note the steady-state value of the filter residual under attacks wouldbe − a ( q ) − N ( q ) F ( q ) f | q =1 . Thus for the multivariate attack with α , the steady-state value of the filterresidual is − a (1) − N (1) F (1) F b α . The proof concludes by noting that N (1) F (1) F b α = ¯ N ¯ F α . (cid:3) ROM STATIC TO DYNAMIC ANOMALY DETECTION 12
In a similar fashion, the robust design perspective in (18) can be modified accordingly as µ (cid:63) := max ¯ N ∈N min α ∈A (cid:110) J ( ¯ N , α ) := | ¯ N ¯ F α | (cid:111) . (23)Notice the relation between the new objective function with the absolute value and the one in (18) withthe infinity-norm. As it appears in the next result, the new setting is in fact a restricted case of the finitereformulation in Theorem 4.3. Theorem 4.7 (Residual long-term behavior: exact convex reformulation and Nash equilibrium) . Considerthe minimax counterpart of the program (18) as defined ϕ (cid:63) := min α ∈A max ¯ N ∈N (cid:110) J ( ¯ N , α ) := | ¯ N ¯ F α | (cid:111) . (24) Each of the program (23) and (24) can be equivalently reformulated through the linear programs µ (cid:63) = max ¯ N, λ b (cid:62) λ s.t. ¯ N ¯ F = λ (cid:62) A ¯ N ∈ N , λ ≥ , (25a) ϕ (cid:63) = min v , v , w, α (cid:62) v + (cid:62) v s.t. ¯ Hw + v − v = ¯ F αv ≥ , v ≥ ,Aα ≥ b . (25b) Moreover, the value of each of these two programs coincide, i.e., µ (cid:63) = ϕ (cid:63) .Proof. See Appendix 1.2. (cid:3)
It is worth noting the difference between the robust perspective of (23) versus the minimax program (24).While in the design perspective of (23) the filter is oblivious to the possible attack scenarios, in the perspectiveof (24) the filter is aware of the attack signal and opts to detect that particular signal in the presence of naturaldisturbances. Obviously, the former setting is the one closer to the reality and, in general, the knowledge ofthe attack signal should help the detection significantly. This observation can indeed be translated throughthe usual weak inequality of µ (cid:63) ≤ ϕ (cid:63) . However, Theorem 4.7 indicates that the filter performance, in viewof the long-term behavior of the worst-case attack scenario, indeed does not depend on the exact knowledgeof the attacker signal and the inequality holds as the equality. We summarize this discussion in the followingremark. Remark 4.8 (Nash equilibrium interpretation) . If the linear programs (25a) (25b) admit a positive optimalvalue ϕ (cid:63) = µ (cid:63) > , then the resulting filter can detect all the admissible multivariate attacks described byDefinition 2.5 along with a non-zero steady-state residual level. On the other hand, if the optimal valuescoincide with ϕ (cid:63) = µ (cid:63) = 0 , it then implies that there is no linear filter being able to decouple the admissibleattack with α (cid:63) , the solution to (25b) , from the natural disturbances in a long-term horizon. Numerical Results
Test system and diagnosis filter description
In order to validate the effectiveness of the diagnosis filter with application to power system cyber security,we employed the IEEE 39-bus system which is well-known as a standard system for testing of new powersystem analysis. As shown in Figure 2, this system consists of 3 areas and 10 generators where 7 of them
ROM STATIC TO DYNAMIC ANOMALY DETECTION 13 time(s)
10 20 30 40 50 60 ( p . u . ) -1-0.500.5 Disturbance ∆ P l FDI Attack f tie , FDI Attack f tie FDI Attack f tie , , f tie (a) Load disturbance and basic attack time(s)
10 20 30 40 50 60 ( p . u . ) -1-0.500.5 Disturbance ∆ P l FDI Attack f tie , FDI Attack f tie , FDI Attack f tie FDI Attack f tie , , f tie (b) Load disturbance and stealthy attack time(s) r S , t i e (c) Residual of static detector under basic attack time(s) r S , t i e (d) Residual of static detector under stealthy attack time(s) r D (e) Residual of dynamic detector under basic attack time(s) r D (f) Residual of dynamic detector under stealthy attack
Figure 3.
Static detector in (2) versus dynamic detector (diagnosis filter) from Corollary 4.4under basic and stealthy attacks.are equipped with AGC for frequency control. All the participating generators in each area are with equalparticipation factors. The total load of the three-area system is 5 .
483 GW for the base of 100 MVA and 60 Hz.The generator specifications and AGC parameters of each area are referred to [4], and the linear frequencydynamics model has been developed in the preceding Section 3. Thus we result in a 19-order model in theform of (4).We apply the diagnosis filter proposed in Section 4 to detect multivariate disruptive attacks on the mea-surements of AGC system. In the following simulations, we set the degree of the dynamic residual generator d N = 3 which is much less than the order of the dynamics model, the sampling time T s = 0 . a ( q ) = ( q − p ) d N / (1 − p ) d N where p is a user-defined variable acting as the pole of the transfer operator R ( q ), and it is normalized insteady-state value for all feasible poles. The pole is set to be p = 0 . Simulation results
To evaluate the performance of the diagnosis filter, the disturbances d i = ∆ P l i are modeled as stochasticload patterns. To capture its uncertainty, as shown in Figure 3a and Figure 3b, we mainly model ∆ P l inArea 1 as random zero-mean Gaussian signals. It should be noted that tie-line power flow measurementsare much more vulnerable to cyber attacks, comparing with frequency measurements (e.g., the anomalies infrequency can be easily detected by comparing the corrupted reading with the normal one.) [5]. Therefore asindicated in Figure 2 we mainly focus on the scenario that there are 5 vulnerable tie-line power measurements,namely ∆ P tie , , ∆ P tie , , ∆ P tie , ∆ P tie , and ∆ P tie . Recalling Definition 2.5 for stealthy attack basis,thus there exist 3 basis vectors in the spanning set and we model them as follows: f = [0 . . T , f = [0 . .
15 0 .
25 0 0] T , f = [0 0 0 0 . . T (all in p . u . ). Here each basis vector lies in the range space ROM STATIC TO DYNAMIC ANOMALY DETECTION 14 of the output matrix that the corrupted measurements still align with an actual physical state, bypassingthe static detector r S [ · ]. Furthermore, without loss of generality we set A = (cid:62) and b = 1 . A and η = 10 in the set N . The design variable ¯ N of the robust residual generator is first derived by solving(18) through (LP i ). The optimal value achieves maximum for i = 2 that γ (cid:63) = 300, which implies a robustdetection during the transient behavior as Corollary 4.4. For the given ¯ N , the multivariate attack coordinates α = [2 . − . (cid:62) are obtained by solving the inner minimization of (18). Next, we look into the steady-statebehavior of the filter with the above sets N and A . For this, following Theorem 4.7 we solve (23) and (24)through the programs (25a) and (25b). It turns out that the derived optimal values satisfy the equality ϕ (cid:63) = µ (cid:63) = 0, indicating that the optimal multivariate attack with α (cid:63) , the optimizer of the program (25b)and an optimal solution to (24), is a stealthy attack in the long-term horizon. We highlight that, thanks tothe fact that the optimal values of the programs (25a) (25b) form a Nash equilibrium, even with the exactinformation of the stealthy attack coefficients α (cid:63) , we still cannot decouple the long-term behavior of theresidual from the natural disturbances; see Remark 4.8.In the first simulation, we begin with a general scenario where the multivariate attack is not carefullycoordinated, i.e., basic attack. Thus as shown in Figure 3a, only 4 of 5 vulnerable measurements are com-promised that f tie , = 0 . p.u. , f tie = 0 . p.u. , f tie , = − . p.u. and f tie = − . p.u. . Note that sincethe injected data on ∆ P tie , and ∆ P tie are inconsistent, the static detector is also expected to be triggered.To test the detectors in a more realistic setup, we also consider the presence of process and measurementsnoises. The process noise term added to the state equation of Area 1 is zero-mean Gaussian noises with thecovariance matrix R X = 0 . × diag([1 1 0 .
03 1 1 1 1] (cid:62) ), i.e., the covariance of the noise to the frequency is0.009 and the covariance of other states’ noise is 0.03 [1]. Similarly, the measurement noise term added tothe measurements of Area 1 is with the covariance matrix R Y = 0 . × diag([1 1 1 0 .
03 1 1 1 1 1] (cid:62) ), i.e., thecovariance of the frequency measurement is 0.009 and the covariance of other measurements’ noise is 0.03 [1].Note the residue r S of BDD in (2) becomes r S [ k ] = ( I − C ( C (cid:62) R − Y C ) − C (cid:62) R − Y ) Y [ k ] under the noisy system.The attacks are launched at k min = 30 sec. In Figure 3c and Figure 3e, results of the static detector in (2) andthe proposed dynamic detector (diagnosis filter) are presented. Both detectors have succeeded to generate adiagnostic signal when attacks occurred, and the diagnosis filter residual r D is significantly decoupled fromstochastic load disturbances, and keeps sensitive to the multivariate attacks for a successful detection undernoisy system settings.In the second simulation, to challenge the detectors, now the multivariate attacks have been launchedon all the 5 vulnerable measurements and the derived attack coefficient α from the optimization results hasbeen used for a more intelligent adversary. Thus in Figure 3b, the corruptions become f tie , = 0 . p.u. , f tie , = 0 . p.u. , f tie = 0 . p.u. , f tie , = − . p.u. and f tie = − . p.u. . This corresponds to the worstcase for the diagnosis filter that the adversary is given the knowledge of the residual generator’s parameter ¯ N that it tries to minimize the payoff function over A . Besides, the noisy system settings have been considered.Figure 3d and Figure 3f demonstrate all the simulation results. In Figure 3d, the static detector becomestotally blind to the occurrence of such an intelligent attack. However, as we can see in Figure 3f, even inthe worst case, the diagnosis filter works perfectly well under the noisy system, generate a residual “alert”for the presence of multivariate attacks. We can also see that the residual output becomes close to zerovalue again after a successful detection during the transient behavior in Figure 3f, which is consistent tothe aforementioned result ϕ (cid:63) = µ (cid:63) = 0 and Remark 4.8. These simulations also prove the effectiveness androbustness of the proposed diagnosis filter design. ROM STATIC TO DYNAMIC ANOMALY DETECTION 15
Further discussions
In this section we elaborate several practical aspects of the proposed filter in the preceding section.5.3.A.
Diagnosis sensitivity to filter poles
While the denominator of the filter a ( q ) in (14) is chosen rather arbitrarily, up to a stability condition, thepoles however has a significant impact on the residual sensitivity. As a general rule, the smaller the poles, thefaster the residual responds, and the more sensitive the residual responds to model imprecision and noises.Simulation results in Figure 4 in Appendix 2.3 numerically illustrate this relation when the filter poles vary.5.3.B. Other types of attacks
In addition to a smart multivariate measurement attacks, the main focus of this study, there are severalother types of attacks that we briefly discuss in the following: • Denial-of-service (DoS) attack : A type of availability attack where the attacker aims to prevent somespecific data from being delivered to the respective destinations. • Replay attack : A two-stage attack where the adversary gathers a sequence of data packets at stage 1,and then replays the recorded data afterwards at stage 2.From a detection point of view, DoS attacks are trivially detectable without any sophisticated mechanismsas the absence of data is not stealthy. In the typical DoS attack modeling, the missing data is typicallyreplaced with the last received ones [31]. In such a mechanism, the DoS can be treated as an “injection”attack. We investigate the performance of our filter in the presence of this class of attacks in Figure 5 inAppendix 2.3. Numerical results confirm that the proposed filter can successfully detect the DoS attacks. Inregard with the replay attack, the articles [22, 14] offer sufficient conditions under which plausible attacksmay remain stealthy irrespective of the detection mechanism providing that the attacker has access all thenecessary data channels and excite attack of stage 2 at a suitable time.5.3.C.
Observer-based diagnosis filters
Another major technique for anomaly detection builds on observer-based techniques. In this view, theestimate of the system states, or in more general setting output observer , is a reference to alert the abnormal-ity [11]. We close this section by a brief summary of the differences between these approaches and the oneproposed in this study. • The observer-based approaches typically yield diagnosis filters with higher dynamical system degreesthan the approach proposed in this study. A low-order diagnosis filter is often more desired due topractical aspects of online implementation particularly for large-scale power systems. • Observer-based diagnosis filters usually rely on a precondition of system observability. An extendedversion of such filters relaxes this condition to the so-called Luenberger-type conditions [2]. Ourdiagnosis filter, however, requires a weaker condition reflected through the feasibility condition of theresulting optimization programs, e.g., when the program (16) in Lemma 4.2 is feasible. • Thanks to the optimization-based framework, unlike the observer-based approaches, we have a sys-tematic approach to incorporate a multivariate attack scenario into the framework.
ROM STATIC TO DYNAMIC ANOMALY DETECTION 16 Conclusion
In this article, we investigated the problem of anomaly detection in the power system cyber security with aparticular focus on exploiting the dynamics information where tempering multiple measurements data may bepossible. Our study showed that a dynamical perspective to the detection task indeed offers powerful diagnosistools to encounter attack scenarios that may remain stealthy from a static point of view. The effectivenessof this result was validated by simulations in the IEEE 39-bus system. Future research directions that weenvision include an extension to nonlinear systems, as well as a setting exposed to the “dynamic” (time-variant) attacks in Remark 4.5, as opposed to the linear models and stationary attack scenarios studied inthis article.
Appendix I: Technical Proofs
Proof of Theorem 4.3
Let us recall that ¯
N V ( α ) = (cid:104) N F F b α N F F b α · · · N d N F F b α (cid:105) , and as such, the payoff function of therobust reformulation (18) is J ( ¯ N , α ) = max i | N i F F b α | where i ∈ { , · · · , d N } . By introducing an auxiliaryvariable β in the simplex set B := { β ∈ R d N +2 | β ≥ , (cid:62) β = 1 } , one can rewrite J as J ( ¯ N , α ) = max β ∈B d N (cid:88) i =0 ( β i − β i +1 ) N i F F b . In this light, the original robust strategy (18) can be equivalently described viamax ¯ N ∈N min α ∈A max β ∈B (cid:40) d N (cid:88) i =0 ( β i − β i +1 ) N i F F b α (cid:41) . Note that given a fixed ¯ N the inner minimax optimization is indeed a bilinear objective in the decisionvariables and the respective feasible sets A and B are convex. Since one of the sets, B , is also compact, thenthe zero-duality gap holds. Therefore, interchanging the optimization over α ∈ A and β ∈ B yields γ (cid:63) = max ¯ N ∈N , β ∈B (cid:40) min α ∈A d N (cid:88) i =0 ( β i − β i +1 ) N i F F b α (cid:41) . (26)The inner minimization of (26) is a (feasible) linear program. We can use the duality again. To this end, letus assume that the decision variables ¯ N and β are fixed and consider the Lagrangian function L ( α ; λ ) = b (cid:62) λ + (cid:16) d N (cid:88) i =0 ( β i − β i +1 ) N i F F b − λ (cid:62) A (cid:17) α, where optimizing over an unconstrained variable α becomesmin α L ( α ; λ ) = b (cid:62) λ if d N (cid:80) i =0 ( β i − β i +1 ) N i F F b = λ (cid:62) Aλ ≥ −∞ otherwise,Using the above characterization as the most inner optimization program in (26) leads tomax λ b (cid:62) λ s.t. d N (cid:80) i =0 ( β i − β i +1 ) N i F F b = λ (cid:62) A,λ ≥ . (27) ROM STATIC TO DYNAMIC ANOMALY DETECTION 17
It then suffices to combine maximizing over the auxiliary variable λ together with the variables ¯ N and β to arrive at the main result in (19). Proof of Theorem 4.7
We first prove the convex reformulation. For a given ¯ N ∈ N , the inner minimization of (23) can betranslated as min α ∈A , r r s.t. ¯ N ¯ F α − r ≤ , − ¯ N ¯ F α − r ≤ . The Lagrangian of the inner minimization reads as L ( α, r ; β, λ ) = b (cid:62) λ + (cid:0) ( β − β ) ¯ N ¯ F − λ (cid:62) A (cid:1) α + (1 − β − β ) r. Optimizing over the variables α, r yieldsmin α, r L ( α, r ; β, λ ) = b (cid:62) λ if ( β − β ) ¯ N ¯ F = λ (cid:62) Aβ + β ≤ β ≥ , β ≥ , λ ≥ −∞ otherwise.Then, combining maximization over the auxiliary variables λ , β , β together with the variable ¯ N arrives atthe optimization program, µ (cid:63) = max ¯ N, β , β , λ b (cid:62) λ s.t. ( β − β ) ¯ N ¯ F = λ (cid:62) A,β + β ≤ , β ≥ , β ≥ , ¯ N ∈ N , λ ∈ R n b , λ ≥ . (28)Note that the actual program (25a) is a restriction of (28) where the variables β and β are restricted to β = 1 and β = 0. Next, we show that this restriction is indeed without loss of generality. To this end,suppose the tuple ( β (cid:63) , β (cid:63) , ¯ N (cid:63) , λ (cid:63) ) is an optimal solution to the program (28). Note that the optimal variables β (cid:63) and β (cid:63) may satisfy one of the following three properties:(i) β (cid:63) = β (cid:63) : In this case, λ (cid:63) = 0, and therefore the optimal value µ (cid:63) = 0. This optimal solution can betrivially achieved in the program (25a) by setting ¯ N = 0.(ii) β (cid:63) > β (cid:63) : Observe that the tuple (cid:0) β (cid:48) = 1, β (cid:48) = 0, ¯ N (cid:48) = ¯ N (cid:63) , λ (cid:48) = λ (cid:63) / ( β (cid:63) − β (cid:63) ) (cid:1) is a feasible solutionwith the objective value b (cid:62) λ (cid:63) / ( β (cid:63) − β (cid:63) ). Since b (cid:62) λ (cid:63) ≥ β (cid:63) − β (cid:63) ∈ (0 , β (cid:63) − β (cid:63) = 1. Thatis, β (cid:63) = 1 and β (cid:63) = 0.(iii) β (cid:63) < β (cid:63) : Following similar steps as the previous case together with the symmetric property of thefeasible set N , one can show that the optimal value of the program (28) also coincides with therestricted version in (25a).This concludes the proof of the convex reformulation from (23) to (25a). In regard with the minimaxproblem (24), let us recall the symmetric property of the feasible set N in the variable ¯ N . With a fixed α ,the inner maximization can be directly formed as max ¯ N ∈N ¯ N ¯ F α whose Lagrangian becomes L ( ¯ N ; v, w ) = − ( (cid:62) v + (cid:62) v ) + (cid:0) w (cid:62) ¯ H (cid:62) + v (cid:62) − v (cid:62) − ( ¯ F α ) (cid:62) (cid:1) ¯ N (cid:62) , ROM STATIC TO DYNAMIC ANOMALY DETECTION 18
Optimizing over the variable ¯ N leads tomin ¯ N L ( ¯ N ; v, w ) = − (cid:62) v − (cid:62) v if (cid:40) ¯ Hw + v − v = ¯ F αv ≥ , v ≥ −∞ otherwise.Thus, combining minimization over the auxiliary variables v , v , w together with the variable α , the minimaxoptimization (24) can be reformulated as the linear program (25b).Finally, we show that the solution to programs (25) indeed forms a Nash equilibrium between the pro-grams (23) and (24). Thus far, we have reformulated maximin and minimax problems as linear programs (25).The idea is to show that these programs have the same optimal values. In fact, we show that the programsare dual of each other, and that the strong duality holds when both programs are feasible. To this end, weresort to the duality of (25a) with the Lagrangian L ( ¯ N , λ ; α, v, w ) = (cid:0) w (cid:62) ¯ H (cid:62) + v (cid:62) − v (cid:62) − ( ¯ F α ) (cid:62) (cid:1) ¯ N (cid:62) + ( α (cid:62) A (cid:62) − b (cid:62) ) λ − ( (cid:62) v + (cid:62) v ) . Optimizing over the variables ¯ N , λ yieldsmin ¯ N,λ L ( ¯ N , λ ; α, v, w ) = − (cid:62) v − (cid:62) v if ¯ Hw + v − v = ¯ F αAα ≥ bv ≥ , v ≥ −∞ otherwise.It is not difficult to see that the above program coincides with the program (25b); this concludes the proof. Appendix II: System Parameters & Added Simulation Results
Dynamic Feedback Controller Modeling
Consider a dynamical system (e.g., the electrical power system studied in Section 3). Suppose the controlsignal is implemented as a dynamic feedback controller described by the discrete-time dynamics (cid:40) X c [ k + 1] = A c X c [ k ] + B c Y [ k ] ,u [ k ] = C c X c [ k ] + D c Y [ k ] , where the input is the dynamical system measurements Y [ · ], the output the control signal u [ · ], and theinternal state of the controller is denoted by X c ∈ R n c . When an attack occurs on the measurements, itaffects the dynamics of the controller and consequently the involved physical system. To study the controldynamics together with the original dynamical system, one can augment the states of the system (3) togetherwith the controller’s as ˆ X := [ X (cid:62) X (cid:62) c ] (cid:62) . Assuming that the control signal can also be measured, one canalso introduce an augmented measurement signals as ˆ Y = [ Y (cid:62) u (cid:62) ] (cid:62) . Following this procedure, the dynamicsof the closed-loop system is described by (cid:40) ˆ X [ k + 1] = ˆ A cl ˆ X [ k ] + ˆ B d d [ k ] + ˆ B f f [ k ] , ˆ Y [ k ] = ˆ C ˆ X [ k ] + ˆ D f f [ k ] . (29)where the involved matrices are defined asˆ A cl := (cid:34) A x + B u D c C B u C c B c C A c (cid:35) , ˆ B d := (cid:34) B d (cid:35) , ˆ B f := (cid:34) B u D c D f B c D f (cid:35) , ˆ C := (cid:34) C D c C C c (cid:35) , ˆ D f := (cid:34) D f D c D f (cid:35) . ROM STATIC TO DYNAMIC ANOMALY DETECTION 19 time(s)
10 20 30 40 50 60 ( p . u . ) -1-0.500.5 Disturbance ∆ P l FDI Attack f tie , FDI Attack f tie FDI Attack f tie , , f tie (a) Load disturbance and basic attack time(s)
10 20 30 40 50 60 ( p . u . ) -1-0.500.5 Disturbance ∆ P l FDI Attack f tie , FDI Attack f tie , FDI Attack f tie FDI Attack f tie , , f tie (b) Load disturbance and stealthy attack time(s) r D (c) Residual signal r D with pole p = 0 . time(s) r D -10123 (d) Residual signal r D with pole p = 0 . time(s) r D (e) Residual signal r D with pole p = 0 . time(s) r D (f) Residual signal r D with pole p = 0 . time(s) r D (g) Residual signal r D with pole p = 0 . time(s) r D (h) Residual signal r D with pole p = 0 . Figure 4.
Results of dynamic detector (diagnosis filter) with different poles ( p =0 . , . , .
98) under basic and stealthy attacks.In this view, the augmented system (29) shares the same structure as (4) studied in the main part of thearticle for the case of static feedback controller.
AGC Parameters of the three-area 39-bus system
In this subsection we provide the involved matrices and parameters of the three-area 39 system. We takethe model description of Area 1 in the three-area system in Figure 2 of Section 3 as an instance, B ,d = (cid:104) − H (cid:105) (cid:62) ,A = T T − H − H − D H H H
00 0 − T ch , S , − T ch , φ , T ch , − T ch , S , − T ch , φ , T ch , − K I − K I − K I B . ROM STATIC TO DYNAMIC ANOMALY DETECTION 20 time(s) ∆ P t i e , ( p . u . ) -0.2-0.100.10.2 (a) ∆ P tie , under DoS attacks from k dos = 30 time(s) r D (b) Residual signal r D under DoS attacks Figure 5.
Results of dynamic detector (diagnosis filter) under DoS attacks on ∆ P tie , ( p = 0 . C i for Area 1 becomes C = (cid:62) . In Area 1, the vulnerable measurements to cyber attacks are the ones of tie-line power flows ∆ P tie , ,∆ P tie , and ∆ P tie . Thus the AGC signal ∆ P agc would be corrupted into∆ ˙ P agc = − k ( B ∆ ω + ∆ P tie , + f tie , + ∆ P tie , + f tie , ) . Then the parameters regarding multivariate attacks are f = (cid:104) f tie , f tie , f tie (cid:105) (cid:62) ,D ,f = (cid:62) , B ,f = − k − k (cid:62) . Additional simulation results
In Figure 4 we present the simulation results of the residual signal r D from the proposed diagnosis filterunder different poles ( p = 0 . , . , .
98, respectively). We also show the simulation results of the residualsignal r D from the proposed diagnosis filter under DoS attacks in Figure 5. References [1]
A. Ameli, A. Hooshyar, E. El-Saadany, and A. Youssef , Attack detection and identification for automatic generationcontrol systems , IEEE Transactions on Power Systems, (2018), p. 1.[2]
V. Andrieu and L. Praly , On the existence of a kazantzis–kravaris/luenberger observer , SIAM Journal on Control andOptimization, 45 (2006), pp. 432–456.[3]
A. Ashok, M. Govindarasu, and V. Ajjarapu , Online detection of stealthy false data injection attacks in power systemstate estimation , IEEE Transactions on Smart Grid, 9 (2018), pp. 1636–1646.[4]
H. Bevrani , Robust Power System Frequency Control , Power Electronics and Power Systems, Springer, 2008.[5]
C. Chen, K. Zhang, K. Yuan, L. Zhu, and M. Qian , Novel detection scheme design considering cyber attacks on loadfrequency control , IEEE Transactions on Industrial Informatics, 14 (2018), pp. 1932–1941.[6]
T. M. Chen and S. Abu-Nimeh , Lessons from stuxnet , Computer, 44 (2011), pp. 91–93.
ROM STATIC TO DYNAMIC ANOMALY DETECTION 21 [7]
C. Cybersecurity , Framework for improving critical infrastructure cybersecurity version 1.1 , tech. rep., National Instituteof Standards and Technology, Apr. 2018.[8]
R. Deng and H. Liang , False data injection attacks with limited susceptance information and new countermeasures insmart grid , IEEE Transactions on Industrial Informatics, (2018), p. 1.[9]
S. X. Ding , Model-based fault diagnosis techniques: design schemes, algorithms, and tools , Springer Science & BusinessMedia, 2008.[10]
X. Gao, X. Liu, and J. Han , Reduced order unknown input observer based distributed fault detection for multi-agentsystems , Journal of the Franklin Institute, 354 (2017), pp. 1464–1483.[11]
W. Ge and C.-Z. FANG , Detection of faulty components via robust observation , International Journal of Control, 47 (1988),pp. 581–599.[12]
A. Giani, E. Bitar, M. Garcia, M. McQueen, P. Khargonekar, and K. Poolla , Smart grid data integrity attacks ,IEEE Transactions on Smart Grid, 4 (2013), pp. 1244–1253.[13]
S. Gorman , Electricity grid in US penetrated by spies , The Wall Street Journal, 8 (2009).[14]
A. Hoehn and P. Zhang , Detection of replay attacks in cyber-physical systems , in American Control Conference, 2016,pp. 290–295.[15]
G. Hug and J. A. Giampapa , Vulnerability assessment of AC state estimation with respect to false data injection cyber-attacks , IEEE Transactions on Smart Grid, 3 (2012), pp. 1362–1370.[16]
S. Li, Y. Yılmaz, and X. Wang , Quickest detection of false data injection attack in wide-area smart grids , IEEE Trans-actions on Smart Grid, 6 (2015), pp. 2725–2735.[17]
G. Liang, S. R. Weller, J. Zhao, F. Luo, and Z. Y. Dong , The 2015 Ukraine blackout: Implications for false datainjection attacks , IEEE Transactions on Power Systems, 32 (2017), pp. 3317–3318.[18]
J. Liang, L. Sankar, and O. Kosut , Vulnerability analysis and consequences of false data injection attack on powersystem state estimation , IEEE Transactions on Power Systems, 31 (2016), pp. 3864–3872.[19]
L. Liu, M. Esmalifalak, Q. Ding, V. A. Emesih, and Z. Han , Detecting false data injection attacks on power grid bysparse optimization , IEEE Transactions on Smart Grid, 5 (2014), pp. 612–621.[20]
Y. Liu, P. Ning, and M. K. Reiter , False data injection attacks against state estimation in electric power grids , in 16thACM Conference on Computer and Communication Security, New York, 2009, pp. 21–32.[21]
M. A. Massoumnia, G. C. Verghese, and A. S. Willsky , Failure detection and identification , IEEE Transactions onAutomatic Control, 34 (1989), pp. 316–321.[22]
Y. Mo and B. Sinopoli , Secure control against replay attacks , in 47th Annual Allerton Conference on Communication,Control, and Computing, 2009, pp. 911–918.[23]
P. Mohajerin Esfahani and J. Lygeros , A tractable fault detection and isolation approach for nonlinear systems withprobabilistic performance , IEEE Transactions on Automatic Control, 61 (2016), pp. 633–647.[24]
P. Mohajerin Esfahani, M. Vrakopoulou, K. Margellos, J. Lygeros, and G. Andersson , Cyber attack in a two-areapower system: Impact identification using reachability , in American Control Conference, 2010, pp. 962–967.[25]
M. Nyberg and E. Frisk , Residual generation for fault diagnosis of systems described by linear differential-algebraicequations , IEEE Transactions on Automatic Control, 51 (2006), pp. 1995–2000.[26]
K. Ogata , Discrete-time Control Systems (2Nd Ed.) , Prentice-Hall, Inc., Upper Saddle River, NJ, USA, 1995.[27]
K. Pan, A. Teixeira, M. Cvetkovic, and P. Palensky , Cyber risk analysis of combined data attacks against power systemstate estimation , IEEE Transactions on Smart Grid, (2018), p. 1.[28]
E. Rakhshani, D. Remon, A. M. Cantarellas, J. M. Garcia, and P. Rodriguez , Virtual synchronous power strategyfor multiple hvdc interconnections of multi-area agc power systems , IEEE Transactions on Power Systems, 32 (2017),pp. 1665–1677.[29]
D. Sahabandu, S. Moothedath, L. Bushnell, R. Poovendran, J. Aller, W. Lee, and A. Clark , A game theo-retic approach for dynamic information flow tracking with conditional branching , in American Control Conference, 2019,pp. 2289–2296.[30]
D. Sahabandu, B. Xiao, A. Clark, S. Lee, W. Lee, and R. Poovendran , Dift games: Dynamic information flowtracking games for advanced persistent threats , in IEEE Conference on Decision and Control, Dec. 2018, pp. 1136–1143.[31]
L. Schenato , To zero or to hold control inputs with lossy links? , IEEE Transactions on Automatic Control, 54 (2009),pp. 1093–1099.[32]
P. Shukla, A. Chakrabortty, and A. Duel-Hallen , A cyber-security investment game for networked control systems ,in American Control Conference, 2019, pp. 2297–2302.[33]
A. Teixeira, S. Amin, H. Sandberg, K. H. Johansson, and S. S. Sastry , Cyber security analysis of state estimators inelectric power systems , in IEEE Conference on Decision and Control, 2010.
ROM STATIC TO DYNAMIC ANOMALY DETECTION 22 [34]
E. E. Tiniou, P. Mohajerin Esfahani, and J. Lygeros , Fault detection with discrete-time measurements: An applicationfor the cyber security of power networks , in IEEE Conference on Decision and Control, 2013.[35]