Guaranteed methods based on constrained zonotopes for set-valued state estimation of nonlinear discrete-time systems

This paper presents new methods for set-valued state estimation of nonlinear discrete-time systems with unknown-but-bounded uncertainties. A single time step involves propagating an enclosure of the system states through the nonlinear dynamics (predi…

Authors: Brenner S. Rego, Guilherme V. Raffo, Joseph K. Scott

Guaranteed methods based on constrained zonotopes for set-valued state   estimation of nonlinear discrete-time systems
Guaranteed methods based on constrained zonotopes for set-v alued state estimation of nonlinear discrete-time systems I Brenner S. Rego a, ∗ , Guilherme V . Ra ff o a,b , Joseph K. Scott c , Davide M. Raimondo d a Graduate Pr ogr am in Electrical Engineering, F ederal University of Minas Ger ais, Belo Horizonte, MG 31270-901, Brazil b Department of Electr onics Engineering, F ederal University of Minas Ger ais, Belo Horizonte, MG 31270-901, Brazil c Department of Chemical and Biomolecular Engineering, Clemson University , Clemson, SC, USA d Department of Electrical, Computer and Biomedical Engineering, University of P avia, Italy Abstract This paper presents new methods for set-valued state estimation of nonlinear discrete-time systems with unkno wn-but-bounded uncertainties. A single time step in volves propagating an enclosure of the system states through the nonlinear dynamics (prediction), and then enclosing the intersection of this set with a bounded-error measurement (update). When these enclosures are represented by simple sets such as interv als, ellipsoids, parallelotopes, and zonotopes, certain set operations can be v ery conservati ve. Y et, using general con vex polytopes is much more computationally demanding. T o address this, this paper presents two ne w methods, a mean value extension and a first-order T aylor e xtension, for e ffi ciently propagating constrained zonotopes through nonlinear mappings. These e xtend existing methods for zonotopes in a consistent way . Examples sho w that these e xtensions yield tighter prediction enclosures than zonotopic estimation methods, while largely retaining the computational benefits of zonotopes. Moreo ver , they enable tighter update enclosures because constrained zonotopes can represent intersections much more accurately than zonotopes. K e ywor ds: Nonlinear state estimation, Set-based computing, Reachability analysis, Con vex polytopes 1. Introduction The importance of state estimation has become evident in many fields of research ov er the years. Examples of recurrent applications are the localization problem [ 1 , 2 ], state-feedback control [ 3 , 4 , 5 ], and fault detection and isolation (FDI) [ 6 , 7 , 8 ]. If stochastic descriptions of the process and measurement uncer- tainties are available, then Bayesian state estimation methods such as Kalman filtering or particle filtering are typically em- ployed. On the other hand, if only bounds on the uncertainties are known, then set-v alued state estimation methods are applied [ 3 , 9 , 10 , 11 , 12 ]. Recent approaches also consider the presence of both kinds of uncertainties in the system [ 13 ]. Set-valued state estimation methods aim to construct com- pact sets that are guaranteed to enclose all possible trajectories of the system subject to unknown-but-bounded uncertainties [ 9 , 14 ]. Using the standard recursive approach for discrete-time systems, this in volv es first bounding the image of the current enclosure under the dynamics (prediction), and then enclosing the intersection of this set with the set of states consistent with a I c  2019. This manuscript version is made av ailable under the CC-BY -NC- ND 4.0 license https://creativecommons.org/licenses/by- nc- nd/4. 0/ . This work was supported in part by the project INCT under the grant CNPq 465755 / 2014-3, F APESP 2014 / 50851-0, and also by the Brazilian agencies CAPES, and F APEMIG. The material in this paper was not presented at any conference. ∗ Corresponding author Email addr esses: brennersr7@ufmg.br (Brenner S. Re go), raffo@ufmg.br (Guilherme V . Ra ff o), jks9@clemson.edu (Joseph K. Scott), davide.raimondo@unipv.it (Da vide M. Raimondo) bounded-error measurement (update). For discrete-time linear systems, if the initial set is a polytope, then e xact enclosures can in principle be computed using conv ex polytopes [ 15 ]. Ho w- ev er, even for linear dynamics, polytope propagation requires demanding operations (e.g., polytope projection, Minkowski sum, or con version between vertex and halfspace representa- tions) whose complexity grows dramatically with time due to the complexity increase of the resulting sets [ 16 , 17 ]. For these reasons, enclosures are often described by simpler sets including ellipsoids [ 18 , 19 , 20 ], parallelotopes [ 9 , 21 ], zonotopes [ 10 , 22 ], or combinations of these [ 23 ]. Howe ver , the mathematical limi- tations of these sets require certain operations to be conserv ativ e, sometimes quite significantly . Notably , this includes set inter- section, which is critical for the update step in set-valued state estimation [ 9 , 10 , 19 ]. The article [ 24 ] proposes the use of zonotope bundles to describe intersections of zonotopes without explicit computation. Ho wev er , the Minko wski sum and linear image (see Section 3 ) are outer -approximated. In the article [ 25 ], constrained zonotopes are introduced to o vercome man y of the limitations of zonotopes. These sets are closed under intersection, Minko wski sum, and linear image, and are capable of describing arbitrary conv ex polytopes if the complexity of the set description is not limited. E ffi cient algorithms for lin- ear set-v alued state estimation and also FDI using constrained zonotopes are described in [ 25 , 8 ]. In contrast to the linear case, e ff ective set-v alued state estimation for nonlinear systems is still an open challenge [ 3 , 14 , 26 , 27 ]. Early approaches in this field used inclusion functions based on interval arithmetic [ 28 ] to propagate bounds Pr eprint submitted to Automatica September 2, 2019 through the nonlinear dynamics, and used interv al-based set in- version techniques to enclose the set of states consistent with the current measurement [ 1 , 3 , 29 ]. Improved accurac y is achiev ed using refinements (i.e., unions of intervals) as well as more advanced interv al methods such as contractor and separator al- gebras [ 1 , 30 ]. Unfortunately , even these methods often pro vide conserv ativ e bounds without e xtensive refinement, which is only tractable for systems with relatively few states [ 3 ]. An inter- esting new interv al method based on discrete-time di ff erential inequalities is proposed in [ 31 ], b ut it only applies to Euler discretized systems with limited step size. A few alternatives for nonlinear set-valued state estimation can be found in the literature. Polytopes are used in [ 32 ] for systems with nonlinear dynamics and linear measurements. The prediction step is performed based on a linearization of the dy- namics, in which conserv ativ e interval enclosures are used to bound the linearization error . Notably , the computed error bound is valid for any v alue of the state, rather than being computed as a function of the current state enclosure, which can lead to a very conserv ative result. Another issue with this method is that the complexity of the polytopic enclosures gro ws rapidly with time, which results in a v ery high computational burden because the complexity of the required set operations scales poorly with increasing polytope complexity . More e ffi cient methods based on zonotopes are proposed in [ 14 ] and [ 26 ]. The propagation step in [ 14 ] is based on the Mean V alue Theorem and is referred to as the mean value extension , while the approach in [ 26 ] uses a first-order T aylor expansion with a rigorous remainder bound, and is referred to as the first-or der T aylor e xtension . Updates are then achiev ed by methods for outer-approximating the intersec- tion of a zonotope with a strip (i.e., a linear measurement with bounded error). An alternativ e zonotope-based prediction step using DC programming is proposed in [ 33 ], but with the same update as in [ 14 ]. Even for linear measurements, the symmetry of zonotopes is kno wn to cause significant errors in the update step [ 25 ]. General conv ex polytopes in halfspace representation are used in [ 27 ] to enable an exact update. Prediction is then done by representing the polytope as an intersection of zono- topes and applying the mean value extension. Unfortunately , con version between these representations is computationally demanding, and the increasing complexity of the zonotope b un- dle with time is not addressed. Constrained zonotopes hav e recently been applied to nonlinear state estimation in [ 34 ] and shown to pro vide much higher accuracy than existing zonotopic methods. Nevertheless, the method in [ 34 ] uses an interval parti- tioning scheme and is therefore intractable for high-dimensional systems. Finally , in an e ff ort to ov ercome the limitations of con vex sets, polynomial zonotopes are introduced in [ 35 ] and used for reachability analysis. Ho wev er , update algorithms for polynomial zonotopes hav e not been dev eloped. In this context, the main contributions of this paper are two new methods for nonlinear set-valued state estimation based on constrained zonotopes. W e follow the standard algorithmic steps typically used for set-valued state estimation (i.e., pre- diction, update, and reduction). For the prediction step, we use new generalizations of the mean value extension and first- order T aylor extension discussed abov e that enable constrained zonotopes, rather than zonotopes, to be e ff ecti vely propagated through nonlinear discrete-time dynamics. Since this class of sets corresponds to an alternativ e representation of con ve x poly- topes, these generalizations can be vie wed as nov el approaches for implicitly propagating con vex polytopes through nonlinear mappings. The generalization of these methods to constrained zonotopes is not straightforward and requires significant mod- ifications to the existing proofs. Ho wev er , it results in much tighter prediction enclosures than e xisting zonotopic methods in some cases, as shown in the numerical examples. More- ov er , using constrained zonotopes for prediction also enables the update step to be done much more e ff ectiv ely using the generalized intersection operation for constrained zonotopes. On the other hand, reduction becomes more complex than for zonotopes, leading to interesting trade-o ff s between accuracy and computational e ffi cienc y . W e inv estigate these trade-o ff s both through numerical examples and by de veloping a detailed computational complexity analysis of the proposed methods and their zonotopic counterparts. T o the best of our knowledge, such an analysis has not been conducted pre viously for either zonotopes or constrained zonotopes in the context of nonlinear set-valued estimation. The remainder of the paper is organized as follows. The nonlinear set-valued state estimation problem is stated in Section 2 . Essential mathematical background is presented in Section 3 , including a discussion of constrained zonotopes and their main properties. Section 4 develops the main results of the paper; namely , the proposed mean v alue and first-order T aylor extensions for constrained zonotopes, heuristics for selecting the point at which the approximation is performed, and the computational comple xity analysis. Numerical examples are presented in Section 5 to demonstrate the e ff ectiv eness of these extensions for set-v alued state estimation of nonlinear discrete- time systems. Finally , Section 6 concludes the paper . 2. Problem f ormulation Consider a class of discrete-time systems with nonlinear dy- namics and linear measurements, described by x k = f ( x k − 1 , u k − 1 , w k − 1 ) , y k = Cx k + D u u k + D v v k , (1) for k ≥ 1, with y 0 = Cx 0 + D u u 0 + D v v 0 , where x k ∈ R n denotes the system state, u k ∈ R n u is a known input, w k ∈ R n w is the process disturbance, y k ∈ R n y is the measured output, and v k ∈ R n v is the measurement uncertainty , with x 0 the initial state. The nonlinear mapping f is assumed to be of class C 2 , and the disturbances and uncertainties are assumed to be bounded, i.e., w k ∈ W k and v k ∈ V k , where W k and V k are known compact sets. This work proposes ne w methods to perform set-valued state estimation for nonlinear systems as in ( 1 ) . The exact characteri- zation of sets X k containing the ev olution of the system states is very di ffi cult in the nonlinear case, if not intractable [ 29 , 36 , 37 ]. Therefore, in the set-membership framework the objecti ve is to enclose such sets as tightly as possible by guaranteed outer bounds ˆ X k on the possible trajectories of the system states x k . 2 Such outer bounds must be consistent with the pre vious estimate ˆ X k − 1 , kno wn inputs u k − 1 , the current measurement y k , and also with the bounds on the disturbances and uncertainties W k − 1 , V k . Giv en an initial condition x 0 ∈ ˆ X 0 , a common approach is to proceed through the well-kno wn prediction-update algorithm, which consists in computing compact sets ¯ X k and ˆ X k such that ¯ X k ⊇ { f ( x , u k − 1 , w ) : x ∈ ˆ X k − 1 , w ∈ W k − 1 } , (2) ˆ X k ⊇ { x ∈ ¯ X k : Cx + D u u k + D v v = y k , v ∈ V k } , (3) in which ( 2 ) is referred to as the pr ediction step , and ( 3 ) as the update step . Our goal is to obtain accurate outer bounds ¯ X k and ˆ X k accord- ing to ( 2 ) and ( 3 ) , respectiv ely . Following these definitions, and considering the initial condition x 0 ∈ ˆ X 0 , the property x k ∈ ˆ X k is guaranteed by construction for all k ≥ 1 [ 9 , 10 , 33 ]. 3. Preliminaries 3.1. Set operations and constr ained zonotopes A few common set operations are often used to compute enclosures satisfying ( 2 ) and ( 3 ) [ 10 , 25 ]. Consider Z , W ⊂ R n , R ∈ R m × n , and Y ⊂ R m . Define the linear mapping, Minkowski sum, and generalized intersection, as R Z , { Rz : z ∈ Z } , (4) Z ⊕ W , { z + w : z ∈ Z , w ∈ W } , (5) Z ∩ R Y , { z ∈ Z : Rz ∈ Y } , (6) respectiv ely . Using ellipsoids or parallelotopes, the linear map- ping ( 4 ) can be computed exactly , but ( 5 ) and ( 6 ) must be ov er- approximated [ 9 , 18 ]. For intervals, the Minkowski sum ( 5 ) is exact, but ( 4 ) and ( 6 ) are conservati ve due to the wrapping e ff ect 1 . In contrast, con vex polytopes are closed under ( 4 ) – ( 6 ) . Moreov er , ( 4 ) and ( 5 ) can be computed e ffi ciently in verte x representation (V -rep), and ( 6 ) can be computed e ffi ciently in half-space representation (H-rep). Howe ver , con version between H-rep and V -rep is computationally e xpensiv e. Zonotopes [ 36 ] allow ( 4 ) and ( 5 ) to be computed exactly and with low computa- tional burden, but ( 6 ) is not a zonotope in general and must be ov er-approximated [ 10 , 38 ]. Constrained zonotopes are an extension of zonotopes, recently proposed in [ 25 ], and are the class of sets of main interest in this work. Definition 1. A set Z ⊂ R n is a constrained zonotope if there exists ( G z , c z , A z , b z ) ∈ R n × n g × R n × R n c × n g × R n c such that Z = { c z + G z ξ : k ξ k ∞ ≤ 1 , A z ξ = b z } . (7) Equation ( 7 ) is called constrained generator r epr esentation (CG-rep), in which each column of G z is a generator , c z is the center , and A z ξ = b z are constraints . In this work, we 1 The generalized intersection in ( 6 ) is not conserv ative when R = I , which corresponds to the standard intersection ∩ . refer to ξ as the generator variables . Let B ∞ ( A z , b z ) = { ξ ∈ R n g : k ξ k ∞ ≤ 1 , A z ξ = b z } and B n g ∞ = { ξ ∈ R n g : k ξ k ∞ ≤ 1 } be, respectiv ely , the constrained unitary hypercube and the n g - dimensional unitary hypercube 2 . Then, a constrained zonotope Z can alternati vely be interpreted as an a ffi ne transformation of B ∞ ( A z , b z ), giv en by Z = c z ⊕ G z B ∞ ( A z , b z ). Note that the linear equality constraints in ( 7 ) allow constrained zonotopes to represent any con vex polytope provided that the comple xity of the CG-rep ( 7 ) is not limited. In fact, Z is a constrained zonotope i ff it is a conv ex polytope [ 25 ]. W e use the compact notation Z = { G z , c z , A z , b z } for constrained zonotopes, and Z = { G z , c z } for zonotopes. In addition to ( 4 ) and ( 5 ) , the intersection ( 6 ) can also be computed exactly with constrained zonotopes. Let Z = { G z , c z , A z , b z } ⊂ R n , W = { G w , c w , A w , b w } ⊂ R n , Y = { G y , c y , A y , b y } ⊂ R m , and R ∈ R n × m . The set operations ( 4 ) – ( 6 ) are computed in CG-rep as R Z = { RG z , Rc z , A z , b z } , (8) Z ⊕ W = ( h G z G w i , c z + c w , " A z 0 0 A w # , " b z b w #) , (9) Z ∩ R Y =          h G z 0 i , c z ,           A z 0 0 A y RG z − G y           ,           b z b y c y − Rc z                    . (10) These operations can be performed e ffi ciently and cause only a moderate increase in the complexity of the CG-rep. Other useful operations with constrained zonotopes are presented in the following. Property 1 provides a simple method for computing the interv al hull of a constrained zonotope by solving 2 n linear programs (LPs), while Proposition 1 provides a way to obtain the closest point in a constrained zonotope to another point in space (in the 1-norm sense) through the solution of a single LP . For simplicity , the subscripts of the v ariables in ( 7 ) will be suppressed henceforth when not necessary . Property 1. [ 25 , 34 ] Let Z = { G , c , A , b } ⊂ R n and let G j denote the j -th row of G . The interval hull [ ζ L , ζ U ] ⊇ Z is obtained by solving the follo wing linear programs for each j = 1 , 2 , . . . , n : ζ L j = min ξ n c j + G j ξ : k ξ k ∞ ≤ 1 , A ξ = b o , ζ U j = max ξ n c j + G j ξ : k ξ k ∞ ≤ 1 , A ξ = b o . Proposition 1. Let Z = { G , c , A , b } ⊂ R n and h ∈ R n . A point z ∈ Z that minimizes k z − h k 1 is giv en by z ∗ = c + G ξ ∗ , where ξ ∗ is a solution to the linear program min ξ k c − h + G ξ k 1 s.t. A ξ = b , k ξ k ∞ ≤ 1 . Pr oof. By definition, k z ∗ − h k 1 = k c − h + G ξ ∗ k 1 2 For simplicity , we drop the use of the superscript n g for B ∞ ( A z , b z ). This dimension can be directly inferred from the number of columns of A z . 3 ≤ k c − h + G ξ k 1 , ∀ ξ ∈ B ∞ ( A , b ) . But, for an y z ∈ Z , there e xists ξ ∈ B ∞ ( A , b ) such that z = c + G ξ . Therefore, k z ∗ − h k 1 ≤ k z − h k 1 , ∀ z ∈ Z .  The presence of the equality constraints in ( 7 ) may result in c < Z (i.e., when 0 < B ∞ ( A , b )). Some of the techniques proposed in this paper require that c ∈ Z (Section 4.3 ). T o accommodate this, Proposition 2 pro vides a simple method for computing an alternati ve (more comple x) CG-rep for Z whose center is any desired point in space. Proposition 2 (Rescaling with desired center) . Let Z = { G , c , A , b } ⊂ R n and let ˜ ξ L , ˜ ξ U ∈ R n g satisfy B ∞ ( A , b ) ⊆ [ ˜ ξ L , ˜ ξ U ]. Choose any desired center h ∈ R n lying in the range of G and let ξ L , ξ U ∈ R n g be solutions to the linear program: min ξ L , ξ U      1 2 ( ξ U − ξ L )      1 s.t. c + 1 2 G ( ξ L + ξ U ) = h , ξ L ≤ ˜ ξ L , ξ U ≥ ˜ ξ U . Letting ξ m = 1 2 ( ξ L + ξ U ) and E r = 1 2 diag ( ξ U − ξ L ), an equi valent CG-rep of Z with center h is giv en by Z =          h GE r 0 i , h ,           AE r 0 0 A GE r − G           ,           b − A ξ m b c − h                    . (11) Pr oof. It is first shown that Z is contained in the set ¯ Z = { GE r , c + G ξ m , AE r , b − A ξ m } . Choose any z ∈ Z . There must e xist ξ ∈ B ∞ ( A , b ) such that z = c + G ξ . Since ξ ∈ [ ξ L , ξ U ], there must exist δ ∈ B n g ∞ such that ξ = ξ m + E r δ . Thus, z ∈ Z = ⇒ ∃ δ ∈ B n g ∞ : z = c + G ( ξ m + E r δ ) , A ( ξ m + E r δ ) = b = ⇒ z ∈ ¯ Z . Therefore, Z ⊆ ¯ Z and it is true that Z = ¯ Z ∩ Z . Since ξ L and ξ U satisfy c + 1 2 G ( ξ L + ξ U ) = c + G ξ m = h , then representing ¯ Z ∩ Z as in ( 10 ) giv es ( 11 ).  Remark 1. The linear program in Proposition 2 does not require [ ξ L , ξ U ] ⊆ B n g ∞ . Therefore, the midpoint 1 2 ( ξ L + ξ U ) can assume any desired v alue, and it is always possible to satisfy c + 1 2 G ( ξ L + ξ U ) = h if h is in the range of G . Thus, the linear program is always feasible. Remark 2. The optimization problems in this section can be readily rewritten as standard form LPs by using additional deci- sion variables and constraints [ 39 ]. 4. Nonlinear state estimation This section presents two new methods for set-v alued state estimation of the class of nonlinear discrete-time systems de- scribed by ( 1 ) . Focusing on the prediction step ( 2 ) , we address the problem of propagating a constrained zonotope ˆ X k − 1 through a nonlinear mapping, with ˆ X 0 , W k , and V k being constrained zonotopes as well. This is done by extending, in a consistent way , two existing approaches for propagating zonotopes through non- linear mappings; namely , the mean value extension in [ 14 , 36 ] and the first-or der T aylor e xtension in [ 26 ]. The methods de- scribed in these works rely , respectiv ely , on the Mean V alue Theorem and T aylor’ s Theorem for the calculation of rigorous outer bounds of the range of the nonlinear mapping in order to obtain zonotopes enclosing the system trajectories. Both methods are based on intersection with strips for performing the update step ( 3 ) . The key adv antage of our new e xtensions is that they allo w the entire state estimation procedure to be done using constrained zonotopes in CG-rep. Therefore, the update ( 3 ) can be done by generalized intersection (with linear measurements), which is kno wn to generate highly asymmetrical sets that cannot be accurately enclosed by ellipsoids, interv als, parallelotopes, or zonotopes. Using the methods dev eloped below , such sets can be directly propagated to the next time step without prior simpli- fication to a symmetric set. This ov ercomes a major source of conservatism in e xisting methods based on the aforementioned enclosures, while largely retaining the e ffi ciency of computa- tions with zonotopes. In addition, our methods e xpand the use of the important tools dev eloped in [ 25 ] to the class of nonlinear discrete-time systems described in Section 2 . In the remainder of the paper , functions with set-valued arguments are consistently used to denote e xact image of the set under the function; e.g., µ ( X , W ) , { µ ( x , w ) : x ∈ X , w ∈ W } . The methods below mak e use of interval arithmetic in sev eral places. Here we recall some of its main concepts. Let IR denote the set of compact intervals in R . Then x ∈ IR is a real compact set defined by x = { a ∈ R : x L ≤ a ≤ x U } , with shorthand notation [ x L , x U ]. The midpoint and radius are defined by mid ( x ) = 1 2 ( x U + x L ) and rad ( x ) = 1 2 ( x U − x L ). The diameter is diam ( x ) = 2 rad ( x ). Interval v ectors and matrices are defined by { a ∈ R n : a L i ≤ a i ≤ a U i } and { A ∈ R n × m : A L i j ≤ A i j ≤ A U i j } , respectiv ely , with midpoint and radius defined component- wise.  ( f ( X ) ) denotes an interv al enclosure of a vector v alued function f ov er X ⊂ R n . The notation  ( f ( X ) ) is used ev en when X is not an interv al. In this case, it is assumed that the interv al hull of X is employed in the operation. Inclusion functions satisfy f ( X ) ⊆  ( f ( X ) ) . See [ 28 ] for definitions of basic interv al arithmetic operations and examples. In addition, the follo wing notations are defined to be used in our proofs. Let κ be a function of class C 2 and z denote its argu- ment. Then, κ q denotes the q -th component of κ , ∇ κ denotes the gradient of κ , and H κ q is an upper triangular matrix describing half of the Hessian of κ q . Specifically , H ii κ q = (1 / 2) ∂ 2 κ q /∂ z 2 i , H i j κ q = ∂ 2 κ q /∂ z i ∂ z j for i < j , and H i j κ q = 0 for i > j . 4.1. Mean value extension This section presents the first of two ne w methods for en- closing the range of a nonlinear function µ ov er a set of inputs described by constrained zonotopes. This method is referred to as the mean value extension of µ (because it relies on the Mean V alue Theorem), and is a consistent generalization of the method for zonotopes in [ 14 ]. Due to significant di ff erences with respect to its zonotopic counterpart, a new theorem (Theorem 2 ) together with a detailed proof is provided for the ne w method. 4 The method in [ 14 ] relies on a zonotope inclusion operator that computes a zonotopic enclosure of the product of an interv al matrix with a unitary box. W e first generalize this operator to constrained zonotopes. Theorem 1. Let X = p ⊕ M B ∞ ( A , b ) ⊂ R m be a constrained zonotope with n g generators and n c constraints, let J ∈ IR n × m be an interval matrix, and consider the set S = J X , { ˆ Jx : ˆ J ∈ J , x ∈ X } ⊂ R n . Let ¯ X = ¯ p ⊕ ¯ M B ¯ n g ∞ be a zonotope satisfying X ⊆ ¯ X , let m be an interval vector such that m ⊇ ( J − mid ( J )) ¯ p and mid ( m ) = 0 , and let P ∈ R n × n be a diagonal matrix defined by P ii = 1 2 diam( m i ) + 1 2 ¯ n g X j = 1 m X k = 1 diam( J ik ) | ¯ M k j | , (12) for all i = 1 , 2 , . . . , n . Then S is contained in the CZ-inclusion S ⊆  ( J , X ) , mid( J ) X ⊕ P B n ∞ . (13) Pr oof. Choose any s ∈ S . It will be shown that s ∈  ( J , X ). By the definition of S , there must exist x ∈ X and ˆ J ∈ J such that s = ˆ Jx . Adding and subtracting mid( J ) x , s = mid( J ) x + ( ˆ J − mid( J )) x . Since x ∈ X ⊆ ¯ X , there exists γ ∈ B ¯ n g ∞ such that x = ¯ p + ¯ M γ . Therefore, s = mid ( J ) x + ( ˆ J − mid ( J ))( ¯ p + ¯ M γ ) . By the choice of m , there must exist ˆ m ∈ m such that s = mid( J ) x + ˆ m + ( ˆ J − mid( J )) ¯ M γ . (14) Let η = ˆ m + ( ˆ J − mid( J )) ¯ M γ . Then η i = ˆ m i + ¯ n g X j = 1 (( ˆ J − mid( J )) ¯ M ) i j γ j , = ˆ m i + ¯ n g X j = 1        m X k = 1 ( ˆ J ik − mid( J ik )) ¯ M k j        γ j . By the triangle inequality and the fact that | γ j | ≤ 1, | η i | ≤ | ˆ m i | + ¯ n g X j = 1        m X k = 1 | ( ˆ J ik − mid( J ik )) || ¯ M k j |        | γ j | , ≤ 1 2 diam( m i ) + 1 2 ¯ n g X j = 1 m X k = 1 diam( J ik ) | ¯ M k j | . Therefore, η ∈ P B n ∞ . From ( 14 ), this implies that s = mid( J ) x + η ∈ mid( J ) X ⊕ P B n ∞ =  ( J , X ) . Thus S ⊆  ( J , X ).  Remark 3. In Theorem 1 , a zonotope satisfying X ⊆ ¯ X can be easily obtained by performing n c iterated constraint eliminations on X using the method in [ 25 ]. Moreover , m can be obtained by simply ev aluating ( J − mid ( J )) ¯ p with interval arithmetic. These methods are used in this work. Finally , the enclosure ( 13 ) has n g + n generators and n c constraints. The following theorem provides the mean v alue extension for constrained zonotopes. Theorem 2. Let µ : R n × R n w → R n be continuously di ff eren- tiable and ∇ x µ denote the gradient of µ with respect to its first argument. Let X ⊂ R n and W ⊂ R n w be constrained zonotopes and choose any h ∈ X . If Z is a constrained zonotope such that µ ( h , W ) ⊆ Z and J ∈ IR n × n is an interv al matrix satisfying ∇ T x µ ( X , W ) ⊆ J , then µ ( X , W ) ⊆ Z ⊕  ( J , X − h ) . Pr oof. Choose any ( x , w ) ∈ X × W . It will be shown that µ ( x , w ) ∈ Z ⊕  ( J , X − h ) . For any i = 1 , 2 , . . . , n , the Mean V alue Theorem ensures that ∃ δ [ i ] ∈ X such that µ i ( x , w ) = µ i ( h , w ) + ∇ T x µ i ( δ [ i ] , w )( x − h ) . But the vector ∇ T x µ i ( δ [ i ] , w ) is contained in the i -th ro w of J by hypothesis, and since this is true for all i = 1 , 2 , . . . , n , there exists a real matrix ˆ J ∈ J such that µ ( x , w ) = µ ( h , w ) + ˆ J ( x − h ) . By Theorem 1 and the choice of Z , it follows that µ ( x , w ) ∈ Z ⊕  ( J , X − h ) , as desired.  Remark 4. The interv al matrix J required by Theorem 2 can be obtained by computing the interval hulls of X and W as in Prop- erty 1 and then bounding ∇ T x µ ( X , W ) using interval arithmetic. Similarly , the constrained zonotope Z ⊇ µ ( h , W ) can be obtained by bounding µ ( h , W ) with interval arithmetic. Alternativ ely , an- other mean v alue extension can be applied around some h w ∈ W to obtain µ ( h , W ) ⊆ Z , µ ( h , h w ) ⊕  ( J w , W − h w ) , where J w is an interv al enclosure of ∇ T w µ ( h , W ). Finally , if µ is a ffi ne in w , i.e, µ ( x , w ) , β x ( x ) + B w ( x ) w , then an exact enclosure of µ ( h , W ) is Z = β x ( h ) ⊕ B w ( h ) W . Since the CG-rep ( 7 ) is an alternati ve representation for con- ve x polytopes [ 25 ], the mean value extension de veloped in Theo- rem 2 provides a ne w method for propagating con ve x polytopes implicitly through nonlinear mappings. A related approach can be found in [ 27 ], where con ve x polytopes are represented by in- tersections of zonotopes (i.e., zonotope bundles [ 24 ]). Howe ver , while e ff ectiv e complexity reduction algorithms are av ailable for constrained zonotopes [ 25 ], e ffi cient methods for comple xity control of zonotope bundles ha ve not yet been proposed. Remark 5. The enclosure obtained in Theorem 2 has at most n g + n g w + 2 n generators and n c + n c w constraints (considering Z computed as in the alternati ves presented in Remark 4 ), with n g and n g w denoting the number of generators of X and W , and n c and n c w the number of constraints, respectiv ely . Thus, the complexity of the resulting set gro ws linearly with respect to the number of constraints and generators. 4.2. F irst-or der T aylor extension This section presents the second ne w method for enclosing the range of a nonlinear function η ov er a set of inputs described by constrained zonotopes. This method is referred to as the first-or der T aylor extension of η (because it relies on a first-order T aylor expansion with a rigorous remainder bound), and is a consistent generalization of the method for zonotopes in [ 26 ]. In contrast to Theorem 2 , for the sake of simplicity of the proof, 5 this function has only one ar gument. Even so, it is possible to consider both states and process uncertainties by concatenating them into a single vector . Due to substantial changes with respect to the zonotopic method, the new approach comes with a ne w theorem (Theorem 3 ) and a detailed proof. In the main result belo w , ( · ) i , : denotes the i -th row of a matrix, and ( · ) i j denotes the element from its i -th row and j -th column. Theorem 3. Let η : R m → R n be of class C 2 and z ∈ R m denote its argument. Let Z = { G , c , A , b } ⊂ R m be a con- strained zonotope with m g generators and m c constraints. For each q = 1 , 2 , . . . , n , let Q [ q ] ∈ IR m × m and ˜ Q [ q ] ∈ IR m g × m g be interval matrices satisfying Q [ q ] ⊇ H η q ( Z ) and ˜ Q [ q ] ⊇ G T Q [ q ] G . Moreov er , define ˜ c q = trace n mid( ˜ Q [ q ] ) o / 2 , ˜ G q , : =  · · · mid( ˜ Q [ q ] ii ) / 2 | {z } ∀ i · · ·  mid( ˜ Q [ q ] i j ) + mid( ˜ Q [ q ] ji )  | {z } ∀ i < j · · ·  , ˜ G d = diag( d ) , d q = X i , j     rad( ˜ Q [ q ] i j )     , ˜ A = h ˜ A ζ ˜ A ξ 0 m c 2 (1 + m c ) × n i , ˜ A ζ =                . . . · · · 1 2 A ri A si · · · . . .                − − − − − − − − − − − − − − − − − → ∀ i                   y ∀ r ≤ s , ˜ A ξ =                . . . · · · A ri A s j + A r j A si · · · . . .                − − − − − − − − − − − − − − − − − − − − − − − − − → ∀ i < j                   y ∀ r ≤ s , ˜ b =                . . . b r b s − 1 2 P i A ri A si . . .                                  y ∀ r ≤ s , with indices i , j = 1 , 2 , . . . , m g and r , s = 1 , 2 , . . . , m c . Finally , choose any h ∈ Z and let L ∈ IR n × m be an interv al matrix satisfying L q , : ⊇ ( c − h ) T Q [ q ] for all q = 1 , . . . , n . Then, η ( Z ) ⊆ η ( h ) ⊕ ∇ T η ( h )( Z − h ) ⊕ R , (15) where R = ˜ c ⊕ h ˜ G ˜ G d i B ∞ ( ˜ A , ˜ b ) ⊕  ( L , ( c − h ) ⊕ 2 G B ∞ ( A , b )). Pr oof. Choose any z ∈ Z and q ∈ { 1 , . . . , n } . By T aylor’ s the- orem applied to η q with reference point h , there must exist Γ [ q ] ∈ H η q ( Z ) ⊆ Q [ q ] such that 3 η q ( z ) = η q ( h ) + ∇ T η q ( h )( z − h ) + ( z − h ) T Γ [ q ] ( z − h ) . 3 Let Υ [ q ] belong to the standard Hessian matrix of η q ( Z ). Then, (1 / 2)( z − h ) T Υ [ q ] ( z − h ) = ( z − h ) T Γ [ q ] ( z − h ) holds. See [ 26 ] for a motivation on this approach. Since z ∈ Z , there must exist ξ ∈ B ∞ ( A , b ) such that z = c + G ξ . Thus, defining p = c − h for brevity , η q ( z ) = η q ( h ) + ∇ T η q ( h )( p + G ξ ) + ( p + G ξ ) T Γ [ q ] ( p + G ξ ) . Expanding the product ( p + G ξ ) T Γ [ q ] ( p + G ξ ) yields p T Γ [ q ] ( p + 2 G ξ ) + ξ T ˜ Γ [ q ] ξ , with ˜ Γ [ q ] = G T Γ [ q ] G ∈ ˜ Q [ q ] . Since ˜ Γ [ q ] ∈ ˜ Q [ q ] , it follows that ˜ Γ [ q ] i j = mid ( ˜ Q [ q ] i j ) + rad ( ˜ Q [ q ] i j ) Λ [ q ] i j for some Λ [ q ] i j ∈ B 1 ∞ . Additionally , ξ i ∈ [ − 1 , 1] implies that ξ 2 i ∈ [0 , 1], and hence ξ 2 i = 1 2 + 1 2 ζ i for some ζ i ∈ [ − 1 , 1]. Considering these two facts, ξ T ˜ Γ [ q ] ξ = 1 2 X i mid( ˜ Q [ q ] ii ) + 1 2 X i mid( ˜ Q [ q ] ii ) ζ i + X i < j (mid( ˜ Q [ q ] i j ) + mid( ˜ Q [ q ] ji )) ξ i ξ j + X i , j rad( ˜ Q [ q ] i j ) ξ i ξ j Λ [ q ] i j , where the third summation results from the fact that ξ i ξ j = ξ j ξ i . Thus, by defining the new generator v ariables ¯ ξ =  · · · ζ i | {z } ∀ i · · · ξ i ξ j | {z } ∀ i < j · · · ξ i ξ j Λ [ q ] i j | {z } ∀ i , j , q · · ·  T , with i , j = 1 , 2 , . . . , m g , q = 1 , 2 , . . . , n , we have that ξ T ˜ Γ [ q ] ξ = ˜ c q + [ ˜ G ¯ G d ] q , : ¯ ξ , where ¯ G d = blkdiag( N [1] , N [2] , . . . , N [ n ] ) 4 , N [ q ] =  · · · rad( ˜ Q [ q ] i j ) | {z } ∀ i , j · · ·  ∈ R 1 × m 2 g . Therefore, we have established that η q ( z ) = η q ( h ) + ∇ T η q ( h )( z − h ) + p T Γ [ q ] ( p + 2 G ξ ) + ˜ c q + [ ˜ G ¯ G d ] q , : ¯ ξ . This holds for e very q = 1 , 2 , . . . , n . Moreover , L satisfies L q , : ⊇ p T Q [ q ] for all q = 1 , 2 , . . . , n by definition, so there must e xist ˆ L ∈ L such that ˆ L q , : = p T Γ [ q ] for all q = 1 , 2 , . . . , n . Therefore, η ( z ) = η ( h ) + ∇ T η ( h )( z − h ) + ˆ L ( p + 2 G ξ ) + ˜ c + [ ˜ G ¯ G d ] ¯ ξ . (16) Furthermore, the equality constraints A ξ = b imply that A ξ ξ T A T = bb T . Thus, considering ξ 2 i = 1 2 + 1 2 ζ i , the r -th row and s -th column of this matrix equality yields 1 2 X i A ri A si ζ i + X i < j ( A ri A s j + A r j A si ) ξ i ξ j = b r b s − 1 2 X i A ri A si , with r , s = 1 , 2 , . . . , m c . Such constraints are linear in ¯ ξ , and non-repeating for r ≤ s , therefore ¯ A ¯ ξ = ˜ b holds, where ¯ A = [ ˜ A ζ ˜ A ξ 0 ˜ m c × nm 2 g ], with ˜ m c = m c 2 (1 + m c ). Hence ξ ∈ B ∞ ( A , b ) = ⇒ ¯ ξ ∈ B ∞ ( ¯ A , ˜ b ). Combining this with ( 16 ) , we hav e prov en the enclosure η ( Z ) ⊆ η ( h ) ⊕ ∇ T η ( h )( Z − h ) ⊕ L ( p ⊕ 2 G B ∞ ( A , b )) ⊕ ˜ c ⊕ [ ˜ G ¯ G d ] B ∞ ( ¯ A , ˜ b ). In fact, this enclosure can be greatly simplified by not- ing that the columns of ¯ A corresponding to the variables [ · · · ξ i ξ j Λ [ q ] i j · · · ] are all zero, and hence [ ˜ G ¯ G d ] B ∞ ( ¯ A , ˜ b ) = ˜ G B ∞ ([ ˜ A ζ ˜ A ξ ] , ˜ b ) ⊕ ¯ G d B nm 2 g ∞ . 4 In this work blkdiag ( A , B , . . . ) denotes a block diagonal matrix with blocks A , B , . . . . 6 Since ¯ G d is block diagonal and each N [ q ] is a ro w v ector , ¯ G d B nm 2 g ∞ is an interval, and is equi valent to ˜ G d B n ∞ , with ˜ G d defined as in the statement of the theorem. Thus, [ ˜ G ¯ G d ] B ∞ ( ¯ A , ˜ b ) = ˜ G B ∞ ([ ˜ A ζ ˜ A ξ ] , ˜ b ) ⊕ ˜ G d B n ∞ , = [ ˜ G ˜ G d ] B ∞ ([ ˜ A ζ ˜ A ξ 0 ˜ m c × n ] , ˜ b ) = [ ˜ G ˜ G d ] B ∞ ( ˜ A , ˜ b ) . Therefore, η ( Z ) ⊆ η ( h ) ⊕ ∇ T η ( h )( Z − h ) ⊕ L ( p ⊕ 2 G B ∞ ( A , b )) ⊕ ˜ c ⊕ [ ˜ G ˜ G d ] B ∞ ( ˜ A , ˜ b ) , and ( 15 ) follows immediately from the definition of R .  Remark 6. Regarding the definitions of ˜ G , ˜ A ζ , ˜ A ξ , and ˜ b in Theorem 3 , the ordering of the indices i < j and r ≤ s is irrelev ant, as long as it is the same for all variables. Remark 7. The interval matrices Q [ q ] required by Theorem 3 can be obtained by computing the interv al hull of Z (Property 1 ) and then bounding H η q ( Z ) using interval arithmetic. Moreov er , ˜ Q [ q ] and L can be obtained by e valuating G T Q [ q ] G and ( c − h ) T Q [ q ] using interval arithmetic. As stated before, process disturbances can be taken into ac- count in ( 15 ) by considering the augmented v ector z = ( x , w ) with Z = X × W ⊂ R n + n w and h = ( h x , h w ) ∈ Z . W ith X = { G x , c x , A x , b x } and W = { G w , c w , A w , b w } , the Cartesian product Z is easily computed by X × W = (" G x 0 0 G w # , " c x c w # , " A x 0 0 A w # , " b x b w #) . Remark 8. In Theorem 3 , ˜ G has P m g j = 1 j = 1 2 m g ( m g + 1) columns, ˜ G d ∈ R n × n , ˜ A has P m c s = 1 s = 1 2 m c ( m c + 1) rows, and  ( L , ( c − h ) ⊕ 2 G B ∞ ( A , b )) has m g + n generators and m c constraints (Remark 3 ). Therefore, the resulting enclosure in ( 15 ) has 1 2 m 2 g + 5 2 m g + 2 n generators and 1 2 m 2 c + 5 2 m c constraints. If Z = X × W , then the enclosure has 1 2 ( n g + n g w ) 2 + 5 2 ( n g + n g w ) + 2 n generators and 1 2 ( n c + n c w ) 2 + 5 2 ( n c + n c w ) constraints, which is a polynomial increase in complexity in terms of both generators and constraints. 4.3. Selection of h The methods proposed in the pre vious sections require a choice of h ∈ X = { G x , c x , A x , b x } in order to compute a con- strained zonotope enclosure for the prediction step ( 2 ) . As shown in Section 5.1 , this choice may drastically a ff ect the accurac y of the obtained enclosure. In the mean value extension for inter- vals and zonotopes, a usual choice of h ∈ X is the center of X [ 14 , 28 ]. Howe ver , with constrained zonotopes, since the center of the CG-rep may not belong to X 5 , a di ff erent point h ∈ X must be chosen. A simple and inexpensi ve choice is the center of the interval hull of X . Unfortunately , ev en this point may not belong to X in some cases 6 . Nevertheless, this choice can 5 From X = c x ⊕ G x B ∞ ( A x , b x ), c x < X as long as @ ξ ∈ B ∞ ( A x , b x ) satisfying G x ξ = 0 . 6 An example is the polytope with vertices (0 , 0 , 0), (1 , 1 , 0), (0 , 1 , 0), (0 , 1 , 1). be applied rigorously by simply checking h ∈ X beforehand by solving an LP [ 25 ]. In the follo wing, we analyze alternati ve choices of h ∈ X . Firstly , we focus on suitable choices v alid for the mean value e x- tension (Theorem 2 ). This extension relies on the CZ-inclusion (Theorem 1 ), and therefore requires the computation of a zono- tope enclosing X − h . In this w ork, we assume that this zonotope is computed through constraint eliminations (see Remark 3 ). Let { G (  ) , c (  ) , A (  ) , b (  ) } denote the constrained zonotope obtained by reducing to  the number of remaining constraints in X − h . Follo wing the constraint elimination algorithm in [ 25 ], for each  = n c , n c − 1 , . . . , 1, the remaining constraints A (  ) ξ = b (  ) are first preconditioned through Gauss-Jordan elimination with full piv oting and then subjected to a rescaling procedure before the next constraint is eliminated. The entire procedure can be repre- sented by the follo wing recursive equations (see Proposition 5 and the Appendix in [ 25 ] for details), where ¯ ( · ) denotes v ariables after preconditioning, ˜ ( · ) denotes v ariables after rescaling, and Λ G , Λ A , ξ m , and ξ r are defined as in [ 25 ]: ˜ c (  ) = c (  ) + ¯ G (  ) ξ (  ) m , c (  − 1) = ˜ c (  ) + Λ (  ) G ˜ b (  ) , ˜ G (  ) = ¯ G (  ) diag( ξ (  ) r ) , G (  − 1) = ˜ G (  ) − Λ (  ) G ˜ A (  ) , ˜ A (  ) = ¯ A (  ) diag( ξ (  ) r ) , A (  − 1) = ˜ A (  ) − Λ (  ) A ˜ A (  ) , ˜ b (  ) = ¯ b (  ) − ¯ A (  ) ξ (  ) m , b (  − 1) = ˜ b (  ) − Λ (  ) A ˜ b (  ) . (17) Careful examination of the algorithm in [ 25 ] re veals that the actions taken during preconditioning, rescaling, and constraint elimination are all independent of the center of the original con- strained zonotope, which in this case is c ( n c ) = c x − h . Therefore, with exception of the center , the variables ( · ) (  ) can be obtained by eliminating the constraints of X prior to choosing h . Consider- ing procedure ( 17 ) , the following corollary pro vides a choice of h that leads to a tight enclosure by reducing the conservati veness of the CZ-inclusion  ( J , X − h ). Corollary 1. Let X = { G x , c x , A x , b x } ⊂ R n , and consider µ , W , and J as defined in Theorem 2 . Assume that ¯ G (  ) , ξ (  ) m , Λ (  ) G , and ˜ b (  ) are obtained by eliminating all n c constraints from X according to ( 17 ) and set h = c x + n c X  = 1  ¯ G (  ) ξ (  ) m + Λ (  ) G ˜ b (  )  . (18) Let ¯ X = { G (0) , c (0) } be obtained by eliminating all n c constraints from X − h according to ( 17 ) , let m ⊇ ( J − mid ( J )) c (0) be computed by standard interv al arithmetic, and suppose that  ( J , X − h ) is computed as in Theorem 1 with this choice of ¯ X and m . Finally , let Z ⊇ µ ( h , W ). If h ∈ X , then µ ( X , W ) ⊆ Z ⊕  ( J , X − h ) . Moreover ,  ( J , X − h ) ⊆  ( J , X − ˆ h ) for any ˆ h ∈ X , ˆ h , h . Pr oof. For h ∈ X , µ ( X , W ) ⊆ Z ⊕  ( J , X − h ) follows directly from Theorem 2 . No w , let us sho w that  ( J , X − h ) ⊆  ( J , X − ˆ h ) holds for any ˆ h ∈ X , ˆ h , h . Recursiv e computation of ( 17 ) leads to c (0) = c ( n c ) + n c X  = 1  ¯ G (  ) ξ (  ) m + Λ (  ) G ˜ b (  )  , (19) 7 where c ( n c ) = c x − h . Therefore, c (0) = 0 i ff h is given by ( 18 ) , thus m = 0 , and diam ( m ) = 0 . Note that in ( 12 ) , ¯ M , G (0) is in variant with respect to h , and since J ⊇ ∇ T x µ ( X , W ), then J is also in variant with respect to h . Consequently , the second term in ( 12 ) is not a function of h . Therefore, P B n ∞ ⊆ ˆ P B n ∞ = (1 / 2) diag ( diam ( ˆ m )) B n ∞ ⊕ P B n ∞ , with ˆ m computed using ˆ h ∈ X . The result then follows from ( 13 ).  By Corollary 1 , the enclosure obtained in Theorem 2 is tight- ened by choosing h such that c (0) is equal to zero. Unfortunately , the h giv en by ( 18 ) may not belong to X , so an alternati ve to obtain tight bounds is to reduce the size of the box m by solving min h {k diam( m ) k 1 : h ∈ X } , (20) with m ⊇ ( J − mid ( J )) ¯ p computed using interval arithmetic, where ¯ p , c (0) . Recall that c (0) is the center of the zonotope obtained by eliminating all the constraints of X − h . Lemma 1. Let X = { G x , c x , A x , b x } ⊂ R n , J ∈ IR n × n . Assume that ¯ G (  ) , ξ (  ) m , Λ (  ) G , and ˜ b (  ) are obtained by eliminating all n c constraints of X according to ( 17 ) . Then h = c x + G x ξ ∗ is the solution to ( 20 ) i ff ξ ∗ is the solution to the linear program min ξ k Θ ¯ p k 1 , s.t. A x ξ = b x , k ξ k ∞ ≤ 1 , (21) with ¯ p = − G x ξ + P n c  = 1  ¯ G (  ) ξ (  ) m + Λ (  ) G ˜ b (  )  , Θ j j = P n i = 1 diam ( J i j ), and Θ i j = 0 for i , j . Pr oof. Each element of ( J − mid ( J )) ∈ IR n × n is a symmet- ric interv al satisfying ( J i j − mid ( J i j )) = (1 / 2) diam ( J i j )[ − 1 , 1], and for every a ∈ R , a [ − 1 , 1] = | a | [ − 1 , 1] holds. Therefore m i = P n j = 1 (1 / 2) diam ( J i j ) | ¯ p j | [ − 1 , 1]. Consequently , diam ( m i ) = P n j = 1 diam( J i j ) | ¯ p j | , and k diam( m ) k 1 = n X i = 1 n X j = 1 diam( J i j ) | ¯ p j | = n X j = 1        n X i = 1 diam( J i j )        | ¯ p j | = n X j = 1 Θ j j | ¯ p j | = k Θ ¯ p k 1 . The equality ¯ p = − G x ξ + P n c  = 1  ¯ G (  ) ξ (  ) m + Λ (  ) G ˜ b (  )  and the con- straints in ( 21 ) follow directly from ( 19 ) and h ∈ X .  Lemma 1 yields an optimal choice of h ∈ X that can be used in Theorem 2 to reduce conservatism in the CZ-inclusion, and requires only the solution of an LP . Note that formulating ( 21 ) requires the knowledge of ¯ G (  ) , ξ (  ) m , Λ (  ) G , and ˜ b (  ) , which are obtained from the iterated constraint elimination process. As stated before, constraint elimination can be performed o ver X to obtain the required data prior to the solution of ( 21 ) . Once the optimal h is obtained, constraint elimination can be repeated, or equiv alently , the zonotope obtained using h = 0 can simply be translated by − h . Remark 9. Note that if the h giv en by Corollary 1 belongs to X , then this coincides with the solution provided by Lemma 1 . W e summarize the proposed choices of h ∈ X for use in Theorem 2 as follows: C1) h is giv en by the center of the interv al hull of X , if it satisfies h ∈ X ; C2) h is obtained by solving ( 21 ). W e now focus on suitable choices valid for the first-order T aylor extension. As with the mean v alue extension, the usual choice of h ∈ X in first-order T aylor e xtensions for intervals and zonotopes is the center of X [ 28 , 26 ]. The next corollary shows that this choice leads to a tight enclosure if it belongs to X . Corollary 2. Let Z = { G , c , A , b } = X × W ⊂ R m , and consider η , ˜ c , ˜ G , ˜ G d , ˜ A , ˜ b , L , and Q [ q ] as defined in Theorem 3 , with q ∈ { 1 , 2 , . . . , n } . If h = c ∈ Z , then η ( Z ) ⊆ η ( h ) ⊕ ∇ T η ( h )( Z − h ) ⊕ ˜ c ⊕ h ˜ G ˜ G d i B ∞ ( ˜ A , ˜ b ). Pr oof. For h = c , L q , : ⊇ ( c − h ) T Q [ q ] = 0 , q = 1 , 2 , . . . , n . Therefore L = 0 holds, and  ( L , ( c − h ) ⊕ 2 G B ∞ ( A , b )) = 0 . The result then follows from ( 15 ).  By inspecting Corollary 2 , it is clear that the enclosure in ( 15 ) is tightened since  ( L , ( c − h ) ⊕ 2 G B ∞ ( A , b )) = 0 . Howe ver , since this point may not belong to X , a good alternative may be to consider the closest point in X to its center, obtained by means of Proposition 1 . By the definition of L , this heuristic leads to smaller values of diam ( L ), and therefore reduces the size of  ( L , ( c − h ) ⊕ 2 G B ∞ ( A , b )) (see ( 12 ) ). A third option is to apply Proposition 2 to obtain an alternative CG-rep of X with any desired center . In this case, the ne w center is chosen as some point in X , ¯ h ∈ X , and then h is chosen as h = ¯ h . A simple choice of new center ¯ h ∈ X for such a procedure is the center of the interval hull of X . The proposed alternati ves for use in Theorem 3 are summarized as follows: C3) h is giv en by the closest point in X to the center of X , computed through Proposition 1 ; C4) h is the center of X , if it satisfies h ∈ X . Otherwise, ¯ h ∈ X is chosen as the center of the interv al hull of X , the center of X is mov ed to ¯ h using Proposition 2 , and h is gi ven by h = ¯ h . 7 4.4. Update step An enclosure for the prediction step ( 2 ) can be obtained in CG-rep using either Theorem 2 or Theorem 3 . Therefore, due to linearity of the measurement in ( 1 ) , an exact bound for the update step ( 3 ) can be directly obtained by computing the generalized intersection of two constrained zonotopes as follo ws. Given the prediction set ¯ X k , a constrained zonotope V describing bounds on measurement errors, the current input u k and measurement y k , an exact enclosure for the update step is obtained using the definition ( 6 ), giv en by ˆ X k = ¯ X k ∩ C (( y k − D u u k ) ⊕ ( − D v V )) . (22) 7 Note that this choice may lead to the same v alue of h provided by C1 , b ut X is described by a di ff erent CG-rep with center h . If this point is not in X , ¯ h can be chosen as the point obtained from Proposition 1 instead. 8 It is well known that the intersection in ( 22 ) can not be computed exactly using zonotopes, and must be over -approximated [ 10 , 14 ]. As a consequence, the enclosures of the system states obtained after many iterations of prediction and update may be quite conservati ve using zonotopes. Howe ver , with constrained zonotopes all operations in ( 22 ) are easily computed through ( 8 ) – ( 10 ) . These lead to an enclosure with n g + n g v generators, and n c + n c v + n y constraints, where n g and n c are the number of generators and constraints of ¯ X k , respectiv ely Remark 10. Iterated computations of the proposed extensions (Theorems 2 and 3 ) and ( 22 ) result in at most a quadratic in- crease in the complexity of the CG-rep ( 7 ) . As with zonotopes, this can be e ff ecti vely addressed using order r eduction algo- rithms [ 40 , 41 ] that o ver-approximate a constrained zonotope by another with lo wer complexity . E ffi cient methods for reducing the number of generators and constraints of the CG-rep ( 7 ) with reasonable conservati veness were proposed in [ 25 ]. 4.5. Complexity analysis T able 1 sho ws the computational complexity 8 of our meth- ods for the prediction and update steps, as well as complexity reduction to the same number of generators and constraints of the set prior to prediction. Specifically , we use the mean v alue extension (Theorem 2 ) and the first-order T aylor extension (The- orem 3 ) for the prediction steps, while the update steps are both giv en by the generalized intersection ( 22 ) . These methods are denoted by CZMV and CZFO, respectiv ely . The computational complexities of their zonotope counterparts are also presented for comparison, denoted analogously by ZMV and ZFO, which use the mean value approach in [ 14 ] and the first-order T ay- lor approach in [ 26 ], respectiv ely , for the prediction step. The update algorithm proposed in [ 38 ] is used for both ZMV and ZFO because it provided the best trade-o ff between accuracy and e ffi ciency in our numerical experiments with zonotopes. Complexity reduction is applied after the update step in all four methods using the reduction methods in [ 25 ] for constrained zonotopes and Method 4 in [ 41 ] for zonotopes. For constrained zonotopes, constraint elimination is performed prior to generator reduction. The complexities in T able 1 take into account the growth of the number of generators and constraints after each step (see Remarks 5 and 8 ). The dimensions in T able 1 are specified by the definitions ˆ X k − 1 = { G x , c x } , W = { G w , c w } , and V = { G v , c v } or ˆ X k − 1 = { G x , c x , A x , b x } , W = { G w , c w , A w , b w } , and V = { G v , c v , A v , b v } with G x ∈ R n × n g , c x ∈ R n , A x ∈ R n c × n g , b x ∈ R n c , G w ∈ R n w × n g w , c w ∈ R n w , A w ∈ R n c w × n g w , b w ∈ R n c w , G v ∈ R n v × n g v , c v ∈ R n v , A v ∈ R n c v × n g v , b v ∈ R n c v , u k ∈ R n u , and y k ∈ R n y . For simplicity , we define m = n + n w , m g = n g + n g w , m c = n c + n c w , δ n = n − n y , δ w = n g w − n c w , δ v = n g v − n c v , and ˜ δ = m 2 g − m 2 c . Moreo ver , we consider that scalar real function and scalar inclusion function ev aluations have complexity O (1). These correspond to ev aluations of the nonlinear dynamics in ( 1 ) and its deriv atives using real and interv al arithmetic, respec- tiv ely . The comple xities of the basic operations on zonotopes 8 W e use the standard O ( · ) notation defined in [ 42 ]. and constrained zonotopes used to deri ve the figures in T able 1 can be found in the Appendix. The dominant terms in the prediction step of CZMV come from the computation of the interval hulls of X and W and the CZ-inclusions  ( J , X − h ) and  ( J w , W − h w ) in Theorem 2 and Remark 4 . In the case of CZFO, the dominant terms come from the computation of the interv al matrices ˜ Q [ q ] , the interv al hull of Z = X × W , and the CZ-inclusion  ( L , ( c − h ) ⊕ 2 G B ∞ ( A , b )) in Theorem 3 . Note that the worst-case complexities of the predic- tion steps of our methods are higher than the zonotope methods, while the update steps are cheaper due to the generalized inter - section ( 22 ) . Even so, the complexity of the proposed methods are still polynomial. For a simplified analysis, assuming that all of the variables in T able 1 increase linearly with n , the total com- plexities for ZMV , ZFO, CZMV and CZFO are O ( n 4 ), O ( n 5 ), O ( n 5 ), and O ( n 8 ), respecti vely . On the other hand, even basic polytope operations are known to be e xponential [ 43 ]. Besides, despite the higher comple xities of CZMV and CZFO in com- parison to the zonotope methods, the y provide more accurate enclosures as shown in the ne xt section. 5. Numerical examples This section presents numerical results for the two new set- valued state estimation methods enabled by the results in the previous section. The imposed limits on the comple xity of the sets used are described separately for each example belo w . 5.1. Example 1 T o demonstrate the e ff ect of the di ff erent choices of h , we first analyze one iteration of the prediction step for the nonlinear system [ 44 ] x 1 , k = 3 x 1 , k − 1 − x 2 1 , k − 1 7 − 4 x 1 , k − 1 x 2 , k − 1 4 + x 1 , k − 1 + w 1 , k − 1 , x 2 , k = − 2 x 2 , k − 1 + 3 x 1 , k − 1 x 2 , k − 1 4 + x 1 , k − 1 + w 2 , k − 1 , (23) with X 0 = (" 0 . 2 0 . 4 0 . 2 0 . 2 0 − 0 . 2 # , " − 1 1 # , h 2 2 2 i , − 3 ) , (24) where w k ∈ R 2 denotes process uncertainties, which are zero in this first scenario. Figure 1 sho ws the constrained zonotope X 0 and the enclo- sures of the one-step reachable set obtained by Theorem 2 using C1 – C4 . Since the complexity of the enclosure for C4 is higher than for the other methods (see Proposition 2 ), the reduction methods in [ 25 ] were used to reduce the number of generators and constraints in this enclosure to match the other methods before comparison. In this example, the choice of h has a mod- erate impact in the enclosure computed by Theorem 2 , with C2 providing the least conserv ative result, as expected. Therefore, C2 is employed in Theorem 2 henceforth. Figure 1 also sho ws the enclosures of the one-step reachable set obtained by Theorem 3 with C1 – C4 . Clearly , the enclosures 9 T able 1: Computational complexity O ( · ) of the state estimators. Step ZMV ZFO Prediction n 2 n g + nn w n g w n ( m 2 m g + mm 2 g ) Update n y ( n 3 ( m g + n ) + n 2 ( m g + n ) 2 + n u + n v n g v ) n y ( n 3 ( m 2 g + n ) + n 2 ( m 2 g + n ) 2 + n u + n v n g v ) Reduction n 2 ( m g + n ) + n ( n g w + n )( m g + n ) n 2 ( m 2 g + n ) + n ( m 2 g + n ) 2 Step CZMV CZFO Prediction n 2 n g + nn w n g w + ( nn g + n c )( n g + n c ) 3 + ( n w n g w + n c w )( n g w + n c w ) 3 n ( m 2 m g + mm 2 g ) + ( mm g + m c )( m g + m c ) 3 Update n y n ( m g + n ) + n y n u + n y n v n g v n y n ( m 2 g + n ) + n y n u + n y n v n g v Reduction ( n c w + n c v + n y )( m g + m c + n g v + n c v + n + n y ) 3 ( m 2 c + n c v + n y )( m 2 g + m 2 c + n g v + n c v + n + n y ) 3 + ( n + n c ) 2 ( n g + δ n + δ w + δ v ) + ( n + n c )( δ n + δ w + δ v )( n g + δ n + δ w + δ v ) + ( n + n c ) 2 ( ˜ δ + δ n + δ v ) + ( n + n c )( ˜ δ + δ n + δ v ) 2 x 1 x 2 x 1 x 2 x 1 x 2 -5 -4 -3 -2 -1 0 -6 -5 -4 -3 -2 -1 -5 -4 -3 -2 -1 0 -6 -5 -4 -3 -2 -1 -1.5 -1.2 -1 0.7 1 1.3 Figure 1: T op: The constrained zonotope X 0 with ‘ × ’ denoting its center . Left: Enclosures obtained by applying Theorem 2 to ( 23 ) . Right: Enclosures obtained by applying Theorem 3 to ( 23 ) . The real vector h is determined by C1 (green), C2 (red), C3 (yellow), and C4 (blue). For the mean value extension (left), C4 is overlapped with C1 . Black dots denote uniform samples from X 0 propagated through ( 23 ). produced by Theorem 3 are strongly a ff ected by the choice of h , with C4 providing the least conserv ativ e result. In addition, note that the enclosures provided by the first-order T aylor ex- tension are more conserv ativ e than those obtained by the mean value e xtension. Ho wev er , experience with zonotopes and inter - vals (see [ 44 ] for detailed examples) suggests that the relativ e merits of these two methods will depend on the dynamics of the system, as well as the shape and size of the set X 0 , and the maximum allo wed number of generators and constraints. This is corroborated by the next results. W e consider now the linear measurement equation " y 1 , k y 2 , k # = " 1 0 − 1 1 # " x 1 , k x 2 , k # + " v 1 , k v 2 , k # , (25) with bounds k w k k ∞ ≤ 0 . 4 and k v k k ∞ ≤ 0 . 4, where v k ∈ R 2 denote measurement uncertainties. The initial states x 0 are bounded by the zonotope 9 X 0 = (" 0 . 1 0 . 2 − 0 . 1 0 . 1 0 . 1 0 # , " 0 . 5 0 . 5 #) . (26) 9 Note that X 0 , W and V are expressed as zonotopes for the purpose of a fair comparison with the zonotope methods. x 1 x 2 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 0 0.5 1 1.5 Figure 2: The initial set X 0 (gray), the initial bounded uncertain measurements (dashed lines), the intersection computed as in [ 38 ] (yello w), and the constrained zonotope (green) computed by ( 22 ). T o generate process measurements, ( 23 ) was simulated with x 0 = (0 . 8 , 0 . 65) ∈ X 0 and process and measurement uncertainties drawn from uniform random distributions. The number of gen- erators and constraints of the constrained zonotopes was limited to 20 and 5, respectiv ely , while the number of generators of the zonotopes was limited to 20. Figure 2 sho ws the results of the initial update step using the intersection algorithm in [ 38 ] and the generalized intersection ( 22 ) computed using ( 10 ) , which yields a constrained zonotope. Clearly , since the generalized intersection is not a symmetric set, it cannot be described by a zonotope. In contrast, the resulting constrained zonotope corre- sponds to the exact intersection, pro viding far less conserv ativ e bounds in the first update step. Figure 3 sho ws the first four time steps of CZMV with h giv en by C2 in a scenario without process uncertainties ( w k = 0 ). For comparison, the zonotopes computed using ZMV are also depicted. CZMV provides much less conservativ e enclosures than ZMV for this example, demonstrating the e ff ecti veness of the proposed nonlinear state estimation strategy . Figure 4 sho ws the radii (half the length of the longest edge of the interv al hull) of the sets provided by CZMV using C2 and ZMV over 100 time steps considering process disturbances. Since ( 23 ) is a ffi ne in w k , the enclosure Z ⊇ µ ( h , W ) in Theorem 2 was computed as described at the end of Remark 4 . CZMV provided less conserv ativ e bounds than the zonotopes computed by ZMV, with a CZMV -to-ZMV average radius ratio (ARR, i.e., the ratio of the radius of the CZMV set at k ov er the ZMV set at k av eraged o ver all time steps k ) of only 51 . 4%. Figure 4 10 x 1 x 2 0 2 4 6 8 10 12 14 16 18 -1.5 -1 -0.5 0 0.5 1 1.5 k = 0 k = 1 k = 2 k = 3 k = 4 Figure 3: Results from the first four time steps of set-valued state estimation (after update) using the constrained zonotopic method CZMV (green) and the zonotopic method ZMV (yellow). Black dots denote uniform samples from X 0 propagated through ( 23 ) that are consistent with the current measurement. k Radius 0 10 20 30 40 50 60 70 80 90 100 0 0.2 0.4 0.6 0.8 1 Figure 4: Radii of the update sets obtained by applying CZMV ( 4 ), ZMV ( ◦ ), CZFO ( + ), and ZFO ( × ) to ( 23 ) with process disturbances. also compares the radii of the update sets computed by ZFO and CZFO with h giv en by C4 . As in the pre vious case, CZFO provides less conservati ve bounds than ZFO, with the CZFO- to-ZFO ARR being only 53 . 66%. The size of the sets provided by CZMV and CZFO were quite similar, with CZFO being less conservati ve (the CZFO-to-CZMV ARR was 98 . 75%). In both experiments, the number of generators and constraints were limited to 20 and 5, respecti vely . The ARR for di ff erent numbers of constraints are shown in T able 2 , with the av erage computed considering in addition simulations with di ff erent numbers of generators. Execution times are sho wn in T able 3 . These were obtained using MA TLAB 9.1 with CPLEX 12.8 and INTLAB 9, in a laptop with 8GB RAM and an Intel Core i7 4510U 3.1 GHz processor . The use of constrained zonotopes in CZMV and CZFO results in sets that are slightly more complex than those generated by ZMV and ZFO (specifically , the set description inv olves fiv e equality constraints that are not present in the zonotopes from ZMV and ZFO). Howe ver , this example shows that this increase in complexity is compensated by greatly impro ved accuracy . T able 2: A verage radius ratio of the estimators with varying numbers of con- straints. Each a verage is tak en ov er 3 separate simulations using n g ∈ { 8 , 12 , 20 } . n c CZMV / ZMV CZFO / ZFO CZFO / CZMV 1 54 . 1% 59 . 8% 104 . 5% 3 51 . 6% 54 . 0% 99 . 0% 5 51 . 6% 53 . 7% 98 . 5% T able 3: A verage total times per iteration of the estimators with v arying numbers of constraints. Each average is taken over 30 separate simulations using n g ∈ { 8 , 12 , 20 } . T imes spent only on complexity reduction are shown in parenthesis. n c ZMV ZFO CZMV CZFO 0 5 . 79 (0 . 24) ms 14 . 03 (2 . 93) ms – – 1 – – 19 . 2 (1 . 6) ms 76 . 4 (57 . 1) ms 3 – – 20 . 7 (1 . 8) ms 4 . 28 (4 . 26) s 5 – – 22 . 5 (1 . 9) ms 8 . 93 (8 . 91) s 5.2. Example 2 Consider the quadrotor unmanned aerial vehi- cle (U A V) described in [ 45 ] with state vector ζ = [ x y z u v w φ θ ψ p q r ] T , where [ x y z ] T is the position of the U A V with respect to the inertial frame, [ u v w ] T is the velocity v ector expressed in the inertial frame, [ φ θ ψ ] T are Euler angles describing the orientation of the U A V , and [ p q r ] T is the angular velocity v ector expressed in the body frame. The equations of motion are [ 45 ]: ˙ ζ =                                                                                                    ˙ x = u , ˙ y = v , ˙ z = w , ˙ u = 1 m (cos ψ sin θ cos φ + sin ψ sin φ ) U 1 + 1 m D x , ˙ v = 1 m (sin ψ sin θ cos φ − cos ψ sin φ ) U 1 + 1 m D y , ˙ w = − g + 1 m (cos θ cos φ ) U 1 + 1 m D z , ˙ φ = p + q sin φ tan θ + r cos φ tan θ , ˙ θ = q cos φ − r sin φ, ˙ ψ = q sin φ sec θ + r cos φ sec θ , ˙ p = I yy − I zz I x x qr + l I x x U 2 , ˙ q = I zz − I x x I yy pr + l I yy U 3 , ˙ r = I x x − I yy I zz pq + 1 I zz U 4 , (27) where m , I x x , I yy , I zz , and l are physical parameters, g is the gravitational acceleration, U 1 is the total thrust generated by the propellers, U 2 is the di ff erence of thrusts between the left and right propellers, U 3 is the di ff erence of thrusts between the front and back propellers, U 4 is the di ff erence of torques between clockwise and counter-clockwise turning propellers, and d = [ D x D y D z ] T are disturbance forces applied to the U A V with k d k ∞ ≤ 1. The experiment consists in obtaining guaranteed bounds on the system states ζ while the quadrotor U A V tracks a vertical helix trajectory defined by the reference 11 values x ref ( t ) = 1 2 cos  t 2  , y ref ( t ) = 1 2 sin  t 2  , z ref ( t ) = 1 + t 10 , ψ ref = π 3 , subject to the disturbance forces described by D x = 1 N for t ∈ [5 , 15) s, D y = 1 N for t ∈ [8 , 15) s, and D z = 1 N for t ∈ [10 , 15) s. These forces are zero otherwise. The dynamic feedback controller in [ 45 ] is used to track the reference trajectory above 10 . The simulation parameters are m = 0 . 7 Kg, l = 0 . 3 m, I x x = I yy = I zz = 1 . 2416 Kg · m 2 , g = 9 . 81 m / s 2 . W e consider a realistic scenario in which the av ailable measurements are provided by sensors located at the quadrotor U A V , which include: (i) a Global Positioning System (GPS); (ii) a barometer; and (iii) an Inertial Measurement Unit (IMU). The measurements are a ff ected by bounded uncertainties as described in T able 4 . The velocity vector [ u v w ] T is not measured. T able 4: Measured variables with error bounds. Sensor V ariables Noise bounds GPS { x , y } ± 0 . 15 m Barometer { z } ± 0 . 51 m IMU { φ, θ, ψ } ± 2 . 618 · 10 − 3 rad { p , q , r } ± 16 . 558 · 10 − 3 rad / s The nonlinear equations ( 27 ) were discretized by Eu- ler approximation with sampling time 0 . 01 s. The ini- tial states ζ 0 are bounded by X 0 = { G 0 , 0 } , where G 0 = diag  2 , 2 , 2 , 1 , 1 , 1 , π 6 , π 6 , π 2 , π 12 , π 12 , π 12  . T o generate process mea- surements, the discrete-time dynamics were simulated with ζ 0 = [0 . 5 0 1 0 1 × 5 π/ 3 0 1 × 3 ] T ∈ X 0 and process and mea- surement noises dra wn from uniform distrib utions. Figure 5 shows the trajectory performed by the quadrotor U A V along with the interv al hulls 11 of the enclosures computed by the meth- ods CZMV and ZMV, projected onto ( x , y , z )-axes. CZMV was implemented with h giv en by C2 , and since ( 27 ) is also a ffi ne in w k , d k , Theorem 3 was implemented with Z ⊇ µ ( h , W ) computed as described at the end of Remark 4 . The number of constraints and generators of the computed constrained zono- topes was limited to 40 and 12, respecti vely , while the number of generators of the computed zonotopes was limited to 40. The interval hulls of the constrained zonotopes obtained by CZMV were smaller than those from ZMV, demonstrating the accuracy of the proposed method. Figure 6 sho ws the radii of the constrained zonotopes and zonotopes computed by CZMV and ZMV, respecti vely . Both algorithms were capable of pro- viding tight bounds on the system states ζ k ∈ R 12 . Nevertheless, CZMV provided less conserv ati ve bounds than ZMV, e ven for a high-order nonlinear dynamical system such as ( 27 ) (the CZMV - to-ZMV ARR was 74 . 41%). Finally , Figure 7 compares the radii of the update sets computed by ZFO and CZFO with h giv en by C3 12 . Once again, CZFO provided less conserv ativ e bounds than ZFO (the CZFO-to-ZFO ARR was 74 . 45%). The results 10 In this experiment, the control action is computed using the real states ζ k . The approach in [ 11 ] can be used for feedback connection using a point that belongs to ˆ X k . 11 The con version from CG-rep to H-rep (see [ 25 ]) for the purposes of exact drawing is intractable for the constrained zonotopes in this example. 12 The increased complexity of the constrained zonotopes provided by C4 proved to be intractable for this e xample. from CZMV and CZFO were again very similar , with CZMV providing mar ginally better results (the CZMV -to-CZFO ARR was 99 . 93%). The ARR for di ff erent numbers of constraints are shown in T able 5 , with the average computed considering in addition simulations with di ff erent numbers of generators. Exe- cution times are shown in T able 6 . Note that most of the CZFO ex ecution times are smaller than the ones presented in T able 3 . This might be counter-intuiti ve since the current example has more state variables. Howe ver the use of C4 in Example 1 results in a relati vely more complex enclosure and therefore requires a much higher e xecution time for generator reduction and constraint elimination. Note that in this example, the com- putational times of the state estimators were greater than the considered sampling time of 0 . 01 s. Nevertheless, this fact does not in validate the obtained results since these were run in a nu- merical simulation. Better times can be achiev ed by optimized implementation of the algorithms and using more po werful hard- ware, for instance. Besides, note that even though the current ex ecution times of CZMV and CZFO would in principle pre vent their use in fast applications, the improved accurac y can signifi- cantly reduce the number of time steps required for guaranteed fault detection and isolation for systems in which execution time is not critical. T able 5: A verage radius ratio of the estimators with varying numbers of con- straints. Each average is tak en o ver 3 separate simulations with n g ∈ { 25 , 30 , 40 } . n c CZMV / ZMV CZFO / ZFO CZFO / CZMV 3 82 . 5% 82 . 2% 99 . 9% 6 76 . 7% 77 . 1% 100 . 7% 12 75 . 0% 74 . 8% 99 . 7% T able 6: A verage total times per iteration of the estimators with v arying numbers of constraints. Each average is taken over 15 separate simulations using n g ∈ { 25 , 30 , 40 } . T imes spent only on complexity reduction are shown in parenthesis. n c ZMV ZFO CZMV CZFO 0 53 . 3 (1 . 0) ms 174 . 2 (54 . 2) ms – – 3 – – 104 . 7 (7 . 3) ms 266 . 9 (128 . 6) ms 6 – – 109 . 5 (8 . 4) ms 615 . 0 (471 . 6) ms 12 – – 127 . 8 (12 . 3) ms 2 . 62 (2 . 46) s 6. Conclusions and future w ork This paper proposed two no vel approaches for set-valued state estimation of nonlinear discrete-time systems with unkno wn-but- bounded process and measurement uncertainties. A mean value extension and a first-order T aylor extension were developed based on constrained zonotopes, a generalization of zonotopes capable of describing strongly asymmetric conv ex sets. In ad- dition, measurement data were e ff ecti vely taken into account by means of generalized intersection, which resulted in far less conservati ve results than existing methods based on zonotopes. The accuracy of the proposed methods was demonstrated by means of three numerical examples, the third one being an ex- periment with a quadrotor U A V , considering a realistic scenario with uncertain measurements provided by sensors located on the aircraft. In the latter, e xecution times were longer than the 12 x (m) y (m) z (m) -0.4 -0.2 0 0.2 0.4 0.6 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.5 1 1.5 2 2.5 x (m) y (m) -0.4 -0.2 0 0.2 0.4 0.6 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 Figure 5: The trajectory performed by the quadrotor U A V (solid line) and the interval hulls of the constrained zonotopes (red boxes) and zonotopes (green boxes) estimated by CZMV and ZMV , respectively , projected onto ( x , y , z ). k Radius 0 500 1000 1500 0 0.5 1 1.5 2 Figure 6: Radii of the update sets computed by CZMV ( − ) and ZMV ( −− ) for the quadrotor U A V experiment. considered sampling time. Ne vertheless, an optimized imple- mentation of the methods, as well as more powerful hardware, implementation in C ++ and parallelization techniques, could be used to achie ve better times. This issue is left as a future work seeking the practical implementation in a real aircraft. k Radius 0 500 1000 1500 0 0.5 1 1.5 2 Figure 7: Radii of the update sets computed by CZFO ( − ) and ZFO ( −− ) for the quadrotor U A V experiment. Appendix. Computational complexity details T able 7 sho ws the computational comple xity of the basic operations used in the zonotope and constrained zonotope meth- ods. These comple xities assume generic inputs with dimensions R ∈ R n r × n in ( 4 ) ; Y = { G y , c y , A y , b y } in ( 6 ) , with G y ∈ R n r × n g r , c y ∈ R n r , A y ∈ R n c r × n g r , and b y ∈ R n c r ; J ∈ IR n s × n in Theo- rem 1 ; X = { G x , c x } or X = { G x , c x , A x , b x } , with G x ∈ R n × n g , c x ∈ R n , A x ∈ R n c × n g , and b x ∈ R n c ; k g and k c are the number of generators and constraints removed in the order reduction process, respecti vely . ‘Set inclusion’ refers to the zonotope in- clusion in [ 14 ] for zonotopes and the CZ-inclusion (Theorem 1 ) for constrained zonotopes. ‘Closest point’ and ‘Change center’ correspond to Propositions 1 and 2 , respectiv ely , which are LPs. For the latter , the bounds ˜ ξ L , ˜ ξ U are obtained using Algorithm 1 in [ 25 ]. Note that the interv al hull of zonotopes does not require the solution of LPs (see Remark 3 in [ 36 ]). In addition, we consider that each LP is solved at least with the performance of the simplex method presented in [ 46 ], which is O ( N d N 3 c ) with N d and N c the number of decision variables and constraints, respectiv ely . Note that these numbers can be inferred for each respectiv e LP directly from T able 7 . F or a detailed deriv ation of the computational complexities in T ables 1 and 7 , please see the supplementary material. T able 7: Computational complexity O ( · ) of basic operations. Operation Zonotopes Constrained zonotopes Linear mapping nn g n r nn g n r Minkowski sum n n Generalized intersection – nn g n r + n r n g r Interval hull nn g nn g ( n g + n c ) 3 Set inclusion nn g nn s n g + n c ( n g + n c ) 3 + nn 2 g n c Closest point – ( n + n g )( n + n g + n c ) 3 Change center – n g ( n + n g ) 3 + n 2 g n c Generator reduction n 2 n g + k g nn g ( n + n c ) 2 n g + k g ( n + n c ) n g Constraint elimination – k c ( n g + n c ) 3 + k c nn 2 g References [1] L. Jaulin, A nonlinear set membership approach for the localization and map building of underw ater robotics, IEEE T ransactions on Robotics 25 (1) (2009) 88–98. [2] S. Saeedi, M. T rentini, M. Seto, H. Li, Multiple-robot simultaneous local- ization and mapping: A revie w , Journal of Field Robotics 33 (1) (2016) 3–46. [3] L. Jaulin, Robust set-membership state estimation; application to underwa- ter robotics, Automatica 45 (1) (2009) 202–206. 13 [4] M. A. Santos, B. S. Rego, G. V . Ra ff o, A. Ferramosca, Suspended load path tracking control strategy using a tilt-rotor UA V, Journal of Advanced T ransportation 2017 (2017) 1–22. [5] F . A. Goodarzi, T . Lee, Global formulation of an e xtended Kalman filter on SE(3) for geometric control of a quadrotor U A V, Journal of Intelligent and Robotic Systems 88 (2) (2017) 395–413. [6] C. Combastel, Q. Zhang, A. Lalami, F ault diagnosis based on the enclosure of parameters estimated with an adaptiv e observer , in: Proc. of the 17th IF A C W orld Congress, 2008, pp. 7314–7319. [7] S. Zhao, B. Huang, F . Liu, Fault detection and diagnosis of multiple-model systems with mismodeled transition probabilities, IEEE T ransactions on Industrial Electronics 62 (8) (2015) 5063–5071. [8] D. M. Raimondo, G. R. Marseglia, R. D. Braatz, J. K. Scott, Closed-loop input design for guaranteed fault diagnosis using set-valued observers, Automatica 74 (2016) 107–117. [9] L. Chisci, A. Garulli, G. Zappa, Recursiv e state bounding by parallelotopes, Automatica 32 (7) (1996) 1049–1055. [10] V . T . H. Le, C. Stoica, T . Alamo, E. F . Camacho, D. Dumur, Zonotopic guaranteed state estimation for uncertain systems, Automatica 49 (11) (2013) 3418–3424. [11] B. S. Rego, D. M. Raimondo, G. V . Ra ff o, Path tracking control with state estimation based on constrained zonotopes for aerial load transportation, in: Proc. of the 57th IEEE Conference on Decision and Control, 2018, pp. 1979–1984. [12] B. S. Rego, G. V . Ra ff o, Suspended load path tracking control using a tilt-rotor U A V based on zonotopic state estimation, Journal of the Franklin Institute 356 (4) (2019) 1695–1729. [13] C. Combastel, Merging Kalman filtering and zonotopic state bounding for robust fault detection under noisy en vironment, in: Proc. of the 9th IF AC Symposium on Fault Detection, Supervision and Safety for T echnical Processes, 2015, pp. 289–295. [14] T . Alamo, J. Brav o, E. Camacho, Guaranteed state estimation by zonotopes, Automatica 41 (6) (2005) 1035–1043. [15] A. Girard, C. L. Guernic, E ffi cient reachability analysis for linear systems using support functions, in: Proc. of the 17th IF AC W orld Congress, 2008, pp. 8966–8971. [16] E. W alter , H. Piet-Lahanier , Exact recursive polyhedral description of the feasible parameter set for bounded-error models, IEEE Transactions on Automatic Control 34 (8) (1989) 911–915. [17] J. S. Shamma, K.-Y . T u, Set-valued observers and optimal disturbance rejection, IEEE T ransactions on Automatic Control 44 (2) (1999) 253–264. [18] F . Schweppe, Recursive state estimation: Unknown but bounded errors and system inputs, IEEE Transactions on Automatic Control 13 (1) (1968) 22–28. [19] C. Durieu, E. W alter , B. Polyak, Multi-input multi-output ellipsoidal state bounding, Journal of Optimization Theory and Applications 111 (2) (2001) 273–303. [20] B. T . Polyak, S. A. Nazin, C. Durieu, E. W alter, Ellipsoidal parameter or state estimation under model uncertainty , Automatica 40 (7) (2004) 1171–1179. [21] A. V icino, G. Zappa, Sequential approximation of feasible parameter sets for identification with set membership uncertainty , IEEE Transactions on Automatic Control 41 (6) (1996) 774–785. [22] C. Combastel, A state bounding observer based on zonotopes, in: Proc. of the 2003 European Control Conference, 2003, pp. 2589–2594. [23] S. B. Chabane, C. S. Maniu, T . Alamo, E. Camacho, D. Dumur, Improv ed set-membership estimation approach based on zonotopes and ellipsoids, in: Proc. of the 2014 European Control Conference, 2014, pp. 993–998. [24] M. Altho ff , B. H. Krogh, Zonotope bundles for the e ffi cient computation of reachable sets, in: Proc. of the 50th IEEE Conference on Decision and Control, 2011, pp. 6814–6821. [25] J. K. Scott, D. M. Raimondo, G. R. Marseglia, R. D. Braatz, Constrained zonotopes: a new tool for set-based estimation and fault detection, Auto- matica 69 (2016) 126–136. [26] C. Combastel, A state bounding observer for uncertain non-linear continuous-time systems based on zonotopes, in: Proc. of the 44th IEEE Conference on Decision and Control, and 2005 European Control Confer- ence, 2005, pp. 7228–7234. [27] J. W an, S. Sharma, R. Sutton, Guaranteed state estimation for nonlinear discrete-time systems via indirectly implemented polytopic set computa- tion, IEEE Transactions on Automatic Control 63 (12) (2018) 4317–4322. [28] R. E. Moore, R. B. Kearfott, M. J. Cloud, Introduction to Interv al Analysis, SIAM, Philadelphia, P A, USA, 2009. [29] M. Kie ff er , L. Jaulin, E. W alter, Guaranteed nonlinear state estimation using interval analysis, in: Proc. of the 37th Conference on Decision and Control, 1998, pp. 3966–3971. [30] L. Jaulin, Inner and outer set-membership state estimation, Reliable Com- puting 22 (2016) 47–55. [31] X. Y ang, J. K. Scott, Accurate set-based state estimation for nonlinear discrete-time systems using di ff erential inequalities with model redun- dancy , in: Proc. of the 57th IEEE Conference on Decision and Control, 2018, pp. 680–685. [32] J. S. Shamma, K.-Y . Tu, Approximate set-v alued observers for nonlinear systems, IEEE Transactions on Automatic Control 42 (5) (1997) 648–658. [33] T . Alamo, J. M. Bravo, M. J. Redondo, E. F . Camacho, A set-membership state estimation algorithm based on DC programming, Automatica 44 (1) (2008) 216–224. [34] B. S. Rego, D. M. Raimondo, G. V . Ra ff o, Set-based state estimation of nonlinear systems using constrained zonotopes and interval arithmetic, in: Proc. of the 2018 European Control Conference, 2018, pp. 1584–1589. [35] M. Altho ff , Reachability analysis of nonlinear systems using conservativ e polynomialization and non-conv ex sets, Hybrid Systems: Computation and Control (2013) 173–182. [36] W . K ¨ uhn, Rigorously computed orbits of dynamical systems without the wrapping e ff ect, Computing 61 (1) (1998) 47–67. [37] A. Platzer , E. M. Clarke, The image computation problem in hybrid sys- tems model checking, in: Proc. of the International W orkshop on Hybrid Systems: Computation and Control, 2007, pp. 473–486. [38] J. M. Bravo, T . Alamo, E. F . Camacho, Bounded error identification of systems with time-varying parameters, IEEE Transactions on Automatic Control 51 (7) (2006) 1144–1150. [39] S. Boyd, L. V andenberghe, Con vex Optimization, Cambridge, 2004. [40] A.-K. K opetzki, B. Sch ¨ urmann, M. Altho ff , Methods for order reduction of zonotopes, in: Proc. of the 56th Conference on Decision and Control, 2017, pp. 5626–5633. [41] X. Y ang, J. K. Scott, A comparison of zonotope order reduction techniques, Automatica 95 (2018) 378–384. [42] T . H. Cormen, C. E. Leiserson, R. L. Ri vest, C. Stein, Introduction to Algorithms, 3rd Edition, Cambridge, 2009. [43] W . Hagemann, E ffi cient geometric operations on conv ex polyhedra, with an application to reachability analysis of hybrid systems, Mathematics in Computer Science 9 (3) (2015) 283–325. [44] D. M. Raimondo, S. Riv erso, S. Summers, C. N. Jones, J. L ygeros, M. Morari, A set-theoretic method for v erifying feasibility of a f ast explicit nonlinear model predictiv e controller, in: R. Johansson, A. Rantzer (Eds.), Distributed Decision Making and Control, Springer London, London, 2012, pp. 289–311. [45] V . Mistler , A. Benallegue, N. K. M’Sirdi, Exact linearization and noninter- acting control of a 4 rotors helicopter via dynamic feedback, in: Proc. of the 10th IEEE International W orkshop on Robot and Human Interacti ve Communication, 2001, pp. 586–593. [46] J. A. Kelner , D. A. Spielman, A randomized polynomial-time simplex algorithm for linear programming, in: Proc. of the 38th Annual A CM Symposium on Theory of Computing, 2006, pp. 51–60. 14 Supplementary material  Brenner S. Rego a, ∗ , Guilherme V . Ra ff o a,b , Joseph K. Scott c , Davide M. Raimondo d a Graduate Progr am in Electrical Engineering, F ederal University of Minas Gerais, Belo Horizonte, MG 31270-901, Brazil b Department of Electronics Engineering, F ederal University of Minas Gerais, Belo Horizonte, MG 31270-901, Brazil c Department of Chemical and Biomolecular Engineering, Clemson University , Clemson, SC, USA d Department of Electrical, Computer and Biomedical Engineering, University of P avia, Italy Abstract This supplement details the deri v ation of the computational complexities of basic operations on zonotopes and constrained zonotopes, and the set-valued state estimators presented in the main paper . K e ywor ds: Nonlinear state estimation, Set-based computing, Reachability analysis, Con ve x polytopes 1. Computational complexity details: basic operations 1.1. Zonotopes 1.1.1. Linear image Let Z = { G z , c z } ⊂ R n , R ∈ R n r × n , with c z ∈ R n , G z ∈ R n × n g , and consider the linear image R Z = { RG z , Rc z } . Then, O ( R Z ) = O ( Rc z ) + O ( RG z ) = O ( nn r ) + O ( nn g n r ) = O ( nn g n r ) . 1.1.2. Minkowski sum Let Z = { G z , c z } ⊂ R n , W = { G w , c w } ⊂ R n , with c z ∈ R n , G z ∈ R n × n g , c w ∈ R n , G z ∈ R n × n g w , and consider the Minko wski sum Z ⊕ W = n h G z G w i , c z + c w o . Then, O ( Z ⊕ W ) = O ( c z + c w ) = O ( n ) . 1.1.3. Generator reduction Let Z = { G , c } ⊂ R n , with c ∈ R n , G ∈ R n × n g . Let k g ≤ ( n g − n ) denote the number of generators to be eliminated. The computational complexity of generator reduction of zonotopes is the same presented in [ 1 , 2 ]: O (eliminate k g generators from Z ) = O ( n 2 n g + k g n g n ) .  This work was supported in part by the project INCT under the grant CNPq 465755 / 2014-3, F APESP 2014 / 50851-0, and also by the Brazilian agencies CAPES, and F APEMIG. The material in this paper was not presented at any conference. ∗ Corresponding author Email addresses: brennersr7@ufmg.br (Brenner S. Rego), raffo@ufmg.br (Guilherme V . Ra ff o), jks9@clemson.edu (Joseph K. Scott), davide.raimondo@unipv.it (Davide M. Raimondo) Supplementary material August 29, 2019 1.1.4. Interval hull Let Z = { G , c } ⊂ R n , with c ∈ R n , G ∈ R n × n g . Consider the interval hull of Z gi ven by  c − ζ , c + ζ  , where [ 3 ] ζ i = n g X j = 1 | G i j | . Then, O (interval hull of Z ) = O ( c − ζ ) + O ( c + ζ ) + nO  n g X j = 1 | G i j |  = O ( c + ζ ) + nO  n g X j = 1 | G i j |  = O ( n ) + nO ( n g ) = O ( nn g ) . 1.1.5. Zonotope inclusion Let Z = p ⊕ M B n g ∞ , with p ∈ R n , M ∈ IR n × n g . Consider the zonotope inclusion  ( Z ) = p ⊕ [ mid ( M ) P ] B n g + n ∞ , where P ii = P n g j = 1 1 2 diam( M i j ), P i j = 0 for i , j [ 4 ]. Then, O (  ( Z )) = O (mid( M )) + O ( P ) = O (mid( M )) + nO ( P ii ) = O ( nn g ) + nO ( n g ) = O ( nn g ) . 1.2. Constrained zonotopes 1.2.1. Linear image Let Z = { G z , c z , A z , b z } ⊂ R n , R ∈ R n r × n , with c z ∈ R n , G z ∈ R n × n g , A z ∈ R n c × n g , b z ∈ R n c , and consider the linear image R Z = { RG z , Rc z , A z , b z } . Then, O ( R Z ) = O ( Rc z ) + O ( RG z ) = O ( nn r ) + O ( nn g n r ) = O ( nn g n r ) . 1.2.2. Minkowski sum Let Z = { G z , c z , A z , b z } ⊂ R n , W = { G w , c w , A w , b w } ⊂ R n , with c z ∈ R n , G z ∈ R n × n g , A z ∈ R n c × n g , b z ∈ R n c , c w ∈ R n , G w ∈ R n × n g w , A w ∈ R n c w × n g w , b w ∈ R n c w , and consider the Minko wski sum Z ⊕ W = ( h G z G w i , c z + c w , " A z 0 0 A w # , " b z b w #) . Then, O ( Z ⊕ W ) = O ( c z + c w ) = O ( n ) . 1.2.3. Generalized intersection Let Z = { G z , c z , A z , b z } ⊂ R n , Y = { G y , c y , A y , b y } , R ∈ R n r × n , with G z ∈ R n × n g , c z ∈ R n , A z ∈ R n c × n g , b z ∈ R n c , G y ∈ R n r × n g r , c y ∈ R n r , A y ∈ R n c r × n g r , b y ∈ R n c r . Consider the generalized intersection Z ∩ R Y =          h G z 0 i , c z ,           A z 0 0 A y RG z − G y           ,           b z b y c y − Rc z                    . Then, O ( Z ∩ R Y ) = O ( RG z ) + O ( − G y ) + O ( c y − Rc z ) = O ( nn g n r ) + O ( n r n g r ) + O ( n r + nn r ) = O ( nn g n r + n r n g r ) . 2 1.2.4. Generator reduction Let Z = { G , c , A , b } ⊂ R n , with c ∈ R n , G ∈ R n × n g , A ∈ R n c × n g , b ∈ R n c . Let k g ≤ n g − ( n + n c ) denote the number of generators to be eliminated. The computational complexity of generator reduction of constrained zonotopes is the same presented in [ 1 , 2 ], with reduction operated ov er the lifted zonotope Z + = (" G A # , " c − b #) . Therefore, O (eliminate k g generators from Z ) = O (( n + n c ) 2 n g + k g n g ( n + n c )) . 1.2.5. Constraint elimination Let Z = { G , c , A , b } ⊂ R n , with c ∈ R n , G ∈ R n × n g , A ∈ R n c × n g , b ∈ R n c . Let k c ≤ n c denote the number of constraints to be eliminated from Z . Follo wing the constraint elimination algorithm in [ 1 ], for each eliminated constraint the remaining constraints are first preconditioned through Gauss-Jordan elimination with full pi voting and then subjected to a rescaling procedure before the next constraint is eliminated. Therefore, O (eliminate 1 constraint from Z ) = O (pre-conditioning) + O (rescaling) + O (eliminate the constraint) . From the Appendix in [ 1 ], we hav e that O (pre-conditioning) + O (rescaling) = O ( n 2 c n g + n c n 2 g ) . T o eliminate the constraint, one of the generators of Z is first chosen by minimizing an approximation of the Hausdor ff distance between Z and the set obtained by constraint elimination. This requires solving (A.8) in [ 1 ] for each j ∈ { 1 , 2 , . . . , n g } , which has the reported complexity O (( n g + n c ) 3 ), and computing the Hausdor ff distance H ∗ j = k Gd ∗ j k 2 2 + k d ∗ j k 2 2 , with d ∗ j the solution of (A.8) for each j . This is done with the following comple xity: O (eliminate the constraint) = O (Hausdor ff distance minimization) + n g O ( H ∗ j ) = O (( n g + n c ) 3 ) + n g O ( nn g ) = O (( n g + n c ) 3 + nn 2 g ) . Then, O (eliminate 1 constraint from Z ) = O ( n 2 c n g + n c n 2 g ) + O (( n g + n c ) 3 + nn 2 g ) = O (( n g + n c ) 3 + nn 2 g ) . Finally , the complexity of eliminating k c constraints from Z can be bounded by O (eliminate k c constraints from Z ) = k c O (eliminate 1 constraint from Z ) = O ( k c ( n g + n c ) 3 + k c nn 2 g ) . 1.2.6. Interval hull (Pr operty 1) Let Z = { G , c , A , b } ⊂ R n , with c ∈ R n , G ∈ R n × n g , A ∈ R n c × n c , b ∈ R n c . Consider the interval hull of Z computed as in Property 1. In the standard form, each LP in Property 1 has N d = n g + 1 decision variables and N c = 2 n g + n c + 1 constraints. Then, O (interval hull of Z ) = 2 nO ( N d N 3 c ) = O ( nN d N 3 c ) = O ( n (1 + n g )(1 + 2 n g + n c ) 3 ) = O ( nn g ( n g + n c ) 3 ) . 1.2.7. CZ-inclusion (Theor em 1) Let X = p ⊕ M B ∞ ( A , b ) ⊂ R m , J ∈ IR n × m , ¯ X = ¯ p ⊕ ¯ M B ¯ n g ∞ ⊇ X , m ⊇ ( J − mid ( J )) ¯ p ,  ( J , X ) = mid ( J ) X ⊕ P B n ∞ , with p ∈ R m , M ∈ R m × n g , A ∈ R n c × n g , b ∈ R n c , P ii = 1 2 diam ( m i ) + 1 2 P ¯ n g j = 1 P m k = 1 diam ( J ik ) | ¯ M k j | , P i j = 0 for i , j . Assume that ¯ X is obtained by eliminating n c constraints from X , and m is computed using interval arithmetic. Then, O (  ( J , X )) = O (mid( J ) X ⊕ P B n ∞ ) + O (( J − mid( J )) ¯ p ) + nO  (1 / 2)diam( m i ) + 1 2 ¯ n g X j = 1 m X k = 1 diam( J ik ) | ¯ M k j |  3 + O (eliminate n c constraints from X ) = O ( nm + nmn g + n ) + O ( nm + nm ) + nO ( ¯ n g m ) + O ( n c ( n g + n c ) 3 + n c mn 2 g ) = O ( nmn g ) + O ( nm ) + O ( nmn g ) + O ( n c ( n g + n c ) 3 + mn c n 2 g ) = O ( nmn g + n c ( n g + n c ) 3 + mn 2 g n c ) . 1.2.8. Closest point (Pr oposition 1) Let Z = { G , c , A , b } ⊂ R n , h ∈ R n , with c ∈ R n , G ∈ R n × n g , A ∈ R n c × n g , b ∈ R n c . In the standard form, the LP in Proposition 1 has N d = n + n g + 1 decision variables and N c = 2 n + 2 n g + n c + 1 constraints. Then, O (Proposition 1) = O ( c + G ξ ∗ ) + O ( N d N 3 c ) = O ( n + nn g ) + O (( n + n g + 1)(2 n + 2 n g + n c + 1) 3 ) = O (( n + n g )( n + n g + n c ) 3 ) . 1.2.9. Change center (Pr oposition 2) Let Z = { G , c , A , b } =          h GE r 0 i , h ,           AE r 0 0 A GE r − G           ,           b − A ξ m b c − h                    ⊂ R n , with c ∈ R n , G ∈ R n × n g , A ∈ R n c × n g , b ∈ R n c , h ∈ R n , E r = 1 2 diag ( ξ U − ξ L ), ξ m = 1 2 ( ξ U + ξ L ). Consider ( ξ L , ξ U ) computed by solving the LP in Proposition 2, in which ( ˜ ξ L , ˜ ξ U ) is obtained by Algorithm 1 in the Appendix in [ 1 ]. In the standard form, the LP in Proposition 2 has N d = 3 n g decision variables and N c = n + 4 n g constraints. Then, O (Proposition 2) = O ( GE r ) + O ( AE r ) + O ( − G ) + O ( b − A ξ m ) + O ( c − h ) + O ((1 / 2)diag( ξ U − ξ L )) + O ((1 / 2)( ξ U + ξ L )) + O ( N d N 3 c ) + O (Algorithm 1) = O ( nn 2 g ) + O ( n c n 2 g ) + O ( nn g ) + O ( n c ) + O ( n c n g ) + O ( n ) + O ( n g ) + O ( n g ) + O ( N d N 3 c ) + O (Algorithm 1) = O ( nn 2 g ) + O ( n c n 2 g ) + O ( N d N 3 c ) + O (Algorithm 1) = O ( nn 2 g ) + O ( n c n 2 g ) + O ((3 n g )( n + 4 n g ) 3 ) + O (Algorithm 1) = O ( nn 2 g ) + O ( n c n 2 g ) + O ( n g ( n + n g ) 3 ) + O (Algorithm 1) = O ( n g ( n + n g ) 3 + n 2 g n c ) + O ( n 2 g n c ) = O ( n g ( n + n g ) 3 + n 2 g n c ) . 2. Computational complexity details: ZMV 2.1. Pr ediction step Let µ : R n × R n w → R n be continuously di ff erentiable and let ∇ x µ denote the gradient of µ with respect to its first argument. Let ˆ X k − 1 = { G x , c x } ⊂ R n , and W = { G w , c w } ⊂ R n w with c x ∈ R n , G x ∈ R n × n g , c w ∈ R n w , and G w ∈ R n w × n g w . Let Z be a zonotope such that µ ( c x , W ) ⊆ Z and let M ∈ IR n × n g be an interval matrix satisfying ∇ T x µ ( ˆ X k − 1 , W ) G x ⊆ M . According to Theorem 4 in [ 4 ], µ ( ˆ X k − 1 , W ) ⊆ ¯ X k = Z ⊕  ( M B n g ∞ ) , in which ¯ X k has n g + n g w + 2 n generators if Z is also computed by the mean value e xtension. Therefore, O (ZMV prediction ) = O ( Z ⊕  ( M B n g ∞ )) + O (  ( M B n g ∞ )) + O ( Z ) + O (interval hull of ˆ X k − 1 , W ) = O ( n ) + O ( nn g + n 2 n g ) + O ( n w n g w + nn w n g w ) + O ( nn g + n w n g w ) = O ( n 2 n g + nn w n g w ) . 4 2.2. Update step Let y k = Cx k + D u u k + D v v k , with y k ∈ R n y , x k ∈ R n , u k ∈ R n u , v k ∈ R n v . Consider the zonotopes ¯ X k = { G x , c x } ⊂ R n , V = { G v , c v } ⊂ R n v with c x ∈ R n , G x ∈ R n × ¯ n g , c v ∈ R n v , and G v ∈ R n v × n g v . For simplicity , define m g = n g + n g w . Follo wing the intersection method proposed in [ 5 ], the update step is computed iteratively using one row of y k at a time, as described below . Algorithm: Update step. (1) Assign ˜ X k = { ˜ G x , ˜ c x } ← ¯ X k , i ← 1. (2) Assign ˜ y ← i -th row of y k , p T ← i -th ro w of C , d T u ← i -th ro w of D u , d T v ← i -th ro w of D v . (3) Compute the strip S ← { x : | p T x − d | ≤ σ } , where d = ˜ y − d T u u k − d v c v , σ = P n g v  = 1 P n v j = 1 | ( D v ) i j ( G v ) j  | . (4) Compute the tight strip ˜ S ← { x : | p T x − ˜ d | ≤ ˜ σ } , where ˜ σ = 1 2 ( ¯ σ U − ¯ σ L ) , ˜ d = d + 1 2 ( ¯ σ U + ¯ σ L ) , with ¯ σ L = max {− σ , p T ˜ c x − k ˜ G T x p k 1 − d } , ¯ σ U = min { σ , p T ˜ c x + k ˜ G T x p k 1 − d } . (5) Assign ˜ X k ← ˜ X k ∩ ˜ S = { G ( j ∗ ) , c ( j ∗ ) } , where j ∗ = arg min j det( G ( j ) G ( j ) T ) , with j ∈ { 0 , 1 , . . . , ¯ n g } , c ( j ) =            ˜ c x +       ˜ d − p T ˜ c x p T ( ˜ G x ) : , j       ( ˜ G x ) : , j if 1 ≤ j ≤ ¯ n g and p T ( ˜ G x ) : , j , 0 ˜ c x otherwise G ( j ) =        [ g j 1 g j 2 . . . g j ¯ n g ] if 1 ≤ j ≤ ¯ n g and p T ( ˜ G x ) : , j , 0 ˜ G x otherwise g j  =                  ( ˜ G x ) : , −       p T ( ˜ G x ) : , p T ( ˜ G x ) : , j       ( ˜ G x ) : , j if  , j       ˜ σ p T ( ˜ G x ) : , j       ( ˜ G x ) : , j if  = j (6) If i < n y , go to Step 2 and assign i ← i + 1. Otherwise, assign ˆ X k ← ˜ X k . The computational complexity is O (ZMV update ) = n y O (update with the i -th row of y k ) = n y  O (Step (3)) + O (Step (4)) + O (Step (5))  . But, O (Step (3)) = O (1) + O ( n u ) + O ( n v ) + O ( n v n g v ) = O ( n u + n v n g v ) , O (Step (4)) = O (1) + O ( n ) + O ( n ¯ n g ) = O ( n ¯ n g ) . O (Step (5)) = ¯ n g O ( n 3 ) + ¯ n g O ( n 2 ¯ n g ) + ¯ n g O ( n ) + ¯ n g O ( n ¯ n g ) = O ( n 3 ¯ n g + n 2 ¯ n 2 g ) . Since ¯ n g = m g + 2 n , hence O (ZMV update ) = n y  O ( n u + n v n g v ) + O ( n ¯ n g ) + O ( n 3 ¯ n g + n 2 ¯ n 2 g + n ¯ n 2 g )  = n y  O ( n 3 ¯ n g + n 2 ¯ n 2 g + n u + n v n g v )  = O ( n y ( n 3 ( m g + n ) + n 2 ( m g + n ) 2 + n u + n v n g v )) . 5 2.3. Or der r eduction Let n g denote the desired number of generators of ˆ X k after order reduction, i.e., ˆ X k must hav e the same complexity as ˆ X k − 1 . Order reduction is performed by eliminating k g generators from ˆ X k . For simplicity , define m g = n g + n g w . The prediction and update steps of the ZMV lead to ˆ X k with m g + 2 n generators. Therefore, k g = n g w + 2 n . Consequently , O (ZMV reduction ) = O (eliminate k g generators from ˆ X k ) = O ( n 2 ( m g + n ) + k g ( m g + n ) n ) = O ( n 2 ( m g + n ) + ( n g w + n )( m g + n ) n ) = O ( n 2 ( m g + n ) + n ( n g w + n )( m g + n )) . 3. Computational complexity details: ZFO Let η : R m → R n be of class C 2 and let z = ( x k − 1 , w k − 1 ) ∈ R m denote its argument, where m = n + n w . Let Z = ˆ X k − 1 × W = { G , c } ⊂ R m be a zonotope with m g = n g + n g w generators. For each q = 1 , 2 , . . . , n , let Q [ q ] ∈ IR m × m and ˜ Q [ q ] ∈ IR m g × m g be interv al matrices satisfying Q [ q ] ⊇ H η q ( Z ) and ˜ Q [ q ] ⊇ G T Q [ q ] G . Moreover , define ˜ c , ˜ G , ˜ G d as in Theorem 3. Then, the method in [ 6 ] ensures that η ( Z ) ⊆ ¯ X k = η ( c ) ⊕ ∇ T η ( c )( G B m g ∞ ) ⊕ ˜ c ⊕ h ˜ G ˜ G d i B ˜ m g ∞ , where ˜ m g = (1 / 2) m 2 g + (3 / 2) m g + n , and therefore ¯ X k has (1 / 2) m 2 g + (5 / 2) m g + n generators. Thus, O (CZFO prediction ) = O ( η ( c ) ⊕ ∇ T η ( c )( G B m g ∞ ) ⊕ ˜ c ⊕ h ˜ G ˜ G d i B ˜ m g ∞ ) = O ( n ) + O ( ˜ c ) + O ( ˜ G ) + O ( ˜ G d ) + nO ( ˜ Q [ q ] ) + O (interval hull of Z ) = O ( n ) + O ( nm g ) + O ( nm 2 g ) + O ( nm 2 g ) + nO ( m 2 m g + mm 2 g ) + O ( mm g ) = O ( n ( m 2 m g + mm 2 g )) . 3.1. Update step Let y k = Cx k + D u u k + D v v k , with y k ∈ R n y , x k ∈ R n , u k ∈ R n u , and v k ∈ R n v . Consider the zonotopes ¯ X k = { G x , c x } ⊂ R n and V = { G v , c v } ⊂ R n v with c x ∈ R n , G x ∈ R n × ¯ n g , c v ∈ R n v , and G v ∈ R n v × n g v . For simplicity , define m g = n g + n g w . The ZFO also follows the intersection method proposed in [ 5 ]. Therefore, the update step is equiv alent to the ZMV update, but with ¯ n g = (1 / 2) m 2 g + (5 / 2) m g + n . Hence O (ZFO update ) = n y  O ( n 3 ¯ n g + n 2 ¯ n 2 g + n u + n v n g v )  = O ( n y ( n 3 ( m 2 g + n ) + n 2 ( m 2 g + n ) 2 + n u + n v n g v )) . 3.2. Or der r eduction Let n g denote the desired number of generators of ˆ X k after order reduction, i.e., ˆ X k must hav e the same complexity as ˆ X k − 1 . Order reduction is performed by eliminating k g generators from ˆ X k . For simplicity , define m g = n g + n g w . The prediction and update steps of ZFO lead to ˆ X k with (1 / 2) m 2 g + (5 / 2) m g + n generators. Therefore, k g = (1 / 2) m 2 g + (5 / 2) m g + n − n g . Consequently , O (ZFO reduction ) = O (eliminate k g generators from ˆ X k ) = O ( n 2 ( m 2 g + m g + n ) + k g ( m 2 g + m g + n ) n ) = O ( n 2 ( m 2 g + n ) + ( m 2 g + n )( m 2 g + n ) n ) = O ( n 2 ( m 2 g + n ) + n ( m 2 g + n ) 2 ) . 6 4. Computational complexity details: CZMV 4.1. Pr ediction step (Theorem 2) Let µ : R n × R n w → R n be continuously di ff erentiable and let ∇ x µ denote the gradient of µ with respect to its first argument. Let ˆ X k − 1 = { G x , c x , A x , b x } ⊂ R n and W = { G w , c w , A w , b w } ⊂ R n w with c x ∈ R n , G x ∈ R n × n g , A x ∈ R n c × n g , b x ∈ R n c , c w ∈ R n w , G w ∈ R n w × n g w , A w ∈ R n c w × n g w , and b w ∈ R n c w , and choose any h ∈ ˆ X k − 1 . Let Z be a constrained zonotope such that µ ( h , W ) ⊆ Z and let J ∈ IR n × n be an interval matrix satisfying ∇ T x µ ( ˆ X k − 1 , W ) ⊆ J . Theorem 2 ensures that µ ( ˆ X k − 1 , W ) ⊆ ¯ X k = Z ⊕  ( J , ˆ X k − 1 − h ) , where ¯ X k has n g + n g w + 2 n generators and n c + n c w constraints if Z is also computed by the mean value extension. Therefore, O (CZMV prediction ) = O ( Z ⊕  ( J , ˆ X k − 1 − h )) + O (  ( J , ˆ X k − 1 − h )) + O ( Z ) + O (interval hull of ˆ X k − 1 , W ) = O ( n ) + O ( n 2 n g + n c ( n g + n c ) 3 + nn 2 g n c ) + O ( nn w n g w + n c w ( n g w + n c w ) 3 + n w n 2 g w n c w ) + O ( nn g ( n g + n c ) 3 + n w n g w ( n g w + n c w ) 3 ) = O ( n 2 n g + nn w n g w + ( nn g + n c )( n g + n c ) 3 + ( n w n g w + n c w )( n g w + n c w ) 3 ) . 4.2. Update step Let ¯ X k = { G x , c x , A x , b x } ⊂ R n , V = { G v , c v , A v , b v } ⊂ R n w , u k ∈ R n u , and y k ∈ R n y with c x ∈ R n , G x ∈ R n × ¯ n g , A x ∈ R ¯ n c × ¯ n g , b x ∈ R ¯ n c , c v ∈ R n v , G v ∈ R n v × n g v , A v ∈ R n c v × n g v , and b w ∈ R n c v . The update step is performed according to ˆ X k = ¯ X k ∩ C (( y k − D u u k ) ⊕ ( − D v V )) , where C ∈ R n y × n , D u ∈ R n y × n u , and D v ∈ R n y × n v . As a result, from the prediction step of CZMV , ¯ X k has ¯ n g = n g + n g w + 2 n generators and ¯ n c = n c + n c w constraints. Therefore, O (CZMV update ) = O ( ¯ X k ∩ C (( y k − D u u k ) ⊕ ( − D v V ))) = O ( n ( n g + n g w + 2 n ) n y + n y n g v ) + O ( n y ) + O ( n y + n y n u ) + O ( n y n v n g v ) = O ( n y n ( m g + n ) + n y n u + n y n v n g v ) , where m g = n g + n g w . 4.3. Or der r eduction Let n g and n c denote the desired number of generators and constraints of ˆ X k after order reduction, i.e., ˆ X k must hav e the same complexity as ˆ X k − 1 . Order reduction is performed by eliminating k c constraints and then k g generators from ˆ X k . For simplicity , define m g = n g + n g w , m c = n c + n c w , δ n = n − n y , δ w = n g w − n c w , and δ v = n g v − n c v . The prediction and update steps of CZMV lead to ˆ X k with n g + n g w + 2 n + n g v generators and n c + n c w + n c v + n y constraints. Therefore, k c = n c w + n c v + n y , and k g = n g w + 2 n + n g v − k c = n g w + 2 n + n g v − n c w − n c v − n y . Consequently , O (eliminate k c constraints from ˆ X k ) = O ( k c ( m g + 2 n + n g v + m c + n c v + n y ) 3 + k c n ( m g + 2 n + n g v ) 2 ) = O (( n c w + n c v + n y )( m g + m c + n g v + n c v + n + n y ) 3 + ( n c w + n c v + n y ) n ( m g + n g v + n ) 2 ) = O (( n c w + n c v + n y )( m g + m c + n g v + n c v + n + n y ) 3 . Note that by eliminating k c constraints from ˆ X k , the same number of generators (i.e., k c ) is also remo ved [ 1 ]. Hence, O (elim. k g gen. from ˆ X k after const. elim.) = O (( n + n c ) 2 ( m g + 2 n + n g v − k c ) + k g ( m g + 2 n + n g v − k c )( n + n c )) 7 = O (( n + n c ) 2 ( δ n + δ w + δ v + n g ) + ( δ n + δ w + δ v )( δ n + δ w + δ v + n g )( n + n c )) . Finally , O (CZMV reduction ) = O (eliminate k c constraints from ˆ X k ) + O (elim. k g gen. from ˆ X k after const. elim.) = O (( n c w + n c v + n y )( m g + m c + n g v + n c v + n + n y ) 3 + ( n + n c ) 2 ( n g + δ n + δ w + δ v ) + ( n + n c )( δ n + δ w + δ v )( n g + δ n + δ w + δ v )) . 5. Computational complexity details: CZFO 5.1. Pr ediction step Let η : R m → R n be of class C 2 and let z = ( x k − 1 , w k − 1 ) ∈ R m denote its argument, where m = n + n w . Let Z = ˆ X k − 1 × W = { G , c , A , b } ⊂ R m be a constrained zonotope with m g = n g + n g w generators and m c = n c + n c w constraints. For each q = 1 , 2 , . . . , n , let Q [ q ] ∈ IR m × m and ˜ Q [ q ] ∈ IR m g × m g be interv al matrices satisfying Q [ q ] ⊇ H η q ( Z ) and ˜ Q [ q ] ⊇ G T Q [ q ] G . Moreover , define ˜ c , ˜ G , ˜ G d , ˜ A , ˜ b , and L as in Theorem 3. Choosing any h ∈ Z , Theorem 3 ensures that η ( Z ) ⊆ ¯ X k = η ( h ) ⊕ ∇ T η ( h )( Z − h ) ⊕ R , where R = ˜ c ⊕ h ˜ G ˜ G d i B ∞ ( ˜ A , ˜ b ) ⊕  ( L , ( c − h ) ⊕ 2 G B ∞ ( A , b )) and ¯ X k has (1 / 2) m 2 g + (5 / 2) m g + 2 n generators and (1 / 2) m 2 c + (5 / 2) m c constraints. Therefore, O (CZFO prediction ) = O ( η ( h ) ⊕ ∇ T η ( h )( Z − h ) ⊕ R ) + O ( R ) = O ( n ) + O ( nmm g ) + O ( R ) = O ( nmm g ) + O ( R ) . But O ( R ) = O ( ˜ c ⊕ h ˜ G ˜ G d i B ∞ ( ˜ A , ˜ b ) ⊕  ( L , ( c − h ) ⊕ 2 G B ∞ ( A , b ))) = O ( n ) + O ( ˜ c ) + O ( ˜ G ) + O ( ˜ G d ) + O ( ˜ A ) + O ( ˜ b ) + O ( L ) + O (  ( L , ( c − h ) ⊕ 2 G B ∞ ( A , b ))) + nO ( ˜ Q [ q ] ) + O (interval hull of Z ) = O ( n ) + O ( nm g ) + O ( nm 2 g ) + O ( nm 2 g ) + O ( m g m 2 c + m 2 g m 2 c ) + O ( m g m 2 c ) + O ( nm 2 ) + O ( n + nmm g + m c ( m g + m c ) 3 + mm 2 g m c ) + nO ( m 2 m g + mm 2 g ) + O ( mm g ( m g + m c ) 3 ) = O ( n ( m 2 m g + mm 2 g ) + ( mm g + m c )( m g + m c ) 3 ) . Therefore, O (CZFO prediction ) = O ( nmm g ) + O ( R ) = O ( n ( m 2 m g + mm 2 g ) + ( mm g + m c )( m g + m c ) 3 ) . 5.2. Update step Let ¯ X k = { G x , c x , A x , b x } ⊂ R n , V = { G v , c v , A v , b v } ⊂ R n w , u k ∈ R n u , and y k ∈ R n y with c x ∈ R n , G x ∈ R n × ¯ n g , A x ∈ R ¯ n c × ¯ n g , b x ∈ R ¯ n c , c v ∈ R n v , G v ∈ R n v × n g v , A v ∈ R n c v × n g v , and b w ∈ R n c v . The update step is performed according to ˆ X k = ¯ X k ∩ C (( y k − D u u k ) ⊕ ( − D v V )) , where C ∈ R n y × n , D u ∈ R n y × n u , and D v ∈ R n y × n v . As a result, from the prediction step of CZFO, ¯ X k has ¯ n g = (1 / 2) m 2 g + (5 / 2) m g + 2 n generators and ¯ n c = (1 / 2) m 2 c + (5 / 2) m c constraints, where m g = n g + n g w , and m c = n c + n c w . Therefore, O (CZFO update ) = O ( ¯ X k ∩ C (( y k − D u u k ) ⊕ ( − D v V ))) = O ( n ((1 / 2) m 2 g + (5 / 2) m g + 2 n ) n y + n y n g v ) + O ( n y ) + O ( n y + n y n u ) + O ( n y n v n g v ) = O ( n y n ( m 2 g + n ) + n y n u + n y n v n g v ) . 8 5.3. Or der r eduction Let n g and n c denote the desired number of generators and constraints of ˆ X k after order reduction, i.e., ˆ X k must hav e the same complexity as ˆ X k − 1 . Order reduction is performed by eliminating k c constraints and then k g generators from ˆ X k . For simplicity , define m g = n g + n g w , m c = n c + n c w , δ n = n − n y , δ v = n g v − n c v , and ˜ δ = m 2 g − m 2 c . The prediction and update steps of CZFO lead to ˆ X k with (1 / 2) m 2 g + (5 / 2) m g + 2 n + n g v generators and (1 / 2) m 2 c + (5 / 2) m c + n c v + n y constraints. Therefore, k c = (1 / 2) m 2 c + (5 / 2) m c + n c v + n y − n c and k g = (1 / 2) m 2 g + (5 / 2) m g + 2 n + n g v − n g − k c = (1 / 2) m 2 g + (5 / 2) m g + 2 n + n g v − n g − (1 / 2) m 2 c − (5 / 2) m c − n c v − n y + n c . Consequently , O (eliminate k c constraints from ˆ X k ) = O ( k c ( m 2 g + m g + n + n g v + m 2 c + m c + n c v + n y ) 3 + k c n ( m 2 g + m g + n ) 2 ) = O (( m 2 c + n c v + n y )( m 2 g + n + n g v + m 2 c + n c v + n y ) 3 + ( m 2 c + n c v + n y ) n ( m 2 g + n ) 2 ) = O (( m 2 c + n c v + n y )( m 2 g + m 2 c + n g v + n c v + n + n y ) 3 . Once again, by eliminating k c constraints from ˆ X k the same number of generators (i.e., k c ) is also remov ed [ 1 ]. Hence, O (elim. k g gen. from ˆ X k after const. elim.) = O (( n + n c ) 2 ( m 2 g + n + n g v − m 2 c − n c v − n y ) + k g ( m 2 g + n + n g v − m 2 c − n c v − n y )( n + n c )) = O (( n + n c ) 2 ( ˜ δ + δ n + δ v ) + ( ˜ δ + δ n + δ v )( ˜ δ + δ n + δ v )( n + n c )) = O (( n + n c ) 2 ( ˜ δ + δ n + δ v ) + ( ˜ δ + δ n + δ v ) 2 ( n + n c )) . Finally , O (CZFO reduction ) = O (eliminate k c constraints from ˆ X k ) + O (elim. k g gen. from ˆ X k after const. elim.) = O (( m 2 c + n c v + n y )( m 2 g + m 2 c + n g v + n c v + n + n y ) 3 + ( n + n c ) 2 ( ˜ δ + δ n + δ v ) + ( ˜ δ + δ n + δ v ) 2 ( n + n c )) . References [1] J. K. Scott, D. M. Raimondo, G. R. Marseglia, R. D. Braatz, Constrained zonotopes: a new tool for set-based estimation and f ault detection, Automatica 69 (2016) 126–136. [2] X. Y ang, J. K. Scott, A comparison of zonotope order reduction techniques, Automatica 95 (2018) 378–384. [3] W . K ¨ uhn, Rigorously computed orbits of dynamical systems without the wrapping e ff ect, Computing 61 (1) (1998) 47–67. [4] T . Alamo, J. Brav o, E. Camacho, Guaranteed state estimation by zonotopes, Automatica 41 (6) (2005) 1035–1043. [5] J. M. Brav o, T . Alamo, E. F . Camacho, Bounded error identification of systems with time-varying parameters, IEEE T ransactions on Automatic Control 51 (7) (2006) 1144–1150. [6] C. Combastel, A state bounding observer for uncertain non-linear continuous-time systems based on zonotopes, in: Proc. of the 44th IEEE Conference on Decision and Control, and 2005 European Control Conference, 2005, pp. 7228–7234. 9

Original Paper

Loading high-quality paper...

Comments & Academic Discussion

Loading comments...

Leave a Comment