Secure State Estimation and Attack Reconstruction in Cyber-Physical Systems: Sliding Mode Observer Approach

Abstract

A cyber-physical system (CPS) is a tight coupling of computational resources, network communication, and physical processes. They are composed of a set of networked components, including sensors, actuators, control processing units, and communication agents that instrument the physical world to make “smarter.” However, cyber components are also the source of new, unprecedented vulnerabilities to malicious attacks. In order to protect a CPS from attacks, three security levels of protection, detection, and identification are considered. In this chapter, we will discuss the identification level, i.e., secure state estimation and attack reconstruction of CPS with corrupted states and measurements. Considering different attack plans that may assault the states, sensors, or both of them, different online attack reconstruction approaches are discussed. Fixed-gain and adaptive-gain finite-time convergent observation algorithms, specifically sliding mode observers, are applied to online reconstruction of sensor and state attacks. Next, the corrupted measurements and states are to be cleaned up online in order to stop the attack propagation to the CPS via the control signal. The proposed methodologies are applied to an electric power network, whose states and sensors are under attack. Simulation results illustrate the efficacy of the proposed observers.

1. Introduction

Cyber-physical systems (CPS) are the integration of the cyber-world of computing and communications with the physical world. In many systems, control of a physical plant is integrated with a wireless communication network, for example, transportation networks, electric power networks, integrated biological systems, industrial automation systems, and economic systems [1, 2]. Since CPSs use open computation and communication platform architectures, they are vulnerable to suffering adversarial physical faults or cyber-attacks. Faults and cyber-attacks are referred to as attacks throughout this chapter.

Recent real-world cyber-attacks, including multiple power blackouts in Brazil [3], and the Stuxnet attack [4] in 2010, showed the importance of providing security to CPSs. Identification and modeling process as [5, 6] which are based on data can be seriously affected by corrupted data. As a result, information security techniques [7] may be not sufficient for protecting systems from sophisticated cyber-attacks. It is suggested in [8] that information security mechanisms have to be complemented by specially designed resilient control systems. Controlling CPS with sensors and actuators, who are hijacked/corrupted remotely or physically by the attackers, is a challenge. The use of novel control/observation algorithms is proposed in this chapter for recovering CPS performance online if an attacker penetrates the information security mechanisms.

Cyber security of CPS must provide three main security goals: availability, confidentiality, and integrity [7]. This means that the CPS is to be accessible and usable upon demand, the information has to be kept secret from unauthorized users, and the trustworthiness of data has to be guaranteed. Lack of availability, confidentiality, and integrity yields denial of service, disclosure, and deception, respectively. A specific kind of deception attack called a replay attack has been investigated when the system model is unknown to the attackers but they have access to the all sensors [9, 10]. Replay attacks are carried out by “hijacking” the sensors, recording the readings for a certain time, and repeating such readings while injecting them together with an exogenous signal into the system’s sensors. It is shown that these attacks can be detected by injecting a random signal, unknown to the attacker, into the system. In the case when the system’s dynamic model is known to the attacker, another kind of deception attack, called a cover attack, has been studied in [11], and the proposed algorithm allows cancelling out the effect of this attack on the system dynamics. In systems with unstable modes, false data injection attacks are applied to make some unstable modes unobservable [12]. Denial of service attacks assaults data availability through blocking information flows between different components of the CPS. The attacker can jam the communication channels, modify devices, and prevent them from sending data, violate the routing protocols, etc. [13]. In a stealth attack, the attacker modifies some sensor readings by physically tampering with the individual meters or by getting access to some communication channels [14, 15]. As a result, detecting and isolating of cyber-attacks in CPSs has received immense attention [16]. However, how to ensure the CPS can continue functioning properly if a cyber-attack has happened is another serious problem that should be investigated; therefore, the focus of this chapter is on resilient control of CPS.

In [17], new adaptive control architectures that can foil malicious sensor and actuator attacks are developed without reconstructing the attacks, by means of feedback control only. A sparse recovery algorithm is applied to reconstruct online the cyber-attacks in [18]. Sliding mode control with advantages of quick response and strong robustness is one of the best approaches to control CPS [19, 20, 21, 22]. In [23], a finite-time convergent higher-order sliding mode (HOSM) observer, based on a HOSM differentiator and a sparse recovery algorithm, are used to reconstruct online the cyber-attack in a nonlinear system. Detection and observation of a scalar attack by a sliding mode observer (SMO) has been accomplished for a linearized differential-algebraic model of an electric power network when plant and sensor attacks do not occur simultaneously [24]. Cyber-attacks against phasor measurement unit (PMU) networks are considered in [25], where a risk mitigation technique determines whether a certain PMU should be kept connected to network or removed. In [26] a sliding mode-based observation algorithm is used to reconstruct the attacks asymptotically. This reconstruction is approximate only, since pseudo-inverse techniques are used.

In this chapter, CPSs controlled by a control input subject to sensor attacks and state/plant attacks are considered. The corrupted measurements propagate the attack signals to the CPS through the control signals causing CPS performance degradation. The main challenge that is addressed in the chapter is online exact reconstruction of the sensor and state attacks with an application to an electric power network. The contribution of this chapter is:

• Novel fixed and adaptive-gain SMO for the linearized/linear CPS under attack are proposed for the online reconstruction of sensor attacks. The time-varying attacks are reconstructed via the proposed SMO that includes a newly designed dynamic filter. Note that the well-known SMO proposed in [27] reconstructs the slow-varying perturbations only.

• A super twisting SMO is applied to reconstruct the state/plant time-varying attacks of the linearized/linear CPS under attack.

• For online state/plant attack reconstruction in nonlinear CPS under attack, a higher-order sliding mode disturbance observer [28] is used.

• An algorithm that use sliding mode differentiation techniques [29] in concert with the finite-time convergent observer for the sparse signal recovery is applied to online reconstruction of time-varying attack in nonlinear CPS under attack when we have limited measurements and more possible sources of attack [30].

2. Motivation example: electric power network under attack

In a real-world power network, only a small group of generator rotor angles and rates is directly measured, and typical attacks aim at injecting disturbance signals that mainly affect the sensorless generators [24].

The small-signal version of the classic structure-preserving power network model is adopted to describe the dynamics of a power network. Consider a connected power network consisting of $n_1$ generators $\left\{g_1,\dots ,g_{n_1}\right\}$ and $n_2$ load buses $\left\{b_{n_1+1},\dots ,b_{n_1+n_2}\right\}$ . The interconnection structure of the power network is encoded by a connected susceptance-weighted graph G. The vertices of G are the generators $g_i$ and the buses $b_i$ . The edges of G are the transmission lines $\left\{b_i,b_j\right\}$ and the connections $\left\{g_i,b_i\right\}$ weighted by their susceptance values. The Laplacian associated with the susceptance-weighted graph is the symmetric susceptance matrix $L\in {\mathbb{R}}^{\left(n_1+n_2\right)\times \left(n_1+n_2\right)}$ defined by $L^{\theta }=\left[\begin{array}{cc}{L}_{g,g}^{\theta }& L_{g,l}^{\theta }\\ L_{l,g}^{\theta }& L_{l,l}^{\theta }\end{array}\right]$ [8].

The CPS that motivates the results presented in this work is the US Western Electricity Coordinating Council (WECC) power system [8] under attack with three generators and six buses, whose electrical schematic is presented in Figure 1. The mathematical model of the power network in Figure 1 under sensor stealth attack and deception attack can be represented as the following descriptor equations that consist of differential and algebraic equations [8]:

Figure 1.

The WECC power system [8].

$$\left[\begin{array}{ccc}I& 0& 0\\ 0& M_g& 0\\ 0& 0& 0\end{array}\right]\left[\begin{array}{c}\dot{\delta }\\ \dot{\omega }\\ \dot{\theta }\end{array}\right]=-\left[\begin{array}{ccc}0& -I& 0\\ L_{g,g}^{\theta }& E_g& L_{g,l}^{\theta }\\ L_{l,g}^{\theta }& 0& L_{l,l}^{\theta }\end{array}\right]\underset{x}{\underbrace{\left[\begin{array}{c}\delta \\ \omega \\ \theta \end{array}\right]}}+\underset{B}{\underbrace{\left[\begin{array}{c}0\\ B_{\omega }\\ B_{\theta }\end{array}\right]}}{d}_x+\left[\begin{array}{c}0\\ P_{\omega }\\ P_{\theta }\end{array}\right],y=\mathit{Cx}+{\mathit{Dd}}_y$$ E1

where the state vector $x={\left[\begin{array}{ccc}{\delta }^T& {\omega }^T& {\theta }^T\end{array}\right]}^T$ includes the vector of rotor angles $\delta \in {\mathbb{R}}^3$ , the vector of generator speed deviations from synchronicity $\omega \in {\mathbb{R}}^3$ , as well as the vector of voltage angles at the buses $\theta \in {\mathbb{R}}^6$ . The $y\in {\mathbb{R}}^p$ is the measurement vector, $d_x\in {\mathbb{R}}^{m_1}$ is the Deception attack corrupting the states, and $d_y\in {\mathbb{R}}^{m-m_1}$ is the stealth attack vector spoofing the measurements. Note that the states of the plant are under attack even if they are not attacked directly but via propagation.

The measurement corruption attacks through an output control feedback. The matrices $E_g,M_g\in {\mathbb{R}}^{3\times 3}$ are diagonal whose nonzero entries consist of the damping coefficients and the normalized inertias of the generators, respectively:

$$M_g=\left[\begin{array}{ccc}0.125& 0& 0\\ 0& 0.034& 0\\ 0& 0& 0.016\end{array}\right],E_g=\left[\begin{array}{ccc}0.125& 0& 0\\ 0& 0.068& 0\\ 0& 0& 0.048\end{array}\right]$$ E2

The inputs $P_{\omega }$ and $P_{\theta }$ are due to known changes in the mechanical input power to the generators and real power demands at the loads. The matrices $B\in {\mathbb{R}}^{12\times m_1}$ and $D\in {\mathbb{R}}^{p\times \left(m-m_1\right)}$ are the attack distribution matrices, and $C\in {\mathbb{R}}^{p\times 12}$ is the output gain matrix. The $L^{\theta }\in {\mathbb{R}}^{9\times 9}$ with $\begin{array}{c}{L}_{g,g}^{\theta }\in {\mathbb{R}}^{3\times 3},L_{g,l}^{\theta }\in {\mathbb{R}}^{3\times 6},L_{l,g}^{\theta }\in {\mathbb{R}}^{6\times 3},L_{l,l}^{\theta }\in {\mathbb{R}}^{6\times 6}\end{array}$ is giving by

$$L^{\theta }=\left[\begin{array}{ccccccccc}0.058& 0& 0& -0.058& 0& 0& 0& 0& 0\\ 0& 0.063& 0& 0& -0.063& 0& 0& 0& 0\\ 0& 0& 0.059& 0& 0& -0.059& 0& 0& 0\\ -0.058& 0& 0& 0.265& 0& 0& -0.085& -0.092& 0\\ 0& -0.063& 0& 0& 0.296& 0& -0.161& 0& -0.072\\ 0& 0& -0.059& 0& 0& 0.330& 0& -0.170& -0.101\\ 0& 0& 0& -0.085& -0.161& 0& 0.246& 0& 0\\ 0& 0& 0& -0.092& 0& -0.170& 0& 0.262& 0\\ 0& 0& 0& 0& -0.072& -0.101& 0& 0& 0.173\end{array}\right]$$ E3

Note that ${\omega }_i\to 0 \forall i=1,2,3$ in a case of the nominal performance of the studied network. Consider the case when the outputs of system, which are the measurement sensors ${\omega }_1,{\omega }_2,{\omega }_3$ , are corrupted by the following stealth attacks.

$$\begin{array}{ccc}{d}_1=-{\omega }_1+2sin\left(\mathit{\pi t}\right),& d_2=-{\omega }_2+cos\left(0.5\mathit{\pi t}\right),& d_3=-{\omega }_3+sin\left(\mathit{\pi t}\right)\end{array}$$ E4

The system (1) was simulated with and without above attacks. Based on the simulation results shown in Figures 2 and 3, the stealth attack in (4) yields inappropriate degradation of the power network performance.

Figure 2.

Comparing corrupted sensor measurements ( ${\omega }_1,{\omega }_2,{\omega }_3$ under attack) and sensor measurements when there is no attack.

Figure 3.

Comparing corrupted states ( ${\delta }_1,{\delta }_2,{\delta }_3$ under attack) and stats when there is no attack.

This motivates why online reconstruction of the attacks followed by cleanup of the measurements prior to using them in control signal is of prime importance for retaining the performance of the power network (as it will be shown in Section VI where the proposed SMO is applied to achieve this goal). The case study of the power network (1) will be further discussed in details in Section 6.

3. Cyber-physical system dynamics

Consider the following completely observable and asymptotically stable system

$$\begin{array}{l}\dot{x}=f\left(x\right)+B\left(x\right)d\left(t\right)\\ y=C\left(x\right)+\mathit{Dd}\left(t\right)\end{array}$$ E5

where $x\in {\mathbb{R}}^n$ is the state vector, $f\left(x\right)\in {\mathbb{R}}^n$ is a smooth vector field, $d\left(t\right)\in {\mathbb{R}}^m$ denotes the attack/fault vector which is additive and matched to the control signal, $y\in {\mathbb{R}}^p$ is the measurement vector, $p\ge m$ , $C\left(x\right)\in {\mathbb{R}}^p$ is the output smooth vector field, $B\left(x\right)\in {\mathbb{R}}^{n\times m}$ and $D\in {\mathbb{R}}^{p\times m}$ denote the attack/fault distribution matrices. For notational convenience, and without affecting generality, the input distribution matrices can be partitioned as

$$B\left(x\right)=\left[\begin{array}{cc}{B}_1\left(x\right)& {\mathbf{0}}_1\end{array}\right],D=\left[\begin{array}{cc}{\mathbf{0}}_2& D_1\end{array}\right]$$ E6

where $B_1\left(x\right)\in {\mathbb{R}}^{n\times m_1},D_1\in {\mathbb{R}}^{p\times \left(m-m_1\right)},{\mathbf{0}}_1\in {\mathbb{R}}^{n\times \left(m-m_1\right)},{\mathbf{0}}_2\in {\mathbb{R}}^{p\times m_1}$ where $m_1\le m$ .

Assumption (A1): $B_1\left(x\right),D_1$ are of full rank.

The attack/fault vector is partitioned accordingly as

$$d=\left[\begin{array}{c}{d}_x\\ d_y\end{array}\right]\begin{array}{ccc}& \mathit{\text{where}}& d_x\in {\mathbb{R}}^{m_1}\begin{array}{cc}\mathit{\text{and}}& d_y\in {\mathbb{R}}^{m-m_1}\end{array}\end{array}$$ E7

Therefore, Eq. (5) can be rewritten as

$$\begin{array}{l}\dot{x}=f\left(x\right)+B_1\left(x\right)d_x\left(t\right)\\ y=C\left(x\right)+D_1{d}_y\left(t\right)\end{array}$$ E8

where $d_x\left(t\right)$ , $d_y\left(t\right)$ represent the state and the sensor attack vectors, respectively. Different attack strategies are shown in Table 1 and discussed in Section 1.

Attack plan	$d_x\left(t\right)\ne 0$	$d_y\left(t\right)\ne 0$	Access to all sensors	Need to know the system model
Stealth attack		√
Deception attack	√
Reply attack	√	√	√
Covert attack	√	√		√
False data injection attack		√		√

Table 1.

Cyber-attack strategies.

Since $p\ge m-m_1$ , the system (8) can be partitioned using a nonsingular transformation $M\in {\mathbb{R}}^{p\times p}$

$$y=M\overline{y}$$ E9

selected so that

$$M^{-1}{D}_1=\left[\begin{array}{c}{\mathbf{0}}_{\left(p-\left(m-m_1\right)\right)\times \left(m-m_1\right)}\\ {\overline{D}}_{1_{\left(m-m_1\right)\times \left(m-m_1\right)}}\end{array}\right]$$ E10

Taking into account (10), system (8) is reduced to

$$\begin{array}{l}\dot{x}=f\left(x\right)+B_1\left(x\right)d_x\left(t\right)\\ {\overline{y}}_1=C_1\left(x\right),{\overline{y}}_2=C_2\left(x\right)+{\overline{D}}_1{d}_y\left(t\right)\end{array}$$ E11

where ${\overline{y}}_1\in {\mathbb{R}}^{p_1}$ with $p_1=p-\left(m-m_1\right)$ and ${\overline{y}}_2\in {\mathbb{R}}^{p_2}$ where $p_2=m-m_1$ . Note that the state attack vector $d_x\left(t\right)$ is additive and matched to the control input that is embedded in system Eq. (11) already.

4. Problem formulation

Assumption (A2): Attacks are detectable, i.e., the invariant zeros of Eq. (11) are stable.

The problem is to protect the closed loop system (11) from the sensor attack $d_y\in {\mathbb{R}}^{m-m_1}$ and state/plant attack $d_x\left(t\right)\in {\mathbb{R}}^{m_1}$ by means of designing fixed-gain and adaptive-gain SMOs that allow: (a) reconstructing online the sensor attack $d_y$ , the state/plant attack $d_x\left(t\right)$ , and the plant states $x$ so that

$$\begin{array}{ccc}{\widehat{d}}_x\left(t\right)\to d_x\left(t\right),& {\widehat{d}}_y\left(t\right)\to d_y\left(t\right),& \widehat{x}\to x\end{array}$$ E12

as time increases and.

(b) “cleanup” of the plant and sensors so that the dynamics of the CPS under attack (11) approaches,

$$\begin{array}{cc}{\dot{x}}_{\mathit{\text{clean}}}=f\left(\widehat{x}\right)+B_1\left(\widehat{x}\right)\left(d_x\left(t\right)-{\widehat{d}}_x\left(t\right)\right),& y_{\mathit{\text{clean}}}=y-D_1{\widehat{d}}_y=C\left(\widehat{x}\right)+D_1\left(d_y\left(t\right)-{\widehat{d}}_y\left(t\right)\right)\end{array}.$$ E13

as time increases, to.

Note that Eq. (13) represents the compensated CPS that converges to CPS without attack as time increases.

5. Results: secure state estimation

In this chapter, for the linearized case of the system in Eq. (5), two SMOs for state estimation and attack reconstruction are discussed. Two other SMO strategies for nonlinear system (5) are also proposed and investigated.

5.1 Attack reconstruction in linear system via filtering by adaptive sliding mode observer

Consider the linearized system in Eq. (5) with $C\left(x\right)=\mathit{Cx}$ and $B\left(x\right)=B$

$$\dot{x}=\mathit{Ax}+\mathit{Bd}\left(t\right),y=\mathit{Cx}+\mathit{Dd}\left(t\right)$$ E14

5.1.1 System’s transformation

Considering system Eq. (14) and assuming assumption (A1) holds, then as show in [29] there exists a matrix $N\in R^{\left(n-p\right)\times n}$ such that the square matrix

$$T_c=\left[\begin{array}{c}N\\ C\end{array}\right]$$ E15

is nonsingular and the change of coordinates $x\mapsto T_{c}x$ creates, without loss of generality, a new state-space representation $\left(A^{\prime },B^{\prime },C^{\prime },D\right)$ where

$$\begin{array}{ccc}{A}^{\prime }=T_c{{\mathit{AT}}_c}^{-1},& B^{\prime }=T_{c}B,& C^{\prime }={{\mathit{CT}}_c}^{-1}=\left[\begin{array}{cc}{\mathbf{0}}_{p\times \left(n-p\right)}& I_{p\times p}\end{array}\right]\end{array}$$ E16

After the linear changing of coordinate, the CPS Eq. (14) is rewritten as

$$\begin{array}{ccc}\begin{array}{l}{\dot{x}}_1=A_{11}{x}_1+A_{12}{x}_2+B_{1}d\\ {\dot{x}}_2=A_{21}{x}_1+A_{22}{x}_2+B_{2}d\\ y=x_2+\mathit{Dd}\end{array}& \mathit{\text{where}}& \begin{array}{cc}{A}^{\prime }=\left[\begin{array}{cc}{A}_{11}& A_{12}\\ A_{21}& A_{22}\end{array}\right],& B^{\prime }=\left[\begin{array}{c}{B}_1\\ B_2\end{array}\right]\end{array}\end{array}$$ E17

with $\begin{array}{cc}{x}_1\in R^{n-p},& x_2\in R^p\end{array}$ , $B_1\in R^{\left(n-p\right)\times m}$ , $B_2\in R^{p\times m}$ , $A_{11}\in R^{\left(n-p\right)\times \left(n-p\right)}$ , $A_{12}\in R^{\left(n-p\right)\times p}$ , $A_{21}\in R^{p\times \left(n-p\right)},A_{22}\in R^{p\times p}$ . It is well known that $\left(A,C\right)$ is observable if and only if $\left(A_{11},A_{21}\right)$ is observable [31].

Defining a further change of coordinates ${\overline{x}}_1=x_1+{\mathit{Lx}}_2$ where $L\in {\mathbb{R}}^{\left(n-p\right)\times p}$ is the design matrix, then the system Eq. (17) can be rewritten as

$$\begin{array}{cc}\begin{array}{l}{{\displaystyle \dot{{\displaystyle \overline{x}}}}}_1={{\displaystyle \tilde{A}}}_{11}{{\displaystyle \overline{x}}}_1+{{\displaystyle \tilde{A}}}_{12}{x}_2+{{\displaystyle \tilde{B}}}_{1}d\\ {{\displaystyle \dot{x}}}_2={{\displaystyle \tilde{A}}}_{21}{{\displaystyle \overline{x}}}_1+{{\displaystyle \tilde{A}}}_{22}{x}_2+{{\displaystyle \tilde{B}}}_{2}d\end{array}& , y=x_2+\mathit{Dd}\end{array}$$ E18

where $\begin{array}{cc}{\tilde{A}}_{11}=A_{11}+{\mathit{LA}}_{21},& {\tilde{A}}_{12}=-A_{11}L+A_{12}-{\mathit{LA}}_{21}L+{\mathit{LA}}_{22}\end{array}$ , ${\tilde{B}}_1=B_1+{\mathit{LB}}_2$ , ${\tilde{A}}_{21}=A_{21}$ , ${\tilde{A}}_{22}=A_{22}-A_{21}L$ , ${\tilde{B}}_2=B_2$ . Since $\left(A_{11},A_{21}\right)$ is observable, there exist choices of the matrix $L$ so that the matrix ${\tilde{A}}_{11}=A_{11}+{\mathit{LA}}_{21}$ is Hurwitz.

Assumption (A3): The attack $d\left(t\right)$ and its derivative are norm bounded, i.e.,

$‖d‖<k_d$ and $‖\dot{d}‖<l_d$ where $k_d,l_d>0$ and are known.

Since $p>m$ , there exists a nonsingular scaling matrix $Q\in R^{p\times p}$ such that

$$\mathit{QD}=\left[\begin{array}{c}{\mathbf{0}}_{\left(p-m\right)\times m}\\ D_2\end{array}\right]$$ E19

where $D_2\in R^{m\times m}$ is nonsingular. Define $\overline{y}$ as the scaling of the measured outputs $y$ according to $\overline{y}=\mathit{Qy}$ . Partition the output of the CPS into unpolluted measurements ${\overline{y}}_1\in {\mathbb{R}}^{p-m}$ and polluted measurements ${\overline{y}}_2\in {\mathbb{R}}^m$ as

$$\overline{y}=\left[\begin{array}{c}{\overline{y}}_1\\ {\overline{y}}_2\end{array}\right]=\left[\begin{array}{c}{Q}_1{x}_2\\ Q_2{x}_2+D_{2}d\end{array}\right]={\mathit{Qx}}_2+\left[\begin{array}{c}{\mathbf{0}}_{\left(p-m\right)\times m}\\ D_2\end{array}\right]d$$ E20

Scale state component $x_2$ and define ${\overline{x}}_2={\mathit{Qx}}_2$ . Then Eq. (18) can be rewritten as

$$\begin{array}{ccc}\begin{array}{l}{{\displaystyle \dot{{\displaystyle \overline{x}}}}}_1={{\displaystyle \overline{A}}}_{11}{{\displaystyle \overline{x}}}_1+{{\displaystyle \overline{A}}}_{12}{{\displaystyle \overline{x}}}_2+{{\displaystyle \overline{B}}}_{1}d\\ {{\displaystyle \dot{{\displaystyle \overline{x}}}}}_2={{\displaystyle \overline{A}}}_{21}{{\displaystyle \overline{x}}}_1+{{\displaystyle \overline{A}}}_{22}{{\displaystyle \overline{x}}}_2+{{\displaystyle \overline{B}}}_{2}d\end{array}& , & {\displaystyle \overline{y}}={{\displaystyle \overline{x}}}_2+\left[\begin{array}{c}\mathbf{0}\\ D_2\end{array}\right]d\end{array}$$ E21

where ${\overline{A}}_{11}={\tilde{A}}_{11}$ , ${\overline{A}}_{12}={\tilde{A}}_{12}{Q}^{-1}$ , ${\overline{B}}_1={\tilde{B}}_1$ , ${\overline{A}}_{21}=Q{\tilde{A}}_{21}$ , ${\overline{A}}_{22}=Q{\tilde{A}}_{22}{Q}^{-1}$ , and ${\overline{B}}_2=Q{\tilde{B}}_2$ . Define ${\overline{x}}_2=\mathit{col}\left({\overline{x}}_{21},{\overline{x}}_{22}\right)$ , where ${\overline{x}}_{21}\in {\mathbb{R}}^{p-m}$ and ${\overline{x}}_{22}\in {\mathbb{R}}^m$ . Consequently the system in Eq. (21) can be written in partitioned form as

$$\begin{array}{l}\begin{array}{l}{\displaystyle \dot{{\displaystyle \overline{x}}}}={\displaystyle \overline{A}}{\displaystyle \overline{x}}+{\displaystyle \overline{B}}d\\ {{\displaystyle \overline{y}}}_1={{\displaystyle \overline{C}}}_1{\displaystyle \overline{x}}, {{\displaystyle \overline{y}}}_2={{\displaystyle \overline{C}}}_2{\displaystyle \overline{x}}+D_{2}d\end{array},{\displaystyle \overline{x}}=\left[\begin{array}{c}{{\displaystyle \overline{x}}}_1\\ {{\displaystyle \overline{x}}}_{21}\\ {{\displaystyle \overline{x}}}_{22}\end{array}\right],{\displaystyle \overline{A}}=\left[\begin{array}{ccc}{{\displaystyle \overline{A}}}_{11}& {{\displaystyle \overline{A}}}_{12a}& {{\displaystyle \overline{A}}}_{12b}\\ {{\displaystyle \overline{A}}}_{21a}& {{\displaystyle \overline{A}}}_{22a}& {{\displaystyle \overline{A}}}_{22b}\\ {{\displaystyle \overline{A}}}_{21b}& {{\displaystyle \overline{A}}}_{22c}& {{\displaystyle \overline{A}}}_{22d}\end{array}\right],{\displaystyle \overline{B}}=\left[\begin{array}{c}{{\displaystyle \overline{B}}}_1\\ {{\displaystyle \overline{B}}}_{21}\\ {{\displaystyle \overline{B}}}_{22}\end{array}\right]\\ {{\displaystyle \overline{C}}}_1=\left[\begin{array}{ccc}{\mathbf{0}}_{\left(p-m\right)\times \left(n-p\right)}& I_{\left(p-m\right)\times \left(p-m\right)}& {\mathbf{0}}_{\left(p-m\right)\times m}\end{array}\right], {{\displaystyle \overline{C}}}_2=\left[\begin{array}{cc}{\mathbf{0}}_{m\times \left(n-m\right)}& I_{m\times m}\end{array}\right]\end{array}$$ E22

where ${\overline{A}}_{11}$ is Hurwitz and the virtual measurement ${\overline{y}}_1$ presents the protected measurements and ${\overline{y}}_2$ shows the attacked/corrupted measurements.

5.1.2 Attack observation

A SMO is proposed to reconstruct the attack in order to clean up the measurements and states and to allow the use of clean measurement in the control signal.

Define a (sliding mode) observer for the system Eq. (22) as

$$\dot{\overline{z}}=\overline{A}\overline{z}+{\overline{G}}_1\left({\overline{y}}_1-{\overline{z}}_{21}\right)+{\overline{G}}_2\left({\overline{y}}_2-{\overline{z}}_{22}\right)-G_n\upsilon$$ E23

where $\overline{z}=\mathit{col}\left({\overline{z}}_1,{\overline{z}}_{21},{\overline{z}}_{22}\right)$ is conformal with the partition of $\overline{x}$ in Eq. (22). In Eq. (23), $\upsilon$ is a nonlinear injection signal that depends on $\left({\overline{y}}_2-{\overline{z}}_{22}\right)$ and is used to induce a sliding motion in the estimation error space, and

$${\overline{G}}_1=\left[\begin{array}{c}{\overline{A}}_{12a}\\ {\overline{A}}_{22a}-A_{22}^s\\ {\mathbf{0}}_{m\times \left(p-m\right)}\end{array}\right],{\overline{G}}_2=\left[\begin{array}{c}{\overline{A}}_{12b}\\ {\overline{A}}_{22b}\\ {\overline{A}}_{22d}-A_{33}^s\end{array}\right],G_n=\left[\begin{array}{c}{\mathbf{0}}_{\left(n-p\right)\times m}\\ {\mathbf{0}}_{\left(p-m\right)\times m}\\ {\mathrm{I}}_{\mathrm{m}\times m}\end{array}\right]$$ E24

are the gain matrices where ${\overline{A}}_{12a}\in {\mathbb{R}}^{\left(n-p\right)\times \left(p-m\right)}$ , ${\overline{A}}_{22a}\in {\mathbb{R}}^{\left(p-m\right)\times \left(p-m\right)}$ , ${\overline{A}}_{12b}\in {\mathbb{R}}^{\left(n-p\right)\times m}$ , ${\overline{A}}_{22b}\in {\mathbb{R}}^{\left(p-m\right)\times m}$ , ${\overline{A}}_{22d}\in {\mathbb{R}}^{m\times m}$ , and the matrices $A_{22}^s\in {\mathbb{R}}^{\left(p-m\right)\times \left(p-m\right)}$ and $A_{33}^s\in {\mathbb{R}}^{m\times m}$ are user-selected Hurwitz matrices, while $A_{33}^s$ is symmetric negative definite. The injection signal $\upsilon \in {\mathbb{R}}^m$ is defined as

$$\upsilon =-\left(\rho +\eta \right)\frac{{\overline{y}}_2-{\overline{z}}_{22}}{‖{\overline{y}}_2-{\overline{z}}_{22}‖},\rho ,\eta >0$$ E25

where scalar gain $\rho$ will be defined in the sequel, and $\eta$ is a positive design scalar.

Assumption (A4): Matrix $\left(\mathit{sI}-A^{\ast }\right)$ is invertible, where $A^{\ast }=\overline{A}-\overline{B}{D}_2^{-1}{\overline{C}}_2-{\overline{G}}_1{\overline{C}}_1$ .

Defining $\overline{e}=\overline{x}-\overline{z}$ , then it follows $\overline{e}=\mathit{col}\left({\overline{e}}_1,{\overline{e}}_{21},{\overline{e}}_{22}\right)$ where ${\overline{e}}_1={\overline{x}}_1-{\overline{z}}_1$ , ${\overline{e}}_{21}={\overline{x}}_{21}-{\overline{z}}_{21}$ , ${\overline{e}}_{22}={\overline{x}}_{22}-{\overline{z}}_{22}$ . It follows

$$e_{y_2}={\overline{y}}_2-{\overline{z}}_{22}={\overline{e}}_{22}+D_{2}d$$ E26

and by direct substitution from Eqs. (22) and (23) that

$$\dot{\overline{e}}=\left[\begin{array}{ccc}{\overline{A}}_{11}& \mathbf{0}& \mathbf{0}\\ {\overline{A}}_{21a}& A_{22}^s& \mathbf{0}\\ {\overline{A}}_{21b}& {\overline{A}}_{22c}& A_{33}^s\end{array}\right]\overline{e}-\left[\begin{array}{c}{\overline{A}}_{12b}\\ {\overline{A}}_{22b}\\ {\overline{A}}_{22d}-A_{33}^s\end{array}\right]D_{2}d+\left[\begin{array}{c}{\overline{B}}_1\\ {\overline{B}}_{21}\\ {\overline{B}}_{22}\end{array}\right]d+\left[\begin{array}{c}\mathbf{0}\\ \mathbf{0}\\ {\mathrm{I}}_{\mathrm{m}}\end{array}\right]\upsilon$$ E27

The idea is to force a sliding motion on

$$e_{y_2}={\overline{y}}_2-{\overline{z}}_{22}=\mathbf{0}$$ E28

The first main results, based on the SMO with the fixed-gain injection term, is formulated in the following theorem.

Theorem 1: Assuming (A3)–(A4) hold and $m_0>0$ satisfies the condition

$$\begin{array}{ccc}‖\varphi \left(t\right)‖\le m_0{k}_d,& \varphi =\left[\begin{array}{cc}{\overline{A}}_{21b}& {\overline{A}}_{22c}\end{array}\right]{\overline{e}}_{11}-\left({\overline{A}}_{22d}-{\overline{B}}_{22}{D}_2^{-1}\right)D_{2}d,& {\overline{e}}_{11}=\mathit{col}\left({\overline{e}}_1,{\overline{e}}_{21}\right)\end{array}$$ E29

Then, as soon as the sliding mode is established in finite time in Eq. (27) on the sliding surface Eq. (28) by means of the injection term Eq. (25) with $\rho =m_0{k}_d+{‖D_2‖}_{\infty }{l}_d$ , the attack $d$ is asymptotically estimated as

$$\begin{array}{ccc}\widehat{d}=G^{\ast }\left(s\right){\upsilon }_{\mathit{eq}}& \mathit{\text{where}}& G^{\ast }\left(s\right)=C^{\ast }{\left(\mathit{sI}-A^{\ast }\right)}^{-1}{B}^{\ast },\end{array}{B}^{\ast }=\left[\begin{array}{c}{\mathbf{0}}_{\left(n-p\right)\times m}\\ {\mathbf{0}}_{\left(p-m\right)\times m}\\ I_{m\times m}\end{array}\right],C^{\ast }=\left[\begin{array}{cc}{\mathbf{0}}_{m\times \left(n-m\right)}& -D_2^{-1}\end{array}\right]$$ E30

where ${\upsilon }_{\mathit{eq}}$ is the equivalent injection term [31] and a close approximation and ${\overline{\upsilon }}_{\mathit{eq}}$ can be obtained in real time by low-pass filtering of the switching signal Eq. (25) [29]. Replacing ${\upsilon }_{\mathit{eq}}$ by ${\overline{\upsilon }}_{\mathit{eq}}$ in Eq. (30) gives

$$\overline{\widehat{d}}=G^{\ast }\left(s\right){\overline{\upsilon }}_{\mathit{eq}}$$ E31

Proof of the Theorem 1 is omitted for brevity.

Remark 1: The SMO (31) is a dynamic filter that allows reconstructing the time-varying attack $d\left(t\right)$ . This filter is the main novel feature of the proposed observer.

5.1.3 Adaptive-gain attack observer design

In Eq. (29), it was assumed that the perturbation term $\phi$ is locally norm-bounded and $\rho >0$ in Eq. (25) is known. In many practical cases, the boundary of attacks is unknown, and the gain of the sliding mode injection term Eq. (25) in the fixed-gain observer in Eq. (23) can be overestimated. The gain overestimation could increase chattering that is difficult to attenuate.

The constant gain $\rho >0$ can be replaced by an adaptive-gain $\rho \left(t\right)$ by applying the dual layer nested adaptive sliding mode observation algorithm [32], i.e.,

$$\upsilon =-\left(\rho \left(t\right)+\eta \right)\frac{{\overline{y}}_2-{\overline{z}}_{22}}{‖{\overline{y}}_2-{\overline{z}}_{22}‖}$$ E32

A sufficient condition to ensure sliding on $e_{y_2}=\mathbf{0}$ in finite time is

$$\rho \left(t\right)>‖A_{33}^s{e}_{y_2}+\varphi +D_2\dot{d}‖$$ E33

An error signal is defined as

$$\sigma \left(t\right)=\rho \left(t\right)-\frac{1}{\alpha }‖{\overline{\upsilon }}_{\mathit{eq}}\left(t\right)‖-\epsilon$$ E34

where the scalars $0<\alpha <1$ , $\epsilon >0$ . The adaptation dynamics of $\rho \left(t\right)$ in Eq. (32) is defined as [32].

$$\dot{\rho }\left(t\right)=-r\left(t\right)\mathit{sign}\left(\sigma \left(t\right)\right)$$ E35

where the time-varying scalar $r\left(t\right)>0$ satisfies an adaptive scheme. It is assumed that $r\left(t\right)$ has the structure

$$r\left(t\right)={\mathrm{\ell }}_0+\mathrm{\ell }\left(t\right)$$ E36

where ${\mathrm{\ell }}_0$ is a fixed positive scalar. The evolution of $\mathrm{\ell }\left(t\right)$ is chosen to satisfy an adaptive law [32]:

$$\dot{\mathrm{\ell }}\left(t\right)=\left\{\begin{array}{l}\begin{array}{cc}\gamma \left|\sigma \left(t\right)\right|& \mathit{\text{if}}\left|\sigma \left(t\right)\right|>{\sigma }_0\end{array}\\ \begin{array}{cc}0\begin{array}{cc}& \end{array}& \mathit{\text{otherwise}}\end{array}\end{array}\right.$$ E37

where $\gamma >0,{\sigma }_0>0$ are design scalars. The second main results are summarized in Theorem 2 as:

Theorem 2: Consider the system in Eq. (27) and

$$a\left(t\right)=A_{33}^s{e}_{y_2}+\varphi +D_2\dot{d}$$ E38

and assume that $\left|a\left(t\right)\right|<a_0,\left|\dot{a}\left(t\right)\right|<a_1$ , where $a_0$ and $a_1$ are finite but unknown. A SMO is designed as in Eq. (23) with the adaptive injection term in Eqs. (32)–(37). If $\epsilon >0$ in (34) is chosen to satisfy

$$\frac{1}{4}{\epsilon }^2>{\sigma }_0^2+\frac{1}{\gamma }{\left(\frac{{\mathit{qa}}_1}{\alpha }\right)}^2$$ E39

for any given ${\sigma }_0$ , $q>1$ , and, $0<\alpha <1$ , then the injection term (32) exploiting the dual layer adaptive scheme given by Eqs. (35)–(37) drives $\sigma \left(t\right)$ to a domain $\left|\sigma \left(t\right)\right|<\epsilon /2$ in finite time and consequently ensures a sliding motion $e_y=0$ can be reached in finite time and sustained thereafter. The gains $r\left(t\right)$ and $\rho \left(t\right)$ remain bounded. The sensor attack signal $d\left(t\right)$ is reconstructed as in Eq. (30) with the equivalent adaptive injection term ${\upsilon }_{\mathit{eq}}$ or ${\overline{\upsilon }}_{\mathit{eq}}$ .

Proof of Theorem 2 is based on the results in [32] and is omitted for brevity.

Remark 2: The proposed unit vector injection gain-adaptation algorithm in Eqs. (32)–(37) does not require the knowledge of the boundaries $k_d,l_d>0$ in $‖d‖<k_d$ and $‖\dot{d}‖<l_d$ .

5.2 State estimation and attack reconstruction in linear systems by using super twisting SMO

Consider the completely observable linearized system Eq. (11) with $C_1\left(x\right)=C_{1}x$ , $C_2\left(x\right)=C_{2}x$ , $B_1\left(x\right)=B$ , that is,

$$\dot{x}=\mathit{Ax}+B_1{d}_x\left(t\right),{\overline{y}}_1=C_{1}x\begin{array}{cc},& {\overline{y}}_2=C_{2}x+{\overline{D}}_1{d}_y\left(t\right)\end{array}$$ E40

where $B_1\in {\mathbb{R}}^{n\times m_1}$ , $C_1\in {\mathbb{R}}^{\left(p-\left(m-m_1\right)\right)\times n}$ , $C_2\in {\mathbb{R}}^{\left(m-m_1\right)\times n}$ .

Assumption (A5): The number of uncorrupted/protected measurements is equal or larger than the number of state/plant attack, i.e., $p_1=p-\left(m-m_1\right)\ge m_1$ .

The system Eq. (40) is assumed to have an input-output vector relative degree $r=\left\{r_1,r_2,\dots ,r_{p_1}\right\}$ , where relative degree $r_i$ for $i=1,2,\dots ,p_1$ is defined as follows:

$$\begin{array}{l}\begin{array}{cc}{C}_{1i}{A}^j{B}_1=0& \begin{array}{ccc}\mathit{\text{for}}& \mathit{all}& j<r_i-1\end{array}\end{array}\\ C_{1i}{A}^{r_i-1}{B}_1\ne 0\end{array}$$ E41

Without loss of generality, it is assumed that $r_1\le \dots \le r_{p_1}$ .

5.2.1 Attack observation

Assumption (A6): there exists a full rank matrix.

$$C_a=\left[\begin{array}{c}\begin{array}{c}{C}_1\\ ⋮\\ C_1{A}^{r_{{\alpha }_1}-1}\\ ⋮\end{array}\\ C_{p_1}\\ ⋮\\ C_{p_1}{A}^{r_{{\alpha }_{p_1}}-1}\end{array}\right]$$ E42

where integers $1\le r_{{\alpha }_i}\le r_i$ are such that $\mathit{rank}\left(C_{a}B\right)=\mathit{rank}\left(B\right)$ and $r_{{\alpha }_i}$ are chosen such that ${\sum }_{i=1}^{p_1}{r}_{{\alpha }_i}$ is minimal.

The following SMO [33] is used to estimate the states of system Eq. (40):

$$\dot{\widehat{x}}=A\widehat{x}+G_l\left(y_a-C_a\widehat{x}\right)+G_n{\upsilon }_c\left(y_a-C_a\widehat{x}\right)$$ E43

where the matrices of appropriate dimensions $G_l$ and $G_n$ are to be designed, and ${\upsilon }_c\left(.\right)$ is an injection vector

$${\upsilon }_c\left(y_a-C_a\widehat{x}\right)=\left\{\begin{array}{c}\begin{array}{cc}-\rho \frac{P\left(y_a-C_a\widehat{x}\right)}{‖P\left(y_a-C_a\widehat{x}\right)‖}& \begin{array}{cc}\mathit{\text{if}}& \left(y_a-C_a\widehat{x}\right)\ne 0\end{array}\end{array}\\ \begin{array}{cc}0& \begin{array}{cccc}\begin{array}{cccc}& & & \end{array}& \mathit{\text{otherwise}}& & \end{array}\end{array}\end{array}\right.$$ E44

where $\rho >0$ is larger than the upper bound of unknown input $d\left(t\right)$ .

The definition of the symmetric positive definite matrix $P$ can be found in [33]. The auxiliary output $y_a$ is defined by

$$y_a=\left[\begin{array}{c}\begin{array}{c}\begin{array}{c}{y}_1\\ \nu \left(y_1-{y^1}_1\right)\end{array}\\ ⋮\\ \nu \left({\tilde{y}}_1^{r_1-1}-{\tilde{y}}_1^{r_1-1}\right)\\ ⋮\end{array}\\ y_{p_1}\\ ⋮\\ \nu \left({\tilde{y}}_{p_1}^{r_{p_1}-1}-y_{p_1}^{r_{p_1}-1}\right)\end{array}\right]$$ E45

where the constituent signals in Eq. (45) are given from the continuous second-order sliding mode observer as

$$\begin{array}{l}{{\displaystyle \dot{y}}}_i^1=\nu \left(y_i-y_i^1\right)\\ {{\displaystyle \dot{y}}}_i^2=E_1\nu \left({{\displaystyle \tilde{y}}}_i^2-y_i^2\right)\\ ⋮\\ {{\displaystyle \dot{y}}}_i^{r_{\mathit{\alpha i}}-1}=E_{r_{\mathit{\alpha i}}-2}\nu \left({{\displaystyle \tilde{y}}}_i^{r_{\mathit{\alpha i}}-1}-y_i^{r_{\mathit{\alpha i}}-1}\right)\end{array}$$ E46

for $1\le i\le p_1$ , with

$$\begin{array}{ccc}{\tilde{y}}_i^1=y_i& ,& {\tilde{y}}_i^j=\nu \left({\tilde{y}}_i^{j-1}-y_i^{j-1}\right),\begin{array}{cc}& 2\le j\le \end{array}{r}_{{\alpha }_i}-1\end{array}$$ E47

The scalar function $E_i$ is defined as

$$E_i=1\mathit{\text{if}}\left|{\tilde{y}}_j^{i+1}-y_j^{i+1}\right|\le \mathit{\text{ε for}}\mathit{all}j\le i,\mathit{\text{else}}{E}_i=0$$ E48

and the continuous injection term $\nu \left(.\right)$ is given by the super twisting algorithm [34]:

$$\begin{array}{l}\nu \left(s\right)=\xi \left(s\right)+{\lambda }_s{\left|s\right|}^{\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}\mathit{sign}\left(s\right)\\ {\displaystyle \dot{\xi }}\left(s\right)={\beta }_s\mathit{sign}\left(s\right)\begin{array}{cc},& {\lambda }_s,{\beta }_s>0\end{array}\end{array}$$ E49

Theorem 3: Assuming the assumptions (A5) and (A6) hold for system Eq. (40), then state/plant attacks are reconstructed as follows:

$${\widehat{d}}_x={\left({\left(C_{a}B\right)}^T{C}_{a}B\right)}^{-1}{\left(C_{a}B\right)}^T{C}_a{G}_n{\left({\upsilon }_c\right)}_{\mathit{eq}}$$ E50

Proof: Defining the state estimation error as $e=x-\widehat{x}$ and the augmented output estimation error $e_y=C_{a}x-\overline{y}$ with

$$\begin{array}{cc}{e}_y={\left[e_1^1,\dots ,e_1^{r_{\mathit{\alpha i}}-1},\dots ,e_{p_1}^1,\dots ,e_{p_1}^{r_{\mathit{\alpha i}}-1}\right]}^T,& y={\left[y_1^1,\dots ,y_1^{r_{\mathit{\alpha i}}-1},\dots ,y_{p_1}^1,\dots ,y_{p_1}^{r_{\mathit{\alpha i}}-1}\right]}^T\end{array}$$ E51

then it follows that

$$\dot{e}=x-\dot{\widehat{x}}=\mathit{Ae}+B_1{d}_x\left(t\right)-G_l\left(y_a-C_a\widehat{x}\right)-G_n{\upsilon }_c\left(y_a-C_a\widehat{x}\right)$$ E52

By choosing suitable gains ${\lambda }_s$ and ${\beta }_s$ in the output injections Eq. (49), then.

$$y_a=C_{a}x$$ E53

for all $t>T$ [33]. Then, the error dynamics Eq. (52) is rewritten as

$$\dot{e}=\left(\overline{A}-G_l{C}_a\right)e+{\overline{B}}_1{d}_x\left(t\right)-G_n{\upsilon }_c\left(C_{a}e\right)$$ E54

Since $\mathit{rank}\left(C_a{\overline{B}}_1\right)=\mathit{rank}\left({\overline{B}}_1\right)$ and by assumption the invariant zeros of the triple $\left(A,B,C_a\right)$ lie in the left half plane, based on the design methodologies in [35], It follows that $e=0$ is an asymptotically stable equilibrium point of Eq. (52) and dynamics are independent of $d_x\left(t\right)$ once a sliding motion on the sliding manifold $s=C_{a}e=0$ has been attained. During the sliding mode $\dot{s}=s=0$ , it is

$$\dot{s}=C_a\dot{e}=C_a\left(\overline{A}-G_l{C}_a\right)e+C_a{\overline{B}}_1{d}_x\left(t\right)-C_a{G}_n{\upsilon }_c\left(C_{a}e\right)=0$$ E55

as $e\to 0$ ; then

$$C_a{G}_n{\left({\upsilon }_c\right)}_{\mathit{eq}}\to C_a{\overline{B}}_1{d}_x\left(t\right)$$ E56

where ${\left({\upsilon }_c\right)}_{\mathit{eq}}$ is the equivalent output error injection required to maintain the system on the sliding manifold. Since $C_a{\overline{B}}_1$ is full rank, the attack reconstruction is obtained as (50).

According to (A1), ${\overline{D}}_1$ is full rank; then sensor attacks in Eq. (40) are reconstructed

$${\widehat{d}}_y\left(t\right)={{\overline{D}}_1}^{-1}\left({\overline{y}}_2-C_2\widehat{x}\right)$$ E57

5.3 The state and disturbance observer for nonlinear systems using higher-order sliding mode differentiator

Consider the locally stable system Eq. (11) where ${\overline{y}}_1$ and $B_1\left(x\right)$ are ${\overline{y}}_1={\left[\begin{array}{cccc}{y}_1& y_2& ,\dots ,& y_{p_1}\end{array}\right]}^T$ , $B=\left[b_1,b_2,\dots ,b_{m_1}\right]\in {\mathbb{R}}^{n\times m_1}$ , $b_i\in {\mathbb{R}}^n,\forall i=1,\dots ,m_1$ are smooth vector fields defined on an open $\mathrm{\Omega }\subset {\mathbb{R}}^n$ . According to (A5), we consider $p_1=m_1$ here. The following properties introduced by Isidori [36] are assumed for $x\in \mathrm{\Omega }$ .

Assumption (A7): The system in Eq. (11) is assumed to have vector relative degree $r=\left\{r_1,r_2,\dots ,r_{m_1}\right\}$ and total relative degree $r_t=\sum _{i=1}^{m_1}{r}_i,r_t\le n$ , i.e.,

$$\begin{array}{l}{L}_{\mathit{bj}}{L}_f^k{y}_i\left(x\right)=0\forall j=1,\dots ,m_1,\forall k<r_i-1,\forall i=1,\dots ,m_1\\ L_{\mathit{bj}}{L}_f^{r_i-1}{y}_i\left(x\right)\ne 0\text{for}\mathrm{at}\text{least}\mathrm{one}1\le j\le m_1\end{array}$$ E58

Assumption (A8): The following Lie derivative matrix is of full rank.

$$L\left(x\right)=\left[\begin{array}{cccc}{L}_{b_1}{L}_f^{r_1-1}{y}_1& L_{b_2}{L}_f^{r_1-1}{y}_1& \cdots & L_{b_{m_1}}{L}_f^{r_1-1}{y}_1\\ L_{b_1}{L}_f^{r_2-1}{y}_2& L_{b_2}{L}_f^{r_2-1}{y}_2& \cdots & L_{b_{m_1}}{L}_f^{r_2-1}{y}_2\\ ⋮& ⋮& ⋮& ⋮\\ L_{b_1}{L}_f^{r_{m_1}-1}{y}_{m_1}& L_{b_2}{L}_f^{r_m-1}{y}_{m_1}& \cdots & L_{b_{m_1}}{L}_f^{r_{m_1}-1}{y}_{m_1}\end{array}\right]$$ E59

Assumption (A9): The distribution $\mathrm{\Gamma }=\mathit{\text{span}}\left\{b_1,b_2,\dots ,b_{m_1}\right\}$ is involutive [36].

The system given by Eq. (11) with the involutive distribution $\mathrm{\Gamma }$ and total relative degree $r_t$ can be rewritten as

$$\begin{array}{l}{\dot{\delta }}_i={\left[\begin{array}{cccc}\begin{array}{cc}0& 1\end{array}& 0& \cdots & 0\\ \begin{array}{cc}0& 0\end{array}& 1& \cdots & 0\\ \begin{array}{cc}⋮& ⋮\end{array}& ⋮& \cdots & ⋮\\ \begin{array}{cc}0& 0\end{array}& 0& 0& 0\end{array}\right]}_{r_i\times r_i}{\delta }_i+\left[\begin{array}{c}0\\ 0\\ ⋮\\ {L_f}^{r_i}{y}_i\left(x\right)\end{array}\right]+\left[\begin{array}{c}0\\ 0\\ ⋮\\ {\sum }_{j=1}^{m_1}{L}_{b_j}{L_f}^{r_i-1}{y}_i\left(x\right)d\left(t\right)\end{array}\right],\forall i=1,\dots ,m_1\\ \dot{\gamma }=g\left(\delta ,\gamma \right)\end{array}$$ E60

where $\delta ={\left[\begin{array}{cccc}{\delta }_1& {\delta }_2& \cdots & {\delta }_{m_1}\end{array}\right]}^T$ and

$${\delta }_i=\left[\begin{array}{c}{\delta }_{i1}\\ {\delta }_{i2}\\ ⋮\\ {\delta }_{{\mathit{ir}}_1}\end{array}\right]=\left[\begin{array}{c}{\eta }_{i1}\left(x\right)\\ {\eta }_{i2}\left(x\right)\\ ⋮\\ {\eta }_{{\mathit{ir}}_1}\left(x\right)\end{array}\right]=\left[\begin{array}{c}{y}_i\left(x\right)\\ L_f{y}_i\left(x\right)\\ ⋮\\ {L_f}^{r_1-1}{y}_i\left(x\right)\end{array}\right]\in {\mathbb{R}}^{r_i}\begin{array}{cc}\forall i=1,\dots ,m_1,& \end{array}\gamma =\left[\begin{array}{c}{\gamma }_1\\ {\gamma }_2\\ ⋮\\ {\gamma }_{n-r}\end{array}\right]=\left[\begin{array}{c}{\eta }_{r+1}\left(x\right)\\ {\eta }_{r+2}\left(x\right)\\ ⋮\\ {\eta }_n\left(x\right)\end{array}\right]$$ E61

With an involutive distribution $\mathrm{\Gamma }$ as defined in (A9), it is always possible to identify the variables ${\eta }_{r+1}\left(x\right),\dots ,{\eta }_n\left(x\right)$ which satisfy

$$L_{b_j}{\eta }_i\left(x\right)=0\forall i=r+1,\dots ,n,\forall j=1,\dots ,m_1$$ E62

Assumption (A10): The norm-bounded solution of the internal dynamics $\dot{\gamma }=g\left(\delta ,\gamma \right)$ is assumed to be locally asymptotically stable [29].

If assumption (A9) is satisfied, then it is always possible to find $n-r$ functions ${\eta }_{r+1}\left(x\right),\dots ,{\eta }_n\left(x\right)$ such that

$$\mathrm{\Psi }\left(x\right)=\mathit{col}\left\{{\eta }_{11}\left(x\right),\dots ,{\eta }_{1r_1}\left(x\right),\dots ,{\eta }_{m_{1}1}\left(x\right),\dots ,{\eta }_{m_1{r}_{m_1}}\left(x\right),\right.\left.{\eta }_{r+1}\left(x\right),\dots ,{\eta }_n\left(x\right)\right\}\in {\mathbb{R}}^n$$ E63

is a local diffeomorphism in a neighborhood of any point $x\in \overline{\mathrm{\Omega }}\subset \mathrm{\Omega }\subset {\mathbb{R}}^n$ , i.e.,

$$x={\mathrm{\Psi }}^{-1}\left(\delta ,\gamma \right)$$ E64

In order to estimate the derivatives ${\delta }_{\mathit{ij}}\left(t\right)\forall i=1,\dots ,m_1,\forall j=1,\dots ,r_i$ of the output.

$y_i$ in finite time, higher-order sliding mode differentiators [28] are used here

$$\begin{array}{l}\begin{array}{ccc}{{\displaystyle \dot{z}}}_0^i=v_0^i,& v_0^i=-{\lambda }_0^i{\left|z_0^i-y_i\left(t\right)\right|}^{\left(r_i/\left(r_i+1\right)\right)}\mathit{sign}\left(z_0^i-y_i\left(t\right)\right)+z_1^i,& {{\displaystyle \dot{z}}}_1^i=v_1^i\end{array}\\ ⋮\\ \begin{array}{ccc}{{\displaystyle \dot{z}}}_{r_i-1}^i=v_{r_i-1}^i,& v_{r_i-1}^i=-{\lambda }_{r_i-1}^i{\left|z_{r_i-1}^i-v_{r_i-2}^i\right|}^{\left(1/2\right)}\mathit{sign}\left(z_{r_i-1}^i-v_{r_i-2}^i\right)+z_{r_i}^i,& {{\displaystyle \dot{z}}}_{r_i}^i=-{\lambda }_{r_i}^i\mathit{sign}\left(z_{r_i}^i-v_{r_i-1}^i\right)\end{array}\end{array}$$ E65

for $i=1,\dots ,m_1$ . By construction,

$$\begin{array}{l}{{\displaystyle \widehat{\delta }}}_1^1={{\displaystyle \widehat{\eta }}}_1^1\left(x\right)=z_0^1,\dots ,{{\displaystyle \widehat{\delta }}}_1^1={{\displaystyle \widehat{\eta }}}_{r_1}^1\left(x\right)=z_{r_1-1}^1,{\displaystyle \dot{{\displaystyle \widehat{\delta }}}}{ }_{r_1}^1={\displaystyle \dot{{\displaystyle \widehat{\eta }}}}{ }_{r_1}^1\left(x\right)=z_{r_1}^1\\ ⋮\\ {{\displaystyle \widehat{\delta }}}_1^{m_1}={{\displaystyle \widehat{\eta }}}_1^{m_1}\left(x\right)=z_0^{m_1},\dots ,{{\displaystyle \widehat{\delta }}}_{r_{m_1}}^{m_1}={{\displaystyle \widehat{\eta }}}_{r_{m_1}}^{m_1}\left(x\right)=z_{r_{m_1}-1}^{m_1},{{\displaystyle \dot{{\displaystyle \widehat{\delta }}}}}_{r_1}^{m_1}={{\displaystyle \dot{{\displaystyle \widehat{\eta }}}}}_{r_{m_1}}^{m_1}\left(x\right)=z_{r_{m_1}}^1\end{array}$$ E66

Therefore, the following exact estimates are available in finite time:

$$\begin{array}{l}{{\displaystyle \widehat{\delta }}}_i={\left({{\displaystyle \widehat{\delta }}}_{i1},{{\displaystyle \widehat{\delta }}}_{i2},\dots ,{{\displaystyle \widehat{\delta }}}_{{\mathit{ir}}_1}\right)}^T={\left({{\displaystyle \widehat{\eta }}}_{i1}\left({\displaystyle \widehat{x}}\right),{{\displaystyle \widehat{\eta }}}_{i2}\left({\displaystyle \widehat{x}}\right),\dots ,{{\displaystyle \widehat{\eta }}}_{{\mathit{ir}}_1}\left({\displaystyle \widehat{x}}\right)\right)}^T\in {\mathbb{R}}^{r_i}\\ \begin{array}{cc}\forall i=1,\dots ,m_1,& {\displaystyle \widehat{\delta }}={\left({{\displaystyle \widehat{\delta }}}^1,{{\displaystyle \widehat{\delta }}}^2,\dots ,{{\displaystyle \widehat{\delta }}}^{m_1}\right)}^T\end{array}\in {\mathbb{R}}^{r_t}\end{array}$$ E67

Next, integrate Eq. (60) with $\delta$ replaced by $\widehat{\delta }$ ; estimate of internal dynamics is

$$\dot{\widehat{\gamma }}=g\left(\widehat{\delta },\widehat{\gamma }\right)$$ E68

and with some initial condition from the stability domain of the internal dynamics, a asymptotic estimate $\widehat{\gamma }$ can be obtained locally

$${\displaystyle \widehat{\gamma }}=\left(\begin{array}{c}{{\displaystyle \widehat{\gamma }}}_1\\ {{\displaystyle \widehat{\gamma }}}_2\\ ⋮\\ {{\displaystyle \widehat{\gamma }}}_{n-r}\end{array}\right)=\left(\begin{array}{c}{{\displaystyle \widehat{\eta }}}_{r+1}\left({\displaystyle \widehat{x}}\right)\\ {{\displaystyle \widehat{\eta }}}_{r+2}\left({\displaystyle \widehat{x}}\right)\\ ⋮\\ {{\displaystyle \widehat{\eta }}}_n\left({\displaystyle \widehat{x}}\right)\end{array}\right)$$ E69

Therefore, the asymptotic estimate for the mapping (63) is identified as

$$\mathrm{\Psi }\left(\widehat{x}\right)=\mathit{col}\left\{{\widehat{\eta }}_{11}\left(\widehat{x}\right),\dots ,{\widehat{\eta }}_{1r_1}\left(\widehat{x}\right),\dots ,{\widehat{\eta }}_{m_{1}1}\left(\widehat{x}\right),\dots ,{\widehat{\eta }}_{m_1{r}_{m_1}}\left(\widehat{x}\right)\left(,{\widehat{\eta }}_{r+1}\left(\widehat{x}\right),\dots ,{\widehat{\eta }}_n\left(\widehat{x}\right)\right\}\right.$$ E70

asymptotic estimate $\widehat{x}$ of the state vector $x$ can be identified via Eqs. (67) and (69)

$$\widehat{x}={\mathrm{\Psi }}^{-1}\left(\widehat{\delta },\widehat{\gamma }\right)$$ E71

Since the finite-time exact estimates ${\widehat{\dot{\delta }}}_{{\mathit{ir}}_i}$ of ${\dot{\delta }}_{{\mathit{ir}}_i}$ , $\forall i=1,\dots ,m_1$ are available via the higher-order sliding mode differentiator, and using the estimates $\begin{array}{c}\widehat{\delta },\widehat{\gamma }\end{array}$ for $\delta ,\gamma$ , an asymptotic estimate $\widehat{d}\left(t\right)$ of disturbance $d\left(t\right)$ in Eq. (11) is identified as [28].

$$\widehat{d}\left(t\right)=L^{-1}\left({\mathrm{\Psi }}^{-1}\left(\widehat{\delta },\widehat{\gamma }\right)\right)\left[\left(\begin{array}{c}{\widehat{\dot{\delta }}}_{1r_1}\\ {\widehat{\dot{\delta }}}_{2r_2}\\ ⋮\\ {\widehat{\dot{\delta }}}_{m_1{r}_{m_1}}\end{array}\right)-\left(\begin{array}{c}{L}_f^{r_1}{y}_{11}\left({\mathrm{\Psi }}^{-1}\left(\widehat{\delta },\widehat{\gamma }\right)\right)\\ L_f^{r_2}{y}_{12}\left({\mathrm{\Psi }}^{-1}\left(\widehat{\delta },\widehat{\gamma }\right)\right)\\ ⋮\\ L_f^{r_{m_1}}{y}_{1m_1}\left({\mathrm{\Psi }}^{-1}\left(\widehat{\delta },\widehat{\gamma }\right)\right)\end{array}\right)\right]$$ E72

where $L\left({\mathrm{\Psi }}^{-1}\left(\widehat{\delta },\widehat{\gamma }\right)\right)={\sum }_{j=1}^{m_1}{L}_{b_j}{L_f}^{r_i-1}{y}_{1i}\left(x\right)$ . Finally, $\widehat{x}\left(t\right)$ and $\widehat{d}\left(t\right)$ are obtained.

from Eqs. (71) and (72).

Remark 3: The convergence $\widehat{d}\to d$ can be achieved only locally and as time increases due to the local asymptotic stability of the norm-bounded solution of the internal dynamics $\dot{\gamma }=g\left(\delta ,\gamma \right)$ . However convergence will be achieved in finite time if the total relative degree $r=n$ and no internal dynamics exist.

Considering Eq. (11) and ${\overline{D}}_1$ is full rank, sensor attack can be reconstructed as

$${\widehat{d}}_y\left(t\right)={{\overline{D}}_1}^{-1}\left({\overline{y}}_2-C_2\left(\widehat{x}\right)\right)$$ E73

5.4 Attack reconstruction in nonlinear system by sparse recovery algorithm

In some applications, there are a limited number of measurements, $p$ , and more sources of attack, $m$ . Previously, we investigated the cases where $p>m$ . Now, consider system (5) with more attacks than measurements, $m>p$ .

Notice that a more general format of (5) is considered here where matrix $D$ is a function of $x$ as well.

Assumption (A11): Assume that the attack vector $d\left(t\right)$ is sparse, meaning that numerous attacks are possible, but the attacks are not coordinated, and only few nonzero attacks happen at the same time.

5.4.1 Sparse recovering algorithm

The problem of recovering an unknown input signal from measurements is well known, as a left invertibility problem, as seen in several works [30, 37], but this problem was only treated in the case where the number of measurements is equal or greater than the number of unknown inputs. The left invertibility problem in the case of fewer measurements than unknown inputs has no solution or more exactly has an infinity of solutions.

In particular, the objective of exact recovery under sparse assumptions denoted for the sake of simplicity as “sparse recovery” (SR) is to find a concise representation of a signal using a few atoms from some specified (over-complete) dictionary,

$$\xi =\mathrm{\Phi }\overline{s}+{\epsilon }_0$$ E74

where $\overline{s}\in {\mathbb{R}}^N$ are the unknown inputs with no more than $j$ nonzero entries, $\xi \in {\mathbb{R}}^M$ are the measurements, ${\epsilon }_0$ is a measurement noise, and $\mathrm{\Phi }\in {\mathbb{R}}^{M\times N}$ is the dictionary where $M\ll N$ .

Definition 1: The Restricted Isometry Property (RIP) condition of $j$ -order with constant ${\varsigma }_j\in \left(0,1\right)$ ( ${\varsigma }_j$ is as small as possible for computational reasons) of the matrix $\mathrm{\Phi }$ yields

$$\left(1-{\varsigma }_{\overline{s}}\right){‖\overline{s}‖}_2^2\le {‖\mathrm{\Phi }\overline{s}‖}_2^2\le \left(1+{\varsigma }_{\overline{s}}\right){‖\overline{s}‖}_2^2$$ E75

for any $j$ sparse of signal $\overline{s}$ . Considering ${\mathrm{\Phi }}_{\mathrm{\Gamma }}$ as the index set of nonzero elements of $\overline{s}$ , then Eq. (75) is equivalent to [23]:

$$1-{\varsigma }_{{\displaystyle \overline{s}}}\le \mathit{eig}\left({\mathrm{\Phi }}_{\mathrm{\Gamma }}^T{\mathrm{\Phi }}_{\mathrm{\Gamma }}\right)\le 1+{\varsigma }_{{\displaystyle \overline{s}}}$$ E76

where ${\mathrm{\Phi }}_{\mathrm{\Gamma }}$ is the sub-matrix of $\mathrm{\Phi }$ with active nodes.

The problem of SR is often cast as an optimization problem that minimizes a cost function constructed by leveraging the observation error term and the sparsity inducing term [37], i.e.,

$${\overline{s}}^{\ast }=\begin{array}{cc}arg& \underset{\overline{s}\in {\mathbb{R}}^N}{min}\frac{1}{2}\end{array}{‖\xi -\mathrm{\Phi }\overline{s}‖}_2^2+\lambda \mathrm{\Theta }\left(\overline{s}\right)$$ E77

In Eq. (77) the original sparsity term is the quasi norm ${\left|\overline{s}\right|}_0$ ; but as long as the RIP conditions hold, it can be replaced by $\mathrm{\Theta }\left(\overline{s}\right)={‖\overline{s}‖}_1\triangleq {\sum }_i\left|{\overline{s}}_i\right|$ . Note that $\lambda >0$ in Eq. (77) is the balancing parameter and ${\overline{s}}^{\ast }$ is the critical point, i.e., the solution of Eq. (74). Typically, for sparse vectors $\overline{s}$ with j-sparsity, where $j$ must be equal or smaller than $\frac{M-1}{2}$ [37], the solution to the SR problem is unique and coincides with the critical point of Eq. (74) providing that RIP condition for $\mathrm{\Phi }$ with order $2j$ is verified. In other words, in order to guarantee the existence of a unique solution to the optimization problem Eq. (74), $\mathrm{\Phi }$ should satisfy restricted isometry property [37].

Under the sparse assumption of $\overline{s}$ and the fulfillment of the j-RIP condition of the matrix $\mathrm{\Phi }$ , the estimation algorithm proposed in [37] is

$$\mu {\displaystyle \dot{v}}\left(t\right)=-{⌈v\left(t\right)+\left({\mathrm{\Phi }}^T\mathrm{\Phi }-I_{N\times N}\right)a\left(t\right)-{\mathrm{\Phi }}^T\xi ⌋}^{\beta } ,\text{and} {\displaystyle \widehat{{\displaystyle \overline{s}}}}\left(t\right)=a\left(t\right)$$ E78

where $v\in {\mathbb{R}}^N$ is the state vector, $\widehat{\overline{s}}\left(t\right)$ represents the estimate of the sparse signal $\overline{s}$ of (74), and $\mu >0$ is a time-constant determined by the physical properties of the implementing system. ${⌈.⌋}^{\beta }=$ ${\left|.\right|}^{\beta }\mathit{sign}\left(.\right)$ and $a\left(t\right)=H_{\lambda }\left(v\right)$ where $H_{\lambda }\left(.\right)$ is a continuous soft thresholding function:

$$H_{\lambda }\left(v\right)=max\left(\left|v\right|-\lambda ,0\right)sgn\left(v\right)$$ E79

where $\lambda >0$ is chosen with respect to the noise and the minimum absolute value of the nonzero terms.

Under Definition 1, the state $v$ of Eq. (78) converges in finite time to its equilibrium point $v^{\ast }$ , and $\widehat{\overline{s}}\left(t\right)$ in Eq.(78) converges in finite time to $\widehat{s}\ast$ of Eq. (77).

5.4.2 Attack reconstruction

The measured output under attack $y$ of the system Eq. (5) is fed to the input of the low-pass filter that facilitates filtering out the possible measurement noise

$$\dot{z}=\frac{1}{\tau }\left(-z+C\left(x\right)+D\left(x\right)d\left(t\right)\right)$$ E80

The filter output $z\in {\mathbb{R}}^p$ is available. Then, system Eq. (5) with filter Eq. (80) is rewritten as

$$\left\{\begin{array}{l}\dot{\xi }=\eta \left(\xi \right)+\mathrm{\Omega }d\left(t\right)\\ \psi =\overline{C}\xi \end{array}\right.$$ E81

where $\psi \in {\mathbb{R}}^p$ , and

$$\begin{array}{l}\xi ={\left[\begin{array}{c}z\\ x\end{array}\right]}_{\left(p+n\right)\times 1},\eta \left(\xi \right)=\left[\begin{array}{cc}-\frac{1}{\tau }{I}_{p\times p}& 0\\ 0& 0\end{array}\right]\left[\begin{array}{c}z\\ x\end{array}\right]+\left[\begin{array}{c}\frac{1}{\tau }C\left(x\right)\\ f\left(x\right)\end{array}\right],\\ C=\left[\begin{array}{cccc}{C}_1& C_2& \dots & C_{p+n}\end{array}\right]=\left[{\begin{array}{cc}{I}_{p\times p}& 0\end{array}}_{p\times n}\right]\end{array}$$ E82

$$\mathrm{\Omega }=\left[\begin{array}{c}\frac{1}{\tau }D\left(x\right)\\ B\left(x\right)\end{array}\right]=\left[\begin{array}{cccc}{\mathrm{\Omega }}_1& {\mathrm{\Omega }}_2& \dots & {\mathrm{\Omega }}_m\end{array}\right]\begin{array}{cc},{\mathrm{\Omega }}_i\in {\mathbb{R}}^{p+n}& \forall i=1,\dots ,m\end{array}\begin{array}{cc}& \end{array}$$

If assumption (A2), (A7), and (A9) hold for system Eq. (81), i.e., the relative degree vector of Eq. (81) is $r=\left\{r_1,r_2,\dots ,r_p\right\}$ , the distribution $\mathrm{\Gamma }=\mathit{\text{span}}\left\{{\mathrm{\Omega }}_1,{\mathrm{\Omega }}_2,\dots ,{\mathrm{\Omega }}_m\right\}$ is involutive, and if zero dynamics exist, they are assumed asymptotically stable and may be left alone. Here it is assumed that there are no zero dynamics in system Eq. (81) and it is presented as

$${\dot{\mathrm{\Upsilon }}}_i=\left[\begin{array}{ccccc}0& 1& 0& \cdots & 0\\ 0& 0& 0& \cdots & 0\\ ⋮& ⋮& ⋮& \cdots & ⋮\\ 0& 0& 0& 0& 0\end{array}\right]{\mathrm{\Upsilon }}_i+\left[\begin{array}{c}0\\ 0\\ ⋮\\ L_f^{r_i}{\psi }_i\left(\xi \right)\end{array}\right]+\left[\begin{array}{c}0\\ 0\\ ⋮\\ \sum _{j=1}^m{L}_{{\mathrm{\Omega }}_j}{L}_f^{r_i-1}{\psi }_i\left(\xi \right)d_j\end{array}\right]\begin{array}{cc},& {\mathrm{\Upsilon }}_i=\left[\begin{array}{c}{\mathrm{\Upsilon }}_1^i\left(\xi \right)\\ {\mathrm{\Upsilon }}_2^i\left(\xi \right)\\ ⋮\\ {\mathrm{\Upsilon }}_{r_i}^i\left(\xi \right)\end{array}\right]=\left[\begin{array}{c}{\psi }_i\left(\xi \right)\\ L_f{\psi }_i\left(\xi \right)\\ \\ L_f^{r_i-1}{\psi }_i\left(\xi \right)\end{array}\right]\end{array}$$ E83

for $i=1,\dots ,p$ , where ${\psi }_i\left(\xi \right)$ is the $i^{\mathit{th}}$ entry of vector $\psi \left(\xi \right)$ and satisfies

$${\dot{\mathrm{\Upsilon }}}_{r_i}^i\left(\xi \right)=L_f^{r_i}{\psi }_i\left(\xi \right)+\sum _{j=1}^m{L}_{{\mathrm{\Omega }}_j}{L}_f^{r_i-1}{\psi }_i{d}_j,i=1,\dots ,p$$ E84

Then, the following algebraic equation is found from Eq. (84):

$$Z_p=F\left(\xi \right)d\left(t\right)$$ E85

where $Z_p\in {\mathbb{R}}^p$ , $F\left(\xi \right)\in {\mathbb{R}}^{p\times m}$ , and

$$\begin{array}{cc}{Z}_p=\left[\begin{array}{c}{\dot{\mathrm{\Upsilon }}}_{r_1}^1\\ ⋮\\ {\dot{\mathrm{\Upsilon }}}_{r_p}^p\end{array}\right]-\left[\begin{array}{c}{L}_f^{r_1}{\psi }_1\left(\xi \right)\\ ⋮\\ L_f^{r_p}{\psi }_p\left(\xi \right)\end{array}\right],& F\left(\xi \right)=\left[\begin{array}{cccc}{L}_{{\mathrm{\Omega }}_1}{L}_f^{r_1-1}{\psi }_1& L_{{\mathrm{\Omega }}_2}{L}_f^{r_1-1}{\psi }_1& \cdots & L_{{\mathrm{\Omega }}_{\alpha }}{L}_f^{r_1-1}{\psi }_1\\ L_{{\mathrm{\Omega }}_1}{L}_f^{r_2-1}{\psi }_2& L_{{\mathrm{\Omega }}_2}{L}_f^{r_2-1}{\psi }_2& & L_{{\mathrm{\Omega }}_{\alpha }}{L}_f^{r_2-1}{\psi }_2\\ ⋮& & & ⋮\\ L_{{\mathrm{\Omega }}_1}{L}_f^{r_p-1}{\psi }_p& L_{{\mathrm{\Omega }}_2}{L}_f^{r_m-1}{\psi }_p& \cdots & L_{{\mathrm{\Omega }}_{\alpha }}{L}_f^{r_p-1}{\psi }_p\end{array}\right]\end{array}$$ E86

Finally, filtered system Eq. (5), as it is rewritten in Eq. (85), is in the same form of Eq. (74). Then, sparse recovery algorithm discussed in Section 5.4.1 is applied to Eq. (85) to reconstruct $d\left(t\right)$ .

Remark 4: The derivatives ${\dot{\mathrm{\Upsilon }}}_{r_1}^1,\dots ,{\dot{\mathrm{\Upsilon }}}_{r_p}^p$ are computed exactly in finite time using higher-order sliding mode differentiators [28] discussed in Eqs. (65) and (66).

6. Case study

Consider the mathematical models (1)–(4) of the US Western Electricity Coordinating Council (WECC) power system [8] with three generators and six buses (Figure 1) when the sensors of the generator speed deviations from synchronicity are under stealth attack and plant is under deception attack.

Assumption (A12): The matrix $L_{l,l}^{\theta }$ in (3) is nonsingular.

If (A12) holds, then the variable $\theta$ can be rewritten as

$$\theta ={\left(L_{l,l}^{\theta }\right)}^{-1}\left(-R_{l,g}^{\theta }\delta +P_{\theta }+B_{\theta }d\right)$$ E87

Substituting (87) into (1), then it follows that

$$\begin{array}{l}\left[\begin{array}{c}\dot{\delta }\\ \dot{\omega }\end{array}\right]=\left[\begin{array}{cc}0& I_{p\times p}\\ M_g^{-1}\left(-L_{g,g}^{\theta }+L_{g,l}^{\theta }{\left(L_{l,l}^{\theta }\right)}^{-1}{L}_{l,g}^{\theta }\right)& -M_g^{-1}{E}_g\end{array}\right]\left[\begin{array}{c}\delta \\ \omega \end{array}\right]+\left[\begin{array}{c}0\\ P_{\mathit{\theta \omega }}\end{array}\right]+\left[\begin{array}{c}{B}_{\delta }\\ B_{\mathit{\theta \omega }}\end{array}\right]d\left(t\right)\begin{array}{cc},& y=C\left[\begin{array}{c}\delta \\ \omega \end{array}\right]+\left[\begin{array}{c}{D}_{\delta }\\ D_{\omega }\end{array}\right]d\left(t\right)\end{array}\\ P_{\mathit{\theta \omega }}=M_g^{-1}\left(P_{\omega }-L_{g,l}^{\theta }{\left(L_{l,l}^{\theta }\right)}^{-1}{P}_{\theta }\right)\begin{array}{cc},& B_{\mathit{\theta \omega }}=M_g^{-1}\left(B_{\omega }-L_{g,l}^{\theta }{\left(L_{l,l}^{\theta }\right)}^{-1}{B}_{\theta }\right)\end{array}\end{array}$$ E88

6.1 Simulation setup

1. The three sensors of rotor angles, $\delta \in {\mathbb{R}}^3$ , are assumed protected from attack, but the three sensors of the generator speed deviations from synchronicity, $\omega \in {\mathbb{R}}^3$ , are assumed to be attacked.

2. The $B_{1\omega }=I_3,B_{1\theta }=0_{6\times 3},D_{\delta }=0_{3\times 6}$ are given, and then Eq. (88) is reduced to

$$\left\{\begin{array}{l}\begin{array}{l}\begin{array}{cc}\dot{\upsilon }={\phi }_{\delta }\left(\delta ,\omega \right),& \begin{array}{cccc}\begin{array}{cc}& \end{array}& & & \end{array}\end{array}\\ \dot{\omega }={\phi }_{\omega }\left(\delta ,\omega \right)+P_{\mathit{\theta \omega }}+M_g^{-1}{d}_x\left(t\right)\end{array}\\ \begin{array}{cc}{y}_1=C_1\upsilon ,& y_2=C_2\omega +D_{1\omega }{d}_y\left(t\right)\end{array}\end{array}\right),\mathit{\text{where}} C_1=C_2=I_{3\times 3},D_{\omega }=\left[\begin{array}{c}\begin{array}{cccc}\begin{array}{ccc}0& 1& 2\end{array}& 0& 1& 1\end{array}\\ \begin{array}{cccc}\begin{array}{ccc}1& 0& 0\end{array}& 2& 1& 0\end{array}\\ \begin{array}{cccc}\begin{array}{ccc}0& 0& 1\end{array}& 0& 1& 0\end{array}\end{array}\right]$$ E89

Remark 5: $D_{1\omega }$ satisfies RIP condition defined in Eq. (75).

In the first step of attack reconstruction, $d_x\left(t\right)$ is estimated by using protected measurements $y_1$ and the SMO described in Section 5.2. It is easy to verify that

$$\begin{array}{cc}\begin{array}{l}{\overline{C}}_{\delta 1}\overline{B}=0\begin{array}{cc},& \end{array}{\overline{C}}_{\delta 1}A\overline{B}\ne 0\\ {\overline{C}}_{\delta 2}\overline{B}=0\begin{array}{cc},& \end{array}{\overline{C}}_{\delta 2}A\overline{B}\ne 0\\ {\overline{C}}_{\delta 3}\overline{B}=0\begin{array}{cc},& \end{array}{\overline{C}}_{\delta 3}A\overline{B}\ne 0\end{array}& C_a=\left[\begin{array}{c}\begin{array}{c}{C}_1\\ C_{1}A\\ C_2\end{array}\\ C_{2}A\\ C_3\\ C_{3}A\end{array}\right]=\left[\begin{array}{cccccc}1& 0& 0& 0& 0& 0\\ 0& 0& 0& 1& 0& 0\\ 0& 1& 0& 0& 0& 0\\ 0& 0& 0& 0& 1& 0\\ 0& 0& 1& 0& 0& 0\\ 0& 0& 0& 0& 0& 1\end{array}\right],y_a=\left[\begin{array}{c}\begin{array}{c}{y}_1\\ \mu \left(y_1-{\widehat{y}}_1\right)\\ y_2\end{array}\\ \mu \left(y_2-{\widehat{y}}_2\right)\\ y_3\\ \mu \left(y_3-{\widehat{y}}_3\right)\end{array}\right]\end{array}$$ E90

where ${\overline{C}}_{\mathit{\delta i}}$ is the ith row of ${\overline{C}}_{\delta }$ . The states of the system, $\widehat{\delta },\widehat{\omega }$ , and plant attacks ${\widehat{d}}_x\left(t\right)$ are reconstructed using Eqs. (43) and (50). Then, $\widehat{\omega }$ is used in Eq. (89) to find

$$D_{\omega }{d}_y\left(t\right)=y_2-\widehat{\omega }$$ E91

There are six sources $d_{y1},\dots ,d_{y6}$ attacking three measurements ${\omega }_1,{\omega }_2,{\omega }_3$ , and at any time, just one out of six attack signals is nonzero. The SR algorithm in Section 5.2 is applied to find ${\widehat{d}}_y\left(t\right)$ . The following attacks are considered for simulation.

$$\begin{array}{l}\left[\begin{array}{c}{d}_{x1}\\ d_{x2}\\ d_{x3}\end{array}\right]=1\left(t-10\right).\left[\begin{array}{c}sin\left(0.5t\right)\\ 1\left(t\right)-1\left(t-4\right)+1\left(t-8.5\right)-1\left(t-13\right)+1\left(t-17.5\right)\\ cos\left(t\right)+0.5sin\left(3t\right)\end{array}\right],\\ d_y\left(t\right)=1\left(t-10\right).{\left[\begin{array}{cccc}\begin{array}{ccc}0& 0& 0\end{array}& 0& sin\left(t\right)& 0\end{array}\right]}^T\end{array}.$$ E92

Deception attacks $d_{x1}$ , $d_{x2}$ , and $d_{x3}$ are reconstructed very accurately as shown in Figures 4–6. The only nonzero sensor attack is detected and accurately estimated by using the SR algorithm as shown in Figure 7. In Figure 8a and 8b, the corrupted system outputs (which are system states in our case) are compared to the “cleaned” outputs that are computed by subtracting the estimated attacks from the corrupted sensors and actuators and to the system outputs when the system is not under attack.

Figure 4.

Plant attack $d_{x_1}$ compared to estimated ${\widehat{d}}_{x_1}$ .

Figure 5.

Plant attack $d_{x_2}$ compared to estimated ${\widehat{d}}_{x_2}$ .

Figure 6.

Plant attack $d_{x_3}$ compared to estimated ${\widehat{d}}_{x_3}$ .

Figure 7.

Sensor attack $d_y$ reconstruction.

Figure 8.

(a) Corrupted output $y_1,y_2,y_3$ compared with compensated and without any attack output and (b) corrupted output $y_4,y_5,y_6$ compared with compensated and without any attack output.

7. Conclusion

The critical infrastructures like power grid, water resources, etc. are large interconnected cyber-physical systems whose reliable operation depends critically on their cyber substructure. In this chapter, cyber-physical systems when their sensors and/or states are under attack or experiencing faults are investigated. The sensor and states/plant attacks are reconstructed online by using a fixed-gain and adaptive-gain sliding mode observers. As soon as the attacks are reconstructed, corrupted measurements and states are cleaned from attacks, and the control signal that uses cleaned measurements provides cyber-physical system performance close to the one without attack. The effectiveness of the proposed approach is shown by simulation results of a real electrical power network with sensors under stealth attack and states under deception attacks.