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

Cyber-Physical Systems are very vulnerable to sparse sensor attacks. But current protection mechanisms employ linear and deterministic models which cannot detect attacks precisely. Therefore, in this paper, we propose a new non-linear generalized model to describe Cyber-Physical Systems. This model includes unknown multivariable discrete and continuous-time functions and different multiplicative noises to represent the evolution of physical processes and random effects in the physical and computational worlds. Besides, the digitalization stage in hardware devices is represented too. Attackers and most critical sparse sensor attacks are described through a stochastic process. The reconstruction and protection mechanisms are based on a weighted stochastic model. Error probability in data samples is estimated through different indicators commonly employed in 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. An experimental validation based on simulation tools and real deployments is also carried out. Both, the new technology performance and scalability are studied. Results prove that the proposed solution protects Cyber-Physical Systems against up to 92% of attacks and perturbations, with a computational delay below 2.5 s. The proposed model shows a linear complexity, as recursive or iterative structures are not employed, just algebraic and probabilistic functions. In conclusion, the new model and reconstruction mechanism can protect successfully Cyber-Physical Systems against sparse sensor attacks, even in dense or pervasive deployments and scenarios


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 highprecision 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.

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 seriesparallel 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 selfmaintained 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.Numerical data filtering (Kalman) [16] Parameter estimation for numerical models and time series [17] Pattern recognition for numerical models and time series [18] Honeypots for attacker capture using the game theory [19] Software-Defined Network technologies to create intrusion detection systems [21,24,23] Decentralized In this paper, we address this challenge, with a continuous-time generic multidimensional nonlinear model, and a protection policy based on probabilistic decision-making schemes.

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.

Secure State Estimation. Model Description
A CPS is supported by a dense sensing platform including N different sensor nodes n m .These nodes monitor and control a catalogue of P different physical variables x i (t) (1).Each variable is monitored in K i different geographical locations g x i s (2), so every node controls a different physical variable in a different location (3).If any node n m monitors more than one variable, we are analyzing it as two independent nodes located in the same geographical position.
Hereinafter we are naming x s i (t) the value of physical variable The information to be finally injected into the computational processes (or system state) is a set of M discrete state variables y j [k], related to the physical variables through a vector unknown function, S (•), 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.
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 F x s i (•) 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 F (•) could be developed as Taylor's series.
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 are unknown coefficients as function 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 F (10) we are introducing because of this truncation is difficult to estimate as function F x s i (•) is unknown, but the Lagrange formula represents its analytical expression.Because this error is not easy to compute, the value for R F 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 F (11).We are assuming function F x s i grows exponentially (the maximum increasing speed) in all directions, so term is the unit.Besides, we choose the maximum value for term which is achieved for The second process represented by the system function S (•) is the transduction phase (ii).Physical variables x s i (t) are transformed into electrical signals v s i (t) through an unknown function T m which is different for each sensor node n m (12).Functions T 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 ε   (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 R 1 are positive integer numbers.
Unknown function T 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 1 r! d r T m d x s i 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 T is introduced (17).An estimation for the maximum value of error E T is proposed too (18), in order to enable the uncertainty management.
Although white noises ε x s i r are affected by functions T 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 ε x s i r 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 χ x s i 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 χ x s i depends on the transformations and the value of R T parameter.Taylor's series for function T m may be finally rewritten (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).
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 m x s i r (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 z x s i r (t) to the data electrical signal v s i (t), according to the previously described distribution (26).R 3 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).
Finally, as every sensor node n m has a different function T m , the whole CPS is represented by a set of N different Eq. ( 27), which can be represented in one vector expression (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 d s i [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 T m (30).This period T m is different for each sensor node n m .
Figure 1: Ideal sampling scheme In this digitalization process only the quantification noise q s i (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).
where m is the quantification step, fixed for every node n m .
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, k e is the last time instant where an event took place, and e m is a parameter, different for each node n m .If signal d s i changes its value more than e m units, a new event is triggered.
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).
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 k d units in the data transmission, while an SSA is represented by injected false signals h s i [k] (35).Then, the received signal by the remote central control platform w s i [k] depends on function φ m [k] and distributions b and c , but also is affected by transmission errors and noises.
Our model considers R w s i multiplicative white Gaussian uncorrelated noises, ϕ w s i r with zero mean and unitary variance, affecting the data transmission.
Finally, information injected into computational processes y j [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 S (•) is the processing stage (v).Processing processes may combine different transmitted signals w s i [k] according to function P y j (•) 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.
Function P y j (•) may be developed as Taylor's series, in a similar way as done for function F x s i .Considering unknown coefficients β r , the final equation of our model may be described as a polynomial (37).Introducing an error E P , whose maximum value may be estimated using the same techniques described before (38), the model may be truncated and limited to R P terms (39) so it can be managed by computational infrastructures.
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).

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.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 B j [p, k] for each one of the M state variables y j considered in the CPS.
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: 1 the probability density function, 2 the Short-Time Fourier transform, 3 the first-return map, 4 the autocorrelation and 5 the first order forward difference.To evaluate all these indicators, the protection algorithm operates with two numerical series.Y j Rr (42) represents the series of the last R r reconstructed states for the j-th state variable, and Y j Run (43) represents the series of the last R un unreconstructed states for the j-th state variable (being k the current time instant).
In order to make feasible the calculation of all these indicators from time series Y j Rr and Y j Run , in this work we propose a specifically tailored definition.The calculation process for every indicator is described below.
Regarding the probability density function 1 , the probability μ 1 j of any unreconstructed state y j [k] for the j-th state variable to happen in a given CPS, may be evaluated considering the previous reconstructed states Y j  Rr achieved by that CPS and the Laplace definition of probability (44), and being δ [•] the Kronecker's delta function.The probability μ j pdf for the entire Y j Run 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).
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) 2 , 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, R sam ], typically the Hann (Hanning) window (48) with a width of R sam samples.Then, the STFT Y j Rr for the reconstructed states Y j Rr (49) may be calculated using a numerical algorithm, as well as the STFT Y j Run for the reconstructed states Y j Run (50).

STFT {y
Then, using the Euclidean definition for distance, we can analyze how different Y j Rr and Y j Run are (51).As the distance μ j Fou gets bigger, the probability of unreconstructed states Y j Run 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 μ j Fou (52).
Another indicator we can use to identify situations where the unreconstructed CPS state is manipulated is the first-return map 3 .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 μ 2 j 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 μ j rt for the entire Y j Run 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 μ j rt (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 4 .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 j Run or Y j Rr , 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 y j Rr 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 μ j au [r] can be now obtained (61).
This autocorrelation μ j au [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 R cor samples in the autocorrelation function μ j au [r] (62).The resulting indicator μ j au 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 μ j au (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 5 .The first order forward difference μ j diff [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 μ j diff of all (R un − 1) samples in the first order forward difference μ j diff [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 μ j diff 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).
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).
On the other hand, SSA-attacked state (p = 2) is characterized by very a low probability μ pdf but very high distances μ Fou and μ rt .As well as a relevant non-null value in the aggregated autocorrelation μ au and the aggregated first order forward differences μ diff .On the contrary, DoS-attacked states (p = 3) are usually associated to moderate values for the probability μ pdf (states are delayed but not manipulated) while still very high distances μ Fou and μ rt (as they are delayed, states are not coherent with the historical series).The aggregated autocorrelation μ au and the aggregated first order forward differences μ diff tend also to be quite reduced.Following a similar philosophy to employed before, we can define the evolution laws for the probabilities γ 2 [k] (65) and γ 3 [k] (66).
In an equivalent manner we may calculate the probability distribution for all stochastic subprocesses B j [p, k] (67) Based on this stochastic process B [p, k], and all subprocesses B j [p, k], we propose a decision function with different thresholds to identify and trigger the proper protection and/or reconstruction mechanism at every time instant k.At this point we are also considering the series Ŷ j Run of predicted states (68), according to the proposed model (see Section 3.1).Time instants are exactly the same to the ones observed in series Y j Run of unreconstructed states.For the calculation of this series of estimated states, probabilities ρ a , ρ b and ρ c are obtained probabilities γ 2 and γ 3 (69).
Fig. 4 shows the proposed decision algorithm.This is an original contribution firstly presented in this work.In the first step it is evaluated if any global probability γ i is θ init units higher than any other probability (70).If that is the case, the situation represented by that probability γ i is considered to be the actual situation of the last received unreconstructed states Y j Run ∀j.If γ 0 is the highest probability, states are unperturbed, and they are added with no modification to the series of reconstructed secure states (71).If γ 1 is the highest probability, states are noisy.The reconstruction action depends on how noisy the unreconstructed states are (72).If the Mean Square Error (MSE) between series Y j Run and Ŷ j Run (73) is lower than threshold θ low noise , noise is negligible and unreconstructed states Y j Run are added with no modification to the series of reconstructed secure states.On the contrary, if the MSE is higher than threshold θ high noise , noise is considered too invasive and next R un reconstructed secure states are taken from the predicted series Ŷ j Run .In any other situation, noise is relevant but not dominant, and reconstructed states are calculated as the average between unreconstructed Y j Run and predicted Ŷ j Run series.For MSE calculation, predicted values Ŷ j Run are obtained considering all possible perturbation sources (for examples, parameters a, b and c takes the most probable value).Finally, if probability γ 2 or probability γ 3 is the highest, the CPS is under an attack (SSA or DoS respectively).In both circumstances, unreconstructed states are not secure and next R un reconstructed secure states are taken from the predicted series Ŷ j Run (74).In this last situation, predicted values are obtained in absent of attacks of any kind (i.e., a, b and c are null).If no global probability γ i is θ init units higher than any other probability, the same algorithm described above is applied to every j −th state variable.If any probability γ j i is θ j init units higher than any other, this is considered to be the actual situation in the CPS for this state variable (75).The next steps in the algorithm are equivalent to the description above, just using specific thresholds θ j,low noise and θ j,high noise for the situation when γ j 2 is the dominant probability (76).The objective, in this case, is to reconstruct the CPS state, variable by variable.This approach is slower and computationally more costly, so it is only triggered when the global analysis is not conclusive.
If no probability γ j i is θ j init units higher than any other for any j − th state variable, the stochastic process B [p, k] is not precise enough.Then, all the algorithm and calculations are repeated for larger values of R r and R un sizes.Then, results may be more precise when operating with more samples.But, if the proper reconstruction actions could not be selected before the maximum values for R r and R un sizes are reached, the global algorithm is run one last time.In this last case, the situation represented by the highest probability γ i (with no restriction) determines the reconstruction action, according to the algorithm described before.In any case, R r and R un sizes are always returned to the initial values.
Sizes for parameters R un and R r are actually very important and sensible.Large values for those sizes avoid most spurious numerical and transitory effects, but they reduce the precision and sensitivity of the protection and reconstruction algorithm to detect short-term attacks and high-frequency noises.While reduced values for parameters R un and R r behave exactly the opposite.The balance cannot be generalized and therefore must be found for every specific application.

Experimental Validation
To validate the proposed mechanisms for the protection of Cyber-Physical Systems against Sparse Sensor and Denial of Service attacks, an experimental validation was conducted.Section 4.1 describes the experimental methodology, while Section 4.2 presents the obtained results.

Experimental Methodology and Environment
The experiments were based on an emulated industrial scenario with real hardware devices (microcontrollers).The experimental works were divided into two different phases.First, we focused on analyzing the precision and attack detection capacity of the proposed technology.The second phase focused on studying the performance and scalability of the proposed model and the reconstruction and protection mechanism.
For all the experiments, the proposed CPS was supported by a collection of ESP-32 microcontrollers.Its number is variable depending on the experiment.ESP-32 microcontrollers are low-cost System-on-Chip provided with Wireless Fidelity (WiFi) and Bluetooth capabilities.It is based on a Tensilica Xtensa LX6 processor, and it includes several peripheral interfaces (Universal Asynchronous Receiver-Transmitter-UART-, Pulse Width Modulation-PWM-, Serial Peripheral Interface-SPI-, etc.), so it can handle a large catalog of different sensors.In our experiments, each ESP-32 node was provided with two sensors, monitoring four physical variables in total.The first sensor was a CCS811 sensor to monitor air quality.It can provide two different variables: carbon dioxide equivalent (eCO2) and organic volatile compounds concentration (TVOC).The second sensor is a DTH-11 device, which generates measurements for the environmental humidity and temperature.The measurement periodicity is variable and depends on the experiment.
All these sensors employed a WiFi connection to send all the collected information to a cloud server, located within the same building.Hypertext Transfer Protocol (HTTP) messages and Representational State Transfer (REST) interfaces were employed to support these communications.The server was a Linux-based machine (Ubuntu 18.04 LTS) with the following hardware characteristics: Dell R540 Rack 2U, 96 GB RAM, two processors Intel Xeon Silver 4114 2.2G, HD 2 TB SATA 7,2K rpm.In this server, both the proposed model and the reconstruction and protection algorithm were hosted and executed.A Node.js server was deployed to collect all data from the sensor nodes and send them to the computational process executing our proposal.A supervisory process was continuously evaluating the evolution and performance of the proposed algorithms and model.The acquired information was employed to carry out a statistical analysis using the MATLAB 2022a software, to validate our hypotheses.All experiments were repeated twelve times to remove possible spurious effects.The results for every measurement are obtained as the average of all these individual twelve realizations.
In the first phase, we performed two different experiments.The first experiment was aimed at analyzing the precision of the proposed model (Section 3.1) by comparing (using the Mean Square Error metric) the information received by the computational processes in the real CPS deployment and the samples predicted by the proposed model.Data were collected for 24 h, and the relative (percentage) Mean Square Error was calculated for all the acquired samples.The experiment was repeated for different values of parameters R P , R T , and R F , which control the complexity of the proposed model.In this experimental phase, these three parameters are considered to have the same value.
The second experiment in this first phase was aimed at analyzing the probability of the proposed reconstruction and protection algorithm to successfully detect the real situation that occurs in the CPS.Some additional ESP32 nodes were deployed to increment the electrical noise in the environment and/or perform Sparse Sensor and Denial of Service attacks.Different situations were generated, with a duration of ten minutes.It was monitored if the proposed algorithm was able to identify them properly.The second experiment had a duration of 24 h too.Results were processed to generate a confusion matrix representing the algorithm's behavior.The experiment was repeated for different values of R un , and R r parameters.During these experiments, both parameters had the same value.
In the second experimental phase, we evaluated the performance and scalability.We measured the computational time needed for the proposed model and the reconstruction and protection algorithm to obtain a final and stable output.The first experiment focused on the mathematical model.The calculation time was analyzed for different values of parameters R P , R T , and R F (all three had the same value) and different quantities of state variables (M).To allow this experiment, the number of sensor nodes in the CPS was increased with each realization.Each configuration was operating for 24 h.The result was obtained as the average of all measurements collected.Finally, the second experiment in this second phase evaluated the computational time required by the reconstruction and protection algorithm.The experiment was repeated for different values of R un , and R r parameters.During these experiments, both parameters had the same value.Different quantities of state variables (M) were also considered.Each configuration was operating for 24 h.The result was obtained as the average of all collected measurements.

Results
To evaluate the behavior of the proposed technology, first, we analyze the precision of our model (Section 3.1), comparing the predicted future CPS states and the actual state finally achieved.Fig. 5 shows the results.As can be seen, the evolution is exponential, as expected from the error in the Taylor series, as the number of terms increases.In general, all configurations show good behavior, although models with only two terms introduce an error of 12% (which may be too high for some applications).The minimum error (2%) may be achieved for models with more than 12 terms.This error is caused by the truncation of the Taylor series, so they can be numerically computed.But, as a counterpart, the resulting finite series does not perfectly represent the original function and we are introducing a numerical error.

Figure 5: Precision of the proposed model
Other limitations in the proposed model (such as the numerical precision of the underlying hardware platform) are also affecting, so, only by increasing the number of terms in Taylor's series cannot reduce the global error as much as desired.But for a very large catalog of applications, an error of 2% is acceptable and can be tolerated.Even, for those scenarios where computationally lightweight solutions are preferred, schemes with four or five terms generate an error of around 6%, which is a standard error for mass non-critical applications.In common applications, errors below 10% can be handled.From these results, we can conclude that the proposed model represents with good precision the physical processes in CPS.
Similarly, we need to analyze the capacity of the proposed reconstruction algorithm to successfully detect the real situation that is happening in the CPS.Fig. 6 shows the results of this experiment.For all possible situations, three regions are identified.First, for low values of R un and R r parameters (below 20), transitory effects are dominant, indicators do not capture properly the CPS behavior, and sensitivity (rate of situation correctly classified) decreases.Random natural variations may be relevant when very short periods are analyzed.To focus on global tendencies, larger collections of data samples are needed.In that way, later, in the central region, R un and R r parameters present good enough values (between 20 and 60), and the proposed algorithm works with a very satisfactory behavior (sensitivity is between 91% and 98%).In this region, short-time transitory effects are not dominant because larger time series are employed in the reconstruction mechanism, and global tendencies are easily detected.But when the values for R un and R r parameters increase beyond a certain limit (60 samples in this case), real fluctuations effects and high-frequency perturbations are ignored when they are aggregated in a large operation.Then, sensitivity decreases.In this region, even natural fluctuations and changes are not significant compared to long-term tendencies.Relevant changes are ignored, because we are integrating too many samples in the same series, and calculation algorithms do not have enough sensitivity.

Figure 6: Precision of the proposed reconstruction and protection
In conclusion, R un and R r parameters must be balanced: small values cause instabilities, while too big values generate a loss of sensitivity.Values between 20 and 60 are the most appropriate region, as shown above (Fig. 5).
On the other hand, the proposed algorithm does not show the same sensitivity when detecting the different situations in a CPS.In general, situations whose probability is calculated using functions with a higher growth rate (such as the exponential) are more sensitive to the quality of indicators representing the CPS (first return maps, STFT, etc.) and then more sensible to the value of R un and R r parameters.This is because small changes in the exponents may generate big variations in the final function.Values must be selected very carefully and according to previous observations.For example, the probability for the SSA-attacked situation is only supported by an exponential function, so it is the one with the most relevant fluctuations.The noisy situation, which includes a linear term, is much flatter.That means SSA-attacked situations are much more difficult to detect, and probably a heuristic calibration process is required in real deployments and scenarios.
Anyway, the sensitivity of the proposed algorithm (in balanced values of R un and R r parameters) is very satisfactory.Noisy situations are detected on 98% of the occasions, while unperturbed and SSA-attacked situations are correctly identified on 92% of the times.The DoS-attacked situation is the one with the worst behavior, but its sensitivity is just slightly lower: 91%.With these results, we can conclude that the proposed algorithm can reconstruct and protect Cyber-Physical Systems against attacks and perturbations.
To go deeper into the analysis of these data, we present the complete confusion matrix (Table 2) for the configuration R un = R r = 40.As can be seen, most errors when identifying the situation in the CPS are false detections of the noisy situation.Probably, that is caused by the linear term in its probability function, which does not reduce its value as much as the exponential function.If this sensitivity needs to be improved, that linear term should be enriched with new indicators and functions.
It is also important to evaluate the performance and scalability of the proposed solution to identify its limitations.Fig. 7 shows the computational time required for the proposed model to operate.As can be seen, the evolution of the computational delay is linear.This is because our model consists of additions and multiplications, without loops or recursive problems.Besides, each new state variable is independent of the others, so the increase is linear.This facilitates the employment of this protection and reconstruction solution in future dense and pervasive scenarios, where up to ten million sensors per square kilometer could be deployed.
Moreover, the delay is always in the range of milliseconds.Additions and multiplications are performed very efficiently on modern computers, and they require a short time to perform millions of operations.In this case, even for a CPS that includes 100 devices (i.e., 400 state variables) and very complex models (with almost 20 terms in the Taylor series), the computational time required to operate the model is below 100 milliseconds (70 milliseconds, to be precise).Considering the most usual Cyber-Physical Systems capture information from the environment every few seconds, this delay is satisfactory.
Finally, the same scalability and performance analysis must be applied to the reconstruction and protection algorithm.Fig. 8 shows the results.Here, again, the evolution is almost linear, because all the proposed computational procedures do not require any loop or recursive processing.In this case, for the largest deployment (one hundred devices) and a typical value for R un and R r parameters, the delay is above one second.This may be slightly above the acceptable maximum for certain critical realtime applications.For smaller deployments and the same configuration, the delay is below one second (between 100 and 600 milliseconds).But even delays above one second are acceptable in mass noncritical applications, where data are acquired every few seconds.Due to linear evolution, scalability is guaranteed (even in future scenarios) as consumed resources grow at the same rate as the number of sensor nodes in the physical platform.In conclusion, considering the limitations that may arise in critical real-time applications, the performance of the proposed reconstruction and protection mechanism is satisfactory.

Conclusions and Future Works
This paper presents a new stochastic model to represent the behavior of Cyber-Physical Systems precisely.This model includes unknown multivariate discrete and continuous-time functions and different multiplicative noises to represent the evolution of physical processes and random effects in the physical and computational worlds.As a novelty, in this model, engineered processes such as the digitalization stage are represented too.Additionally, and contrary to the commonly employed deterministic attackers, in this new model attackers are described through a stochastic process.Standard error sources are estimated through different indicators and non-linear techniques (such as the Fourier transform, first-return maps, or the probability density function).Finally, the reconstruction mechanism consists of a weighted stochastic model combining all error sources.The actual reconstructed value is generated as the output from a decision algorithm.
Experimental results show that the precision of the proposed model is above 90%, with a residual error between 6% and 2% for the most common configurations.Additionally, the sensitivity of the proposed reconstruction and protection algorithm is up to 92%.Considering all this, the proposed solution is a valid security scheme for CPS.Future works will analyze new indicators and probability functions to improve the sensitivity, especially in noisy situation.In addition, the solution will be deployed in real industrial scenarios with legacy systems, to study the impact of second-order effects such as reduced connectivity or human accidents and manipulations.
as thermal noise or environmental radiation) are mixed with real physical variables x s i (t) in the transformation function T m .On the other hand, multiplicative electrical noises ξ

r
(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[•]

Figure 2 :
Figure 2: Protection and reconstruction mechanism

Figure 4 :
Figure 4: Reconstruction and protection algorithm

Figure 8 :
Figure 8: Computational delay and scalability (reconstruction and protection algorithm)

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

Table 1 :
State-of-the-art

Table 2 :
Confusion matrix for the configuration R un = R r= 40