Ellipsoidal constrained state estimation in presence of bounded disturbances
aa r X i v : . [ ee ss . S Y ] D ec Ellipsoidal constrained state estimation inpresence of bounded disturbances
Yasmina BECIS-AUBRY ∗ December 8, 2020
Abstract
This contribution proposes a recursive, input-to-state stable and ready-to-use online algorithm for the state estimation of linear discrete-time sys-tems with unknown but bounded disturbances corrupting both the stateand the sporadic measurement vectors and subject to linear inequalityand equality constraints on the state vector.Two set representation techniques are used: the state vector is char-acterized by an ellipsoid whereas the disturbances vectors are bounded bypossibly degenerate zonotopes. The proposed algorithm is decomposedinto two steps: time and observation updating that uses a switching esti-mation gain.
There is no more need to praise the interests of the set-membership state es-timation techniques neither is there a necessity to recall how interesting alter-native they offer to conventional state estimation methods where the statisticalassumptions on the disturbances can not be satisfied in certain practical situa-tions, nor how increasing attention they are currently receiving since the noisesof their models are assumed only to be bounded. Nevertheless, the stabilityquestion is rarely addressed in this kind of estimation approach.On the other hand, the constrained state filtering has been widely studiedin stochastic context [Sim10], [DL13], [JZ13]. In [AIBS19], constrained Kalmanfilter variations were reexamined and an alternative derivation of the optimalconstrained Kalman filter for time variant systems was proposed. The litera-ture is less abundant on this subject when it comes to bounded error framework.LMI techniques were employed for ellipsoidal set-membership constrained statefiltering with linear [YL09a] and linearized nonlinear [YL09b] equalities. In[NBH15], the authors used a combined stochastic and set-membership uncer-tainty representation by integrating, into the Kalman filter structure, ellipsoidalconstraints on the state vector as a relaxation of equality constraints. All theseworks faced a same challenge, not arising here, in inverting the estimation er-ror covariance matrix, becoming inevitably singular, when dealing with equalityconstraints. ∗ *The author is with Universit´e d’Orl´eans, Laboratoire PRISME EA 4229 (Univ. Orl´eans- INSA CVL). 63 av. de Lattre de Tassigny, 18020 Bourges Cedex, FRANCE. Tel. +33 2 4823 84 78 [email protected] .
1n our early paper [BABD08], we presented a state bounding estimation al-gorithm for linear discrete-time systems, where the state and all the (processand measurement) disturbances were characterized by multidimensional ellip-soids. In order to guarantee the input-to-state stability of the estimation errorand the size decrease of the state bounding ellipsoid at the measurement updat-ing stage, a polynomial equation had to be solved, at each time step, involvingthe computation of the eigenvalues and eigenvectors of a matrix having thesame dimension as the state vector. This issue was partially solved in [SLZ + − norm of somequantities of interest, resulting typically in ellipsoidal bounding sets that arenothing else that balls or hyperspheres undergoing rotations and scalings and2. those bounding the ∞− norm of such quantities, leading to intervals (whichare boxes or hypercubes), parallelotopes (skewed, stretched, or shrunken imagesof such boxes) or zonotopes (flattened images of boxes, which are also general-ization of parallelotopes). The drawbacks of using exclusively the ellipsoid asbounding set for both state and disturbances vectors were highlighted above.Now, when using intervals, the recourse to interval computing softwares includ-ing time costly operations such as subpavings and contractors are inevitableto overcome the conservatism of such aligned with the coordinates axes boxes[JKDW12], [RJ15]. And when dealing with zonotopes to characterize the outerbound of the set of all possible values of the state vector, it is necessary to usesome tools such as LMI to manage the growth of the number of generators,inherent to the zonotopes summing, during the time update, and to their inter-section, during the measurement correction ( cf. [Com15] and references within).This is why we chose two different bounding techniques: the ellipsoids (basedon the 2 − norm) to characterize the set of all possible values of the state vec-tor at each time step and the zonotope (based on the ∞− norm) to bound thedisturbances. Moreover, we are interested here in the state estimation of lineardiscrete-time systems subject not only to bounded process and measurementdisturbances but also to all kinds of linear constraints applied to the state vec-tor, i.e. , equalities (modelled by hyperplanes) and inequalities (represented bypolyhedrons and zonotopes); all, noises and constraints, manifesting themselvessporadically, not at all time steps.The paper is organized as follows. After this introduction, which is com-pleted by some notations and definitions, in the second section, the constrainedset-membership state estimation problem with sporadic measurements is formu-lated. The third section concerns the time-prediction stage, while the correctionstage of the estimation algorithm is detailed in forth one. Its properties and sta-bility are studied in the fifth section. Numerical simulations are presented inthe sixth and finally, a brief conclusion terminates the paper.2 otations and definitions
1. The symbol := ( resp . =: ) means that the Left Hand Side ( resp . RHS) isdefined to be equal to the Right Hand Side ( resp . LHS). Normal lowercaseletters are used for scalars, capital letters for matrices, bold lowercaseletters for vectors and calligraphic capital letters for sets. IR, IR ∗ , R + denote the sets of real, non-zero and nonnegative numbers resp . IN andIN ∗ are the sets of nonnegative and positive integers resp . l, m, n, p ∈ INdesignate vectors and matrices dimensions. The subscript k ∈ IN is thediscrete time step and i, j ∈ IN ∗ are vector and matrix component indices.2. x i is the i th component of the vector x . a ij is the i th row and j th columnelement of A := (cid:2) a ij (cid:3) = [ a j ] ∈ IR n × m and a j ∈ IR n is its j th column vector(if n = 0 or m = 0, A is an empty matrix).3. n ∈ IR n and 0 m × n ∈ IR n × m are zero vector and zero matrix resp. and I n ∈ IR n × n is the identity matrix.4. A T , A † , rank( A ), Ker ( A ) and R ( A ) stand resp . for the transpose, Moore-Penrose inverse, rank, kernel and range of the matrix A . If A is square,tr( A ), | A | = det( A ) and A − , are its trace, determinant and inverse (ifany) resp .5. Diag( x i ) i ∈{ , ··· ,k } is a diagonal matrix where x , . . . , x k are its diagonalelements.6. A Symmetric matrix M is Positive Definite, denoted by SPD or M > resp . Positive Semi-Definite or non-negative definite, denoted by SPSDor M ≥
0) if and only if ∀ x ∈ IR n – { } , x T M x > resp . x T M x ≥ resp . non-negative). The matrix inequality M > N ( resp . M ≥ N ) means that M − N > resp . M − N ≥ k x k := k x k := √ x T x is the 2-norm of the vector x ; k A k , := X j k a j k and k A k := k A k := σ max ( A ).8. B np := { z ∈ IR n | k z k p ≤ } is a unit ball in IR n for the p − norm. B n and B n ∞ := [ − , n are the centred unit hypersphere and hypercube/box resp .9. S ⊕ S := { x ∈ IR n | x = x + x , x ∈ S , x ∈ S } is the Minkowski sumof the sets S , S ⊂ IR n and ⊕ mi =1 S i := S ⊕ · · · ⊕ S m .10. E ( c , P ) := { x ∈ IR n | ( x − c ) T P − ( x − c ) ≤ } is an ellipsoid in IR n , where c ∈ IR n is its center and P ∈ IR n × n is a SPD matrix that defines its shape,size and orientation in the IR n space. It can be also viewed as an affinetransformation of matrix M (where M T M = P ) of the unit Euclideanball B n : E ( c , M T M ) = { x ∈ IR n | x = c + M z , z ∈ B p } . If M is not SPDbut only SPSD , then the ellipsoid is degenerate. It has an empty interiorin the case where p < n .11. H ( d , a ) := { x ∈ IR n | x T d = a } is a hyperplane in IR n of normal vec-tor d ∈ IR n and whose signed distance from the origin is a k d k . Letalso G ( d , a ) := { x : x T d ≤ a } be one of the two halfspaces into which3he hyperplane divides the IR n space and G ( − d , − a ) is the other one.Now let D ( d , a ) := G ( d , a ) ∩ G ( − d , − a ), i.e. , D ( d , a ) := { x ∈ IR n : (cid:12)(cid:12) x T d − a (cid:12)(cid:12) ≤ } which is the strip of IR n , of width 2 k d k − , that can also beseen as a degenerate unbounded ellipsoid or zonotope centred at H ( d , a ). P ( C, d ) = T mi =1 G ( c i , d i ) is a polyhedron.12. Z ( c , L ) := { x ∈ IR n | x = c + L z , z ∈ B m ∞ } = ⊕ qj =1 { t j l j , | t j | ≤ } ⊕ { c } is azonotope of center c , obtained by affine transformation, of shape matrix L ∈ IR n × m , of the unit box B m ∞ , where m can be smaller, equal to orgreater than n . A zonotope is also a convex polyhedron with centrallysymmetric faces in all dimensions. Besides its vertex representation suit-able for geometrical sum, there is a halfspace representation, suitable forintersection: Z H ( D, a ) = T i D ( d i , a i ).13. The support function of a set S ⊂ IR n is ρ S : IR n → IR, x ρ S ( x ) := sup u ∈S u T x . H (cid:0) x , ρ S ( x ) (cid:1) is the supporting hyperplane of S and S ⊂ G (cid:0) x , ρ S ( x ) (cid:1) . ρ E ( c ,P ) ( x ) = c T x + √ x T P x cf. [Che94]. Consider the following linear discrete time system x k = A k − x k − + B k − τ k − + R k − ν k − , (1a)where x ∈ E ( ˆ x , ς P ) =: E ⊂ IR n and ν k ∈ B m ∞ , (1b)where x k ∈ IR n , τ k ∈ IR l and ν k ∈ IR m are the unknown state vector to be es-timated, a known and bounded control vector and an unobservable bounded pro-cess noise vector with unknown statistical characteristics, resp ., E ( ˆ x , ς P ) =: E is a known ellipsoid ( cf. § ˆ x ∈ IR n is the initial estimate of x k at k = 0, P ∈ IR n × n is a SPD matrix, ς ∈ IR ∗ + is a scaling positive scalar,the product ς P is chosen as large as the confidence in ˆ x is poor; B m ∞ is theunit ball for the ∞− norm in IR m ( cf. § A k ∈ IR n × n and B k ∈ IR n × l areknown state and input matrices, resp . and R k ∈ IR n × m is the generator ma-trix defining the shape of the zonotope bounding the unknown input vector: η k ∈ Z ( n , R k ) , where η k := R k ν k . Now consider the following measurements:¯ y k i ≤ f Tk i x k ≤ ¯ y k i , i ∈ P k := { , . . . , p k } ⊂ IN (2)The output matrix F k := [ f k . . . , f k pk ] ∈ IR n × p is time varying and so is thenumber of its columns, p k ∈ IN, which can be zero sometimes. Indeed, the mea-surements are available in varying amounts, at not all but only some sporadic,not a priori known, time steps k . Three cases can be exhaustively enumer-ated: 1) for some i ∈ D k ⊂ P k , both (finite and distinct) bounds are available:¯ y k i < ¯ y k i ; 2) for some i ∈ G k := (¯ G k ∩ ¯ G k ) ⊂ P k , only one bound, either ¯ y k i (if i ∈ ¯ G k ) or ¯ y k i (if i ∈ ¯ G k ) is available, in this case, the other (unavailable)bound is considered as ∓∞ . 3) for some other i ∈ H k ⊂ P k , the boundsare equal: ¯ y k i = ¯ y k i . The sets D k , ¯ G k , ¯ G k and H k form a partition for P k : P k = D k ∪ ¯ G k ∪ ¯ G k ∪ H k . The measurement inequalities (2) can be rewrittenas f Tk i x k ≤ ¯ y k i and ¯ y k i → −∞ , if i ∈ ¯ G k , (3a) − f Tk i x k ≤ − ¯ y k i and ¯ y k i → + ∞ , if i ∈ ¯ G k , (3b) f Tk i x k = ¯ y k i , and ¯ y k i = ¯ y k i if i ∈ H k , (3c)4 (cid:12)(cid:12) γ ki f Tk i x k − y k i (cid:12)(cid:12)(cid:12) ≤ , otherwise ( i ∈ D k ) (3d)where γ k i := ¯ y ki − ¯ y ki and y k i := ¯ y ki +¯ y ki γ ki (3e)(3a) ⇔ x k ∈ ¯ G k i := G ( f k i , ¯ y k i ) , ∀ i ∈ ¯ G k , (4a)(3b) ⇔ x k ∈ ¯ G k i := G ( − f k i , − ¯ y k i ) , ∀ i ∈ ¯ G k , (4b)(3c) ⇔ x k ∈ H k i := H ( f k i , ¯ y k i ) , ∀ i ∈ H k , (4c)(3e) ⇔ x k ∈ D k i := D (cid:16) γ ki f k i , y k i (cid:17) , ∀ i ∈ D k , (4d)where G , H and D are a halfspace, a hyperplane and a strip resp. ( cf. § Assumptions 2.1
From now on, we assume that1. all known matrices and vectors intervening in (1) and (3) , as well as theSPD P and ς ∈ IR ∗ + are bounded;2. all the columns of all the matrices intervening in (1) and those of F k arenonzero;3. the matrix F k ( H k ) := [ f k i ] i ∈ H k , intervening in (3c) or (4c) , has full col-umn rank (in order to avoid contradictory constraints leading to an emptyset); Aims 2.2
We are intending here to design an estimator ˆ x k for the state vector x k of the system (1) - (2) , such that,1. a set (ellipsoid E k of center ˆ x k ) containing all possible values of the truestate vector x k is quantified, at each time step k ∈ IN ∗ (standard require-ment for a set-membership approach);2. the state estimate vector ˆ x k is acceptable , i.e., it belongs to all the setsdefined in (3) .3. under some conditions, the estimator ˆ x k is ISS, (Input-to-State Stable, cf.Theorem A.5). This is one of the distinguishing features of the algorithmdesigned here. The other distinguishing feature is that, unlike the other set-membership tech-niques, such as those using exclusively intervals, zonotopes or polytopes, the onedetailed here delivers an optimal ( w.r.t. some chosen criterions) set, withoutany conservatism. Since the only measured information about the true statevector x k consists in its belonging to the sets defined in (3), there is no betterestimate than the one that belongs to these sets. But such an estimator is notunique and is not necessarily stable so the most suitable one will be chosenamong the set of all possible estimators by optimizing a given criterion.Let E k := E (ˆ x k , ς k P k ) be the ellipsoid containing all possible values of thetrue state vector x k . Note that the singular values of the shape matrix ς k P k correspond to the semi-lengths of its axes, whose directions are defined by theassociated–orthogonal since P k is symmetric–eigenvectors. In what follows, wehave to determine the progression law for the ellipsoid E k (and thence for thestate estimate vector ˆ x k ) such that the aims i .– iii . are fulfilled.5 Time update (prediction stage)
Let E k/k − := E (ˆ x k/k − , ς k − P k/k − ) be the ellipsoid including the “reachableset” of every possible value of x k − ∈ E k − that evolves according to theplant dynamics eq. (1a), subject to (1b). The following theorem gives theparametrized family of ellipsoids E k/k − that contains the sum of the ellipsoidresulting of the endomorphism of matrix A k − applied to E k − on one hand andthe zonotope Z ( n , R k − ), on the other. Theorem 3.1 If x k − ∈ E k − := E (ˆ x k − , ς k − P k − ) and x k obeys to (1) , then ∀ µ := ( µ , . . . , µ m ) T ∈ (cid:0) IR ∗ + (cid:1) m , x k ∈ E (ˆ x k/k − , ς k − P k/k − ) =: E k/k − := E k/k − m ⊇ E k/k − m − ⊇ . . . ⊇ E k/k − where E k/k − i := E (ˆ x k/k − , ς k − P k/k − i ) and ˆ x k/k − := A k − ˆ x k − + B k − τ k − , (5a) P k/k − := P k/k − m , (5b) P k/k − := A k − P k − A Tk − ; (5c) and, ∀ i ∈ { , . . . , m } , P k/k − i := (1 + µ i ) P k/k − i − + 1 + µ i µ i ς k − r k − i r Tk − i , (5d) r k i (the i th column of R k ) being the generator vector of the zonotope containingall possible values of the process noise η k := R k ν k . Proof . Conforming to (1a), the set containing every possible value of x k canbe schematized by (cid:0) A k − E k − ⊕ (cid:8) B k − τ k − (cid:9)(cid:1) ⊕ Z ( n , R k ) . (6)First, A k − E k − = E ( A k − ˆ x k − , ς k − A k − P k − A Tk − ) is the image of the ellip-soid E k − by the endomorphism of matrix A k − and E k/k − is its translation bythe known vector B k − τ k − . Secondly, E k/k − is the outer-bounding ellipsoidof the Minkowski sum ( cf. § .1.9.) of E k/k − and the zonotope Z ( n , R k − ): E k/k − ⊃ (cid:0) E k/k − ⊕ Z ( n , R k − ) (cid:1) . Thirdly, the zonotope Z ( n , R k ), where R k = (cid:2) r k , · · · , r k m (cid:3) , can be represented as the sum of m degenerate ellipsoids[KV14]: Z ( n , R k ) = m ⊕ i =1 E ( n , r k i r Tk i ) . Now, the Minkowski sum of two ellipsoids E ( c , P ) and E ( c , P ) is not anellipsoid, in general, yet can be bounded by a parametrized ellipsoid [MN96]: E ( c , P ( µ )) ⊃ E ( c , P ) ⊕ E ( c , P ) , ∀ µ ∈ IR ∗ + , (7)where c = c + c and P ( µ ) = (1 + µ ) P + (1 + 1 /µ ) P . (8)Applying this result sequentially to E k/k − ⊃ E k/k − ⊕ (cid:18) m ⊕ i =1 E ( n , r k − i r Tk − i ) (cid:19) , (9)eventuates in (5). ❑ The parameters µ i := µ k i are positive scalars, chosen in such a way as tominimize the size of E k/k − . The most telling measures of the size of an ellipsoidare surely the volume and the sum of the squared semi-axes lengths. Sincethe eigenvalues of ς k − P k/k − are the squared semi-axes lengths of E k/k − , theformer is proportional to their product, i.e. , to ς nk − det( P k/k − ) and the latteris equal to ς k − tr( P k/k − ). 6 heorem 3.2 Let E k +1 /k defined in Thm 3.1.1. If P k +1 /k is SPD, then E k +1 /k has the minimum volume if µ = µ vol k ,where, ∀ i ∈ { , · · · , m } , µ vol k i := n q ( n − h k i + 4 nh k i − n − n h k i , (10a) where h k i := ς − k r Tk i P − k +1 /k i − r k i . (10b) E k +1 /k has the minimum sum of the squared axes lengths, if µ = µ tr k ,where µ tr k i := r r Tki r ki ς k tr( P k +1 /ki − ) , ∀ i ∈ { , · · · , m } . (11) and the recursive formula (5b) - (5d) becomes: P k +1 /k = (cid:0) ¯ µ k ¯ µ k (cid:1)(cid:0) P k +1 /k + ¯ µ k ς k ¯ R k (cid:1) , (12a) where ¯ µ k := q ς k tr( P k +1 /k ) , where P k +1 /k given in (5c) , (12b)¯ µ k := m X i =1 k r k i k = k R k k , , (12c)¯ R k := m X i =1 k r k i k − r k i r Tk i . (12d) Proof .1. If P k/k − is SPD then P k/k − i is so, ∀ i ∈ { , . . . , m } . The volume of an el-lipsoid being proportional to the determinant of its shape matrix and ς k − being considered as constant at time step k , µ vol k i = arg min µ i det( P k +1 /k i ),(10) can be deduced from [Che94] p. 84.2. The sum of the squared semi-axes lengths of an ellipsoid being the traceof its shape matrix, µ tr k i = arg min µ i tr( P k +1 /k i ); (11) can be derived from[MN96]. As for (12), it is a direct consequence of the result [DWP01]saying that the minimum trace ellipsoid containing the Minkowski sum of m ellipsoids is : E ( c , P ) := ⊕ mi =1 E ( c i , P i ) , (13)where c := m X i =1 c i and P := (cid:16) m X i =1 p tr( P i ) (cid:17)(cid:16) m X i =1 (cid:0)p tr( P i ) (cid:1) − P i (cid:17) . (14)Then, after noticing thattr( r k − i r Tk − i ) = r Tk − i r k − i = k r k − i k (15)and that m X i =1 r k i r Tk i q tr( r k − i r Tk − i ) = R k M − k R Tk , (16)(14) applied to E k/k − ⊕ (cid:0) ⊕ mi =1 E ( n , r k − i r Tk − i (cid:1) , leads clearly to (12).It is also stated in [DWP01] that such an ellipsoid is the same that theone obtained sequentially in (11). ❑ emark 3.1 It is worth noting that the volume minimisation problem arg min µ i det( P k/k − i ) has an explicit solution here. If the unknown input vec-tor was bounded by an ellipsoid, as was the case in [MN96, DWP01, BABD08,SLZ + µ vol k i would bethe unique positive solution of an n − order polynomial to be solved at each timestep k . Remark 3.2
Because of the equality constraints introduced by the measure-ments i ∈ H k , the matrix P k looses rank during the correction stage. Hence,the E k/k − ’s volume minimization, (10) , should be avoided in favor of (11) be-cause of the inversion of P k/k − i involved in the former, in (10b) . Remark 3.3
When minimizing the sum of squared axes lengths of E k/k − , itis not necessary to compute the m intermediate values of P k/k − i , given by therecursive formula (5b) - (5d) . P k/k − can be computed directly using (12) instead.Note also that the matrix M k is invertible thanks to the assumption 2. Remark 3.4
It is possible to minimize the weighted sum of the squared axeslengths of E k/k − : tr( CP k/k − C T ) , for any C ∈ IR n C × n , n C ∈ IN ∗ . In this case,the optimal value for µ would be (cf. [Che99]) µ tr k i := s r Tk − i C T C r k − i ς k − tr( CP k/k − i − C T ) , i ∈ { , · · · , m } , (17)¯ µ k := (cid:13)(cid:13) CR k R Tk C T (cid:13)(cid:13) , = m X j k C r k i k , (18) M k := Diag( k C r k i k ) i ∈{ , ··· ,m } . (19)Given the ellipsoid at the previous time step E k − , Thm 3.1 provides theellipsoid E k/k − whose center is given by (5a) and whose shape matrix is given,up to the factor ς k − , by the recursive formula (5b)-(5d) which depends on µ ;Thm 3.2 offers the optimal values for this parameter according to two criterions,the choice of which is let to the user, in the absence of equality constraints.Otherwise, the shape matrix is calculated directly by (12). The dynamic state evolution equation (1) allowed to compute the predictedellipsoid E k/k − which contains all possible values of the state vector x k takinginto account all the measurements up to time step k − x k , obtained from the measurements:(3) ⇔ x k ∈ \ i ∈ G k G k ∩ \ i ∈ D k D k ∩ \ i ∈ H k H k j , if p k = 0 (20)It is interesting to note that the intersection of half-spaces can be consideredas a possibly unbounded polyhedron and that the intersection of strips is azonotope: \ i ∈ G k G k =: P k := P (cid:0) F k ( G k ) , ˇ y k ( G k ) (cid:1) , (21) \ i ∈ D k D k =: Z k := Z H (cid:0) F k ( D k ) , ˇ y k ( D k ) (cid:1) . (22)8he correction stage consists in performing the intersection between E k/k − andthe set (20), allowing to find E k ⊃ S k in light of the current measurements,where S k := (cid:16)(cid:0) E k/k − ∩ \ i ∈ G k G k i (cid:1) ∩ \ i ∈ D k D k i (cid:17) ∩ \ i ∈ H k H k i (23)= (cid:16)(cid:0) E k/k − ∩ P k (cid:1) ∩ Z k (cid:17) ∩ \ i ∈ H k H k i . (24)This intersection does not result in an ellipsoid in general and has to be circum-scribed by such a set. We shall begin by working on the intersection E k/k − ∩G k i in § § \ i ∈ D k D k i ; § § \ i ∈ H k H k i . Finally, all these results will be compiled in a unique state estimationalgorithm in § The intersection between the ellipsoid E k/k − obtained in § P k can be reformulated as the intersection of E k/k − and a series of strips D . Tograsp this idea, take any closed convex set S and a hyperplane H intersectingit. The intersection of S with a halfspace G delimited by H is nothing else thatits intersection with the strip formed between H and a support hyperplane of S , parallel to H and contained in G . Now, if H doesn’t intersect S , the latteris either a subset of G or lies outside of it, and if H is tangent to S (being itssupport hyperplane), then S is either again a subset of G or it has only onepoint in common with it. This is shown in the theorem below, that providesthe parameters of the intersecting strip, in the case where S is an ellipsoid. Theorem 4.1 (ellipsoid-halfspace intersec.)
Let c ∈ IR n , f ∈ IR n – { n } , P ∈ IR n × n SPSD, ς ∈ IR + and ¯ y ∈ IR .If ¯ y < − ¯ ρ , (case 1) E ( c , ςP ) ∩ G ( f , ¯ y ) = ∅ ; (25a) else if ¯ y ≥ ¯ ρ , (case 2) E ( c , ςP ) ∩ G ( f , ¯ y ) = E ( c , ςP ); (25b) else if ¯ y = − ¯ ρ , (case 3) E ( c , ςP ) ∩ G ( f , ¯ y ) = E ( c , ςP ) ∩ H ( f , − ¯ ρ ) = { c − ς ( f T P f ) − P f } ; (25c) else ( − ¯ ρ < ¯ y < ¯ ρ ), (case 4) E ( c , ςP ) ∩ G ( f , ¯ y ) = E ( c , ςP ) ∩ D ( γ f , y ) , (25d) where γ := (¯ y + ¯ ρ ) and y := γ (¯ y − ¯ ρ ) (25e)¯ ρ := ρ E ( c ,ςP ) ( − f ) = − c T f + p ς f T P f (cf. § (25f)¯ ρ := ρ E ( c ,ςP ) ( f ) = c T f + p ς f T P f . (25g)9 roof . Definition 4.1
The signed distance from a set
S ⊂ IR n to a vector x ∈ IR n is ζ ( S , x ) := max k u k =1 u T x − ρ S ( u ) . Proposition 4.2 ([KV06])
The signed distance from an ellipsoid to a hyper-plane is given by: ζ (cid:0) E ( c , P ) , H ( d , a ) (cid:1) := k d k − (cid:16)(cid:12)(cid:12) a − c T d (cid:12)(cid:12) − √ d T P d (cid:17) . (26)Let E := E ( c , ςP ). The signed distance from E to H ( f , ¯ y ) is ζ := k f k − (cid:16)(cid:12)(cid:12) ¯ y − c T f (cid:12)(cid:12) − p ς f T P f (cid:17) . (27) • ζ ≥ H ( f , ¯ y ) does not intersect E in more than one point:1. if c T f > ¯ y , then E ⊂ G ( − f , − ¯ y ) and E ∩ G ( f , ¯ y ) = ∅ ⇔ (25a);2. if c T f ≤ ¯ y , then E ⊂ G ( f , ¯ y ) ⇒ E ∩ G ( f , ¯ y ) = E ⇔ (25b);3. if c T f − ¯ y = p ς f T P f , then H ( f , ¯ y ) is tangent to E and E ∩G ( f , ¯ y ) = E ∩ H ( f , ¯ y ) = { ˇ c } , where ˇ c is calculated using (43d), letting δ ← ¯ y − c T f , with δ = ς f T P f , f ← f , P ← ς − P . • If ζ ≤
0, then H ( f , ¯ y ) intersects E and H ( − f , ρ ( − f )) is the ellipsoid’ssupporting hyperplane of normal vector − f which is contained in G ( f , ¯ y ).Indeed, x ∈ H ( − f , ρ ( − f )) ⇔ x T f − c T f = − p ς f T P f ≤ ¯ y − c T f ⇒ x T f ≤ ¯ y ⇔ x ∈ G ( f , ¯ y ) . (28)Thence, G (cid:0) − f , ρ ( − f ) (cid:1) is its supporting halfspace and E ⊂ G (cid:0) − f , ρ ( − f ) (cid:1) .Therefore, E ∩ G ( f , ¯ y ) = (cid:16) G (cid:0) − f , ρ ( − f ) (cid:1) ∩ G ( f , ¯ y ) (cid:17) ∩ E (29) (cid:12)(cid:12) ¯ y − c T f (cid:12)(cid:12) < p ς f T P f means that 0 < ¯ y + ρ ( − f ) < p ς f T P f , entailing,on one hand, G ( f , ¯ y ) = (cid:8) x (cid:12)(cid:12) x T f ≤ ¯ y (cid:9) = (cid:26) x (cid:12)(cid:12)(cid:12)(cid:12) y + ρ ( − f ) x T f ≤ y ¯ y + ρ ( − f ) (cid:27) = G ( γ − f , y + 1) . (30)and G ( − f , ρ ( − f )) = (cid:26) x (cid:12)(cid:12)(cid:12)(cid:12) − y + ρ ( − f ) x T f ≤ ρ ( − f )¯ y + ρ ( − f ) (cid:27) = G ( − γ − f , − y + 1) , (31)on the other hand. Finally, the proof (25d)–(25e) is achieved thusly: E ∩ G ( f , ¯ y ) = E ∩ (cid:16) G (cid:0) γ − f , y + 1 (cid:1) ∩ G (cid:0) − γ − f , − y + 1 (cid:1)(cid:17) (32)= E ∩ D (cid:0) γ − f , y (cid:1) . ❑ In the previous paragraph, we showed that the incorporation of the measure-ments i ∈ G k result, as for those i ∈ D k , from the intersection of the predictedellipsoid with a zonotope, formulated as an intersection of several strips. We10eed now to overbound this intersection by an ellipsoid. To begin with, thetheorem below presents a family of parametrized ellipsoids that contain an el-lipsoidal layer, coming out of the intersection of E ( c , ςP ) with the strip D ( f , y ),which can be considered–interestingly enough–as an ellipsoid unbounded in alldirections orthogonal to f . Theorem 4.3 (ellips./strip inters.)
Let c ∈ IR n , ς ∈ IR ∗ + , y ∈ IR , P ∈ IR n × n SPSD and f ∈ IR n – { n } ,if y − > ¯ ρ or y + 1 < − ¯ ρ , (case 1) D ( f , y ) ∩ E ( c , ςP ) = ∅ ; (33a) else if y + 1 ≥ ¯ ρ and y − ≤ − ¯ ρ , (case 2) D ( f , y ) ∩ E ( c , ςP ) = E ( c , ςP ); (33b) else if y = − ¯ ρ − , (case 3.a) D ( f , y ) ∩ E ( c , ςP ) = { c − ς ( f T P f ) − } ; (33c) else if y = ¯ ρ + 1 , (case 3.b) D ( f , y ) ∩ E ( c , ςP ) = { c + ς ( f T P f ) − } ; (33d) else ( − ¯ ρ < y + 1 < ¯ ρ or − ¯ ρ < y − < ¯ ρ ) , ∀ β ∈ ]0 , , (case 4) D ( f , y ) ∩ E ( c , ςP ) = D (˘ f , ˘ y ) ∩ E ( c , ςP ) ⊂ E (˘ c ( β ) , ˘ ς ( β ) ˘ P ( β )) =: E ( β ) , (34a) where ˘ f := γ f and ˘ y := γ ( f T c + δ ) , (34b)˘ P ( β ) := P − αβP f f T P, (34c)˘ c ( β ) := c + αβδP f , (34d)˘ ς ( β ) := ς + αβ (cid:0) γ (1 − β ) − − δ (cid:1) , (34e) α := (cid:0) f T P f (cid:1) − , (34f) δ := (¯ y + ¯ y − f T c ) = (¯ y + ¯ y − ¯ ρ + ¯ ρ ) , (34g) γ := (¯ y − ¯ y ) , (34h)¯ y := min( y + 1 , ¯ ρ ) and ¯ y := max( y − , ¯ ρ ) (34i) and ¯ ρ and ¯ ρ are defined in (25f) – (25g) . Proof . Let E := E ( c , ςP ). The signed distance from the ellipsoid E to each ofthe two hyperplanes H ( f , y ∓ D ( f , y ), is ζ := k f k − (cid:0)(cid:12)(cid:12) y ∓ − f T c (cid:12)(cid:12) − p ς f T P f (cid:1) . (35)When ζ >
0, the ellipsoid doesn’t intersect any of both hyperplanes meaningeither that it is situated between them i.e. , contained in the strip (case 2)or that the ellipsoid is located outside the strip, in which case (case 1), theintersection is empty. In the case 3, the interior of the ellipsoid is outside thestrip touching it in only one point and the case 3 of Thm 4.1 is then applicable: D ( f , y ) ∩ E = G ( f , y + 1) ∩ E (case 3.a) and D ( f , y ) ∩ E = G ( − f , − y + 1) ∩ E (case 3.b). In the case 4, where ζ ≤
0, the intersection is not empty. It is thenpossible to introduce the following lemma, based on the results of [FH82] and[TWS97]: 11 emma 4.4 ∀ y ∈ IR , c ∈ IR n , f ∈ IR n , σ ∈ IR ∗ + and SPD P ∈ IR n × n , if D ( f , y ) ∩ E ( c , σP ) = ∅ , then ∀ ω ∈ IR ∗ + , E (˜ c ( ω ) , ˜ σ ( ω ) ˜ P ( ω )) ⊃ D ( f , y ) ∩ E ( c , σP ) , where ˜ P ( ω ) := P − ω ( ωα + 1) − P f f T P, ˜ c ( ω ) := c + ω ( ωα + 1) − δP f = c + ω ˜ P ( ω ) f − δ, ˜ σ ( ω ) := σ + ω (1 − α ( ω + α ) − δ ) ,δ := y − f T c and α is given in (34f) . This lemma is also the mono-output case of the “observation update” part ofThm 1, [BABD08]. (34c), (34d) and (34e) are obtained by setting ω := αβ (1 − β ) − , thus β = ω ( α + ω ) − . But before applying the lemma above,it is suitable to reduce the strip D ( f , y ) in case where one of the two hyperplanesdoes not intersect the ellipsoid E , i.e. , when either y + 1 > ¯ ρ or y − < − ¯ ρ ,by translating the aforementioned hyperplane so that it becomes tangent to theellipsoid, as proposed in [BBC90]. The new strip so obtained is D ( γ − f , ˘ y ),where γ and ˘ y := y are given in (25e) and obtained by applying (case 4) of Thm4.1 to E ∩ G ( − f , − y + 1) and to E ∩ G ( f , y + 1). ❑ β Now, the optimal value of the weighting parameter β with respect to a judi-ciously chosen criterion is derived. In their well-known paper [FH82], Fogel andHuang give two optimal values of ω := αβ − β : the first minimizing the determi-nant of ˘ ς ˘ P and the second, its trace, thus optimizing the volume and the sum, resp ., of the squared semi-axes lengths of the ellipsoid E (cid:0) ωω + α (cid:1) , defined in (34).Contrary to all such algorithms in the literature, [MN96, KV97, DWP01, Che05],that minimize the size of the ellipsoid E ( β ), the optimal value of β chosen hereis the one that fulfills some stability criterion of the estimation algorithm tobe derived, in the manner of [TWS97, BABD08, SLZ + Theorem 4.5
Let E ( β ) given by (34) , where c ∈ IR n , ς ∈ IR ∗ + , y ∈ IR , P ∈ IR n × n SPSD and f ∈ IR n – { n } meet case 4 of Thm 4.3, then ˘ ς ( β ) definedin (34e) satisfies ˘ ς ( β ) = max x ∈D (˘ f , ˘ y ) ∩E ( c ,ςP ) V β ( x ) (36a) where V β ( x ) := (cid:0) x − ˘ c ( β ) (cid:1) T ˘ P ( β ) † (cid:0) x − ˘ c ( β ) (cid:1) ; (36b) and its minimum is given by β ∗ := arg min β ∈ ]0 , ˘ ς ( β ) = ( − γ | δ | − if | δ | > γ , otherwise ; (36c) where ˘ f , ˘ y , α and δ are defined in (34b) , (34f) and (34g) . Proof . Applying the generalization of the Sherman-Morrison formula to thepseudo-inverse of the matrix (34c) ( cf.
Corollary 3.5 [Xu17]), we can write˘ P ( β ) † = P † + αβ − β P † P f f T P P † . (37)12ince c ∈ R ( P ) (being the center of the ellipsoid of shape matrix ςP ) and αβP f ∈ R ( P ), by the use of (34d), it is clear that ˘ c ( β ) ∈ R ( P ). Now, noticingthat ( P P † ) T = P † P , and recalling that, for all x ∈ R ( P ), P P † x = x , thenreplacing (37) in (36b) leads to V β ( x ) := (cid:0) x − ˘ c ( β ) (cid:1) T (cid:16) P † + αβ − β P † P f f T P P † (cid:17)(cid:0) x − ˘ c ( β ) (cid:1) = (cid:0) x − ˘ c ( β ) (cid:1) T (cid:16) P † + αβ − β f f T (cid:17)(cid:0) x − ˘ c ( β ) (cid:1) (38)Inserting (34d) in (38), we can show, by the mean of some standard algebraicmanipulations, that , ∀ x ∈ R ( P ), V β ( x ) = αβγ − β (˘ y − ˘ f T x ) − αβδ + ( x − c ) T P † ( x − c ) , (39)max x ∈D (˘ f , ˘ y ) V β ( x ) = αβγ − β − αβδ + ( x − c ) T P † ( x − c ) , (40)max x ∈D (˘ f , ˘ y ) ∩E ( c ,ςP ) V β ( x ) = αβγ − β − αβδ + ς = ˘ ς ( β ) . (41)The optimal value of β is obtained by zeroing the derivative of ˘ ς :d˘ ς d β ( β ∗ ) = 0 ⇔ γ (1 − β ∗ ) − − δ = 0 ⇔ β ∗ = 1 − γ | δ | − . (42)Since β ∗ ≥
0, this solution is conditioned by | δ | > γ ; if | δ | ≤ γ , the solution tothe above minimization problem would be β ∗ = 0. ❑ Remark 4.1
The representation of the output noise vector’s bounding set asan intersection of strips, rather than as an ellipsoid, enables this optimizationproblem to have an analytical solution.
Now let us examine the intersection of an ellipsoid with a hyperplane. Thisintersection is the projection of the ellipsoid on the subspace represented bythis hyperplane and leads to a degenerate ellipsoid of lesser dimension, whoseshape matrix loses one rank with each intersecting (not parallel) hyperplane.The following theorem gives the expression of thusly obtained ellipsoid.
Theorem 4.6 (ellips./hyperplane inters.)
Let c ∈ IR n , P ∈ IR n × n SPSD, ς ∈ IR ∗ + , f ∈ IR n – { n } and y ∈ IR ,if y > ¯ ρ or y < − ¯ ρ , (case 1) E ( c , ςP ) ∩ H ( f , y ) = ∅ ; (43a) else if y = ¯ ρ = − ¯ ρ , (case 2) E ( c , ςP ) ∩ H ( f , y ) = E ( c , ςP ); (43b) otherwise ( if − ¯ ρ ≤ y ≤ ¯ ρ ) (case 3) E ( c , ςP ) ∩ H ( f , y ) = E (ˇ c , ˇ ς ˇ P ) , (43c) where ˇ c := c + αδP f , (43d)ˇ P := P − αP f f T P , (43e) V β is optimized on D (˘ f , ˘ y ) ∩E ( c , ςP ), it is then obvious that x ∈ R ( P ), since x ∈ E ( c , ςP ). ς := ς − αδ , (43f) δ := y − f T c (43g) and where α , ¯ ρ and ¯ ρ are defined in (34f) , (25f) and (25g) resp . Proof . To start with, recall that an affine map F : IR n → IR n , x L x + a turns an ellipsoid E ( c , P ) into another one E ( L c + a , L T P L ) and the hyperplane H ( f , y ) into H ( L † f , y + f T L † a ). Throughout this proof, we’ll be changing coor-dinate systems but dealing with one and the same hyperplane H := H ( f , y ) andone and the same ellipsoid E := E ( c , ςP ). Consider the vector. f ∈ IR n − { n } .Two cases (different from those of the Thm) will be distinguished dependingon whether f ∈ Ker ( P ) (1) or f / ∈ Ker ( P ) (2).1. f ∈ Ker ( P ). This means that the matrix P is not SPD but only SPSD(having at least one zero eigenvalue) and f T P f = 0. In this case E ⊂ H ′ , where H ′ := { x ∈ IR n | f T x = c } is the hyperplane of normal vector f and containingthe center c of E . If c / ∈ H , i.e. , f T c = y (corresponding to case 1 of the Thm,with f T P f = 0), E is a subset of the hyperplane H ′ parallel to H and E ∩H = ∅ ,as in (43a). Otherwise (case 2), H ′ = H , meaning that E ⊂ H and
E ∩ H = E ,as in (43b).2. Consider now f / ∈ Ker ( P ) and let ¯ n := rank( P ) ≤ n . We shall define theaffine transformation that maps the unit hypersphere or ball into the ellipsoid E ( c , ςP ): B n F −−→ E ( c , ςP ) , i.e. , F : ¯ x x = ( ςP ) (¯ x + c ) . (44)Now consider its (pseudo-)inverse transform F † that maps the ellipsoid into (apossibly degenerate) unit ball: E ( c , ςP ) F † −−→ E (¯ c , ¯ P ) , where¯ c = n and ¯ P = ¯ I n, ¯ n where ¯ I n, ¯ n := (cid:20) I ¯ n ¯ n × ¯ n ¯ n × ¯ n ¯ n × ¯ n (cid:21) , ¯ n := n − ¯ n (45)and H ( f , y ) F † −−→ H (¯ f , ¯ y ), i.e. , ¯ x ∈ H ⇔ ¯ x T ¯ f = ¯ y . In the new coordinatessystem transformed thusly, the unit normal vector to the hyperplane H and itsminimum signed distance from origin are resp .¯ f := P f p f T P f , with (cid:13)(cid:13) ¯ f (cid:13)(cid:13) = 1 and ¯ y := (cid:0) y − f T c (cid:1)p ς f T P f . (46)Let e := [1 0 . . . T the first vector of the identity matrix and H := I n − k ¯ f − e k (¯ f − e )(¯ f − e ) T (47)is the Householder symmetric ( H = H T ) and unitary ( HH T = I n ) matrix thattransforms ¯ f into e : H ¯ f = e ⇔ ¯ f = H T e = h .Next, let F : ¯ x ¯¯ x = H (¯ x − ¯ y ¯ f ) that transforms the former (second)coordinate system into the third one, in which the considered hyperplane isorthogonal to e and contains the origin: ¯¯ x ∈ H ⇔ ¯¯ x T e = 0, i.e. , H (¯ f , ¯ y ) F −−→ H (¯¯ f , ¯¯ y ) where ¯¯ f := e and ¯¯ y := 0 . (48)The (possibly degenerate) unit ball E (¯ c , ¯ P ) is transformed, by F , into the(possibly degenerate) hypersphere E (¯¯ c , ¯¯ P ), where¯¯ c := H ¯ c − ¯ yH ¯ f = − ¯ y e and ¯¯ P := H ¯ P H T = ¯ I n, ¯ n (49)Now, the distance between the center of the ellipsoid E (¯¯ c , ¯¯ P ) and the hyperplane H (¯¯ f , ¯¯ y ), (cid:12)(cid:12)(cid:12) ¯¯ c T ¯¯ f − ¯¯ y (cid:12)(cid:12)(cid:12) = (cid:12)(cid:12) − ¯ y e T e − (cid:12)(cid:12) = | ¯ y | is compared to the projection of theradius of the former onto the normal vector to the latter: q ¯¯ f T ¯¯ P ¯¯ f = q e T ¯ I n, ¯ n e = 1 . (50)14f | ¯ y | > f T P f = 0), then E ∩ H = ∅ . Otherwise (case 3), thespheroid resulting from the intersection of the (possibly degenerate) hyper-sphere E ( − ¯ y e , ¯ I n, ¯ n ) and the hyperplane H ( e ,
0) is E (¯¯ˇ c , ¯¯ˇ P ) where¯¯ˇ c := ¯¯ c + e T ¯¯ ce = − ¯ y e + ( e T e )¯ y e = n , (51a)¯¯ˇ P := (cid:0) − ( e T ¯¯ c ) (cid:1)(cid:0) ¯ I n, ¯ n − ee T (cid:1) = (cid:0) − ¯ y (cid:1)(cid:0) ¯ I n, ¯ n − ee T (cid:1) . (51b)This ellipsoid is expressed in the third coordinate system. Well, we have tofind its expression in the orignal one and for this purpose, the inverse formertransformations will be applied in reverse order: E (¯¯ˇ c , ¯¯ˇ P ) F ◦ F − −−−−−−→ E (ˇ c , ˇ σ ˇ P ). Tostart with, we’ll apply the inverse transformation F to the spheroid: E (¯¯ˇ c , ¯¯ˇ P ) F − −−−→ E (¯ˇ c , ¯ˇ P ), to obtain¯ˇ c := H T ¯¯ˇ c + ¯ y ¯ f = ¯ y ¯ f (52a)¯ˇ P := H T ¯¯ˇ P H = (cid:0) − ¯ y (cid:1)(cid:0) H T ¯ I n, ¯ n H − H T ee T H (cid:1) = (cid:0) − ¯ y (cid:1)(cid:0) ¯ I n, ¯ n − ¯ f ¯ f T (cid:1) . Then, applying F : E (¯ˇ c , ¯ˇ P ) F −−→ E (ˇ c , ˇ σ ˇ P ), yields toˇ c := ( ςP ) ¯ˇ c + c = c + ( ςP ) ¯ y ¯ f (53a)ˇ σ ˇ P = ( ςP ) T ¯ˇ P ( ςP ) = ς (cid:0) − ¯ y (cid:1)(cid:0) P − P ¯ f ¯ f T P (cid:1) . (53b)Lastly, choosing ˇ σ := ς (cid:0) − ¯ y (cid:1) and ˇ P := P − P ¯ f ¯ f T P and replacing afterwards¯ y , ¯ f , ˇ c , ˇ P and ˇ σ by their respective expressions, (46) and (53), we get to(43e) − (43f). ❑ Theorem 4.7
Let us set the following assignments ¯ y k i ← min(¯ y k i , ¯ ρ k i ) and ¯ y k i ← max(¯ y k i , − ¯ ρ k i ) . (54a) If x k satisfies (1) meeting (3) , then x k ∈ S k ⊆ E k := E (ˆ x k , ς k P k ) , (54b) where S k is defined in (24) and ∀ k ∈ IN ∗ , ς k := ς k pk , P k := P k pk and ˆ x k := ˆ x k pk ; (54c) ς k := ς k − , P k := P k/k − and ˆ x k := ˆ x k/k − ; (54d)ˆ x k/k − and P k/k − are given in (5a) and (12) ; and for i ∈ { , . . . , p k } , P k i := P k i − − α k i β k i ϕ k i ϕ Tk i , (54e)ˆ x k i := ˆ x k i − + α k i β k i δ k i ϕ k i , (54f) ς k i := ς k i − − α k i β k i δ k i ; (54g) where α k i := ( θ − k i , if p k = 0 and λ k i = 0 , , otherwise ; (54h) β k i := , if ¯ y k i = ¯ y k i and − ¯ ρ k i = ¯ ρ k i , − γ k i | δ k i | − , else if | δ k i | > γ k i , and ( − ¯ ρ k i < ¯ y k i or ¯ y k i < ¯ ρ k i ) , , otherwise; (54i) δ k i := (¯ y k i + ¯ y k i − ¯ ρ k i + ¯ ρ k i ) , (54j) A spheroid is a possibly degenerate hypersphere. k i := 12 (¯ y k i − ¯ y k i ) , (54k) θ k i := f Tk i ϕ k i , (54l) ϕ k i := P k i − f k i , (54m)¯ ρ k i := λ k i + f Tk i ˆ x k i − and ¯ ρ k i := 2 λ k i − ¯ ρ k i , (54n) λ k i := ( ς k i − θ k i ) . (54o) Proof . Direct application of Thms 4.1, 4.3, 4.5 and 4.6 to S k given in (24). ❑ The time prediction stage given by Thms 3.1 and (12) and the measurementcorrection phase, given by Thm 4.7 are concatenated to form the hole state esti-mation algorithm presented in Algorithm 1, where N is the number of samples. Remark 4.2
In the case where H k = ∅ , the matrix P k i loses rank with eachintersecting hyperplane H k i , i ∈ H k , thusly entailing the progressive flatteningof the ellipsoid E k i . Depending on the rank of the matrix R k (on which noassumption is made), the rank of P k +1 /k can be recovered at the time-updatephase. Remark 4.3
Setting either α k i = 0 or β k i = 0 results in freezing E k i − , mean-ing that the corresponding measurements f k i , ¯ y k i , ¯ y k i do not bring any usefulinformation. Remark 4.4
The cases 1 of Thms 4.1, 4.3 and 4.6 are not explicitly treatedin this theorem assuming that they can not occur since the intervening mea-surements should be consistent with the system model; yet the case where themeasurement f k i , ¯ y k i , ¯ y k i is aberrant is implicitly considered, setting again ei-ther α k i = 0 or β k i = 0 , preventing so the updating of the ellipsoid E k i − . Remark 4.5
This algorithm is of low computational complexity. Indeed, all theoperations are simple sums and products: they were optimized and are thencesuitable for high dimensional systems with many measurements. The intermedi-ate variables α k i , θ k i , λ k i , ϕ k i were added on to perform redundant vector andmatrix operations only once and noticing that f Tk i ˆ x k i − = (¯ ρ k i − ¯ ρ k i ) , allows todetermine δ k i := (¯ y k i + ¯ y k i − f Tk i ˆ x k i − ) and ¯ ρ k i := λ k i − f Tk i ˆ x k i − using additionof scalars, in (54j) and (54n) resp ., rather than multiplication of possibly highdimensional vectors. Remark 4.6
For more numerical stability and in order to avoid the explosionof the matrix P k , caused by the set summations at the prediction step, the as-signments (54d) can be replaced by P k ← ς kpk ς P k pk , ¯ ς k ← ς ς kpk ¯ ς k − and ς k ← ς .Then P k would, by itself, represent the shape of the ellipsoid E k up to a constantfactor ς − and the new variable ¯ ς k is introduced to keep track of the decreasingparameter ς k , with ¯ ς = ς . The proposed algorithm is designed in such a way as to fulfill the requirements1. - 3., expressed in the § lgorithm 1 Computation of the ellipsoid E (ˆ x k , ς k P k ) Input: ˆ x , ς , P , N Output: ˆ x k , ς k , P k n ← size of ˆ x for k = 1 , , . . . , N doInput: A k − , B k − , R k − , τ k − , intervening in (1), F k = [ f k i ] i ∈{ ,...,p k } , ¯ y k = [¯ y k i ] i ∈{ ,...,p k } , ¯ y k = [¯ y k i ] i ∈{ ,...,p k } , as in (3); Initialisation: P k/k − ← A Tk − P k − A k − ; { // Time prediction // } ¯ µ k − := q ς k − tr( P k/k − ); ¯ µ k − := m X i =1 k r k − i k as in (12b) and (12c);¯ R k − := m X i =1 k r k − i k − r k − i r Tk − i , as in (12d); P k/k − := (cid:0) ¯ µ k − ¯ µ k − (cid:1)(cid:0) P k/k − + ¯ µ k − ς k − ¯ R k − (cid:1) as stated in (12);ˆ x k/k − := A k − ˆ x k − + B k − τ k − conforming to (5a); p k ← number of columns of F k ; { // Measurement correction // } if p k = 0 then ˆ x k ← ˆ x k/k − ; P k ← P k/k − ; ς k ← ς k − ; elseInitialisation: ς k ← ς k − ; P k ← P k/k − ; ˆ x k ← ˆ x k/k − ; for i = 1 , . . . p k do ϕ k i := P k i − f k i ; θ k i := f Tk i ϕ k i as in (54m) and (54l); if θ k i = 0 then ˆ x k ← ˆ x k/k − ; P k ← P k/k − ; ς k ← ς k − ; else α k i := θ − k i ; λ k i := ( ς k i − θ k i ) as in (54h), (54o);¯ ρ k i := λ k i + f Tk i ˆ x k i − ; ¯ ρ k i := 2 λ k i − ¯ ρ k i , as in (54n);¯ y k i ← min(¯ y k i , ¯ ρ k i ); ¯ y k i ← max(¯ y k i , − ¯ ρ k i ) as in (54a); δ k i := (¯ y k i + ¯ y k i − ¯ ρ k i + ¯ ρ k i ); as in (54j); γ k i := (¯ y k i − ¯ y k i ), as in (54k); if ¯ y k i = ¯ y k i and − ¯ ρ k i = ¯ ρ k i then β k i = 1; else if | δ k i | > γ k i then β k i = 1 − γ k i | δ k i | − ; else β k i = 0; end if P k i := P k i − − α k i β k i ϕ k i ϕ Tk i , as in (54e);ˆ x k i := ˆ x k i − + α k i β k i δ k i ϕ k i , as in (54f); ς k i := ς k i − − α k i β k i δ k i as in (54g); end ifend for ˆ x k ← x k pk ; P k ← P k pk ; ς k ← ς k pk ; end ifend for -input system should be globally asymptotically stable. More formaldefinitions and results are given in the Appendix A.3. Theorem 5.1
Consider the system (1) subject to (3) and its state estimationalgorithm given by Thms 3.1 and 4.7.1. If x ∈ E ( ˆ x , σ P ) , then ∀ k ∈ IN ∗ , x k ∈ E (ˆ x k , ς k P k ) ;2. The vector ˆ x k is acceptable , i.e., it satisfies (4a) – (4d) : ∀ k ∈ K , ˆ x k ∈ \ i ∈ G k G k i ∩ \ i ∈ D k D k i ∩ \ i ∈ H k H k i , (55) where K := { k ∈ IN | p k = 0 } .3. The sequence ( ς k ) k ∈ IN ∗ is decreasing and convergent on IR + , where ς k isdefined in (54g) and it satisfies ς k = max x ∈S k V k ( x ) , where S k is given by (24) and V k ( x ) := (cid:0) x − ˆ x k (cid:1) T ˘ P † k (cid:0) x − ˆ x k (cid:1) . (56a) Proof .1. This point is satisfied by construction. Indeed, from (24), Thms 3.1 and4.7, we have x ∈ E ⇒ (cid:0) x ∈ E / (cid:1) ∧ (cid:18) x ∈ (cid:0) P ∩ Z ∩ \ i ∈ H H i (cid:1)(cid:19) ⇒ x ∈ S ⇒ x ∈ E ⇒ · · · ⇒ x k − ∈ E k − ⇒ (cid:0) x k ∈ E k/k − (cid:1) ∧ (cid:18) x k ∈ P k ∩ Z k ∩ \ i ∈ H k H k j (cid:19) ⇒ x k ∈ S k ⇒ x k ∈ E k , ∀ k ∈ IN ∗ . (57)2. This point is also granted by construction. To check it, considerˆ x k := ˆ x k/k − . From (54f), f Tk i ˆ x k i = f Tk i ˆ x k i − + α k i β k i δ k i f Tk i ϕ k i . If f Tk i ϕ k i = 0 or | δ k i | ≤ γ k i , it means that ˆ x k i is already in D k i or H k i .Else, f Tk i ˆ x k i = f Tk i ˆ x k i − + β k i δ k i . (58)Now, if i ∈ H k , β = 1 according to (54i), then inserting (54j) in (58),results in f Tk i ˆ x k i = y k i meaning that ˆ x k i ∈ H k i . Otherwise, β = 1 − γ k i | δ k i | − and f Tk i ˆ x k i − y k i = −
1, if δ k i < − γ k i and f Tk i ˆ x k i − y k i = 1,if δ k i > γ k i ; this means that ˆ x k i ∈ D k i . Combining these results for i ∈ G k ∪ D k , leads to ˆ x k ∈ S k , where S k is defined in (24) and considering(3), the proof of the point 2. is achieved.3. From (54g) of Thm 4.5, ς k i − ς k i − = − α k i β k i δ k i . Since α k i , defined in(54h), is a quadratic form when it is non-zero, it is obvious that ς k i − ς k i − ≤
0. From (54c) and (54d), ς k := ς k pk and ς k − =: ς k , then ς k − ς k − = ς k pk − ς k = − p k X i =0 α k i β k i δ k i ≤ . (59)18he sequence (cid:0) ς k (cid:1) k ∈ IN is decreasing, bounded above by ς and hence con-vergent. ❑ Theorem 5.2
Let ¯ F k := F k ( G k ∪ D k ) := [ f k i ] i ∈ G k ∩ D k ∈ IR n × ¯ p k , (60a)¯ p k := Card( G k ∩ D k ) and ¯ K = { i ∈ IN ∗ | ¯ p k = 0 } . (60b) If the pairs { A k , ¯ F Tk } and { A k , R k } are sporadically observable and completelyuniformly controllable resp ., then1. the volume of E k and all its axes lengths are bounded;2. V k , given in (56a) , is an ISS-Lyapunov function for the estimation error ˜ x k := x k − ˆ x k , which is ISS . Proof . The proof is detailed in the Appendix A. ❑ The matrices of the system model (1) and (3) are generated randomly for twovalues of the state dimension: n = 10 and n = 100, with m = n , ̺ = q = r = s = n , π = n and µ = µ tr k ; the input vector B k − τ k − contains sine entriesof random magnitude and frequency. P = 100 I n , ς = 1 and ˆ x k randomlychosen on the boundary of E . The measurements are available at all time steps K = { , . . . , N } in the case 1; at some randomly chosen time steps, in case 2and K = ∅ in the case 3 (where only prediction stage is performed without anymeasurement correction). For each case, the simulations are run 25 times underMATLAB R2018b on Intel Core i7 (2.3GHz, 8G RAM), each one for a differentsystem model and containing N = 100 time steps. The results are summarizedin Table 1. Let ς tr( P ) := mean ς k tr( P k ) k ∈{ ,...,N } : the average sum of E k ’s squared axeslengths, k ˜ x k := mean k ˜ x k k k ∈{ ,...,N } : the mean estimation error vector norm and T : theaverage computational time for the simulation horizon of N time steps. It is n K ς tr( P N ) ς N tr( P ) ς tr( P ) k ˜ x N kk ˜ x k k ˜ x k T ( ms )10 case 1 0.010 1 141 0.027 3.02 18case 2 0.026 2 206 0.053 4.60 15case 3 0.058 7 540 0.066 7.25 10100 case 1 0.712 7 · · · cf. Definition A.3 in the Appendix A.2 cf. Definitions A.4, A.5 and the Lemma A.5. CONCLUSION
We have proposed an ellipsoidal state characterization for discrete-time lin-ear dynamic models with linear-in-state sporadic measurements, which are cor-rupted by additive unknown process and measurement disturbances, enclosedby zonotopes, on one hand and subject to linear equality and inequality con-straints on the other hand. Here is a turnkey, ready to use, easily implementablealgorithm, without any parameter to tune.A particular attention was accorded first to the stability of the estimationalgorithm, which is ISS despite of the irregularity of the measurements; then toits computational efficiency. Indeed, the proposed algorithm is composed of onlylow demanding, optimized in this sense operations (matrix sums and products),no costly tools nor heavy operations such as interval arithmetic or LMI, not evenmatrix inversion have to be performed, what makes this algorithm suitable forhigh dimensional systems. Furthermore, the challenge faced in other Kalman-like algorithms, inherent to the inversion of the state error covariance matrix–which is inevitably singular in presence of equality constraints–is circumventedhere since the matrix P k is actually never inverted. A Appendix: Stability analysis
A.1 Kalman filter analogy
To prove Thm 5.2, we’ll be using the observability and controllability propertiesof the Kalman filter. For this purpose, we have to show the analogy of thelatter with the proposed algorithm. Consider the following linear time-varyingstochastic system with some bounded matrix ¯ A k ∈ IR n × n : ξ k = ¯ A k − ξ k − + B k − τ k − + R k − w k − , k ∈ IN ∗ (61a)ˇ y k = ¯ F Tk ξ k + v k , ∀ k ∈ ¯ K , ( cf. (60b)) (61b)where w k ∼ N ( n , W k ) (61c) v k ∼ N ( ¯ p k , V k ) (61d)where ξ k ∈ IR n is the unknown state vector, ˇ y k := [ y k i ] i ∈ G k ∩ D k ∈ IR ¯ p k , y k := (¯ y k − ¯ y k ) and ¯ F k ∈ IR n × ¯ p k , defined in (60), are the output vectorand the observation matrix resp .; B k − τ k − is the known input intervening in(1a); and w k − and v k are gaussian centred noise vectors of covariance matrices W k − and V k resp. Now consider the Kalman filter, designed for the system (61):ˆ ξ k = ˆ ξ k/k − + K k δ k (62a)¯ P k = ( I n − K k ¯ F Tk ) ¯ P k/k − (62b) δ k := ˇ y k − ¯ F Tk ˆ ξ k/k − (62c) K k := ( ¯ P k/k − ¯ F k ( ¯ F Tk ¯ P k − ¯ F k + V k ) − , if k ∈ ¯ K n × ¯ p k , otherwise; (62d)ˆ ξ k/k − = ¯ A k − ˆ ξ k − + B k − τ k − (62e)¯ P k/k − = ¯ A k − ¯ P k − ¯ A Tk − + R k − W k − R Tk − . (62f)20he time prediction stage, (5), of Thm 3.1 can be seen as the prediction stage ofthe Kalman filter (62e)-(62f) and the measurement correction stage (54), givenin Thm 4.7 is nothing else than (62a)-(62d). This is stated in Proposition A.1.Forasmuch as the Kalman filter undergoes numerical stability issues when thesystem (61) is subject to equality constraints (the matrix ¯ F Tk P k − ¯ F k + V k in theKalman gain, (62d), becoming ill-conditioned), (3c) are not considered for themoment. Proposition A.1 If ˆ x k is computed in line with the algorithm composed of (5a) , (12) and (54) and if ˆ ξ k is the Kalman estimator (62) designed for thesystem (61) , such that ¯ A k := λ k A k , λ k := q ¯ µ k ¯ µ k , (63a) W k := ¯ µ k +¯ µ k ς k Diag(¯ µ − k i ) i ∈{ , ··· ,m } , (63b) V k := Diag (cid:0) ω − k i (cid:1) i ∈ G k ∩ D k , where ω k i := α ki β ki − β ki (63c) and where ¯ µ k , ¯ µ k , α k i , β k i are defined in (12b) , (12c) , (54h) , (54i) resp . ς k := ς k − − δ Tk Υ k δ k , where Υ k := Diag (cid:0) α k i β k i (cid:1) i ∈ G k ∩ D k , for a fixed ς k ; (64) and if ˆ x = ˆ ξ and ¯ P = P , then ∀ k ∈ IN ∗ , ˆ x k = ˆ ξ k and ¯ P k = P k . (65) Proof . Replacing β k i = ( ω k i + α k i ) − ω k i and α k i from (54h) in (54e), the lattercan be rewritten P k i = P k i − − P k i − f k i (cid:0) f Tk i P k i − f k i + ω − k i (cid:1) − f Tk i P k i − . Then, using the inversion lemma, it comes that P − k = P − k ¯ pk = P − k + ¯ p k X i =1 ω k i f k i f Tk i . (66)Recalling that P k = P k ¯ pk and that P k = P k/k − and noticing that ¯ p k X i =1 ω k i f k i f Tk i = ¯ F Tk V k ¯ F k (67)we have P − k = P − k/k − + ¯ F Tk V − k ¯ F k . (68)Applying the inversion lemma again to (68), the algorithm (54) can be rewrittenas (61b), (62) and (63). Finally, using P k +1 /k defined in (12) and considering(63a) and (63b), we obtain (62e)-(62f), which completes the proof. ❑ A.2 Observability and controllability
Before examining the observability and the controllability of the studied sys-tem (1)-(2), we need to define the controllability and observability gramians, oflength l ∈ IN: ¯ C k + l,k := k + l − X i = k ¯ λ − i +1 ,k Φ k,i +1 R i W i R Ti Φ Tk,i +1 (69a)¯ O k + l,k := k + l X i = k ¯ λ i,k Φ Ti,k ¯ F i V − i ¯ F Ti Φ i,k (69b)21here Φ k + l,k := A k + l − . . . A k , with Φ k,k + l = Φ − k + l,k , is the state transition ma-trix associated to A k ∈ IR n × n which is assumed to be invertible ;¯ λ k + l,k := λ k + l − . . . λ k , where λ k is defined in (63a); V k ∈ IR ¯ p k × ¯ p k is an SPDmatrix, R k ∈ IR n × m and ¯ F k ∈ IR n × ¯ p k . Definition A.1 (uniform complete controllability)
The matrix pair { ¯ A k , R k W k } is uniformly completely controllable , if there exist positive con-stants ¯ ̺ and ¯ ̺ and a positive integer h , such that, for all k ≥ h , ¯ ̺ I n ≤ ¯ C k,k − h ≤ ¯ ̺ I n . (70) Definition A.2 (uniform complete observability)
The matrix pair { ¯ A k , ¯ F Tk } is uniformly completely observable , if there exist positive constants ¯ ̺ and ¯ ̺ and a positive integer h , such that, for all k ≥ h , ¯ ̺ I n ≤ ¯ O k,k − h ≤ ¯ ̺ I n . (71) Proposition A.2 ([SG95, BGEGG09])
Let ¯ K = IN ∗ (cf. (60b) ). If the ma-trix pairs { ¯ A k , R k W k } and { ¯ A k , ¯ F Tk } are uniformly completely controllable andobservable resp ., the estimation covariance matrix of the Kalman filter (62) ,designed for the system (61) , satisfies the following inequalities, for all k ≥ l : ( ¯ O k,k − l + ¯ C − k,k − l ) − ≤ ¯ P k ≤ ¯ O − k,k − l + ¯ C k,k − l . Proposition A.3
The pairs { ¯ A k , R k W k } and { ¯ A k , ¯ F Tk } are uniformly com-pletely controllable and observable resp ., if and only if { A k , R k } and { A k , ¯ F Tk } have the respective properties. Proof . Since λ k and W k , given in (63a) and (63b), are a both bounded andpositive ( resp . SPD), the observability and controllability gramians, associ-ated to the matrices A k , R k and ¯ F k : O k,k − l := P k + li = k Φ Ti,k ¯ F i V − i ¯ F Ti Φ i,k and C k,k − l := P k + l − i = k Φ k,i +1 R i R Ti Φ Tk,i +1 , are SPD bounded matrices if and only if ¯ O k,k − l and ¯ C k,k − l , given by (69), associated to λ k A k , R k W k and ¯ F k are alsobounded SPD matrices. ❑ It is needless to say that it is difficult to ensure the full rank for the matrixsum ¯ O k,k − h (69b) on a time window of constant length h when dealing withsporadic measurements. The system (1)-(3) can therefore not be uniformlycompletely observable. Let us then introduce a new observability criterion forthis kind of systems by using the observabilty gramian with a variable length: Definition A.3 (sporadic observability)
Assuming that ¯ F k := 0 n × ¯ p k , ∀ k / ∈ ¯ K (cf. (60b) ), the pair { ¯ A k , ¯ F Tk } is said sporadically observable , if there existpositive constants ¯ ̺ and ¯ ̺ and a positive integer h , such that, for all k ≥ κ k ( h ) , ¯ ̺ I n ≤ ¯ O k,k − κ k ( h ) ≤ ¯ ̺ I n , where κ k ( h ) ∈ IN is s.t.Card (cid:0) { i ∈ ¯ K| k − κ k ( h ) ≤ i ≤ k } (cid:1) = h (72) and where Card( S ) stands for the cardinality (number of elements) of the set S . The direct consequence of Proposition A.3 applied to the system with all i ∈ G k ∩ D k ∩ H k measurements (3) including equality constraints (3c) can be statedas follows: 22 orollary A.4 Consider the system (1) with (3) and the matrix P k computedin line with (5) , (12) and (54) . Let H k ∈ IR n × ¯ n k whose columns form or-thonormal basis for R (cid:0) P k (cid:1) where ¯ n k := rank( P k ) and let ¯ P k := H k P k H Tk . Ifthe pair { A k , ¯ F Tk } is sporadically observable and { A k , R k } is completely uni-formly controllable, then there exist positive finite numbers ̺ k and ̺ k , s.t. , forall k ≥ κ k ( l ) , ̺ k I ¯ n k ≤ ¯ P k ≤ ̺ k I ¯ n k , (73) κ k being defined in (72) . A.3 Input-to-State stability
In this paragraph, we shall examine the ISS concept. Before doing so, let usrecall some comparison functions, widely used in stability analysis. A continuousfunction ψ : IR + → IR + is called positive definite if it satisfies ψ (0) = 0 and ψ ( t ) > ∀ t >
0. A positive definite function is of class K if it is strictlyincreasing and of class K ∞ if it is of class K and unbounded. A continuousfunction ψ : IR + → IR + is of class L if ψ ( t ) is strictly decreasing to 0 as t → ∞ and a continuous function ψ : IR + × IR + → IR + is of class K L if it isof class K in the first argument and of class L in the second argument. Definition A.4 ([JW01])
The system z ( k + 1) = f ( z ( k ) , u ( k )) , (74a) where f ( , ) = and z (0) := z , (74b) is globally input-to-state stable (ISS), if there exists a K L -function β and a K -function ψ such that, for each bounded input sequence u [0 ,k ] := { u , . . . , u k } and each z ∈ IR n , (cid:13)(cid:13) z ( k, z , u [0 ,k − ) (cid:13)(cid:13) ≤ β ( k x k , k ) + ψ ( sup i ≤ k − k u i k ) , where z ( k, z , u [0 ,k − ) is the trajectory of the system (74) , for the initial state z ∈ IR n and the input sequence u [0 ,k − . Definition A.5 ([JW01])
A continuous function V : IR n → IR + is an ISS-Lyapunov function for the system (74), if there exists K ∞ -functions ψ and ψ such that for all z ∈ IR n , ψ ( k z k ) ≤ V ( z ) ≤ ψ ( k z k ) and there exists a K ∞ -function ψ and a K -function χ such that for all z ∈ IR n and all u ∈ IR m , V (cid:0) f ( z , u ) (cid:1) − V ( z ) ≤ − ψ ( k z k ) + χ ( k u k ). Lemma A.5 ([JW01])
The system (74) is ISS if it admits a continuous ISS-Lyapunov function.
A.4 Stability of the proposed estimation algorithm
This section is dedicated to the proof of Theorem 5.2.1. Let ¯ n k := n − ¯ n k = n − rank( n ) and H k ∈ IR n × ¯ n k whose columns formorthonormal basis for Ker (cid:0) P k (cid:1) , s.t. , H k := [ H k | H k ] is a unitary ma-trix. Hence, we have H Tk P k H k = Bdiag( ¯ P k , ¯ n k × ¯ n k ). The E k ’s semi-axeslengths are the singular values of the matrix ς k P k , which are those of ς k H Tk P k H k = ς k ¯ P k . On one hand, it is shown, at the point 3. of Thm5.1, that the sequence (cid:0) ς k (cid:1) k ∈ IN is decreasing, bounded above by ς and23onvergent. On the other hand, as stated in (73) of Corollary A.4, thesingular values of ¯ P k are bounded and so are the ellipsoid’s axes lengths,as well as their product representing the ellipsoid’s volume.2. First, for any possible value of the true state vector x k , we have x k ∈ E (ˆ x k , ς k P k ) ⇔ ˜ x k := x k − ˆ x k ∈ E ( n , ς k P k ) ⇔ ˜ x k = ( ς k i P k ) u k , u k ∈ B n , ( cf. §
1. 8.). (75)It means that ˜ x k ∈ R (cid:0) P k (cid:1) , which is a subspace of IR n of dimension ¯ n k ≤ n : H Tk ˜ x k = ς k H Tk P k H k H Tk u k = ς k h ¯ P k ¯ n k × ¯ n k ¯ n k × ¯ n k ¯ n k × ¯ n k i H Tk u k = h ς k ¯ P k ¯ u k ¯ n k i , where ¯ u k := H Tk u k ∈ B ¯ n , meaning that ∀ x k ∈ E k , ˜ x k = H Tk [˜ x Tk T ¯ n k ] T , where ˜ x k := H Tk ˜ x k and H Tk ˜ x k = ¯ n k . Now we shall show that V k is an ISS-Lyapunov function for all possiblevalues of ˜ x k ∈ R (cid:0) P k (cid:1) . First, V k := ˜ x Tk H k H Tk P † k H k H Tk ˜ x k = (cid:2) ˜ x Tk | ˜ x Tk (cid:3)h ¯ P − k i(cid:2) ˜ x Tk | ˜ x Tk (cid:3) T = ˜ x Tk ¯ P − k ˜ x k noticing that k ˜ x k k = (cid:13)(cid:13) H Tk ˜ x k (cid:13)(cid:13) = k ˜ x k k and by virtue of (73), it can bededuced that ψ k ( k ˜ x k k ) ≤ V k ≤ ψ k ( k ˜ x k k ) , (76)where ψ k i : IR + → IR + , i ∈ { , } ,t ψ k i ( t ) = ̺ − k i t , are K ∞ functions. (77)Second, from the point 3. of Thm 5.1, we have V k − V k/k − ≤ − p k X i =0 α k i β k i δ k i ≤ V k/k − := ˜ x Tk/k − P † k/k − ˜ x k/k − (79)and where ˜ x k/k − := A k − ˜ x k − + η k − (80)Basing on the same reasoning as done in Lemma 3 in [SLZ + a , b ∈ IR n and any matrices A, B ∈ IR n × n ,( a + b ) T ( A + B ) † ( a + b ) ≤ a T A † a + b T B † b (81)and considering P k +1 /k given by (12), we have V k/k − ≤ ¯ µ k − ¯ µ k − +¯ µ k − (cid:16) ˜ x Tk − A Tk − (cid:0) A k − P k − A Tk − (cid:1) † A k − ˜ x k − + ς k − ¯ µ k − η Tk − ¯ R † k − η k − (cid:17) ≤ ¯ µ k − ¯ µ k − +¯ µ k − ˜ x Tk − P † k − ˜ x k − + ς k − ¯ µ k − +¯ µ k − (cid:13)(cid:13)(cid:13) ¯ R † k − (cid:13)(cid:13)(cid:13) k η k − k . (82)Now, from (78), we have V k − V k − ≤ V k/k − − V k − ; (83)and consequently, (82) becomes V k − V k − ≤ V k/k − − V k − ≤ − ̺ k V k − + ψ k ( k η k − k ) , (84)where ̺ k := ¯ µ k ¯ µ k +¯ µ k > φ k : IR + → IR + , t φ k ( t ) = ς k σ k ¯ µ k +¯ µ k t , is a K ∞ function, where σ k = (cid:13)(cid:13)(cid:13)(cid:0) R k M − k R Tk (cid:1) † (cid:13)(cid:13)(cid:13) >
0. This means that V k is anISS-Lyapunov function for the system of state vector ˜ x k . Thus applyingLemma A.5 completes the proof of the theorem.The cases where k / ∈ ¯ K can be viewed as measurements i for which α k i = 0or β k i = 0. 24 emark A.1 H k := [ H k | H k ] ∈ IR n × n is a unitary matrix which rotates P k into a basis where it has two-bloc-diagonal form ( H k can be obtained by QRdecomposition of P k or by SVD of P k ). H Tk ˜ x k = [˜ x Tk | ˜ x Tk ] T gives the compo-nents of the state estimation error vector, ˜ x k , in this new rotated basis, wherethe ¯ n k first components are ISS (point 2. of Thm 5.1) and the last ¯ n k ones arezero, meaning that the corresponding estimations are equal to their true values,in this rotated basis. References [AIBS19] Leif Erik Andersson, Lars Imsland, Edmund F. Brekke, andFrancesco Scibilia. On kalman filtering with linear state equal-ity constraints.
Automatica , 101:467 – 470, 2019.[BABD08] Y. Becis-Aubry, M. Boutayeb, and M. Darouach. State estimationin the presence of bounded disturbances.
Automatica , 44:1867–1873, 2008.[BBC90] G. Belforte, B. Bona, and V. Cerone. Parameter Estimation Al-gorithm for a Set-Membership Description of Uncertainty.
Auto-matica , 26(5):887–898, September 1990.[BGEGG09] Vibhor L. Bageshwar, Demoz Gebre-Egziabher, William L. Gar-rard, and Tryphon T. Georgiou. Stochastic observability test fordiscrete-time kalman filters.
Journal of Guidance Control and Dy-namics , 32(4):1356–1370, 2009.[Che94] F. L Chernousko.
State estimation for dynamic systems . BocaRaton: CRC Press, 1994. Includes bibliographical references (p.293-299) and index.[Che99] F. L. Chernousko. What is ellipsoidal modelling and how to use itfor control and state estimation? In Isaac Elishakoff, editor,
Whysand Hows in Uncertainty Modelling , pages 127–188, Vienna, 1999.Springer Vienna.[Che05] F. L. Chernousko. Ellipsoidal state estimation for dynamical sys-tems.
Nonlinear Analysis , 63:872–879, 2005.[Com15] Christophe Combastel. Zonotopes and kalman observers: Gainoptimality under distinct uncertainty paradigms and robust con-vergence.
Automatica , 55:265 – 273, 2015.[DL13] Z. Duan and X. R. Li. The role of pseudo measurements inequality-constrained state estimation.
IEEE Transactions onAerospace and Electronic Systems , 49(3):1654–1666, July 2013.[DWP01] C. Durieu, E. Walter, and B. Polyak. Multi-input multi-outputellipsoidal state bounding.
Journal of Optimization Theory andApplications , 111(2):273–303, 2001.25FH82] E. Fogel and Y. F. Huang. On the value of information in systemidentification - bounded noise case.
Automatica , 18(2):229–238,1982.[JKDW12] Luc Jaulin, Michel Kieffer, Olivier Didrit, and Eric Walter.
Ap-plied Interval Analysis: with Examples in Parameter and StateEstimation, Robust Control and Robotics . Springer London Ltd,2012.[JW01] Zhong-Ping Jiang and Yuan Wang. Input-to-state stability fordiscrete-time nonlinear systems.
Automatica , 37:857–869, 2001.[JZ13] Chaoyang Jiang and Yong-An Zhang. Some results on linear equal-ity constrained state filtering.
International Journal of Control , 86,12 2013.[KV97] A. Kurzhanskiy and I. V´alyi.
Ellipsoidal Calculus for Estimationand Control . Systems & Control: Foundations & Applications.Birkhauser, Boston, Basel, Berlin, 1997.[KV06] A. A. Kurzhanskiy and P. Varaiya. Ellipsoidal toolbox. Techni-cal Report UCB/EECS-2006-46, EECS Department, University ofCalifornia, Berkeley, May 2006.[KV14] Alexander A. Kurzhanskiy and Pravin Varaiya.
Dynamics andControl of Trajectory Tubes . Birkh¨auser, 1 edition edition, October27 2014.[MN96] D. Maksarov and J. P. Norton. State bounding with ellipsoidal setdescription of the uncertainty.
International Journal of Control ,65(5):847–866, 1996.[NBH15] Benjamin Noack, Marcus Baum, and Uwe Hanebeck. State es-timation for ellipsoidally constrained dynamic systems with set-membership pseudo measurements. In , pages 297–302, 09 2015.[RJ15] Nacim Ramdani and Luc Jaulin.
Interval Methods and Applica-tions , volume 8. Mathematics in Computer Science, June 2015.[SG95] Y. Song and J. W. Grizzle. The Extended Kalman Filter as aLocal Asymptotic Observer for Discrete-time Nonlinear Systems.
Journal of Mathematical Systems Estimation and Control , 5(1):59–78, 1995.[Sim10] Dan Simon. Kalman filtering with state constraints: A survey oflinear and nonlinear algorithms.
Control Theory and Applications,IET , 4:1303 – 1318, 09 2010.[SLZ +
18] Qiang Shen, Jieyu Liu, Xiaogang Zhou, Qian Zhao, and WangQi. Low-complexity iss state estimation approach with boundeddisturbances.
International Journal of Adaptive Control and SignalProcessing , 32:1473–1488, July 2018.26TWS97] Guojie Tan, Changyun Wen, and Yeng Chai Soh. Identificationfor systems with bounded noise.
IEEE Transactions on AutomaticControl , 42(7):996–1001, 1997.[Xu17] Xuefeng Xu. Generalization of the sherman–morrison–woodburyformula involving the schur complement.
Applied Mathematics andComputation , 309:183 – 191, 2017.[YL09a] F. Yang and Y. Li. Set-membership filtering with state con-straints.
IEEE Transactions on Aerospace and Electronic Systems ,45(4):1619–1629, Oct 2009.[YL09b] Fuwen Yang and Yongmin Li. Set-membership filtering fordiscrete-time systems with nonlinear equality constraints.