# An information theoretic vulnerability metric for data integrity attacks on
smart grids

Xiuzhen Ye1∗ Iñaki Esnaola1,2 Samir M. Perlaza 3,2 and Robert F. Harrison1
<EMAIL_ADDRESS>1Dept. of Automatic Control and Systems Engineering,
University of Sheffield, Sheffield S1 3JD, UK 2Dept. of Electrical
Engineering, Princeton University, Princeton, NJ 08544, USA 3INRIA, Sophia
Antipolis 06902, France

###### Abstract

A novel metric that describes the vulnerability of the measurements in power
systems to data integrity attacks is proposed. The new metric, coined
vulnerability index (VuIx), leverages information theoretic measures to assess
the attack effect on the fundamental limits of the disruption and detection
tradeoff. The result of computing the VuIx of the measurements in the system
yields an ordering of their vulnerability based on the level of exposure to
data integrity attacks. This new framework is used to assess the measurement
vulnerability of IEEE 9-bus and 30-bus test systems and it is observed that
power injection measurements are overwhelmingly more vulnerable to data
integrity attacks than power flow measurements. A detailed numerical
evaluation of the VuIx values for IEEE test systems is provided.

##  1 Introduction

Supervisory Control and Data Acquisition (SCADA) systems and more recently
advanced communication systems facilitate efficient, economic and reliable
operation of power systems GJ_PSanalysis_1994 . For instance, the
communication system transmits the measurements to a state estimator that
evaluates the operational status of the system accurately
AA_PSstateestimation_04 . However, the integration between the physical layer
and the cyber layer exposes the system to cybersecurity threats. Cyber
incidents highlight the vulnerability of power systems to sophisticated
attacks. To ensure the security and reliability of power system operation, it
is essential to quantitatively characterize the vulnerabilities of the system
in order to set up appropriate security mechanisms WZ_CN_13 . To that end,
security metrics provide operationally meaningful vulnerability descriptors
and identify the impact that security threats pose to the system. Moreover,
security metrics enable operators to assess the defence mechanisms
requirements to be embedded into cybersecurity policies, processes, and
technology JA_PearsonEducation_07 . For example, the Common Vulnerability
Scoring System (CVSS) analysis Information Technology (IT) system
MS_IEEESecurityandPrivacy_06 . Typical security metrics for power systems
focus on integrity, availability, and confidentiality as envisioned by the
cybersecurity working group in the NIST Smart Grid interoperability panel
NISTIR_securityguideline_10 . System security objectives are categorized into
system vulnerability, defence power, attack severity, and situations to
develop security metrics in a systematic manner PM_ACM_16 . A cyberphysical
security assessment metric (CP-SAM) based on quantitative factors is proposed
to assess the specific security challenges of microgrid systems VV_TSG_19 .

This fragmented landscape showcases a wide variety of metrics available that
depend on the security services, threat characteristics, and system
parameters. Remarkably, there is a lack of general data integrity
vulnerability metrics for power systems. For instance, the impact of data
injection attacks (DIAs) LY_TISSEC_11  can be assessed with a wide variety of
criteria that depend on the objectives of the attackers CKKPT_SPM_12 ;
MO_JSAC_13 ; IE_TSG_16 . A large body of literature addresses DIAs that
compromise both the confidentiality and integrity of the information contained
by the system measurements NG_Stateestimation_21 . With the unprecedented data
acquisition capabilities available to cyberphysical systems, the attackers can
even learn the statistical structure of the system and incorporate the
underlying stochastic process to launch the attacks SE_TSG_19 ; YE_SGC_20 .
DIAs that operate within a Bayesian framework by leveraging stochastic models
of the system are studied in SE_TSG_19  and YE_SGC_20 . From the perspective
of the operator, the introduction of stochastic descriptors opens the door to
information theoretic quantifications of the measurement vulnerability.

In this paper, we propose a novel information theoretic metric to assess the
vulnerability of measurements in power systems to data integrity attacks.
Specifically, we characterize the fundamental information loss induced by data
integrity attacks via mutual information and the stealthiness of the attack
via Kullback-Leibler divergence. Our aim is to provide a metric that is
grounded on fundamental principles, and therefore, inform the vulnerabilities
of the measurements in the system to a wide range of threats. This is enabled
by the use of information theoretic measures which characterize the amount of
information acquired by the measurements in the system in fundamental terms.

The rest of the paper is organized as follows: In Section 2, we introduce a
Bayesian framework with linearized dynamics for DIAs. Information theoretic
attacks are presented in Section 3. The vulnerability metric on information
theoretic attacks is proposed in Section 4. In Section 5, we characterize the
vulnerability of measurements in uncompromised systems and propose an
algorithm to evaluate the vulnerability of measurements. The vulnerability of
measurements of the IEEE test systems is presented in Section 6. The paper
concludes in Section 7.

The main contributions of this paper follow: (1) A notion of vulnerability for
the measurements in the system is proposed. The proposed notion is
characterized by the information theoretic cost induced by random attacks.
Specifically, mutual information and KL divergence are used to construct a
quantitative measure of vulnerability. (2) The vulnerability assessment of the
measurements is posed as a minimization problem and closed-form expressions
are obtained for the case in which the initial state of the system is
uncompromised. (3) An algorithm that computes the proposed vulnerability
indices for general state estimators in power systems is proposed. (4) The
proposed framework is numerically evaluated in IEEE 9-bus and 30-bus test
systems to obtain qualitative characterizations of the vulnerability of the
measurements in the systems.

Notation: We denote the number of state variables on a given system by $n$ and
the number of the measurements by $m$. The set of positive semidefinite
matrices of size $n\times n$ is denoted by $S_{+}^{n}$. The $n$-dimensional
identity matrix is denoted as $\textbf{I}_{n}$. For a matrix ${\bf
A}\in\mathbb{R}^{m\times n}$, then $({\bf A})_{ij}$ denotes the entry in row
$i$ and column $j$ and $\textnormal{diag}({\bf A})$ denotes the vector formed
by the diagonal entries of ${\bf A}$. The elementary vector ${\bf
e}_{i}\in\mathbb{R}^{n}$ is a vector of zeros with a one in the $i$-th entry.
Random variables are denoted by capital letters and their realizations by the
corresponding lower case, e.g., $x$ is a realization of the random variable
$X$. Vectors of $n$ random variables are denoted by a superscript, e.g.,
$X^{n}=(X_{1},\ldots,X_{n})^{\sf{T}}$ with corresponding realizations denoted
by ${\bf x}$. Given an $n$-dimensional vector
$\hbox{\boldmath$\mu$}\in\mathbb{R}^{n}$ and a matrix
$\hbox{\boldmath$\Sigma$}\in S_{+}^{n}$, we denote by
$\mathcal{N}(\hbox{\boldmath$\mu$},\hbox{\boldmath$\Sigma$})$ the multivariate
Gaussian distribution of dimension $n$ with mean $\mu$ and covariance matrix
$\Sigma$. The mutual information between random variables $X$ and $Y$ is
denoted by $I(X;Y)$ and the Kullback-Leibler (KL) divergence between the
distributions $P$ and $Q$ is denoted by $D(P\|Q)$.

##  2 System model

###  2.1 Observation Model

In a power system the state vector ${\bf x}\in{\mathbb{R}^{n}}$ that contains
the voltages and phase angles at all the buses describes the operational state
of the system. State vector ${\bf x}$ is observed by the acquisition function
$F:{\mathbb{R}^{n}}\rightarrow{\mathbb{R}^{m}}$. A linearized observation
model is considered for state estimation, which yields the observation model
as

$Y^{m}={\bf H}{\bf x}+Z^{m},$ (1)

where ${\bf H}\in{\mathbb{R}^{m\times n}}$ is the Jacobian of the function $F$
at a given operating point and is determined by the parameters and topology of
the system. The vector of measurements $Y^{m}$ is corrupted by additive white
Gaussian noise introduced by the sensors GJ_PSanalysis_1994 ,
AA_PSstateestimation_04 . The noise vector $Z^{m}$ follows a multivariate
Gaussian distribution, that is,

$Z^{m}\sim\mathcal{N}(\textbf{0},\sigma^{2}\textrm{{I}}_{m}),$ (2)

where $\sigma^{2}$ is the noise variance.

In a Bayesian estimation framework, the state variables are described by a
vector of random variables $X^{n}$ with a given distribution. In this study,
we assume $X^{n}$ follows a multivariable Gaussian distribution GE_TSG_18
with zero mean and covariance matrix $\hbox{\boldmath$\Sigma$}_{X\\!X}\in
S_{+}^{n}$, that is,

$X^{n}\sim\mathcal{N}(\textbf{0},\hbox{\boldmath$\Sigma$}_{X\\!X}).$ (3)

From (1), it follows that the vector of measurements is with zero mean and
covariance matrix $\hbox{\boldmath$\Sigma$}_{Y\\!Y}\in S_{+}^{m}$, that is,

$Y^{m}\sim\mathcal{N}(\textbf{0},\hbox{\boldmath$\Sigma$}_{Y\\!Y}),$ (4)

where

$\hbox{\boldmath$\Sigma$}_{Y\\!Y}\stackrel{{\scriptstyle\Delta}}{{=}}{\bf
H}\hbox{\boldmath$\Sigma$}_{X\\!X}{\bf
H}^{\sf{T}}+\sigma^{2}\textrm{{I}}_{m}.$ (5)

###  2.2 Attack Setting

Let us denote the measurements corrupted by the malicious attack vector
$A^{m}\in\mathbb{R}^{m}$ as $Y_{A}^{m}$, that is,

$Y_{A}^{m}={\bf H}X^{n}+Z^{m}+A^{m},$ (6)

where $A^{m}\sim P_{A^{m}}$ is the distribution of the random attack vector
$A^{m}$. With a fixed covariance matrix $\hbox{\boldmath$\Sigma$}_{A\\!A}$,
when the additive disturbance to the system, that is, $Z^{m}+A^{m}$ follows a
multivariate Gaussian distribution, the mutual information between the state
variables $X^{n}$ and the compromised measurements $Y_{A}^{m}$ denoted by
$I(X^{n};Y^{m}_{A})$ is minimized SI_TIT_13 . Hence, from the Lévy-Cramér
decomposition theorem Levy_levycramertheorem_35 ; cramer_levycramertheorem_35
, it holds that for the sum $Z^{m}+A^{m}$ to be Gaussian, given that $Z^{m}$
satisfies (2), then, $A^{m}$ must be Gaussian. Therefore, in the following, we
assume that

$A^{m}\sim\mathcal{N}(\textbf{0},\hbox{\boldmath$\Sigma$}_{A\\!A}),$ (7)

where $\textbf{0}=(0,0,\ldots,0)$ and
$\hbox{\boldmath$\Sigma$}_{A\\!A}\in\mathcal{S}_{+}^{m}$ are the mean vector
and the covariance matrix of the random attack vector $A^{m}$. The assumption
in (7) is further discussed in Section 3. Consequently, the vector of
compromised measurements $Y_{A}^{m}$ follows a multivariate Gaussian
distribution with zero mean and covariance matrix
$\hbox{\boldmath$\Sigma$}_{Y_{A}\\!Y_{A}}\in\mathcal{S}_{+}^{m}$, that is,

$Y^{m}_{A}\sim\mathcal{N}(\textbf{0},\hbox{\boldmath$\Sigma$}_{Y_{A}\\!Y_{A}}),$
(8)

with

$\hbox{\boldmath$\Sigma$}_{Y_{A}\\!Y_{A}}\stackrel{{\scriptstyle\Delta}}{{=}}{\bf
H}\hbox{\boldmath$\Sigma$}_{X\\!X}{\bf
H}^{\sf{T}}+\sigma^{2}\textrm{{I}}_{m}+\hbox{\boldmath$\Sigma$}_{A\\!A}.$ (9)

##  3 Information Theoretic Attacks

The aim of the attack is two fold. Firstly, the attack aims to disrupt the
state estimation procedure. Secondly, it aims to stay undetected. For the
first objective, we minimize the mutual information between the vector of
state variables $X^{n}$ in (3) and the vector of compromised measurements
$Y_{A}^{m}$ in (6), that is, $I(X^{n};Y_{A}^{m})$. In the other words, the
attack yields less information about the state variables contained by the
compromised measurements. The stealthy constraint in the second objective is
captured by the Kullback Leibler (KL) divergence between the distribution
$P_{Y^{m}_{A}}$ in (6) and the distribution $P_{Y^{m}}$ in (1), that is,
$D(P_{Y_{A}^{m}}\|P_{Y^{m}})$. The Chernoff-Stein Lemma book_EIT  states that
the minimization of KL divergence leads to the minimization of the asymptotic
detection probability.

The following propositions characterize mutual information and KL divergence
with Gaussian state variables and attacks, respectively [SE_TSG_19, , Prop. 1,
2].

######  Proposition 1.

The mutual information between the vectors of random variables $X^{n}$ in (3)
and $Y_{A}^{m}$ in (8) is

$\displaystyle
I(X^{n};Y^{m}_{A})=\dfrac{1}{2}\log\dfrac{|\hbox{\boldmath$\Sigma$}_{X\\!X}||\hbox{\boldmath$\Sigma$}_{Y_{A}\\!Y_{A}}|}{|\hbox{\boldmath$\Sigma$}|},$
(10)

where the matrices $\hbox{\boldmath$\Sigma$}_{X\\!X}$ and
$\hbox{\boldmath$\Sigma$}_{Y_{A}\\!Y_{A}}$ are in (3) and (9), respectively;
and the matrix $\Sigma$ is the covariance matrix of the joint distribution of
$X^{n}$ and $Y_{A}^{m}$, that is,
$(X^{n};Y_{A}^{m})\sim\mathcal{N}(\textbf{0},\hbox{\boldmath$\Sigma$})$ with

$\displaystyle\hbox{\boldmath$\Sigma$}=\left(\begin{matrix}\hbox{\boldmath$\Sigma$}_{X\\!X}&\hbox{\boldmath$\Sigma$}_{X\\!X}{\bf
H}^{\sf T}\\\ {\bf H}\hbox{\boldmath$\Sigma$}_{X\\!X}&{\bf
H}\hbox{\boldmath$\Sigma$}_{X\\!X}{\bf H}^{\sf
T}+\sigma^{2}\textbf{I}_{m}+\hbox{\boldmath$\Sigma$}_{A\\!A}\end{matrix}\right),$
(11)

where the real $\sigma\in\mathbb{R}_{+}$ is in (2); the matrices ${\bf H}$ and
$\hbox{\boldmath$\Sigma$}_{A\\!A}$ are in (1) and (7), respectively.

######  Proposition 2.

The KL divergence between the distribution of random vector $Y_{A}^{m}$ in (8)
and the distribution of random vector $Y^{m}$ in (4) is

$\displaystyle
D(P_{Y_{A}^{m}}\|P_{Y^{m}})=\frac{1}{2}\left(\log\frac{\left|\hbox{\boldmath$\Sigma$}_{Y\\!Y}\right|}{\left|\hbox{\boldmath$\Sigma$}_{Y_{A}\\!Y_{A}}\right|}-m+\textnormal{tr}\left(\hbox{\boldmath$\Sigma$}_{Y\\!Y}^{-1}\hbox{\boldmath$\Sigma$}_{Y_{A}\\!Y_{A}}\right)\right),$
(12)

where the matrices $\hbox{\boldmath$\Sigma$}_{Y\\!Y}$ and
$\hbox{\boldmath$\Sigma$}_{A\\!A}$ are in (5) and (7), respectively.

Particularly, the attack construction is proposed in the following
optimization problem SE_TSG_19 ; YE_SGC_20 :

$\min_{P_{A^{m}}}I(X^{n};Y^{m}_{A})+\lambda D(P_{Y_{A}^{m}}\|P_{Y^{m}}),$ (13)

where $\lambda\in\mathbb{R}_{+}$ is the weighting parameter that determines
the tradeoff between mutual information and KL divergence. Note that the
optimization domain in (13) is the set of $m$-dimensional Gaussian
multivariate distributions. The optimal Gaussian attack for $\lambda\geq 1$ as
a solution to (13) is given by SE_TSG_19

$P_{A^{m}}^{*}\sim\mathcal{N}(\mathbf{0},{\lambda^{-1/2}}{\bf
H}\hbox{\boldmath$\Sigma$}_{X\\!X}{\bf H}^{\sf T}).$ (14)

Note that the attack realizations from (14) are nonzero with probability one,
that is, $\mathbb{P}[\left|\textnormal{supp}({A^{m}})\right|=m]=1$, where

$\textnormal{supp}({A^{m}})\stackrel{{\scriptstyle\Delta}}{{=}}\left\\{i:\mathbb{P}\left[A_{i}=0\right]=0\right\\}.$
(15)

The attack implementation requires access to the sensing infrastructure of the
industrial control system (ICS) operating the power system. For that reason,
the attack construction incorporates the sparsity constraint by limiting the
optimization domain over the attack vector $A^{m}$ in (6) to the distributions
with cardinality of the support satisfying $|\textnormal{supp}({A^{m}})|=k\leq
m$, that is,

$\mathcal{P}_{k}\stackrel{{\scriptstyle\Delta}}{{=}}\bigcup_{i=1}^{k}\left\\{{A^{m}}\sim\mathcal{N}(\mathbf{0},\bar{\hbox{\boldmath$\Sigma$}}):\left|\textnormal{supp}({A^{m}})\right|=i\right\\}.$
(16)

The resulting attack construction with the sparsity constraints is

$\min_{\mathcal{P}_{k}}I(X^{n};Y^{m}_{A})+\lambda
D(P_{Y_{A}^{m}}\|P_{Y^{m}}).$ (17)

The following theorem provides the optimal single sensor attack construction.

######  Theorem 1.

[YE_SGC_20, , Th. 1] The solution to the sparse stealth attack construction
problem in (17) for the case $k=1$ is

$\bar{\hbox{\boldmath$\Sigma$}}^{*}=v{\bf e}_{i}{\bf e}_{i}^{\sf T},$ (18)

where

$\displaystyle i$
$\displaystyle=\operatorname*{arg\,min}_{j\in\\{1,2,\ldots,m\\}}\left\\{\left(\hbox{\boldmath$\Sigma$}_{Y\\!Y}^{-1}\right)_{jj}\right\\}$
(19) $\displaystyle v$
$\displaystyle=-\frac{\sigma^{2}}{2}+\frac{1}{2}\left(\sigma^{4}-\frac{4(\underline{$w$}\sigma^{2}-1)}{\lambda\underline{$w$}^{2}}\right)^{\frac{1}{2}}$
(20)

with
$\underline{$w$}\stackrel{{\scriptstyle\Delta}}{{=}}(\hbox{\boldmath$\Sigma$}_{Y\\!Y}^{-1})_{ii}$.

##  4 Vulnerability Metric on Information Theoretic Attacks

###  4.1 Attack Structure with Sequential Measurement Selection

To assess the impact of the attacks to different measurements, we model the
entries of the random attack vector $A^{m}$ with independency, that is,

$P_{A^{m}}=\prod_{i=1}^{m}P_{A_{i}},$ (21)

where $A_{i}$ is the $i$-th entry of $A^{m}$ and for all
$i\in\\{1,2,\ldots,m\\}$, the distribution $P_{A_{i}}$ is Gaussian with zero
mean and variance $v\in\mathbb{R}_{+}$, that is, $A_{i}\sim\mathcal{N}(0,v)$.
Consider that $k$ sensors have been attacked with $k\in\\{0,1,2,\ldots,m-1\\}$
and let the covariance matrix of the corresponding attack vector $A^{m}$ in
(6) be

$\hbox{\boldmath$\Sigma$}\in\mathcal{S}_{k},$ (22)

where $\mathcal{S}_{k}$ is the set of $m$-dimensional positive semidefinite
matrix with $k$ positive entries in the diagonal, that is,

$\mathcal{S}_{k}\stackrel{{\scriptstyle\Delta}}{{=}}\\{{\bf
S}\in\mathcal{S}_{+}^{m}:\|\textnormal{diag}({\bf S})\|_{0}=k\\}.$ (23)

Let the set of measurements that have not been compromised be

$\displaystyle\mathcal{K}_{o}\stackrel{{\scriptstyle\Delta}}{{=}}\\{i\in\\{1,2,\ldots,m\\}:(\hbox{\boldmath$\Sigma$})_{ii}=0\\},$
(24)

where $(\hbox{\boldmath$\Sigma$})_{ii}$ is the entry of $\Sigma$ in row $i$
and column $i$. The sequential measurement selection imposes the following
structure in the covariance matrix of the attack vector in (7):

$\hbox{\boldmath$\Sigma$}_{A\\!A}=\hbox{\boldmath$\Sigma$}+v{\bf e}_{i}{\bf
e}_{i}^{\sf T},$ (25)

where $i\in\mathcal{K}_{o}$ and $v\in\mathbb{R}_{+}$. From equation (25), the
cost function
$f:\mathcal{S}_{k}\times\mathbb{R}_{+}\times\mathbb{R}_{+}\times\mathcal{K}_{o}\rightarrow\mathbb{R}_{+}$
defined by adding (10) and (12) is as follows:

$\displaystyle f(\hbox{\boldmath$\Sigma$},\lambda,v,i)$ (26)
$\displaystyle\stackrel{{\scriptstyle\Delta}}{{=}}$ $\displaystyle
I(X^{n};Y^{m}_{A})+\lambda D(P_{Y_{A}^{m}}\|P_{Y^{m}})$ (27) $\displaystyle=$
$\displaystyle\dfrac{1}{2}\log\dfrac{|\hbox{\boldmath$\Sigma$}_{X\\!X}||\hbox{\boldmath$\Sigma$}_{Y_{A}\\!Y_{A}}|}{|\hbox{\boldmath$\Sigma$}|}$
(28)
$\displaystyle+\frac{1}{2}\lambda\left(\log\frac{\left|\hbox{\boldmath$\Sigma$}_{Y\\!Y}\right|}{\left|\hbox{\boldmath$\Sigma$}_{Y_{A}\\!Y_{A}}\right|}-m+\textnormal{tr}\left(\hbox{\boldmath$\Sigma$}_{Y\\!Y}^{-1}\hbox{\boldmath$\Sigma$}_{Y_{A}\\!Y_{A}}\right)\right)$
$\displaystyle=$
$\displaystyle\dfrac{1}{2}\log\dfrac{|\hbox{\boldmath$\Sigma$}_{Y_{A}\\!Y_{A}}|}{|\sigma^{2}\textbf{I}_{m}+\hbox{\boldmath$\Sigma$}_{A\\!A}|}+\frac{1}{2}\lambda\left(\log\frac{\left|\hbox{\boldmath$\Sigma$}_{Y\\!Y}\right|}{\left|\hbox{\boldmath$\Sigma$}_{Y_{A}\\!Y_{A}}\right|}+\textnormal{tr}\left(\hbox{\boldmath$\Sigma$}_{Y\\!Y}^{-1}\hbox{\boldmath$\Sigma$}_{A\\!A}\right)\right),$
(29) $\displaystyle=$
$\displaystyle\frac{1}{2}(1-\lambda)\textnormal{log}\left|\hbox{\boldmath$\Sigma$}_{Y\\!Y}+\hbox{\boldmath$\Sigma$}+v{\bf
e}_{i}{\bf e}_{i}^{\sf
T}\right|-\frac{1}{2}\textnormal{log}\left|\hbox{\boldmath$\Sigma$}+v{\bf
e}_{i}{\bf e}_{i}^{\sf T}+\sigma^{2}\textbf{I}_{m}\right|$
$\displaystyle+\frac{1}{2}\lambda\left(\textnormal{tr}\left(\hbox{\boldmath$\Sigma$}_{Y\\!Y}^{-1}\left(\hbox{\boldmath$\Sigma$}+v{\bf
e}_{i}{\bf e}_{i}^{\sf
T}\right)\right)+\log\left|\hbox{\boldmath$\Sigma$}_{Y\\!Y}\right|\right),$
(30)

where the inequality in (28) holds from plugging in (10) and (12) into (27);
the equality in (29) follows from cancelling
$|\hbox{\boldmath$\Sigma$}_{X\\!X}|$ in the first term [Seber_book, , Sec.
14.17] and noting that
$\hbox{\boldmath$\Sigma$}_{Y_{A}\\!Y_{A}}=\hbox{\boldmath$\Sigma$}_{Y\\!Y}+\hbox{\boldmath$\Sigma$}_{A\\!A}$
in (9); and the equality in (30) holds from plugging in (25) into (29).

###  4.2 Information theoretic vulnerability of a measurement

We propose a notion of vulnerability that is linked to the information
theoretic cost function proposed in SE_TSG_19  to characterize the disruption
and detection tradeoff incurred by the attacks. Taking the state of the system
with $k$ compromised measurements as the baseline, we quantify the
vulnerability of measurement $i$ in terms of the cost decrease induced by
attacking the correspond sensor $i$ that obtains the measurement, with
$i\in\mathcal{K}_{o}$. In the following, we define the vulnerability of a
measurement.

######  Definition 1.

The function
$\Delta:\mathcal{S}^{m}_{+}\times\mathbb{R}_{+}\times\mathbb{R}_{+}\times\mathcal{K}_{o}\rightarrow\mathbb{R}_{+}$,
where $\mathcal{K}_{o}$ is in (24), defines the vulnerability of measurement
$i$ in the following form:

$\Delta(\hbox{\boldmath$\Sigma$},\lambda,v,i)\stackrel{{\scriptstyle\Delta}}{{=}}f(\hbox{\boldmath$\Sigma$},\lambda,v,i)-f(\hbox{\boldmath$\Sigma$},\lambda,0,i),$
(31)

where the function $f$ is defined in (26).

Note that the attacker aims to minimize (26) by choosing an index $i$ and a
variance $v$, and therefore, the definition above implies that given that $k$
sensors in $\\{1,2,\ldots,m\\}\setminus\mathcal{K}_{o}$ are already attacked
in the system, the most vulnerable measurement is obtained by solving the
following minimization problem

$\min_{i\in\mathcal{K}_{o}}\Delta(\hbox{\boldmath$\Sigma$},\lambda,v,i),$ (32)

where $\mathcal{K}_{o}$ is defined in (24).

##  5 Vulnerability of Measurements

###  5.1 Vulnerability analysis of uncompromised systems

We first consider the case in which no sensors are under attacks, that is,
$k=0$ and the following holds

$\displaystyle\hbox{\boldmath$\Sigma$}=$ $\displaystyle\ \textbf{0},$ (33)
$\displaystyle\mathcal{K}_{o}=$ $\displaystyle\ \\{1,2,\ldots,m\\}.$ (34)

The attacker selects a single sensor and corrupts the corresponding
measurement with a given budget $v\leq v_{0}$. We quantify the vulnerability
of measurement $i$ in terms of $\Delta(\hbox{\boldmath$\Sigma$},\lambda,v,i)$
defined in (31).

For the uncompromised system case, the optimization problem in (32) can be
solved in closed form expression. The following theorem provides the solution.

######  Theorem 2.

The solution to the problem in (32), with
$\mathcal{K}_{o}~{}=~{}\\{1,2,\ldots,m\\}$, is

$i=\operatorname*{arg\,min}_{j\in\\{1,2,\ldots,m\\}}\left\\{\left(\hbox{\boldmath$\Sigma$}_{Y\\!Y}^{-1}\right)_{jj}\right\\},$
(35)

where $\hbox{\boldmath$\Sigma$}_{Y\\!Y}$ is in (5).

###### Proof.

We start by noting that from (33), it yields that the vulnerability of
measurement $i$ in (31) is $\Delta(\textbf{0},\lambda,v,i)$. From the equality
in (30), the function $f(\textbf{0},\lambda,0,i)$ is a constant with respect
to $i$. Hence, for $\hbox{\boldmath$\Sigma$}=\textbf{0}$, the optimization
problem in (32) is equivalent to

$\min_{i\in\mathcal{K}_{o}}f(\textbf{0},\lambda,v,i),$ (36)

where $\mathcal{K}_{o}$ is defined in (34). Let $\lambda\in\mathbb{R}_{+}$ and
$v\in\mathbb{R}_{+}$. From (30), the resulting problem in (36) is equivalent
to the following optimization problem:

$\displaystyle\min_{i\in\\{1,2,\ldots,m\\}}$
$\displaystyle(1-\lambda)\textnormal{log}\left|\hbox{\boldmath$\Sigma$}_{Y\\!Y}+v{\bf
e}_{i}{\bf e}_{i}^{\sf T}\right|-\textnormal{log}\left|v{\bf e}_{i}{\bf
e}_{i}^{\sf T}+\sigma^{2}\textbf{I}_{m}\right|$ $\displaystyle+\lambda
v\textnormal{tr}\left(\hbox{\boldmath$\Sigma$}_{Y\\!Y}^{-1}{\bf e}_{i}{\bf
e}_{i}^{\sf T}\right)$ (37) $\displaystyle=\min_{i\in\\{1,2,\ldots,m\\}}$
$\displaystyle(1-\lambda)\textnormal{log}\left|\textbf{I}_{m}+v\hbox{\boldmath$\Sigma$}_{Y\\!Y}^{-1}{\bf
e}_{i}{\bf e}_{i}^{\sf T}\right|-\textnormal{log}(v+\sigma^{2})$
$\displaystyle+\lambda
v\textnormal{tr}\left(\hbox{\boldmath$\Sigma$}_{Y\\!Y}^{-1}{\bf e}_{i}{\bf
e}_{i}^{\sf T}\right)$ (38) $\displaystyle=\min_{i\in\\{1,2,\ldots,m\\}}$
$\displaystyle(1-\lambda)\textnormal{log}\left(1+v\textnormal{tr}\left(\hbox{\boldmath$\Sigma$}_{Y\\!Y}^{-1}{\bf
e}_{i}{\bf e}_{i}^{\sf T}\right)\right)$ $\displaystyle+\lambda
v\textnormal{tr}\left(\hbox{\boldmath$\Sigma$}_{Y\\!Y}^{-1}{\bf e}_{i}{\bf
e}_{i}^{\sf T}\right),$ (39)

where the equivalence in (37) holds from plugging in
$\hbox{\boldmath$\Sigma$}=\textbf{0}$ into the equality in (30); the equality
in (38) follows from removing a constant
$(1-\lambda)\textnormal{log}\left|\hbox{\boldmath$\Sigma$}_{Y\\!Y}\right|$
from the first term; and the equality in (39) follows from the fact that
$\hbox{\boldmath$\Sigma$}_{Y\\!Y}^{-1}{\bf e}_{i}{\bf e}_{i}^{\sf T}$ is a
matrix with nonzero entries in the $i$-th column and all the other entries are
zeros.

We now proceed by defining
$t\stackrel{{\scriptstyle\Delta}}{{=}}v\textnormal{tr}\left(\hbox{\boldmath$\Sigma$}_{Y\\!Y}^{-1}{\bf
e}_{i}{\bf e}_{i}^{\sf T}\right)$, with $t\in\mathbb{R}_{+}$. Hence, the
equality in (39) is equivalent to

$\displaystyle\min_{t\in\mathbb{R}_{+}}$
$\displaystyle(1-\lambda)\textnormal{log}\left(1+t\right)+\lambda t.$ (40)

Note that (40) is monotonically increasing with respect to $t$. Therefore, the
cost function in (39) is monotonically increasing with respect to
$\textnormal{tr}\left(\hbox{\boldmath$\Sigma$}_{Y\\!Y}^{-1}{\bf e}_{i}{\bf
e}_{i}^{\sf T}\right)$. This completes the proof. ∎

From Theorem 2, it follows that the identification of the most vulnerable
measurement is independent of $\lambda$, introduced in (26), and the variance
$v$. That is, it exclusively depends on the system topology and parameters
denoted by $\hbox{\boldmath$\Sigma$}_{Y\\!Y}$ defined in (5). This result
coincides with Theorem 1 in the sense that in the attack construction for
$k=1$, the most vulnerable measurement is characterized in (19), which is
independent of the value of $\lambda$. The following corollary formalizes this
observation.

######  Corollary 1.

Consider the parameters $\hbox{\boldmath$\Sigma$}=\textbf{0}$,
$v\in\mathbb{R}_{+}$ and $\lambda\in\mathbb{R}_{+}$. The vulnerability ranking
for all measurements

${\bf s}\stackrel{{\scriptstyle\Delta}}{{=}}(s_{1},s_{2},\ldots,s_{m})$ (41)

is such that for all measurements $i$, with $i~{}\in~{}\\{1,2,\ldots,m\\}$,
$s_{i}\in\\{1,2,\ldots,m\\}$ and

$\textnormal{tr}\left(\hbox{\boldmath$\Sigma$}_{Y\\!Y}^{-1}{\bf e}_{s_{1}}{\bf
e}_{s_{1}}^{\sf
T}\right)\leq\textnormal{tr}\left(\hbox{\boldmath$\Sigma$}_{Y\\!Y}^{-1}{\bf
e}_{s_{2}}{\bf e}_{s_{2}}^{\sf
T}\right)\leq\ldots\leq\textnormal{tr}\left(\hbox{\boldmath$\Sigma$}_{Y\\!Y}^{-1}{\bf
e}_{s_{m}}{\bf e}_{s_{m}}^{\sf T}\right).$ (42)

For all $i\in\\{1,2,\ldots,m\\}$, the $i$-th most vulnerable measurement is
$s_{i}$.

###  5.2 Vulnerability index (VuIx)

The vulnerability analysis of uncompromised systems in Section 5.1 is
constrained to $k=0$. To generalize the vulnerability analysis to compromised
systems when $k>0$, in the following we propose a novel metric, that is,
vulnerability index, for all $i\in\mathcal{K}_{o}$.

######  Definition 2.

For $k\in\\{1,2,\ldots,m-1\\}$ and $\mathcal{S}_{k}$ in (23), consider the
parameters $\hbox{\boldmath$\Sigma$}\in\mathcal{S}_{k}$, $v\in\mathbb{R}_{+}$,
$\lambda\in\mathbb{R}_{+}$. Consider also the set
$\\{(i,\Delta):i\in\mathcal{K}_{o}\\}$, with $\mathcal{K}_{o}$ in (24) and

$\Delta_{i}\stackrel{{\scriptstyle\Delta}}{{=}}\Delta(\hbox{\boldmath$\Sigma$},\lambda,v,i).$
(43)

Let the vulnerability ranking

${\bf r}=(r_{1},r_{2},\ldots,r_{|\mathcal{K}_{o}|})$ (44)

be such that for all $i\in\\{1,2,\ldots,|\mathcal{K}_{o}|\\}$,
$r_{i}\in\mathcal{K}_{o}$ and moreover,

$\Delta_{r_{1}}\leq\Delta_{r_{2}}\leq\ldots\leq\Delta_{r_{|\mathcal{K}_{o}|}}.$
(45)

The vulnerability index (VuIx) of measurement $r_{j}\in\mathcal{K}_{o}$ is
$j$, that is, $\textnormal{VuIx}(r_{j})=j$.

Figure 1: Vulnerability index (VuIx) when $k=1$, SNR = 10 dB, $\lambda=2$ and
$\rho=0.1$ on the IEEE 9-bus system.

Figure 2: Vulnerability index (VuIx) when $k=2$, SNR = 10 dB, $\lambda=2$ and
$\rho=0.1$ on the IEEE 9-bus system.

Figure 3: Vulnerability index (VuIx) when $k=1$, SNR = 30 dB, $\lambda=2$ and
$\rho=0.1$ on the IEEE 9-bus system.

Figure 4: Vulnerability index (VuIx) when $k=2$, SNR = 30 dB, $\lambda=2$ and
$\rho=0.1$ on the IEEE 9-bus system.

Note that the measurement with the smallest VuIx is the most vulnerable
measurement that corresponds to the solution to the optimization problem in
(32). The proposed VuIx for $i\in\mathcal{K}_{o}$ is obtained from Algorithm
1.

Algorithm 1 Computation of Vulnerability Index (VuIx)

1:${\bf H}$ in (1);
$\sigma^{2}$ in (2);
$\hbox{\boldmath$\Sigma$}_{X\\!X}$ in (3);
$\hbox{\boldmath$\Sigma$}\in\mathcal{S}_{k}$ in (22);
$\lambda\in\mathbb{R}_{+}$ and $v\in\mathbb{R}_{+}$.

2:the VuIx for all $i\in\mathcal{K}_{o}$.

3:Set $\mathcal{K}_{o}$ in (24)

4:for $i\in\mathcal{K}_{o}$ do

5: Compute $\Delta(\hbox{\boldmath$\Sigma$},\lambda,v,i)$ in (31)

6:end for

7:Sort $\Delta(\hbox{\boldmath$\Sigma$},\lambda,v,i)$ in ascending order

8:Set ${\bf r}=(r_{1},r_{2},\ldots,r_{|\mathcal{K}_{o}|})$

9:Set the VuIx of measurement $r_{j}\in\mathcal{K}_{o}$ as $j$.

Figure 5: Vulnerability index (VuIx) when $k=1$, SNR = 10 dB, $\lambda=2$ and
$\rho=0.1$ on the IEEE 30-bus system.

Figure 6: Vulnerability index (VuIx) when $k=2$, SNR = 10 dB, $\lambda=2$ and
$\rho=0.1$ on the IEEE 30-bus system.

Figure 7: Vulnerability index (VuIx) when $k=1$, SNR = 30 dB, $\lambda=2$ and
$\rho=0.1$ on the IEEE 30-bus system.

Figure 8: Vulnerability index (VuIx) when $k=2$, SNR = 30 dB, $\lambda=2$ and
$\rho=0.1$ on the IEEE 30-bus system.

##  6 Numerical results

In this section, we numerically evaluate the VuIx of the measurements on a
direct current (DC) setting for the IEEE Test systems UoW_ITC_99 . The voltage
magnitudes are set to 1.0 per unit, that is, the measurements of the systems
are active power flow between the buses that are physically connected and
active power injection to all the buses. The Jacobian matrix ${\bf H}$ in (1)
determined by the topology of the system and the physical parameters of the
branches is generated by MATPOWER matpower . We adopt a Toeplitz model for the
covariance matrix $\hbox{\boldmath$\Sigma$}_{X\\!X}$ that arises in a wide
range of practical settings, such as autoregressive stationary processes.
Specifically, we model the correlation between state variable $X_{i}$ and
$X_{j}$ with an exponential decay parameter $\rho\in\mathbb{R}_{+}$, which
results in the entries of the matrix
$(\hbox{\boldmath$\Sigma$}_{X\\!X})_{ij}=\rho^{|i-j|}$ with
$(i,j)\in\\{1,2,\ldots,n\\}\times\\{1,2,\ldots,n\\}$. In this setting, the
VuIx of the measurements is also a function of the correlation parameter
$\rho$, the noise variance $\sigma^{2}$, and the Jacobian matrix ${\bf H}$.
The noise regime in the observation model is characterized by the signal to
noise ratio (SNR) defined as

$\textrm{SNR}\stackrel{{\scriptstyle\Delta}}{{=}}10\log_{10}\left(\frac{\textnormal{tr}({\bf
H}\hbox{\boldmath$\Sigma$}_{X\\!X}{\bf H}^{\sf T})}{m\sigma^{2}}\right).$ (46)

For all $\lambda\in\mathbb{R}_{+}$ and $v\in\mathbb{R}_{+}$, we generate a
realization of $k$ attacked indices
$\mathcal{K}_{a}\subseteq\\{1,2,\ldots,m\\}$ that is uniformly sampled from
the set of sets given by

$\tilde{\mathcal{K}}=\left\\{\mathcal{A}\subseteq\\{1,2,\ldots,m\\}:|\mathcal{A}|=k\right\\}.$
(47)

We then construct a random covariance matrix describing the existing attacks
on the system as

$\tilde{\hbox{\boldmath$\Sigma$}}=\sum_{i\in{\mathcal{K}_{a}}}{\bf e}_{i}{\bf
e}_{i}^{\sf T},$ (48)

with $\mathcal{K}_{a}\in\tilde{\mathcal{K}}$. In the numerical simulation, we
obtain the vulnerability of measurement $i$ by computing

$\Delta(\tilde{\hbox{\boldmath$\Sigma$}},\lambda,1,i),$ (49)

where $i\in\mathcal{K}_{o}$ is in (24) and $\Delta$ is defined in (31).

###  6.1 Assessment of vulnerability index (VuIx)

Fig. 2 and Fig. 2 depict the mean and variance of the VuIx obtained from
Algorithm 1 for all the measurements with SNR = 10 dB, $\lambda=2$ and
$\rho=0.1$ on the IEEE 9-bus system when $k=1$ and $k=2$, respectively. It is
observed that power injection measurements yield higher priority vulnerability
indices, which indicates that power injection measurements are more vulnerable
to data integrity attacks. Most power injection measurements correspond to
higher ranked vulnerability indices but there are instances of power flow
measurements with a higher ranked VuIx than that of some power injection
measurements. Interestingly, the power injection measurements with lower
vulnerability indices correspond to the buses that are more isolated in the
system, that is, the buses with a lower number of connections. On the other
hand, the power flow measurements with higher ranked vulnerability indices
correspond to the branches with higher admittance. The VuIx for $k=0$ obtained
in Corollary 1 is depicted for the purpose of serving as a reference to assess
the deviation when $k>0$. Interestingly, the VuIx of most measurements does
not change substantially for different values of $k$, which suggests that the
VuIx is insensitive to the state of the system.

Fig. 4 and Fig. 4 depict the mean and variance of the VuIx from Algorithm 1
for all the measurements with SNR = 30 dB, $\lambda=2$ and $\rho=0.1$ on the
IEEE 9-bus system when $k=1$ and $k=2$, respectively. Interestingly, the mean
of the VuIx for most of the measurements does not deviate significantly from
the case when $k=0$. Instead, most of the variances deviate significantly in
comparison with the cases in Fig. 2 and Fig. 2 with SNR = 10 dB. Fig. 6 and
Fig. 6 depict the results on IEEE 30-bus systems with the same setting as in
Fig. 2 and Fig. 2, respectively. Fig. 8 and Fig. 8 depict the results on IEEE
30-bus systems with the same setting as in Fig. 4 and Fig. 4, respectively.
Surprisingly, the mean of the VuIx in larger systems coincides with that
obtained for the case $k=~{}0$, which shows that the VuIx is a robust security
metric for large systems. Interestingly, the power injection measurements
corresponding to the least connected buses decrease in the VuIx when SNR = 10
dB.

###  6.2 Comparative vulnerability assessment of power flow and power
injection measurements

In Section 6.1 we have established that power injection measurements and power
flow measurements are qualitatively different in terms of the VuIx. To provide
a quantitative description of this difference, Fig. 9 depicts the probability
of a given VuIx $i\in\left\\{1,2,\ldots,m-|\mathcal{K}_{a}|\right\\}$ being
taken by a power injection measurement or a power flow measurement for the
IEEE 9-bus and 30-bus systems when $\lambda=2$, $k=2$, SNR = 30 dB and
$\rho=0.1$. Specifically, Fig. 9 depicts the probability of the following
events:

$\displaystyle{\sf Flow}_{i}\\!$ $\displaystyle\\!:\textnormal{VuIx $i$
corresponds to a power flow measurement},$ $\displaystyle{\sf Inj}_{i}\\!$
$\displaystyle\\!:\textnormal{VuIx $i$ corresponds to a power injection
measurement}.$

It is observed that in both systems, small VuIx are more likely to correspond
to power injection measurements than to power flow measurements, that is,
$\mathbb{P}[{\sf Inj}_{i}]>\mathbb{P}[{\sf Flow}_{i}]$ for small values of
$i$. Conversely, it holds that $\mathbb{P}[{\sf Inj}_{i}]<\mathbb{P}[{\sf
Flow}_{i}]$ for large values of $i$. In fact, small VuIx corresponding to
power injection measurements is with probability one, which shows that the
most vulnerable measurements in the system are always power injection
measurements. Conversely, the larger VuIx values corresponding to power flow
measurements is with probability one, which indicates that the least
vulnerable measurements are always power flow measurements.

Figure 9: Probability mass function of Vulnerability index (VuIx) for power
injection measurements and power flow measurements when $\lambda=2$, $k=2$,
SNR = 30 dB and $\rho=0.1$ on IEEE 9-bus and 30-bus systems, respectively.

Interestingly, there is a clear demarcation for each system for which
$\mathbb{P}[{\sf Inj}_{i}]$ and $\mathbb{P}[{\sf Flow}_{i}]$ change rapidly
with the VuIx value, which suggests a phase transition type phenomenon for
measurement vulnerability.

The probability of VuIx taken by power injection measurements has high
probability mass in higher priority vulnerability indices. One the other hand,
power flow measurements are with higher probability mass in low ranked VuIx.
Precisely, the probability of the vulnerability indices with higher priority
taken by power injection measurements is 1 in both IEEE 9-bus and 30-bus
systems. Meanwhile, the probability of the lower ranked vulnerability indices
taken by power flow measurements is 1. Note that the probability of medium
ranked vulnerability indices taken by power injection measurements drops
significantly, which indicates that there are some power flow measurements
that are equally as vulnerable as power injection measurements. We observe
that these power flow measurements correspond to the branches with higher
admittance. The power injection measurements with lower vulnerability indices
correspond with the buses that are isolated in the systems.

Fig. 10 depicts the distribution of VuIx for power injection measurements and
power flow measurements on the IEEE 9-bus and 30-bus systems when $\lambda=2$,
$k=2$, SNR = 30 dB and $\rho=0.1$. Specifically, Fig. 10 depicts the
probability mass function of the following events:

$\displaystyle{\sf VuIx(Flow)}=i\\!$ $\displaystyle\\!:\textnormal{VuIx for
power flow measurements is $i$},$ $\displaystyle{\sf VuIx(Inj)}=i\\!$
$\displaystyle\\!:\textnormal{VuIx for power injection measurements is $i$}.$

Power injection measurements have much higher probability with high ranked
VuIx. VuIx, whereas power flow measurements have much higher probability with
low ranked VuIx. It is worth noting that the probability mass functions are
close to uniform for high and low vulnerability index ranges. This suggests
that the most vulnerable measurements in the system are contained with high
probability in a subset of the power injection measurements. Conversely, the
least vulnerable measurements comprise the majority of the power flow
measurements with no apparent preference over the majority. Interestingly, in
the 30-bus system, the probability of lowest ranked VuIx for power flow
measurements experiences a sharp increase.

##  7 Conclusion

In this paper, we have proposed, from a fundamental perspective, a novel
security metric referred to as vulnerability index (VuIx) that characterizes
the vulnerability of power system measurements to data integrity attacks. We
have achieved this by embedding information theoretic measures into the metric
definition. The resulting VuIx framework evaluates the vulnerability of all
the measurements in the systems and enables the operator to identify those
that are more exposed to data integrity threats. We have tested the framework
for IEEE test systems and concluded that power injection measurements are more
vulnerable to data integrity attacks than power flow measurements.

Figure 10: Probability of Vulnerability index (VuIx) corresponds to power
injection measurements and power flow measurements when $\lambda=2$, $k=2$,
SNR = 30 dB and $\rho=0.1$ on IEEE 9-bus and 30-bus systems, respectively.

## References

  * [1] Grainger JJ, Stevenson WD.  Power system analysis.  McGraw-Hill; 1994.
  * [2] Abur A, Exposito AG.  Power system state estimation: Theory and implementation.  CRC press; 2004.
  * [3] Wang W, Lu Z.  Cyber security in the smart grid: Survey and challenges.  Computer networks. 2013 Jan;57(5):1344-71.
  * [4] Jaquith A.  Security metrics: replacing fear, uncertainty, and doubt.  Pearson Education; 2007.
  * [5] Mell P, Scarfone K, Romanosky S.  Common vulnerability scoring system.  IEEE Security & Privacy. 2006;4(6):85-9.
  * [6] Pallitteri VY, Brewer TL.  Guidelines for Smart Grid Cybersecurity.  NIST Interagency/Internal Report (NISTIR), National Institute of Standards and Technology; 2014.  Available from: https://doi.org/10.6028/NIST.IR.7628r1.
  * [7] Pendleton M, Garcia-Lebron R, Cho JH, Xu S.  A survey on systems security metrics.  ACM Computing Surveys. 2017 Dec;49(4):1-35.
  * [8] Venkataramanan V, Hahn A, Srivastava A.  CP-SAM: Cyber-physical security assessment metric for monitoring microgrid resiliency.  IEEE Trans Smart Grid. 2022 Mar;11(2):1055-65.
  * [9] Liu Y, Ning P, Reiter MK.  False data injection attacks against state estimation in electric power grids.  ACM Trans Info Syst Sec. 2011 May;14(1):1-33.
  * [10] Cui S, Han Z, Kar S, Kim TT, Poor HV, Tajer A.  Coordinated data-injection attack and detection in the smart grid: A detailed look at enriching detection solutions.  IEEE Signal Process Mag. 2012 Aug;29(5):106-15.
  * [11] Ozay M, Esnaola I, Vural FTY, Kulkarni SR, Poor HV.  Sparse attack construction and state estimation in the smart grid: Centralized and distributed models.  IEEE J Sel Areas Commun. 2013 Jul;31(7):1306-18.
  * [12] Esnaola I, Perlaza SM, Poor HV, Kosut O.  Maximum distortion attacks in electricity grids.  IEEE Trans Smart Grid. 2016 Jul;7(4):2007-15.
  * [13] Bretas A, Bretas N, London Jr JB, Carvalho B.  Cyber-physical power systems state estimation.  Elsevier; 2021.
  * [14] Sun K, Esnaola I, Perlaza SM, Poor HV.  Stealth attacks on the smart grid.  IEEE Trans Smart Grid. 2019 Aug;11(2):1276-85.
  * [15] Ye X, Esnaola I, Perlaza SM, Harrison RF.  Information theoretic data injection attacks with sparsity constraints.  In: Proc. IEEE Int. Conf. on Smart Grid Comm. Tempe, AZ, USA: IEEE; 2020\. .
  * [16] Genes C, Esnaola I, Perlaza SM, Ochoa LF, Coca D.  Robust recovery of missing data in electricity distribution systems.  IEEE Trans Smart Grid. 2018 Jun;10(4):4057-67.
  * [17] Shomorony I, Avestimehr AS.  Worst-case additive noise in wireless networks.  IEEE Trans Inf Theory. 2013 Jun;59(6):3833-47.
  * [18] Lévy P.  Propriétés asymptotiques des sommes de variables aléatoires enchaînées.  J Math Pures Appl. 1935;14:109-28.
  * [19] Cramér H.  Über eine Eigenschaft der normalen Verteilungsfunktion.  Math Z. 1936;41:405-14.
  * [20] Cover TM, Thomas JA.  Elements of information theory.  John Wiley & Sons; 1999.
  * [21] Seber GA.  A matrix handbook for statisticians. vol. 15.  John Wiley & Sons; 2008.
  * [22] UW.  Power Systems Test Case Archive.  University of Washington; 1999.  Available from: https://labs.ece.uw.edu/pstca/.
  * [23] Zimmerman RD, Murillo-Sánchez CE, Thomas RJ.  MATPOWER: Steady-state operations, planning, and analysis tools for power systems research and education.  IEEE Trans Power Syst. 2010 Feb;26(1):12-9.


# Solving the Conflict-Free Electric Vehicle Routing Problem Using SMT Solvers

Sabino Francesco Roselli1 and Martin Fabian1 and Knut Åkesson1 We gratefully
acknowledge financial support from Chalmers AI Research Centre (CHAIR) and AB
Volvo (Project ViMCoR). 1Department Electrical Engineering, Chalmers
University of Technology, Göteborg, Sweden {rsabino, fabian<EMAIL_ADDRESS>
###### Abstract

The Vehicle Routing Problem (VRP) is the combinatorial optimization problem of
computing routes to serve customers while minimizing a cost function,
typically the travelled distance or the number of vehicles required for a
given performance. Industrial applications of the problem in manufacturing
plants is the scheduling and routing of Automated Guided Vehicles (AGVs) to
deliver material between storage areas and assembly stations. For large fleets
of vehicles it is necessary to take the limited space of plant floor into
account during scheduling and routing in order to limit the number of AGVs
that are at certain locations at a given time. In addition, AGVs are most
often powered by batteries and therefore have limited operating range and non-
negligible charging time that will also affect the scheduling and routing
decisions. In this paper we provide a model formulation for the scheduling and
routing of AGVs with given time-windows for delivering material, restricted by
capacity constraints on the path network, and with the need for periodic
recharge. The problem is modelled and solved using optimizing Satisfiability
Modulo Theory (SMT) solvers. The approach is evaluated on a set of generated
problem instances, showing that the solver can handle medium size instances in
a reasonable amount of time.

##  I Introduction

Automation is increasingly used material handling and delivery [1]. Many
companies are implementing _just-in-time_ delivery to workstations [2]. This
can be implemented by a fleet of Automated Guided Vehicles (AGVs) that pick up
and deliver the material and components to the due workstation just-in-time
for being assembled.

Just-in-time delivery benefits from using optimization to solve the scheduling
and routing problems in order to minimize the travelled-distance and using as
few vehicles as necessary while still delivering all material within the given
time-windows. Thus, an optimization problem can be formulated, typically
variants of the Vehicle Routing Problem (VRP) [3], a combinatorial
optimization problem of designing routes for vehicles to serve customers so
that a cost function. Typically either the number of vehicles or the total
travelled distance, or a combination of them, is minimized. As the delivery of
components to the assembly line has to be done just-in-time, vehicles cannot
arrive too late or too early to the workstation. There exists a variation of
the standard VRP that accounts for this additional constraint; the vehicle
routing problem with time windows (VRPTW) [4]. AGVs are often powered by
batteries and may need to recharge during service. There exists a variation of
the standard VRP that accounts for limited operating range of the vehicles and
charging stations; the electric vehicle routing problem (E-VRPTW) [5].
Finally, another important aspect of the system is the capacity of the road
segments in terms of the number of vehicles a segment can accommodate.
Basically, if a segment is too narrow, two vehicles cannot travel in opposite
directions through it. As for the previous cases, there exists a variation of
the standard VRP that accounts for limited capacity of the road segments; the
dispatch and conflict-free routing problem (DCFRP) [6].

In this paper we are going to design a model to address all these features at
once, the henceforth called the _Conflict-Free Electric Vehicle Routing
Problem_ , or CF-EVRP. The cost function we use to evaluate the goodness of a
solution is the total travelled distance.

For all the variations mentioned in the previous paragraph, there exist Mixed
Integer Linear Programming (MILP) models. Using MILP to handle the problem
offers high flexibility in terms of extending the model as new constraints are
required. On the other hand, the computational complexity increases as the
model grows (in terms of number of customer and/or vehicles and, therefore,
variables and constraints), especially when the new requirements involve
combinatorial logical constraints, since the linearization of logical
constraints involves a large number of additional variables and _auxiliary
constraints_. For this reason MILP solvers are questionable for large scale
problems and specific-purpose algorithms involving local search [7], column
generation [8] or stochastic methods [9, 10] are used instead.

Our previous results [11] showed that Satisfiability Modulo Theory (SMT, [12,
13]) solvers outperformed MILP solvers on combinatorial problems involving
many logical constraints. A model to address the system previously discussed
involves many of these logical constraints and, therefore, it is reasonable to
assume that an SMT solver could handle large problem instances of the CF-EVRP.

Moreover, there is recent work focusing on fleets of electric vehicles [14,
15], as well as conflict-free routing [16, 17]. However, to the best of our
knowledge, there is no work focusing on both aspects at the same time, and SMT
has not been tested on this particular problem. We have therefore designed
problem instances of the CF-EVRP to test the performance of the state-of-the-
art SMT solver Z3 [18].

The contributions of this paper are: (i) Designing a model for the CF-EVRP in
order to represent a just-in-time material handling system based on AGVs that
takes into account time windows for delivery, limited operating range of the
vehicles and distributed charging stations, as well as limitations in the road
segments capacity in terms of number of vehicles that can transit at the same
time. (ii) Design of problem instances to test the model’s correctness and
performance. (iii) Benchmark of a state-of-the-art SMT solver over the
generated problem instances.

The paper is organized as follows: in Section II a formal description of the
problem is provided, together with a mathematical model and detailed
descriptions of its constraints. In Section III, the results of the analysis
over the set of generated problem instances is presented. Finally, conclusions
are drawn in Section IV.

##  II Problem Formulation

In the variations of the VRP mentioned in Section I, the real-world map is
abstracted and it is possible to travel directly from any customer to any
other customer. For this reason it is possible to represent the problem as a
complete undirected graph, where each customer is a node and two nodes are
connected with each other by a weighted edge representing their distance.
Customers are located at specific coordinates and the weight of the edge
connecting two customers is their Euclidean or Manhattan distance.

However, in the CF-EVRP the capacity constraints on the road segments are
taken into account and therefore, it is necessary to know where every vehicle
is located at each time. Moreover, there can be several ways to travel from
one customer to another, each with a different length. For these reasons it is
not possible to make an abstraction of the real world map. Instead, a
weighted, directed graph represents the plant layout; the nodes are the
customers or the depot stations (as in the classical VRP), but they can also
be intersections between road segments, while the edges represent the road
segments and the edges’ bi-dimensional weight represents their length and
their capacity in terms of vehicles that can travel through it at the same
time.

The following definitions are provided:

  * •

Node: a location in the plant.

  * •

Edge: a road segment that connects two nodes.

  * •

Hub: a node that can accommodate an unbounded number of vehicles at the same
time; nodes that are not hubs can only accommodate a single vehicle.

  * •

Vehicle: a transporter, e.g. an AGV, that is able to move between locations in
the plant and perform pickup and delivery operations.

  * •

Task: either a pickup or a delivery operation. A task is always associated
with a location in the plant where the task is executed.

  * •

Job: one or more pickup tasks and one delivery task. Pickup may have
precedence constraints among them, while the delivery task for a job always
happens after all pickups for that job are completed.

  * •

Time Window: a pair of integer values indicating the earliest and latest time
at which a vehicle can execute a task. All tasks have a time window which,
unless explicitly specified, is the entire time horizon.

  * •

Depot: node at which one or more vehicles start and must return to after
completing the assigned jobs.

The requirements of the problem are summarized as follows:

  * •

all jobs have to be completed; for a job to be completed a vehicle has to be
assigned to it and visit the locations of the job’s tasks within their
respective time windows.

  * •

when a vehicle is assigned to more than one job, it has to execute the job’s
tasks sequentially; all tasks of a job must be completed before it can execute
any task of another job.

  * •

there can be an arbitrary number of depots; vehicles have to return to the
depot they were dispatched from.

  * •

vehicles are powered by batteries with limited capacity but with the ability
to recharge at a charging station. It is assumed that charging and discharging
of the batteries is linear with respect to distance.

  * •

different road segments in the plant have different capacities in terms of
number of vehicles they can accommodate.

  * •

pickup and delivery duration (the time when the pickup and delivery tasks
respectively are executed) are considered to be zero, as these times can be
considered negligible compared to the travelling time.

  * •

not all vehicles are eligible to execute all jobs.

###  II-A Notation

In order to model the VRP described above, we define the following notation
that is used to define constraints:

$\displaystyle N\textrm{: the set of nodes}$ $\displaystyle N_{H}\subseteq
N\textrm{: set of \emph{hub} nodes}$ $\displaystyle E\subseteq N\times
N\textrm{: the set of edges}$ $\displaystyle d_{nn^{\prime}},\ \forall
n,n^{\prime}\in N\textrm{: the length of the edge connecting}$ nodes $n$ and
$n^{\prime}$ 111The need to increase the value _B_ is discussed on page II-C.
$\displaystyle g_{nn^{\prime}},\forall n,n^{\prime}\in N\textrm{: the capacity
of the edge connecting}$ nodes $n$ and $n^{\prime}$ $\displaystyle
A_{n}\subset N,\ \forall n\in N\textrm{: the set of nodes having an incoming}$
edge from node $n$ $\displaystyle h=max_{n,n^{\prime}\in
N}{(d_{nn^{\prime}})}\textrm{: the value of the longest road-segment}$
$\displaystyle B\textrm{: an upper bound representing the time when all jobs}$
should be finished $\displaystyle T=B+h\textrm{: the time horizon}$
$\displaystyle J:\textrm{the set of jobs}$ $\displaystyle K_{j},\ \forall j\in
J\textrm{: the tasks of job $j$}$ $\displaystyle L_{jk},\in N\ \forall j\in
J,k\in K_{j}\textrm{: the location of task $k$ of job $j$}$ $\displaystyle
V\textrm{: the set of all vehicles}$ $\displaystyle V_{j}\subseteq V,\ \forall
j\in J\textrm{: the set of vehicles eligible for job $j$}$ $\displaystyle
s_{i}\in N,\ \forall i\in V\textrm{: starting node of vehicle $i$}$
$\displaystyle\mathit{OR}\textrm{: the maximum operating range of the
vehicles}$ $\displaystyle C\textrm{: the charging coefficient}$ $\displaystyle
D\textrm{: the discharging coefficient}$ $\displaystyle F\textrm{: the set of
charging stations}$ $\displaystyle P_{jk}\subset K_{j},\ \forall j\in J,k\in
K_{j}\textrm{: the set of tasks to execute}$ before task $k$ of job $j$
$\displaystyle l_{jk},\ \forall j\in J,k\in K_{j}\textrm{: the time window's
lower bound for}$ task $k$ of job $j$ $\displaystyle u_{jk},\ \forall j\in
J,k\in K_{j}\textrm{: the time window's upper bound for}$ task $k$ of job $j$

The set of variables used to build the model are:

  * $x_{ij}$: Boolean variables that evaluate to true if the $i$-th vehicle is assigned to the $j$-th job

  * $y_{ijk}$: Boolean variables that evaluate to true if the $i$-th vehicle is assigned to the $k$-th task of the $j$-th job

  * $z_{ijkt}$: Boolean variables that evaluate to true if the $i$-th vehicle serves the $k$-th task of the $j$-th job at time $t$

  * $rc_{it}$: Integer positive variables that keep track of the remaining operating range of vehicle _i_ at time _t_

  * $at_{int}$: Boolean variables that evaluate to true if the $i$-th vehicle is at the $n$-th node at time _t_

  * $move_{int}$: Boolean variables that evaluate to true if the $i$-th vehicle decides to move to the $n$-th node at time _t_

The following logical operators are used to express cardinality constraints
[19]:

  * $\textrm{{AMO}}(a):$ at most one variable of the set $a$ evaluates to true;

  * $\textrm{{AMN}}(a,n):$ at most $n$ variables of the set $a$ evaluate to true;

  * $\textrm{If}(c,o_{1},o_{2}):$ if condition $c$ is _true_ returns $o_{1}$, else returns $o_{2}$.

For further clarity, the constraint sets have been divided into sections,
depending on the features of the system they model. The _job assignment_ is
concerned with the allocation of the vehicles to the jobs; the _vehicle
movements_ defines the rules for the vehicles to move through the plant; the
_conflict-free routing_ deals with the capacity constraints on the road
segments; the _battery management_ defines how the vehicles’ batteries
discharge and recharge according to their state.

###  II-B Jobs Assignment

The following constraints are about assignment of jobs to vehicles, making
sure that all jobs are assigned one and only one vehicle according to the
specifications, and that such vehicle performs service to the job’s tasks:

$\displaystyle z_{ijkt}\implies$ $\displaystyle at_{iL_{jk}t}$
$\displaystyle\forall i\in V,j\in J,k\in K_{j},t=0,\dots,T$ (1) $\displaystyle
x_{ij}\implies$ $\displaystyle\bigwedge_{\begin{subarray}{c}i^{\prime}\in V\\\
i\neq i^{\prime}\end{subarray}}{\neg{x_{i^{\prime}j}}}\qquad\forall i\in
V,j\in J$ (2) $\displaystyle y_{ijk}\implies$
$\displaystyle\bigwedge_{\begin{subarray}{c}i^{\prime}\in V\\\ i\neq
i^{\prime}\end{subarray}}{\neg{y_{i^{\prime}jk}}}\qquad\forall i\in V,j\in
J,k\in K_{j}$ (3) $\displaystyle y_{ijk}\implies$
$\displaystyle\bigvee_{t=0}^{T}{z_{ijkt}}\qquad\forall i\in V,j\in J,k\in
K_{j}$ (4) $\displaystyle z_{ijk_{1}t}\implies$
$\displaystyle\bigwedge_{\begin{subarray}{c}k_{2}\in P_{jk_{1}}\\\
t^{\prime}=t,\dots,T\end{subarray}}{\neg{z_{ijk_{2}t^{\prime}}}}$
$\displaystyle\forall i\in V,j\in J,k\in K_{j},t=0,\dots,T$ (5) $\displaystyle
x_{ij}\implies$ $\displaystyle\bigwedge_{k\in K_{j}}{y_{ijk}}\qquad\forall
i\in V,j\in J$ (6) $\displaystyle y_{ijk}\implies$
$\displaystyle\bigvee_{t=l_{jk},\dots,u_{jk}}{z_{ijkt}}\qquad\forall i\in
V,j\in J,k\in K_{j}$ (7) $\displaystyle\neg{x_{ij}}$
$\displaystyle\qquad\forall i\notin V_{j},j\in J$ (8)

(1) states that in order to serve a task at time $t$, a vehicle has to be at
the task location at time $t$; (2) states that only one vehicle can be
assigned to a job; similarly, (3) states that only one vehicle can be assigned
to a task; (4) states that, in order to execute a task, a vehicle must be at
the task location at some point in time; (5) enforces the precedence
constraint among tasks of the same job; (6) states that, when assigned to a
job, a vehicle must execute all tasks of the job; (7) states that a vehicle
must execute tasks within their time windows; (8) states that only eligible
vehicles can execute the corresponding job.

###  II-C Vehicles’ Movements

The following constraints are related to the vehicles’ movements, specifying
where vehicles can go given their current location, how long it will take to
reach the next location, and how being at a location relates to performing a
service at that location:

$\displaystyle at_{is_{i}0}\quad$ $\displaystyle\forall i\in V$ (9)
$\displaystyle at_{is_{i}B}\quad$ $\displaystyle\forall i\in V$ (10)
$\displaystyle\textrm{AMO}_{n\in N}{(at_{int})}\quad$ $\displaystyle\forall
i\in V,t=0,\dots,T$ (11) $\displaystyle\textrm{AMO}_{i\in V}{(at_{int})}\quad$
$\displaystyle\forall n\in N\setminus N_{H},t=0,\dots,T$ (12)
$\displaystyle\textrm{AMO}_{n\in N}{(move_{int})}\quad$ $\displaystyle\forall
i\in V,t=0,\dots,T$ (13) $\displaystyle at_{int}\implies\neg{move_{int}}\quad$
$\displaystyle\forall i\in V,n\in N,t=0,\dots,T$ (14) $\displaystyle
at_{int}\implies\bigwedge_{n^{\prime}\notin A_{n}}$
$\displaystyle{\neg{move_{in^{\prime}t}}}$ $\displaystyle\forall i\in V,n\in
N,t=0,\dots,T$ (15) $\displaystyle(at_{int}\wedge\bigwedge_{n^{\prime}\in N}$
$\displaystyle{\neg{move_{in^{\prime}t}}})\implies at_{int+1}$
$\displaystyle\forall i\in V,$ $\displaystyle n\in N,t=0,\dots,T-1$ (16)
$\displaystyle(at_{int}\wedge move_{in^{\prime}t})\implies$
$\displaystyle(\bigwedge_{\begin{subarray}{c}n^{\prime\prime}\in N\\\
t^{\prime}=t+1,\dots,t+d_{nn^{\prime}}-1\end{subarray}}$
$\displaystyle{\neg{at_{in^{\prime\prime}t^{\prime}}}}\wedge
at_{in^{\prime}t+d_{nn^{\prime}}})$ $\displaystyle\forall i\in
V,(n,n^{\prime})$ $\displaystyle\in E,t=0,\dots,T-d_{nn^{\prime}}$ (17)
$\displaystyle(z_{ij_{1}k_{1}t_{1}}\wedge z_{ij_{1}k_{2}t_{2}})\implies$
$\displaystyle\bigwedge_{\begin{subarray}{c}j_{2}\in J,j_{1}\neq j_{2}\\\
k^{\prime}\in K_{j_{2}}\\\
t^{\prime}=t_{1}+1,\dots,t_{2}\end{subarray}}{\neg{z_{ij_{2}k^{\prime}t^{\prime}}}}$
$\displaystyle\forall i\in V,j_{2}\in J,k_{1},k_{2}$ $\displaystyle\in
K_{j_{1}},k_{1}\neq k_{2},$ $\displaystyle t_{1}=0,$
$\displaystyle\dots,T,\,t_{2}=t_{1},\dots,T$ (18)

(9) states that all vehicles are at their respective depots at time zero; (10)
states that all vehicles are at the depot at time $B$, these constraints force
the vehicles to return to the depot after they execute the jobs they are
assigned to; (11) states that a vehicle can be at most in one location at a
time. Normally, it would make sense to state that a vehicle should be
_exactly_ in one place at a time. However, given the way the vehicle movement
is modeled, this constraint needs to be relaxed because moving can take more
than one time-step; (13) states that a vehicle can move to at most one
location at a time; (14) forbids a vehicle to move to the node where it
already is, and (15) states that if vehicle $i$ at time-step $t$ is at node
$n$, it cannot move to any non-adjacent node; (12) states that at most one
vehicle can, at the same time-step, be at a node that is not a _hub_ ; (16)
states that, if a vehicle is not moving at a certain time-step, it will be at
the same location at the next time-step; (17) states that, if a vehicle moves
to a new location, it will be at no node while it is transiting on the edge
connecting the start and arrival nodes, and it will be at the arrival node
after as many steps as it takes to traverse the edge, defined by the edge’s
weight. It is for this reason that constraint (11) needs to be relaxed. Also,
in (17) $t^{\prime}=t+1,\dots,t+d_{nn^{\prime}}-1$, therefore
$t=0,\dots,T-d_{nn^{\prime}}$. This means that there are no constraints
applying to events taking place in the interval $[T-d_{nn^{\prime}},T]$. For
this reason it is necessary to extend the time horizon to $T=B+h$, the steps
$B,\dots,h$ are not actually used in the model but this way it is guaranteed
that all constraints apply for $T=0,\dots,B$; (18) states that if a vehicle is
executing two tasks $k_{1}$ and $k_{2}$ of job $j$, then the vehicle cannot
execute a task of any other job different from $j$ in between. This
constraints is used to make sure that if a vehicle is selected for multiple
jobs, it has to execute them in sequence (i.e., finish all tasks of a job
before executing any other task).

###  II-D Conflict-Free Routing

The following constraints are related to the conflict-free routing, specifying
how one or more vehicle’s location affects the other vehicles’ movements:

$\displaystyle\bigwedge_{i\in\bar{V}}{at_{int}}$
$\displaystyle\implies\bigwedge_{n^{\prime}\in
A_{n}}{\textrm{AMN}_{i\in\bar{V}}}(move_{in^{\prime}t},g_{nn^{\prime}})$
$\displaystyle\\{\bar{V}\in$ $\displaystyle 2^{V}\ \big{|}\
|\bar{V}|>g_{nn^{\prime}}\\},n\in N,t=0,\cdots,T$ (19)
$\displaystyle\bigwedge_{i_{1}\in\bar{V_{1}}}{at_{i_{1}nt}}$
$\displaystyle\wedge\bigwedge_{i_{2}\in\bar{V_{2}}}{at_{i_{2}n^{\prime}t}}\wedge\bigwedge_{i_{1}\in\bar{V_{1}}}{move_{i_{1}n^{\prime}t}}\implies$
$\displaystyle\bigwedge_{t^{\prime}=t}^{t+d_{nn^{\prime}}}$
$\displaystyle\textrm{AMN}_{i_{2}\in\bar{V_{2}}}{\left(move_{i_{2}nt^{\prime}},g_{nn^{\prime}}-|\bar{V_{1}}|\right)}$
$\displaystyle\\{\bar{V_{1}}$ $\displaystyle\in 2^{V}\ \big{|}\
|\bar{V_{1}}|\leq g_{nn^{\prime}}\\},t=0,\cdots,T-d_{nn^{\prime}}$
$\displaystyle\\{\bar{V_{2}}$ $\displaystyle\in 2^{V}\ \big{|}\
|\bar{V_{2}}|\geq g_{nn^{\prime}}-|\bar{V_{1}}|\\},(n,n^{\prime})\in E$ (20)

(19) states that if _n_ vehicles are on the same node, they cannot transit on
the same edge having capacity $n-1$ at the same time-step; (20) states that if
_n_ vehicles are on two adjacent nodes and decide to traverse the edge
connecting them and whose capacity is $n+m$, at most _m_ vehicles can traverse
the edge from the other node before the _n_ vehicles are done transiting.

###  II-E Battery Management

The following constraints are related to the battery management, i.e., where
the vehicles can charge, how fast, and how far they can travel before needing
to recharge:

$\displaystyle rc_{it}\geq 0\wedge rc_{it}\leq\mathit{OR}\qquad\quad$
$\displaystyle\forall i\in V,t=0,\dots,T$ (21) $\displaystyle(at_{int}\wedge
move_{in^{\prime}t})\implies\qquad$
$\displaystyle\bigwedge_{t^{\prime}=t+1,\dots,t+d_{nn^{\prime}}+1}$
$\displaystyle{rc_{it^{\prime}}=rc_{it^{\prime}-1}-D}$ $\displaystyle\forall
i\in V,(n,n^{\prime})\in E,$ $\displaystyle t=0,\dots,T-d_{nn^{\prime}}$ (22)
$\displaystyle(at_{int}\wedge\bigwedge_{n^{\prime}\in
N}{\neg{move_{in^{\prime}t}}})\implies$ $\displaystyle rc_{it+1}=rc_{it}$
$\displaystyle\forall i\in V,n\in N$ $\displaystyle\setminus F,t=0,\dots,T-1$
(23) $\displaystyle(at_{int}\wedge\bigwedge_{n^{\prime}\in
N}{\neg{move_{in^{\prime}t}}})\implies$ $\displaystyle rc_{it+1}=rc_{it}+C$
$\displaystyle\forall i\in V,n\in$ $\displaystyle F,t=0,\dots,T-1$ (24)

(21) defines the domain of the $rc_{it}$ variables; (22) states that, if a
vehicle is moving, its energy decreases at each time-step according to the
discharge coefficient; (23) states that if a vehicle is not moving, its energy
level does not change. It is necessary to specify that the vehicle is at some
node because if the vehicle is at no node, it is moving; (24) states that if a
vehicle is at a charging station, its energy level increases at each time-step
according to the charging coefficient.

###  II-F Objective Function

Finally, the cost function to minimize (25) is the total travelled distance,
calculated as the cumulative distance travelled by each and every vehicle:

$\displaystyle\sum{\textrm{If}((at_{int}\wedge
move_{in^{\prime}t}),d_{nn^{\prime}},0)}$ $\displaystyle\forall i\in V,n\in N$
$\displaystyle,t=0,\dots,T-1$ (25)

##  III Computational Analysis

Figure 1: Test Case with five vehicles _a_ to _e_ that have to execute six
jobs _1_ to _5_ and _X1_ , each consisting of a pickup with time window [0,T]
and a delivery with time window [20,25].

In order to evaluate the model presented in Section II we generated a set of
problem instances, defined by a weighted, directed graph that represents a
hypothetical plant layout; The number of nodes is in this case a parameter
that affects the problem size, as well as the number of edges. In this work we
restricted the capacity of all road-segments to be _one_. For the graph it is
specified which nodes are _hubs_. Another parameter of the instance is the
number of vehicles available; to each vehicle a starting node is assigned;
more vehicles can start from the same node (in this case the node must be a
_hub_). A list of jobs to execute is given; each job is characterized by a
number of tasks and the vehicle(s) eligible to execute it. In turn, each task
is characterized by a location (node), a precedence list, and a time window.
Other parameters of the instance are the charging and discharging
coefficients, as well as the operating range of the vehicles and the time
horizon, by when all jobs should be finished and the vehicles should be back
at their starting locations.

TABLE I: Evaluation of the CF-EVRP over a set of generated problem instances.
Instances are sorted by the parameters N-V-J (nodes, vehicles, jobs), value of
edge reduction, and time horizon. For each resulting class, the number of
instances solved within the time limit of 10800 seconds (Sol) is reported, as
well as the average generation time (GT, in seconds), and the amount (SC) and
average solving time (ST, also in seconds) for the instances that turned out
to be satisfiable (on the left), and unsatisfiable (on the right). The symbol
“-” means that no instance for that category was either satisfiable or
unsatisfiable, depending on where the symbol appears.

Time Horizon |  | N-V-J
---|---|---
| 15-3-5 | 25-4-7 | 35-6-8
| Edge Reduction | Edge Reduction | Edge Reduction
| 0 | 25 | 50 | 0 | 25 | 50 | 0 | 25 | 50
15 | Sol | 5/5 | 5/5 | 5/5 | 5/5 | 5/5 | 5/5 | 5/5 | 5/5 | 5/5
GT | 12 | 14 | 17 | 80 | 96 | 59 | 87 | 197 | 137
ST | 2 | 3 | 3 | 3 | - | 5 | - | 34 | - | 35 | - | 15 | - | 33 | - | 108 | - | 218
| SC | 3 | 2 | 2 | 3 | 0 | 5 | 0 | 5 | 0 | 5 | 0 | 5 | 0 | 5 | 0 | 5 | 0 | 5
20 | Sol | 5/5 | 5/5 | 5/5 | 5/5 | 5/5 | 5/5 | 5/5 | 5/5 | 5/5
GT | 22 | 31 | 34 | 117 | 148 | 133 | 149 | 241 | 246
ST | 12 | 2 | 22 | 4 | - | 169 | 281 | 34 | 1517 | 107 | - | 1669 | 1619 | 37 | 2375 | 144 | - | 818
| SC | 4 | 1 | 4 | 1 | 0 | 5 | 3 | 2 | 2 | 3 | 0 | 5 | 2 | 3 | 1 | 4 | 0 | 5
25 | Sol | 5/5 | 5/5 | 4/5 | 5/5 | 4/5 | 2/5 | 2/5 | 4/5 | 5/5
GT | 38 | 47 | 55 | 192 | 222 | 233 | 257 | 295 | 391
ST | 19 | 29 | 69 | 61 | 204 | 39 | 1380 | 292 | 3327 | 246 | 1046 | 100 | 1911 | - | 2801 | 2551 | 7754 | 2826
| SC | 4 | 1 | 4 | 1 | 2 | 2 | 4 | 1 | 3 | 1 | 1 | 1 | 2 | 0 | 1 | 3 | 1 | 4
30 | Sol | 5/5 | 5/5 | 3/5 | 4/5 | 3/5 | 2/5 | 1/5 | 2/5 | 2/5
GT | 64 | 71 | 97 | 371 | 339 | 414 | 568 | 591 | 607
ST | 557 | - | 549 | - | 987 | 1152 | 5995 | - | 3916 | - | 9325 | 3414 | 6052 | - | 9812 | 703 | 9503 | 629
| SC | 5 | 0 | 5 | 0 | 2 | 1 | 4 | 0 | 3 | 0 | 1 | 1 | 1 | 0 | 1 | 1 | 1 | 1

Figure 1 shows what a problem instance may look like; there are five depot
stations located at nodes A, E, M, U, and Y, which are also the charging
stations in the problem. There is no _hub_. Vehicles are initially located as
indicated in the figure and they have to return to their initial location
before the last time-step. All edges in this problem have weight _one_ , the
charging and discharging coefficients are both _one_ , the time horizon _T_ is
equal to _35_ , and the vehicles operating range is _15_. Tasks are reported
in the figure next to their location; for instance, _task 2_ of _job 1_ is at
node _I_. Time windows, if they exist, are reported below the task name, in
between square brackets; for instance, the time window for _task 1_ of _job 2_
is _[20,25]_. In this problem jobs _1_ to _5_ can only be executed by vehicles
_a_ to _e_ respectively, while _Job x1_ can be executed by any vehicle. For
this particular instance the solver took about 90 seconds to generate the
model and 540 seconds to find the first feasible solution, which correspond to
a total travelled distance of _102_. Jobs _1_ to _5_ are executed by vehicles
_a_ to _e_ respectively, while _Job x1_ is executed by _Vehicle d_ after
executing _Job 4_ and recharging at _Node A_.

In the computational analysis we evaluated 160 problem instances structured as
explained in the previous paragraphs. The parameters we used for the problems
generation are:

  * •

_N-V-J_ (nodes, vehicles, and jobs, respectively). _N-V-J_ can be either
_15-3-5_ , or _25-4-7_ , or _35-6-8_. The corresponding layouts are grids of
size $5\times 3$, $5\times 5$ (as in Figure 1), and $5\times 7$. Each job is
composed by one pickup and one delivery;

  * •

_Edge Reduction_ can be either _0_ , or _25_ or _50_. The value _0_ correspond
to a complete grid, as in Figure 1; the value _25_ corresponds to a grid where
3, 6, and 9 pairs of edges (i.e. $(A,B)$ and $(B,A)$) are removed for _N_
equal to 15, 25, and 35 respectively; the value _50_ corresponds to a grid
where 6, 12, and 18 pairs of edges are removed for _N_ equal to 15, 25, and 35
respectively.

  * •

Time Horizon can be either _15_ , _20_ , _25_ , or _30_.

  * •

All edges have weight equal to _one_.

For each combination of the above-mentioned parameters, five different
problems are generated by randomly assigning the starting node for each type
of vehicle, the operating range, the charging coefficient, the number of
vehicles available per each type, the type of vehicle required for each job,
and the location of the jobs’ tasks.

Problem instances are designed to be on average satisfiable if the capacity
constraints are not included. This can be achieved by acting on the operating
range of the vehicles and on the time windows of the jobs. The goal is to
produce both satisfiable and unsatisfiable problem instances in order to
evaluate the performance of the solver in both cases. Moreover,
unsatisfiability can be caused by different factors, such as too tight time
windows and/or time horizon, not sufficient operating range, or too low
connectivity of the graph (in this case there are not enough paths to reach
the desired locations without getting into another vehicle’s way). By tuning
the previously defined parameters it is possible to trigger different type of
unsatisfiability and evaluate the solver performance accordingly.

For the analysis we used Z3 4.8.9. The time limit is set to 10800 seconds
(three hours); the model generation time is measured separately, since it is
implementation-dependent and can be dealt with using more efficient
formulations, as discussed in our previous work [20]. All the experiments were
performed on an _Intel Core i7 6700K, 4.0 GHZ, 32GB RAM_ running _Ubuntu-18.04
LTS_. Though Z3 allows for optimization of the objective function, the size of
the problems evaluated is such that no optimum is expected to be found in any
reasonable time. Therefore Z3 is set to find satisfiable (hence sub-optimal)
solutions 333The implementation of the model presented in Section II, as well
as the problem instances are available at
https://github.com/sabinoroselli/VRP.git..

In Table I the results of the computational analysis are summarized. For each
combination of the values introduced in the previous paragraph, it is reported
the number of solved problems (out of _five_), the average model generation
time, and the average solving time, both for the instances that turned out to
be satisfiable and for the ones that did not; the amount of satisfied and
unsatisfied instances is also reported. As expected, all the parameters are
directly proportional to the model size, since the model generation time
increases as they increase and to the solving time. The solving time is also
directly affected by the increase of the parameters. In general, given
problems of the same class, proving unsatisfiability seems to be faster than
finding a feasible solution; there are exceptions, especially for values of
edge reduction of 25 and 50. Even when the problem is satisfiable, a smaller
number of edges makes the problem harder, since the number of possible
conflict-free routes decreases. For _N-V-J_ equal to 25-4-7 and 36-6-8 and a
time horizon of 15 time units all problems turn out to be unsatisfiable; this
result does not come unexpected, since the graph and the gap between number of
vehicles and number of jobs are too large to execute all jobs and go back to
the starting point in time. We assumed that an obvious unsatisfiability could
be faster to prove, but the results seem to contradict our hypothesis.

##  IV Conclusion

In this paper we formulated a model to represent just-in-time material
handling system based on AGVs that takes into account limited operated range
and road segments capacity. We also designed a set of problem instances to
test the model, implemented for the state-of-the-art SMT solver Z3. The solver
turned out to be able to handle medium and large models; however, given their
complexity, these models corresponds to relatively small systems. The approach
presented can nonetheless be used to solve abstracted problems, i.e. a
simplification of the plant, as part of a compositional algorithm. Also, it
can be used to test the correctness of alternative approaches because, given
enough running time, it can provide the true optimum. Finally, the benchmark
set we generated can be used to compare performance of different approaches.

## References

  * [1] K Azadeh, MBM deKoster and D Roy  “Robotized Warehouse Systems: Developments and Research Opportunities”  In _ERIM Report Series Research in Management_ , 2017
  * [2] Kee-hung Lai and TC Edwin Cheng  “Just-in-time logistics”  CRC Press, 2016
  * [3] Kris Braekers, Katrien Ramaekers and Inneke Van Nieuwenhuyse  “The vehicle routing problem: State of the art classification and review”  In _Computers & Industrial Engineering_ 99 Elsevier, 2016, pp. 300–313
  * [4] Martin Desrochers, Jacques Desrosiers and Marius Solomon  “A new optimization algorithm for the vehicle routing problem with time windows”  In _Operations research_ 40.2 INFORMS, 1992, pp. 342–354
  * [5] Michael Schneider, Andreas Stenger and Dominik Goeke  “The electric vehicle-routing problem with time windows and recharging stations”  In _Transportation Science_ 48.4 INFORMS, 2014, pp. 500–520
  * [6] Nirup N Krishnamurthy, Rajan Batta and Mark H Karwan  “Developing conflict-free routes for automated guided vehicles”  In _Operations Research_ 41.6 INFORMS, 1993, pp. 1077–1090
  * [7] Olli Bräysy and Michel Gendreau  “Vehicle routing problem with time windows, Part I: Route construction and local search algorithms”  In _Transportation science_ 39.1 INFORMS, 2005, pp. 104–118
  * [8] Sarmad Riazi et al.  “Scheduling and routing of AGVs for large-scale flexible manufacturing systems”  In _2019 IEEE 15th International Conference on Automation Science and Engineering (CASE)_ , 2019, pp. 891–896  IEEE
  * [9] Barrie M Baker and MA Ayechew  “A genetic algorithm for the vehicle routing problem”  In _Computers & Operations Research_ 30.5 Elsevier, 2003, pp. 787–800
  * [10] Yue-Jiao Gong et al.  “Optimizing the vehicle routing problem with time windows: a discrete particle swarm optimization approach”  In _IEEE Transactions on Systems, Man, and Cybernetics, Part C (Applications and Reviews)_ 42.2 IEEE, 2011, pp. 254–267
  * [11] Sabino Roselli, Kristofer Bengtsson and Knut Åkesson  “SMT Solvers for Job-Shop Scheduling Problems: Models Comparison and Performance Evaluation”  In _Congress of CASE, the Portuguese Operational Research Society_ , 2017
  * [12] Clark W Barrett, Roberto Sebastiani, Sanjit A Seshia and Cesare Tinelli  “Satisfiability modulo theories”  In _Handbook of satisfiability_ 185, 2009, pp. 825–885
  * [13] Leonardo De Moura and Nikolaj Bjørner  “Satisfiability Modulo Theories: Introduction and Applications”  In _Commun. ACM_ 54.9 New York, NY, USA: ACM, 2011, pp. 69–77  DOI: 10.1145/1995376.1995394
  * [14] Gerhard Hiermann, Jakob Puchinger, Stefan Ropke and Richard F Hartl  “The electric fleet size and mix vehicle routing problem with time windows and recharging stations”  In _European Journal of Operational Research_ 252.3 Elsevier, 2016, pp. 995–1018
  * [15] Federico Rossi, Ramon Iglesias, Mahnoosh Alizadeh and Marco Pavone  “On the interaction between Autonomous Mobility-on-Demand systems and the power network: Models and coordination algorithms”  In _IEEE Transactions on Control of Network Systems_ 7.1 IEEE, 2019, pp. 384–397
  * [16] Emmanouil Thanos, Tony Wauters and Greet Vanden Berghe  “Dispatch and conflict-free routing of capacitated vehicles with storage stack allocation”  In _Journal of the Operational Research Society_ Taylor & Francis, 2019, pp. 1–14
  * [17] Keisuke Murakami  “Time-space network model and MILP formulation of the conflict-free routing problem of a capacitated AGV system”  In _Computers & Industrial Engineering_ 141 Elsevier, 2020, pp. 106270
  * [18] Leonardo De Moura and Nikolaj Bjørner  “Z3: An efficient SMT solver”  In _International conference on Tools and Algorithms for the Construction and Analysis of Systems_ , 2008, pp. 337–340  Springer
  * [19] Carsten Sinz  “Towards an optimal CNF encoding of boolean cardinality constraints”  In _International conference on principles and practice of constraint programming_ , 2005, pp. 827–831  Springer
  * [20] Sabino Francesco Roselli, Kristofer Bengtsson and Knut Åkesson  “Compact Representation of Time-Index Job Shop Problems Using a Bit-Vector Formulation”  In _2020 IEEE 16th International Conference on Automation Science and Engineering (CASE)_ , 2020, pp. 1590–1595  IEEE