Zonotopes and Kalman observers: Gain optimality under distinct uncertainty paradigms and robust convergence

Automatica 55 (2015) 265–273

Contents lists available at ScienceDirect

Automatica journal homepage: www.elsevier.com/locate/automatica

Brief paper

Zonotopes and Kalman observers: Gain optimality under distinct uncertainty paradigms and robust convergence✩ Christophe Combastel 1 ECS-Lab (Electronics and Control Systems Laboratory, EA 3649) - ENSEA, 6, avenue du Ponceau, 95014 Cergy-Pontoise, France

article

info

Article history: Received 15 May 2014 Received in revised form 15 December 2014 Accepted 26 February 2015

Keywords: Observers Uncertainty Dynamic systems Kalman filters Intervals Robust stability State estimation Zonotopes

abstract State bounding observation based on zonotopes is the subject of this paper. Dealing with zonotopes is motivated by set operations resulting in simple matrix calculations with regard to the often huge number of facets and vertices of the equivalent polytopes. Discrete-time LTV/LPV systems with state and measurement uncertainties are considered. Based on a new zonotope size criterion called FW -radius, and by merging optimal and robust observer gain designs, a Zonotopic Kalman Filter (ZKF) is proposed with a proof of robust convergence. The notion of covariation is introduced and results in an explicit bridge between the zonotopic set-membership and the stochastic paradigms for Kalman Filtering. No intersection is used and the influence of the reduction operator limiting to a tunable maximum the size of the matrices involved in the zonotopic set computations is fully taken into account in the LMI-based robust stability analysis. A numerical example illustrates the effectiveness of the proposed ZKF. © 2015 Elsevier Ltd. All rights reserved.

1. Introduction State observation is often a key step for the satisfaction of advanced monitoring and/or control requirements in many engineering applications. Explicitly characterizing how the model and measurement uncertainties can influence the possible state values is also useful, not only when an automated decision-making is further required (e.g. fault diagnosis Ding, 2008), but also in some control frameworks (Efimov, Raïssi, & Zolghadri, 2013). Two paradigms can be used to model uncertainties: The stochastic one relies on probability theory and mainly deals with random variables. Usually, assumptions about their probability distribution are required. For instance, state estimators based on Kalman filtering (Kalman, 1960; Lewis, 1986; Maybeck, 1979; Sorenson, 1983) rely on covariance matrices to model usually Gaussian state and measurement random perturbations. Though well suited to take the

✩ Part of this work was done with the project MAGIC-SPS (Guaranteed Methods and Algorithms for Integrity Control and Preventive Monitoring of Systems) funded by the French National Research Agency (ANR) in accordance with the decision no ANR 2011 INS 006 02. The material in this paper was not presented at any conference. This paper was recommended for publication in revised form by Associate Editor Andrey V. Savkin under the direction of Editor Ian R. Petersen. E-mail address: [email protected]. 1 Tel.: +33 1 30 73 66 64.

http://dx.doi.org/10.1016/j.automatica.2015.03.008 0005-1098/© 2015 Elsevier Ltd. All rights reserved.

distribution of random noises into account, this modeling of uncertainty may be less representative when dealing with large disturbances mostly related to not well-known deterministic behaviors (e.g. load torque of a motor under incompletely specified operating conditions). The set-membership paradigm relying on unknown but bounded uncertainties (Schweppe, 1968) can lead to descriptions requiring no assumption about the probability distributions. Variants of ellipsoidal state sets have been proposed (Bertsekas & Rhodes, 1971; Kurzhanskiy & Valyi, 1997) as well as solutions addressing robustness issues (Petersen & Savkin, 1999). Interval analysis (Jaulin, Kieffer, Didrit, & Walter, 2001; Moore, 1966) has also led to state bounding algorithms, either based on set predictions/intersections (Chisci, Garulli, & Zappa, 1996; Combastel, 2003; Jaulin et al., 2001; Raïssi, Ramdani, & Candau, 2004) resembling Kalman filtering, or based on interval observers (Combastel, 2013; Gouzé, Rappaport, & Hadj-Sadok, 2000; Mazenc & Bernard, 2011; Raïssi, Efimov, & Zolghadri, 2012). The latter usually provide computationally efficient set-membership estimations with proven stability, but essentially deal with interval hulls which may be a rough enclosure of the consistent state sets. Though better taking domain shapes into account, single ellipsoids may hardly give tight enclosures of interval vectors and general polytopes (Ziegler, 1994) often suffer from the complexity of vertices/facets enumeration wrt the space dimension. This motivates the focus of this paper on zonotopes (Kühn, 1998; Ziegler, 1994), a class of polytopes

266

C. Combastel / Automatica 55 (2015) 265–273

whose shape is implicitly represented by a rectangular matrix. Basic operations over zonotopes often reducing to simple matrix calculations, zonotopic state bounding observers have already been proposed (Alamo, Bravo, & Camacho, 2005; Combastel, 2003, 2005; Le, Stoica, Alamo, Camacho, & Dumur, 2013; Montes de Oca, Puig, & Blesa, 2012). In Le et al. (2013), the so-called P-radius is introduced as a zonotope size criterion and permits to prove that zonotopic state estimates have non-increasing sizes in time, for timeinvariant systems under the assumption that the required complexity reduction operator (Kühn, 1998) has a negligible influence. By minimizing a matrix norm serving as a zonotope size criterion named F -radius (FW -radius if the norm is weighted) and by using the covariation of a zonotope which is a set-membership analog to covariance, a Zonotopic Kalman Filter (ZKF) is here derived from a discrete time-varying model subject to unknown but bounded uncertainties. ZKF computes minimal zonotopic sets enclosing all the admissible states. Explicit links between the set-membership and the stochastic paradigms for Kalman filtering are given. A robust stability analysis using Linear Matrix Inequalities (LMI) and fully taking the reduction operator into account is proposed. Compared to previous zonotopic observers using a singular value decomposition (svd), the complexity order of ZKF is significantly improved. The paper organization begins with preliminaries in Section 2 and the problem formulation in Section 3. The optimal observer gain minimizing the FW -radius of the predicted zonotopic state set is studied in Section 4. The ZKF algorithm is given in Section 5. After comparison with the stochastic Kalman Filter in Section 6, the robust stability is analyzed in Section 7. A numerical example illustrates the efficiency of ZKF in Section 8. 2. Preliminaries

⟨R⟩ = ⟨0, R⟩ is called centered zonotope. Any permutation of the columns of R leaves it invariant. The Minkowski sum of two sets S1 and S2 is S1 ⊕ S2 = {s1 + s2 , (s1 , s2 ) ∈ S1 × S2 }. The linear image of the set S ⊂ Rn by L ∈ Rq×n is L ⊙ S = {Ls, s ∈ S }. Zonotopes form a class of polytopic sets implicitly represented by matrices and leading to efficient set computations (Kühn, 1998; Ziegler, 1994). This class is closed under Minkowski sum ⊕ (computed from a matrix concatenation like [R1 , R2 ] in (5)) and linear image ⊙ (computed from a matrix product like LR in (6)):

⟨c1 , R1 ⟩ ⊕ ⟨c2 , R2 ⟩ = ⟨c1 + c2 , [R1 , R2 ]⟩, L ⊙ ⟨c , R⟩ = ⟨Lc , LR⟩. ⟨c , R⟩ ⊂ ⟨c , b(R)⟩, b(R) = diag (|R|1),

R = [r1 , . . . , rj , . . . , rp ],

↓q,W R = [R> , b(R< )] ∈ Rn×q ,

∂X f (X ) is a short notation for ∂ f (X )/∂ X . If the function f returns scalar values and X = [Xij ] is a matrix, then ∂X f (X ) = [∂Xji f (X )].  tr (.) denoting the trace of a square matrix (tr (A) = i Aii ), and X , A, B, C being matrices of appropriate size, it comes:

3. Problem formulation

tr (A) = tr (AT ),

(1)

xk+1 = Ak xk + Bk uk + Ek vk ,

(2)

yk = Ck xk + Dk uk + Fk wk ,

(3)

x0 ∈ ⟨c0 , R0 ⟩ ⊂ R .

tr (AB) = tr (BA),

∂X tr (AX B) = A B , T

∂X tr (AXBX C ) = BX CA + B X A C . T

∥rj ∥2W ≥ ∥rj+1 ∥2W ,

T

T

T

T

T

The matrix M is spd if it is symmetric (M = M T ) and positive definite (M ≻ 0 i.e. ∀x ̸= 0, xT Mx > 0). Then, M −1 exists and is spd. Similar definitions and properties hold for the negative case. The Schur complement of C in L = [A, B; BT , C ] is S = A − BC −1 BT and L ≻ 0 ⇔ A ≻ 0 ∧ S ≻ 0. Let R = [r1 , . . . , rp ] ∈ Rn×p be a matrix with p non zero columns and M ∈ Rn×n . For any such R, if M is spd, then RT MR is spd and tr (RT MR) > 0. Similarly, M ≺ 0 ⇒ tr (RT MR) < 0. Let W ∈ Rn×n be a spd matrix and R ∈ Rn×p , ∥R∥F ,W =

tr (RT WR) is the weighted Frobenius norm of R. ∥R∥F , obtained for W = In (identity: no weight), is invariant under orthogonal trans-



formations. Also: ∥R∥2F ,W =

p

i =1 n

∥ri ∥2W , where ∥ri ∥W =

riT Wri

is a weighted vector norm in R . 2.2. Zonotopes

A zonotope ⟨c , R⟩ ⊂ Rn with center c ∈ Rn and generator matrix R ∈ Rn×p is a polytopic set defined as the linear image of the unit hypercube [−1, +1]p ⊂ Rp by R:

⟨c , R⟩ = {c + Rs, ∥s∥∞ ≤ 1}.

(4)

(8)

If p ≤ q, Then ↓q,W R = R, Else R> = [r1 , . . . , rq−n ],

T

(6) (7)

(7) shows how a zonotope ⟨c , R⟩ can be enclosed within an aligned box (or interval hull) defined by b(R) ∈ Rn×n , where | · | is the element-by-element absolute value operator, 1 is a column vector of ones, and diag (.) returns a diagonal matrix from a vector of diagonal elements. Such a box enclosure usually being too conservative, a reduction operator can be used: Kühn (1998). In this work, a weighted version of the reduction operator first introduced in Combastel (2003) and improved in Combastel (2005) will be denoted ↓q,W , where q ≥ n specifies the maximum number of columns of matrix ↓q,W R satisfying the inclusion property ⟨R⟩ ⊂ ⟨↓q,W R⟩. The reduction operator ↓q,W first consists in sorting the columns of R ∈ Rn×p on decreasing weighted vector norm ∥ · ∥W (8) and enclosing the set ⟨R< ⟩ generated by the p − q + n smaller columns into a box:

2.1. Matrix calculus (Boyd, El Ghaoui, Feron, & Balakrishnan, 1994; Golub & Van Loan, 1996)

T

(5)

R< = [rq−n+1 , . . . , rp ].

(9) (10)

The zonotopic state observation of systems modeled by uncertain time-varying discrete-time dynamics as in (11)–(12) is the main subject of this paper:

nx

vk ∈ ⟨0, Inv ⟩,

(11)

wk ∈ ⟨0, Inw ⟩,

(12) (13)

nM (resp. pM ) denoting the number of lines (resp. columns) of any matrix M, xk ∈ Rnx , uk ∈ Rnu , yk ∈ Rny , vk ∈ Rnv , wk ∈ Rnw respectively refer to the state, the known input, the known output (measurement), the state evolution uncertainty, and the measurement uncertainty vectors. ∀k ≥ 0, Ak , Bk , Ek , Ck , Dk , Fk denote non empty matrices with appropriately fixed size. All the related matrix sequences, e.g. A: = (A0 , . . . , Ak , . . .), k ∈ N, are not necessarily known a priori, but each matrix with time index k is assumed to be known at time k, possibly expressed as a matrix function, e.g. A(.) : Rnα → Rnx ×nx , of a known (measured) scheduling parameter vector αk ∈ Rnα so that Ak = A(αk ), Bk = B(αk ), etc. Notice that a robust detectability assumption further stated in Section 7 will be required to prove the convergence of the proposed set-based observer. The initial state x0 is assumed to belong to a zonotope ⟨c0 , R0 ⟩ specified by its center c0 ∈ Rnx and non empty generator n ×p matrix R0 ∈ R x R0 (13), (4). ∀k ≥ 0, vk ∈ [−1, +1]nw = ⟨0, Inv ⟩ where Inv is the identity matrix of size nv × nv i.e. vk is bounded by a unit hypercube expressed as a centered zonotope in (11), and idem for wk in (12). For the sake of simplified notations, the discrete time index k can be omitted and the index + replaces k + 1 in the following.

C. Combastel / Automatica 55 (2015) 265–273

4. Optimal gain of zonotopic observer 4.1. Observer structure Intuitively, the proposed set-based observer can be viewed as a Kalman Filter where the usually Gaussian probability density functions (pdf) would have been replaced by zonotopic sets. More precisely, the considered zonotopic observer structure is obtained by induction after first assuming that the state vector x satisfies the inclusion property x ∈ ⟨c , R⟩ at a given time k ≥ 0, which is true at k = 0 (13). The introduction of a time-varying matrix gain G in (11)–(12) gives: x+ = Ax + Bu + E v + G(y − Cx − Du − F w).

(14)

Though leaving x+ unchanged, G changes the modeling of uncertainty propagation by weighting the relative influences of the state prediction uncertainty vk (bounded state noise) and the measurement uncertainty wk (bounded measurement noise). G can thus be viewed as an observer gain giving more freedom degrees than prediction alone (G = 0) to improve the state estimation i.e. to reduce the state estimation error. In the next paragraph 4.2, a set size criterion will be used to quantify the latter. Proposition 1 (Observer Structure). Given a system modeled as in (11)–(12) with x0 satisfying (13), let us recursively define c (state center estimate) and R (shape matrix of state bounding zonotope) as in (15)–(16): c+ = (A − GC )c + (B − GD)u + Gy, R+ = [(A − GC )R¯ , E , −GF ],

R¯ = ↓q,W R.

∥R1 ∥2F ,W + ∥R2 ∥2F ,W .

Definition 3 (F -Radius). The F -radius of the zonotope ⟨c , R⟩ ⊂ Rn is the Frobenius norm of R: ∥⟨c , R⟩∥F = ∥R∥F . It equals the FW radius for W = In . Definition 4 (Covariation). The covariation of the zonotope ⟨c , R⟩ is cov(⟨c , R⟩) = RRT . So, minimizing the F -radius of a zonotope like ⟨c+ , R+ ⟩ (15)–(16) is equivalent to minimizing the trace of its covariation P+ = R+ RT+ . Indeed, J = tr (P+ ) = ∥R+ ∥2F (1). By analogy, in the weighted case, minimizing the FW -radius of the zonotope ⟨c+ , R+ ⟩ is equivalent to minimizing the criterion JW = tr (WP+ ) = ∥R+ ∥2F ,W also expressed as a (weighted) function of the covariation. As for the FW -radius, JW = J for W = In . 4.3. Optimal observer gain Some covariation matrices are first introduced: P = RRT ,

(16)

P¯ = R¯ R¯ T = cov(↓q,W R).

(17)

Proof. Assuming x ∈ ⟨c , R⟩ at time k (true at k = 0) and ↓q,W

preserving inclusion, x ∈ ⟨c , R¯ ⟩. Since v ∈ ⟨0, Inv ⟩ and w ∈ ⟨0, Inw ⟩ (11)–(13), from (14), it comes: x ∈ ⟨c , R⟩ ⇒ x+ ∈ ⟨c+ , R+ ⟩ = ((A − GC ) ⊙ ⟨c , R¯ ⟩) ⊕

((B − GD) ⊙ ⟨u, 0⟩) ⊕ (G ⊙ ⟨y, 0⟩) ⊕ ((−GF ) ⊙ ⟨0, Inw ⟩) ⊕ (E ⊙ ⟨0, Inv ⟩).

2 Since ∥R∥2F ,W = i=1 ∥ri ∥W (see end of Section 2.1), the FW -radius is a zonotope size criterion involving the norm of all the generator segments i.e. all the columns ri of R as in (8). The weighting spd matrix W makes the FW -radius benefit from similar advantages as the P-radius (Le et al., 2013) in terms of LMI-based robust stability analysis. Also, the Minkowski sum of zonotopes (5) is closely related to the FW -radii as: ∥⟨c1 , R1 ⟩⊕⟨c2 , R2 ⟩∥2F ,W = ∥[R1 , R2 ]∥2F ,W =

(15)

Then, the state inclusion property (17) holds ∀k ≥ 0: x ∈ ⟨c , R⟩.

267

p

(18)

Applying zonotope properties (5)–(6) show that c+ and R+ in (18) are the same as in (15) and (16). Therefore, x+ ∈ ⟨c+ , R+ ⟩. This gives the proof by induction.  (15)–(16) define the structure of the considered zonotopic observer: given a time-varying matrix gain G, their explicit computation provides state estimates with guaranteed zonotopic bounds: the observation error is x − c ∈ ⟨0, R⟩. When several iterations are considered, the otherwise increasing number of columns pR in (16) is solved by using the reduction operator ↓q,W (8)–(10). Indeed, ∀k > 0, pR ≤ q + pE + pF . 4.2. Optimality criterion The observer structure in Proposition 1 raises the issue of recursively computing an optimal observer gain G∗ at time k. A precision criterion related to the size of the next state bounding zonotope ⟨c+ , R+ ⟩ = ⟨c+ (G), R+ (G)⟩ (15)–(16) can be used to that purpose. The FW -radius, a zonotope size criterion consistent with the norm used in the reduction operator ↓q,W , is thus considered: Definition 2 (FW -Radius). Let W ∈ Rn×n be a spd matrix: W = W T ≻ 0. The weighted Frobenius radius (FW -radius) of the zonotope ⟨c , R⟩ ⊂ Rn is the weighted Frobenius norm of R: ∥⟨c , R⟩∥F ,W = ∥R∥F ,W .

Qv = EE T ,

Qw = FF T ,

(19) (20)

Given x ∈ ⟨c , R⟩ (enclosure at time k) and given a weighting spd matrix W , an optimal observer gain G∗ minimizing the uncertainty about x+ can be obtained by minimizing the FW -radius of ⟨c+ , R+ ⟩: Theorem 5 (Optimal Observer Gain). Given R+ as in (16) and W = W T ≻ 0, the optimal observer gain G∗ = arg minG ∥R+ ∥2F ,W is given by (21)–(22). It does not depend on W . Therefore, G∗ minimizes not only the FW -radius but also the F -radius of ⟨c+ , R+ ⟩. In other words, G∗ also minimizes J = tr (P+ ) i.e. the trace of the covariation of ⟨c+ , R+ ⟩ ∋ x+ . G∗ = AK ∗ ,

¯ , L = PC T

K ∗ = LS −1 ,

¯ + Qw . S = C PC T

(21) (22)

Proof. JW = ∥R+ ∥2F ,W being convex wrt G, G∗ is the value of G

such that ∂G JW = 0. From (1), JW = tr (WP+ ) with P+ = R+ RT+ =

(A − GC )P¯ (A − GC )T + Qv + GQw GT according to (16), (19), (20). So, ∂G JW = ∂G tr (W (A − GC )P¯ (A − GC )T + WQv + WGQw GT ). Using (1), the symmetry of P¯ and W , and removing the additive terms independent on G in the argument of tr (.), ∂G JW = 0 gives ∂G tr (WGSGT ) = 2∂G tr (WALGT ) where L and S are as in (22). Replacing X by G in (3) (resp. (2)) and thanks to the symmetry of S, the equality rewrites as (WAL)T = S T GT W which yields (21) after transposition and left multiplication by W −1 which exists since W is spd.  Theorem 5 shows the independence of the optimal observer gain G∗ wrt to the weighting matrix W . Therefore, W can be viewed as a spd matrix parameter that can vary freely on the discrete time scale k without influencing the computation of G∗ . The robust stability analysis (Section 7) heavily relies on this. 5. Zonotopic Kalman Filter (ZKF) The analogy between K ∗ (21)–(22) and the Kalman gain is striking. In the same time, the observer structure (15)–(16) explicitly

268

C. Combastel / Automatica 55 (2015) 265–273

gives a zonotope enclosing the state (17) and, equivalently, the observation error x − c ∈ ⟨0, R⟩. The proposed Zonotopic Kalman Filter (ZKF) is obtained by choosing G∗ (21)–(22) as optimal observer gain G. Theorem 5 shows that both the FW -radius ∥R+ ∥F ,W and the F -radius ∥R+ ∥F of ⟨c+ , R+ ⟩ ∋ x+ (or ⟨0, R+ ⟩ ∋ (x+ − c+ )) are then minimized. As a byproduct of such minimization of uncertainty terms, c+ stands for the best point estimate of the state x+ (at time k + 1) using the measurements y until time k. The explicit computation of c+ , R+ and the covariation matrix P+ , which prevents from computing the matrix product RRT at each iteration, results in the ZKF algorithm (24)–(31):

(P˜ = (I − K ∗ C )P¯ ),

(27)

Let us consider stochastic (Lewis, 1986; Maybeck, 1979) rather than bounded-error uncertainties in the problem formulation stated in Section 3. Denoting x ∼ N (c , P ) a Gaussian random vector x centered on c with covariance P, this can be done by substituting vk ∼ N (0, Inv ), wk ∼ N (0, Inw ), x0 ∼ N (c0 , P0 ) for vk ∈ ⟨0, Inv ⟩, wk ∈ ⟨0, Inw ⟩, x0 ∈ ⟨c0 , R0 ⟩ in (11), (12), (13), respectively. Taking Inv (resp. Inw ) as covariance for v (resp. w ) is not restrictive since it is always possible to find E (resp. F ) for other covariance values than identity. The random vectors x0 as well as v (state noise) and w (measurement noise) at each time step are all assumed to be mutually independent. The involved Gaussian processes then follow the usual Kalman Filter (KF) assumptions. By introducing an observer gain G as in (14), it comes: x ∼ N (c , P ) ⇒ x+ ∼ N (c+ , P+ ) where c+ is as in (15) and P+ = (A − GC )P (A − GC ) + Qv + GQw GT with Qv = EE T and Qw = FF T . Searching for the optimal observer gain minimizing J = tr (P+ ) (here, trace of covariance) to reduce the imprecision around c+ leads to the discretetime KF (35)–(38) in delayed form (or one-step ahead predictor form) which is easy to embed in digital control loops:

E˘ = [E , −G F ],

(28)

ε = y − (Cc + Du),

ZKF iteration: Inputs : c , R, P (= RRT ), u, y, α, A, B, C , D, E , F , q, (W , G) = RobDetect(α, . . .), (G : not used) R¯ = ↓q,W R, P¯ = P − R< RT< + b(R< )2 ,

(23)

Qv = EE T ,

(25)

Qw = FF T ,

ε = y − (Cc + Du), ∗

K = LS

−1

,

G = AK , ∗

∗

¯ T, L = PC

S = CL + Qw ,

c˜ = c + K ε, ∗

A˘ = A − G C , ∗

∗

c+ = Ac˜ + Bu, R+ = [A˘ R¯ , E˘ ],

˜ + Qv ), (or P+ = APA T Outputs : c+ , R+ , P+ (= R+ R+ ). T

P+ = A˘ P¯ A˘ + E˘ E˘

T

T

(24)

(26)

r = |R|1.

˜ + Qv , P+ = APA T

P˜ = (I − K C )P¯ . ∗

−1

c˜ = c + K ε,

(31)

c+ = Ac˜ + Bu,

(32)

(33) (34)

˜ T + Qv with Proof. Substituting G∗ (21) in (33) yields P+ = APA ˜P = (I − K ∗ C )P¯ (I − K ∗ C )T + K ∗ Qw K ∗T = (I − K ∗ C )P¯ + ∆K ∗T and ¯ T + K ∗ Qw = −L + K ∗ S = 0 (21).  ∆ = −(I − K ∗ C )PC

S = CPC T + Qw ,

,

(30)

Corollary 6. The optimal gain (Kalman-like gain) given in Theorem 5 yields two equivalent expressions of P+ : A˘ = A − G∗ C ,

T

K = PC S

The zonotopic shape matrix R however captures more information about guaranteed state bounds than a single interval vector or ellipsoid would usually do (e.g. see Fig. 3). Under constant W (then, no need for (23)), the complexity of the ZKF iteration (Appendix D) is mainly that of a zonotopic prediction plus that of a Kalman Filter for discrete-time LTV systems. In (26), ε is an innovation term. The first expression of P+ in (31), where E˘ E˘ T = Qv + G∗ Qw G∗T , does not require P˜ (27). The second expression of P+ is valid only for G = G∗ (optimal gain value) whereas the first one remains consistent with R+ (30) for any observer gain G. Indeed,

P+ = A˘ P¯ A˘ T + Qv + G∗ Qw G∗T ,

∗

(29)

In (24), ↓q,W R results from (8)–(10) and P¯ from (20). RobDetect (23) aims at computing a (possibly constant) weighting matrix W influencing the sorting of generator segments in ↓q,W R. W interprets as a Lyapunov matrix relaxing the conditions required to prove the stability of ZKF compared to W = Inx (non-weighted case). Numerical methods for such design of W and RobDetect will be discussed in Sections 7 and 8. ↓q,W permits to fix the number of columns of R to a maximum along the iterations in spite of the concatenation in (30), which is the classical scheme of zonotopic prediction (Combastel, 2003; Kühn, 1998). A careful analysis shows that the ZKF iteration (24)–(31) computes c+ and R+ as in Proposition 1 (15)–(16) with the observer gain G taken equal to the optimal value G∗ given in Theorem 5. So, the inclusion property (17) is satisfied and (7) gives a state bounding interval vector: x ∈ c ± r = [c − r , c + r ] ⊂ Rnx ,

6. Comparison with stochastic Kalman Filter

∗

(35) (36)

P˜ = (I − K C )P , ∗

˜ T + Qv . P+ = APA

(37) (38)

(35)–(37) and (38) respectively stand for the update/correction step and the prediction step. c = xˆ k/k−1 is the state estimate at time k using measurements until k − 1, and P = Pˆ k/k−1 is its covariance. ε is the innovation term and S is its covariance (35). K ∗ is the optimal Kalman gain (36). c˜ = xˆ k/k is the updated state estimate using yk , and P˜ = Pˆ k/k is its covariance (37). c+ = xˆ k+1/k is the state prediction at the next time step and P+ = Pˆ k+1/k is its covariance. Theorem 7 (ZKF vs. KF). Based on the bounded-error problem formulation (Section 3) leading to ZKF (23)–(31) and its stochastic counterpart leading to KF (35)–(38), If c0 is chosen equal for ZKF and KF, if the initial covariation P0 = R0 RT0 in ZKF equals the initial covariance P0 in KF, and if q → ∞ (which yields no over-approximation induced by ↓q,W i.e. R¯ = R,

P¯ = P and no need for W in ZKF), Then, ∀k ≥ 0, in spite of distinct uncertainty paradigms (boundederror vs stochastic), the iterations of ZKF and KF (i) give the same Kalman gain K ∗ , (ii) return the same punctual state estimate c, (iii) compute equal covariations (P in ZKF) and covariances (P in KF). Proof. By induction using (9), (23)–(31) to simplify ZKF, and (35)– (38) to compare with KF.  Additional facts related to the comparison between ZKF and KF can be stated as follows: (i) With ZKF, the probability distribution of uncertainties within the bounds can be arbitrary (not specified). At the additional computational cost of a zonotopic prediction (24), (30), ZKF provides zonotopic state sets enclosing all the possible values consistent with the bounded-error model and the past measurements. (ii) With KF, the probability distribution of uncertainties is specified as Gaussian. No worst-case bounding domains can be obtained because of the infinite support of Gaussian distributions but, given some probability level, ellipsoidal confidence domains can be computed.

C. Combastel / Automatica 55 (2015) 265–273

269

(iii) Theorem 7 underlines the analogy between covariation with ZKF (zonotopic bounds) and covariance with KF (Gaussian probability distributions). (iv) Infinitely many distinct centered zonotopes share2 a same covariation matrix whereas, given some confidence level, a centered confidence ellipsoid is uniquely determined by a covariance matrix.

Theorem 10. Let R ∈ Rn×p be the generator matrix of a zonotope ⟨R⟩. Let W ∈ Rn×n be a spd matrix with all its eigenvalues in [λ, λ] ⊂ R, λ > 0. Let n ≤ q < p and R¯ = ↓q,W R = [R> , b(R< )] ∈ Rn×q as

7. Robust stability analysis of ZKF

∥R¯ ∥2F ,W − ∥R∥2F ,W ≤

7.1. Robust detectability assumption

Proof. See Appendix B.

Lemma 8. The time-varying discrete-time linear system xk+1 = Ak xk is γ -stable (i.e. stable with decay rate γ ) if there exists a bounded sequence of spd matrices Wk and γ ∈ ]0, 1] ⊂ R+ such that, ∀k ≥ 0,

i.e. R¯ k has never more than q columns. With no loss of generality, let us consider pR¯ k = q (otherwise, ↓q,W induces no over-

 γ Wk ⋆

ATk Wk+1 Wk+1



≻ 0.

(39)

Proof. In (39), ⋆ denotes a symmetric term (here, WkT+1 Ak ). The Schur complement of (39) gives γ Wk − ATk Wk+1 Wk−+11 WkT+1 Ak ≻ 0. So, ∀k ≥ 0, ∀xk ̸= 0, γ xTk Wk xk − xTk ATk Wk+1 Ak xk > 0. Choosing Vk = xTk Wk xk as Lyapunov function, Vk+1 < γ Vk .  Remark. In the time-invariant case and γ = 1 (no requirement for a decay rate), Lemma 8 reduces to the Schur stability of A (verifiable through a LMI in W ).  The system xk+1 = Ak xk is said robustly γ -stable wrt the polytopic matrix sequence set [Ak ] if it is γ -stable for all state matrix sequence consistent with [Ak ] i.e. ∀Ak ∈ [Ak ], ∀k ≥ 0. Robust stability then refers to robust 1-stability (limit case of γ = 1). LMI-based numerical methods as outlined in Appendix A give constructive solutions to design observer gain sequences Gk and Lyapunov matrix sequences Wk satisfying the robust detectability Assumption 9, where (1 − γ ) > 0 can be viewed as a detectability margin: Assumption 9. Given any two bounded sequences of matrices Ak ∈ [Ak ] and Ck ∈ [Ck ], k ≥ 0, where [Ak ] and [Ck ] refer to two sets of bounded matrix sequences, it is assumed that a bounded sequence of observer gains Gk ∈ [Gk ] and a bounded sequence of spd matrices Wk ∈ [Wk ] (both possibly depending on the known scheduling parameter vector αk ) have been designed such that xk+1 = (Ak − Gk Ck )xk is robustly γ -stable wrt [Ak ] and [Ck ], with γ < 1, and with Wk as sequence of Lyapunov matrices having all their eigenvalues within some constant real interval [λ, λ], λ > 0. The robust design in Assumption 9 allows the implementation of function RobDetect (23) returning Wk and Gk based on the available data at time k. Such data possibly include the value of a scheduling vector αk as in (23) to efficiently compute online Wk and, optionally, Gk from the results of an offline robust design. Indeed, Gk is not used online, but further required to prove convergence.

in (9), i.e. ⟨R¯ ⟩ is a reduced zonotope with at most q generator segments satisfying ⟨R⟩ ⊂ ⟨R¯ ⟩ ⊂ Rn . Let d = p − q ∈ N and µ = ((λ/λ) (d + n) − 1)(d + n). Then,



 µ ∥R∥2F ,W . d+q

(40)

The reduction operator ↓q,W (8)–(10) in (24) ensures ∀k, pR¯ k ≤ q

approximation since pR¯ k < q ⇒ R¯ = R according to (9)). (30) and

E˘ k in (28) give ∀k, (nR+ = nx ) ∧ (pR+ = q + pE + pF ). Then, Theorem 10 can be applied to Rk , ∀k, under fixed n = nRk , p = pRk and d = p − q = pE + pF , which also fixes µ. Thus, based on (40), ∀k, the relative error on the FW -radius of ⟨ck , Rk ⟩ induced by ↓q,W in (24) can be made arbitrarily small simply by choosing a sufficiently large value for the complexity tuning parameter q ∈ N. 7.3. Convergence result about the computed zonotopes by merging robust and optimal designs Theorem 11. ZKF (23)–(31) is considered with (11)–(13) and under the robust detectability Assumption 9. Let Ek ∈ Rnx ×pE and Fk ∈ Rnx ×pF be two bounded sequences of matrices with fixed size and k ∈ N. Then,

∀k ≥ 0 ,

xk ∈ ⟨ck , Rk ⟩,

(41)

∃q0 , ∀q ≥ q0 , ∃β ≥ 0, ∀k ≥ 0, ∀k > 0,

Rk 2F ,Wk

∥ ∥

≤ β,

pRk ≤ q + pE + pF .

(42) (43)

In words, ZKF (23)–(31) ensures the inclusion of the unknown state n ×p xk within the computed zonotope ⟨ck , Rk ⟩ (41), where Rk ∈ R x Rk . All the computed zonotopes are generated by a bounded number of line segments, pRk , which can be tuned through the complexity reduction parameter q ∈ N (43). Moreover, a value q0 ∈ N exists such that choosing q ≥ q0 ensures ⟨ck , Rk ⟩ is a sequence of zonotopes with bounded shape (i.e. with bounded FW -radius) as k → ∞ (42), for any possible time-varying pair (Ak , Ck ) as specified in Assumption 9. The proof of Theorem 11 and, in particular, the proof of (42) heavily relies on a second theorem giving more insight into the convergence properties of ZKF: Theorem 12. Under the same assumptions as in Theorem 11, a ¯ be an upper bound bounded sequence ψk satisfying (44) exists. Let ψ of ψk as k → ∞, and let µ be defined as in (40) for Rk ∈ Rnx ×(q+d) with d = pE + pF .

∀k ≥ 0 ,

(∥Ek ∥2F ,Wk+1 + ∥Gk Fk ∥2F ,Wk+1 ) ≤ ψk ,

γ¯ , γ (1 + µ/(d + q)),

ϕk , ∥Rk ∥2F ,Wk ,

(44) (45)

7.2. Influence of the reduction operator

ϕk+1 ≤ γ¯ ϕk + ψk ,

In order to further study the robust convergence of ZKF, the zonotopic outer approximation provided by the reduction operator ↓q,W is first quantified in Theorem 10.

Provided q is large enough, γ¯ < 1 can always be achieved since γ < 1 in (45). Then, (46) both ensures that the FW -radius ϕk of ⟨ck , Rk ⟩ is bounded as k → ∞, and gives an evaluation of its upper bound when k → ∞.

¯ 1 − γ¯ ). ϕ∞ ≤ ψ/(

(46)

Proof. See Appendix C. 2 This is a consequence of the additional freedom degrees granted by the usually rectangular zonotope shape matrices.

The proof of Theorem 11 concludes this paragraph. It gathers the three main properties of ZKF (23)–(31):

270

C. Combastel / Automatica 55 (2015) 265–273

(i) The inclusion property (41) finds its roots in (15)–(17). (ii) The tunable bounded complexity of the computed zonotopes (43) results from ↓q,W , (30), and E˘ k in (28). (iii) Under the robust detectability Assumption 9, the ability of ZKF to compute zonotopic state estimates having a bounded shape (42) for all k (and, in particular, as k → ∞) directly follows from Theorem 12. 8. Numerical example To illustrate ZKF (23)–(31), a time-varying discrete-time system modeled as in (11)–(13) is considered with: Ak = A(αk ) = A(0) + A(1) α1,k + A(2) α2,k , (0)

(1)

(A , A , A ) = 0.01 · . . .    67 −21 −76 24 10 21 67 −24 −76  3 10 −3 105 −33 ,  0 3

10



33

−7

0 0 0 0

0 0



0 0 Bk =   , 0 1

0

 −6 19  , 0 

   −0.2 1 0 0.4  1 0 + ξ, −0.5 1 0 k 0.1 1 0    0 0 , Dk = ,

1 1 Ek =  1 1

0.1 0

−0.2

0.1 Fk = 0

0 , 0.1

0 0

0



10 0 0

19 6 0 0

0

8 0  , −11 0

 



105

(48)

−3



0 0 0 0

−1

Ck =

(47)

(2)

1

0

(50)



⟨c0 , R0 ⟩ = [−40, +40]4 :

(49)

c0 = 04 , R0 = 40I4 .

(51)

ξk = 5(1 + sin(k/5)) for 100 < k ≤ 150 and ξk = 0 otherwise (49). Observer gain scheduling under constant Lyapunov matrix W is considered here to satisfy the robust detectability Assumption 9. LMI constraints derived from Proposition 13 by using a polytopic enclosure of the known scheduling vector αk ∈ [α, α] ⊂ R2 are like:



γW ⋆

A(α)T Y − C T Z (α) ≻ 0, Y + YT − W



where γ = 0.7 (< 1), α = −α = [0.5, 2]T , α successively takes the value of each of the 4 vertices of the polytopic enclosure [α, α], A(α) = A(0) + A(1) α1 + A(2) α2 and C = Ck as in (48), (50). The unknowns are: the spd matrix W ∈ Rnx ×nx , Y ∈ Rnx ×nx , Z (i) ∈ Rny ×nx for i = 0, 1, 2 with Z (α) = Z (0) + Z (1) α1 + Z (2) α2 . In addition to these four LMIs, two LMIs like Inx ≺ W ≺ λInx can be used to balance the eigenvalues of W . An implementation under Matlab making use of Yalmip (Löfberg, 2004) and MPT (Herceg, Kvasnica, Jones, & Morari, 2013) returned (up to rounding errors): 6.2 −5.2 W = 2.8 2.0

−5.2 7.8 0.8 −3.5



5.0 −4.5 Y = 1.8 1.7



Z (α) = Z

(0)

G(α) = (Y

−3.4 6.7 1.4 −3.1 +Z

(1)

2.8 0.8 12.7 −5.9 2.8 0.2 10.1 −5.8

α1 + Z

(2)

T



(52)

0.6 −0.3 , 0.1  8.1



α2 ,

(0)

) Z (α) = G

−1 T

2.0 −3.5 , −5.9 9.1

matrix sequences are bounded. Assumption 9 is so ensured with [λ, λ] = [1, 18]. Moreover, W being chosen constant, (23) reduces to an off-line assignment (52) in the implementation of ZKF. The main iteration is then given by (24)–(31). The ‘‘true’’ system to be observed is simulated using (11)–(12) with x0 = [10, 0, 10, −10]T and wk computed as a uniform random noise in [−1, +1]2 . Two scenarios, namely v = 0 and v ̸= 0 are reported in Figs. 1–3 and Table 1. The former (v = 0) is obtained from vk = 0 for all k. The latter (v ̸= 0) results from vk = [v1,k , v2,k ]T ∈ R2 where v1,k is a periodic square signal with its first period (12 samples) defined as −1 for k = 0, . . . , 5, and +1 for k = 6, . . . , 11. Also, v2,k = −v1,k . Thus, ∀k ≥ 0, vk ∈ [−1, +1]2 . Both scenarios are simulated with q = 100 and with the same values of uk , αk and ξk . αk = [0, 0]T (resp. α, α ) for k ∈ [0, 50] ⊂ N (resp. k ∈ ]50, 100], k ∈ ]150, 300]). For k ∈ ]100, 150], αk varies linearly from α to α . Notice ξk (49) is also varying (non linearly) in time for k ∈ ]100, 150]. Scalar intervals bounding the ith state xi,k of xk are computed ∀k as [xi,k , xi,k ] = ci,k ± ri,k (32) for i = 1, . . . , nx (Figs. 1 and 2). Table 1 reports the rms values of several sequences N 1 1 2 1/2 placed in columns. For instance, rms(e: ) = ( N + k=0 nx ∥ek ∥2 ) 1 nx with N = 300, ek = xk − ck ∈ R , nx = 4, evaluates the overall observation error for each scenario i.e. v = 0/v ̸= 0 and for each correction scheme among the next four ones: 0: null observer gain (prediction only, leading to diverging envelops), G: robust observer gain (54)–(55), G∗ : optimal gain of ZKF (28), svd: intersection-based correction using a singular value decomposition as in Combastel (2003). Except for rms(e: ), no distinction is made between v = 0/v ̸= 0, since rms values are identical in both cases. The state estimation is significantly improved (both in terms of observation error and bounding domain size) with the optimal gain G∗ , compared to the robust gain G designed to prove the convergence. This interpretation of Table 1 is confirmed by Fig. 1 for v = 0, and Fig. 2 for v ̸= 0. Indeed, G∗ quickly contracts the initial state uncertainty and leads to tight state bounding intervals (Fig. 2) and zonotopes (Fig. 3), so showing a good control of the wrapping effect by ZKF. The size improvement of the estimated state sets granted by zonotopic domain shapes compared to interval hulls (boxes) is quantified by the lower values of ∥R: ∥F compared to ∥b(R: )∥F in Table 1. In this study, the correction based on svd shows very similar state bounding domains compared to G∗ , but at the price of a significant increase in runtime (more than 4 times) coming from a complexity in O (q3 ) instead of O (q log q) (Appendix D). The average runtime of an iteration of ZKF was ≃ 0.26 ms using an Intel(R) Core(TM) i7-3517U CPU @ 1.90 GHz after setting the Matlab task affinity (Win8) to a single thread of execution. The first numerical experiments give encouraging results to further implement ZKF on real-time targets. 9. Conclusion

such that : (1)

(G(0) , G(1) , G(2) ) = · · · (55)     0.209 0.006 −1.693 −0.701 −5.477 −1.346 −0.639 0.172  , ,  1.038 0.302  −0.012 −0.003 −0.004 −0.001 −0.403 1.207   0.471 0.108 0.003   0.010  0.060 −0.107 . −0.001 0.000 Then, ∀k ≥ 0, provided αk ∈ [α, α], taking Ak = A(αk ), Zk = Z (αk ), Ck = C , Wk = Wk+1 = W and Yk = Y , Proposition 13 ensures that the system x˜ k+1 = (Ak − Gk Ck )˜xk with observer gain Gk = (Zk Yk−1 )T is γ -stable with Wk as Lyapunov matrix sequence. Thanks to the affine dependence of A(.) and Z (.) on the bounded scheduling vector αk and the invertibility of Y , all the involved

(2)

+ G α1 + G α2 ,

(53) (54)

By merging optimal and robust observer gain designs, a Zonotopic Kalman Filter (ZKF) is proposed to compute optimal

C. Combastel / Automatica 55 (2015) 265–273

271

Table 1 Comparison of different correction schemes (rms values).

0 G G∗ svd

e : = x : − c:

r:

∥b(R: )∥F

∥R: ∥F

∥R: ∥F ,W

3.223 / 50.02 1.407 / 6.239 0.654 / 3.295 0.685 / 3.312

250 14.9 5.98 6.07

500 29.8 11.97 12.13

79.5 13.2 8.46 8.50

218.4 25.37 20.24 20.30

Runtime for 300 iterations: ≃ 78 ms (except svd: ≃ 350 ms).

Fig. 2. Scenario: v ̸= 0. Interval estimate of third state x3,k given by [x3,k , x3,k ], either computed with the optimal gain G∗k (black) or the robust gain Gk (magenta). Also, ‘‘true’’ state x3,k (blue) and its punctual estimate i.e. center c3,k (red) obtained with G∗k . Top: whole time range. Bottom: zoom at beginning (left) and at end (right) of simulation. The state bounds computed with G∗k (black) tightly enclose the influence of unknown but bounded vk on x3,k (blue). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

foundations to jointly deal with epistemic and random uncertainties. Another connection is the joint use of optimality and robust stability: this paper addresses set-membership state observation, but the dual(?) problem of constraint finite-time optimal control has also been the subject of recent works (Besselmann, Löfberg, & Morari, 2012). Much remains to be done to exploit these connections. Appendix A. Outline of some numerical methods to satisfy the robust detectability assumption Proposition 13. Given a bounded sequence of matrix pairs (Ak , Ck ) ∈ Rn×n × Rm×n for all k ≥ 0, if there exists a bounded sequence of spd matrices Wk ∈ Rn×n , bounded sequences of matrices Yk ∈ Rn×n , Zk ∈ Rm×n , and γ ∈ ]0, 1] ⊂ R+ such that, ∀k ≥ 0,

Fig. 1. Top: input uk and first output y1,k for scenarios v = 0 (blue) and v ̸= 0 (green). Bottom (v = 0): Interval estimate of first state x1,k given by [x1,k , x1,k ], either computed with the optimal gain G∗k (black) or the robust gain Gk (magenta). Also, ‘‘true’’ state x1,k (blue) and its punctual estimate i.e. center c1,k (red) obtained with G∗k . (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

zonotopic state estimates while giving a proof of robust convergence for LPV systems. The robustness wrt the influence of the complexity reduction operator is taken into account. Natural extensions include dealing with some classes of non-linear systems which could result in an Extended ZKF. Notice that the current formulation readily allows the enclosure of bounded non-linearities since the set-membership paradigm is not subject to any assumption about the distribution of values within the uncertainty domains. By comparison, an explicit characterization of the evolution (prediction) of covariance seems more difficult in a probabilistic context with non-linearities/non-Gaussian distributions. Merging the strengths of both paradigms is however appealing in practical applications and the explicit connection made in this paper gives

 γ Wk ⋆

ATk Yk − CkT Zk Yk + YkT − Wk+1



≻ 0.

(A.1)

Then, the (closed-loop) system xk+1 = (Ak − Gk Ck )xk with observer gain Gk = (Zk Yk−1 )T is γ -stable with Wk as bounded Lyapunov matrix sequence. Proof. State feedback/observer gain design duality and steps analog to the proof of Theorem 3 in de Oliveira, Bernussou, and Geromel (1999).  Numerical methods for the design of sequences of observer gains and Lyapunov matrices satisfying Assumption 9 can rely on Proposition 13. For instance, based on a polytopic representation of system matrices, the analogy between the left term in (A.1) and the terms of the sum in (19) in Oliveira and Peres (2009) gives a constructive way to compute a gain scheduled observer gain G(αk ) ensuring the robust γ -stability of A(αk ) − G(αk )C (αk ) with W (αk ) as returned Lyapunov matrix sequence, possibly under fixed bounds for the variations of the scheduling vector αk . Proposition 13 can also serve as a basis for simpler LMI-based robust observer gain designs (Section 8).

272

C. Combastel / Automatica 55 (2015) 265–273

Fig. 3. Scenario: v ̸= 0. 2D projections of the 4D zonotopic state estimate ⟨c300 , R300 ⟩ ∋ x300 at time k = 300. State estimation with optimal gain G∗k (resp. robust gain Gk ): blue (resp. magenta) sets. Yellow stars: 2D projections of the ‘‘true’’ state xk which is close to the zonotope boundaries because of the significant excitation vk ̸= 0 in this scenario. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Appendix B. Proof of Theorem 10

• After permutation of its columns to sort them on decreasing order, ∥rj ∥2W ≥ ∥rj+1 ∥2W for j = 1, . . . , (p − 1), R rewrites as [R> , R< ] ∈ Rn×p , R> = [r1 , . . . , rq−n ] ∈ Rn×(q−n) and R< = [rq−n+1 , . . . , rp ] ∈ Rn×(p−q+n) , leaving ⟨R⟩ invariant. ∥R∥2F ,W =  ∥R> ∥2F ,W + ∥R< ∥2F ,W , ∥R> ∥2F ,W = qj=−1n ∥rj ∥2W ≥ (q − n)∥rq−n ∥2W  p 2 2 and ∥R< ∥2F ,W = j=q−n+1 ∥rj ∥W ≤ (p − q + n)∥rq−n ∥W , the inequalities following from the sorting order. Then, (q − n)∥R< ∥2F ,W ≤ (p − q + n)∥R> ∥2F ,W , (q − n)∥R< ∥2F ,W ≤ (p − q + n)(∥R∥2F ,W − ∥R< ∥2F ,W ), ∥R< ∥2F ,W ≤ ((p − q + n)/p)∥R∥2F ,W , and (B.1) comes: ∥R< ∥2F ,W ≤ ((d + n)/(d + q))∥R∥2F ,W . ∥b(R< )∥

2 F ,W

≤ (λ/λ)(d + n)∥

R< 2F ,W

∥

(B.1)

.

(B.2)

• Proof of (B.2). ∀x ∈ Rn , λxT x ≤ xT Wx ≤ λxT x i.e. λ∥x∥22 ≤ ∥x∥2W ≤ λ∥x∥22 . Thus, ∥b(R< )∥2F ,W ≤ λ∥b(R< )∥2F , ∥R< ∥2F ≤ (1/λ)∥R< ∥2F ,W , and (B.2) follows after proving that ∥b(R< )∥2F ≤ (d + n)∥R< ∥2F . Given any matrix M, let us denote Mi: (resp. M:j ) the ith line (resp. jth column) of M. From b(.) in (7) and the equivalence between the 1 and 2 vector norms giving ∀x ∈ Rm , ∥x∥21 ≤ m∥x∥22 , n n 2 2 it comes ∥b(R< )∥2F = i=1 ∥(R< )i: ∥1 ≤ (d + n) i=1 ∥(R< )i: ∥2 since (R< )i: is a (d + n)-dimensional (line) vector with d = p − q. n 2 2 Then, noticing that i=1 ∥(R< )i: ∥2 = ∥R< ∥F completes the proof of (B.2).

• The third step consists in jointly using (B.1) and (B.2) to prove (40) in Theorem 10. Since R rewrites as [R> , R< ] (after columns permutation) and R¯ = [R> , b(R< )], ∥R∥2F ,W = ∥R> ∥2F ,W + ∥R< ∥2F ,W and ∥R¯ ∥2F ,W = ∥R> ∥2F ,W +∥b(R< )∥2F ,W . From the last two equalities, ∥R¯ ∥2F ,W − ∥R∥2F ,W = ∥b(R< )∥2F ,W − ∥R< ∥2F ,W . Then, by successively applying (B.2) and (B.1), (40) follows.



Appendix C. Proof of Theorem 12

• E: , F: (resp. G: , W: ) are bounded sequences of time-varying matrices: Theorem 11 (resp. Assumption 9) and ∀x ∈ Rn , ∥x∥2W+ ≤

λ∥x∥22 . So, ∀k, ∥E ∥2F ,W+ + ∥GF ∥2F ,W+ ≤ λ(∥E ∥2F + ∥GF ∥2F ), and an explicit (not depending on k + 1) bounded sequence ψ: satisfying (44) exists thanks to the boundedness of E: , F: and G: . • Let us focus on (46). The FW+ -radius of ⟨c+ , R+ ⟩ in ZKF is: ∥R+ ∥2F ,W+ = ∥[(A − G∗ C )R¯ , E , −G∗ F ]∥2F ,W+ (30) and the optimal gain G∗ (28) minimizes ∥R+ ∥2F ,W+ whatever the spd matrix W+ is

(Theorem 5). Thus, taking any gain G instead of G∗ , like, for instance, G given by the robust design in Assumption 9 (2nd output of RobDetect (23)) gives: ∥R+ ∥2F ,W+ ≤ ∥[(A − GC )R¯ , E , −GF ]∥2F ,W+

i.e. ϕ+ ≤ ∥(A − GC )R¯ ∥2F ,W+ + ψ , so merging optimal and robust gain designs. The γ -stability of (A−GC ) with γ < 1 (Assumption 9) and the Schur complement of (39) give: (A − GC )T (W+ )(A − GC ) − γ W ≺ 0. Then, using M ≺ 0 ⇒ tr (R¯ T M R¯ ) < 0 (see 2.1) and linearity of tr, tr (R¯ T (A − GC )T (W+ )(A − GC )R¯ ) < γ tr (R¯ T W R¯ ), ∥(A − GC )R¯ ∥2F ,W+ < γ ∥R¯ ∥2F ,W , and ϕ+ ≤ γ ∥R¯ ∥2F ,W +ψ . Substituting ϕ for

∥R∥2F ,W in (40) gives ϕ+ ≤ γ¯ ϕ + ψ as in (45)–(46) or, equivalently:   µ ∥R+ ∥2F ,W+ ≤ γ 1 + ∥R∥2F ,W + ψ. d+q

¯ 1 − γ¯ ) in (46), with no loss of generality, • To prove ϕ∞ ≤ ψ/( ¯ as k → ∞ is considered. Also, a sequence ψ: such that ψk → ψ let ϕ¯ be recursively defined as: ϕ¯ + = γ¯ ϕ¯ + ψ, ϕ¯ 0 = ϕ0 . Then, the positive dynamic of (ϕ¯ − ϕ) ensures ∀k, ϕ ≤ ϕ¯ . Moreover, ¯ 1 − γ¯ ) (steady-state of ϕ¯ ), which completes the proof ϕ¯ ∞ = ψ/( ¯ 1 − γ¯ ) as k → ∞.  that ϕ is upper bounded by ψ/( Appendix D. Complexity analysis Using (Hunger, 2007) to count the number of flops (floating point operations) of elementary matrix operations and K ∗ involving the inverse of spd matrix S, the number of flops in an iteration of ZKF under constant W , (24)–(31), is: ny − nw − nx − nv − q + 2nu nx + 2nu ny + 4nv nx

+ 3nw nx + nx q + 6nv n2x + 4nw n2x + 2nw n2y + 6nx n2y + 8n2x ny + 4n2x q + 3n2x + 6n3x + n3y + 2nw nx ny ,

(D.1)

C. Combastel / Automatica 55 (2015) 265–273

with a possible decrease of nx (2nx − 1)(nv + nw + q) flops if W = Inx which simplifies the sorting criterion in ↓q,W . (D.1) is linear in q, so a sorting algorithm with worst-case complexity in O (q log q) in (8) imposes the overall worst-case complexity order in q. However, profiling the execution shows that sorting usually requires less than 30% of runtime of ↓q,W , even for large q. The cubic terms n3x , n3y in (D.1) come from covariation computations and would also appear with covariances in KF. Notice also that the complexity of svd(R¯ ) (full svd required in Combastel, 2003) involves a term in O (q3 ) (Golub & Van Loan, 1996), so explaining the runtime improvement with ZKF in Table 1. For the numerical example in Section 8, the flops count (D.1) is 67q + 1222 with nx = 4, nu = 1, ny = nv = nw = 2, while ⟨R⟩ ⊂ R4 is a zonotope with up to 364 208 facets.3 References Alamo, T., Bravo, J. M., & Camacho, E. F. (2005). Guaranteed state estimation by zonotopes. Automatica, 41(6), 1035–1043. Bertsekas, D. P., & Rhodes, I. B. (1971). Recursive state estimation for a setmembership description of uncertainty. IEEE Transactions on Automatic Control, 16(2), 117–128. Besselmann, T., Löfberg, J., & Morari, M. (2012). Explicit MPC for LPV systems: stability and optimality. IEEE Transactions on Automatic Control, 57(9), 2322–2332. Boyd, S., El Ghaoui, L., Feron, E., & Balakrishnan, V. (1994). SIAM studies in applied mathematics, Linear matrix inequalities in system and control theory. Philadelphia: SIAM. Chisci, L., Garulli, A., & Zappa, G. (1996). Recursive state bounding by parallelotopes. Automatica, 32(7), 1049–1055. Combastel, C. (2003). A state bounding observer based on zonotopes. In Proc. of European control conference, ECC. Combastel, C. (2005). A state bounding observer for uncertain non-linear continuous-time systems based on zonotopes. In 44th IEEE conf. on decision and control, and European control conf., CDC-ECC (pp. 7228–7234). Combastel, C. (2013). Stable interval observers in C for linear systems with timevarying input bounds. IEEE Transactions on Automatic Control, 58(2), 481–487. de Oliveira, M. C., Bernussou, J., & Geromel, J. C. (1999). A new discrete-time robust stability condition. Systems & Control Letters, 37(4), 261–265. Ding, S. (2008). Model-based fault diagnosis techniques: design schemes, algorithms, and tools. Springer. Efimov, D., Raïssi, T., & Zolghadri, A. (2013). Control of nonlinear and LPV systems: interval observer-based framework. IEEE Transactions on Automatic Control, 58(3), 773–778. Golub, G. H., & Van Loan, C. F. (1996). Matrix computations (3rd ed.). Johns Hopkins University Press. Gouzé, J. L., Rappaport, A., & Hadj-Sadok, Z. (2000). Interval observers for uncertain biological systems. Ecological Modelling, 133, 45–56. Herceg, M., Kvasnica, M., Jones, C. N., & Morari, M. (2013). Multi-parametric toolbox 3.0. In Proc. of the European control conference, Zürich, Switzerland, July 17–19 (pp. 502–510). http://control.ee.ethz.ch/~mpt.





3 = 2 q + nv + nw . nx − 1

273

Hunger, R. (2007). Floating point operations in matrix–vector calculus. Technical report. Technische Universität München, September. https://mediatum.ub.tum.de/doc/625604. Jaulin, L., Kieffer, M., Didrit, O., & Walter, E. (2001). Applied interval analysis. Springer. Kalman, R. E. (1960). A new approach to linear filtering and prediction problems. Transactions of the ASME—Journal of Basic Engineering, 82, 35–45. Kühn, W. (1998). Rigorously computed orbits of dynamical systems without the wrapping effect. Computing, 61(1), 47–67. Kurzhanskiy, A., & Valyi, I. (1997). Ellipsoidal calculus for estimation and control. Boston, MA: Birkhäuser. Le, V. T. H., Stoica, C., Alamo, T., Camacho, E. F., & Dumur, D. (2013). Zonotopic guaranteed state estimation for uncertain systems. Automatica, 49(11), 3418–3424. Lewis, R. (1986). Optimal estimation with an introduction to stochastic control theory. John Wiley & Sons, Inc.. Löfberg, J. (2004). YALMIP: a toolbox for modeling and optimization in MATLAB. In Proceedings of the CACSD conference, Taipei, Taiwan. Maybeck, P. S. (1979). Stochastic models, estimation and control. Academic Press. Mazenc, F., & Bernard, O. (2011). Interval observers for linear time-invariant systems with disturbances. Automatica, 47, 140–147. Montes de Oca, S., Puig, V., & Blesa, J. (2012). Robust fault detection based on adaptive threshold generation using interval LPV observers. International Journal of Adaptive Control and Signal Processing, 26(3), 258–283. Moore, R. E. (1966). Interval analysis. Englewood Cliffs, NJ: Prentice-Hall. Oliveira, R. C. L. F., & Peres, P. L. D. (2009). Time-varying discrete-time linear systems with bounded rates of variation: stability analysis and control design. Automatica, 45(11), 2620–2626. Petersen, I. R., & Savkin, A. V. (1999). Control engineering, Robust Kalman filtering for signals and systems with large uncertainties. Boston: Birkhäuser. Raïssi, T., Efimov, D., & Zolghadri, A. (2012). Interval state estimation for a class of nonlinear systems. IEEE Transactions on Automatic Control, 57(1), 260–265. Raïssi, T., Ramdani, N., & Candau, Y. (2004). Set membership state and parameter estimation for systems described by nonlinear differential equations. Automatica, 40(10), 1771–1777. Schweppe, F. C. (1968). Recursive state estimation: unknown but bounded errors and system inputs. IEEE Transactions on Automatic Control, 13(1), 22–28. Sorenson, H. (1983). Special issue on applications of Kalman filtering. IEEE Transactions on Automatic Control, 28(3), 253–434. Ziegler, G. M. (1994). Graduate texts in mathematics: Vol. 152. Lectures on polytopes. Springer-Verlag.

Christophe Combastel received his M.Sc. degree in Engineering (1997) and his Ph.D. in Control Systems (2000), both from the National Polytechnic Institute of Grenoble (G-INP), France. Since 2001, he is Associate Professor at ENSEA (Ecole Nationale Supérieure de l’Electronique et de ses Applications, Cergy, France) where he was responsible of the control system department from 2002 to 2005. He is also with the Electronics and Control Systems Laboratory (ECS-Lab, EA3649) where his current research interests include interval and set-membership algorithms for integrity control applications ranging from on-line fault diagnosis to verified model-based design.