Let A ∈ R m × n ( m ≪ n ) and b ∈ R m be a sensing matrix and measurement vector. Consider the following optimization problem: min x | | x | | 0 subject to b = A x . (1)

Notice that   ||x||₀ equals the number of non-zero entries of x. Lemma 1 provides a sufficient condition for obtaining a unique solution to Eq. 1.

Lemma 1. (Candes et al., 2006) If the sparsest solution to Eq. 1 has  ||x||₀ = q and m ≥ 2q and all subsets of 2q columns of A are full rank, then the solution is unique.

Consider a full rank coding matrix C ∈ R ℓ × n ( ℓ > n ) and a set of attacked measurements y = Cx + e where e is an arbitrary and unknown sparse error vector. The classical error correction problem aims at recovering the vector x ∈ R n from the attacked measurements y. To achieve this goal, we need to compute the error vector e since given y = Cx + e and e, we can compute Cx. Then, x can be computed using the full rank matrix C (Candes et al., 2006). Candes et al. (Candes et al., 2006) construct a matrix F such that FCx = 0 for all x. Then, by applying F to y, they obtain Y ̄ = F ( C x + e ) = F e . (2)

Hence, the classical error correction techniques can be transformed into reconstructing a sparse error vector e from the measurement vector Y ̄ = F e . According to Lemma 1, any error vector e that satisfies ‖e‖₀ ≤ q, can be reconstructed if all the subsets of 2q columns of matrix F have full rank.

Consider a discrete-time linear system as follows: x [ k + 1 ] = A x [ k ] , y [ k ] = C x [ k ] + e [ k ] , (3) where x [ k ] ∈ R n and y [ k ] ∈ R p denote the states and outputs of the system at time slot k, respectively. e [ k ] ∈ R p denotes attack signals.

Let Y ∈ R p ⋅ K be the collection of the measurements over K time slots, and let E _q,K denote the set of error vectors e [ 0 ] ; … ; e [ K − 1 ] ∈ R p ⋅ K where each error vector e[k] satisfies ‖e[k]‖₀ ≤ q ≤ p. We have: Y ≜ y [ 0 ] y [ 1 ] ⋮ y [ K − 1 ] = C x [ 0 ] + e [ 0 ] C A x [ 0 ] + e [ 1 ] ⋮ C A K − 1 x [ 0 ] + e [ K − 1 ] = C C A ⋮ C A K − 1 x [ 0 ] + E q , K ≜ Φ x [ 0 ] + E q , K , (4) where Φ ∈ R p ⋅ K × n is the K-step observability matrix of the system at hand. Note that (A, C) represents the system dynamics matrix and if the system is observable, then this can be easily achievable. If the system is not observable, in order to guarantee full rank of Φ, we need to add/allocate measurement units to achieve full rank condition. We assume that rank(Φ) = n; otherwise, we cannot determine x [0] even if there was no attack E _q,K = 0. For instance, if the system is not observable, there is no way to infer any attack signal.

By following a two-step procedure Candes et al. (2006) and Hayden et al. (2016) we can now reconstruct the initial state x [0] from y [k]’s, where k = 0, … , K−1 First, we compute the error vector E _q,K , and then solve for x [0]. To compute E _q,K , we pre-multiply Eq. 4 by Q 1 Q 2 ⊤ , where Q 1 Q 2 ∈ R p ⋅ K × p ⋅ K is orthogonal, Q 1 ∈ R p ⋅ K × n , Q 2 ∈ R p ⋅ K × ( p ⋅ K − n ) , and R 1 ∈ R n × n is a rank-n upper triangular matrix by using QR decomposition of F and then we obtain: Q 1 ⊤ Q 2 ⊤ Y = R 1 0 x [ 0 ] + Q 1 ⊤ Q 2 ⊤ E q , K . (5)

Using the second block row, we have: Y ̄ ≜ Q 2 ⊤ Y = Q 2 ⊤ E q , K , (6) where Q 2 ⊤ ∈ R ( p ⋅ K − n ) × p ⋅ K . Using Lemma 1, the above optimization has a unique, s-sparse solution (where s ≤ q ⋅ K) if all subsets of 2s columns (at most 2q ⋅ K columns) of Q 2 ⊤ are full rank. Using this observation, we consider solving the following l ₁-minimization problem: E ̂ q , K = arg min E norm E l 1  subject to  Y ̄ = Q 2 ⊤ E . (7)

Here, we approximate the l ₀-minimization problem with an l ₁-minimization problem to obtain a convex decoder (Candes and Tao, 2005) (Tropp, 2004).

Using the first block row of Eq. 5, we can now compute x [0] as follows: x [ 0 ] = R 1 − 1 Q 1 ⊤ ( Y − E ̂ q , K ) . (8)

The following result provides the conditions under which exists a unique solution to Eq. 8. The proof follows by using Lemma 1 and the fact that the null space of Q 2 ⊤ is equal to the column space of Φ.

Lemma 2. Chang et al. (2018). There exists a unique solution x [0] if all subsets of 2s columns of Q 2 ⊤ are linearly independent and Φ is full column rank.

The role of Lemma 2 guarantees accurate decoding when the assumption holds.

Consider a microgrid with m + n + l buses. m of the m + n + l buses, indexed by N ( S ) = { 1 , … , m } , have synchronous generators connected to them, and n of the microgrid buses, indexed by N ( I ) = { m + 1 , … , m + n } , have inverter-interfaced power supplies attached to them. The remaining buses, indexed by N ( L ) = { m + n + 1 , … , m + n + l } , have electric loads and no generation3. We further assume that the network interconnecting the buses is linear, i.e., the network can be represented by the nodal admittance matrix Y = G + − 1 B where G and B denote the conductance and susceptance matrices, respectively.

We use the standard structure-preserving model to describe the microgrid’s dynamics (Bergen and Hill, 1981). This model that incorporates the dynamics of generators’ rotor angle and response of load power output to frequency deviations, allows us to preserve the structure of the grid, i.e., no load bus elimination is made. We introduce fictitious buses representing the internal generation voltages for synchronous generators and inverter-interfaced power supplies. Each of these buses is connected to either a synchronous generator or inverter-interfaced power supply bus via reactances that account for transient reactances and connecting lines. These reactances can be considered as transmission lines (Sauer and Pai, 1998). Therefore, in the augmented network, the number of buses is 2(m + n) + l. We further denote the set of fictitious buses for synchronous generators and inverter-interfaced power supplies by N f ( S ) and N f ( I ) , respectively. The augmented network can be modeled by graph, G ( N , E ) , where N = N ( S ) ∪ N ( I ) ∪ N ( L ) ∪ N f ( S ) ∪ N f ( I ) ; and where each element in the edge set E corresponds to an electrical line between a pair of buses in the augmented network. We assume that the network topology is fixed and known to the microgrid operator.

Let V _i and θ _i denote the voltage magnitude and phase angle of bus i ∈ N , respectively. We use a structure-preserving model with constant complex voltage behind reactance [34, Sec. 7.9.2], to represent the dynamics of each synchronous generator. For a synchronous generator at fictitious bus i ∈ N f ( S ) , let θ _i denote the angle of its voltage as measured with respect to a synchronous reference rotating at the nominal system electrical frequency ω ₀. Further, let ω _i denote its rotor electrical angular speed, and let P i m be the turbine’s mechanical power. For each synchronous generator i ∈ N f ( S ) , we have θ ̇ i = ω i − ω 0 , (9) M i ω ̇ i = P i m − D i ( ω i − ω 0 ) − ∑ j ∈ N i V i V j | y i j | sin ( θ i − θ j + ϕ i j ) , (10) τ i P i m ̇ = − P i m + P i s − R i ( ω i , meas − ω 0 ) , (11) where y i j = g i j + − 1 b i j with g _ij and b _ij being the elements of the conductance and susceptance matrices, respectively, and ϕ _ij equals arctan (g _ij /b _ij ). N i is the set of neighbors of node i in graph G ( N , E ) , D _i (in s/rad) is the generator damping coefficient, and M _i (in s ²/rad) is the scaled inertia parameter. Further, R _i is the frequency-power speed-droop characteristic constant, τ _i is the generator governor time constant, and P i s denotes the generator’s power set-point. Notice that ω _i,meas is the measured value of state ω _i for all i ∈ N f ( S ) .

The dynamics of inverter-based power supplies can similarly be represented by a constant voltage behind reactance model, augmented to include a frequency-droop controller. For an inverter-based power supply at fictitious bus i ∈ N f ( I ) , let θ _i denote the angle of the bus voltage measured with respect to the nominal system electrical frequency ω ₀, and let P i s denote the generation set-point. For each inverter-based power supply i ∈ N f ( I ) , we have D i θ ̇ i = P i s − ∑ j ∈ N i V i V j | y i j | sin ( θ i − θ j + ϕ i j ) , (12) where D _i (in s/rad) is the speed-droop characteristic slope of the power supply.

In general, the relationship between the real power drawn at load bus i ∈ N ( S ) ∪ N ( I ) ∪ N ( L ) and the bus voltage and frequency is nonlinear. However, for constant voltages and small frequency deviations around the nominal frequency ω ₀, it is reasonable to assume that the real power drawn by the load equals P i d + D i θ i ̇ where D _i > 0 and P i d > 0 are the constant frequency coefficient of load and the nominal load at bus i ∈ N ( S ) ∪ N ( I ) ∪ N ( L ) , respectively. Therefore, for constant voltages and small frequency deviations, the dynamics of the phase angle at load bus i ∈ N ( S ) ∪ N ( I ) ∪ N ( L ) can be modeled by D i θ i ̇ = − P i d − ∑ j ∈ N i V i V j | y i j | sin ( θ i − θ j + ϕ i j ) . (13)

Here, we ignore reactive powers for generators and loads. Notice that reactive power analysis is unnecessary in this study.

Each bus i is equipped with a controller, transceiver, and measurement unit that allow bus i to communicate with the central controller (i.e., the microgrid operator) as well as other bus controllers. As mentioned earlier, fictitious buses represent the internal generation voltages for synchronous generators and inverter-interfaced power supplies. Each inverter-interfaced power supply is equipped with a controller, transceiver, and measurement unit so that the inverter-interfaced power supply at fictitious bus i can measure and communicate its internal voltage phase angle.

The rotor angle of synchronous generators has an electro-mechanic nature, and hence it cannot be measured via electric measurement units. However, we can estimate rotor angle by using other electric parameters of synchronous generators. In the literature, several studies are addressing the problem of estimating the rotor angle of synchronous generators (Tripathy et al., 2010) (Venkatasubramanian and Kavasseri, 2004). Each synchronous generator is equipped with a controller, transceiver, and measurement unit, and that the measurement unit at bus i is using a rotor angle estimator to compute the internal voltage phase angle for the synchronous generator at bus i.

A communication network connects the local controllers and the microgrid operator. This network allows the local controllers to send their feedback information, including rotors’ speeds and voltage phase angles, to the central controller. It also allows the microgrid operator to send different control commands to the local controllers. In this study, we assume that the communication paths from the microgrid operator to the local controllers are secure while other paths are not secure. Note that typically, central controllers are protected against cyber attacks strongly, and that local controllers are more vulnerable to cyber attacks (CIP Standards, 2020). In particular, we consider the following attack types:

• Attack 1: an attack that corrupts the communication paths from the transceivers to the central controller.

We further assume that the set of attacked measurements can change with time. To illustrate different attack types, we refer the reader to Figure 1.

Next, we propose a secure dynamic state estimator for AC microgrids.

We obtain a discrete-time approximation of the microgrid dynamics using the forward Euler discretization scheme. Let us begin with the synchronous generators, and use the same approach for inverter-interfaced power supplies and loads.

Synchronous Generators: For each synchronous generator at fictitious bus i ∈ N f ( S ) , we have θ i [ k + 1 ] = θ i [ k ] + δ ( ω i [ k ] − ω 0 ) , ω i [ k + 1 ] = α i ω i [ k ] + η i P i m [ k ] + β i + ∑ j ∈ N i f i j ( θ i [ k ] , θ j [ k ] ) , P i m [ k + 1 ] = κ i P i m [ k ] + ζ i [ k ] − ν i ( ω i [ k ] + 0,1 e i [ k ] ) , where f _ij (θ _i [k], θ _j [k]) = γ _ij [k] sin (θ _i [k] −θ _j [k] + ϕ _ij ), η i = δ / M i , ν i = R i δ / τ i , κ i = ( 1 − δ / τ i ) , α i = ( M i − δ D i ) / M i , β i = δ D i ω 0 / M i , ζ i [ k ] = δ ( P i s [ k ] + R i ω 0 ) / τ i , γ i j [ k ] = − δ V i [ k ] V j [ k ] | y i j | / M i . (16)

Here we define some auxiliary variables for presentation purposes. There is no meaning to these parameters/variables. We defined the main variables and parameter in system model.

Notice that ω _i,meas [k] equals 0,1 y i [ k ] which is equal to ( ω i [ k ] + 0,1 e i [ k ] ) . Using Eq. 14–15, we have f i j ( θ i [ k ] , θ j [ k ] ) = γ i j [ k ] sin ( θ i [ k ] − θ j [ k ] + ϕ i j ) = γ i j [ k ] sin ( y i [ k ] − e i [ k ] − y j [ k ] + e j [ k ] + ϕ i j ) = γ i j s [ k ] − γ i j s [ k ] e i j c [ k ] − γ i j c [ k ] e i j s [ k ] , (17) where γ i j s [ k ] = γ i j [ k ] sin ( y i [ k ] − y j [ k ] + ϕ i j ) , γ i j c [ k ] = γ i j [ k ] cos ( y i [ k ] − y j [ k ] + ϕ i j ) , e i j c [ k ] ≜ 1 − cos ( e i [ k ] − e j [ k ] ) , e i j s [ k ] ≜ sin ( e i [ k ] − e j [ k ] ) . (18)

Notice that the coefficients of γ i j c [ k ] and γ i j s [ k ] can be calculated using the measurement received from the local controllers.

Now, the state space model of the synchronous generator at fictitious bus i ∈ N f ( S ) can be represented by x i [ k + 1 ] = A i x i [ k ] + u i [ k ] − 0 H i [ k ] ϵ i [ k ] 0 − B i e i [ k ] y i [ k ] = C i x i [ k ] + e i [ k ] (19) where the state vector x i [ k ] = θ i [ k ] , ω i [ k ] , P i m [ k ] ⊤ , and A i = 1 δ 0 0 α i η i 0 − ν i κ i , B i = 0 0 0 0 0 ν i , C i = 1 0 0 0 1 0 , u i [ k ] = − δ ω 0 , β i + ∑ j ∈ N i γ i j s [ k ] , ζ i [ k ] ⊤ , ϵ i [ k ] = e i j 1 c [ k ] , … , e i j n ( i ) c [ k ] , e i j 1 s [ k ] , … , e i j n ( i ) s [ k ] ⊤ , H i [ k ] = γ i j 1 s [ k ] , … , γ i j n ( i ) s [ k ] , γ i j 1 c [ k ] , … , γ i j n ( i ) c [ k ] .

Notice that n(i) denotes the number of neighbors of node i in graph G ( N , E ) , i.e., n ( i ) = | N i | , and that nodes j ₁, … , j _n(i) represent the neighbors of node i in G ( N , E ) . More precisely, we have H i [ k ] ∈ R 1 × 2 n ( i ) and ϵ i [ k ] ∈ R 2 n ( i ) . Next, by using the same approach, we obtain the state space model of inverter-interfaced power supplies and loads.

Inverter-interfaced Power Supplies: The state space model of the inverter-interfaced power supply at fictitious bus i ∈ N f ( I ) can be described by x i [ k + 1 ] = x i [ k ] + u i [ k ] − H i [ k ] ϵ i [ k ] , y i [ k ] = x i [ k ] + e i [ k ] , (20) where x _i [k] = θ _i [k] and u i [ k ] = δ P i s [ k ] / D i + ∑ j ∈ N i γ i j s [ k ] . H _i [k] and ϵ _i [k] are defined above.

Loads: The state space model of the load at bus i ∈ N ( S ) ∪ N ( I ) ∪ N ( L ) can be described by x i [ k + 1 ] = x i [ k ] + u i [ k ] − H i [ k ] ϵ i [ k ] , y i [ k ] = x i [ k ] + e i [ k ] , (21) where x _i [k] = θ _i [k] and u i [ k ] = − δ P i d [ k ] / D i + ∑ j ∈ N i γ i j s [ k ] . H _i [k] and ϵ _i [k] are defined above.

Consider an enlarged system composed of all the system dynamics: X [ k + 1 ] = A X [ k ] + U [ k ] − H [ k ] ϵ [ k ] − B E [ k ] , Y [ k ] = y 1 [ k ] y 2 [ k ] ⋮ y | N | [ k ] = C X [ k ] + E [ k ] , (22) where A ≜ blkdiag { A 1 , … , A m , I m + 2 n + l } ∈ R ( | N | + 2 m ) × ( | N | + 2 m ) , B ≜ blkdiag { B 1 , … , B m , 0 m + 2 n + l } ∈ R ( | N | + 2 m ) × ( | N | + m ) , C ≜ blkdiag { C 1 , … , C m , I m + 2 n + l } ∈ R ( | N | + m ) × ( | N | + 2 m ) , X [ k ] ≜ x 1 [ k ] ⊤ , … , x | N | [ k ] ⊤ ⊤ ∈ R ( | N | + 2 m ) , U [ k ] ≜ u 1 [ k ] ⊤ , … , u | N | [ k ] ⊤ ⊤ ∈ R ( | N | + 2 m ) , ϵ [ k ] ≜ ϵ 1 [ k ] ⊤ , … , ϵ | N | [ k ] ⊤ ⊤ ∈ R 2 ∑ i ∈ N n ( i ) , E [ k ] ≜ e 1 [ k ] ⊤ , … , e | N | [ k ] ⊤ ⊤ ∈ R | N | + m , H [ k ] ≜ blkdiag 0 1 × 2 n ( 1 ) H 1 [ k ] 0 1 × 2 n ( 1 ) , … , 0 1 × 2 n ( m ) H m [ k ] 0 1 × 2 n ( m ) , , H m + 1 [ k ] , … , H | N | [ k ] ∈ R ( | N | + 2 m ) × 2 ∑ i ∈ N n ( i ) .

Consider K time slots k = 0, … , K−1, and collect the measurements corresponding the K time slots in vector Y ̄ as follows: Y ̄ = Y [ 0 ] Y [ 1 ] − C U [ 0 ] Y [ 2 ] − C A U [ 0 ] − C U [ 1 ] ⋮ Y [ K − 1 ] − C ∑ k = 0 K − 2 A K − 2 − k U [ k ] . (23)

Y ̄ ∈ R K ( | N | + m ) can be rewritten in the following form: Y ̄ = Φ X [ 0 ] + Ψ E ̄ , (24) where Ψ = Ψ 1 Ψ 2 with Ψ 1 ∈ R K ( | N | + m ) × K ( | N | + m ) and Ψ 2 ∈ R K ( | N | + m ) × 2 K ∑ i ∈ N n ( i ) , and Ψ 1 = I | N | + m 0 ⋯ ⋯ − C B I | N | + m ⋯ ⋯ − C A B − C B I | N | + m ⋯ ⋮ ⋮ ⋱ ⋮ − C A K − 2 B ⋯ − C B I | N | + m , Ψ 2 = 0 0 ⋯ ⋯ − C H [ 0 ] 0 ⋯ ⋯ − C A H [ 0 ] − C H [ 1 ] 0 ⋯ ⋮ ⋮ ⋱ ⋮ − C A K − 2 H [ 0 ] ⋯ − C H [ K − 2 ] 0 , Φ = C ; C A ; C A 2 ; ⋯ ; C A K − 1 , E ̄ = E [ 0 ] ; ⋯ ; E [ K − 1 ] ; ϵ [ 0 ] ; ⋯ ; ϵ [ K − 1 ] .

The number of columns of the matrix H [k] and the number of rows of the column vector ϵ [k] can be reduced from 2 ∑ i ∈ N n ( i ) to 2 | E | by using the following properties: e i j c [ k ] = 1 − cos ( e i [ k ] − e j [ k ] ) = e j i c [ k ] , e i j s [ k ] = sin ( e i [ k ] − e j [ k ] ) = − e j i s [ k ] . (25)

Notice that we have γ _ij [k] = γ _ji [k] and ϕ _ij = ϕ _ji . Hence, the dimension of Ψ₂ and E ̄ can also be reduced as follows: dim ( Ψ 2 ) → K ( | N | + m ) × 2 K | E | , dim ( E ̄ ) → K | N | + m + 2 | E | . (26)

We now choose Ω ∈ R ( K | N | − 2 | N | ) × K ( | N | + m ) which annihilates Φ, i.e., ΩΦ = 0. We then have: Y ̄ = Ω Y ̄ = Ω Ψ ̃ E ̃ , (27) where Ψ ̃ and E ̃ are the reduced Ψ and E ̄ .

In error correction, accurate decoding is guaranteed when the coding matrix (i.e., Ω Φ ̃ ) satisfies the Restricted Isometric Properties (RIP). Notice that RIP can be obtained by randomly selecting a coding matrix a priori. In (Chang et al., 2018), Theorem 1 provides a sufficient condition for perfect reconstruction of the states against sensor attacks and describes estimator design by using a state feedback controller. However, in the current setting, since there is a limitation to manipulate the coding matrix, we combine our secure estimation scheme with a Kalman Filter (KF) to enhance its performance in practice.

The microgrid operator measures rotor angle from the three synchronous generator buses, and rotor angle and speed from the three fictitious synchronous generator buses. Hence, the operator has access to nine measurements that are subject to cyber attacks. We assume that from t = 1.1 s onwards, the attacker randomly chooses a set of five measurements out of the nine measurements and corrupts them with random signals at each time step. Figure 3 shows the simulation results for the three scenarios: 1) there is no attack (No Attack), 2) the measurements are attacked and they are not protected with any secure estimator (SE), and 3) the measurements are attacked and the microgrid operator uses the proposed SE. Notice that in Figure 3, we only show phase angles and rotor speeds of the three fictitious synchronous generator buses.

In the microgrid without cyber attacks, the rotor speeds converge to 60 Hz after an initial transient period. As mentioned earlier, no secure estimation is used in Scenario 2 while the micorgrid operator uses the proposed secure state estimation to recover the system states in Scenario 3. Therefore, both rotor angles and speeds cannot be estimated correctly in Scenario 2, as shown in Figure 3. However, the operator can perfectly estimate both rotor angles and speeds in Scenario 3. Figure 4 shows the attack signal, secure estimator’s estimated attack signal, and the estimation error. The results show that the proposed dynamic state estimator accurately estimates the system state. Hence, the system’s dynamics are identical to the system without any attacks.

The microgrid operator measures phase angles from the inverter-interfaced power supply buses. Hence, the operator has access to 50 measurements that are subject to cyber attack: 25 correspond to phase angle measurements of the inverter-interfaced power supply buses and 25 correspond to the phase angle measurements of the corresponding fictitious buses. We assume that from t = 1.1 s onwards, the attacker randomly chooses a set of five measurements out of the 50 measurements and corrupts them with random signals at each time step. The simulation results are shown in Figure 5 and Figure 6. The results show that both phase angles and rotor speeds are severely affected when the measurements are attacked and no secure estimator is used by the operator. However, the microgrid operator can perfectly recover the attack signals and restore the system’s normal dynamics as if there was no attack when secure state estimation is used.