OCT 2017

1

Characterization of Model-Based Detectors for CPS
Sensor Faults/Attacks

arXiv:1710.04407v1 [] 12 Oct 2017

Carlos Murguia and Justin Ruths

Abstract—A vector-valued model-based cumulative sum
(CUSUM) procedure is proposed for identifying faulty/falsified
sensor measurements. First, given the system dynamics, we derive
tools for tuning the CUSUM procedure in the fault/attack free
case to fulfill a desired detection performance (in terms of false
alarm rate). We use the widely-used chi-squared fault/attack detection procedure as a benchmark to compare the performance of
the CUSUM. In particular, we characterize the state degradation
that a class of attacks can induce to the system while enforcing
that the detectors (CUSUM and chi-squared) do not raise alarms.
In doing so, we find the upper bound of state degradation that is
possible by an undetected attacker. We quantify the advantage of
using a dynamic detector (CUSUM), which leverages the history
of the state, over a static detector (chi-squared) which uses a
single measurement at a time. Simulations of a chemical reactor
with heat exchanger are presented to illustrate the performance
of our tools.
Index Terms—Cyber Physical Systems, Model-based
fault/attack detection, Security, CUSUM, Chi-squared.

I. I NTRODUCTION
During the past half-century, scientific and technological
advances have greatly improved the performance of control systems. From heating/cooling devices in our homes,
to cruise-control in our cars, to robotics in manufacturing
centers. However, these new technologies have also led to
vulnerabilities of some our most critical infrastructures–e.g.,
power, water, transportation. Advances in communication and
computing power have given rise to adversaries with enhanced
and adaptive capabilities. Depending on attacker resources and
system defenses, attackers may deteriorate the functionality of
systems even while remaining undetected. Therefore, designing efficient fault/attack detection schemes and attack-robust
control systems is of key importance for guaranteeing the
safety and proper operation of critical systems. Tools from
sequential analysis and fault detection have to be adapted to
deal with the systematic, strategic, and persistent nature of
attacks. These new challenges have attracted the attention of
many researchers in the control and computer science communities [1]-[10]. Lately, there has been increasing interest in
studying systems performance degradation induced by attacks
Carlos Murguia is with the Center for Research in Cyber Security (iTrust),
Engineering Systems and Design (ESD) Pillar, Singapore University of
Technology and Design, Singapore. E-mail: murguia_rendon@sutd.edu.sg.
Justin Ruths is with the Departments of Mechanical and Systems Engineering, University of Texas at Dallas, USA, e-mail: jruths@utdallas.edu.
This work was supported by the National Research Foundation (NRF),
Prime Minister’s Office, Singapore, under its National Cybersecurity R&D
Programme (Award No. NRF2014NCR-NCR001-40) and administered by the
National Cybersecurity R&D Directorate.
Manuscript received May 5, 2017; revised May 5, 2017.

that remain hidden or undetected by detection procedures [1][2],[5]-[6],[10]. Quantifying the system degradation provides a
measure of impact to assess the performance of control structures, estimation schemes, and detection procedures against
this class of intelligent attacks. For instance, in [5]-[6], for
arbitrary detection procedures, the authors quantify how much
the attacker can deviate the estimate of the state from its
attack-free values while remaining stealthy. They characterize
stealthiness of attacked sequences using the Kullback-Leibler
Divergence [11] between the attack-free and the attacked
sequence. In the same spirit, the authors in [1]-[2] study how
attacks propagate through the control structure to degrade the
system dynamics while remaining undetected by the detection mechanism. In particular, the authors in [1] characterize
undetectability (for a class of deterministic LTI systems) as
the ability of attackers to excite only the zero dynamics of
the system [12] (making its effect undetectable from output
measurements). In [2],[10],[13], the authors propose a notion
of stealthiness by attacks that do not change the alarm rate
of the detector by more than a small amount (thus making
it hard for the operator to distinguish between an attack-free
and an attacked system, i.e., these attacks remain hidden from
the detector). As a measure of impact, they characterize the
reachable sets that these hidden attacks can induce to the
system.
Most of the current work on security of control systems
has focused on static detection procedures (either bad-data
or chi-squared detectors), which identify anomalies based on
a single measurement at a time [1]-[4],[10]. There is only
a small amount of literature considering the use of dynamic
change detection procedures such as the Sequential Probability
Ratio Test (SPRT) or the Cumulative Sum (CUSUM) [14],
which employ measurement history, in the context of security of Cyber-Physical Systems (CPS) [7],[15]-[16]. Dynamic
detectors present an appealing alternative to the aforementioned static procedures. Using measurement history provides
extra degrees of freedom for improving the performance of
our fault/attack detection strategies; in particular, against low
amplitude persistent attacks [16].
This paper addresses, for Linear Time-Invariant (LTI) systems subject to sensor/actuator noise, the problem of characterizing CUSUM dynamic and chi-squared static detectors
in terms of false alarm rates and performance degradation
under a class of attacks. Standard Kalman filters are proposed
to estimate the state of the physical process. Both detectors
employ a distance measure that is a quadratic function of
the residual (the error between sensor measurements and the
estimated outputs). In the chi-squared procedure, at each time

OCT 2017

instant if the distance measure is larger than a threshold, an
alarm is raised, indicating a possible compromised sensor.
In the CUSUM procedure, the distance measure values are
accumulated over time, and if this accumulated value is greater
than expected an alarm is triggered. Fundamentally a detector
aims to properly raise an alarm when a fault/attack happens
and not raise an alarm when there is no fault/attack. Deviation
from this ideal performance is captured by false positives
(alarms are raised when there are no faults/attacks), also called
false alarms, and false negatives (a fault/attack happens, but
no alarm is raised). Although minimizing both false positives
and false negatives is best, often they must be traded off based
on which is more tolerable. In this context, detectors with high
sensitivity would have high rates of false positives in favor of
low rates of false negatives (and vice versa).
In order to provide an equitable comparison between detectors (e.g., static versus dynamic), we first require the ability to
tune each type of detector to a similar level of sensitivity. Todate, however, there is not a complete characterization of how
features of the system (e.g., system matrices, control/estimator
gains, noise, sampling) affect the selection of the CUSUM
parameters to achieve a desired sensitivity, quantified by the
rate of false alarms. Our first contribution in this paper is
to provide systematic tools to tune the CUSUM detector,
and for completeness, also the chi-squared (static) detector,
in the fault/attack free case based on the system dynamics,
the Kalman filter, the stochastic properties of the distance
measure, and a desired false alarm rate. In particular, sufficient
conditions for mean square boundedness of the CUSUM
sequence are derived when it is driven by a quadratic form
of the residual. Then, using a Markov chain approximation
of the CUSUM sequence, we give a procedure for selecting
the decision threshold such that a desired false alarm rate is
satisfied.
Second, for a class of zero-alarm attacks (attacks that
prevent the detector from raising alarms), we characterize
the impact of the attack sequence on the system dynamics
when the vector-valued CUSUM and chi-squared detectors
are deployed for attack detection. From an empirical point of
view, zero-alarm attacks have been actively used to assess the
resilience of dynamical systems against attacks [7],[16],[17].
zero-alarm attacks provide a simple, deterministic, yet representative class of attacks that can be easily scaled to systems
with different dynamics and properties; thus making them a
good choice for assessing the performance of attack detectors
in terms of state degradation.
In our preliminary work [8], we have started analyzing these
ideas. The contributions of this manuscript with respect to [8]
are the following: A comprehensive and complete exposition
of all the results and methodologies; the main results have been
revised and improved and the corresponding proofs (which
are not given in our preliminary work) are included in this
paper; we formulate an all-new measure for attack degradation
centered around the concept of input-to-state stability; and a
benchmark simulation experiment used in the fault-detection
literature [18],[19] (a chemical reactor with heat exchanger)
is presented to illustrate the performance of our tools.

2

A. Notation
Throughout this paper, the following notation is used: the
symbol R stands for the real numbers, R>0 (R≥0 ) denotes
the set of positive (non-negative) real numbers. The symbol
N stands for the set of natural numbers. The Euclidian norm
in Rn is denoted by kxk, kxk2 = xT x, where T denotes
transposition. The induced norm of a matrix A ∈ Rn×n ,
denoted by kAk, is defined as kAk = maxx∈Rn ,kxk=1 kAxk.
The n × n identity matrix is denoted by In or simply I if no
confusion can arise. Similarly, n × m matrices composed of
only ones and only zeros are denoted by 1n×m and 0n×m ,
respectively, or simply 1 and 0 when their dimensions are
clear. If a quadratic form xT P x with a symmetric matrix
P = P T is positive definite (semidefinite), then P is called
positive definite (semidefinite). For positive definite (semidefinite) matrices, we use the notation P > 0 (P ≥ 0); moreover,
P > Q (P ≥ Q) means that the matrix P − Q is positive
definite (semidefinite). The spectrum of a matrix A is denoted
by spec[A], tr[A] denotes its trace, and ρ[A] is its spectral
radius. The notation λmin [A] (λmax [A]) stands for the smallest
(largest) eigenvalue of the square matrix A. The notation
E[x] stands for the expected value of x and Ey [x] denotes
the expected value of x conditional to y. The variance of a
random variable x is denoted by var[x]. The notation pr[·]
denotes probability and x ∼ N (µ, Σ) means that x ∈ Rn is a
vector-valued normally distributed random variable with mean
µ ∈ Rn and covariance matrix Σ ∈ Rn×n . For simplicity of
notation, we often suppress the explicit dependence of time t.
II. S YSTEM D ESCRIPTION & ATTACK D ETECTION
We study
 LTI stochastic systems of the form:
x(tk+1 ) = F x(tk ) + Gu(tk ) + v(tk ),
(1)
y(tk ) = Cx(tk ) + η(tk ),
with sampling time-instants tk , k ∈ N, state x ∈ Rn , measured
output y ∈ Rm , control input u ∈ Rl , matrices F , G, and C
of appropriate dimensions, and i.i.d. multivariate zero-mean
Gaussian noises v ∈ Rn and η ∈ Rm with covariance matrices
R1 ∈ Rn×n , R1 ≥ 0 and R2 ∈ Rm×m , R2 ≥ 0, respectively.
The initial state x(t1 ) is assumed to be a Gaussian random
vector with covariance matrix R0 ∈ Rn×n , R0 ≥ 0. The
processes v(tk ), k ∈ N and η(tk ), k ∈ N and the initial
condition x(t1 ) are mutually independent. It is assumed that
(F, G) is stabilizable and (F, C) is detectable. At the timeinstants tk , k ∈ N, the output of the process y(tk ) is sampled
and transmitted over a communication network. The received
output ȳ(tk ) is used to compute control actions u(tk ) which
are sent back to the process, see Fig. 1. The complete controlloop is assumed to be performed instantaneously, i.e., the
sampling, transmission, and arrival time-instants are equal. In
this paper, we focus on attacks on sensor measurements. That
is, in between transmission and reception of sensor data, an
attacker may replace the signals coming from the sensors to the
controller1, see Fig. 1. After each transmission and reception,
1 Such an attack can also be accomplished by installing malware on the
controller equipment, in which case true measurements reach the controller,
but are manipulated before they are used.

OCT 2017

3

assumption that guarantees that the Kalman filter converges),
the estimator gain:

−1
Lk = L := F P C T R2 + CP C T
,
(7)
leads to the minimal steady state covariance matrix P , with
P given by the solution of the algebraic Riccati equation:
F P F T − P + R1 = F P C T (R2 + CP C T )−1 CP F T . (8)
The reconstruction method given by (4)-(8) is referred to as
the steady state Kalman filter, cf. [22].

Fig. 1: Cyber-physical system under sensor attacks.

the attacked output ȳ takes the form:
ȳ(tk ) := y(tk ) + δ(tk ) = Cx(tk ) + η(tk ) + δ(tk ),
(2)
m
where δ(tk ) ∈ R denotes additive sensor attacks/faults.
Denote xk := x(tk ), uk := u(tk ), vk := v(tk ), ȳk := ȳ(tk ),
ηk := η(tk ), and δk := δ(tk ). Using this new notation, the
attacked system is written in the following compact form:


xk+1 = F xk + Guk + vk ,
ȳk = Cxk + ηk + δk .

(3)

Remark 1 If the stochastic processes vc (tk ) and η(tk ) are
non-Gaussian, using spectral factorization [20], [21], we could
rewrite them as output signals coming from linear filters,
say G1 (q) and G2 (q), with Gaussian stochastic processes as
inputs, say w1 (tk ) and w2 (tk ); that is, vc (tk ) = G1 (q)w1 (tk )
and η(tk ) = G2 (q)w2 (tk ), where q denotes the forwardshift operator. Then, by extending the system dynamics with
the filters and considering the non-Gaussian noises vc (t) and
η(tk ) as new states, the extended system is written as a LTI
system perturbed by Gaussian noise, see, for instance, [20],
[21] for details.

A. Steady state Kalman filter (attack/fault free case)
To estimate the state of the process, a one step-ahead
estimator with the following structure is proposed:

(4)
x̂k+1 = F x̂k + Guk + Lk ȳk − C x̂k ,
with estimated state x̂k ∈ Rn , x̂1 = E[x(t1 )], and gain
matrix Lk ∈ Rn×m . Define the estimation error ek := xk − x̂k .
The matrix Lk is designed to minimize the covariance matrix
Pk := E[ek eTk ] in the absence of attacks. Given the discretetime dynamics (3) and the estimator (4), the estimation error
is governed by the difference equation:

ek+1 = F − Lk C ek − Lk ηk − Lk δk + vk .
(5)
If the pair (F, C) is detectable, the covariance matrix converges to steady state in the sense that, in the attack-free case,
limk→∞ Pk = P exists [22]. Let δk = 0; then, from (5), the
mean value of ek is given by

E[ek+1 ] = F − Lk C E[ek ].
(6)
Because x̂1 = E[x(t1 )], the mean value of the estimation error
equals 0n×1 independent of Lk . We assume that the system
has reached steady state before an attack occurs. Then, the
estimation of the random sequence xk , k ∈ N can be obtained
by the estimator (4) with Pk and Lk in steady state. It can be
verified that, if CP C T + R2 is positive definite (a standard

Remark 2 It is well known that, if the noise sequences vk
and ηk are Gaussian, the Kalman filter (4)-(8) gives the best
estimate x̂k+1 of the state xk+1 (in terms of minimum-meansquare estimation error) from noisy measurements. Moreover,
if the noise is not Gaussian, the Kalman filter is the best linear
estimator; although there may exist nonlinear estimators with
better performance, cf. [23].
B. Residuals and hypothesis testing
Attacks can be regarded as induced faults in the system.
Then, it is reasonable to use existing fault detection techniques
to identify sensor attacks. The main idea behind fault detection
theory is the use of an estimator to forecast the evolution
of the system in the absence of faults. This prediction is
compared with the actual measurements from the sensors. If
the difference between what it is measured and the estimation
(often referred to as residual) is larger than expected, there
might be a fault in the system. Although the notion of residuals
and model-based detectors is now routine in the fault detection literature, the primary focus has been on detecting and
isolating faults with specific structures (e.g., constant biases in
sensor measurements or random faults in sensors and actuators
following specific distributions). Now, in the context of an
intelligent adversarial attacker, new challenges arise to understand the effect that an intruder can have on the system given
the dynamics, the estimator, and the detector structure. In this
work, we use the steady state Kalman filter introduced in the
previous section as our estimator.
Consider the discrete-time process dynamics (3), the steady
state Kalman filter (4)-(8), and the corresponding error difference equation (5). Define the residual sequence rk , k ∈ N
as
rk := ȳk − C x̂k = Cek + ηk + δk .
(9)
Then, 
rk evolves according to the difference equation:
ek+1 = F − LC ek − Lηk + vk − Lδk ,
(10)
rk = Cek + ηk + δk .
If there are no faults/attacks, the mean of the residual is
E[rk+1 ] = CE[ek+1 ] + E[ηk+1 ] = 0m×1 ,
(11)
and the covariance matrix is given by
T
E[rk+1 rk+1
] = CP C T + R2 =: Σ ∈ Rm×m .
(12)
For this residual, we identify two hypothesis to be tested: H0
the normal mode (no faults/attacks) and H1 the faulty mode
(with faults/attacks). Under the normal mode, the statistics of
the residual are:

E[rk ] = 0m×1 ,
H0 :
(13)
E[rk rkT ] = Σ.

OCT 2017

Therefore, when an fault/attack occurs in the system (H1 ), we
expect that the statistics of the residual are different from the
normal mode, i.e., 
E[rk ] 6= 0m×1 , or
H1 :
(14)
E[rk rkT ] 6= Σ.
There exist many well-known hypothesis testing techniques
which may be used to examine the residual and subsequently
detect faults/attacks. For instance, Sequential Probability Ratio
Testing (SPRT) [24], [25], Cumulative Sum (CUSUM) [14],
[26], Generalized Likelihood Ratio (GLR) testing [27], Compound Scalar Testing (CST) [28], etc. Each of these techniques
has its own advantages and disadvantages depending on the
scenario. The most utilized and powerful one is, arguably
the SPRT, which minimizes the time to reach a decision for
given probabilities of false detection (i.e., declaring H1 when
it is actually H0 ). In this manuscript, we mainly focus on the
CUSUM procedure which is a version of SPRT that permits
repeated detection [26]. However, for comparison, we also
present results about a particular case of CST, namely the
so-called chi-squared change detection procedure.
C. Distance measures and CUSUM procedure
Change detection theory was founded by Wald in 1947 when
his book "Sequential Analysis" was published and the SPRT
was first introduced. Subsequently, the CUSUM procedure
[26] was proposed by Page to detect changes in the mean
of random variables by testing a weighted sum of the last few
observations, i.e., a moving average. As Page pointed out, the
CUSUM is equivalent to a repeated SPRT in which the test
is restarted once a change has been detected. The input to the
CUSUM procedure is a distance measure zk ∈ R, k ∈ N, i.e.,
a measure of how deviated the estimator is from the sensor
measurements. We propose the quadratic distance measure
zk := rkT Σ−1 rk ,
(15)
where rk and Σ are the residual sequence and its covariance
matrix defined in (9) and (12), respectively. If there are no
attacks, E[rk ] = 0 and E[rk rkT ] = Σ; it follows that

E[zk ] = tr[Σ−1 Σ] + E[rk ]T Σ−1 E[rk ]



= m,
var[z
]
=
2tr[Σ−1 ΣΣ−1 Σ] + 4E[rk ]T Σ−1 ΣΣ−1 E[rk ]

k


= 2m,
(16)
see, e.g., [11] for details. Moreover, since rk ∼ N (0, Σ),
then zk = rkT Σ−1 rk follows a chi-squared distribution with
m degrees of freedom, cf. [11]. Other options are based on
likelihood ratios. In this case, instead of directly using the
sequence zk to drive the CUSUM procedure, the log-likelihood
ratio Λk (zk ) between the two hypotheses is employed:
f 1 (z|H1 )
Λk (zk ) := log z0k
,
(17)
fzk (z|H0 )
where fzjk (z|Hj ) denotes the Probability Density Function
(PDF) of the distance measure zk , k ∈ N under Hj ,
j = {0, 1}. A problem to address when using log-likelihood
ratios for detecting attacks or unstructured faults is the fact
that the PDF of the faulty sequence f1 (zk |H1 ) is unknown.

4

Actually, in the case of attacks, the adversary may induce any
arbitrary (and possibly) non-stationary sequence zk . Assuming
the statistical properties of the attack sequences may limit our
ability to detect a wide range of attacks [7].
The CUSUM procedure of Page driven by the distance
measure zk is defined as follows.
CUSUM:

 S1 = 0,
Sk = max(0, Sk−1 + zk − b), if Sk−1 ≤ τ,
(18)

Sk = 0 and k̃ = k − 1,
if Sk−1 > τ.
Design parameters: bias b ∈ R>0 and threshold τ ∈ R>0 .
Output: alarm time(s) k̃.
The idea is that the test sequence Sk accumulates the
distance measure zk and alarms are triggered when Sk exceeds
the threshold τ . The test is reset to zero each time Sk becomes
negative or larger than τ . If zk is an independent non-negative
sequence (which is our case) and b is not sufficiently large, the
CUSUM sequence Sk grows unbounded until the threshold τ
is reached, no matter how large τ is set. In order to prevent
these drifts, inevitably leading to false alarms, the bias b must
be selected properly based on the statistical properties of the
distance measure. Once the the bias is chosen, the threshold τ
must be selected to fulfill a required false alarm rate A∗ (see
Section III-B).
III. CUSUM-T UNING
To enhance the performance of the CUSUM procedure, the
bias b and the threshold τ must be selected appropriately.
We have already mentioned that too small a bias can lead to
inevitable growth of the CUSUM test sequence. At the same
time, too large a bias may hide the effect of faults/attacks. In
what follows, we provide tools for selecting these parameters
given the statistical properties of the distance measure zk introduced in (16). In particular, we provide sufficient conditions
on the bias b such that, in the absence of faults/attacks, the
sequence Sk of the CUSUM remains bounded (independent of
the reset due to τ ) in mean-squared sense. This is important
to avoid false alarms due to the inherent divergence of Sk .
Subsequently, we characterize the false alarm rate of the
CUSUM in terms of b and τ given a desired false alarm rate.
A. Boundedness
First, we introduce the following concept of boundedness
of stochastic processes, cf. [29],[30], followed by sufficient
conditions for boundedness of the CUSUM sequence.
Definition 1 The sequence Sk , k ∈ N is said to be bounded
in mean square, if
 
sup ES1 Sk2 < ∞,
k∈N

is satisfied, i.e., the second moment of Sk is finite.

Theorem 1 Consider the discrete-time process (3) and the
steady state Kalman filter (4)-(8). Assume that there are no
attacks to the system, i.e., δk = 0. Let the CUSUM (18)

OCT 2017

with bias b ∈ R>0 and threshold τ ∈ R>0 be driven by
the distance measure zk = rkT Σ−1 rk , k ∈ N with residual
sequence rk ∼ N (0, Σ), k ∈ N. Then, if the bias is set
larger than the number of measurements, b > b̄ := m, the
CUSUM sequence Sk , k ∈ N is bounded in mean square
sense independent of the threshold τ .
The proof of Theorem 1 is presented in the appendix.
Remark 3 Notice that boundedness of the first moment
follows from boundedness of the second moment and Jensen’s
inequality [11]. Then, b > b̄ = m implies that the expected
value ES1 [Sk ], k ∈ N is finite.
The result stated in Theorem 1 implies that for b > b̄, the
second moment (and hence the first) of the sequence Sk ,
k ∈ N does not diverge. Consequently, we avoid false alarms
due to intrinsic growth of the CUSUM sequence. Note that if
the bias b is selected greater than but close to b̄, small changes
in the distance measure zk would lead to divergence of Sk .
Therefore, the smaller the bias, the higher the sensitivity
against changes in (or uncertain characterization of) the
residual signals.

B. False Alarms
Once the bias is selected such that boundedness of the second moment E[Sk2 ] is guaranteed, the next step is to select the
threshold τ to fulfill a desired false alarm rate. The occurrence
of an alarm in the CUSUM when there are no faults/attacks
to the CPS is referred to as a false alarm. Operators need to
tune this false alarm rate depending on the application. To do
this, the threshold τ must be selected to fulfill a desired false
alarm rate A∗ . Let A ∈ [0, 1] denote the false alarm rate of the
procedure defined as the expected proportion of observations
which are false alarms, i.e., for the CUSUM procedure, A :=
pr[Sk ≥ τ ], see [31] and [32]. Define the run length K of
the CUSUM (18) as the number of iterations needed such that
SK > τ (without attacks):
K := min{k ≥ 1 : Sk > τ }.
(19)
The expected value E[K] of K is known in the literature as the
Average Run Length (ARL). The ARL is inversely proportional
to the false alarm rate A [32], [31], i.e.,
A = 1/ARL.
(20)
Then, for a given b > b̄, the problem of selecting τ to satisfy a
desired false alarm rate A∗ can be reformulated as the problem
of selecting τ such that
ARL = 1/A∗ .
(21)
To determine a pair (b, τ ) satisfying (21), an expression for the
ARL = E[K] is required but, in general, its exact evaluation
is analytically intractable [33]. The problem of approximating
the ARL for CUSUM procedures has been addressed by many
authors during the last decades. For instance, the authors
in [33]-[35] propose Wiener process approximations of the
ARL using analogies between the CUSUM and the SPRT

5

for normally distributed distance measures. Although these
techniques lead to explicit formulas for evaluating the ARL,
the obtained approximations are often too conservative, see
[34]-[35]. Accurate numerical methods have been proposed
by, for instance, [36]-[39]. These methods rely on two main
techniques, namely Markov chain and integral equation approaches. Both methods give accurate predictions of the ARL
(see [36] for a comparison); however, we find the Markov
chain approach more constructive and easier to implement. In
this work, we use the result of Evans and Brook [37]. With
this result, we outline a procedure for selecting the threshold
τ given the bias b and a required false alarm rate A∗ .

For given b > b̄ and some τ ∈ R>0 , consider the sequence
Sk generated by the CUSUM procedure (18) driven by the
distance measure zk = rkT Σ−1 rk , k ∈ N. Given the recursive
nature of the CUSUM procedure and independence of vk and
ηk , k ∈ N, the sequence Sk forms a Markov chain taking
values on the non-negative real line [40]. By discretizing
the probability distribution of the distance measure, it is
possible to subdivide the CUSUM sequence Sk into a finite
set of partitions. The idea is to approximate the continuous
scheme by a Markov chain having N + 1 states labeled
as {E0 , E1 , . . . , EN }, where EN is absorbing. Then, the
probability that the chain remains in the same state at the
next step should correspond to the case when Sk does not
change in value by more than a small amount, say 12 ∆S ,
i.e., the next distance measure zk does not differ from the
bias b by more than 21 ∆S . The constant ∆S determines the
width of the grouping interval involved in the discretization
of the probability distribution of zk . The interval width 12 ∆S
must be selected such that the probability of jumping from Ej ,
j ∈ {0, . . . , N −1} to the absorbing state EN is approximately
equal to the probability that the CUSUM sequence Sk jumps
beyond the threshold τ from a position Sk−1 ∈ (0, τ ) which
corresponds approximately to the state Ej . This requirement
is satisfied by taking
2τ
,
(22)
∆S :=
2N − 1
see [37] for details. Then, the transition probabilities from a
starting state Ej , j = 0, . . . , N − 1, can be determined from
the probability distribution of zk − b = rkT Σ−1 rk − b, as:
pr(Ej → E0 ) = pr(zk − b ≤ −j∆S + 12 ∆S ),
pr(Ej → EN ) = pr((N − j)∆S − 21 ∆S < zk − b),

pr(Ej → Eν ) = pr(zk − b ≤ (ν − j)∆S + 12 ∆S )

− pr(zk − b < (ν − j)∆S − 12 ∆S ).
Note that pr(E0 → EN ) = pr(zk − b > τ ). For given
b and τ , the states {E0 , . . . , EN } and the above transition
probabilities forms a Markov chain whose transition matrix
can be constructed from the probability distribution of zk − b.
Denote Tχ := pr(zk − b ≤ χ∆S + 21 ∆S ) and pχ :=
pr(χ∆S − 12 ∆S < zk − b ≤ χ∆S + 21 ∆S ). Then, the Markov

OCT 2017

6

transition matrix P ∈ R(N +1)×(N +1) is given by:


T0
p1
p2 . . . pN −1 1 − TN −1
 T−1 p0
p1 . . . pN −2 1 − TN −2 



 ..
..
..
..
..

 .
.
.
.
.



T
p
p
.
.
.
p
1
−
T
P := 
N −1−j
N −1−j  . (23)
 −j 1−j 2−j

 ..
..
..
..
..

 .
.
.
.
.


T1−N p2−N p3−N . . . p0
1 − T0 
0
0
0 ...
0
1
Since the state EN is absorbing, the last row consists of zeros
except for the last entry. To compute the transition probabilities
Tχ and pχ of P, we need the Cumulative Distribution Function
(CDF) of the shifted distance measure zk − b = rkT Σ−1 rk − b.
If there are no attacks, rk ∼ N (0, Σ); therefore, zk −b follows
a shifted chi-squared
 distribution
with CDF:
x+b
P m
,
, for x ≥ −b,
2
2
Fzk −b (x) :=
(24)
0, for x < −b,
where P(·, ·) denotes the regularized lower incomplete gamma
function [11]. Then, the entries of the transition matrix are
given by


1

 pχ = Fzk −b χ∆S + 2 ∆S
−Fzk −b χ∆S − 21 ∆S ,
(25)



1
Tχ = Fzk −b χ∆S + 2 ∆S .

Define the transformation T := (IN 0N ×1 ) ∈ R
the matrix:
R := T PT T ∈ RN ×N .

N ×(N +1)

and

(26)

The matrix R is known as the fundamental matrix associated
with the Markov transitionmatrix P.Note that
R
∗
P=
.
01×N 1
Then, all entries of R are non-negative and its row sums
are less than one. Therefore, by Gershgorin circle theorem,
the eigenvalues of R satisfy: 1 > |λN | ≥ . . . ≥ |λ1 |. It
follows that ρ[R] < 1, where ρ[·] denotes spectral radius;
therefore, the matrix (IN − R) is invertible [41]. Next, having
introduced the transition matrix P of the approximated
Markov chain and the fundamental matrix R, we can compute
an approximation Ã of the false alarm rate A based on the
result in [37], equation (20), and (22)-(26).
Theorem 2 Assume that there are no attacks on the system
and let the CUSUM (18) with bias b > b̄ = m and
threshold τ ∈ R>0 be driven by the distance measure
zk = rkT Σ−1 rk with residual sequence rk ∼ N (0, Σ), k ∈ N.
For a finite number of partitions N ∈ N, consider the
fundamental matrix R, defined in (26), obtained from the
transition matrix P (22)-(25), and define
µ := (IN − R)−1 1N ×1 = [µ1 , . . . , µN ]T .
(27)
Then, the false alarm rate A = 1/ARL, is approximately
given by Ã := µ−1
1 . Moreover, as N → ∞, Ã → A, i.e.,
limN →∞ Ã = A.

Proof : Consider the Markov chain M given by the states
{E0 , E1 , . . . , EN } and the transition matrix P (22)-(25).
Let K̃ ∈ N denote the number of iterations needed to
reach the absorbing state EN from E0 . The random variable
K̃ follows a discrete phase-type distribution with E[K̃] =

µ1 and µ1 as defined in (27), see [42]. By construction,
M is a finite state approximation of the continuous Markov
chain formed by the CUSUM sequence Sk ∈ R≥0 , k ∈ N
driven by zk = rkT Σ−1 rk , k ∈ N. It follows that E[K̃] =
µ1 ≈ E[K] = ARL where K denotes the run length of
the CUSUM defined in (19). Then, from (20), we have that
A = ARL−1 ≈ µ−1
1 . Next, increasing the number of partitions
N would reduce the width of the grouping interval ∆S (22),
such that, as N → ∞, the Markov chain M retrieves the
continuous scheme given by the CUSUM sequence Sk ∈ R≥0 ,
k ∈ N, (18); therefore, A = 1/ limN →∞ E[K̃].

Remark 4 Theorem 2 provides a tool for approximating the
false alarm rate A of the CUSUM procedure for given bias
b and threshold τ . In particular, for a given b > b̄, it
provides a map S : R>0 → (0, 1) from the threshold τ
to the approximated false alarm rate Ã, i.e., τ 7→ S(τ ),
Ã = S(τ ). Given that Fzk −b (z) is a continuous function
for all z ∈ R, it can be proved that S(τ ) is continuous for
all τ ∈ R>0 ; then, simple bisection methods can be used to
determine the threshold τ = τ ∗ ∈ R>0 required to satisfy
Ã = S(τ ∗ ) = A∗ for given b > b̄.
IV. C HI - SQUARED T UNING
The CUSUM approach to fault/attack detection offers an
compelling alternative to the more popular chi-squared detector. Here, we use the chi-squared approach as a benchmark
to compare the performance of the CUSUM. Consider again
the residual sequence rk , (10), and its covariance matrix
Σ, (12). The chi-squared procedure is defined as follows:
Chi-squared procedure:
If zk = rkT Σ−1 rk > α, k̃ = k.
Design parameter: threshold α ∈ R>0 .
Output: alarm time(s) k̃.

(28)

The idea is that alarms are triggered if zk exceeds the
threshold α. Similar to the CUSUM procedure, the parameter
α is selected to satisfy a required false alarm rate A∗ .
Theorem 3 Assume that there are no attacks on the system
and consider the chi-squared procedure (28) with threshold
∗
α ∈ R>0 , rk ∼ N (0, Σ). Let α = α∗ := 2P−1 ( m
2 , 1 − A ),
−1
where P (·, ·) denotes the inverse regularized lower incomplete gamma function, then A = A∗ .
Proof : Let K̃ denote the run length of the chi-squared procedure (28) defined as the number of iterations needed such
that k = K̃ implies zk > α when there are no attacks.
As with the CUSUM procedure, the Average Run Length
is given by ARL = E[K̃] and the false alarm rate satisfies
A = 1/ARL. The random variable K̃ follows a geometric
distribution [11]; therefore, ARL = E[K̃] = 1/pr(zk > α)
and A = pr(zk > α). Each element of the sequence zk ,
k ∈ N is an i.i.d. random variable with CDF given by
x
m x
F̃zk (x) = P( m
2 , 2 ). Then, A = pr(zk > α) = 1 − P( 2 , 2 )
and the result follows.


OCT 2017

V. D ETECTOR P ERFORMANCE UNDER Z ERO -A LARM
ATTACKS
In this section, we assess the performance of the CUSUM
procedure by quantifying the effect of the attack sequence δk
on the estimation error when the CUSUM procedure is used
to identify anomalies. To maintain an equitable comparison
between detectors in this section, some assumption must be
made about their false positive rate (false alarm rate) and
false negative rate (the rate at which true attacks are not
detected). Using the tools introduced in prior sections, we can
calibrate the CUSUM and chi-squared detectors to have the
same false alarm rate. Here, we consider a class of zero-alarm
attacks, i.e., attack sequences that keep the detector from
raising alarms. This implies that the entire attacked distribution
(of zk or Sk ) is at or below the decision threshold, effectively
maximizing the false negative rate (since the true positive rate
is zero, i.e., the true attack is never detected). Zero-alarm
attacks provide a concise quantification of attacker impact
on the system performance. In particular, we characterize
the estimation error deviation due to zero-alarm attacks. This
serves as a useful proxy for the capabilities of the attacker
due to the detection mechanism. To do this, using the notion
of input to state stability [43]-[44], we derive upper bounds
on the trajectories of the estimation error given the system
dynamics, the attack sequence, and the CUSUM parameters.
Furthermore, we compare the performance of the CUSUM
against the chi-squared detector. In this paper, we have now
characterized the rate of false alarms based on the type of
detector and the system and detector parameters. We now
close the loop on this analysis by identifying the impact that
attackers can have while taking advantage of the detector
structure.
A. Zero-alarm Attacks
Here, we quantify the damage that attacks may induce to
the estimation error dynamics while enforcing that alarms are
not raised by the detector. We assume that the attacker has
perfect knowledge of the system dynamics, the Kalman filter,
control inputs, measurements, and detection procedure (either
CUSUM or chi-squared). It is further assumed that all the
sensors can be compromised by the attacker at each time step
(a worst-case scenario).
First, consider the chi-squared procedure (28) and write zk
in terms of the estimation error ek :
zk = (Cek + ηk + δk )T Σ−1 (Cek + ηk + δk ).
(29)
Because ek and ηk have infinite support, to prevent zk from
going beyond the threshold α, the attack sequence δk must
compensate for the term Cek +ηk . By assumption, the attacker
has access to yk = Cxk +ηk (real-time sensor measurements).
Moreover, given its perfect knowledge of the Kalman filter,
the adversary can compute the estimated output C x̂k and then
construct yk −C x̂k = Cek +ηk . For a given chi-squared threshα
old α, define the sequence
:= {δ̄kα ∈p
Rm |(δ̄kα )T δ̄kα ≤ α};
p αδ̄k p
α
α T
, m
,..., m
] and δ̄kα =
for instance, δ̄kα = [ m
√
[ α, 0, . . . , 0]T . Let k = k ∗ denote the starting attack instant
for some k ∗ ≥ 1. Then, for k ≥ k ∗ , it follows that
1
(30)
δk = −Cek − ηk + Σ 2 δ̄kα → zk ≤ α,

7

1

where Σ 2 denotes the symmetric square root matrix of Σ, is a
feasible attack sequence given the capabilities of the attacker.
Sequences δk of the form (30) define a class of attacks that
can be launched by the opponent while preventing the chisquared detector from raising alarms, i.e., zero-alarm attacks.
The estimation error dynamics under the attack (30) is given
by
1

ek+1 = F ek − LΣ 2 δ̄kα + vk , k ≥ k ∗ .

(31)

Remark 5 Note that if ρ[F ] > 1, then kE[ek ]k diverges to
infinity as k grows for any nonstabilizing δ̄kα [22]. That is, zeroalarm attacks of the form (30) may destabilize the system if
ρ[F ] > 1. If ρ[F ] ≤ 1, then kE[ek ]k may or may not diverge
to infinity depending on algebraic and geometric multiplicities
of the eigenvalues with unit modulus of F (a known fact from
stability of LTI systems [22]).
Using the superposition principle of linear systems, the
estimation error ek can be written as ek = evk + eδk , where
evk denotes the part of ek driven by noise and eδk is the part
driven by attacks. Using this new notation, we can write the
dynamics (31) as follows:
(32)
evk+1 = F evk + vk ,
1

(33)
eδk+1 = F eδk − LΣ 2 δ̄kα , k ≥ k ∗ ,
v
δ
with ek∗ = ek∗ and ek∗ = 0. Therefore, the contribution of
zero-alarm attacks to ek is solely determined by eδk generated
by (33). For a sequence sk ∈ Rn , k ∈ N, let s[k∗ ,k] denote the
truncation of sk from k ∗ to k, i.e., s[k∗ ,k] := {sk∗ , . . . , sk }
and ||s[k∗ ,k] || := supk∗ ≤N ≤k ksN k. For any matrix An×n such
that ρ[A] < 1, let k·k∗ denote some matrix norm satisfying
kAk∗ < 1 (such a norm always exists if ρ[A] < 1 [41]).
Proposition 1 Consider the process (3), the Kalman filter
(4)-(8), and the chi-squared procedure (28) with threshold
α ∈ R>0 . Assume ρ[F ] < 1 and let c ∈ R>0 be some
constant satisfying F k ≤ c F k ∗ for all k ∈ N. Let the
α
sensors be attacked by the sequence (30); then, for all δ̄[k
∗ ,k] ,
k > k∗ ∈ N, the trajectories of (33) satisfy the inequalities:
k−k∗

 eδ ≤ γ χ2 := √α c||LΣ 21 || 1 − kF k∗

,

k
k
1 − kF k∗
(34)
√
1

α c||LΣ 2 ||

χ2
δ

 limk→∞ ek ≤ γ̄ :=
.
1 − kF k∗
Proof : The solution of (33), for k > k ∗ , is given by
∗
k−1−k
X
1
α
F i LΣ 2 δ̄k−i−1
eδk = −
,
i=0

it follows that
Pk−1−k∗
1
α
∗
eδk ≤ ||LΣ 2 || i=0
F i ||δ̄[k
∗ ,k−1] ||, k > k .
Because ρ[F ] < 1, there exists a matrix norm, say k·k∗ , such
that kF k∗ < 1 (see Lemma 5.6.10 in [41]). Moreover, because
all norms are equivalent in finite dimensional vector spaces
[41], there exists a constant c ∈ R>0 satisfying kDk ≤ c kDk∗
for all D ∈ Rn×n [45]. It follows that
Pk−1−k∗
1
i
α
kF k∗ ||δ̄[k
eδk ≤ c||LΣ 2 || i=0
∗ ,k−1] ||
∗

1 − kF k∗k−k
α
≤ c||LΣ ||
||δ̄[k
∗ ,k−1] ||,
1 − kF k∗
1
2

OCT 2017

Pn
i
because
i=0 kF k∗ is a geometric series. By construction,
α
for all k ≥ k ∗ , (δ̄kα )T δ̄kα ≤ α which implies ||δ̄[k
∗ ,k−1] || =
√
α
supk∗ ≤N ≤k−1 δ̄N ≤ α; therefore, the estimation error
driven by attacks eδk satisfies the inequalities in (34).

The effect of the attack sequence (30) on the upper bound
2
of the estimation error (34) is determined by the sequence γkχ .
2
The sequence γkχ quantifies the impact of zero-alarm attacks
when the chi-squared detector is used to detect anomalies, i.e.,
2
γkχ gives a measure of the detector performance for the class
of attacks in2 (30) in terms of estimation error deviation. The
sequence γkχ depends on the norm k·k∗ which could be any
matrix norm satisfying kF k∗ < 1. If usual matrix norms (e.g.,
k·k1 , k·k2 , k·k∞ , etc.) do not satisfy this condition, in the
proof of Lemma 5.6.10 in [41], the authors give a procedure
for constructing such a norm provided that ρ(F ) < 1. For
given norm k·k∗ satisfying kF k∗ < 1, the constant c can be
taken as c = inf{c ∈ R>0 : F k − c F k ∗ ≤ 0, ∀ k ∈ N}
which can be obtained numerically.
Next, consider the CUSUM procedure and write (18) in
terms of the estimation error ek :
1
Sk = max(0, Sk−1 + ||Σ− 2 (Cek + ηk + δk )||2 − b), (35)
if Sk−1 ≤ τ ; and Sk = 0, if Sk−1 > τ . As with the chisquared procedure, we look for attack sequences that maintain
the CUSUM statistic below the threshold τ preventing alarms
to be raised. Let the attack start at some k = k ∗ ≥ 2 and
Sk∗ −1 ≤ τ , i.e., the attack does not start immediately after a
false alarm. Define τ̄k := {τ̄k ∈ Rm |τ̄kT τ̄k ≤ τ + b − Sk−1 }
and δ̄kb := {δ̄kb ∈ Rm |(δ̄kb )T δ̄kb ≤ b} for given threshold τ and
bias b. Consider the attack sequence:

1
−Cek − ηk + Σ 2 τ̄k , k = k ∗ ,
(36)
δk =
1
−Cek − ηk + Σ 2 δ̄kb , k > k ∗ .
It follows that Sk∗ = max(0, Sk∗ −1 + τ̄kT τ̄k − b) ≤ τ ,
Sk∗ +1 = max(0, Sk∗ + (δ̄kb ∗ +1 )T δ̄kb ∗ +1 − b) ≤ τ , and
Sk∗ +N = max(0, Sk∗ +N −1 + (δ̄kb ∗ +1 )T δ̄kb ∗ +1 − b) ≤ τ for
all N ∈ N. That is, the class of attack sequences in (36)
prevents the CUSUM procedure from raising alarms. Note
that the attacker can only induce this sequence by exactly
knowing Sk∗ −1 , i.e., the value of the CUSUM sequence one
step before the attack. This is a strong assumption since it
represents a real-time quantity that is not communicated over
the communication network. Even if the opponent has access
to the parameters of the CUSUM (b, τ ) given the stochastic
nature of the residual, the attacker would need to know the
complete history of observations (from when the CUSUM was
started) to be able to reconstruct Sk∗ −1 from data. This is
an inherent security advantage in favor of the CUSUM over
static detectors like the bad-data or chi-squared. Nevertheless,
for evaluating the worst case scenario, we assume that the
attacker has access to Sk∗ −1 . By construction, the estimation
error dynamics under the attack sequence (36) is written as:
1
ek∗ +1 = F e∗k − LΣ 2 τ̄k∗ + vk∗ , and, for k > k ∗ ,
1
(37)
ek+1 = F ek − LΣ 2 δ̄kb + vk .
Note that (31) and (37) have the same dynamics but different
initial condition. Therefore, we may expect upper bounds on
kek k similar to (34) obtained for the chi-squared. Again, we
write ek as ek = evk + eδk , where evk denotes the part of ek

8

driven by noise and eδk is the part driven by attacks of (37).
Using this new notation, we can write the dynamics (37) as:
evk+1 = F evk + vk ,
(38)
1

eδk∗

with
= 0,
F ek∗ + vk∗ .

eδk+1 = F eδk − LΣ 2 δ̄kb , k > k ∗ ,
eδk∗ +1

1
2

= −LΣ τ̄k∗ ,

evk∗

(39)

= ek∗ , and

evk∗ +1

=

Proposition 2 Consider the process (3), the Kalman filter
(4)-(8), and the CUSUM procedure (18) with threshold
τ ∈ R>0 and bias b > b̄ = m ∈ N>0 . Assume ρ[F ] < 1
and let c ∈ R>0 be some constant satisfy F k ≤ c F k ∗
for all k ∈ N. Let the sensors be attacked by the sequence
b
∗
∈ N, the trajectories
(36); then, for all τ̄k∗ and δ̄[k
∗ ,k] , k > k
of (39) satisfy the inequalities:

k−k∗
√
1 − kF k∗

1

CS
δ
2

ek ≤ γk := b c||LΣ ||


1 − kF k∗


k−k∗ −1

1







 limk→∞ eδk

+ c||LΣ 2 τ̄k∗ || kF k∗
√
1
b c||LΣ 2 ||
CS
.
≤ γ̄ :=
1 − kF k∗

,

(40)

Proof : The solution of (39), for k > k ∗ + 1, is given by
∗
k−2−k
X
1
1
δ
k−k∗ −1
b
ek = −F
LΣ 2 τ̄k∗ −
F i LΣ 2 δ̄k−i−1
.
i=0

Then, following the same lines as in the proof of Proposition
1, we can write the following
∗
1
−1
eδk ≤ c||LΣ 2 τ̄k∗ || kF kk−k
∗
∗

1 − kF k∗k−k
b
||δ̄[k
+ c||LΣ ||
∗ +1,k−1] ||,
1 − kF k∗
By construction, for all k > k ∗ , (δ̄kb )T δ̄kb ≤√b which implies
b
b
≤ b; therefore, the
||δ̄[k
δ̄N
∗ +1,k−1] || = supk∗ +1≤N ≤k−1
trajectories of the estimation error dynamics (39) satisfy the
inequalities in (40).

1
2

When considering the CUSUM, the effect of attacks of
the form (36) on the upper bound of the estimation error is
determined by the sequence γkCS in (40). Note that, in steady
state, the (τ̄k∗ )-dependent term in γkCS decreases exponentially
2
to zero. Then, the sequences γkχ and γkCS take constant values
asymptotically. It follows that we can directly quantify the
performance ratio (in terms of steady state deviation of eδk )
between the two detectors. This is stated in the following
corollary of Proposition 1 and Proposition 2.
Corollary 1 In steady state:
2

2

γkχ
γ̄ χ
= CS =
CS
k→∞ γ
γ̄
k
lim

r

α
.
b

(41)

B. Detector Comparison
For the case of zero-alarm attacks, we have derived upper
bounds on eδk for both the chi-squared and CUSUM procedures provided that ρ[F ] < 1 (otherwise eδk diverges
under the class of zero-alarm attacks considered). To compare
2
these results, we use the ratio of sequences (γkχ /γkCS ). From
p
2
Corollary 1, limk→∞ (γkχ /γkCS ) = α/b. That is, in steady
state, the (τ̄k∗ )-dependent term in γkCS exponentially decreases

OCT 2017

9

2500

Sk

2000
1500

b = 0.85b̄
b = 0.95b̄
b = 1.05b̄

1000
500
0
0

1000

2000

3000

4000

5000

k
2

Fig. 2: CUSUM evolution for different values of bias b.

Fig. 4: Asymptotic ratio (γ̄ χ /γ̄ CS ) versus the false alarm rate
A for different values of CUSUM bias b.
use this system to demonstrate our results. The state, input,
and 
output vectors
of the considered reactor are:
 

 

C0
C0
Cu
 

 

 T0 
x(t) :=   , u(t) :=  Tu  , y(t) :=  T0  ,

 Tw 


Tw
Tw,u

Tm





where

2

Fig. 3: Upper bounds γkχ and γkCS , and deviation of eδk due
to the zero-alarm attacks (δ̄kα,1 , δ̄kb,1 ) and (δ̄kα,2 , δ̄kb,2 ). Attacks
are induced at k = 2 × 103 .
to zero and, if b < α, under the same class of zero-alarm
attacks, the CUSUM procedure leads to smaller steady state
deviations on kek k than the chi-squared procedure. In general,
to increase the chances of attack detection, it is desired to
select b as close as possible to b̄ in Theorem 1. It follows
that b ≈ b̄ = m. On the other hand, according to Theorem
∗
3, α must be selected as α = α∗ = 2P−1 ( m
2 , 1 − A ) to
∗
fulfill a desired false alarm rate A . In this case, we want
to select A∗ close to zero, such that there are only a few
false alarms. Let A∗ ∈ {0.01, 0.1} and m = 2, i.e., false
alarms between 1% and 10% and two dimensional outputs;
∗
then, α = 2P−1 ( m
2 , 1−A ) ∈ [4.60, 9.21] and b ≈ b̄ = 2. This
implies that, in steady state, for the same class of attacks and
A∗ ∈ [0.01, 0.1], the chi-squared procedure leads to at least
two times larger upper bounds than the CUSUM. Actually,
2
for having α = b (which implies limk→∞ (γkχ /γkCS ) = 1),
it is necessary to allow for a rate of A∗ = 0.63, which is
high for practical purposes. For the CUSUM procedure, the
threshold τ is selected to fulfill the desired A∗ . Given that
there are no exact closed-form expressions to relate τ and
A∗ (we provided a numeric approximation), it is not possible
to exactly tell how large the τ would need be to satisfy
A∗ . However, as mentioned, the contribution of τ̄k∗ to ek
vanishes exponentially, i.e., independent of how large τ is,
its contribution to γkCS is zero in steady state.
VI. S IMULATION E XPERIMENTS
The authors in [18],[19] study the fault detection problem
for a well stirred chemical reactor with heat exchanger. We

C0




T0





 Tw
Tm



Cu




 Tu


Tw,u

: Concentration of the chemical product,
: Temperature of the product,
: Temperature of the jacket water of heat exchanger,
: Coolant temperature,
: Inlet concentration of reactant,
: Inlet temperature,
: Coolant water inlet temperature.

We linearize the nonlinear model introduced in [18] about
the origin x(t) = 04×1 and then discretize it with sampling
time h = 0.05. The resulting discrete-time linear system is
given by (3)-(8) with matrices as given in (42). The original
model in [18] does not consider sensor/actuator noise, we have
included noise to increase the complexity of our simulation
experiments. First, assume no attacks, i.e., δk = 0, and
consider the CUSUM procedure (18) with distance measure
zk = rkT Σ−1 rk and residual sequence (10). According to
Theorem 1, the bias b must be selected larger than b̄ = m = 3
to ensure mean square boundedness of Sk independent of the
threshold τ . Figure 2 depicts the evolution of the CUSUM
for b ∈ {0.85b̄, 0.95b̄, 1.05b̄1} and k ∈ [1, 5000]. For the
purpose of illustrating this unbounded growth, we have omitted
the reset procedure of the CUSUM. Note that the bound
for b is tight, small deviations from b̄ lead to (boundedness)
unboundedness of Sk . Next, for desired false alarm rates A∗ ∈
{0.25, 0.10, 0.02}, we compute the corresponding thresholds
τ = τ ∗ using Theorem 2 and Remark 4. For these thresholds,
in Table 1, we present the actual false alarm rate A (obtained
by simulation) and the desired A∗ . Note that the difference
between A and A∗ is less that 0.05 in all cases.
In Figure 3, we present the evolution of eδk when both the
chi-squared and the CUSUM are deployed for attack detection
and the attack sequences δk are zero-alarm attacks of the form
introduced in (30) and (36), respectively.
We consider two
p
α,1
α
attack
sequences,
first,
δ̄
=
δ̄
=
α/m
δ̄1 , δ̄kb = δ̄kb,1 =
k
k
p
p
b/mδ̄1 , τ̄k∗ = τ̄k1∗ =
(τ + b − Sk−1 )/mδ̄1 , and δ̄1 =
√
1
.
The
second
attack
is
δ̄kα = δ̄kα,2 = αδ̄2 , δ̄kb = δ̄kb,2 =
m×1
√
p
bδ̄2 , and τ̄k∗ = τ̄k2∗ = τ + b − Sk−1 δ̄2 , where δ̄2 denotes

OCT 2017

10

A∗ = 0.25
τ∗

b/b̄

τ =

1.05
1.15
2.00

1.0282
0.6872
−

A∗ = 0.10

A (Simul.)
0.2041
0.2010
−

τ =

τ∗

A∗ = 0.02

A (Simul.)

τ = τ∗

A (Simul.)

0.0899
0.0885
0.0953

12.3208
10.0327
4.1002

0.0196
0.0184
0.0202

3.9602
3.3699
0.2528

Table 1. Simulation Experiments. Results from Theorem 2 and Remark 1.

dure for identifying compromised sensors. In particular, steady
state Kalman filters have been proposed to estimate the state
of the physical process; then, these estimates have been used
to construct residual variables (between sensor measurements
and estimations) which drive the CUSUM procedure. Using
stability results for stochastic systems and Markov chain
approximations of the CUSUM sequence, we have derived
systematic tools for tuning the CUSUM procedure such that
mean square boundedness of the CUSUM sequence is guaranteed and the desired false alarm rate is fulfilled. For a class
of zero-alarm attacks, we have characterized the performance
of the proposed CUSUM procedure in terms of the effect
that the attack sequence can induce on the system dynamics.
Then, we have compared this performance against the one
obtained using chi-squared procedures. For the linearized
model of the chemical reactor considered in [19], [18], by
means of a simulation study, we have showed that our tools are
useful for tuning the CUSUM procedure and provide accurate
predictions about the performance of the detection scheme.

the unitary singular vector corresponding to the largest singular
1
value of (I − F )−1 LΣ 2 . It can be proved that this selection of
δ̄2 maximizes the steady state value of eδk . For the CUSUM,
we select b = 2b̄ = 6 and τ = τ ∗ = 4.1002 such that A ≈
A∗ = 0.02 (see Table 1). Likewise, we select α = α∗ =
2P−1 ( 22 , 1 − 0.02) = 9.83 such that, according to Theorem 3,
A = A∗ = 0.02. The attacks are induced at k = k ∗ = 2×103.
Note that, as stated in Proposition 1 and Proposition 2, given
that ρ[F ] < 1, inequalities (34) and (40) are satisfied for both
attacks. Moreover, as mentioned in Section V-B, we expect
that the CUSUM leads to smaller steady state deviation on
kek k because b < α and δ̄kα,i = aδ̄kb,i , i = 1, 2 for some
a ∈ R>0 . This is exactly what we see in Figure 3. Note that
the ratio γ̄χ2 /γ̄CS = 1.28 is fixed by our choice of false alarm
rate and bias. In Figure 4, we depict the evolution of the ratio
γ̄χ2 /γ̄CS versus the false alarm rate A for different values of
CUSUM bias b.
VII. C ONCLUSIONS
In this paper, for a class of stochastic linear time-invariant
systems, we have characterized a model-based CUSUM proce








F




















L












=






=



0.8353

0

0

0

0.8324

0

0
0

0.0001
0.0280

0.1633
0.0172

0.8271

0

0

0

0.8243

0
0

0.0002
0.0481

0








0.0031
G=


0 

0.9320


0.0458

0

0

0.0457

0
0

0
0.0007

0




1


0 

,
C
=
0

0.0231

0
0.0006


1.0169

0.0002
 , R2 = 0.01 × I3 , R0 = R1 = I4 , Σ = 
 0

0.1619
0
0.0543

R EFERENCES
[1] F. Pasqualetti, F. DÃűrfler, and F. Bullo, “Attack detection and identification in cyber-physical systems,” IEEE Transactions on Automatic
Control, vol. 58, pp. 2715–2729, 2013.
[2] Y. Mo, E. Garone, A. Casavola, and B. Sinopoli, “False data injection
attacks against state estimation in wireless sensor networks,” in Decision
and Control (CDC), 2010 49th IEEE Conference on, 2010, pp. 5967–
5972.
[3] C. Kwon, W. Liu, and I. Hwang, “Security analysis for cyber-physical
systems against stealthy deception attacks,” in American Control Conference (ACC), 2013, 2013, pp. 3344–3349.
[4] F. Miao, Q. Zhu, M. Pajic, and G. J. Pappas, “Coding sensor outputs
for injection attacks detection,” in Decision and Control (CDC), 2014
IEEE 53rd Annual Conference on, 2014, pp. 5776–5781.
[5] C.-Z. Bai and V. Gupta, “On kalman filtering in the presence of a
compromised sensor: Fundamental performance bounds,” in American
Control Conference (ACC), 2014, 2014, pp. 3029–3034.
[6] C. Z. Bai, F. Pasqualetti, and V. Gupta, “Security in stochastic control
systems: Fundamental limitations and performance bounds,” in American Control Conference (ACC), 2015, 2015, pp. 195–200.

0
1.0169
0.0001

0

0

1
0

0
1


0

0 ,
0

(42)


0

0.0001 .
1.0105

[7] A. Cárdenas, S. Amin, Z. Lin, Y. Huang, C. Huang, and S. Sastry, “Attacks against process control systems: Risk assessment, detection, and
response,” in Proceedings of the 6th ACM Symposium on Information,
Computer and Communications Security, 2011, pp. 355–366.
[8] C. Murguia and J. Ruths, “Cusum and chi-squared attack detection of
compromised sensors,” in proceedings of the IEEE Multi-Conference on
Systems and Control (MSC), 2016.
[9] ——, “Characterization of a cusum model-based sensor attack detector,”
in proceedings of the 55th IEEE Conference Decision and Control
(CDC),, 2016.
[10] C. Murguia, N. van de Wouw, and J. Ruths, “Cusum and chi-squared
attack detection of compromised sensors,” in proceedings of the 20th
World Congress of the International Federation of Automatic Control
(IFAC), accepted, 2017.
[11] M. Ross, Introduction to Probability Models, Ninth Edition. Orlando,
FL, USA: Academic Press, Inc., 2006.
[12] H. Nijmeijer and A. van der Schaft, Nonlinear dynamical control
systems. New York: Springer, 1990.
[13] Z. Guo, D. Shi, K. H. Johansson, and L. Shi, “Optimal linear cyberattack on remote state estimation,” IEEE Transactions on Control of
Network Systems, vol. 4, no. 1, pp. 4–13, 2017.

OCT 2017

[14] F. Gustafsson, Adaptive Filtering and Change Detection. West Sussex,
Chichester, England: John Wiley and Sons, LTD, 2000.
[15] K. Cheolhyeon, Y. Scott, and H. Inseok, “Real-time safety assessment
of unmanned aircraft systems against stealthy cyber attacks,” Journal of
Aerospace Information Systems, vol. 13, pp. 27–45, 2015.
[16] D. I. Urbina, J. A. Giraldo, A. A. Cardenas, N. O. Tippenhauer,
J. Valente, M. Faisal, J. Ruths, R. Candell, and H. Sandberg, “Limiting the impact of stealthy attacks on industrial control systems,” in
Proceedings of the 2016 ACM SIGSAC Conference on Computer and
Communications Security, ser. CCS ’16, 2016, pp. 1092–1105.
[17] D. I. Urbina, J. Giraldo, A. A. Cardenas, J. Valente, M. Faisal, N. O.
Tippenhauer, J. Ruths, and H. Sandberg, “Survey and new directions
for physics-based attack detection in control systems,” in NIST Standard
Reference Simulation Website, NIST Standard Reference Database Number 16-010, National Institute of Standards and Technology, Gaithersburg MD, 20899, 2016.
[18] J. Chen and R. Patton, Robust Model-Based Fault Diagnosis for Dynamic
Systems. Springer Publishing Company, Incorporated, 2012.
[19] K. Watanabe and D. M. Himmelblau, “Fault diagnosis in nonlinear
chemical processes. part ii. application to a chemical reactor,” AIChE
Journal, vol. 29, 1983.
[20] T. T. Ha, Theory and Design of Digital Communication Systems.
Cambridge University Press, 2010.
[21] L. Ljung and T. Glad, Modeling of Dynamic Systems. Englewood Cliffs:
PTR Prentice Hall, 1994.
[22] K. J. Aström and B. Wittenmark, Computer-controlled Systems (3rd Ed.).
Upper Saddle River, NJ, USA: Prentice-Hall, Inc., 1997.
[23] B. Anderson and J. Moore, Optimal Filtering. Englewood Cliffs, NJ:
Prentice-Hall, 1979.
[24] A. Wald, “Sequential tests of statistical hypotheses,” Ann. Math. Statist.,
vol. 16, pp. 117–186, 1945.
[25] A. Willsky, “A survey of design methods for failure detection in dynamic
systems,” Automatica, vol. 12, pp. 601 – 611, 1976.
[26] E. Page, “Continuous inspection schemes,” Biometrika, vol. 41, pp. 100–
115, 1954.
[27] M. Basseville, “Detecting changes in signals and systems - a survey,”
Automatica, vol. 24, pp. 309 – 326, 1988.
[28] J. Gertler, “Survey of model-based failure detection and isolation in
complex plants,” Control Systems Magazine, IEEE, vol. 8, pp. 3–11,
1988.
[29] R. Agniel and E. Jury, “Almost sure boundedness of randomly sampled
systems,” SIAM Journal on Control, vol. 9, pp. 372–384, 1971.
[30] T. Tarn and R. Yona, “Observers for nonlinear stochastic systems,”
Automatic Control, IEEE Transactions on, vol. 21, pp. 441–448, 1976.
[31] C. van Dobben de Bruyn, Cumulative sum tests : theory and practice.
London : Griffin, 1968.
[32] B. Adams, W. Woodall, and C. Lowry, “The use (and misuse) of
false alarm probabilities in control chart design,” Frontiers in Statistical
Quality Control 4, pp. 155–168, 1992.
[33] R. A. Khan, “Wald’s approximations to the average run length in cusum
procedures,” Journal of Statistical Planning and Inference, vol. 2, pp.
63 – 77, 1978.
[34] C. Park, “A corrected wiener process approximation for cusum arls,”
Sequential Analysis, vol. 6, pp. 257–265, 1987.
[35] M. Reynolds, “Approximations to the average run length in cumulative
sum control charts,” Technometrics, vol. 17, pp. 65–71, 1975.
[36] C. Champ and S. Rigdon, “A a comparison of the markov chain and the
integral equation approaches for evaluating the run length distribution of
quality control charts,” Communications in Statistics - Simulation and
Computation, vol. 20, pp. 191–204, 1991.
[37] D. A. E. D. Brook, “An approach to the probability distribution of cusum
run length,” Biometrika, vol. 59, no. 3, pp. 539–549, 1972.
[38] A. Luceno and J. Puig-pey, “Evaluation of the run-length probability
distribution for cusum charts: Assessing chart performance,” Technometrics, vol. 42, pp. 411–416, 2000.
[39] W. Woodall, “The distribution of the run length of one-sided cusum
procedures for continuous random variables,” Technometrics, vol. 25,
pp. 295–301, 1983.
[40] S. Meyn and R. Tweedie, Markov Chains and Stochastic Stability.
Springer-Verlag, 1993.
[41] R. A. Horn and C. R. Johnson, Matrix Analysis, 2nd ed. New York,
NY, USA: Cambridge University Press, 2012.
[42] T. Dayar, “On moments of discrete phase-type distributions,” in Formal
Techniques for Computer Systems and Business Processes, M. Bravetti,
L. Kloul, and G. Zavattaro, Eds. Versailles, France: Springer, 2005,
pp. 51–63.

11

[43] “On characterizations of the input-to-state stability property,” Systems
and Control Letters, vol. 24, pp. 351 – 359, 1995.
[44] “Input-to-state stability for discrete-time nonlinear systems,” Automatica, vol. 37, pp. 857 – 869, 2001.
[45] G. Stewart, Matrix Algorithms. Society for Industrial and Applied
Mathematics, 2001.

A PPENDIX A
P ROOF OF T HEOREM 1


Define the functions Vk := Sk2 and ∆Vk := ESk Vk+1 −
Vk . Along (18) with distance measure zk and independent of
τ , we have that


(43)
∆Vk = ESk (Sk + zk+1 − b)+2 − Sk2 ,
+2
2
where ρ := max(0, ρ) . Note that, by construction, Sk , zk ∈
R≥0 for all k ∈ N. First, consider Sk ∈ [b, ∞), it follows that
(Sk + z − b)+2 = (Sk + z − b)2 ,
for all z ∈ R≥0 ; therefore



(44)
ESk Vk+1 = ESk (Sk + zk+1 − b)2 ],
for Sk ∈ [b, ∞) and zk+1 ∈ R≥0 . Next, consider Sk ∈ [0, b)
which implies Sk − b < 0, then
+2
z ∈ [0, b − Sk ] ⇒ Sk + z − b
= 0,
+2
2
z ∈ (b − Sk , ∞] ⇒ Sk + z − b
= Sk + z − b .
It follows that



(45)
ESk Vk+1 ≤ ESk (Sk + zk+1 − b)2 ],
for Sk ∈ [0, b) and zk+1 ∈ R≥0 . Using (44), (45), and
independence between Sk and zk+1 , we have that



(46)
ESk Vk+1 ≤ ESk (Sk + zk+1 − b)2 ]

2
= (Sk − b)2 + 2(Sk − b)E[zk+1 ] + E[zk+1
],
for Sk , zk+1 ∈ R≥0 . Using (16) and the relation:
var[zk ] = E[zk2 ] − E[zk ]2 ,
inequality (46) amounts to


ESk Vk+1 ≤ (Sk − b)2 + 2(Sk − b + 1)m + m2 , (47)
and, therefore,


∆Vk = ESk Vk+1 − Sk2 ,

≤ −2(b − m)Sk + (b − m)2 + 2m,
(48)
for Sk , zk+1 ∈ R≥0 . From (48), given that b ∈ R>0 and
Sk ∈ R>0 by construction, it is easy to verify that
∆Vk < 0 ⇔ b > b̄ := m and Sk > S̄, S̄ :=

(b−m)2 +2m
2(b−m) .

(49)

Therefore, b ∈ (b̄, ∞) implies:

∆Vk < 0 for Sk ∈ (S̄, ∞),
(50)
∆Vk ≥ 0 for Sk ∈ [0, S̄].


Recall that ∆Vk := ESk Vk+1 − Vk . Assume that for some
k = k ∗ ∈ N, Sk∗ ∈ (S̄, ∞); then, from (50), ∆Vk∗ < 0 and
consequently
(51)
ESk∗ [Vk∗ +1 ] < Vk∗ .
Next, for k = k ∗ + 1, let Sk∗ +1 ∈ (S̄, ∞), then
ESk∗ +1 [Vk∗ +2 ] < Vk∗ +1 .
Using (51), (52), and the property
ESk∗ [Vk∗ +2 ] = ESk∗ [ESk∗ +1 [Vk∗ +2 ]],
we have
ESk∗ [Vk∗ +2 ] < ESk∗ [Vk∗ +1 ] < Vk∗ .

(52)
(53)
(54)

OCT 2017

12

Continuing this way, we obtain
(55)
Vk∗ > ESk∗ [Vk∗ +1 ] > . . . > ESk∗ [Vk∗ +n ],
∗ ∗
∗
for n ≥ 2 and k = {k , k + 1, . . . , k + n} such that
Sk ∈ (S̄, ∞). Therefore, from (55), it can be concluded that
the second moment ESk∗ [Vk ] = ESk∗ [Sk2 ] decreases monotonically for Sk ∈ (S̄, ∞) and b ∈ (b̄, ∞); and consequently,
ESk∗ [Sk2 ] < ∞. Next, assume that for some k = k ∗ ∈ N,
Sk∗ ∈ [0, S̄] and b ∈ (b̄, ∞); then, from (50), ∆Vk∗ ≥ 0 and,
by (48), it follows that

(56)
∆Vk∗ ≤ −2 b − m Sk∗ + (b − m)2 + 2m ≥ 0.
Since Sk∗ ∈ [0, S̄], there always exist constants a ∈ (0, 1) and
β ∈ R>0 satisfying

∆Vk∗ ≤ −2 b − m Sk∗ + (b − m)2 + 2m
(57)
≤ −aVk∗ + β, for Sk∗ ∈ [0, S̄].


Given that ∆Vk = ESk Vk+1 − Vk , by (57), we have
ESk∗ [Vk∗ +1 ] ≤ (1 − a)Vk∗ + β, for Sk∗ ∈ [0, S̄]. (58)
Next, for k = k ∗ + 1, let Sk∗ +1 ∈ [0, S̄], then
(59)
ESk∗ +1 [Vk∗ +2 ] ≤ (1 − a)Vk∗ +1 + β.
By (58), (59), and the property
(60)
ESk∗ [Vk∗ +2 ] = ESk∗ [ESk∗ +1 [Vk∗ +2 ]],
we have
ESk∗ [Vk∗ +2 ] ≤ ESk∗ [(1 − a)Vk∗ +1 + β]
= (1 − a)ESk∗ [Vk∗ +1 ] + β
≤ (1 − a)2 Vk∗ + (1 − a)β + β. (61)
Continuing this way, we obtain
n−1
X
(1 − a)i ,
(62)
ESk∗ [Vk∗ +n ] ≤ (1 − a)n Vk∗ + β
i=0

for n ≥ 1 and k = {k ∗ , k ∗ + 1, .P
. . , k ∗ + n} such P
that Sk ∈
k−1
∞
[0, S̄]. Given that a ∈ (0, 1) and i=0 (1 − a)i ≤ i=0 (1 −
a)i = a1 , then, as n → ∞, it is satisfied that
β
2
(63)
ESk∗ [S∞
]≤ .
a
Therefore, by (62) and (63), the second moment ESk∗ [Sk2 ]
does not grow unbounded for Sk ∈ [0, S̄] and b ∈ (b̄, ∞),
i.e., ESk∗ [Sk2 ] < ∞. So far, combining the preliminary results
presented above, we have proved that for all Sk ∈ [0, b] ∪
(b, ∞), the second moment is finite and either decreasing or
uniformly bounded in k ∈ N provided that the conditions
of Theorem 1 are satisfied. Hence, for the residual sequence

rk ∼ N (0, Σ), b > b̄ → ES1 [Sk2 ] < ∞ for all k ∈ N.

