Network Topology-aware Mitigation of Undetectable PMU Time Synchronization Attacks

—Time Synchronization attacks constitute a major threat to PMU-based smart grid applications, their cost-efficient detection and mitigation is thus of utmost importance. In this paper we propose a mitigation approach based on authenticated network-based time synchronization. Our approach relies on the observation that a time synchronization attack can be undetectable if and only if it targets at least three time references in the power system, and such attacks need to be mitigated through appropriate security controls. We first provide a formal proof of this result, including a characterization of the degrees of freedom of the attacker in constructing an attack. We then formulate the problem of mitigating undetectable attacks at minimum cost as an integer linear program, and prove that it is NP-hard. To solve the problem, we propose two approximation algorithms based on (1) computing shortest paths, and (2) solving a linear relaxation of the problem. Extensive simulations suggest the superiority of the proposed algorithms on IEEE benchmark power system graphs compared to baseline solutions. We report mitigation cost savings of at least 76% compared to a naive approach for mitigation and at least 30% compared the state-of-the-art.


I. INTRODUCTION
Wide Area Monitoring Protection And Control (WAMPAC) systems rely primarily on Phasor Measurement Units (PMUs).The high-frequency, timestamped measurements of voltage and current phasors generated by PMUs enable a variety of applications, such as phase angle monitoring [1], oscillation monitoring [2], and fault localization [3].Synchrophasor measurements require, however, precise time synchronization between PMUs, and thus secure and reliable time synchronization is of utmost importance to modern power system operators.
PMUs have traditionally relied on space-based time synchronization, provided by multiple satellites (e.g., GPS) equipped with atomic clocks.Recent years have seen the emergence of Network-Based Time Synchronization (NBTS) as an alternative to space based solutions [4], enabling operators to reduce the dependence on external systems.In the case of NBTS, each PMU clock acts as a slave clock and receives time information through a packet switched network from a master clock that is connected to an accurate time reference.The most prominent protocol for NBTS is the Precision Time Protocol (PTP) [5].
Unfortunately both space-based and network-based time synchronization are vulnerable to cyber-attacks, referred to as Time Synchronization Attacks (TSAs) [6] [7].Unauthenticated civilian GPS signals are vulnerable to spoofing attacks, and despite optional security features in PTPv2.1 [8], PTP messages are not authenticated by default and are thus also vulnerable to spoofing attacks.A TSA would induce a phase angle shift in synchrophasor measurements and can have a serious impact on WAMPAC applications [6].Intuitively, one would expect that Linear State Estimation (LSE) combined with commonly used Bad Data Detection (BDD) algorithms could serve for the detection of TSAs [9].Nonetheless, as shown recently, a careful attacker can compute TSAs that bypass state-of-the-art BDDs [10].The detection and mitigation of TSAs thus require additional security controls.
In this paper we formulate the problem of cost-efficient mitigation of TSAs by combining PTPv2.1 message authentication [8] and BDD.Our problem formulation is based on the observation that undetectable TSAs [10] are feasible against only a subset of PMUs, and thus it is possible to minimize the system level cost of deploying PTP authentication (e.g., in terms of network equipment to upgrade and associated key management) and the number of devices containing keying material by protecting an appropriately chosen subset of PMUs.Even though PTPv2.1 message authentication cannot eliminate all types of TSAs, e.g., delay attacks through inserting a delay-box, it does mitigate cyber attacks that involve message spoofing or injection.Therefore, in the rest of the paper, TSA mitigation refers to the mitigation of these classes of TSAs.
The main contributions of the paper are threefold.First, we prove a necessary and sufficient condition for the existence of undetectable TSAs and for the attackers degrees of freedom in performing an attack.Second, we use these findings to formulate the problem of mitigating undetectable TSAs at minimum cost as an Integer Linear Program (ILP).Third, we propose approximation algorithms to solve the problem and show through extensive simulations on IEEE benchmark systems that the proposed cyber-physical mitigation allows to make TSAs detectable at up to 95% cost reduction.
The rest of the paper is organized as follows.Section II accounts for previous research on the detection and mitigation of TSAs.In Section III we present a model for PTP time synchronization in power systems and introduce the theory of undetectable TSAs.Necessary and sufficient conditions for the existence of undetectable TSAs are shown in Section IV.Section V formulates the problem of mitigating undetectable TSAs at minimum cost and presents the proposed approximation algorithms.Section VI evaluates the performance of the proposed algorithms through extensive simulations.Finally, Section VII concludes the paper.

II. RELATED WORK
Undetectable TSAs are akin to false data injection attacks (FDIA) [11; 12], but typical FDIA mitigation approaches that protect the integrity of the measured data are not efficient against TSAs [12; 13; 14].Solutions have also been proposed for detecting FDIAs based on historical data and meta-data, e.g., authors in [15] analyzed correlations in consecutive estimated system states using wavelet transform and deep neural networks.Moreover, [16] proposed using a Kalman filter for the detection of FDIAs.More recently, graph signal processing has been utilized to detect FDIAs and identify the attacked PMUs [17; 18].A recent survey for FDIA detection methods is presented in [19].
A number of recent works proposed anomaly detectors specifically designed for detecting TSAs, mostly relying on k k th source vertex and k th terminal vertex information from the time synchronization system (GPS or PTP).For detecting GPS spoofing against PMUs, [20] proposed utilizing the carrier-to-noise-ratio of the GPS signal.For PTP, [21] proposed introducing guard clocks in order to detect PTP delay attacks.More recent work proposed combining data from both the power system and the time synchronization system to increase the detection accuracy [22].Nevertheless, anomaly detectors are prone to false negatives, and thus they may not detect skillful adversaries.
Works focusing on the mitigation of TSAs proposed authentication of GPS messages [23] or PTP messages [24].One disadvantage of such works is that they do not consider the system level cost of TSA mitigation, (e.g., cost of upgrading network switches in a PTP network, overhead due to key management, and compatibility issues) and are thus cost-inefficient for power systems.Similar to our work, [25] aims to reduce the TSA mitigation cost by deploying PTP authentication only on the subset of the network equipment that is vulnerable to undetectable TSAs.While [25] motivated the approach by a conjecture that only TSAs against 3 time references or more can be undetectable [26], this paper provides a proof of this conjecture in Section IV, and provides a different problem formulation in Section V, which allows to achieve significant cost savings compared to that in [25], by considering a wider range of TSA mitigation strategies.

III. SYSTEM MODEL A. Power System and Network Model
Table I summarizes the most important notations used in the paper.We consider a power system with N buses and M ≥ N measurements taken by PMUs to ensure system observability.
We denote the vector of PMU measurements by z ∈ C M , z m = |z m |e jθm , where j = √ −1 and θ m is the phase angle.We further denote the system state (voltage phasors at the buses) by x∈C N .The linear measurement model is then where H ∈ C M×N denotes the measurement matrix, and e ∈ C M is white Gaussian measurement noise.The linear model ( 1) is justified by that PMUs measure current and voltage phasors, and it allows for efficient computation of the Least-Squares (LS) state estimate x = (H † H) −1 H † z, where H † denotes the conjugate transpose of H. Alternative approaches using non-linear least squares based on power flow measurements or using dynamic state estimation (DSE) are more computationally intensive as they involve iterative algorithms [27].The verification matrix is defined as , and allows us to express the measurement residual r =F z ∈C M .The residual r is typically used for BDD, e.g., using the LNR test [9] for identifying measurements with suspiciously high residuals.Bad measurement data identified by the BDD are typically discarded, as long as system observability is maintained.The measurements are taken by a set T ={τ 1 ,...,τ T } of PMUs, T ≤M, which depend on precise time synchronization.We denote by M i ≥ 1 the number of measurements taken by PMU τ i .We consider securing the time references of a subset of these PMUs using PTP, and thus we model the communication infrastructure used for time synchronization in the WAMPAC by an undirected graph G = (V,E), where V = T ∪ N is the set of vertices, N is the set of network switches and routers, and E is the set of edges (communication links) in the graph.We consider that the edges E of G form a tree that spans the PMUs T .This assumption is motivated by that the active communication topology in power systems typically is a tree, and PTP networks use a tree topology as well for disseminating time references in a network [5].Furthermore, we denote by v r ∈V the root vertex of the tree, which corresponds to where the PTP master clock is deployed.

B. Attack Model
We consider an attacker that is able to spoof PTP messages [7] or the GPS signals [6].This way, the attacker can manipulate the time references of a subset T a ={τ a 1 ,...,τ a P }⊆T of PMUs, where P is the number of manipulated time references (and PMUs).To capture the dependence of measurements on the attacked time references we define the attack-measurement matrix Ψ∈{0,1} M×P as for all m∈{1,...,M} and p∈{1,...,P } [10].Due to the attack, the phasors measured by PMU τ a p will be rotated by α p .Thus, using the notation u p =e jαp , the m th measurement taken by PMU τ a p ,p∈ {1,...,P } can be written as z a p,m = |z p,m |e j(θp,m+αp) = z p,m u p , where |z p,m | is the measured magnitude, θ p,m is the unattacked phase angle, and α p is the phase angle shift introduced by the attack.

C. Undetectable TSAs
Let z a be the attacked measurement vector.If linear state estimation based on (1) is employed in the WAMPAC then a TSA should ideally yield high residuals and consequently it should be detected by BDD.It is thus natural to introduce the notion of undetectability as follows.
Definition 1.A TSA u = (u 1 , ... , u P ) against PMUs T a is undetectable if and only if it does not change the measurement residual, i.e., F z =F z a .
A necessary and sufficient condition for a TSA to be undetectable according to Defintion 1 was formulated in [10], as follows.
Lemma 1.Consider a TSA against PMUs T a .The TSA is undetectable if and only if the vector u ∈ C P , s.t., u p =e jαp ,p∈{1,...,P } satisfies where W =Ψ T diag(z) † F † F diag(z)Ψ is the complex attack angle matrix, W ∈C P ×P , and is Hermitian.
Observe that an undetectable TSA corresponds to any solution to (2) other than (u i = 1,∀i ∈ {1,...,P }), which corresponds to not attacking any time reference.Interestingly, when rank(W )=1, (2) can be simplified to where W 1i is the i th element of the row with highest 2 norm of the W matrix, and |u i | = 1,∀i ∈ {1,...,P }.The row with highest 2 norm in W is chosen in order to avoid choosing rows that are equal to the zero vector, and in order to provide better computational stability of the computed solutions.Moving the term related to u P from the LHS to the RHS of the above equation, the undetectability condition can be written as A geometric interpretation of the above in the complex plane, defined by the real and imaginary parts of the equation, is that the solution to (3) are the intersection points between an annular region (LHS of (3)), and a circle (RHS) [26].In what follows we denote by AR 1,P −1 = (c 1,P −1 ,r o 1,P −1 ,r i 1,P −1 ) the annular region with center c 1,P −1 = −  Lemma 2. Consider the case when rank(W )=1.
• If P =1 then there is no undetectable TSA.The equivalence classes are calculated by computing a metric called the minimum index of separation (IoS*), as proposed in [26], for all T 2 pairs of PMUs.When IoS*=1 for a pair of PMUs, this means that the rank of the W matrix corresponding to attacking this pair is always equal to 1, independent on the measurement values z.Vulnerable PMU pairs that have pairwise IoS*=1 were shown to form equivalence classes [26], i.e., IoS * (a,b)=1 and IoS * (b,c)=1 implies IoS * (a,c)=1.Consequently, solving (3) for a set of PMUs T a that is a subset of an equivalence class can result in undetectable TSAs according to Definition 1, while solving the same equation for a set of PMUs that is not a subset of an equivalence class does not result in undetectable TSAs.
TSAs against P = 2 time references were shown to be computable using a closed form expression in [10].For P ≥3, [26] proposed a recursive algorithm that selects u i ,i ∈ {3,...,P } (i.e., P −2 variables) from unions of connected sets of feasible values based on the geometric interpretation of (3), and thus conjectured that the solution to (3) has P − 2 degrees of freedom.The next section provides a formal definition of degrees of freedom and a proof of the conjecture, which is the first contribution of this paper.

IV. FEASIBILITY ANALAYSIS OF TSAS
In this section we investigate the feasibility of computing undetectable TSAs by solving (3).We first provide necessary and sufficient conditions under which an attack does not exist for P ≥2 and rank(W ) = 1.Next, we provide conditions for the set of attacks not to be a singleton and we characterize the degrees of freedom of the attack solution space.

A. Infeasibility of Attacks
We first formulate a lemma that we will use later for deriving conditions for the infeasibility of attacks.where W 1i ,i∈{1,•••,q}, u i ,i∈{1,•••,q}, and s are complex numbers, and |u i |=1.Then there exists exactly one solution if and only if one or both of the two conditions is satisfied, where Proof.The proof is given in the appendix.
We can now formulate the following result regarding the infeasibility of attacks.
Proposition 1.Consider an attack on P ≥ 2 PMUs such that rank(W )=1.Then an undetectable TSA does not exist if and only if at least one of the following conditions is satisfied. where Proof.Recall that every attack corresponds to a solution to (3).We know that (3) has at least one solution (u i = 1,∀i ∈ {1,•••,P }), i.e., no attack.If this solution is the only feasible solution to (3), then there does not exist a feasible attack.By comparing equations ( 3) and ( 4) we observe that they are identical when q = P and s = 0. Therefore, we can use Lemma 3 with q = P and s = 0 to obtain the necessary and sufficient conditions (7) and (8).
Note that (7) holds if and only if W 1i are co-directed.Since W is Hermitian, i.e., W ii is real, therefore W 1i must be real, and thus W is a real symmetric matrix.Similarly, (8) holds if and only if all W 1i ,i = i * point in the opposite direction of W 1i * , which implies also that all W 1i are real, and hence W is real.This is extremely improbable in practice due to the noisy nature of the complex-valued measurements affecting W .

B. Conditions for finding a non-empty attack solution set
Based on the previous discussion, we can now formulate the following result for the existence of undetectable TSAs.
Proposition 2. Consider that rank(W ) = 1 and P ≥ 2. Then an undetectable TSA exists if and only if where Proof.Observe that it is impossible to have w(P,0)> We know from Proposition 1 that equality on at least one side of ( 10) is a necessary and sufficient condition for the non-existence of attacks.Therefore, having strict inequality on both sides, i.e., (9), is a necessary and sufficient condition for the existence of undetectable TSAs, which proves the proposition.

C. Degrees of freedom of the set of undetectable TSAs
In order to characterize the set of undetectable TSAs, we rely on the notion of the degrees of freedom of the solution of a system of equations.In our work, the system of equations consists of equation ( 3) as well as the constraint |u p | = 1, p ∈ {1,...,P }, in the variables (u 1 ,...,u P ).Let us denote by U ⊆ Π P p=1 S 1 the set of feasible solutions, where S 1 is the circle group.Definition 2. We say that U has n degrees of freedom around u ∈ U if there is a nonempty open set Θ ⊂ R n and an injective mapping ψ :Θ→U such u=ψ(α) for some α∈Θ.
As an example, S 1 has one degree of freedom around any point u ∈ S 1 .Without loss off generality, we can consider that the attacker starts to compute a solution u ∈ U by choosing u P ∈ U P , where U P ⊆ S 1 is the set of all feasible values of u P (c.f., the bold arcs in Fig. 1).For the chosen u P , the attacker then chooses u P −1 ∈U P −1 (u p ), etc.In general, when choosing u i ,i∈{1,...,P }, let us denote by u + i = (u i+1 ,...,u P ) the values already chosen, by u − i = (u 1 ,...,u i−1 ) the values that remain to be chosen, and by U i (u + i )⊆S 1 the set of all feasible values of u i given earlier choices u + i such that there is Next we define the notion of free variables and leading variables, tightly related to the degrees of freedom of the solution set.
) is a free variable w.r.t. the solution of (3) if its value could be arbitrarily chosen from a set of feasible values.The degrees of freedom of the solution of (3) is the number of arbitrarily chosen variables whose values can be fixed (i.e., free variables) so that the values of the remaining variables (i.e., leading variables) can be uniquely determined to satisfy (3).In other words, the solution set U has n degrees of freedom around every interior point u∈U if there are n free variables in (u 1 ,•••,u P ).That is, given an arbitrary ordering of variable choice, the first n chosen variables (i.e., (u P −n+1 ,...,u P )) will be free variables, and the remaining variables (i.e., (u 1 ,...,u P −n )) will be leading variables.
In what follows we characterize the degrees of freedom of the attack solution set for rank(W )=1.
Theorem 1.Consider that rank(W ) = 1 and a TSA is feasible.Then the solution set U has P −2 degrees of freedom around every interior point u∈U.
In order to prove this result, we first introduce two lemmas.Lemma 4. If an undetectable TSA is feasible for P ≥ 3 then u P is a free variable.
Proof.Let us denote by U P the set of feasible values for u P , excluding the boundary values corresponding to intersection points in Figure 1.Furthermore, let us define the set of feasible attack angles Θ P = {α P : u P = e jα P ,u P ∈ U P }.Note that the mapping ψ P :Θ P →U P is injective, as required by Definition 3. By Proposition 1 we know that it is with negligible probability that AR 1,P −1 and O P intersect in only one point.Furthermore, Θ P is either a non-empty connected bounded open interval (e.g., if O P =O 1 P in Figure 1), or the union of two such intervals (e.g., if O P = O 2 P in Figure 1), both cases satisfying the requirement on Θ i in Definition 3. Therefore, u P is a free variable according to Definition 3.This establishes that for P ≥3 the attack solution has at least one degree of freedom corresponding to u P being a free variable.Next, let us consider that u + q = (u q+1 ,•••,u P ), q ≥ 3 are already chosen.The next lemma establishes the conditions for u q ∈ U q (u + q ) to be a free variable.
Then u q is a free variable if and only if where s j =−W 1j (u j −1).
Proof.The proof is given in the appendix.
We can now proceed to prove Theorem 1.
Proof of Theorem 1.By Lemma 4, the first chosen variable (i.e., u P ) is always a free variable.Next, Lemma 5 establishes that u P −1 ∈U(u P ) is also a free variable given that u P is chosen such that the intersection point between AR 1,P −1 and O P does not lie on the inner or outer bounding circle of the annular region (i.e., the end points of the bold arcs in Figure 1).Similarly, Lemma 5 states that the intersection point between AR 1,q and O q+1 (corresponding to u q+1 ) must not lie on the inner or outer bounding circle of the annular region Therefore, choosing suitable values for u + q (i.e., values that do not lie in the inner or outer circles) ensures that u q , q ∈{3,•••,P −1} are free variables.Note however that u 2 and u 1 cannot be free variables because AR 1,1 becomes a circle, and thus it can intersect with O 2 in at most two points.Since u P is also a free variable by Lemma 4, the number of free variables is P −2.
Observe that the order in which the terms are moved from the left hand side to the right hand side of (3) determines which P −2 variables in {u 1 ,...,u P } will be free variables and which 2 variables will be leading variables.

V. MINIMUM COST RANK-1 TSA MITIGATION PROBLEM
In this section we formulate the problem of mitigating undetectable TSAs at minimum cost, making use of Theorem 1. Recall that our approach mitigates TSAs involving message spoofing or injection against sets of PMUs with rank(W ) = 1.In principle, undetectable TSAs could still be possible even when rank(W ) > 1.However, from our experience, identifying sets of PMUs for which rank(W ) > 1 and computing undetectable TSAs for those usually requires the attacker to manipulate a larger number of PMUs.Hence attacks against PMUs with rank(W )=1 are easiest to perform.Moreover, the sets of PMUs for which rank(W ) > 1 in most cases contain subsets of PMUs that have a rank-1 W matrix.Therefore, mitigating rank-1 TSAs would typically mitigate rank-k attacks, for k >1.
Together with Lemma 2, Theorem 1 characterizes the minimum number of PMUs that an attacker needs to manipulate in order to launch an undetectable TSA.The theorem implies that for an equivalence class C i , C i = |C i |, any subset of 3 ≤ P ≤ C i time references in C i can be attacked in an undetectable manner (i.e., 1 or more degrees of freedom).Furthermore, an attack against any subset of P =2 time references (i.e., zero degrees of freedom) can be constructed, but due to the lack of degrees of freedom in the solution, the attack can be detected by BDD in practice, as shown in [26].Therefore, it is sufficient to only consider the mitigation of TSAs against P ≥3 time references if BDD is used in the system.

A. Securing Time References
We consider a power system operator that wants to upgrade its network infrastructure to mitigate TSAs.We assume that LSE based on (1) is used with BDD, and hence the objective is to mitigate attacks against the collection {C 1 ,...,C C } of equivalence classes of time references vulnerable to undetectable TSAs, as defined in Section III.We consider mitigation through authenticated time synchronization, e.g., using PTPv2.1.Upgrading a device (vertex) to PTPv2.1 means that the device will be able to use authenticated PTP messages.An edge is then secured if both incident vertices are upgraded.Therefore, we call a path in G to be secured if and only if all edges along the path are secured, and hence all vertices along the path are upgraded.Based of the previous discussion, the time reference of a PMU τ t can be either absolutely or relatively secured, defined as follows.
Definition 5. A time reference τ t ∈ C i is absolutely secured if a path in G from the root vertex v r ∈V (PTP master) to τ t is secured, including all intermediate vertices.Definition 6.Two time references τ t ,τ c ∈C i are relatively secured (denoted by τ t ↔ τ c ) if a path in G between τ t and τ c is secured, including all intermediate vertices.
In practice, securing a path means that the network equipment on the path have to be upgraded to support authenticated PTP messages, and related key management.Absolutely securing a time reference τ t mitigates any TSA including τ t .On the contrary, having τ t ↔τ c makes that a TSA must impose the same phase angle shift on the corresponding PMU measurements.Observe that absolutely securing τ t can be thought of as relatively securing τ t and v r .We can thus define the augmented equivalence classes C + i =C i ∪{v r },∀i∈{1,...,C}.In what follows we show the effect of relatively securing time references in an equivalence class C i on the rank of the resulting W matrix, and hence the vulnerability of C i to undetectable TSAs.Recall that attacking time references in C i would always result in a rank-1 W matrix.However, relatively securing two time references would make them dependent, thus effectively reducing the number of independent time references in C i by one.However, in order to make sure the results in Theorem 1 apply to the reduced set of time of references (i.e., that they form a vulnerable equivalence class), we have to ensure that attacking any subset yields a rank-1 W matrix, which is proven by the following Proposition.In the proposition, we consider an attack against P time references affecting m ≥ P measurements.Furthermore, we say that the matrix Ψ is targeting the set of measurements M a (denoted by Ψ → M a ) if and only if Ψ i,j =0,∀i / ∈M a , and ∀i∈M a ∃j ∈{1,...,P } s.t.Ψ i,j =1.Proof.
then rank(V ) = 1 as well, since W = V † V and rank(A) = rank(A † A) for any matrix A. Now observe that V can be also written as V =G 0 Ψ 0 , where G 0 is defined as G 0 =F diag(z) 0 s.t.diag(z) 0 ∈C M×m is the sub-matrix of diag(z) including only the column indexes in M a , and Ψ 0 ∈ {0,1} m×P is the sub-matrix of Ψ including only the rows with indexes in M a .
Since Ψ corresponds to attacking one measurement per time reference, then w.l.o.g.we can assume that Ψ 0 =I P , i.e., the P ×P identity matrix.This implies that rank(V )=rank(G 0 )=1.Given that rank(AB) ≤ min(rank(A),rank(B)) for any two matrices A and B, the rank of any V = G 0 Ψ0 = G Ψ s.t.Ψ0 ∈ {0,1} m× P , 1 ≤ P ≤ m, and Ψ0 → M a will be at most 1.Then it follows that rank( W ) ≤ 1, since W = V † V .Note that rank( W ) = 0 if and only if the sum of each row in G 0 is 0, which is extremely unlikely due to the noisy nature of z.Therefore rank( W )=1.
The above result illustrates the effect of relative securing of time references of the same equivalence class.Since relative securing of two PMUs forces a TSA to have the same time shift for both PMUs, the time reference of these PMUs will not be independent.Therefore, this corresponds to having a modified Ψ matrix targeting the same set of measurements.The proposition then shows that if the rank of the original W matrix is equal to 1 then the rank of the W matrix obtained using the modified Ψ matrix will also be equal to 1. Therefore, for an equivalence class C i with m measurements, rank(W ) = 1 holds independent of the number of time references in C i that are relatively secured.Moreover, given that the number of degrees of freedom of a rank-1 TSA against C i is P −2 (Theorem 1), and that a rank-1 TSA with P =2 is detectable by BDD, relatively securing two time references in C i decrements P by 1, and similarly decrements the number of degrees of freedom by 1.As a consequence, relatively securing time references in C i , s.t., there are only P = 2 independent time references will render any attack against C i detectable by BDD.
The optimal solution is highlighted in red, while the optimal solution when considering only absolutely securing is highlighted in blue.

B. TSA Mitigation Problem Formulation
In what follows we formulate the problem of mitigating TSAs at minimum cost.We define the cost as the number of network equipment and PMUs (vertices in V) that need to be upgraded to support authenticated PTP messages.In the formulation we make use of the collection k :τ 1 ↔τ 2 ,∀C +4 k ∈C +4 .In the above problem formulation, observe that mitigation relies on identifying the collection C +4 of vulnerable quadruplets, which is based on the collection C of equivalence classes.In turn, equivalence classes in C include only PMUs for which the pairwise IoS* metric is equal to 1, which holds only when the corresponding W matrix is rank-1.Therefore, by definition, MIN-TM attempts only to mitigate those attacks corresponding to rank(W )=1.Observe also that the above problem formulation is based on that in order to mitigate all undetectable TSAs with rank(W )=1, each equivalence class can contain at most two independent time references.
In the following we present an Integer Linear Programming (ILP) formulation of the problem, generalizing an ILP formulation proposed for the Group Steiner Tree (GST) problem [28].The formulation is based on converting G into a directed graph and then solving a max-flow problem, as demonstrated in Fig. 2. The first step in the conversion is to replace each undirected edge e∈E with two directed edges, one in each direction.The next step is to add two sets V s = {v s,1 ,...,v s,K } (source vertices) and V t = {v t,1 ,...,v t,K } (terminal vertices), each consisting of K additional vertices, to G. Next, we add a directed edge from each vertex v s,k to each of its corresponding four vertices in C +4 k , and another directed edge from each vertex in C +4 k to its corresponding terminal vertex v t,k , as shown in Fig. 2. We let the augmented directed graph be G =(V ,E ) where (i,l)∈δ + (i) where y and y are decision variables that indicate whether a vertex or an edge is included G * , respectively.f is a decision variable indicating the flow from each source vertex in V s to its corresponding terminal vertex in V t on each directed edge in E , δ + (i) is the set of directed edges (i,l),∀l ∈ V originating from vertex i, δ − (i) is the set of directed edges (l,i),∀l ∈V terminating at vertex i, and d i,k is defined as A special case of the problem arises when considering absolute securing of time references only.The mitigation problem can then be solved by absolutely securing at least C i − 2 time references per equivalence class C i [25].However, solving this problem is not guaranteed to minimize the mitigation cost, since the entire paths from vertices in C i to v r have to be secured, resulting a cost at least as high as when securing paths between pairs of vertices in each C i is allowed.This difference is demonstrated in Fig. 2, by showing the optimal solution of [25] (blue) and that of MIN-TM (red).Proof.We prove NP-hardness through reduction from the case when considering absolute securing only (we call it MIN-TM-A), which was shown to be NP-hard [25].Given an instance of MIN-TM-A with G =(V,E) and C +4 , we construct an instance of MIN-TM with G = (V,E) and k,l is a leaf vertex, and v k,l is the parent of v  Select (12) with constr.(13) and V k ← vertices on the path from v s * k and v t * k 12: V * ←V * ∪V k 13: end for

C. Proposed Mitigation Algorithms
In what follows we propose two approximation algorithms for solving MIN-TM.Both algorithms solve the problem greedily for each set of vertices C +4 k ,k ∈ {1,...,K}, and then combine the respective results.We denote by c * the cost of the optimal solution for MIN-TM.
Shortest Path Greedy (SP-Greedy): For each C +4 k ∈ C +4 choose the shortest path among all 4  2 = 6 paths between pairs of vertices in C +4 k .The output of the algorithm is the union of all vertices in all K chosen shortest paths.We denote by c SP the cost of the solution achieved by SP-Greedy.Proposition 5. SP-Greedy is a K-approximation of MIN-TM.i.e., c SP ≤Kc * Proof.Consider that the cost (number of vertices) of the shortest path chosen by SP-Greedy for C +4 k is β k , k ∈ {1, ... , K}.Now consider the worst case that there exists a path (that was not chosen by SP-Greedy) that has a cost of max k∈{1,...,K} β k , and that this path by itself solves MIN-TM.This is indeed the worst case, since SP-greedy would have chosen this path if it had lower cost.The optimal solution for MIN-TM in this case has a cost of c * = max k∈{1,...,K} β k .The solution by SP-Greedy would have a cost of c SP = K k=1 β k , which can be upper-bounded by Kmax k∈{1,...,K} β k =Kc * .
The second proposed algorithm includes solving a linear relaxation of ILP (12).In order to avoid trivial solutions (solutions with zero cost, e.g., using only edges (v s ,v 6 ) and (v 6 ,v t ) in Fig. 2) of the resulting LP, we complement the formulation with the extra constraint where Linear Programming Greedy (LP-Greedy): Algorithm 1 shows the pseudo code of LP-Greedy.Among all 4  2 paths between pairs of vertices in each C +4 k ∈C +4 , the algorithm chooses the path with the highest flow in the fractional solution of the LP relaxation of ILP (12), along with the extra constraint (13).The output of the algorithm is the union of all vertices in the K paths chosen.We denote by c LP the cost of the solution achieved by LP-Greedy.Proposition 6. LP-Greedy is a 4K-approximation of MIN-TM.i.e., c LP ≤4Kc * Proof.Solving the LP for the first time yields a cost that is at most i by a factor of 4 in the worse case.Since this constraint is added K times for each set in C +4 , then the cost of the resulting solution c LP will be at most 4Kc * .
The computations needed to implement the proposed mitigation scheme are shown in Algorithm 2.

VI. NUMERICAL RESULTS
In what follows we evaluate the performance of both proposed approximation algorithms through extensive simulations on synthetic graph topologies and on IEEE benchmark power systems.To quantify how well the proposed algorithms solve MIN-TM, we compared their performance in terms of the mitigation cost to the optimal solution obtained by brute-force search (denoted by Brute-Force) for small graphs.To evaluate the performance on larger problem instances, we compared the algorithms to the optimal fractional solution of the linear relaxation of ILP (12) along with the constraint (13).Note that this approach only provides a lower bound on the optimal solution, and is thus denoted by Optimal-LB.We also compare the performance of the proposed algorithms to the greedy algorithm proposed in [25], which works as SP-Greedy but considering only absolute securing of time references.Since v r will be included in each shortest path, the output will always be a tree, and thus we refer to this algorithm as SP-Greedy-T.In our simulations, we further optimize SP-Greedy, LP-Greedy, and SP-Greedy-T by removing all shortest paths for quadruplets that were already secured by earlier paths.All simulations were carried out on a notebook with Intel Core i7-8550 CPU @ 1.8 GHz and 16 GB of RAM.

A. Performance on Synthetic graphs
To generate random tree graphs with |V| vertices, we randomly choose |V| − 1 edges from all possible |V| 2 edges, such that the graph connectivity is ensured.We then choose v r to be the vertex with the highest betweenness-centrality in G.To simulate equivalence classes, we randomly choose C disjoint subsets of V to form the collection C, considering two scenarios for the cardinality of the equivalence classes: (I) C = 3, C i ∼ U (3,4), and (II) C =5,C i ∼U (3,7), where U(a,b) is a discrete uniform distribution defined by the interval [a, b].Fig. 3 shows the mitigation cost (Fig. 3a) and the execution time (Fig. 3b) for the considered methods to solve MIN-TM, as a function of the number of vertices |V| in G.Each point represents an averaged value from 200 simulations, along with 95% confidence intervals.From Fig 3a we first observe that Optimal-LB is a loose lower bound on the optimal solution (Brute-force) for smaller graphs in Scenario I, suggesting that the proposed algorithms actually perform much better than indicated by their distance to Optimal-LB.Compared to a naive approach of securing all vertices (i.e., mitigation cost = |V|, assuming T =V), the proposed algorithms (SP-Greedy and LP-Greedy) achieve cost savings ranging between 50% and 89% for Scenario I, and between 44% and 70% for scenario II.Compared to SP-Greedy-T [25], we observe that the proposed algorithms perform similarly.In fact, SP-Greedy-T performs slightly better, especially for Scenario II when C and C i (and hence K) were larger.SP-Greedy-T favors shortest paths that are intersecting, which is expected to reduce its cost when vertices in C are placed randomly in the graph.
Fig. 3b confirms that the average execution time of SP-Greedy, LP-Greedy, and SP-Greedy-T scale polynomially with |V|, in contrast to Brute-Force which has an exponential execution time.Moreover, the execution time of SP-Greedy-T was consistently one order of magnitude lower than SP-Greedy for both Scenario I and II.The execution time of LP-Greedy was the highest among the polynomial time algorithms, one to two orders of magnitude higher than SP-Greedy.

B. Performance on IEEE Benchmark Power Systems
We used power system topology information in the MATPOWER package [29] in Matlab to evaluate the proposed mitigation scheme on the 118-bus, the 145-bus, and the 300-bus IEEE benchmark power systems.We used the effective rank ratio metric introduced in [30] for computing the equivalence classes C (and hence, C and C i ).For each of the considered benchmark power systems, we chose V to be the set of buses in the system.The set E of edges was chosen to be a random subset of connections between buses, such that G =(V,E) forms a tree.Similar to synthetic graphs, v r was chosen as the vertex with the highest betweenness-centrality in G. Furthermore, we set the set T of time references to be the buses (vertices) that have a PMU installed.In a practical deployment, each bus with a PMU may correspond to multiple vertices, e.g., one for the PTP switch in the corresponding substation, and one for the PMU itself, but this simplifcation does not change the solution to the MIN-TM problem.
Fig. 4 shows the mitigation cost (Fig. 4a) and the execution time (Fig. 4b) for the considered methods, as a function of the number of measurements M deployed in the system.Each point represents an averaged value from 100 deployments of M voltage and current injection phasor measurements, along with 95% confidence intervals.The M measurements were located at random, given that they ensure the observability of the power system and the absence of critical measurements.That is, the rank of the measurement matrix H and all its possible M − 1 × N sub-matrices is N (i.e., they have full column rank).Note that we do not deploy PMUs on zero-injection buses (i.e., buses that do not have either generators or loads, as  specified in the MATPOWER package).The number of measurable buses N * (i.e., not zero-injection buses) for the IEEE 118-bus, 145bus, and 300-bus systems was 108, 61, and 224, respectively.Fig. 4a shows higher cost savings achieved by the proposed algorithms on IEEE benchmark systems compared to synthetic graphs.A naive approach that secures all time references will have an average mitigation cost of |T | = 0.75N * when M = N (sparsest possible PMU deployment).Compared to that, SP-Greedy and LP-Greedy can reduce the cost by at least 87%, 95%, and 76% for the IEEE 118-bus, 145-bus, and 300-bus systems, respectively.Moreover, SP-Greedy and LP-Greedy significantly outperform SP-Greedy-T for all three IEEE benchmark systems, reducing the cost by at least 40%, 53%, and 30% for the IEEE 118-bus, 145-bus, and 300-bus systems, respectively.In this realistic setting, an equivalence class C i typically consists of vertices that are closer to each other than they are to v r , and thus relatively securing time references in C i is more efficient than absolute securing.Moreover, even though the performance of the proposed algorithms is not close to Optimal-LB, Fig. 3a suggests that the distance to the optimal solution achieved by e.g., Brute-Force would be much smaller.Overall, Fig. 4a demonstrates an important trade-off between the cost of deploying extra measurements in the system and that of securing time references.For the execution time, Fig. 4b shows that SP-Greedy-T is slightly more time efficient than SP-Greedy, while LP-Greedy takes signifcantly higher execution time (upto 3 orders of magnitudes) than that of SP-Greedy.Overall, the results suggest that SP-Greedy performs well in terms of both mitigation cost and execution time compared to its competitors.When time complexity is not an issue, the operator can simply choose the algorithm that yields the smallest cost.

VII. CONCLUSION
We considered the problem of mitigating undetectable TSAs against power systems, by using a combination of LSE and PTP authentication.We provided necessary and sufficient conditions for the existence of undetectable TSAs that implies that a TSA is undetectable if and only if at least 3 time references are manipulated.Based on that, we then formulated the problem of mitigating undetectable TSAs at minimum cost, and proposed two approximation algorithms for solving it.Beside being computationally efficient, the proposed algorithms were shown to result in huge cost savings compared to the naive mitigation approach, and significant savings compared to previous work, both for synthetic graph topologies and for realistic IEEE benchmark systems.
Observe that the existence of only one solution to (4) implies that there can be only one intersection point between the circle and the annular region.Proving the inverse, i.e., that the existence of only one intersection point implies that only one solution to (4) exists, will be shown later.
The annular region and the circle intersect in only one point in three cases, depending on the relative locations of O q and AR 1,q−1 .
1) Case A: In the first case, the circle O q is outside the annular region, thus the distance between the center of the annular region and the center of the circle equals the sum of the outer radius of the annular region and the radius of the circle, i.e., |c 1,q−1 −c q |=r o 1,q−1 +r q .
Substituting the definitions for O q and AR 1,q−1 into ( 14) we obtain which yields (5).
2) Case B: The second case arises when the radius of the circle O q is equal to the sum of (1) the distance between the centers of the annular region and the circle, and (2) the outer radius of the annular region, i.e., Substituting into (15) we obtain |W 1i | = w(q,s).( 16) 3) Case C: The third case arises if the inner radius of the annular region is equal to the sum of (1) the distance between the centers of the annular region and the circle, and (2) the radius of the circle, i.e., r i 1,q−1 = |c 1,q−1 −c q |+r q r i 1,q−1 −r q = |c 1,q−1 −c q |.
Recall that from case A, (5) is a sufficient condition for having one intersection point between O q and AR 1,q−1 .Furthermore, observe that ( 16) and ( 18) can be combined as |W 1i | = w(q,s), where i * = argmax i∈{1,•••,q} |W 1i |, which proves that ( 6) is also a sufficient condition for having one intersection point.Therefore, only one intersection point can be found if either of ( 5) or ( 6) holds.Furthermore, the three cases A, B and C are the only cases where the circle has only one intersection point with AR 1,p−1 .Therefore, if O q and AR 1,q−1 intersect in only one point, then either (5) or ( 6) must hold.
To prove that ( 5) and ( 6) are necessary and sufficient conditions for having only one solution to (4), we need to show that there exists only one solution to (4) if and only if there exists only one intersection point between AR 1,q−1 and O q .Recall that a unique solution to (4) implies that only one intersection point exists.It remains to prove the opposite, i.e., that the existence of only one intersection point between AR 1,q−1 and O q implies the existence of only one solution to (4).To prove that, that we consider the conditions ( 5) and ( 6) for having only one intersection point.
Notice that since s = q i=1 W 1i (u i −1), we can rewrite ( 5) and ( 6) as Observe that since |u i | = 1, then ( 19) cannot be satisfied unless all W 1i u i ,i ∈ {1,••• ,q} point in the same direction, which is the direction of their sum.In other words, they point in the direction of ( q i=1 W 1i )+s.Therefore for all j ∈ {1,••• ,q}.It is clear that (21) yields only one solution (u 1 ,•••,u q ), and hence there exist only one solution to (4).
On the other hand, (20) cannot be satisfied unless all W 1i u i ,i ∈ {1,•••,q}\{i * } point in the same direction, which is the opposite direction to W 1i * u i * and ( q i=1 W 1i )+s.Therefore Again, it is clear that (22) yields only one solution (u 1 ,•••,u q ), and hence there exists only one solution to (4).
Proof of Lemma 5.For computing u q we can rearrange (3) to q−1 i=1 W 1i (u i −1)=−W 1,q (u q −1)+ where s j = W 1j (u j − 1).The solution of (23) corresponds to the intersection between an annular region AR 1,q−1 (the LHS) and a circle O q , with the same definitions as in Lemma 3, for s = P j=q+1 s j .Similar to Lemma 4, having more than one intersection point is equivalent to that the feasible set of attack angles Θ q (θ + q ) = {θ q : u q = e jθq ,u q ∈ U q (u + q )} is a non-empty connected bounded open interval or the union of two such intervals, and thus u q is a free variable.
Therefore, u q will not be a free variable if (23) has one solution or no solution.To identify the conditions for which (23) has one Furthermore, we denote by O P = (c P ,r P ) the circle with center c P = W 1P and radius r P = |W 1P |.

Figure 1
shows an illustration of the geometrical problem and the solution to be computed (shown in red) for two different cases of O P .Based on the previous discussion, let us recall the following result from[26].
This article has been accepted for publication in IEEE Transactions on Control of Network Systems.This is the author's version which has not been fully edited and content may change prior to final publication.Citation information: DOI 10.1109/TCNS.2024.3368316This work is licensed under a Creative Commons Attribution 4.0 License.For more information, see https://creativecommons.org/licenses/by/4.0/

Fig. 2 :
Fig. 2: An illustration of converting the original graph G (to the left) into the extended directed graph G (to the right), for a simple case of K = 1 vulnerable quadruplets, withC +4 1 = {v r ,v 2 ,v 4 ,v 6 }.The optimal solution is highlighted in red, while the optimal solution when considering only absolutely securing is highlighted in blue.
vulnerable quadruplets of time references, where k ∈ {1,...,K} and K = |C +4 | = C i=1 Ci 3 .The problem of mitigating undetectable TSAs at minimum cost can then be formulated as follows.Minimum-Cost TSA Mitigation (MIN-TM): Given the communication infrastructure graph G = (V,E) of the WAMPAC and a collection C +4 of vulnerable quadruplets of time references.Find a subgraph G * = (V * , E * ) of G with minimum |V * | s.t.∃τ 1 ,τ 2 ∈C +4 This article has been accepted for publication in IEEE Transactions on Control of Network This is the author's version which has not been fully edited and content may change prior to final publication.Citation information: DOI 10.1109/TCNS.2024.3368316This work is licensed under a Creative Commons Attribution 4.0 License.For more information, see https://creativecommons.org/licenses/by/4.0/V=V ∪V s ∪V t , and |E |=2|E|+4K +4K.The objective is to find a sub-graph with minimum number of vertices that includes a path from each source vertex v s,k to its corresponding terminal vertex v t,k (excluding trivial solutions using only two edges), yielding the ILP min

Proposition 4 .
The MIN-TM problem is NP-hard.

( 1 )
k,l .Let the set of extra edges be denoted by E k,l .Furthermore, let and C +4 = {C +4 1 , ..., C +4 K }.Solving this instance of MIN-TM is guaranteed to solve MIN-TM-A since paths from any vertex in C +3 k to v r will be shorter than paths between any two vertices in C +3 k .The solution for MIN-TM-A can be reconstructed by removing the added vertices V k,l and edges E k,l .Algorithm 1 LP-Greedy input: Collection of vulnerable quadruplets C +4 output: Set of secured vertices V * 1: V * ←φ 2: C ←φ (set of extra constraints) 3: f ← Solve LP relaxation of (12) with constraint (13) 4: for C +4 k ∈C +4 do 5: This article has been accepted for publication in IEEE Transactions on Control of Network Systems.This is the author's version which has not been fully edited and content may change prior to final publication.Citation information: DOI 10.1109/TCNS.2024.3368316This work is licensed under a Creative Commons Attribution 4.0 License.For more information, see https://creativecommons.org/licenses/by/4.0/
This article has been accepted for publication in IEEE Transactions on Control of Network Systems.This is the author's version which has not been fully edited and content may change prior to final publication.Citation information: DOI 10.1109/TCNS.2024.3368316This work is licensed under a Creative Commons Attribution 4.0 License.For more information, see https://creativecommons.org/licenses/by/4.0/

TABLE I :
Table of Notation i Set of time references in equivalence class i, |C i |= C The pairs of time references that are vulnerable to undetectable TSAs for P = 2 form equivalence classes C = {C 1 ,...,C C }, where C =|C| is the number of equivalence classes.
• If P =2 then there is 1 undetectable TSA.• • For P ≥ 3, undetectable TSAs against T a exist if ∃C i ∈ C s.t.T a ⊆ C i .In addition, the set of undetectable TSAs (i.e., This article has been accepted for publication in IEEE on Control of Network Systems.This is the author's version which has not been fully edited and content may change prior to final publication.Citation information: DOI 10.1109/TCNS.2024.3368316Proposed Mitigation Scheme input: Measurement matrix H, network graph G =(V,E), and root vertex v r output: Set of secured vertices V This work is licensed under a Creative Commons Attribution 4.0 License.For more information, see https://creativecommons.org/licenses/by/4.0/Algorithm 2 * 1: C ← Equivalence classes computed using IoS* (Lemma 2) 2: