Alexandre Kircher
email: alexandre.kircher@ec-lyon.fr

Laurent Bako
email: laurent.bako@ec-lyon.fr

Eric Blanco
email: eric.blanco@ec-lyon.fr

Mohamed Benallouch
email: mohamed.benallouch@ecam.fr

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

Keywords: Secure state estimation, resilient estimators, optimal estimation, Cyber-physical systems. 1

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 proposed general estimation framework. 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.

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) [START_REF] Cardenas | Secure control: Towards survivable cyber-physical systems[END_REF]. 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 [START_REF] Huber | Robust Statistics[END_REF], hybrid system identification [START_REF] Bako | Identification of switched linear systems via sparse optimization[END_REF], 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 [START_REF] Mishra | Secure state estimation against sensor attacks in the presence of noise[END_REF], [START_REF] Pasqualetti | Attack detection and identification in cyber-physical systems[END_REF]. 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 [START_REF] Sinopoli | Kalman filtering with intermittent observations[END_REF], [START_REF] Mo | Secure estimation in the presence of integrity attacks[END_REF], [START_REF] Ren | Secure state estimation with byzantine sensors: A probabilistic approach[END_REF]. Another category of approaches, which are inspired by some recent results in compressive sampling [START_REF] Candès | An introduction to compressive sampling[END_REF], [START_REF] Foucart | A mathematical introduction to compressive sensing[END_REF], 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 [START_REF] Fawzi | Secure estimation and control for cyber-physical systems under adversarial attacks[END_REF]. 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 [START_REF] Fawzi | Secure estimation and control for cyber-physical systems under adversarial attacks[END_REF] 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, [START_REF] Shoukry | Event-triggered state observers for sparse sensor noise/attacks[END_REF] suggests a reformulation which is argued to reduce computational cost by using the concept of event-triggered update ; [START_REF] Pajic | Attack-resilient state estimation for noisy dynamical systems[END_REF] considers an observation model which includes dense noise along with the sparse attack signal. In [START_REF] Chang | Secure estimation based kalman filter for cyber-physical systems against sensor attacks[END_REF], the assumption of a fixed number of attacked sensors is relaxed. Finally, the recent paper [START_REF] Han | Convex optimization based state estimation against sparse integrity attacks[END_REF] 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 [START_REF] Bako | On a class of optimization-based robust estimators[END_REF], [START_REF] Candès | Highly robust error correction by convex programming[END_REF].

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 [START_REF] Kircher | Analysis of resilience for a state estimator for time-discrete linear systems[END_REF], 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 [START_REF] Han | Convex optimization based state estimation against sparse integrity attacks[END_REF] 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 optimizationbased 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. R ≥0 (respectively R>0) is the set of nonnegative (respectively positive) reals. R * designates the set of real numbers excluding zero. We note R a the set of (column) vectors with a real elements and R a×b , the set of real matrices with a rows and b columns. If M ∈ R a×b , then M will designate the transposed matrix of M . I will refer to the (square) identity matrix of appropriate dimension. The notation • will denote a norm over a given set (which will be specified when necessary). • p denotes the p norm (for p ≥ 1) or the p quasi-norm (for 0 < p < 1)

defined for z = z1 • • • za in R a by z p = (|z1| p + • • • + |za| p ) 1/p
. The limit of this when p ↓ 0 gives the so-called 0-norm • 0 of z, i.e., the number of nonzero entries in z. Its limit when p → +∞ gives the infinity norm denoted z ∞ and returning the maximum value of the |zi|. For x ∈ R, e x refers to the exponential function applied to x. If S is a set, then P(S) is the power set of S. If S is a finite set, the notation |S| refers to the cardinality of S. K∞ functions [START_REF] Kellett | A compendium of comparison function results[END_REF]. We name K∞ the set of functions f : R ≥0 → R ≥0 which are continuous, zero at zero, strictly increasing and satisfy lim λ→+∞ f (λ) = +∞. If f ∈ K∞, then it admits an inverse, denoted here f -1 , which also lies in K∞. Similarly, we use the notation Ksat,a to denote the set of saturated functions f : R ≥0 → R ≥0 which are continuous, zero at zero, strictly increasing on [0, a] and such that f (λ) = f (a) for all λ ≥ a. Supremum. Given a function f over R a and a subset S of R a , the notation sup z∈S f (z) < b, with b ∈ R, will mean that for all z in S, f (z) < b. This notation includes the case where the supremum is b but is not attained by any element of S.

II. THE RESILIENT ESTIMATION PROBLEM

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

Σ : xt+1 = Atxt + wt yt = Ctxt + ft (1)
where xt ∈ R n is the state vector of the system at time t and yt ∈ R ny is the output vector at time t; {At} and {Ct} are families of matrices with appropriate dimensions; {wt} is an unobserved bounded noise sequence. As to {ft}, it is regarded here as an (unobserved) arbitrary noise sequence affecting the measurements.

For clarity of the exposition, it may be convenient to view ft as a combination of two types of sequences: a bounded sequence {vt} and a sparse sequence {st} (this decomposition is indeed always possible for an arbitrary noise signal). Hence, we may write

ft = vt + st, (2)
where vt is a sensor noise of moderate amplitude and st 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 {st} may account for adversarial attacks in the same spirit as in [START_REF] Fawzi | Secure estimation and control for cyber-physical systems under adversarial attacks[END_REF], [START_REF] Han | Convex optimization based state estimation against sparse integrity attacks[END_REF], intermittent sensor faults, or data losses, in particular when a communication network is involved in the data acquisitiontransmission chain. In the sequel, we may also refer to {wt} and {vt} in ( 1) and ( 2) as dense noises and to the largest elements of {st} as outliers.

For the sake of simplicity, the sparse (and potentially arbitrary large) noise is assumed here to affect only the measurement equation of the system (1). 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 x0, . . . , xT -1 of the system (1) on a time period T = {0, . . . , T -1} given T measurements y0, ..., yT -1 of the system output. We shall seek an accurate estimate of the state despite the uncertainties in the system equations (1) modeled by wt and ft 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 those instances of {ft} having the largest amplitudes. Such an estimator will be called resilient, see Definition 2 for a formal characterization of this property.

Denote with

XT

-1 = x0 x1 . . . xT -1 (3)
the actual state trajectory of the system Σ on a finite time horizon of length T . Similarly, we use the notation

YT -1 = y0 y1 • • • yT -1 (4)
to refer to the collection of measurements on the same time horizon.

The state estimation problem is approached here from an offline perspective, that is, 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 T = {0, . . . , T -1}. To simplify further the formulas, we also pose T = {0, . . . , T -2} while S = {1, . . . , ny} will be a set indexing the sensors (or the rows of the matrices Ct in (1)).

III. OPTIMIZATION-BASED APPROACH TO RESILIENT STATE ESTIMATION

A. The state estimator

In this section we present an optimization-based framework for solving the state estimation problem defined above. To define for-mally the proposed state estimator, let us first introduce the to-beminimized objective function. Given the matrices {(At, Ct)} of the system (1) and T output measurements

Y = y0 • • • yT -1 , we consider a performance function VΣ : R ny ×T × R n×T → R ≥0 defined by VΣ(Y, Z) = λ t∈T φt(zt+1 -Atzt) + t∈T ψt(yt -Ctzt) (5)
where Z = z0 • • • zT -1 ∈ R n×T is a hypothetical trajectory matrix with zi denoting the i-th column of Z; λ > 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Σ(Y, Z). {φt} and {ψt} are two families of positive functions (called here loss functions) defined on R n and R ny respectively. For the sake of simplicity, we will assume throughout the paper that for all t in T, φt and ψt can be expressed by

φt(z) = φ(Wtz) ∀z ∈ R n (6) ψt(e) = ψ(Vte) ∀e ∈ R ny , (7)
where φ : R n → R ≥0 and ψ : R ny → R ≥0 are two fixed loss functions and {Wt} and {Vt} are two families of nonsingular weighting matrices with appropriate dimensions.

Definition 1. Given a system Σ such as the one in (1) and given an output measurement matrix Y ∈ R ny ×T , we define a state estimator to be a set-valued map E : R ny ×T → P(R n×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

E(Y ) = arg min Z∈R n×T VΣ(Y, Z). (8)
As such the estimator E is well-defined if for any fixed Y , VΣ(Y, Z) admits a non empty minimizing set, that is, if there exists at least one Z such that VΣ(Y, Z) ≥ VΣ(Y, Z ) for all Z ∈ R n×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 φ and ψ entering in the definition of the objective function VΣ.

B. Well-definedness of the estimator

Let us start by stating the properties required for the loss functions involved in the definition of VΣ. 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 ξ : R a×b → R ≥0 which will be required to satisfy a subset (depending of the specific usage) of the following properties: (P1) Positive definiteness: ξ(0) = 0 and ξ(Z) > 0 for all Z = 0 (P2) Continuity: ξ is continuous (P3) Symmetry: ξ(-Z) = ξ(Z) for all Z ∈ R a×b (P4) Generalized Homogeneity (GH): There exists a K∞ function q : R ≥0 → R ≥0 such that for all λ ∈ R * and for all Z ∈ R a×b ,

ξ(Z) ≥ q 1 |λ| ξ(λZ). (9)
(P5) Generalized Triangle Inequality (GTI): There exists a positive real number γ ξ such that for all Z1, Z2 in R a×b

ξ(Z1 -Z2) ≥ γ ξ ξ(Z1) -ξ(Z2). (10)
It can be usefully observed, for the future developments, that [START_REF] Fawzi | Secure estimation and control for cyber-physical systems under adversarial attacks[END_REF] can be equivalently written as ξ(Z1 + Z2) ≤ γ -1 ξ ξ(Z1) + γ -1 ξ ξ(Z2). Examples of loss functions. Note that norms on R a×b satisfy naturally the properties (P1)-(P5) with q : λ → λ and γ ξ = 1, hence yielding the classic homogeneity property and triangle inequality. It can also be checked that functions ξ of the form ξ(Z) = Z p with p > 0, fully qualify as loss functions in the sense that they fulfill all the properties (P1)-(P5). In this case, γ ξ in [START_REF] Fawzi | Secure estimation and control for cyber-physical systems under adversarial attacks[END_REF] can be taken equal to 2 1-1/p if 0 < p ≤ 1 and 2 1-p otherwise. Lastly we note that if : R a×b → R ≥0 satisfies (P1)-(P3) and (P5), then so does the function ξ defined by ξ(Z) = 1e -(Z) (see Lemma 9 in the appendix). Similarly, saturated functions of the form ξ(Z) = min( (Z), R0) for some R0 > 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 ξ : R a×b → R ≥0 is convex and satisfies property (P4) with a K∞ function q, then it also satisfies (P5) with γ ξ = 2q(1/2).

Observe that quadratic functions ξ : R a×b → R ≥0 of the form ξ(Z) = Tr(Z QZ) with Q ∈ R a×a being a positive definite matrix and Tr referring to the trace of a matrix, satisfy properties (P1)-(P4) with a K∞ function q : λ → λ 2 . Since such functions are convex, it follows from Lemma 1 above that they also verify (P5) for γ ξ = 2q(1/2) = 1/2. Remark 1. In virtue of (6)-( 7), the families of functions {φt} and {ψt} satisfy (P1)-(P5) whenever φ and ψ satisfy (P1)-(P5).

We now recall from [START_REF] Kircher | Analysis of resilience for a state estimator for time-discrete linear systems[END_REF] a technical lemma which will play a fundamental role in analyzing the properties of the estimator [START_REF] Cardenas | Secure control: Towards survivable cyber-physical systems[END_REF]. In particular, our proof of well-definedness relies on this lemma.

Lemma 2 (Lower Bound of a loss function). Let ξ : R a×b → R ≥0 be a function which has properties (P1)-(P2) and (P4) with a K∞ function q. Then, for all norm

• on R a×b , ξ(Z) ≥ Dq( Z ) ∀Z ∈ R a×b (11)
where

D = min Z =1 ξ(Z) > 0. (12)
Proof: We start by observing that the unit hypersphere S = Z ∈ R a×b : Z = 1 is a compact set in the topology induced by the norm • . By the extreme value theorem, ξ being continuous, admits necessarily a minimum value on S, i.e., there is Z ∈ S such that ξ(Z) ≥ D ξ(Z ) > 0 for all Z ∈ S. For any nonzero

Z ∈ R a×b , Z Z ∈ S so that ξ( Z Z ) ≥ D.
On the other hand, by the relaxed homogeneity of ξ,

ξ(Z) ≥ q( Z )ξ( Z Z ) ≥ Dq( Z ).
Moreover, this inequality holds for Z = 0. It therefore holds true for any Z ∈ R a×b .

Proposition 1 (Well-definedness of the estimator). Let the loss functions φ and ψ in (6)- [START_REF] Candès | An introduction to compressive sampling[END_REF] satisfy properties (P1)-(P5) and assume that the LTV system (1) is observable on [0, T -1] in the sense that the observability matrix

O0,T -1 (C0) (C1A0) • • • (CT -1AT -2 . . . A1A0) (13
) has full column rank. Then the estimator (8) is well-defined, i.e., the objective function VΣ(Y, •) attains its minimum for any fixed Y .

Hence, the condition of the proposition guarantees that E(Y ) is non empty for all Y ∈ R ny ×T . Before proving this result, we first make the following observation.

Lemma 3 (Equivalent condition of Observability). Consider the objective function VΣ defined in [START_REF] Candès | Highly robust error correction by convex programming[END_REF] where {(φt, ψt)} are defined as in (6)-( 7) with φ and ψ satisfying (P1)-(P4). Then the following two statements are equivalent:

(i) The system is observable on the time interval [0, T -1].

(ii) There exists a K∞ function q such that for all

Z = z0 z1 . . . zT -1 in R n×T , VΣ(0, Z) ≥ q( z0 ) (14)
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 {φt} and {ψ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Σ(Y, •) is coercive (i.e., continuous and radially unbounded) for any given Y and then apply a result1 in [22, Thm 1.9] to conclude on the attainability of the infimum (which certainly exists since VΣ(Y, •) is a positive function). Clearly, VΣ(Y, •) is continuous as a consequence of φ and ψ being continuous by assumption (see property (P2)). We then just need to prove the radial unboundedness of VΣ(Y, •), i.e., lim Z →+∞ VΣ(Y, Z) = +∞ for an arbitrary norm • on the Zspace and for all fixed Y . Since ψ satisfies property (P5), there exists a constant γ ψ > 0 such that ψt(yt -Ctzt) ≥ γ ψ ψt(Ctzt)ψt(yt). Applying this property leads naturally to

VΣ(Y, Z) ≥ F (Z) - t∈T ψt(yt),
where

F (Z) = λ t∈T φt(zt+1 -Atzt) + γ ψ t∈T ψt(Ctzt). (15)
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 • on R n×T , there exists a K∞ function q such that F (Z) ≥ q( Z ).

Combining this with the inequality above, we obtain that

VΣ(Y, Z) ≥ q( Z ) - t∈T ψt(yt)
which implies the radial unboundedness of VΣ(Y, •) 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 φ and ψ ensure that E(Y ) is a non empty set for any given Y . We then call any member X = x0 x1 . . . xT -1 of E(Y ), an estimate of the state trajectory X of system (1) on the time interval T. In particular, xt is called an estimate of the state xt at time t ∈ T.

To conclude this section, note that the definition of the estimator in [START_REF] Cardenas | Secure control: Towards survivable cyber-physical systems[END_REF] does not require any convexity assumption on the objective function VΣ. 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 φ and ψ. In effect, the elements of 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Σ(Y, Z) is a convex function of Z [START_REF] Boyd | Convex optimization[END_REF], [START_REF] Grant | CVX: Matlab software for disciplined convex programming[END_REF]. Nevertheless it is fair to recognize that nonconvex optimization methods can be successfully implemented as well though with less

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 [START_REF] Cardenas | Secure control: Towards survivable cyber-physical systems[END_REF]. For this purpose, let Y = y 0 . . . , y T -1 denote the noise-free output matrix of (1), i.e., the output defined by y t = Ctx t , t = 0, . . . , T -1, with x t+1 = Atx t and x 0 = x0 (x0 being the true initial state of ( 1)). Let Fr, a subset of R ny ×T containing 0, denote a matrix of measurement noise components.

Definition 2 (Resilience of an estimator). The set-valued estimator E defined in (8) is called resilient against the set Fr of measurement noise if there exists a K∞ function g such that, when the process noise {wt} is zero, it holds that for any measurement noise matrix

F ∈ R ny ×T , X -X ≤ g inf Ω∈Fr d(F -Ω) ∀ X ∈ E(Y + F ) (16)
with • denoting some norm, d : R ny ×T → R ≥0 being a function subject to (P1)-(P5) and Y referring to the noise-free output as defined above. Hence infΩ∈F r d(F -Ω) denotes some pseudodistance from F to the set Fr.

Since 0 ∈ Fr, a consequence of property [START_REF] Kircher | Analysis of resilience for a state estimator for time-discrete linear systems[END_REF] is that E(Y ) = {X}, which follows from [START_REF] Kircher | Analysis of resilience for a state estimator for time-discrete linear systems[END_REF] 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 Σ 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 Fr, that is, E(Y + F ) = {X} for all measurement noise F ∈ Fr. In case the matrix F does not lie in Fr but is situated at a bounded distance from Fr, then the estimator error will no longer be equal to zero but will just be bounded. Throughout this paper, we consider a set Fr defined as follows. For

F = f0 • • • fT -1 ∈ R ny ×T , let T c 0 (F ) = {t ∈ T : ψt(ft) > 0} and T0(F ) = {t ∈ T : ψt(ft) = 0}.
For r a positive integer, define Fr to be the set of matrices in R ny ×T having at most r nonzero columns, i.e.,

Fr = {F : |T c 0 (F )| ≤ r} . ( 17
)
For the need of making explicit the resilience property (in the results to be presented) with respect to the set Fr, we will need the following lemma.

Lemma 4. Consider the set Fr of measurement noise matrix defined in [START_REF] Mishra | Secure state estimation against sensor attacks in the presence of noise[END_REF] and select a (pseudo distance) function d : R ny ×T → R ≥0 defined by d(F ) = t∈T ψt(ft) with ψt a function defined as in [START_REF] Candès | An introduction to compressive sampling[END_REF] and having the properties (P1)-(P5). Then infΩ∈F r d(F -Ω) is equal to the sum of the Tr smallest terms in {ψt(ft) : t ∈ T}.

Proof: Let I c r (F ) denote the index set of the r largest entries of the vector ψ0(f0) . . . ψT -1(fT -1) and Ir(F ) denote the index set of its Tr smallest entries. Then, with the notation

Ω = ω0 • • • ωT -1 , inf Ω∈Fr d(F -Ω) = inf Ω∈Fr t∈T ψt(ωt -ft) = inf Ω∈Fr t∈Ir (F ) ψt(ωt -ft) + t∈I c r (F ) ψt(ωt -ft) = inf Ω∈Fr T 0 (Ω)=Ir (F ) t∈Ir (F ) ψt(ft) + t∈I c r (F ) ψt(ωt -ft) = t∈Ir (F ) ψt(ft)
where the notation T0(Ω) is defined in the lines preceding Eq. ( 17). The infimum is reached here for Ω ∈ Fr such that ωt = 0 ∀t ∈ Ir(F ) and ωt = ft ∀t ∈ I c r (F ). Hence infΩ∈F r d(F -Ω) is, as claimed, the sum of the Tr smallest values among {ψt(ft) : t ∈ T}.

We will also introduce in the sequel a notion of approximate resilience of E. This terminology refers to Definition 2 when the right hand side of ( 16) is modified as g infΩ∈F r d(F -Ω) + δ with δ some nonnegative real number.

B. Some notational conventions for the analysis

For convenience, let us introduce a few more notations. Let Φ :

R n×T → R ≥0 and Ψ T : R n×T → R ≥0 be defined by

Φ(Z) = t∈T φt(zt+1 -Atzt) (18)
Ψ T (Z) = t∈T ψt(Ctzt) (19)
We also introduce the partial cost function Ψ I defined for any I ⊂ T by Ψ I (Z) = t∈I ψt(Ctzt). We will assume throughout the paper that the loss functions φ and ψ 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 γ φ and γ ψ respectively. Finally, let us pose

HΣ(Z) = λγ φ Φ(Z) + γ ψ Ψ T (Z). (20)
We will organize the resilience analysis for the estimator (8) along two cases: first, the scenario where the gross error vector sequence {st} 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 ∈ R ny ×T the matrix formed from the sequence {st : t ∈ T}, 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).

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Σ and on the loss functions φ and ψ. Definition 3. Let r be a nonnegative integer. Assume that the system Σ in (1) is observable on [0, T -1]. We then define the r-Resilience index of the estimator E in (8) (when applied to Σ) to be the real number pr given by pr = sup

Z∈R n×T Z =0 sup I⊂T |I|=r Ψ I (Z) HΣ(Z) ( 21
)
where HΣ is as defined in [START_REF] Pasqualetti | Attack detection and identification in cyber-physical systems[END_REF]. The supremum is taken here over all nonzero Z in R n×T and over all subsets I of T with cardinality equal to r.

The index pr can be interpreted as a quantitative measure of the observability of the system Σ. The observability is needed here to ensure that the denominator HΣ(Z) of ( 21) is different from zero whenever Z = 0. Furthermore, it should be remarked that Ψ I (Z) ≤ HΣ(Z) for any I ⊂ T, which implies that the defining suprema of pr are well-defined. Note that pr is an increasing function of r and satisfies p0 = 0 and pT = 1. More discussions on the numerical evaluation of pr are deferred to Section VI.

In order to state the resilience result for the estimator (8) when applied to system Σ, let us introduce a last notation to be used in the analysis. Let ε ≥ 0 be a given number. For any admissible sequence {ft} t∈T in (1), we can split the time index set T into two disjoint label sets,

Tε = {t ∈ T : ψt(ft) ≤ ε} , (22)
indexing those ft which are upper bounded by ε and T c ε = {t ∈ T : ψt(ft) > ε} indexing those ft which are possibly unbounded. It is important to keep in mind that ε 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 {ψ(ft)}). For example, taking ε = 0 would be appropriate for analyzing the properties of the estimator when ft 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 Σ defined by [START_REF] Bako | Identification of switched linear systems via sparse optimization[END_REF] with output Y together with the state estimator [START_REF] Cardenas | Secure control: Towards survivable cyber-physical systems[END_REF] in which the loss functions φ and ψ are assumed to obey (P1)-(P5). Denote with γ φ and γ ψ the constants associated with the GTI (P5) and q φ and q ψ the K∞ functions associated with the GH (P4) for φ and ψ respectively. Let ε ≥ 0 and set r = |T c ε |.

If the system is observable on [0, T -1] and pr < 1/(1 + γ ψ ), then for any norm • on R n×T , X -X ≤ h 2βΣ(ε) D 1 -(1 + γ ψ )pr + δ(ε) ∀ X ∈ E(Y ) (23)
with X denoting the true state matrix from (1), D = min Z =1 HΣ(Z) > 0 and βΣ(ε), δ(ε) and h being defined by (28) Similarly, by making use of (1), observe that ψt(yt-Ct xt) = ψt(ft-Ctet). We now apply the GTI and symmetry of ψt in two different ways depending on whether t belongs to Tε or T c ε :

βΣ(ε) = λ t∈T φt(wt) + t∈Tε ψt(ft), (24)
δ(ε) = 1 -γ ψ D [1 -(1 + γ ψ )pr] t∈T c ε ψt(ft) (25)
h(α) = max q -1 φ (α), q -1 ψ (α) , α ∈ R ≥0 (26) Proof: Let X in E(Y ). By definition of E in (8), we have VΣ(Y, X) ≤ VΣ(Y,
∀t ∈ Tε, ψt(yt -Ct xt) ≥ γ ψ ψt(Ctet) -ψt(ft) ∀t ∈ T c ε , ψt(yt -Ct xt) ≥ γ ψ ψt(ft) -ψt(Ctet)
These inequalities imply that the second term on the left hand side of ( 27) is lower bounded as follows

t∈Tε [γ ψ ψt(Ctet) -ψt(ft)] + t∈T c ε [γ ψ ψt(ft) -ψt(Ctet)] ≤ t∈T ψt(yt -Ct xt) (29)
Combining ( 27), ( 28) and (29) gives

λγ φ t∈T φt(et+1 -Atet) + γ ψ t∈T ψt(Ctet) -(1 + γ ψ ) t∈T c ε ψt(Ctet) ≤ 2 λ t∈T φt(wt) + t∈Tε ψt(ft) + t∈T c ε (1 -γ ψ )ψt(ft)
which, by using ( 19), ( 20), [START_REF] Shoukry | Event-triggered state observers for sparse sensor noise/attacks[END_REF], can be written as

HΣ(E) -(1 + γ ψ )Ψ T c ε (E) ≤ 2βΣ(ε) + t∈T c ε (1 -γ ψ )ψt(ft) with E = e0 e1 • • • eT -1 .
As T c ε has r elements, applying the definition of pr gives

Ψ T c ε (E) ≤ prHΣ(E) (30)
By the assumption that pr < 1/(1 + γ ψ ) we have that 1 -(1 + γ ψ )pr > 0, and consequently, that

HΣ(E) ≤ 1 1 -(1 + γ ψ )pr 2βΣ(ε) + (1 -γ ψ ) t∈T c ε ψt(ft) (31)
Given that the system is observable on [0, T -1], it can be shown, thanks to Lemma 7 in the Appendix, that HΣ 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

• HΣ(E) ≥ Dq ( E ) (32)
with D defined by D = min Z =1 HΣ(Z) and q (α) = min{q φ (α), q ψ (α)}. Finally, the result follows by selecting h to be h = q -1 with q -1 denoting the inverse of q , 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 ψ satisfies the triangle inequality.

Corollary 1 (Resilience property). Let the conditions of Theorem 1 hold with the additional requirement that γ ψ = 1. Then

X -X ≤ h 2βΣ(ε) D 1 -2pr ∀ X ∈ E(Y ). (33)
Proof: The proof is immediate by considering the bound in [START_REF] Sharon | Minimum sum of distances estimator: Robustness and stability[END_REF] and observing that δ(ε) expressed in [START_REF] Sinopoli | Kalman filtering with intermittent observations[END_REF] vanishes when γ ψ = 1, hence eliminating completely the contribution of the extreme values of {ft} 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(βΣ(ε)) with g ∈ K∞ defined by g(α) = h 2α/ D(1 -2pr) . Moreover, βΣ(ε) reduces to t∈Tε ψt(ft) = infΩ∈F r d(F -Ω) when the process noise is zero (see [START_REF] Shoukry | Event-triggered state observers for sparse sensor noise/attacks[END_REF] and Lemma 4). Hence, E qualifies, in the sense of Definition 2, as an estimator which is resilient to the set Fr of measurement noise defined in [START_REF] Mishra | Secure state estimation against sensor attacks in the presence of noise[END_REF].

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 {ft} t∈T . Considering in particular the function βΣ(ε), we remark that it can be overestimated as follows

βΣ(ε) ≤ λ t∈T φt(wt) + |Tε|ε. ( 34
)
The first term on the left hand side of (34) represents the uncertainty brought by the dense noise {wt} over the whole state trajectory. It is bounded since {wt} 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 ft whose magnitude is smaller that ε.

Because βΣ is a function of ε, the bound in (33) represents indeed a family of bounds parameterized by ε. Since ε 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 ε , satisfy

ε ∈ arg min ε≥0 h 2βΣ(ε) D(1 -2pr) : r = |T c ε |, pr < 1/2 .
Another interesting point is that the inequality stated by Theorem 1 holds for any norm • on R n×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 pr is, the tighter the error bound will be, which suggests that the estimator is more resilient when pr is lower. A similar reasoning applies to the number D which is desired to be large here. These two parameters (i.e., pr 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 ψ is fundamental for achieving strict resilience. When ψ does not satisfy this property (i.e., when γ ψ = 1), the term δ(ε) in ( 23) is unlikely to vanish completely. However we can prevent it from growing excessively by an appropriate choice of ψ in [START_REF] Candès | An introduction to compressive sampling[END_REF]. To see this, assume for example that ψ is defined by ψ(y) = 1e -(y) . Then since ψ(y) ≤ 1 for all y, δ(ε) saturates to a constant value regardless of how large the ft are for t ∈ T c ε . On the other hand, this choice introduces a new technical challenge related to the fact that the function q in (32) is no longer a K∞ 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 ψ, we obtain the following approximate resilience result.

Corollary 2 (Case where γ ψ = 1). Let the conditions of Theorem 1 hold. Assume further that the loss function ψ in (7) is defined by ψ(y) = 1e -(y) where : R ny → R ≥0 satisfies (P1)-(P5). In particular, assume that property (P4) is satisfied by with a K∞ function q such that (9) is an equality relation. Also, let ε ≥ 0 be such that b(ε)

2βΣ(ε) + r(1 -γ ψ )r o (ε) D [1 -(1 + γ ψ )pr] < 1, (35)
where

r = r(ε) = |T c ε | and r o (ε) = max t∈T c ε ψt(ft) ≤ 1.
Then there exists a continuous and strictly increasing function hsat : [0, 1] → [0, 1] (obeying hsat(0) = 0 and hsat(1) = 1) such that for any norm

• on R n×T , X -X ≤ h -1 sat b(ε) ∀ X ∈ E(Y ). ( 36
)
with D in (35) defined as in the proof of Theorem 1 using the norm • .

Proof: That the particular function ψ 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 is just a saturated function in Ksat,1 (as defined in the notation section) with bounded range [0, 1]. This results in fact from Lemma 9 and the proof of Lemma 2. We can therefore write

q ( E ) ≤ 1 D (1 -(1 + γ ψ )pr) 2βΣ(ε) + (1 -γ ψ ) t∈T c ε ψt(ft) ≤ b(ε) < 1
with q ∈ Ksat,1. Note from the definition of the class Ksat,1, that q ( E ) < 1 implies that E < 1 (since otherwise we would have q ( E ) = 1). Letting hsat be the restriction of such a function q on [0, 1], we have q ( E ) = hsat( E ) ≤ b(ε) with hsat being invertible. We can now apply -1 sat to each member of this inequality to reach the desired result since b(ε) lies in the range of hsat.

D. Resilience to attacks on the individual sensors

We now consider the situation where the matrix S ∈ R ny ×T formed from {st} in (2) may be sparse entrywise i.e., a relatively important fraction of the entries of S are equal to zero2 . 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 ψt to have a separable structure. To be more specific, let ψt be such that for any e = [e1

• • • en y ] ∈ R ny ψt(e) = ny i=1 ψti(ei) (37)
where, consistently with [START_REF] Candès | An introduction to compressive sampling[END_REF], ψti(ei) = ψ • i (Vtiei) with Vti ∈ R>0 and ψ • i : R → R+, i = 1, . . . , ny, being some loss functions on R enjoying the properties (P1)-(P5). As in the statement of Corollary 1, we shall require that γ ψ with Vt being a diagonal matrix having the Vti, i = 1, . . . , ny, on its diagonal.

To state the resilience property in this particular setting, we partition the index set T × S of the entries of S as

Λε = {(t, i) ∈ T × S : ψti(fti) ≤ ε} Λ c ε = {(t, i) ∈ T × S : ψti(fti) > ε} (39)
with fti denoting the i-th entry of the vector ft ∈ R ny . Also, in order to account for the specificity of the new scenario, let us refine slightly the r-Resilience index (21) to be pr = sup

Z∈R n×T Z =0 sup I⊂T×S |I|=r (t,i)∈I ψti(c ti zt) HΣ(Z) ( 40
)
where HΣ is defined as in [START_REF] Pasqualetti | Attack detection and identification in cyber-physical systems[END_REF] from ψt in (38) and c ti is i-th row of the observation matrix Ct. The difference between pr in [START_REF] Ren | Secure state estimation with byzantine sensors: A probabilistic approach[END_REF] and pr in (40) resides in the index sets for counting possible error occurrences which are T and T × 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 Σ defined by (1) with output Y together with the state estimator (8) in which φ is assumed to obey (P1)-(P5) and ψ is defined as in (38). Let ε ≥ 0 and set r = |Λ c ε | with Λ c ε defined in (39).

If the system is observable on [0, T -1] and if pr < 1/2, then there exists a K∞ function h such that for all norm To some extent, Theorem 2 can be viewed as a special case of Theorem 1 in which the function ψ is taken to be the 1 norm and the data set is modified to be T×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 Fr in (17) of measurement noise in the sense of Definition 2 (see the proof of Corollary 1).

• on R n×T , X -X ≤ h 2 βΣ(ε) D(1 -2pr) ∀ X ∈ E(Y ) (41
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 Σ defined by (1) and let r = |Λ c 0 | (which means that we consider every nonzero occurrence of fit as an outlier by taking ε = 0 in (39)). If the conditions of Theorem 2 hold, pr < 1/2, and if wt = 0 in (1) for all t, then the estimator defined by (8) retrieves exactly the state trajectory of the system, i.e., E(Y ) = {X}.

Proof: This follows directly from the fact that βΣ(0) = 0 in the case where there is no dense noise wt and ε = 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 = |Λ c 0 | of nonzero in the sequence {fti} (t,i)∈T×S is small enough for pr 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 max r : pr < 1/2 .

(42)

V. A SPECIAL VARIANT OF THE ESTIMATOR E

In this section, we consider a constrained reformulation of the estimator E defined in [START_REF] Cardenas | Secure control: Towards survivable cyber-physical systems[END_REF]. 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 wt in (1) is identically equal to zero and the sequence {ft} is sparse in the sense that its dense component vt 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 ZΣ ⊂ R n×T defined by

ZΣ = Z = z0 A0z0 • • • AT -1 • • • A1A0z0 : z0 ∈ R n
of possible trajectories starting in any initial state z0 ∈ R n . Following this idea, we consider the estimator E • defined by

E • (Y ) = arg min Z∈Z Σ VΣ(Y, Z).
Then E • (Y ) can be rewritten more simply in the form

E • (Y ) = Z = z0 A0z0 • • • AT -1 • • • A1A0z0 : z0 ∈ arg min z∈R n V • Σ (Y, z) (43)
where

V • Σ (Y, z) = t∈T ψt(yt -Mtz) (44)
with

Mt = CtAt-1 • • • A1A0 (45)
for all t ≥ 1 and M0 = C0. Hence the estimation of the state trajectory reduces to estimating the initial state x0. This can be viewed as a robust regression problem, like the ones discussed in [START_REF] Han | Convex optimization based state estimation against sparse integrity attacks[END_REF], [START_REF] Bako | On a class of optimization-based robust estimators[END_REF]. Generalizing a result in [START_REF] Bako | On a class of optimization-based robust estimators[END_REF], we derive next a necessary and sufficient condition for exact recovery of the true state, which holds if and only if arg min z∈R n V • Σ (Y, z) = {x0} with x0 being the exact initial state of the system Σ. 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 I of T, let

Ψ • I (z) = t∈I ψt(Mtz). ( 46
)
Definition 4 (r-th concentration ratio). Let {ψt} be a family of loss functions defined by [START_REF] Candès | An introduction to compressive sampling[END_REF] in which ψ is assumed to satisfy (P1), (P3) and (P5) with constant γ ψ = 1. Let M = {Mt} t∈T be a sequence of matrices such that the function Ψ • T defined in (46) is positive definite. We call r-th concentration ratio of M , the number defined by

νr(M ) = sup z∈R n z =0 sup I⊂T |I|=r Ψ • I (z) Ψ • T (z) (47)
At a fixed r, νr(M ) quantifies a genericity property for the sequence M = {Mt} t∈T . In view of the particular structure of the collection M in (45), note that Ψ • T is positive definite whenever the system Σ is observable on T. Furthermore, ν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 [0, T -1], then for all Z ∈ ZΣ initiated from z in R n , we have VΣ(0, Z) = Ψ • T (z) ≥ q( z ) for some K∞ function q. It follows that

νr(M ) ≤ sup z∈R n z =0 sup I⊂T |I|=r Ψ • I (z) q( z ) (48)
Hence the more observable (i.e., the larger the gain function q), the smaller νr(M ). For all (Y, z0) ∈ R ny ×T × R n with Y expressed columnwise in the form Y = y0 • • • yT -1 , consider now the following notations:

I 0 (Y, z0) = {t ∈ T : yt -Mtz0 = 0} I c (Y, z0) = {t ∈ T : yt -Mtz0 = 0}.
Theorem 3 (Exact Recoverability Condition). Consider the cost function (44) where M = {Mt} is assumed to have been constructed as in (45) from the matrices of system [START_REF] Bako | Identification of switched linear systems via sparse optimization[END_REF]. Assume that the loss functions {ψt} involved in (44) are defined by [START_REF] Candès | An introduction to compressive sampling[END_REF] in which ψ is assumed to satisfy (P1), (P3) and (P5) with constant γ ψ = 1. Let r be a positive integer. If the system (1) is observable on [0, T -1], then the two following propositions are equivalent:

(i) For all Y in R ny ×T and all z0 in R n , |I c (Y, z0)| ≤ r ⇒ arg min z∈R n V • Σ (Y, z) = z0 (49)
(ii) The index νr(M ) satisfies

νr(M ) < 1/2 (50)
Proof: (i) ⇒ (ii): Assume that (i) holds. Consider an arbitrary subset I of T such that |I| ≤ r. Let z0 = 0 be a vector in

R n . Construct a sequence Y in R ny ×T such that yt = 0 if t ∈ I and yt = Mtz0 otherwise. Then I c (Y, z0) ⊂ I, so that |I c (Y, z0)| ≤ r. It then follows from (i) that arg min z∈R n V • Σ (Y, z) = z0 which means that V • Σ (Y, z0) < V • Σ (Y, z) for all z ∈ R n , z = z0. In particular, V • Σ (Y, z0) < V • Σ (Y, 0) which, by taking into account the definition of Y , gives Ψ • I (z0) < Ψ • I c (z0), where I c = T \ I. Since Ψ • T (z0) = Ψ • I (z0) + Ψ • I c (z0), we see that Ψ • I (z0) Ψ • T (z0) < 1/2
This reasoning works for every nonzero z0 and for every subset I of T. We can hence conclude that νr(M

) < 1/2. (ii) ⇒ (i): Assume that (ii) holds. Let (Y, z0) ∈ R ny ×T × R n be such that |I c (Y, z0)| ≤ r.
We then need to prove that arg min z∈R n V • Σ (Y, z) = z0 . Since the assertion (ii) is assumed true, it follows from (47) and (50) that

2Ψ • I c (z 0 ) < Ψ • T (z 0 ) ∀z 0 ∈ R n , z 0 = 0 (51)
where, for simplicity, we have posed I c = I c (Y, z0). In the derivation of (51), we have used the obvious fact that r1 ≤ r2 ⇒ νr 1 (M ) ≤ νr 2 (M ). If we pose I = I 0 (Y, z0) = T \ I c , then the inequality ( 51) is equivalent to

t∈I c ψt(Mtz 0 ) < t∈I ψt(Mtz 0 ) (52)
Now we observe that for all t in I = I 0 (Y, z0), yt = Mtz0, so that ψt(Mtz 0 ) = ψt yt -Mt(z0 + z 0 ) . On the other hand, for t ∈ I c = I c (Y, z0), if we apply the GTI [START_REF] Fawzi | Secure estimation and control for cyber-physical systems under adversarial attacks[END_REF] with γ ψ = 1, we obtain

ψt(Mtz 0 ) = ψt yt -Mtz0 -(yt -Mt(z0 + z 0 ) ≥ ψt(yt -Mtz0) -ψt(yt -Mt(z0 + z 0 ))
Combining all these remarks with (52) yields

t∈I c ψt(yt -Mtz0) -ψt(yt -Mt(z0 + z 0 ) < t∈I ψt(yt -Mt(z0 + z 0 )) Rearranging this gives V • Σ (Y, z0) < V • Σ (Y, z0 + z 0 ) for all z 0 ∈ R n with z 0 = z0. This is equivalent to arg min z∈R n V • Σ (Y, z) = z0 . Hence (ii) holds as claimed.
From the statement of Theorem 3, we infer that under the assumption that only the sparse noise {st} is active (i.e., there is no dense noise (wt, vt)) in the system equations ( 1), E • (Y ) = {X} whenever νr(M ) < 1/2, i.e, the estimator E • returns exactly the true state. For a given system, if one can evaluate numerically the index νr(M ), then it becomes possible to assess the number rmax max {r : νr(M ) < 1/2} of gross errors that can be corrected by the estimator E • when applied to that system. We will get back to the computational matter in Section VI.

A. Special case of 0-norm loss based estimator

Consider the special case where the loss functions {ψt} are defined, for all t ∈ T, by

∀e ∈ R ny , ψt(e) = 1 if e = 0 0 otherwise (53)
This corresponds to the block 0-norm. Note that such functions satisfy the assumptions (P1), (P3) and (P5) requested in the definition 4 of νr(M ) and in the statement of Theorem 3. Hence νr(M ) is well-defined in this case.

Corollary 4. Consider system (1) under the assumption that wt = 0 for all t. Assume observability of the system on [0, T -1]. Consider the estimator (43) in which the cost V • Σ is defined from the family of loss functions {ψt} expressed in (53). Then for all (Y, z0 Proof: Let us start by observing that with the particular loss functions invoked in the statement of the corollary, Ψ • T (z) denotes the number of t ∈ T for which ψt(Mtz) = 0. It follows from the definition of µ(M ) that Ψ • T (z) ≥ Tµ(M ) + 1 for all z = 0. The reason for this is that if ψt(Mtz) was to be equal to zero more than µ(M ) -1 times, then z would be necessarily equal to zero. As a result we get

) ∈ R ny ×T × R n , |I c (Y, z0)| < T -µ(M ) + 1 2 ⇒ E • (Y ) = {X} ,
νr(M ) ≤ r T -µ(M ) + 1 Hence by applying Theorem 3, |I c (Y, z0)|/(T -µ(M ) + 1) < 1/2
is a sufficient condition for exact recovery by the 0-norm based estimator.

Remark 2. Assume that ψt is defined to be the counting norm, i.e., ψt(e) = e 0 (55)

Then ψ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

νr(M ) = sup z∈R n z =0 sup I⊂T×S |I|=r (t,i)∈I Mtiz 0 Ψ • T (z) (56)
where Mti refers to the i-th row of Mt. Further, let

μ(M ) = min k : ∀I ⊂ T × S, |I| = k ⇒ rank( MI ) = n
with MI ∈ R |I|×n denoting the matrix obtained by stacking the row vectors {Mti : (t, i) ∈ I}. Then a result similar to Corollary 4 is obtainable: if the number the measurements corrupted by a nonzero error (among the nyT available) is strictly less than (nyTμ(M )+ 1)/2, then the estimator E • expressed in (43) (with ψt being the 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 < ny sensors may be compromised by attackers, then exact recovery is achieved if

k < ny 2 - μ(M ) -1 2T . ( 57
)
Taking into consideration the fact that μ(M ) -1 < T , it can then be seen that ( 57) is equivalent to k ≤ ny/2 -1 where the notation r , for r ∈ R, refers to the smallest integer larger or equal to r.

To sum up, when the ψt are defined as in (55), the estimator (43) is able to return the true state matrix even when ny/2 -1 sensors get faulty over the entire observation horizon. This is reminiscent of a result stated in [START_REF] Fawzi | Secure estimation and control for cyber-physical systems under adversarial attacks[END_REF] which therefore appears to be a consequence of Theorem 3.

B. Stability of the class of estimators E • with respect to dense noise

We have argued that the class of estimators E • in ( 43) is able to obtain exactly the true state matrix when there is no dense noises (wt, vt) in the system equations and only the sparse noise {st} 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

xt+1 = At xt, x0 = x0 yt = Ct xt + st + (vt + ṽt), (58)
produces the same output as system (1). Here, ṽt = Ct wt, with wt =

t-1 k=0 At-1 • • • A k+1 w k a definition which uses the convention that the product At-1 • • • A k+1 = I if k = t -1.
Then the idea is to apply the estimator E • to (58) by neglecting the dense component (vt + ṽt) of the output equation. To state the resilience result for E • , consider for a given ε ≥ 0, a partition ( Tε, Tc ε ) of T defined as in [START_REF] Rockafellar | Variational analysis[END_REF] with ft replaced by ft st + (vt + ṽt) = ft + ṽt. Theorem 4. Consider the estimator (43) for the system (1). Assume that the loss functions {ψt} involved in (44) are defined by [START_REF] Candès | An introduction to compressive sampling[END_REF] in which ψ is assumed to satisfy (P1)-(P5) with constant γ ψ = 1. Let ε ≥ 0 and set r = | Tc ε |. Denote with N a norm on R n×T defined by N (Z) = max t∈T zt with zt being the t-th column of Z and • being a norm on R n . If the system (1) is observable on [0, T -1] and νr(M ) < 1/2, then there exists a K∞ function α such that for all norm

• on R n×T , N ( X -X) ≤ RΣα -1 (ρ) + max t∈T wt ∀ X ∈ E • (Y ), (59)
where RΣ is some constant depending on the system Σ and

ρ = 2 D1 1 -2νr(M ) t∈ Tε ψt( ft) with ft = ft + ṽt, and D1 = min z =1 Ψ • T (z). Proof: Let x0 ∈ arg min z∈R n V • Σ (Y, z).
We first provide a bound on the error e0 = x0-x0 with x0 denoting the true initial state of system [START_REF] Bako | Identification of switched linear systems via sparse optimization[END_REF]. By exploiting the fact that

V • Σ (Y, x0) ≤ V • Σ (Y, x0
) and noting that yt = Mtx0 + ft, we reach the inequality

t∈T ψt( ft -Mte0) ≤ t∈T ψt( ft).
By then reasoning quite similarly as in the proof of Theorem 1, we get

Ψ • T (e0) -2Ψ • T c ε (e0) ≤ t∈ Tε ψt( ft)
which, by exploiting (47) and the assumption that νr(M ) < 1/2, leads to

Ψ • T (e0) ≤ 2 1 -2νr(M ) t∈ Tε ψt( ft)
Applying now Lemma 2, we conclude that for any norm • on R n , there exists a K∞ function α such that e0 ≤ α -1 (ρ). Now by observing that for any

X ∈ E • (Y ), X -X = e0 A0e0 • • • AT -1 • • • A0e0 -0 w0 • • • wT -1
the result follows by posing 3

RΣ = max t∈T At-1 • • • A0 ind with
• ind being the matrix norm induced by the vector norm • on R n .

The interest in Theorem 4 is that it provides a condition of resilience for the estimator E • 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 ν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 pr, pr or ν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 [START_REF] Candès | The restricted isometry property and its implications for compressed sensing[END_REF] is comparatively as hard to evaluate. We note however that when the dimension of the state is small enough, νr(M ) can be exactly computed by taking inspiration from a method presented in [START_REF] Sharon | Minimum sum of distances estimator: Robustness and stability[END_REF] 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 [START_REF] Bako | On a class of optimization-based robust estimators[END_REF]. The next lemma provides such an overestimate for νr(M ).

Lemma 5 (An estimate of νr). Assuming all quantities are welldefined (see the conditions in Definitions 3 and 4), the following statements hold:

(a) νr ≤ pr (b) If µ(M ) ≤ T -1 then νr(M ) ≤ rν o 1 + ν o , (60)
where

ν o = max t∈T min λ t ∈R T λt ∞ : VtMt = k∈T λ tk V k M k , λtt = 0
(61) 3 We use here the convention that

A t-1 • • • A 0 = I if t = 0.
In (61), the λ tk denote the entries of the vector λt ∈ R T and {Vt} refers to the sequence of nonsingular weighting matrices involved in [START_REF] Candès | An introduction to compressive sampling[END_REF].

The proof of statement (a) is straightforward by noticing that (47) follows from [START_REF] Ren | Secure state estimation with byzantine sensors: A probabilistic approach[END_REF] by constraining the variable Z to be in ZΣ. As to the proof of statement (b), it follows a similar reasoning as the proof of Theorem 2 in [START_REF] Bako | On a class of optimization-based robust estimators[END_REF]. The interest of this lemma is twofold. First it suggests that the resilience condition of E • is weaker (in the sense that it is easier to achieve) than that of E. Second, it provides an upper bound on νr(M ) which can be computed by solving a convex optimization problem (see Eqs (60)-( 61)). More specifically, given ν o in (61), checking numerically whether |T c ε | < 1/2(1 + 1/ν o ) provides a sufficient condition for ν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 pr defined in (40) can also be overestimated via convex optimization if the loss functions φt and ψt in ( 6)-( 7) are both taken to be norms. Lemma 6 (An estimate of pr). Consider the resilience parameter pr defined in (40) where we assume that ψti(e) = |e| for all (t, i, e) ∈ T × S × R and φt is an arbitrary norm. Then

pr ≤ r b1 (62)
where b1 = inf

(t,i)∈T×S inf Z∈R n×T HΣ(Z) : c ti zτ = 1 = 1 p1 . (63)
Proof: (See Appendix D) Note, under the assumptions of Lemma 6, that inf Z∈R n×T HΣ(Z) : c ti zτ = 1 is a convex optimization problem for any given (t, i). Hence, solving for b1 in (63) requires solving T ny convex problems and picking the smallest value among all. The interest of the lemma is that it provides an overestimate of pr which is numerically computable. Based on the so obtained overestimate of pr, we see from Theorem 1 that the estimator ( 8) is resilient to r outliers if r < b1/2. Moreover, we can deduce an underestimate of the number of outliers that the estimator ( 8) is able to handle as rmax = max {r : r < b1/2}. As a last remark in this paragraph, let us observe that Lemma 6 is also applicable to overestimate pr defined in [START_REF] Ren | Secure state estimation with byzantine sensors: A probabilistic approach[END_REF] in the case of a single-output system, i.e., when ny = 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 4 , an example of Linear Time-Invariant system in the form [START_REF] Bako | Identification of switched linear systems via sparse optimization[END_REF]. We select a single-input single-output example where the pair (A, C) is given by

A = 0.7 0.45 -0.5 1 , C = 1 2 . ( 64
)
We instantiate the loss functions in (5) as follows: For all t in T and for all (z, e) ∈ R n × R ny , φt(z) = ψ(Wtz) and ψt(e) = ψ(Vte) where the weighting matrices Wt and Vt and the functions φ and ψ will be specified below for each experiment.

A. Numerical certificate of exact recoverability Suppose in this section that the process noise wt and the dense component vt of ft (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 {st}. The times of occurrence of the nonzeros values in the sequence {st} 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 E • is implemented by directly solving the optimization problem defined in (43) through the CVX interface [START_REF] Grant | CVX: Matlab software for disciplined convex programming[END_REF]. Note that the implementation of the estimator (43)-( 44) requires computing the matrices Mt 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 Mt which might not be of the same order of magnitude. To preserve the contribution of each term of (44), we introduce special weighting matrices {Vt} (in the loss function ψ selected as the 1 norm) to normalize the rows of these matrices so that they all have unit 2-norm. Vt is therefore selected to be a diagonal matrix of the form Vt = diag(Vt1, • • • , Vtn y ), where

Vti = 1/ c i A t 2 if c i A t = 0 1 otherwise . (65)
Here, c i , i = 1, . . . , ny, denote the i-th row of the matrix C. Indeed the effect of the weighting function in (44) is equivalent to changing yt and Mt respectively to ỹt = Vtyt and Mt = VtCA t . Posing M = Mt , it can be checked using the methods discussed in Section VI (See Eq. ( 60)) that at least rmax = 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 |Λ c 0 |/T of nonzero values in the sequence {st}. 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%.

B. Performances in the presence of dense noise

We consider now the more realistic scenario where the process noise {wt} and the measurement noise {vt} 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 [-a, a]. For comparison purpose, we conduct the estimation with several estimators:

• an instance of the estimator [START_REF] Cardenas | Secure control: Towards survivable cyber-physical systems[END_REF], denoted E 2 2 , 1 in the sequel, in which the loss function φt is quadratic and ψt is the 1-norm and λ = 1000 (see ( 5) and ( 6)-( 7))

• an instance of the estimator (8) denoted E 1 , 1 in which both loss functions φt and ψt are the 1-norm with λ = 10

• the estimator E • defined in (43) In addition we implement oracle versions 5

of E 2 2 , 1 and of E 2 2 , 2 2
(the latter corresponding to an instance of E where both φt and ψt are instantiated as quadratic functions).

Experiment 1: Resilience test. Keeping the level of both dense noises (i.e., wt and vt) 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 E 2 2 , 1 , E 1 , 1 and E • 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 {st} ranging from 0 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 E 2 2 , 1 and those obtained by a standard least squares estimator E 2 2 , 2 2 (i.e. with φt and ψ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 {st}. 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 {st} 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 rmax = 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 pr proposed in Lemma 6.

Experiment 2: Stability with respect to dense noise. Now, we fix the sparsity level of the time sequence {st} to 0.2 and let the powers of the dense noise {(wt, vt)} 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 φt and ψ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 λ in E. To assess the influence of the regularization parameter λ on the performance of the estimators E 2 2 , 1 and E 1 , 1 , we fix the amplitudes of both dense noises {wt} and {vt} at the same level as in Experiment 1 (i.e. SNR equal to 30dB and ratio of non-zeros entries in {st} equal to 30%). In this setting, we consider a set of values of λ 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 when λ 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 E • in Experiment 1 in the same configuration, hence suggesting that E tends indeed to E • in behavior as λ becomes large. Finally, it is interesting to observe that both performance curves exhibit minima located around λ = 750 and λ = 10 for E 2 2 , 1 and E 1 , 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 10 -3 10 -2 10 -1 10 0 10 1 10 2 10 3 10 4 10 5 10 6 , 1 E 1 , 1 for the system (64), SNR= 30dB and 30% of non-zero entries in {st} for different values of the regularization parameter λ. may be convex (e.g., p-norms) or nonconvex (e.g., 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. q2(λ2). It follows that q(λ1) < min i∈{1,2} qi(λ2) = q(λ2) and hence q is strictly increasing. We now show that q(λ) tends to infinity when λ → +∞. Let M > 0 be an arbitrary positive number. Since q1 and q2 tend to infinity, there exist η1 and η2 such that λ ≥ η1 ⇒ q1(λ) ≥ M and λ ≥ η2 ⇒ q2(λ) ≥ M . By taking η = max i∈{1,2} ηi, it holds that q(λ) ≥ M whenever λ ≥ η, or equivalently that, lim λ→+∞ q(λ) = +∞.

Proof of Lemma 7: The sum ξ1 + ξ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 ξ1 • satisfies properties (P1)-(P3) too.

Proof of (j): Assume that ξ1 and ξ2 satisfy (P4) with K∞ functions q1 and q2 respectively. For all λ = 0 and all Z ∈ R a×b , (9) yields ξi(Z) ≥ min

If we define q so that for all λ ∈ R ≥0 , q(λ) = min i∈{1,2} qi(λ), then q is a K∞ function (see Lemma 8 above) such that for all λ = 0 and Z ∈ R a×b , ξ1(Z) + ξ2(Z) ≥ q 1 |λ| (ξ1(λZ) + ξ2(λZ)) given the linearity of . We can then conclude that ξ1 • also verifies property (P4).

Proof of (jj): Assume that ξ1 and ξ2 satisfy (P5) for γ1 and γ2 respectively. Let γ = min i∈{1,2} γi. Similarly to the first case, for all Z1, Z2 in R a×b and i in {1, 2}, (10) yields ξi(Z1 -Z2) ≥ γξi(Z1)ξi(Z2) (70) which gives ξ1(Z1-Z2)+ξ2(Z1-Z2) ≥ γ (ξ1(Z1) + ξ2(Z1))-(ξ1(Z2) + ξ2(Z2)) (71) therefore ξ1 + ξ2 satisfies property (P5). Moreover, for all Z1 and Z2 in R c×d , ξ1( (Z1 -Z2)) = ξ1( (Z1) -(Z2)) ≥ γξ1( (Z1)) -ξ1( (Z2))

(72) so ξ1 • satisfies (P5) too.

B. Proof of Lemma 3

(i) ⇒ (ii): Assuming that the system is observable on the interval [0, T -1], we need to prove that there exists a K∞ function q which verifies [START_REF] Huber | Robust Statistics[END_REF]. The idea of the proof is to apply Lemma 7 to the function F of R n×T defined by F (Z) = VΣ(0, Z) with VΣ defined as in [START_REF] Candès | Highly robust error correction by convex programming[END_REF]. To begin with, we note that F can be decomposed as F = ξ • where ξ : R n×(T -1) ×R n×T → R ≥0 is a loss function such that for Z = z0 To apply Lemma 7 to F , we need to check that F fulfills the properties (P1)-(P3). In virtue of the assumptions on φt and ψt agreed in the statement of the lemma, the first two properties are obviously satisfied. The third will be satisfied if is injective, a propriety which we now check. Let Z be such that (Z) = 0. Then Therefore, thanks to the recursive relation (73), we can conclude that Z = 0, and so, the linear mapping is injective. We can therefore apply Lemma 7 to conclude that F satisfy indeed (P1)-(P4). Now, consider a matrix norm • ind on R n×T induced by two vector norms • T and • defined respectively on R T and R n in the sense that

Z ind = sup η∈R T η =0 Zη η T
Applying Lemma 2 to F with the so-defined induced norm, we infer that there exists D > 0 defined as in [START_REF] Grant | CVX: Matlab software for disciplined convex programming[END_REF] and a K∞ function q , such that for all Z in R n×T ,

F (Z) ≥ Dq ( Z ind) (75)
If we denote with e1 the canonical vector of R T with all entries equal to zero except the first one which is equal to 1, then Ze1 = z0. However, by definition of the induced norm, we know that Ze1 / e1 T ≤ Z ind. Therefore, as q is an increasing function, we get that q ( z0 / e1 T ) ≤ q ( Z ind). By posing q : λ → Dq (λ/ e1 T ), it is easy to see that q is a K∞ function so that for all Z in R n×T , V (0, Z) = F (Z) ≥ q( z0 ).

(ii) ⇒ (i): Assume that there exists q in K∞ such that for all Z = z0 z1 . . . zT -1 in R n×T such that (14) holds. We want to prove that the matrix O0,T -1 defined in ( 13) is of full column rank, which is equivalent to showing that for z in R n , O0,T -1z = 0 implies z = 0. For all z ∈ R n , construct a sequence

Z * = z * 0 • • • z * T -1
as follows: z * 0 = z and z * t+1 = Atz * t for all t ∈ {0, . . . , 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

V (0, Z * ) = T -1 t=0 ψt(Ctz * t ) ≥ q( z * 0 ) (76)
Now, observe that if O0,T -1z = 0, then it follows from the recursive relation z * t+1 = Atz * t that for all t in {0, . . . , T -1}, Ctz * t = 0. Injecting this in (76) imposes that q( z * 0 ) ≤ 0 which necessarily implies that z = 0 as q is a K∞ function. Therefore, the matrix O0,T -1 is injective and the system is observable on the interval

  X), which gives explicitly λ t∈T φt(xt+1 -At xt) + t∈T ψt(yt -Ct xt) ≤ λ t∈T φt(wt) + t∈T ψt(ft) (27) Using the fact that xt+1 = Atxt + wt from (1) and applying the GTI and the symmetry properties of φt, we can write φt(xt+1 -At xt) = φt(xt+1 -xt+1 -At(xtxt) + wt) ≥ γ φ φt(et+1 -Atet)φt(wt) with et = xtxt. It follows that the first term on the left hand side of (27) is lower bounded as follows λ t∈T [γ φ φt(et+1 -Atet)φt(wt)] ≤ λ t∈T φt(xt+1 -At xt).

• i = 1 .

 1 It follows that one can set ψ • i to be the absolute value without loss of generality. Let therefore set ψ • i (ei) = |ei| so that ψti(ei) = |Vtiei| and ψt(e) = Vte 1 (38)

  ) with X denoting the true state matrix from (1) and βΣ(ε) defined by βΣ(ε) = λ t∈T φt(wt) + (t,i)∈Λε ψti(fti) D and h are defined as in the statement of Theorem 1 with HΣ being constructed from ψ in (38).

  where µ(M ) defined by µ(M ) = min k : ∀I ⊂ T, |I| = k ⇒ rank(MI ) = n (54) is the minimum number k such that any matrix MI ∈ R |I|ny ×n formed by stacking vertically the matrices of the collection {Mt : t ∈ I} indexed by I ⊂ T with |I| = k, has full column rank.

Figure 1 :

 1 Figure 1: Probability of exact recovery (expressed in percentage) by the estimator (43) in the presence of only sparse measurement noise {st}. The level of sparsity of the noise is expressed in terms of a fraction of nonzero values in the sequence {st : t ∈ T} with |T| = T = 100.

Figure 2 :

 2 Figure 2: Average relative estimation error (in logarithm scale) induced by different estimators versus sparsity level of the sparse noise {st}. The relative error is expressed here as X -X 2 / X 2 where X and X denote the true and estimated state matrices respectively. Parameters of the estimator E in (8): λ = 1000, Wt = I2 and Vt = 1 for all t.

Figure 3 :

 3 Figure 3: Average relative estimation error (in log scale) induced by different estimators for different levels of both dense noises wt and vt. Parameters of the estimator E in (8): λ = 1000, Wt = I2 and Vt = 1 for all t.

1 Estimator E 1 , 1 Figure 4 :

 1114 Figure 4: Average relative estimation error induced by E 2, 1 E 1 , 1 for the system (64), SNR= 30dB and 30% of non-zero entries in {st} for different values of the regularization parameter λ.

  ) therefore ξ1 + ξ2 verifies property (P4). Besides, for all λ = 0 andZ in R c×d , ξ1( (Z)) ≥ q1 1 |λ| ξ1(λ (Z)) = q1 1 |λ| ξ1( (λZ))(69)

2 ,

 2 • • • zT -2 in R n×(T -1) , Y = y0 • • • yT -1 in R ny ×T , and : R n×T → R n×(T -1) × R n×T a linear mapping such that for all Z = z0 • • • zT -1 in R n×T , (Z) = z1 -A0z0 • • • zT -1 -AT -2zT -C0z0 • • • CT -1zT -1 .

  ∀t ∈ {0, . . . , T -2}, zt+1 -Atzt = 0 (73) ∀t ∈ {0, . . . , T -1}, Ctzt = 0 (74) An immediate consequence of (73)-(74) is that O0,T -1z0 = 0 which yields z0 = 0 because the system is observable on [0, T -1].

Lemma 9 .

 9 [0, T -1]. C. Technical results for proving Corollary 2 This section contains some technical steps of the proof of Corollary 2. If : R ny → R ≥0 satisfies (P1)-(P3) and (P5), then so does the function ψ defined by ψ(y) = 1e -(y) . Moreover if

  , tends to suggest that low values of λ yield quite poor results. Conversely,

	Relative estimation error	10 -1 10 0 10 1	Estimator E 2 2 , 1 Estimator E 1 , 1 Estimator E • Oracle E 2 2 , 1 Oracle E 2 2 , 2 2
		0	0.2	0.4	0.6	0.8	1
				Sparsity level of sparse noise

Note that radial unboundedness is equivalent to level-boundedness in the terminology of[START_REF] Rockafellar | Variational analysis[END_REF]. theoretical guarantees of reaching global optimality with general purpose solvers.IV. THE RESILIENCE PROPERTY OF THE PROPOSED CLASS OFESTIMATORSIn this section, we prove that the state estimator proposed in[START_REF] Cardenas | Secure control: Towards survivable cyber-physical systems[END_REF] possesses the resilience property under some conditions. More specifically, our main result states that the estimation error 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 {ft} provided that the number of such outliers is below some threshold.

In 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.

While the theory presented in the previous sections applies to a general LTV system, the choice of an LTI example for the numerical illustrations is made here just for the sake of simplicity.

By oracle version of an estimator, we refer here to an implementation of this estimator which is aware of the sparse noise sequence {st}.

APPENDIX

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

A. A useful technical lemma

Lemma 7. Let ξ1, ξ2 : R a×b → R ≥0 be two functions which satisfy properties (P1)-(P3) and let : R c×d → R a×b be an injective linear mapping. Then ξ1 + ξ2 and ξ1 • verify (P1)-(P3). In addition, the following holds:

(j) If ξ1, ξ2 verify (P4), then ξ1 + ξ2 and ξ1 • verify (P4) .

(jj) If ξ1, ξ2 verify (P5), then ξ1 + ξ2 and ξ1 • 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 γ, their sum still verifies those properties.

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

Lemma 8 (Minimum function of two K∞ functions). If q1 and q2 are two K∞ functions, then so is the function q defined by

Proof: We have to prove that q is continuous, strictly increasing and satisfies q(0) = 0 and lim λ→+∞ q(λ) = +∞. First of all, it is clear that q(0) = 0. Also, continuity of q is immediate from that of q1 and q2 by noting that q = (q1 + q2 -|q1 -q2|)/2. To see the strict increasingness of q, consider λ1 and λ2 in R ≥0 such that λ1 < λ2. Then q(λ1) ≤ q1(λ1) < q1(λ2) and q(λ1) ≤ q2(λ1) < fulfills (P4), then ψ satisfies the same property but with a function q in Ksat,a for a = 1.

Proof: It is straightforward to check that ψ obeys (P1)-(P3). By assumption, obeys (P5). Denote therefore the associated constant with γ (which, by [START_REF] Fawzi | Secure estimation and control for cyber-physical systems under adversarial attacks[END_REF], is necessarily less than or equal to 1). To see then that (P5) is also satisfied by ψ, we just need to check that b) . From this it follows that for (77) to hold, it is enough that

which can indeed be checked to be true by applying the identity 1 + γ αγ ≤ α γ , see e.g., [3, Fact 1.9.2]. In effect, it follows from this identity that

In conclusion, (77) holds and therefore ψ satisfies (P5). It remains now to check (P4). This follows directly from Lemma 10 below, from which we know that ψ(y) ≥ q (1/λ)ψ(λy) with q is a saturated function in Ksat,1.

Lemma 10. Let : R ny → R ≥0 be a function satisfying properties (P1)-(P2) and (P4). In particular, assume that property (P4) is satisfied by with a K∞ function q such that (9) is an equality relation. Let

1e -(y/λ) for λ = 0 and y = 0. Then the function q : R ≥0 → [0, 1] defined by q (λ) = inf y =0 g(y, λ) for λ > 0 and q (0) = 0, is well-defined, continuous and strictly increasing on [0, 1]. Moreover we have

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

with q(1/λ) ≥ 1 for 0 < λ ≤ 1 and q(1/λ) < 1 for λ > 1. We therefore obtain

The so obtained q is clearly continuous wherever it is well defined. Moreover, since lim λ→0 q (λ) = q (0) = 0, we conclude that q is continuous on its entire domain. From the properties of q, we deduce that q is strictly increasing on [0, 1]. Lastly, we observe that the inequality in the statement of the lemma is a direct consequence of the definition of q .

D. Proof of Lemma 6

The starting point of the proof is the observation that for every integer r in {1, . . . , T }, pr ≤ r p1. Hence it suffices to show that p1 = 1/b1 and is as expressed in (63). Recall that by definition, p1 = sup

Without loss of generality, assume that c ti = 0 for all (t, i) ∈ T × S Then for any (t, i),

Recalling that HΣ(Z) is a norm under the conditions of the lemma, the second equality in the expression of βti above follows from the (strict) homogeneity property of norms. As to the last equality, it follows from the fact that c ti zt is a scalar which induces the possibility to replace the constraint |c ti zt| = 1 indifferently either by c ti zt = 1 or by c ti zt = -1. Now by invoking the definition of p1, it can be seen that p1 = sup (t,i)∈T×S