1

Colonel Blotto Game for Secure State Estimation in
Interdependent Critical Infrastructure
Aidin Ferdowsi∗ , Walid Saad∗ , Narayan B. Mandayam†
Wireless@VT, Bradley Department of Electrical and Computer Engineering,
Virginia Tech, Blacksburg, VA, USA, Emails: {aidin,walids}@vt.edu
WINLAB, Dept. of ECE, Rutgers University, New Brunswick, NJ, USA, Email: narayan@winlab.rutgers.edu
∗

arXiv:1709.09768v1 [cs.SY] 28 Sep 2017

†

Abstract—Securing the physical components of a city’s interdependent critical infrastructure (ICI) such as power, natural
gas, and water systems is a challenging task due to their
interdependence and large number of involved sensors. Using
a novel integrated state-space model that captures the interdependence, a two-stage cyber attack on ICI is studied in which the
attacker first compromises the ICI’s sensors by decoding their
messages, and, subsequently, it alters the compromised sensors’
data to cause state estimation errors. To thwart such attacks,
the administrator of the CIs must assign protection levels to
the sensors based on their importance in the state estimation
process. To capture the interdependence between the attacker and
the ICI administrator’s actions and analyze their interactions, a
Colonel Blotto game framework is proposed. The mixed-strategy
Nash equilibrium of this game is derived analytically. At this
equilibrium, it is shown that the administrator can strategically
randomize between the protection levels of the sensors to deceive
the attacker. Simulation results coupled with theoretical analysis
show that, using the proposed game, the administrator can
reduce the state estimation error by at least 50% compared to
any non-strategic action. The results also show that the ICI’s
administrator must consider the CIs inside a city as a unified ICI
for security analysis instead of assigning independent protection
levels to each individual CI, as is conventionally done.

I. I NTRODUCTION
The services delivered by a smart city’s critical infrastructure (CI) such as power, natural gas, and water will be highly
interdependent [1]–[4]. CIs are cyber-physical systems (CPSs)
that encompass physical infrastructure whose performance
is monitored and controlled by a cyber system, typically
consisting of a massive number of sensors. These CPSs
exhibit close interactions between their cyber and physical
components [4]–[6]. The state estimation of the CIs, which
uses cyber elements to monitor the physical elements, is a
crucial stage for controlling their functionality. However, the
interdependence between CIs and the high synergy between
their physical and cyber components make them vulnerable to
attacks and failures [7]–[10].
Numerous solutions have been presented for securing state
estimation of CPSs as well as for CI failure detection and identification [11]–[14]. In [11], the authors presented a controltheoretic approach for attack detection and identification in
noiseless environments using centralized and distributed attack
detection filters. The works in [12]–[14] considered the estimation of a CPS under stealthy deception and replay cyber-attacks
using a Kalman filter. Moreover, the security of interdependent
critical infrastructure (ICI) has been studied in recent works
such as [15]–[17]. In [15], the authors assessed the security of
interdependent power and natural gas systems under multiple
hazards, considering the ICI’s performance as a measurement
for security. In [16], the authors proposed an agent-based
This research was supported by the U.S. National Science Foundation under
Grants OAC-1541105, OAC-1541069, and EAGER-1745829.

model to capture the effects of interdependencies and quantify
the coupling strength between ICI. Also, the impact of natural
and human-included disasters been studied in [17].
Furthermore, the security and protection of sensor networks,
which collect data from CIs, has been studied in [18]–[20]. In
[18], the authors proposed a novel method for physical attack
protection with human virtualization in the context of data
centers using sensors that detect an impending physical/human
attack and, then, alarms to mitigate the attack. The authors in
[19], introduced a learning algorithm to extract features from
the sensor messages to detect the cyber attacks. The work in
[20] proposed a distributed observer for state estimation of
CIs in lossy sensor networks with cyber attacks. The works in
[21]–[23] used game-theoretic tools to study the interactions
between the defender of a single CPS and an attacker that
seeks to compromise the various nodes of a CPS.
However, the works in [11]–[20] do not consider the limitations of the available security resources for the protection,
detection, and identification of CI attacks. For instance, in
practical smart cities, resource limitations may substantially
affect the security of the CIs. In a city, to prevent an attacker
from breaking into the sensors of an ICI, numerous methods
can be adopted such as encryption of sensor data, implementation of attack detection filters, or periodic monitoring
algorithms. However, because of massive data transmission
from sensors to the central processing unit, such security
solutions will require a large number of computation, a high
communication bandwidth, or a considerable amount of financial resources, all of which constitute limited resources
for the ICI’s administrator. Therefore, unlike the idealized
security solutions in [11]–[20], due to resource limitations,
the administrator of an ICI has to prioritize between the
protection of the cyber components of ICI based on their
importance in the state estimation process [24]. Another key
limitation in the current literature is that the majority of
the existing works, such as [11]–[14] and [21]–[23] do not
take into account the interdependence between the CPSs.
Meanwhile, those that account for interdependencies such as in
[15]–[17] are mostly based on graph-theoretic constructs that
abstract much of the functionalities of the CIs. In practice, the
functionalities of CIs in a city are largely interdependent and
cannot be simply captured by a graph. For instance, power
generation in generators which are supplied by natural gas
requires instantaneous natural gas transmission from natural
gas CI. Therefore, the sensors of each CI will collect the
physical measurements of interdependent CIs and, to protect
them from cyber attacks, the ICI administrator must consider
realistic models for the CI interdependencies.

2

The main contribution of this paper is a novel gametheoretic framework for analyzing and optimizing the security
of a large-scale ICI’s state estimation. To build this unified
security framework, this paper makes several contributions:
•

•

•

•

We first introduce a novel integrated state-space model
that captures the dynamics of an ICI consisting of power,
natural gas, and water distribution systems. For enabling
state estimation of the proposed ICI dynamics model, we
implement a centralized Kalman filter that uses the sensor
measurements collected to estimate the ICI’s state.
We consider a two-stage cyber attack that targets the
sensors of the ICI so as to manipulate the state estimation.
In the first stage, the attacker aims to compromise the
ICI sensors by breaking their protection algorithm (e.g.,
end-to-end encryption or sensor attack detection filter). In
the second stage, the attacker manipulates the ICI’s state
estimation by altering the compromised sensors’ data to
induce state estimation errors. To defend against such
attacks and protect the sensors, the administrator of the
ICI can assign different protection levels to each group
of sensors. The protection levels can vary across sensors
for two reasons: a) the available resources which are used
in the protection procedure are limited and b) the state
estimation sensors have different importance levels.
Since the actions of the attacker and the defender are
interdependent, we propose a Colonel Blotto game framework [25] to analyze the interactions between the attacker
and the administrator. In this game, the attacker chooses
the set of sensors to compromise while the administrator
assigns protection levels to the sensors. In contrast to
existing works on Colonel Blotto for CPS security [21]–
[23], our game considers the interdependence between
multiple CPSs. For this game, we derive the mixedstrategy Nash equilibrium for the administrator and the
attacker as a function of their available resources and the
maximum state estimation error due to the attack.
We simulate a cyber attack to a power-gas-water ICI in
which a sensor network collects data from the physical
components of the ICI. Simulation results show that the
mixed strategy for the administrator increases the protection of the cyber system of a large-scale ICI and reduces
the estimation error of the ICI by at least a factor of 50%
compared to a baseline case in which the administrator
of the ICI does not randomize between the assignment
of protection levels. Theoretical and simulation results
also show that the ICI administrator must consider all the
CIs as a unified system while assigning the protection
levels to the sensors, rather than separately protecting
them without accounting for their interdependence.

The rest of the paper is organized as follows. Section
II presents the dynamic ICI model and studies the state
estimation using Kalman filtering. In Section III, we analyze
the maximum reachability of the estimation error. The gametheoretic framework is discussed in Section IV. In Section
V, we present and analyze the simulation results. Finally,
conclusions are drawn in Section VI.

II. I NTERDEPENDENT C RITICAL I NFRASTRUCTURE AND
ATTACK M ODEL
Consider an ICI as a CPS whose physical system consists of
three interdependent power, natural gas, and water distribution
CIs and whose cyber system is a network of sensors that collect
data from the physical components of the CIs and transmit it
to a central processing unit. We first derive a state-space model
for the physical system of each CI separately and then present
the general ICI model. Finally, we discuss the associated cyber
system and its vulnerability to attacks.
A. Physical System
The power system can be modeled as a linear dynamic
system whose inputs are the electrical power demands from
the load buses [26] and [27]. We focus only on generators
that are supplied by natural gas [28]and we consider water as
a requirement for the vapor condensation and cooling down
in some of the generators [29]. Then, the dynamic model for
a power infrastructure can be given by:
ẋe (t) = Ae xe (t) + B e ue (t),

(1)
y eg (t) = C eg xe (t), y ew (t) = C ew xe (t),
where Ae is the power CI state matrix, B e is the power
demand matrix, and C eg and C ew are the matrices of natural
gas demand and water demand of power CI, respectively.
Moreover, xe , ue , y eg , and y ew are vectors that capture power
state variables, power demand, natural gas demand of power
CI, and water demand of power CI.
The natural gas and water CIs are designed to supply natural
gas and water to consumers in a city. Due to the pressure loss
at the junctions of these two CIs, gas compressors and water
pumps are used to compensate the pressure loss. We can thus
capture the operation of the natural gas CI using the following
dynamic model [30]:
ẋg (t) = Ag xg (t) + B g ug (t),

(2)
y ge (t) = C ge xg (t) + D g ug (t),
where Ag is the matrix of natural gas CI state variables, B g is
the natural gas demand matrix, C g is the power demand matrix
of the natural gas CI, and D g is the matrix of relationship
between the natural gas output and power demand of natural
gas CI. xg , ug , and y ge are the vectors of natural gas CI state
variables, natural gas demand, and power demand of natural
gas CI. Similarly, for the water CI, we have [31]:
ẋw (t) = Aw xw (t) + B w uw (t),

(3)
y we (t) = C we xw (t) + D w uw (t),
where xw , uw , and y we are, respectively, the vectors of water
CI state variables, water demand, and power demand of water
CI. Aw is the state matrix of water CI, B g is the water demand
matrix, C g is the power demand matrix of water CI, and D g is
the matrix of relationship between the water output and power
demand of water CI. Note that, we derived all the matrices in
(1), (2), and (3) and the interdependence of CIs as summarized
in Appendix A.
Considering the interdependence between power-gas-water
CIs, we can derive a unified model for each CI as follows:
ẋe (t) = Ae xe (t) + B e (ue (t) + T ge y ge (t) + T we y we (t)) ,

3

ẋg (t) = Ag xg (t) + B g (ug (t) + T eg y eg (t)) ,
ẋw (t) = Aw xw (t) + B w (uw (t) + T ew y ew (t)) ,

(4)

where
y eg (t) = C eg xe (t), y ew (t) = C ew xe (t),
y ge (t) = C ge xg (t) + D g (ug (t) + T eg y eg (t)) ,

(5)
y we (t) = C we xw (t) + D w (uw (t) + T ew y ew (t)) .
Here, T ge ,T we ,T ew , and T eg are matrices connecting the
inputs and outputs of the three CIs whose elements are equal
to one if the output of one CI is connected to the input of
another CI or is equal to zero otherwise. By substituting (5)
into (4), we will have the following state-space model for the
interdependent critical gas-power-water infrastructure:
ẋ(t) = Ax(t) + Bu(t),

(6)

where

 e

xe (t)
u (t)
x(t) =  xg (t)  , u(t) =  ug (t)  ,
(7)
xw (t)
uw (t)
and A and B are defined in (8). (6) captures the dynamics
of an ICI. In this model, the state variables of three CIs are
mutually interdependent, and changes in one CI can affect the
other two CIs.
In Fig. 1, we present an illustrative example of an ICI that
can show the interdependence between the state variables of
three CIs. In this example, we consider 10 generators out of
which 6 are supplied by natural gas and 4 require water flow to
operate. Here, 11 natural gas pipelines, and 11 water pipelines
distribute natural gas and water to the demand junctions. Based
on this example we find the matrices A and B in (6) and
simulate the ICI. To illustrate how the changes in one CI can
affect the state variables of other CIs, we increase the power
demand in generator 5, ue5 , at time t = 20.5. Fig. 2 shows the
change of state variables of the natural gas pipeline between
junctions 2 and 3 and state variables of water pipeline between
the junctions 2 and 3. The reason is that any increase in power
demand results in increase of electric power generation, and
due to the interdependence between electric power generation
and the consumption of the natural gas and water, the state
variables of the natural gas and water CIs change.


B. Cyber System
To monitor the state variables in (6), a cyber system is
needed. For the considered ICI , the cyber system will consist
of a number of sensors spread around the ICI and collecting
different measurements from the ICI’s components. Sensors
and meters in the power infrastructure measure the instantaneous frequency of the generator, the mechanical input power
to the generator, and the line powers between the generators.
In the natural gas and water CI, sensors collect the outlet
pressure, and inlet flow rate of each pipeline. As shown in Fig.
1, we consider a sensor network that is used to collect data
from the ICI and send it to a central server. Each area in Fig.
1 corresponds to the set of all neighboring. The sensor data
collected from each CI can be expressed as a linear equation
of the states of the ICI, as follows:
y(t) = Cx(t),

(9)

Fig. 1: An illustrative example of an ICI
where x(t) is given in (6), and y is a vector of all the
sensor data at each time instant, and C is a matrix relating
the states to the sensor data. However, due to the inaccuracy
in measurements and the process noise in the infrastructure,
the owner of each CI must estimate the system state at each
time instant. Due to the interdependence between the CIs, their
owners have to share the collected data from the components
with a single administrator who has access to the ICI model
[32]. Note that a lack of cooperation between the owners of
the CIs can yield estimation error since the administrator will
not be able to capture the interdependencies. Considering the
process and measurement noise, and also the discrete sensor
data, we can rewrite the state-space model equations and the
sensor outputs as a discrete linear dynamic system:
x(k + 1) = Ax(k) + Bu(k) + w(k),
(10)
y(k) = Cx(k) + l(k),
where A, B, and C are defined in (6) and (9), x(k) ∈ Rn
is the vector of state variables of the ICI at time step k,

4



Ae +B e (T ge D g T eg C eg +T we D w T ew C ew )

B e T ge C ge

B e T we C we

B g T eg C eg

Ag

0

0

Aw

A,

w

B T

ew

C

ew

Be
,B ,  0
0




B e T ge D g
Bg
0


B e T we D w
 . (8)
0
w
B

ence between the state x(k) and its estimate x̂:
200

Magnitude of state variable

150

xw
231
xg231
xg232
xg234
ue5

Start of change in
the power demand

100

e(k) , x(k) − x̂(k).
(15)
Using (14) and (15), we have:
e(k + 1) = (A − KCA)e(k) + (I − KC)w(k)
(16)
−Kv(k) + KCB (u(k) − (A − KCA)u(k − 1)) ,
We also define the residue of the Kalman filter:
z(k) , y(k) − CAx̂(k).

(17)

50

Because of process and measurement noise, we need to
validate the estimation of the states and detect the failure of
the estimation filter. We use X 2 failure detector allowing the
detector computes the following value at each time step [34]:

0

g(k) = z T (k)Z −1 z(k),

-50
20

20.1

20.2

20.3

20.4 20.5 20.6
Time (second)

20.7

20.8

20.9

21

Fig. 2: Interdependence between the state variables of power,
gas, and water CI.
u(k) ∈ Rp is the vector of external inputs of the ICI at
time step k, w(k) ∈ Rn is the process noise at time k,
and l(k) ∈ Rp is the measurement noise at time k and p
is the number of sensors. Due to the discrete sensor data,
hereinafter, we use (10) in our analysis which is the discrete
model for the ICI. For notational simplicity, we use (A, B, C)
for the discrete state-space model, however, we transform these
matrices to the discrete form using available methods such as
bilinear transformation. Note that x(0) is the initial state of
the ICI, and w(k), l(k) and x(0) are independent Gaussian
random variables with x(0) ∼ N (0, Ψ), w(k) ∼ N (0, Φ),
and l(k) ∼ N (0, Ω).
The ICI administrator seeks to estimate the state of the
ICI using (10). However, due to sensor error and operation
noise, a noise-resilient method is needed to estimate the state
variables. To this end, in [33], the authors showed that by
using a Kalman filter, one can compute the state estimation
x̂(k) from observations y(k). Since the initial time of the ICI
is considered −∞, the Kalman filter converges to a fixed gain
linear estimator. To find the state estimate of the system, we
first compute the Kalman state probability matrix P :

−1
P = AP AT + Φ − AP C T CP C T + Ω CP AT , (11)
then, we compute the Kalman fixed gain as follows:

−1
K = P C T CP C T + Ω
.

(12)

Next, we find the state estimation vector at time k dependent
to knowing the state estimation at time k − 1, x̂(k|k − 1) as:
x̂(k|k − 1) = Ax̂(k − 1|k − 1) + Bu(k − 1).

(13)

Finally, we compute the state estimate x̂(k|k) using a Kalman
filter:
x̂(k|k) = (I − KC)x̂(k|k − 1) + Ky(k),
(14)
where the initial state is defined as x̂(0|0) = x(0).
We can define the estimation error at time k as the differ-

(18)

where Z is the diagonal covariance matrix of the residue
z(k), which implies that the residues are independent of each
other. If g(k) exceeds the threshold level, then the detector
will trigger an alarm.
C. Attack Model
Consider the cyber system of the ICI in Fig. 1, where
sensors collect measurements from the physical components of
the ICI and transmit the measurement data to a central node in
their proximity. Then different central nodes will transmit the
data to a central server that will calculate the estimation of ICI
state variables using the presented Kalman filter in (14). We
refer to the group of sensors which connect to a single central
node, as a sensor cluster (SC). We consider a two-stage attack
model to the cyber system of our ICI. In the first stage, the
attacker aims to compromise the ICI’s state estimation sensors
by collecting their data and sending it to its central processor.
After compromising some of the sensors, in the second stage,
the attacker alters the sensor data to increase the estimation
error of the ICI. In the first stage, the attacker has to break
the security solution that is implemented by the administrator
of the ICI (referred to as the defender, hereinafter).
Our model can be used to capture any ICI security solution
that can include an end-to-end encryption of the sensor data
[35], an implementation of attack detection filter [11], or a
physical protection of the sensors [18]. Therefore, to compromise any SC inside the ICI, the attacker has to collect the
broadcast data from the sensors to the central nodes and compromise the implemented security solution. However, this requires processing of the collected data from the sensors across
the ICI, physical presence of the attacker in the proximity of
the central nodes to collect data, or communication resources
for transmission of the collected data to the attacker’s central
processing unit. Since processing, communication and human
resources are limited, the attacker needs to prioritize between
the sensors based on their importance in the state estimation
of the ICI. From the defender’s point of view, implementing
the aforementioned security solutions, requires computational
resources, communication bandwidth, or financial resources

5

which are restricted in availability for the defender. Therefore,
the defender must also prioritize between the sensors of the
ICI that it seeks to protect.
In summary, the attacker aims to maximize the number of
compromised sensors and the defender seeks to protect the
SCs of the ICI from this cyber attack, under strict resource
limitations at both sides. To analyze this interactions between
the attacker and the defender, first, we analyze the second stage
of attack to find the maximum estimation error caused by the
cyber attack and quantify the importance each SC in the ICI,
then using these values we can formally analyze the attackerdefender interaction and derive optimal defense strategies.

k→∞

k→∞

= lim M Λk−1 M ∆e(1) = 0,

In this section, we analyze the impact of the second stage of
the cyber attack in order to quantify the ability of an attacker
to increase the estimation error by altering the sensor data. We
assume a worst-case scenario for security analysis in which the
attacker has complete knowledge about the system as done in
[11] and was able to compromise some of the SCs in the first
stage. We assume that the attacker can change the data of the
compromised sensors to a desired value in order to disturb the
ICI’s state estimation. Given the set of all compromised SCs,
A, we define attack vector y a , [y a1 T , . . . , y aN T ]T where N
is the number of SCs, and y ai is the 1 × Ni attack vector on
SC i where Ni is the number of sensors in SC i. Also, y i = 0
if i ∈
/ A. Therefore, the linear relationship between the state
variables of the ICI and the sensor data under attack will be:
ȳ(k) = Cx(k) + l(k) + y a (k),
(19)
where ȳ(k) is the vector of sensor measurements under attack,
and y a (k) is independent from w(k), v(k), and x(0). Here,
we assume that the attack to the sensors starts from k = 1.
When the ICI’s cyber system is under attack, the Kalman state
estimation filter of the ICI in (14) changes as follows:
x̄(k|k − 1) = Ax̄(k − 1|k − 1) + Bu(k − 1),
(20)
x̄(k|k) = (I − KC)x̄(k|k − 1) + K ȳ(k),
where x̄(k|k) is the estimate of the states under attack. The
new residue and estimation error are defined as follows:
z̄(k) , ȳ(k) − CAx̄(k − 1),
(21)
ē(k) , x(k) − x̄(k).
We can define the error difference between the ICI under attack
and in absence of attack as follows:
(22)

Using (14) and (20), we can find the following model for the
difference in error and residue:
∆e(k + 1) = (A − KCA)∆e(k) − Ky a (k + 1), (23)
∆z(k + 1) = CA∆e(k) + y a (k + 1).

Proof. Since y a (1) is a vector with y a (1)[i] = 0 for
i ∈
/ A and the attack is an impulse input, then, we have
y a (k) = 0 for k > 1. Therefore, using (23), we have
∆e(k) = Q(k−1) Ky a (1). For a stable system, the impulse
response of the state vector ∆e(k) returns to the origin at
k → ∞ for an arbitrary error deviation in k = 1, ∆e(1), or:
lim ∆e(k) = lim Q(k−1) ∆e(1)

III. M AXIMUM S TATE E STIMATION E RROR IN THE
C OMPROMISED S ENSORS

∆e(k) , ē(k) − e(k), ∆z(k) , z̄(k) − z(k).

Proposition 1. If Q is a Hurwitz stable matrix, then the
maximum CED caused by an impulse attack to a set of sensors
n T
A is:
q m (y a ) , max y a (1)K T EKy a (1),
o
(26)
T
y a (1)K T QT EQKy a (1) .

(24)

We define cumulative error difference (CED) at time step k:
q(k) , ∆eT (k)E∆e(k),
(25)
where E ∈ Rn×n is a positive semidefinite diagonal matrix
which is the relative cost of each state error. For convenience,
we define Q , A − KCA which is a Hurwitz stable matrix
since the ICI model is a stable system [33]. Next, we derive
the maximum CED caused by an impulse attack to an SC.

(27)

k→∞

for any ∆e(1), where the columns of M are the eigenvectors
and Λ is the diagonal matrix of eigenvalues of Q [33]. Assume
that λi where i = 1, . . . , n are the eigenvalues of Q, then, all
the elements of ∆e(k) are a linear combination of components
λk−1
. And since in a Hurwitz stable matrix |λi | < 1 then
i
(∆e(k)[i])2 < (∆e(k − 1)[i])2 for k > 2. Therefore, we have:
n
n
X
X
E[i, i](∆e(k − 1)[i])2
E[i, i](∆e(k)[i])2 <
q(k) =
i=1

i=1

= q(k − 1),

k>2

(28)

where E[i, i] is the element in the i-th column and i-th row of
matrix E. Note that based on the definition of diagonal semi
definite matrix, E[i, i] has a non-negative value. (28) states
that the maximum value of q(k) for k > 1 occurs in q(2).
Thus, the maximum value of q(k) is max{q(1), q(2)}. Since
the initial condition of ∆e(k) is ∆e(1) = −Ky a (1), this
results in:
T
q(1) = y a (1)K T EKy a (1),
(29)
T
q(2) = y a (1)K T QT EQKy a (1).
This proves (26).

Proposition 1 shows that the maximum CED caused by an
impulse attack occurs during first two time instants, q(1), q(2),
after the initiation of the impulse attack to the sensors. In the
presence of the attack, the failure detector in (18) computes
the following value in each time step:
ḡ(k) = z̄ T (k)Z −1 z̄(k)

(30)

Now, we define two new parameters for the analysis of
probability of failure in the system as follows:
β(k) , Pr(g(k) > g t ),
t

β̄(k) , Pr(ḡ(k) > g ),

(31)
(32)

where β(k) and β̄(k) capture probabilities of failure in absence
and existence of attack, respectively, and g t is the failure
trigger threshold.
Definition 1. An impulse attack to set A is α-feasible if:
q
D(z(k)||z̄(k)) = ||∆z(k)||S = ∆z T (k)S∆z(k) ≤ α.
(33)
−1
for all k = 1, . . . , ∞, where S = Z /2 and D(z(k)||z̄(k))
is the Kullback-Leibler (KL) distance between z(k) and z̄(k).

6

Using [14, Theorem 1], we can directly prove the convergence of β̄(k) to β(k) as D(z(k)||z̄(k)) → 0, as follows.
Lemma 1. For any  > 0, there exists α > 0, such that if
D(z(k)||z̄(k)) ≤ α,

(34)

for k = 1, . . . , ∞, then

following quadratic program with quadratic constraints:
o
n T
a
a
aT
a
(40)
qαm (A, A, B, C), max
max
y
R
y
,
y
R
y
1
2
ya
o
n T
T
T
(41)
s.t. max y a P 1 y a , y a P 2 y a , y a P 3 y a ≤ α2 ,
y a [i] = 0,

β̄(k) ≤ β(k) + 

(35)

i∈
/ A,

(42)

where

for all k = 1, . . . , ∞.

R1 = K T EK, R2 = K T QT EQK,

Lemma 1 shows that, if the probability of alarm triggering
β(k) at time k increases by a value of  in presence of attack,
β̄(k) = β(k) + , then, there exists a value for α such that an
impulse attack can be designed with a KL distance lower than
α. Now, if the attacker wants to design an α-feasible impulse
attack then it should change the sensor data such that the KL
distance never exceeds α. Next, we find the maximum KL
distance caused by an impulse attack to a set of sensors.

P 1 = S, P 2 = K T AT C T SCAK,

Lemma 2. If Q is a Hurwitz stable matrix, then the maximum
KL distance caused by an impulse attack to a set of sensors
A will be:
(q
T
m
a
(CAQKy a (1)) SCAQKy a (1)
D (y ) , max
q

y

aT

)
q
T
(1)Sy a (1), (CAKy a (1)) SCAKy a (1)

P 3 = K T QT AT C T SCAQK.
Proof. From Proposition 1, we know that the maximum CED
caused by an impulse attack which we define it as vector y a
in time step k = 1 to a sensor set A is:
o
n T
T
q m (y a )=max y a K T EKy a , y a K T QT EQKy a ,
(43)
where y a [i] = 0 for i ∈
/ A and y a [i] is the i-th entity of vector
y a . From Lemmas 1 and 2, we know that the maximum KL
distance caused by an α-feasible attack to the sensor set A
can not exceed α and therefore we have:
2

(36)

(Dm (S a )) < α2 ,
(
T

T

max y a Sy a , (CAKy a ) SCAKy a ,
Proof. From (23) and (24) we have:

)

a

∆z(1) = y (1),

(37)
(k−2)

∆z(k) = −CAQ

a

Ky (1),

k>1

(38)

therefore, all the elements of vector ∆z(k) will be a linear
combination of eigenvalues of Q powered by k − 2, λk−2
, for
i
k > 2 and thus (∆z(k)[i])2 < (∆z(k − 1)[i])2 for k > 3.
Then, we have:
v
u n
uX
D(z(k)||z̄(k)) = t
S[i, i](∆z(k)[i])2 ,
i=1

v
u n
uX
<t
S[i, i](∆z(k − 1)[i])2 ,

(39)

i=1

= D(z(k − 1)||z̄(k − 1)), k > 3
where S[i, i] is the element in the i-th row and i-th column of S. In addition, since Z is a diagonal matrix with
positive entities, then its inverse has positive values, hence
S = 12 Z −1  0. (39) implies that the KL distance is
decreasing for k > 2, and hence, the maximum KL distance
will occur in k = 1, 2, or 3 and this proves (36).

Lemma 2 finds the maximum KL distance caused by an
impulse attack to a set of sensors. We use the maximum error
caused by an impulse attack and maximum KL distance to
find the maximum CED caused by an α-feasible attack in
the following theorem. This theorem essentially quantifies the
maximum CED that the attacker can cause without triggering
the alarm to a set of sensors.
Theorem 1. The maximum CED caused by an impulse αfeasible attack to a set of sensors A is the solution of the

a T

(CAQKy ) SCAQKy a

< α2 ,

(44)

then, y a should maximize (43) with constraints in (44), and
considering y a [i] = 0 for i ∈
/ A. Also, since E and S are
positive definite matrices then R1 , R2 , P 1 , P 2 , and P 3 are
all semi positive definite matrices. Because of positive semi
definitiveness, (40) is a convex function and (41) is a convex
bounded constraint and the solution of the maximization
problem will lie at the boundaries [36].

Theorem 1 provides a method for the attacker to find the
maximum CED caused by altering a set of sensors without
triggering failure alarm. To solve the optimization problem in
Theorem 1, known techniques such as quadratic programming
can be used [36]. Using theorem 1, we can assign a value to
quantify the maximum CED for each of the ICI’s SCs. To do
so, for each of the ICI’s SCs we calculate the following value:
ϕi (A, B, C) , qαm (Ni , A, B, C),

(45)

where ϕi is the value of SC i in the state estimation and
Ni is the set of sensors inside SC i. This value captures the
importance of each SC for the attacker and the defender in
the first stage of attack, because the attacker can increase
the estimation error by ϕi in the second stage of attack after
compromising the SC i in the first stage. Based on this value
both the attacker and the defender can prioritize between their
actions in the first stage. Since we can now quantify the
importance of different SCs under attack, next, we study how
the ICI can defend against the first stage of attack during which
the attacker and the defender should allocate their available
resources on all the SCs based on their values.

7

IV. ICI S ECURITY R ESOURCE A LLOCATION AS A
C OLONEL B LOTTO G AME
In this section, we analyze the resource allocation of the
attacker and the defender in the first stage of the cyber
attack. In the considered model, the available resources for the
defender and attacker are denoted by Rd , and Ra , respectively.
Consequently, the defender and the attacker must simultaneously allocate their resources across a finite number of SCs,
N . Moreover, each SC in ICI i has a value, ϕi (A, B, C)
given by (45) which quantifies the maximum CED caused by
compromising SC i. This value captures both the cyber and
physical nature of the ICI as per (45). Hereinafter, we use
subscripts a and d to denote the attacker and the defender,
respectively.
j T
Also, r j = [r1j , . . . , rN
] denotes player j’s allocation
vector across N SCs. In each SC i, the defender assigns
a protection level which requires rid resources. In contrast,
the attacker spends some effort to break the sensor’s security
mechanism, which requires ria resources in SC i. For instance,
in end-to-end encryption, for any encryption algorithm the
defender must consider a number of computations in the
decryption of each SC’s messages in the central server [35]. To
break this encryption, the attacker must collect the messages
of each SC and break the encryption using a large number of
computation, which requires the attacker to assign a portion
of its available computational resources for each SC. Such a
resource limitation is not restricted to cases of encryption as it
can also be applied to other protection methods such as attack
detection filters [11].
Therefore, for any protection method, in each SC, if the
defender allocates more resources than the attacker then the
defender prevents that SC from being compromised. In this
case, we assign the normalized value of SC i to the defender
and zero to attacker if the defender wins SC i. In contrast,
if the attacker allocates a higher number of resources in each
SC, then the attacker can compromise that SC. In this case,
we assign the normalized value of SC i to the attacker and
zero to the defender if the attacker wins SC i (i.e., in this
case, the CED is zero, and the defender perfectly protects its
SC). Also, in case of equal allocation of resources, which has
the probability of zero due to the continuous action space of
the attacker and the defender, we share the normalized value
of each SC equally between the attacker and the defender.
Therefore, in each SC i, the normalized payoff for the attacker
and defender is given by:

j
−j

φi (A, B, C), if ri > ri ,
j j −j
(46)
vi (ri , ri ) = φi (A,B,C)
,
if rij = ri−j ,
2


0,
if rij < ri−j ,
where −j is the opponent of j and
ϕi (A, B, C)
φi (A, B, C) = PN
.
(47)
m=1 ϕm (A, B, C)
The total payoff of the defender and the attacker resulting
from allocations across all N SCs is the sum of the individual
payoffs in (46) received from each individual SC:
N
X
uj (r j , r −j )=
vij (rij , ri−j ).
(48)
i=1

Here, we define the total CED caused by the allocation vectors
r a and r d as follows:
N
X
π(r a , r d ) , ua (r a , r d )
ϕm (A, B, C),
(49)
m=1

since ua (r a , r d ) captures summation of the estimation errors
from all the SCs. The attacker aims to increase its utility
function in (48) by maximizing the number of compromised
SCs which results in maximizing the total state estimation
error. Also, the defender seeks to increase its utility function
in (48) by maximizing the number of protected SCs from the
cyber attack to minimize the state estimation error. Moreover,
the payoff for each player depends on the actions of both
players and, thus, we can use a game-theoretic approach to
solve this problem [37]. In particular, next, we first model the
problem as a two-player Colonel Blotto game [25] between
the attacker and the defender, and then present the solution for
the game. The Colonel Blotto game framework is particularly
suitable for the considered ICI security problem since, in
this game, two colonels simultaneously allocate their available
military resources on N battlefields, where the winner of
each battlefield is the colonel with a more allocated resources
and both the colonels aim to maximize the number of won
battlefields. This is similar to the problem in (48), in which
SCs are the battlefields and the defender maximizes the
number of protected SCs while the attacker seeks to maximize
the number of compromised SCs.
A. Game Formulation and Pure Strategy Nash Equilibrium
To model the interdependent decision making
processes of the attacker and defender, we introduce
a
noncooperative
Colonel
Blotto
game
[25]
j
P, {Qj }j∈P , {Rj }j∈P , N, {φai , φdi }N
,
{u
}
defined
j∈P
i=1
by six components: a) the players which are the attacker a
and the defender d in the set P , {a, d}, b) the strategy
spaces Qj for j ∈ P, c) available resource Rj for j ∈ P,
d) number of the SCs N , e) normalized value of each SC
i for j ∈ P, φji (A, B, C), and f) the utility function, uj ,
for each player. For both players, the set of pure strategies
Qj corresponds to the different possible resource allocations
(
)
across the SCs:
N
X
j
j
j
j j
Q = r
ri ≤ R , ri ≥ 0 .
(50)
i=1

Also, the utility function of each player, uj , can be defined as
in (48). The utility function in (48) is a symmetric case for the
Colonel Blotto game where φdi (A, B, C) = φai (A, B, C) =
φi (A, B, C), which indicates that the values of SCs are equal
for the defender and the attacker. In the following subsection,
first we present the solution of the Colonel Blotto game for a
general case of φdi (A, B, C) 6= φai (A, B, C), then we derive
the solution of symmetric case. For notational simplicity,
hereinafter, we drop the arguments (A, B, C) in the notation
of variables ϕi and φji .
One of the most important solution concepts for noncooperative games is that of the Nash equilibrium (NE). The NE
characterizes a state at which no player j can improve its
utility by changing its own strategy, given the strategy of the

8

other player is fixed. For a noncooperative game, the NE in
pure (deterministic) strategies can be defined as follows:
Definition 2. A pure-strategy Nash equilibrium of a noncoop∗
erative game is a vector of strategies [r a ∗ , r d ] ∈ Qa × Qd
such that ∀j ∈ P, the following holds true:
∗

∗

∗

uj (r j , r −j ) ≥ uj (r j , r −j ), ∀r j ∈ Qj .

(51)

The NE characterize a stable game state at which the
defender cannot improve the protection of the ICI’s SCs by
unilaterally changing its action r d given that the action of
the attacker is fixed at r a ∗ . At the NE, the attacker cannot
increase the state estimation error of the ICI by changing
its action, r a , when the defender keeps its action fixed at
r d∗ . However, the NE is not guaranteed to exist in pure
strategies. In particular, for a Colonel Blotto game, without
loss of generality, if Rd > Ra , then it can be proven that, for
N Ra > Rd there exist no pure-strategy NE [25]. However, it
is proven that there exists at least one NE in mixed strategies
[38] for noncooperative games. When using mixed strategies,
each player will assign a probability for playing each one of
its pure strategies. For an ICI security problem, the use of
mixed strategies is motivated by two facts: a) both players must
randomize over their strategies in order to make it nontrivial
for the opponent to guess their potential action, and b) the
allocation of resources can be repeated over an infinite time
duration and mixed strategies can capture the frequency of
choosing certain strategies for both players. A mixed strategy,
which can be termed as a distribution of resources, for player
j is an N -variate distribution function Gj : RN
+ → [0, 1]
with support contained in player j’s set of feasible allocations,
Qj . We also define univariate marginal distribution functions
(MDFs) {Fij }N
i=1 : R+ → [0, 1] for each SC i and can be
called as distribution of resources on each SC i.

Therefore, player j’s optimization problem considering its
constraint on the available resource is:

N Z ∞ h
i
X
j −j j
j
j j
φi Fi (ri ) − ζ ri dFi + ζ j rj , (54)
max
j
{Fi }N
i=1 i=1
j

0

where ζ is a multiplier for player j’s expected resource allocation constraint. For each i = 1, . . . , N , the corresponding
first-order condition for maximizing (54) is given by:
i
d h j −j j
φi Fi (ri ) − ζ j = 0,
j
dri
φji d −j j
F (ri ) = 1,
(55)
ζ j drij i
where (55) is equivalent to the necessary condition for a single
all-pay auction game where player j’s value for the prize in
φj
φ−j
φj
i
auction is ζ ji [40]. In such an all-pay auction, if ζ ji ≥ ζ −j
the solution of (55) is described as follows:

 j
"
#
φi
φ−j
i
−j
−
j
−j
r
φ
ζ
ζ
−j
i
 + j , r ∈ 0,
, (56)
Fi (r) = 
φji
φi
ζ −j
j
j
ζ
ζ
"
#
r
φ−j
j
i
Fi (r) = −j ,
r ∈ 0, −j . (57)
φi
ζ
ζ −j

Now, to find the multipliers (ζ a , ζ d ), let µ ,
φa
i
φd
i

In a game-theoretic setting, each player chooses its own
mixed-strategy distribution to maximize its expected utility.
We first derive the solution for a special case of our problem
in which the attacker and the defender consider the expected
allocation of their resources on each SC instead of exact
allocation. This is a special case of the Colonel Blotto game
known as the General Lotto game [39]. In a Colonel Blotto
game, the sum of allocated resources cannot exceed the limited
resources for the players as in (50). In contrast, in a in General
Lotto game, the sum of expected allocated resource on SCs
cannot exceed the(restricted resource of players:
)
N
X
j
j
j
j j
Q = r
Ei (r) ≤ R , ri ≥ 0 .
(52)
i=1

where Eji (r) is the expected value of resources allocated by
player j on SC i. In this case, the utility of each player j ∈ P
is defined as the expected value over its mixed strategies:
−j N
U j (Gj , G−j ) = U j ({Fij }N
i=1 , {Fi }i=1 )

N Z ∞
X
(53)
=
φji Fi−j (rij )dFij .
i=1

0

and assume

that Ωa (µ) is the set of SCs in which
> µ. Then using
(54), (56), and (57), we have:
 a 2
φi
X φd
X
ζa
i

 = Ra ,
(58)
+
φd
2ζ d
2 i
i∈Ω
/ a (µ)

i∈Ωa (µ)


X
i∈Ωa (µ)

B. Mixed-Strategy Nash Equilibrium Solution

ζa
ζd

2

φd
i
ζd



ζd

2

φa
i
ζa

+

X
i∈Ω
/ a (µ)

φai
= Rd .
2ζ a

(59)

From [39, Propostion 1] we know that there exists at least one
solution to system of equations in (58) and (59).
Now that we characterized the functions that maximize the
expected utility of players in (53), we first define the solution
concept of mixed strategy Nash equilibrium (MSNE) and then,
finalize the solution of Lotto game by deriving its MSNE. The
MSNE is defined as follows:
Definition 3. A mixed strategy profile G∗ constitutes a mixed
strategy Nash equilibrium if for player j we have:
∗

∗

∗

U j (Gj , G−j ) ≥ U j (Gj , G−j ) ∀Gj ∈ G j
(60)
j
where G is the set of all probability distributions for player
j over its action space Qj .
The MSNE for this game characterizes a state of the system
at which the defender has chosen its optimal randomization
over the allocation of resources on SCs and, thus, cannot
improve the protection of ICI’s SCs by changing this choice.
Also, the MSNE for the attacker is a probability distribution
that captures the allocation of its resources over the SCs
in a way to maximize the state estimation error when the
defender chooses its MSNE strategies. Using the definition of
the MSNE, we define the expected CED at MSNE as follows:
∗

∗

d N
Π(Ga ∗ , Gd ) = Π({F a ∗i }N
i=1 , {F i }i=1 )

9

∗

, U a (Ga ∗ , Gd )

N
X

ϕi .

(61)

m=1

It is proven in [39, Theorem 1] that for each solution (ζ d , ζ a )
for system of equations in (58) and (59), each player in
a General Lotto game has a unique MSNE with univariate
marginal distributions in (56) and (57). In the following
proposition, we characterize the solution for our problem when
the values of the ICI’s SCs for both attacker and defender are
equal and, then, we find the expected state estimation error.
Proposition 2. For the problem of resource allocation over
SCs having equal values for the attacker and the defender
, at the MSNE, the MDFs for the
φai = φdi , PNϕi (v)
i=1 ϕi (v)
attacker and the defender, when the defender’s resources are
greater than the attacker’s resources, Rd ≥ Ra , will be given
by:


r Ra
Ra
a∗
, r ∈ [0, 2φi Rd ], (62)
Fi (r)= 1 − d +
R
2φi Rd Rd
r
∗
Fid (r)=
,
r ∈ [0, 2φi Rd ], (63)
2φi Rd
PN
Ra
and the expected CED at MSNE will be 2R
d
i=1 ϕi .
Proof. Considering φai = φdi and Rd ≥ Ra the solution for
the system of equations in (59) and (58) is ζ a = 2R1 d , ζ d =
Ra
. Then, substituting ζ a and ζ d into (56) and (57) results
2(Rd )2
in marginal distributions in (62) and (63). To find the expected
utility of the attacker, we have:



N
N
U a {Fia ∗ }i=1 , Fid∗ i=1 =
"
# N
N Z 2φi Rd
X
X Ra
Ra
Ra
r
φi
dr
=
φ
=
, (64)
i
2φi Rd 2φi (Rd )2
2Rd
2Rd
0
i=1
i=1
Also, the expected utility of the defender can be derived as:


Ra
N
N
,
(65)
U d Fid∗ i=1 , {Fia ∗ }i=1 = 1 −
2Rd
and from (61) the expected
PN CED at MSNE is:
Ra
d∗ N
Π({F a ∗i }N
,
{F
}
)
=

i i=1
i=1
i=1 ϕi .
2Rd
PN
Ra
The value 2R
d
i=1 ϕi captures the expected estimation
error of ICI for the case in which, in the first stage of
the attack, the available resource for the attacker and the
defender are Ra and Rd , respectively. From Proposition 2,
we can conclude two important points: a) the probability of
allocation
of resources
greater than 2φi Rd to SC i is zero,


j
d
Pr ri > 2φi R = 0, b) as much the ratio of the attacker’s
a

,
available resource to the defender’s available resource, R
Rd
increases, the expected state estimation error increases. Next,
we prove that, if the defender concentrates only on one of
the CIs without considering their interdependence, then the
expected estimation error increases.
Theorem 2. Suppose the sets N e , N g , N w contain the SCs
inside the power, natural gas, and water CIs. Then, the ICI’s
expected estimation error increases if the defender does not
consider the interdependence between the CIs.
Proof. Without loss of generality, suppose that the defender
only protects the water CI while the attacker considers all of
the ICI’s SCs. In this case, we assume that the value of the

SCs in the natural gas and power CI for the defender is ε → 0.
The water CI’s SCs, however, will have new values as follows:
ϕi
φdi = P
, i ∈ N w,
(66)
ϕ
w
m
m∈N
while the values of SCs for the attacker are similar to (47).
Then, the ratio of the values for the attacker and defender is
as follows:
P
(
w ϕm
κ , Pm∈N
< 1, if i ∈ N w ,
φai
N
m=1 ϕm
=
(67)
ϕi
φdi
→ ∞,
if i ∈
/ N w.
ε

To find the values of ζ a and ζ d we consider two cases for µ:
If µ < κ, then Ωa (µ) = N , and, thus, from (58) and (59) we
κRd
κRd
have ζ d = 2R1 a , ζ a = 2(R
a )2 which results in µ = Ra where
it violates the condition µ <Sκ.
If µ ≥ κ, then Ωa (µ) = N g N e , and from (58) and (59) we
Ra
κRd
κ
d
= 2(R
find ζ a = 2R
d, ζ
d )2 which results in µ = Ra ≥ κ.
Therefore, the MDFs for the attacker and the defender are:
(
r =h 0, i ∈
/ Niw
1,
a
(68)
F̄i (r) =
a
a
2φi Rd
rR κ
1− R
+ 2(R
,i ∈ N w
d )2 φ ,r∈ 0,
κ
Rd
i
(
1,
r = h0, i ∈
/ N wi
d
d
F̄i (r) =
(69)
rκ
, r ∈ 0, 2φκi R , i ∈ N w ,
2Rd φi
Using above distribution functions and (61), we can find the
expected CED at the MSNE as follows:
 a
 N
R κ 1−κ X
d N
a N
ϕi . (70)
Π({F̄i }i=1 , {F̄i }i=1 ) =
+
2Rd
2
i=1
To reduce the expected CED in Proposition 2, we must have:
 a
 N
N
R κ 1−κ X
Ra X
ϕ
<
ϕi ,
+
i
2Rd
2
2Rd i=1
i=1
 a

R
(κ − 1)
− 1 < 0,
(71)
Rd
a
and since κ < 1 then we need to have R
> 1 which results in
Rda
a contradiction because we know that R
≤ 1. Therefore, the
Rd
defender is never better off if it does not allocate resources to
the water infrastructure. Hence, when the defender allocates
its resources only on a subset of CIs, although the expected
allocated resources on these CIs increases, the estimation error
of their states will also increase due to the interdependence of
the state variables of the ICI.

Theorem 2 illustrates the role of the interdependence between power, natural gas, and water CIs in the state estimation.
The defender must consider all the CIs and their interdependence in the security analysis. Otherwise, if the defender only
protects one of the CIs, then the attacker can cause higher
estimation error on all the state variables of the ICI.
The solutions presented in Proposition 2 and Theorem 2
are for the General Lotto game in which the constraints on
resources hold true in expectation as in (52). Next, we analyze
a special case for the values of the SCs in which the solution
of the General Lotto game in Proposition 2 can be applied
to our original Colonel Blotto game. Using the results in [41,
Corollary 1] and [41, Corollary 2] we can show directly:
Lemma 3. Consider ψ = [ψ1 , . . . , ψN ] is a random permutation vector of SC values φ1 , . . . , φN , then if there exists a

10

Ra /Rd =0.2

if N is odd, (72)

k=0
N
X

(−1)k ψk = 0,

if N is even, (73)

k=1

Then, we can apply the solution of the General Lotto game in
Proposition 2 to the Colonel Blotto game
Finding ψ that satisfies the conditions in Proposition 3 is
computationally hard for large values of N . Thus, we now
prove that, for special values of SCs when the number of SCs
is sufficiently large, we can apply the solution of the General
Lotto game in Proposition 2 to the Colonel Blotto game.
a

R
Theorem 3. For N sufficiently large, if φi < 2R
d and the least
integer
common
multiple
of
all
the
SC
values
is
smaller than
F
N , (φ1 , . . . , φN ) < 2N , then the solution of the General
Lotto game in Proposition 2 constitutes an MSNE for the
Colonel Blotto game.

Proof. From Proposition 2, we have rij < 2φi Rd , and, hence,
to satisfy (50), we must have ria < Ra and rid < Rd which
Ra
results in φi < 2R
d . Now, suppose that the least integer
F
common multiple of all the SC values (φ1 , . . . , φN ) = Φ,
then from Lemma 3 we know that to use the solution of
General Lotto game in Proposition 2 for the Colonel Blotto
game, we have to have a ψ such that:
N
−1
X
0
(−1)2+k ψi+k
> 0, ∀i = 1, . . . , N, if N is odd, (74)
k=0
N
X

Percentage of compromised SCs in simulations
Average percentage of compromised SCs in simulations
Theoretical expected percentage of compromised SCs

30
Percentage of compromised SCs

permutation ψ such that we have:
N
−1
X
(−1)2+k ψi+k > 0, ∀i = 1, . . . , N,

25
20
15
10
5
0
5

10

15

20
25
30
Index of simulation

35

40

45

50

Fig. 3: Percentage of compromised SCs in 50 simulations.
as maxi {ιi } < 2N , then the solution of the General Lotto
game in Proposition 2 constitutes an MSNE for the Colonel
Blotto game under the resource constraints in (52).
In this section, we derived two key theorems. Theorem 2
highlights the importance of interdependence in the security
analysis and proves that the defender cannot minimize the
state estimation error by allocating resources to a single CI.
Moreover, we characterized the solution of our Colonel Blotto
game for the ICI model in Theorem 3 and we derived the
expected estimation error at the MSNE as a function of the
attacker’s and the defender’s available resources.
V. S IMULATION R ESULTS AND A NALYSIS

(−1)

k

ψk0

= 0,

if N is even, (75)

k=1
ψi0 =

where
ψi Φ. It is proven in [42] that, for i = 1, . . . , N , if
ψi0 ∈ {1, . . . , 2m }, then for lim m
N = 1 there exists a vecm,N →∞

N
tor ψ satisfying (74) and (75). Therefore, if we have
 Φ <N2
0
N
m
0
then ψi < Φ < 2 ' 2 and, hence, ψi ∈ 1, . . . , 2 .
This means that there exists a Blotto game solution using
the MDFs in Proposition 2 and the N -copula which maps
marginal univariate distributions in Proposition 2 to an N variate distribution can be expressed as a uniform distribution
[25, Theorem 4]:
 a NY
N
N
R
a a
a
a
G (r1 , . . . , rN )=
φ
,
r
∈
[0, 2φi Rd ]. (76)
i
2(Rd )2 i=1
i=1

N Y
N
N
1
d
d
φ
,
r
∈
[0, 2φi Rd ]. (77)
Gd (r1d , . . . , rN
)=
i
2Rd
i=1
i=1

×

×


Theorem 3 provides a condition for the SCs to apply the
solution of Proposition 2 to our Colonel Blotto game, under
the more general constraint on resources in (52). The following
corollary proves that, in an ICI with a large number of SCs,
the normalized values of SCs, φi , can be approximated such
that the condition in Theorem 3 can be satisfied.
Corollary 1. If there exists 0 < φ < 1 such that φi ' ιi φ
where ιi is an integer number for i = 1, . . . , N , then as long

For our simulations, we study the ICI example in Fig. 1
which captures a real-world ICI scenario such as in [16], [27],
[30] and [31]. In this ICI, each power generator and each
junction in the natural gas and water systems has a demand
profile, which specifies power, gas, and, water demand, at
each time. We design A and B such the dynamic ICI model
constitutes an asymptotically stable system. Also, we consider
32 SCs which collect sensor data from different physical
components. Based on the sensor network architecture in Fig.
1, we generate the matrix C. In our simulations, we consider
0.5-feasible attacks that we use to compute the values of the
SCs, φi (A, B, C). Moreover, due to the large number of the
SCs and small values of the SCs, we can approximate the
values of the SCs as φi ' ιi 10−4 and maxi {ιi } = 661 < 232
which satisfies the condition in Corollary 1 and, therefore, the
solution of Blotto game in Theorem 3 will hold.
In Fig. 3, we simulate the presented ICI in Fig. 1 for the
a
= 0.2.
case in which the ratio of the available resources is R
Rd
Over 50 simulation runs, both players, empirically, play their
MSNE based on the MDFs in Proposition 2. The average
percentage of compromised SCs in this case is 11% which
is very close to the theoretical expected percentage of the
Ra
compromised SCs which is 2R
d = 10%. Fig. 3 shows that
since the attacker and the defender randomize between their
strategies in each simulation, the percentage of compromised
SCs of each simulation might differ from the theoretical expected percentage of compromised SCs. However, the average

11

Theoretical expected % of compromised SAs both playing MSNE
% of Compromised SCs both playing MSNE

100

% of Compromised SCs, the att. best response, the def. φ portion of Rd

Percentage of compromised SCs

90

i

% of Compromised SCs, the def. best response, the att. φ i portion of R

80

a

70
60
50
40
30
20
10
0
0.1

0.2

0.3

0.4
0.5
0.6
0.7
Ratio of available resources (Ra /Rd )

0.8

0.9

1

Fig. 4: Average percentage of compromised SCs.

Fig. 5: The average estimation error caused by the attack to
a
1 10
= 20
the ICI at MSNE for two cases: R
, 20 .
Rd
percentage of compromised SCs for 50 simulations conforms
to the theoretical prediction.
To analyze the efficiency of the MSNE , we consider three
baseline approaches: a) both players play their MSNE (which
is the solution of the game), b) the attacker allocates a portion
φi of its resource on SC i while the defender plays best
response, and c) the defender allocates a portion φi of its
resource on SC i while the attacker plays a best response.
The best response of player j is the pure strategy which
maximizes its utility function. Fig. 4 shows that, when the
defender allocates a portion φi of its resource on SC i without
mixing the allocation of resources, the attacker can predict the
defender’s strategy and compromise more SCs. In contrast,
if the attacker allocates a portion φi of its resources on SC
i without playing a mixed strategy, then the defender can

protect almost all the SCs since it knows the exact amount
of attacker’s allocated resource on each SC. Thus, the attacker
is clearly better off randomizing using an MSNE. Therefore,
Fig. 4 shows that if the defender plays MSNE, it can protect
the SCs at least 50% better than the case in which it allocates
a portion φi of its resource to every SC i. Moreover, Fig. 4
shows that, as much the portion of the attacker’s resources
to the defender’s resources increases, the percentage of the
compromised SCs by the attacker increases.
In Fig. 5, we simulate the use of Kalman filtering for the
estimation of one of the ICI’s state variables, ω6 , in presence
a
.
of the attack and in absence of the attack, for two values of R
Rd
Fig. 5 shows that the estimation of the state
variable
ω
causes
6
a
= 10
10 times higher error in the case of R
20 compared to
Rd 1
a
R
1
the case of Rd 2 = 20 when both players play their MSNE.
To explain this result, from Proposition 2, we knowPthat the
N
Ra
expected estimation error of the attack is Πa = 2R
d
i=1 ϕi .
Hence,a the error caused by the attacker increases by a factor
(R
)
Rd 1
of R
= 10.
( Rad )2
To analyze the interdependence between the three CIs, we
simulate our model for cases in which the defender protects
only one of the CIs. Fig. 6a shows the average percentage of
compromised SCs. From Fig. 6a, we can see that the defender
loses more SCs to the attacker by concentrating on the security
of only one CI. Moreover, Fig. 6b shows the average CED
caused by the cyber attack on the ICI. Fig. 6b shows that the
estimation of state variables of each CI depends on the data
of the SCs in other CIs. Therefore, when the defender focuses
only on one CI to protect, the attacker can disturb this CI’s
state estimation by attacking to the SCs in other CIs. The
simulation results in Fig. 6 corroborate the theoretical results
in Theorem 2, where we proved that the defender cannot
protect the estimation error by allocating its resources only
on one CI without considering the interdependence between
the CIs.
In Fig. 7, we analyze the estimation error of one of the
state variables of the pipeline between junction 2 and 3 of
natural gas CI. In this simulation, we examine two games
where in first case, we consider that the defender can allocate
resources to all the SCs of ICI, while in the second case,
we consider that the defender can only protect the SCs of
the natural gas CI. Note that the ratio between the players’
a
= 0.25 in two cases. Fig. 7 shows that, at
resources is R
Rd
MSNE although the defender protects only the SCs of the
natural gas CI, the average estimation error of a natural gas
CI state variable in this case is larger than the estimation error
in the case of protecting all the SCs of ICI. Moreover, although
the allocated resources on natural gas CI is increased but from
Fig. 7 we see that the average estimation error in natural
gas state variable is increased with the factor of 0.47
0.18 = 2.6.
This difference between the error estimation of two cases is
very close to the theoretical expected utility of the attacker in
Theorem 2. From Theorem 2, we know that the ratio
of the
a
( RRdκ +1−κ)
expected utility of the attacker in two cases is
.
Ra
P

g

ϕ

Rd

i
= 0.38. Therefore, the ratio of the
Here, κ = Pi∈N
N
i=1 ϕi
expected utility of the attacker in two cases is 2.86. Fig. 7

12

0.22
80

Def. allocates on all SCs
Def. allocates only on power SCs
Def. allocates only on gas SCs
Def. allocates only on water SCs

0.18
0.16

60

0.14
50

CED of ICI

Percentage of compromised SCs

70

Def. allocates on all SCs
Def. allocates only on power SCs
Def. allocates only on gas SCs
Def. allocates only on water SCs

0.2

40

0.12
0.1
0.08

30

0.06
20
0.04
10

0.02

0
0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

Ratio of available resources (Ra /Rd )

(a) The average percentage of compromised SCs in ICI when the
defender allocates resources to only one SC.

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

Ratio of available resources (Ra /Rd )

(b) The average CED caused by the cyber attack to ICI when the
defender allocates resources to only one SC.

Fig. 6: Interdependence analysis of the ICI.
general model of protection of SCs of an ICI when the
available resources of the defender and the attacker are limited.
In particular, we have formulated the problem of allocating the
limited resources of the defender and the attacker as a Colonel
Blotto game. We have derived the MSNE of the defender and
the attacker in closed-form as a function of the values of SCs
and the available resources for the attacker and the defender.
The derived MSNE gives insights on the allocation of the
resources on each SC and also underlines the interdependence
of the three infrastructure. The defender’s best strategy when
trying to protect only one of the infrastructure is to allocate
the available resource on all the SCs of the ICI not only on
the infrastructure which it aims to protect. Simulation results
verify that the derived MSNE is the best strategy for the
defender and due to the interdependence of three CIs, the
defender must consider ICI as a unified system in the security
analysis.
R EFERENCES
Fig. 7: Comparison between the average estimation error of
a natural gas state caused by the attack to the ICI when the
defender protects all the ICI VS. only the natural gas CI.
illustrates the interdependence of three CIs. To protect the
states of only one of the CIs, the defender has to consider
the interdependence of CIs, and allocate its resources on all
the SCs of ICI not only the CI which it desires to protect.
VI. CONCLUSION
In this paper, we have analyzed the problem of allocating
limited protection resources on sensor areas (SCs) of an ICI
using a game-theoretic approach. We have modeled the dynamic system of interdependent power, natural gas, and water
infrastructure. We have analyzed the state estimation of the
states of an ICI and the maximum reachability of estimation
errors of different SCs of the ICI. We have considered a

[1] H. Farhangi, “The path of the smart grid,” IEEE power and energy
magazine, vol. 8, no. 1, 2010.
[2] J. Gao, S. V. Buldyrev, H. E. Stanley, and S. Havlin, “Networks formed
from interdependent networks,” Nature physics, vol. 8, no. 1, p. 40,
2012.
[3] S. M. Rinaldi, J. P. Peerenboom, and T. K. Kelly, “Identifying, understanding, and analyzing critical infrastructure interdependencies,” IEEE
Control Systems, vol. 21, no. 6, pp. 11–25, 2001.
[4] T. Nam and T. A. Pardo, “Smart city as urban innovation: Focusing
on management, policy, and context,” in Proceedings of the 5th International Conference on Theory and Practice of Electronic Governance,
ser. ICEGOV ’11. New York, NY, USA: ACM, Sept 2011, pp. 185–194.
[5] J. Shi, J. Wan, H. Yan, and H. Suo, “A survey of cyber-physical systems,”
in 2011 International Conference on Wireless Communications and
Signal Processing (WCSP), Nanjing, China, Nov 2011, pp. 1–6.
[6] S. K. Khaitan and J. D. McCalley, “Design techniques and applications
of cyberphysical systems: A survey,” IEEE Systems Journal, vol. 9, no. 2,
pp. 350–365, 2015.
[7] L. Xie, Y. Mo, and B. Sinopoli, “Integrity data attacks in power market
operations,” IEEE Transactions on Smart Grid, vol. 2, no. 4, pp. 659–
666, Dec 2011.
[8] A. Sanjab and W. Saad, “Data injection attacks on smart grids with
multiple adversaries: A game-theoretic perspective,” IEEE Transactions
on Smart Grid, vol. 7, no. 4, pp. 2038–2049, July 2016.

13

[9] O. Yagan, D. Qian, J. Zhang, and D. Cochran, “Optimal allocation
of interconnecting links in cyber-physical systems: Interdependence,
cascading failures, and robustness,” IEEE Transactions on Parallel and
Distributed Systems, vol. 23, no. 9, pp. 1708–1720, Sept 2012.
[10] A. Ferdowsi, A. Sanjab, W. Saad, and N. B. Mandayam, “Game theory for secure critical interdependent gas-power-water infrastruct-water
infrastructure,” in Proceedings of 10th Resilience Week. Wilmington,
DE, USA: IEEE, Sept 2017.
[11] F. Pasqualetti, F. Dörfler, and F. Bullo, “Attack detection and identification in cyber-physical systems,” IEEE Transactions on Automatic
Control, vol. 58, no. 11, pp. 2715–2729, 2013.
[12] H. Fawzi, P. Tabuada, and S. Diggavi, “Secure estimation and control for
cyber-physical systems under adversarial attacks,” IEEE Transactions on
Automatic Control, vol. 59, no. 6, pp. 1454–1467, 2014.
[13] C. Kwon, W. Liu, and I. Hwang, “Security analysis for cyber-physical
systems against stealthy deception attacks,” in 2013 American Control
Conference, Washington, DC, USA, June 2013, pp. 3344–3349.
[14] Y. Mo, E. Garone, A. Casavola, and B. Sinopoli, “False data injection
attacks against state estimation in wireless sensor networks,” in Proceedings of 49th IEEE Conference on Decision and Control (CDC), Atlanta,
GA, USA, Dec 2010, pp. 5967–5972.
[15] M. Ouyang and Z. Wang, “Resilience assessment of interdependent
infrastructure systems: With a focus on joint restoration modeling and
analysis,” Reliability Engineering & System Safety, vol. 141, pp. 74–82,
2015.
[16] C. Nan and G. Sansavini, “A quantitative method for assessing resilience
of interdependent infrastructures,” Reliability Engineering & System
Safety, vol. 157, pp. 35–53, 2017.
[17] S. E. Chang, T. McDaniels, J. Fox, R. Dhariwal, and H. Longstaff, “Toward disaster-resilient cities: characterizing resilience of infrastructure
systems with expert judgments,” Risk Analysis, vol. 34, no. 3, pp. 416–
434, 2014.
[18] J. Szefer, P. Jamkhedkar, Y.-Y. Chen, and R. B. Lee, “Physical attack protection with human-secure virtualization in data centers,” in IEEE/IFIP
International Conference on Dependable Systems and Networks Workshops (DSN 2012), Boston, MA, USA, June 2012, pp. 1–6.
[19] P. Düssel, C. Gehl, P. Laskov, J.-U. Bußer, C. Störmann, and J. Kästner,
Cyber-Critical Infrastructure Protection Using Real-Time PayloadBased Anomaly Detection. Berlin, Heidelberg: Springer Berlin Heidelberg, 2010, pp. 85–97.
[20] A. Ashok, A. Hahn, and M. Govindarasu, “Cyber-physical security of
wide-area monitoring, protection and control in a smart grid environment,” Journal of advanced research, vol. 5, no. 4, pp. 481–489, 2014.
[21] A. Gupta, G. Schwartz, C. Langbort, S. S. Sastry, and T. Baar, “A threestage colonel blotto game with applications to cyberphysical security,”
in 2014 American Control Conference, Portland, OR, USA, June 2014,
pp. 3820–3825.
[22] A. Ferdowsi, W. Saad, B. Maham, and N. B. Mandayam, “A Colonel
Blotto game for interdependence-aware cyber-physical systems security
in smart cities,” in Proceedings of the 2nd International Workshop on
Science of Smart City Operations and Platforms Engineering. Pittsburgh, PA, USA: ACM, 2017, pp. 7–12.
[23] G. Schwartz, P. Loiseau, and S. S. Sastry, “The heterogeneous Colonel
Blotto game,” in 2014 7th International Conference on NETwork Games,
COntrol and OPtimization (NetGCoop). Trento, Italy: IEEE, 2014, pp.
232–238.
[24] A. Rullo, D. Midi, E. Serra, and E. Bertino, “A game of things:
Strategic allocation of security resources for IoT,” in Proceedings of
the Second International Conference on Internet-of-Things Design and
Implementation, ser. IoTDI ’17. New York, NY, USA: ACM, 2017,
pp. 185–190.
[25] B. Roberson, “The Colonel Blotto game,” Economic Theory, vol. 29,
no. 1, pp. 1–24, Sep 2006.
[26] P. Kundur, N. J. Balu, and M. G. Lauby, Power system stability and
control. McGraw-hill New York, 1994, vol. 7.
[27] G. Andersson, “Dynamics and control of electric power systems,”
Lecture notes, pp. 227–0528, 2012.
[28] W. N. Lubega and A. M. Farid, “Quantitative engineering systems
modeling and analysis of the energy–water nexus,” Applied Energy, vol.
135, pp. 142–157, 2014.
[29] C. Maffezzoni, “Boiler-turbine dynamics in power-plant control,” Control Engineering Practice, vol. 5, no. 3, pp. 301–312, 1997.
[30] R. Alamian, M. Behbahani-Nejad, and A. Ghanbarzadeh, “A state space
model for transient flow simulation in natural gas pipelines,” Journal of
Natural Gas Science and Engineering, vol. 9, pp. 51–59, 2012.

[31] J. Burgschweiger, B. Gnädig, and M. C. Steinbach, “Optimization
models for operative planning in drinking water networks,” Optimization
and Engineering, vol. 10, no. 1, pp. 43–73, Mar 2009.
[32] A. Zanella, N. Bui, A. Castellani, L. Vangelista, and M. Zorzi, “Internet
of things for smart cities,” IEEE Internet of Things Journal, vol. 1, no. 1,
pp. 22–32, 2014.
[33] P. R. Kumar and P. Varaiya, Stochastic systems: Estimation, identification, and adaptive control. SIAM, 2015.
[34] J. Chen and R. J. Patton, Robust model-based fault diagnosis for dynamic
systems. Springer Science & Business Media, 2012, vol. 3.
[35] H. Chan and A. Perrig, “Security and privacy in sensor networks,”
computer, vol. 36, no. 10, pp. 103–105, 2003.
[36] D. P. Baron, “Quadratic programming with quadratic constraints,” Naval
Research Logistics Quarterly, vol. 19, no. 2, pp. 253–260, 6 1972.
[37] G. Bacci, S. Lasaulce, W. Saad, and L. Sanguinetti, “Game theory
for networks: A tutorial on game-theoretic tools for emerging signal
processing applications,” IEEE Signal Processing Magazine, vol. 33,
no. 1, pp. 94–119, Jan 2016.
[38] T. Başar and G. J. Olsder, Dynamic noncooperative game theory. SIAM,
1998.
[39] D. Kovenock and B. Roberson, “Generalizations of the general lotto and
colonel blotto games,” 2015.
[40] M. R. Baye, D. Kovenock, and C. G. De Vries, “The all-pay auction
with complete information,” Economic Theory, vol. 8, no. 2, pp. 291–
305, 1996.
[41] C. Thomas, “N-dimensional blotto game with heterogeneous battlefield
values,” Economic Theory, pp. 1–36, 2017.
[42] C. Borgs, J. T. Chayes, S. Mertens, and B. Pittel, “Phase diagram
for the constrained integer partitioning problem,” Random Structures
& Algorithms, vol. 24, no. 3, pp. 315–380, 2004.
[43] K. A. Pambour, R. Bolado-Lavin, and G. P. Dijkema, “An integrated
transient model for simulating the operation of natural gas transport
systems,” Journal of Natural Gas Science and Engineering, vol. 28, pp.
672–690, 2016.

A PPENDIX A
C RITICAL I NFRASTRUCTURE DYNAMIC M ODEL
In this appendix, we derive a dynamic system model for the
power, natural gas, and water CI. Moreover, we analyze the
dynamic model of interdependence between these three CIs.
A. Power Infrastructure
To analyze the power system, we consider a synchronous
generator connected to a transmission line as the study system.
For a large-scale power system, consisting of ne generators
interconnected through a transmission network, the model
derived in [26] and [27] is used. In this model, each generator
is considered as a subsystem, with the input to each subsystem
i being the power demand from the connected bus, Pei . Any
changes in the power demand, Pei , will result in a change in
the frequency of the generator and the mechanical input to the
generator. The block diagram of two connected generators is
shown in Fig. 8. In this model, the dynamics of each subsystem
i can be written as follows:
Mechanical dynamics:


ne
X
Di
1
Pmi (t) − Pei (t) −
Pij .
ω̇i (t) = − ωi (t) +
Ji
Ji
j=1,j6=i

(78)
Mechanical feedback dynamics:
1
1
Ṗmi (t) = − Pmi (t) − o ωi ,
Tti
Pi Tti
Line dynamics:
Ṗij = Pijo (ωi − ωj ),

(79)

(80)

14


−1
l

 Ji , x [m] = Plij for j = 1, . . . , ne , j 6= i
e
H i [m] = J1i , xl [m] = Plji for j = 1, . . . , ne , j 6= i


0,
otherwise,
(85)
where xl [m] is the m-th entity of vector xl .

Fig. 8: Block diagram of two connected generators
Parameter
δio
ωi
Ji , D i
Pmi
P ei
Tti
xlij
Pij
o
Pij
Ei
eg
ew
neg
new
ηig
nl

Description
the constant operational power angle of the i-th generator,
in rad
relative speed of the i-th generator, in rad/s
rotor inertia and damping factor of generator i
mechanical input power of the i-th generator
active power demand in the subsystem i
direct-axis transient time constant of mechanical power of
i-th generator
equivalent reactance of the bus between the generators i and
j.
Power flow in the line between generators i and j
linearized power flow in the constant voltage
constant voltage of generator i
vector of indexes of all the generators supplied by natural
gas
vector of indexes of all the generators which require water
number of generators supplied by gas
number of generators require water
the efficiency of transformation of gas fuel to mechanical
power in the generator i
number of lines connected between the subsystems

Table I: Electrical system parameter description
where Pijo is the linearized power flow in the constant voltage
as follows:
Ei Ej
Pijo =
cos(δio − δjo ).
(81)
xlij
All the parameters in the above equations are defined in Table
I.
Using (78), (79), and (80) we can summarize the state space
model for each power subsystem as follows:
ẋei (t)

=

Aei xei (t)

+

B ei uei (t)



xei1
xei2



ωi
Pmi
#

+

H ei xl ,

i = 1, . . . , ne
(82)

where
xei

=
"

Aei

=



i
−D
Ji

−1
Pio Tti

=
1
Ji
−1
Tti



, uei =

, B ei



uei1

−1
Ji





=



Pe i



(83)

In (82), the input to a subsystem, uei , is the demand power
from each generator, Pei . Any change in power demand results
in the deviation of state variables. Since one of the state variables in (82) is the mechanical power input to the generator,
Pmi , then any change in the demand power will result in the
deviation of mechanical power input to the generator. However
mechanical input to the generator is excited by an external
input such as nuclear energy, coal, natural gas, wind, or water
flow. Here, we focus on dynamics of the generators excited by
the natural gas and for the other types of generators we assume
a known mechanical input without any dynamics. We also
assume that the mechanical input to the generators supplied
by natural gas is proportional to the flow rate of arriving gas
to the generator from the natural gas network. Therefore, we
g
out
have Pmi = ηig ϑout
i , where ηi and ϑi are defined in Table II.
Power generation is also dependent on some form of water
input. For instance, thermal power generation requires large
volumes of water for cooling purposes. Hydroelectric power
requires flowing water to drive generating turbines [28]. We
model the dependence of the power generation on water
using the boiler turbine dynamics presented in [29]. Since
control of the generator’s temperature requires water flow
to the generator, the relationship between the input water
and the generator’s mechanical power can be written as:
Pmi = ηiw ξiout , where ηiW and ξiout are defined in Table III.
(82) considers each generator separately. However, to analyze the power infrastructure as a unified system we need
to define a centralized model for the power system. In this
regard, we present the state space model in (1), where
 e 


x1
u1
 .. 




xe =  .  , ue =  ...  ,
(86)
 xen 
e
u
ne
xl

 e
A1
0
0
H e1
.. 

..

.
0
. 
Ae(2ne +nl )×(2ne +nl ) =  0
 , (87)
 0
0 Aene H eng 
Ael

0


P o , if 2j = 1, . . . , ne , & xl [i] = P j m ,

 2j m
2
Ael
−P jo m , if 2j = 1, . . . , ne , & xl [i] = Pm j ,
nl ×2ne [i, j] =
2

2

0,
otherwise,
(88)

=



.

(84)

and xl is a vector with size nl consisting all the line powers
Pij . Also, H ei is a row vector with entities as follows:


e

B e1


B = 0
0

0
..
.
0

0




0 ,
B ene

(89)

15

Parameter
in
gij
out
gij
comp

gi
gidef
ϑdi
ϑin
ij

ϑout
ij

Fig. 9: Block diagram of
gas pipeline model
( natural
1
if j/2 = eg [i],
g
ηj/2
eg
C neg ×(2ne +nl ) [i, j] =
(90)
0,
otherwise,
( 1
if j/2 = ew [i],
w
ηj/2
(91)
C ew
new ×(2ne +nl ) [i, j] =
0,
otherwise,
and y eg and y ew are the vectors of natural gas and water
demand from power system, respectively. Also eg and ew are
the vectors whose entities are the indexes of generators require
natural gas and water.

γij
nout
i
pipe
ng

Description
the gas pressure at the inlet of the pipeline connecting
junction i to j
the gas pressure at the outlet of the pipeline connecting
junction i to j
the pressure produced by the compressor in the junction i
the predefined pressure in the junction i
the gas flow demand in junction i
the gas flow rate at the inlet of the pipeline connecting
junction i to j
the gas flow rate at the outlet of the pipeline connecting
junction i to j
The sector area of a pipeline between the junction i and j
Number of pipelines suppling junction i
Total number of the pipelines in the gas system

Table II: Natural gas system parameter description
of inlet gas flow rate in the junctions has to be equal to the
summation of the outlet gas flow rate. Therefore, the inletoutlet gas rate relationship in each junction i is as follows:
X
X
d
ϑout
ϑin
mi =
ij + ϑi ,
{m|m→i}

{j|i→j}


ϑout
mi

=P

B. Natural Gas Infrastructure
To analyze the performance of a natural gas system, we
need to model the transient flow of the gas pipelines, as done
in [30]. In this model, the inlet pressure and outlet flow rate
at each pipeline are the inputs to the pipeline system. Also,
the outlet gas pressure and the inlet flow rate are the outputs
of the pipeline. The reason that outlet flow rate is the input
of the pipeline system is that the outlet flow rate determines
the demand gas flow rate in the output of the system. Also,
the pressure in the inlet of the pipelines is controlled by the
components such as the compressor. Therefore, we consider
the inlet pressure of the pipeline as another input to the
pipeline system. The block diagram of a pipeline is shown
in Figure 9 and the state space equations for a pipeline are
derived as follows:
−1 g τ1 in
ẋg1 =
x + g ,
ρ1 1 ρ1
−1 g
1
x2 + g in ,
ẋg2 =
ρ̂1
ρ̂1
(92)
−1 g τ2 out
g
ẋ3 =
x3 − ϑ ,
ρ̂1
ρ̂1
−1 g
1
g
ẋ4 =
x + ϑout ,
τ3 4 τ3
and the outputs of the pipeline are:


τ̂1
τ̂1 τ2 out
g
out
g = x1 + 1 −
xg3 −
ϑ ,
ρ̂1
ρ̂1
(93)
−τ̂2 g
τ̂2
ϑin =
x2 + xg4 − g in .
ρ̂1
ρ̂1
The parameters τ1 , τ2 , τ3 , τ̂1 , ρ1 , and ρ̂1 are computed from
the linearization of the pressure and gas flow rate inside the
pipeline. For detailed analysis, the reader can refer to [30,
Annex A]. Other variables in the natural gas network are
described in Table II. In the network of gas distribution
system, we need to take into account the connection between
the pipelines. To do so, first we consider that the summation



γmi1
{q|q→i}

X
γqi1

ϑin
ij



+

(94)

ϑdi  ,

{j|i→j}

where the notation j|j → i means all the junctions supplying
junction i and j|i → j means all the junctions supplied by
junction i. Also, the inlet-outlet pressure relationship in each
junction can be written as follows:
1 X out
giout , out
gji ,
ni
j|j→i
(95)
giout + gicomp = gidef ,
in
gij
= gidef , ∀j|i → j.
A compressor in each junction compensates the pressure loss
[43]. The power required at each compressor i is a function of
the pressure produced by the compressor and can be expressed
as follows:

Pei = ηic (gidef − giout ),

(96)

gidef

where
is defined by the designer of the system based on
the characteristics of the junction. Therefore, the gas pressure
in the inlet of the pipelines is fixed.
Now, we can write the state space model of a pipeline
between junction i and j in the gas CI as follows:


n
X
ẋgij = Agij xgij + B gij ugij +
y gjm  , i = 1, . . . , ng
m|j→m


y gij

=

C gij xgij

+

D gij

ugij

n
X

+


y gjm  ,

(97)

m|j→m


ge  ge
y ge
ui − D ge
i = Ci
i

n
X


y gji  ,

j|j→i

where


xgij 1
 g   def 
 xgij  g
uij 1
gi
g
2 

xij =  g  , uij =
=
,
ϑdj
xij 3
ugij 2
g
xij 4


(98)

16

y gij =






Agij = 



g
yij
1
g
yij
2
−1
ρij 1

0
0
0




=
1
−1
ρ̂ij 1

0
0

τi1
ρi1
1
ρ̂i1

out
gij
ϑij in

0
0
−1
ρ̂ij 1

0



, y ge
i =
0
0
0
−1
τij 1



Pe i



,

(99)




,


(100)
Fig. 10: Block diagram of water pipeline model



0
0





2
,
−τji
B gij = 
1
 0

P
ρ̂ij 1

{j|j→i} τji1 
τji1
P
0
τij 3
{j|j→i} τji1




τ̂ij 1
1
0
1
−
0
ρ̂ij 1
,
C gij = 
−τ̂
0
1
0 ρ̂ijij2
1
"
τ̂ij 1 τij 2 #
0
ρ̂ij 1
D gij =
,
τ̂ij 1
0
ρ̂ij 1
i
h
 c 
 def 
1
0 , uge
ηi , D ge
gi
C ge
.
nout
i =
i =
i =
i

(101)

{m|m→i}

{j|i→j}



(102)

(103)
(104)

and
is defined the power demand in the compressor of
each junction i.
The relationship in (97) considers the dynamics of a single
pipeline system. To analyze the natural gas system as a
unified system, we propose the centralized state space model
consisting of all the state variables of all the pipelines in (2),
where

 ge 



ug1
xg1
y1
 ..  g  ..  ge  .. 
g
x =  . ,u =  . ,y =  . .
ug pipe
xg pipe
y ge
ng
ng

(105)
In (105), each vector xgm corresponds to one of the pipeline
state vectors, xgij . Moreover, each vector ugm in (105) corresponds to one of the pipeline input vectors, ugij . To find
matrices Ag , B g , C g , D g , we need to start from the last
pipelines in the natural gas infrastructure and find the outputs
of these pipelines and use them as the inputs of previous
pipelines. By proceeding this method until reaching to the
source junction of natural gas we can derive the entities of
mentioned matrices.

γmi

ξmi = P

γqi



(107)

ξij + ξid  .

{j|i→j}

Moreover, due to the pressure loss in the pipelines, there exist
a water pump in each junction which controls the pressure of
water in the junctions. Then, the water pressure relationship
in each junction is as follows:
1 X out
wji ,
wiout , out
ni
j|j→i
(108)
wiout + wipump = widef ,
in
wij
= widef , ∀j|i → j.
The power needed for the pump in each junction i is a function
of pressure provided by the pump and can be written as
follows:

Pei = ηip (widef − wiout ),

(109)

widef

where
is a predefined value by the designer of the water
system which indicates the fixed inlet water pressure in each
pipe.
Similar to the natural gas system, we can present the state
space model for each pipeline in the water system as follows:


n
X
w  w
w w
 , i = 1, . . . , nw
uij +
yw
ẋw
jm
ij = Aij xij + B ij
m|j→m



n
X

w w
w  w
yw
uij +
ij = C ij xij + D ij


,
yw
jm

(110)

m|j→m


we  we
y we
ui − D we
i = Ci
i

C. Water Infrastructure
To analyze the water infrastructure we use the analytical
framework presented in [31]. In this model the flow rate is
considered constant in the outlet and inlet of each pipeline.
Moreover, the relationship between the flow, inlet and outlet
pressure of the pipeline can be expressed as follows:
1
1
ṙ = −r + ξ − win ,
ν
θ
(106)
wout = −θr + win .
where all the parameters are described in Table III. Figure
10 shows the block diagram model of the pipeline in water
system. Also at each junction of water system the summation
of inlet flow has to be equal to the summation of the outlet


X

{q|q→i}

y ge
i

ng

flow of water. Therefore, at each junction i we have:
X
X
ξmi =
ξij + ξid ,

n
X


,
yw
ji

j|j→i

where


widef
ξjd





rij ,
, yw
(111)
ij =
h
i


−1 , B gij = ν1ij −1
Aw
,
(112)
θij
ij =






−θij
1 0
Pe i ,
Cw
, Dw
, y we
ij =
ij =
i =
0
0 0
(113)
h
i




1
0 , uwe
ηip , D we
widef .
C we
nout
i =
i =
i =
i
(114)
Vector y we
identifies
the
power
demand
from
the
pumps
in
i
each junction.
xw
ij

=



rij



, uw
ij

=

17

Parameter
ηiw
wiout
wiin
ξi
ri
ν
θ

Description
a constant value based on the inner characteristics of the
turbine
the outlet pressure of water at the pipeline i
the inlet pressure of water at the pipeline i
the flow rate of the water in pipeline i
friction of the water at pipeline i
the kinematic viscosity of water
friction coefficient of water

Table III: Water system parameter description
The system presented in (110) is a decentralized model
which considers only one pipeline. Using (107) (108), and
(110), for each pipeline in the water system we can find
the centralized model for the interconnected network of the
pipelines in (3), where
 we 




uw
xw
y1
1
1
 ..  w  ..  we  .. 
w
x =  . ,u =  . ,y =  . ,
uw
xw
y we
nw
npipe
npipe
w
w
(115)
where each xw
corresponds
to
one
of
the
pipeline
state
m
w
corresponds
to
one
of
the
pipeline
vectors, xw
,
and
each
u
m
ij
w
input vectors, uw
ij . Also to find the entities of the matrices A ,
w
w
w
B , C , D , the same procedure in the natural gas system
model can be proceeded. We start from the junctions that only
have the demand outlet flow without any pipelines connected
to their outlet. Then, we introduce the outputs of each pipeline
as inputs to the pipelines of previous step. We continue this
procedure until we reach to the source junctions.
In summary, we proposed the state space modeling for each
of the CIs and analyzed their interdependence. Moreover, we
derived the matrices of dynamic model of each CI.

