extfieldsofstudy:
Computer Science
Mathematics

text:
Second-Order Asymptotics of Hoeffding-Like Hypothesis Tests

We consider a binary statistical hypothesis testing problem, where $n$ independent and identically distributed random variables $Z^n$ are either distributed according to the null hypothesis $P$ or the alternate hypothesis $Q$, and only $P$ is known. For this problem, a well-known test is the Hoeffding test, which accepts $P$ if the Kullback-Leibler (KL) divergence between the empirical distribution of $Z^n$ and $P$ is below some threshold. In this paper, we consider Hoeffding-like tests, where the KL divergence is replaced by other divergences, and characterize, for a large class of divergences, the first and second-order terms of the type-II error for a fixed type-I error. Since the considered class includes the KL divergence, we obtain the second-order term of the Hoeffiding test as a special case.

I. INTRODUCTION
Statistical hypothesis testing is known to have applications in areas such as information theory, signal processing, and machine learning. The most simple form of hypothesis testing is binary hypothesis testing, where the goal is to determine the distribution of a random variable Z between a null hypothesis P and an alternate hypothesis Q. There can be two types of errors in a binary hypothesis testing: The type-I error is the probability of declaring the hypothesis as Q when the true distribution is P . The type-II error is the probability of declaring the hypothesis as P when the true distribution is Q. In general, we are interested in analyzing the trade-off between these two types of errors for a given test. This test may have full or only partial access to the distributions P and Q.
When both P and Q are known, the likelihood ratio test (also known as Neyman-Pearson test [1]) achieves the optimal error trade-off between type-I and type-II errors. The Neyman-Pearson test is also investigated in an asymptotic setting, where one observes n independent copies of Z and the errors are analyzed asymptotically as n tends to infinity. In this case, we are often interested in the behaviour of the type-II error β n for a fixed type-I error α n . It is known that for α n ≤ ǫ, ǫ ∈ (0, 1), the corresponding type-II error satisfies [2,Prop. 2.3] − ln β n = nD(P Q) − nV (P Q)Q −1 (ǫ) + o( √ n), as n → ∞ where D(P Q) is the Kullback-Leibler (KL) divergence between P and Q, defined as Q −1 (·) denotes the inverse of the tail probability of the standard Normal distribution, and o( √ n) denotes a term that grows more slowly than √ n as n → ∞. Thus the first-order term of − ln β n , defined as lim n→∞ − ln βn n and sometimes referred to as the error exponent, is given by D(P Q). The second-order term, defined as − ln βn−nD(P Q) √ n , is given by − V (P Q)Q −1 (ǫ). The case where the test has only partial access to the distributions P and Q is generally studied under the name of composite hypothesis testing. Here, the assumption is that P and Q belong to certain classes of distributions denoted by P and Q, respectively. Several special cases of composite hypothesis testing have been investigated in the literature. Two such cases are: 1) P is known and Q can be any distributions but P ; 2) P is known and Q belongs to a special class Q. A test proposed by Hoeffding [3], known as the Hoeffding test, is suitable for these two cases. In the Hoeffding test, the null hypothesis P is accepted if the KL divergence between the type t Z n (the empirical distribution) of the observations Z n = (Z 1 , · · · , Z n ) and P is below some threshold c. Otherwise, the alternative hypothesis is accepted. Mathematically, if D(t Z n P ) < c, then accept H 0 ; otherwise accept H 1 .
In [4,Th. III.2], it was shown that the Hoeffding test achieves the same first-order term of − ln β n as the Neyman-Pearson test. Consequently, not having access to the distribution of the alternate hypothesis does not affect the first-order term. To analyze the second-order term of − ln β n in the above two cases of composite hypothesis testing, Watanabe proposed a test that is second-order optimal in some sense [5].
The Hoeffding test does not require side information about the special class Q to which Q belongs. It is therefore perhaps not surprising that tests that have knowledge of Q (case 2 of composite hypothesis testing) may outperform the Hoeffding test. Such tests include the generalized likelihood ratio test (GLRT) [6] and test via mismatched divergence [7]. For example, it has been observed through simulations that the test via mismatched divergence performs better than the Hoeffding test [7]. However, an analytical comparison between the Hoeffding test and other tests is missing since refined asymptotics of the Hoeffding are not available.
In this paper, we wish to understand the behavior of the second-order term of the Hoeffding test as well as of Hoeffding-like tests where the KL divergence is replaced by other divergences, in the following referred to as divergence test. For a large class of divergences, we demonstrate that the second-order term is given by− V (P Q)Q −1 χ 2 ,k−1 (ǫ), where Q −1 χ 2 ,k−1 (·) denotes the inverse of the tail probability of the chi-squared distribution with k − 1 degrees of freedom. This class of divergences includes the KL divergence, so we obtain the second-order term of the Hoeffding test, as a special case.
The rest of this paper is organized as follows. Section II presents notations and the problem formulation. Section III discusses our main results. Section IV presents the proof of our results. Section V concludes the paper with a summary and discussion of our results.

II. NOTATIONS AND PROBLEM FORMULATION A. Notations
Let f (x) and g(x) be two real-valued functions.
Let {a n } be a sequence of positive real numbers. We say that a sequence of random variables is X n = o p (a n ) if for every ǫ > 0, we have lim n→∞ P X n a n > ǫ = 0.
Similarly, we say that X n = O p (a n ) if for every ǫ > 0 there exist M ǫ > 0 and N ǫ ∈ N such that P X n a n > M ǫ ≤ ǫ, n ≥ N ǫ .

B. Divergence and Divergence Test
We consider a divergence test on a random variable Z that takes value in a discrete set Z = {a 1 , · · · , a k }, where k ≥ 2. We denote the probability distribution function of Z by a k-length vector P = (P 1 , · · · , P k ) T , where be the set of all such probability distributions. The set of all positive probability distributions with P i > 0, i = 1, · · · , k on Z is denoted by P(Z). Let Z n = (Z 1 , · · · , Z n ) denote an n-length i.i.d. sequence of Z. By denoting the realization of Z n by z n , we define the type distribution t z n = (t z n (a 1 ), · · · , t z n (a k )) T of the sequence z n as (t z n (a 1 ), · · · , t z n (a k )) T N a 1 |z n n , · · · , N a k |z n n T where N a i |z n , i = 1, · · · , k denotes the number of times a i occurs in the sequence z n . Since k i=1 P i = 1, a probability vector P = (P 1 , · · · , P k ) T ∈ P(Z) can be represented by any (k − 1) components of P . With this representation, the set P(Z) can be identified as a (k − 1)-dimensional manifold; see [8]. In this paper, we choose the coordinate system of the probability distribution P = (P 1 , · · · , P k ) T as the first (k − 1) components denoted by P = (P 1 , · · · , P k−1 ) T .
Given any two probability distributions P, Q ∈ P(Z), one can define a non-negative binary functioñ D(P Q) called a divergence which represents a measure of discrepancy between them. A divergence is not necessarily symmetric in its arguments and also need not satisfy the triangle inequality; see [8], [9] for more details. Mathematically, a divergence is defined as follows.
Definition 1 Consider two points P and Q in P(Z) whose coordinates are P = (P 1 , · · · , P k−1 ) T and Q = (Q 1 , · · · , Q k−1 ) T . A divergenceD(P Q) between P and Q is a smooth function of P and Q (we may writeD(P Q) =D(P Q)) satisfying the following conditions: 1)D(P Q) ≥ 0 for any P, Q ∈ P(Z).
There are many classes of divergences that are widely used in various fields of science and engineering; see, e.g., [10] for more details. A well-known example of a divergence is the f -divergence defined as follows: Let f : (0, ∈ ∞) → R be a convex function with f (1) = 0. Then, the f -divergence between P and Q is defined as When f (u) = u log u, D f is the Kullback-Leibler divergence. For the f -divergence is the α-divergence given by We next define a divergence test for case 1 of composite hypothesis testing (see Section I), That is, it is known that under hypothesis H 0 the distribution is P , while under hypothesis H 1 the distribution is anything but P . A Hoeffding-like test or divergence test TD n (r) for testing H 0 :P = P against the alternative H 1 :P = Q is defined as follows: Upon observing z n , ifD(t z n P ) < r for some r > 0, then the null hypothesis P is accepted; else P is rejected. For r > 0, define the acceptance region for P as AD n (r) z n |D(t z n P ) < r .
Then, the type-I and the type-II errors are given by α n (TD n (r)) P n AD n (r) c β n (TD n (r)) Q n (AD n (r)).
Our goal is to analyze the asymptotics of the type-II error when the type-I error satisfies α n ≤ ǫ, 0 < ǫ < 1.
Following the discussion in Section I, we define the first-order term β ′ and the second-order term β ′′ for the divergence test as follows: if the limits exist. For the Hoeffding test, it is known that β ′ = D(P Q). In this paper, we generalize this result to the divergence test for a class of divergences defined in Definition 2. We further obtain β ′′ for the divergence test.
III. SECOND-ORDER ASYMPTOTICS OF THE DIVERGENCE TEST We shall consider the following class of divergences.
Definition 2 Let Ξ denote the class of divergences satisfying the following conditions: 1) For any P, T ∈ P(Z), the second-order Taylor approximation ofD(T P ) around T = P is given byD for some constant η > 0 (that possibly depends on P ). In (7), T = (T 1 , · · · , T k−1 ) T , P = (P 1 , · · · , P k−1 ) T , and d χ 2 is the χ 2 -divergence defined as 2) For any P, T ∈ P(Z), the divergenceD satisfies the Pinsker-type inequalitỹ 3) The CDF of the statistic η −1 nD(t z n P ) can be approximated as for some positive sequence δ n satisfying lim n→∞ δ n = 0, where F χ 2 ,k−1 is the CDF of the chi-squared random variable χ 2 k−1 with k − 1 degrees of freedom. The following theorem characterizes β ′ and β ′′ of the divergence test for any divergenceD in Ξ.
Theorem 1 LetD ∈ Ξ and 0 < ǫ < 1. Consider the divergence test TD n (r) for testing H 0 : P against the alternative H 1 : Q, where P, Q ∈ P(Z) and P = Q. Then, the type-II error satisfies the following: Part 1: There exists a threshold value r n satisfying α n (TD n (r n )) ≤ ǫ (11) such that, as n → ∞, Part 2: For all r n > 0 satisfying (11), we have as n → ∞ where Q −1 Since the sequence δ n in (12) and (13) tends to zero as n → ∞, we have that O(max{δ n √ n, ln n}) = o( √ n). Consequently, (12) and (13) imply that, for the optimal threshold value r n , − ln β n (TD n (r n )) = nD(P Q) − nV (P Q)Q −1 Thus, Theorem 1 characterizes the first and second-order terms of the divergence test for anyD ∈ Ξ. It can be verified that Q −1 Consequently, the second-order term β ′′ of the divergence test is strictly smaller than the second-order term of the Neyman-Pearson test.
For certain divergences in Ξ, we can obtain more precise asymptotics than the one in (14) by tightening the o( √ n) term. To this end, we shall first discuss the three conditions stated in Definition 2, which we shall do one by one.

Condition 1:
WhenD is the f -divergence D f , it follows from a Taylor-series expansion of D f (T P ) around T = P that D f satisfies (7) with η = f ′′ (1) 2 [11,Th. 4.1]. Since both the α-divergence and the KL divergence D belong to the f -divergence class with f ′′ (1) = 1, we obtain thatD = {D α , D} satisfy (7) with η = 1 2 [12]. Furthermore, for the Renyi divergence of order α > 0, defined as There are many more divergences that satisfy (7), such as the Battacharya divergence, the Tsallis divergence, and the Sharma-Mittal divergence; see [10] for more details. In information geometry, such divergences are referred to as invariant divergences [10, Eq. 2.29].
Condition 2: It is well-known that the KL divergence D satisfies Pinsker's inequality [13, Lemma 11.6.1]. From [14,Th. 3], it further follows that Pinsker's inequality can be generalized to many fdivergences under some conditions on f . In particular, it follows from [14,Cor. 6] that the α-divergence satisfies a Pinsker-type inequality when 3 ≤ α ≤ 3, α = ±1: It follows from the same corollary that the Renyi divergence I α also satisfies (9). Furthermore, for any for some positive sequence a n such that a n → 0, then T −P 1 = O( √ a n ), which is sufficient to prove Theorem 1. In fact, condition 2 can be replaced by the condition that, ifD(T, P ) = O(a n ) for some positive sequence a n such that a n → 0, then T − P 1 = O( √ a n ). Condition 3: Note that this condition is a Berry-Esseen-type convergence of the statistic η −1 nD(t Z n P ) to the chi-squared distribution with k − 1 degrees of freedom. For both the α-divergence and the KL divergence, we have η = 1 2 . The following result implies that, for these divergences, (10) See [18] for more details.
From the above discussion, we conclude thatD = {D, D α } with −3 ≤ α ≤ 3, α = ±1 satisfy all three conditions in Definition 2 with δ n = 1 √ n and thus belong to the class Ξ. We thus have the following corollary to Theorem 1, which sharpens the third-order term in (14).
for testing H 0 : P against the alternative H 1 : Q, where P, Q ∈ P(Z) and P = Q. Then, the type-II error of TD n (r n ) minimized over all the threshold values r n satisfying α n (TD n (r n )) ≤ ǫ can be characterized as Note that, whenD is the KL divergence, the divergence test TD n (r) is the Hoeffding test. Hence, the above corollary recovers the second-order asymptotics of the Hoeffding test as a special case.
IV. PROOF OF THEOREM 1 A. Proof of Part 1 of Theorem 1 From (10), it follows that there exist M 0 > 0 and N 0 ∈ N such that For 0 < ǫ < 1, let Then, from (16), it follows that the type-I error is upper-bounded by ǫ: We next consider the type-II error for r n in (17). To this end, define Then, where the last step follows from [13,Th. 11.1.4]. In (18), P n denotes the set of types with denominator n and T (P ) is the type class ofP , defined as We next derive a lower bound on D(P Q) forP ∈ P n ∩AD(r n ). To this end, we shall use the following auxiliary results.
Lemma 3 For any two probability distributions Q = (Q 1 , · · · , Q k ) T and T = (T 1 , · · · , T k ) T , the secondorder Taylor approximation of the function T → D(T Q) around the null hypothesis T = P is given by Proof: See Appendix A. By the Pinsker-type inequality, we obtain that, for every T ∈ AD(r), we have T − P 1 ≤k √ r for some constantk. This implies that T − P 1 ≤k 1 √ r for some constantk 1 . Since any two norms on a finite-dimensional Euclidean space are equivalent, there exists a constant C 0 > 0 such that the Euclidean distance is bounded by Thus, for any T ∈ AD(r n ) with r n in (17), it follows that Then, for any T ∈ AD(r n ) with r n in (17), (7) and (20) and can be written as The expansion (23) implies that there exist M > 0 and N 1 ≥ N 0 such that, for all n ≥ N 1 and T ∈ AD(r n ), Next define We have the following lemma.
Lemma 4 Let r n be given in (17). Then, we have Proof: It follows from (25) that, for T ∈ AD(r n ) and n ≥ N 1 , ηd χ 2 (T, P ) − M n 3/2 ≤D(T P ) < r n . This implies that We next use (24) to lower-bound the term D(P Q) forP ∈ P n ∩ A D (r). Indeed, it follows from (24) that there exist M 2 > 0 andN 2 ∈ N such that, for all n ≥N 2 and T ∈ AD(r n ), Consequently, for all n ≥N 2 and T ∈ AD(r n ), where the last inequality follows since d χ 2 (T, P ) is non-negative.
To proceed further, let us consider the following function, defined for a probability distribution Γ = (Γ 1 , · · · , Γ k ) T ∈ P(Z) where For n ≥ N 1 andP ∈ P n ∩ AD(r n ), we have that where the last inequality follows from (26).
To compute the right-most minimum in (30), we define the following quantities: We then have the following result.
Remark 1 Since Γ * and P are probability distributions, by taking summation on both sides of (34), we obtain By the assumption of Theorem 1, we have that P = Q. Hence, D(P Q) − α i cannot be zero for all i = 1, · · · , k. This implies that there exist two distinct indices l, l ′ ∈ {1, · · · , k}, such that D(P Q)−α l < 0 and D(P Q) − α l ′ > 0. It follows that l ∈ I, which in turn yields that B is a non-empty set. Let Since r ′ n is vanishing as n → ∞, we can chooseN such that, for n ≥N , r ′ n satisfies √ r ′ n < √ V (P Q) τ with τ defined in (33). Then, (30) and (35) yield that, forP ∈ P n ∩ AD(r n ) and n ≥ N 2 max{N 1 ,N 2 ,N}, It follows that, for anyP ∈ P n ∩ AD(r n ), (27) can be further lower-bounded as Applying (39) to (18), the type-II error can be upper-bounded as β n (TD n (r n )) ≤ P ∈Pn∩AD(rn) where the last inequality follows since |P n | ≤ (n + 1) |Z| [13, Th. 11.1.1]. By taking logarithms on both sides of (40), we obtain ln β n (TD n (r n )) ≤ |Z| ln(n + 1) − nD(P Q) + n V (P Q)r ′ n + By (17) and (37), A Taylor-series expansion of √ x + δ around x > 0 yields then that where the last inequality holds since the function Q −1 χ 2 ,k−1 (·) is monotonically decreasing. It follows that (41) can be further upper-bounded as ln β n (TD n (r n )) ≤ |Z| ln(n + 1) − nD(P Q) + n V (P Q)r ′ n + Using a Taylor-series expansion of Q −1 Consequently, (44) can be expanded as ln β n (TD n (r n )) ≤ |Z| ln(n + 1) − nD(P Q) + nV (P Q) Q −1 This proves Part 1 of Theorem 1.

B. Proof of Part 2 of Theorem 1
Define for 0 < ǫ < 1 and n ∈ N, Lemma 6 For M 0 and δ n given in (16), we have Proof: For any r > 0, we can write P n D (t Z n P ) ≥ r = P n n ηD (t Z n P ) ≥ nr η .
Then, it follows from (16) that, for n ≥ N 0 and any r > 0, Substituting the value ofṙ n in (49), we get for n ≥ N 0 , Thus, by the definition of r ǫ n r ǫ n ≤ṙ n = η n Q −1 Furthermore, (50) implies that r ǫ n exists. Indeed, the set R ǫ n is lower-bounded by 0. Then, it follows from (50) that, for n ≥ N 0 , we haveṙ n ∈ R ǫ n . This implies that R ǫ n is a non-empty subset of the real numbers R that has a lower bound. Hence, by the completeness property of R, r ǫ n exists. We next derive a lower bound on r ǫ n . Leẗ Substituting the value ofr n in (49), we get for n ≥ N 0 , P n (D(t Z n P ) ≥r n ) ≥ ǫ + M 0 δ n > ǫ.
The above equation implies thatr n does not belong to the set R ǫ n . It follows from the property of the infimum that r ǫ n ≥r n = η n Q −1 Thus, from (50) and (52) we obtain η n Q −1 Taylor-series expansions of Q −1 χ 2 ,k−1 (ǫ + 2M 0 δ n ) and Q −1 For any T ∈ AD(r ǫ n ) with r ǫ n in (54), it follows from (21) that 13 Then, (7) and (20) can be written as The first equation (56) implies that there existM ′ 0 > 0 andN ′ 0 ∈ N such that, for all n ≥N ′ 0 and for T ∈ AD(r ǫ n ), Furthermore, by the definition of the χ 2 -divergence in (8) and the fact that P i < 1, i = 1, · · · , k, we have that d χ 2 (T, P ) ≥ T − P 2 2 . Since r ǫ n = O 1 n , these two inequalities imply that T − P 2 = O 1 √ n and, hence, also T − Then we have the following lemma.
Then, using (59), we have for T ∈Ā χ 2 r ǫ n η − n . Hence, we obtain T ∈ AD(r ǫ n ). From (57), we obtain that, for T ∈ AD(r ǫ n ), there exist M ′ 2 > 0 and N ′ 2 ∈ N such that If r < r ǫ n , then, by definition of r ǫ n , the type-I error α n (TD) exceeds ǫ. Hence such a threshold violates (11). We thus assume without loss of optimality that r ≥ r ǫ n . In this case, Together with Lemma 7, this implies that, for r ǫ n and n ≥ N ′ 1 , the type-II error of the divergence test TD n (r) can be lower-bounded as β n (TD n (r)) ≥ Q n (D(t Z n P ) < r ǫ n ) = z n : t z n ∈AD(r ǫ n ) By (54),r n vanishes as n → ∞. We can therefore find anÑ ′ such that, for n ≥Ñ ′ , we have that √r n < √ V (P Q) τ with τ defined in (33). To lower-bound the type-II error further, let us consider the following lemma.
V. CONCLUSIONS For the divergence test and for a large class of divergences, we established the first-order and the second-order terms of the type-II error given that the type-I error is upper-bounded by a given value. The divergence test does not require the knowledge of the distribution of the alternative hypothesis, hence, it is suitable for problems where we only have access to the distribution of the null hypothesis. The class of divergences considered in this paper includes well-known divergences such as the KL divergence and the α-divergence, and our divergence test specializes to the Hoeffding test if the chosen divergence is the KL divergence. It is well-known that the Hoeffding test is first-order optimal in the sense that its first-order term is equal to the first-order term of the Neyman-Pearson test. However, our results demonstrate that the second-order term of the divergence test for the class of divergences considered in this paper, and therefore also of the Hoeffding test, is strictly smaller than the second-order term of the optimal Neyman-Pearson test. The question whether there exists a test that only requires the knowledge of the distribution of the null hypothesis and that achieves the second-order term of the Hoeffding test for the some Q's while it achieves a higher second-order term for the remaining Q's is yet to be explored.

APPENDIX A PROOF OF LEMMA 3
For the probability distribution T = (T 1 , · · · , T k ) T and Q ∈ P(Z), the function f Q (T ) D(T P ) is a smooth function of T = (T 1 , · · · , T k−1 ) T given by The partial derivatives of f Q (·) with respect to T i , i = 1, · · · , k − 1 up to third orders are given by Consequently, the second-order Taylor approximation of f Q (T) around P is given by where the remainder term R 3 (T) in Lagrange's form is given by for Π = P + t(T − P), t ∈ (0, 1). The absolute value of this term can be upper-bounded by When T − P 2 → 0, we have that P i − δ < T i < P i + δ for an arbitrary δ > 0 and i = 1, · · · , k − 1. and Using (83), (84), and (79), one can show that, for any T satisfying 0 < T − P 2 < δ, Hence, from (82), where C ′ 1 3! C 1 (k −1) 3 and the last inequality follows since |T i −P i | ≤ T−P 2 for all i = 1, · · · , k −1. and Thus, by (88), (89), and the definition of the χ 2 -divergence, (80) can be written as as T − P 2 → 0.


Automatic Rate and Rate-Integrating Mode-Switchable Axisymmetric Gyroscope

This paper describes a novel rate and rate-integrating mode-switchable axisymmetric gyroscope. A precession angle tracking algorithm is developed to enable the gyro to switch automatically between rate and rate-integrating modes at preset rate points through a digital control system within the gyro. We also propose a vibrating amplitude control method for the rate-integrating mode that directly extracts the angular rate output to ensure switching stability. In rate mode, the bias instability and angle random walk of the gyro reach 0.106°/h and 0.011°/√h, respectively. Additionally, an input range of over ±5000°/s is measured in rate-integrating mode. The scale factor nonlinearity reaches approximately 116 ppm over the full-scale range. The control system implements effective steering control of the gyroscope, with a switching delay of 10 ms from rate mode to rate-integrating mode and 100 ms from rate-integrating to rate mode. The proposed system actualizes a new type of gyroscope with high accuracy and a wide input range, which combines the benefits of rate and rate-integrating modes.

Introduction
Silicon micro gyroscopes are widely used in the automotive, consumer electronics, industrial, and aerospace fields due to their small size, low power consumption, and low cost [1]. Axisymmetric gyroscopes can be categorized as either rate gyroscopes (RGs) or rate-integrating gyroscopes (RIGs), depending on whether they measure angular rates or rotation angles [2][3][4][5][6]. RG operating in force-to-rebalance mode applies a feedback force that cancels the Coriolis force and detects the angular rate by measuring the magnitude of the feedback force. This approach is highly accurate, but the size of the supply voltage limits the range of this mode. RIG mode, also called whole-angle mode, is an operational mode in which the Coriolis force causes the oscillation pattern of the gyro to precess freely, and the input angle is detected by measuring the change in the precession angle of the resonant modes. Theoretically, RIG mode has unlimited bandwidth and range as well as high stability of the angular gain [7][8][9]. However, manufacturing errors can result in damping asymmetry and stiffness asymmetry [10], leading to a minimum rate threshold for RIGs [11,12]. The RIG mode cannot operate at rates below this threshold, but it is suitable for higher rates.
The research focus of RGs and RIGs is on improving accuracy and reducing the threshold, respectively. The accuracy of RGs can be improved by using electrostatic tuning and quadrature closed-loop control to reduce the frequency split (∆f ) and quadrature errors [13][14][15]. However, the magnitude of the feedback force in the circuit is limited by the operating voltage, making it difficult to achieve a large range for RGs;in the abovementioned references, the range does not exceed ±300 • /s. For RIGs, the coarse compensation for the quadrature electrostatic tuning and frequency electrostatic tuning as well as quadrature closed-loop control can only eliminate quadrature errors, so the residual error still has a 2. Theory

Dynamic Equation
The dynamics of axisymmetric gyroscopes are generally in a non-ideal state due to processing and material factors, which cause problems such as damping and stiffness asymmetries. The dynamic equations of axisymmetric gyroscopes can be written as [22]: ..
y are the velocities of x, y, and ..
y are the accelerations of x, y; ∆ω and ω are, respectively, associated with the difference and the average of the resonant frequencies of the two modes; 1/τ is the damping coefficient; ∆(1/τ) is the damping error; θ τ , θ ω are the included angles between the damping principal axis and the stiffness principal axis, respectively; k is the angular gain; n is the order of vibration modes [23]; Ω is the input angular rate; M is the effective mass of the gyro; and F x , F y are the forces in the xand y-directions. From the point of view of vibration equations, the particle motion trajectory of the gyro can be regarded as an ellipse, as shown in Figure 1b.
where x and y are the vibration amplitudes of two resonant modes, x  , y  are the ities of x, y, and x  , y  are the accelerations of x, y; Δω and ω are, respectively, asso with the difference and the average of the resonant frequencies of the two modes the damping coefficient; Δ(1/τ) is the damping error; θτ, θω are the included ang tween the damping principal axis and the stiffness principal axis, respectively; k angular gain; n is the order of vibration modes [23]; Ω is the input angular rate; M effective mass of the gyro; and Fx, Fy are the forces in the x-and y-directions. Fr point of view of vibration equations, the particle motion trajectory of the gyro can garded as an ellipse, as shown in Figure 1b. The vibration amplitudes can be written as: cos cos sin sin sin cos cos sin where a is the amplitude of oscillation, q is the quadrature motion, θ is the prec angle of the ellipse motion of vibration, and φ is the phase of the oscillation. The control forces along the x-and y-axes can be obtained as follows:  The vibration amplitudes can be written as: where a is the amplitude of oscillation, q is the quadrature motion, θ is the precession angle of the ellipse motion of vibration, and ϕ is the phase of the oscillation. The control forces along the xand y-axes can be obtained as follows: in which F a = F ac cos ωt + F as sin ωt F q = F qc cos ωt + F qs sin ωt (4) In (4), F as and F ac are the sine and cosine components of F a , respectively, F qs and F qc are the sine and cosine components of F q , respectively. We can obtain the ellipse parameters by averaging the equations of motion in Equation (1) over a period of vibration [24]: In (5), a, q, θ can be controlled by F as , F qc , F qs , respectively. The motion parameters are all affected by damping asymmetry and stiffness asymmetry.

Control Algorithm
According to Lynch's algorithm [25], the vibration amplitudes can be demodulated by multiplication to obtain the following four vibration pattern parameters: x c = a cos θ cos ϕ + q sin θ sin ϕ x s = a cos θ sin ϕ − q sin θ cos ϕ y c = a sin θ cos ϕ − q cos θ sin ϕ y s = a sin θ sin ϕ + q cos θ cos ϕ (6) The control variables for the energy (E), quadrature (Q), precession angle (θ), and phase reference error (ϕ) can be obtained by applying a combination operation: where θ can be used for angle control, which allows the oscillation pattern to be fixed at the reference value θ 0 in RG operation mode; ϕ is controlled by a phase-locked loop (PLL) to track the resonant frequency of oscillation. The variables E, Q, and θ form energy, quadrature, and angle control loops, respectively, through proportional-integral (PI) controllers, which are used to maintain the vibration amplitude, reduce the quadrature error, and fix the precession angle, respectively. The corresponding outputs of PI controllers F E , F Q , and F θ can be described as: (8) where K P and K I are PI parameters and E 0 is the reference value for the energy loop. The feedback forces F x and F y can then be derived as: The open-loop angle control enables the gyro to operate in RIG mode with a freelymoving oscillation pattern regulated by the Coriolis force, that is F θ = 0. The closed-loop angle control forces the gyro to operate in RG mode by fixing the precession angle, the PI controller output F θ reflects the detected angular rate. Therefore, the control algorithm enables the gyro to operate in RIG or RG mode, and it realizes the switching between the two modes through the angle open-loop/closed-loop control, providing a theoretical basis for the switching system.

Angular Rate Output in RIG Mode
Considering a detected system cannot identify the high or the low rotation rates with angle information, it is necessary to obtain the angular rate information of the RIG mode in real time for switching from RIG to RG mode. In RIG mode, the force applied in the angle control open-loop scheme is used to control the direction of the oscillation [26]. Additionally, the Coriolis force arising from the rotation causes the oscillation pattern to freely precess, indicating that the input angular rate can be detected by measuring the Coriolis force in the angle control loop. As mentioned in [27], the force containing angular rate information can be extracted through the demodulation of vibration amplitudes without considering quadrature errors, but the underlying expressions for this operation have not yet been derived, due to their work mainly being focused on using other forces to compensate for the relevant errors. In RIG mode, in addition to the change in the precession angle (from which we obtain the angle output), the angular rate can be obtained through a demodulation and a combinatorial calculation of amplitudes of two resonance modes. Therefore, we independently derived the theoretical equation for the rate output in RIG mode, and we performed a successful experimental verification of the equation.
In the absence of the quadrature error and input angular rates, Equation (2) can be described as: x = a cos θ 0 cos ωt The Coriolis forces of the xand y-axes F cx , F cy in RIG mode can be written as: where F c is the amplitude of the Coriolis force, which is proportional to the input angular rate Ω. When Ω is input, the forces from (11) produce the phase-shifted displacements after passing through the gyro, and the displacement changes accordingly: where k fd is the force-to-displacement coefficient. Equation (12) indicates that Ω influences the amplitude of two resonate modes, and it has a considerable effect at larger rates. The detected angular rate V Ω can be calculated as: (13) where k e is the gain of the hardware circuit. The rate signal V Ω is proportional to Ω, and it is related to M, k fd , and k e . Based on Equation (13), the calculation block diagram can be described in Figure 2.
Considering a detected system cannot identify the high or the low rotation rates with angle information, it is necessary to obtain the angular rate information of the RIG mode in real time for switching from RIG to RG mode. In RIG mode, the force applied in the angle control open-loop scheme is used to control the direction of the oscillation [26]. Additionally, the Coriolis force arising from the rotation causes the oscillation pattern to freely precess, indicating that the input angular rate can be detected by measuring the Coriolis force in the angle control loop. As mentioned in [27], the force containing angular rate information can be extracted through the demodulation of vibration amplitudes without considering quadrature errors, but the underlying expressions for this operation have not yet been derived, due to their work mainly being focused on using other forces to compensate for the relevant errors. In RIG mode, in addition to the change in the precession angle (from which we obtain the angle output), the angular rate can be obtained through a demodulation and a combinatorial calculation of amplitudes of two resonance modes. Therefore, we independently derived the theoretical equation for the rate output in RIG mode, and we performed a successful experimental verification of the equation.
In the absence of the quadrature error and input angular rates, Equation (2) can be described as: The Coriolis forces of the x-and y-axes Fcx, Fcy in RIG mode can be written as: where Fc is the amplitude of the Coriolis force, which is proportional to the input angular rate Ω. When Ω is input, the forces from (11) produce the phase-shifted displacements after passing through the gyro, and the displacement changes accordingly: cos cos sin cos sin cos cos cos where kfd is the force-to-displacement coefficient. Equation (12) indicates that Ω influences the amplitude of two resonate modes, and it has a considerable effect at larger rates. The detected angular rate VΩ can be calculated as: where ke is the gain of the hardware circuit. The rate signal VΩ is proportional to Ω, and it is related to M, kfd, and ke. Based on equation (13), the calculation block diagram can be described in Figure 2.   Figure 3a shows the simulation results for the angular rate directly output in RIG mode by a vibrating amplitude control method, where the input angular rates are ±100 • /s, ±200 • /s, and ±300 • /s. The simulation parameters are listed in Table 1. Figure 3b shows that the average output is proportional to the input angular rate and that the scale factor is approximately 2.46 mV/ • /s. The simulation and the measurement results prove the correctness of (13) and provide a basis for improving the switching stability. Figure 3a shows the simulation results for the angular rate directly output in RIG mode by a vibrating amplitude control method, where the input angular rates are ±100°/s, ±200°/s, and ±300°/s. The simulation parameters are listed in Table 1. Figure 3b shows that the average output is proportional to the input angular rate and that the scale factor is approximately 2.46 mV/°/s. The simulation and the measurement results prove the correctness of (13) and provide a basis for improving the switching stability.

Parameters
Value Units fc 1 5045

Switching Control
To ensure that mode switching is carried out reliably and automatically, we choose different preset rate points for the two operational modes, where A and B are the preset RG-to-RIG and RIG-to-RG thresholds, respectively. ΩRG and ΩRIG are the rate outputs of RG and RIG modes, respectively, as shown in Figure 4. We can set switching thresholds as required for the condition A > B, as the system can misjudge or make repeated switches when A ≤ B, resulting in unstable or failed switches. The switch module consists of three functions for automatic comparison, switching, and precession angle tracking ( Figure 4). The gyro operates in RG mode when powered on and switches to RIG mode with open-loop angle control when the output satisfies RG A Ω > . In RIG mode, the system tracks the precession angle information in real time. When RIG B Ω < is satisfied, the information stops refreshing, and it is used as the reference value for the closed-loop angle control to switch back to RG mode.

Parameters
Value Units

Switching Control
To ensure that mode switching is carried out reliably and automatically, we choose different preset rate points for the two operational modes, where A and B are the preset RG-to-RIG and RIG-to-RG thresholds, respectively. Ω RG and Ω RIG are the rate outputs of RG and RIG modes, respectively, as shown in Figure 4. We can set switching thresholds as required for the condition A > B, as the system can misjudge or make repeated switches when A ≤ B, resulting in unstable or failed switches. The switch module consists of three functions for automatic comparison, switching, and precession angle tracking (Figure 4). The gyro operates in RG mode when powered on and switches to RIG mode with open-loop angle control when the output satisfies |Ω RG | > A. In RIG mode, the system tracks the precession angle information in real time. When |Ω RIG | < B is satisfied, the information stops refreshing, and it is used as the reference value for the closed-loop angle control to switch back to RG mode.   Figure 5 shows the simulated rate and the angle outputs of the automatic mode switching. The output signals in Figure 4 can be divided into three stages. In the first stage, the gyro operates in RG mode with a fixed precession angle. As the input rate accelerates in the second stage, the gyro switches from RG mode to RIG mode and the   Figure 4 can be divided into three stages. In the first stage, the gyro operates in RG mode with a fixed precession angle. As the input rate accelerates in the second stage, the gyro switches from RG mode to RIG mode and the precession angle begins to change when the input rate is above 200 • /s. In the third stage, the gyro switches back to RG mode and the precession angle becomes constant again when the input rate is less than 100 • /s. The RG and RIG outputs are expressed in terms of angular rate by dividing the output voltage by the scale factor. The simulation implements mode switching for the gyro operated in RG mode when under low rates and in RIG mode when under high rates.  Figure 5 shows the simulated rate and the angle outputs of the automatic mode switching. The output signals in Figure 4 can be divided into three stages. In the first stage, the gyro operates in RG mode with a fixed precession angle. As the input rate accelerates in the second stage, the gyro switches from RG mode to RIG mode and the precession angle begins to change when the input rate is above 200°/s. In the third stage, the gyro switches back to RG mode and the precession angle becomes constant again when the input rate is less than 100°/s. The RG and RIG outputs are expressed in terms of angular rate by dividing the output voltage by the scale factor. The simulation implements mode switching for the gyro operated in RG mode when under low rates and in RIG mode when under high rates.

Device and Test Platform
A cobweb-like disk resonator gyroscope (CDRG) [28] was used in the experiments. Figure 6a,b shows the structure and the mode shapes of the CDRG. The CDRG is composed of 10 concentric polygonal spider web rings connected by eight spokes to a single central anchor. Each ring is connected by 16 identical rectangular beams. Table 2 presents the main structural parameters of the CDRG. Figure 6c shows the test platform and control circuit. Figure 7 shows a block diagram of the mode switching control system.

Device and Test Platform
A cobweb-like disk resonator gyroscope (CDRG) [28] was used in the experiments. Figure 6a,b shows the structure and the mode shapes of the CDRG. The CDRG is composed of 10 concentric polygonal spider web rings connected by eight spokes to a single central anchor. Each ring is connected by 16 identical rectangular beams. Table 2 presents the main structural parameters of the CDRG. Figure 6c shows the test platform and control circuit. Figure 7 shows a block diagram of the mode switching control system.

Scale Factor
The DC voltages VQ and VT are applied in corresponding electrodes to the aniso stiffness compensation for a CDRG. The raw output of PI controller Fθ is used as the rat data in RG mode. The RIG's rate output is measured by (13). . Block diagram of mode switching control system. The 'demod & calculator' module in the FPGA is used to demodulate and calculate the required control signals; the PLL module is used to track the resonant frequency of the gyro; the switch controller module is used to compare the angular rate and the preset rate and send the switching control flag; and the feedback controller and 'mod' modules are used to generate and to modulate the feedback signals to the gyro.

Scale Factor
The DC voltages V Q and V T are applied in corresponding electrodes to the anisostiffness compensation for a CDRG. The raw output of PI controller F θ is used as the rate data in RG mode. The RIG's rate output is measured by (13).
±30°/s, ±50°/s, ±100°/s, ±200°/s in RG mode and ±200°/s, ±300°/s, ±500°/s, ±1000°/s, ±3000°/s, ±5000°/s in RIG mode. Figure 8a,b shows the scale factor fitting curves in RG and RIG modes, respectively. The measured value is 6.52 mV/°/s in RG mode and 250.1 μV/°/s in RIG mode. The scale factor nonlinearity of the RG and RIG modes in their measurement ranges reaches approximately 116 ppm and 56 ppm, respectively, and the scale factor stability is 334 ppm and 102 ppm, respectively. The RIG operation mode achieves better scale factor performance at larger rates.

Zero-Rate Output Performance in RG Mode
After the open-loop mode-matching operation, zero-rate output (ZRO) data from the FPGA were collected through the serial port at room temperature over an acquisition time of more than 30 min and a sampling frequency of 1 Hz. Figure 9 shows the Allan variance curve of the ZRO data. The figure shows that the bias instability (BI) and angle random walk (ARW) are 0.106°/h and 0.011°/√h, respectively.

Zero-Rate Output Performance in RG Mode
After the open-loop mode-matching operation, zero-rate output (ZRO) data from the FPGA were collected through the serial port at room temperature over an acquisition time of more than 30 min and a sampling frequency of 1 Hz. Figure 9 shows the Allan variance curve of the ZRO data. The figure shows that the bias instability (BI) and angle random walk (ARW) are 0.106 • /h and 0.011 • / √ h, respectively.
±3000°/s, ±5000°/s in RIG mode. Figure 8a,b shows the scale factor fitting curves in RG and RIG modes, respectively. The measured value is 6.52 mV/°/s in RG mode and 250.1 μV/°/s in RIG mode. The scale factor nonlinearity of the RG and RIG modes in their measurement ranges reaches approximately 116 ppm and 56 ppm, respectively, and the scale factor stability is 334 ppm and 102 ppm, respectively. The RIG operation mode achieves better scale factor performance at larger rates.

Zero-Rate Output Performance in RG Mode
After the open-loop mode-matching operation, zero-rate output (ZRO) data from the FPGA were collected through the serial port at room temperature over an acquisition time of more than 30 min and a sampling frequency of 1 Hz. Figure 9 shows the Allan variance curve of the ZRO data. The figure shows that the bias instability (BI) and angle random walk (ARW) are 0.106°/h and 0.011°/√h, respectively.

Mode Switching
According to the experimental scale factor results, the maximum nonlinearity of the RG's scale factor occurs at 200 • /s. Thus, we set the thresholds to A = 150 • /s and B = 75 • /s based on the A > B condition. Figure 10 shows the experimental results of the mode switching and the corresponding precession angle output. In Figure 10a, the red lines represent the rate output in RG mode, whereas the blue lines represent the measured rate in RIG mode. The gyro starts in RG mode and the vibration pattern is fixed at 0 • . When the input angular rate accelerates beyond 150 • /s, the switching control input flag is set to the high level, the gyro switches to RIG mode, where the vibration pattern starts to precess freely. As the input decelerates below 75 • /s, the flag is reset to the low level, the gyro switches to RG mode, and the vibration pattern becomes fixed again. As the scale factors of the RG and RIG are not equal, their outputs are expressed as the measured voltage divided by their scale factors in degrees per second. The switching error from RG mode to RIG mode is about ±1 • /s with a switching delay of 10 ms, compared with about ±4 • /s and 100 ms from RIG mode to RG mode. The switching delay is the time from when the output of the former mode stops to when the output of the latter mode becomes stable. The switching system has a measurement range of ±0-150 • /s in RG mode and ±75-5000 • /s in RIG mode, although the two operation modes do not work simultaneously.
switching and the corresponding precession angle output. In Figure 10a, the red lines represent the rate output in RG mode, whereas the blue lines represent the measured rate in RIG mode. The gyro starts in RG mode and the vibration pattern is fixed at 0°. When the input angular rate accelerates beyond 150°/s, the switching control input flag is set to the high level, the gyro switches to RIG mode, where the vibration pattern starts to precess freely. As the input decelerates below 75°/s, the flag is reset to the low level, the gyro switches to RG mode, and the vibration pattern becomes fixed again. As the scale factors of the RG and RIG are not equal, their outputs are expressed as the measured voltage divided by their scale factors in degrees per second. The switching error from RG mode to RIG mode is about ±1°/s with a switching delay of 10 ms, compared with about ±4°/s and 100 ms from RIG mode to RG mode. The switching delay is the time from when the output of the former mode stops to when the output of the latter mode becomes stable.
The switching system has a measurement range of ±0-150°/s in RG mode and ±75-5000°/s in RIG mode, although the two operation modes do not work simultaneously. The gyro starts in RG mode and switches to RIG mode after a delay of about 10 ms when the input is greater than 150°/s. It switches back to RG mode after a delay of about 100 ms when the input is less than 75°/s. The angle remains fixed in RG mode and precesses freely in RIG mode. Figure 10a also shows that the control system switches between RG and RIG modes successfully and automatically. The rate output of the RIG (blue lines) still has an angle-dependent bias error during rotation due to residual asymmetric errors as clearly demonstrated in (5) and (11). However, in this paper, the gyro only operates in RIG mode at high rates, so the angle-dependent bias has less impact on the performance of the RIG in this situation.
Because the RG and RIG modes output the angular rate in different ways, there will be breakpoints in the rate data after switching, and the number of breakpoints is related to the switching delay. The RG-to-RIG delay can be further reduced by program opti- Figure 10. Process of mode switching. (a) Angular rate outputs of the gyro and (b) corresponding precession angle. The gyro starts in RG mode and switches to RIG mode after a delay of about 10 ms when the input is greater than 150 • /s. It switches back to RG mode after a delay of about 100 ms when the input is less than 75 • /s. The angle remains fixed in RG mode and precesses freely in RIG mode. Figure 10a also shows that the control system switches between RG and RIG modes successfully and automatically. The rate output of the RIG (blue lines) still has an angledependent bias error during rotation due to residual asymmetric errors as clearly demonstrated in (5) and (11). However, in this paper, the gyro only operates in RIG mode at high rates, so the angle-dependent bias has less impact on the performance of the RIG in this situation.
Because the RG and RIG modes output the angular rate in different ways, there will be breakpoints in the rate data after switching, and the number of breakpoints is related to the switching delay. The RG-to-RIG delay can be further reduced by program optimization in FPGA. The RIG-to-RG delay can be further improved by optimizing the PI parameters. The bandwidth of the CDRG will be presented as the bandwidth of RG mode if the input angular rate is less than the switching threshold A, whereas the bandwidth of the CDRG will be limited by the switching delay if the input angular rate is greater than A. Figure 10b shows that the flag matches the characteristics of angle output of each operating mode. RIG mode can output both an angular rate and an angle measurement, allowing for full-scale angle and rate detection.

Conclusions
With the continuous improvement of axisymmetric MEMS gyroscopes, their applications are advancing towards the navigation level and a wide measurement range. This paper has described a switching control system that takes advantage of the high accuracy of RGs and the wide input range of RIGs. In RIG mode, we present the successful simulating and the experimental verification of the algorithm to extract the output angular rate by demodulating and calculating the vibration amplitudes, thus providing a new method of rate detection without the angle derivation operation. The automatic switching between RG/RIG modes is achieved by presetting the switching rate point and controlling the switching flag. The RIG-to-RG switching error is ±4 • /s and the switching delay is 100 ms. The RG-to-RIG switching error is ±1 • /s and the switching delay is 10 ms. The gyro with the proposed system has a BI of 0.106 • /h, an ARW of 0.011 • / √ h, and an input range of over ±5000 • /s, achieving both a high accuracy and a large range of MEMS axisymmetric gyros.
On this basis, through further optimization of the error model, control program, and self-compensation techniques, it should be possible to improve the output accuracy of the RG and to reduce the threshold of the RIG and the delays of switching, laying the foundation for broadening the application range of silicon MEMS axisymmetric gyros. Moreover, the proposed system and method can be applied to other axially symmetric gyros, such as hemispherical gyros.


text:
A Trilevel Model for Segmentation of the Power Transmission Grid Cyber Network

Network segmentation of a power grid's communication system is one way to make the grid more resilient to cyber attacks. We develop a novel trilevel programming model to optimally segment a grid communication system, taking into account the actions of an information technology (IT) administrator, attacker, and grid operator. The IT administrator is given an allowance to segment existing networks, and the attacker is given a fixed budget to attack the segmented communication system in an attempt to inflict damage on the grid. Finally, the grid operator is allowed to redispatch the grid after the attack in order to minimize damage. The resulting problem is a trilevel interdiction problem, which we solve by leveraging current research in bilevel branch and bound. We demonstrate the benefits of optimal network segmentation through case studies on the 9-bus WSCC system and the 30-bus IEEE system. These examples illustrate that network segmentation can significantly reduce the threat.

E
Security enclaves B Balancing authority entities x e,r Indicates whether or not enclave e ∈ E(S) communicates with relay r y e,f Indicates whether or not enclave e communicates with enclave f q n,e Indicates whether or not enclave e is in entity n t e,n Indicates whether or not enclave e communicates with entity n (in the tier below) Attacker decisions: z e Indicates whether or not the attacker infiltrates enclave e δ r Indicates whether or not relay r is compromised v k Indicates whether or not line k is operational u d Indicates whether or not the load d is online w g Indicates whether or not generator g is operational T HE electric grid, like many infrastructure systems, historically relied on air-gapped and specialized control networks for security. For convenience and cost savings, grid/utility owners are increasingly integrating traditional information technologies (IT) into their control systems. This shift has increased the attack surface for electric grids and has made them more vulnerable to cyber attacks. This was illustrated in the December 2015 cyber attack against the Ukrainian power grid [1]. According to [2], nation-state actors in particular pose a significant threat to cyber-physical systems and also possess the ability to gather the necessary information in advance of attacks. Criminal organizations, terrorists, hackers, and hackivists may also pose a threat to these systems [3].
Network segmentation is one strategy that has been proposed to improve cyber network security [4]. Dividing networks into small segments and restricting communications between segments can limit the scope of attacks and the attackers' ability to pivot within a network. Several governmental and regulatory agencies have provided guidance that motivates segmentation and offer strategies for implementing it [5], [6], [7], [8]. However, as noted in [9], this guidance tends to provide general design principles and architectures, but not specific recommendations or tools for cyber network designers.
Designing segmentation strategies for cyber-physical systems presents a particularly difficult challenge as it requires consideration of two complex interconnected systems. In an IT setting, segmentation decisions can be made based on business functions and the expected consequences of disruptions to these functions (e.g. interruption to payment systems, loss of customer data). In a grid setting, the cyber and physical layers are connected and distributed over wide areas; coordinated attacks create the potential for consequences that are greater than the sum of their parts. A key recommendation in [3] is the need to "evaluate the potential risk of a coordinated cyberattack on geographically distributed targets." In this paper, we explore the use of optimization to identify strategies for segmenting cyber networks. We use an explicit representation of the cyber and grid network layers to capture the interconnected nature of these systems and the resulting impacts of attacks. We also employ trilevel optimization techniques to capture the adaptive nature of and interactions between cyber defenders, attackers, and infrastructure operators. Finally, given the potential capabilities of the attacker, we assume that the attacker has full knowledge of the system when crafting attacks.

A. Literature Review
Several researchers have previously considered applying optimization to cyber security problems that are related to network segmentation. In [10], a heuristic optimization algorithm is used to place the minimum number of intrusion detection systems (IDS) on the nodes of a cyber network while ensuring that no nodes are more than n hops from the nearest detection system. This heuristic approach is also applied by [11] to place trust systems on the smart grid cyber network.
In [12], trust system placement is combined with network segmentation. A mixed-integer linear program (MILP) is used to partition a power system supervisory control and data acquisition (SCADA) network into a collection of enclaves such that a trust node is placed on at least one end of all arcs connecting different enclaves. This ensures that all communications between segments are scanned. The authors also develop a heuristic algorithm for partitioning and placing trust nodes on larger networks.
In [9], simulated annealing is used to identify segmentation strategies for a computer network in an air operations center that supports 40 military missions. An emulation test bed is used to simulate the performance of each segmentation strategy against a "dumb" worm attack that propagates in a manner similar to a disease in a Susceptible-Exposed-Infectious-Recovered (SEIR) model. The quality of each segmentation strategy is based on mission delay and the availability of devices on the network. In [13], this approach is extended to an algorithm that co-evolves the parameters available to the attacker and defender in order to optimize both of their strategies. The attacker still employs a "dumb" worm attack and cannot control how it propagates, but they tune the intensity of the attack to balance the number of components compromised against the probability of detection.
In [14], a heuristic method that uses minimum spanning trees is used to improve performance and scalability for the trust node placement problem on smart grid SCADA systems. Variants of this model are considered in [15] and [16]. These variants account for node centrality and communication latency when deciding where to place trust nodes. In [17], this approach is further expanded to consider both link coverage (the number of links adjacent to a trust node) and path tolerance (the longest path not passing through a trust node) when making placement decisions.
In all of the previously mentioned references that focus on cyber-physical security, only the cyber network is considered, while effects on the underlying physical system it controls are not modeled. There are several examples where segmentation decisions are made based on the underlying physical system managed by the cyber network. In [18], an emulation test bed is used to compare two segmentation strategies developed by subject matter experts (SME's) for SCADA systems in a chemical plant. They demonstrate that in process-flow operations, security and safety gains can be made by separating SCADA devices for in-flow and outflow values. This approach is extended in [19] to include an SME-informed heuristic approach that creates segmentation strategies to separate control of in-flows and out-flows. In [20], network segmentation is used to prevent domino effects in storage tank facilities (e.g. chained explosions of nearby tanks) by an attacker who can attack at most one segment. Belief networks are used to estimate the probability of chained reactions and a graph centrality measure is used to find the segmentation strategy. These works use testbeds to study the effects of network segmentation on the underlying physical system. However, they do not use optimization to make their network segmentation decisions.
Prior research has used bilevel programming to consider the effects on the power grid from an intentional attack, though these model only the physical grid, not the cyber communication system. For example, [21] and [22] formulate and solve a bilevel program that models an attacker with perfect knowledge of the grid and the response of a grid operator. The attacker directly de-energizes power grid components such as transmission lines and generators. The grid operator then responds through generator redispatch and load shed to minimize unmet demand. This attacker-defender model is adapted for a cyber attacker in [23], though still without modeling the cyber network explicitly.
Finally, several authors have employed trilevel programming to preemptively protect the power grid against malicious attacks in [24], [25], [26], and [27]. These works focus on hardening power grid components such as transmission lines, generators, and buses against optimal physical attacks on the grid. However, they do not consider the grid's communication system in modeling these attacks.

B. Contributions and Paper Organization
In this work, we propose segmenting the grid's communication network optimally by explicitly modeling both the cyber network and the physical network. Specifically, we use trilevel optimization to model the interconnected decisions made by an IT administrator segmenting a grid's communication network, a cyber attacker maliciously trying to damage the grid by attacking the segmented network, and a grid operator mitigating damage from the attack. The main contributions of this paper are: • The formulation of a trilevel programming model to identify provably optimal segmentation strategies for a power transmission cyber-physical system which reduce the severity of the worst attack available to a malicious cyber attacker. • The use of an infrastructure model (i.e., a DC optimal power flow (DCOPF) grid model) to inform network segmentation decisions. • A description of how to reformulate the trilevel model as a bilevel program, which can be solved using algorithms from the literature (e.g. bilevel branch-and-bound). • A case study showing network segmentation results for both the 9-bus WSCC and 30-bus IEEE test cases, where segmentation reduces the severity of the worst-case attack by 31% and 56%, respectively. In the remainder of the paper, Sections II and III give a description of the model. Section IV details our solution methodology. Section V presents our case studies on the 9bus and 30-bus IEEE systems. We conclude in Sections VI and VII with ideas for future work and conclusions.

II. CYBER-PHYSICAL NETWORK SEGMENTATION MODEL
We first describe the trilevel network segmentation model at a high level. We give the mathematical formulation in the next section. When we talk about the cyber-physical system, we draw a distinction between entities, the bodies which control the grid, enclaves, the networks that make up an entity's IT infrastructure, and physical components such as the lines, generators, and loads on the physical grid. While our model is more general, in this paper we consider three types of entities: balancing authorities, control centers, and substations.

A. Modeling the Communication Network
For a description of how a power transmission grid communication system network is structured, we refer the reader to [28]. We simplify the communication network to a 3-tier forest where each substation enclave, control center enclave, and balancing authority enclave is represented as a node, and parent-child relationships represent that the parent enclave communicates with the child enclave (i.e., data flows between the parent and child enclave). Every substation enclave is the child of exactly one control center enclave. In turn, every control center enclave is the child of exactly one balancing authority enclave. We show an example of the forest representation of the 9-bus WSCC test system communication network in Fig. 1a. In this example, every entity has only one enclave. The balancing authority enclave communicates with two different control center enclaves. The first of these communicates with three substation enclaves, and the second with six substation enclaves.
The physical power grid is connected to its communication system through relays, which are drawn as the leaves of the tree in Fig. 1a. Individual relays are controlled by exactly one substation enclave and can be used by an attacker to deenergize transmission lines, generators, and loads. In Fig. 1c, each shaded block of grid components is colored to indicate which control center enclave from Fig. 1a ultimately controls it.

B. Modeling Network Segmentation
In the outermost problem of the trilevel model, we model the decisions of an IT administrator segmenting the network in order to minimize the load shed from a worst-case cyber attack. We model network segmentation by allowing the IT administrator to partition any security enclave into two or more security enclaves while still respecting the underlying control structure among the entities. By only subdividing existing assignments, we aim to generate new segmentation strategies that are compatible with the original topology. In the graph representation, segmentation corresponds to expanding the graph by dividing nodes. If a node is divided, the new subnetwork must be assigned to the same entity as the original node. For example, if a substation enclave is divided, both new enclaves are subnetworks of the same substation entity. Edges are added to the expanded graph according to the following rules: 1) A new enclave may only communicate with (i.e., be a child of) a parent enclave if the entities containing the new enclave and the parent enclave communicated in the original graph. As an example, a substation enclave may only communicate with a control center enclave if the corresponding substation and control center originally communicated before segmentation.
(a) The 9-bus system communication network before segmentation: Arcs indicate that the parent component communicates with the child component. The balancing authority, control center, and substation enclaves comprise the communication system, and its control over the relays is shown as the leaves of the tree. (Note that in this example, we include fewer relays than in the 9-bus system used in the case study in Section V-B.) (b) The 9-bus system communication network after segmentation. The balancing authority enclave has been subdivided into two enclaves, each of which communicates with one of the control centers. The contol centers also have two enclaves each, as does Substation 4.  2) After partitioning a substation enclave, that substation's relays must be reassigned so that each relay is controlled by exactly one of the new substation enclaves. The result of this partitioning is that a substation's grid connections are separated by the substation's enclaves, making that substation less vulnerable overall. 3) Every substation enclave must communicate with at least one relay.
Continuing our example, see Fig. 1b for a possible segmentation of the 9-bus network. We see that, after network segmentation, the balancing authority has two enclaves, one for each control center. Both control centers gain a new enclave.
For each control center, the substations it originally controlled are divided between its two new enclaves. Finally, Substation 4 is segmented into two enclaves, each communicating with a separate enclave within Control Center 2. Note that each of the two relays controlled by Substation 4 is assigned to a different enclave. The corresponding diagram of the physical grid is shown in Fig. 1d.

C. Attacker-Defender Model
In the inner two problems of the trilevel model, we model a cyber attacker who finds the highest-load-shed attack on the segmented network given that the grid operator can redispatch after their attack in order to minimize load shed. We assume that attacks begin by infiltrating balancing authority enclaves, then progress by gaining access to the control center enclaves followed by substation enclaves, until the attack reaches the relays, where it can produce an effect on the grid. At each step, the attacker can never infiltrate a node in the graph without first infiltrating its parent. We assume the attacker has an upper limit on the number of enclaves that can be compromised, and every time the attacker gains access to a security enclave, exactly one unit of this budget is consumed. In other words, one unit of budget is used for every node in the graph the attacker reaches. When the attacker gains access to a substation enclave, we assume that the attacker will use all relays controlled by that substation enclave to de-energize every accessible grid component. To gain some intuition about the benefits of network segmentation, note for example that originally in Fig. 1, an attacker would need to access three enclaves to gain control of the two relays controlled by Substation 4. After segmentation, the attacker would need to access five enclaves to gain control of the same relays.
Finally, given an attack, the grid operator solves a DCOPF on the remaining components to minimize total load shed.

III. TRILEVEL FORMULATION
We now give a mathematical formulation of the problem described in Section II. Because many constraints involving (parent, child) entity pairs hold at multiple levels -both (balancing authority, control center) and (control center, substation) -we use the set T to denote the set of all three communication network entities and the set Z = {(B, C), (C, S)} to denote adjacent tiers of entities. The trilevel model is: where D is the feasible region of the network designer, A(x, y) is the feasible region of the attacker, given the design decisions x and y, and O(u, v, w) is the DCOPF feasible region given the attack defined by u, v, and w. We model D with the following constraints: q s,e ≤ r∈Rs x e,r , ∀s ∈ S, e ∈ E 1 (S) q s,e ≥ x e,r , ∀s ∈ S, r ∈ R s , e ∈ E 1 (S) Q s,e ≤ r∈Rs x e,r , ∀s ∈ S, e ∈ E 0 (S) Q s,e ≥ x e,r , ∀s ∈ S, r ∈ R s , e ∈ E 0 (S)

e∈E(A)
y e,f = 1, n∈T q n,e = 1, ∀T ∈ T , e ∈ E 1 (T ) y e,f q n,f , t e,n ≥ y e,f Q n,f , t e,n ≥ y e,f q n,f , ∀(A, B) ∈ Z, e ∈ E(A), t e,n ≤ Q m,e , ∀(A, B) ∈ Z, e ∈ E 0 (A), t e,n ≤ q m,e , ∀(A, B) ∈ Z, e ∈ E 1 (A), y e,f ∈ {0, 1}, ∀(e, f ) ∈ (E(C) × E(S)) ∪ (E(B) × E(C)) (15) x e,r ∈ {0, 1}, e ∈ E(S), ∀r ∈ R t e,n ∈ {0, 1}, ∀(e, n) ∈ (E(C) × S) ∪ (E(B) × C) (18) In (2), we require that every substation enclave controls at least one relay, and in (3) we enforce that every relay is controlled by exactly one substation enclave. Referring back to the graph in Fig. 1, the x r,e variables represent whether or not there is an arc from enclave e in the substation tier to relay r. Constraints (4)-(7) enforce that an enclave e is assigned to a substation if and only if some relay controlled by that substation is controlled via e.
In (8), we require that, between any two adjacent tiers, every enclave in the tier below is controlled by exactly one enclave in the tier above. Since y e,f determines whether or not there is an arc from enclave e to enclave f , these constraints enforce that substation enclaves and control center enclaves have exactly one parent.
Constraints (9) guarantee that, for each of the tiers of the communication network, every new enclave is assigned to exactly one entity (substation, control center, or balancing authority depending on the tier). Constraints (10)-(12) enforce that enclave e communicates with entity n in the layer below e if and only if entity n has an enclave f that communicates with e. Returning to the graphs in Fig. 1, t e,n is an indicator of whether enclave e is a parent of some enclave assigned to entity n. Constraints (13)- (14) enforce that if enclave e communicates with entity n in the network layer below it, then e must be assigned to an entity m which controls n. Note that B m is an instance of the set T n with T = B and n = m. Finally, (15)- (18) give variable domains.
Note that (10) and (12) contain bilinear terms. Since all variables in these products are binaries, they can be reformulated linearly with the introduction of one new binary variable per bilinear term. That is, we introduce β e,f,n ∈ {0, 1} for all (A, B) ∈ Z, e ∈ E(A), f ∈ E 1 (B), n ∈ B and the constraints β e,f,n ≤ y e,f , ∀(A, B) ∈ Z, e ∈ E(A), β e,f,n ≤ q n,f , ∀(A, B) ∈ Z, e ∈ E(A), β e,f,n ≥ y e,f + q n,f − 1, ∀(A, B) ∈ Z, e ∈ E(A), and replace the y e,f q n,f terms in (10) and (12) with β e,f,n For fixed network design decisions x and y, A(x, y) is defined by: x e,r z e , ∀r ∈ R Constraint (22) enforces the attacker's budget. Constraint (23) requires that the attacker can only access enclaves which are controlled by already-accessed enclaves. That is, with respect to the graphs in Fig. 1, the attack must begin at a balancing authority enclave, and no child can be accessed if its parent was not. Constraint (24) enforces that a relay is compromised when the enclave that controls it has been accessed by the attacker. Constraints (23) and (24) contain bilinear terms that can be linearized using the same technique that linearizes (10) and (12).

Constraints (25)-(26) enforce that the line k is opened if and only if some relay that controls it is compromised. Constraints (27)-(30) do the same for generators and loads.
Finally, given the attack decision, we define O(u, v, w), the operator's DCOPF feasible region: Equation (31) enforces flow balance at each substation, while the line power flow approximation is enforced in (32). We linearize (32) by replacing it with: Constraints (33) and (34) enforce thermal limits and maximum generation capacity, forcing flow or generation to 0 if the component is compromised. Note that we can drop (33) after adding (39). Also note that in (34), we assume the minimum generation capacity for all generators is 0. While this is not accurate, it is necessary in order to make (31)-(36) feasible for all possible values of u, v, and w, the importance of which we describe in Section IV. In (35), we upper bound load shed by the total load and require that we shed all compromised loads. Bounds on the phase angles are enforced in (36).

IV. SOLUTION METHODOLOGY
Because the grid operator's problem is linear and has an objective opposite the attacker's objective, we can take its dual and linearize the bilinearities that appear in the dual objective function in order to reformulate the inner two levels as one mixed integer linear program (MILP). For details on this procedure, see [22]. Note that the linearization requires upper bounds on the dual variables. The authors of [29] show that verifying the correctness of these dual bounds is as hard as solving the original bilevel problem. Thus, we adopt the heuristic commonly used in prior literature, which is to assume that the duals are bounded above by the power capacity of the component corresponding to their index [21]. After this reformulation, where we take the dual of the linear program which minimizes d∈L l d over the set {(θ, f, p, l) : (31), (34)-(36), (37)-(39)} and combine the resulting maximization problem with the attacker problem, we have a bilevel problem with integer leader and mixed integer follower: where X = (δ, z, u, v, w, f, p, d, θ, µ, ξ + , ξ − , and A D (x, y) is defined by: where µ is the dual of the balance constraint (31), ξ − and ξ + are the duals of (37) and (38) respectively, λ + and λ − are the duals of (33), γ is the dual of the upper bound in (34), α + and α − are the duals of (35), and β + and β − are the duals of (36 where we are assuming that α + d is bounded above by the total demand at load d. Note that for every outer-problem solution (x, y, q, t, β) ∈ D, problem (40) has a finite objective value. That is, A D (x, y) = ∅. This is because, regardless of the network segmentation and the attack, it is always possible to prevent infeasibility by shedding load or turning off generators. We refer to this property as relatively complete recourse. Note that it is because the generator dispatch lower bounds are 0 that we have this property.
Problem (40) is a mixed-integer bilevel model, which we solve using the bilevel branch and bound algorithm from [30]. In [31], [32], and [30], callbacks within the IBM CPLEX solver [33] are used to create a bilevel branch-and-cut solver. The authors make their software available for academic use at [34]. To easily generate the required MPS file for the bilevel branch and bound solver, we formulated our model with Pyomo [35] and [36].

A. An Equivalent Solution Methodology
We found that, when solving (40) with the solver from [34], the solver could not close the gap within a time limit of 72 hours. We suspect this may be a symptom of numerical issues in the cuts. To overcome this, we solved a slightly modified problem, in which the variable β and the constraints it appears in are moved to the inner problem of (40). Formally, the formulation is as follows: Note that, in (49), we have lost the relatively complete recourse property: Fixing the outer level problem's choice of y and q, it is possible to choose β in order to make one or more of (19)- (21) infeasible. However, because the solver from [34] uses the high-point relaxation of (49) to get a lower bound, it assumes the outer problem will not choose to make the inner problem infeasible. That is, we are actually solving the problem: where Essentially, there is an implicit constraint added to the outer problem that its solution is in the projection of the inner problem's feasible region onto the outer problem variables. We therefore show that (50) is equivalent to (40). We first show that the optimal objective value of (50) is a lower bound to that of (40). Note first that P equals the projection of D onto the (x, y, q, t) variables. This follows from the relatively complete recourse property of (40). We can rewrite P as Next, note that the projection of A D (x, y, q) onto the space of the X variables is a restriction of A D (x, y) (since constraints (19)- (21) have been moved to the inner problem). Since, in projected space, we have an equivalent outer problem region and a restricted inner problem region, we have a lower bound. Next, we show that this bound is achieved. To see this, consider a solution (x, y, q, t, β, X) to (50). Then, by our argument above, (x, y, q, t, β) ∈ D, and by the definition of A D (x, y, q), X ∈ A D (x, y). Thus, this solution is feasible in (40). This completes the argument: When using the branch and bound solver from [34] to solve (49), we are solving (40) exactly.

V. RESULTS
We present a case study on the 9-bus WSCC and 30-bus IEEE systems [37] using the methodology given in Section IV-A. Note that, as we have formulated the network segmentation problem, the network designer is forced to use their entire budget, regardless of if every enclave in the budget is necessary to reduce that load shed the attacker can cause. While it would be possible to reformulate the problem to avoid this, it would be necessary to use the network designer's objective to motivate them to be frugal with their budget. This is no longer an interdiction problem, and is in essence a multiobjective variant. For simplicity, we leave this as future work, and instead suggest solving the proposed model for a variety of network designer budgets in order to find a minimal one which achieves a satisfactory reduction in load shed. We will show an example of this process for the 9-bus system.

A. Hardware and Software Specification
We ran our experiments on a Linux server with two Intel® Xeon® Silver 4210 CPU @ 2.20GHz processors and 256 GB RAM. For the CPLEX-based solver from [34], we used CPLEX 12.9.0 with two threads.

B. 9-bus WSCC System
To create intuitive results that are easy to validate, we first use the 3-generator 9-bus WSCC test system [37]. We use a simple communication network consisting of a single balancing authority and two control centers. The three substations with generators are assigned to one control center and the remaining six substations are assigned to the other control center. There is exactly one security enclave at each substation and each grid component at that substation is controlled through a relay hosted by its substation's security enclave. See Fig. 2a for a depiction of this cyber-physical system before segmentation. Note that this starting topology does not correspond exactly to the communication network we showed in Fig. 1. In Fig. 1, we reduced the number of relays in the system for the sake of readability.
We do two case studies on this network: The first defends against an attacker with a budget of 5 and the second defends against an attacker with a budget of 8. We begin with the case where the budget is 5. Before network segmentation, the attacker can cause a complete blackout by pivoting from (a) The 9-bus system before segmentation. The worst-case attack with budget 5 is marked with a skull-and-crossbones on each enclave and component which is compromised. The attacker infiltrates the Balancing Authority enclave, the Control Center 2 enclave, and the three substation enclaves whose relays control the loads and lines at Buses 5, 6, and 8.
(b) The 9-bus system after segmentation, where the network designer was allowed a budget of two control center enclaves and the attacker had a budget of five enclaves. After segmentation, the attacker spends three units of budgets to infiltrate the balancing authority enclave and both the Control Center 2 enclaves. They only have two units of budget remaining, so they infiltrate buses 4 and 7, which islands the two largest generators, resulting in 225 MW of load shed. Note that the attack on the second yellow Control Center 2 enclave does not cause any load shed; it only allows the attacker to use the full attack budget.  the balancing authority, through the second control center, to Substations 5, 6, and 8. These substations contain all of the system's load of 315 MW. See Fig. 2a for a visualization of this attack. Note that the attacker has multiple optimal solutions: For example, they could alternatively pivot from the balancing authority to Control Center 1, and then compromise Substations 1, 2, and 3. This allows them to shut down all the generators, also resulting in 315 MW of load shed.
In Table I, we present results on a partial parameter sweep of network designer budgets for this problem. We see that two enclaves at either the balancing authority or control center level are necessary to prevent any load shed. Furthermore, two enclaves at the balancing authority level prevent 35 MW more load shed than two enclaves at the control center level. Figure 2b displays the optimal network segmentation for a budget of two control center enclaves. The model segments both the control center networks into two enclaves, and reallocates the substations so that it is no longer possible to compromise all three generators or all three loads. After network segmentation, one worst-case attack infiltrates both enclaves in Control Center 2. The attacker can then only access two substations. They choose Substations 4 and 7 to cause a load shed of 225 MW, reducing the load shed possible in the unsegmented network by 28.6%. Note that the attack on the solid yellow right enclave of Control Center 2 is unnecessary. However, due to the segmentation, there is nothing productive for the attacker to do with this last unit of budget.
In the second case study, we decide an optimal network segmentation to protect against an attacker with a budget of 8 enclaves. Clearly, since the attacker with budget 5 was able to shed all the load (see Figure 2a), in the unsegmented case with budget 8, the attacker sheds 315 MW of load also. Using a process similar to what we showed in Table I, we select a network designer budget of two balancing authority enclaves and four control center enclaves. The resulting segmentation is visualized in Fig. 3. The balancing authority is divided into three enclaves, each of which communicates with two control center enclaves. No load bus or generator bus is on the same control center enclave, so, to remove all the loads or all the generators, the attacker would need to spend three units of budget at the substation level and three at the control center level. However, this leaves them only two at the balancing authority level. Because each balancing authority enclave is mapped to one Control Center 1 enclave and one Control Center 2 enclave, they would need nine or more units of budget to attack enough enclaves to cause a complete blackout. Thus, they cannot shed all the load, and 90 MW (the smallest load) is served. Again, network segmentation reduced the load shed by 28.6%. Note in the attack shown in Fig. 3, the attacker wastes two units of budget: Similarly to the budget-5 case, the attacks on the blue balancing authority enclave and the purple control center enclaves are unproductive, but they cannot be reassigned so that the attack increases the load shed.
Both of these examples illustrate that optimal network segmentation isolates the most critical infrastructure in a way that drastically reduces the damage caused by the worst-case attack.

C. 30-bus System
The 30-bus IEEE system provides a larger example with a more complex topology. For this system, we use the voltage levels to create a fictitious communication network. All highvoltage grid components communicate with a single control center and balancing authority. The lower-voltage grid components are split between two control centers. Both these control centers communicate with a second balancing authority. This network is depicted in Fig. 4a. Balancing Authority 1 and Control Center 1 control the high-voltage grid components, and the remaining two control centers and Balancing Authority 2 control the low-voltage components.
In this case study, we consider an attack budget of six enclaves. We again search the results for a range of network designer budgets, trying to find a minimal one with satisfactory reduction of load shed. We present results where we allow two extra balancing authority enclaves, one extra control center enclave, and one extra substation enclave. Our model segments both balancing authorities into two security enclaves each and divides Control Center 1 into two enclaves. Because substations controlled by Control Center 1 have the largest loads and generation, the model assigns substations to the two enclaves in Control Center 1 so that the largest loads communicate with only one control center enclave. Control Centers 2 and 3 are each assigned to their own enclave in Balancing Authority 2. Finally, the substation enclave is used (a) The 30-bus system before segmentation: The network is controlled by two balancing authorities and 3 control centers, each with only one security enclave. The worst-case attack with budget 6 compromises Balancing Authority 1, Control Center 1, and four substation enclaves controlling buses Blaine, Claytor, Kumis, and Reusens.
(b) The 30-bus system after segmentation: Balancing Authority 1 and Control Center 1 have two enclaves each, meaning the attacker would require 4 units of budget to gain access to the same substations. In lieu of this, they attack only the green control center enclave, compromising Blaine, Fieldale, and Reusens. While it is possible for them to also attack Glen Lyn or Cloverdl, the resulting load shed is not higher. Before network segmentation, the attacker can cause a load shed of 119.2 MW, or 63.1% of the total load, by pivoting from the first balancing authority, through its control center, to substations Kumis, Claytor, Blaine, and Reusens. See Fig. 4a. After network segmentation, the worst-case attack uses two units of budget to give the attacker access to the substations controlled by the green control center enclave. This allows them to attack substations Blaine, Fieldale, and Reusens, causing a total load shed of 52.8 MW, or 27.9% of the total demand. This means that network segmentation reduced the worst-case load shed by 55.7%. See Fig 4b for a depiction. Note that the attack on the blue balancing authority enclave is wasted. In fact, the attack on Fieldale is also not necessary: All of the load shed stems from directly compromising the loads at Blaine and Reusens, and the generators at Glen Lyn and Fieldale are not necessary for serving the remaining load.
This case study illustrates an intuitive network segmentation strategy. The largest loads are connected to substations Blaine and Reusens, so an attacker will expend all resources to infiltrate those substations. However, this network segmentation strategy ensures that the attacker will not be able to shed any further load. While it may seem intuitive to put Blaine and Reusens on separate control center enclaves, note that the attacker still has enough budget to access both, expending four units of budget to attack both Balancing Authority 1 enclaves and both Control Center 1 enclaves, then using the remaining two units to attack Blaine and Reusens. The segmentation therefore limits the attacker's ability to compromise additional enclaves.
Finally, we note that the segmentation of Balancing Authority 2 and of the Hancock substation are indeed necessary in this example. Without segmenting Balancing Authority 2, there is an attack on 6 enclaves which sheds 104.7 MW of load by compromising the Balancing Authority 2 enclave, the enclaves in Control Centers 1 and 2, the Roanoke substation with the 6 MW load, and the Cloverdl substation under Control Center 3. By segmenting Balancing Authority 2, this attack is prevented since it takes an additional unit of budget to gain access to both control centers. With only this segmentation and not the extra substation enclave at Hancock, there still exists an attack which shed 53 MW of load, targeting the yellow Balancing Authority 2 enclave, Control Center 2, the Hancock substation with the 11 MW load, Bus 15, Bus 17, and Bus 20. Including the new Hancock substation enclave, the worst attack on the low-voltage side of the network is reduced to 51.7 MW, compromising the red Balancing Authority 2 enclave, Control Center 3, the Roanoke substation with the 6 MW load, Bus 23, and the accessible Cloverdl substation.
This case study shows that even moderate budgets of additional enclaves can greatly decrease the possible load shed achievable by a cyber attacker, in this case reducing it by more than half.
VI. FUTURE WORK The network segmentation model we develop can be improved and extended in several ways. More details from the communication network presented in [28] can be included to result in a more realistic network. The attacker model can be made more realistic by adding more complex movement throughout the network, using enclave access fees based on real data, and adding context to the attacker's budget. For example, an IT administrator could perform analysis to determine the maximum number of enclaves that could likely be comprised before the IDS systems on the network would generate an alert, and use that as the basis for making segmentation decisions. We use the simplest grid operator model available. Future work can focus on using a higher-fidelity power flow model and allowing grid operator transmission switching. In addition, solution strategies that would allow solving the model on realistically-sized networks would be of interest.
VII. CONCLUSION Segmentation is a cyber defense strategy that has been proposed for improving network security. However, specific strategies for optimally implementing these recommendations, while accounting for the underlying system being protected, are lacking. To address this, we develop a trilevel cyber-physical power-transmission-system networksegmentation model that is the first of its kind. It models an IT administrator who must decide how to preemptively segment the power transmission system communication network in preparation for an attacker with perfect information of the whole system. We use an attacker model where a malicious actor must move from balancing authorities through control centers to substations in order to disable grid components through relays. The attacker anticipates how the grid operator will use a DCOPF to redispatch the generators after the attack is executed. We reformulate our model to a mixed-integer bilevel problem and use bilevel branch-and-cut to solve it. Finally, we present network segmentation results on the 9bus WSCC and 30-bus IEEE test systems, demonstrating the benefit of optimal network segmentation in terms of reducing the worst-case load shed the attacker can achieve.

VIII. ACKNOWLEDGMENTS
We would like to thank Santanu Dey for his help in showing the correctness of the methodology in Section IV-A.