Over-approximating reachable tubes of linear time-varying systems
aa r X i v : . [ m a t h . O C ] F e b Serry, Reissig Computing reachable tubes 1
Over-approximating reachable tubes of lineartime-varying systems
Mohamed Serry and Gunther Reissig
Abstract
We present a method to over-approximate reachable tubes over compact time-intervals, for linearcontinuous-time, time-varying control systems whose initial states and inputs are subject to compact convexuncertainty. The method uses numerical approximations of transition matrices, is convergent of first order,and assumes the ability to compute with compact convex sets in finite dimension. We also present a variantthat applies to the case of zonotopic uncertainties, uses only linear algebraic operations, and yields zonotopicover-approximations. The performance of the latter variant is demonstrated on an example.
Index Terms
Reachability, linear time-varying systems, MSC: Primary, 93B03; Secondary, 34A60
I. Introduction
Reachable (or attainable ) sets and tubes are central concepts in systems and control theory, withmyriads of applications. See, e.g. [1]–[7] and the references given therein. The efficient computationof accurate approximations of these sets is a challenging problem whose diverse variants have beenattracting research attention for decades. In this paper, we focus on over-approximating reachabletubes of linear time-varying control systems of the form ˙ x ( t ) = A ( t ) x ( t ) + B ( t ) u ( t ) (1)over compact time-intervals [ t , t f ] , where A : [ t , t f ] → R n × n and B : [ t , t f ] → R n × m are time-varyingmatrices. Both the initial state x ( t ) and the input signal u are subject to compact convex uncertainty.The problem is mathematically formalized in Section III.In numerous applications it is critical to formally verify that all solutions of the system (1) alwaysavoid certain predefined unsafe regions ; see, e.g. [4, Sect. 3] and the references given therein. That is,these applications require proof that the reachable tube over the time-interval [ t , t f ] (and not onlythe reachable set at some time t ∈ [ t , t f ] ) of the system (1) does not intersect any unsafe region. Astubes cannot, in general, be determined exactly, intersection tests need to rely on over-approximations(and not on mere approximations) in place of the actual tubes. The over-approximations should beas precise as possible to avoid excessive conservatism of the verification, and need to be representedin a form that facilitates to reliably and efficiently verify disjointness from unsafe regions.One of the earliest techniques of reachability analysis, the hyperplane method , approximates reach-able sets by intersections of supporting halfspaces and by convex hulls of the respective support points[1], [8]. More recent techniques rely on a variety of additional classes of sets including, e.g. ellipsoids,hyper-intervals, and zonotopes [5]–[7], [9]–[25]. As for reachable tubes, the standard approach todayis to apply the method proposed in [9] or one of its extensions, e.g. [7], [10]–[13], which compute over-approximations in the form of finite unions of zonotopes [7], [9], [10] and of more general convex sets This work has been supported by the German Research Foundation (DFG) under grant no. RE 1249/4-1.This work has beenaccepted for publication in the
IEEE Trans. Automatic Control
Serry, Reissig Computing reachable tubes [11]–[13]. As a result of such representation, disjointness from a polyhedral (or convex) unsafe regioncan be verified by solving a linear (or convex) feasibility problem. While the method is particularlyefficient and converges, i.e., it is capable of producing arbitrarily precise over-approximations, itsapplication is limited to the time-invariant special case of (1). Its extension in [14] additionally allowsfor uncertain coefficients A and B in (1), but does not converge even if A and B are precisely known,in which case (1) is again required to be time-invariant.Another prominent class of methods, ellipsoidal techniques [5], solve the more general problem offeedback synthesis for linear time-varying plants with two competing inputs. When applied to thesystem (1), these methods yield a set-valued function E defined on the interval [ t , t f ] whose valueat any time is a finite intersection of ellipsoids containing the reachable set at that time as a subset.While arbitrarily precise over-approximations are obtained when a sufficient amount of ellipsoids iscomputed, the approach suffers from two shortcomings. Firstly, the ellipsoids result from numericallysolving linear-quadratic optimal control problems derived from (1), yet numerical errors incurredin the course of the solution are not taken into account. Hence, mere approximations rather thanover-approximations might actually be computed. Secondly, approximations of reachable tubes areobtained only implicitly, as the union over t ∈ [ t , t f ] of E ( t ) , and so they are disjoint from an unsaferegion R if and only if the graph of the set-valued map E is disjoint from the set [ t , t f ] × R . Verifyingthe latter condition is a great challenge since the graph of E is not, in general, convex. The issue hasso far been resolved only for the time-invariant special case of (1); see [15]. Moreover, while ellipsoidaltechniques have been generalized to handle nonlinear dynamics, the extensions still suffer from boththe aforementioned shortcomings, e.g. [16].Other approaches use differential inequalities, comparison principles, interval arithmetic, and com-binations thereof, and compute interval over-approximations [6], [17]–[20]. While these techniquesmay allow for uncertain coefficients A and B in (1) [17] or even for nonlinear dynamics [6], [18]–[20],they are all conservative, i.e., arbitrarily precise over-approximations of reachable tubes cannot beobtained, and the methods in [6], [17]–[19] additionally suffer from both shortcomings mentioned inour discussion of ellipsoidal techniques. Finally, the reachable tube can also be characterized as asublevel set of the viscosity solution of a partial differential equation called Hamilton-Jacobi-Bellmanequation [5]. However, solving the latter numerically is avoided in practice as this would requirediscretizing the state space and so the computational effort would scale exponentially with the statespace dimension.To conclude, efficient methods to compute arbitrarily precise over-approximations of reachabletubes of the system (1), that are additionally represented in a form suitable for formal verificationpurposes, are currently limited to the time-invariant special case of (1). This is in stark contrast tothe importance of the general time-varying case of (1) in several fields of application, e.g. [26].In Section IV-A of this paper, we present a method that produces over-approximations that areconvergent of first order, does not require discretization of either the input or the state space, usesnumerical approximations of transition matrices rather than closed-form solutions, and assumes theability to compute with compact convex sets in finite dimension. A variant that applies to thecase of zonotopic uncertainties, uses only linear algebraic operations, and yields zonotopic over-approximations, is subsequently presented in Section IV-B. In Section V, we demonstrate the perfor-mance of the latter variant on an example.
II. Preliminaries
A. Notation
Given two sets A and B and a positive integer p , B \ A and A × B denotes the relative complementof the set A in the set B , and the product of A and B , respectively, and A p = A × · · · × A ( p factors). R , R + , Z and Z + denote the sets of real numbers, non-negative real numbers, integers andnon-negative integers, respectively, and N = Z + \ { } . [ a, b ] , ] a, b [ , [ a, b [ , and ] a, b ] denote closed, open erry, Reissig Computing reachable tubes 3 and half-open, respectively, intervals with end points a and b , e.g. [0 , ∞ [ = R + . [ a ; b ] , ] a ; b [ , [ a ; b [ , and ] a ; b ] stand for discrete intervals, e.g. [ a ; b ] = [ a, b ] ∩ Z , [1; 4[ = { , , } , and [0; 0[ = ∅ .Given any map f : A → B , the image of a subset C ⊆ A under f is denoted f ( C ) , f ( C ) = { f ( c ) | c ∈ C } . We denote the identity map X → X : x x by id , where the domain of definition X will always be clear form the context.Arithmetic operations involving subsets of a linear space X are defined pointwise, e.g. αM := { αy | y ∈ M } and the Minkowski sum M + N := { y + z | y ∈ M, z ∈ N } , if α ∈ R and M, N ⊆ X .The convex hull of M is denoted conv( M ) . By k· k we denote any norm on X , B ⊆ X is the closed unitball w.r.t. k · k , and the norm of a non-empty subset M ⊆ X is defined by k M k := sup x ∈ M k x k . Themaximum norm on R n is denoted k · k ∞ , k x k ∞ = max {| x i | | i ∈ [1; n ] } for all x ∈ R n . The Hausdorffdistance d H is defined in the Appendix.We say that a map is of class C k if it is continuous and k times continuously differentiable, k ∈ Z + .Given a non-empty set X ⊆ R n and a compact interval [ a, b ] ⊆ R , X [ a,b ] denotes the set of allmeasurable maps [ a, b ] → X . Integration is always understood in the sense of Lebesgue. Given normson R n and R m , the linear space R n × m of n × m matrices is endowed with the usual matrix norm, k A k = sup k x k≤ k Ax k for A ∈ R n × m .We use the asymptotic notation O ( · ) in the usual way [27]. In particular, let X ⊆ R n , f : F ⊆ R × X → R + , g : G ⊆ R → R + , H : F × R + → R + and a ∈ R be such that a = lim i →∞ s i for somesequence ( s i , x i ) i ∈ N in F , and suppose that s ∈ G whenever ( s, x ) ∈ F . Then f ( s, x ) ≤ H ( s, x, O ( g ( s ))) as s → a , uniformly in x , if there exist k : G → R + and a neighborhood U ⊆ R of a such that k ( s ) = O ( g ( s )) as s → a and f ( s, x ) ≤ H ( s, x, k ( s )) whenever ( s, x ) ∈ F ∩ ( U × X ) , and similarly for a ∈ {∞ , −∞} . B. Linear Time-Varying Control Systems
Given u : [ t , t f ] → R m , a map x : [ t , t f ] → R n is a solution of the system (1) (generated by u ) if x isabsolutely continuous and (1) holds for (Lebesgue) almost every t ∈ [ t , t f ] . We shall always assumethat A and B are continuous and that u is integrable, which implies both existence and uniquenessof solutions [28]. The general solution of the system (1) is the map ϕ defined by the requirement thatfor all p ∈ R n , s ∈ [ t , t f ] and integrable u , ϕ ( · , s, p, u ) is the unique solution of (1) defined on [ t , t f ] and satisfying ϕ ( s, s, p, u ) = p . The map ϕ ( t, s, · , , which is linear, is called the transition matrix at ( t, s ) of the system and is denoted by φ ( t, s ) . The map φ : [ t , t f ] × [ t , t f ] → R n × n is of class C , andthe identities ϕ ( t, s, p, u ) = φ ( t, s ) p + Z ts φ ( t, τ ) B ( τ ) u ( τ ) dτ,φ ( s, s ) = id , and φ ( t, s ) φ ( s, T ) = φ ( t, T ) hold for all s, t, T ∈ [ t , t f ] , all p ∈ R n , and all integrable u ;see, e.g. [28]. Moreover, D φ ( t, s ) = A ( t ) φ ( t, s ) and D φ ( s, t ) = − φ ( s, t ) A ( t ) hold for all s, t ∈ [ t , t f ] ,where D i φ denotes the partial derivative of φ with respect to (w.r.t.) the i th argument. If A isadditionally of class C k , k ≥ , then φ is of class C k +1 . Finally, Gronwall’s lemma implies k φ ( t, s ) k ≤ e | t − s | M and k φ ( t, s ) − id k ≤ e | t − s | M − (2)for all s, t ∈ [ t , t f ] , provided that k A ( t ) k ≤ M for all t ∈ [ t , t f ] . C. Reachable Sets and Tubes
Given non-empty, compact, convex subsets X ⊆ R n and U ⊆ R m , and a, b, t ∈ [ t , t f ] satisfying a ≤ b , the sets R ( t ) = (cid:8) ϕ ( t, t , x , u ) (cid:12)(cid:12) x ∈ X , u ∈ U [ t ,t ] (cid:9) , R ([ a, b ]) = [ s ∈ [ a,b ] R ( s ) Serry, Reissig Computing reachable tubes are the reachable set at time t and the reachable tube over the time interval [ a, b ] , respectively, of thesystem (1). Both R ( t ) and R ([ a, b ]) are non-empty and compact, and R ( t ) is additionally convexand is conveniently written using a set-valued integral, R ( t ) = φ ( t, t ) X + R tt φ ( t, s ) B ( s ) U ds.
See,e.g. [29]. Moreover, the well-known semi-group property of reachable sets [30] yields the identity R ( b ) = φ ( b, a ) R ( a ) + Z ba φ ( b, s ) B ( s ) U ds. (3)
III. Problem Statement
We consider the system (1), where both the initial state x ( t ) ∈ X and the input u ( t ) ∈ U aresubject to uncertainty, represented by the set X and U , respectively. We assume the following.( H ) n ∈ N , t , t f ∈ R and t < t f .( H ) X ⊆ R n and U ⊆ R m are non-empty, compact, and convex.( H ) A and B are of class C , and k A ( t ) k ≤ M A , k ˙ A ( t ) k ≤ M ˙ A , k B ( t ) k ≤ M B , and k ˙ B ( t ) k ≤ M ˙ B for all t ∈ [ t , t f ] , where M A , M ˙ A , M B , M ˙ B ∈ R and M A > . Here, ˙ A denotes the derivative of themap A : [ t , t f ] → R n × n , and similarly for ˙ B .( H ) Denote D = { ( t, s ) ∈ [ t , t f ] × [ t , t f ] | t ≥ s } . Then k φ ( t, s ) − e φ ( t, s ) k ≤ θ ( t − s ) for all ( t, s ) ∈ D, (4) θ ( h ) = O ( h ) as h → , (5)where e φ : D → R n × n approximates the transition matrix φ of (1) and θ : R + → R + is monotonicallyincreasing.We note that ( H ) is the requirement that the approximation e φ of φ has consistency order [31, Def. 4.7]. Under assumptions ( H )-( H ), this requirement is satisfied by the vast majority ofnumerical methods to solve initial value problems. See, e.g. [31, Example 4.8], as well as Lemma A.4in the Appendix.The problem data t , t f , A , B , X and U are fixed throughout the paper, and so are the constants M A , M ˙ A , M B , and M ˙ B , as well as the functions φ , e φ and θ and the set D . Throughout the paper, allthat data is subject to the standing hypotheses ( H )-( H ). III.1 Problem.
Devise a convergent method that over-approximates R ([ t , t f ]) , in the sense thatgiven the problem data and a time discretization parameter N , a superset b R N of R ([ t , t f ]) is obtained,satisfying b R N → R ([ t , t f ]) in Hausdorff distance as N → ∞ . IV. Proposed method
In order to solve Problem III.1 for any given value of the time discretization parameter N , weshall over-approximate reachable sets R ( t i ) and reachable tubes R ([ t i , t i +1 ]) of the control system(1) for equidistant points of time t i ∈ [ t , t f ] , i ∈ [0; N ] . The approximation will be convergent offirst order [31], meaning that the Hausdorff distance between R ( t i ) and its approximation is of order O (1 /N ) , and similarly for tubes. The respective method of over-approximation, presented in SectionIV-A, applies to any uncertainty sets X and U satisfying Hypothesis ( H ) and assumes the abilityto compute with compact convex sets in finite dimension. Our algorithmic solution subsequentlypresented in Section IV-B applies to the case of zonotopic uncertainties, uses only linear algebraicoperations, and involves an additional approximation step that yields zonotopic over-approximationsof reachable tubes retaining first order accuracy. erry, Reissig Computing reachable tubes 5 A. Over-approximation of Reachable Sets and Reachable Tubes
We consider the system (1) under our standing hypotheses ( H )-( H ). Given a time discretizationparameter N ∈ N , we define sequences (Ω i ) i ∈ [0; N ] and (Γ i ) i ∈ [1; N ] of subsets Ω i , Γ i ⊆ R n by the followingrequirements for all i ∈ [1; N ] . h = ( t f − t ) /N and t i = t + ih, (6a) Ω = X , (6b) Ω i = e φ ( t i , t i − )Ω i − + hB ( t i ) U + ( α h + θ h k Ω i − k ) B , (6c) Γ i = conv (Ω i − ∪ (Ω i + ( β h + γ h k Ω i − k ) B )) . (6d)Here, k · k denotes any norm and the maps α, β, γ : R + → R + are defined by r ( s ) = exp( sM A ) − − sM A , β ( s ) = s M ˙ B k U k , (6e) α ( s ) = r ( s ) k U k M ˙ B + M A M B M A , γ ( s ) = r ( s ) (cid:18) M ˙ A M A (cid:19) . (6f)For convenience, here and in the sequel we often use α h in place of α ( h ) , and similarly for β , γ and θ .By (6a), we define an equidistant grid with step size h , of points t , . . . , t N , spanning the timeinterval [ t , t f ] . The requirements (6b)-(6c) iteratively define sets Ω i , which are supposed to approx-imate the reachable sets R ( t i ) , and in turn, (6d) uses these approximations as well as their inflatedversions to define sets Γ i , which are supposed to approximate reachable tubes R ([ t i − , t i ]) . As we shallshow, due to our careful definition of the maps α , β and γ depending on the time-varying problemdata, both Ω i and Γ i actually are over-approximations, with approximation error of order O (1 /N ) .We now set out to state formally and to prove what we have just described in informal terms.In doing so, we shall use the superscript N to indicate that, e.g. the sequence (Ω Ni ) i ∈ [0; N ] has beencomputed by our method (6) for a specific value of the time discretization parameter, and similarlyfor h , t i and Γ i . IV.1 Proposition (Reachable Sets).
For each N ∈ N , let sequences ( t Ni ) i ∈ [0; N ] and (Ω Ni ) i ∈ [0; N ] bedefined by (6a)-(6c) and (6e)-(6f ).Then R ( t Ni ) ⊆ Ω Ni for all N ∈ N and all i ∈ [0; N ] , and d H ( R ( t Ni ) , Ω Ni ) ≤ O (1 /N ) as N → ∞ ,uniformly w.r.t. i . For our proof, we need the following auxiliary results.
IV.2 Lemma.
We have the estimate d H ( I ( a, b ) , J ( a, b )) ≤ α ( b − a ) whenever t ≤ a ≤ b ≤ t f , where I ( a, b ) = R ba φ ( b, s ) B ( s ) U ds , J ( a, b ) = ( b − a ) B ( b ) U , and α is defined in (6f ).Proof. Let a, b ∈ [ t , t f ] , a ≤ b . The assumption ( H ) on U implies J ( a, b ) = R ba B ( b ) U ds , and usingFilippov’s Lemma [29], we obtain d H ( I ( a, b ) , J ( a, b )) ≤ k U k Z ba k φ ( b, s ) B ( s ) − B ( b ) k ds. (7)Next, using (2) and the identity φ ( b, s ) B ( s ) − B ( b ) = Z sb φ ( b, z )( ˙ B ( z ) − A ( z ) B ( z )) dz for all s ∈ [ t , t f ] , we see that the integrand in (7) is bounded by ( M ˙ B + M A M B )(e ( b − s ) M A − /M A ,which proves the lemma. IV.3 Lemma.
Let a ∈ R + , b ∈ R , K : N → N , and for each N ∈ N , let ( x Ni ) i ∈ [0; K ( N )] be a sequencein R + . Suppose that K ( N ) = O ( N a ) , x N = O ( N a + b ) and x Ni ≤ (1 + O ( N − a )) x Ni − + O ( N b ) hold as Serry, Reissig Computing reachable tubes N → ∞ , uniformly w.r.t. i .Then x Ni ≤ O ( N a + b ) as N → ∞ , uniformly w.r.t. i .Proof. By our hypotheses, there exist maps p, q, r : N → R + satisfying p ( N ) = O ( N a + b ) , q ( N ) = O ( N − a ) and r ( N ) = O ( N b ) as N → ∞ , and x N ≤ p ( N ) and x Ni ≤ (1 + q ( N )) x Ni − + r ( N ) (8)for all sufficiently large N ∈ N and all i ∈ [1; K ( N )] . Define f ( N, i ) = (1 + q ( N )) i for all N ∈ N and all i ∈ [0; K ( N )] , to arrive at f ( N, i ) ≤ exp( q ( N ) K ( N )) . Then, by our assumptions on q and K , the map f is bounded. In view of (8) and the variation-of-constants formula we conclude that x Ni ≤ O ( N a + b ) as claimed. Proof of Proposition IV.1.
For the sake of simplicity, throughout this proof we drop the superscript N from our notation. Let h be defined by (6a).The first claim holds for i = 0 and all N ∈ N as R ( t ) = X = Ω . Assume that R ( t i ) ⊆ Ω i holds for some N ∈ N and some i ∈ [0; N [ . Then, using the identity (3) and Lemma IV.2 as wellas Lemma A.2(v), we obtain R ( t i +1 ) ⊆ φ ( t i +1 , t i )Ω i + hB ( t i +1 ) U + α ( h ) B . Moreover, φ ( t i +1 , t i )Ω i ⊆ e φ ( t i +1 , t i )Ω i + θ ( h ) k Ω i k by the estimate (4) and Lemma A.2(iii)(v), and so R ( t i +1 ) ⊆ Ω i +1 .To prove the second claim, we use the triangle inequality, assumption ( H ) and estimates (2), (4)and (5) to obtain the bound k e φ ( t i , t i − ) k ≤ O (1 /N ) as N → ∞ , uniformly w.r.t. i . In turn, (5),(6c), ( H ) and the fact that α ( s ) = O ( s ) as s → together imply k Ω i k ≤ (1 + O (1 /N )) k Ω i − k + O (1 /N ) , and so k Ω i k ≤ O (1) as N → ∞ , uniformly w.r.t. i , by Lemma IV.3. It follows that α ( h ) + θ ( h ) k Ω i − k ≤ O (1 /N ) and d H ( e φ ( t i , t i − )Ω i − , φ ( t i , t i − ) R ( t i − )) ≤ (1 + O (1 /N )) e i − + O (1 /N ) ,where e i = d H ( R ( t i ) , Ω i ) . Moreover, d H ( hB ( t i ) U, R t i t i − φ ( t i , s ) B ( s ) U ds ) ≤ O (1 /N ) by Lemma IV.2,and so e i ≤ (1 + O (1 /N )) e i − + O (1 /N ) as N → ∞ , uniformly w.r.t. i . Then e i ≤ O (1 /N ) by LemmaIV.3, as claimed.The following theorem, and its corollary immediately obtained using Lemma A.2(vii), provide a firstsolution to Problem III.1. IV.4 Theorem (Reachable Tubes).
For each N ∈ N , let sequences ( t Ni ) i ∈ [0; N ] and (Γ Ni ) i ∈ [1; N ] bedefined by (6).Then R ([ t i − , t i ]) ⊆ Γ Ni for all N ∈ N and all i ∈ [1; N ] , and d H ( R ([ t i − , t i ]) , Γ Ni ) ≤ O (1 /N ) as N → ∞ , uniformly w.r.t. i . IV.5 Corollary.
Under the hypotheses and in the notation of Theorem IV.4, denote b R N = ∪ i ∈ [1; N ] Γ Ni .Then R ([ t , t f ]) ⊆ b R N for all N ∈ N , and d H ( R ([ t , t f ]) , b R N ) = O (1 /N ) as N → ∞ . Our proof of Theorem IV.4 uses the following auxiliary result.
IV.6 Lemma.
Let γ be defined by (6f ), and let a, b ∈ [ t , t f ] with a < b . Then k φ ( t, a ) − ψ ( t, a, b ) k ≤ ( t − a )( b − a ) − γ ( b − a ) for all t ∈ [ a, b ] , where ψ ( t, a, b ) = id +( φ ( b, a ) − id)( t − a ) / ( b − a ) .Proof. The claim is obvious for t ∈ { a, b } , so we suppose that t ≤ a < t < b ≤ t f . Then, by a changeof variable, t − ab − a Z ba A ( τ ( s )) φ ( τ ( s ) , a ) ds = Z ta A ( s ) φ ( s, a ) ds, where τ ( s ) = ( t − a )( s − a ) / ( b − a ) + a ∈ [ a, s ] , and so the difference φ ( t, a ) − ψ ( t, a, b ) can be writtenas t − ab − a Z ba A ( τ ( s )) φ ( τ ( s ) , a ) − A ( s ) φ ( s, a ) ds. (9)As the map s A ( s ) φ ( s, a ) is smooth, the integrand in (9) takes the form R τ ( s ) s ( ˙ A ( z )+ A ( z ) ) φ ( z, a ) dz .The claim then follows from the estimate (2) and assumption ( H ). erry, Reissig Computing reachable tubes 7 We mention in passing that, in the time-invariant case of (1) with B ( t ) = id , our estimates inLemmas IV.2 and IV.6 reduce to those in [11, Lemma 2] and [11, p. 260, last inequ.], respectively.Another related but less precise estimate is given in [15, Lemma 1]. The mathematical tools we haveused to treat the general time-varying case are quite different from the ones used in [11], [15]. Proof of Theorem IV.4.
For the sake of simplicity, throughout this proof we drop the superscript N from our notation. Moreover, we do not mention the domains N , [1; N ] and [ t i − , t i ] of N , i and t , andasymptotic estimates are always meant to hold for N → ∞ , uniformly w.r.t. i and t . The map ψ isdefined in Lemma IV.6, and h is defined in (6a).We claim that R ( t ) ⊆ E i,t for all N , i and t , where E t,i = ψ ( t, t i − , t i ) R ( t i − ) + ( t − t i − ) h − M i ,M i = hB ( t i ) U + ( α h + β h + γ h k Ω i − k ) B . Indeed, using Lemmata IV.6 and A.2(iii),(v), the compactness of reachable tubes, and PropositionIV.1, we see that φ ( t, t i − ) R ( t i − ) ⊆ ψ ( t, t i − , t i ) R ( t i − ) + t − t i − h γ h k Ω i − k B . (10)Moreover, we obviously have k B ( t i ) − B ( t ) k ≤ hM ˙ B , and in turn, d H (( t − t i − ) B ( t ) U, ( t − t i − ) B ( t i ) U ) ≤ ( t − t i − ) β h /h by Lemma A.2(iii). Then Lemmata IV.2 and A.2(v), the compactness of U , and thefact that α ( s ) /s is monotonically increasing in s , imply Z tt i − φ ( t, s ) B ( s ) U ds ⊆ t − t i − h ( hB ( t i ) U + ( α h + β h ) B ) (11)for all N , i and t . Our claim then follows from the identity (3). Moreover, the estimates from whichthe inclusions (10) and (11) have been obtained also show that d H ( R ( t ) , E i,t ) ≤ O (1 /N ) .Next observe that the set E i,t takes the form { (1 − λ ) x + λφ ( t i , t i − ) x + λm | x ∈ R ( t i − ) , m ∈ M i } ,where λ = ( t − t i − ) /h , and so E i,t ⊆ F i,t for all N , i and t , where F i,t is defined to be the set { (1 − λ ) x + λφ ( t i , t i − ) y + λm | x, y ∈ R ( t i − ) , m ∈ M i } . Moreover, if z ∈ F i,t , then there exist x, y ∈ R ( t i − ) and m ∈ M i satisfying z = (1 − λ ) x + λφ ( t i , t i − ) y + λm . We define x ′ = (1 − λ ) x + λy ∈ R ( t i − ) and z ′ = (1 − λ ) x ′ + λφ ( t i , t i − ) x ′ + λm to obtain z − z ′ = λ (1 − λ )( φ ( t i , t i − ) − id)( y − x ) . The estimate (2) then implies k z − z ′ k ≤ (exp( hM A ) − k Ω i − k / ,and as k Ω i k ≤ O (1) by Proposition IV.1 and the compactness of reachable tubes, we arrive at d H ( E i,t , F i,t ) ≤ O (1 /N ) .So far, we have shown that R ( t ) ⊆ F i,t for all N , i and t , and that d H ( R ( t ) , F i,t ) ≤ O (1 /N ) . Itfollows that R ([ t i − , t i ]) ⊆ ∪ t ∈ [ t i − ,t i ] F i,t , and by Lemma A.2(vii), the Hausdorff distance between thetwo sets does not exceed O (1 /N ) . Next observe that ∪ t ∈ [ t i − ,t i ] F i,t = conv ( R ( t i − ) ∪ ( φ ( t i , t i − ) R ( t i − ) + M i )) by Lemma A.1(ii), and that φ ( t i , t i − ) R ( t i − ) ⊆ e φ ( t i , t i − )Ω i − + θ ( h ) k Ω i − k B by Lemma A.2(iii),(v),the estimate (4), and Proposition IV.1. Thus, ∪ t ∈ [ t i − ,t i ] F i,t ⊆ Γ i for all N and i , and the aforementionedresults also show that the distance of the two sets does not exceed O (1 /N ) , which completes ourproof.So far, we have demonstrated that our method (6) yields over-approximations of reachable setsand tubes, for any uncertainty sets X and U satisfying Hypothesis ( H ), assuming the ability tocompute with compact convex sets in finite dimension. By suitably representing these sets and the setoperations in (6), thereby possibly specializing to a subclass of sets, the method can be implementedon a computer. See e.g. [13], [32] for a discussion of the merits of several classes of sets and theirrepresentations in reachability analysis. Serry, Reissig Computing reachable tubes
B. Zonotopic Over-approximation
In this section, we present a variant of our method (6) for the class of zonotopes, i.e., for sets ofthe form ¯ Z ( c, G ) = c + G [ − , q (12)for some c ∈ R n , G ∈ R n × q , and q ∈ Z + , where c is the center and the columns of G are the generators of ¯ Z ( c, G ) . In particular, we assume that the uncertainty of the system (1) is given as zonotopes, X = ¯ Z ( a, E ) and U = ¯ Z ( c, G ) , (13)where a ∈ R n , c ∈ R m , E ∈ R n × p , G ∈ R m × q , and p, q ∈ Z + .A problem with zonotopic implementations of (6) is that zonotopes are not closed under convexhulls, and so the sets Γ i defined in (6d) are not, in general, zonotopes. We here follow an idea by Girard [9] and replace Γ i by a zonotope obtained using the enclosure operator Enc : ( R n × R n × p ) → R n × R n × (2 p +1) given by Enc(( b, F ) , ( c, G )) = (cid:18) b + c , (cid:18) F + G , b − c , F − G (cid:19)(cid:19) (14)for all b, c ∈ R n , F, G ∈ R n × p and p ∈ Z + . Specifically, we propose the following variant of our method(6) for the case of zonotopic uncertainties (13). Given a time discretization parameter N ∈ N , we shallcompute sequences ( b i ) i ∈ [0; N ] , ( F i ) i ∈ [0; N ] , ( d i ) i ∈ [1; N ] and ( H i ) i ∈ [1; N ] satisfying the following conditionsfor all i ∈ [1; N ] . b = a and F = E, (15a) m i − = k ( b i − , F i − ) k ∞ and K i = hB ( t i ) G, (15b) b i = e φ ( t i , t i − ) b i − + hB ( t i ) c (15c) F i = (cid:16) e φ ( t i , t i − ) F i − , K i , ( α h + θ h m i − ) id (cid:17) , (15d) ( d i , J i ) = Enc (cid:16) ( b i − , F i − ) , ( b i , e φ ( t i , t i − ) F i − ) (cid:17) , (15e) H i = ( J i , K i , ( α h + β h + ( γ h + θ h ) m i − ) id) , (15f)where h , t i , α , β and γ are given by (6a), (6e) and (6f) and the norm k · k in (6e) and (6f) is themaximum norm.We note that the norm in (15b) is straightforward to compute. See Lemma A.3. Moreover, (15a)-(15d) is a straightforward implementation of the set operations in (6b)-(6c) into linear algebraicoperations on centers and generators, and using induction it easily follows that Ω i = ¯ Z ( b i , F i ) for all i ∈ [0; N ] , (16)provided that the norm k · k in (6) is the maximum norm and B is the respective closed unit ball.Thus, by Proposition IV.1, the pairs ( b i , F i ) produced by algorithm (15) represent zonotopic over-approximations of reachable sets R ( t i ) with first order approximation error.The case of reachable tubes is more involved and is the subject of Theorem IV.7 and its CorollaryIV.8 below. We shall demonstrate that the pairs ( d i , H i ) produced by the algorithm (15) representzonotopes ¯ Z ( d i , H i ) over-approximating the sets Γ i defined in (6d). Then, by Theorem IV.4, thesezonotopes over-approximate the reachable tubes R ([ t i − , t i ]) , and we shall also show that first orderconvergence is retained. This way, we obtain a solution to Problem III.1 which applies in the casethat the uncertainty of the system (1) is given as zonotopes, and, in contrast to the more generalalgorithm (6), this solution can be directly implemented on a computer. As before, we shall use thesuperscript N to indicate that, e.g. the sequence ( F Ni ) i ∈ [0; N ] has been computed by our method (15)for a specific value of the time discretization parameter, and similarly for b i , d i and H i . erry, Reissig Computing reachable tubes 9 IV.7 Theorem (Zonotopic Over-approximation of Reachable Tubes).
Assume (13), and foreach N ∈ N , let sequences ( t Ni ) i ∈ [0; N ] , ( d Ni ) i ∈ [1; N ] and ( H Ni ) i ∈ [1; N ] be defined by (6a), (6e), (6f ) and(15), where the norm k · k in (6e) and (6f ) is the maximum norm, and denote Λ Ni = ¯ Z ( d Ni , H Ni ) .Then R ([ t i − , t i ]) ⊆ Λ Ni for all N ∈ N and all i ∈ [1; N ] , and d H ( R ([ t i − , t i ]) , Λ Ni ) ≤ O (1 /N ) as N → ∞ , uniformly w.r.t. i . IV.8 Corollary.
Under the hypotheses and in the notation of Theorem IV.7, denote b R N = ∪ i ∈ [1; N ] Λ Ni .Then R ([ t , t f ]) ⊆ b R N for all N ∈ N , and d H ( R ([ t , t f ]) , b R N ) = O (1 /N ) as N → ∞ . Our proof of Theorem IV.7 uses the following auxiliary result.
IV.9 Lemma.
Let Ω , Γ , W ⊆ R n be non-empty, compact and convex, and suppose that ∈ W . Then conv(Ω ∪ (Γ + W )) ⊆ W + conv(Ω ∪ Γ) , (17) and the Hausdorff distance between the two sets does not exceed k W k .Proof. Let r ∈ conv(Ω ∪ (Γ + W )) . Then, by Lemma A.1(ii), there exist λ ∈ [0 , , x ∈ Ω , y ∈ Γ and z ∈ W such that r = λx + (1 − λ )( y + z ) = λx + (1 − λ ) y + (1 − λ ) z. Notice that (1 − λ ) z ∈ W as , z ∈ W . Hence, r ∈ W + conv(Ω ∪ Γ) which implies (17). Let s ∈ W + conv(Ω ∪ Γ) , then there exist λ ∈ [0 , , x ∈ Ω , y ∈ Γ , z ∈ W such that s = λx + (1 − λ ) y + z .Define t = λx + (1 − λ )( y + z ) ∈ conv(Ω ∪ (Γ + W )) . Then we have s − t = λz , and so k s − t k ≤ k W k ,which proves the bound. Proof of Theorem IV.7.
For each N ∈ N , let h N and sequences (Ω Ni ) i ∈ [0; N ] , (Γ Ni ) i ∈ [1; N ] , ( b Ni ) i ∈ [0; N ] , ( F Ni ) i ∈ [0; N ] , ( J Ni ) i ∈ [1; N ] , ( K Ni ) i ∈ [1; N ] and ( m Ni ) i ∈ [0; N [ be defined by (6) and (15). In the sequel, we dropthe superscript N from our notation, and often we do not mention the domains N and [1; N ] of N and i . Everything is w.r.t. the maximum norm, here denoted by k · k . This applies, in particular, tothe norm k · k and to the unit ball B in (6).We claim that Γ i ⊆ Λ i for all N and i , and that d H (Γ i , Λ i ) ≤ O (1 /N ) , where asymptotic estimatesare always meant to hold for N → ∞ , uniformly w.r.t. i . The theorem then follows from an applicationof Theorem IV.4.To prove the claim, let N ∈ N and i ∈ [1; N ] , and denote P = Ω i − , L = e φ ( t i , t i − ) , w = hB ( t i ) c , M = LP + w , and W = hB ( t i )( U − c ) + ( α h + β h + ( γ h + θ h ) k P k ) B . Then Γ i = conv( P ∪ ( M + W )) by(6c) and (6d), and so Γ i ⊆ W +conv( P ∪ ( LP + w )) by Lemma IV.9, and in turn, Lemma A.3(iv), (15c)and (15e) imply Γ i ⊆ W + ¯ Z ( d i , J i ) . Then W + ¯ Z ( d i , J i ) = Λ i by (15b), (15f) and Lemma A.3(i), whichproves the first part of our claim. From Lemmata IV.9 and A.3(iv) we additionally obtain the bound d H (Γ i , Λ i ) ≤ k W k + k L − id kk F i − k . Next, Lemma A.3(iii) shows that k F i − k ≤ k Ω i − k , and thetriangular inequality, Proposition IV.1, and the estimate (2) yield k L − id k ≤ O (1 /N ) . Finally, by theboundedness of B from assumption ( H ), and by the fact that k Ω i k ≤ O (1) by Proposition IV.1 andthe compactness of reachable tubes, we obtain k W k = O (1 /N ) , which implies d H (Γ i , Λ i ) = O (1 /N ) as claimed.To close this section, we discuss the complexity of the proposed method. It is easily seen that thememory requirement of algorithm (15) is determined by the need to store the computed zonotopicover-approximation. The zonotope Λ Ni , i ∈ [1; N ] , obtained in Theorem IV.7, has p + 1 + (2 i − q + n ) generators, and consequently, the over-approximation b R N obtained in Corollary IV.8 consists of N zonotopes in R n with a total of ( q + n ) N + (2 p + 1) N generators. In particular, the memoryrequired by the zonotopic over-approximation, and in turn, the memory required by our method,is of order O ( n N ) as one of the variables n or N tends to infinity and the other one is fixed,where we have assumed both p = O ( n ) and q = O ( n ) . Regarding the run time of algorithm (15),we additionally assume m = O ( n ) and consider only arithmetic operations. Then the computational Figure 1. Zonotopic over-approximations computed by the proposed method for the example in Section V, for selected values ofthe discretization parameters N and N d : Reachable tube R ([ t , t f ]) for N d = 4 (left), and maximum bridge displacement, uponall nodal points, obtained from reachable tube R ([ t , t f ]) , for N d ∈ { , , , , } (right). effort is dominated by the multiplication of an n × n matrix by an n × ( p + ( i − q + n )) matrix instep (15d). Hence, the run time of algorithm (15) is of order O ( n N ) as one of the variables n or N tends to infinity and the other one is fixed. V. Numerical Example
In this section, we demonstrate the performance of the proposed method on a comprehensivenumerical example. Specifically, we consider several instances of a reduced order model, obtained byfinite difference approximation, of an infinite dimensional pedestrians-footbridge system.
A. Footbridge Model
Analyzing the dynamic response of footbridges is crucial for their structural integrity and thesafety of pedestrians [33], [34]. In this example, we consider a continuum model of a pedestrians-footbridge system. The pedestrians’ stroll on the footbridge generates dynamic load which triggersdeformation of the footbridge. The lateral displacement of the footbridge is given by the function q : I × [0 , L ] → R , where L denotes the length of the bridge, and I = [ t , t f ] is the compact timeinterval under consideration. Let S = I × [0 , L ] . The function q satisfies (see, e.g. [33, eq. 29] [34,eqs. 2.91, 2.92] [35, eq. 2]) mD q + EID q + cD q = f p ( t, y ) , ( t, y ) ∈ S , (18a) f p ( t, y ) = f cos( ωt ) q ( t, y ) + w ( t, y ) , ( t, y ) ∈ S , (18b) | w ( t, y ) | ≤ ¯ w, ( t, y ) ∈ S , (18c) q ( t,
0) = q ( t, L ) = D q ( t,
0) = D q ( t, L ) = 0 , t ∈ I , (18d) q ( t , y ) = D q ( t , y ) = 0 , y ∈ [0 , L ] , (18e)where D ki q denotes the k th order partial derivative of the map q w.r.t. its i th argument, D := D , m denotes the mass per unit length, EI is the bending stiffness, c is a damping coefficient, f p ( t, y ) is erry, Reissig Computing reachable tubes 11 Figure 2. Performance of the proposed method on the example in Section V for the time discretization parameter N = 100 ,depending on the dimension n of the state space of the system (1): Total number of components of the generators of the computedzonotopic over-approximation (left), and run time (right). the load per unit length, f and ω are model parameters, w ( t, y ) is a bounded uncertain term, and ¯ w ∈ R + is the bound on the uncertainty. B. Reduced Order Model
Now, we deduce a reduced order model from system (18) by means of finite difference. Let N d be a spatial discretization parameter with N d ≥ , and h y := L/N d , y i := ih y , i ∈ [0; N d ] . Byreplacing the th order spatial derivative in (18a) with second order centered difference (see, e.g. [36]),and the second order spacial derivatives in (18d) with first order forward and backward differences,respectively, and considering the homogeneous boundary and initial conditions in (18d) and (18e), weobtain the approximating model m ¨ z + K ( t ) z + c ˙ z = v ( t ) , where z (0) = ˙ z (0) = 0 , v ( t ) ∈ [ − ¯ w, ¯ w ] N d − , t ∈ [ t , t f ] , and K ( · ) is obtained from the finite difference approximation. By setting x = ( z, ˙ z ) , wearrive at a problem for system (1), where A ( t ) = (cid:18) − m K ( t ) − cm id (cid:19) and B = (cid:18) (cid:19) , (19a) U = [ − ¯ w/m, ¯ w/m ] N d − , X = { } . (19b)The dimension of the system (1) is n = 2( N d − . We note that system (19), for the case N d = 4 ,corresponds to a damped and perturbed version of the well-known Mathieu equation which modelsvarious physical phenomena and engineering systems; see, e.g. [26]. We also note that the approachfollowed above to obtain the reduced order linear time-varying (LTV) problem (19) from (18) hasbeen followed in the literature to construct benchmark problems for reachability analysis of lineartime-invariant (LTI) systems [37]. In this example, we set L = 10 , m = 2 , c = f = EI = ω = 1 , ¯ w = 0 . , and [ t , t f ] = [0 , . We aim at over-approximating the reachable tube R ([ t , t f ]) of (1)for several instances of N d , using the proposed method. The computed over-approximations will beused to obtain bounds on the bridge displacements and will be analyzed in terms of accuracy andcomputational costs. Figure 3. Comparison between the zonotopic variant, with N = 800 , and the ET both applied to the example in Section V fordifferent instances of the discretization parameter N d : Estimation of maximum bridge displacement at time t f , and run time(right). C. Implementation
To address the problem described above, we employ the zonotopic variant of the proposed methodas given in (15). Moreover, we take advantage of the smoothness of the matrix-valued function A ( · ) in (19a) and use a second order approximation ˜ φ of the transition matrix as given in Lemma A.4.The zonotopic variant is implemented in MATLAB (2019a), and MATLAB is run on an AMD Ryzen5 2500U/ 2GHz processor. Plots of zonotopes are produced with the help of software CORA [38]. D. Results
First, we demonstrate the over-approximations obtained by the proposed method. Fig. 1 illustratesseveral over-approximations of R ([ t , t f ]) , with N d = 4 . As seen in the mentioned figure, the accuracyof the over-approximations increases as the value of N increases, which matches with the findingsof this work. Next, bounds on the displacement of the bridge are obtained based on several over-approximations of R ([ t , t f ]) for several instances of N d . Fig. 1 also illustrates bounds on bridgedisplacement obtained for several values of N d and N . The quality of the bounds improves as N increases due to the increased accuracy of the computed over-approximations.Finally, the scalability of the proposed method is illustrated by considering a fixed value of thetime discretization parameter N and selected values of the spatial discretization parameter N d . Fig. 2indicates a memory requirement of order O ( n ) , as predicted by the discussion at the end of SectionIV, and a run time of order O ( n . ) approximately, which is less than the predicted O ( n ) . Thedifference is due to the fact that MATLAB takes advantage of the sparse structure of the matrix e Φ in Lemma A.4, inherited from the matrix A in (19a). E. Comparison: Ellipsoidal Techniques
In this subsection, we illustrate the performance of the zonotopic variant of the proposed methodin comparison with ellipsoidal techniques [5] implemented in the ellipsoidal toolbox (ET) [39]. erry, Reissig Computing reachable tubes 13
As we have discussed in the Introduction, ellipsoidal techniques yield reachable tubes only in implicitform, and hence we restrict the scope of the comparison to reachable sets only. Here, we considerdifferent instances of system (19) with n ∈ [2; 18] , or equivalently, N d ∈ [4; 12] . The ellipsoidal set ˜ U := ( ¯ w/m ) B , where B is the -norm closed unit ball in R N d − , is considered as the input setwhen applying the ET, as the ET is directly applicable to ellipsoidal initial and input sets only.Since ˜ U ⊆ U , the ET is given the advantage of using a smaller input set. As the sets B ˜ U and X are degenerate, the ET requires defining a regularization parameter, which we have set to be − ,which introduces full dimensional conservative substitutes. Moreover, we arbitrarily use the directionvector e = (1 , , , . . . ) in our computations of ellipsoidal approximations. The zonotopic variant isimplemented with N = 800 and restricted to compute reachable sets only (computations in equations(15e) and (15f) are omitted). Both techniques are set to obtain over-approximations of R ( t f ) whichare subsequently used to estimate the maximum bridge displacement, upon all nodal points, at time t f .Fig. 3 (left) shows that for instances of system (19), with N d ∈ [4; 12] , the zonotopic variantperforms very well in comparison with the ET in terms of estimating maximum bridge displacementdespite the inherent disadvantage of using a larger input set and of accounting for approximationerrors which are not considered by the ET. As seen from Fig. 3 (left), the effect of these errors is morepronounced for increasing state space dimension, which is due to rapidly growing estimates of matrixnorms of the system matrices and their derivatives (growth is of order O ( n ) as a result of the finitedifference approximation of the th order derivative in equation (18a)). Fig. 3 (right) illustrates thatthe zonotopic variant outperforms the ET in terms of computational time for the instances of system(19), with N d ∈ [4; 12] . We note, however, that the ET computes additionally under-approximationsof reachable sets, which might contribute to the relatively higher computational time. VI. Conclusion
We have proposed a method to compute over-approximations of reachable tubes for LTV systemsthat are additionally represented in a form suitable for formal verification purposes. The method hasbeen inspired by existing techniques for LTI systems, and, when applied to that special case, it isalmost equivalent to those in [9], [11], except that it additionally requires to repeatedly compute convexhulls. We have also presented a zonotopic variant of the method and demonstrated its performance onan example, which indicates that the computational effort is comparable to that of existing methodsapproximating reachable sets rather than tubes. The accuracy of our method could be improved byimplementing component-wise estimates as in e.g. [12], in place of matrix and vector norms.
Appendix
A.1 Lemma (Convex Sets).
Let Ω , Γ ⊆ R n be convex, α, β ∈ R , and L ∈ R m × n . Then the followingholds.(i) The sets α Ω , Ω + Γ and L Ω are convex. They are additionally compact if Ω and Γ are so.(ii) If Ω and Γ are additionally non-empty, then conv(Ω ∪ Γ) = { λx + (1 − λ ) y | x ∈ Ω , y ∈ Γ , λ ∈ [0 , } .Proof. For (i), see [40, Sect. 3] and note that images of compact sets under continuous maps arecompact. For (ii), see [40, Th. 3.3].The
Hausdorff distance d H (Ω , Γ) of two non-empty bounded subsets Ω , Γ ⊆ R n w.r.t. k · k is definedby d H (Ω , Γ) = inf { ε > | Ω ⊆ Γ + ε B , Γ ⊆ Ω + ε B } , and is used to measure the extent by which the two sets Ω and Γ differ from each other. This distancesatisfies the triangle inequality, it is a metric when restricted to non-empty compact subsets of R n [41], and it additionally enjoys the properties in the following lemma. A.2 Lemma (Hausdorff Distance).
Let Ω , Ω ′ , Γ , Γ ′ ⊆ R n be non-empty and bounded, and let A, B ∈ R m × n and δ, ε ∈ R + . Then the following holds:(i) d H (Ω + Γ , Ω ′ + Γ ′ ) ≤ d H (Ω , Ω ′ ) + d H (Γ , Γ ′ ) .(ii) d H ( A Ω , A Γ) ≤ k A k d H (Ω , Γ) .(iii) d H ( A Ω , B Ω) ≤ k A − B kk Ω k .(iv) k Ω k = d H (Ω , { } ) .(v) If Ω and Γ are additionally closed, then d H (Ω , Γ) ≤ ε iff Ω ⊆ Γ + ε B and Γ ⊆ Ω + ε B .(vi) If d H (Ω , Γ) ≤ ε , then d H (Ω , Γ + δ B ) ≤ ε + δ .(vii) Let (Ω i ) i ∈ I and (Γ i ) i ∈ I be families of non-empty subsets of R n . Then d H ( ∪ i ∈ I Ω i , ∪ i ∈ I Γ i ) ≤ sup i ∈ I d H (Ω i , Γ i ) .Proof. For (i), (iii) and (v), see [41, Lemma 2.2], [25, Lemma 0.1.2.7], and the discussion in [42,p. 48]. The definition of d H directly implies (ii) and (iv), and (i) and (iv) imply (vi). To prove (vii),we may assume that sup i ∈ I d H (Ω i , Γ i ) < ε for some real ε . Then d H (Ω i , Γ i ) < ε for every i ∈ I , andin turn Ω i ⊆ ε B + Γ i . It follows that ∪ i ∈ I Ω i ⊆ ε B + ∪ i ∈ I Γ i , and similarly with the roles of Ω i and Γ i interchanged. Hence, d H ( ∪ i ∈ I Ω i , ∪ i ∈ I Γ i ) ≤ ε , and since the bound ε was arbitrary, we are done.In the following result, given a norm k · k on R n , it is assumed that the norm of any matrix A ∈ R n × p is w.r.t. the maximum norm on R p , i.e., k A k = sup {k Ax k | k x k ∞ ≤ } . A.3 Lemma (Zonotopes).
Let b, c ∈ R n , F ∈ R n × p , G ∈ R n × q , and L ∈ R m × n , and denote Ω = ¯ Z ( b, F ) , Γ = ¯ Z ( c, G ) . Then the following holds.(i) Ω + Γ = ¯ Z ( b + c, ( F, G )) , where ( F, G ) ∈ R n × ( p + q ) .(ii) L Γ = ¯ Z ( Lc, LG ) .(iii) k Γ k = k ( c, G ) k , where ( c, G ) ∈ R n × ( q +1) . In particular, k Γ k ∞ = max n | c i | + P qj =1 | G i,j | (cid:12)(cid:12)(cid:12) i ∈ [1; n ] o .(iv) If p = q , then conv(Ω ∪ Γ) ⊆ ¯ Z (Enc(( b, F ) , ( c, G ))) , where the operator Enc is defined in (14).Moreover, the Hausdorff distance between the two sets does not exceed k F − G k . In particular, for theHausdorff distance w.r.t. the maximum norm, the bound equals max nP qj =1 | F i,j − G i,j | (cid:12)(cid:12)(cid:12) i ∈ [1; n ] o .Proof. For (i) and (ii), see [43, Prop. 1.4, 1.5]. To prove (iii), note that k ¯ Z (0 , ( c, G )) k = k ( c, G ) k bythe very definition (12) of ¯ Z , and that Γ ⊆ ¯ Z (0 , ( c, G )) . It remains to show that k ¯ Z (0 , ( c, G )) k ≤ k Γ k .To this end, let p = αc + y for some α ∈ [ − , and y ∈ G [ − , q . Then p = λ ( c + y ) + (1 − λ )( y − c ) for λ = (1 + α ) / ∈ [0 , , and so k p k ≤ max {k c + y k , k c − y k} as k · k is convex. It follows that k p k ≤ k Γ k since c + y, c − y ∈ Γ . The set inclusion claim in (iv) is known [9]; we sketch a proof: Denote E = ¯ Z (Enc(( b, F ) , ( c, G ))) and let x = b + F λ for some λ ∈ [ − , p . Then x = ( b + c ) / b − c ) / F + G ) λ/ F − G ) λ/ ∈ E . This shows that Ω ⊆ E , and similarly we obtain Γ ⊆ E . The claimfollows as E is convex. To prove the estimate, which improves Girard’s result [9], let x ∈ E . Then x = b + c + α ( b − c ) + ( F + G ) µ + ( F − G ) ν for some α ∈ [ − , and some µ, ν ∈ [ − , p . Define y = λ ( c + Gµ ) + (1 − λ )( b + F µ ) for λ = (1 − α ) / ∈ [0 , . Then y ∈ conv(Ω ∪ Γ) by Lemma A.1(ii),and x − y = ( F − G )( ν − αµ ) / . Since k ν − αµ k ∞ ≤ , we arrive at k x − y k ≤ k F − G k , which provesthe bound. A.4 Lemma (Taylor’s method of order ). Suppose that hypotheses ( H )-( H ) in Section IIIhold. Additionally assume that A is of class C and that k ¨ A ( t ) k ≤ M ¨ A holds for all t ∈ [ t , t f ] . Define e φ : D → R n × n by e φ ( t, s ) = id +( t − s ) A ( s ) + ( t − s ) ( ˙ A ( s ) + A ( s ) ) / . Then condition ( H ) in SectionIII holds with θ given by θ ( h ) = (1 + 3 M ˙ A /M A + M ¨ A /M A )(exp( hM A ) − h M A / − hM A − .Proof. The map θ is clearly monotonically increasing, and θ ( h ) = O ( h ) as h → , which implies (5).Given s ∈ [ t , t f ] , Taylor’s formula for φ ( · , s ) about the point s reads φ ( t, s ) = e φ ( t, s ) + ( t − s ) R (1 − z ) D φ ( s + z ( t − s ) , s ) dz . Using the identity D φ ( t, s ) = ( ¨ A ( t ) + 2 A ( t ) ˙ A ( t ) + ˙ A ( t ) A ( t ) + A ( t ) ) φ ( t, s ) as well as the estimate (2) we obtain the estimate (4). erry, Reissig Computing reachable tubes 15 References [1] G. Basile and G. Marro,
Controlled and conditioned invariants in linear system theory . Englewood Cliffs, NJ: PrenticeHall Inc., 1992.[2] F. Blanchini and S. Miani,
Set-theoretic methods in control , ser. Systems & Control: Foundations & Applications. Boston,MA: Birkhäuser Boston Inc., 2008.[3] P. Tabuada,
Verification and control of hybrid systems . Springer, 2009.[4] M. Althoff, S. Bak, M. Forets, G. Frehse, N. Kochdumper, R. Ray, C. Schilling, and S. Schupp, “ARCH-COMP19 CategoryReport: Continuous and hybrid systems with linear continuous dynamics,” in
Proc. 6th Intl. Worksh. Appl. Verificationfor Continuous and Hybrid Systems (ARCH) , ser. EPiC Series in Computing, G. Frehse and M. Althoff, Eds., vol. 61.EasyChair, 2019, pp. 14–40.[5] A. B. Kurzhanski and P. Varaiya,
Dynamics and control of trajectory tubes , ser. Systems & Control: Foundations &Applications. Birkhäuser/Springer, 2014, theory and computation.[6] G. Reissig, A. Weber, and M. Rungger, “Feedback refinement relations for the synthesis of symbolic controllers,”
IEEETrans. Automat. Control , vol. 62, no. 4, pp. 1781–1796, Apr. 2017. doi:10.1109/TAC.2016.2593947[7] M. Althoff, “Reachability analysis of large linear systems with uncertain inputs in the Krylov subspace,”
IEEE Trans.Automat. Control , vol. 65, no. 2, pp. 477–492, 2020.[8] T. Pecsvaradi and K. S. Narendra, “Reachable sets for linear dynamical systems,”
Inform. and Control , vol. 19, pp. 319–344,1971.[9] A. Girard, “Reachability of uncertain linear systems using zonotopes,” in
Proc. 8th Intl. Workshop on Hybrid Systems:Computation and Control (HSCC), Zürich, Switzerland, Mar. 9-11, 2005 , ser. Lect. Notes Comput. Sci., M. Morari andL. Thiele, Eds., vol. 3414. Springer, 2005, pp. 291–305.[10] A. Girard, C. Le Guernic, and O. Maler, “Efficient computation of reachable sets of linear time-invariant systems withinputs,” in
Proc. 9th Intl. Workshop on Hybrid Systems: Computation and Control (HSCC), Santa Barbara, CA, U.S.A.,Mar. 29-31, 2006 , ser. Lect. Notes Comput. Sci., J. Hespanha and A. Tiwari, Eds. Berlin: Springer, 2006, vol. 3927, pp.257–271.[11] C. Le Guernic and A. Girard, “Reachability analysis of linear systems using support functions,”
Nonlinear Anal. HybridSyst. , vol. 4, no. 2, pp. 250–262, 2010.[12] G. Frehse, C. Le Guernic, A. Donzé, S. Cotton, R. Ray, O. Lebeltel, R. Ripado, A. Girard, T. Dang, and O. Maler, “SpaceEx:scalable verification of hybrid systems,” in
Proc. 23rd Intl. Conf. Computer Aided Verification (CAV), Snowbird, UT, USA,Jul.14-20, 2011 , ser. Lect. Notes Comput. Sci., vol. 6806. Springer, Heidelberg, 2011, pp. 379–395.[13] M. Althoff and G. Frehse, “Combining zonotopes and support functions for efficient reachability analysis of linear systems,”in
Proc. IEEE Conf. Decision and Control (CDC), Las Vegas, U.S.A., 12-14 Dec. 2016 , 2016, pp. 7439–7446.[14] M. Althoff, B. H. Krogh, and O. Stursberg, “Analyzing reachability of linear dynamic systems with parametric uncertainties,”in
Modeling, design, and simulation of systems with uncertainties , ser. Math. Eng. Springer, Heidelberg, 2011, pp. 69–94.[15] O. Botchkarev and S. Tripakis, “Verification of hybrid systems with linear differential inclusions using ellipsoidalapproximations,” in
Proc. 3rd Intl. Workshop on Hybrid Systems: Computation and Control (HSCC), Pittsburgh, U.S.A.,Mar. 23-25, 2000 , ser. Lect. Notes Comput. Sci., N. A. Lynch and B. H. Krogh, Eds., vol. 1790. Springer, 2000, pp. 73–88.[16] M. E. Villanueva, R. Quirynen, M. Diehl, B. Chachuat, and B. Houska, “Robust MPC via min-max differential inequalities,”
Automatica J. IFAC , vol. 77, pp. 311–321, 2017.[17] M. Serry and G. Reissig, “Hyper-rectangular over-approximations of reachable sets for linear uncertain systems,”in
Proc. 57th IEEE Conf. Decision and Control (CDC), Miami, FL, USA, 17-19 Dec. 2018 , 2018, pp. 6275–6282.doi:10.1109/CDC.2018.8619276[18] K. Shen and J. K. Scott, “Exploiting nonlinear invariants and path constraints to achieve tighter reachable set enclosuresusing differential inequalities,”
Math. Control Signals Systems , vol. 32, no. 1, pp. 101–127, 2020.[19] M. Zamani, G. Pola, M. Mazo, Jr., and P. Tabuada, “Symbolic models for nonlinear control systems without stabilityassumptions,”
IEEE Trans. Automat. Control , vol. 57, no. 7, pp. 1804–1809, 2012.[20] N. S. Nedialkov, K. R. Jackson, and G. F. Corliss, “Validated solutions of initial value problems for ordinary differentialequations,”
Appl. Math. Comput. , vol. 105, no. 1, pp. 21–68, 1999.[21] V. Veliov, “Second-order discrete approximation to linear differential inclusions,”
SIAM Journal on Numererical Analysis ,vol. 29, no. 2, pp. 439–451, 1992.[22] W.-J. Beyn and J. Rieger, “Numerical fixed grid methods for differential inclusions,”
Computing , vol. 81, no. 1, pp. 91–106,2007.[23] G. Reissig and M. Rungger, “Symbolic optimal control,”
IEEE Trans. Automat. Control , vol. 64, no. 6, pp. 2224–2239, Jun.2019. doi:10.1109/TAC.2018.2863178[24] S. Bak and P. S. Duggirala, “Simulation-equivalent reachability of large linear systems with inputs,” in
Computer aidedverification. Part I , ser. Lect. Notes Comput. Sci. Springer, 2017, vol. 10426, pp. 401–420.[25] R. Baier, “Set-valued integration and the discrete approximation of reachable sets,” doctoral thesis, Univ. Bayreuth,Bayreuth, Germany, 1995, (“Mengenwertige Integration und die diskrete Approximation erreichbarer Mengen”, in German).[26] T. I. Fossen and H. Nijmeijer, Eds.,
Parametric resonance in dynamical systems . Springer, New York, 2012, proc. Workshopheld in Longyearbyen, Jun. 22–26, 2011.[27] N. G. de Bruijn,
Asymptotic methods in analysis , 3rd ed. Dover, 1981, corr. reprint.[28] D. L. Lukes,
Differential equations , ser. Mathematics in Science and Engineering. London: Academic Press Inc., 1982, vol.162.[29] H. Hermes, “The generalized differential equation ˙ x ∈ R ( t, x ) ,” Advances in Math. , vol. 4, pp. 149–169, 1970.[30] F. L. Chernous ′ ko, State Estimation of Dynamic Systems . CRC Press, 1994. [31] P. Deuflhard and F. Bornemann,
Scientific computing with ordinary differential equations , ser. Texts in Applied Mathematics.Springer-Verlag, New York, 2002, vol. 42.[32] C. Le Guernic, “Reachability analysis of hybrid systems with linear continuous dynamics,” Ph.D. dissertation, UniversitéGrenoble I - Joseph Fourier, 28 Oct. 2009.[33] G. Piccardo and F. Tubino, “Parametric resonance of flexible footbridges under crowd-induced lateral excitation,”
J. SoundVibration , vol. 311, no. 1-2, pp. 353–371, 2008.[34] F. Gazzola,
Mathematical models for suspension bridges . Springer, 2015.[35] J. Bodgi, S. Erlicher, and P. Argoul, “Lateral vibration of footbridges under crowd-loading: Continuous crowd modelingapproach,” in
Key Engineering Materials , vol. 347. Trans Tech Publ, 2007, pp. 685–690.[36] B. Fornberg, “Generation of finite difference formulas on arbitrarily spaced grids,”
Math. Comp. , vol. 51, no. 184, pp.699–706, 1988.[37] H.-D. Tran, L. V. Nguyen, and T. T. Johnson, “Large-scale linear systems from order-reduction,” in
Proc. 3rd Intl.Worksh. Appl. Verification for Continuous and Hybrid Systems (ARCH) , ser. EPiC Series in Computing, G. Frehse andM. Althoff, Eds., vol. 43. EasyChair, 2017, pp. 60–67.[38] M. Althoff, “An introduction to CORA 2015,” in
Proc. 1st and 2nd Intl. Worksh. Appl. Verification for Continuous andHybrid Systems (ARCH) , ser. EPiC Series in Computing, G. Frehse and M. Althoff, Eds., vol. 34. EasyChair, 2015, pp.120–151.[39] P. Gagarinov and A. A. Kurzhanskiy, “Ellipsoidal toolbox release 2.0.1,” UC Berkeley, Dept EECS, Tech. Rep., May 2014.http://systemanalysisdpt-cmc-msu.github.io/ellipsoids/[40] R. T. Rockafellar,
Convex analysis , ser. Princeton Mathematical Series, No. 28. Princeton, N.J., U.S.A.: PrincetonUniversity Press, 1970.[41] F. S. De Blasi, “On the differentiability of multifunctions,”
Pacific J. Math. , vol. 66, no. 1, pp. 67–81, 1976.[42] R. Schneider,
Convex bodies: the Brunn-Minkowski theory , ser. Encyclopedia of Mathematics and its Applications.Cambridge University Press, Cambridge, 1993, vol. 44.[43] V. T. H. Le, C. Stoica, T. Alamo, E. F. Camacho, and D. Dumur,