Igor Ryadchikov

Stanislav Aranovskiy
email: stanislav.aranovskiy@centralesupelec.fr

Evgeny Nikulchev

Jian Wang

Dmitry Sokolov

Differentiator-based velocity observer with sensor bias estimation: an inverted pendulum case study

Keywords: Velocity observer, differentiator, sensor bias, inverted pendulum, reaction wheel

In this paper we consider the problem of velocity estimation and stabilization for balancing an inverted pendulum equipped with a reaction wheel. A homogeneous differentiator is proposed for velocity estimation, and it is shown that a bias in sensor readings yields steadystate estimation error. The proposed observer is augmented with a reduced-order bias estimator and local asymptotic stability of the coupled observers is shown. The proposed solution is tested and compared with another approach on an experimental setup.

INTRODUCTION

This research is motivated by the walking robot with a non-anthropomorphic dynamic stabilization system that we have recently reported in [START_REF] Ryadchikov | Stabilization system of a bipedal non-anthropomorphic robot anywalker[END_REF]. For this robot, the idea was to provide auxiliary means of vertical stabilization by installing inside the robot's body two flywheels (shown in red in the Figure 1, left), improving vertical stability without affecting the rest of the mechanical system. The flywheels are orthogonal to each other and the problem of vertical stabilization of the robot can be considered for each axis separately. Thus, we consider a stabilization of a one-dimensional reactionwheel based inverted pendulum as a basic approximation for the robot motion control system.

Inverted pendulums, including reaction-wheel setups, have been extensively studied. [START_REF] Spong | Nonlinear control of the reaction wheel pendulum[END_REF] proposed a switching controller to swing the pendulum upwards from its downright equilibrium. The monograph [START_REF] Block | The reaction wheel pendulum[END_REF] deeply covers the pendulum motion, including linear and nonlinear models, and models of the sensors and actuators that are used for feedback control. Finally, for our paper, the most relevant work is the Cubli by [START_REF] Gajamohan | The cubli: A cube that can jump up and balance[END_REF], a very interesting project created at ETH Zürich.

The key element for the stabilization of an inverted pendulum is velocity estimation both for the pendulum and for the reaction wheel. While in the Cubli device, a gyroscope sensor is used to measure the velocity, in our case we use optical angular encoders, and the velocities are to be estimated from the available measurements. Velocity (or momenta) observer design for a mechanical system is a well-known and long-standing problem, see, for example, the recent work by [START_REF] Aranovskiy | A globally exponentially stable speed observer for a class of mechanical systems: experimental and simulation comparison with high-gain and sliding mode designs[END_REF] where several solutions are considered and compared. While nonlinear Luenberger-like observers and Kalman-based filters operate with the whole vector of measurements and can be computationally demanding, a common engineering approach is to consider each degree of freedom separately and estimate velocity as a derivative of a scalar position (angle) signal; in such a case velocity estimation is considered as a differentiator design problem. There are multiple known solutions for numerical differentiation, such as linear approximation based on Taylor series expansion, e.g. a firstorder difference used by [START_REF] Spong | Nonlinear control of the reaction wheel pendulum[END_REF], sliding-mode exact differentiators first proposed by [START_REF] Levant | Robust exact differentiation via sliding mode technique[END_REF], or high-gain differentiators as described by [START_REF] Vasiljevic | Error bounds in differentiation of noisy signals by high-gain observers[END_REF]. For this paper, we are particularly interested in the homogeneous differentiator (HOMD) proposed by [START_REF] Perruquetti | Finite-time observers: Application to secure communication[END_REF] that ensures finite-time con-vergence and is robust with respect to possible model uncertainties.

An important issue for pendulum stabilization is a sensor bias potentially present in the position measures. This bias can be caused by imprecisions in sensor placement; however, there are other possible sources. For example, in the one-dimensional prototype considered in this paper the motor cable creates a disturbance that shifts the equilibrium point from zero readings of the position sensor. As it is shown by [START_REF] Block | The reaction wheel pendulum[END_REF] and illustrated in §4, such a bias may shift the stable equilibrium of the closedloop system away from the origin, that is undesirable behavior for robot stabilization. Thus, the bias should be estimated and taken into account by the control design. Unfortunately, differentiator-based observers cannot be used here since the bias is not a derivative of any measured signal. One possible solution is to use an extended fullstate observer instead of differentiators. Another solution was used in the Cubli project by [START_REF] Gajamohan | The cubli: A cube that can jump up and balance[END_REF], where the sensor offset was estimated at a slower time scale from the observed equilibrium, see §6 for more details.

In this paper, we propose a novel solution. More precisely, a contribution of this paper is twofold. First, we design a finite-time HOMD-based velocity observer for a reaction-wheel inverted pendulum and use this observer for pendulum stabilization. Second, we propose to combine this observer with a reduced-order sensor bias estimator. Our analysis shows that the resulting coupled nonlinear observer does not preserve the finite-time convergence property, however, under a proper choice of tuning gains the coupled system is locally asymptotically stable and estimation errors converge to zero with a certain domain of attraction. Moreover, it is also shown that the gains can be tuned by the means of linear matrix inequalities (LMIs). The proposed solution is implemented in an experimental setup and compared with the solution proposed by [START_REF] Gajamohan | The cubli: A cube that can jump up and balance[END_REF].

The rest of the paper is organized as follows. In §2 we present a model of the considered system and a stabilizing LQR control. Next, in §3 we design HOMD-based observers to estimate pendulum and reaction wheel velocities. The behavior of the closed-loop system with biased measurements is analyzed in §4 and a reduced-order bias observer is proposed and analyzed in §5. Finally, experimental results and comparison with the bias observer used for the Cubli project are given in §6.

PENDULUM MODEL AND LQR CONTROL

Model description and used notations. In this section we follow the notations from [START_REF] Block | The reaction wheel pendulum[END_REF]: our main variables are θ and θ r , where θ is the angle between the pendulum and the vertical, and θ r is the angle of the reaction wheel w.r.t the pendulum. Refer to Figure 1 for an illustration; Table 1 provides all necessary measurements and notations used in the paper.

Neglecting the friction, we derive the equations of motion (the derivations are based on standard arguments for mechanical systems and are omitted here due to the lack of space): Resolution of the pendulum angle measurement, rad -6.28•10 -4

J r θr + J r θ = kI,

(J + J r ) θ + J r θr = mlg sin θ, (1)
where the symbols are as defined in Table 1, ml := m p l p + m r l r , J := J p + m p l 2 p + m r l 2 r , and I is the current in the motor windings. Assuming an internal fast-time-scale current loop, we consider I as our input control signal.

The control goal is to locally stabilize the pendulum in the upper equilibrium, that is to drive the variables θ, θ, and θr to zero, while θ r ∈ R. Let us further denote the upper equilibrium as the set Ω 0 := θ = θ = θr = 0, θ r ∈ R ⊂ R 4 , which is an invariant set of the system (1). In what follows we say that a control law (locally) stabilizes the system (1) if under this control law the set Ω 0 is (locally) attractive.

Define the state variable vector x := θ θ θr . Then the control goal is to construct the control input I such that the state x locally converges to the origin.

LQR control design. To achieve the goal, the model ( 1) is linearized around the equilibrium x eq := [0 0 0] yielding the linearized model ẋ = Ax + BI, θ = Cx, (2) where θ is the measured output and

A :=   0 1 0 mlg J 0 0 -mlg J 0 0   , B :=   0 -k J (J+Jr)k J Jr   , C := 1 0 0 .
To stabilize the pendulum the LQR is designed:

u = -Kx, (3) where
K := [k 1 k 2 k 3 ] is the gain vector minimizing the cost function ∞ 0
x (τ )Qx(τ ) + Ru 2 (τ ) dτ, where the matrix Q > 0 and the scalar R > 0 are the design parameters.

Under the control law (3), the upper equilibrium of the pendulum is locally asymptotically stable with a domain of attraction depending on the design parameters. However, to implement the law (3) measurements of the velocities θ and θr are required; if these variables are not measured directly, then a velocity observer should be used to generate the estimates θ and θr , and the law (3) takes the form

u = -K θ θ θr . (4)
The estimates θ, θr can be obtained with a linear Luenberger-like observer. Being widely used and rather simple in implementation, this approach is based on the linearized dynamics (2). Hence, the linear observer converges only in a neighborhood of the equilibrium and, moreover, the transient performance can significantly degrades for initial conditions being far from the origin. This drawback can be partially alleviated using linear observers with time-varying gains, e.g. an observer with gains scheduling, but such systems are typically harder to implement. In what follows, we propose to use nonlinear velocity observers, namely homogeneous differentiators, to obtain the estimates with a finite time of convergence.

FINITE TIME VELOCITY OBSERVER

By design, the homogeneous differentiator proposed in [START_REF] Perruquetti | Homogeneous finite time observer for nonlinear systems with linearizable error dynamics[END_REF]; [START_REF] Perruquetti | Finite-time observers: Application to secure communication[END_REF] is applied to a chain of integrators with known inputoutput injections and with a scalar output signal. Thus, to use it for velocity estimation for the pendulum system considered in this paper, we apply it to each degree of freedom separately, i.e. we construct two observers, one for the pendulum itself, and another for the reaction wheel.

Pendulum velocity observer.

From (1) we obtain the dynamics of the pendulum: θ = -k J I + mlg J sin(θ). Let xp ∈ R 2 be an estimate of the vector x p := θ θ and define the estimation error as e p := xp -x p , where we recall that e p,1 = θ -θ is measured. Denote for any real numbers x and α > 0

x α := |x| α sgn(x), where sgn(•) is the sign function. Following [START_REF] Perruquetti | Homogeneous finite time observer for nonlinear systems with linearizable error dynamics[END_REF], we construct the homogeneous velocity observer as a differentiator of x p,1 :

ẋp

,1 = xp,2 -k p,1 e p,1 αp , ẋp,2 = - k J I + mlg J sin(θ) -k p,2 e p,1 2αp-1 , θ = xp,2 , (5)
where α p , k p,1 , k p,2 are the design parameters. Then the error dynamics is given by ėp

,1 = e p,2 -k p,1 e p,1 αp , ėp,2 = -k p,2 e p,1 2αp-1 . ( 6
)
Note that e p = 0 is the only equilibrium of the system (6). Proposition 1. [START_REF] Perruquetti | Finite-time observers: Application to secure communication[END_REF]). Consider the system (6), where α p ∈ 1 2 , 1 and k p,1 , k p,2 are chosen such that the polynomial s 2 + k p,1 s + k p,2 is Hurwitz. Let e p (t) be the solution with the initial conditions e p (0) ∈ R 2 \ {0}. Then the origin is finite-time stable, i.e. there exists T = T (e p (0)) > 0 such that e p (t) is defined and unique on [0, T ), bounded, and lim t→T e p (t) = 0. T is called the settling-time function of the system (6). The settling-time function can be extended at the origin by T (0) = 0.

The proof of Proposition 1 follows applying Theorem 10 from [START_REF] Perruquetti | Finite-time observers: Application to secure communication[END_REF] to the system (6).

Reaction wheel velocity observer. For the reaction wheel dynamics we have from (1): θr = (J+Jr)k

JJr I -mlg J sin(θ). Let xr ∈ R 2 be an estimate of the vector x r := θ r θr and define the estimation error as e r := xr -x r , where e r,1 is measured. Similarly to (5), we construct the homogeneous velocity observer as ẋr,1 = xr,2 -k r,1 e r,1 αr ,

ẋr,2 = (J + J r )k JJ r I - mlg J sin(θ) -k r,2 e r,1 2αr-1 , θr = xr,2 , (7)
where α r , k r,1 , k r,2 are the design parameters. Then the error dynamics is given by ėr,1 = e r,2 -k r,1 e r,1 αr , ėr,2 = -k r,2 e r,1 2αr-1 .

(8)

The system ( 8) is similar to (6) and the finite-time stability of the equilibrium e r = 0 follows from Proposition 1 choosing 1 2 < α r < 1 and k r,1 , k r,2 such that the polynomial s 2 + k r,1 s + k r,2 is Hurwitz. Closed-loop system behavior. Given the finite-time observers, it follows that the control law (4), ( 5), ( 7) locally stabilizes the system (1). Indeed, it is known that the control law (3) that stabilizes the linearized system (2) is locally stabilizing for the nonlinear system (1). Let us denote the region of attraction under the control law (3) as Ω α . For the observers ( 5) and ( 7) there exists a common settling time T = T (e p (0), e r (0)), which is the maximum of the settling times of these observers, and the control laws (3) and ( 4) are equivalent for t ≥ T . Since the trajectories xp and xr are bounded, there exist constants ēp > 0, ēr > 0, and a set Ω β ⊂ Ω α , such that for all initial conditions satisfying x(0) ∈ Ω β , e p (0) < ēp , e r (0) < ēr trajectories x(t) stay in Ω α for t ∈ [0, T ].

Then starting from t = T we can consider the system (1) under the control law (4), ( 5), (7) as the system (1) under the control law (3) with the initial condition x(T ) ∈ Ω α ensuring x(t) → 0.

The set Ω β , which is the domain of attraction under the control (4), ( 5), ( 7), depends on the design parameters and initial conditions, and is typically smaller than the set Ω α , which is the domain of attraction under the control (3). However, given e p (0) = e r (0) = 0 we have T = 0 and the sets Ω β and Ω α coincide.

BEHAVIOR UNDER BIASED MEASUREMENT

The proposed control law (4), ( 5), ( 7) was successfully implemented with the experimental testbed, see Section 6 for details on the equipment. Given that the optical encoder measuring θ is perfectly adjusted, stabilization of the pendulum in the upper equilibrium and convergence x(t) → 0 were achieved.

However, it was found that if the zero reading of the sensor does not coincide with the equilibrium position and the angle θ is measured with a certain constant offset, then the control goal is not achieved. This problem is well-known, see e.g. [START_REF] Block | The reaction wheel pendulum[END_REF]; [START_REF] Gajamohan | The cubli: A cube that can jump up and balance[END_REF]. Denote the constant offset as d and let the measured signal be y = θ + d.

Using y in place of the real value θ in (3) yields u = -K y θ θr = -Kx-K [d 0 0] . Then the equilibrium of the linearized system (2) under this control law can be found as

x eq,d = (A -BK) -1 BK [d 0 0] = 0 0 -d k1 k3
.

It means that when measurements are biased, then the pendulum is stabilized at the upper position, however the reaction wheel does not stop and θr converges to a nonzero constant value.

Consider now behavior of the velocity observers ( 5), ( 7) under biased measurements, where we impose the following assumption.

Assumption 1. The sensor offset d is sufficiently small and sin(d) ≈ d, cos(d) ≈ 1.

Using y as a measurement of θ in ( 5) and noting that θ -y = e p,1 -d, the error dynamics (6) takes the form

ėp,1 = e p,2 -k p,1 e p,1 -d αp , ėp,2 = a 1 cos(θ)d -k p,2 e p,1 -d 2αp-1 , (10)
where we define a 1 := mlg J . Obviously, e p = 0 is not an equilibrium of this system. Considering the pendulum in the upper half-plane where cos(θ) > 0, at the equilibrium of the system (10) we have the steady-state velocity estimation error

e 0 p,2 := k p,1 a 1 cos(θ) k p,2 |d| αp 2αp -1 sgn(d). (11)
Analyzing the error dynamics of the reaction wheel observer, it can be shown that the steady-state velocity estimation error is e 0 r,2 = -e 0 p,2 . Applying these results to the control law (4), it follows

x eq,d = (A -BK) -1 BK   d e 0 p,2 e 0 r,2   =   0 0 e 0 p,2 k3-k2 k3 -d k1 k3   . ( 12
)
The conclusion is that when the measurements are biased, the velocity estimates do not converge to the real values but are biased. The pendulum stabilizes in the upper equilibrium, but the reaction wheel keeps rotating. To deal with the biased sensor we propose to construct a reducedorder bias observer, described in the following section. Remark 1. It is worth noting that if the reaction wheel position θ r is measured with a constant bias d r while the pendulum position θ is measured without a distortion, then the error dynamics of the observer (7) becomes ėr,1 = e r,2 -k r,1 e r,1 -d r αr , ėr,2 = -k r,2 e r,1 -d r 2αr-1 .

It follows that at the equilibrium we have e r,2 = 0, i.e. a bias of the reaction wheel sensor does not lead to a steady-state bias of the estimate xr,2 of the reaction wheel velocity.

REDUCED-ORDER BIAS OBSERVER

Reduced-order observer design

Define the extended state vector z := x d = θ θ d . Then the measured signal is given by y = [1 0 1] z. Let ẑ be an estimate of z and define e := ẑz. Then e 3 is the bias estimation error. Motivated by Assumption 1, we expect that the estimate of the bias is also small and the following approximation holds: sin(z 1 ) = sin(y -z 3 ) ≈ sin(y -ẑ3 ) + cos(y -ẑ3 )e 3 . Then the dynamics of z can be written as

ż =   z 2 -k J I + a 1 sin(z 1 ) 0   = z 2 a 1 cos(y -ẑ3 )e 3 0 + 0 β z 0 , ( 13
)
where the signal β z is available, β z := -k J I +a 1 sin(y -ẑ3 ). If the velocity z 2 is measured, then the reduced-order Luenberger-like linear observer of d can be constructed as

v = β z , ẑ3 = L(v -z 2 ), d = ẑ3 , ( 14
)
where L is the design parameter. Then ė3 = -La 1 cos(y -ẑ3 )e 3 .

(15) As we consider the pendulum around the upper equilibrium, it is reasonable to impose the following assumption. Assumption 2. For trajectories of the system there exists 1 c 0 > 0 such that min cos(y -d), cos(θ) > c 0 along these trajectories.

Under this Assumption, it is obvious that choosing L > 0 in (15) yields to exponential convergence of e 3 to zero. With the estimator ( 14), the stabilizing control law (4) can be rewritten as

u = -K y - d θ θr . (16)

Closed-loop convergence

Let us now consider what happens when the reduced-order bias observer is in the loop with the homogeneous velocity observer. In such a case the state z 2 in the observer ( 14) is not measured directly but generated by the observer (5). At the same time, the estimate ẑ3 is used to compensate the bias and we substitute y -ẑ3 in place of θ in (5). The joint observers dynamics is now given by ż1 =ẑ

2 -k p,1 ẑ1 + ẑ3 -y αp , ż2 = - k J I + a 1 sin(y -ẑ3 ) -k p,2 ẑ1 + ẑ3 -y 2αp-1 , v = - k J I + a 1 sin(y -ẑ3 ), ẑ3 =L(v -ẑ2 ). (17)
To proceed with the dynamics analysis note that ẑ1 + ẑ3y = e 1 + e 3 and that under the imposed assumptions the following approximation holds: sin(y-ẑ 3 ) = sin(z 1 -e 3 ) ≈ sin(z 1 )-cos(z 1 )e 3 . Denote s := [e 1 + e 3 e 2 e 3 ] . Then the error dynamics of the observer yields ṡ1

= s 2 -k p,1 s 1 αp + Lk p,2 s 1 2αp-1 , ṡ2 = -a 1 cos(z 1 )s 3 -k p,2 s 1 2αp-1 , ṡ3 = Lk p,2 s 1 2αp-1 . (18)
Since cos(z 1 ) > c 0 , the only equilibrium of ( 18) is the origin s = e = 0. However, the error dynamics ( 18) is not homogeneous for α p < 1 and thus it does not have the finite-time stability property. An important observation is that for α p = 1 the system (18) becomes a linear timevarying (due to cos(z 1 )) system. The following proposition establishes convergence of the extended observer for α p in a vicinity of 1. Proposition 2. Consider the observer (17) for the system (13) under Assumption 2 with the error dynamics (18). Define

A m := -k p,1 + Lk p,2 1 0 -k p,2 0 -a 1 c 0 Lk p,2 0 0 , A M := -k p,1 + Lk p,2 1 0 -k p,2 0 -a 1 Lk p,2 0 0
, and

A 1 :=   Lk p,2 -1 2 k p,1 0 0 -k p,2 0 0 Lk p,2 0 0   ,
where c 0 is defined in Assumption 2. If the parameters of the observer ( 17) are chosen such that the there exists a positive-definite symetric matrix P satisfying the LMIs P A m + A m P + γP ≤ 0,

P A M + A M P + γP ≤ 0, P A 1 + A 1 P ≤ 0, (19)
for some γ > 0, then there exist ε α > 0 and a compact set Ω α , such that for α ∈ (1 -ε α , 1] and all initial conditions s(0) ∈ Ω α it holds s → 0 and lim t→∞ |ẑ -z| = 0.

The proof of the proposition is omitted due to the lack of space.

Remark 2. The LMIs ( 19) are used to check if the given set of gains k p,1 , k p,2 , and L is stabilizing. Using standard LMIs arguments, see [START_REF] Boyd | Linear matrix inequalities in system and control theory[END_REF], this feasibility check can be reformulated in order to find (if any) such gains, that (19) are feasible.

Estimation of the domain of attraction Ω. Suppose we are given P being a solution of (19) for some c 0 and γ and ∆ M being the maximum value ∆ α . First, compute Q and find λ Q as the smallest generalized eigenvalue of γP and Q, i.e. the smallest value such that det (γP

-λ Q Q) = 0. Then V is negative definite if ln s 2 1 ∆ α < λ Q . Next, compute s M as s M = exp λ Q ∆ M .
To estimate Ω we have to find such C that s P s < C ⇒ s 2 1 < s M , i.e. we are looking for the ellipsoid s P s = C that touches the plane s 1 = √ s M . Since P is positive definite, there exists matrix R = R such that P = RR. Then s P s = v v where the new coordinates are v = Rs. Now s P s = C

n v = R -1 n s .
Then the distance from the plane in v coordinates to the origin is given by

√ C = n v pv |nv| = s M n s R -2 ns , hence C = s M (P -1 ) 1,1
, where P -1 1,1 is the (1, 1) element of the matrix P -1 . Remark 3. It is worth noting that the set Ω is obviously a conservative estimate of the real domain of attraction.

EXPERIMENTS

The hardware for the tests (shown in Figure 1) is assembled from off-the-shelf components. The reaction wheel (a bicycle brake rotor) is driven by a 70W Maxon EC 45 flat brushless DC motor, which is controlled using Maxon EPOS2 50/5 in torque mode. A STM32F407 discovery board was chosen as the main computing unit.

For our experiments, we consider the problem of the upright pendulum stabilization, where the initial position is in a neighborhood of the equilibrium. We use the homogeneous velocity observers ( 5) and ( 7) with the feedback control law (4). At it is mentioned in Section 4, if the optical encoder is perfectly adjusted and d = 0 in (9), then the proposed controller achieves the goal and stabilizes the pendulum. In what follows we consider the problem of stabilization under biased measurements, where the sensor offset d in ( 9) is d ≈ 5 • ≈ 0.09. Table 2 lists the parameters used in the experiments.

First, we illustrate that the same control law (4), (5), and (7) cannot stabilize the pendulum for d = 0. As it is shown in ( 12), the pendulum arrives to the upper position, but the reaction wheel does not stop and maintains a nonzero constant velocity. The experimental results are depicted in Figure 2. Here the measured angle y converges to a nonzero value, such that y ≈ d, and the physical angle θ ≈ 0. Velocity estimates xp,2 and xr,2 converge to nonzero values, where the estimation errors are predicted by (11). The physical velocities converge to θ ≈ 0 and θr ≈ 156 radians per second, that corresponds to (12).

Next we apply the approach by [START_REF] Gajamohan | The cubli: A cube that can jump up and balance[END_REF].

The key observation of this approach is that under the control law (4) the pendulum angle θ reaches zero, and is steady state it holds y ≈ d. Thus, the measurement y can be used as an estimate d, which is then injected in (16). To avoid stability issues, estimation of d is performed in a slower time scale than stabilization that is ensured by a low-pass linear filter: where τ d > 0. The drawback of this approach is that if τ d is large enough, then the stabilization goal is achieved but with a long transient time. If the value τ d is reduced to accelerate the transients, then the closed-loop stability can be compromised. Results of the experiments with the velocity observers ( 5), ( 7), bias observer (20), and control law ( 16) are depicted in Figure 3 for τ d = 20 and τ d = 2 and illustrate the trade-off between the transient time and the closed-loop stability.

τ d d dt d(t) + d(t) = y(t), (20)
Finally, we perform an experiment with the reduced-order bias observer (14). Results of the experiment are shown in Figure 4 and illustrate stabilization of the pendulum with fast transients.

CONCLUSION

The problem of velocity observer design for an inverted reaction-wheel pendulum under biased measurements was considered in the paper. The proposed solution combines a HOMD-based velocity observer with a reduced-order bias estimator. It is shown that the stability of the proposed observer is ensured under a proper choice of tuning pa- rameters. Experimental results illustrate the performance of the proposed approach.

As a further research direction, we consider robustness analysis for time-varying parameters and model uncertainties and application of the proposed nonlinear observer for walking robot stabilization, see [START_REF] Ryadchikov | Stabilization system of a bipedal non-anthropomorphic robot anywalker[END_REF]

  This work was supported by the Russian Ministry of Education and Science, through the grantNo. 8.2321No. 8. .2017. .

Fig. 1 .

 1 Fig. 1. Left: AnyWalker robot uses reaction wheels (shown in red) as an auxiliary stabilization system. Middle and right: 1D pendulum hardware and corresponding notations.

  is a sphere v v = C with the radius √ C in the new coordinates. The plane s 1 = √ s M in the s coordinates is defined by the normal vector n s := [1 0 0] and the point p s = [ √ s M 0 0] . In the v coordinates we have the point p v = Rp s and the normal vector

Fig. 2 .Fig. 3 .

 23 Fig. 2. The control law (4), (5), and (7) for d ≈ 0.09.

Fig. 4 .

 4 Fig. 4. The control (5), (7), and (16) with the reduced order observer (14) for d ≈ 0.09.

Table 1 .

 1 Parameters of the experimental setup.

	Description	Symbol	Value
	Mass of the pendulum, kg	mp	0.58
	Pivot -pendulum's center of mass dis-	lp	0.10
	tance, m
	Moment of inertia of the pendulum	Jp	3.8•10 -3
	around its center of mass, kg•m 2
	Mass of the reaction wheel, kg	mr	0.35
	Pivot -reaction wheel axis distance, m	lr	0.22
	Moment of inertia of the reaction	Jr	12.48 • 10 -4
	wheel, kg•m 2
	Current-to-torque gain, N/A	k	3.69 • 10 -2
	Sampling frequency, Hz	-	500
	Resolution of the reaction wheel angle	-	6.54•10 -2
	measurement, rad

Table 2 .

 2 Controller parameters used in experiments.

	Description	Symbol	Value
	LQR gains in (4) and (16)	K	-581.6 -83.4 -1.2
	Parameters in (5)	[k p,1 , k p,2 , αp]	[7, 12, 0.75]
	Parameters in (7)	[k r,1 , k r,2 , αr]	[7, 12, 0.75]
	Gain of the bias observer (14)	L	0.033

In practical applications this assumption can be ensured restricting, e.g. by a projection operator, possible variations of ẑ3 subject to known bounds of the displacement d.


Alexandre Kircher

Laurent Bako

Éric Blanco

Mohamed Benallouch

On the resilience of a class of Correntropy-based state estimators

Keywords: Secure state estimation, Maximum Correntropy Criterion (MCC), optimal estimation, Cyber-physical systems

This paper deals with the analysis of a class of offline state estimators for LTI discrete-time systems in the presence of an arbitrary measurement noise which can potentially take any value. The considered class of estimators is defined as the solution of an optimization problem involving a performance function which can be interpreted as a generalization of cost functions used in the Maximum Correntropy Criterion. The conclusion of the analysis is that if the system is observable enough, then the considered class of estimators is resilient, which means that the obtained estimation error is independent from the highest values of the measurement noise. In the case of systems with a bounded process noise, the considered class of estimators provides a bounded estimation error under the appropriate conditions despite not being designed for this scenario.

INTRODUCTION

The problem considered in this paper is the estimation of the (hidden) state of a system despite the presence of an arbitrary measurement noise which can take any value. This kind of problem typically occurs when designing an estimator to cope with faulty sensors [START_REF] Niedfeldt | Robust estimation with faulty measurements using recursive-RANSAC[END_REF]. Arbitrary noises are also a suitable way of modeling attacks of diverse natures, like replay attacks [START_REF] Mo | Secure control against replay attacks[END_REF], which pose some challenges in estimation scenarios where the data are collected via a communication network. This is case for example for Cyber-Physical Systems [START_REF] Cardenas | Secure Control: Towards Survivable Cyber-Physical Systems[END_REF].

To address the estimation problem in the face of such challenging uncertainties, many approaches have been developed in the literature among which faulty measurements isolation [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], compressive sampling-inspired methods [START_REF] Fawzi | Secure Estimation and Control for Cyber-Physical Systems Under Adversarial Attacks[END_REF][START_REF] Pajic | Attack-Resilient State Estimation for Noisy Dynamical Systems[END_REF], or event-triggered resilient estimation [START_REF] Shoukry | Event-Triggered State Observers for Sparse Sensor Noise/Attacks[END_REF][START_REF] Liu | Eventtriggered resilient filtering with missing measurements[END_REF].

In parallel to those methods, estimators based on the Maximum Correntropy Criterion (MCC) were developed to tackle the resilient state estimation problem. In order to design new filters more robust to non-Gaussian noises, the core idea is to use the correntropy as a similarity measure to design the cost function. However, there are several types of similarity measures and also several methods to derive the estimator. For instance, (Chen et al., 2017a;[START_REF] Liu | State space maximum correntropy filter[END_REF] use Gaussian kernel functions as cost functions and derives a fixed-point algorithm in order to update the posterior estimation in a Kalman Filter. The scope of estimators inspired by the MCC then grew, including for example optimization-based estimators with cost functions under the form of a sum of exponential absolute value functions [START_REF] Chen | Robust State Estimator Based on Maximum Exponential Absolute Value[END_REF], or online estimators which derive an MCC-based cost function under the form of a Kalman Filter with approximated weights based on the prior estimation at each step [START_REF] Kulikova | Square-root algorithms for maximum correntropy estimation of linear discrete-time systems in presence of non-Gaussian noise[END_REF].

The MCC framework has already given very promising results, such as in machine learning with robust facial recognition [START_REF] He | Maximum Correntropy Criterion for Robust Face Recognition[END_REF] or in robust channel estimation for wireless communications [START_REF] Ma | Maximum correntropy criterion based sparse adaptive filtering algorithms for robust channel estimation under non-Gaussian environments[END_REF]. However, the question of the resilience of such algorithms is an open question which however gets more and more interest: in particular, [START_REF] Bako | Robustness analysis of a maximum correntropy framework for linear regression[END_REF] and [START_REF] Chen | Effects of Outliers on the Maximum Correntropy Estimation: A Robustness Analysis[END_REF] discussed this question under the scope of regression.

The goal of this paper is to give theoretical guarantees of resilience for a class of state estimators inspired by the kernel functions used in MCC-based approaches. The studied system model is given in the form of a LTI discretetime state-space representation affected by an arbitrary measurement noise. There is indeed no assumption on its statistical properties, on its values or on the distribution of its largest values over time and among sensors. The main theoretical result of this paper states that under a certain sufficient condition on the observability of the system, the estimation error induced by those estimators is upper bounded by a value which does not depend on the largest values of the measurement noise. The obtained bound is therefore valuable in order to understand how the characteristics of the system impact the performances of the estimator as it is directly linked to the parameters of the system itself. It is however non-computable. We nonetheless provide simulation results for the discussed class of estimators which show good performances in the presence of arbitrary noise and even when a process noise is added in the state equation, for which it was not designed in the first place.

Outline. The rest of the paper is organized as follows. The system to be estimated will be introduced in Section 2 while the class of estimators will be defined in Section 3. Section 4 will be focused on deriving theoretical guarantees on the boundedness of the estimation error. Finally, Section 5 will provide simulation results to assess the performances of the class of estimators, in the case without process noise in subsection 5.1 and in the case with process noise in subsection 5.2; Section 6 provides some concluding remarks.

Notations. Throughout this paper, R ≥0 (respectively R >0 ) designates the set of nonnegative (respectively positive) reals. We note R a the set of (column) vectors with a real elements and for any vector z in R a , z i with i in {1, ..., a} is the i-th component of z. Moreover, R a×b is 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 . • 2 is the Euclidean norm, defined by z 2 = √ z z for all z in R a . When applied to a matrix, • 2 will designate the matrix norm induced by the Euclidean norm. For a finite set S, the notation |S| will refer to the cardinality of S, S c the complementary set of S and P(S) will be the power set of S, i.e. the collection of all subsets of S.

STUDIED SYSTEM TO BE ESTIMATED

In this paper, we consider a Linear Time-Invariant discrete-time system Σ :

x t+1 = Ax t y t = Cx t + f t (1)
with A ∈ R n×n the state transition matrix of the system and C ∈ R ny×n the observation matrix of the system.

x t ∈ R n is the state of the system at time t ∈ Z + ; in particular x 0 ∈ R n is called the initial state, while y t ∈ R ny designates the output of the system. {f t } ⊂ R ny is a noise sequence which can take potentially any value. A convenient way of describing it is as follows:

f t = v t + s t ,
(2) where {v t } is a bounded white noise sequence whereas {s t } is a sparse noise sequence, meaning that most of its values are equal to zero but its non-zero values can be arbitrarily large. It is however important to understand that s t is sparse with regards to both time and sensors indexes, which implies that {s t } can present several nonzero values in a row on a specific sensor. The bounded white noise v t usually represents measurement noise in the output equation, while the sparse noise can represent many different things such as sensor failures, intermittent network outage or false data injection attacks.

Problem. In this paper, we consider the problem of estimating the state trajectory of system (1) on a fixed time horizon T = {0, 1, ..., T -1}, which means estimating the matrix X (x0

x 1 • • • x T -1) ∈ R n×T with the measurement matrix Y (y0 y 1 • • • y T -1)
containing T measurements and the model (1) of the system. The resulting estimated trajectory must be accurate even though the estimation process is undermined by the presence of f t whose characteristics are already described above.

THE CLASS OF STATE ESTIMATORS

To address the problem stated above, we propose an optimization-based solution. To do so, we define the following cost function:

V Σ (Y, z 0 ) = (t,j)∈T×J
e -γψ(yt,j -θ t,j z0)

(3) with J = {0, 1, ..., n y }, γ ∈ R >0 a user-defined parameter and z 0 ∈ R n . For every (t, j) ∈ T×J, θ t,j = c j A t where c j is the j-th row of the observation matrix C: every row vector θ t,j is a row of the observability matrix of system Σ over the time horizon

T, O = C (CA) • • • (CA T -1
) .

ψ : R → R is a real function which is assumed to verify the following properties:

(P1) Positive definiteness: ψ(0) = 0 and ψ(a) > 0 for all non-zero real a.

(P2) Symmetry: ψ(-a) = ψ(a) for all a ∈ R. (P3) Non-decreasingness: for any a 1 , a 2 in R, |a 1 | < |a 2 | implies ψ(a 1 ) ≤ ψ(a 2 ). (P4) Generalized Triangle Inequality (GTI): there exists α ∈ R >0 such that for all a 1 , a 2 in R, ψ(a 1 -a 2 ) ≥ αψ(a 1 ) -ψ(a 2 ). ( 4
)
Many functions verify this set of properties, such as every absolute value power function a → |a| p with p ∈ N (see [START_REF] Bako | Robustness analysis of a maximum correntropy framework for linear regression[END_REF]). It aims at generalizing the kernels used in Maximum Correntropy approaches (see (Chen et al., 2017a) for Gaussian kernel and [START_REF] Chen | Robust State Estimator Based on Maximum Exponential Absolute Value[END_REF] for exponential absolute value kernels). The estimator can therefore be defined as the set-valued map Ψ Σ : R ny×T → P(R n ) such that for any possible measurement matrix

Y in R ny×T , Ψ Σ (Y ) = arg max z0∈R n V Σ (Y, z 0 ). ( 5
)
We can note that this rather defines a class of estimators which has as many members as there are ψ functions verifying (P1)-(P4).

To obtain an estimate of the whole trajectory, we then need to simulate system (1) with any estimated initial state x0 ∈ Ψ Σ (Y ). From now on, the expression "estimated trajectories" will designate any matrix X(x 0 ) such that 6) With this estimation framework defined, our first goal is to provide a theoretical analysis of our estimator in order to exhibit its resilience properties.

X(x 0 ) = x0 Ax 0 • • • A T -1 x0 , x0 ∈ Ψ Σ (Y ). (

THEORETICAL RESULTS

In this section, we will conduct an analysis which will result in an upper bound on the estimation error independent from the largest values of the arbitrary noise sequence {f t }. This analysis is inspired by the one presented in [START_REF] Bako | Robustness analysis of a maximum correntropy framework for linear regression[END_REF] but in the state estimation framework. A discussion about the estimator performances on a system with process disturbances will also be provided.

Preliminaries

Before stating our main theoretical result, we will first introduce a few notations. For convenience, and without any loss of generality, we can assume that there is no (t, j) in T × J such that θ t,j 2 = 0 (1) . We then define the following variable:

σ Σ = min (t,j)∈T×J θ t,j 2 > 0 (7)
Moreover, given λ in [0; 1] and z 0 ∈ R n , we also define the following set of indexes (t, j):

J Σ (z 0 , λ) = {(t, j) ∈ T × J : |θ t,j z 0 | ≥ λ θ t,j 2 z 0 2 } (8)
This set selects the rows in the observability matrix which are almost in the same direction as z 0 : indeed, if (t, j) is in J Σ (z 0 , λ), then the vector θ t,j is within the cone of direction z 0 and of half-top angle arccos(λ). Thus, λ = 1 is the case where θ t,j and z 0 are colinear and λ = 0 is the case where θ t,j can be in any direction. As a result, |J Σ (z 0 , λ)| can be considered as a local observability measure with regards to z 0 given that it represents how many θ t,j are almost in the same direction as z 0 with a tolerance λ.

We also define a global observability parameter:

R Σ (λ) = inf z0∈R n |J Σ (z 0 , λ)|. (9)
For any λ in [0; 1], R Σ (λ) will be between 0 and n y T : in particular, if λ = 0, then R Σ (λ) = n y T . The higher R Σ (λ) will be for high values of λ, the more the system will be considered observable.

The following lemma provides more context about the link between R Σ (λ) and the observability of the system:

Lemma 1. Consider the system Σ defined in (1). This system is (i) observable, i.e. rank(O) = n, if and only if (ii) there exists λ ∈]0; 1] such that R Σ (λ) = 0.

Proof. (i) ⇒ (ii): if the system is observable, let's assume that for all λ in ]0; 1], R Σ (λ) = 0. As the set {|J Σ (z 0 , λ)|} z0∈R n is a subset of N, its infimum is necessarily attained, which entails that there exists

a λ ∈ R n such that J Σ (a λ , λ) = ∅ for any λ ∈]0; 1]. For a fixed λ, J Σ (a λ , λ) = ∅ implies that for all (t, j) in T × J, |θ t,j a λ | < λ θ t,j 2 a λ 2
. By squaring the two sides of the inequalities and adding them for every (t, j), we obtain

a λ O Oa λ < λ a λ 2 2 (t,j)∈T×J θ t,j 2 2 (10)
From [START_REF] Bernstein | Matrix Mathematics: Theory, Facts, and Formulas -Second Edition[END_REF] where λ e = σ min / (t,j)∈T×J θ t,j 2 2 > 0. This is a clear contradiction for any λ in ]0; λ e ]∩]0; 1], which proves the implication.

(ii) ⇒ (i): by contraposition, if the system is not observable, then there exists a ∈ R n = 0 such that Oa = 0. Given that Oa = θ 0,1 a • • • θ 0,ny a θ 1,1 a • • • θ T -1,ny a , we have for every (t, j) in T × J, θ t,j a = 0. As a result, for any λ ∈]0; 1], J Σ (a, λ) = ∅ since |θ t,j a| = 0 cannot be greater than λ θ t,j 2 a 2 : this leads to R Σ (λ) = 0 for any λ ∈]0; 1], which proves the implication. This lemma shows that the observability is an equivalent condition to the existence of a λ ∈]0; 1] such that R Σ (λ) = 0: Σ will thus be supposed to be observable from now on.

Finally, we need to define a last notation. Given ε a positive real number, for any noise sequence {f t,j }, it is possible to split the set T × J into two disjoint subsets,

I ε = {(t, j) ∈ T × J : |f t,j | ≤ ε},
(11) which gathers the indexes of f t,j such that their absolute value is smaller than ε, and I c ε = {(t, j) ∈ T×J : |f t,j | > ε} which consists of the indexes of outliers in {f t,j } with regards to ε. ε therefore acts like a threshold we can choose and tune in order to conduct our analysis.

Main result

In the following theorem, we are going to express under which circumstances the norm of the estimation error on the initial state, i.e. e 0 = x0x 0 , and consequently the estimation error over the whole trajectory E = X -X, can be bounded by a value which does not depend on the largest values of the noise sequence {f t }: Theorem 2. Consider the state estimator (5) for system (1) under the assumption that its loss function ψ defined in (3) verifies properties (P1)-(P4). Let ε > 0. For any noise sequence {f t } and initial state x 0 in (1), generating a measurement matrix Y , such that

1 1 + e -γψ(ε) R Σ (λ) + e -γψ(ε) |I ε | > n y T, ( 12
)
is verified for some λ in ]0; 1], the following holds true :

∀x 0 ∈ Ψ Σ (Y ), ψ(λσ Σ e 0 2 ) ≤ 1 γα ln (1/µ) ( 13
)
with e 0 = x0x 0 and

µ = 1 + e -γψ(ε) |I ε | + R Σ (λ) -n y T 1 1 + e -γψ(ε) R Σ (λ) +e -γψ(ε) |I ε | -n y T (14)
Moreover, if ψ is (strictly) increasing on R ≥0 , then it admits an invert function ψ (-1) and

N ( X -X) ≤ M Σ λσ Σ ψ (-1) 1 γα ln(1/µ) ( 15
)
with X as defined in (6), N ( X -X) = max t∈T xtx t 2 and M Σ a constant depending on the system dynamics.

Proof. By definition (5) of the estimator Ψ Σ , for any measurement matrix Y in R ny×T , we have

V Σ (Y, x 0 ) ≤ V Σ (Y, x0 ) (16) for every x0 in Ψ Σ (Y ). This yields (t,j)∈T×J e -γψ(ft,j ) ≤ (t,j)∈T×J e -γψ(yt,j -θ t,j x0) (17)
Left side sum in (17) will now be decomposed as T ×

J = I ε ∪ I c ε by definition (11) of I ε . If (t, j) is in I ε , then |f t,j | < ε, which implies ψ(f t,j ) ≤ ψ(ε) because of (P3),
and subsequently e -γψ(ft,j ) ≥ e -γψ(ε) . For any (t, j) in I c ε , we also have e -γψ(ft,j ) > 0, which yields

|I ε |e -ψ(ε) ≤ (t,j)∈T×J e -γψ(yt,j -θ t,j x0) (18)
Moreover, for all (t, j) in T × J, we always have ψ(y t,jθ t,j x0 ) = ψ(θ t,j x 0 + f t,jθ t,j x0 ) = ψ(f t,jθ t,j e 0 ) with e 0 = x0x 0 . For (t, j) in I ε , we apply the GTI (4) to ψ, leading to ψ(f t,jθ t,j e 0 ) ≥ αψ(θ t,j e 0 )ψ(f t,j ). For (t, j) in I c ε , as ψ verifies (P1), we have ψ(f t,jθ t,j e 0 ) ≥ 0, so for any (t, j) in I c ε , e -γψ(ft,j -θ t,j e0) ≤ 1. As a result, we obtain

|I ε |e -ψ(ε) ≤ e γψ(ft,j ) (t,j)∈Iε e -γαψ(θ t,j e0) + |I c ε | ⇔ e -γψ(ε) |I ε |(1 + e -ψ(ε) ) -n y T ≤ (t,j)∈Iε e -γαψ(θ t,j e0) as |I c ε | = n y T -|I ε
| and e γψ(ft,j ) < e γψ(ε) for (t, j) ∈ I ε . Given λ between 0 and 1, we will now decompose I ε depending on whether (t, j) belongs to J Σ (e 0 , λ) or not. If (t, j) is in I ε ∩J Σ (e 0 , λ), then ψ(θ t,j e 0 ) ≥ ψ(λ σ Σ 2 e 0 2 ). As a → e -a is decreasing, it yields

e -γψ(ε) |I ε |(1 + e -ψ(ε) ) -n y T ≤ |I ε ∩ J Σ (e 0 , λ)| e -γαψ(λσ Σ e0 2 ) -1 + |I ε | Moreover, we have |I ε ∩ J Σ (e 0 , λ)| = |I ε | + |J Σ (e 0 , λ)| -|I ε ∪ J Σ (e 0 , λ)| ≥ |I ε | + R Σ (λ)
n y T By choosing λ so that condition ( 12) is met, we obtain

|I ε | + R Σ (λ) -n y T > 0 given that |I ε | + R Σ (λ) ≥ 1 1 + e -γψ(ε) R Σ (λ) + e -γψ(ε) |I ε | > n y T.
(19) Since e -γψ(λσ Σ e0 2 ) -1 ≤ 0, we obtain

e -γψ(ε) |I ε |(1 + e -ψ(ε) ) -n y T -|I ε | ≤ (|I ε | + R Σ (λ) -n y T ) e -γαψ(λσ Σ e0 2) -1 (20) Dividing both sides by |I ε | + R Σ (λ)
n y T and simplifying the left hand side of the inequality then gives us µ ≤ e -γαψ(λσ Σ e0 2 )

(21) with µ as defined in ( 14). When condition ( 12) is met, µ is positive, so we can apply ln to both sides, yielding

ψ(λσ Σ e 0 2 ) ≤ 1 γα ln(1/µ) ( 22
)
If ψ is increasing on R ≥0 , then it is obviously invertible on that interval, which entails

e 0 2 ≤ 1 λσ Σ ψ (-1) 1 γα ln(1/µ) (23)
Finally, with X as defined in ( 6), xt = A t x0 and x t = A t x 0 , so we have xtx t = A t e 0 for all t ∈ T. Consequently,

N ( X -X) = max t∈T A t e 0 2 ≤ (max t∈T A t 2 ) e 0 2 . This eventually yields N ( X -X) ≤ max t∈T A t 2 λσ Σ ψ (-1) 1 γα ln(1/µ) (24)
which is the desired result with

M Σ = max t∈T A t 2 .
What this theorem states is that the estimator Ψ Σ yields estimates which are bounded despite the presence of an arbitrary noise potentially taking any value. Unfortunately, condition ( 12) cannot be computed due to the fact that obtaining R Σ (λ) is a combinatorial problem and that having access to |I ε | requires knowing the proportion of outliers in {f t,j } with regards to ε.

Nevertheless, it gives relevant information about what is important to ensure the resilience of the estimator. Indeed, the inequality is composed of three terms: on the left-hand side, there are two terms, one which depends on a quantitative observability of the system through R Σ (λ), one which depends on the number of reasonable measurements with regards to ε. The right-hand side there is composed of a constant term equal to the total number of measurements on the time horizon. To promote this condition, there are two important things:

• We need the observability criterion R Σ (λ) to be as large as possible • The number of outliers must be somehow limited The interpretation of condition ( 12) is that the estimation error is bounded if the system is observable enough, and the more it is observable, the more outliers the estimator is able to handle.

In addition, the bound itself in (15) gives information on what impacts the quality of the estimation: besides the conclusions obtained through the study of condition ( 12), we can see that the actual values of λ and ε play a role in the tightness of the bound. Given that it exists for any (λ, ε) which verify condition (12), the best bound obtained through Theorem 1 would be

N (E) ≤ M Σ min (λ,ε) µ>0 1 λσ Σ ψ (-1) 1 γα ln(1/µ) (25)

Case of systems with dynamic bounded noise

We now consider a system Σ w of the form Σ w :

x t+1 = Ax t + w t y t = Cx t + f t (26)
where {w t } is a bounded noise sequence. Adding a noise component w t in the state equation is a common way of modeling process disturbances which shifts the state from the dynamic induced by A. Even though our estimator is not designed for such systems, it can still be applied to Σ w . Indeed, for any system Σ w , there exists a system Σ such that

Σ : xt+1 = Ax t y t = C xt + ft (27)
with x0 = x 0 , ft = C wt + f t and wt = t-1 k=0 A k w k . This new system verifies the structure defined in (1) and gives the exact same output as Σ w . In addition, the gap between x t and xt is equal to x t -xt = A t x 0 + wt -A t x 0 = wt . As a result, we can draw the following corollary from Theorem 2: Corollary 3. Consider the state estimator (5) for system (1) under the assumption that its loss function ψ defined in (3) is invertible and verifies properties (P1)-(P4). Let ε > 0. For any noise sequence {f t } and initial state x 0 in (1), generating a measurement matrix Y , such that condition ( 12) is verified for some λ in ]0; 1], the following holds true :

∀x 0 ∈ Ψ Σ (Y ), N ( X(x 0 ) -X) ≤ M Σ λσ Σ ψ (-1) 1 γα ln(1/µ) + max t∈T wt 2 (28)
with X(x 0 ) as defined in (6), w as defined in ( 27) and µ, N and M Σ as defined in Theorem 2.

Proof. By applying Theorem 2 to system Σ defined in (27), we obtain

e 0 2 ≤ 1 λσ Σ ψ (-1) 1 γα ln(1/µ) (29)
In the case of the systems Σ w , we have for all t ∈ T, e t = xtx t = A t e 0 -wt , so by considering the norm N , we have

N (E) ≤ M Σ X -X 2 + max t∈T wt 2 (30)
which yields the desired result.

DISCUSSIONS ON THE IMPLEMENTATION OF THE ESTIMATOR

For our numerical tests, we considered a system (1) with the following parameters

A = 0.7 0.45 -0.5 1 , C = (1 2) , x 0 = 1 2 . ( 31
)
In this system, f t was decomposed as in (2): v t is a uniformly distributed white noise, while s t is a sparse vector. It is generated in two steps: (1) the indexes (t, j)

for which s t,j = 0 are uniformly selected over T × J, and then (2) for each (t, j) such that s t,j = 0, the value of s t,j is set as the realization of a Gaussian process of variance 100 and mean 0.

The optimization problem which defines the estimator is differentiable, but non-convex, which makes its numerical implementation difficult. To assess its resilience properties, we designed an Iterative Reweighted Least Squares (IRLS) algorithm (see Appendix A) to try to solve the optimization problem when ψ is the square function, i.e. ∀a ∈ R, ψ(a) = a 2 . We also chose γ = 1. Note however that there is no theoretical evidence as to the convergence of this iterative process towards the true solution of the (nonconvex) optimization problem in (5).

To perform the following tests, each setting was realized 100 times and the obtained relative estimation errors X -X 2 / X 2 were averaged. The IRLS estimator, implementing the estimator Ψ Σ discussed in this paper, and an Oracle Least Squares (OLS) estimator were implemented.

The mention "Oracle" means it works on a version of y t unaffected by sparse noise s t , namely y wt = y ts t . The OLS was chosen as a reference for comparison given that it is the best estimator with regards to the covariance of the estimation error in presence of Gaussian noises [START_REF] Geer | Least squares estimation[END_REF](Geer, , p. 1041) and should therefore provide good results with bounded noises. Additionally, to give a reference when the IRLS is presented in the case with process noise, we also implemented a 1 -norm based estimator designed for the presence of disturbance noise and for which the resilience was discussed in our previous work [START_REF] Kircher | Analysis of resilience for a State Estimator for Linear Systems[END_REF].

Sparsity test

First, we conduct a test in absence of w t to assess how the proportion of outliers were affecting the performances of the estimator with a fixed signal-to-noise ratio (SNR) of 30dB for v t . The obtained results are presented on Figure 1. As expected, the IRLS estimator maintains acceptable performances as the ratio of outliers in the sparse noise increases: its average relative estimation error is almost identical to the one of the OLS estimator until 70% of non-zero values in s t .

Performances of the estimator in presence of process noise

For the tests presented in this subsection we add a uniformly distributed white noise w t in the state equation of the system. We then conduct two tests, one to see how the presence of a process noise w t of 30dB degrades the performances of Ψ Σ with increasing ratio of outliers in f t (First test), and a second one to assess the impact of increasing SNR on the estimator (Second Test).

First test. The results obtained for the IRLS, 1 -norm based and OLS estimators are gathered in Figure 2. The performance of the IRLS estimator are logically worse than in the previous subsection as a process noise w t with a SNR of 40dB was introduced, and we see that the 1 -norm based estimator is better on ]0; 60] given that its relative estimation error is closer to the one of the OLS estimator. However, as the analysis suggested, the performances of the IRLS estimator are still acceptable, and its estimation error seems bounded until 80% of outliers. In particular, we have proven a resilience property for this class of state estimators so that the resulting estimation error is bounded by a value which is independent of the extreme values of the measurement noise provided that a certain condition linked to the observability of the system is met. Moreover, we proposed simulation results to assess the performance of this class, observing that it has good performances with and without the presence of a process noise. In future works, we would like to investigate the resilience property for an estimator defined by a cost function which would take into account the presence of process noise. Derivating theoretical guarantees for the designed IRLS algorithm would also be an interesting development.

Appendix = e -γ(yt,j -θ t,j x (i) 0 ) 2 (A.2) where x (i) 0 refers to the solution of (A.1) obtain at step i. At each step, Problem (A.1) is solved through the CVX Solver in MATLAB [START_REF] Grant | The CVX Users' Guide[END_REF]. When the relative difference between two consecutive estimated initial states is smaller than the user-defined threshold ε, the algorithm returns the state trajectory obtained from the last estimated state. The whole algorithm can be stated as follows:

Algorithm 1 Iterative Reweighted Least Squares (IRLS) Algorithm 1: Inputs: γ, Y , Σ, ε 2: Initialization: 3: ∀(j, t) ∈ J × T, k

(1) t,j ← e -γy 2 t,j 4:

x(1) 0 ← arg min z0∈R n V

(1)

Σ,c (Y, z 0 )

5: ∀(t, j) ∈ T × J, k

(2) t,j ← e -γ(yt,j -θ t,j

x(1) 0 ) 2 6: η ← 10 8 7: i ← 1 8: End of Initialization. 9: while η > ε do 10:

i ← i + 1 11:

x(i) 0 ← arg min z0∈R n V (i)

Σ,c (Y, z 0 )

12:

∀(t, j) ∈ T × J, k (i+1) t,j
← e -γ(yt,j -θ t,j

x(i) 0 ) 2 13:

η ←

x(i) 0 -x (i-1) 0

x(i-1) 0 14: end while 15: return X(x

(i) 0 )

Fig. 1 .

 1 Fig. 1. Relative Estimation Error with different ratios of corrupted values for s t and a SNR of 30dB for v t

  Second Test. In this test, the ratio of outliers is fixed at 10%, while the SNR of both w t and v t simultaneously vary from 10dB to 100dB. Figure3displays the relative estimation error obtained with those different SNR for the three estimators. Despite the presence of 10% of outliers in s t , the estimation error of both resilient estimators follow the same trend as the estimation error of the OLS one: it decreases as the SNR increases, which is logical given that noises are smaller compared to other quantities in the system equations. Nonetheless, once again, we observe that the IRLS estimator gives overall worse results than the 1norm based one, which was to be expected considering it is outside of the scope of its design. Its performances are however still in the same range as the 1 -norm based estimator. Invariant discrete-time systems in the face of uncertainties modeled as measurement noise in the system output equation. This noise sequence assumes values of possibly arbitrarily large amplitude which occur intermittently in time and across the available sensors. For this problem we have considered a class of state estimators through a cost function which generalizes the ones used in MCC-based approaches.

	It is even better than the 1 -norm based estimator on the interval [60; 90].
	Relative estimation error (log)	-2 -1 0 1	0	20 IRLS estimator 40 1 -norm based estimator 60 OLS estimator	80	100
				Ratio of sparse noise (in %)
	Fig. 2. Relative Estimation Error with different ratios of corrupted values for s t and a SNR of 30dB for both w t and v t
	Relative estimation error (log)	0 -1 -2 -3 -4 -5		OLS estimator IRLS estimator 1 -norm based estimator
			20	40	60	80	100
		Signal-to-noise ratio of bounded noises (dB)
	Fig. 3. Relative Estimation Error with varying SNR and 10% of corrupted values
				6. CONCLUSION
	In this paper, we have considered the problem of esti-mating the state of Linear Time-

(1) If O contains null rows, the following analysis can be conducted with a matrix collecting all the non-zero rows of O.