Data-Driven Set-Based Estimation using Matrix Zonotopes with Set Containment Guarantees
Alexander Berndt, Amr Alanwar, Karl Henrik Johansson, Henrik Sandberg
DData-Driven Set-Based Estimation using Matrix Zonotopeswith Set Containment Guarantees
Alexander Berndt, Amr Alanwar, Karl Henrik Johansson, and Henrik Sandberg
Abstract — We propose a method to perform set-based stateestimation of an unknown dynamical system using a data-driven set propagation function. Our method comes with set-containment guarantees, making it applicable to the estimationof safety-critical systems. The method consists of two phases:(1) an offline learning phase where we collect noisy state-input data to determine a function to propagate the state-setahead in time; and (2) an online estimation phase consistingof a time update and a measurement update. It is assumedthat sets bound measurement noise and disturbances, but weassume no knowledge of their statistical properties. These setsare described using zonotopes, allowing efficient propagationand intersection operations. We propose two approaches toperform the measurement update. The method is extended toconstrained zonotopes. Simulations show that the proposed esti-mator yields state sets comparable in volume to the confidencebounds obtained by a Kalman filter approach, but with theaddition of state set-containment guarantees. We observe thatusing constrained zonotopes yields smaller sets, but with highercomputational cost compared to unconstrained zonotopes.
I. I
NTRODUCTION
Set-based estimation involves the computation of a time-updated set, which is guaranteed to contain the system’s truestate at each time step given bounded uncertainties. Existingset-based observers require a system model to propagatethe state set at each time step. We address the problem ofpropagating the state set using only noisy offline state andinput data and merging this with online measurements toobtain a time-varying state set which is guaranteed to containthe true system’s state at each time-step. This problem isessential in safety-critical applications [1].Two popular set-based estimators are interval observersand set-membership observers. Interval-based observers gen-erally generate state estimates by utilizing an observer gainto fuse a model-based time update of the state with currentmeasurements. For example, the authors in [2] proposean exponentially stable interval-based observer for time-invariant linear systems. Set-membership observers follow ageometrical approach by intersecting the state-space regionsconsistent with the model with those from the measurementsto obtain the current state set [3]. This approach has beenextended to sensor networks with event-based communi-cation in [4] and multi-rate systems in [5]. Various setrepresentations have been used for set-membership observers
The authors are with the Division of Decision and Control Sys-tems at KTH Royal Institute of Technology. {alberndt, alanwar, kallej,hsan}@kth.seThis work was supported by the Swedish Research Council, the Knutand Alice Wallenberg Foundation, and the Democritus project on Decision-making in Critical Societal Infrastructures by Digital Futures. such as ellipsoids [6], polytopes [7] and zonotopes [8].Zonotopes are a special class of polytopes for which onecan efficiently compute linear maps and Minkowski sums –both frequent operations performed by set-based observers.All the aforementioned observers use a model of theunderlying system to propagate the state set. However, iden-tifying a system model is often time-consuming, and theidentified model is not necessarily well-suited for estimationor control. Recent works based on Willems’ fundamentallemma [9] have shown that system trajectories can be useddirectly to synthesize controllers. The authors in [10] presentan extended Kalman filter and model predictive control(MPC) scheme computed directly from system trajectories.Stability and robustness guarantees for such a data-drivencontrol scheme are presented in [11], and for an MPC schemein [12]. An alternative approach is to find a set of modelsthat is consistent with data and use this set of models topropagate a state set [13]. Up to our knowledge, this paperpresents the first data-driven set-based observer.Our contribution is a novel method to perform set-based state estimation with set-containment guarantees givenbounded, noisy measurements and known inputs. Themethod consists of an offline learning phase to determinea state-propagation function f ( · ) directly from data, and an online estimation phase to perform a time update using f ( · ) and measurements iteratively to track the system state. Wepresent two approaches to perform the measurement updateutilizing either the singular value decomposition (SVD) ofthe observation matrix or an optimization formulation. Wecompare the approaches in simulation. Our method is shownto yield set-based state estimates similar in size to σ confi-dence bounds of an approach based on system identificationand a Kalman filter, but with the addition of set-containmentguarantees.The rest of this paper is outlined as follows. Sec. IIintroduces preliminary concepts and the problem statement.We present our method in Sec. III and evaluate it in Sec. IV.Finally, Sec. V concludes the paper.II. P RELIMINARIES AND P ROBLEM S TATEMENT
We introduce first zonotopes and matrix zonotopes.
Definition 1. (Zonotope [14])
Given a center c ∈ R n anda number ξ ∈ N of generator vectors in a generator matrix G = [ g (1) , ..., g ( ξ ) ] ∈ R n × ξ , a zonotope is a set Z = (cid:110) x ∈ R n (cid:12)(cid:12)(cid:12) x = c + ξ (cid:88) i =1 β i g ( i ) , − ≤ β i ≤ (cid:111) . (1) a r X i v : . [ ee ss . S Y ] J a n e use the shorthand notation Z = (cid:104) c, G (cid:105) . Given two zonotopes Z and Z , we use the notation + for the Minkowski sum, and Z − Z to denote Z + ( −Z ) . Definition 2. (Matrix zonotope [1, p.52])
Given a centermatrix C ∈ R n × k and ξ ∈ N generator matrices G ( i ) ∈ R n × k where i ∈ { , . . . , ξ } , a matrix zonotope is the set M = (cid:110) X ∈ R n × k (cid:12)(cid:12)(cid:12) X = C + ξ (cid:88) i =1 β i G ( i ) , − ≤ β i ≤ (cid:111) . We use the notation M = (cid:104) C, G (1: ξ ) (cid:105) , where G (1: ξ ) =[ G (1) , . . . , G ( ξ ) ] . We consider estimating the set of all possible system statesusing an array of q sensors. Our system is described as x ( k + 1) = Ax ( k ) + Bu ( k ) + w ( k ) , (2) y i ( k ) = C i x ( k ) + v i ( k ) , i ∈ { , . . . , q } , (3)where x ( k ) ∈ R n is the system state, u ( k ) ∈ R m the input, y i ( k ) ∈ R p i the measurement, x (0) ∈ X the initial condi-tion where X is the initial bounding zonotope. Furthermore,the system matrices A ∈ R n × n and B ∈ R n × m are unknown whereas C i ∈ R p i × n is known for all i ∈ { , . . . , q } .The noise w ( k ) ∈ Z w and v i ( k ) ∈ Z v,i are assumed tobelong to the bounding zonotopes Z v,i = (cid:104) c v,i , G v,i (cid:105) for i ∈ { , . . . , q } and Z w = (cid:104) c w , G w (cid:105) ⊂ R n , respectively. Wedenote the Frobenius norm by (cid:107) . (cid:107) F and the null space of amatrix A by ker ( A ) .Let R k denote a set containing x ( k ) given the exact system model and bounded, but unknown , process andmeasurement noise. The problem addressed in this paperis to develop an algorithm that returns a set ˆ R k ⊃ R k ,which is guaranteed to contain the true state x ( k ) at eachtime instance k , i.e., x ( k ) ∈ ˆ R k for all k , given modeluncertainties and measurement noise.III. D ATA - DRIVEN S ET - BASED E STIMATION
Our proposed data-driven set estimator consists of twophases: an offline learning phase and an online estimationphase . The online phase consists of iteratively performinga time update and a measurement update. In the offlinephase, we compute the function to perform the time update. ˜ R k ⊂ R n and ˆ R k ⊂ R n are the time and measurementupdated sets at k , respectively. A. Offline Learning Phase
The objective of this phase is to compute a function f : R n × R m → R n , such that ˜ R k +1 = f ( ˆ R k , U k ) , i.e., f returns ˜ R k +1 given a known input zonotope U k and themeasurement updated set ˆ R k at time-step k such that we canguarantee x ( k + 1) ∈ ˜ R k +1 for all k . During this phase, weassume that we have access to an input sequence u ( k ) andnoisy state measurements z ( k ) such that z ( k ) = x ( k ) + γ ( k ) , (4)where the noise γ ( k ) is bounded by the zonotope Z γ = (cid:104) c γ , G γ (cid:105) , i.e., γ ( k ) ∈ Z γ . Given an experiment yielding a sequence of noisy data of length T , we can construct thefollowing sequences Z T = (cid:2) z (1) . . . z ( T ) (cid:3) ,Z T − = (cid:2) z (0) . . . z ( T − (cid:3) ,U T − = (cid:2) u (0) . . . u ( T − (cid:3) . (5)Furthermore, we denote the sequence of unknown processnoise w ( k ) as W T − := (cid:2) w (0) . . . w ( T − (cid:3) . Here, W T − ∈ M w where M w = (cid:104) C M w , G (1: αT ) M w (cid:105) is the matrixzonotope resulting from the concatenation of noise zonotopes Z w = (cid:104) c Z w , [ g (1) Z w , . . . , g ( α ) Z w ] (cid:105) as C M w = (cid:2) c Z w . . . c Z w (cid:3) ,G (1+( i − T ) M w = (cid:104) g ( i ) Z w n × ( T − (cid:105) ,G ( j +( i − T ) M w = (cid:104) n × ( j − g ( i ) Z w n × ( T − j ) (cid:105) ,G ( T +( i − T ) M w = (cid:104) n × ( T − g ( i ) Z w (cid:105) , for all i = { , . . . , α } , j = { , . . . , T − } [13]. In a similarfashion, we describe the matrix zonotope of γ ( k ) as M γ = (cid:104) C M γ , G (1: αT ) M γ (cid:105) . We substitute the sequences (5) into (4) and(2) and bound the noises by their matrix zonotopes yielding Z T − M γ = A (cid:0) Z T − − M γ (cid:1) + BU T − + M w . (6)Assuming a bound on the noise term Aγ ( k ) , we can write Aγ ( k ) ∈ Z Aγ , which we combine into a matrix zonotope M Aγ as above. With this, we rearrange (6) to obtain Z T − M γ + M Aγ − M w = (cid:2) A B (cid:3) (cid:20) Z T − U T − (cid:21) . We assume that the data matrix D := (cid:2) Z (cid:62) T − U (cid:62) T − (cid:3) (cid:62) has full row rank n + m , implying that the input sequencesufficiently excites the system [13]. As a result, the set M Σ is guaranteed to contain the true system matrices, i.e., (cid:2) A B (cid:3) ∈ M Σ = (cid:0) Z T − M γ + M Aγ − M w (cid:1) (cid:20) Z T − U T − (cid:21) † , where D † denotes the Moore-Penrose pseudo-inverse of D ,which is a right inverse. Hence, ˜ R k +1 = M Σ (cid:0) ˆ R k × U k (cid:1) + Z w , (7)which defines f ( · ) introduced above. Remark 1.
Note that in our computation of f ( · ) we assumea known bound on Aγ ( k ) . Dealing with measurement noisewithout this assumption remains an open problem in thefield of data-driven control. Notable other works such as[11] require a similar assumption to derive controllers inthe noisy-measurement setting.B. Online Estimation Phase using Zonotopes In this subsection, we present the online estimation phase .We are now considering the system (2) with observations(3). This phase consists of a time update and a measurementupdate. In Sec. III-A, we derived the function f ( · ) for thetime update. We next present two approaches to perform themeasurement update.
1) Approach 1 - Reverse-Mapping:
For this approach, weaim to determine the mapping of an observation y i ( k ) to theorresponding state-space region. Specifically, we constructa zonotope Z x | y i ( k ) ⊂ R n that contains all possible x ∈ R n given y i ( k ) , C i and bounded noise v i ( k ) ∈ Z v,i satisfying(3), for each i . This can be written as Z x | y i ( k ) = (cid:110) x ∈ R n (cid:12)(cid:12)(cid:12) C i x = y i ( k ) − Z v,i , (cid:111) . (8)In Prop. 1, we show how Z x | y i ( k ) can be determined usingthe SVD of the observation matrix C i . We omit the timeindex ( k ) when possible for clarity. Proposition 1.
Assume (cid:107) x (cid:107) ≤ K . Given a measurement y i ( k ) with noise v i ( k ) ∈ Z v,i = (cid:104) c v,i , G v,i (cid:105) satisfying (3),the possible states x that correspond to this measurement arecontained within the zonotope Z x | y i = (cid:104) c x | y i , G x | y i (cid:105) , where c x | y i = V Σ − r i × r i U (cid:62) (cid:0) y i ( k ) − c v,i (cid:1) ,G x | y i = (cid:2) V Σ − r i × r i U (cid:62) G v,i V M (cid:3) , (9) for all M ≥ K , with U , V , Σ and V obtained from theSVD of C i . Assuming C i has rank r i , then C i = (cid:2) U U (cid:3) (cid:20) Σ r i × r i r i × ( n − r i ) ( p i − r i ) × r i ( p i − r i ) × ( n − r i ) (cid:21) (cid:20) V (cid:62) V (cid:62) (cid:21) , (10) where a matrix with non-positive index is an empty matrix.Proof. From (10), we rewrite (3) as U Σ V (cid:62) x = y i − v i , so x = V Σ − U (cid:62) ( y i − v i ) . Since v i is bounded by Z v,i , wecan write x = V Σ − U (cid:62) (cid:0) y i − c v,i (cid:1)(cid:124) (cid:123)(cid:122) (cid:125) c x | yi − V Σ − U (cid:62) G v,i (cid:124) (cid:123)(cid:122) (cid:125) G (cid:48) x | yi β, | β | ≤ . This set corresponds to all possible x values within therange space of C i satisfying (3). By definition, if r i = n ,then V = ∅ , V spans the domain of x , and (cid:104) c x | y i , G (cid:48) x | y i (cid:105) sufficiently defines all possible x satisfying (3). However,if r i < n , V only spans a subset of the domain of x . Toensure Z x | y i contains all possible x we include a basis for ker ( C i ) in G x | y i by appending the generator V M to G x | y i ,and ensuring M ≥ K such that V M includes all x valuesin the directions of V such that (cid:107) x (cid:107) ≤ K . In both casesfor r i , the generator matrix can be written as G x | y i = (cid:104) G (cid:48) x | y i V M (cid:105) = (cid:2) V Σ − U (cid:62) G v,i V M (cid:3) , and the set Z x | y i = (cid:104) c x | y i , G x | y i (cid:105) . This result extends tothe case when r i < p i using similar argumentation in therespective cases r i = n and r i < n . Remark 2.
In our case, Z x | y i ( k ) will eventually be in-tersected with ˜ R k = (cid:104) ˜ c k , ˜ G k (cid:105) . It is therefore sufficient toset M ≥ radius ( ˜ R k ) + (cid:107) V (cid:62) ˜ c k (cid:107) instead of the moreconservative M ≥ K , where radius ( ˜ R k ) returns theradius of a minimal hypersphere containing ˜ R k [15]. Having determined the sets Z x | y i ( k ) for all i ∈ { , . . . , q } ,we can compute the measurement updated set ˆ R k given thepredicted set ˜ R k and each measurement set Z x | y i ( k ) as ˆ R k = ˜ R k ∩ qi =1 Z x | y i ( k ) , (11)which can be performed using the standard intersectionoperations presented in [15].
2) Approach 2 - Implicit Intersection:
Contrary to Ap-proach 1, here, we do not explicitly determine the sets Z x | y i ( k ) . Instead, ˆ R k is determined directly from the set ˜ R k , the measurements y i ( k ) and some weights λ ik for i ∈{ , . . . , q } . We then optimize over the weights to minimizethe volume of ˆ R k . Proposition 2.
The intersection of ˜ R k = (cid:104) ˜ c k , ˜ G k (cid:105) and the q regions for x corresponding to y i ( k ) as specified in (3) canbe over-approximated by the zonotope ˆ R k = (cid:104) ˆ c k , ˆ G k (cid:105) , ˆ c k = ˜ c k + q (cid:88) i =1 λ ik (cid:16) y i ( k ) − C i ˜ c k + c v,i (cid:17) , (12) ˆ G k = (cid:20) ( I − q (cid:80) i =1 λ ik C i ) ˜ G k λ k G v, . . . λ qk G v,q (cid:21) , (13) where λ ik ∈ R n × p i for i ∈ { , . . . , q } are weights.Proof. The proof is based on [16, Prop.1] but with zonotopesas measurements instead of strips. Let x ∈ ˜ R k ∩Z x | y ∩· · ·∩Z x | y q . Then there exists a z such that x = ˜ c k + ˜ G k z . Addingand subtracting (cid:80) qi =1 λ ik C i ˜ G k z yields x = ˜ c k + q (cid:88) i =1 λ ik C i ˜ G k z + ( I − q (cid:88) i =1 λ ik C i ) ˜ G k z. (14)From (3), we obtain C i x = y i + c v,i + G v,i d i . Using x =˜ c k + ˜ G k z yields C i ˜ G k z = y i ( k ) − C i ˜ c k + c v,i + G v,i d i ,which we insert into (14) to obtain x = ˜ c k + q (cid:88) i =1 λ ik (cid:16) y i ( k ) − C i ˜ c k + c v,i + G v,i d i (cid:17) + (cid:16) I − q (cid:88) i =1 λ ik C i (cid:17) ˜ G k z, = (cid:20) ( I − q (cid:80) i =1 λ ik C i ) ˜ G k λ k G v, . . . λ qk G v,q (cid:21)(cid:124) (cid:123)(cid:122) (cid:125) ˆ G k zd ... d q (cid:124) (cid:123)(cid:122) (cid:125) z b + ˜ c k + q (cid:88) i =1 λ ik ( y i ( k ) − C i ˜ c k + c v,i ) (cid:124) (cid:123)(cid:122) (cid:125) ˆ c k = ˆ G k z b + ˆ c k . Note that z b ∈ [ − , since b i ∈ [ − , and z ∈ [ − , . ˆ R k adheres to Def. 1 with center ˆ c k and generators ˆ G k .As in [8], we find the optimal weights λ ik ∈ R n × p i from ¯ λ ∗ = arg min ¯ λ || ˆ G k || F , (15)where ¯ λ = [ λ k . . . λ qk ] .The online estimation phase is illustrated in the blockdiagram of Fig. 1. The detailed estimation phase is presentedin Algorithm 1. The function measZon() executes Prop. 1,and optZon() Prop. 2. The function reduce ( ˜ R k +1 ) reducesthe order of ˜ R k +1 using the method proposed in [17],which ensures the number of generators in ˜ R k +1 remainsrelatively low, avoiding potential tractability issues aftermultiple iterations.ig. 1: The proposed method showing the offline learning phase yielding f ( · ) and the online estimation phase which utilizes f ( · ) to perform the time update, followed by a measurement update yielding the set ˆ R k at time-step k . Algorithm 1
Online Estimation Phase ˆ R = X k = 1 while True do ˜ R k = f ( ˆ R k − , (cid:104) u ( k − , (cid:105) ) using (7) if Approach 1 thenforeach i ∈ { , . . . , q } do Z x | y i ( k ) = measZon (cid:0) y i ( k ) , Z v,i , C i (cid:1) using (9) end ˆ R k = ˜ R k (cid:84) qi =1 Z x | y i ( k ) if Approach 2 then (cid:104) ˆ c k , ˆ G k (cid:105) = optZon ( ˜ R k , y ( k ) , C, Z v )ˆ G ∗ k , ¯ λ ∗ ← Solve (15) ˆ R k = (cid:104) ˆ c k , ˆ G ∗ k (cid:105) ˜ R k = reduce ( ˜ R k ) using [17] k ← k + 1 end C. Online Estimation Phase using Constrained Zonotopes
When intersecting zonotopes, the result is an over-approximation of the true intersection. However, it is possibleto determine the exact intersection of constrained zonotopes.
Definition 3. (Constrained Zonotope [18]) An n -dimensional constrained zonotope is C = { x ∈ R n | x = c C + G C β, A C β = b C , (cid:107) β (cid:107) ∞ ≤ } , (16) where c C ∈ R n is the center, G C ∈ R n × n g the generatormatrix and A C ∈ R n c × n g and b C ∈ R n c the constraints. Inshort, we write C = (cid:104) c C , G C , A C , b C (cid:105) . When using constrained zonotopes, we replace the timeand measurement updated sets ˜ R k and ˆ R k by the constrainedzonotopes ˜ C k and ˆ C k , respectively.
1) Approach 1 - Reverse-Mapping:
This approach worksdirectly with constrained zonotopes. The sets Z x | y i ( k ) ofProp. 1 are constrained zonotopes with no A C , b C constraints.The intersection in (11) becomes ˆ C k = ˜ C k ∩ qi =1 Z x | y i ( k ) which can be performed as described in [15].
2) Approach 2 - Implicit Intersection:
We adapt Prop. 2to use constrained zonotopes.
Proposition 3.
The intersection of ˜ C k = (cid:104) ˜ c k , ˜ G k , ˜ A k , ˜ b k (cid:105) and q regions for x corresponding to y i ( k ) as in (3) can bedescribed by the constrained zonotope ˆ C k = (cid:104) ˆ c k , ˆ G k , ˆ A k , ˆ b k (cid:105) with weights λ ik ∈ R n × p i for i ∈ { , . . . , q } where ˆ c k = ˜ c k + q (cid:88) i =1 λ ik (cid:0) y i ( k ) − C i ˜ c k + c v,i (cid:1) , ˆ G k = (cid:20) ( I − q (cid:80) i =1 λ ik C i ) ˜ G k λ k G v, . . . λ qk G v,q (cid:21) , (17) ˆ A k = ˜ A k . . . C ˜ G k − G v, . . . ... . . . C q ˜ G k . . . − G v,q , (18) ˆ b k = ˜ b k y ( k ) − C c k + c v, ... y q ( k ) − C q c k + c v,q . (19) Proof.
We follow a similar approach to [19, Thm. 6.3], butextend the proof by defining measurement sets as zonotopesinstead of strips. Z x | y i refers to Z x | y i ( k ) unless specifiedotherwise. Let x k ∈ ˜ C k ∩ Z x | y ∩ · · · ∩ Z x | y q , then thereexists a z k ∈ (cid:2) − n g × , n g × (cid:3) such that x k = ˜ c k + ˜ G k z k , ˜ A k z k = ˜ b k . (20)Using (3) and the measurement noise (cid:104) c v,i , G v,i (cid:105) , we write C i x = y i ( k ) + c v,i + G v,i d i , (21)where d i ∈ [ − , . Inserting (20) into (21) yields C i ˜ G k z k = y i ( k ) − C i ˜ c k + c v,i + G v,i d i , (22)which, combined with (20), yields ˜ A k . . . C G k − G v, . . . ... . . . C q G k . . . − G v,q (cid:124) (cid:123)(cid:122) (cid:125) ˆ A k z k d ... d q (cid:124) (cid:123)(cid:122) (cid:125) γ = ˜ b k y ( k ) − C c k + c v, ... y q ( k ) − C q c k + c v,q (cid:124) (cid:123)(cid:122) (cid:125) ˆ b k . (23)dding and subtracting (cid:80) qi =1 λ i,k C i ˜ G k z k to (20) yields x k = ˜ c k + q (cid:88) i =1 λ ik C i ˜ G k z k + ( I − q (cid:88) i =1 λ ik C i ) ˜ G k z k . (24)If we now insert (23) into (24), we obtain x = (cid:20) ( I − q (cid:80) i =1 λ ik C i ) ˜ G k λ k G v, . . . λ m i k G v,q (cid:21)(cid:124) (cid:123)(cid:122) (cid:125) ˆ G k γ + ˆ c k − + q (cid:88) i =1 λ jk (cid:0) y i ( k ) − C i ˜ c k + c v,i (cid:1)(cid:124) (cid:123)(cid:122) (cid:125) ˆ c k = ˆ G k γ + ˆ c k . Hence, x ( k ) ∈ ˆ C k and ( ˜ C ∩ Z x | y ∩ · · · ∩ Z x | y q ) ⊂ ˆ C k .Conversely, let x ( k ) ∈ ˆ C k . Then, there exists a γ suchthat (16) in Def. 3 is satisfied. Partitioning γ into γ =( z k , d . . . , d q ) , it follows that we can construct a constrainedzonotope ˜ C k = { ˜ c k , ˜ G k , ˜ A k , ˜ b k } given that (cid:107) z k (cid:107) ∞ ≤ .Thus, x ( k ) ∈ ˜ C . Similarly, we can get the constraints in (21).Inserting (20) in (22) results in obtaining all the equationsin (21). Therefore, x ( k ) ∈ Z x | y i ( k ) , ∀ i ∈ { , . . . , q } . Thus, x ( k ) ∈ ( ˜ C k ∩ Z x | y ∩ · · · ∩ Z x | y q ) and ˆ C k ⊂ ( ˜ C k ∩ Z x | y ∩· · · ∩ Z x | y q ) , which concludes the proof.IV. E VALUATION
We evaluate our method by considering an input-drivenvariant of the rotating target described in [8]. We set A = (cid:20) . − . . . (cid:21) , B = (cid:20) . (cid:21) (25)with q = 3 measurements parameterized as follows C = (cid:2) . (cid:3) , C = (cid:2) . − . (cid:3) , C = (cid:20) − . .
20 0 . (cid:21) , Z v, = (cid:104) , (cid:105) , Z v, = (cid:104) , (cid:105) , Z v, = (cid:104) [0 0] (cid:62) , I (cid:105) . The noise signals are characterized by the zonotopes Z γ = (cid:104) [0 0] (cid:62) , . I (cid:105) and Z w = (cid:104) [0 0] (cid:62) , . I (cid:105) . We run the offline learning phase with T = 500 and inputs sampleduniformly from the set U = (cid:104) , (cid:105) . The noise signals v v,i ( k ) , w ( k ) and γ ( k ) are sampled uniformly from theirrespective zonotope sets using the command randPoint ( Z ) as described in [15].After learning f ( · ) , we run the online estimation phase .The initial state set is X = (cid:104) [0 0] (cid:62) , I (cid:105) and the trueinitial state is x (0) = (cid:2) −
10 10 (cid:3) (cid:62) . Once again, we samplethe inputs uniformly from U . We evaluate both the zonotopeand constrained zonotope methods, each time using either ofthe two proposed measurement update approaches. Fig. 2ashows the bounds of ˆ R k in the x state dimension forboth approaches. Fig. 2b shows the equivalent results whenour method uses constrained zonotopes. As expected, x ( k ) is always contained within ˆ R k (or ˆ C k ) at each time step.Although both measurement update approaches yield similarset sizes on average, the set evolution of Approach 2 iscomparatively smoother. We observed that the intersectionoperation from [15] used to compute (11) in Approach 1yielded sets with varying degrees of over-approximation,resulting in the jagged set evolutions seen in Fig. 2. By (a) Using zonotopes showing bounds of ˆ R k in x (b) Using constrained zonotopes showing bounds of ˆ C k in x Fig. 2: Bounds of the set ˆ R k in (a), and ˆ C k in (b), projectedonto the first state dimension x of x ( k ) using measurementupdate approaches 1 and 2. Fig. 3: Sets ˆ R k using measurement update approaches 1 and2, and the equivalent sets ˆ C k using constrained zonotopes( CZ ), compared to the KF’s σ confidence bounds.improving this intersection method, we expect Approach 1to yield a smoother and tighter state set evolution.Furthermore, we compare our results with N4SID subspaceidentification [20] combined with a Kalman filter (KF). InFig. 3, we show the sets ˆ R k and ˆ C k , using either mea-surement update approach, using zonotopes or constrainedzonotopes. We also show the ellipse corresponding to the σ uncertainty bound of the KF estimate, indicating that ourestimator provides state sets comparable in size to that of theKF.Referring to both Fig. 2 and Fig. 3, it is clear that the con-strained zonotopes yield smaller state sets at each time step.However, this comes at the cost of increased computationalload. Running our simulations on a Dell laptop with an 8-core i5-8365U processor at 1.6GHz, the average computationtime per iteration for Approach 1 increased from . s to . s. when using constrained zonotopes; for Approach 2,the corresponding times were . s and . s respectively.For all our approaches, we observed that reducing the orderof the sets, which reduces the number of generators in ˆ R (or ˆ C ), was critical to keep the computational load low. -20-15-10-5051015 Fig. 4: ˆ R projected onto the x axis for varying conditionnumbers κ = σ max σ min of the data matrix D .We also investigate how the condition number of the datamatrix D , used during the offline learning phase , affects thesize of ˆ R and ˆ C . We denote the condition number of D by κ = σ max σ min , where σ max and σ min denote the max. and min.singular values of D respectively. This ratio gives an insightinto the degree to which the full row rank condition on D , issatisfied. A higher condition number is expected to yield alarger state set obtained from f ( · ) for a fixed initial state setand input. To evaluate this, we run identical simulations,but with a decreasing amplitude of the input u ( k ) (whichnaturally increases κ ). We then plot the bounds of the set ˆ R (and ˆ C ) for varying κ as shown in Fig. 4. Indeed, we observethat higher condition numbers yield large state sets. However,we note that for κ values higher than approximately , ˆ R remains the same. This is because ˆ R is the result of theintersection of ˜ R and the measurement-sets Z x | y i (3) , andwhen ˜ R is very large, ˆ R is essentially only the intersectionthe measurement sets. We conclude that in order to ensureadequately small state sets, we must ensure that the conditionnumber of D is sufficiently small. For the simulations shownin Fig. 2 and Fig. 3, the condition number was κ = 1 . .V. C ONCLUSIONS AND R ECOMMENDATIONS
In this paper, we introduced a novel zonotope-basedmethod to perform set-based state estimation with set con-tainment guarantees using a data-driven set propagationfunction. We presented two approaches to perform the mea-surement update which merges the time updated state setwith the observed measurements. We extended our methodto use constrained zonotopes, which yielded smaller state setsat the cost of increased computational load. Our results showstate sets comparable in size to the σ uncertainty boundsobtained when running N4SID subspace identification and aKalman filter, but with the added feature of set-containmentguarantees and without requiring any knowledge of thestatistical properties of the noise.Future work includes evaluating our proposed estimator onreal-world examples as well as gaining more insight into thelimitations of our method when applied to more complexdynamical systems. Additionally, improving the zonotopeintersection operation to lessen the degree of overapproxi-mation of the resultant state set would yield tighter state setestimates at each time step. R
EFERENCES[1] M. Althoff,
Reachability analysis and its application to the safetyassessment of autonomous cars . PhD thesis, Technische UniversitätMünchen, 2010.[2] F. Mazenc and O. Bernard, “Interval observers for linear time-invariantsystems with disturbances,”
Automatica , vol. 47, no. 1, pp. 140–147,2011.[3] G. Belforte, B. Bona, and V. Cerone, “Parameter estimation algorithmsfor a set-membership description of uncertainty,”
Automatica , vol. 26,no. 5, pp. 887–898, 1990.[4] L. Ma, Z. Wang, H.-K. Lam, and N. Kyriakoulis, “Distributed event-based set-membership filtering for a class of nonlinear systems withsensor saturations over sensor networks,”
IEEE Transactions on Cy-bernetics , vol. 47, no. 11, pp. 3772–3783, 2016.[5] L. Orihuela, S. Roshany-Yamchi, R. A. García, and P. Millán,“Distributed set-membership observers for interconnected multi-ratesystems,”
Automatica , vol. 85, pp. 221–226, 2017.[6] C. Durieu, E. Walter, and B. Polyak, “Multi-input multi-output ellip-soidal state bounding,”
Journal of Optimization Theory and Applica-tions , vol. 111, no. 2, pp. 273–303, 2001.[7] J. Blesa, V. Puig, and J. Saludes, “Robust fault detection usingpolytope-based set-membership consistency test,”
IET Control Theory& Applications , vol. 6, no. 12, pp. 1767–1777, 2012.[8] A. Alanwar, J. J. Rath, H. Said, and M. Althoff, “Distributed set-basedobservers using diffusion strategy,” arXiv:2003.10347 , 2020.[9] J. C. Willems, P. Rapisarda, I. Markovsky, and B. L. De Moor, “Anote on persistency of excitation,”
Systems & Control Letters , vol. 54,no. 4, pp. 325–329, 2005.[10] D. Alpago, F. Dörfler, and J. Lygeros, “An Extended Kalman Filterfor Data-Enabled Predictive Control,”
IEEE Control Systems Letters ,vol. 4, no. 4, pp. 994–999, 2020.[11] C. De Persis and P. Tesi, “Formulas for data-driven control: Stabi-lization, optimality, and robustness,”
IEEE Transactions on AutomaticControl , vol. 65, no. 3, pp. 909–924, 2019.[12] J. Berberich, J. Köhler, M. A. Müller, and F. Allgöwer, “Data-drivenmodel predictive control with stability and robustness guarantees,”
IEEE Transactions on Automatic Control , 2020.[13] A. Alanwar, A. Koch, F. Allgöwer, and K. H. Johansson, “Data-drivenreachability analysis using matrix zonotopes,” arXiv:2011.08472 ,2020.[14] W. Kühn, “Rigorously computed orbits of dynamical systems withoutthe wrapping effect,”
Computing , vol. 61, no. 1, pp. 47–67, 1998.[15] M. Althoff, “An introduction to CORA 2015,” in
Proceedings of theWorkshop on Applied Verification for Continuous and Hybrid Systems ,2015.[16] V. T. H. Le, C. Stoica, T. Alamo, E. F. Camacho, and D. Dumur,“Zonotope-based set-membership estimation for multi-output uncer-tain systems,” in
IEEE International Symposium on Intelligent Control ,pp. 212–217, 2013.[17] A. Girard, “Reachability of uncertain linear systems using zonotopes,”in
Hybrid Systems: Computation and Control , LNCS 3414, pp. 291–305, 2005.[18] J. K. Scott, D. M. Raimondo, G. R. Marseglia, and R. D. Braatz,“Constrained zonotopes: A new tool for set-based estimation and faultdetection,” vol. 69, pp. 126–136, 2016.[19] A. Alanwar, V. Gassmann, X. He, H. Said, H. Sandberg, K. H.Johansson, and M. Althoff, “Privacy preserving set-based estimationusing partially homomorphic encryption,” arXiv:2010.11097 , 2020.[20] P. Van Overschee and B. De Moor, “N4SID: Subspace algorithmsfor the identification of combined deterministic-stochastic systems,”