Structural-Constrained Methods for the Identification of False Data Injection Attacks in Power Systems

Power system functionality is determined on the basis of power system state estimation (PSSE). Thus, corruption of the PSSE may lead to severe consequences, such as disruptions in electricity distribution, maintenance damage, and financial losses. Classical bad data detection (BDD) methods, developed to ensure PSSE reliability, are unable to detect well-designed attacks, named unobservable false data injection (FDI) attacks. In this paper, we develop novel structural-constrained methods for the detection of unobservable FDI attacks and the identification of the attacked buses. The proposed methods are based on formulating structural, sparse constraints on both the attack and the system loads. First, we exploit these constraints in order to compose an appropriate model selection problem. Then, we develop the associated generalized information criterion (GIC) for this problem. However, the GIC method’s computational complexity grows exponentially with the network size, which may be prohibitive for large networks. Therefore, based on the proposed structural and sparse constraints, we develop two novel low-complexity methods for the practical identification of unobservable FDI attacks: 1) a modification of the state-of-the-art orthogonal matching pursuit (OMP) method; and 2) a method that utilizes the graph Markovian property in power systems, i.e. the second-neighbor relationship between the power data at the system’s buses. In order to analyze the performance of the proposed methods, the appropriate oracle and clairvoyant detectors are also derived. The performance of the proposed methods is evaluated on the IEEE-30 bus test case.


I. INTRODUCTION
M ODERN electrical grids are monitored by energy management systems (EMSs). The EMS evaluates the power system state estimation (PSSE) for multiple monitoring purposes, including stability assessment, control, and security [1], [2]. To ensure the reliability of the measurements, residualbased bad data detection (BDD) methods have been integrated into the EMS [1]. However, BDD methods cannot detect well-designed attacks, named unobservable false data injection (FDI) attacks [3], [4]. Unobservable FDI attacks are achieved by manipulating measurements based on the power network topology [3], where the topology matrix is either known or can be estimated from historical data [5]- [8]. Unobservable FDI attacks may inflict severe damage by influencing the PSSE [4], [9], [10]. Therefore, appropriate tools for the detection, identification, and estimation of these attacks, that ensure the reliability of the PSSE, are essential for obtaining high power quality and maintaining stable power system operation.
The problem of detecting unobservable FDI attacks has been considered in the last decade. Methods that strategically protect a basic set of measurements were proposed in [11]- [13]. In the same context, using synchronized phasor measurement units (PMUs) has been suggested [13]- [15]. However, these approaches require integrating additional hardware into the power system, which results in high cost, a long installation period, and additional security vulnerabilities [16]. In addition, they may result in an unobservable system [1], [17]. Detection by the moving target defense technique, where the system configuration is actively changed, has been proposed in [18], [19]. However, in these methods, effective attack detection carries the cost of operating far from the optimal state and, thereby, incurs economic losses [19]. Detection methods based on machine learning and data mining have been suggested in [20]- [23], but these methods require a large set of historic power system data, which is usually unavailable. Detection based on load forecasting has been suggested in [24], yet obtaining a reliable forecast is not ensured and may require extensive resources [25]. In conclusion, all the above methods cannot serve as practical solutions for the detection of unobservable FDI attacks in power systems. Moreover, these methods do not provide attack identification, i.e. localization of the attacked buses.
In this paper, we suggest a novel compressive sensing (CS) parametric approach, based on structural constraints, for the detection of unobservable FDI attacks and the identification of the attacked buses. Sparse recovery CS techniques have become a foremost research area in the last two decades (see, e.g. [26]- [29] and references therein). In power systems, CS algorithms have been used for multi-line outage identification [30], [31], gross error identification [32], identification of imbalances [33], and BDD [34]. CS algorithms designed against unobservable FDI attack have been provided in [35]- [37], where these methods exploit temporal correlations of the power measurements. However, the methods in [35]- [37] require multi-time measurements and do not consider the difference between the structures of the unobservable FDI attack and of the system measurements.
In this paper, we formulate structural and sparse constraints for both the FDI attack and the change in the power system loads between two consecutive time samples. Accordingly, we formulate a model selection problem, where in each model the attack is assumed to be represented by a different set of dictionary elements, i.e. the topology matrix columns. Based on the model selection formulation and the generalized information criterion (GIC) selection rule [38], we develop a structural-constrained method for the detection, identification, and estimation of unobservable FDI attacks. However, the proposed GIC-based method requires an exhaustive search, for which the required computational complexity increases exponentially with the network size, and is therefore practically infeasible. Thus, based on the proposed structural and sparse constraints, we develop two novel low-complexity methods for unobservable FDI attack identification: 1) a modification of the state-of-the-art orthogonal matching pursuit (OMP) method [39]; and 2) a method that is based on the graph Markovian property in power systems, i.e. second-neighbor relationship [40]. Finally, we show by numerical simulation that the proposed GIC, OMP, and Graph Markovian GIC (GM-GIC) methods outperform the detection and identification performance of existing methods. In particular, the proposed parametric methods that are based on information criteria outperform our previous non-parametric detection methods in [41], [42] that are heuristic in nature and are based on the Graph Fourier Transform (GFT).
The main contribution of this paper is threefold. First, we present a new model that takes into account the physical, statistical, and structural properties of unobservable FDI attacks and the typical load changes in power systems. Second, by leveraging the proposed model, we derive a model-selectionbased approach for unobservable FDI attack identification. Finally, we propose an OMP-based method and the novel lowcomplexity GM-GIC method that utilizes both the sparsity and the graphical properties of the problem in order to reduce complexity, while preserving high capabilities for the identification of unobservable FDI attacks. Further, we demonstrate that the GM-GIC method can be used in the general context of sparse recovery of a graph signal from the outputs of a graph filter, and is not limited to power system applications.
The remainder of this paper is organized as follows. In Section II we introduce the necessary background on the power flow equations, BDD, and unobservable FDI attacks. The proposed structural-constrained model for unobservable FDI attacks is presented in Section III. This model is then used in Section IV to develop the GIC-based approach for the detection and identification of unobservable FDI attacks, which is the basis for the two low-complexity methods developed in Section V. Next, a simulation study is presented in Section VI and conclusions in Section VII.
In this paper vectors are denoted by boldface lowercase letters and matrices are denoted by boldface uppercase letters. The operators ||·|| and ||·|| 0 denote the Euclidean norm and the zero seminorm, respectively, where the latter specifies the number of non-zero elements in the vector. The operators (·) T and (·) −1 are the transpose and inverse operators, respectively. The linear space spanned by the A matrix columns is denoted by col(A). For an index set, Λ ⊂ {1, . . . , L}, θ Λ is the |Λ|dimensional subvector of θ containing the elements indexed by Λ, where |Λ| denotes the set's cardinality. For any index sets, Λ 1 and Λ 2 , A Λ 1 ,Λ 2 is the submatrix composed by the columns and rows of A associated with Λ 1 and Λ 2 , respectively.

A. System model
A power system can be represented as an undirected weighted graph, G(V, E), where the set of vertices, V, is the set of buses (substations), and the edge set, E, is the set of transmission lines between these buses. The weight in each line is defined according to the π-model of transmission lines [1]. Hence, the weight over line (n, k) ∈ E is given by the admittance of the transmission line, Y n,k . Specifically, in the direct current (DC) model, where branches are without resistance loss, only the imaginary part of the admittance (the susceptance) is considered as the line weight.
The direct current (DC) power flow model is a linearized representation of the power measurements. This model defines the relation between the active power measurements in the buses and power flows in the transmission lines, z = [z 1 , · · · , z M ] T ∈ R M , and the voltage angles ("states") at the N buses, θ = [θ 1 , . . . , θ N ] T , in the power system network [1]. Based on the DC model, the standard observation model for data injection attacks is given by [43]: The topology matrix, H, is a constant M × N Jacobian matrix, N < M, which is composed of the susceptance elements (as described, for example, in [1], [44]). In addition, e ∈ R M is a zero-mean Gaussian additive noise vector with covariance matrix R, e ∼ N(0, R). The attack is denoted by a ∈ R M . It is assumed that H is a full-rank matrix, i.e. rank(H) = N − 1, and that the system is fully observable.

B. Power system state estimation (PSSE)
The PSSE is commonly used for multiple monitoring purposes, where under the DC model, the system states are the voltage angles, θ . The PSSE is commonly computed using the weighted least squares (WLS) estimator [1]: where In order to robustify the PSSE against errors, classical BDD methods are implemented [1]. These BDD methods are based on the residual error: where the last equality is obtained by substituting the estimated state vector as defined in (2). The residual error from (4) is used, for example, in the Largest Normalized Residual r N maxtest and the χ 2 -test [1], which is given by where H 1 is the hypothesis that the measurements are corrupted by bad data, which may be the result of an attack, and H 0 is the null hypothesis. The threshold, γ BDD , is determined to obtain a desired false alarm probability. However, both the Largest Normalized Residual, which is an identification method, and the χ 2 -test, which is a detection method, can neither detect nor identify the presence of unobservable FDI attacks, as described in the next subsection.

C. Unobservable FDI attacks
The unobservable FDI attack satisfies where c ∈ R N is an arbitrary constant vector. By substituting the well-designed attack from (6) in the observation model from (1) we obtain From a comparison between the models in (1) and in (7), it can be verified that by observing z, the state vector in (1), θ , cannot be distinguished from its corrupted (attacked) version in (7), θ + c, since both θ and c are unknown vectors. As a result, the residual error, obtained by substituting the WLS estimation (which is based on the unobservable FDI model from (7)) in (4) is and therefore cannot be utilized in order to indicate the presence of an unobservable FDI attack. Consequently, all residual-based methods, as well as all the other methods that are based solely on the model in (7), are expected to fail in detecting unobservable FDI attacks. Therefore, methods for the identification of manipulated measurements that integrate additional constraints are required.

III. STRUCTURAL-CONSTRAINED MODELING FOR UNOBSERVABLE FDI ATTACK IDENTIFICATION
Due to the unobservable FDI attack formation and further sparsity restrictions, the additional power caused by the attack is constrained to a small subset of buses in the network. Hence, in this section, a new framework that is based on observing power changes in small subnetworks of the system is proposed. First, the assumptions behind this framework are presented in Subsection III-A. Then, the new structural-constrained model for power measurements in the presence of unobservable FDI attacks is developed in Subsection III-B.

A. Assumptions
The proposed framework, which facilitates structural constraints on the unobservable FDI attack and the typical load changes, is constructed by applying the following definitions and assumptions: A.1 Difference-based model: Similarly to the models in [24] and [45], we assume that two consecutive time samples are observed from the model in (7) at times t and t + 1. The first, z t , is free from malicious attacks, i.e. c t = 0, while the second, z t+1 , may contain an unobservable FDI attack, i.e. c t+1 = c. Thus, we obtain where c = 0 if there is no attack. It should be noted that this assumption is not restrictive since, in practice, the proposed detection method can be integrated into an adaptive scheme in similar manner to the detection method outlined in [45]. The adaptive scheme is initialized by a secure (free of an unobservable FDI attack) measurement, and then the system is constantly monitored by observing the difference between two consecutive measurements. Consequently, only a single measurement used for initializing the adaptive scheme is required to be free of an attack. A.2 Restricted measurements: Since generator buses are heavily secured and obtain direct communication to the control center [46], [47], we assume that generator bus measurements cannot be manipulated. Additionally, it is assumed that 'zero load' measurements cannot be manipulated. Mathematically, this assumption implies that where the set {V \ L} includes the generator and 'zero load' buses, i.e. all buses except the load buses stored in L. A.3 State sparsity: The number of manipulated state variables is bounded by the sparsity term K c , which is considered to be significantly smaller than the cardinality of the node set, |V|, i.e. K c ≪ |V|. Following this sparsity restriction, we define the set that includes all possible supports of the state attack c. As a result, there exists a node subset Λ i ∈ G K c , where i = 1, 2, · · · , |G K c |, that fully contains the attack and satisfies H : It is shown in [3], that this assumption (also used in [37]) stems directly from the commonly-used sparsity restriction on the number of manipulated meters, which states that the attack vector, a, is sparse (see, e.g. [44]). A.4 Typical load changes (quasi-static system): Power systems under normal conditions are quasi-static systems that only change slightly over a short period of time [35], [48]. Therefore, it is considered that typical load changes satisfy where η is a relatively small tuning parameter that can be obtained from the system statistics. A.5 Typical load changes (structural properties): The typical load changes w.r.t. the actual load measurements at a specific moment are determined by the consumption demand [25]. Thus, representing the typical load changes in a matrix form will output a non-sparse vector that is not related to the system topology matrix. As a result, composing the typical load changes requires all of the topology matrix columns. Moreover, a close representation can be achieved only by using the vast majority of the topology matrix columns, which is equivalent to excluding only a small subset of topology matrix columns. Mathematically, we assume that for any Λ i ∈ G K c , there exists a small non-negative parameter ε i such that the orthogonal projection matrix, P ⊥ Λ i , is an (1 − ε i ) ℓ 2 -subspace embedding of H L,V (see Definition 1 in [49], Subsection 2.1): where is the projection matrix onto the space col(H L,Λ i ) and It can be seen that Assumptions A.1-A.3 define constraints on the unobservable FDI attack while Assumptions A.4-A.5, define constraints on the typical load changes.

B. Structural-constrained model
In the following, we use Assumptions A.1-A.5 to construct the new model. First, following Assumption A.1, we consider the difference-based model by taking the difference between two consecutive time samples that satisfy the model in (7), with c = 0 for the first time sample and an arbitrary c in the second time sample. Thus, the difference-based observation model is given by where ∆θ is the change in the state vector, θ , between the two consecutive time samples. Furthermore, by assuming that the noises of the two consecutive time samples are independent, we obtain that the difference-based measurement noise in (16) is distributed according to ∆e ∼ N(0, 2R), where, for the sake of simplicity, it is assumed that 2R = σ 2 e I. According to Assumption A.2, only the changes in the load buses' measurements are relevant for the attack. Thus, the relevant measurement model for the detection of attacks is obtained by taking only the load measurements from (16), which results in where H L,V is the associated submatrix of the topology matrix H, and ∆e L is the corresponding noise. Assumption A.3 implies that the state attack vector support, Λ, is in the set G K c . Hence, as shown in (12), there exists a subset of nodes, Λ i ∈ G K c , that fully contains the attack. By substituting (12) in (17) we obtain Thus, the projection of the measurement vector in (18) onto col(H L,Λ i ) satisfies where (15). In addition, by substituting the following property of projection matrices (see, e.g. p. 46 in [50]): (20) in (14) from Assumption A.5 and rearranging the equation, we obtain Hence, by substituting (21) in (13), we obtain that the projection of H L,V ∆θ onto the column space of the submatrix, H L,Λ i , is bounded by where, according to Assumptions A.4-A.5, ε i η is a small parameter for any i = 1, . . . , |G K c |.
Based on (18) and (22), identifying the subset of attacked buses, Λ i , under the considered model can be formulated as the following multiple hypothesis testing problem: where ε i η is small for any Since εη i is small, we adopt standard practices from the CS literature [26]- [29], where it is common to exclude low amplitude samples from the sparse approximation in order to develop tractable methods. That is, instead of solving (23) one can solve the following modified hypothesis testing problem: where i = 1, . . . , |G K c |. The multiple hypothesis testing in (24) provides a new framework for identifying the sparse state attack vector c from measurements contaminated by additive noise and the nuisance parameter ∆θ . In the modified hypothesis testing in (24), the state attack vector c and the load changes ∆θ are not in the same subspace. Hence, in contrast to the use of (7), under the framework in (24) the attack is observable. Consequently, the formulation in (24) is appropriate for the development of unobservable FDI attack identification algorithms such as the GIC-based identification algorithm developed in Section IV and the low-complexity practical algorithms described in Section V. The performance of the proposed methods is examined w.r.t. the values η and ε i , i = 1, . . . |G K c |, from the hypothesis testing in (23) in Subsection IV-B and empirically in Section VI (see Fig. 3).
A possible interpretation of the problem in (23) and (24) is as a special case of Matched Subspace Detection (MSD) [51], in which the measurements are composed of two deterministic signals and additive noise. In a similar manner to in our framework, one of the two deterministic signals is considered as the target signal (here, the attack) and the other is considered as the background signal (here, the load changes). However, our framework deviates from the standard MSD problem by: 1) the spanning subspace of the target signal is unknown; 2) sparsity restrictions are assumed on the target signal; and 3) the subspace of the target is contained within the linear space that spans the background signal. Thus, standard MSD techniques cannot be applied to solve (23) or (24). A different perspective, presented in Subsection V-B5, is the applicativity of (24) in the context of graph signal processing (GSP), where H is a general graph filter.

A. GIC approach
In the following, we derive the identification of unobservable FDI attacks by selecting the most likely hypothesis in (24), i.e. the most suitable choice of attack support, Λ i , from the set of candidate supports, {G K c ∪ / 0}, given the measurement vector, ∆z L . In order to solve (24), we implement the GIC selection rule [38], which is widely employed in signal and array processing. The GIC method chooses the hypothesis H i which maximizes the sum of two terms: the likelihood term for data encoding, L(·), and a penalty function, τ(·), that inhibits the number of free parameters of the model from becoming very large. For the considered hypothesis testing in (24) and a given difference-based state vector ∆θ , the GIC statistic is given by is the log-likelihood function of ∆z under the ith hypothesis, which is associated with the support, Λ i , andĉ ML|i is the ML estimation of the state attack vector, c. Therefore, in the general case the GIC statistic is a function of ∆θ , which is an unknown deterministic vector. However, as presented in the following, the GIC selection rule in this case is independent of ∆θ and, therefore, is a valid rule. The term τ(|Λ i |, |L|) is a penalty function, which increases with the number of free unknown parameters that is determined in this case by the number of manipulated variables, |Λ i |, and the number of load buses, |L|. In particular, a special case of the GIC family is the Akaike information criterion (AIC), for which the penalty term is Further discussion on the penalty term is provided in the simulation study in Section VI. (25) satisfies

Proposition 1. The GIC statistic in
where const is a constant term that does not depend on the Proof. The proof is provided in Appendix A.
It should be noted that for the null hypothesis in Proposition 1, in which Λ 0 = / 0, we use the convention that P Λ 0 ∆z L = 0. It can be seen from the GIC statistic in (27) that the selected hypothesis is the one that maximizes the sum of the two terms: 1) the projection of the load measurements onto the associated "attack subspace", col(H L,Λ i ), by computing P Λ i ∆z L ; and 2) a penalty function, −τ(|Λ i |, |L|), that encourages a sparse solution. As a result, the representation of the GIC selection rule from (25) by (27) shows that it is not a function of the unknown states, ∆θ , and therefore, it is a valid selection rule. An intuition regarding the GIC selection rule in Proposition 1 can be drawn from a comparison of the GIC statistic in (27) with the bound in (22). This comparison implies that a high energy level for the projected signal in (19), ||P Λ i ∆z L || 2 , cannot be associated with conventional states, and thus provides an indication of the presence of an unobservable FDI attack in the subspace associated with The proposed structural-constrained GIC algorithm for the identification of unobservable FDI attacks is provided in Algorithm 1 and is based on the GIC statistic in Proposition 1. Detection of unobservable FDI attacks is obtained as a by-product of the identification solution of Algorithm 1, s, where the proposed structural-constrained GIC-based detector decides that there is no attack if s = 0 and that an unobservable FDI attack is present for the case where s = 0.

Input:
-difference-based measurements: ∆z -network parameters: H, L -set of candidate state attack supports: G K c -GIC penalty function: τ(·, ·) Output: selected hypothesis: (27) The computational complexity of Algorithm 1 is dominated by Step 2 in the for loop between Steps 1-3, in which the GIC statistic in (27) and K c is the maximal sparsity level from (11). For any k = 1, 2, . . . , K c , the binomial coefficient, |V| k , grows by O(|V| k ) (p. 36, in [28]). Thus, the number of times the GIC statistic is evaluated in Algorithm 1, which is the number of elements in (28), is in the order of O(∑ K c k=0 |V| k ). The number of matrixvector multiplications required for the computation of each GIC statistic in (27) Thus, for large-scale power systems, where |V| and |L| are significantly large compared with K c , the complexity is in the order of O(|V| K c (K c + 1)|L| 2 ), and thus it may be infeasible to use the structural-constrained GIC method. In response to this issue, in Section V, we provide two low-complexity approximations.

B. Performance analysis of the GIC method
In this section, we provide a theoretical performance analysis of the GIC method. In particular, we derive: 1) the oracle GLRT detector [28] that has access to the support of the sparse vector, Λ i ∈ G k c ; and 2) the clairvoyant detector [52], which is a GLRT detector that has full access to the nuisance parameters. While both the oracle and the clairvoyant detectors are infeasible since they are based on unavailable information, they provide insights insights with regard to the influence of the approximation error, ε i η, in the general case, where the support is unknown.
In [53] it is shown that the GIC test provides the same decision rule as the generalized likelihood ratio test (GLRT) for a binary hypothesis testing problem with a fixed threshold. Hence, by using the GIC test in (27), the GLRT for the associated binary (modified) hypothesis testing problem in (24) with only H 0 and a single H i is given by where γ GLRT is the oracle detector threshold, which is set to achieve a desired false alarm probability. The following proposition states upper and lower bounds on the probability of false alarm, P f a , in detecting unobservable FDI attacks with a given support, Λ i , for the state attack vector, c.
Proposition 2. The probability of false alarm of the GLRT in (29) satisfies the following inequality: where Q v (a, b) is the generalized Marcum Q-function of order v [54].
Proof. The proof appears in Appendix B.
The upper and lower bounds in (30) describe the influence of η and γ GLRT on the probability of false alarm of the GLRT detector in (29) and are dictated by the property 0 ≤ ε i ≤ 1 shown in Appendix B. Since Q v (a, b) decreases as b decreases [54], we obtain the expected result that both the upper and lower bounds of P f a decrease as the threshold γ GLRT decreases. In contrast, since the function Q v (a, b) increases as a increases [54], the upper bound on the r.h.s. of (30) increases as η increases. This result indicates that when the bound on the nuisance parameter, H∆θ , in (13) is tighter, we can guarantee a tighter upper bound on the probability of false alarm, P f a . Moreover, under Assumption A.5, ε i from (21) is expected to be small for all i = 1, 2, . . . , |G K c |. Thus, following the derivations in Appendix B, the probability of false alarm is expected to be close to the lower bound in (30). In particular, the extreme case of ε i = 0 occurs when the load changes are outside of the column space col( For this case, the probability of false alarm achieves the lower bound in (30), i.e.
and the probability of detection (see (73) in Appendix B) is Hence, following the discussion on the properties of the generalized Marcum Q function in Appendix B, the probability of detection in (31) increases as the attack energy, ||H L,Λ i c Λ i || 2 , increases.
In the following, we analyze the performance of the associated clairvoyant detector. The clairvoyant detector in the general case is the detector obtained by using the otherwise unknown signal parameter values of the composite hypothesis testing (see, e.g., Chapter 6.5, [52] and [55]). Therefore, the clairvoyant detector can serve as a reference for the achievable performance of practical detectors. In the considered model, the clairvoyant detector assumes perfect knowledge regarding the nuisance parameter, H L,V ∆θ . Thus, similar to the GLRT in (29), the clairvoyant GLRT is the GLRT detector for the binary (modified) hypothesis testing problem in (24) with only H 0 and a single H i , but now when H∆θ is assumed to be known. That is, the clairvoyant GLRT is given by where γ cl is the clairvoyant detector threshold. It should be noted that the clairvoyant detector is a theoretical detector that provides an upper bound on the performance of the probability of detection of practical detectors.
The following proposition presents a special case in which the proposed GLRT in (29) achieves the clairvoyant GLRT. Proposition 3. If ε i = 0, then the GLRT in (29) coincides with the clairvoyant GLRT in (32).
Proof. In Appendix B it is shown that ε i = 0 implies that ||P Λ i H L,V ∆θ || 2 = 0. Hence, as a result of known norm properties we obtain that P Λ i H L,V ∆θ = 0. Thus, by substituting this result in the clairvoyant GLRT in (32) we obtain the proposed GLRT in (29).
As a result of Proposition 3, if ε i = 0, then the detection of the attack is not influenced by the presence of the nuisance parameter, H∆θ . In other words, the value of H∆θ does not affect the performance of the proposed GLRT detector in (29) for the observed case.
Similar to the clairvoyant GLRT, we define the clairvoyant GIC as the GIC selection rule for the hypothesis testing in (23) that assumes perfect knowledge of the nuisance parameter, H L,V ∆θ . Therefore, similar to the derivation of (27), it can be shown that the clairvoyant GIC is given by The performance of the clairvoyant GIC is used as a benchmark on the performance of the proposed methods, as demonstrated in the simulations.

V. LOW-COMPLEXITY IDENTIFICATION METHODS
As shown at the end of Subsection IV, the computational complexity of the structural-constrained GIC method makes it impractical for large power system networks, where |V| and |L| are large. In this section, we develop two low-complexity methods that rely on Assumptions A.1-A.5: 1) an OMP-based method [39]; and 2) a novel method that exploits the graph Markovian property of order two between the nodes (buses) in the graph representation of the power system [40].

A. OMP method
The OMP algorithm in [39] is an efficient method for the recovery of sparse signals. The basic principle behind the OMP algorithm is to iteratively find the support set of the sparse vector. The OMP method proceeds by finding the column of the CS matrix that correlates most strongly with the signal residual. The residual is constructed in each iteration by projecting the measurements onto the linear space spanned by the remaining columns that were not selected in previous iterations.
In this subsection, we apply the OMP algorithm for the sparse recovery of the state attack vector, c, which is a sparse signal as described in Assumption A.3, from the measurements in (17). It should be noted that the measurement model in (17) contains a nuisance parameter vector, H L,V ∆θ , which is not a part of the conventional sparse recovery model. Nevertheless, similarly to in the derivations in Subsection III-B, it can be shown that under Assumptions A.1-A.5 the nuisance parameter has a negligible effect on the OMP selection criterion. Thus, the conventional OMP method is valid for this setting.
The main iteration of the OMP algorithm, which is performed on the load measurements model from (17), is given as follows. Suppose Λ ( j) is the estimated support set of c in the jth iterative step. In the j + 1th iteration, we compute where P˜k is obtained by replacing Λ i =k in (15), and is the signal residual at the jth iteration. By substituting (35) and P ⊥ Λ ( j) = I − P Λ ( j) on the r.h.s. of (34) we obtain By limiting the number of iterations to the maximal sparsity level, K c , we obtain that Λ ( j) ∈ G K c . In addition, any single bus, k ∈ V, is an element in G K c , i.e. k ∈ G K c , as well. Therefore, under Assumptions A.1-A.5, we obtain from (22) that the nuisance parameter, H L,V ∆θ , has a minor effect on both the terms P {k} ∆z L and P Λ ( j) ∆z L . Consequently, the nuisance parameter has a minor effect on (36) and, thus, its influence on the OMP selection procedure in (34)  update: r OMP|( j) = P ⊥ Λ ( j) ∆z L 8: end for 9: returnΛ OMP =Λ ( j) structural-constrained OMP algorithm for the identification of unobservable FDI attacks is provided in Algorithm 2.
The computational complexity of Algorithm 2 is dominated by Step 3 (computing (34)), which requires O(2|L| 2 + 3|L|) matrix-vector multiplications for each k = 1, 2, . . . , |V|. The loop between Steps 2-8 is performed at most K c times. As a result, the total complexity of Algorithm 2 is in the order of O(|V|K c |L| 2 ), which is significantly lower than the complexity of the GIC method (see discussion after (28)).
In general, the OMP method is used in a variety of applications due to its low computational complexity. However, the OMP method is a greedy algorithm with no optimal recovery guarantees, and usually requires an incoherent dictionary in order to provide high performance [27]. In our case, the CS matrix is the topology matrix H L,V , which may be highly correlated, and thus, with large mutual coherence. In order to address this issue, in the following subsection, we develop a novel low-complexity method that uses the power system graph representation Markovian properties.

B. Graph Markovian GIC (GM-GIC)
In this section, we develop the low complexity GM-GIC method for the measurement model in (17) and under Assumptions A.1-A.5, which considers the graphical representation of the power system. Accordingly, in this subsection the system buses are referred to as the graph nodes. The power system graphical representation is utilized in Subsection V-B1 to analyze the affect of an unobservable FDI attack on the measurements. Based on this analysis, we derive the GM-GIC method, which consists of the following four stages: 1) lowscale pre-screening (Subsection V-B2); 2) node partitioning (Subsection V-B3); 3) local GIC, on the partitioned subsets; and 4) sparsity correction. The GM-GIC method, including Stages 3 and 4, is summarized in Subsection V-B4. Furthermore, in Subsection V-B5 it is demonstrated that the GM-GIC method can be applied in GSP applications with sparse signals, in addition to its application for FDI attack identification in power systems.

1) Attack analysis:
In the following, we analyze the effect of an unobservable FDI attack on the system measurements. Specifically, we focus on how an attack, with a support Λ, affects the single-node measurement subspace, col(H L,m ), for any node m ∈ V. Considering the measurement model in (17), the unobservable FDI attack can be linearly decomposed as follows: Based on the decomposition in (37), we can analyze the influence of an attack on a single-node measurement subspace by summing over the individual influences on each node, k ∈ Λ. The following proposition evaluates the effect of a single-node attack on a single-node measurement subspace.
Proposition 4. The projection of a single-node attack on node k ∈ Λ, H L,k c k , onto the measurement subspace associated with node m ∈ V, col(H L,m ), satisfies and where P {m} is obtained by substituting Λ = {m} in (15).
Proof. The proof appears in Appendix C.
This proposition demonstrates that the single-node attack obtains a 'local' effect. That is, a single-node attack on node k does not affect the single-node measurement subspace of node m if the geodesic distance (hop distance) [56] between nodes k and m, d(k, m), is greater than two, d(k, m) > 2. By multiplying (37) by P {m} and substituting (38), we obtain As a result of (40), the single-node measurement subspaces affected by the unobservable FDI attack are only those which are associated with nodes in the set where, according to (41), the set A includes the attacked nodes, those in Λ, and first-or second-order neighbors of attacked nodes. This observation can be interpreted as a second-order graph Markov property [40]. In general, the power system network generates a sparse graph where nodes are connected to 2-5 neighbors [57]. Thus, the nodes in the power system have at most 5 first-order and 25 second-order neighbor nodes. As a result, the number of nodes affected by the attack is bounded: where K c is the sparsity term defined in Assumption A.3 and A is defined in (41). Moreover, the power system network diplays a local behavior for the connectivity pattern, where only substations geographically close are likely to be connected. Thus, the bound in (42) is not a tight bound, and |A| is expected to be significantly lower than the r.h.s of (42). In conclusion, it is implied that for large networks |A| ≪ |V|.
2) Pre-screening: The first stage of the GM-GIC method is a pre-screening stage. In order to evaluate if the dictionary matrix column associated with node m ∈ V, H L,m , is correlated with the attack, Hc, we derive the GLRT (in a similar manner to in (29)) for the associated binary hypothesis testing problem in (24) with only H 0 and a single H i , selected as H i = {m}, which results in where ρ is set to determine a desired false alarm rate. Applying the GLRT detector in (43) for any node in the system yields the following set of suspicious nodes: Thus, the node m will be included in the set S if abnormal energy is detected in the single-node measurement subspace, col(H :,m ). By considering the hypothesis testing in (24) in which Λ i = {m}, the projection of the nuisance parameter, P {m} H L,V ∆θ , is negligible. Therefore, the inclusion of the node m in the set S indicates that the measurement subspace, col(H :,m ), is affected by the attack, i.e. that P {m} H L,Λ c Λ = 0. Thus, from (40), the set of suspicious nodes in (44) can be interpreted as an estimator of the set A in (41). From (39), it can be seen that the norm of the projected single-node attack on single-node measurement subspaces of attacked nodes is higher than the norm of the projection onto measurement subspaces associated with other nodes. Thus, considering the discussion after (37), the effect of the attack, H L,Λ c Λ , on attacked nodes is expected to be higher than its effect on their first-and second-order neighbor nodes. Therefore, even if the estimator S does not cover all the nodes in A, the attacked nodes are still expected to be included.

3) Node partitioning:
The second stage of the GM-GIC method is the partitioning of the suspicious set, S, in (44), into disjoint subsets. This partitioning provides the condition that an attack on nodes located in one of the subsets does not affect the measurement subspaces associated with nodes in the other subsets. According to the proposed partitioning, a general set S is partitioned into Q disjoint subsets, {S q } Q q=1 , if the following is satisfied: for any two different subsets, S q and S p , p = q, selected from (40) and (45), an attack on any node subset in S q , Λ q ∈ S q , does not affect the column space associated with the nodes in any of the other subsets, e.g. col(H L,S p ), i.e.
Therefore, a node partitioning that satisfies (45) enables identification of the attacked buses in each subset separately, and, thus, reduces the problem dimension. In contrast, the GIC from Algorithm 1 considers the entire node set V.
In practice, finding a node partitioning that satisfies (45) is performed as follows. First, based on the graph representation of the given power system, we generate the unweighted undirected graphG = (S, EG), where S are the nodes and the set of edges is defined as (47) in which d(·, ·) refers to the geodesic distance measured on the original graph that represents the power system network. Then, we find the connected components ofG (e.g. by using the Matlab command conncomp). According to the definition in (47), we obtain that d(k, m) > 2 for any k and m that belong to different connected components ofG. Therefore, selecting the partition subsets {S} Q q=1 to be the node sets of the connected components ofG satisfies (45).

4) Summary: GM-GIC method:
In this subsection we summarize the proposed GM-GIC method. In the first stage of this method, a reduced set, composed of suspicious nodes, S, is extracted from the node set V, according to (44). In the second stage, the set S is partitioned into Q disjoint subsets, {S q } Q q=1 , as described in Subsection V-B3. In the third stage, the GIC method is applied on each subset separately, by replacing the set of candidate state attack supports, G K c , in Algorithm 1 with As a result, for each subset S q we obtain a partial estimation of the support set, which is denoted byΛ q . The total support set of the state attack vector is the union of all the partial estimatesΛ The estimated support in (49) may exceed the sparsity limit |Λ GM-GIC temp | ≤ K c . Hence,Λ GM-GIC temp may include node-elements which are not attacked. It should be noted that in this stage, Assumption A.2 is relaxed and the sparsity restriction is only imposed on the partial (separated) estimates, i.e. Λ q ≤ K c . In order to reduce the identification errors that may be induced by this relaxation, in the fourth stage, the estimated support in (49) is corrected to satisfy |Λ GM-GIC | ≤ K c . This stage is performed by: 1) evaluating the state attack ML estimation, c ML|î Λ GM-GIC , by replacing Λ i with (49) in (65) and 3) preserving only the K c elements with the highest absolute values,Λ The proposed GM-GIC algorithm for the identification of unobservable FDI attacks is provided in Algorithm 3. The computational complexity of the different stages of the GM-GIC method in Algorithm 3 are as follows: 1) the pre-screening stage in Step 1 requires O(|V|(2|L| 2 + 3|L|)) matrix-vector multiplications; 2) the node partitioning stage in Step 3 is dominated by the graph partitioning, which is implemented by the conncomp(·) Matlab command and requires O(|S| + |S| 2 ) computations (considering that the number of edges inG cannot be anticipated in advance, this analysis is returnΛ GM-GIC = / 0 4: end if 5: Node-partitioning: 1) generateG = (S, EG) by computing (47); 2) partitionG into its connected components (e.g. by conncomp in Matlab); and 3) set {S q } Q q=1 to be these components 6: for q = 1 . . . Q do 7: generate the graph G  (65) 13: correctΛ GM-GIC by (50) and (51) 14: end if 15: returnΛ GM-GIC based on the worst case in whichG is a fully connected graph); 3) in Steps 4-7 the GIC method is applied on each of the subsets {S q } Q q=1 , so that, from Subsection IV, the complexity is in the order of O((∑ Q q=1 |S q | K c )(K c + 1)|L| 2 ); and 4) the sparsity correction stage, in Steps 9-12, is dominated by Step 10, which requires O(K c (K c + 1)|L|) matrix-vector multiplications. In power systems, the number of load nodes (buses) |L| is, in general, in scale with the number of nodes |V|. Thus, from Subsections V-B1-V-B2, we obtain that |S| ≪ |L| and the total computational complexity of Algorithm 3 is O((∑ Q q=1 |S q | K c )(K c + 1)|L| 2 ). The computational complexity of the different methods is summarized in Table I. It can be seen that the computational complexity of the GM-GIC method is significantly lower than that of the complexity of the GIC method, but may be higher than that of the complexity of the OMP method. Finally, the worst complexity for the GM-GIC method is obtained when S cannot be partitioned. In this case, the complexity is in the order of O(|S| K c (K c + 1)|L| 2 ), which is still significantly lower than the computational complexity of the GIC method.
In contrast to the GIC and the OMP methods, the computational complexity of the GM-GIC method also depends on the underlying structure of the power system (and not only the dimensions), which is represented by the topology matrix H. In particular, the results of the pre-screening and node partitioning stages, S and {S q } q , which are detailed in Subsections V-B2 and V-B3, respectively, depend on the formation of the topology matrix H. The method will perform effectively when the properties detailed at the end of Subsection V-B1 are satisfied.

5) Sparse signal recovery in general GSP applications:
The proposed GM-GIC heuristic is designed to exploit the graphical properties of power system networks, which are detailted in Subsection V-B1 as part of the attack analysis, and include: 1) graph Markovity; 2) graph sparsity; and 3) local graph connectivity behavior. Thus, the GM-GIC method can be utilized for sparse recovery in GSP frameworks, where, if Properties 1)-3) are satisfied, the GM-GIC method has advantages, compared with conventional sparse recovery algorithms, in terms of computational complexity and accuracy of the sparse recovery.
The emerging field of GSP provides new methodologies for the analysis of signals in applications with underlying relations that could be modeled by a graph [58]- [60]. In particular, the propagation of a process that originates from a sparse input through the graph vertex domain is commonly modeled in the GSP literature as an output of a graph filter [61]- [63]; this approach has various applications, such as locating the source of a disease [64], [65] or monitoring anomalies in sensor networks [66]. The proposed GM-GIC method can be applied to this problem of sparse recovery in GSP models as follows. Let us assume the graph G = (V, E) with the set of nodes (vertices) V and the set of edges E. We consider the recovery of a sparse graph signal, x ∈ R |V| , s.t. ||x|| 0 ≪ |V|, from the noisy output of a graph filter, F ∈ R |V|×|V| : where e ∼ N(0, σ 2 n ) is the system noise, and the graph filter F is linear and shift-invariant [58]- [61]. Hence, F is a polynomial in a graph shift operator (GSO), S ∈ R |V|×|V| , as follows [58]- [60]: where h 0 , . . . , h Ψ are the filter's coefficients and Ψ ≥ 1 is the filter's order. The GSO, S, is defined as a matrix that satisfies ∀k, m ∈ V, where d(k, m) is the geodesic distance between nodes k and m. In particular, it can be verified that (54) implies ∀k, m ∈ V, which leads to It can be seen that the problem in (24) is reduced to the problem in (52) when the topology matrix H is replaced with the graph filter F and when the nuisance parameter is absent, ∆θ = 0. It can be seen that for a graph filter of order Ψ = 1, (56) is identical to the property in (78)-(79) from Appendix C, where in this case F = H V,L . The property in (78)-(79) is the basis for Proposition 4 that enables the derivation of the GM-GIC method. Thus, for this case, the GM-GIC method can be implemented (without any change) on (52) for the recovery of the sparse graph signal, x, which replaces the state attack vector, c, and without nuisance parameters. For the general case, where Ψ > 1, the GM-GIC method can be applied on (52) after replacing the constraint in (47) (that is part of the node partitioning stage in Algorithm 3) with the constraint 1 ≤ d(n, k) ≤ 2Ψ. For this case, in which (56) is considered instead of (79), the condition in (38) from Proposition 4 is changed to d(k, m) > 2Ψ.
VI. SIMULATIONS In this section, the performance of the proposed methods is demonstrated for the tasks of detection and identification of unobservable FDI attacks. The simulations are conducted on the IEEE-30 bus test case, where the topology matrix and measurement data are extracted using the Matpower toolbox for Matlab [67]. The simulation set-up is described in Subsection VI-A, while the proposed methods' performance is demonstrated in Subsection VI-B.

A. Set-up 1) Measurements:
The simulation study is conducted on the difference-based model in (16). For the sake of simplicity of implementation, we ensure Assumption A.2 by defining the setṼ to include only state variables that are related to load measurements, and then restricting the support of the state attack vector to Λ ⊆Ṽ. In particular, for the IEEE-30 bus test case it can be verified thatṼ = {14, 16, 17, 18, 19, 20}. In each simulation, we set the cardinality of the state attack support vector to be |Λ| = K a and draw the support vector, Λ, uniformly from the set {Λ ∈Ṽ : |Λ| = K a }. The attack values on the chosen support are randomly drawn from the uniform distribution over the interval [−1, 1]. Then, the values are scaled to obtain a desired value of attack norm, ||a||.
The difference-based state vector is obtained by first generating the state vectors at times t and t + 1, θ t and θ t+1 , respectively, and then subtracting θ t from θ t+1 to obtain ∆θ = θ t+1 − θ t . The state vector at time t, θ t , is set to the values in the IEEE test case [67]. Based on θ t and considering Assumption A.5, the state vector at time t + 1, θ t+1 , is simulated by first updating the load measurements with random scaling: and then computing θ t+1 by the rundcpf(·) Matpower command in [67]. In addition, throughout the simulations the noise is modeled as in (7) with σ 2 e = 0.01. 2) Methods: The proposed methods include the structuralconstrained GIC method from Section III, and the lowcomplexity methods in Section V: the structural-constrained OMP method and the GM-GIC method. The penalty function for the GIC and GM-GIC methods defined in (27) is set by where ζ and γ GIC are user-defined regularization parameters and δ [·] is the Kronecker delta function. The term ζ |Λ i | is set to encourage sparse solutions for the identification problem, while the term γ GIC δ (|Λ i |) is set to maintain a desired false alarm rate for the detection problem. In particular, the probability of false alarm is set to P FA = 0.05 unless stated otherwise. Furthermore, by selecting the AIC in (26), the GIC tuning parameter in (58) is set to ζ = 2. The maximal sparsity rate is set to K c = 6. The performance of the proposed methods is compared with the following methods that were all modified according to the difference-based model in (16): The sparse optimization technique in [37], which considers an ℓ 2 relaxation for the attack sparsity restriction (denoted as ℓ 2 ). M. 2 The sparse optimization technique in [35], which considers an ℓ 1 relaxation for the attack sparsity restriction (denoted as ℓ 1 ).
All numerical results in this section were obtained using at least 500 Monte Carlo simulations. The detection thresholds were computed from simulated historic data obtained by 500 off-line simulations of (16) (57) is set to σ 2 s = 0.05. This figure shows that the detection performance of the proposed low-complexity methods, GM-GIC, and OMP, converges to the performance of the optimal GIC method. Furthermore, the proposed GIC, GM-GIC, and OMP methods provide a higher detection probability for any chosen false alarm rate, when compared to the previous methods. It should be noted that the sparse ℓ 2 and ℓ 1 methods were developed based on models with multiple time measurements and where the sparsity pattern of the attack is defined over both the time and the bus domains. This is in contrast with the considered model, in which only two time samples are provided and where the attack is only sparse in the bus domain. Finally, as expected from Subsection II-C, the BDD method cannot detect unobservable FDI attacks and is no better than flipping a coin.
2) Identification: The identification performance is measured by the capability to classify each state variable as manipulated or not. Therefore, we evaluate the identification performance by the F-score classification metric [68]: where Λ is the true support of the state attack vector, c, and Λ is the estimated support by a given method. The terms t p , f p , and f n , are the true-positive, false-positive, and falsenegative probabilities, respectively. The F-score metric takes values between 0 and 1, where 1 means perfect identification. In Fig. 2, the identification performance, measured by the Fscore metric, of the GIC, GM-GIC, OMP, ℓ 2 , and ℓ 1 methods is examined w.r.t the attack characteristics. Methods M.3-M.5 are not included in the comparison since they provide solely detection and do not identify the attacked buses' locations. In Fig. 2.a the F-score is presented versus the normalized attack norm, ||a|| K a , for K a = 4 and σ 2 s = 0.1. In addition, in Fig. 2.b the F-score is presented versus the number of attacked elements, K a , for ||a|| = 1.2 and σ 2 s = 0.05. It can be observed from both figures that the proposed methods have a significantly higher F-score than those of the previous sparsity-based methods, ℓ 2 and ℓ 1 . In addition, it can be seen from Fig. 2.a that the performance of all the methods improves with the increase of the attack measured by ||a|| K a . Similarly, Fig. 2.b shows that the performance of all the methods degrades as K a increases. Nonetheless, the F-score of the proposed methods is above 0.8 even when a considerable portion (> 0.2) of the system is attacked. Finally, it can be seen that both the GM-GIC and OMP methods show relatively close results to those of the GIC method, where the GM-GIC method provides higher identification rates than the OMP method.
3) Identification under mismatch: The use of the GIC approach for identifying and detecting unobservable FDI attacks is based on the assumption that the nuisance parameter H∆θ is bounded and small, as described in (13). In Fig. 3, we illustrate the robustness of the performance to this assumption by evaluating the influence of the norm of the nuisance parameter, ||H∆θ ||, on the identification performance of the different methods. To this end, we consider the worst case, in which the equality in (13) holds, i.e. ||H L,V ∆θ || 2 = η. In particular, the F-score of the different methods, as well as of  (33), is compared for different rates of η for K a = 4 and ||a|| = 0.05K a . Fig. 4 shows that the F-score of the different methods, excluding the clairvoyant GIC from (33), decreases as η increases, due to the model mismatch. The F-score of the clairvoyant GIC is independent of η since it uses the true value of ∆θ and is not based on the approximated model, defined by η. It can be seen that this degradation in the identification performance is significantly greater for the sparsity-based methods, ℓ 2 and ℓ 1 , than for the proposed methods. Therefore, the proposed methods are more robust to a mismatched model in comparison to the existing methods. Finally, for small values of η the proposed GIC method converges to the clairvoyant GIC.
4) Run-time: In Fig. 4, the averaged run-time of the identification methods is presented versus the ratio of node elements attacked, K a |V| , for σ 2 s = 0.05 and ||a|| = 1.2. It can be seen that the GIC method requires a run-time which is significantly higher than those of the other methods. This makes the GIC method impractical for large networks. The CS methods [35], [37] have a higher averaged run-time than the proposed lowcomplexity methods, GM-GIC and OMP. As expected from the discussion in Subsection V-B4, the GM-GIC method requires a  Fig. 3. Unobservable FDI attack identification: The F-score of the different methods with K a = 4 and ||a|| = 0.05K a . significantly lower run-time than the GIC method, but a higher run-time than the OMP method.

VII. CONCLUSIONS
In this paper, we introduce novel methods for the identification of unobservable FDI attacks in power systems. Unlike classical BDD methods that depend on measurement residuals and fail to detect unobservable FDI attacks, the proposed methods successfully detect and identify these attacks by utilizing power system intrinsic structural information. Specifically, the proposed methods are based on defining structural constraints on both the attack and the typical load changes, and then formulating the identification problem as a model selection problem. We develop the GIC selection rule for the identification task, as well as two low-complexity methods: 1) the structural-constrained OMP method, which is a modification to the standard OMP method that accounts for the proposed structural and sparse constraints; and 2) the novel GM-GIC method that exploits the graph sparsity, graph Markovity, and the local behavior of graph connectivity in power systems. As by-products of the identification process, the proposed methods also enable the detection of unobservable FDI attacks and PSSE in the presence of these attacks. The low-complexity methods also enable fast implementation that can be integrated into an adaptive scheme that performs detection continuously to see if an attack is present at any given time. In addition, we show the relations between the assumed problem and the problems of sparse recovery of graph signals and MSD. In particular, the proposed GM-GIC method can be applied for denoising an output of a general low-order GSP filter with a sparse input. This problem arises in a variety of GSP applications, such as epidemiology and anomaly detection in sensor networks. In addition, performance analysis of the proposed approach is investigated by introducing a clairvoyant and an oracle detector of the considered FDI detection problem. Our numerical simulations show that the proposed methods outperform existing methods for detection and identification, and are robust to the model assumptions. Moreover, the GM-GIC method, that integrates the structural information regarding the underling graph behind the data, obtains better performance than that of the OMP method and the existing sparsity-based methods. The OMP and GM-GIC methods require shorter run-time than the GIC method. Finally, the performance of the proposed approaches is compared with that of derived clairvoyant and oracle detectors. (60) Thus, the ML estimator of the state attack vector, c Λ i , is obtained by maximizing (60) w.r.t. c Λ i . Since (60) is a concave function w.r.t. c Λ i , the ML estimator can be computed by equating the derivative of (60) w.r.t. c Λ i to zero, as follows: which implies that the ML estimator iŝ By using the projection matrix from (15), P Λ i , and its orthogonal projection, P ⊥ Λ i , we can decompose H L,V ∆θ as follows: By substituting the constraint from (24), P Λ i H L,V ∆θ = 0, in (63), one obtains Substitution of (64) in (62) results in where the last equality is obtained from the properties of projection matrices. The generalized log likelihood is obtained by substituting (64) and (65) in (60): Additionally, from the properties of projection matrices, it can be verified that (see, e.g. p. 46 in [50]) By substituting the constraint from (24), P Λ i H L,V ∆θ = 0, in (67) we obtain ∀Λ i ∈ G K c . By substituting (68) in (66)  is independent of the hypothesis, i.e. is not a function of Λ i , maximizing the GIC rule from (70) w.r.t. i = 0, 1, . . . , |G K c |, is equivalent to maximizing the r.h.s. of (27) w.r.t. the same candidate sets, i = 0, 1, . . . , |G K c |.

APPENDIX B ORACLE GLRT
While the GLRT in (29) was developed for the binary modified hypothesis testing problem in (24), , it should be analyzed based on the binary version of the original hypothesis testing in (23) with only H 0 and a single H i . In this case, the GLRT detector is distributed as follows (see, e.g. Appendix 7B in [52]): where χ 2 r (λ ) denotes a non-central χ-square distribution with r degrees of freedom and a non-centrality parameter λ . From (71) it can be verified that the probability of false alarm and the probability of detection of the GLRT from (29) are given by (see Subsection 2.2.3 in [52]): and respectively, where Q v (a, b) is the generalized Marcum Qfunction of order v [54]. From [54], Q v (a, b) strictly increases as a increases. That is, for all a 1 ≥ 0 and a 2 , v, b > 0. Thus, the probability of false alarm in (72) strictly increases as ||P Λ i H L,V ∆θ || increases. Consequently, by considering the worst case in (22) of ||P Λ i H L,V ∆θ )|| 2 = ε i η, we obtain that the probability of false alarm in (72) strictly increases as ε i η increases.
The following proposition defines upper and lower bounds on the approximation error ε i . Proposition 5. The approximation error from (21) can be bounded by 0 ≤ ε i ≤ 1.

APPENDIX C SINGLE-NODE ATTACK
In Proposition 1 in [40] it was shown that under the assumption that the injected powers in the nodes are Gaussian distributed and mutually independent, which holds under the model in (16), the columns of the topology matrix satisfy Thus, by substituting (79) in (80) we obtain (38). In addition, by using the properties of projection matrices it can be verified (see Theorem 2.22 in [50]) that ||P {m} H L,k c k || ≤ ||H L,k c k ||, ∀k, m ∈ V, and P {k} H L,k c k = H L,k c k .