


\section{Introduction}
Cyber-Physical System (CPS) and Internet of Things (IoT) are playing an increasingly important role in critical infrastructures and everyday life. 
The related area of CPS and IoT continues to emerge and expand as costs drop and the confluence of sensors, platforms and networks increases\cite{2018DHS_report}.
Simultaneously, the cyber-security risks and attack surfaces are also increasing \cite{cardenas2008research} since CPS relies on remote sensing devices, communication channels, and spatially distributed processors, which are prone to failures under cyber attacks on the data acquisition and communication channels.
The unintentional faults or malicious attacks could cause severe system damage, economic loss, and environmental degradation, e.g., the Stuxnet launched on Iran’s nuclear facilities\cite{STUXNET}, power blackouts in Brazil \cite{samba_stop} and Ukraine \cite{Ukraine_Blackout}, etc. The research community has recognized the importance of CPS security, especially the design of secure detection, estimation, and control strategy\cite{cardenas2009challenges}. 


Recently, substantial research efforts have been devoted to secure state estimation against malicious sensors.
One of the research paths is the finite horizon approach, where only measurements in a finite time-window are considered when recovering the system state. 
The problem is combinatorial in nature since searching for the corrupted sensors involves minimizing the $\ell_0$ norm \cite{FawziTAC2014}, which is an NP problem.
Fawzi et al.~\cite{FawziTAC2014} propose an algorithm using its convex relaxation, i.e., $\ell_1$ optimization, to solve the problem in the absence of noise. 
Similarly, Shoukry and Tabuada \cite{ShoukryTAC2016} adopt a $2$-norm batch optimization approach to solve the state estimation problem.
Shoukry et al.~\cite{Shoukry2017} propose an algorithm in the presence of noise by searching for reliable sensors using consistency check, and the searching complexity is reduced based on Satisfiability Modulo Theory.
However, in these works, the sensory data out of the time window are discarded, which may cause performance degradation and estimation delay.


Another solution is the switch estimator\cite{Mishra2017TCNS}\cite{yorie}\cite{luTAC2019} where multiple estimates are maintained based on measurements from different subsets of sensors, and the system operator switches between these estimates based on the evaluation of their reliability by consistency checking or malicious detection algorithms. However, the combinatorial nature of estimation candidates or sensor subsets may incur heavy computation and storage burden of the devices.
In view of these problems, Liu et al. \cite{liuxinghua-TAC2020} design local estimators whose weighted sum coincides with the Kalman estimate with a certain probability in the absence of attack. The local estimates are fused securely by a quadratic programming problem with an $\ell_1$ term to handle the sparse outliers. 


Even though there are efficient algorithms computing the estimation by introducing $\ell_1$ relaxation, the design of the secure state estimation is still computational challenging. For secure state recovery problem in the presence of $p$ malicious sensors, it is required that the system needs to be $2p$-sparse observable\cite{ShoukryTAC2016}, and calculating the sparse observability index for general systems is proved to be NP-hard\cite{mao2021computational}. Besides sparse observability, other results \cite{FawziTAC2014}\cite{liuxinghua-TAC2020}\cite{sandberg_TAC2014} impose stronger conditions for system resiliency, whose validations are also NP problems.
However, it has been observed by Mao et al. \cite{mao2021computational} that when the eigenvalues of the system matrix $A$ have unitary multiplicity ($A$ is non-derogatory), the computational complexity of secure state reconstruction without noise can be significantly reduced.
Leveraging upon this assumption, we propose our secure estimation scheme in the formulation of LTI system with Gaussian noise, and it has the following merits:
\begin{itemize}[left=0pt]
	\item In the presence of $p$ compromised sensors, the proposed estimation is secure if the system is $2p$-sparse observable. In the absence of attack, the proposed estimation coincides with Kalman estimation with certain probability, which can be adjusted.
	\item During the designing phase, the sparse observability index can be computed with low complexity. During the algorithm operating phase, the proposed estimation is formulated as the solution of a convex optimization problem which can be computed efficiently.
\end{itemize}

\textit{Organization:} We introduce the problem formulation and preliminary results in Section \ref{sec:problem}. The main results are provided in Section \ref{sec:main_result} and collaborated by numerical simulation in Section \ref{sec:sim}. Section \ref{sec:conclusion} finally concludes the paper.

\textit{Notations:}
Cardinality of a set ${\mathcal{S}}$ is denoted as $|{\mathcal{S}}|$. Conjugate transpose of a matrix $A$ is denoted as $A'$.
The determinant of a matrix is represented by $\det(\cdot)$. 
Diagonal matrix with diagonal elements $A_1,\cdots,A_k$ is denoted as $\text{diag}(A_1,\cdots,A_k)$.
Denote the span of row vectors of matrix $A$ as $\rs(A)$.
All-one vector with size $m\times 1$ is denoted as $\mathbf{1}_{m}$. 
${\mathbb{R}}^n$ is the space of real vectors with $n$ entries. 
The set of complex matrices with $m$ rows and $n$ columns is denoted as ${\mathbb{C}}^{m\times n}$.

\section{Problem Formulation and Preliminary Results}\label{sec:problem}	
\subsection{Secure dynamic state estimation}
In this paper, we consider the linear time-invariant system with Gaussian noise:
\begin{align}
x(k+1)&=A x(k)+w(k) , \label{eq:system} \\
y(k)&=C x(k)+v(k)+a(k) ,\label{eq:y_i_def}
\end{align}
where $x(k) \in \mathbb{R}^{n}$ is the system state, $w(k) \sim {N}(0, Q)$ and $v(k) \sim {N}(0, R)$ are i.i.d. Gaussian process noise and measurement noise with zero mean and covariance matrix $Q$ and $R$.  
The vector $y(k)\in \mathbb{R}^{m}$ is the collection of measurements from all $m$ sensors, and $i$-th entry $y_i(k)$ is the measurement from sensor $i$.
The vector $a(k)$ denotes the bias injected by an adversary and $a_i(k)$ is the attack on sensor $i$. Define $$z(k)=C x(k)+v(k)$$ as the true measurements without the attack.
The initial state $x(0) \sim {N}(0, \Sigma)$ is assumed to be zero mean Gaussian random variable and is independent from the process noise $\{w(k)\}$.

The secure dynamic estimation problem aims at recovering system state $x(k)$ at every time $k$ based on all history observations $\{y(t),0\leq t\leq k\}$ which have been partly manipulated by the malicious attacker.
It is conventional in the literature \cite{FawziTAC2014}\cite{Shoukry2017} that the attacker can only compromise a fixed subset of sensors with known maximum cardinality. 
Denote the index set of all sensors as ${\mathcal{R}} \triangleq\{1,2, \ldots, m\}$. 
Define the support of vector $a\in{\mathbb{R}}^{n}$ as $\supp(a)\triangleq \left\{i| 1\leq i\leq n , a_i\neq0 \right\}$ where $a_i$ is the $i$-th entry of vector $a$.
We have the following assumptions on the malicious adversary. 

\begin{definition}[Sparse Attack]\label{def:attack}
	The attack is called a $p$-sparse attack if the vector sequence $a(k)$ satisfies that,
	there exists a time invariant index set ${\mathcal{I}}\subseteq {\mathcal{R}} $ with $|{\mathcal{I}}| = p$ such that $\bigcup_{k=1}^{\infty} \supp\left\{a(k)\right\} = {\mathcal{I}}$.
\end{definition}

Closely related to the sparse attack, we introduce the notion of sparse observability that characterizes the system observability in the presence of attack.

\begin{definition}[Sparse observability]\label{df:sparse_obs}
	The sparse observability index of system \eqref{eq:system}-\eqref{eq:y_i_def} is the largest integer $s$ such that system\footnote{The matrix $C_{{\mathcal{R}}\setminus{\mathcal{I}}}$ represents the matrix composed of rows of $C$ with row index in ${\mathcal{R}}\setminus{\mathcal{I}}$.} $(A,C_{{\mathcal{R}}\setminus{\mathcal{I}}})$ is observable for any set of sensors ${\mathcal{I}}\subset{\mathcal{R}}$ with cardinality $|{\mathcal{I}}| = s$. When the sparse observability index is $s$, we say that the system with pair $(A,C)$ is $s$-sparse observable.
\end{definition}
Define $y(k_1:k_2)$ as the sequence $\{y(k_1),y(k_1+1),\cdots,y(k_2)\}$. Similar notation is also applied on $z(k)$.	
For linear Gaussian noise system, the estimation is secure if the estimation error is bounded by a constant term irrelevant to the attack. 

\begin{definition}[Secure estimator]\label{def:resi}
	An estimator is an infinite sequence of mappings $g=\{g_k\}_{k=1}^{{\infty}}$ where $g_k$ is a mapping from all the history observations to an estimation at time $k$:
	$$g_k\left(y(0:k)\right)=\hat{x}(k).$$
	Define the estimation difference introduced by attack as
	$$d_k\triangleq \left\|g_k\left(z(0:k)\right)-g_k\left(y(0:k)\right) \right\|_2 .
	$$
	The estimator is said to be secure against the $p$-sparse attack if the following holds:

	$$\sup_{k\in{\mathbb{Z}}^+} {\mathbb{E}} \left[d_k^2\right] < {\infty} ,$$
	where ${\mathbb{E}}$ is the expectation with respect to the probability measure generated by the Gaussian noise $\{w(k)\}$ and $\{v(k)\}$.
\end{definition}
If all sensors are benign, i.e., $a(k)=\mathbf{0}$ for all $k$, the optimal state estimator is the classical Kalman filter:
\begin{align*}
\hat{x}(k)&=\hat{x}(k | k-1)+K(k)\left[y(k)-C \hat{x}(k  | k-1)\right] ,\\
P(k)&=P(k  | k-1)-K(k) C P(k  | k-1),
\end{align*}
where
\begin{align*}
&\hat{x}(k  | k-1)=A \hat{x}(k-1), P(k  | k-1)=A P(k-1) A{'}+Q ,\\	
&K(k)=P(k  | k-1) C{'}\left(C P(k  | k-1) C{'}+R\right)^{-1},
\end{align*}
with initial condition $\hat{x}(0  |-1)=0,\ P(0  |-1)=\Sigma $.
It is well-known that for observable system, the estimation error covariance matrices $P(k)$ and the gain $K(k)$ will converge to the following matrices $P$ and $K$:
\begin{align*}
&P \triangleq \lim _{k \rightarrow \infty} P(k),\ P_{+}=A P A{'}+Q ,\\
&K \triangleq P_{+} C{'}\left(C P_{+} C{'}+R\right)^{-1}.
\end{align*}
Since typically the control system will be running for an extended period of time, we focus on the case where the Kalman filter is in steady state, and thus the Kalman filter reduces to the following fixed-gain linear estimator:
\begin{equation}\label{eq:fix_gain_kalman}
\hat{x}(k+1)=(A-K C A) \hat{x}(k)+K y(k+1) .
\end{equation}
Before introducing our work, we first recall some results in the previous work that decomposes the fix gain Kalman filter to local estimates and recovers it by an optimization problem. 


\subsection{Preliminary Results}\label{sec:preli}
The following preliminary results in this subsection are from \cite{liuxinghua-TAC2020}.
We introduce the following assumption:
\begin{assumption}\label{as:distinct_eigvalue}
	$A-K C A$ has $n$ distinct eigenvalues. Moreover, $A-K C A$ and $A$ do not share any eigenvalue.
\end{assumption}
Since $A-K C A$ has distinct eigenvalues, it can be diagonalized as:
\begin{equation}\label{eq:VLambda}
A-K C A=V \Pi V^{-1},
\end{equation}
where $\Pi$ is a diagonal matrix with the eigenvalues of $A-KCA$ as diagonal entries. Denote these $n$ eigenvalues (diagonal entries) as $\pi_{1},\cdots,\pi_{n}$. 
We design local estimation $\zeta_{i}(k)$ with the following dynamic:
\begin{equation}\label{eq:def_zeta}
\zeta_{i}(k+1)=\Pi \zeta_{i}(k)+\mathbf{1}_{n} y_{i}(k+1) .
\end{equation}
In other words, $\zeta_{i}(k)$ is the state of a estimator with dynamics similar as $A-KCA$ while only takes observations from sensor $i$.

Define $G_i$ as
\begin{equation}\label{eq:def_Gi}
G_{i} \triangleq\left[\begin{array}{c}
C_{i} A\left(A-\pi_{1} I\right)^{-1} \\
\vdots \\
C_{i} A\left(A-\pi_{n} I\right)^{-1}
\end{array}\right].
\end{equation}
In the following, we claim that local estimate $\zeta_i(k)$ is actually a stable estimate of state $x(k)$ transformed by the constant matrix $G_i$, i.e., their difference $\epsilon_i(k)\triangleq\zeta_{i}(k)-G_ix(k)$ is a stationary Gaussian process.
The following result has been proved in \cite{liuxinghua-TAC2020}. 
\begin{lemma}[\hspace{-0.0001pt}\cite{liuxinghua-TAC2020}]\label{lm:epsilon}
	Let $\epsilon_i(k)\triangleq\zeta_{i}(k)-G_ix(k)$, then 
	\begin{align}
	\epsilon_{i}(k+1)= \Pi \epsilon_{i}(k)+\left(G_{i}-\mathbf{1}_{n} C_{i}\right) w(k) \notag \\
	-\mathbf{1}_{n} v_{i}(k+1)&-\mathbf{1}_{n} a_{i}(k+1) . \label{eq:epsilon}
	\end{align}
\end{lemma}
Noticing that $\Pi$ is a strictly stable matrix (all eigenvalues are within the open unit disk) and $w(k),v(k)$ are Guaissain zero-mean noise, when the attack is absence, i.e., $a_i(k)=0$, the residue $\epsilon_i(k)$ is a stationary Gaussian process.
Define $\tilde{Q} \in {\mathbb{C}}^{m n \times m n}$ as the covariance of noise term $\left(G_{i}-\mathbf{1}_{n} C_{i}\right) w(k) -\mathbf{1}_{n} v_{i}(k+1)$ for all $i$, i.e.,
\begin{align*}
\tilde{Q} \triangleq&
\cov\left(	\begin{bmatrix}
G_1-\mathbf{1}_{n} C_1 \\
\vdots \\
G_m-\mathbf{1}_{n} C_m
\end{bmatrix} w(k)\right)+
\cov\left(	\begin{bmatrix}
\mathbf{1}_{n} v_{1}(k+1)\\
\vdots \\
\mathbf{1}_{n} v_{m}(k+1)
\end{bmatrix} \right)
\\
=&
\begin{bmatrix}
G_1-\mathbf{1}_{n} C_1 \\
\vdots \\
G_m-\mathbf{1}_{n} C_m
\end{bmatrix}
Q\begin{bmatrix}
G_1-\mathbf{1}_{n} C_1 \\
\vdots \\
G_m-\mathbf{1}_{n} C_m
\end{bmatrix}^{'}
+ R\otimes \mathbf{1}_{n\times n},
\end{align*}
where $\otimes$ is the Kronecker product.
Define $\tilde{\Pi}\in {\mathbb{C}}^{m n \times m n}$ as
$$
\tilde{\Pi} \triangleq\left[\begin{array}{ccc}
\Pi & & \\
& \ddots & \\
& & \Pi
\end{array}\right].
$$
As $k{\rightarrow}{\infty}$, the stable covariance of $\epsilon(k)\triangleq\left[\epsilon_1(k){'},\cdots,\epsilon_m(k){'}\right]{'}$ is the solution $\tilde{W}$ of the following Lyapunov equation:
$$
\tilde{W}=\tilde{\Pi} \tilde{W} \tilde{\Pi}{'}+\tilde{Q}.
$$ 
The matrix $\tilde{W}$ is well-defined since $\Pi$ is strictly stable. As a result, the secure estimation can be recovered by the solution of the following optimization problem where $\zeta(k)\triangleq\left[\zeta_1(k){'},\cdots,\zeta_m(k){'}\right]{'} $ and $G\triangleq\left[G{'}_1,\cdots,G{'}_m\right]{'} $.
\begin{subequations}\label{pb:old_lasso}
	\begin{align}
	&\underset{\check{x}(k), \mu(k), \nu(k)}{\operatorname{minimize}}\quad \frac{1}{2} \mu(k){'} \tilde{W}^{-1} \mu(k) + \gamma \|\nu(k)\|_1 \label{min:old_lasso} \\
	& \text{ subject to}\quad
	\zeta(k)=
	G \check{x}(k)+\mu(k)+\nu(k). \label{eq:old_lasso}
	\end{align}
\end{subequations}
The parameter $\gamma$ is a non-negative constant chosen by the system operator.
According to \cite{liuxinghua-TAC2020}, the solution $\check{x}(k)$ to problem (\ref{pb:old_lasso}) is a secure estimation and has the following properties.
\begin{theorem}[\hspace{-0.01pt}\cite{liuxinghua-TAC2020} Theorem 4]\label{th:TAC}
	In the presence of $p$-sparse attack, the state estimation $\check{x}(k)$ solved from problem \eqref{pb:old_lasso} is secure if the following inequality holds for all $u \neq \mathbf{0}$, $u\in{\mathbb{R}}^n$:
	\begin{equation}\label{eq:cond}
	\sum_{i \in \mathcal{I}}\left\|G_{i} u\right\|_{1}<\sum_{i \in \mathcal{I}^{c}}\left\|G_{i} u\right\|_{1}, \quad \forall\ {\mathcal{I}}\subset {\mathcal{R}}, |\mathcal{I}|\leq p .
	\end{equation}
\end{theorem}
Even though Theorem \ref{th:TAC} establishes the sufficient condition of the estimation to be secure, validating \eqref{eq:cond} is NP-hard.
In the following section, we reduce the complexity of checking condition \eqref{eq:cond} by transforming $G_i$ to its canonical form under the assumption that the geometric multiplicities of all the eigenvalues of $A$ are $1$ ($A$ is non-derogatory). 


\section{Main Results}\label{sec:main_result}
In this section, under the assumption on $A$, we first prove that the span of rows of $G_i$ coincides with the observability space of $(A,C_i)$, which implies that the matrix $G_i$ has a canonical form under row operations. Based on the canonical form, a new secure estimation scheme is proposed.  


\subsection{Canonical form of $G_i$}\label{subsec:transform}
In order to prevent degeneration problem, we introduce the following assumption.
\begin{assumption}\label{as:geo_mul_1}
	System matrix $A$ is non-singular, and all the eigenvalues of $A$ have geometric multiplicity $1$.
\end{assumption}
Without loss of generality, based on Assumption \ref{as:geo_mul_1}, we assume that $A$ is in the following Jordan canonical form:
\begin{align*}
A=\begin{bmatrix}
J_{1} & \mathbf{0} & \cdots & \mathbf{0} \\
\mathbf{0} & J_{2} & \cdots & \mathbf{0} \\
\vdots & \vdots & \ddots & \vdots \\
\mathbf{0} & \mathbf{0} & \cdots & J_{l} 
\end{bmatrix},\
J_k=
\begin{bmatrix}
\lambda_{k} & 1 & 0 &  \cdots & {0} \\
{0} & \lambda_{k} & 1 & \cdots & {0} \\
\vdots & \ddots & \ddots   & \ddots & \vdots \\
\vdots & \ddots & \ddots  & \ddots & 1 \\
{0} & \cdots & \cdots & 0 & \lambda_{k} 
\end{bmatrix},
\end{align*}
and the (complex) eigenvalues with respect to distinct blocks are different, i.e., $k_1\neq k_2$ implies $\lambda_{k_1}\neq\lambda_{k_2}$.
Define the observability matrix of system $(A,C_i)$ as 
\begin{equation}\label{eq:def_O}
O_{i} \triangleq\left[\begin{array}{c|c|c|c}
C_{i}{'} &
\left(C_{i} A\right){'} &
\cdots &
\left(C_{i} A^{n-1}\right){'}
\end{array}\right]{'}.
\end{equation}
Before continuing on, we need the following notation of state observability. 
Define ${\mathcal{S}}_j$ as the index set of sensors that can observe state $j$, i.e.
\begin{equation}\label{eq:def_Oc}
{\mathcal{S}}_j\triangleq \{i\in{\mathcal{R}}\ |\ O_i{'} e_j\neq \mathbf{0} \},
\end{equation}
where ${\mathcal{R}}\triangleq \{1,2,\cdots,m\}$ is the index set of all sensors and $e_j$ is the $n$-dimensional canonical basis vector with 1 on the $j$-th entry and 0 on the other entries.
We have the following theorem characterizing the structure of $G_i$. 
\begin{theorem}\label{th:span}
	Assume system matrix $A$ satisfies Assumption \ref{as:geo_mul_1}, then the following equation holds:
	\begin{align}
	\rs(G_i)=\rs(O_i)=\rs(H_i) ,
	\end{align}
	where $H_i$ is the following diagonal matrix
	\begin{equation*}
	H_i\triangleq \begin{bmatrix}
	{\mathbb{I}}_{i\in{\mathcal{S}}_1} & & & \\
	&{\mathbb{I}}_{i\in{\mathcal{S}}_2} & &  \\
	& & \ddots &  \\
	& & & {\mathbb{I}}_{i\in{\mathcal{S}}_n}
	\end{bmatrix} ,
	\end{equation*}
	and $\mathbb{I}_{i\in{\mathcal{S}}}$ is the indicator function that takes the value 1 when ${i\in{\mathcal{S}}}$ and value 0 when ${i\notin{\mathcal{S}}}$.
\end{theorem}

\begin{proof}
	Define the characteristic polynomial of $A$ as $p(x)=a_n x^n +\cdots+a_1 x +a_0$.
	Define polynomial fraction  $q_\pi(x)$ with respect to constant $\pi$ as
	$q_\pi(x)\triangleq\frac{p(x)-p(\pi)}{x-\pi}$ where $x\neq \pi$.
	For polynomial fraction $q_\pi(x)$, we have
	\begin{align*}
	&\frac{p(x)-p(\pi)}{x-\pi}\\
	=&\frac{1}{x-\pi}\left(a_n(x^n-\pi^n)+a_{n-1}(x^{n-1}-\pi^{n-1})+\cdots+a_1(x-\pi)\right)\\
	=&a_n\left(x^{n-1}+\pi x^{n-2}+\pi^2x^{n-3}+\cdots+\pi^{n-1}\right)\\
	&+a_{n-1}\left(x^{n-2}+\pi x^{n-3}+\pi^2x^{n-4}+\cdots+\pi^{n-2}\right)\\
	&+\cdots \\
	&+a_{2}\left(x+\pi\right)\\
	&+a_1,
	\end{align*}
	where the last inequality comes from the fact that 
	$$\frac{x^n-\pi^n}{x-\pi}=\pi^0x^{n-1}+\pi^1 x^{n-2}+\pi^2x^{n-3}+\cdots+\pi^{n-1}x^0.$$
	The above interpretation of $q_\pi(x)$ is also valid when $x$ is a square matrix.
	As a result, by rearranging the terms of $A$ with the same power, one obtains
	\begin{align*}
	q_{\pi_j}(A)=&a_n\left(A^{n-1}+\pi_j A^{n-2}+\cdots+\pi_j^{n-1}I\right)\\
	&+a_{n-1}\left(A^{n-2}+\pi_j A^{n-3}+\cdots+\pi_j^{n-2}I\right)\\
	&+\cdots+a_{2}\left(A+\pi_j I\right)+\cdots+a_1 I\\
	=&a_nA^{n-1}+\left(a_n\pi_j+a_{n-1}\right)A^{n-2}+\cdots\\
	&+\left(a_n\pi_j^{n-1}+a_{n-1}\pi_j^{n-2}+\cdots+a_2\pi_j+a_1\right)I.
	\end{align*}
	Define $b_{j,k}$ as the constant scalar coefficient of $A^k$ in $q_{\pi_j}(A)$:
	\begin{equation}\label{eq:bjk}
	b_{j,k}\triangleq\sum_{i=0}^{n-k-1} a_{i+k+1} \pi_j^i.
	\end{equation}
	Thus, $q_{\pi_j}(A)=\sum_{k=0}^{n-1}b_{j,k}A^{k}$.
	The $j$-th row of matrix $G_i$ can be reformulated as
	\begin{align*}
	&C_{i} A\left(A-\pi_{j} I\right)^{-1}\notag \\
	=&-\frac{1}{p(\pi_j)} C_i A q_{\pi_j}(A)\notag \\
	=&-\frac{1}{p(\pi_j)} C_i \left(\sum_{k=0}^{n-1}b_{j,k}A^{k+1}\right)\notag \\
	=&-\frac{1}{p(\pi_j)}\begin{bmatrix}
	b_{j,0} & b_{j,1} & \cdots  & b_{j,n-1} 
	\end{bmatrix} 
	\begin{bmatrix}
	C_i A \\
	C_i A^2 \\
	\vdots \\
	C_i A^n
	\end{bmatrix}\notag \\
	=&-\frac{1}{p(\pi_j)}\begin{bmatrix}
	b_{j,0} & b_{j,1} & \cdots  & b_{j,n-1} 
	\end{bmatrix} O_i A.
	\end{align*}
	
	
	Therefore, $G_i$ can be interpreted as follows
	\begin{align}\label{eq:Gi=DOA}
	G_i = \mathcal{D}_1\begin{bmatrix}
	b_{1,0} & b_{1,1} & \cdots  & b_{1,n-1} \\
	b_{2,0} & b_{2,1} & \cdots  & b_{2,n-1} \\
	\vdots & \vdots & \ddots  & \vdots \\
	b_{n,0} & b_{n,1} & \cdots  & b_{n,n-1} 
	\end{bmatrix}
	O_i A 
	= & \mathcal{D}_1\mathcal{D}_2\mathcal{D}_3 O_i A ,
	\end{align}
	where $\mathcal{D}_1\triangleq\text{diag}\left(-\frac{1}{p(\pi_1)},-\frac{1}{p(\pi_2)},\cdots,-\frac{1}{p(\pi_n)}\right)$,
	\begin{align*}
	\mathcal{D}_2\triangleq\hspace{-2pt}
	\begin{bmatrix}
	\pi_1^{n-1} & \pi_1^{n-2} & \cdots  & 1 \\
	\pi_2^{n-1} & \pi_2^{n-2} & \cdots  & 1 \\
	\vdots & \vdots & \cdots  & \vdots \\
	\pi_n^{n-1} & \pi_n^{n-2} & \cdots  & 1
	\end{bmatrix}, 
	\mathcal{D}_3\triangleq\hspace{-2pt}
	\begin{bmatrix}
	a_n & 0 & \cdots &   0 \\
	a_{n-1} & a_n & \cdots &   0 \\
	\vdots & \vdots & \ddots  & \vdots \\
	a_1 & a_2 & \cdots  & a_n 
	\end{bmatrix}.
	\end{align*}
	According to Assumption \ref{as:distinct_eigvalue}, all $\pi_j$ are distinct eigenvalues and they are not the eigenvalues of $A$, i.e. the diagonal matrix $\mathcal{D}_1$ and the Vandermonde matrix $\mathcal{D}_2$ are invertible. Moreover, $a_n=1$. Therefore, the lower triangular Toeplitz matrix $\mathcal{D}_3$ is invertible and thus $\rs(G_i)=\rs(O_i A)$ from equation \eqref{eq:Gi=DOA}. 		
	We continue to prove $\rs(O_i)=\rs(O_i A)$. Considering that $A^n=-a_{n-1}A^{n-1}-\cdots-a_0 I$, one obtains the following equation \eqref{eq:O_OA}.
	\begin{equation}
	\label{eq:O_OA}
	O_i A=
	\begin{bmatrix}
	0 & 1 & 0 &  \cdots & 0 \\
	0 & 0 & 1 &  \cdots & 0 \\
	\vdots & \vdots & \vdots & \ddots & \vdots \\
	0 & 0 & 0 &  \cdots & 1 \\
	-a_0 & -a_1 & -a_2 & \cdots &  -a_{n-1}
	\end{bmatrix}
	O_i .
	\end{equation}
	According to Assumption \ref{as:geo_mul_1}, $A$ is invertible and $a_0=(-1)^n\det(A)\neq 0$, which leads to the equation that $\rs(O_i)=\rs(O_i A)$.
	Since $A$ is assumed to be in the Jordan canonical form and all eigenvalues have geometric multiplicity 1, one can verify that nonzero columns of $O_i$ are linear independent, and equivalently nonzero columns of $G_i$ are linear independent. Therefore, $i\in{\mathcal{S}}_j$ is equivalent to that $j$-th column of $G_i$ is non-zero, i.e., $G_i$ has the same row-span with the canonical form $H_i$.
\end{proof}


Recall that $\rs(\cdot)$ represents the span of row vectors of a matrix.
According to Theorem \ref{th:span} and definition of $\rs(\cdot)$, one directly obtains the following corollary:
\begin{corollary}\label{co:P_i}
	For every sensor index $i\in{\mathcal{R}}$, there exists an invertible matrix $P_i$ such that $P_iG_i=H_i$.
\end{corollary}

After transformation $P_i$, matrix $G_i$ is transformed into canonical form $H_{i}$ whose rows are either canonical basis vectors or zero vectors. 
The non-zero entries of $H_i$ records the state observability of sensor $i$. Therefore, the sparse observability index can be directly obtained from diagonal matrices $H_i,i\in\{1,2,\cdots,m\}$ or, equivalently, from sets ${\mathcal{S}}_j,j\in\{1,2,\cdots,n\}$.
\begin{corollary}\label{co:sparse_obs}
	Denote $s$ as the sparse observability index of system \eqref{eq:system}-\eqref{eq:y_i_def}. Then 
	\begin{equation*}
	s=\min_{j\in\{1,2,\cdots,n\}} |{\mathcal{S}}_{j}| - 1 .
	\end{equation*}
\end{corollary} 
\begin{proof}
	According to the definition of $s$, for arbitrary $\overline{s}$ that satisfies $\overline{s}\geq s+1$, there exists a state index $j^*$ and a sensor index set ${\mathcal{I}}^*$ with $|{\mathcal{I}}^*|=\overline{s}$ such that ${\mathcal{S}}_{j^*} \cap \left({\mathcal{R}}\setminus{\mathcal{I}}^*\right)=\varnothing$.
	As a result, state $j^*$ can not be observed by any sensor in ${\mathcal{R}}\setminus {\mathcal{I}}^*$, i.e.,
	\begin{equation*}
	e_{j^*}\notin \rs(O_i),\ \forall i\in{\mathcal{R}}\setminus {\mathcal{I}}^*.
	\end{equation*}
	and thus system $(A,C_{{\mathcal{R}}\setminus{\mathcal{I}}^*})$ is not observable.
	For arbitrary $\underline{s}$ that satisfies $\underline{s}\leq s$, arbitrary $j$ and arbitrary ${\mathcal{I}}$ with $|{\mathcal{I}}|=\underline{s}$, one obtains ${\mathcal{S}}_{j} \cap \left({\mathcal{R}}\setminus{\mathcal{I}}\right)\neq\varnothing$, which means for all $j$, there exists $i^*\in{\mathcal{R}}\setminus {\mathcal{I}}$ such that: $e_{j}\in \rs(O_{i^*})$. Therefore, system $(A,C_{{\mathcal{R}}\setminus{\mathcal{I}}})$ is observable. According to Definition \ref{df:sparse_obs}, the system is $s$-sparse observable.
\end{proof}


In view of Corollary \ref{co:sparse_obs}, the system sparse observability index can be obtained with low computational complexity since calculating cardinality of ${\mathcal{S}}_j$ can be done in \emph{polynomial time} with respect to $m$ and $n$.
In the following subsection, we will propose an information fusion scheme that is secure in the presence of $p$-sparse attack as long as the system is $2p$-sparse observable.


\subsection{Secure Information Fusion}
Recalling the transformation $P_i$ introduced in Corollary \ref{co:P_i}, define $\tilde{P} \triangleq \text{diag}\left(P_1,\cdots,P_m\right),\ \tilde{M}\triangleq\tilde{P}\tilde{W}\tilde{P}{'}$ and
\begin{align}
{\mathcal{Y}} (k)\triangleq
\begin{bmatrix}
P_1\zeta_{1}(k) \\
\vdots \\
P_m\zeta_{m}(k)
\end{bmatrix}, \
H\triangleq\begin{bmatrix}
H_{1} \\
\vdots \\
H_{m}
\end{bmatrix} .	
\end{align}
We present the following optimization problem whose solution is our proposed estimation. The constant scalar $\gamma$ is an adjustable parameter.
\begin{subequations}\label{pb:lasso}
	\begin{align}
	&\underset{\tilde{x}(k), \mu(k),\nu(k)}{\operatorname{minimize}}\quad \frac{1}{2} \mu(k){'} \tilde{M}^{-1} \mu(k) + \gamma\|\nu(k)\|_1  \\
	&\ \text{subject to }\quad
	{\mathcal{Y}}(k)=H \tilde{x}(k) +\mu(k)+\nu(k). 
	\end{align}
\end{subequations}

We first prove that, in the absence of attack, the estimation  $\tilde{x}(k)$ is optimal for certain probability that can be adjusted.
We need to define the following matrices to facilitate the understanding of the result. 
Define
$$F_i\triangleq V\diag(V^{-1}K_i),\ F=\begin{bmatrix} F_1&F_2&\cdots & F_m \end{bmatrix},$$
where $V$ is defined in \eqref{eq:VLambda} and $K_i$ is the $i$-th column of Kalman gain $K$. Moreover, $\diag(V^{-1}K_i)$ is an $n\times n$ diagonal matrix with the $j$-th diagonal entry equals to $j$-th element of vector $V^{-1}K_i$.
Recalling that $\epsilon(k)\triangleq\left[\epsilon_1(k)',\cdots,\epsilon_m(k)'\right]'$ and $\epsilon_i(k)=\zeta_i(k)-G_ix(k)$ from Lemma \ref{lm:epsilon},
the following theorem characterizes the estimation $\tilde{x}(k)$ in the absence of attack. 



\begin{theorem}\label{th:no_attack}
	In the absence of attack, if the parameter $\gamma$ in problem (\ref{pb:lasso}) satisfies
	\begin{align}\label{eq:kalman_cond}
	\left\|\tilde{M}^{-1}\left(I-GF\right)\epsilon(k)\right\|_\infty \leq \gamma,
	\end{align}
	then our proposed estimation $\tilde{x}(k)$ solved from problem \eqref{pb:lasso} is equivalent to the Kalman estimation $\hat{x}(k)$ defined in (\ref{eq:fix_gain_kalman}), i.e.,
	\begin{equation}\label{eq:eq_to_kalman}
	\tilde{x}(k)=\hat{x}(k).
	\end{equation}
\end{theorem}
\begin{remark}
	Noticing that $\epsilon(k)$ is a stationary Gaussian process from Lemma \ref{lm:epsilon}, the probability that \eqref{eq:kalman_cond} holds can be calculated given the parameters $A,C,Q,R$ of the linear Gaussian system \eqref{eq:system}-\eqref{eq:y_i_def}.
	By tuning design parameter $\gamma$, the probability of recovering the Kalman estimation can be adjusted.
\end{remark}
\begin{proof}
	Consider the following least square problem by suppressing $\nu(k)=\mathbf{0}$ in problem \eqref{pb:lasso}.
	\begin{subequations}\label{pb:least_square}
		\begin{align}
		&\underset{\tilde{x}_\text{ls} (k), \varphi(k)}{\operatorname{minimize}}\quad\ \frac{1}{2} \varphi(k)' \tilde{M}^{-1} \varphi(k)   \\
		&\text{subject to }\quad
		{\mathcal{Y}}(k)=H \tilde{x}_\text{ls}(k) +\varphi(k). 
		\end{align}
	\end{subequations}
	It can be verified that the solution is equivalent to problem (30) in \cite{liuxinghua-TAC2020}.
	Thus, the solution $\varphi(k)$ of problem \eqref{pb:least_square} satisfies the following equation due to Theorem 1 in \cite{liuxinghua-TAC2020}:
	\begin{equation}\label{eq:ls_recover}
	\tilde{x}_\text{ls}(k)=\hat{x}(k),\ \varphi(k)=(I-GF)\epsilon(k),
	\end{equation}
	where $\hat{x}(k)$ is the fixed gain Kalman estimation defined in \eqref{eq:fix_gain_kalman}.
	Considering the KKT condition of problem \eqref{pb:lasso}, one obtains that if $\left\|\tilde{M}^{-1} \varphi(k)\right\|_{\infty}\leq\gamma$, then the solution $\tilde{x}(k),\mu(k),\nu(k)$ satisfy
	\begin{equation}\label{eq:lasso_recover}
	\tilde{x}(k)=\tilde{x}_\text{ls}(k)=\hat{x}(k),\ \mu(k)=\varphi(k),\ \nu(k)=\mathbf{0}.
	\end{equation}
	Combining \eqref{eq:ls_recover} and \eqref{eq:lasso_recover}, result in Theorem \ref{th:no_attack} is obtained.
\end{proof}

In the presence of $p$-sparse attack, 
we have the following theorem demonstrating that the estimation $\tilde{x}(k)$ solved from problem \eqref{pb:lasso} is secure. 

\begin{theorem}\label{th:main}
	In the presence of arbitrary $p$-sparse attack, if the system $(A,C)$ is $2p$-sparse observable, the proposed estimation $\tilde{x}(k)$ solved from (\ref{pb:lasso}) is secure.
\end{theorem}


The proof of Theorem \ref{th:main} is provided in Appendix \ref{ap:th_main} for legibility.
Since sparse observability index only requires simple computation according to Corollary \ref{co:sparse_obs}, this work reduces the complexity of evaluating system vulnerability significantly under the assumption of geometric multiplicity.
For general $A$ that has eigenvalues with geometric multiplicity larger than 1 ($A$ is derogatory), computing sparse observability index is an NP-hard problem \cite{ShoukryTAC2016}\cite{mao2021computational}\cite{yanwen_CDC19}, and there is no computational efficient solution unless P$=$NP. Simultaneously, for algorithm online operation, the computing of our proposed estimation involves solving a convex optimization problem based on LASSO \cite{LASSOTibshirani}, which can be done efficiently.

Moreover, the condition of $2p$-sparse observable is necessary for secure state recovering in the presence of $p$ compromised sensors \cite{ShoukryTAC2016}, which establishes the fundamental limit. Our proposed estimator achieves this limit and also holds optimality in the absence of attack for certain probability, which can be adjusted by parameter $\gamma$, according to Theorem \ref{th:no_attack}. 
By properly tuning parameter $\gamma$, the system operator can achieve the trade-off between the performance in normal operation and performance under attack. 




\section{Simulation}\label{sec:sim}
We use an inverted pendulum for the numerical simulation. The physical parameters are illustrated in Fig. \ref{fig:invpen}. 
The control input $u(k)$ is the force applied on the cart, and we assume that there are no frictions of any form.
The states $x_1,x_2,x_3,x_4$ represent cart position coordinate, cart velocity, pendulum angle from vertical and pendulum angle velocity respectively. 
\begin{figure}[htpb]
	\centering\vspace{-10pt}
	\input{inv_pen.tex}\vspace{-10pt}
	\caption{Illustration of the inverted pendulum.}\label{fig:invpen}
\end{figure}

Consider the system linearized at $x_3=x_4=0$, and we sample the continuous-time linear system periodically every $0.1$ seconds. The system matrix $A$ can be transformed into the following Jordan canonical form $$
\begin{bmatrix}
1  & 1 &  0 & 0\\
0  & 1 & 0 &  0\\
0  & 0 & 0.642 & 0  \\
0  & 0 & 0 &  1.557 
\end{bmatrix},
$$ 
and there is a Jordan block with size $2\times 2$ and the geometric multiplicity of all the eigenvalues are 1.


The output matrix $C$ is defined as $C_1=C_2=C_3=[1\ 0\ 0\ 0]$,  $C_4=[0\ 0\ 1\ 0]$, i.e., first three sensors are monitoring the cart position $x_1$ and the fourth sensor is monitoring pendulum angle $x_3$.
In this case, the first 3 sensors can observe all the states, and the fourth sensor can only observe state $x_3$ and $x_4$.
According to Corollary \ref{co:sparse_obs}, the sparse observability index is 2, and the system can at most withstand 1 corrupted sensor.

The noise covariances of the system satisfy $w(k) \sim {N}(0, 0.001\cdot I_4)$ and $v(k) \sim {N}(0, 0.001\cdot I_4)$, where $I_n$ is the identity matrix with size $n\times n$. 
The controller of the system is designed as a Linear-Quadratic Regulator (LQR), and the feedback matrix is chosen as 
$$K_{\rm lqr}=\begin{bmatrix}
-0.604 & -1.678 & -39.514 & -9.721
\end{bmatrix}.
$$
We demonstrate the performance of different estimators of close-loop system where $u(k)=-K_{\rm lqr} x(k)$.

\begin{figure}[htpb]
	\input{states_under_attack_closeloop.tex}
	\caption{Value of system states and their estimation in the presence of attack. The black solid line represents the real state. The red dashed line and blue dotted line represent Kalman estimation and our proposed secure estimation respectively. The initial state is $x(0)=[0,\ 1,\ 0,\ 1]{'}$.} \label{fig:close_loop}

\end{figure}

Fig. \ref{fig:close_loop} presents the Kalman estimation and our proposed secure estimation in the presence of attack on sensor 4. The attack $a_4(k)$ is a time-independent random value uniformly distributed on open interval $\left(-\frac{\pi}{2},\frac{\pi}{2}\right)$. In other words, the measurement of the pendulum angle is randomly shifted by the malicious adversary in magnitude of $\pm 90 ^{\circ}$. The attack is launched since the beginning time. As shown in Fig \ref{fig:close_loop}, Kalman estimation (denoted as red dashed line) has large estimation error under the attack while our estimation (denoted as blue dotted line) converges to the real state quickly.


\begin{figure}[htpb]
	\centering
	\input{MSE.tex}
	\caption{Estimation mean square error (MSE) in the absence and in the presence of attack with varying tuning parameter $\gamma$.} \label{fig:MSE}

\end{figure}

Fig. \ref{fig:MSE} illustrates the estimation mean square error (MSE) of our proposed estimator with varying $\gamma$ in the absence and in the presence of attack on sensor 4. 
The attack launched on the system is the same as in Fig. \ref{fig:close_loop}.
The MSE of the oracle Kalman estimation, .i.e., Kalman estimation not affected by the attack, is illustrated by the red dashed line in Fig. \ref{fig:MSE}.
As shown in the figure, by properly choosing $\gamma$, the MSE of our proposed estimator is smaller than that of Kalman estimation (the blue line is below the red horizontal line), with the cost that MSE without attack is slightly larger. 
Moreover, the MSE of our proposed estimator without attack coincides with that of Oracle Kalman filter when $\gamma$ is large (e.g., larger than 50), which validates the performance of our estimator in normal operation. 

\begin{figure}[htpb]
	\centering
	\input{MSE_mag.tex}
	\caption{Estimation mean square error (MSE) in the absence and in the presence of attack with varying attack magnitude $\|a\|_{\infty}$. The parameter $\gamma$ is set as 5.} \label{fig:MSE_mag}

\end{figure}

Fig. \ref{fig:MSE_mag} illustrates the estimation mean square error (MSE) of our proposed estimator and Kalman filter with varying attack magnitude $\|a\|_{\infty}$. The
attack signal $a_4(k)$ is uniformly distributed in interval $(-\|a\|_{\infty}, \|a\|_{\infty})$. As the magnitude of attack signal increases, MSE of Kalman estimator increases significantly while MSE of our proposed secure estimator remains low. This result demonstrates that our proposed estimator holds low estimation error in the presence of malicious bias data in large magnitude. At the same time, in the absence of attack (when $\|a\|_{\infty}=0$), the estimation error is only slightly larger than that of the optimal Kalman filter.






\section{Conclusion}\label{sec:conclusion}
This paper considers LTI system with Gaussian noise against malicious attack on a subset of sensors. 
We improve upon the previous work by reducing the computational complexity in the designing phase while remaining low computational complexity for online estimation, under the assumption that all the geometric multiplicities of eigenvalues of $A$ are $1$ ($A$ is non-derogatory).
To achieve this, we prove that the span of the rows of $G_i$ is equivalent to the observable space associated to sensor $i$, based on which the canonical form $H_i$ is designed. The proposed estimator is formulated as a convex optimization problem based on $H_i$. We further prove that the proposed estimation is secure if the system is $2p$-sparse observable, which is easily checked due to the simple form of $H_i$. Moreover, in the absence of attack, the proposed estimation coincides with Kalman estimation for certain probability, which can be adjusted by tuning parameter $\gamma$ to balance between the performance with and without attack.




