Observability Robustness under Sensor Failures: a Computational Perspective
aa r X i v : . [ m a t h . O C ] J a n Observability Robustness under Sensor Failures:Complexities and algorithms
Yuan Zhang, Yuanqing Xia, and Kun Liu
Abstract —The problem of determining the minimal number ofsensors whose removal destroys observability of a linear timeinvariant system is studied. This problem is closely related to theability of unique state reconstruction of a system under adversarialsensor attacks, and the dual of it is the inverse to the recentlystudied minimal controllability problems. It is proven that thisproblem is NP-hard both for a numerically specific system, andfor a generic system whose nonzero entries of its system matricesare unknown but can take values freely (also called structuredsystem). Two polynomial time algorithms are provided to solvethis problem, respectively, on a numerical system with boundedmaximum geometric multiplicities, and on a structured systemwith bounded matching deficiencies, which are often met bypractical engineering systems. The proposed algorithms can beeasily extended to the case where each sensor has a non-negativecost. Numerical experiments show that the structured system basedalgorithm could be alternative when the exact values of systemmatrices are not accessible.
Index Terms —Observability robustness, secure estimation, opti-mization, computational complexity, recursive tree
I. INTRODUCTIONModern control systems are often equipped with computer-based components, which could be vulnerable to cyber attacks[1]. As such, security has become a more and more importantissue in control and estimation of cyber-physical systems, suchas chemical processes, power grids and transportation networks[2–5]. Typical attacks could be imposed on actuators, sensors orcontrollers of a control system. Among the related issues, oneproblem has attracted researchers’ interest, which is that, underwhat condition it is possible to uniquely reconstruct the statesof a system by observations of system outputs in the presenceof adversarial sensor attacks [2, 3, 6]. To be specific, considerthe following linear time invariant (LTI) system under attack: ˙ x ( t ) = Ax ( t ) , y ( t ) = Cx ( t ) + e ( t ) , (1)where A ∈ R n × n , C ∈ R r × n , x ( t ) ∈ R n , y ( t ) ∈ R r arerespectively the state and output vectors, and e ( t ) ∈ R r is theattack signal injected to sensors. The vector e ( t ) is sparse in thesense that only the components corresponding to the attackedsensors can be nonzero, t ≥ , but identities of the attackedsensors are not known.It is found in [2, 6] that, unique state reconstruction ofsystem (1) under s ∈ N attacks (i.e., the maximum number ofnonzero entries in e ( t ) over all t ) is possible, if and only system(1) remains observable after the removal of arbitrary sets of s sensors. Such condition is called s -sparse observability in [3], and has been extended to nonlinear systems and thedistributed secure state estimation scenarios [4, 5, 7]. For agiven system and s ∈ N , verifying the aforementioned condition The authors are with the School of Automation, Bei-jing Institute of Technology, Beijing, China. (email: { zhangyuan14,xia_yuanqing,kunliubit } @bit.edu.cn ). seems combinatorial in nature. To the best of our knowledge,however, no literature has systematically discussed its compu-tational complexity or provided the associated (approximation)algorithms. This motivates our present work on observabilityrobustness under sensor failures.More precisely, we will consider the following problem: givenan LTI system, determine the minimal number of sensors whoseremoval makes the resulting system unobservable. We will callthis problem minimum sensor removal observability problem(minSRO) . Here, sensor removal could also be regarded as theresult of denial-of-service attacks on a specific set of sensors [1],which means that the attacked sensors cannot output any systeminformation.On the other hand, the dual of the minSRO, i.e., determiningthe minimal number of actuators (inputs) whose removal destroyssystem controllability, is the inverse to the minimal controllabil-ity problem (MCP) , which can be formulated as selecting theminimal number of inputs from a given set of inputs to ensuresystem controllability [8, 9]. The MCP is proven to be NP-hard, and can be approximated gracefully using simple greedyalgorithms [8, 9]. An extension of the MCP is to determine thesparsest input matrices such that the system remains controllableagainst a desired level of actuator failures [10]. Notice that thisproblem is different from the minSRO, as the former is a designproblem where the input matrices are to be determined, whilethe latter is to measure robustness of a given system.Robustness of controllability and observability of anetworked system under structural perturbations (such asedge/node/actuator/sensor removals or deletions) have beenstudied in several works [11–16]. For example, [11] hasdiscussed observability preservation under sensor removals, and[12] has studied controllability preservation under simultaneousfailures in both the edges and nodes. Controllability robustnessis also measured by the number of additional inputs neededfor controllability under edge/node removals in [15, 16].The above works mainly focus on the classifications ofedges and nodes according to effects of their failures onsystem controllability/observability, rather than the associatedoptimization problems. A similar problem to the minSRO understructured system framework with ununiform actuator costshas been proven to be NP-hard in [14]. However, no efficientalgorithm was provided therein.In this paper, it is proven that the minSRO is generally NP-hard, thus confirming the conjecture implicitly made in [2, 3] thatverifying the s -sparse observability condition is computation-ally intractable. Nevertheless, a polynomial time algorithm basedon recursive trees for solving this problem is provided on systemswith bounded maximum geometric multiplicities (often met bypractical engineering systems). Particularly, this result indicatesthat, in contrast to the MCP, the computational intractability ofthe minSRO is essentially caused by the increase of maximum geometric multiplicities of system matrices, rather than thatof numbers of system states or sensors. Considering that thisalgorithm requires accurate eigenspace decompositions of systemmatrices, we then focus on the minSRO under the structuredsystem framework. It is proven that such problem is still NP-hard for the structured systems. Afterwards, a polynomial timealgorithm is provided on structured systems with bounded match-ing deficiencies. This algorithm is further accelerated by theDulmage-Mendelsohn decomposition. Numerical experimentsshow that for randomly generated systems, the structured systembased algorithms could be good alternatives for observabilityrobustness evaluation when the exact values of system matricesare not available. All the above results could be directly extendedto similar problems on controllability robustness under actuatorfailures by duality between controllability and observability.The rest of this paper is organized as follows. Section IIprovides the problem formulation, and Section III presents somepreliminaries. Sections IV and V give the computational com-plexity and algorithms for the considered problem, respectively.Section VI studies the structure counterpart of the minSRO, withSection VII providing some numerical results. The last sectionconcludes this paper. Notations:
Symbols R , C and N denote the sets of real,complex and integer numbers, respectively. For an m × n matrix M , letting J ⊆ { , ..., m } and J ⊆ { , ..., n } , M J ,J denotesthe submatrix of M whose rows are indexed by J and columnsare indexed by J , M J the submatrix consisting of rows indexedby J , and M [ J ] the submatrix consisting of columns indexedby J . By col { X i | Ni =1 } (resp. diag { X i | Ni =1 } ) we denote thecomposite (resp. diagonal) matrix with its i th row block (resp.diagonal block) being X i , and span( X ) the space spanned by row vectors of X .II. P ROBLEM F ORMULATION
Consider the following LTI system ˙ x ( t ) = Ax ( t ) + Bu ( t ) , y ( t ) = Cx ( t ) (2)where x ( t ) ∈ R n , u ( t ) ∈ R m and y ( t ) ∈ R r are respectively thestate, input and output vectors, A ∈ R n × n , B = [ b , ..., b m ] ∈ R n × m and C = col { c i | ri =1 } ∈ R r × n are respectively the statetransition, input and output matrices, with b i ∈ R n , c i ∈ R × n .Without loss of generality, assume that every row of C isnonzero . Let V . = { , ..., r } .In the field of cyber-physical security, it is of both theoreticaland practical interest to know whether a given system like (2)can preserve its observability under the removal of a cardinality-constrained set of sensors [2–5]. As mentioned earlier, thisis not only related to the ability of unique reconstruction ofsystem states under adversarial sensor attacks, but also significantto measure observability robustness/resilience under denial-of-service attacks on sensors. We say System (2) is s -robustobservable, s ∈ N , if under the removal of arbitrary s sensors, theresulting system is still observable. To measure such observabil-ity robustness/resilience, the following minSRO is considered. Problem Given ( A, C ) in (2) min S ⊆ V | S | ( A, C V \ S ) is unobservable Denote the cardinality of the optimal solution to Problem 1by r min . Then, it is easy to see that, System (2) is s -robustobservable for any integer s ≤ r min − . Moreover, accordingto [2, 6], it is possible to uniquely reconstruct the system statesunder s ′ sensor attacks with s ′ ≤ (cid:4) r min − (cid:5) (see Section I), where ⌊•⌋ denotes the floor of a scalar.Next, we will also consider the structure counterpart ofProblem 1. This is motivated by the fact that, the exact valuesof system matrices may sometimes be hard to know in practisedue to modeling errors or parameter uncertainties. Instead, thesparsity patterns of system matrices are often easily accessible.Under this circumstance, one may be interested in the genericproperties of a system, i.e., properties that hold almost every-where in the corresponding parameter space [17]. Controllabilityand observability are such generic properties.To be specific, a structured matrix is a matrix whose entriesare either fixed zero or free parameters. Denote the set of all n × m structured matrices by { , ∗} n × m , where ∗ denotes thefree parameters. A structured system is a system whose systemmatrices are represented by structured matrices. Such system issaid to be structurally observable , if there exists one numericalrealization of it (i.e., a pair of numerical matrix that has the samesparsity pattern as the system matrices) that is observable. Let ¯ A and ¯ C be two structured matrices consistent with the sparsitypatterns of A and C . We consider the following problem. Problem Given ( ¯ A, ¯ C ) of ( A, C ) in (2) min S ⊆ V | S | ( ¯ A, ¯ C V \ S ) is structurally unobservable By definition of structural observability, the optimal valueof Problem 2 is an upper bound of that of Problem 1. Asobservability is a generic property, it is expected that, the optimalvalues of these two problems tend to coincide with each otherin randomly generated systems, which will also be demonstratedby the numerical experiments in Section VII.III. P
RELIMINARIES
This section presents some necessary preliminaries on observ-ability, structural observability and graph theories.
Lemma [PBH test,[18]] System (2) is observable, if andonly if for each λ ∈ C , there exists no x ∈ C n , x = 0 , such that Ax = λx and Cx = 0 .Given System (2), suppose A has p ≤ n distinct eigenvalues,and denote the i th one by λ i . Let k i be the dimension of thenull space of λ i I − A , i.e., k i is the geometric multiplicity of λ i . Let X i be an eigenbasis of A associated with λ i , that is, X i consists of k i vectors which are linearly independent spanningthe null space of λ i I − A . With these definitions, the followinglemma is immediate from the PBH test. Lemma
System (2) is observable, if and only if CX i is of full column rank for i = 1 , ..., p .Criteria for structural observability often have explicit graph-ical presentations. To this end, for System (2), let X = { x , ..., x n } , Y = { y , ..., y m } denote the sets of state ver-tices and output vertices respectively. Denote the edges by E X , X = { ( x i , x j ) : ¯ A ji = 0 } , E X , Y = { ( x j , y i ) : ¯ C ij = 0 } .Let D ( ¯ A, ¯ C ) = ( X ∪ Y , E X , X ∪ E X , Y ) be the system digraphassociated with ( ¯ A, ¯ C ) , and D ( ¯ A ) = ( X , E X , X ) . A state vertex x ∈ X is said to be output-reachable, if there exists a path from x to an output vertex y ∈ Y in D ( ¯ A, ¯ C ) .The generic rank of a structured matrix M , denoted by grank( M ) , is the maximum rank M can achieve as the functionof its free parameters. The bipartite graph associated with amatrix M is given by B ( M ) = ( R ( M ) , C ( M ) , E ( M )) , wherethe left vertex set R ( M ) (resp. right vertex set C ( M ) ) cor-responds to the row (resp. column) index set of M , and theedge set corresponds to the set of nonzero entries of M , i.e., E ( M ) = { ( i, j ) : i ∈ R ( M ) , j ∈ C ( M ) , M ij = 0 } . A matchingof a bipartite graph is a subset of its edges among which anytwo do not share a common vertex. A vertex is said to beunmatched associated with a matching, if it is not in edges ofthis matching. The maximum matching is the matching withthe largest number of edges among all possible matchings. It isknown that, grank( M ) equals the cardinality of the maximummatching of B ( M ) [19]. Lemma
The pair ( ¯ A, ¯ C ) is structurally observable, ifand only if (a) every state vertex is output-reachable in D ( ¯ A, ¯ C ) ,and (b) grank([ ¯ A ⊺ , ¯ C ⊺ ]) = n .We shall call (a) of Lemma 3 the output-reachability condition ,and (b) the matching condition .IV. C OMPUTATIONAL C OMPLEXITY OF P ROBLEM
1A natural question is whether Problem 1 is solvable inpolynomial time. In this section we will give a negative answer.
Theorem Problem 1 is NP-hard, even when C = I n .We will give a proof based on the linear degeneracy problem ,defined as follows. Definition
The linear degeneracy problem is to deter-mine whether a given d × p rational matrix F = [ f , ..., f p ] con-tains a degenerate submatrix of order d , i.e., det[ f i , ..., f i d ] = 0 ,for some { i , ..., i d } ⊆ { , ..., p } .In [21], it is proven that the linear degeneracy problem isNP-complete, and there are infinitely many integral matricesassociated with which this problem is NP-complete. Proof of Theorem 1:
We give a reduction from the lineardegeneracy problem . Let X be an arbitrary k × n integral matrix( k < n ) with full row rank, and α = max {| X ij |} . Let X ⊥ bethe basis matrix of the null space of X , i.e., X ⊥ is an n × ( n − k ) matrix spanning the null space of X . X ⊥ can be constructed viathe Gaussian elimination method in polynomial time O ( n ) [22].Moreover, let α = max {| X ⊥ ij |} , and α max = max { α , α } .Notice that the entries of X ⊥ are rational, and X ⊥ multipliedby any nonzero scalar is still a basis matrix of the null spaceof X . Hence, the encoding length of α (i.e., log α ) can bepolynomially bounded by k and n ; so is α max . Define H ( η ) as H ( η ) = (cid:2) X ⊺ , X ⊥ + η n × ( n − k ) (cid:3) , where n × ( n − k ) denotes the n × ( n − k ) matrix whose en-tries are all one. Then, clearly, det H (0) = 0 . We are toshow that one can find a scalar η ∗ satisfying η ∗ > α and det H ( η ∗ ) = 0 in polynomial time. Notice that det H ( η ) isa polynomial of η with degree at most one as the coefficientmatrix of η in H ( η ) is rank-one. Hence, H ( η ) can be expressedas det H ( η ) = [det H (1) − det H (0)] η + det H (0) . There-fore, in arbitrary set consisting of distinct rational numberswhich are bigger than α , there must exist one η ∗ , such that det H (0)(1 − η ∗ ) + det H (1) η ∗ = 0 , leading to det H ( η ∗ ) = 0 .Moreover, it holds that det H (0) ≤ α n max n n by noting that det H (0) consists of the summations of n ! ≤ n n signed productsof precisely one entry per row and column of H (0) . Similarly, det H (1) ≤ ( α max +1) n n n . Hence, both det H (0) and det H (1) have encoding lengths (which are n log α max + n log n and n log ( α max + 1) + n log n respectively) polynomially boundedby n .Afterwards, let matrices P = H ( η ∗ ) , Γ = diag { I k , , , ..., n − k + 1 } , and construct the system as A = P Γ P − , C = I n . Since the encoding lengths of entries of P are polynomiallybounded by n , its inversion P − can be computed in polynomialtime and has polynomially bounded encoding lengths too. Hence, A can be computed in polynomial time.We claim that the optimal value of Problem 1 associated with ( A, C ) is no more than n − k , if and only if there exists an k × k submatrix of X which has zero determinant.Indeed, from the construction of A , the ( k + i ) th column of P , denoted by P [ k + i ] , is the eigenvector of A associated withthe eigenvalue i + 1 , i = 1 , ..., n − k . Notice that, all entries of P [ k + i ] are nonzero as η ∗ > α . Hence, all n rows of C need tobe removed, such that the resulting C V \ S P [ k + i ] fails to be of fullcolumn rank, where S ⊆ V . From Lemma 2, the optimal valueof Problem 1 associated with ( A, C ) then equals the minimalnumber of rows (resp. columns) whose removal from X ⊺ (, X )makes the resulting matrix fail to be of full column (row) rank.If such value is no more than n − k , there must exist an k × k + submatrix of X which is not of full row rank for some k + ≥ k .Then, the linear degeneracy problem associated with X is yes.Conversely, suppose that there is an k × k submatrix of X withzero determinant, and denote it by X [ ¯ S ] , ¯ S ⊆ V . Then, clearlyone just needs to remove the columns indexed by V \ ¯ S from X ,such that the resulting X [ ¯ S ] fails to be of full row rank. Hence,the optimal value of Problem 1 is no more than | V \ ¯ S | = n − k .Combining the fact that the linear degeneracy problem associatedwith X is NP-complete, we prove the NP-hardness of Problem 1. (cid:4) The above result indicates that, there is in general no efficientway to verify whether a given system is s -robust observable fora given s ∈ N unless P = N P . This confirms the conjectureimplicitly made in [2, 3, 6] where they tend to think thatcomputing the maximal number of sensor attacks that can betolerated is computationally intractable (see Section I).Notice that in the proof of Theorem 1, C = I n which meansthat each sensor measures exactly one state variable. Theorem 1thus indicates that, given ( A, C ) , it is NP-hard to determine theminimal number of state variables that need to be prevented frombeing directly measured (i.e., the minimal number of nonzerocolumns of C that need to be removed), such that the resultingsystem becomes unobservable.By duality between controllability and observability, it is easyto obtain the following corollary on controllability robustnessunder actuator failures/removals. Corollary For System (2), it is NP-hard to determine theminimal number of actuators whose removal makes the resultingsystem uncontrollable.
Taking the input and output into consideration together, wehave the following conclusion.
Corollary For System (2), it is NP-hard to determine theminimal total number of sensors and actuators whose removalmakes the resulting system neither controllable nor observable.
Proof:
Construct ( A, C ) as suggested in the proof of Theorem1, i.e., A = P Γ P − , C = I n . Construct B such that B ∈ R n × k and ( A, B ) is controllable. Since the eigenbases P of A areavailable, the matrix B can be constructed in polynomial timeas suggested in [23]. Notice that k is the maximum geometricmultiplicity of eigenvalues of A . According to [23], the minimalnumber of inputs (actuators) that ensures controllability of theassociated system equals k . Hence, removing any one of these k actuators can make the system uncontrollable. As a result, theminimal total number of actuators and sensors whose removalcauses uncontrollability and unobservability equals r min + 1 ,where r min is the optimum of the minSRO with ( A, C ) . Thelatter problem is shown to be NP-hard in Theorem 1. Therequired result follows. (cid:4) V. A
LGORITHMS
In this section, we first give a polynomial time algorithmfor Problem 1 on systems with bounded eigenvalue geometricmultiplicities. Then, we extend this algorithm to the case whereeach sensor has a non-negative cost. We shall assume that acollection of eigenbases of A are computationally available, anddenote them by X i | pi =1 . The symbols p , k i | pi =1 are defined inSection III. A. Simple Dynamics Case
At the beginning, we consider the simple spectrum case where A has no repeated eigenvalues. In this case, { X i | ni =1 } becomesa collection of eigenvectors of A . For the i th eigenvalue of A ,define F i as the collection of rows of C that are not orthogonal to X i . Since X i is a vector, c j X i becomes a scalar, j ∈ { , ..., r } .One has F i = { j : c j X i = 0 , j ∈ { , ..., r }} . Then, from Lemma2, it is obvious that r min = min i ∈{ ,...,n } | F i | . Obviously, finding r min can be done in polynomial time. Theabove analysis may seem trivial. In the following we will extendit to the general case where the geometric multiplicities ofeigenvalues of A can be greater than one. B. General Case
Now we assume that the geometric multiplicities of eigen-values of A are bounded by some fixed constant ¯ k ∈ N .That is, k i ≤ ¯ k , i = 1 , ..., p , as n and r increase. For mostpractical engineering systems, this assumption is reasonable. Infact, it is found in [24] that, with probability tending to ,the symmetric matrices with random entries have no repeatedeigenvalues, and Erdos-Renyi (ER) random graph has simplespectrum asymptotically almost surely.Let us focus on an individual eigenvalue λ i , i ∈ { , ..., p } .Let r i be the minimum number of rows whose removal from Y ( i ) . = col { c j X i | rj =1 } makes the remaining matrix fail to befull of column rank. To determine r i , an exhaustive combina-torial search needs to compute the ranks of at most (cid:18) r (cid:19) + (cid:18) r (cid:19) + · · · + (cid:18) rr − k i (cid:19) → O (2 r ) submatrices, which increasesexponentially with r even when k i is bounded. Hence, the directcombinatorial search is not computationally efficient. Algorithm 1 : Algorithm for Problem 1
Input:
System parameters ( A, C ) Output:
The optimal solution F min to Problem 1 Calculate the eigenbases { X i | pi =1 } of A . for κ = 1 to p do Initialize τ =0 , T (0)1 = ∅ , L − = 1 . while τ < k κ do for i = 1 to L τ − do ¯ T ( τ ) i = T ( τ ) i ∪ (cid:26) j : span( Y ( κ ) { j } ) ∈ span( Y ( κ ) T ( τ ) i ) , j ∈ V \ T ( τ ) i (cid:27) . Ω ( τ ) i = V \ ¯ T ( τ ) i . Let c ( τ ) i = P il =1 | Ω ( τ ) l | , c ( τ )0 =0 . Rewrite Ω ( τ ) i . = { j , ..., j | Ω ( τ ) i | } and let T ( τ +1) c ( τ ) i − + q = ¯ T ( τ ) i S { j q } for q = 1 , ..., | Ω ( τ ) i | . end for L τ = c ( τ ) L τ − , τ + 1 ← τ . end while T [ κ ]max = arg max n | ¯ T ( k κ − | , · · · , | ¯ T ( k κ − L ( kκ − | o , F κ = V \ T [ κ ]max . end for Return F min = arg min {| F | , ..., | F p |} and r min = | F min | .In what follows, a polynomial time algorithm based ontraversals over a recursive tree is presented. Here, ‘recursive’comes from the fact that this tree is constructed recursively .The pseudo code of this algorithm is given in Algorithm 1.The intuition behind Algorithm 1 is that, for the κ th eigenvalueof A , κ ∈ { , ..., p } , instead of directly determining r κ , wetry to determine r − r κ , which is the maximum number ofrows of Y ( κ ) that fail to have full column rank. Algorithm 1first builds a recursive tree and then searches the maximumreturn value among the leafs of this tree. An illustration ofsuch recursive tree is given in Fig. 1. To build the recursivetree for the κ th eigenvalue, the root is ¯ T (0)1 = ∅ . In the τ thlayer, τ = 0 , ..., k κ − , the set T ( τ ) i ⊆ V has the propertythat Y ( κ ) T ( τ ) i has rank τ , i = 1 , ..., L τ − , where L τ − is thenumber of nodes in the τ th layer. Then, Y ( κ )¯ T ( τ ) i is obtained byadding the maximal rows of Y ( κ ) V \ T ( τ ) i to Y ( κ ) T ( τ ) i while its rank ispreserved. The collection of sets { ¯ T ( τ ) i | L τ − i =1 } form the τ th layerof the tree. Next, each set ¯ T ( τ ) i generates | Ω ( τ ) i | sets, constituting { T ( τ +1) j | L τ j =1 } for the ( τ + 1) th layer, where Ω ( τ ) i is such setwhose arbitrary element j makes matrix Y ( κ )¯ T ( τ ) i ∪{ j } have rank τ + 1 . After τ = k κ − , the complete recursive tree is obtained,and one just needs to search the set with maximum cardinalityamong the leafs of this tree, i.e., T [ κ ]max , whose cardinality equalsexactly r − r κ . The key principle of Algorithm 1 is the rank-oneupdate rule expressed by (3). Theorem Given a system ( A, C ) with geometric multiplic-ities of eigenvalues of A bounded by ¯ k , A ∈ R n × n , C ∈ R r × n , In this paper, the root of a tree is indexed as the th layer. Fig. 1: Illustration of the recursive tree in Algorithm 1 associatedwith Y (1) . In this example, ¯ T (0)1 = ∅ , ¯ T (1) i = { i } for i = 1 , , , ¯ T (1)2 = ¯ T (1)4 = { , } , and ¯ T (2) j = { , , , } for j =1 , , , , , , , , , , ¯ T (2) l = { , , } for l = 7 , , , , ¯ T (2)4 = ¯ T (2)15 = { , } , ¯ T (2)11 = ¯ T (2)17 = { , } .Algorithm 1 can determine the optimal solution to Problem 1 intime complexity at most O (¯ k r ¯ k +1 n + rn ) . Proof:
As analyzed above, the main procedure of Algorithm 1is to determine the maximum number of rows of Y ( κ ) that failto have full column rank for each eigenvalue λ κ , κ = 1 , ..., p .To justify the recursive procedure, observe that for any S ⊆ V , j ∈ V , rank Y ( κ ) S S { j } = ( rank Y ( κ ) S , if span( Y ( κ ) { j } ) ⊆ span( Y ( κ ) S )rank Y ( κ ) S + 1 , if span( Y ( κ ) { j } ) / ∈ span( Y ( κ ) S ) (3)Hence, for the κ th eigenvalue, the recursive tree has depthexactly k κ . Moreover, { ¯ T ( k κ − , · · · , ¯ T ( k κ − L kκ − } index all thesubmatrices formed by rows of Y ( κ ) with the property that: ithas rank k κ − , and adding any rest rows from Y ( κ ) can makeit have full column rank. The optimality of the solution returnedfrom Algorithm 1 then follows immediately.In the recursive tree for the κ th eigenvalue, each parent nodehas at most r child nodes. Thus, the τ th layer has at most r τ nodes, τ = 0 , ..., k κ − . Hence, the total nodes in that treeis at most P k κ − τ =0 r τ = ( r k κ − / ( r − → O ( r k κ − ) In practical systems, different sensors may incur differentcosts to be removed. Here, cost could be used to measure thebudget of a sensor to be destroyed, or the difficulty/fragility tobe attacked. In such case, an attacker/protector may be moreinterested in the minimal cost set of sensors whose removalcauses unobservability. Let c ( i ) ≥ denote the cost of the i thsensor, i ∈ V . Given S ⊆ V , let c ( S ) = P i ∈ S c ( i ) . Then,Problem 1 can be generalized to the following Problem 3 Problem Given ( A, C ) in (2) min S ⊆ V c ( S )( A, C V \ S ) is unobservable The following theorem reveals that Problem 3 is polynomialtime solvable for systems with bounded maximum geometricmultiplicities of A . Theorem Given a system ( A, C ) with geometric multiplic-ities of eigenvalues of A bounded by ¯ k , A ∈ R n × n , C ∈ R r × n ,Algorithm 2 can determine the optimal solution to Problem 3with complexity at most O (¯ k r ¯ k +1 n + rn ) . Proof: Note that the distinction between Algorithm 2 andAlgorithm 1 lies in Line 4 of Algorithm 2. To be specific, inAlgorithm 2, the cardinality of a set of sensors is replaced by itstotal cost. With this observation, the proof of Theorem 3 followssimilar arguments to that of Theorem 2. (cid:4) Algorithm 2 : Algorithm for Problem 3 Input: System parameters ( A, C ) , sensor costs { c ( i ) } i =1 ,...,r Output: The optimal solution F c min to Problem 3 Calculate the eigenbases { X i | pi =1 } of A . for κ = 1 to p do Run the same procedure in Lines 3-11 of Algorithm 1 . T [ κ ]max = arg max { ¯ T ( kκ − i } n c ( ¯ T ( k κ − ) , · · · , c ( ¯ T ( k κ − L ( kκ − ) o , F κ = V \ T [ κ ]max . end for Return F c min = arg min { F i } {| c ( F ) | , ..., | c ( F p ) |} and r c min = c ( F c min ) .VI. C OMPLEXITY AND A LGORITHMS FOR P ROBLEM A. Computational Complexity Theorem Problem 2 is NP-hard.We first give some definitions and lemmas needed to proceedwith our proof. Let ¯1 m × n be an m × n structured matrix whoseentries are all free parameters. A set of columns of a structuredmatrix are said to be linearly dependent, if the generic rank of the matrix consisting of them is less than the number of theircolumns. A circuit of an m × n matrix M , is a set J ⊆ { , ..., n } such that rank( M [ J ] ) = | J | − . The set of indices of anylinearly dependent columns of M contains at least one circuitof M . A clique of an undirected graph is a subgraph such thatany two of its vertices are adjacent. The k -clique problem is todetermine whether there is a clique with k vertices, which isNP-complete [22]. Lemma Let matrix ¯ M ∈ { , ∗} n × l . If grank([ ¯1 n × m , ¯ M ]) < n , where m < n , then grank( ¯ M ) < n − m . Proof: Suppose grank( ¯ M ) = p ≥ n − m . Then, there exist twosets J ⊆ { , ..., n } and J ⊆ { , ..., l } with | J | = | J | = p ,such that grank( ¯ M J ,J ) ≥ n − m . Then, grank([ ¯1 n × m , ¯ M ]) ≥ grank( (cid:20) ¯1 n − p,m ¯1 p,m ¯ M J ,J (cid:21) ) = min { n − p, m } + p ≥ n , whichcauses a contradiction. (cid:4) Proof of Theorem 4: Our proof is a reduction from the NP-complete k -clique problem to one instance of Problem 2. Wedivide it into three steps.Step 1). Let G = ( V , E ) be an undirected graph, and |V| = p , |E| = r . The incidence matrix of G , denoted by In ( G ) , isan r × p matrix so that [ In ( G )] ij = 1 if v j ∈ V ( e i ) and otherwise, where v j ∈ V , e i ∈ E . Construct the output matrix ¯ C as ¯ C = (cid:20) In ( G ) , ¯1 r × ( (cid:16) k (cid:17) − k − (cid:21) , where k ≥ (so that (cid:0) k (cid:1) − k − ≥ ), and ‘ ’ in In ( G ) represents a free entry.Let q . = (cid:0) k (cid:1) . Without any losing of generality, assume that p > k + 1 , which does not affect the NP-completeness of the k -clique problem associated with G . It demonstrates in Page 53of [25] that: 1) Every circuit J of ¯ C ⊺ must satisfy | J | ≥ q . 2) If G has an k -clique, then ¯ C ⊺ has a circuit with cardinality q . 3)If ¯ C ⊺ has a circuit with cardinality q , then G has an k -clique.From the aforementioned results, it can be declaimed that, ¯ C ⊺ has q linearly dependent columns, if and only if the digraph G has an k -clique. Indeed, from Fact 1), if ¯ C ⊺ has q linearlydependent columns, then such columns must form a circuit of ¯ C ⊺ . Otherwise, ¯ C ⊺ has a circuit with cardinality less than q .With Fact 3), the ‘only if’ direction is obtained. The ‘if’ directioncomes directly from Fact 2).Step 2). Construct the state transition matrix ¯ A as ¯ A = (cid:20) q × q q × ( n − q ) ¯1 ( n − q ) × q ¯1 ( n − q ) × ( n − q ) (cid:21) , where n . = p + q − k − > q . Construct the systemdigraph D ( ¯ A, ¯ C ) , and let state vertices X = X ∪ X with X = { x , ..., x q } , X = { x q +1 , ..., x n } , output vertices Y = { y , ..., y r } . From the construction, it is clear that every statevertex in X \{ x n } is an in-neighbor of x n , and x n is reachableto every output in Y . Hence, all the r outputs must be deletedso that at least one state vertex is output-unreachable. In otherwords, the minimal number of outputs (sensors) needed to bedeleted to destroy the output-reachability of at least one statevertex equals r (Lemma 3).Step 3). We declaim that the optimal value of Problem 2associated with ( ¯ A, ¯ C ) is no more than r − q , if and only ifthere exist q linearly dependent columns in ¯ C ⊺ . Let ¯ B . = ¯ C ⊺ .Denote the optimal solution to Problem 2 by S opt ⊆ V .For the one direction, suppose there are q linearly dependent columns in ¯ C ⊺ , and denote the set of its indices by S ld . Let S re = V \ S ld . It holds that grank([ ¯ A ⊺ , ¯ B [ V \ S re ] ]) = grank([ ¯ A ⊺ , ¯ B [ S ld ] ]) ≤ grank( ¯ A ) + grank( ¯ B [ S ld ] ) < n − q + q = n. As a result, S re is a feasible solution to Problem 2 associatedwith ( ¯ A, ¯ C ) . Since | S re | = r − q , we have | S opt | ≤ r − q .For the other direction, suppose that | S opt | ≤ r − q . Then,since the minimal number of sensors whose removal destroysthe output-reachability condition equals r , it must hold that grank([ ¯ A ⊺ , ¯ B [ V \ S opt ] ]) < n . As grank( ¯ A ) = n − q , fromLemma 4, it is valid that grank( ¯ B [ V \ S opt ] ) < q . Noting that | V \ S opt | ≥ q , we have that the columns of ¯ C ⊺ indexed by V \ S opt are linearly dependent.Since we have shown that verifying whether there are q linearly dependent columns in ¯ C ⊺ is equivalent to verifyingwhether the digraph G has an k -clique, which is NP-complete,and the above-mentioned reduction is in polynomial time, it isconcluded that Problem 2 is NP-hard. (cid:4) The key challenge to prove Theorem 4 lies in constructing asystem associate with which there is an explicit relationship insize between the number of sensors whose removal destroys theoutput-reachability condition, and that of sensors whose removaldestroys the matching condition. While [14] has proven theNP-hardness of determining the minimal cost of inputs whoseremoval causes structural uncontrollability, it cannot be obtainedby duality that the minSRO for structured systems is NP-hard,as in their proof different inputs may have different costs. B. Algorithms We will call the number of vertices of R ( ¯ A ) thatare not matched in any maximum matching of B ( ¯ A ) =( R ( ¯ A ) , C ( ¯ A ) , E ( ¯ A )) , as the matching deficiency , which is equalto n − grank( ¯ A ) . Assumption The matching deficiency of ¯ A is bounded bysome fixed constant ¯ k s as n and r increase. Remark Assumption 1 is satisfied by many practical sys-tems. For example, in consensus networks or multi-agent systems[26], every state vertex often has a self-loop, under whichcircumstance B ( ¯ A ) has a zero matching deficiency. Moreover,connected ER random networks can be controlled using onlyone dedicated input with probability [27], which means thatthe associated matching deficiencies are no more than [16]with probability .To determine the minimal number of sensors whose removaldestroys the output-reachability condition, we resort to thestrongly connected component (SCC) decomposition. A digraphis said to be strongly connected, if there is a path from everyone of its vertices to the other one. An SCC of G is a subgraphthat is strongly connected and no additional edges or verticesfrom G can be included in this subgraph without breaking itsproperty of being strongly connected. We call an SCC the sinkSCC, if there is no outgoing edge from this SCC to any otherSCCs. Suppose there are l sink SCCs in D ( ¯ A ) = ( X , E X , X ) . Foreach sink SCC N i , denote by N ( N i ) the out-neighbors of N i in D ( ¯ A, ¯ C ) (hence N ( N i ) ⊆ Y ). Define J yi = { j : y j ∈ N ( N i ) } ;that is, J yi is the set of sensors that are out-neighbors of N i .Then, it is easy to verify that, the minimal number of sensors whose removal destroys the output-reachability condition equals min ≤ i ≤ l | J yi | .To determine the minimal sensors whose removal destroysthe matching condition (denoted by J yma ), the similar ideato Algorithm 1, by constructing a recursive tree, is utilized.This tree has a similar structure to Fig. 1; see the pseudocode in Algorithm 3. The root of this tree is an empty set.In the τ th layer, each node ¯ T ( τ ) i indexes a subset of sensorssuch that grank( col { ¯ A, ¯ C ¯ T τi } ) = grank( ¯ A ) + τ , and no otherelement from V \ ¯ T ( τ ) i can be added to ¯ T τi without increasing grank( col { ¯ A, ¯ C ¯ T τi } ) . Hence, in the τ m . = ( n − grank( ¯ A ) − thlayer, the leaf with the maximum cardinality is the comple-mentary set to the minimal rows of ¯ C whose removal from col { ¯ A, ¯ C } makes the resulting matrix have generic rank n − .This means the total layers of the recursive tree is exactly n − grank( ¯ A ) . Theorem Given a system ( ¯ A, ¯ C ) , where ¯ A ∈ { , ∗} n × n , ¯ C ∈ { , ∗} r × n , under Assumption 1, Algorithm 3 provides anoptimal solution to Problem 2 with time complexity at most O (( r + n ) . r ¯ k s ) . Proof : The proof for optimality of Algorithm 3 follows theabove arguments, and is similar to that of Theorem 2. The SCCdecomposition along with the step to obtain J yre has complexityat most O ( n + nr ) [22]. In building the recursive tree, the depthis exactly n − grank( ¯ A ) , which is no more than ¯ k s . In the i thlayer, there is at most r i nodes, i = 0 , ..., n − grank( ¯ A ) − .Hence, the number of total nodes in the recursive tree is at most P ¯ k s − τ =0 r τ = ( r ¯ k s − / ( r − → O ( r ¯ k s − ) ( r ≥ ). In eachnode, the operation updating ¯ T ( τ ) i can be done by calling themaximum matching algorithm on bipartite graphs for at most r times, which incurs time complexity O (( r + n ) . r ) in theworst case [22]. Hence, the total time complexity of Algorithm3 is O ( n + rn + ( r + n ) . r ¯ k s ) , i.e, O (( r + n ) . r ¯ k s ) . (cid:4) C. Acceleration of Algorithm 3 The main procedure dominating the complexity of Algo-rithm 3 is to obtain the minimal sensors J yma whose removaldestroys the matching condition. A method to accelerate this pro-cedure is first to do Dulmage-Mendelsohn decomposition (DM-decomposition) on ¯ A . Then implement Lines 2-11 of Algorithm3 on the reduced system, which usually has a dimension muchless than that of the original system. Definition Let ¯ M ∈ { , ∗} m × n .There exist two permutation matrices P r ∈ R m × m and P c ∈ R n × n , such that P r ¯ MP c = ¯ M ¯ M ¯ M ¯ M M ¯ M M M , where ¯ M , ¯ M , and ¯ M are square with zero-free diagonals.The columns of ¯ M (or the rows of ¯ M ) correspond to theunmatched vertices of at least one maximum matching in B ( ¯ M ) ,and the rows (or columns) of ¯ M correspond to the vertices thatappear in every maximum matching in B ( ¯ M ) . The generic rankof ¯ M equals m minus the number of rows of ¯ M .Using the DM-decomposition, suppose there are two permu- tation matrices P r and P c with compatible dimensions, such that P r [ ¯ A ⊺ , ¯ C ⊺ ] (cid:20) P c I (cid:21) = ¯ A ¯ A ¯ A ¯ A ¯ B A ¯ A ¯ B A ¯ B A ¯ B , where the left block of the right-hand side matrix of the above-mentioned formula is the DM-decomposition of ¯ A ⊺ . Since (cid:20) ¯ A 12 ¯ A 130 ¯ A (cid:21) has full generic rank, it can be validated that grank([ ¯ A ⊺ , ¯ C ⊺ ]) = (cid:12)(cid:12)(cid:12)(cid:12) C ( (cid:20) ¯ A ¯ A A (cid:21) ) (cid:12)(cid:12)(cid:12)(cid:12) +grank( (cid:20) ¯ A ¯ B ¯ A ¯ B (cid:21) ) where C ( • ) denotes the set of column indices of a matrix. Hence,to determine J yma , we just need to identify the minimal numberof columns of col { ¯ B , ¯ B } whose deletion makes (cid:2) ¯ A 34 ¯ B A 44 ¯ B (cid:3) row generic rank deficient. As col { ¯ B , ¯ B } is obtained by per-mutating columns of ¯ C , the set of indices of the aforementionedcolumns exactly correspond to J yma . Hence, the reduced system (cid:18)(cid:20) ¯ A ⊺ ¯ A ⊺ (cid:21) , [ ¯ B ⊺ , ¯ B ⊺ ] (cid:19) has the same solution as that of the original system ( ¯ A, ¯ C ) on J yma , where zero matrix is added to [ ¯ A ⊺ , ¯ A ⊺ ] to square it.For a sparse matrix ¯ A , |C ( col { ¯ A , ¯ A } ) | is usually much lessthan |C ( ¯ A ) | . Hence, it is perhaps safe to expect that the DM-decomposition can accelerate Algorithm 3. Algorithm 3 : Algorithm for Problem 2 Input: System parameters ( ¯ A, ¯ C ) Output: The optimal solution J opt to Problem 2 Find all the sink SCCs {N i | li =1 } of D ( ¯ A ) . Let J yi = { j : y j ∈ N ( N i ) } . J yre = arg min ≤ i ≤ l | J yi | . Initialize τ =0 , T (0)1 = φ , L − = 1 . while τ < n − grank( ¯ A ) do for i = 1 to L τ − do ¯ T ( τ ) i = T ( τ ) i ∪ { j : grank( col { ¯ A, ¯ C T ( τ ) i ∪{ j } } )= grank( col { ¯ A, ¯ C T ( τ ) i } ) , j ∈ V \ T ( τ ) i } . Ω ( τ ) i = V \ ¯ T ( τ ) i . Let c ( τ ) i = P il =1 | Ω ( τ ) l | , c ( τ )0 =0 . Rewrite Ω ( τ ) i . = { j , ..., j | Ω ( τ ) i | } and let T ( τ +1) c ( τ ) i − + q = ¯ T ( τ ) i S j q for q = 1 , ..., | Ω ( τ ) i | . end for L τ = c ( τ ) L τ − , τ + 1 ← τ . end while T max = arg max n | ¯ T ( τ m − | , · · · , | ¯ T ( τ m − L τm − | o , J yma = V \ T max , with τ m . = n − grank( ¯ A ) − . Return J opt = arg min {| J yma | , | J yre |} , with the optimal value | J opt | = min {| J yma | , | J yre |} .VII. N UMERICAL R ESULTS We implement some numerical experiments to validate theeffectiveness of our proposed algorithms. In our simulations,we generate ER random networks with n nodes and use theiradjacent matrices as the state transition matrices, n rangingfrom to . For each n , the probability of the existenceof a directed edge between two nodes is set to be cn , which 20 40 60 80 100 12000.10.20.30.40.50.60.70.80.91 Number of Nodes (n) R a t i o o f R e m o v ed S en s o r s Degree basedRandomAlgorithm 1Algorithm 3Algorithm 3+DM (a) c = 3 20 40 60 80 100 12000.10.20.30.40.50.60.70.80.91 Number of Nodes (n) R a t i o o f R e m o v ed S en s o r s Degree basedRandomAlgorithm 1Algorithm 3Algorithm 3+DM (b) c = 4 Fig. 2: Experiments on ER random networks. Note that thevertical axis is the ratio between the number of removed sensorsand that of the total sensors.means that the average degree (the sum of in-degree and out-degree) is c [22]. The weight of each directed edge is randomlygenerated in the internal [ − , . For each n , randomly generate observable systems (to save experiment time, all these systemshave maximum geometric multiplicities less than ), where of the nodes are randomly selected to be measured bydedicated sensors (i.e., each sensor measures exactly one node).Five methods are adopted to identify a set of sensors whoseremoval destroys (structural) observability, the first of which isselecting the measured nodes in the descending order of their in-degrees, the second of which is randomly selecting the measurednodes, and the rest three of which are respectively based onAlgorithm 1, Algorithm 3 and Algorithm 3 accelerated by theDM-decomposition. The results are shown in Fig. 2.Fig. 2 shows that, for each fixed n , the ratio of removedsensors returned by Algorithm 1 is much less than that ofthe degree based and the random selection methods. This isreasonable, as Algorithm 1 returns the optimal solutions. More-over, Algorithms 1, 3 and the DM-decomposition acceleratedAlgorithm 3 return almost the same values for each fixed n . Thisindicates that, when we are not accessible to the exact values ofsystem matrices, the structured system based Algorithm 3 couldbe an alternative for the observability robustness evaluation.Comparing Fig. 2(a) and (b), it seems that networks with strongerconnectivity (i.e., bigger average degrees) tend to have betterobservability robustness under sensor attacks.VIII. C ONCLUSIONS The problem of determining the minimal number of sensorswhose removal destroys system observability is considered. Thisproblem is closely related to secure state estimation underadversarial sensor attacks, and its dual is also the inverse tothe MCP [8, 9]. It is shown that this problem is NP-hard,both for a numerically specific system, and for a structuredsystem, thus confirming the conjecture that verifying the s -sparse observability condition in [2, 3] is NP-hard. Neverthe-less, polynomial time algorithms are provided to solve thisproblem, respectively, on a numerical system with boundedmaximum geometric multiplicities, and on a structured systemwith bounded matching deficiencies, which are often met by A deep insight shows that for some n , Algorithm 1 returns values slightlysmaller than those of Algorithm 3. This is consistent with the fact that for agiven ( A, C ) , the optimal value of Problem 2 is a theoretical upper bound ofthat of Problem 1. practical engineering systems. Numerical experiments show thatwhen we have no access to the exact values of system matrices,the structured system based algorithms could be alternative forobservability robustness evaluation.Rpractical engineering systems. Numerical experiments show thatwhen we have no access to the exact values of system matrices,the structured system based algorithms could be alternative forobservability robustness evaluation.R