Reachability Analysis of Large Linear Systems with Uncertain Inputs in the Krylov Subspace

1 Reachabili ty Analysis of Lar ge Linear Systems with Uncertain Inp u ts in the Krylo v Subsp ace Matthias Althof f Abstract —One often wishes f or the ability to form ally analyze large-sca le systems—typically , ho wev er , one can e ither f ormally analyze a rather small system or informally analyze a large-sca le system. This work tries to further close this perf ormance gap fo r reachability analysis of li near systems. Reachability analysis can capture the whole set of possible solutions of a dynamic system and is thus used to pr ov e that unsafe states ar e neve r reached; this r equires full consideration of arbitrarily v arying uncertain in p uts, since sensor n oise or dist u rbances usually do not f ollow any patter ns. W e use Krylo v methods in this w ork to compute reachable sets for large-sca le li n ear systems. While Krylov methods hav e been used bef ore in r eachability analysis, we ov erco me the previous l i mitation that inp uts must be (p iecewise) constant. As a result, we can compute rea chable sets of systems with se vera l thousand state variables for b ounded, but arbitrarily vary ing input s. Index T erms —Reachability analysis, linear systems, Krylov subspace, u ncertain inp uts, large-scale systems. I . I N T RO D U C T I O N Reachability analysis compu tes th e set of possible solutions of dynam ic systems subject to uncer tain initial states and inputs. The av ailab ility of all possible solutions can be u sed for many purpo ses: form al veriﬁcation o f dyn amic sy stems with discrete a n d/or co ntinuou s dynamics [ 1], [2]; c omputatio n o f in variance sets [3], [4]; compu tation o f the region of attraction [5]; optimizatio n of constraine d systems with uncertain ties [6], [7]; set-based observers [8], [9] ; and con f ormance check- ing [10]. The theory of efﬁciently co mputing reachable sets is adv ancing rap idly an d the ir usefulness has already been demonstra te d for many application s, such a s automa te d driving [11], [12], robotics [13], [14], p ower systems [15], [16], and analog/mixed- signal c ir cuits [17], [18]. a) Reachab ility analysis of linear systems: The reach- ability analysis in par ticular of linear continu ous systems has been in tensely r e sear ched. While the main ap proach for computing r eachable sets o f linear systems by using the superp osition principle has no t signiﬁcantly ch anged in recent year s [19, Eq . 4.1 4], much p rogress has b een made by experimen ting with different set representations: ellipsoids [20], polytopes [21], zonotopes [22], zonotope bundles [2 3], support f unctions [ 2 4], lev el sets [25], and comb ination of support functions and zonotop es [2 6]. In particular, when using zonoto p es, suppo r t function s, or the comb ination thereof , one can efﬁciently com pute systems with several hun dred continuo us state variables. Recently , a new tech n ique ha s been pr oposed which comb ines simulation results by using Matthi as Althof f is wit h t he De partment of Computer Science, T echnisch e Univ ersit ¨ at M ¨ unchen, 85748 Garching, Germany , email: althoff@tum.d e the superp osition pr inciple to represen t reac hable states via generalized star sets for certain inp uts [27]– [29] an d un certain, piecewise-constant inpu ts [30]. Although this techniq ue can compute very large sy stems—in some ca ses up to a billion states variables [ 29]—it can not consider unce rtain, ar bitrarily time-varying in puts and requires a for mally veriﬁed solver for linear systems. Even mo r e recently , de c o mposition techn iques have bee n developed to speed up reach ability analysis of linear systems [31]; howe ver, arbitrarily -varying inputs can not be considered in [ 3 1]. Since o ne cann ot exactly compute the reachable set o f lin ear systems, excep t wh en all eigenv alues are real or im aginary [3 2], one typ ically deman ds tight over- approx imations of reach able sets. b) F urther use for nonlin ear systems: Reachable set computatio ns of linear system s are often embed ded in al- gorithms for compu tin g reachable sets of nonlinear systems, either by partitio ning the state space into conservati vely lin- earized regions [ 3 3] or b y conservativ ely and continu ously linearizing a non linear system alo ng its cen ter traje c tory [34], [35]. W ith conservative linearization, we refer to techniques that comp ensate linearization er rors by ad ding unc e rtainty , e.g., in the for m of additive un certain in puts. c) Or d er r ed uction techniques: For large-scale linear systems with u ncertain time-varying in puts beyond 1000 continuo us state variables, howe ver, e ven the most efﬁcient reachability a lgorithms becom e too slow f or practical use. On e of the main reaso ns is that the exponen tial m a trix e At of the system ma trix A in the linea r system dyn amics ˙ x = Ax + B u may be unbearab ly time- consumin g to co mpute [36]. I t sho uld be n oted that the dimension is not the on ly critical param eter; sparsity o f A and the sensitivity of the matrix exponential [36, Sec. 2] are also impo rtant p a rameters inﬂuencin g the computatio n time. For this r eason, method s have been devel- oped to reduce the o rder o f the in vestigated system, w h ich are generally r eferred to as or der redu ction techn iques. Ther e exist a v ast nu mber of different techniqu es surveyed in e.g., [37]–[39]. Mo st or der red uction techn iques aim at achieving a similar inp ut/outpu t b ehavior when the sy stem is initially in a steady state; they r ely on th e fact that large- scale systems often h av e a large numb er of state variables but a r ather small num ber o f inp ut and ou tp ut variables. A typical example is that of inﬁnite-dim ensional systems, which are spatially discretized, c ontrolled by few actu ators, and sensed by f ew sensors. M any application s of reachab ility analysis, howe ver, such as forma l veriﬁcation , require that all or m a ny state variables are accurately app roximated as well. For instance, in a large power ne twork , all voltages and f r equencies have to stay within certain bounds. Subsequ ently , we revie w p revious work o n combinin g reachab ility an a lysis with ord er r eduction 2 technique s. d) Or der r ed uction for reac hability analysis: T o th e best knowledge of the author, the ﬁr st work combinin g reacha b il- ity analysis with order reduction techniq ues is [40]. Ther e, simulating th e solutions of all vertices o f an in itial set was required to guar a n tee an erro r boun d f o r a set of initial states— this appr o ach is e xponen tial in the nu mber of state v ariables when each variable is un certain w ith in an interval. At the time the approach in [40] was p r oposed, too ls f or reach ability analysis had an exponen tial c o mplexity a s well, so that overall the com putation tim e c ould be signiﬁcantly reduced . Howev er , modern too ls such as Spa c eEx [41], Flow* [4 2], HyLAA [28], XSpeed [43], or CORA [44] ha ve a polynomial com p lexity , as d emonstrated in [45], [46], and thus would m ost cer tain ly outperf orm the techniqu e p roposed in [40], even with out any order red u ction. The same author s later combine d reach a b ility analysis with Kr ylov subspace appro ximation methods [47]. The advantage of this tech n ique is that it does not scale exponentially in the nu mber of state variables; h owe ver, it can only handle linear systems with ﬁxed inp u t, wh ich are also referred to as af ﬁne systems. Recently , the work from [40] w as co ntinued in [48], which conside rs order reduction in th e input/ou tput sense for stable linear system s. In contrast, this work c an also handle unstable systems and r e construct the whole set of states, not just th e o utputs. For n onlinear systems, non-r igorou s redu ction techn iques h av e been pr esented in [4 9]; howe ver, un like in this work, the results are not formal. e) Appr oximate bisimulation: Somewhat related to ord er reduction tech n iques is (ap proxim ate) bisimulatio n [50], [5 1], which basically sho ws that t wo systems ha ve similar ou tput behavior for the same in puts, wh ile at the same tim e a relation S exists so that ( x ( t ) , ˜ x ( t )) ∈ S . Howe ver , even bisimulatio n technique s for linear systems do n ot scale to the system dimensions handled in this work, since they r equire solving linear ma tr ix ineq u alities [52] to ﬁnd a bisimulation relation [53]. f) Contribution: This work e xtends the state of the art in reachab ility analysis of linear systems by prop osing the ﬁrst metho d with unc e r tain, tim e - varying inputs in the reduce d Krylov subspace. This red uction makes it possible to c ompute reachable sets of systems with a number of co ntinuou s state variables that was previously infeasible. Please note that in- put/outp ut or d er re duction techniqu e s such as transform ation, truncation , and projec tion (com bination of transfo rmation and truncation ) cann o t recon struct the state, in c ontrast to the reduction p resented in this work. Howe ver , we also improve computatio nal efﬁciency wh e n only the in put/outp ut b e h avior is of interest. Ou r pr oposed tech nique rigor o usly consid ers reduction err ors. I n contrast to [47], ou r appro ach a) does not rely on computin g errors by the norm of the sy stem matrix k A k (see (8 )), which would q uickly accumulate error s wh en k A k is gre a ter than 1 (m ostly true in p ractice; fo r the presented example o f a bridg e k A k = 4 . 6 347 · 10 8 ), b) does not requ ire to enlarging the reach able set equally in all dim ensions to accoun t for erro rs, which can cau se large over-approximation s, and c) can be applied to zo notop e s, which use a compact gen erator representatio n [54, Th eorem 7]. Our a pproach is imp le m ented in CORA [44] and will be pub licly released with th e new CORA version. g) Organization: Th e paper is organized as follows: Sec. II presents pre liminaries fro m the are a s of Kry lov subspace approxim a tion, set representation, and r eachability analysis. The c o mputation of re a chable sets in th e Kr ylov subspace is p resented step by step: The hom ogeneo u s solution is described in Sec. II I and the input solution in Sec. I V. Then, these are c ombined in Sec. V to d emonstrate the overall algorithm . W e close with nu merical exam p les in Sec. VI and the conclusions in Sec. VII. I I . P R E L I M I N A R I E S Let us ﬁrst recall some important b asics required for this work: Krylov subspace appr oximation , rep resentation o f con- tinuous s ets in high-d imensional spaces, and com p utation of reachable sets of linear systems with unce rtain in p uts. A. Krylov Subspace Appr o ximation The main obstacle towards reachab ility analysis o f large- scale linear sy stems is the ev aluation of e C v , wh ere C ∈ R n × n and v ∈ R n . T o co mpute e C v more efﬁciently , we intr o duce the Krylov sub space K ξ = span ( v, C v, . . . , C ξ − 1 v ) , where span ( · ) r eturns the lin ear span o f a set of vectors and ξ d enotes the dimension of the subspace. Several po ssibilities have b een d ev eloped for ap proxim ating e C v in th e Krylov subspace [37], [55], [56 ]. In this work, we use the simplest approa c h , wh ich is also one of the m ost popular: the Arnoldi algorithm , as presented in Alg. 1 ( see [55, Sec. 2 . 1] an d [57, Alg. 1] ). Please n o te that w ∗ in Alg. 1 d enotes the co mplex conjuga te of w and k . k return s the Eu clidean norm. A fur ther reason for ch o osing the Arno ldi algorithm is tha t tight a- posteriori and a-prior error boun ds exist [58], [59]. Algorithm 1 Arnoldi iteration Require: C , v , m a x order ξ , tolerance tol Ensure: H , V 1: v (1) = v / k v k 2: for k = 1 . . . ξ do 3: w = C v ( k ) 4: for j = 1 . . . k do 5: h j,k = w ∗ v ( j ) 6: w := w − h j,k v ( j ) 7: end for 8: h k +1 ,k = k w k 9: if h k +1 ,k ≤ tol k C k then 10: happy-breakdown 11: end if 12: v ( k +1) = w /h k +1 ,k 13: end for 14: V = [ v (1) , v (2) , . . . , v ( ξ ) ] The results of th e Arnoldi iteration (usin g the stabi- lized Gram-Sch midt process) are an ortho gonal basis V = [ v (1) , v (2) , . . . , v ( ξ ) ] of th e Krylov subspace K ξ and the u p per 3 ξ × ξ He ssen b erg matrix H consisting of the elements h ij from Alg. 1. Please note that we are nu mbering vectors with superscripted numbers in parentheses in order to a void c o n- fusion with powers. Using H , V , and e 1 = [1 , 0 , 0 , . . . , 0] T , the ev aluation of e C v can be app roximated as [55, eq. 3-4]: e C v ≈k v k V e H e 1 , (1) e C t v ≈k v k V e H t e 1 . (2) The following rigoro us a-p riori bound for th e appr o ximation error exists [ 5 5, eq. 29]:   e C t v − k v k V e H t e 1   ≤ 2 k v k k C t k ξ e k C t k ξ ! . A tighter a-pr iori bo u nd is proposed in [59] fo r d ifferent types of matrices (ske w-Hermitian matrices, p ositiv e deﬁnite matrices, etc.). Even tighter r esults are obtain ed by usin g a-posterior i b o und. Although they are computationally more expensiv e, they make it possible to obtain error b ounds for long p rediction h orizons s o th at o ne typ ic a lly does no t have to recom p ute the Arn oldi iter ation; this discussion is detailed in Sec. V -B. In th is work , we u se the a-p osteriori bo und f rom [58, Eq. 4.1 ] so tha t ∀ t ∈ [0 , t f ] :   e C t v − k v k V e H t e 1   ≤ k v k ǫ norm t, (3) where the comp utation of ǫ norm is d escribed in detail in Append ix B. B. Set Repr esentation by Zonoto pes As alread y summarized in the in troductio n, se veral set representatio ns for reachab ility ana ly sis of lin ear systems hav e been propo sed. When considering u ncertain in puts, zon o topes [22] a n d su pport function s [ 24] d emonstrate g ood p erfor- mance. Recently , [2 6] sh owed th at combinin g zonotope s and support fun ctions provides e ven b etter beneﬁts: zo notopes can compute the solution mo re efﬁciently , while su pport fun ctions can r e p resent m ore g eneral initial sets. Since most initial sets are multid imensional intervals (special c a se of zonoto pes), we use zo notopes an d neglect their co mbination with supp o rt function s to fo c us o n the novel aspects of this work. Deﬁnition II.1 (Zonotope) Give n a c e n ter c ∈ R n and so- called generators g ( i ) ∈ R n , a zonotop e is deﬁn ed a s Z := n x ∈ R n    x = c + p X i =1 β i g ( i ) , β i ∈ [ − 1 , 1] o . W e write in sho rt Z = ( c, g (1) , . . . , g ( p ) ) and deﬁne the order of a zon otope as o := p n , where p is the numb er of generato r s. A z o notop e can be seen as the Minkowski addition of line segments [ − 1 , 1] g ( i ) , which pr ovid es an idea of how a zonotop e is con structed; see Fig. 1. The iterative compu ta tio n o f r eachable sets for linear sys- tems req u ires set-b ased ad dition, which is often r eferred to as Minkowski ad dition ( X ⊕ Y := { x + y | x ∈ X , y ∈ Y } ), and set-based multiplication ( X ⊗ Y := { x y | x ∈ X , y ∈ Y } ). No te that t he symbol fo r set-based mu ltiplication is often omitted for simplicity of notation , an d that one or bo th op erands can be singletons. The m u ltiplication with a m atrix M ∈ R o × n and P S f r a g r e p l a c e m e n t s c g (1) g (2) g (3) c ⊕ g (1) c ⊕ g (1) ⊕ g (2) c ⊕ g (1) ⊕ g (2) ⊕ g (3) construct ion direct ion ” ⊕ ” ” ⊕ ” Fig. 1. Step-by-step construct ion of a two-dimension al zonotope . the Minkowski addition of two zono topes Z 1 = ( c , g (1) , . . . , g ( p 1 ) ) and Z 2 = ( d , h (1) , . . . , h ( p 2 ) ) are a direct conseq u ence of the zono tope deﬁn itio n (see [ 60]): Z 1 ⊕ Z 2 = ( c + d, g (1) , . . . , g ( p 1 ) , h (1) , . . . , h ( p 2 ) ) , M ⊗ Z 1 = ( M c, M g (1) , . . . , M g ( p 1 ) ) . (4) Also, th e conve x hull o f Z 1 and Z 2 (both h aving equal order) is required (see [22]): conv ( Z 1 , Z 2 ) ⊆ 1 2 ( c + d, g (1) + h (1) , . . . , g ( p 1 ) + h ( p 1 ) , c − d, g (1) − h (1) , . . . , g ( p 1 ) − h ( p 1 ) ) . (5) For th e multiplicatio n o f an interval m atrix M with a zo no- tope, th e matrix M is split into a real-valued matrix M ∈ R n × n and an inter val matrix with bou nd S ∈ R n × n , such that M = M ⊕ [ − S, S ] . After intr oducing S j as the j th row of S , the result is over - approx imated as shown in [61, Theorem 3.3] by MZ 1 ⊆ ( M Z 1 ⊕ [ − S , S ] Z 1 ) ⊆ ( M c 1 , M g (1) , . . . , M g ( p 1 ) , h (1) , . . . , h ( n ) ) , h ( i ) j = ( S j ( | c | + P p 1 k =1 | g | ( k ) ) , fo r i = j 0 , for i 6 = j . Another importan t oper ation is the en closure o f a zonoto pe by an axis-aligned box (see [61, Prop osition 2 .2]): box ( Z ) :=( c, ˆ h (1) , . . . , ˆ h ( n ) ) , ˆ h ( i ) j = ( ( P p 1 k =1 | g | ( k ) ) j , f or i = j 0 , fo r i 6 = j . (6) W e also require some ne w o perations on zono topes in th e Krylov sub space, wh ich are introd uced later in Sec. III. C. R eachability A nalysis Reachable set computations are typically performed iter a - ti vely f or sh o rt time intervals τ k := [ t k , t k +1 ] . (7) In this work, constant-size tim e intervals t k := k δ are used to focus on the m a in inn ovations, where k ∈ N is the time step and δ ∈ R + is refer red to as the time in c rement. An extension to variable time increm ents is d escribed in [41]. 4 W e rec apitulate the reac h ability analy sis of a linear differ- ential inclusion ˙ x ∈ Ax ( t ) ⊕ U , (8) where x ∈ R n , A ∈ R n × n , an d U ⊂ R n is a set of uncertain inpu ts. Please note th a t this also includ es the form ˙ x = A x ( t ) + B u ( t ) , u ( t ) ∈ ˜ U often used in control theory , since on e could e a sily choo se U = B ˜ U . Let χ ( t ; x 0 , u ( · )) denote the solution to (8) f or an initial state x (0) = x 0 and the inpu t trajector y u ( · ) . For a set of initial states X 0 ⊂ R n and a set o f po ssible inpu t values U ⊂ R m , the set of r eachable states is R e ([0 , t f ]) := n χ ( t ; x 0 , u ( · ))    x 0 ∈ R (0) , t ∈ [0 , t f ] , ∀ τ ∈ [0 , t ] u ( τ ) ∈ U o . (9) The sup erscript e o n R e ([0 , t f ]) den otes the exact r e achable set, which can not be compu ted for g eneral lin ear systems [6 2]. Therefo re, an over -ap proxim a tion R ([0 , t f ]) ⊇ R e ([0 , t f ]) is computed as accurately as possible, while at the same time ensuring that the comp utations are efﬁcient and scale well with the system dimension n . For lin ear systems, the main task is to co m pute the r eachable set of the ﬁrst time interval [0 , δ ] . Most of this pap er will deal with computin g the reach able set for the initial tim e interval since the p ropagatio n for la ter time in tervals is rather simple, as shown later in Alg. 2. W e take adv antage of the superpo sition princip le of lin e ar systems by computin g the following reachab le sets separ ately and later joinin g them together: the r eachable set of the h omogen eous solution R h ( t ) ( u ( τ ) = 0 in (9)), th e reachab le set of the pa r ticular solutio n R p ( t ) due to the uncer tain inp ut U ( x 0 = 0 in (9)), an d the reac h able set R ǫ correcting the initial assumption that trajectories are straight lines w ith in [0 , δ ] . According to [61, Sec. 3.2], the reachable set fo r [0 , δ ] is computed as shown in Fig. 2: 1) Starting f r om X 0 , compute the set of all h o mogen eous solutions R h ( δ ) . 2) Obtain the conv ex hu ll o f X 0 and R h ( δ ) to app roximate the reachable set for the time inte r val [0 , δ ] . 3) Compu te R ([0 , δ ]) := S t ∈ [0 ,δ ] R ( t ) by considering uncertain inputs by add ing R p ( δ ) an d acco unting for the curvature of trajecto ries by adding R ǫ . P S f r a g r e p l a c e m e n t s R h ( δ ) X 0 con vex hull of X 0 , R h ( δ ) R ([0 , δ ]) ➀ ➁ ➂ enlar ge- ment Fig. 2. Steps for computing the reachable set of a linea r system for the ﬁrst time interv al. I I I . H O M O G E N E O U S S O L U T I O N I N T H E K RY L OV S U B S P AC E In this section, the b asic idea of compu ting reach able sets as pr e sented in Sec. I I-C is exten d ed so that r eachable sets can be computed in th e Krylov subspa ce. W e ﬁrst present new technique s fo r the homog eneous solu tion of points in time and time inter vals. Subsequ ently , we consider f or the ﬁrst time how reachable sets can be computed f or ar bitrarily ch a nging inpu t trajectories within the Krylov su b space. A. Solution for a P oint in T ime The we ll-known homogen e ous so lution of a lin ear time- in variant system with initial state x 0 is x h ( t ) = e At x 0 . W e can bound the exact solution using the lemma b elow . For that lemma and subsequent derivations we intro duce [ − 1 , 1 ] n as an n - d imensional vector whose entrie s are intervals [ − 1 , 1] . Analogou sly , we write [ − 1 , 1 ] n × m to represent an n × m matrix whose entries are intervals [ − 1 , 1 ] . Lemma III.1 (Single sta te homo geneous so lut io n) After obtainin g V and H fr om Alg. 1 with in puts C = A , v = x 0 , we can boun d the ho m ogeneo us solution x h ( t ) = e At x 0 by x h ( t ) ∈ k x 0 k V e H t e 1 ⊕ [ − 1 , 1 ] n k x 0 k ǫ norm t = ˆ x h ( t ) ⊕ E red ( t, x 0 ) , (10) wher e ˆ x h ( t ) := k x 0 k V e H t e 1 E red ( t, x 0 ) = [ − 1 , 1 ] n k x 0 k ǫ norm t (11) and ǫ norm is computed as described in Ap pendix B. Pr oof: The lem ma directly follows from enlarging the approx imate so lu tion in (1) b y the erro r k x 0 k ǫ norm t f rom (3 ). Since we ha ve for a vector a ∈ R n and a scalar b ∈ R that k a k < b ⇒ a ∈ [ − 1 , 1 ] n b , one obtains from (3) that e A t x 0 − k x 0 k V e H t e 1 ∈ [ − 1 , 1 ] n k x 0 k ǫ norm t, which proves the lemma. For reachability analysis, on e has to comp ute the homo ge- neous solu tio n for a set of in itial states. Replacing the single initial state in Lemma III .1 by a set of initial states X 0 is not tr ivial since the ma tr ices V an d H depend on each initial state x 0 ∈ X 0 . T o indic a te this dep endency , we write V ( x 0 ) and H ( x 0 ) from now on to stress th at tho se ma tr ices have been ob tained fro m the state x 0 . Since X 0 is re p resented as a zonotop e in th is work , the homog eneous r eachable set can be computed by the following th eorem: Theorem III.1 (Homoge neo us solution for a point in time) The r eachable set of th e homogeneou s so lu tion R h ( t ) := e A t X 0 for the initial zonotope X 0 = ( c, g (1) , . . . , g ( p ) ) can 5 be over -appr o ximated by th e zonoto p e R h ( t ) ⊆ ( ˆ c, ˆ g (1) , . . . , ˆ g ( p ) ) ⊕ R h,er r , ˆ c = k c k V ( c ) e H ( c ) t e 1 , ˆ g ( i ) = k g ( i ) k V ( g ( i ) ) e H ( g ( i ) ) t e 1 , R h,er r = [ − 1 , 1 ] n  k c k + p X i =1 k g ( i ) k  ǫ norm t. Pr oof: Inserting th e d eﬁnition o f a zonotop e (Def. II.1) into R h ( t ) = e A t X 0 and using (4), we obtain R h ( t ) = e A t c ⊕ p M i =1 [ − 1 , 1] e A t g ( i ) . Using Lemma III.1 yields R h ( t ) ⊆k c k V ( c ) e H ( c ) t e 1 ⊕ [ − 1 , 1 ] n k c k ǫ norm t ⊕ p M i =1  [ − 1 , 1] k g ( i ) k V ( g ( i ) ) e H ( g ( i ) ) t e 1 ⊕ [ − 1 , 1 ] n k g ( i ) k ǫ norm t  = k c k V ( c ) e H ( c ) t e 1 ⊕ p M i =1 [ − 1 , 1] k g ( i ) k V ( g ( i ) ) e H ( g ( i ) ) t e 1 ⊕ [ − 1 , 1 ] n  k c k + p X i =1 k g ( i ) k  ǫ norm t | {z } = R h,err . This results in the zonotop e of th e theor em. Since R h,er r is a m ultidimension al interval, an d thu s also a zono tope, the addition of the two zo notopes R h,er r and ( ˆ c, ˆ g (1) , . . . , ˆ g ( p ) ) results in the zono tope R h ( t ) . W e p ropose two different representation s of R h,er r . Th e ﬁrst one uses the equ iv alen c e [ − 1 , 1 ] n  k c k + p X i =1 k g ( i ) k  ǫ norm t | {z } =:ˆ ǫ norm t = n M i =1 [ − 1 , 1] e i ˆ ǫ norm t, (12) where e i is a unit vector with the i th element being 1 and all others 0 . Let u s also introd uce 0 n and 1 n as an n-d imensional vector of zeros and one s, respectively . Using (12), we can write R h,er r = ( 0 n , ˆ h (1) , . . . , ˆ h ( n ) ) with ˆ h ( i ) j = ( ˆ ǫ norm t, f or i = j 0 , for i 6 = j , which makes it possible to obtain (see Th e o rem III .1) R h ( t ) ⊆ (ˆ c, ˆ g (1) , . . . , ˆ g ( p ) ) ⊕ R h,er r =(ˆ c, ˆ g (1) , . . . , ˆ g ( p ) , ˆ h (1) , . . . , ˆ h ( n ) ) . T o a void add ing new generato rs ˆ h (1) , . . . , ˆ h ( n ) , one can also bound R h,er r by an in te r val vector as presented in the next lemma. Lemma III.2 (Interval vector multiplicat ion) Using [ b , b ] = b c ⊕ [ − b ∆ , b ∆ ] ( b ∆ ≥ 0 ) we can pr ovide the bo und e A t [ b , b ] ⊆k b c k V ( b c ) e H ( b c ) t e 1 ⊕ [ − µ , µ ] , wher e µ = k b ∆ k ¯ V ( b ∆ ) e ¯ H ( b ∆ ) t e 1 + 1 n ( k b ∆ k ¯ ǫ norm t + k b c k ǫ norm t ) and ¯ V , ¯ H ar e obtained fr om Alg. 1 and ¯ ǫ norm t fr om App endix B using C = | A | . Please n ote that the ab solute values are computed elementwise, i.e., | A | i,j = | A i,j | . Pr oof: W e ﬁrst sho w that e A t [ b , b ] can be o ver- approx imated by e A t [ b , b ] = e A t ( b c ⊕ [ − b ∆ , b ∆ ]) ⊆ e A t b c ⊕ e A t [ − b ∆ , b ∆ ] ⊆ e A t b c ⊕ [ − | e A t | b ∆ , | e A t | b ∆ ] ⊆ e A t b c ⊕ [ − e | A | t b ∆ , e | A | t b ∆ ] , where the over -approx imation ach iev ed in the last line directly follows from the T aylor series of e A t and A i,j ≤ C i,j . Using the above result and Lem ma II I.1, we obtain e A t b c ⊕ [ − e | A | t b ∆ , e | A | t b ∆ ] = k b c k V ( b c ) e H ( b c ) t e 1 ⊕ [ − 1 , 1 ] n k b c k ǫ norm t ⊕ [ − γ , γ ] , (13) where γ = k b ∆ k ¯ V ( b ∆ ) e ¯ H ( b ∆ ) t e 1 + 1 n k b ∆ k ¯ ǫ norm t. The result in (13) can be simpliﬁed to k b c k V ( b c ) e H ( b c ) t e 1 ⊕ [ − µ , µ ] , with µ as in the lemm a . B. Solution for a T ime I n terval In the previous subsection , we ov er-approximated the ho- mogene o us solution fo r points in time. This subsection over- approx imates x h ( t ) for a time interval [0 , δ ] . Since x h ( t ) ∈ ˆ x h ( t ) ⊕ E red ( t, x 0 ) as shown in (10), where E red ( t, x 0 ) is monoto nically in creasing according to (11), we have that ∀ t ∈ [0 , δ ] : E red ( t, x 0 ) ⊆ E red ( δ, x 0 ) . It rem a ins to app roximate ˆ x h ( t ) within time inter vals: ∀ t ∈ [0 , δ ] : ˆ x h ( t ) ≈ x 0 + t δ ( ˆ x h ( δ ) − x 0 ) . T o con sider the error of th is appro x imation, we use a ﬁnite T aylo r expan sion of the exponen tial matrix of η th order w ith error matrix E (see [ 63, Prop. 2]): e H ( x 0 ) t = η X i =0 1 i ! ( H ( x 0 ) t ) i + E ( t, x 0 ) , (14) where E is enclo sed by an interval matrix : E ( t, x 0 ) ∈ E ( t, x 0 ) = [ − W ( t, x 0 ) , W ( t, x 0 )] , (15) W ( t, x 0 ) = ∞ X i = η +1 t i i ! | H ( x 0 ) | i = e | H ( x 0 ) | t − η X i =0 t i i ! | H ( x 0 ) | i . 6 W e a r e prop o sing an e nclosure of the error x h ( t ) − ˆ x h ( t ) , which is m ultiplicative with the initial state, since large initial states result in larger errors for a given time horizon . Lemma III.3 (Co rrection interval matrix F ) The in terval matrix F ( x 0 ) = η X i =2 [( i − i i − 1 − i − 1 i − 1 ) δ i , 0] H i ( x 0 ) i ! ⊕ E ( δ, x 0 ) , with E ( δ, x 0 ) an d η accor d ing to ( 1 5) en sur es the enclosure of the exact solution: ∀ t ∈ [0 , δ ] : ˆ x h ( t ) ∈ x 0 + t δ ( ˆ x h ( δ ) − x 0 ) ⊕ k x 0 k V ( x 0 ) F ( x 0 ) e 1 . (16) Pr oof: Let u s start by rearran ging (1 6): ˆ x h ( t ) − x 0 − t δ ( ˆ x h ( δ ) − x 0 ) ∈ k x 0 k V ( x 0 ) F ( x 0 ) e 1 . After replacing ˆ x h ( t ) by k x 0 k V ( x 0 ) e H ( x 0 ) t e 1 from (11) and using k x 0 k V ( x 0 ) e 1 = x 0 (this fo llows fro m the fact tha t the ﬁrst column of V ( x 0 ) is x 0 / k x 0 k according to Alg. 1 ), the above inclusion becom e s k x 0 k V ( x 0 ) e H ( x 0 ) t e 1 − k x 0 k V ( x 0 ) e 1 − t δ  k x 0 k V ( x 0 ) e H ( x 0 ) δ e 1 − k x 0 k V ( x 0 ) e 1  ∈k x 0 k V ( x 0 ) F ( x 0 ) e 1 . Next, w e use the result A ( B + C ) D = AB D + AC D from linear algebra to simplify the above eq uation to k x 0 k V ( x 0 )  e H ( x 0 ) t − I − t δ  e H ( x 0 ) δ − I   e 1 ∈k x 0 k V ( x 0 ) F ( x 0 ) e 1 . By comp aring the inner p art  o f k x 0 k V ( x 0 )  e 1 , one o btains ∀ t ∈ [0 , δ ] : e H ( x 0 ) t − I − t δ ( e H ( x 0 ) δ − I ) ∈ F ( x 0 ) . Substituting e H ( x 0 ) t by (1 4) and canc e ling linear terms yields η X i =2 ( t i − t δ i − 1 ) 1 i ! H i ( x 0 ) + E ( t, x 0 ) − t δ E ( δ , x 0 ) ∈ F ( x 0 ) . The interval of t i − t δ i − 1 for t ∈ [0 , δ ] is o btained exactly b y computin g th e minimum a nd maxim um o f t i − t δ i − 1 for which only one extreme value exists since th e second d e riv ative is nonnegative: d dt ( t i − t δ i − 1 ) = 0 ⇒ t min = i − 1 i − 1 δ . This means that the maximum values are to be fo u nd a t the bo rders of t ∈ [0 , δ ] , which are bo th 0 for t = 0 an d t = δ . Thu s, { t i − t δ i − 1 | t ∈ [0 , δ ] } = [ t i min − t min δ i − 1 , 0] = [( i − i i − 1 − i − 1 i − 1 ) δ i , 0] . It remains to b ound E ( t, x 0 ) − t δ E ( δ , x 0 ) f or t ∈ [0 , δ ] . Using (15) we have that E ( t, x 0 ) − t δ E ( δ , x 0 ) ≤     E ( t, x 0 ) − t δ E ( δ , x 0 )     (15) =     ∞ X i = η +1 t i i ! | H ( x 0 ) | i − t δ ∞ X i = η +1 δ i i ! | H ( x 0 ) | i     =     ∞ X i = η +1 ( t i − tδ i − 1 ) | H ( x 0 ) | i i !     = ∞ X i = η +1    t i − tδ i − 1    | H ( x 0 ) | i i ! ≤ ∞ X i = η +1 δ i | H ( x 0 ) | i i ! = E ( δ , x 0 ) . Thus, ∀ t ∈ [0 , δ ] : E ( t, x 0 ) − t δ E ( δ , x 0 ) ∈ E ( δ, x 0 ) , which completes the proo f . When the initial state x 0 is substituted by the set of initial states X 0 , Lemma III.3 can be gener alized as describ e d in the subsequen t theorem . Theorem III.2 (Solution for time intervals) The r eachable set for t ∈ [0 , δ ] is over -appr oximated by R h ([0 , δ ]) := [ t ∈ [0 ,δ ] R h ( t ) = conv  X 0 , R h ( δ )  ⊕ N , wher e N = k c k V ( c ) F ( c ) e 1 ⊕ p M i =1 k g ( i ) k V ( g ( i ) ) F ( g ( i ) ) e 1 . Pr oof: Let us start b y for mulating the reachable set for t ∈ [0 , δ ] usin g x h ( t ) ∈ ˆ x h ( t ) ⊕ E red ( t, x 0 ) from (1 0) and ˆ x h ( t ) from (16); we also use ∀ t ∈ [0 , δ ] : E red ( t, x 0 ) ⊆ E red ( δ, x 0 ) : R h ( t ) ⊆ n x 0 + t δ ( ˆ x h ( δ ) − x 0 )    x 0 ∈ X 0 o ⊕ [ x 0 ∈X 0 k x 0 k V ( x 0 ) F ( x 0 ) e 1 ⊕ t δ [ x 0 ∈X 0 E red ( δ, x 0 ) . (17) Since ∀ t ∈ [0 , δ ] : n x 0 + t δ ( ˆ x h ( δ ) − x 0 )    x 0 ∈ X 0 o ⊕ t δ [ x 0 ∈X 0 E red ( δ, x 0 ) ⊆ co nv ( X 0 , R h ( δ )) we can simplify (17) to R h ([0 , δ ]) ⊆ conv ( X 0 , R h ( δ )) ⊕ [ x 0 ∈X 0 k x 0 k V ( x 0 ) F ( x 0 ) e 1 ⊆ conv ( X 0 , R h ( δ )) ⊕ ˆ N , 7 where ˆ N := S x 0 ∈X 0 k x 0 k V ( x 0 ) F ( x 0 ) e 1 . It remain s to over - approx imate ˆ N wh en the initial set is a zon otope: [ x 0 ∈X 0 k x 0 k V ( x 0 ) F ( x 0 ) e 1 ⊆ k c k V ( c ) F ( c ) e 1 | {z } =: N (0) ⊕ p M i =1 [ − 1 , 1] k g ( i ) k V ( g ( i ) ) F ( g ( i ) ) e 1 | {z } = k g ( i ) k V ( g ( i ) ) F ( g ( i ) ) e 1 =: N ( i ) =: N , and th e in terval vectors N ( i ) are added using standar d interval arithmetic [64]. Thu s, N is an interval vector , which is conv erted to a zono tope and add ed to the c o n vex hull. The scalar interval [ − 1 , 1] can be moved rig ht in f ront of F ( g ( i ) ) and since all entries of F ( g ( i ) ) ha ve symmetric bounds, we have that [ − 1 , 1] F ( g ( i ) ) = F ( g ( i ) ) so tha t [ − 1 , 1 ] can be removed. Next, we derive the set of solu tions due to uncer ta in inpu ts. I V . I N P U T S O L U T I O N I N T H E K RY L OV S U B S PAC E In this section, we ob ta in for the ﬁrst time the set of input solutio ns for uncertain , time-varying inp uts in the Kr ylov subspace. W e ﬁrst co nsider over -a p prox im ations fo r solutio ns of a single con stant in put. Next, we g eneralize this re sult to uncertain b u t constant inputs. Finally , we deri ve an over- approx imation for arbitra r ily-varying an d uncertain inp uts. Th e ﬁrst lemma pr ovides the over-approxim ation fo r a con stant, known inpu t. Lemma IV .1 (Krylov er ror o f constant input so lution) After obtain ing ˜ V ( u ) , ˜ H ( u ) fr om Alg. 1 a nd ˜ ǫ norm δ fr o m Append ix B with in p uts C =  A u 0 1 × n 0  , v =  0 n 1  , we can b ound the p articular so lu tion ( aka in put solution) fo r constant inputs x p, const ( δ ) = Z δ 0 e A ( δ − t ) d t u = Z δ 0 e At d t u (18) by x p, const ( δ ) ∈ P ˜ V ( u ) e ˜ H ( u ) δ e 1 ⊕ [ − 1 , 1 ] n ˜ ǫ norm δ = ˜ x p ( δ ) ⊕ ˜ E red ( δ ) , wher e P = [ I , 0 n ] , ( I is the n × n identity matrix) ˜ x p ( δ ) = P ˜ V ( u ) e ˜ H ( u ) δ e 1 ˜ E red ( δ ) = [ − 1 , 1 ] n ˜ ǫ norm δ. (19) Pr oof: W e ﬁrst rewrite the solu tion of (18) to the differ- ential equatio n be low as also performe d in [57, Sec. 5]:  ˙ x p, const ˙ ˇ x  | {z } ˙ ˜ x =  A u 0 1 × n 0  | {z } =: ˜ A ( u )  x p, const ˇ x  | {z } ˜ x , ˜ x (0) = ˜ x 0 =  0 n 1  . After inserting the projection matrix P = [ I , 0 n ] (see (19)), we can write x p, const ( t ) = P e ˜ A ( u ) t ˜ x 0 . This makes it p ossible to reform ulate the Krylov error of the input solution: k P e ˜ A ( u ) t ˜ x 0 | {z } x p, const ( t ) − P k ˜ x 0 k ˜ V ( u ) e ˜ H ( u ) t e 1 | {z } Krylov approximation of x p, const ( t ) k ≤k P k k e ˜ A ( u ) t ˜ x 0 − k ˜ x 0 k ˜ V ( u ) e ˜ H ( u ) t e 1 k , ( k P k = k ˜ x 0 k = 1) = k e ˜ A ( u ) t ˜ x 0 − ˜ V ( u ) e ˜ H ( u ) t e 1 k ≤ ˜ ǫ norm t and ˜ ǫ norm t is obta in ed as in App endix B. In orde r to gen eralize the pr evious results to arbitrar ily- varying inp uts, we requ ir e the following co rollary o n the input solution for general time boun ds: Corollary IV .1 (Input solution for genera l time bounds) The partial input solution for x p, const ([ t 0 , t e ]) := Z t e t 0 e A ( δ − t ) d t u , wher e 0 ≤ t 0 ≤ t e ≤ δ , can be over-appr oxima ted as x p, const ([ t 0 , t e ]) ∈ P ˜ V ( u )( e ˜ H ( u ) ( δ − t 0 ) − e ˜ H ( u ) ( δ − t e ) ) e 1 ⊕ [ − 1 , 1 ] n ˜ ǫ norm ( t e − t 0 ) . Pr oof: Let u s rewrite x p, const ([ t 0 , t e ]) as x p, const ([ t 0 , t e ]) := Z t e t 0 e A ( δ − t ) d t u = Z δ − t 0 δ − t e e At d t u =  Z δ − t 0 0 e At d t − Z δ − t e 0 e At d t  u. After applying Lemma IV .1, we obtain x p, const ([ t 0 , t e ]) = P ˜ V ( u ) e ˜ H ( u ) ( δ − t 0 ) e 1 − P ˜ V ( u ) e ˜ H ( u ) ( δ − t e ) e 1 ⊕ [ − 1 , 1 ] n ˜ ǫ norm ( δ − t 0 − ( δ − t e )) = P ˜ V ( u )( e ˜ H ( u ) ( δ − t 0 ) − e ˜ H ( u ) ( δ − t e ) ) e 1 ⊕ [ − 1 , 1 ] n ˜ ǫ norm ( t e − t 0 ) . Next, we over-approximate the solution for con stant inpu ts which are u ncertain with in the set U . In many cases, it is desired to have the solution f o r con stant inputs, e.g., if a control system with zero-o r der hold is co n sidered. Theorem IV .1 (Rea chable set for co nstant inputs) The r ea chable set of the in put solution ˆ R p ( δ ) := Z δ 0 e A ( δ − t ) d t U = Z δ 0 e At d t U , wher e U = ( c u , g (1) u , . . . , g ( q ) u ) is a zonotope, c an be over- appr o ximated by the zonotope ˆ R p ( δ ) ⊆ ( ˜ c u , ˜ g (1) u , . . . , ˜ g ( q ) u ) ⊕ ˆ R p,err , ˜ c u = P ˜ V ( c u ) e ˜ H ( c u ) δ e 1 , ˜ g ( i ) u = P ˜ V ( g ( i ) u ) e ˜ H ( g ( i ) u ) δ e 1 , ˆ R p,err = [ − 1 , 1 ] n  ˜ ǫ norm ( c u ) + q X i =1 ˜ ǫ norm ( g ( i ) u )  δ. 8 Pr oof: After inserting th e d e ﬁn ition of a zon otope (Def. II.1) into ˆ R p ( δ ) := R δ 0 e A ( δ − t ) d t U an d using (4), we obtain ˆ R p ( δ ) = Z δ 0 e A ( δ − t ) d t c u ⊕ q M i =1 [ − 1 , 1] Z δ 0 e A ( δ − t ) d t g ( i ) u . Using Lemma IV .1, one receives ˆ R p ( δ ) ⊆ P ˜ V ( c u ) e ˜ H ( c u ) δ e 1 ⊕ [ − 1 , 1 ] n ˜ ǫ norm ( c u ) δ ⊕ q M i =1  [ − 1 , 1] P ˜ V ( g ( i ) u ) e ˜ H ( g ( i ) u ) δ e 1 ⊕ [ − 1 , 1 ] n ˜ ǫ norm ( g ( i ) u ) δ  = P ˜ V ( c u ) e ˜ H ( c u ) δ e 1 ⊕ q M i =1 [ − 1 , 1] P ˜ V ( g ( i ) u ) e ˜ H ( g ( i ) u ) δ e 1 ⊕ [ − 1 , 1 ] n  ˜ ǫ norm ( c u ) + q X i =1 ˜ ǫ norm ( g ( i ) u )  δ | {z } = ˆ R p,err . This results in the zon otope of th e theo r em. Since ˆ R p,err is a m ultidimension al interval, and thus also a zo notope, the addition o f th e two zonoto pes ˆ R p,err and (˜ c u , ˜ g (1) u , . . . , ˜ g ( q ) u ) results in the zono tope ˆ R p ( δ ) . The above deri vations on ly hold for constant inputs. W e generalize the previous r esults to arbitrar y inp ut trajectories in the following theor em. Theorem IV .2 (Rea chable set for vary ing inputs) The r ea chable set du e to uncertain inputs R p ( δ ) = n x p ( δ )    x p ( δ ) = Z δ 0 e A ( δ − t ) u ( t ) dt, ∀ t ∈ [0 , δ ] : u ( t ) ∈ U o can be over -appr oximated by R p ( δ ) ⊆ P  η M j =1 n ˜ V ( u ) ˜ H j ( u ) j !    u ∈ U o δ j ⊕ n ˜ E ( δ, u )    u ∈ U o e 1 ⊕ [ − 1 , 1 ] n n ˜ ǫ norm ( u )    u ∈ U o δ, wher e ˜ E ( δ, u ) is obtained fr om (15) by r ep lacing H ( x 0 ) with ˜ H ( u ) . The proof can be fo u nd in Append ix A. Obviously , ˆ R p ( δ ) ⊆ R p ( δ ) , since ˆ R p ( δ ) is exact except for adding ˆ R p,err , which is id e ntically added to R p ( δ ) , an d R p ( δ ) also co ntains all solutions for non -constant inp uts. Next, we con sider a particu lar solution of Theo rem IV .2 when U is a zonotope . Corollary IV .2 (V a rying inputs with zonot opic bounds) The r eachab le set du e to uncertain inputs R p ( δ ) = n x p ( δ )    x p ( δ ) = Z δ 0 e A ( δ − t ) u ( t ) dt, ∀ t ∈ [0 , δ ] : u ( t ) ∈ U o when U = ( c u , g (1) u , . . . , g ( q ) u ) is a zon otope, can b e over- appr o ximated by R p ( δ ) =( ˜ c u , ˜ g (1 , 1) u , . . . , ˜ g ( q, 1) u , ˜ g (1 , 2) u , . . . , ˜ g ( q, 2) u , . . . , ˜ g ( q,η ) u ) ⊕ ˆ R p,err , ˜ c u = P  η X j =1 ˜ V ( c u ) ˜ H j ( c u ) j ! δ j  e 1 , ˜ g ( i,j ) u = P ˜ V ( g ( i ) u ) ˜ H j ( g ( i ) u ) j ! δ j e 1 , R p,err = P  ˜ E ( δ, c u ) ⊕ q M i =1 ˜ E ( δ, g ( i ) u )  e 1 ⊕ [ − 1 , 1 ] n  ˜ ǫ norm ( c u ) + q X i =1 ˜ ǫ norm ( g ( i ) u )  δ. Pr oof: After inserting the deﬁnition of a zo notope (Def. II.1) into the result of Theo rem I V .2, we obtain R p ( δ ) ⊆ P  η M j =1  ˜ V ( c u ) ˜ H j ( c u ) j ! ⊕ q M i =1 [ − 1 , 1] ˜ V ( g ( i ) u ) ˜ H j ( g ( i ) u ) j !  δ j ⊕ ˜ E ( δ, c u ) ⊕ q M i =1 [ − 1 , 1] ˜ E ( δ, g ( i ) u ) | {z } = ˜ E ( δ,c u ) ⊕ L q i =1 ˜ E ( δ,g ( i ) u )  e 1 ⊕ [ − 1 , 1 ] n  ˜ ǫ norm ( c u ) ⊕ q M i =1 [ − 1 , 1]˜ ǫ norm ( g ( i ) u )  δ | {z } =[ − 1 , 1 ] n  ˜ ǫ norm ( c u )+ P q i =1 ˜ ǫ norm ( g ( i ) u )  δ , which can be simpliﬁed to R p ( δ ) ⊆ P  η X j =1 ˜ V ( c u ) ˜ H j ( c u ) j ! δ j  e 1 ⊕ η M j =1 q M i =1 [ − 1 , 1]  P ˜ V ( g ( i ) u ) ˜ H j ( g ( i ) u ) j ! δ j  e 1 ⊕ R p,err . This results in the zon otope of th e theo r em. Since R p,err is a multidimensional interval, and thus also a zono - tope, the addition of the two zo notop e s R p,err and (˜ c u , ˜ g (1 , 1) u , . . . , ˜ g ( q, 1) u , ˜ g (1 , 2) u , . . . , ˜ g ( q, 2) u , . . . , ˜ g ( q,η ) u ) results in the zonotop e R p ( δ ) . V . P RO P A G AT I O N So far , we h av e only described the comp utation of a reachable set fo r the ﬁrst time inter val. This section presents how the initial re sults are prop agated fo r further consecu tiv e time intervals τ i . 9 A. Overall Alg orithm T o g rasp th e implemen ted prop agation scheme more easily , we start with one of the simplest p r opagatio n procedu res (see [22, Alg. 1]): R h ( τ k +1 ) = e Aδ R h ( τ k ) , (20) R p ( τ k +1 ) = e Aδ R p ( τ k ) ⊕ R p ( τ 0 ) . (21) In or der to keep the number of g enerator s of zonotopes that are mu ltiplied by th e matrix exponential small, the wrap ping- free appro ach from [65] is used in th is work as a basis; we later in troduce mod iﬁcatio ns to make the best u se of th e propo sed computation in the Kry lov subspace. The wrappin g- free app roach modiﬁes the above procedure by introdu cing the auxiliary reachable set R b , which is enclosed by an axis- aligned b ox: R p ( τ 0 ) = bo x ( R b ( τ 0 )) ( see line 3 of Alg. 2 and (6)). This has the effect that the re p resentation of R p ( τ k +1 ) in (21) d oes not gr ow i n comp lexity due to the Minkowski addition, but stays a simple ax is-aligned box . The mo d iﬁed computatio n of R p ( τ k +1 ) accordin g to [65] is: R b ( τ k +1 ) = e Aδ R b ( τ k ) , (22) R p ( τ k +1 ) = R p ( τ k ) ⊕ bo x ( R b ( τ k +1 )) . (23) In o r der to take advantage o f the fact that the Arnoldi iteration does not hav e to b e recomp uted for different times as shown in (2), we furth er mo dify the c omputatio n of R b ( τ k +1 ) to R b ( τ k +1 ) = e At k R b ( τ 0 ) . (24) Another m odiﬁcation, which is no t pro posed in [65], is that we change (20): instead o f comp uting the h omoge n eous solution for the ﬁrst time interval o nly and then pr o pagating it, we apply Theo r em II I.2 in each time step. Th e reaso n is that for large systems, one saves much comp u tation time when the matrix expo nential multiplication is p erforme d with R h ( t k ) instead of R h ( τ k ) due to a fewer n umber of generators, while the com putation in Theorem III.2 is negligib le whe n using z onotop es. Since the sets in Alg. 2 are in dexed withou t explicitly indicating tim e , we d istinguish sets representin g points in time by an asterisk ( see e.g., R ∗ h, 1 ) from th e ones representin g time intervals. T o sum up, in the Krylov subspace, the computatio n of (20) and (2 4) is c a r ried ou t by Theorem II I.1 (see line 7 and line 8 in Alg . 2 ). T he propagation of the pa r ticulate solutio n (23) is r ealized in line 9 o f Alg. 2 and the ho mogen eous solu tion for a time interval is o b tained in line 10. The aggr egatio n of the homogene o us and the particulate solu tion are p erforme d in line 5 and line 11 of Alg. 2. When only o utput values y ( t ) = C x ( t ) + D u are o f interest, one can save sub stantial computation time by ﬁrst multiplyin g the ortho gonal b ases V ( c ) , V ( g ( i ) ) , ˜ V ( c u ) , and ˜ V ( g ( i ) u ) with C to obtain V ( c ) := C V ( c ) , V ( g ( i ) ) := C V ( g ( i ) ) , ˜ V ( c u ) := C ˜ V ( c u ) , an d ˜ V ( g ( i ) u ) := C ˜ V ( g ( i ) u ) as also presented in e.g., [29]. Otherwise, on e would ﬁr st co mpute the full-dim ensional reachable set an d afterwards compute the projec tio n instead of directly comp u ting the r eachable set of the outpu t. When only constant inputs are considered, coroll ary _IV.2 ( A, B , U , δ ) in line 2 is replaced b y theore m _IV.1 ( A, B , U , δ ) . Also, when the center u c of the Algorithm 2 Compute R ([0 , t f ]) Require: State m a trix A , input matrix B , initial set X 0 , input set U , time step δ , time horizon t f Ensure: R ([0 , t f ]) 1: R ∗ h, 1 = theore m _III.1 ( X 0 , A, δ ) 2: R ∗ b, 1 = coroll ary _ IV.2 ( U , A, B , δ ) 3: R p, 1 = bo x  conv ( 0 n , R ∗ b, 1 )  4: R h, 1 = theore m _III.2 ( X 0 , R ∗ h, 1 , A, δ ) 5: R 1 = R h, 1 ⊕ R p, 1 6: for k = 1 . . . ⌈ t f /δ − 1 ⌉ do 7: R ∗ h,k +1 = theore m _III.1 ( X 0 , A, ( k + 1 ) δ ) 8: R ∗ b,k +1 = th eorem _III.1 ( R ∗ b, 1 , A, ( k + 1) δ ) 9: R p,k +1 = R p,k ⊕ bo x ( R ∗ b,k +1 ) 10: R h,k +1 = theore m _III.2 ( R ∗ h,k , R ∗ h,k +1 , A, δ ) 11: R k +1 = R h,k +1 ⊕ R p,k +1 12: end for 13: R ([0 , t f ]) = S t f /δ k =1 R k set of uncer tain inpu ts U is lar ge com p ared to the d eviation from the cen ter U ∆ := U ⊕ ( − u c ) , one should move the solution of u c inside th e con vex hu ll co mputation per formed in the orem _III.2 ( A, X 0 , δ ) in line 1 0 (see e .g., [66, Alg. 1]). This, howe ver, is independ ent from the extension to compute in the Krylov subspace and thus not d iscussed in this work. B. Computatio nal Complexity Let us ﬁrst con sid er the computatio nal complexity when applying Alg. 2 without utilizing the Krylov subspace: For each time step, we have to map zonotopes using the matrix exponential, com pute the box enclosu re as well as the conv ex hull of zonoto p es, an d add two zonotop e s. W e intro duce the number of gen erators o f R ∗ h, 1 in Alg. 2 as p h and the number of gen erators of R ∗ b, 1 as p b . Thus, the o rders of those zonotop es are o h = p h n and o b = p b n . The number o f required binary ope rations for each of the previously m entioned high - lev el operands are listed in T ab. I. W e are interested in th e complexity with respect to the dimension n , since the n umber of time steps is typ ically ﬁxed or g i ven by r eaching a ﬁxed- point. As can be seen from T ab . I , the c o mplexity with resp ect to n is cubic for linear maps ( when M is quadratic as for e Aδ ), lin ear for addition, qu adratic fo r over-approximating the conv ex hull, and quadratic fo r the box enclosure. Thus, the overall complexity is dom inated by the linear map, which has complexity O (( o h + o b ) n 3 ) . Note that according to Alg. 2, the ord e r o f the inv olved zon otopes does n o t grow compare d to other propa gation algorithm s ( see e. g., [22], [61]). The only difference wh en comp uting in the Krylov subspace is that the complexity of com puting e Aδ Z chang es. As dis- cussed above, this is also the operation wh ich determines the overall compu tational complexity . W e are d iscussing the case when the Krylov order is chosen such that the expon ential matrix is compu ted up to machine pre c ision so that the results have the same accuracy and the co mparison is fair . The Arnoldi iteration req uires ξ matrix- vector mu ltiplications ( ξ is the dimension of the Krylov sub space; see Alg. 1). Howe ver, 10 since the in volved m atrices are spar se, we can add a spar sity constant s << 1 deﬁned as the fraction o f non-zero entries so that we obtain O ( sξ n 2 ) opera tions. Further, we r equire ξ 2 / 2 inner prod ucts ( O ( ξ 2 n ) ) an d ξ 2 / 2 other oper a tio ns with O ( n ) resulting in O ( ξ 2 n ) , so that the overall co mplexity is O ( ξ 2 n ) + O ( sξ n 2 ) . The Arnold i iteration has to b e computed for ( o h + o b ) n generato r s, so tha t the comp u tational complexity amou nts to O ( ξ 2 ( o h + o b ) n 2 + sξ ( o h + o b ) n 3 ) concerning the Arnoldi iteration. The required compu tation of each expon ential m atrix e H δ in the Krylov su b space o nly depe n ds on ξ and n ot on n . Finally , we r equire matrix/m atrix co mputation s between V and the exponen tial matrix e H δ with comp lexity O ( nξ 2 ) . Thus, the overall com plexity using the Krylov techn ique is O ( ξ 2 ( o h + o b ) n 2 + sξ ( o h + o b ) n 3 ) + O ( nξ 2 ) = O ( ξ 2 ( o h + o b ) n 2 + sξ ( o h + o b ) n 3 ) . A s a result, b oth the classical appr oach and the Kr y lov app roach ar e cubic in th e number of continuo us state variables. Howe ver , the Krylov method can be mu ch faster (as d e m onstrated in the next section ) d ependin g on the sp arsity of the system matr ix A , since ξ << n a n d typically sξ << 1 (compar e O (( o h + o b ) n 3 ) for th e classical approa c h with O ( ξ 2 ( o h + o b ) n 2 + sξ ( o h + o b ) n 3 ) for th e Krylov approa c h ). Please no te that the spa rsity typ ically increa ses the larger a system is ( see e.g ., [6 7, Fig. 3]). Also, it sho uld b e noted th a t although A is o ften very sparse, all entries of e Aδ are ty pically non-zero . Fur th er , observations show that ξ and n are n ot linearly re la te d for a co nstant error ǫ , b ut that n ξ grows con siderably with increasing n [59]. V I . N U M E R I C A L E X P E R I M E N T S T o demonstrate the u sefulness of the presented appro ach, we compu te reachable sets of ben chmark models in [6 8] and perfor m a for m al analysis of a bridge as a represen ta tive of a safety-critical stru cture. T o the b est knowledge of th e auth or , no work on form ally b o unding values of a safety-c r itical struc- ture exists; the most prominent example of a failing safety - critical structure is the collapsed T acoma Narrows Bridge [69]. A. Models a) Benchmark Models: Ben c hmarks for large linear sys- tems hav e b een p roposed in [6 8]. W e have selected the FOM model and the MNA mode ls; the MN A models are inspired T ABLE I R E QU I R E D O P E R A T I O N S . operands Z h = ( c h , g (1) h , . . . , g ( p h ) h ) , ˆ Z h = (ˆ c h , ˆ g (1) h , . . . , ˆ g ( p h ) h ) , Z p = ( c p , g (1) p , . . . , g ( p p ) p ) ⊂ R n , M ∈ R m × n operatio n nr . of binary operations M Z h , see (4) 2 mn ( o h n + 1) Z h ⊕ Z p , see (4) n conv ( Z h , ˆ Z h ) , see (5) 2( o h n + 1) n box ( Z h ) , see (6) o h n 2 by linear differential alge b raic systems E ˙ x = Ax ( t ) + B u ( t ) in [ 70]. In [68] o nly MNA-1 and MN A-5 are presented and the E matrix is ch o sen as the id entity matrix to obtain a ordinar y differential equ ations. W e h av e additionally included the mod e ls MN A- 2 up to MNA-4 b y also cho osing E as the identity matrix. b) B ridge Mo del: The model of the Roosevelt Lake Bridge (Ar izona) is taken from a student work [71], wh ich in vestigated its structu ral dynam ic s. A picture 1 of the b r idge and the correspond ing ﬁnite element mode l are shown in Fig. 3. The bridge is mostly made ou t of steel, except f or the roadway deck. The dam struc tu re in fro nt of th e bridg e (see Fig. 3 (a)) creates a special aerody namic situation, which can excite certain natur al fre quencies of the b ridge [7 2]. Th e bridge model co nsists of 445 nod e s connecting bars and beam s. The exact m odel is provided as an example in CORA [4 4]. Using the pr o posed reachability an alysis, we can in vestigate the e n tire r ange o f possible time responses for all k in ds o f frequen cies o f the exciting forces at once . Thus, we do not miss a single possible freq uency that could cause a pr oblem. P S f r a g r e p l a c e m e n t s r e f e r e n c e t r a j e c t o r y v e h i c l e I (a) Picture of the bridge with the dam accele rating winds in the background. 150 100 50 0 -50 0 20 -100 40 60 0 -150 P S f r a g r e p l a c e m e n t s r e f e r e n c e t r a j e c t o r y v e h i c l e I e g o v e h i c l e p r e v i o u s l y v e r i ﬁ e d p a r t e g o v e h i c l e n e w i n t e n d e d p a r t (b) Finite-eleme nt model of the bridge (distances are in meters). Fig. 3. Picture an d ﬁnite -element model of the R oose velt Lake Bri dge (Arizona) . W e consider two excitations of the bridge: (1) The b ridge is excited by lateral forces simulating the win ds acting on the b ridge and the street deck and (2) vertical excitation is generated b y the street deck c aused by movin g trafﬁc o n the bridge. It should be noted that we consider all k inds of frequen cies of these excitation s so that n o p o ssible solutio n is missed. Due to the large r epresentation size o f the inpu t set X 0 and th e input vector u ( 50 4 0 dimension s), we provide those values within the u p loaded examp le in CORA. B. Results T o b e as pr ecise as po ssible, we have ch osen th e ord er reduction so th a t k v k ǫ norm is belo w the precision o f ﬂo ating 1 The picture is taken from commons.wikimedia .org 11 point numbe r s in MA TLAB, which is 2 . 220 4 · 10 − 16 . The results below are obtained using a stand ard laptop with an Intel i7-352 0M CPU ru nning at 2 .90GHz and 8GB of RAM. c) Reachability Analysis: A summ ary of the results for reachability analysis is presented in T ab . I I. The r e, the com - putation times for r eachability analy sis in the un tr ansforme d space and in the Krylov subspa c e is sh own. The c omputatio n times inclu de all processes, in c luding determ ining the pr o per order of th e K r ylov subsp a ce. W e also show th e different computatio n times depending on whether the f u ll state is required o r only the ou tputs of the system. T o better judge the results, we hav e also ad d ed the state d imension, the input dimen sion, the output dim ension, the time step size, and the maximum Krylov or der . Since the cen ter and each generato r might h av e d ifferent Krylov o r ders, we just present the maxim u m value. In order to determin e the pro per Kr ylov or d er fulﬁlling the error k v k ǫ norm in th is w ork, we iteratively in crement the Krylov order by 20 starting f rom 1 . Please note that the step sizes are rather small since we co nsider the full models a n d do not reduce them by removing higher frequencies. For all models we have ch osen th e initial values of the ﬁrst ten state variables to b e u ncertain within the in terval ζ [ − 1 , 1 ] an d all inputs to be uncer tain within the interval ϑ [ − 1 , 1] . Selected plots o f reachab le sets over time of the b enchmar k examples are p resented in Fig. 4. The results of the brid ge model are illustrated using selected two-dimensional pr ojec- tions of the reachab le set in Fig . 5. 0 0.05 0.1 t -500 0 500 1000 1500 y (a) FOM model. 0 5 10 t -10 -5 0 5 y 9 (b) MNA-5 model. Fig. 4. Reachab le set of outputs for selecte d benchmark models plotted over time. Black lines sho w random simulation s, the gray area shows the reachable set. d) V eriﬁcation: While th e focu s of this work is on reachability analysis, we also brieﬂy demo nstrate a possible veriﬁcation of the b ridge fo r which we ch eck wh ether all states stay within certain axis-aligned bo unds. Since th e bounds consist of 2 · 5040 = 100 80 values, we will only provid e them within CORA. I n th e examp le, the reachable set respects the given bounds; checking the boun ds for all time inte r val takes aroun d 2 0 s, which is small co mpared to 38 35 s fo r the reachability analysis. e) Discussion: Overall, the compu ta tio nal b e n eﬁts o f th e Krylov metho ds typica lly impr ove with the size of the model as shown in T ab. II. When only the outp ut is o f imp ortance, the Krylov method is particular ly efﬁcient since one ca n direc tly include the outp ut matr ix C in the K r ylov meth od by u sing the orthon ormal basis V ( v ) := C V ( v ) as d iscu ssed in Sec. V - A (see also [29]) . Th is av oids c omputin g th e reachable set of all states followed by projecting the result onto the outputs as requ ired by standard metho ds. For th e FOM mo d el and the bridg e mod el, the Kr ylov meth od is aroun d 100 times faster , wh ile the improvements for the MN A mo dels are also substantial—only the sma ll MNA-1 mod el shows similar computatio n times for bo th mo dels. In g eneral, the com putation times of th e K ylov method are less pred ic ta b le com pared to the standard appro ach g iv en the dimension o f the system. While the compu tation tim es o f the standard approa c h are ty p ically monoto nically increasin g with the system dimen sion, the Kry lov meth ods som etimes pr ovides faster resu lts for larger systems. V I I . C O N C L U S I O N S W e hav e pr esented th e ﬁrst w o rk fo r c omputin g t he o ver - approx imative reach able set of linear systems in the Krylov subspace fo r arbitrar ily-varying, boun ded inputs. Unlike most other work o n applying redu c tio n techniq ues fo r re achability analysis, which comp ute an output abstraction, we perfo rm a state ab straction and can f u lly r econstruct the c omplete reachable sets. When using outp ut abstractions, only the reachable outputs can be over-approxim ated—if those ou tput- abstraction techn iq ues would b e used to reconstru ct th e whole reachable set, no r eduction would b e achieved. Du e to the strict con sideration of error bou nds, our a p proach can be used f or formal veriﬁcation an d o th er f o rmal techn iques, such as comp utation o f in variance sets, com putation of the region of attraction, op timization of constrain ed system s with uncertainties, set-based observers, and conf ormance c h ecking. In our nu m erical example we have seen a speed up of two orders of magnitu de althou gh the appr oximation was accurate up to ﬂoating point precision . A C K N OW L E D G M E N T This work was partly supported by the German Research Foundation ( DFG) un der grant number AL 1185/5 -1. A P P E N D I X A P R O O F O F T H E O R E M I V . 2 Pr oof: W e ﬁrst compu te a n over-approx im ativ e r eachable set resulting fro m inputs when assuming constant i nputs for time intervals [ ˆ t k − 1 , ˆ t k [ , wh ere 0 = ˆ t 0 < ˆ t 1 < . . . < ˆ t l − 1 < ˆ t l = δ , of equ al d uration ρ = ˆ t k − ˆ t k − 1 , so that u ( t ) =            u (0) for t ∈ [0 , ˆ t 1 [ u ( ˆ t 1 ) for t ∈ [ ˆ t 1 , ˆ t 2 [ . . . u ( ˆ t l − 1 ) for t ∈ [ ˆ t l − 1 , δ ] . In a seco nd step, we let ρ → 0 to ob tain an over -approx imation for arbitrarily varying inputs. The so lution of p iecewise con- stant inputs is obtained from x p , pw ( δ ) = Z t 1 0 e A ( δ − t ) dt u (0 ) + Z t 2 t 1 e A ( δ − t ) dt u ( t 1 ) + . . . + Z δ t l − 1 e A ( δ − t ) dt u ( t l − 1 ) . 12 T ABLE II C O M P U TA T I O N T I M E S A N D P A R A M E T E R S ( O O M : O U T O F M E M O RY ) . Computation times [s] Model and computation parameters States Outputs Dimension Model Standard Krylov Standard Krylov State Input Output Krylov δ t f ζ ϑ FOM 3070 308.7 3306 21.19 1006 1 1 120 10 − 4 10 − 1 10 0.1 MNA-1 8 3.26 94.83 85.30 82.48 578 9 9 240 10 − 5 10 − 3 100 0.1 MNA-2 OoM OoM OoM 1234 9223 18 18 700 10 − 5 10 − 3 10 0.1 MNA-3 401 1.99 1601 3987.14 867.1 4863 22 22 780 10 − 5 10 − 3 10 0.1 MNA-4 3 44.1 44.11 333.7 21.36 980 4 4 260 10 − 5 10 − 3 10 0.1 MNA-5 OoM 2850 OoM 1 80.6 10913 9 9 440 10 − 1 10 +1 10 0.1 Bridge 490912 3835 490912 3835 5040 2 5040 380 10 − 7 10 − 4 0.001 0.1 Fig. 5. T wo-dimensional projections of reachable sets for select ed projections of the bridge model. Black lines sho w random simulations, the gray area sho ws the reach able set. When the set of inpu ts is u ncertain within U , we obtain ˆ R p,l ( δ ) = Z t 1 0 e A ( δ − t ) dt U ⊕ . . . ⊕ Z δ t l − 1 e A ( δ − t ) dt U . (25) In ord er to ob tain not o nly an appro ximation, b ut an over - approx imation, the solution for a time interval [ t k − 1 , t k [ is further abstracted. From Corrollar y IV .1 we have that ˆ R p ([ t 0 , t e ]) := Z t e t 0 e A ( δ − t ) d t U ⊆ P n ˜ V ( u )  e ˜ H ( u ) ( δ − t 0 ) − e ˜ H ( u ) ( δ − t e )     u ∈ U o e 1 ⊕ [ − 1 , 1 ] n n ˜ ǫ norm ( u )    u ∈ U o ( t e − t 0 ) . The a bove over-approximation is rewritten using a ﬁn ite T aylor 13 series (see (14)): ˆ R p ([ t 0 , t e ]) ⊆ P n ˜ V ( u )  I + ˜ H 1 ( u ) 1! ( δ − t 0 ) 1 + . . . + ˜ H η ( u ) η ! ( δ − t 0 ) η + ˜ E ( δ − t 0 , u ) −  I + ˜ H 1 ( u ) 1! ( δ − t e ) 1 + . . . + ˜ H η ( u ) η ! ( δ − t e ) η + ˜ E ( δ − t e , u )     u ∈ U , ˜ E ( t, u ) ∈ ˜ E ( t, u ) o e 1 ⊕ [ − 1 , 1 ] n n ˜ ǫ norm ( u )    u ∈ U o ( t e − t 0 ) = P n ˜ V ( u ) ˜ H 1 ( u ) 1! [( δ − t 0 ) 1 − ( δ − t e ) 1 ] + . . . + ˜ V ( u ) ˜ H η ( u ) η ! [( δ − t 0 ) η − ( δ − t e ) η ] + ˜ E ( δ − t 0 , u ) − ˜ E ( δ − t e , u )    u ∈ U , ˜ E ( t, u ) ∈ ˜ E ( t, u ) o e 1 ⊕ [ − 1 , 1 ] n n ˜ ǫ norm ( u )    u ∈ U o ( t e − t 0 ) . Finally , th e un certainty o f the inpu t is moved inwards, which results in a furthe r over-approximation : ˆ R p ([ t 0 , t e ]) ⊆ P  n ˜ V ( u ) ˜ H 1 ( u ) 1!    u ∈ U o [( δ − t 0 ) 1 − ( δ − t e ) 1 ] ⊕ . . . ⊕ n ˜ V ( u ) ˜ H η ( u ) η !    u ∈ U o [( δ − t 0 ) η − ( δ − t e ) η ] ⊕ n ˜ E ( δ − t 0 , u ) − ˜ E ( δ − t e , u )    u ∈ U , ˜ E ( t, u ) ∈ ˜ E ( t, u ) o  e 1 ⊕ [ − 1 , 1 ] n n ˜ ǫ norm ( u )    u ∈ U o ( t e − t 0 ) . (26) After introdu cing D ( j ) := n ˜ V ( u ) ˜ H j ( u ) j !    u ∈ U o , one can rewrite (26) a s ˆ R p ([ t 0 , t e ]) = Z t e t 0 e A ( δ − t ) d t U ⊆ P  D (1) [( δ − t 0 ) 1 − ( δ − t e ) 1 ] ⊕ . . . ⊕ D ( η ) [( δ − t 0 ) η − ( δ − t e ) η ] ⊕ n ˜ E ( δ − t 0 , u ) − ˜ E ( δ − t e , u )    u ∈ U , ˜ E ( t, u ) ∈ ˜ E ( t, u ) o  e 1 ⊕ [ − 1 , 1 ] n n ˜ ǫ norm ( u )    u ∈ U o ( t e − t 0 ) . (27) By repeatedly inserting (27) into ( 25), rearrang ing time in- tervals ( δ − t k ) = t l − k , an d using ˜ E ( t, u ) ∈ ˜ E ( t, u ) = [ − ˜ W ( t, u ) , ˜ W ( t, u )] , where ˜ W ( t, u ) = P ∞ i = η +1 t i i ! | ˜ H ( u ) | i from (15), we obtain ˆ R p,l ( δ ) = P  ( D (1) ( t l − t l − 1 ) ⊕ . . . ⊕D ( η ) ( t η l − t η l − 1 ) ⊕ ( D (1) ( t l − 1 − t l − 2 ) ⊕ . . . ⊕D ( η ) ( t η l − 1 − t η l − 2 ) ⊕ . . . ⊕ ( D (1) ( t 1 − t 0 ) | {z } L ⊕ . . . ⊕ D ( η ) ( t η 1 − t η 0 ) | {z } L ⊕ n L ∞ i = η +1 t i l − t i l − 1 i !  − | ˜ H ( u ) | i , | ˜ H ( u ) | i  ⊕ L ∞ i = η +1 t i l − 1 − t i l − 2 i !  − | ˜ H ( u ) | i , | ˜ H ( u ) | i  ⊕ . . . ⊕ L ∞ i = η +1 t i 1 − t i 0 i !  − | ˜ H ( u ) | i , | ˜ H ( u ) | i  | {z } L    u ∈ U o  e 1 ⊕ [ − 1 , 1 ] n n ˜ ǫ norm ( u )    u ∈ U o ( t l − t l − 1 ) ⊕ [ − 1 , 1 ] n n ˜ ǫ norm ( u )    u ∈ U o ( t l − 1 − t l − 2 ) ⊕ . . . ⊕ [ − 1 , 1 ] n n ˜ ǫ norm ( u )    u ∈ U o ( t 1 − t 0 ) | {z } P . (28) The summ ation symbols indicate th at the terms written in on e column are summed up. Since the expressions t j k − t j k − 1 are positive scalars, the fo llowing statement can b e used for the summation: For any two positi ve s calars a, b ∈ R + and the conv ex set S , one can state that a S ⊕ b S = ( a + b ) S . From this follows that D ( j ) ( t j 1 − t j 0 ) ⊕ D ( j ) ( t j 2 − t j 1 ) ⊕ . . . ⊕ D ( j ) ( t j l − t j l − 1 ) = D ( j ) ( t j l − t j 0 ) = D ( j ) t j l = D ( j ) δ j (29) and similarly for the T ay lor remaind e r terms ∞ M i = η +1 t i 1 − t i 0 i !  − | ˜ H ( u ) | i , | ˜ H ( u ) | i  ⊕ . . . ⊕ ∞ M i = η +1 t i l − t i l − 1 i !  − | ˜ H ( u ) | i , | ˜ H ( u ) | i  = ∞ M i = η +1 t i l − t i 0 i !  − | ˜ H ( u ) | i , | ˜ H ( u ) | i  = ˜ E ( δ, u ) (30) and the Krylov er ror terms [ − 1 , 1 ] n n ˜ ǫ norm ( u )    u ∈ U o ( t 1 − t 0 ) ⊕ . . . ⊕ [ − 1 , 1 ] n n ˜ ǫ norm ( u )    u ∈ U o ( t l − t l − 1 ) = [ − 1 , 1 ] n n ˜ ǫ norm ( u )    u ∈ U o δ. (31) 14 Inserting (29) - (31) into (2 8) y ields ˆ R p,l ( δ ) = P  D (1) δ ⊕ . . . ⊕ D ( η ) δ η ⊕ n ˜ E ( δ, u )    u ∈ U o  e 1 ⊕ [ − 1 , 1 ] n n ˜ ǫ norm ( u )    u ∈ U o δ. Since the above resu lt is independen t of the number of intermediate time steps l , we can choose l → ∞ , meanin g that lim l →∞ ˆ R p,l ( δ ) = R p ( δ ) . A P P E N D I X B E R RO R B O U N D F O R T H E A P P ROX I M A T I O N O F M A T R I X E X P O N E N T I A L S In this work, we u se the a-po steriori bound f r om [58, Eq. 4.1]. A method f or ef ﬁcient compu tation of this b ound is presented in the subseq uent p roposition . Proposition B.1 The err o r bound ∀ t ∈ [0 , t f ] :   e C t v − k v k V e H t e 1   ≤ k v k ǫ norm t in (3) hold s for ǫ norm = h ξ +1 , ξ ω ϕ, ϕ = ( e ν ( C ) t f − 1 ν ( C ) t f , for ν ( C ) > 0 , 1 , otherwise . ω = max k ∈{ 1 ,... , ⌈ t f /δ ⌉} ( ω k ) , ω k = e T ξ e H t k e H [0 ,δ ] e 1 , H [0 , δ ] = { 0 . 5 H δ + β 0 . 5 H δ | β ∈ [ − 1 , 1] } , wher e H [0 , δ ] can be r epres ented as a matrix zonotope, i.e., a zonoto pe who se center and generators are replaced by matrices; th e exponential of a matrix zonoto p e e H [0 ,δ ] can b e compu te d a s pr esented in [73]. Furthe r , ν ( C ) := λ max (0 . 5( C + C ∗ )) r eturn s the lar gest eigen va lu e and C ∗ is the co njugate transpo se o f C . Please note tha t the co mputation of ν ( C ) ca n b e done efﬁciently since 0 . 5( C + C ∗ ) is a r ea l, spa rse , symmetric matrix for which ef ﬁcient eigenvalue computatio ns [74], [75] exis t. Results can be o btained up to machine precision [76] a nd err or bou n ds are d erived in [77], [78]. F inally , e ξ is the unit vector whose ξ th value is one and the values of h ξ +1 , ξ and H ar e taken fr om Alg. 1. Pr oof: The a-p osteriori bound from [58, Eq. 4.1 ] is   e C t f v − k v k V e H t f e 1   ≤k v k h ξ +1 , ξ max t ∈ [0 ,t f ] | e T ξ e H t e 1 | e ν ( C ) t f − 1 ν ( C ) . (32) Please n ote that some sign s are chang e d compa r ed to [5 8, Eq. 4.1 ] since in that work, the system is speciﬁed as ˙ x = e − C t instead of ˙ x = e C t . a) E rr or bound is linear in ti me: The error boun d can be en closed by a fu nction that is linear in time. Fr om (32) we can derive th e following bo und: ∀ t ∈ [0 , t f ] :   e C t v − k v k V e H t e 1   ≤k v k h ξ +1 , ξ max ˆ t ∈ [0 ,t ]  | e T ξ e H ˆ t e 1 |  e ν ( C ) t − 1 ν ( C ) = t k v k h ξ +1 , ξ max ˆ t ∈ [0 ,t ]  | e T ξ e H ˆ t e 1 |  e ν ( C ) t − 1 ν ( C ) t ≤ t k v k h ξ +1 , ξ max ˆ t ∈ [0 ,t f ]  | e T ξ e H ˆ t e 1 |  max ˆ t ∈ [0 ,t f ]  e ν ( C ) ˆ t − 1 ν ( C ) ˆ t | {z } =: ϕ ( ˆ t )  (33) Dependin g on ν ( C ) , the f u nction ϕ ( t ) is monoton ically in- creasing or decre a sin g. T o show this, we insert the T aylor series e ν ( C ) t = P ∞ i =0 ( ν ( C ) t ) n n ! in the deriv ativ e ˙ ϕ ( t ) = ν ( C ) te ν ( C ) t − e ν ( C ) t + 1 ν ( C ) t 2 so that we obtain ˙ ϕ ( t ) = ν ( C ) t (1 + ν ( C ) t + ( ν ( C ) t ) 2 2! + . . . ) ν ( C ) t 2 − (1 + ν ( C ) t + ( ν ( C ) t ) 2 2! + . . . ) − 1 ν ( C ) t 2 = ( ν ( C ) t ) 2 + ( ν ( C ) t ) 3 2! + . . . ) − ( ν ( C ) t ) 2 2! − ( ν ( C ) t ) 3 3! − . . . ν ( C ) t 2 = ∞ X i =0  1 ( n + 1 )! − 1 ( n + 2)!  ( ν ( C ) t ) n +2 ν ( C ) t 2 = ν ( C ) ∞ X i =0  1 n + 1 − 1 ( n + 2 )( n + 1)  ( ν ( C ) t ) n n ! ∈ ν ( C ) ∞ X i =0 ]0 , 0 . 5] ( ν ( C ) t ) n n ! = ν ( C )]0 , 0 . 5 ] ∞ X i =0 ( ν ( C ) t ) n n ! = ν ( C ) ]0 , 0 . 5 ] e ( ν ( C ) t ) | {z } > 0 . Thus, ϕ ( t ) is mon otonically increasin g if ν ( C ) > 0 , monoto n- ically decreasing for ν ( C ) < 0 , a nd co n stant for ν ( C ) = 0 . As a consequ ence, the m a ximum is obtained at t = t f for ν ( C ) > 0 an d at t = 0 f o r ν ( C ) ≤ 0 . For the evaluation of t = 0 , we r equire the rule of L ’H ˆ opital: lim t → 0 e ν ( C ) t − 1 ν ( C ) t = lim t → 0 ν ( C ) e ν ( C ) t ν ( C ) = lim t → 0 e ν ( C ) t = 1 . The maximum is then obtain e d by ϕ = max ˆ t ∈ [0 ,t f ] ϕ ( ˆ t ) = ( e ν ( C ) t f − 1 ν ( C ) t f , for ν ( C ) > 0 , 1 , otherwise . (34) b) Com p utation of ma x ˆ t ∈ [0 ,t f ] | e T ξ e H ˆ t e 1 | : T o com p ute max ˆ t ∈ [0 ,t f ] | e T ξ e H ˆ t e 1 | , we compute e H τ k = e H t k e H [0 ,δ ] 15 as presented in [73, Eq. 4.10] u sing the matrix zo notope (see [73, Eq. 4.8 ]) H [0 , δ ] = { 0 . 5 H δ + β 0 . 5 H δ | β ∈ [ − 1 , 1] } . Let us also introdu ce ω k = e T ξ e H τ k e 1 . Since ω k is o ne-dime nsional, the result is the in terval max t ∈ τ k ( e T ξ e H t e 1 ) = sup( ω k ) . For th e entire time horizon, the maxim u m is simply ob tained as ω = max t ∈ [0 ,t f ] | e T ξ e H t e 1 | = max k ∈{ 1 ,..., ⌈ t f /δ ⌉} ( ω k ) . (35) c) F ina l Re sult: using (35) and ( 34) in (33) results in the propo sition: ∀ t ∈ [0 , t f ] :   e C t v − k v k V e H t e 1   ≤ t k v k h ξ +1 , ξ ω ϕ. As a short ﬁnal remark, we would like to m ention th at we have used the MA TLAB function eigs to ob ta in λ max . R E F E R E N C E S [1] A. Banerjee, K. K. V enkatasubra manian, T . Mukherje e, and S. K. S. Gupta, “Ensuring s afety , security , and sustainabilit y of mission-criti cal cybe r-physi cal systems, ” Proce edings of the IEEE , vol. 100, no. 1, pp. 283–299, 2012. [2] E. Asarin, T . Dang, G. Frehse, A. Girard, C. Le Guernic, and O. Maler , “Rece nt progress i n c ontinuou s and hybrid reachab ility analysis, ” in Pr oc. of the IEEE Confere nce on Computer Aided Contr ol Systems Design , 2006, pp. 1582–1587. [3] I. Kol manovsk y and E. G. Gilber t, “Theory and computat ion of dis- turbanc e in vari ant sets for discrete-t ime linear systems, ” Mathemati cal Pr oblems in Engineering , vol. 4, pp. 317–367, 1998. [4] F . Blanchini, “Set in vari ance in control, ” Automat ica , vol . 35, no. 11, pp. 1747 – 1767, 1999. [5] A. El-Guindy , D. Han, and M. Althof f, “Estimating the re gion of attrac tion via forward reacha ble sets, ” in Pr oc. of the American Contr ol Confer ence , 2017, pp. 1263–1270. [6] D. Limon, I. Alvarad o, T . Alamo, and E. F .Camacho, “Robust tube- based MPC for tracking of constrained linear syste ms with additi ve disturban ces, ” J ournal of Pro cess Contr ol , vol. 20, no. 3, pp. 248–260, 2010. [7] B. Sch ¨ urmann and M. Althof f, “Con vex interpol ation control with formal guarant ees for disturbed and constrai ned nonlinear systems, ” in Proc . of Hybrid Systems: Computation and Contr ol , 2017, pp. 121–130. [8] C. Combastel, “ A state boun ding observe r for uncerta in non-line ar continu ous-time systems based on zonotopes, ” in P roc. of the 44th IEEE Confer ence on Decision and Contr ol, and the Europea n Contro l Confer ence , 2005, pp. 7228–7234. [9] V . T . H. Le, C. St oica, T . Alamo, E. F . Camacho, and D. Du mur , Zonotopes: F r om Guaranteed State-estimat ion to Contr ol . W iley , 2013. [10] H. Roehm, J. Oehlerking, M. W oehrle , and M. Althof f, “Reachset conformanc e testin g of hybrid automata, ” in Proc. of Hybrid Systems: Computati on and Contr ol , 2016, pp. 277–286. [11] P . Falco ne, M. Ali, and J. Sj ¨ oberg , “Predic ti ve threat assessment via reacha bilit y analysis and set in varian ce theory , ” IEEE T ransact ions on Intell ige nt T ransportati on Systems , vol. 12, no. 4, pp. 1352–1361, 2011. [12] M. Althof f and J. M. Dolan, “Online veriﬁca tion o f automated road vehi cles using reachabi lity analysis, ” IEEE T ransact ions on Robotics , vol. 30, no. 4, pp. 903–918, 2014. [13] H. T ¨ aubig, U. Frese, C. Hertzber g, C. L ¨ uth, S. Mohr, E. V orobe v , and D. W alter , “Guarantee ing functi onal safety: design for prov ability and computer -aided ve riﬁcati on, ” Au tonomous R obots , vol. 32, no. 3, pp. 303–331, 2012. [14] A. Pereira and M. Althof f, “Overapp roximati ve arm occupanc y predic- tion for human-robot intera ction bui lt from archetypa l move ments, ” in Internati onal Confere nce on R obotics and Automati on , 2016, pp. 1394– 1401. [15] H. N. V ille gas Pico a nd D. C. Aliprantis, “ V oltage ride-t hrough ca- pabili ty veriﬁc ation of wind turbines with fully- rated con verters using reacha bilit y analysis, ” IEEE T ransac tions on Energy Con version , vol. 29, no. 2, pp. 392–405, 2014. [16] A. El-Guindy , D. Han, and M. Althoff, “Formal analysis of drum-boiler units to maximize the load-foll o wing capabili ties of powe r plants, ” IEEE T ransaction s on P ower Systems , vol. 31, no. 6, pp. 4691–4702, 2016. [17] G. Frehse, B. H. Krogh, and R. A. Rutenbar , “V erifying analog oscillat or circui ts using forw ard/bac kward abstrac tion reﬁne ment, ” in Pr oc. of Design, Automation and T est in Eur ope , 2006. [18] M. Althof f, A. Rajhans, B. H. Krogh, S. Y aldiz , X. Li, and L. Pile ggi, “Formal veriﬁcati on of phase-lock ed loops using reachabil ity analysis and continuiza tion, ” Communicati ons of the ACM , vol. 56, no. 10, pp. 97–104, 2013. [19] T . Dang, “V ´ eriﬁcati on et synth ` ese des syst ` emes hybrides, ” Ph. D. disser- tatio n, Institut National Polytechniq ue de Grenoble , 2000. [20] A. A. Kurzhanski y and P . V araiya, “Ellipsoidal techniques for reacha- bility analysi s of discrete-ti me linear systems, ” IEEE T ransactions on Automat ic Contr ol , vol. 52, no. 1, pp. 26–38, 2007. [21] A. Chutinan and B. H. Krogh, “Computational techniques for hybrid system veri ﬁcati on, ” IEEE T ransactio ns on Automati c Contr ol , vol. 48, no. 1, pp. 64–75, 2003. [22] A. Girard, “Reac habili ty of uncert ain linear systems using zonotopes, ” in Hybrid Systems: Computati on and Contr ol , ser . L NCS 3414. Springer , 2005, pp. 291–305. [23] M. Althof f and B. H. Krogh, “Zonotope bundle s for the ef ﬁcient computat ion of reacha ble sets, ” in Proc. of the 50th IE EE Confer ence on Decision and Contr ol , 2011, pp. 6814–6821. [24] A. Girard and C. Le Guernic, “Efﬁcien t reachabilit y analysis for linear systems using support functio ns, ” in Pr oc. of the 17th IF AC W orld Congr ess , 2008, pp. 8966–8971. [25] I. M. Mitche ll, A. M. Bayen, and C. J. T om lin, “ A time-dep endent Hamilton– Jacobi formulation of reacha ble s ets for continuous dynamic games, ” IEEE T ransac tions on Automatic Contr ol , vol. 50, pp. 947–957, 2005. [26] M. Althof f and G. Frehse, “Combining zonotopes and support functions for efﬁc ient reacha bilit y analysis of linear systems, ” in Proc. of the 55th IEEE Conferen ce on Decision and Contr ol , 2016, pp. 7439–7446. [27] P . S. Duggirala and M. V iswanat han, “Parsimonious, simulation based veri ﬁcatio n of linear systems, ” in Proc. of Internationa l Confer ence on Computer Aided V eriﬁcat ion , 2016, pp. 477–494. [28] S. Bak and P . S . Duggirala , “HyLAA: A tool for computing s imulati on- equi vale nt reachabi lity for linea r systems, ” in Pr oc. of the 20th Interna- tional Confer ence on Hybrid Systems: Computati on and Contr ol , 2017, pp. 173–178. [29] S. Bak, H.-D. Tran, and T . T . J ohnson. (2018) Numerical veriﬁcati on of af ﬁne systems with up to a billio n dimensions. arXiv: 1804.01583v2. [30] S. Bak and P . S. Duggirala , “Simulation -equi valent reachabili ty of large linea r syste ms with i nputs, ” in Pr oc. of I nternati onal Confe re nce o n Computer Aided V eriﬁcat ion , 2017, pp. 401–420. [31] S. Bogomolov , M. Forets, G. Frehse, F . V iry , A. Podelski, and C. Schilling, “Reach set approximat ion through decompositio n with low- dimensiona l sets and high-di mensional matrice s, ” in Proc . of the 21st Internati onal Confere nce on Hybrid Systems: Computation and Contr ol , 2018, pp. 41–50. [32] G. Lafferr iere, G. J. Pappas, and S. Y ovine, “ A ne w class of decidable hybrid systems, ” in Hybrid Systems: Computation and Contr ol , ser . LNCS 1569. Springer , 1999, pp. 137–151. [33] E. Asarin, T . Dang, and A. Girard, “Hybridi zatio n metho ds for the analysi s of nonlinear systems, ” A cta Informatic a , vol. 43, pp. 451–476, 2007. [34] M. Althof f, O. Stursberg , and M. Buss, “Reacha bilit y analysis of nonlin- ear systems with uncerta in parameters using conserv ativ e linearizat ion, ” in Proc. of the 47th IEEE Confer ence on Decision and Contr ol , 2008, pp. 4042–4048. [35] T . Dang, O. Maler , and R. T estyli er , “ Accurate hybridiza tion of nonlinear systems, ” in Hybrid Systems: Computation and Contr ol , 2010, pp. 11– 19. [36] C. Mole r and C. V . Loan, “Nin eteen dubious w ays to compute the expo nentia l of a matrix, twenty-ﬁ ve years late r , ” SIAM Revie w , vol. 45, no. 1, pp. 3–49, 2003. [37] Y . Saa d, Iterati ve M ethods for Sparse Linear Systems . Society for Industria l and Applied Mathemati cs, 2003. [38] A. C. Antoula s, A pproxi mation of Larg e-Scal e Dynamical Systems . Society for Industrial and Applied Mathema tics, 2005. [39] W . H. A. Schilders, H. A. va n der V orst, and J. Rommes, Model Or der Reductio n: Theory , Researc h As pects and Applicat ions . Springer , 2008. 16 [40] Z. Han and B. Krogh, “Reachabi lity analysis of hybrid control systems using reduced-orde r models, ” in Pr oc. of the American Contr ol Confer- ence , 2004, pp. 1183–1189. [41] G. Frehse, C. L. Guernic, A. Donz ´ e, S. Cotton, R. Ray , O. Lebeltel, R. Ripado, A. Girard, T . Dang, and O. Maler , “SpaceEx: Scalable veri ﬁcatio n of hybrid systems, ” in Pro c. of the 23r d Internationa l Confer ence on Computer A ided V eriﬁc ation , ser . LNCS 6806. Springer , 2011, pp. 379–395. [42] X. Chen, E. ´ Abrah ´ am, and S. Sankaranar ayanan, “Flo w*: An analyzer for non-linea r hybrid systems, ” in P r oc. of Computer-Aided V eriﬁcation , ser . LNCS 8044. Spri nger , 2013, pp. 258–263. [43] A. Gurung, A. Deka, E. Bartocci, S. Bogomolov , R. Grosu, and R. Ray , “Para llel reachabilit y analysis for hybrid systems, ” in Pr oc. of the ACM /IEEE Internationa l Confer ence on F ormal Methods and Models for System Design , 2016, pp. 12–22. [44] M. Althoff, “ An introduction to CORA 2015, ” in Pr oc. of the W orkshop on A pplied V eriﬁca tion for Continuo us and Hybrid Systems , 2015, pp. 120–151. [45] M. Althof f, S. Bak, D. Catta ruzza, X. Chen, G . Frehse, R. Ray , and S. Schupp, “ARCH-COMP17 cate gory report: Continuous and hybrid systems with linea r continuous dynamics, ” in Pro c. of the 4th Inter- national W orkshop on Applied V eriﬁcat ion for Continuous and Hybrid Systems , 2017, pp. 143–159. [46] M. Althof f, S. Bak, X. Chen, C. Fan , M. Fore ts, G. Frehse, N. Kochdumper , Y . Li, S. Mitra, R. Ray , C. Schilling, and S. Schupp, “ARCH-COMP18 catego ry report: Continuou s and hybrid systems with linea r continuous dynamics, ” in Proc . of the 5th International W orkshop on A pplied V eriﬁca tion for Continuo us and Hybrid Systems , 2018, pp. 23–52. [47] Z. Han and B. H. Krogh, “Reachabi lity analysis of large -scale af ﬁne systems using lo w -dimensional polytopes, ” in Hybrid Systems: Compu- tation and Contr ol , ser . LNCS 3927. Spri nger , 2006, pp. 287–301. [48] H.-D. T ran, L . V . Nguyen, W . Xiang, and T . T . Johnson, “Order - reducti on abstractions for safety veri ﬁcation of high-dimensiona l linear systems, ” Discrete Event Dynamic Systems , vol. 27, no. 2, pp. 443–461, 2017. [49] Y . Chou, X. Chen, and S. Sankarana rayana n, “ A s tudy of m odel-orde r reducti on techniques for veriﬁc ation, ” in Pro c. of Numerical Softwar e V eriﬁcation , 2017, pp. 98–113. [50] A. J. v an der Schaft, “Equi valen ce of dynamical systems by bisimu- latio n, ” IEEE Tr ansacti ons on Automatic Contr ol , vol. 49, no. 12, pp. 2160–2172, 2004. [51] A. Girard and G. J. Pappas, “ Approximate bisimulation : A bridge betwee n computer sc ience and control theory , ” Eur opean J ournal of Contr ol , vol. 17, no. 5-6, pp. 568–578, 2011. [52] J. G. V anAntwe rp and R. D. Braatz, “ A tutorial on linear and bilinear matrix inequali ties, ” Journal of Proc ess Contr ol , vol. 10, pp. 363–385, 2000. [53] A. Girard and G. J. Pappas, “ Approximate bisimulation relations for constrai ned linear systems, ” Automatic a , vol. 43, pp. 1307–1317, 2007. [54] M. Althof f, O. Stursber g, and M. Buss, “Computing reachable sets of hybrid systems using a combination of zonoto pes and polytope s, ” Nonline ar Analysis: Hybrid Systems , vol. 4, no. 2, pp. 233–249, 2010. [55] Y . Saad, “ Analysis of some Kr ylov subspa ce approxi mations to the matrix exponent ial operat or , ” SIAM Journal on Numerical Analysis , vol. 29, no. 1, pp. 209–228, 1992. [56] M. Hochbruc k and C. Lubich, “On Krylo v subspace approximatio ns to the matrix exponent ial operat or , ” SIAM Journal on Numerical Analysis , vol. 34, no. 5, pp. 1911–1925, 1997. [57] R. B. Sidje, “Expokit: A software package for computin g matrix expo nentia ls, ” A CM T ransactions on Mathemat ical Softwar e , v ol. 24, no. 1, pp. 130–156, 1998. [58] Z. Jia a nd H. L v , “ A posteriori erro r estimates of Krylov subspace approximat ions to matrix functions, ” Numerical Algori thms , vol . 69, no. 1, pp. 1–28, 2015. [59] H. W ang and Q. Y e, “Error bounds for the Krylov subspace m ethods for computat ions of m atrix exponent ials, ” SIAM Journal on Matrix Analysis and Applicat ions , vol. 38, no. 1, pp. 155–187, 2017. [60] W . K ¨ uhn, “Rigorously compute d orbits of dynamical systems without the wrapping ef fect, ” Computing , vol. 61, pp. 47–67, 1998. [61] M. Althof f, “Reac habili ty analysis and its a pplic ation to the safety assessment of autonomous cars, ” Disserta- tion, T echnisc he Univ ersit ¨ at M ¨ un chen, 2010, http ://nbn- resolving .de/urn/re solve r .pl?urn:nbn:de:bvb:91-diss-20100715-963752- 1-4. [62] G. Laf ferriere, G. J. Pa ppas, and S. Y ovi ne, “Symbolic re achabi lity computat ion for famili es of linear vec tor ﬁelds, ” Symbolic Computation , vol. 32, pp. 231–253, 2001. [63] M. Althoff, C. Le Guernic, and B. H. Krogh, “Reac hable set compu- tatio n for uncertain time-v arying linear systems, ” in Hybrid Sy stems: Computati on and Contr ol , 2011, pp. 93–102. [64] L. Jaulin, M. Kief fer , and O. Didrit, Applied Interval Analysis . Springer , 2006. [65] A. Girard , C. Le Guern ic, and O. Ma ler , “Ef ﬁ cient computation of reacha ble sets of line ar time-in varia nt systems w ith inputs, ” in Hybrid Systems: Computation and Contr ol , ser . LNCS 3927. Springer , 2006, pp. 257–271. [66] M. Althof f and B. H. Krogh, “Reacha bility analysis of nonlinea r dif ferent ial-a lgebraic syste ms, ” IEEE T ransactions on Automatic Con- tr ol , vol. 59, no. 2, pp. 371–383, 2014. [67] S. Will iams, L. Oliker , R. V uduc, J. Shalf, K. Y elic k, and J. Demmel, “Optimiz ation of sparse matrix-vect or multiplicat ion on emerging mul- ticore platforms, ” P arall el Computing , vol. 35, pp. 178–194, 2009. [68] H.-D. Tran, L. V . Nguyen, and T . T . Johnson, “Large -scale linear systems from order-reduct ion, ” in Proc. of ARCH16. 3rd International W orkshop on A pplied V eriﬁca tion for Continuo us and Hybrid Systems , 2017, pp. 60–67. [69] O. H. Amman, T . von K ´ arm ´ an, and G. B. W oodruff , “The failure of the T acoma Narrows Bridge, ” Federal W orks Agency , T ech. Rep., 1941. [70] Y . Chahlaoui and P . V an Dooren, “ A collect ion of benchmark exa m- ples for model reduction of linear time in v arian t dynamical systems, ” Uni ve rsity of Manchester , T ech. Rep., 2002. [71] M. Pisaroni, “Dynamic response of a 3D s tructure: Roose velt lake bridge, ” June 2012, T echnic al Univ ersity of Munich. [72] A. G. Dav enport and J. P . C. King, “The inﬂuence of topography on the dynamic wind loading of long s pan bridges, ” J ournal of W ind Engineerin g and Industrial Aero dynamics , vol. 36, no. 2, pp. 1373–1382, 1990. [73] M. Al thof f, B. H. K rogh, and O. Stursberg, M odelin g, Design, and Simulati on of Systems with Uncertainti es . Springer , 2011, ch. Analyzing Reacha bilit y of Linear Dynamic Systems with Parametric Uncertain ties, pp. 69–94. [74] R. B. Lehoucq and D. C. Sorensen, “Deﬂation techni ques for an implici tly re-started Arnoldi iteration, ” SIAM Journal on Matrix Analysis and Applicat ions , vol. 17, no. 4, pp. 789–821, 1996. [75] D. Sorensen, “Implicit appli cation of polynomia l ﬁlters in a k-step Arnoldi method, ” SIAM Journal on Matrix A nalysis and Applicati ons , vol. 13, no. 1, pp. 357–385, 1992. [76] H. D. Simon, “ A nalysis of the symmetric Lanczos algorithm with re- orthogona lirat ion methods, ” Linear Algebra and its Applicat ions , vol. 61, pp. 101–131, 1984. [77] C. C. Paige, “ Accurac y and effe cti veness of the Lanczos al gorithm for the symmetric eigenp roblem, ” Linear Algebr a and its Applications , vol. 34, pp. 235–258, 1980. [78] ——, “Error analysis of the Lanczos algorit hm for tridiagonal izing a symmetric matrix, ” IMA Journal of A pplied Mathematic s , vol. 18, no. 3, pp. 341–349, 1976. Matthias Althoff is an assistant professor in com- puter sc ience at T echni sche Uni versi t ¨ at M ¨ unchen , Germany . He recei ved his diploma engineering de- gree in Mech anical Engineeri ng in 2005, and his Ph.D. degree in Electric al E ngineer ing in 2010, both from T echni sche Univ ersit ¨ at M ¨ unchen, Germany . From 2010 to 2012 he was a postdoctora l researche r at Carnegie Mellon Unive rsity , Pittsbur gh, USA, and from 2012 to 2013 an assistant professor at T ech- nische Univ ersit ¨ at Ilmenau, Germany . His research intere sts incl ude formal veriﬁca tion of continuous and hybrid systems, reachabili ty analysi s, planning algorithms, nonlinea r control , automated vehicle s, and power systems.

Reachability Analysis of Large Linear Systems with Uncertain Inputs in the Krylov Subspace

Original Paper

Comments & Academic Discussion

Leave a Comment

Original Paper

Related Papers

Comments & Academic Discussion

Leave a Comment