An optimization framework for resilient batch estimation in Cyber-Physical Systems

TL;DR

This paper introduces a resilient set-valued estimation framework for linear time-varying systems affected by bounded process noise and potentially unbounded impulsive measurement noise, ensuring bounded estimation errors.

Contribution

It develops a novel class of set-valued estimators that maintain bounded errors despite unbounded impulsive measurement noise, enhancing robustness in cyber-physical system state estimation.

Findings

Estimator achieves resilience against unbounded impulsive noise

Error bounds depend on system observability and performance functions

Simulation confirms robustness and bounded error performance

Abstract

This paper proposes a class of resilient state estimators for LTV discrete-time systems. The dynamic equation of the system is assumed to be affected by a bounded process noise. As to the available measurements, they are potentially corrupted by a noise of both dense and impulsive natures. The latter in addition to being arbitrary in its form, need not be strictly bounded. In this setting, we construct the estimator as the set-valued map which associates to the measurements, the minimizing set of some appropriate performance functions. We consider a family of such performance functions each of which yielding a specific instance of the general estimator. It is then shown that the proposed class of estimators enjoys the property of resilience, that is, it induces an estimation error which, under certain conditions, is independent of the extreme values of the (impulsive) measurement noise. Hence, the estimation error may be bounded while the measurement noise is virtually unbounded. Moreover, we provide several error bounds (in different configurations) whose expressions depend explicitly on the degree of observability of the system being observed and on the considered performance function. Finally, a few simulation results are provided to illustrate the resilience property.

Equations

$$\Sigma:\left\{\begin{array}[]{r l}x_{t+1}&=A_{t}x_{t}+w_{t}\\ y_{t}&=C_{t}x_{t}+f_{t}\end{array}\right.$$

$$f_t=v_t+s_t,$$

$$X_{T-1}=\left(\begin{array}{cccc}{x}_0& x_1& \dots & x_{T-1}\end{array}\right)$$

$$Y_{T-1}=\left(\begin{array}{cccc}{y}_0& y_1& \cdots & y_{T-1}\end{array}\right)$$

$$V_{\Sigma }(Y,Z)=\lambda \sum _{t\in {\mathbb{T}}^{\prime }}{\varphi }_t(z_{t+1}-A_t{z}_t)+\sum _{t\in \mathbb{T}}{\psi }_t(y_t-C_t{z}_t)$$

$${\varphi }_t(z)=\varphi (W_{t}z)\forall z\in {\mathbb{R}}^n$$

$${\psi }_t(e)=\psi (V_{t}e)\forall e\in {\mathbb{R}}^{n_y},$$

$$\mathcal{E}(Y)=\underset{Z\in {\mathbb{R}}^{n\times T}}{\mathrm{argmin}}{V}_{\Sigma }(Y,Z).$$

$$\xi (Z)\ge q\left(\frac{1}{|\lambda |}\right)\xi (\lambda Z).$$

$$\xi (Z_1-Z_2)\ge {\gamma }_{\xi }\xi (Z_1)-\xi (Z_2).$$

$$\xi (Z)\ge Dq(\Vert Z\Vert )\forall Z\in {\mathbb{R}}^{a\times b}$$

$$D=\underset{\Vert Z\Vert =1}{\min }\xi (Z)>0.$$

$$\xi (Z)\ge q(\Vert Z\Vert )\xi (\frac{Z}{\Vert Z\Vert })\ge Dq(\Vert Z\Vert ).$$

$${\mathcal{O}}_{0,T-1}\triangleq {\left(\begin{array}{cccc}(C_0{)}^{\top }& (C_1{A}_0{)}^{\top }& \cdots & (C_{T-1}{A}_{T-2}\dots A_1{A}_0{)}^{\top }\end{array}\right)}^{\top }$$

$$V_{\Sigma }(0,Z)\ge q(\Vert z_0\Vert )$$

$$V_{\Sigma }(Y,Z)\ge F(Z)-\sum _{t\in \mathbb{T}}{\psi }_t(y_t),$$

$$F(Z)=\lambda \sum _{t\in {\mathbb{T}}^{\prime }}{\varphi }_t(z_{t+1}-A_t{z}_t)+{\gamma }_{\psi }\sum _{t\in \mathbb{T}}{\psi }_t(C_t{z}_t).$$

$$F(Z)\ge q(\Vert Z\Vert ).$$

$$V_{\Sigma }(Y,Z)\ge q(\Vert Z\Vert )-\sum _{t\in \mathbb{T}}{\psi }_t(y_t)$$

$$\|\hat{X}-X\|\leq g\big{(}\inf_{\Omega\in\mathcal{F}_{r}}d(F-\Omega)\big{)}\quad\forall\hat{X}\in\mathcal{E}(Y^{\star}+F)$$

$${\mathcal{F}}_r=\left\{F:|{\mathbb{T}}_0^c(F)|\le r\right\}.$$

$$\underset{\Omega \in {\mathcal{F}}_r}{\inf }d(F-\Omega )$$

$$\displaystyle=\inf_{\Omega\in\mathcal{F}_{r}}\Big{[}\!\!\sum_{t\in I_{r}(F)}\psi_{t}(\omega_{t}-f_{t})+\!\!\sum_{t\in I_{r}^{c}(F)}\psi_{t}(\omega_{t}-f_{t})\Big{]}$$

$$\displaystyle=\inf_{\begin{subarray}{c}\Omega\in\mathcal{F}_{r}\\ \mathbb{T}_{0}(\Omega)=I_{r}(F)\end{subarray}}\Big{[}\!\sum_{t\in I_{r}(F)}\psi_{t}(f_{t})+\!\!\!\sum_{t\in I_{r}^{c}(F)}\psi_{t}(\omega_{t}-f_{t})\Big{]}$$

$$=\sum _{t\in I_r(F)}{\psi }_t(f_t)$$

$$\Phi (Z)=\sum _{t\in {\mathbb{T}}^{\prime }}{\varphi }_t(z_{t+1}-A_t{z}_t)$$

$${\Psi }_{\mathbb{T}}(Z)=\sum _{t\in \mathbb{T}}{\psi }_t(C_t{z}_t)$$

$$H_{\Sigma }(Z)=\lambda {\gamma }_{\varphi }\Phi (Z)+{\gamma }_{\psi }{\Psi }_{\mathbb{T}}(Z).$$

$$p_r=\underset{\begin{array}{c}Z\in {\mathbb{R}}^{n\times T}\\ Z\ne 0\end{array}}{\sup }\underset{\begin{array}{c}\mathbb{I}\subset \mathbb{T}\\ |\mathbb{I}|=r\end{array}}{\sup }\frac{{\Psi }_{\mathbb{I}}(Z)}{H_{\Sigma }(Z)}$$

$${\mathbb{T}}_{\epsilon }=\left\{t\in \mathbb{T}:{\psi }_t(f_t)\le \epsilon \right\},$$

$$\lVert\hat{X}-X\rVert\leq h\Big{(}\dfrac{2\beta_{\Sigma}(\varepsilon)}{D\big{[}1-(1+\gamma_{\psi})p_{r}\big{]}}+\delta(\varepsilon)\Big{)}\quad\forall\hat{X}\in\mathcal{E}(Y)$$

$${\beta }_{\Sigma }(\epsilon )=\lambda \sum _{t\in {\mathbb{T}}^{\prime }}{\varphi }_t(w_t)+\sum _{t\in {\mathbb{T}}_{\epsilon }}{\psi }_t(f_t),$$

$$\delta (\epsilon )=\frac{1-{\gamma }_{\psi }}{D\left[1-(1+{\gamma }_{\psi })p_r\right]}\sum _{t\in {\mathbb{T}}_{\epsilon }^c}{\psi }_t(f_t)$$

Full text

I Introduction

Context

We consider in this work the problem of designing state estimators which would be resilient against an (unknown) sparse noise sequence affecting the measurements. By sparse noise we refer here to a signal sequence which is of impulsive nature, that is, a sequence which is most of the time equal to zero, except at a few instants where it can take on arbitrarily large values. The problem is relevant for example, in the supervision of Cyber-Physical Systems (CPS) [8]. In this application, the supervisory data may be collected by spatially distributed sensors and then sent to a distant processing unit through some communication network. During the transmission, the data may incur intermittent packet losses or adversarial attacks consisting in e.g., the injection of arbitrary signals. Beyond CPS, there are many other applications where the sparse noise model of uncertainty is relevant: robust statistics [14], hybrid system identification [1], intermittent sensor fault detection, etc.

Related works

The problem of estimating the state of CPS under attacks has been investigated recently through many different approaches. Since the measurements are assumed to be affected by a sequence of outliers which is sparse in time, a natural scheme of solution to the state estimation problem may be to first process the data so as to detect the occurrences of the nonzero instances of that sparse noise, remove the corrupted data and then proceed with classical estimation methods such as the Kalman filter or the Luenberger type of observer [17, 20]. While this approach sounds a priori reasonable, the main challenge remains to achieve an efficient detection and isolation of the outliers. Regarding the scenarios where the sporadic noise is modeled in a probabilistic setting, there exists a body of interesting results providing performance limits of estimation schemes [25, 18, 21].

Another category of approaches, which are inspired by some recent results in compressive sampling [7, 11], rely on sparsity-inducing optimization techniques. A striking feature of these methods is that they do not treat separately the tasks of detection, data cleaning and estimation. Instead, an implicit discrimination of the wrong data is induced by some specific properties of the to-be-minimized cost function. One of the first works that puts forward this approach for the resilient state estimation problem is the one reported in [10]. There, it is assumed that only a fixed number of sensors are subject to attacks (sparse over time but otherwise arbitrary disturbances). The challenge then resides in the fact that at each time instant, one does not know which sensor is compromised. Note however that the assumptions in [10] were quite restrictive as no dense process noise or measurement noise (other than the sparse attack signal) was considered. These limitations open ways for later extensions in many directions. For example, [24] suggests a reformulation which is argued to reduce computational cost by using the concept of event-triggered update ; [19] considers an observation model which includes dense noise along with the sparse attack signal. In [9], the assumption of a fixed number of attacked sensors is relaxed. Finally, the recent paper [13] proposes a unified framework for analyzing resilience capabilities of most of these (convex) optimization-based estimators. Although a bound on the estimation error was derived in this paper, it is not quantitatively related to the properties (e.g., observability) of the dynamic system being observed. The state estimation problem treated there is rather viewed as a linear regression problem similarly to [2, 5].

Contributions

The contributions of the current paper consist in the design and the analysis of a class of optimization-based resilient estimators for Linear Time-Varying (LTV) discrete-time systems. The available model of the system assumes bounded noise in both the dynamics and the observation equation with the latter being possibly affected, additionally, by an unknown but sparse attack signal. Contrary to the settings considered in some existing works, we did not impose here any restriction on the number of sensors which are subject to attacks, that is, any sensor can be compromised at any time. Note also that no statistical significance is attached to the uncertainties modeled by noise. In this setting, by generalizing our previous work reported in [16], the current paper proposes a general (robust) estimation framework for the state of LTV systems. We propose a class of state estimators constructed as the set-valued maps which associate to the output measurements, the minimizing set of some appropriate performance functions. A variety of performance functions are considered for the design of the estimator and handled in a unified analysis framework: convex nonsmooth/smooth loss functions and nonconvex saturated ones. Our main theoretical results concern the resilience analysis of the proposed class of estimators. We show that the estimation error associated with the new class of estimators can be made, under certain conditions, insensitive to the (possibly very large) amplitude of the sparse attack signal. The proposed analysis relies on new quantitative characterizations of the observability property of the system whose state is being observed. Although the derived error bounds may be conservative, they have the important advantage of being explicitly expressible in function of the properties of the considered dynamic system and those of the optimized loss functions. This makes them valuable qualitative tools for assessing the impact of the estimator’s design parameters and that of the system matrices on the quality of the estimation. For example, the proposed error bounds reflect the intuition that the more observable the system is with respect to the new criteria, the larger the number of instances of gross values (of the output noise) it can handle and the smaller the values of the bounds. Finally the paper shows that for some choice of the design functions (loss functions), some instances of the proposed family of estimators enjoy the exact recoverability property in the particular situation where the measurements are corrupted only by sparse noise. We present a condition for this property that can be numerically verified by means of convex optimization. Overall, in comparison with [13] which also considers resilient estimation though in a linear regression setting, we (i) introduce here an alternative definition of resilience, (ii) characterize quantitatively the impact of intrinsic properties (observability) of the system being observed on the quality of the estimation (iii) derive an explicit expression of a bound on the estimation error.

Outline

The rest of the paper is structured as follows. We start by introducing in Section II, the settings for the resilient state estimation problem. We then define in Section III the new class of optimization-based estimators proposed here to address this problem. The analysis of this new framework is presented in Section IV. In Section V, we further discuss the properties of a special constrained version of the initial class of estimators. In Section VI, we comment on the numerical verification of the conditions derived in the analysis part. Some numerical results are described in Section  VII and finally, concluding remarks are given in Section VIII.

Notation

$\mathbb{R}_{\geq 0}$ (respectively $\mathbb{R}_{>0}$) is the set of nonnegative (respectively positive) reals. $\mathbb{R}_{*}$ designates the set of real numbers excluding zero. We note $\mathbb{R}^{a}$ the set of (column) vectors with $a$ real elements and $\mathbb{R}^{a\times b}$, the set of real matrices with $a$ rows and $b$ columns. If $M\in\mathbb{R}^{a\times b}$, then $M^{\top}$ will designate the transposed matrix of $M$. $I$ will refer to the (square) identity matrix of appropriate dimension. The notation $\lVert\cdot\rVert$ will denote a norm over a given set (which will be specified when necessary). $\lVert\cdot\rVert_{p}$ denotes the $\ell_{p}$ norm (for $p\geq 1$) or the $\ell_{p}$ quasi-norm (for $0<p<1$) defined for $z=\begin{pmatrix}z_{1}&\cdots&z_{a}\end{pmatrix}$ in $\mathbb{R}^{a}$ by $\left\|z\right\|_{p}=\left(|z_{1}|^{p}+\cdots+|z_{a}|^{p}\right)^{1/p}$. The limit of this when $p\downarrow 0$ gives the so-called $\ell_{0}$-norm $\left\|\cdot\right\|_{0}$ of $z$, i.e., the number of nonzero entries in $z$. Its limit when $p\uparrow+\infty$ gives the infinity norm denoted $\left\|z\right\|_{\infty}$ and returning the maximum value of the $\left|z_{i}\right|$. For $x\in\mathbb{R}$, $e^{x}$ refers to the exponential function applied to $x$.

If $\mathcal{S}$ is a set, then $\mathcal{P}(\mathcal{S})$ is the power set of $\mathcal{S}$. If $\mathcal{S}$ is a finite set, the notation $|\mathcal{S}|$ refers to the cardinality of $\mathcal{S}$.

$\mathcal{K}_{\infty}$* functions [15].* We name $\mathcal{K}_{\infty}$ the set of functions $f:\mathbb{R}_{\geq 0}\rightarrow\mathbb{R}_{\geq 0}$ which are continuous, zero at zero, strictly increasing and satisfy $\lim_{\lambda\rightarrow+\infty}f(\lambda)=+\infty$. If $f\in\mathcal{K}_{\infty}$, then it admits an inverse, denoted here $f^{-1}$, which also lies in $\mathcal{K}_{\infty}$. Similarly, we use the notation $\mathcal{K}_{\operatorname*{sat},a}$ to denote the set of saturated functions $f:\mathbb{R}_{\geq 0}\rightarrow\interval{0}{a}$ which are continuous, zero at zero, strictly increasing on $\interval{0}{a}$ and such that $f(\lambda)=f(a)$ for all $\lambda\geq a$. Supremum. Given a function $f$ over $\mathbb{R}^{a}$ and a subset $\mathcal{S}$ of $\mathbb{R}^{a}$, the notation $\sup_{z\in\mathcal{S}}f(z)<b$, with $b\in\mathbb{R}$, will mean that for all $z$ in $\mathcal{S}$, $f(z)<b$. This notation includes the case where the supremum is $b$ but is not attained by any element of $\mathcal{S}$.

II The Resilient Estimation Problem

Consider a discrete-time Linear Time-Varying (LTV) system described by

[TABLE]

where $x_{t}\in\mathbb{R}^{n}$ is the state vector of the system at time $t$ and $y_{t}\in\mathbb{R}^{n_{y}}$ is the output vector at time $t$; $\{A_{t}\}$ and $\{C_{t}\}$ are families of matrices with appropriate dimensions; $\{w_{t}\}$ is an unobserved bounded noise sequence. As to $\{f_{t}\}$, it is regarded here as an (unobserved) arbitrary noise sequence affecting the measurements. For clarity of the exposition, it may be convenient to view $f_{t}$ as a combination of two types of sequences: a bounded sequence $\left\{v_{t}\right\}$ and a sparse sequence $\left\{s_{t}\right\}$ (this decomposition is indeed always possible for an arbitrary noise signal). Hence, we may write

[TABLE]

where $v_{t}$ is a sensor noise of moderate amplitude and $s_{t}$ is a sparse noise sequence in the sense that its (entrywise and/or timewise) components are mostly equal to zero but its nonzero elements can take on (possibly) arbitrarily large values. Such a sparse sequence $\left\{s_{t}\right\}$ may account for adversarial attacks in the same spirit as in [10, 13], intermittent sensor faults, or data losses, in particular when a communication network is involved in the data acquisition-transmission chain. In the sequel, we may also refer to $\left\{w_{t}\right\}$ and $\left\{v_{t}\right\}$ in (1) and (2) as dense noises and to the largest elements of $\left\{s_{t}\right\}$ as outliers.

For the sake of simplicity, the sparse (and potentially arbitrary large) noise is assumed here to affect only the measurement equation. Note however that the proposed analysis method can be extended to the more general scenario where the sparse noises may affect both the dynamics and the measurements.

Problem

The problem considered in this paper is the one of estimating the states $x_{0},\ldots,x_{T-1}$ of the system (1) on a time period $\mathbb{T}=\left\{0,\ldots,T-1\right\}$ given $T$ measurements $y_{0},...,y_{T-1}$ of the system output. We shall seek an accurate estimate of the state despite the uncertainties in the system equations (1) modeled by $w_{t}$ and $f_{t}$ the characteristics of which are informally described above. In particular, we would like the to-be-designed estimator to produce an estimate such that the estimation error is, when possible, independent of the maximum amplitude of $\left\{f_{t}\right\}$. Such an estimator will be called resilient, see Definition 2 for a formal characterization of this property.

Denote with

[TABLE]

the actual state trajectory of the system $\Sigma$ on a finite time horizon of length $T$. Similarly, we use the notation

[TABLE]

to refer to the collection of measurements on the same time horizon. The state estimation problem is approached here from an offline perspective, therefore $T$ is fixed. For the sake of simplicity, the $T$ index will be dropped from the variable names and it will be assumed that signal matrices without an index represent values on the period $\mathbb{T}=\{0,\ldots,T-1\}$. To simplify further the formulas, we also pose $\mathbb{T}^{\prime}=\{0,\ldots,T-2\}$ while $\mathbb{S}=\left\{1,\ldots,n_{y}\right\}$ will be a set indexing the sensors (or the rows of the matrices $C_{t}$ in (1)).

III Optimization-based approach to Resilient state estimation

III-A The state estimator

In this section we present an optimization-based framework for solving the state estimation problem defined above. To define formally the proposed state estimator, let us first introduce the to-be-minimized objective function. Given the matrices $\left\{(A_{t},C_{t})\right\}$ of the system (1) and $T$ output measurements $Y=\begin{pmatrix}y_{0}&\cdots&y_{T-1}\end{pmatrix}$, we consider a performance function $V_{\Sigma}:\mathbb{R}^{n_{y}\times T}\times\mathbb{R}^{n\times T}\rightarrow\mathbb{R}_{\geq 0}$ defined by

[TABLE]

where $Z=\begin{pmatrix}z_{0}&\cdots&z_{T-1}\end{pmatrix}\in\mathbb{R}^{n\times T}$ is a hypothetical trajectory matrix with $z_{i}$ denoting the $i$-th column of $Z$; $\lambda>0$ is a user-defined parameter which aims at balancing the contributions of the two terms involved in the expression of the performance index $V_{\Sigma}(Y,Z)$. $\{\phi_{t}\}$ and $\{\psi_{t}\}$ are two families of positive functions (called here loss functions) defined on $\mathbb{R}^{n}$ and $\mathbb{R}^{n_{y}}$ respectively. For the sake of simplicity, we will assume throughout the paper that for all $t$ in $\mathbb{T}$, $\phi_{t}$ and $\psi_{t}$ can be expressed by

[TABLE]

where $\phi:\mathbb{R}^{n}\rightarrow\mathbb{R}_{\geq 0}$ and $\psi:\mathbb{R}^{n_{y}}\rightarrow\mathbb{R}_{\geq 0}$ are two fixed loss functions and $\{W_{t}\}$ and $\{V_{t}\}$ are two families of nonsingular weighting matrices with appropriate dimensions.

Definition 1**.**

Given a system $\Sigma$ such as the one in (1) and given an output measurement matrix $Y\in\mathbb{R}^{n_{y}\times T}$, we define a state estimator to be a set-valued map $\mathcal{E}:\mathbb{R}^{n_{y}\times T}\rightarrow\mathcal{P}(\mathbb{R}^{n\times T})$ which maps $Y$ to a subset of the space of possible trajectories of the system.

We consider a class of state estimators defined by

[TABLE]

As such the estimator $\mathcal{E}$ is well-defined if for any fixed $Y$, $V_{\Sigma}(Y,Z)$ admits a non empty minimizing set, that is, if there exists at least one $Z^{\star}$ such that $V_{\Sigma}(Y,Z)\geq V_{\Sigma}(Y,Z^{\star})$ for all $Z\in\mathbb{R}^{n\times T}$. To ensure this property we will need to put an observability assumption on the system whose state is being estimated and require some further properties on the loss functions $\phi$ and $\psi$ entering in the definition of the objective function $V_{\Sigma}$.

III-B Well-definedness of the estimator

Let us start by stating the properties required for the loss functions involved in the definition of $V_{\Sigma}$. Due to the multiple usages that will be made of these properties, it is convenient to state them for a generic loss function defined on a set of matrices (of which vectors constitute a special case). Throughout this paper, a loss function is a positive function $\xi:\mathbb{R}^{a\times b}\rightarrow\mathbb{R}_{\geq 0}$ which will be required to satisfy a subset (depending of the specific usage) of the following properties:

1. (P1)

Positive definiteness: $\xi(0)=0$ and $\xi(Z)>0$ for all $Z\neq 0$ 2. (P2)

Continuity: $\xi$ is continuous 3. (P3)

Symmetry: $\xi(-Z)=\xi(Z)$ for all $Z\in\mathbb{R}^{a\times b}$ 4. (P4)

Generalized Homogeneity (GH): There exists a $\mathcal{K}_{\infty}$ function $q:\mathbb{R}_{\geq 0}\rightarrow\mathbb{R}_{\geq 0}$ such that for all $\lambda\in\mathbb{R}_{*}$ and for all $Z\in\mathbb{R}^{a\times b}$,

[TABLE]

1. (P5)

Generalized Triangle Inequality (GTI): There exists a positive real number $\gamma_{\xi}$ such that for all $Z_{1}$, $Z_{2}$ in $\mathbb{R}^{a\times b}$

[TABLE]

It can be usefully observed, for the future developments, that (10) can be equivalently written as $\xi(Z_{1}+Z_{2})\leq\gamma_{\xi}^{-1}\xi(Z_{1})+\gamma_{\xi}^{-1}\xi(Z_{2}).$

Examples of loss functions

Note that norms on $\mathbb{R}^{a\times b}$ satisfy naturally the properties (P1)–(P5) with $q:\lambda\mapsto\lambda$ and $\gamma_{\xi}=1$, hence yielding the classic homogeneity property and triangle inequality. It can also be checked that functions $\xi$ of the form $\xi(Z)=\left\|Z\right\|^{p}$ with $p>0$, fully qualify as loss functions in the sense that they fulfill all the properties (P1)–(P5). In this case, $\gamma_{\xi}$ in (10) can be taken equal to $2^{1-1/p}$ if $0<p\leq 1$ and $2^{1-p}$ otherwise. Lastly we note that if $\ell:\mathbb{R}^{a\times b}\rightarrow\mathbb{R}_{\geq 0}$ satisfies (P1)–(P3) and (P5), then so does the function $\xi$ defined by $\xi(Z)=1-e^{-\ell(Z)}$ (see Lemma 9 in the appendix). Similarly, saturated functions of the form $\xi(Z)=\min(\ell(Z),R_{0})$ for some $R_{0}>0$ satisfy (P1)–(P3) and (P5). In the case of convex functions, a link can be established between (P4) and (P5).

Lemma 1** ([16]).**

If $\xi:\mathbb{R}^{a\times b}\rightarrow\mathbb{R}_{\geq 0}$ is convex and satisfies property (P4) with a $\mathcal{K}_{\infty}$ function $q$, then it also satisfies (P5) with $\gamma_{\xi}=2q(1/2)$.

Observe that quadratic functions $\xi:\mathbb{R}^{a\times b}\rightarrow\mathbb{R}_{\geq 0}$ of the form $\xi(Z)=\operatorname*{Tr}(Z^{\top}QZ)$ with $Q\in\mathbb{R}^{a\times a}$ being a positive definite matrix and $\operatorname*{Tr}$ referring to the trace of a matrix, satisfy properties (P1)–(P4) with a $\mathcal{K}_{\infty}$ function $q:\lambda\mapsto\lambda^{2}$. Since such functions are convex, it follows from Lemma 1 above that they also verify (P5) for $\gamma_{\xi}=2q(1/2)=1/2$.

Remark 1**.**

In virtue of (6)-(7), the families of functions $\left\{\phi_{t}\right\}$ and $\left\{\psi_{t}\right\}$ satisfy (P1)–(P5) whenever $\phi$ and $\psi$ satisfy (P1)–(P5).

We now recall from [16] a technical lemma which will play a fundamental role in analyzing the properties of the estimator (8). In particular, our proof of well-definedness relies on this lemma.

Lemma 2** (Lower Bound of a loss function).**

Let $\xi:\mathbb{R}^{a\times b}\rightarrow\mathbb{R}_{\geq 0}$ be a function which has properties (P1)–(P2) and (P4) with a $\mathcal{K}_{\infty}$ function $q$. Then, for all norm $\lVert\cdot\rVert$ on $\mathbb{R}^{a\times b}$,

[TABLE]

where

[TABLE]

Proof.

We start by observing that the unit hypersphere $\mathcal{S}=\left\{Z\in\mathbb{R}^{a\times b}:\left\|Z\right\|=1\right\}$ is a compact set in the topology induced by the norm $\left\|\cdot\right\|$. By the extreme value theorem, $\xi$ being continuous, admits necessarily a minimum value on $\mathcal{S}$, i.e., there is $Z^{\star}\in\mathcal{S}$ such that $\xi(Z)\geq D\triangleq\xi(Z^{\star})>0$ for all $Z\in\mathcal{S}$. For any nonzero $Z\in\mathbb{R}^{a\times b}$, $\dfrac{Z}{\left\|Z\right\|}\in\mathcal{S}$ so that $\xi(\dfrac{Z}{\left\|Z\right\|})\geq D$. On the other hand, by the relaxed homogeneity of $\xi$,

[TABLE]

Moreover, this inequality holds for $Z=0$. It therefore holds true for any $Z\in\mathbb{R}^{a\times b}$. ∎

Proposition 1** (Well-definedness of the estimator).**

Let the loss functions $\phi$ and $\psi$ in (6)-(7) satisfy properties (P1)–(P5) and assume that the LTV system (1) is observable on $\interval{0}{T-1}$ in the sense that the observability matrix

[TABLE]

has full column rank. Then the estimator (8) is well-defined, i.e., the objective function $V_{\Sigma}(Y,\cdot)$ attains its minimum for any fixed $Y$.

Hence, the condition of the proposition guarantees that $\mathcal{E}(Y)$ is non empty for all $Y\in\mathbb{R}^{n_{y}\times T}$. Before proving this result, we first make the following observation.

Lemma 3** (Equivalent condition of Observability).**

Consider the objective function $V_{\Sigma}$ defined in (5) where $\left\{(\phi_{t},\psi_{t})\right\}$ are defined as in (6)-(7) with $\phi$ and $\psi$ satisfying (P1)–(P4). Then the following two statements are equivalent:

1. (i)

The system is observable on the time interval $\interval{0}{T-1}$. 2. (ii)

There exists a $\mathcal{K}_{\infty}$ function $q$ such that for all $Z\leavevmode\nobreak\ =\leavevmode\nobreak\ \begin{pmatrix}z_{0}&z_{1}&\ldots&z_{T-1}\end{pmatrix}$ in $\mathbb{R}^{n\times T}$,

[TABLE]

A proof of this lemma is reported in Appendix -B. The function $q$ can be interpreted here as a gain function which measures how much the system is observable with regards to the two families $\{\phi_{t}\}$ and $\{\psi_{t}\}$: the more the system is observable, the more $q$ amplifies its argument magnitude, making different trajectories more discernible.

Proof of Proposition 1: The idea of the proof is to show that $V_{\Sigma}(Y,\cdot)$ is coercive (i.e., continuous and radially unbounded) for any given $Y$ and then apply a result111Note that radial unboundedness is equivalent to level-boundedness in the terminology of [22]. in [22, Thm 1.9] to conclude on the attainability of the infimum (which certainly exists since $V_{\Sigma}(Y,\cdot)$ is a positive function). Clearly, $V_{\Sigma}(Y,\cdot)$ is continuous as a consequence of $\phi$ and $\psi$ being continuous by assumption (see property (P2)). We then just need to prove the radial unboundedness of $V_{\Sigma}(Y,\cdot)$, i.e., $\lim_{\left\|Z\right\|\rightarrow+\infty}V_{\Sigma}(Y,Z)=+\infty$ for an arbitrary norm $\left\|\cdot\right\|$ on the $Z$-space and for all fixed $Y$. Since $\psi$ satisfies property (P5), there exists a constant $\gamma_{\psi}>0$ such that $\psi_{t}(y_{t}-C_{t}z_{t})\geq\gamma_{\psi}\psi_{t}(C_{t}z_{t})-\psi_{t}(y_{t})$. Applying this property leads naturally to

[TABLE]

where

[TABLE]

It can then be shown (following a similar reasoning as in Appendix -B), under the observability assumption, that $F$ satisfies the conditions of Lemma 2. It follows that for any norm $\left\|\cdot\right\|$ on $\mathbb{R}^{n\times T}$, there exists a $\mathcal{K}_{\infty}$ function $q$ such that

[TABLE]

Combining this with the inequality above, we obtain that

[TABLE]

which implies the radial unboundedness of $V_{\Sigma}(Y,\cdot)$ for any fixed $Y$. Hence the estimator (8) is well-defined as stated. ∎

As it turns out from Proposition 1, observability of system (1) and properties (P1)–(P4) imposed on $\phi$ and $\psi$ ensure that $\mathcal{E}(Y)$ is a non empty set for any given $Y$. We then call any member $\hat{X}=\begin{pmatrix}\hat{x}_{0}&\hat{x}_{1}&\dots&\hat{x}_{T-1}\end{pmatrix}$ of $\mathcal{E}(Y)$, an estimate of the state trajectory $X$ of system (1) on the time interval $\mathbb{T}$. In particular, $\hat{x}_{t}$ is called an estimate of the state $x_{t}$ at time $t\in\mathbb{T}$.

To conclude this section, note that the definition of the estimator in (8) does not require any convexity assumption on the objective function $V_{\Sigma}$. Hence the theoretical analysis to be presented in the next sections does not make use of convexity either. However, we may prefer in practice to select convex loss functions $\phi$ and $\psi$. In effect, the elements of $\mathcal{E}(Y)$ are not necessarily expressible through an explicit formula. So, in practice one would resort instead to numerical solvers to approach the solution of the underlying optimization problem. And the numerical search process is known to be more efficient when $V_{\Sigma}(Y,Z)$ is a convex function of $Z$ [4, 12]. Nevertheless it is fair to recognize that nonconvex optimization methods can be successfully implemented as well though with less theoretical guarantees of reaching global optimality with general purpose solvers.

IV The resilience property of the proposed class of estimators

In this section, we prove that the state estimator proposed in (8) possesses the resilience property under some conditions. More specifically, our main result states that the estimation error $\hat{X}-X$, i.e., the difference between the real state trajectory and the estimated one, is upper bounded by a bound which does not depend on the amplitude of the outliers contained in $\left\{f_{t}\right\}$ provided that the number of such outliers is below some threshold.

IV-A Definition of the resilience of an estimator

Let us start with a formal definition of the resilience property for a state estimator of the form (8). For this purpose, let $Y^{\star}=\begin{pmatrix}y_{0}^{\star}&\ldots,y_{T-1}^{\star}\end{pmatrix}$ denote the noise-free output matrix of (1), i.e., the output defined by $y_{t}^{\star}=C_{t}x_{t}^{\star}$, $t=0,\ldots,T-1$, with $x_{t+1}^{\star}=A_{t}x_{t}^{\star}$ and $x_{0}^{\star}=x_{0}$ ($x_{0}$ being the true initial state of (1)). Let $\mathcal{F}_{r}$, a subset of $\mathbb{R}^{n_{y}\times T}$ containing [math], denote a matrix of measurement noise components.

Definition 2** (Resilience of an estimator).**

The set-valued estimator $\mathcal{E}$ defined in (8) is called resilient against the set $\mathcal{F}_{r}$ of measurement noise if there exists a $\mathcal{K}_{\infty}$ function $g$ such that, when the process noise $\left\{w_{t}\right\}$ is zero, it holds that for any measurement noise matrix $F\in\mathbb{R}^{n_{y}\times T}$,

[TABLE]

*with $\|\cdot\|$ denoting some norm, $d:\mathbb{R}^{n_{y}\times T}\rightarrow\mathbb{R}_{\geq 0}$ being a function subject to (P1)-(P5). Hence $\inf_{\Omega\in\mathcal{F}_{r}}d(F-\Omega)$ denotes some pseudo-distance from $F$ to the set $\mathcal{F}_{r}$. *

Since $0\in\mathcal{F}_{r}$, a consequence of property (16) is that $\mathcal{E}(Y^{\star})=\left\{X\right\}$ which follows from (16) for $F=0$. This fact expresses correctness of the estimator in a nominal situation, i.e., its ability to recover the true state matrix $X$ in the absence of any uncertainty in the a priori known model. Indeed this condition is guaranteed to hold if the system $\Sigma$ is observable over the considered observation time horizon $T$. Another key implication of condition (16) is that the estimation error associated with a resilient estimator is totally insensitive to any measurement noise matrix $F$ which lies in $\mathcal{F}_{r}$, that is, $\mathcal{E}(Y^{\star}+F)=\left\{X\right\}$ for all measurement noise $F\in\mathcal{F}_{r}$. Throughout this paper, we consider a set $\mathcal{F}_{r}$ defined as follows. For $F=\begin{pmatrix}f_{0}&\cdots&f_{T-1}\end{pmatrix}\in\mathbb{R}^{n_{y}\times T}$, let $\mathbb{T}_{0}^{c}(F)=\left\{t\in\mathbb{T}:\psi_{t}(f_{t})>0\right\}$ and $\mathbb{T}_{0}(F)=\left\{t\in\mathbb{T}:\psi_{t}(f_{t})=0\right\}$. For $r$ a positive integer, define $\mathcal{F}_{r}$ to be the set of matrices in $\mathbb{R}^{n_{y}\times T}$ having at most $r$ nonzero columns, i.e.,

[TABLE]

For the need of making explicit the resilience property (in the results to be presented) with respect to the set $\mathcal{F}_{r}$, we will need the following lemma.

Lemma 4**.**

Consider the set $\mathcal{F}_{r}$ of measurement noise matrix defined in (17) and select a (pseudo distance) function $d:\mathbb{R}^{n_{y}\times T}\rightarrow\mathbb{R}_{\geq 0}$ defined by $d(F)=\sum_{t\in\mathbb{T}}\psi_{t}(f_{t})$ with $\psi_{t}$ a function defined as in (7) and having the properties (P1)-(P5). Then $\inf_{\Omega\in\mathcal{F}_{r}}d(F-\Omega)$ is equal to the sum of the $T-r$ smallest terms in $\left\{\psi_{t}(f_{t}):t\in\mathbb{T}\right\}$.

Proof.

Let $I_{r}^{c}(F)$ denote the index set of the $r$ largest entries of the vector $\begin{pmatrix}\psi_{0}(f_{0})&\ldots&\psi_{T-1}(f_{T-1})\end{pmatrix}$ and $I_{r}(F)$ denote the index set of its $T-r$ smallest entries. Then, with the notation $\Omega=\begin{pmatrix}\omega_{0}&\cdots&\omega_{T-1}\end{pmatrix}$,

[TABLE]

where the notation $\mathbb{T}_{0}(\Omega)$ is defined in the lines preceding Eq. (17). The infimum is reached here for $\Omega\in\mathcal{F}_{r}$ such that $\omega_{t}=0$ $\forall t\in I_{r}(F)$ and $\omega_{t}=f_{t}$ $\forall t\in I_{r}^{c}(F)$. Hence $\inf_{\Omega\in\mathcal{F}_{r}}d(F-\Omega)$ is, as claimed, the sum of the $T-r$ smallest values among $\left\{\psi_{t}(f_{t}):t\in\mathbb{T}\right\}$. ∎

We will also introduce in the sequel a notion of approximate resilience of $\mathcal{E}$. This terminology refers to Definition 2 when the right hand side of (16) is modified as g\big{(}\inf_{\Omega\in\mathcal{F}_{r}}d(F-\Omega)+\delta\big{)} with $\delta$ some nonnegative real number.

IV-B Some notational conventions for the analysis

For convenience, let us introduce a few more notations. Let $\Phi:\mathbb{R}^{n\times T}\rightarrow\mathbb{R}_{\geq 0}$ and $\Psi_{\mathbb{T}}:\mathbb{R}^{n_{y}\times T}\rightarrow\mathbb{R}_{\geq 0}$ be defined by

[TABLE]

We also introduce the partial cost function $\Psi_{\mathbb{I}}$ defined for any $\mathbb{I}\subset\mathbb{T}$ by $\Psi_{\mathbb{I}}(Z)=\sum_{t\in\mathbb{I}}\psi_{t}(C_{t}z_{t})$. We will assume throughout the paper that the loss functions $\phi$ and $\psi$ satisfy a subset of the properties (P1)–(P5) and in particular, when they are required to satisfy the GTI (P5), we will denote the associated positive constants with $\gamma_{\phi}$ and $\gamma_{\psi}$ respectively. Finally, let us pose

[TABLE]

We will organize the resilience analysis for the estimator (8) along two cases: first, the scenario where the gross error vector sequence $\left\{s_{t}\right\}$ in (2) is block-sparse in time and then the situation where it is both componentwise and temporally sparse. To be more precise, if we denote with $S\in\mathbb{R}^{n_{y}\times T}$ the matrix formed from the sequence $\left\{s_{t}:t\in\mathbb{T}\right\}$, then the first case refers to columnwise block-sparsity of $S$ while the second is related to an entrywise sparsity. Note that the two cases coincide when the system of interest is single-input single-output (SISO).

IV-C Resilience to intermittent timewise block-sparse errors

We start by introducing the concept of $r$-Resilience index of an estimator such as the one in (8), a measure which depends of the system matrices, the structure of the performance function $V_{\Sigma}$ and on the loss functions $\phi$ and $\psi$.

Definition 3**.**

Let $r$ be a nonnegative integer. Assume that the system $\Sigma$ in (1) is observable on $\interval{0}{T-1}$. We then define the $r$-Resilience index of the estimator $\mathcal{E}$ in (8) (when applied to $\Sigma$) to be the real number $p_{r}$ given by

[TABLE]

where $H_{\Sigma}$ is as defined in (20). The supremum is taken here over all nonzero $Z$ in $\mathbb{R}^{n\times T}$ and over all subsets $\mathbb{I}$ of $\mathbb{T}$ with cardinality equal to $r$.

The index $p_{r}$ can be interpreted as a quantitative measure of the observability of the system $\Sigma$. The observability is needed here to ensure that the denominator $H_{\Sigma}(Z)$ of (21) is different from zero whenever $Z\neq 0$. Furthermore, it should be remarked that $\Psi_{\mathbb{I}}(Z)\leq H_{\Sigma}(Z)$ for any $\mathbb{I}\subset\mathbb{T}$, which implies that the defining suprema of $p_{r}$ are well-defined. Note that $p_{r}$ is an increasing function of $r$ and satisfies $p_{0}=0$ and $p_{T}=1$. More discussions on the numerical evaluation of $p_{r}$ are deferred to Section VI.

In order to state the resilience result for the estimator (8) when applied to system $\Sigma$, let us introduce a last notation to be used in the analysis. Let $\varepsilon\geq 0$ be a given number. For any admissible sequence $\left\{f_{t}\right\}_{t\in\mathbb{T}}$ in (1), we can split the time index set $\mathbb{T}$ into two disjoint label sets,

[TABLE]

indexing those $f_{t}$ which are upper bounded by $\varepsilon$ and $\mathbb{T}_{\varepsilon}^{c}=\left\{t\in\mathbb{T}:\psi_{t}(f_{t})>\varepsilon\right\}$ indexing those $f_{t}$ which are possibly unbounded. It is important to keep in mind that $\varepsilon$ is just a parameter for decomposing the noise sequence in two parts in view of the analysis (and not necessarily a bound on elements of the sequence $\left\{\psi(f_{t})\right\}$). For example, taking $\varepsilon=0$ would be appropriate for analyzing the properties of the estimator when $f_{t}$ is strictly sparse and each of its nonzero elements is treated as an outlier.

Theorem 1** (Upper bound on the estimation error).**

*Consider the system $\Sigma$ defined by (1) with output $Y$ together with the state estimator (8) in which the loss functions $\phi$ and $\psi$ are assumed to obey (P1)–(P5). Denote with $\gamma_{\phi}$ and $\gamma_{\psi}$ the constants associated with the GTI (P5) and $q_{\phi}$ and $q_{\psi}$ the $\mathcal{K}_{\infty}$ functions associated with the GH (P4) for $\phi$ and $\psi$ respectively. Let $\varepsilon\geq 0$ and set $r=|\mathbb{T}_{\varepsilon}^{c}|$.

If the system is observable on $\interval{0}{T-1}$ and $p_{r}<1/(1+\gamma_{\psi})$, then for any norm $\lVert\cdot\rVert$ on $\mathbb{R}^{n\times T}$,*

[TABLE]

with $X$ denoting the true state matrix from (1), $D=\min_{\left\|Z\right\|=1}H_{\Sigma}(Z)>0$ and $\beta_{\Sigma}(\varepsilon)$, $\delta(\varepsilon)$ and $h$ being defined by

[TABLE]

Proof.

Let $\hat{X}$ in $\mathcal{E}(Y)$. By definition of $\mathcal{E}$ in (8), we have $V_{\Sigma}(Y,\hat{X})\leq V_{\Sigma}(Y,X)$, which gives explicitly

[TABLE]

Using the fact that $x_{t+1}=A_{t}x_{t}+w_{t}$ from (1) and applying the GTI and the symmetry properties of $\phi_{t}$, we can write

[TABLE]

with $e_{t}=\hat{x}_{t}-{x}_{t}$. It follows that the first term on the left hand side of (27) is lower bounded as follows

[TABLE]

Similarly, by making use of (1), observe that $\psi_{t}(y_{t}-C_{t}\hat{x}_{t})=\psi_{t}(f_{t}-C_{t}e_{t})$. We now apply the GTI and symmetry of $\psi_{t}$ in two different ways depending on whether $t$ belongs to $\mathbb{T}_{\varepsilon}$ or $\mathbb{T}_{\varepsilon}^{c}$:

[TABLE]

These inequalities imply that the second term on the left hand side of (27) is lower bounded as follows

[TABLE]

Combining (27), (28) and (29) gives

[TABLE]

which, by using (19), (20), (24), can be written as

[TABLE]

with $E=\begin{pmatrix}e_{0}&e_{1}&\cdots&e_{T-1}\end{pmatrix}$. As $\mathbb{T}_{\varepsilon}^{c}$ has $r$ elements, applying the definition of $p_{r}$ gives

[TABLE]

By the assumption that $p_{r}<1/(1+\gamma_{\psi})$ we have that

$1-(1+\gamma_{\psi})p_{r}>0$, and consequently, that

[TABLE]

Given that the system is observable on $\interval{0}{T-1}$, it can be shown, thanks to Lemma 7 in the Appendix, that $H_{\Sigma}$ satisfies properties (P1)–(P4) (the proof of this is quite similar to that of Lemma 3 in Appendix -B). We can therefore apply Lemma 2 to conclude that for any norm $\lVert\cdot\rVert$

[TABLE]

with $D$ defined by $D=\min_{\lVert Z\rVert=1}H_{\Sigma}(Z)$ and $q^{\prime}(\alpha)=\min\{q_{\phi}(\alpha),q_{\psi}(\alpha)\}$. Finally, the result follows by selecting $h$ to be $h=q^{\prime-1}$ with $q^{\prime-1}$ denoting the inverse of $q^{\prime}$, which can be simplified to match its definition in (26). ∎

Strict resilience

Now we state our (strict) resilience result as a consequence of Theorem 1 when the output error-measuring loss function $\psi$ satisfies the triangle inequality.

Corollary 1** (Resilience property).**

Let the conditions of Theorem 1 hold with the additional requirement that $\gamma_{\psi}=1$. Then

[TABLE]

Proof.

The proof is immediate by considering the bound in (23) and observing that $\delta(\varepsilon)$ expressed in (25) vanishes when $\gamma_{\psi}=1$, hence eliminating completely the contribution of the extreme values of $\left\{f_{t}\right\}$ to the error bound. This gives immediately (33). It remains now to make clear that (33) is consistent with the requirement (16) of Definition 2. For this purpose note that the bound in (33) can be written as $g(\beta_{\Sigma}(\varepsilon))$ with $g\in\mathcal{K}_{\infty}$ defined by g(\alpha)=h\big{(}2\alpha/\big{(}D(1-2p_{r})\big{)}\big{)}. Moreover, $\beta_{\Sigma}(\varepsilon)$ reduces to $\sum_{t\in\mathbb{T}_{\varepsilon}}\psi_{t}(f_{t})=\inf_{\Omega\in\mathcal{F}_{r}}d(F-\Omega)$ when the process noise is zero (see (24) and Lemma 4). Hence, $\mathcal{E}$ qualifies, in the sense of Definition 2, as an estimator which is resilient to the set $\mathcal{F}_{r}$ of measurement noise defined in (17). ∎

The resilience property of the estimator (8) lies here in the fact that, under the conditions of Theorem 1 and Corollary 1, the bound in (33) on the estimation error does not depend on the magnitudes of the extreme values of the noise sequence $\left\{f_{t}\right\}_{t\in\mathbb{T}}$. Considering in particular the function $\beta_{\Sigma}(\varepsilon)$, we remark that it can be overestimated as follows

[TABLE]

The first term on the left hand side of (34) represents the uncertainty brought by the dense noise $\left\{w_{t}\right\}$ over the whole state trajectory. It is bounded since $\left\{w_{t}\right\}$ is bounded by assumption (see the description of the system in Section II). The second term is a bound on the sum of those instances of $f_{t}$ whose magnitude is smaller that $\varepsilon$.

Because $\beta_{\Sigma}$ is a function of $\varepsilon$, the bound in (33) represents indeed a family of bounds parameterized by $\varepsilon$. Since $\varepsilon$ is a mere analysis device, a question would be how to select it for the analysis to achieve the smallest bound. Such favorable values, say $\varepsilon^{\star}$, satisfy

[TABLE]

Another interesting point is that the inequality stated by Theorem 1 holds for any norm $\left\|\cdot\right\|$ on $\mathbb{R}^{n\times T}$. Note though that the value of the bound depends (through the parameter $D$) on the specific norm used to measure the estimation error. Moreover, different choices of the performance-measuring norm will result in different geometric forms for the uncertain set, that is, the ball (in the chosen norm) centered at the true state with radius equal to the upper bound displayed in (33).

We also observe that the smaller the parameter $p_{r}$ is, the tighter the error bound will be, which suggests that the estimator is more resilient when $p_{r}$ is lower. A similar reasoning applies to the number $D$ which is desired to be large here. These two parameters (i.e., $p_{r}$ and $D$) reflect properties of the system whose state is being estimated. They can be interpreted, to some extent, as measures of the degree of observability of the system. In conclusion, the estimator inherits partially its resilience property from characteristics of the system being observed. This is consistent with the well-known fact that the more observable a system is, the more robustly its state can be estimated from output measurements.

Approximate resilience

As discussed above, the triangle inequality property of the loss function $\psi$ is fundamental for achieving strict resilience. When $\psi$ does not satisfy this property (i.e., when $\gamma_{\psi}\neq 1$), the term $\delta(\varepsilon)$ in (23) is unlikely to vanish completely. However we can prevent it from growing excessively by an appropriate choice of $\psi$ in (7). To see this, assume for example that $\psi$ is defined by $\psi(y)=1-e^{-\ell(y)}$. Then since $\psi(y)\leq 1$ for all $y$, $\delta(\varepsilon)$ saturates to a constant value regardless of how large the $f_{t}$ are for $t\in\mathbb{T}_{\varepsilon}^{c}$. On the other hand, this choice introduces a new technical challenge related to the fact that the function $q^{\prime}$ in (32) is no longer a $\mathcal{K}_{\infty}$ function but a bounded (saturated) function. Handling this situation will require some additional condition on the upper bound in (31). To sum up, by selecting a saturated loss function for $\psi$, we obtain the following approximate resilience result.

Corollary 2** (Case where $\gamma_{\psi}\neq 1$).**

Let the conditions of Theorem 1 hold. Assume further that the loss function $\psi$ in (7) is defined by $\psi(y)=1-e^{-\ell(y)}$ where $\ell:\mathbb{R}^{n_{y}}\rightarrow\mathbb{R}_{\geq 0}$ satisfies (P1)–(P5). In particular, assume that property (P4) is satisfied by $\ell$ with a $\mathcal{K}_{\infty}$ function $q$ such that (9) is an equality relation. Also, let $\varepsilon\geq 0$ be such that

[TABLE]

where $r=r(\varepsilon)=|\mathbb{T}_{\varepsilon}^{c}|$ and $r^{o}(\varepsilon)=\max_{t\in\mathbb{T}_{\varepsilon}^{c}}\psi_{t}(f_{t})\leq 1$. Then there exists a continuous and strictly increasing function $h_{\operatorname*{sat}}:\interval{0}{1}\rightarrow\interval{0}{1}$ (obeying $h_{\operatorname*{sat}}(0)=0$ and $h_{\operatorname*{sat}}(1)=1$) such that for any norm $\lVert\cdot\rVert$ on $\mathbb{R}^{n\times T}$,

[TABLE]

with $D$ in (35) defined as in the proof of Theorem 1 using the norm $\lVert\cdot\rVert$.

Proof.

That the particular function $\psi$ specified in the statement of the corollary satisfies the properties (P1)–(P3) and (P5) is a question which is fully answered by Lemma 9 in Section -C of the appendix. Consequently, let us observe that the inequality (31) arising in the proof of Theorem 1 still holds true here. As to (32), it also holds as well but with the slight difference that $q^{\prime}$ is just a saturated function in $\mathcal{K}_{\operatorname*{sat},1}$ (as defined in the notation section) with bounded range $\interval{0}{1}$. This results in fact from Lemma 9 and the proof of Lemma 2. We can therefore write

[TABLE]

with $q^{\prime}\in\mathcal{K}_{\operatorname*{sat},1}$. Note from the definition of the class $\mathcal{K}_{\operatorname*{sat},1}$, that $q^{\prime}(\left\|E\right\|)<1$ implies that $\left\|E\right\|<1$ (since otherwise we would have $q^{\prime}(\left\|E\right\|)=1$). Letting $h_{\operatorname*{sat}}$ be the restriction of such a function $q^{\prime}$ on $\interval{0}{1}$, we have $q^{\prime}(\left\|E\right\|)=h_{\operatorname*{sat}}(\left\|E\right\|)\leq b(\varepsilon)$ with $h_{\operatorname*{sat}}$ being invertible. We can now apply $h_{\operatorname*{sat}}^{-1}$ to each member of this inequality to reach the desired result since $b(\varepsilon)$ lies in the range of $h_{\operatorname*{sat}}$. ∎

IV-D Resilience to attacks on the individual sensors

We now consider the situation where the matrix $S\in\mathbb{R}^{n_{y}\times T}$ formed from $\left\{s_{t}\right\}$ in (2) may be sparse entrywise i.e., a relatively important fraction of the entries of $S$ are equal to zero222In this case, the sparsity is expressed in term of fraction of nonzero entries in the matrix $S$ whereas in the timewise block-sparsity case, the sparsity level is measured in term of the fraction of nonzero columns in $S$.. This case is relevant when any individual sensor may be faulty (or compromised by an attacker) at any time. To address the resilient state estimation problem in this scenario, we select the loss functions $\psi_{t}$ to have a separable structure. To be more specific, let $\psi_{t}$ be such that for any $e=[\begin{matrix}e_{1}&\cdots&e_{n_{y}}\end{matrix}]\in\mathbb{R}^{n_{y}}$

[TABLE]

where, consistently with (7), $\psi_{ti}(e_{i})=\psi_{i}^{\circ}(V_{ti}e_{i})$ with $V_{ti}\in\mathbb{R}_{>0}$ and $\psi_{i}^{\circ}:\mathbb{R}\rightarrow\mathbb{R}_{+}$, $i=1,\ldots,n_{y}$, being some loss functions on $\mathbb{R}$ enjoying the properties (P1)–(P5). As in the statement of Corollary 1, we shall require that $\gamma_{\psi_{i}^{\circ}}=1$. It follows that one can set $\psi_{i}^{\circ}$ to be the absolute value without loss of generality. Let therefore set $\psi_{i}^{\circ}(e_{i})=|e_{i}|$ so that $\psi_{ti}(e_{i})=\left|V_{ti}e_{i}\right|$ and

[TABLE]

with $V_{t}$ being a diagonal matrix having the $V_{ti}$, $i=1,\ldots,n_{y}$, on its diagonal.

To state the resilience property in this particular setting, we partition the index set $\mathbb{T}\times\mathbb{S}$ of the entries of $S$ as

[TABLE]

with $f_{ti}$ denoting the $i$-th entry of the vector $f_{t}\in\mathbb{R}^{n_{y}}$. Also, in order to account for the specificity of the new scenario, let us refine slightly the $r$-Resilience index (21) to be

[TABLE]

where $H_{\Sigma}$ is defined as in (20) from $\psi_{t}$ in (38) and $c_{ti}^{\top}$ is $i$-th row of the observation matrix $C_{t}$. The difference between $p_{r}$ in (21) and $\tilde{p}_{r}$ in (40) resides in the index sets for counting possible error occurrences which are $\mathbb{T}$ and $\mathbb{T}\times\mathbb{S}$, respectively.

With these notations, we can provide the following theorem which is the analog of Corollary 1 in the case where the disturbance matrix $S$ (see Eq. (2)) is entrywise sparse.

Theorem 2** (Upper bound on the estimation error: Separable case).**

*Consider the system $\Sigma$ defined by (1) with output $Y$ together with the state estimator (8) in which $\phi$ is assumed to obey (P1)–(P5) and $\psi$ is defined as in (38). Let $\varepsilon\geq 0$ and set $r=|\Lambda_{\varepsilon}^{c}|$ with $\Lambda_{\varepsilon}^{c}$ defined in (39).

If the system is observable on $\interval{0}{T-1}$ and if $\tilde{p}_{r}<1/2$, then there exists a $\mathcal{K}_{\infty}$ function $\tilde{h}$ such that for all norm $\lVert\cdot\rVert$ on $\mathbb{R}^{n\times T}$,*

[TABLE]

with $X$ denoting the true state matrix from (1) and $\tilde{\beta}_{\Sigma}(\varepsilon)$ defined by

[TABLE]

$\tilde{D}$* and $\tilde{h}$ are defined as in the statement of Theorem 1 with $H_{\Sigma}$ being constructed from $\psi$ in (38).*

To some extent, Theorem 2 can be viewed as a special case of Theorem 1 in which the function $\psi$ is taken to be the $\ell_{1}$ norm and the data set is modified to be $\mathbb{T}\times\mathbb{S}$. Hence the proof follows a similar line of arguments as that of Theorem 2. Again it is not hard to see that the result of Theorem 2 implies the property of resilience with respect to the set $\mathcal{F}_{r}$ in (17) of measurement noise in the sense of Definition 2 (see the proof of Corollary 1).

An interesting property of the estimator can be stated in the absence of dense noise, i.e., when only the sparse noise is active:

Corollary 3**.**

Consider the system $\Sigma$ defined by (1) and let $r=|\Lambda_{0}^{c}|$ (which means that we consider every nonzero occurrence of $f_{it}$ as an outlier by taking $\varepsilon=0$ in (39)). If the conditions of Theorem 2 hold, $\tilde{p}_{r}<1/2$, and if $w_{t}=0$ in (1) for all $t$, then the estimator defined by (8) retrieves exactly the state trajectory of the system, i.e., $\mathcal{E}(Y)=\left\{X\right\}$.

Proof.

This follows directly from the fact that $\tilde{\beta}_{\Sigma}(0)=0$ in the case where there is no dense noise $w_{t}$ and $\varepsilon=0$. ∎

Therefore, the estimator (8) has the exact recoverability property, that is, it is able to recover exactly the true state of the system (1) when only the sparse noise is active in the measurement equation provided that the number $r=|\Lambda_{0}^{c}|$ of nonzero in the sequence $\left\{f_{ti}\right\}_{(t,i)\in\mathbb{T}\times\mathbb{S}}$ is small enough for $\tilde{p}_{r}$ to be less than $1/2$. According to our analysis, the number of outliers that can be handled by the estimator in this case can be underestimated by

[TABLE]

V A special variant of the estimator $\mathcal{E}$

In this section, we consider a constrained reformulation of the estimator $\mathcal{E}$ defined in (8). As will be shown shortly, this reformulation also enjoys the resilience property but under a condition which is more easily verifiable from a numerical perspective.

We start by considering the simple scenario where the process noise $w_{t}$ in (1) is identically equal to zero and the sequence $\left\{f_{t}\right\}$ is sparse in the sense that its dense component $v_{t}$ displayed in (2) does not exist. In this setting we can obtain a more resilient (to sparse noise in the measurement) estimator than (8) by making it aware of the absence of dense process noise. This can be achieved by contraining the searched state matrix to be in the set $\mathcal{Z}_{\Sigma}\subset\mathbb{R}^{n\times T}$ defined by

[TABLE]

of possible trajectories starting in any initial state $z_{0}\in\mathbb{R}^{n}$. Following this idea, we consider the estimator $\mathcal{E}^{\circ}$ defined by

[TABLE]

Then $\mathcal{E}^{\circ}(Y)$ can be rewritten more simply in the form

[TABLE]

where

[TABLE]

with

[TABLE]

for all $t\geq 1$ and $M_{0}=C_{0}$. Hence the estimation of the state trajectory reduces to estimating the initial state $x_{0}$. This can be viewed as a robust regression problem, like the ones discussed in [13, 2]. Generalizing a result in [2], we derive next a necessary and sufficient condition for exact recovery of the true state, which holds if and only if $\operatorname*{arg\,min}_{z\in\mathbb{R}^{n}}V_{\Sigma}^{\circ}(Y,z)=\{x_{0}\}$ with $x_{0}$ being the exact initial state of the system $\Sigma$. To this end, we first introduce the concept of concentration ratio of a collection of matrices with respect to a loss function. A notational convention will be necessary for the statement of this property: for any subset $\mathbb{I}$ of $\mathbb{T}$, let

[TABLE]

Definition 4** ($r$-th concentration ratio).**

Let $\{\psi_{t}\}$ be a family of loss functions defined by (7) in which $\psi$ is assumed to satisfy (P1), (P3) and (P5) with constant $\gamma_{\psi}=1$. Let $M=\left\{M_{t}\right\}_{t\in\mathbb{T}}$ be a sequence of matrices such that the function $\Psi_{\mathbb{T}}^{\circ}$ defined in (46) is positive definite. We call $r$-th concentration ratio of $M$, the number defined by

[TABLE]

At a fixed $r$, $\nu_{r}(M)$ quantifies a genericity property for the sequence $M=\left\{M_{t}\right\}_{t\in\mathbb{T}}$. In view of the particular structure of the collection $M$ in (45), note that $\Psi_{\mathbb{T}}^{\circ}$ is positive definite whenever the system $\Sigma$ is observable on $\mathbb{T}$. Furthermore, $\nu_{r}(M)$ can be interpreted to some extent, as a quantitative measure of observability. It is indeed all the smaller as the system is strongly observable. To see this, recall from Lemma 3 that if the system is observable on $\interval{0}{T-1}$, then for all $Z\in\mathcal{Z}_{\Sigma}$ initiated from $z$ in $\mathbb{R}^{n}$, we have $V_{\Sigma}(0,Z)=\Psi_{\mathbb{T}}^{\circ}(z)\geq q(\lVert z\rVert)$ for some $\mathcal{K}_{\infty}$ function $q$. It follows that

[TABLE]

Hence the more observable (i.e., the larger the gain function $q$), the smaller $\nu_{r}(M)$.

For all $(Y,z_{0})\in\mathbb{R}^{n_{y}\times T}\times\mathbb{R}^{n}$ with $Y$ expressed columnwise in the form $Y=\begin{pmatrix}y_{0}&\cdots&y_{T-1}\end{pmatrix}$, consider now the following notations:

[TABLE]

Theorem 3** (Exact Recoverability Condition).**

Consider the cost function (44) where $M=\left\{M_{t}\right\}$ is assumed to have been constructed as in (45) from the matrices of system (1). Assume that the loss functions $\{\psi_{t}\}$ involved in (44) are defined by (7) in which $\psi$ is assumed to satisfy (P1), (P3) and (P5) with constant $\gamma_{\psi}=1$. Let $r$ be a positive integer. If the system (1) is observable on $\interval{0}{T-1}$, then the two following propositions are equivalent:

1. (i)

For all $Y$ in $\mathbb{R}^{n_{y}\times T}$ and all $z_{0}$ in $\mathbb{R}^{n}$,

[TABLE] 2. (ii)

The index $\nu_{r}(M)$ satisfies

[TABLE]

Proof.

(i) $\Rightarrow$ (ii): Assume that (i) holds. Consider an arbitrary subset $\mathbb{I}$ of $\mathbb{T}$ such that $|\mathbb{I}|\leq r$. Let $z_{0}\neq 0$ be a vector in $\mathbb{R}^{n}$. Construct a sequence $Y$ in $\mathbb{R}^{n_{y}\times T}$ such that $y_{t}=0$ if $t\in\mathbb{I}$ and $y_{t}=M_{t}z_{0}$ otherwise. Then $\mathcal{I}^{c}(Y,z_{0})\subset\mathbb{I}$, so that $|\mathcal{I}^{c}(Y,z_{0})|\leq r$. It then follows from (i) that \operatorname*{arg\,min}_{z\in\mathbb{R}^{n}}V_{\Sigma}^{\circ}(Y,z)=\big{\{}z_{0}\big{\}} which means that $V_{\Sigma}^{\circ}(Y,z_{0})<V_{\Sigma}^{\circ}(Y,z)$ for all $z\in\mathbb{R}^{n}$, $z\neq z_{0}$. In particular, $V_{\Sigma}^{\circ}(Y,z_{0})<V_{\Sigma}^{\circ}(Y,0)$ which, by taking into account the definition of $Y$, gives

[TABLE]

where $\mathbb{I}^{c}=\mathbb{T}\setminus\mathbb{I}$. Since $\Psi_{\mathbb{T}}^{\circ}(z_{0})=\Psi_{\mathbb{I}}^{\circ}(z_{0})+\Psi_{\mathbb{I}^{c}}^{\circ}(z_{0})$, we see that

[TABLE]

This reasoning works for every nonzero $z_{0}$ and for every subset $\mathbb{I}$ of $\mathbb{T}$. We can hence conclude that $\nu_{r}(M)<1/2$.

(ii) $\Rightarrow$ (i): Assume that (ii) holds. Let $(Y,z_{0})\in\mathbb{R}^{n_{y}\times T}\times\mathbb{R}^{n}$ be such that $|\mathcal{I}^{c}(Y,z_{0})|\leq r$. We then need to prove that \operatorname*{arg\,min}_{z\in\mathbb{R}^{n}}V_{\Sigma}^{\circ}(Y,z)=\big{\{}z_{0}\big{\}}. Since the assertion (ii) is assumed true, it follows from (47) and (50) that

[TABLE]

where, for simplicity, we have posed $\mathbb{I}^{c}=\mathcal{I}^{c}(Y,z_{0})$. In the derivation of (51), we have used the obvious fact that $r_{1}\leq r_{2}$ $\Rightarrow$ $\nu_{r_{1}}(M)\leq\nu_{r_{2}}(M)$. If we pose $\mathbb{I}=\mathcal{I}^{0}(Y,z_{0})=\mathbb{T}\setminus\mathbb{I}^{c}$, then the inequality (51) is equivalent to

[TABLE]

Now we observe that for all $t$ in $\mathbb{I}=\mathcal{I}^{0}(Y,z_{0})$, $y_{t}=M_{t}z_{0}$, so that \psi_{t}(M_{t}z_{0}^{\prime})=\psi_{t}\big{(}y_{t}-M_{t}(z_{0}+z_{0}^{\prime})\big{)}. On the other hand, for $t\in\mathbb{I}^{c}=\mathcal{I}^{c}(Y,z_{0})$, if we apply the GTI (10) with $\gamma_{\psi}=1$, we obtain

[TABLE]

Combining all these remarks with (52) yields

[TABLE]

Rearranging this gives $V_{\Sigma}^{\circ}(Y,z_{0})<V_{\Sigma}^{\circ}(Y,z_{0}+z_{0}^{\prime})$ for all $z_{0}^{\prime}\in\mathbb{R}^{n}$ with $z_{0}^{\prime}\neq z_{0}$. This is equivalent to \operatorname*{arg\,min}_{z\in\mathbb{R}^{n}}V_{\Sigma}^{\circ}(Y,z)=\big{\{}z_{0}\big{\}}. Hence (ii) holds as claimed. ∎

From the statement of Theorem 3, we infer that under the assumption that only the sparse noise $\left\{s_{t}\right\}$ is active (i.e., there is no dense noise $(w_{t},v_{t})$) in the system equations (1), $\mathcal{E}^{\circ}(Y)=\left\{X\right\}$ whenever $\nu_{r}(M)<1/2$, i.e, the estimator $\mathcal{E}^{\circ}$ returns exactly the true state. For a given system, if one can evaluate numerically the index $\nu_{r}(M)$, then it becomes possible to assess the number $r_{\max}\triangleq\max\left\{r:\nu_{r}(M)<1/2\right\}$ of gross errors that can be corrected by the estimator $\mathcal{E}^{\circ}$ when applied to that system. We will get back to the computational matter in Section VI.

V-A Special case of $\ell_{0}$-norm loss based estimator

Consider the special case where the loss functions $\left\{\psi_{t}\right\}$ are defined, for all $t\in\mathbb{T}$, by

[TABLE]

This corresponds to the block $\ell_{0}$-norm. Note that such functions satisfy the assumptions (P1), (P3) and (P5) requested in the definition 4 of $\nu_{r}(M)$ and in the statement of Theorem 3. Hence $\nu_{r}(M)$ is well-defined in this case.

Corollary 4**.**

Consider system (1) under the assumption that $w_{t}=0$ for all $t$. Assume observability of the system on $\interval{0}{T-1}$. Consider the estimator (43) in which the cost $V_{\Sigma}^{\circ}$ is defined from the family of loss functions $\left\{\psi_{t}\right\}$ expressed in (53). Then for all $(Y,z_{0})\in\mathbb{R}^{n_{y}\times T}\times\mathbb{R}^{n}$,

[TABLE]

where $\mu(M)$ defined by

[TABLE]

is the minimum number $k$ such that any matrix $M_{I}\in\mathbb{R}^{|I|n_{y}\times n}$ formed by stacking vertically the matrices of the collection $\left\{M_{t}:t\in I\right\}$ indexed by $I\subset\mathbb{T}$ with $|I|=k$, has full column rank.

Proof.

Let us start by observing that with the particular loss functions invoked in the statement of the corollary, $\Psi_{\mathbb{T}}^{\circ}(z)$ denotes the number of $t\in\mathbb{T}$ for which $\psi_{t}(M_{t}z)\neq 0$. It follows from the definition of $\mu(M)$ that $\Psi_{\mathbb{T}}^{\circ}(z)\geq T-\mu(M)+1$ for all $z\neq 0$. The reason for this is that if $\psi_{t}(M_{t}z)$ was to be equal to zero more than $\mu(M)-1$ times, then $z$ would be necessarily equal to zero. As a result we get

[TABLE]

Hence by applying Theorem 3, $|\mathcal{I}^{c}(Y,z_{0})|/(T-\mu(M)+1)<1/2$ is a sufficient condition for exact recovery by the $\ell_{0}$-norm based estimator. ∎

Remark 2**.**

Assume that $\psi_{t}$ is defined to be the counting norm, i.e.,

[TABLE]

Then $\psi_{t}$ has a separable structure as illustrated in (37). Consider then defining, still under the observability assumption, an entrywise version of the concentration ratio by

[TABLE]

where $M_{ti}$ refers to the $i$-th row of $M_{t}$. Further, let

[TABLE]

with $\tilde{M}_{I}\in\mathbb{R}^{|I|\times n}$ denoting the matrix obtained by stacking the row vectors $\left\{M_{ti}:(t,i)\in I\right\}$. Then a result similar to Corollary 4 is obtainable: if the number the measurements corrupted by a nonzero error (among the $n_{y}T$ available) is strictly less than $(n_{y}T-\tilde{\mu}(M)+1)/2$, then the estimator $\mathcal{E}^{\circ}$ expressed in (43) (with $\psi_{t}$ being the $\ell_{0}$ norm as in (55)) recovers exactly the true state.

Remark 3**.**

Under the condition of Remark 2, if we consider the scenario where only a set of $k<n_{y}$ sensors may be compromised by attackers, then exact recovery is achieved if

[TABLE]

Taking into consideration the fact that $\tilde{\mu}(M)-1<T$, it can then be seen that (57) is equivalent to $k\leq\left\lceil n_{y}/2-1\right\rceil$ where the notation $\left\lceil r\right\rceil$, for $r\in\mathbb{R}$, refers to the smallest integer larger or equal to $r$. To sum up, when the $\psi_{t}$ are defined as in (55), the estimator (43) is able to return the true state matrix even when $\left\lceil n_{y}/2-1\right\rceil$ sensors get faulty over the entire observation horizon. This is reminiscent of a result stated in [10] which therefore appears to be a consequence of Theorem 3.

V-B Stability of the class of estimators $\mathcal{E}^{\circ}$ with respect to dense noise

We have argued that the class of estimators $\mathcal{E}^{\circ}$ in (43) is able to obtain exactly the true state matrix when there is no dense noises $(w_{t},v_{t})$ in the system equations and only the sparse noise $\left\{s_{t}\right\}$ is active. The question we ask now is whether this set of estimators can, in addition to sparse noise, handle dense process and output noises and to what extent this is possible. The starting point of our reflection is the observation that the dynamical system defined by

[TABLE]

produces the same output as system (1). Here, $\tilde{v}_{t}=C_{t}\tilde{w}_{t}$, with $\tilde{w}_{t}=\sum_{k=0}^{t-1}A_{t-1}\cdots A_{k+1}w_{k}$ a definition which uses the convention that the product $A_{t-1}\cdots A_{k+1}=I$ if $k=t-1$. Then the idea is to apply the estimator $\mathcal{E}^{\circ}$ to (58) by neglecting the dense component $(v_{t}+\tilde{v}_{t})$ of the output equation. To state the resilience result for $\mathcal{E}^{\circ}$, consider for a given $\varepsilon\geq 0$, a partition $(\tilde{\mathbb{T}}_{\varepsilon},\tilde{\mathbb{T}}_{\varepsilon}^{c})$ of $\mathbb{T}$ defined as in (22) with $f_{t}$ replaced by $\tilde{f}_{t}\triangleq s_{t}+(v_{t}+\tilde{v}_{t})=f_{t}+\tilde{v}_{t}$.

Theorem 4**.**

*Consider the estimator (43) for the system (1). Assume that the loss functions $\{\psi_{t}\}$ involved in (44) are defined by (7) in which $\psi$ is assumed to satisfy (P1)–(P5) with constant $\gamma_{\psi}=1$. Let $\varepsilon\geq 0$ and set $r=|\tilde{\mathbb{T}}_{\varepsilon}^{c}|$. Denote with $N$ a norm on $\mathbb{R}^{n\times T}$ defined by $N(Z)=\max_{t\in\mathbb{T}}\left\|z_{t}\right\|$ with $z_{t}$ being the $t$-th column of $Z$ and $\left\|\cdot\right\|$ being a norm on $\mathbb{R}^{n}$.

If the system (1) is observable on $\interval{0}{T-1}$ and $\nu_{r}(M)<1/2$, then there exists a $\mathcal{K}_{\infty}$ function $\alpha$ such that for all norm $\lVert\cdot\rVert$ on $\mathbb{R}^{n\times T}$,*

[TABLE]

where $R_{\Sigma}$ is some constant depending on the system $\Sigma$ and

[TABLE]

with $\tilde{f}_{t}=f_{t}+\tilde{v}_{t}$, and $D_{1}=\min_{\left\|z\right\|=1}\Psi_{\mathbb{T}}^{\circ}(z)$.

Proof.

Let $\hat{x}_{0}\in\operatorname*{arg\,min}_{z\in\mathbb{R}^{n}}V_{\Sigma}^{\circ}(Y,z)$. We first provide a bound on the error $e_{0}=\hat{x}_{0}-x_{0}$ with $x_{0}$ denoting the true initial state of system (1). By exploiting the fact that $V_{\Sigma}^{\circ}(Y,\hat{x_{0}})\leq V_{\Sigma}^{\circ}(Y,x_{0})$ and noting that $y_{t}=M_{t}x_{0}+\tilde{f}_{t}$, we reach the inequality

[TABLE]

By then reasoning quite similarly as in the proof of Theorem 1, we get

[TABLE]

which, by exploiting (47) and the assumption that $\nu_{r}(M)<1/2$, leads to

[TABLE]

Applying now Lemma 2, we conclude that for any norm $\lVert\cdot\rVert$ on $\mathbb{R}^{n}$, there exists a $\mathcal{K}_{\infty}$ function $\alpha$ such that $\left\|e_{0}\right\|\leq\alpha^{-1}(\rho)$. Now by observing that for any $\hat{X}\in\mathcal{E}^{\circ}(Y)$,

[TABLE]

the result follows by posing333We use here the convention that $A_{t-1}\cdots A_{0}=I$ if $t=0$. $R_{\Sigma}=\max_{t\in\mathbb{T}}\left\|A_{t-1}\cdots A_{0}\right\|_{\text{ind}}$ with $\left\|\cdot\right\|_{\text{ind}}$ being the matrix norm induced by the vector norm $\left\|\cdot\right\|$ on $\mathbb{R}^{n}$. ∎

The interest in Theorem 4 is that it provides a condition of resilience for the estimator $\mathcal{E}^{\circ}$ which can be checked numerically as will be discussed in the next section.

VI On the numerical evaluation of the resilience conditions

The analysis results presented in Sections IV and V rely on some functions (resilience index, concentration ratio, …) which characterize quantitatively some properties of the system being observed. A question we ask now is whether it would be possible to evaluate numerically these measures. In effect, computing the $r$-resilience index in (21) would help testing for example the resilience condition in Theorem 1. Similarly, evaluating the concentration ratio $\nu_{r}(M)$ introduced in (47) is the way to assess whether a given estimator is able to return the true state of a given system if we make an hypothesis on the number of potential nonzero errors in the measurements.

Unfortunately, obtaining numerically the numbers $p_{r}$, $\tilde{p}_{r}$ or $\nu_{r}$ require solving some hard nonconvex and combinatorial optimization problems. This is indeed a common characteristic of the concepts which are usually used to assess resilience; for example, the popular Restricted Isometry Property (RIP) constant  [6] is comparatively as hard to evaluate. We note however that when the dimension of the state is small enough, $\nu_{r}(M)$ can be exactly computed by taking inspiration from a method presented in [23] even though at the price of a huge (but affordable) computational cost. Alternatively, a cheaper overestimation can be obtained by means of convex optimization as suggested in [2]. The next lemma provides such an overestimate for $\nu_{r}(M)$.

Lemma 5** (An estimate of $\nu_{r}$).**

Assuming all quantities are well-defined (see the conditions in Definitions 3 and 4), the following statements hold:

• (a)

$\nu_{r}\leq p_{r}$**

• (b)

If $\mu(M)\leq T-1$ then

[TABLE]

where

[TABLE]

In (61), the $\lambda_{tk}$ denote the entries of the vector $\lambda_{t}\in\mathbb{R}^{T}$ and $\left\{V_{t}\right\}$ refers to the sequence of nonsingular weighting matrices involved in (7).

The proof of statement (a) is straightforward by noticing that (47) follows from (21) by constraining the variable $Z$ to be in $\mathcal{Z}_{\Sigma}$. As to the proof of statement (b), it follows a similar reasoning as the proof of Theorem 2 in [2].

The interest of this lemma is twofold. First it suggests that the resilience condition of $\mathcal{E}^{\circ}$ is weaker (in the sense that it is easier to achieve) than that of $\mathcal{E}$. Second, it provides an upper bound on $\nu_{r}(M)$ which can be computed by solving a convex optimization problem (see Eqs (60)-(61)). More specifically, given $\nu^{o}$ in (61), checking numerically whether $|\mathbb{T}_{\varepsilon}^{c}|<1/2(1+1/\nu^{o})$ provides a sufficient condition for $\nu_{r}(M)<1/2$ and so, for the resilience of the estimator (43). In a similar spirit as in Lemma 5, we now show that the parameter $\tilde{p}_{r}$ defined in (40) can also be overestimated via convex optimization if the loss functions $\phi_{t}$ and $\psi_{t}$ in (6)–(7) are both taken to be norms.

Lemma 6** (An estimate of $\tilde{p}_{r}$).**

Consider the resilience parameter $\tilde{p}_{r}$ defined in (40) where we assume that $\psi_{ti}(e)=|e|$ for all $(t,i,e)\in\mathbb{T}\times\mathbb{S}\times\mathbb{R}$ and $\phi_{t}$ is an arbitrary norm. Then

[TABLE]

where

[TABLE]

Proof.

(See Appendix -D) ∎

Note, under the assumptions of Lemma 6, that $\inf_{\begin{subarray}{c}Z\in\mathbb{R}^{n\times T}\end{subarray}}\left\{H_{\Sigma}(Z):c_{ti}^{\top}z_{\tau}=1\right\}$ is a convex optimization problem for any given $(t,i)$. Hence, solving for $b_{1}$ in (63) requires solving $Tn_{y}$ convex problems and picking the smallest value among all. The interest of the lemma is that it provides an overestimate of $\tilde{p}_{r}$ which is numerically computable. Based on the so obtained overestimate of $\tilde{p}_{r}$, we see from Theorem 1 that the estimator (8) is resilient to $r$ outliers if $r<b_{1}/2$. Moreover, we can deduce an underestimate of the number of outliers that the estimator (8) is able to handle as $r_{\max}=\max\left\{r:r<b_{1}/2\right\}$.

As a last remark in this paragraph, let us observe that Lemma 6 is also applicable to overestimate $p_{r}$ defined in (21) in the case of a single-output system, i.e., when $n_{y}=1$.

VII Simulation Results

In this part we will illustrate numerically the resilience properties of the proposed class of estimators. For this purpose, we consider for simplicity, an example of linear time-invariant system in the form (1). We select a single-input single-output example where the pair $(A,C)$ is given by

[TABLE]

We instantiate the loss functions in (5) as follows: For all $t$ in $\mathbb{T}$ and for all $(z,e)\in\mathbb{R}^{n}\times\mathbb{R}^{n_{y}}$, $\phi_{t}(z)=\psi(W_{t}z)$ and $\psi_{t}(e)=\psi(V_{t}e)$ where the weighting matrices $W_{t}$ and $V_{t}$ and the functions $\phi$ and $\psi$ will be specified below for each experiment.

VII-A Numerical certificate of exact recoverability

Suppose in this section that the process noise $w_{t}$ and the dense component $v_{t}$ of $f_{t}$ (see Eq. (2)) are both identically equal to zero. We then focus on testing the exact recoverability property of the estimator (43) in the presence only of the sparse noise $\left\{s_{t}\right\}$. The times of occurrence of the nonzeros values in the sequence $\left\{s_{t}\right\}$ are picked at random. As to its values there are also randomly generated from a zero-mean normal distribution with variance $100^{2}$. Given $T=100$ output measurements and the system matrices in (64), the estimator $\mathcal{E}^{\circ}$ is implemented by directly solving the optimization problem defined in (43) through the CVX interface [12]. Note that the implementation of the estimator (43)-(44) requires computing the matrices $M_{t}$ expressed in (45), which take the form $CA^{t}$ in the LTI case. A problem that may occur however is that if $A$ is Schur stable as is the case here (or unstable), taking successive powers of $A$ produces matrices $M_{t}$ which might not be of the same order of magnitude. To preserve the contribution of each term of (44), we introduce special weighting matrices $\left\{V_{t}\right\}$ (in the loss function $\psi$ selected as the $\ell_{1}$ norm) to normalize the rows of these matrices so that they all have unit $2$-norm. $V_{t}$ is therefore selected to be a diagonal matrix of the form $V_{t}=\operatorname*{diag}(V_{t1},\cdots,V_{tn_{y}})$, where

[TABLE]

Here, $c_{i}^{\top}$, $i=1,\ldots,n_{y}$, denote the $i$-th row of the matrix $C$. Indeed the effect of the weighting function in (44) is equivalent to changing $y_{t}$ and $M_{t}$ respectively to $\tilde{y}_{t}=V_{t}y_{t}$ and $\tilde{M}_{t}=V_{t}CA^{t}$. Posing M=\big{\{}\tilde{M}_{t}\big{\}}, it can be checked using the methods discussed in Section VI (See Eq. (60)) that at least $r_{\max}=30$ erroneous data (out of $T=100$ measurements) can be accommodated by the estimator while still returning exactly the true state.

To investigate empirical performance, we consider different ratios $|\Lambda_{0}^{c}|/T$ of nonzero values in the sequence $\left\{s_{t}\right\}$. For each fixed proportion of nonzero values, we run the estimator over $100$ different realizations of the output measurements. The results, depicted in Figure 1, tend to show that the estimator can still find the true state even for proportions of gross errors as large as $60\%$.

VII-B Performances in the presence of dense noise

We consider now the more realistic scenario where the process noise $\left\{w_{t}\right\}$ and the measurement noise $\left\{v_{t}\right\}$ are nonzero. We further assume them to be bounded, white and uniformly distributed. For the numerical experiments these signals are sampled from an interval of the form $\interval{-a}{a}$. For comparison purpose, we conduct the estimation with several estimators:

• •

an instance of the estimator (8), denoted $\mathcal{E}_{\ell_{2}^{2},\ell_{1}}$ in the sequel, in which the loss function $\phi_{t}$ is quadratic and $\psi_{t}$ is the $\ell_{1}$-norm and $\lambda=1000$ (see (5) and (6)-(7))

• •

an instance of the estimator (8) denoted $\mathcal{E}_{\ell_{1},\ell_{1}}$ in which both loss functions $\phi_{t}$ and $\psi_{t}$ are the $\ell_{1}$-norm with $\lambda=10$

• •

the estimator $\mathcal{E}^{\circ}$ defined in (43)

In addition we implement oracle versions444By oracle version of an estimator, we refer here to an implementation of this estimator which is aware of the sparse noise sequence $\left\{s_{t}\right\}$. of $\mathcal{E}_{\ell_{2}^{2},\ell_{1}}$ and of $\mathcal{E}_{\ell_{2}^{2},\ell_{2}^{2}}$ (the latter corresponding to an instance of $\mathcal{E}$ where both $\phi_{t}$ and $\psi_{t}$ are instantiated as quadratic functions).

Experiment 1: Resilience test

Keeping the level of both dense noises (i.e., $w_{t}$ and $v_{t}$) fixed with amplitude $a=0.03$ for the entries of the former and $a=0.1$ for the latter (yielding a Signal to Noise Ratio (SNR) of about 30 dB in each case), we apply the estimators $\mathcal{E}_{\ell_{2}^{2},\ell_{1}}$, $\mathcal{E}_{\ell_{1},\ell_{1}}$ and $\mathcal{E}^{\circ}$ as defined above) to $100$ different realizations of the output data and we compute the average of the corresponding relative estimation errors. This process is repeated for different fractions of nonzeros in the sparse noise $\left\{s_{t}\right\}$ ranging from [math] to $0.8$. The estimates obtained by these estimators are displayed in Figure 2 in log scale. For the sake of comparison, we also display the oracle estimates given by $\mathcal{E}_{\ell_{2}^{2},\ell_{1}}$ and those obtained by a standard least squares estimator $\mathcal{E}_{\ell_{2}^{2},\ell_{2}^{2}}$ (i.e. with $\phi_{t}$ and $\psi_{t}$ taken to be both quadratic in (5)). By oracle of an estimator, we mean here a version of that estimator which is aware of the true values of the sparse noise sequence $\left\{s_{t}\right\}$. The results tend to show that the estimator (8) remains stable until the (empirical) resilience condition is violated (an event that happens when the sparsity level for the sparse noise is around $60\%$). This is consistent with the resilience property characterized in Theorem 1 and the empirical observations made in Section VII-A according to which the estimator is insensitive to the sparse noise sequence $\left\{s_{t}\right\}$ as long as the number of nonzero values in it (whose magnitudes are possibly arbitrarily large) is less than a certain threshold determined by the properties of the system. While Lemma 6 provides an underestimate of the number of correctable outliers as $r_{\max}=8$ (out of $100$), we can observe that the empirical breakpoint in the current example seems to be indeed around $40\%$. The discrepancy between the two values is partly explained by the pessimism of the upper bound of $p_{r}$ proposed in Lemma 6.

Experiment 2: Stability with respect to dense noise

Now, we fix the sparsity level of the time sequence $\left\{s_{t}\right\}$ to $0.2$ and let the powers of the dense noise $\left\{(w_{t},v_{t})\right\}$ vary jointly from 5 dB to 100 dB in term of SNR. The estimates obtained by the estimators (8) and (43) with the choices of $\phi_{t}$ and $\psi_{t}$ agreed in the beginning of Section VII are displayed in Figure 3 in term of $\log_{10}$ of estimation errors. What this illustrates is that whenever the number of faulty data is reasonable (here $20\%$ of the available measurements), the estimator discussed in this section behaves almost in the same way as when there is no faulty data at all.

Experiment 3: Impact of the regularization parameter $\lambda$ in $\mathcal{E}$

To assess the influence of the regularization parameter $\lambda$ on the performance of the estimators $\mathcal{E}_{\ell_{2}^{2},\ell_{1}}$ and $\mathcal{E}_{\ell_{1},\ell_{1}}$, we fix the amplitudes of both dense noises $\left\{w_{t}\right\}$ and $\left\{v_{t}\right\}$ at the same level as in Experiment 1 (i.e. SNR equal to $30$dB and ratio of non-zeros entries in $\{s_{t}\}$ equal to 30$\%$). In this setting, we consider a set of values of $\lambda$ ranging from $10^{-3}$ to $10^{6}$. For each of these values we perform an estimation over a hundred realizations of the output data and compute the average of the corresponding relative estimation errors. The outcome of this test, depicted in Figure 4, tends to suggest that low values of $\lambda$ yield quite poor results. Conversely, when $\lambda$ tends towards infinity, both estimators’ performance measures saturate at the same value, namely $0.1$. It turns out that this limit value corresponds to the relative error obtained for $\mathcal{E}^{\circ}$ in Experiment 1 in the same configuration, hence suggesting that $\mathcal{E}$ tends indeed to $\mathcal{E}^{\circ}$ in behavior as $\lambda$ becomes large. Finally, it is interesting to observe that both performance curves exhibit minima located around $\lambda=750$ and $\lambda=10$ for $\mathcal{E}_{\ell_{2}^{2},\ell_{1}}$and $\mathcal{E}_{\ell_{1},\ell_{1}}$ respectively.

VIII Conclusion

In this paper, we have considered the problem of estimating the state of linear time-varying systems in the face of uncertainties modeled as process and measurement noises in the system equations. The measurement noise sequence assumes values of possibly arbitrarily large amplitude which occur intermittently in time and accross the available sensors. For this problem we have proposed a class of estimators based on the resolution of a family of parameterizable optimization problems. The discussed family is rich enough to include optimization-based estimators based on various loss functions which may be convex (e.g., $\ell_{p}$-norms) or nonconvex (e.g., $\ell_{p}$ quasi-norms or saturated functions), smooth or nonsmooth. In particular, we have proved a resilience property for the proposed class of state estimators, that is, the resulting estimation error is bounded by a bound which is independent of the extreme values of the measurement noise provided that the number of occurrences (over time and over the whole set of sensors) of such extreme values is limited. Note however that the estimators studied here operate in batch mode, that is, they apply to a finite collection of measurements. In future works we intend to investigate efficient and low cost adaptive versions of the proposed optimization framework. Another interesting research avenue would be to study the level of performance which is achievable if one uses the discussed framework as a method to detect bad data prior to a refinement with standard least squares estimation.

In this appendix, we provide some technical results used in the paper and the associated proofs.

-A A useful technical lemma

Lemma 7**.**

Let $\xi_{1},\xi_{2}:\mathbb{R}^{a\times b}\rightarrow\mathbb{R}_{\geq 0}$ be two functions which satisfy properties (P1)–(P3) and let $\ell:\mathbb{R}^{c\times d}\rightarrow\mathbb{R}^{a\times b}$ be an injective linear mapping. Then $\xi_{1}+\xi_{2}$ and $\xi_{1}\circ\ell$ verify (P1)–(P3). In addition, the following holds:

1. (j)

If $\xi_{1}$, $\xi_{2}$ verify (P4), then $\xi_{1}+\xi_{2}$ and $\xi_{1}\circ\ell$ verify (P4) . 2. (jj)

If $\xi_{1}$, $\xi_{2}$ verify (P5), then $\xi_{1}+\xi_{2}$ and $\xi_{1}\circ\ell$ verify (P5) .

The main point of interest of this lemma is that even if there are functions which satisfy properties (P4) and (P5) with different values of $q$ and $\gamma$, their sum still verifies those properties.

To prove Lemma 7, we will need the following result.

Lemma 8** (Minimum function of two $\mathcal{K}_{\infty}$ functions).**

If $q_{1}$ and $q_{2}$ are two $\mathcal{K}_{\infty}$ functions, then so is the function $q$ defined by

[TABLE]

Proof.

We have to prove that $q$ is continuous, strictly increasing and satisfies $q(0)=0$ and $\lim_{\lambda\rightarrow+\infty}q(\lambda)=+\infty$.

First of all, it is clear that $q(0)=0$. Also, continuity of $q$ is immediate from that of $q_{1}$ and $q_{2}$ by noting that $q=(q_{1}+q_{2}-|q_{1}-q_{2}|)/2$. To see the strict increasingness of $q$, consider $\lambda_{1}$ and $\lambda_{2}$ in $\mathbb{R}_{\geq 0}$ such that $\lambda_{1}<\lambda_{2}$. Then $q(\lambda_{1})\leq q_{1}(\lambda_{1})<q_{1}(\lambda_{2})$ and $q(\lambda_{1})\leq q_{2}(\lambda_{1})<q_{2}(\lambda_{2})$. It follows that $q(\lambda_{1})<\min_{i\in\{1,2\}}q_{i}(\lambda_{2})=q(\lambda_{2})$ and hence $q$ is strictly increasing. We now show that $q(\lambda)$ tends to infinity when $\lambda\rightarrow+\infty$. Let $M>0$ be an arbitrary positive number. Since $q_{1}$ and $q_{2}$ tend to infinity, there exist $\eta_{1}$ and $\eta_{2}$ such that $\lambda\geq\eta_{1}\Rightarrow q_{1}(\lambda)\geq M$ and $\lambda\geq\eta_{2}\Rightarrow q_{2}(\lambda)\geq M$. By taking $\eta=\max_{i\in\{1,2\}}\eta_{i}$, it holds that $q(\lambda)\geq M$ whenever $\lambda\geq\eta$, or equivalently that, $\lim_{\lambda\rightarrow+\infty}q(\lambda)=+\infty$. ∎

Proof of Lemma 7: The sum $\xi_{1}+\xi_{2}$ has clearly the properties (P1)–(P3) as a sum of continuous, even, positive definite functions. Moreover, the composition of a continuous, even, convex positive definite function with an injective linear mapping yields a continuous, even, positive definite function, so $\xi_{1}\circ\ell$ satisfies properties (P1)–(P3) too.

*Proof of (j): * Assume that $\xi_{1}$ and $\xi_{2}$ satisfy (P4) with $\mathcal{K}_{\infty}$ functions $q_{1}$ and $q_{2}$ respectively. For all $\lambda\neq 0$ and all $Z\in\mathbb{R}^{a\times b}$, (9) yields

[TABLE]

If we define $q$ so that for all $\lambda\in\mathbb{R}_{\geq 0}$, $q(\lambda)=\min_{i\in\{1,2\}}q_{i}(\lambda)$, then $q$ is a $\mathcal{K}_{\infty}$ function (see Lemma 8 above) such that for all $\lambda\neq 0$ and $Z\in\mathbb{R}^{a\times b}$,

[TABLE]

therefore $\xi_{1}+\xi_{2}$ verifies property (P4). Besides, for all $\lambda\neq 0$ and $Z$ in $\mathbb{R}^{c\times d}$,

[TABLE]

given the linearity of $\ell$. We can then conclude that $\xi_{1}\circ\ell$ also verifies property (P4).

Proof of (jj): Assume that $\xi_{1}$ and $\xi_{2}$ satisfy (P5) for $\gamma_{1}$ and $\gamma_{2}$ respectively. Let $\gamma=\min_{i\in\{1,2\}}\gamma_{i}$. Similarly to the first case, for all $Z_{1},Z_{2}$ in $\mathbb{R}^{a\times b}$ and $i$ in $\{1,2\}$, (10) yields

[TABLE]

which gives

[TABLE]

therefore $\xi_{1}+\xi_{2}$ satisfies property (P5). Moreover, for all $Z_{1}$ and $Z_{2}$ in $\mathbb{R}^{c\times d}$,

[TABLE]

so $\xi_{1}\circ\ell$ satisfies (P5) too.

∎

-B Proof of Lemma 3

(i) $\Rightarrow$ (ii): Assuming that the system is observable on the interval $\interval{0}{T-1}$, we need to prove that there exists a $\mathcal{K}_{\infty}$ function $q$ which verifies (14). The idea of the proof is to apply Lemma 7 to the function $F$ of $\mathbb{R}^{n\times T}$ defined by $F(Z)=V_{\Sigma}(0,Z)$ with $V_{\Sigma}$ defined as in (5). To begin with, we note that $F$ can be decomposed as $F=\xi\circ\ell$ where $\xi:\mathbb{R}^{n\times(T-1)}\times\mathbb{R}^{n\times T}\rightarrow\mathbb{R}_{\geq 0}$ is a loss function such that for $Z=\begin{pmatrix}z_{0}&\cdots&z_{T-2}\end{pmatrix}$ in $\mathbb{R}^{n\times(T-1)}$, $Y=\begin{pmatrix}y_{0}&\cdots&y_{T-1}\end{pmatrix}$ in $\mathbb{R}^{n_{y}\times T}$,

[TABLE]

and $\ell:\mathbb{R}^{n\times T}\rightarrow\mathbb{R}^{n\times(T-1)}\times\mathbb{R}^{n\times T}$ a linear mapping such that for all $Z=\begin{pmatrix}z_{0}&\cdots&z_{T-1}\end{pmatrix}$ in $\mathbb{R}^{n\times T}$,

[TABLE]

To apply Lemma 7 to $F$, we need to check that $F$ fulfills the properties (P1)–(P3). In virtue of the assumptions on $\phi_{t}$ and $\psi_{t}$ agreed in the statement of the lemma, the first two properties are obviously satisfied. The third will be satisfied if $\ell$ is injective, a propriety which we now check. Let $Z$ be such that $\ell(Z)=0$. Then

[TABLE]

An immediate consequence of (73)–(74) is that $\mathcal{O}_{0,T-1}z_{0}=0$ which yields $z_{0}=0$ because the system is observable on $\interval{0}{T-1}$. Therefore, thanks to the recursive relation (73), we can conclude that $Z=0$, and so, the linear mapping $\ell$ is injective.

We can therefore apply Lemma 7 to conclude that $F$ satisfy indeed  (P1)–(P4). Now, consider a matrix norm $\lVert\cdot\rVert_{\text{ind}}$ on $\mathbb{R}^{n\times T}$ induced by two vector norms $\lVert\cdot\rVert_{T}$ and $\lVert\cdot\rVert$ defined respectively on $\mathbb{R}^{T}$ and $\mathbb{R}^{n}$ in the sense that

[TABLE]

Applying Lemma 2 to $F$ with the so-defined induced norm, we infer that there exists $D>0$ defined as in (12) and a $\mathcal{K}_{\infty}$ function $q^{\prime}$, such that for all $Z$ in $\mathbb{R}^{n\times T}$,

[TABLE]

If we denote with $e_{1}$ the canonical vector of $\mathbb{R}^{T}$ with all entries equal to zero except the first one which is equal to $1$, then $Ze_{1}=z_{0}$. However, by definition of the induced norm, we know that $\lVert Ze_{1}\rVert/\lVert e_{1}\rVert_{T}\leq\lVert Z\rVert_{\text{ind}}$. Therefore, as $q^{\prime}$ is an increasing function, we get that $q^{\prime}(\lVert z_{0}\rVert/\lVert e_{1}\rVert_{T})\leq q^{\prime}(\lVert Z\rVert_{\text{ind}})$. By posing $q:\lambda\mapsto Dq^{\prime}(\lambda/\lVert e_{1}\rVert_{T})$, it is easy to see that $q$ is a $\mathcal{K}_{\infty}$ function so that for all $Z$ in $\mathbb{R}^{n\times T}$, $V(0,Z)=F(Z)\geq q(\lVert z_{0}\rVert)$. (ii) $\Rightarrow$ (i): Assume that there exists $q$ in $\mathcal{K}_{\infty}$ such that for all $Z\leavevmode\nobreak\ =\leavevmode\nobreak\ \begin{pmatrix}z_{0}&z_{1}&\ldots&z_{T-1}\end{pmatrix}$ in $\mathbb{R}^{n\times T}$ such that (14) holds. We want to prove that the matrix $\mathcal{O}_{0,T-1}$ defined in (13) is of full column rank, which is equivalent to showing that for $z$ in $\mathbb{R}^{n}$, $\mathcal{O}_{0,T-1}z=0$ implies $z=0$. For all $z\in\mathbb{R}^{n}$, construct a sequence $Z^{*}=\begin{pmatrix}z_{0}^{*}&\cdots&z_{T-1}^{*}\end{pmatrix}$ as follows: $z_{0}^{*}=z$ and $z_{t+1}^{*}=A_{t}z_{t}^{*}$ for all $t\in\{0,\ldots,T-2\}$. Since the inequality (14) is supposed to be true for any sequence, so it is for the particular sequence $\{z_{t}^{*}\}$ defined above. Applying this inequality to $Z^{*}$ yields

[TABLE]

Now, observe that if $\mathcal{O}_{0,T-1}z=0$, then it follows from the recursive relation $z_{t+1}^{*}=A_{t}z_{t}^{*}$ that for all $t$ in $\{0,\ldots,T-1\}$, $C_{t}z_{t}^{*}=0$. Injecting this in (76) imposes that $q(\lVert z_{0}^{*}\rVert)\leq 0$ which necessarily implies that $z=0$ as $q$ is a $\mathcal{K}_{\infty}$ function. Therefore, the matrix $\mathcal{O}_{0,T-1}$ is injective and the system is observable on the interval $\interval{0}{T-1}$.

-C Technical results for proving Corollary 2

This section contains some technical steps of the proof of Corollary 2.

Lemma 9**.**

If $\ell:\mathbb{R}^{n_{y}}\rightarrow\mathbb{R}_{\geq 0}$ satisfies (P1)–(P3) and (P5), then so does the function $\psi$ defined by $\psi(y)=1-e^{-\ell(y)}$. Moreover if $\ell$ fulfills (P4), then $\psi$ satisfies the same property but with a function $q$ in $\mathcal{K}_{\operatorname*{sat},a}$ for $a=1$.

Proof.

It is straightforward to check that $\psi$ obeys (P1)-(P3). By assumption, $\ell$ obeys (P5). Denote therefore the associated constant with $\gamma_{\ell}$ (which, by (10), is necessarily less than or equal to $1$). To see then that (P5) is also satisfied by $\psi$, we just need to check that

[TABLE]

with $\bar{\gamma}_{\ell}=\gamma_{\ell}^{-1}\geq 1$, which is equivalent to

[TABLE]

Noting that $\ell(a+b)\leq\bar{\gamma}_{\ell}\ell(a)+\bar{\gamma}_{\ell}\ell(b)$, we have $-e^{-\ell(a+b)}\leq-e^{-\bar{\gamma}_{\ell}\ell(a)-\bar{\gamma}_{\ell}\ell(b)}$. From this it follows that for (77) to hold, it is enough that

[TABLE]

Posing $\alpha=e^{-\ell(a)}$ and $\beta=e^{-\ell(b)}$, it suffices that

[TABLE]

which can indeed be checked to be true by applying the identity $1+\bar{\gamma}_{\ell}\alpha-\bar{\gamma}_{\ell}\leq\alpha^{\bar{\gamma}_{\ell}}$, see e.g., [3, Fact 1.9.2]. In effect, it follows from this identity that

[TABLE]

In conclusion, (77) holds and therefore $\psi$ satisfies (P5).

It remains now to check (P4). This follows directly from Lemma 10 below, from which we know that $\psi(y)\geq q^{\star}(1/\lambda)\psi(\lambda y)$ with $q^{\star}$ is a saturated function in $\mathcal{K}_{\operatorname*{sat},1}$. ∎

Lemma 10**.**

Let $\ell:\mathbb{R}^{n_{y}}\rightarrow\mathbb{R}_{\geq 0}$ be a function satisfying properties (P1)–(P2) and (P4). In particular, assume that property (P4) is satisfied by $\ell$ with a $\mathcal{K}_{\infty}$ function $q$ such that (9) is an equality relation. Let

[TABLE]

for $\lambda\neq 0$ and $y\neq 0$. Then the function $q^{\star}:\mathbb{R}_{\geq 0}\rightarrow\interval{0}{1}$ defined by $q^{\star}(\lambda)=\inf_{y\neq 0}g(y,\lambda)$ for $\lambda>0$ and $q^{\star}(0)=0$, is well-defined, continuous and strictly increasing on $\interval{0}{1}$. Moreover we have

[TABLE]

Proof.

Since $g$ is positive on its domain (hence lower-bounded), the defining infimum of $q^{\star}$ is well-defined. Pose $a=e^{-\ell(y)}$. Then by using the continuity property of $\ell$ and its radial unboundedness (see Lemma 2), we see that the range of $a$ when $y$ lives in $\mathbb{R}^{n_{y}}\setminus\{0\}$ is $\interval[open]{0}{1}$. From the assumptions of the lemma, $\ell(y/\lambda)=q(1/\lambda)\ell(y)$ for all $y$ and all $\lambda>0$ and so, $q(1)=1$ and $e^{-\ell(y/\lambda)}=a^{q(1/\lambda)}$. For all $\lambda>0$ we can write

[TABLE]

with $q(1/\lambda)\geq 1$ for $0<\lambda\leq 1$ and $q(1/\lambda)<1$ for $\lambda>1$. We therefore obtain

[TABLE]

The so obtained $q^{\star}$ is clearly continuous wherever it is well defined. Moreover, since $\lim_{\lambda\rightarrow 0}q^{\star}(\lambda)=q^{\star}(0)=0$, we conclude that $q^{\star}$ is continuous on its entire domain. From the properties of $q$, we deduce that $q^{\star}$ is strictly increasing on $\interval{0}{1}$. Lastly, we observe that the inequality in the statement of the lemma is a direct consequence of the definition of $q^{\star}$. ∎

-D Proof of Lemma 6

The starting point of the proof is the observation that for every integer $r$ in $\left\{1,\ldots,T\right\}$, $\tilde{p}_{r}\leq r\tilde{p}_{1}$. Hence it suffices to show that $\tilde{p}_{1}=1/b_{1}$ and is as expressed in (63). Recall that by definition,

[TABLE]

Without loss of generality, assume that $c_{ti}^{\top}\neq 0$ for all $(t,i)\in\mathbb{T}\times\mathbb{S}$ Then for any $(t,i)$,

[TABLE]

where

[TABLE]

Recalling that $H_{\Sigma}(Z)$ is a norm under the conditions of the lemma, the second equality in the expression of $\beta_{ti}$ above follows from the (strict) homogeneity property of norms. As to the last equality, it follows from the fact that $c_{ti}^{\top}z_{t}$ is a scalar which induces the possibility to replace the constraint $|c_{ti}^{\top}z_{t}|=1$ indifferently either by $c_{ti}^{\top}z_{t}=1$ or by $c_{ti}^{\top}z_{t}=-1$.

Now by invoking the definition of $\tilde{p}_{1}$, it can be seen that

[TABLE]

∎