Stochastic Models to Mitigate Sparse Sensor Attacks in Continuous-Time Non-Linear Cyber-Physical Systems

1  Introduction

Cyber-Physical Systems (CPS) are seamless integrations of physical and computational processes [1]. Many different architectures and approaches to support these unions have been reported, from schemes based on the control theory [2] to feedback loops in computational systems [3]. But all proposed CPS implementations include a sensing platform to monitor the evolution of the physical world [1]. That platform is dense, including thousands of networked sensor nodes capable of capturing information through several different physical parameters [4]. Those data must be injected into computational processes, to ensure that the cybernetic and physical worlds evolve together in a feedback control loop [5].

Therefore, precise information about physical processes is essential to ensure that the loop is convergent and follows the expected evolution [6]. However, it is hardly possible to obtain precise information in real applications [4]. Many random effects have an impact on the behavior of CPS, such as noise, transmission errors, measurement, digitalization, and discretization processes [4]. In that way, information finally injected into computational processes is not the raw or authentic information acquired from the physical world, but a non-deterministic transformation of it. And this transformed information prevents the CPS from integrating the computational and physical processes with the expected synchronicity and showing the required behavior [7].

Furthermore, as Cyber-Physical Systems are used in more scenarios and applications, including critical infrastructures, they are more exposed to new risks. Eventual and unexpected cyberattacks are the main ones. Although innovative attack strategies have been reported to exploit specific vulnerabilities of CPS [8], nowadays the greatest risks for CPS are still associated with classic cyberattacks such as the Sparse Sensor Attack (SSA). In the SSA [9], attackers introduce false information and/or cause delays in the sensing platform monitoring the physical world at a low level, so the CPS behavior is altered or denied. It is the most common attack in control solutions, and new uncertainty about the physical information injected into the computational processes is to be handled.

In this context, reconstruction mechanisms to recover the original and real information extracted from the physical world are essential [10]. The state of any CPS may be described as a multidimensional vector, where each position represents a physical parameter. By establishing the analytical law that describes the trajectory of all those state variables in the phase space, the transformed information received may be corrected through a theoretically predicted CPS state [11]. However, in the general case, all physical parameters are not independent, but they are interrelated through complex physical laws [12]. The appearance of complex non-linear laws, together with the need for stochastic terms to describe random effects such as sparse sensor cyberattacks, turns quite difficult to find a general high-precision model. Thus, traditional reconstruction schemes are based on some basic assumptions, so the mathematical expressions describing the evolution of CPS are easier to manipulate and implement [13].

Our work is motivated by limitations and vulnerabilities caused by these simple assumptions, which make CPS weaker against cyberattacks than other state-of-the-art technological systems. Namely:

•   First, Cyber-Physical Systems are assumed to evolve according to a linear law.

•   Second, all terms are considered deterministic, including noise and attacking signals.

•   Third, all physical variables are assumed to be fully independent of each other.

•   And fourth, physical processes are assumed to be discrete, so digitalization and transmission processes do not have to be explicitly considered. Although those linear deterministic models present important advantages (for example, they can be manipulated to find analytical expressions for the detection and identification of SSA), their applicability is very limited [14].

•   Only closed CPS based on a reduced number of physical variables with a smooth and invariant behavior (such as the temperature in a climatized space) are governed and can be secured and protected by such a simple model.

Therefore, more complex and general models are required to protect and mitigate SSA in multidimensional CPS with a continuous-time non-linear behavior. In this paper, we address this challenge.

Three innovative contributions are introduced in this paper:

•   A complex non-linear model to describe the CPS behavior in a general situation.

•   New signals and models for SSA and digitalization processes.

•   The third and final contribution is an innovative reconstruction scheme.

The proposed model describes the behavior of CPS using unknown generic functions, which are developed as Taylor series. This model also includes stochastic terms to represent physical, transmission, and measurement noises. Besides, SS attacks are described as a new signal whose value follows a probabilistic behavior according to a given discrete random variable. Physical processes are represented by continuous-time signals that are discretized using an event-based scheme. The resulting multidimensional model injects discrete data into computational processes, but is too complex to generate analytic expressions to mitigate SSA in CPS. Finally, the proposed reconstruction scheme is supported by a weighted stochastic model where the error probability is estimated through different indicators commonly employed to describe non-linear dynamics (such as the Fourier transform, first-return maps, or the probability density function). A decision algorithm calculates the final reconstructed value considering the previous error probability.

The rest of the paper is organized as follows. Section 2 analyzes the state-of-the-art on cyberattacks and countermeasures in CPS. Section 3 describes the proposed solution, including the mathematical model to describe the behavior of the CPS and the reconstruction and protection scheme to mitigate SSA. Finally, Section 4 describes the experimental validation and the results obtained. Section 5 concludes the paper.

2  Related Works

Cyber-Physical Systems are one of the most promising technological revolutions nowadays. They are expected to govern all production, domestic, and critical digital systems. Due to this relevance, many authors have investigated how to protect CPS against various well-known and innovative attacks. In general, we can distinguish two different protection approaches: those based on control theory and those supported by Information Technologies (IT).

IT protection mechanisms for CPS are usually data processing and filtering modules to remove and correct malicious or corrupted data packets. Stochastic techniques and models [4], advanced filtering algorithms such as the Kalman filter [15], hardware-enabled algorithms such as parameter estimation [16], and pattern recognition techniques to identify unusual information [17] are the most common technologies. As well as game-theory and other common technologies for CPS protection, such as honeypots [18] or Software-Defined Networks [19]. However, a limited number of works supporting this vision may be found, as information theory techniques are high-level and agnostic concerning the underlying hardware platform [20]. And the most critical cyber risks in CPS nowadays are associated with sensor and actuator nodes [8]. Different authors have identified new attack vectors and strategies [8], so feedback loops in CPS can be used to magnify cyberattacks starting in a single hardware node and spreading throughout the entire system. Furthermore, these IT protection technologies are computationally heavy and require long processing times, so they are not effective against fast cyberattacks. Other low-level lightweight techniques are required.

Physical infrastructure protection is, then, a priority in CPS. And most works on CPS security employ control theory to design new hardware protection schemes. Globally, all these technologies follow the same strategy [21]: they estimate or predict a secure future state for the physical platform and/or control loop, which is used to mitigate different types of attacks. Although this paradigm could fully protect CPS [22], it is very difficult to implement in practice and the reported implementation presents different weaknesses. Techniques may be local (or decentralized), distributed, or centralized.

Decentralized state estimation techniques are handled by independent sensing nodes. They are sparse as individual sensors have very limited information and actuation capabilities, so the achieved protection level is poor. Continuous bidimensional linear models are employed to detect perturbations and attacks (typically Denial of Service attacks) and modify the behavior of nodes by, for example, increasing their computational resources [23]. The objective is to guarantee the local stability of the control loops by mitigating all perturbances [21]. In contrast, other decentralized CPS protection schemes use variance-based strategies (also known as ‘secure control’ [14]). This approach is more general and can be applied against a generic cyberattack. Using discrete bidimensional models, tuned filters and tuned control loops can be varied to reduce system errors, even while a cyberattack is running [24]. However, even if local control loops can operate normally, with variance-based techniques the global system is handling corrupted data, and that impacts the later global behavior. Some authors have shown that global system protection requires cooperation and information sharing among all agents [21]. Distributed techniques fill this gap.

Distributed secure state estimation is useful against systemic attacks such as Byzantine attacks [25]. System states are deducted through an optimization process where linear models represent the sensors’ outputs and graphs [26], Markov chains [27], binary decision trees [28], and other mathematical paradigms (such as the Lipschitz continuity) [29] are used to represent the interconnections and transmissions among nodes. Custom quasi-linear models for specific applications, such as series-parallel systems, have been also reported [30,31]. However, these protection mechanisms are passive and cannot deploy countermeasures to mitigate the impact of cyberattacks. Then, they must be complemented with specific controllers [32,33] to apply active protection policies on the CPS. Anyway, the final performance of distributed protection techniques is highly dependent on the number of trusted nodes, not affected by the attack [21,34]. Furthermore, linear and quasi-linear models cannot represent the output of most complex sensing platforms [35]. Thus, reported schemes can only be applied to a reduced number of application scenarios, excluding critical risks such as massive or viral attacks and common nonlinear algorithms.

On the other hand, recently distributed artificial intelligent solutions, such as federated learning [36], Support Vector Machines [37], feature selection [38] or eXplainable Artificial Intelligence (XAI) [39], have also been applied to CPS securitization and intrusion detection. But performance must be enhanced through additional techniques such as reinforcement learning [40]. Intelligent solutions must be designed for very specific attacks, as they are usually focused on Denial-of-Service attacks. Although the final results are promising, the balance between cost and performance is still worse than the one observed in other distributed techniques, and they are preferred to be used for privacy preservation [41].

The main disadvantage of distributed protection mechanisms is the increase in system congestion, due to the large number of transmissions required to run the distributed algorithms. On the contrary, centralized approaches may handle global stability and attacks (as distributed techniques) but with a lower system overload. Most reported works follow this paradigm.

Centralized control is usual in CPS, as it is the traditional approach in legacy Supervisory Control And Data Acquisition (SCADA) systems. Different kinds of multi-dimensional models are employed to represent the state of every single node on the platform. These models can be analytically manipulated to define protection algorithms based on Orthogonal-Triangular (QR) decomposition [42] or Linear–Quadratic (LQ) control [43], mitigating the impact of attacks. Models can be deterministic [44] or include some stochastic terms to represent noises [45]. Besides, continuous [44] and discrete [46] models may be found. However, most of these models are linear and only consider the self-maintained evolution of the node output and the measurement errors (in line with traditional control theory models). While other relevant effects, such as the digitalization process or the transmission protocols, are not considered, although they can be relevant. On the other hand, nonlinear models are very rare [47] and they are only developed for specific use cases. This centralized approach is successful against false information attacks (also known as sparse sensor attacks or deception attacks [14]), as it handles a full picture of the CPS. However, current models are very limited, and analytical protection algorithms have a reduced impact in real applications.

Table 1 summarizes the main current approaches and their associated open challenges.

In this paper, we address this challenge, with a continuous-time generic multidimensional non-linear model, and a protection policy based on probabilistic decision-making schemes.

3  Proposed Scheme

The proposed solution includes two different phases. First, a multidimensional stochastic model is employed to estimate or predict the future state of the CPS. Later, the obtained secure state estimation is compared to the real state produced by the physical platform. Both values are compared using a probabilistic model, where several indicators are considered. The system state may be replaced or corrected using the predicted secure state if the decision-making algorithm indicates the information is false (corrupted or caused by an SSA). This section describes in detail the entire scheme. Section 3.1 introduces the stochastic model, while Section 3.2 presents the decision-making and protection algorithms.

3.1 Secure State Estimation. Model Description

A CPS is supported by a dense sensing platform including N different sensor nodes nm. These nodes monitor and control a catalogue of P different physical variables xi(t) (1). Each variable is monitored in Ki different geographical locations gsxi (2), so every node controls a different physical variable in a different location (3). If any node nm monitors more than one variable, we are analyzing it as two independent nodes located in the same geographical position.

Hereinafter we are naming xis (t) the value of physical variable xi(t) in location gsxi.

{xi(t) i=1,…,P}(1)

{gsxi s=1,…,Ki}(2)

N=∑i=1PKi(3)

The information to be finally injected into the computational processes (or system state) is a set of M discrete state variables yj[k], related to the physical variables through a vector unknown function, 𝒮(⋅), named as “system function” (4). This system function integrates five different processes: (i) the physical world’s evolution, (ii) the transduction phase, (iii) the measurement scheme, (iv) the data transmission, and (v) the final processing stage.

{yj[k] j=1,…,M}=𝒮(xis (t) i=1,…,P s=1,…,Ki)(4)

The physical world (i) is considered a closed autonomous system, with no external intervention, so the future evolution of the physical variables is only determined by the past values of those same variables (5). The function relating the past and future values of the physical variables ℱxis(⋅) is unknown and, in the general case, non-linear. For clarity, we are using vector X→ to represent the full ordered collection of physical variables (6). Although it is unknown, vector function ℱ(⋅)  could be developed as Taylor’s series.

xis (t) t ≥t0=ℱxis({xis (t) t<t0 i=1,..,P s=1,…,Ki}) ∀ i,s(5)

X→ (t)={xis (t) i=1,..,P s=1,…,Ki} ∀ t(6)

Any multidimensional function may be developed as Taylor’s series using the partial derivation (7). For simplicity, we are using a McLaurin development around the origin. In this expression terms (1k1!⋅…⋅kN! ∂k1+…+kN∂(x11)k1⋅…⋅∂(xPKP)kN ℱxis(0→)) are unknown coefficients as function ℱxi(⋅) is unknown too. We are representing them as λr (8). The purpose of our model is to estimate the future CPS state, in order to mitigate any potential SSA. Then, the model must be numerically implementable. Infinite series are not, and they must be limited to the R first terms (9). The error Eℱ (10) we are introducing because of this truncation is difficult to estimate as function ℱxis(⋅) is unknown, but the Lagrange formula represents its analytical expression. Because this error is not easy to compute, the value for Rℱ parameter must be experimentally chosen, so the numerical model is precise enough to represent the behavior of CPS. However, in order to handle uncertainties in our model, we are proposing an estimation (maximum value) for error Eℱ (11). We are assuming function ℱxis grows exponentially (the maximum increasing speed) in all directions, so term (∂k1+…+kN∂(x11)k1⋅…⋅∂(xPKP)kN ℱxis(0→))  is the unit. Besides, we choose the maximum value for term (1k1!⋅…⋅kN!) which is achieved for k1=…=kN−1=1.

ℱxis(X→ (t))=∑k1=0∞…∑kN=0∞(1k1!⋅…⋅kN! ∂k1+…+kN∂(x11)k1⋅…⋅∂(xPKP)kN ℱxis(0→))⋅(x11)k1⋅…⋅(xPKP)kN(7)

ℱxis(X→ (t))=∑r=0k1+…+kN=r∞λr⋅(x11)k1⋅…⋅(xPKP)kN(8)

ℱxis(X→(t))≈∑r=0k1+…+kN=rRℱλr⋅(x11)k1⋅…⋅(xPKP)kN (9)

|Eℱ|≤max{λRℱ+1}=max{1k1!⋅…⋅kN! ∂k1+…+kN∂(x11)k1⋅…⋅∂(xPKP)kN ℱxis(0→) k1+…+kN=Rℱ+1}(10)

|Eℱ|≤1(Rℱ−N+2)!(11)

The second process represented by the system function 𝒮(⋅) is the transduction phase (ii). Physical variables xis(t) are transformed into electrical signals vis(t) through an unknown function 𝒯m which is different for each sensor node nm (12). Functions 𝒯m are unidimensional (scalar) as the transduction process must be bijective to preserve the information. Besides, two different kinds of multiplicative noises affect the transduction phase. On the one hand, multiplicative physical noises εrxis(t) (such as thermal noise or environmental radiation) are mixed with real physical variables xis(t) in the transformation function 𝒯m. On the other hand, multiplicative electrical noises ξrxis(t) are added to the obtained electrical signals, because of the impact of electronic circuits. Each noise (electrical or physical) has a similar probability distribution f[⋅] (13). Noises are white (Gaussian), mutually uncorrelated with zero mean and unitary variances [42]. We are considering all these stochastic processes are stationary, and their probability distribution remains stable along time (14). Parameters Rvis1 and Rvis2 are positive integer numbers.

Unknown function 𝒯m can be developed as Taylor’s series as well, but in this case, expressions are simpler are unidimensional techniques can be applied (15). As done before, unknown coefficients (1r! dr𝒯md(x~is)r(0)) are represented by variables αr. Besides, Taylor’s series must be truncated too, to make our model computationally handleable (16), although an error E𝒯 is introduced (17). An estimation for the maximum value of error E𝒯 is proposed too (18), in order to enable the uncertainty management.

vis(t)=𝒯m(xis(t)+∑r=1Rvis1εrxis(t))+∑r=1Rvis2ξrxis(t)+∑r=1Rvis3mrxis(t)(12)

P(a≤ε≤b)=∫abfε(u) du(13)

P(ξrxis(t)=u) ∼ P(εrxis(t)=u)∼ 12π e−u22 ∀t(14)

𝒯m(x~is)=∑r=0∞1r! dr𝒯md(x~is)r(0)⋅(x~is)rbeing x~is(t)=xis(t)+∑r=1Rvis1εrxis(t)(15)

𝒯m(x~is)≈ ∑r=0R𝒯αr⋅(x~is)r(16)

E𝒯≤αR𝒯+1=1(R𝒯+1)! dR𝒯+1𝒯md(x~is)R𝒯+1(0)(17)

E𝒯 ≤ 1(R𝒯+1)!(18)

Although white noises εrxis are affected by functions 𝒯m and then, they are part of the Taylor’s series (15)–(16) because of the fact they follow a Gaussian distribution, these expressions can be simplified, so only the physical variables are part of the Taylor’s polynomial. Every term in the Taylor’s series where a noise εrxis is included may be considered as a non-monotonous transformation T(⋅) of a Gaussian random variable. The transformation theorem (19) shows that any transformed Gaussian distribution is a new Gaussian distribution with mean μt and variance σt (20). Later, all the transformed Gaussian distributions may be aggregated, and, because of the central limit theorem, the resulting random variable χxis is a Gaussian distribution too. However, mean μtt and variance σtt (21) are unknown, as the final value for the mean and variance of the global distribution χxis depends on the transformations and the value of R𝒯 parameter. Taylor’s series for function 𝒯m may be finally rewritten (22).

f𝒯m(εrxis)(u)=∑ufεrxis(u) |ddεrxis T(εrxis)|ε=u|being {u} the roots of T(εrxis) (19)

P(𝒯m(εrxis)=u)∼ 1σt2π e−(u−μt)22σt2(20)

P(χxis(t)=u)∼ 1σtt2π e−(u−μtt)22σtt2 ∀t(21)

𝒯m(x~is)≈ 𝒯m(xis)=∑r=0R𝒯αr⋅(xis)r+χxis(22)

Finally, it is necessary to estimate the value for mean μtt and variance σtt, so error in the proposed model may be properly handled. In general, errors are bigger as values for the mean μt and variance σt go up. Then, a superior limit is a good approximation for both parameters (23), which may be easily calculated considering the reproducibility of the Gaussian random variables. In order to get the final values for the mean μtt and variance σtt it is enough to apply the same reproducibility law a second time (24).

f𝒯m(εrxis)(u)=∑ufεrxis(u) |ddεrxis T(εrxis)|ε=u| ≤ ∑ufεrxis(u) ∼12π⋅R𝒯 e−(u−∑r=1R𝒯(xis)r)22R𝒯(23)

σtt=R𝒯⋅Rvis1μtt=Rvis1⋅∑r=1R𝒯(xis)r(24)

The transduction phase is open, so it can be affected by SSA and malicious signals. In order to represent this risk, we are considering a set of additive malicious signals mrxis(t) whose value is determined by a stochastic process. This stochastic process is characterized by a Bernoulli distribution Γa, representing the existence of a running SSA. Parameter a is equal to the unit if the attack is running, or zero in the opposite case (25). The attack probability ρa varies with time (as Γa is a stochastic process). The estimation scheme for this probability is part of the attack detection algorithm (see Section 3.2). In our model, a SSA consists of adding a false data signal zrxis(t) to the data electrical signal vis(t), according to the previously described distribution (26). Rvis3 different uncorrelated attackers may be operating over the CPS at the same time. False data signals, in our model, are understood as unreported and unexpected perturbations. This is relevant in order to define a precise attack detection and mitigation strategy (see Section 3.2).

Γa(t) : P(a;t)={1−ρa(t) if a=0 ρa(t) if a=1(25)

mrxis=arxis⋅zrxis(t)(26)

Finally, as every sensor node nm has a different function 𝒯m, the whole CPS is represented by a set of N different Eq. (27), which can be represented in one vector expression (28).

{v11(t)=∑r=0RJαrx11⋅(xis)r+𝒳x11+∑r=1Rv112ξrx11(t)+∑r=1Rv113mrx11(t)⋯vpkp(t)=∑r=0RJarxpkp⋅(xis)+𝒳xpkp+∑r=1Rvp2kpξrxpkp(t)+∑r=1Rvp3kpmrxpkp(t)(27)

V→(t)=∑r=0R𝒯αr→⋅(X→ (t))r+χ→+∑r=1Rv112ξr→(t)+∑r=1Rv113mr→(t)(28)

The third subprocess to be represented in our model is the measurement scheme (iii). This, basically, is a digitalization scheme, developed internally by sensor nodes (see Fig. 1). Discrete signal dis[k] are obtained through an ideal sampling scheme (29), where electrical signals are multiplied by a Dirac comb or impulse train ωTm(t) with period Tm (30). This period Tm  is different for each sensor node nm.

dis(t)=vis(t)⋅ωTm(t)+qis(t)=∑k=−∞∞vis(k⋅Tm)⋅δ(t−k⋅Tm)+qis(k⋅Tm)dis[k]=vis(k⋅Tm)+qis(k⋅Tm) k∈N(29)

ωTm(t)=∑k=−∞∞δ(t−k⋅Tm)(30)

Figure 1: Ideal sampling scheme

In this digitalization process only the quantification noise qis(t) is relevant. This noise, as the digitalization scheme is invariant in time, is also time-invariant, and characterized by a uniform random variable (31).

P(qis(t)=u )∼{1Δm u ∈[−Δm, Δm] 0 otherwise(31)

where Δm is the quantification step, fixed for every node nm.

The fourth process to be represented is the data transmission (iv). In general, hardware platforms in CPS are low-energy, and they sleep most of the time. Being event-based, they only activate the transmission subsystem when an event is detected in the physical world. We are defining function ϕm[k] as event-triggering function (32). This function takes as value the unit in the discrete time instant a new event must be generated. Its value is zero otherwise. Function u(⋅) is the Heaviside step function, ke is the last time instant where an event took place, and em is a parameter, different for each node nm. If signal dis changes its value more than em units, a new event is triggered.

ϕm[k]=u(dis[k]dis[ke]−(em)2)(32)

Data transmission is, once again, an open process, so it is vulnerable to attacks. Denial-of-Service (DoS) attacks in this case. But SSA too (as the transduction phase). Bernoulli distribution Γb represents the probability of a DoS attack to be running. Parameter b is equal to the unit if an attack is being performed, or zero if not (33). The attack probability ρb varies with time (as Γb is a stochastic process). The estimation scheme for this probability is part of the attack detection algorithm (see Section 3.2). Similarly, Bernoulli distribution Γc represents the probability of a SSA to be running at the data transmission stage (34).

Γb[k] : P[b;k]={1−ρb[k] if b=0 ρb[k] if b=1(33)

Γc[k] : P[c;k]={1−ρc[k] if c=0 ρc[k] if c=1(34)

All parameters and their meaning are equivalent to distributions Γa and Γb.

In our model, a DoS attack is represented by an arbitrary delay of kd units in the data transmission, while an SSA is represented by injected false signals his[k] (35). Then, the received signal by the remote central control platform wis[k] depends on function ϕm[k] and distributions Γb and Γc, but also is affected by transmission errors and noises.

wis[k]=c⋅his[k]+(1−c)⋅(b⋅dis[k−kd]+(1−b)⋅dis[k])+∑r=1Rwisφrwis[k]being k ⋮ ϕm[k]=1(35)

Our model considers Rwis multiplicative white Gaussian uncorrelated noises, φrwis with zero mean and unitary variance, affecting the data transmission.

Finally, information injected into computational processes yj[k] may not be the raw physical information, but a transformation of it (mean, minimum or maximum values, for example). Then, the final step in the system function 𝒮(⋅) is the processing stage (v). Processing processes may combine different transmitted signals wis[k] according to function 𝒫yj(⋅) which, in general, is unknown and, even, may change with time (36). This is an internal process where only numerical errors may affect the final result. However, central control systems are usually computationally powerful, and numerical error are negligible.

yj[k]=𝒫yj(W [k]→)=𝒫yj(wis[k] i=1,…,P s=1,…,Ki)(36)

Function 𝒫yj(⋅) may be developed as Taylor’s series, in a similar way as done for function ℱxis. Considering unknown coefficients βr, the final equation of our model may be described as a polynomial (37). Introducing an error E𝒫, whose maximum value may be estimated using the same techniques described before (38), the model may be truncated and limited to R𝒫 terms (39) so it can be managed by computational infrastructures.

𝒫yj(W [k]→)=∑r=0k1+⋯+kN=r∞βr⋅(w11)k1⋅…⋅(wPKP)kNbeing βr=(1k1!⋅…⋅kN! ∂k1+…+kN∂(x11)k1⋅…⋅∂(xPKP)kN 𝒫yj(0→))(37)

|Eℱ|≤max{βR𝒫+1} ≤1(R𝒫−N+2)!(38)

𝒫yj(W [k]→)≈ ∑r=0k1+…+kN=rR𝒫βr⋅(w11)k1⋅…⋅(wPKP)kN(39)

Then, the final analytical model to describe the behavior of CPS includes five different Eq. (40). All parameters and coefficients are known (or may be estimated) but λr, αr and βr which must be calculated. The value for those parameters is obtained from an initial calibration process and an optimization algorithm based on the minimization of the Mean Square Error (MSE).

{xis(t)=∑r=0k1+⋯+kN=rRℱλr⋅(x11)k1⋅…⋅(xPKP)kN     t ≥t0V→(t)=∑r=0R𝒯αr→⋅(X→ (t))r+χ→+∑r=1Rv112ξr→(t)+∑r=1Rv113mr→(t)dis[k]=vis(k⋅Tm)+qis(k⋅Tm)    k∈Nwis[k]=c⋅his[k]+(1−c)⋅(b⋅dis[k−kd]+(1−b)⋅dis[k])+∑r=1Rwisφrwis[k]yj[k]=∑r=0k1+…+kN=rR𝒫βr⋅(w11)k1⋅…⋅(wPKP)kN(40)

3.2 Reconstruction and Protection Mechanisms

The proposed model (see Section 3.1) considers seven sources of perturbations. On the one hand, errors may be caused by four different phenomena: erratic behaviors in the physical variables, electrical noises, quantification noise, and transmission perturbations. And, on the other hand, three different potential attacks affect CPS in the general case: SSA at the transduction phase, and SSA and Denial-of-Service attacks at the transmission phase. Fig. 2 represents the proposed reconstruction and protection mechanisms.

Figure 2: Protection and reconstruction mechanism

As a novelty, the proposed reconstruction and protection mechanism evaluates all these potential perturbations to build a global stochastic process (contrary to traditional deterministic models). This stochastic process B[p, k] (41) is discrete. p is a discrete variable representing the four possible situations a CPS state may achieve: unperturbed (p=0), noisy (p=1), SSA-attacked (p=2) and DoS-attacked (p=3). While k is a variable representing the discrete time. Similarly, we can define M stochastic sub-processes Bj[p, k] for each one of the M state variables yj considered in the CPS.

B[p, k] → {γ0[k]=B[0, k] γ1[k]=B[1, k]γ2[k]=B[2, k]γ3[k]=B[3, k](41)

In the proposed protection mechanism, five indicators are employed to evaluate the probability distribution of the stochastic process B[p, k] at every time instant k: ① the probability density function, ② the Short-Time Fourier transform, ③ the first-return map, ④ the autocorrelation and ⑤ the first order forward difference. To evaluate all these indicators, the protection algorithm operates with two numerical series. YRrj (42) represents the series of the last Rr reconstructed states for the j-th state variable, and Y~Runj (43) represents the series of the last Run unreconstructed states for the j-th state variable (being k the current time instant).

YRrj={yj[k−r−Run] r=1,…,Rr}(42)

Y~Runj={yj~[k−r] r=1,…,Run}(43)

In order to make feasible the calculation of all these indicators from time series YRrj and Y~Runj, in this work we propose a specifically tailored definition. The calculation process for every indicator is described below.

Regarding the probability density function ①, the probability μj1 of any unreconstructed state yj~[k] for the j-th state variable to happen in a given CPS, may be evaluated considering the previous reconstructed states YRrj achieved by that CPS and the Laplace definition of probability (44), and being δ[⋅] the Kronecker’s delta function. The probability μpdfj for the entire Y~Runj series of unreconstructed states may be obtained as the mean value of all the individual probabilities (45). And, finally, the global probability μpdf for all the M state variables may be calculated as the average value (46).

μj1=∑r=1Rrδ[yj[k−r−Run]−yj~[k]]Rr(44)

μpdfj=∑m=1Run∑r=1Rrδ[yj[k−r−Run]−yj~[k−m]]Run⋅Rr(45)

μpdf=1M∑j=1Mμpdfj(46)

But even if the unreconstructed CPS state has a relevant probability, it can still be manipulated and not be coherent with the system evolution. This situation may be detected through two different indicators: the Short-Time Fourier Transform (STFT) and the first-return map. Considering the Short-Time Fourier Transform (STFT) ②, the Fourier spectrum tends to be stable in a CPS, so any abrupt change may indicate an attack. The STFT (47) is equivalent to the traditional Fourier transform, but only considering a limited number of samples (instead of the usual infinite sum) through a window function Ω[k, Rsam], typically the Hann (Hanning) window (48) with a width of Rsam samples. Then, the STFT 𝒴Rrj for the reconstructed states YRrj (49) may be calculated using a numerical algorithm, as well as the STFT 𝒴~Runj for the reconstructed states Y~Runj (50).

STFT{y[k]}=𝒴(m,ν)=∑k=−∞∞y[k]⋅Ω[k−m, Rsam]⋅e−jkν(47)

Ω[k, Rsam]=12−12cos(2πkRsam)(48)

𝒴Rrj=STFT{YRrj}=∑k=−∞∞yj[k]⋅Ω[k−Rr2−Run, Rr]⋅e−jkν(49)

𝒴~Runj=STFT{Y~Runj}=∑k=−∞∞y~j[k]⋅Ω[k−Run2, Run]⋅e−jkν(50)

Then, using the Euclidean definition for distance, we can analyze how different 𝒴Rrj and 𝒴~Runj are (51). As the distance μFouj gets bigger, the probability of unreconstructed states Y~Runj to be manipulated increases. As done before, the global distance μFou for all the M state variables may be calculated as the average value of all partial distances μFouj (52).

μFouj=∥𝒴Rrj−𝒴~Runj∥=(𝒴Rrj−𝒴~Runj)⋅(𝒴Rrj−𝒴~Runj)(51)

μFou=1M∑j=1MμFouj(52)

Another indicator we can use to identify situations where the unreconstructed CPS state is manipulated is the first-return map ③. The first return map is a function Π(⋅), which can be obtained numerically, and shows the relation between consecutive (reconstructed) CPS states (53). The minimum Euclidean distance μj2 between every ordered pair of unreconstructed states π[m] (54) and the first return map Π(⋅) represents how close the unreconstructed states are to the expected behavior (55). The global distance μrtj for the entire Y~Runj series may be obtained as the mean value (56), and the global distance μrt for all the M state variables may be calculated as the average value of all partial distances μrtj (57).

yj[k+1]=Π(yj[k])(53)

π[m]=(y~j[k−m], y~j[k−m+1])(54)

μj2=minη(r) ∈ Π(⋅)|| π[m]−η(r)||=minr ∈ (2, Rr)||π[m]−(yj[k−Run−r], yj[k−Run−r+1])||(55)

μrtj=1Run∑m=2Run(minr ∈ (2, Rr)||π[m]−(yj[k−Run−r], yj[k−Run−r+1])||)(56)

μrt=1M∑j=1Mμrtj(57)

But in some situations, very noisy states are difficult to distinguish from attacks. To clarify and separate these two situations we use the autocorrelation ④. Noise is a random effect, so autocorrelation tend to the null value very quickly. While planned attacks follow a certain structure, and autocorrelation oscillates but not disappears because of these patterns. But autocorrelation cannot be directly applied to series Y~Runj or YRrj, as they contain actual information, and it would be always non-null. Then, before calculating the autocorrelation we are using a stop-band filter to remove the legitimate information (see Fig. 3).

Figure 3: Stop-band filtering for autocorrelation calculation

From the STFT yRrj we can obtain the central frequency Ω0 and the bandwidth Ωc of the j-th information signal (state variable). And then, the stop-band filter in the Laplace domain may be described as a quotient function (58). And the filtering process as a product (59). The resulting filtered signal y¯¯j[k] (60), thus, only contains information about the noise and/or attacks affecting the CPS. The autocorrelation μauj[r] can be now obtained (61).

H(s)=s+Ω02s2+Ωc⋅s+Ω02(58)

𝒴¯¯Runj=𝒴~Runj⋅H(s)(59)

y¯¯J[k]=STFT−1{𝒴¯¯Runj}(60)

μauj[r]=∑m=1Run(y¯¯j[k−m]−ϖj)⋅(y¯¯j[k−m+r]−ϖj)∑m=1Run(y¯¯j[k−m]−ϖj)2 r∈[0,Run2]ϖj=1Run∑r=1Runy¯¯j[k−r] (61)

This autocorrelation μauj[r] should disappear as r parameter increases if the CPS is just noisy. To get that confirmation but avoid possible transitory effects, we are aggregating the last Rcor samples in the autocorrelation function μauj[r] (62). The resulting indicator μauj will be lower as the perturbations in the unreconstructed state are more similar to Gaussian white noise. As in all the previous indicators, the global autocorrelation μau for all the M state variables may be calculated as the average value of all partial distances μauj (63).

μauj=∑r=Run2−RcorRun2μauj[r](62)

μau=1M∑j=1Mμauj(63)

However, some attacks may use perturbations within the information signals’ bandwidth, and autocorrelation may not generate a conclusive result. To analyze this situation, we are using our last indicator, the first order forward difference ⑤. The first order forward difference μdiffj[k] (60) represents the tendency, evolution or growing of the j-th unreconstructed state variable. In general, CPS states fluctuate but do not increase or decrease in a monotonous manner. Even less if the evolution is divergent (for example, exponential). Then, the sum μdiffj  of all (Run−1) samples in the first order forward difference μdiffj[m] is typically very small (61), as growing periods are cancelled by decreasing phases and vice versa. But if the CPS state is manipulated and it does not oscillate but increases or decreases monotonously and diverges, the sum μdiffj  will take very extreme values (positive or negative). The global tendency (aggregated first order differences) μdiff for all the M state variables may be calculated as the average value (62).

μdiffj[m]=y~j[m+1]−y~j[m] m ∈[k−2, k−Run](64)

μdiffj=∑m=k−2k−Runμdiffj[m](65)

μdiff=1M∑j=1Mμdiffj(66)

Using these five indicators, we can now estimate the probability distribution of the stochastic process B[p, k]. Mathematical models for all this probability distribution are a genuine contribution of this work. The unperturbed state (p=0) is only probable when the probability μpdf is very high (close to the unit), and distances μFou and μrt are very low. Any other value is indicating a noisy state (p=1), which is still probable even for smaller values of probability μpdf and bigger values of distances μFou and μrt. But noisy states require a low value for autocorrelation μau to be probable (on the contrary, the CPS may be under a cyberattack). Because of this sensitivity, exponential and power laws are the most adequate ones to represent the probability of the unperturbed state γ0[k] (63), while linear evolutions and slower exponential laws fit the more tolerant behavior of the probability law γ1[k] of the noisy state (64).

γ0[k]=(1−(μpdf−1)2τ01)⋅exp(−μFouτ02)⋅exp(−μrtτ03)(67)

γ1[k]=μpdf⋅exp(−μFouτ12)⋅exp(−μrtτ13)⋅exp(−μauτ14)(68)