High-Confidence Attack Detection via Wasserstein-Metric Computations

Cyber-Physical Systems (CPS) are physical processes that are tightly integrated with computation and communication systems for monitoring and control. Examples include critical infrastructure, such as transportation networks, energy distribution systems, and the Internet. These systems are usually complex, large-scale and insufficiently supervised, making them vulnerable to attacks [1, 2]. A significant literature has studied various denial of service [3], false data-injection [4, 5], replay [6, 7], sensor, and integrity attacks [8, 9, 10, 11]. The majority of these works study attack-detection problems in a control-theoretical framework. This approach essentially employs detectors to identify abnormal behaviors by comparing estimation and measurements under some predefined metrics. However, attacks could be stealthy, and exploit knowledge of the system structure, uncertainty and noise information to inflict significant damage on the physical system while avoiding detection. This motivates the characterization of the impact of stealthy attacks via e.g. reachability set analysis [9, 12, 13]. To ensure computational tractability, these works assume either Gaussian or bounded system noise. However, these assumptions fall short in modeling all natural disturbances that can affect a system. Such systems would be vulnerable to stealthy attacks that disguise themselves via an intentionally selected, unbounded and non-Gaussian distribution. When designing detectors, an added difficulty is in obtaining tractable computations that can handle these more general distributions. More recently, novel measure of concentration has opened the way for online tractable and robust attack detection with probability guarantees under uncertainty. A first attempt in this direction is [14], where a Chebyshev inequality is used to design a detector and, an assessment of the impact of stealthy attacks is given, under the assumption that system noises were bounded. With the aim of obtaining a less conservative detection mechanism, we leverage an alternative measure-concentration result via Wasserstein metric. This metric is built from data gathered on the system, and can provide concentration result significantly sharper than that via Chebyshev inequality. In particular, we address the following question for linear CPSs: How to design an online attack-detection mechanism that is robust to light-tailed distributions of system noise while remaining sensitive to attacks and limiting the impact of the stealthy attack?

To answer the question, we consider a sensor-attack detection problem on a linear dynamical system. The linear system models a remotely-observed process that is subject to an additive noise described by an unknown, not-necessarily bounded, light-tailed distribution. To identify abnormal behavior, we employ a steady-state Kalman filter as well as a benchmark distribution of an offline sequence corresponding to the normal system operation. In this framework, we propose a novel detection mechanism that achieves online and robust attack detection of stealthy attacks in high confidence.

Statement of Contributions: 1) We propose a novel detection measure, which employs the Wasserstein distance between the benchmark and a distribution of the residual sequence obtained online. 2) We propose a novel threshold-detection mechanism, which exploits measure-of-concentration results to guarantee the robust detection of an attack with high confidence using a finite set of data, and which further enables the robust tuning of the false alarm rate. The proposed detector can effectively identify real-time attacks when its behavior differs from that of the system noise. In addition, the detector can handle systems noises that are not necessarily distributed as Gaussian. 3) We propose a quantifiable, probabilistic state-reachable set, which reveals the impact of the stealthy attacker and system noise in high probability. 4) To implement and analyze the proposed mechanism, we formulate a linear optimization problem and a semidefinite problem for the computation of the detection measure as well as the reachable set, respectively. We illustrate our methods in a two-dimensional linear system with irregular noise distributions and stealthy sensor attacks.

This section presents cyber-physical system (CPS) model, and underlying assumptions on the system and attacks.

A remotely-observed, cyber-physical system subject to sensor-measurement attacks, as in Fig. 1, is described as a discrete-time, stochastic, linear, and time-invariant system

	$\displaystyle\boldsymbol{x}(t+1)$	$\displaystyle=\;A\boldsymbol{x}(t)+B\boldsymbol{u}(t)+\boldsymbol{w}(t),$		(1)
	$\displaystyle\boldsymbol{y}(t)$	$\displaystyle=\;C\boldsymbol{x}(t)+\boldsymbol{v}(t)+\boldsymbol{\gamma}(t),$

where $\boldsymbol{x}(t)\in^{n}$, $\boldsymbol{u}(t)\in^{m}$ and $\boldsymbol{y}(t)\in^{p}$ denote the system state, input and output at time $t\in\mathbb{N}$, respectively. The state matrix $A$, input matrix $B$ and output matrix $C$ are assumed to be known in advance. In particular, we assume that the pair $(A,B)$ is stabilizable, and $(A,C)$ is detectable. The process noise $\boldsymbol{w}(t)\in^{n}$ and output noise $\boldsymbol{v}(t)\in^{p}$ are independent zero-mean random vectors. We assume that each $\boldsymbol{w}(t)$ and $\boldsymbol{v}(t)$ are independent and identically distributed (i.i.d.) over time. We denote their (unknown, not-necessarily equal) distributions by $\mathbb{P}_{\boldsymbol{w}}$ and $\mathbb{P}_{\boldsymbol{v}}$, respectively. In addition, we assume that $\mathbb{P}_{\boldsymbol{w}}$ and $\mathbb{P}_{\boldsymbol{v}}$ are light-tailed¹¹1 For a random vector $\boldsymbol{w}$ such that $\boldsymbol{w}\sim\mathbb{P}$, we say $\mathbb{P}$ is $q$-light-tailed, $q=1,2,\ldots$, if $c:=\mathbb{E}_{\mathbb{P}}[\exp(b\|\boldsymbol{w}\|^{a})]<\infty$ for some $a>q$ and $b>0$. All examples listed have a moment generating function, so their exponential moment can be constructed for at least $q=1$., excluding scenarios of systems operating under extreme events, or subject to large delays. In fact, Gaussian, Sub-Gaussian, Exponential distributions, and any distribution with a compact support set are admissible. This class of distributions is sufficient to characterize the uncertainty or noise of many practical problems.

An additive sensor-measurement attack is implemented via $\boldsymbol{\gamma}(t)\in^{p}$ in (1), on which we assume the following {assumption}[Attack model] It holds that 1) $\boldsymbol{\gamma}(t)=\boldsymbol{0}$ whenever there is no attack; 2) the attacker can modulate any component of $\boldsymbol{\gamma}(t)$ at any time; 3) the attacker has unlimited computational resources and access to system information, e.g., $A$, $B$, $C$, $\boldsymbol{u}$, $\mathbb{P}_{\boldsymbol{w}}$ and $\mathbb{P}_{\boldsymbol{v}}$ to decide on $\boldsymbol{\gamma}(t)$, $t\in\mathbb{N}$.

This section presents our online detection procedure, and a threshold-based detector with high-confidence performance guarantees. Then, we propose a tractable computation of the detection measure used for online detection. We finish the section by introducing a class of stealthy attacks.

The distribution $\mathbb{P}_{\boldsymbol{r},\operatorname{D}}$ uses a small number $T$ of samples to ensure the online computational tractability of $z(t)$, so $\mathbb{P}_{\boldsymbol{r},\operatorname{D}}$ is highly dependent on instantaneous samples. Thus, $\mathbb{P}_{\boldsymbol{r},\operatorname{D}}$ may significantly deviate from the true $\mathbb{P}_{\boldsymbol{r}}$, and from $\mathbb{P}_{\boldsymbol{r},\operatorname{B}}$. Therefore, even if there is no attack, the attack detector is expected to generate false alarms due to the system noise as well as an improper selection of the threshold $\alpha$. In the following, we discuss how to select an $\alpha$ that is robust to the system noise and which results in a desired false alarm rate. Note that the value $\alpha$ should be small to be able to distinguish attacks from noise, as discussed later.

Due to space limit, please see ArXiv version [19] for proofs.

Lemma III.1 ensures that the false alarm rate is no higher than $\Delta$ when there is no attack, i.e.,

$$\operatorname{Prob}(\operatorname{Alarm}(t)=1\;|\;\operatorname{no\;attack})\leq\Delta,\quad\forall\;t.$$

Note that the rate $\Delta$ can be selected by properly choosing the threshold $\alpha$. Intuitively, if we fix all the other parameters, then the smaller the rate $\Delta$, the larger the threshold $\alpha$. Also, large values of $N$, $T$, $1-\beta$ contribute to small $\alpha$.

The Wasserstein distance $d_{W,q}(\mathbb{P}_{\boldsymbol{r},\operatorname{B}},\mathbb{P}_{\boldsymbol{r},\operatorname{D}})$, originally solving the Kantorovich optimal transport problem [18], can be interpreted as the minimal work needed to move a mass of points described via a probability distribution $\mathbb{P}_{\boldsymbol{r},\operatorname{B}}(\boldsymbol{r})$, on the space $\mathcal{Z}\subset^{p}$, to a mass of points described by the probability distribution $\mathbb{P}_{\boldsymbol{r},\operatorname{D}}(\boldsymbol{r})$ on the same space, with some transportation cost $\ell$. The minimization that defines $d_{W,q}$ is taken over the space of all the joint distributions $\Pi$ on $\mathcal{Z}\times\mathcal{Z}$ whose marginals are $\mathbb{P}_{\boldsymbol{r},\operatorname{B}}$ and $\mathbb{P}_{\boldsymbol{r},\operatorname{D}}$, respectively.

Assuming that both $\mathbb{P}_{\boldsymbol{r},\operatorname{B}}$ and $\mathbb{P}_{\boldsymbol{r},\operatorname{D}}$ are discrete, we can equivalently characterize the joint distribution $\Pi$ as a discrete mass optimal transportation plan [18]. To do this, let us consider two sets $\mathcal{N}:=\{{1},\dots,{N}\}$ and $\mathcal{T}:=\{{0},\dots,{T-1}\}$. Then, $\Pi$ can be parameterized (by $\lambda$) as follows

	$\displaystyle\Pi_{\lambda}$	$\displaystyle(\xi_{1},\xi_{2}):=\sum_{i\in\mathcal{N}}\sum_{j\in\mathcal{T}}\lambda_{ij}\delta_{\{\boldsymbol{r}^{(i)}\}}(\xi_{1})\delta_{\{\boldsymbol{r}(t-j)\}}(\xi_{2}),$
	$\displaystyle\operatorname{s.t.}$	$\displaystyle\hskip 0.0pt\sum_{i\in\mathcal{N}}\lambda_{ij}=\frac{1}{T},\;\forall\;j\in\mathcal{T},\quad\sum_{j\in\mathcal{T}}\lambda_{ij}=\frac{1}{N},\;\forall\;i\in\mathcal{N},$		(6)
		$\displaystyle\hskip 8.5359pt\lambda_{ij}\geq 0,\;\forall\;i\in\mathcal{N},\;j\in\mathcal{T}.$		(7)

Note that this characterizes all the joint distributions with marginals $\mathbb{P}_{\boldsymbol{r},\operatorname{B}}$ and $\mathbb{P}_{\boldsymbol{r},\operatorname{D}}$, where $\lambda$ is the allocation of the mass from $\mathbb{P}_{\boldsymbol{r},\operatorname{B}}$ to $\mathbb{P}_{\boldsymbol{r},\operatorname{D}}$. Then, the proposed detection measure $z(t)$ in (4) reduces to the following

	$\displaystyle(z(t))^{q}=\min\limits_{\lambda}$	$\displaystyle\;\sum_{i\in\mathcal{N}}\sum_{j\in\mathcal{T}}\lambda_{ij}\|\boldsymbol{r}^{(i)}-\boldsymbol{r}(t-j)\|^{q},$		(P)
	$\displaystyle\operatorname{s.t.}$	$\displaystyle\;\eqref{eq:lambda1},\;\eqref{eq:lambda2},$

which is a linear program over a compact polyhedral set. Therefore, the solution exists and (P) can be solved to global optimal in polynomial time by e.g., a CPLEX solver.

In this section, we address the second question in Problem 2 via reachable-set analysis. In particular, we achieve an outer-approximation of the finite-step probabilistic reachable set, quantifying the effect of the stealthy attacks and the system noise in probability.

Consider an attack sequence $\boldsymbol{\gamma}(t)$ as in (8), where $\bar{\gamma}(t)\in^{p}$ is any vector such that the attack remains stealthy. That is, $\bar{\gamma}(t)$ results in the detected distribution $\mathbb{P}_{\boldsymbol{r},\operatorname{D}}$, which is close to $\mathbb{P}_{\boldsymbol{r},\operatorname{B}}$ as prescribed by $\alpha$. This results in the following representation of the system (2)

$$\small\boldsymbol{\xi}(t+1)=\underbrace{\begin{bmatrix}A+BK&-BK\\ 0&A\end{bmatrix}}_{H}\boldsymbol{\xi}(t)+\underbrace{\begin{bmatrix}I&0\\ I&-L\end{bmatrix}}_{G}\begin{bmatrix}\boldsymbol{w}(t)\\ \bar{\gamma}(t)\end{bmatrix}.$$

(9)

We provide in the following remark an intuition of how restrictive the stealthy attacker’s action $\bar{\gamma}(t)$ has to be.

To quantify the reachable set of the system under attacks, prior information on the process noise $\boldsymbol{w}(t)$ is needed. To characterize $\boldsymbol{w}(t)$, let us assume that, similarly to the benchmark $\mathbb{P}_{\boldsymbol{r},\operatorname{B}}$, we are able to construct a noise benchmark distribution, denoted by $\mathbb{P}_{\boldsymbol{w},\operatorname{B}}$. As before,

$${\operatorname{Prob}}\left(d_{W,q}(\mathbb{P}_{\boldsymbol{w},\operatorname{B}},\mathbb{P}_{\boldsymbol{w}})\leq\epsilon_{\boldsymbol{w},\operatorname{B}}\right)\geq 1-\beta,$$

for some $\epsilon_{\boldsymbol{w},\operatorname{B}}$. Given certain time, we are interested in where, with high probability, the state of the system can evolve from some $\boldsymbol{\xi}_{0}$. To do this, we consider the $M$-step probabilistic reachable set defined as follows

$\displaystyle\mathcal{R}_{\boldsymbol{x},M}=\left\{\hskip-0.77498pt\boldsymbol{x}(M)\in^{n}\;\rule[-11.38092pt]{0.56917pt}{28.45274pt}\;\begin{array}[]{l}\textrm{system }\eqref{eq:attack}\textrm{ with }\boldsymbol{\xi}(0)=\boldsymbol{\xi}_{0},\\ \exists\;\mathbb{P}_{\boldsymbol{w}}\ni d_{W,q}(\mathbb{P}_{\boldsymbol{w}},\mathbb{P}_{\boldsymbol{w},\operatorname{B}})\leq\epsilon_{\boldsymbol{w},\operatorname{B}},\\ \exists\;\mathbb{P}_{\bar{\gamma}}\ni d_{W,q}(\mathbb{P}_{\bar{\gamma}},\mathbb{P}_{\boldsymbol{r},\operatorname{B}})\leq\alpha,\\ \end{array}\right\},$

then the true system state $\boldsymbol{x}(t)$ at time $M$, $\boldsymbol{x}(M)$, satisfies

$$\operatorname{Prob}\left(\boldsymbol{x}(M)\in\mathcal{R}_{\boldsymbol{x},M}\right)\geq 1-\beta,$$

where $1-\beta$ accounts for the independent restriction of the unknown distributions $\mathbb{P}_{\boldsymbol{w}}$ to be “close” to its benchmark.

The exact computation of $\mathcal{R}_{\boldsymbol{x},M}$ is intractable due to the unbounded support of the unknown distributions $\mathbb{P}_{\boldsymbol{w}}$ and $\mathbb{P}_{\bar{\gamma}}$, even if they are close to their benchmark. To ensure a tractable approximation, we follow a two-step procedure. First, we equivalently characterize the constraints on $\mathbb{P}$ by its probabilistic support set. Then, we outer-approximate the probabilistic support by ellipsoids, and then the reachable set by an ellipsoidal bound.

Let us now consider any transport map $f$ of $d_{M,q}$, and define a partition of the support of $\mathbb{P}_{\bar{\gamma}}$ by

$$W_{i}:=\{\boldsymbol{r}\in\mathcal{Z}\;|\;f(\boldsymbol{r})=\boldsymbol{r}^{(i)}\},\;i\in\mathcal{N},$$

where $\boldsymbol{r}^{(i)}$ are samples in $\Xi_{\operatorname{B}}$, which generate $\mathbb{P}_{\boldsymbol{r},\operatorname{B}}$. Then, we equivalently rewrite the value of the objective function defined in $d_{M,q}$, as

$\displaystyle\mathcal{H}(\mathbb{P}_{\bar{\gamma}},W_{1},\ldots,W_{N})=\sum\limits_{i=1}^{N}\int_{W_{i}}\ell^{q}(\xi,\boldsymbol{r}^{(i)})\mathbb{P}_{\bar{\gamma}}(\xi)d\xi,$

$\displaystyle\small\hskip 28.45274pt\operatorname{s.t.}\;\int_{W_{i}}\mathbb{P}_{\bar{\gamma}}(\xi)d\xi=\frac{1}{N},\;\forall\;i\in\mathcal{N},$

(10)

where the $i^{\textup{th}}$ constraints come from the fact that a transport map $f$ should lead to consistent calculation of set volumes under $\mathbb{P}_{\boldsymbol{r},\operatorname{B}}$ and $\mathbb{P}_{\bar{\gamma}}$, respectively. This results in the following equivalent problem to $d_{M,q}$ as

	$\displaystyle(d_{M,q}(\mathbb{P}_{\bar{\gamma}},\mathbb{P}_{\boldsymbol{r},\operatorname{B}}))^{q}=\min\limits_{W_{i},i\in\mathcal{N}}$	$\displaystyle\;\mathcal{H}(\mathbb{P}_{\bar{\gamma}},W_{1},\ldots,W_{N}),$		(P1)
	$\displaystyle\operatorname{s.t.}$	$\displaystyle\;\eqref{eq:mass}.$

Given the distribution $\mathbb{P}_{\bar{\gamma}}$, Problem (P1) reduces to a load-balancing problem as in [22]. Let us define the Generalized Voronoi Partition (GVP) of $\mathcal{Z}$ associated to the sample set $\Xi_{\operatorname{B}}$ and weight $\omega\in^{N}$, for all $i\in\mathcal{N}$, as

		$\displaystyle\mathcal{V}_{i}(\Xi_{\operatorname{B}},\omega)=$
		$\displaystyle\hskip 22.76228pt\{\xi\in\mathcal{Z}\;|\;\|\xi-\boldsymbol{r}^{(i)}\|^{q}-\omega_{i}\leq\|\xi-\boldsymbol{r}^{(j)}\|^{q}-\omega_{j},\;\forall j\in\mathcal{N}\}.$

It has been established that the optimal Partition of (P1) is the GVP [22, Proposition V.1]. Further, the standard Voronoi partition, i.e., the GVP with equal weights $\bar{\omega}:=\boldsymbol{0}$, results in a lower boundof (P1), when constraints (10) are removed [21], and therefore that of $d_{M,q}$. We denote this lower bound as $L(\mathbb{P}_{\bar{\gamma}},\mathcal{V}(\Xi_{\operatorname{B}}))$, and use this to quantify a probabilistic support of $\mathbb{P}_{\bar{\gamma}}$. Let us consider the support set

$$\Omega(\Xi_{\operatorname{B}},\epsilon):=\cup_{i\in\mathcal{N}}\left(\mathcal{V}_{i}(\Xi_{\operatorname{B}})\cap\mathbb{B}_{\epsilon}(\boldsymbol{r}^{(i)})\right),$$

where $\mathbb{B}_{\epsilon}(\boldsymbol{r}^{(i)}):=\{\boldsymbol{r}\in^{p}\;|\;\|\boldsymbol{r}-\boldsymbol{r}^{(i)}\|\leq\epsilon\}$. Then we have the following lemma.

In this way, the support $\Omega(\Xi_{\operatorname{B}},2\alpha)$ contains at least $1-1/2^{q}$ of the mass of all the distributions $\mathbb{P}_{\bar{\gamma}}$ such that $d_{W,q}(\mathbb{P}_{\bar{\gamma}},\mathbb{P}_{\boldsymbol{r},\operatorname{B}})\leq\alpha$. Equivalently, we have

$$\operatorname{Prob}\left(\bar{\gamma}\in\Omega(\Xi_{\operatorname{B}},2\alpha)\right)\geq 1-1/2^{q},$$

where the random variable $\bar{\gamma}$ has a distribution $\mathbb{P}_{\bar{\gamma}}$ such that $d_{W,q}(\mathbb{P}_{\bar{\gamma}},\mathbb{P}_{\boldsymbol{r},\operatorname{B}})\leq\alpha$. We characterize $\mathbb{P}_{\boldsymbol{w}}$ similarly.

Note that in practice, one can choose ball radius factor $s$ large in order to generate support which contains higher portion of the mass of unknown $\mathbb{P}$. However, it comes at a cost on the approximation of the reachable set.

A tight reachable set bound can be now derived by solving

	$\displaystyle\min\limits_{Q,a_{1},a_{2}}$	$\displaystyle\;-\log\det(Q),$		(P2)
	$\displaystyle\operatorname{s.t.}$	$\displaystyle\;\eqref{eq:Q},\;\eqref{eq:W},$

which is a convex semidefinite program, solvable via e.g., SeDuMi [23]. Note that the probabilistic reachable set is

$$\mathcal{R}_{\boldsymbol{x}}:=\cup_{M=1}^{\infty}\mathcal{R}_{\boldsymbol{x},M},$$

which again can be approximated via $Q^{\star}$ solving (P2) for⁹⁹9The set $\mathcal{R}_{\boldsymbol{x}}$ is in fact contained in the projection of $\mathcal{E}(Q^{\star})$ onto the state subspace, i.e., $\mathcal{R}_{\boldsymbol{x}}\subset\{\boldsymbol{x}\;|\;{\boldsymbol{x}}^{\top}\big{(}Q_{xx}-Q_{xe}Q_{ee}^{-1}{Q_{xe}}^{\top}\big{)}\boldsymbol{x}\leq\frac{(2-a)}{1-a}\}$ with $Q^{\star}:=\begin{bmatrix}Q_{xx}&Q_{xe}\\ {Q_{xe}}^{\top}&Q_{ee}\end{bmatrix}$. See, e.g., [13] for details.

$$\mathcal{R}_{\boldsymbol{x}}\subset\mathcal{E}(Q^{\star})=\{\boldsymbol{x}\in^{n}\;|\;{\boldsymbol{\xi}}^{\top}Q^{\star}\boldsymbol{\xi}\leq\frac{(2-a)}{1-a}\}.$$