Uncertainty Quantification in Hybrid Dynamical Systems

Uncertain t y Quan tification in Hybrid Dynamical Systems T uhin Sahai † Jos ´ e Miguel P asini † † United T ec hnologies Researc h Cen ter, 41 1 Silv er Lane, East Hartford, CT 06108, USA Abstract Uncertaint y quantification (U Q ) techn iques are frequently used to ascertain output v ariability in systems with parametr ic uncerta int y . T ra ditional algorithms for UQ are either system-agnostic and slow (such as Monte Car lo) o r fast with str ingent assumptions on smo oth- ness (suc h as po lynomial chaos and Quasi- Mo nt e Carlo). In this work, we develop a fast UQ appro ach for hybrid dyna mical sys tems by ex- tending the po lynomial chaos metho dolog y to these systems. T o cap- ture disco nt inu ities, w e use a wav elet-based Wiener-Haar expansion. W e develop a b oundar y lay er approach to pro pagate uncertaint y through sep ar able reset co nditions. W e also in tr o duce a transp ort theory based approach for propaga ting uncertaint y throug h hybrid dynamical sys- tems. Here the expansion yields a set o f hyp e rb olic equations that are solved by integrating along characteristics. The so lution of the par- tial differential equation along the characteristics allows one to quan- tify uncertaint y in h y brid or switching dynamical systems. The ab ove metho ds are demonstrated on example pro blems. 1 In tro d u ction Uncertain t y Quan tification (UQ) is an area of mathematics that is used to quan tify output distr ibutions giv en parametric un certain t y . T raditio n al ap- proac hes include Mon te Carlo and Quasi-Mon te Carlo m etho ds [1], resp onse surface metho ds [2, 3] as we ll as p olynomial c h aos and probabilistic colloca- tion based approac hes [4]. The p olynomial c haos a p proac h f or u ncertain ty quan tification w as originally pr op osed b y Norb ert Wiener [5]. Assuming that one is given inpu t un certain t y in the form of distributions a ss o ciated 1 with v arious parameters of th e system, p olynomial chaos/ p robabilistic col- lo cation metho ds provide an approac h for fast uncertain ty quantifica tion under the assumption of smo oth dyn amics. In particular, p olynomial c haos pro vid es exp onential con v ergence for smo oth systems and pro cesses with finite v ariance [4]. P olynomial c haos based metho ds ha ve b een used for a m ultitud e of applications, see [6 , 7, 8, 9, 10, 11, 12, 13] for examples. Note that, dep en ding on the ap p lication, one can com b ine v arious UQ approac hes. F or example, a com bin ation of p olynomial c haos and the resp onse su rface metho dology h as b een used to dev elop probabilistic collo cation metho ds for discrete d istributions in [10]. In this w ork, we f o cus on dev eloping UQ tec hniques for h ybr id dynamical systems. Hybrid d y n amical systems theory is used to mod el systems with b oth discrete and con tin u ous d ynamics [14]. E x amp les include the b ouncing ball automaton [15], biologica l net works [16, 17], air traffic managemen t systems [18], comm u nication net works [19], elev ators, and rob otics, to name a few. These systems frequently display ric h d ynamics not seen in contin uous systems. F or example, Zeno b ehavio r in hybrid systems is characte r ized by an infi nite num b er of discrete switc h es in finite time [15, 20]. Hybrid s ystems can b e particularly c hallenging from an analysis standp oint sin ce traditional tec hniques, suc h as p olynomial c haos based m etho ds, assume smoothness, rendering them inapp licable. In this w ork, w e deve lop p olynomial chao s and transp ort theory based metho ds for propagating u n certain t y through h ybr id systems. W e assume that the domains asso ciated with different mo des of op eration of the hybrid system do not o verlap. W e demons trate that, by in tegrating o ver appr opri- ate time-v aryin g regio n s, one can e xtend the p olynomial c haos framew ork to hybrid dynamical systems. W e r esolv e t he issue of state r esets [14] in the separable case b y using b oun d ary la y er appr o ximations. T o capture the disco n tinuities in the p robabilit y distributions of the outpu t v a r iables, w e us e a Haar-w a velet expansion [21]. This expansion has previously b een used in th e p olynomial chaos setting to pr opagate u ncertain ty through dy- namical systems close to bifurcation points [9, 2 2, 23]. Here w e dev elop a metho dology to p ropagate uncertain ty through systems with discontin uities in d ynamics and output along with state r esets. W e also dev elop a trans- p ort th eory based approac h that allo ws one to propagate th e uncertaintie s through th e v arious mo des of the h ybrid d y n amical system. Our pap er is organized as follo ws: in s ection 2 w e define hybrid d ynam- ical sy s tems and the pr ob lem of uncertain ty qu an tification. In section 3 w e first construct the framework for p olynomial chaos in the h y b rid dynamical system setting (3.1). W e then d emons trate th e Haar wa v elet expansion for 2 h yb rid p olynomial chaos in 3.2. The handlin g of sta te resets is considered in section 3.3. Finally , the results on hybrid p olynomial c h aos are pr esented in 3.4. The transp ort operator theo r y based method for propagati ng un- certain t y through hybrid dyn amical systems is d ev elop ed in se ction 4 and conclusions are d r a wn in section 5. 2 Problem definition Let S = ( q , X, f , x (0) , D , E , G, R ) denote a hybrid system S , where q Set of discrete v ariables X Set of con tinuous v ariables f : X × q → T X V ector field x (0) Set of initial conditions D : q → P ( X ) Domain E Set of discrete transitions G : E → P ( X ) Guard cond itions R : E × X → P ( X ) Reset map. In the ab o ve table, T X is the tangen t bun dle of X and P ( X ) is the p o wer set of X . F or more details on the definition of hybrid systems, s ee [14]. W e can use the follo wing r epresen tation for hybrid systems, ˙ x = f ( x, λ, q ) , (1) where x ∈ X is a v ector of state v ariables and the form of f ( x, λ, q ) is dic- tated by q , which represents the mo d e of op eration of the hybrid d y n amical system. The discrete state q is d etermined by the guard conditions that dictate transitions b et we en mo des (see Fig. 1). The reset functions h i ( x ) are a p art of the reset map R . In Eqn . 1 let λ denote the v ector of system parameters an d x (0) the initial condition for the system. If the system parameters λ in Eqn. 1 are un certain (i.e., eac h λ i has an asso ciated distribution) then one t yp ically desires to q u an tify the time- v arying moments (suc h as mean and v ariance) of x ( t ) (note that x (0) may also b e uncertain). As men tioned earlier, although one can use Mon te Carlo based samplin g metho ds [1], they are plagued b y slo w conv ergence. I n par- ticular, the mean is exp ected to con verge as 1 / √ N , where N is t he n um- b er of samples. Quasi-M on te Carlo based sampling method s are exp ected to gi v e a con vergence rate of l og d ( N ) / N , where d is th e d imensionalit y of the random space [24], making these metho d s attractiv e for problems in lo w dimensions. P olynomial c haos based metho ds p r o vide an alternativ e framew ork for uncertain t y quanti fi cation with exp onenti al con vergence f or 3 ˙ x = f ( x, λ ) ˙ x = g ( x, λ ) x + = h 1 ( x − ) x + = h 2 ( x − ) q ′ x ∈ G ( q , q ′ ) x ∈ G ( q ′ , q ) q Figure 1: Sc h ematic for hybrid (switc hin g) systems. pro cesses with fi nite v ariance, b ut they to o su ffer from the curse of dimen- sionalit y [4]. In the next section we extend the p olynomial c haos framew ork to hybrid systems. 3 P olynomial c haos for h ybrid dyn amical sy s tems Starting with a complete p robabilit y space Γ giv en by (Ω , F , P ), where Ω is the sample space, F is the σ -algebra on Ω and P is a pr obabilit y measure, let L 2 (Γ , X ) denote the Hilb ert sp ace of square-integ r able, Γ-measurable, X -v alued random element s . Then one can, in general, define a p olynomial c haos basis { H α ( λ ( ω )) } , wh ere λ ( ω ) is a random vec tor and α = ( α 1 , α 2 , . . . ) is a vect or of non-negativ e indices. W e denote th e probability d en sit y f unc- tion of the random v ector λ b y ρ ( λ ). Generalized p olynomial chaos (gPC) [25], pro vides a framew ork for r ep- resen ting second-order stoc hastic pro cesses r ∈ L 2 (Γ , X ) for a r bitrary d is- tributions of λ by the follo w in g expansion: r ( λ ) = ∞ X | α | =0 a α H α ( λ ) , (2) where | α | = P i α i is th e sum of the indices of α and H α ( λ ) are orthogonal p olynomials on Γ with r esp ect to ρ ( λ ), i.e. Z Γ ρ ( λ ) H α ( λ ) H β ( λ ) dλ = δ αβ , (3) where δ αβ is the K ronec ke r delta p ro duct. Dep ending on ρ ( λ ) one can gen- erate an app ropriate orthogo n al basis for representing r ( λ ). F or example, 4 if ρ is Gaussian, then the appropriate p olynomial chaos basis is the set of Hermite p olynomials; if ρ is the uniform distribution, then the basis is th e set of Legendre p olynomials. F or details on the corresp ondence b et ween d is- tributions and p olynomials see [4, 26]. A framew ork to generate p olynomials for arbitrary d istributions has b een dev elop ed in [25]. In practice, t h e expansion in Eqn. 2 is tru ncated at a particular order, sa y , p . One can then use Galerkin pr o jections to obtain a set of differentia l equations for th e co efficien ts a α in Eqn . 2 [4]. W e no w extend the standard p olynomial c h aos framewo r k to hybrid dy- namical s ystems d espite the presence of switc h ing an d state resets. T o the b est of our knowle d ge, it is the fi rst attempt to d ev elop to ols for f ast uncer- tain t y quant ifi cation f or this class of systems. 3.1 Hybrid P olynomial Chaos Without loss of generalit y , consider the f ollo w ing tw o-mo de hybrid d ynam- ical sys tem as repr esen tativ e of sy s tems in which th e different op erating mo des are asso ciated with non-o verlapping regions: ˙ x = ( f ( x, λ ) if x ≥ 0 g ( x, λ ) otherwise . (4) Here on e d esires to quantify x ( t ; λ ), i.e., d etermine x as a function of time t and parameters λ . The sys tem ab o ve has tw o mo des of op eration determined b y its state. One can parameterize these mo des in the follo w ing wa y: 1 R 1 ( x ) = ( 1 if x ≥ 0 0 otherwise (5) 1 R 2 ( x ) = 1 − 1 R 1 ( x ) . (6) When 1 R 1 ( x ) = 1 (corresp onding to x ≥ 0) the go v erning differential equa- tions are f ( x, λ ), and when 1 R 2 ( x ) = 1 ( x < 0) the go ve r ning differenti al equations are g ( x, λ ). T h u s, one can rewrite Eqn. 4 as, ˙ x = 1 R 1 ( x ) f ( x, λ ) + 1 R 2 ( x ) g ( x, λ ) . (7) This equation extends easily to k mod es of op eration b y constructing indi- cator fun ctions for eac h mode of op eration 1 R 1 , 1 R 2 , . . . , 1 R k of the hybrid system. W e no w expand x ( t ; λ ) in the appropriate o r thogonal p olynomial c haos b asis, x ( t ; λ ) = p X | α | =0 a α ( t ) H α ( λ ) . (8) 5 Dropping the argum en ts of a α ( t ) and H α ( λ ) for s implicit y and us ing the ab o ve relation with Eqn. 7, one gets p X | α | =0 ˙ a α H α = 1 R 1 ( p X | α | =0 a α H α ) f ( p X | α | =0 a α H α , λ ) + 1 R 2 ( p X | α | =0 a α H α ) g ( p X | α | =0 a α H α , λ ) . By multiplying the ab o ve relation b y ρ ( λ ) H k ( λ ), in tegrating o ver Γ, and using orthogonalit y conditions, w e get ˙ a k ( t ) = Z R 1 ( t ) f ( p X | α | =0 a α H α , λ ) ρ ( λ ) H k ( λ ) dλ + Z R 2 ( t ) g ( p X | α | =0 a α H α , λ ) ρ ( λ ) H k ( λ ) dλ. (9) Note that R 1 ( t ) = { λ : P p | α | =0 a α H α ≥ 0 } and R 2 ( t ) = Γ − R 1 ( t ). Thus b y ev aluating th e t wo in tegrals on e can ev olv e a k ( t ) for any index vecto r k . Note, ho w ever, that the regions of integrati on R 1 and R 2 are time-dep end en t quan tities and must b e ev aluated at ev ery in s tant in time. F or a pictorial depiction of R 1 and R 2 , see Fig. 2. 3.2 Hybrid P olynomial Chaos and w a velet expans ions Hybrid systems can disp la y d iscon tinuous b eha vior as a fun ction of the un- certain parameters. In view of this, a smo oth p olynomial chaos expans ion is exp ected to degrade as the discon tinuities b ecome more sev ere. In Ref. [22] the authors dev elop a w av elet-based Wiener-Haar expansion to treat b if u r- cating (bu t smo oth) dynamical s ystems with uncertain in itial conditions that resu lt in disconti nuous b ehavio r . I n this section we adapt the Wiener- Haar expansion to hybrid dynamical sy s tems. In [22, 23], output v ariables a r e expanded i n terms o f Wiener-Haar w av elets expressed as fun ctions of the Cu mulativ e Distribu tion F u nction (CDF) of the uncertain parameters. F or simplicit y , co n sider the univ ariate case. Here w e denote the CDF of the u n certain p arameter λ as u ( λ ) and 6 −4 −2 0 2 4 −5 0 5 10 15 20 λ P S f r a g r e p l a c e m e n t s P p | α | =0 a α ( t ) H α ( λ ) R 1 ( t ) R 1 ( t ) R 1 ( t ) R 2 ( t ) R 2 ( t ) Figure 2: Pictorial r epresen tation of the regions o f in tegration for the tw o in tegrals in Eqn. 9. I n this example, the regions are not simply connected. expand th e state v ector x as 1 x ( t ; λ ) = x 0 ( t ) + P X j =0 2 j − 1 X k =0 x j k ( t ) ψ j k ( u ( λ )) , (10) where, ψ j k ( u ) = 2 j / 2 ψ (2 j u − k ) , (with j = 0 , 1 , . . . and k = 0 , . . . , 2 j − 1) is a family of Haar wa v elets [21], defin ed in terms of the mother wavelet : ψ ( u ) =      1 0 ≤ u < 1 / 2 − 1 1 / 2 ≤ u < 1 0 otherwise. The in dex j determin es the scale of th e w av elet and k its d isplacemen t. Note that { ψ j k } is a family of orthonormal functions on the in terv al [0 , 1] with resp ect to th e un iform density . This makes the family { ψ j k ◦ u } automatically orthonormal w ith resp ect to the p robabilit y density of λ : δ j l δ k m = Z 1 0 ψ j k ( u ) ψ lm ( u ) du = Z ( ψ j k ◦ u )( λ ) ( ψ lm ◦ u )( λ ) ρ ( λ ) dλ. 1 F o r the multiv ariate case, see Ref. [22]. 7 Additionally , all ψ j k ’s are orthogonal to the constan t fun ction on [0 , 1], w hic h implies that the mean of x is giv en by the fir st term in the expansion x 0 ( t ) = Z x ( t ; λ ) ρ ( λ ) dλ and that the v ariance is σ 2 ( t ) = P X j =0 2 j − 1 X k =0 x 2 j k ( t ) . W e no w use this expansion on a switc hin g oscillator example: ¨ x + c ˙ x + x + λ = 0 if x ≥ 0 ¨ x + c ˙ x + x − λ = 0 if x < 0 , (11) whic h can b e rewritten as, ˙ x = y ˙ y = − cy − x − λ 1 R 1 ( x ) + λ 1 R 2 ( x ) . The expans ion in this case is x ( t ; λ ) = x 0 ( t ) + P X j =0 2 j − 1 X k =0 x j k ( t ) ψ j k ( u ( λ )) y ( t ; λ ) = y 0 ( t ) + P X j =0 2 j − 1 X k =0 y j k ( t ) ψ j k ( u ( λ )) . Pro jecting these equations on to the basis fu nctions yields, ˙ x 0 = y 0 ˙ y 0 = − cy 0 − x 0 − Z 1 0 λ ( u ) 1 R 1 ( x ( u )) du + Z 1 0 λ ( u ) 1 R 2 ( x ( u )) du ˙ x j k = y j k ˙ y j k = − cy j k − x j k − Z 1 0 λ ( u ) 1 R 1 ( x ( u )) ψ j k ( u ) du + Z 1 0 λ ( u ) 1 R 2 ( x ( u )) ψ j k ( u ) du Note that to compute λ ( u ) we inv ert the CDF u ( λ ). T o compute the in tegrals needed to ev olv e these equations numerically , w e tak e adv an tage of the fact that Haar w av elets are piecewise constant. 8 Namely , for a giv en truncation order P , if w e divide [0 , 1] into 2 P +1 equal subinterv als, b oth ψ j k (with j ≤ P ) and the truncated expansion for x (and therefore the in d icator functions) are constan t in eac h subinterv al . This implies that in eac h sub in terv al w e only need to calculate the integ r al of λ ( u ), whic h is kno wn a priori. F or the case of a Gaussian λ ∼ N ( µ, σ 2 ), we ha ve λ ( u ) = µ + σ √ 2 erf − 1 (2 u − 1) , whic h has a p rimitiv e, Z λ ( u ) du = µu − σ 1 √ 2 π exp  − [erf − 1 (2 u − 1)] 2  . Therefore the contribution to ˙ y j k of the inte grals in eac h subinterv al l = 0 , . . . , 2 P +1 − 1 is either zero or ± 2 j / 2 times the pr ecomputed v alue Z ( l +1) / 2 P +1 l/ 2 P +1 λ ( u ) du. Section 3.4 present s the r esults obtained with the Wiener-Haar wa v elet ex- pansion in Eq. 10 for hybrid dynamical systems. 3.3 Mo deling state resets A signifi cant c hallenge that hybrid dynamical sys tems present is the p ossi- bilit y of state resets [14]. When a h yb r id sy s tem switc hes from one mod e to another, the state of the system can, in general, b e r eset d iscontin uously . F or exa m p le, in th e case o f the b oun cing ball [14], the ve lo city of the ball c hanges discon tinuously after every im p act. Wh en the hybrid system tran- sitions f r om mo d e q to q ′ the state resets are typically repr esen ted as, x + = h ( x − ) , (12) where x − and x + are the states of the system b efore and after the r eset. Such discon tinuities cann ot b e easily accommod ated within the hybrid p olynomial c haos fr amew ork as d escrib ed in the previous sections. T o circum ven t this problem, one can construct a b oundary la ye r in the vicinit y of the guard condition. W e also in tro d uce a du m m y v ector ( z ) that trac ks the state x outside the b ound ary la yer and is set to x − within the b ound ary la yer. Note that we assu me sep ar ability of the states, i.e. the guard conditions (which d etermine the switc hing b et w een mo des of op era- tion) can b e w ritten indep endently of the state reset conditions in Eqn. 12. 9 In other words, the states th at determine the gu ard conditions do not par- ticipate in th e state reset. L et the reset condition b e in term s of ve ctor x in Eqn. 12, and the guard co n ditions b e in terms of v ector y (giv en by { y : g ( y ) = 0 } ). Note that, [ x y ] T represent s the entire state v ector with the f ollo wing go verning equ ation, ˙ x = f 1 ( x, y ) ˙ y = f 2 ( x, y ) . (13) W e no w construct a b ound ary la y er around th e guard condition for v ector y as follo w s :   ˙ x ˙ y ˙ z   =                             f 1 ( x, y ) f 2 ( x, y ) ( x − z ) /ǫ    if : | g ( y ) | ≥ ǫ    [ h ( z ) − x ] /ǫ ǫf 2 ( x, y ) 0    otherwise. (14) The ab o ve dynamical system is constructed such that x − ev olv es to h ( x − ) in ∆ t ≈ ǫ , where ǫ is a sm all parameter. By replacing eac h r eset condition with an equation of the form giv en by Eqn. 14, one obtains a n ew dyn amical system without resets that app ro xi- mates the original. On this new dynamical system one can use the expan- sion from Sec. 3.1 and evolv e it using Eqn . 9. In other w ord s, the fr amew ork generalizing p olynomial c h aos to hybrid systems can b e augmented using Eqn. 14 to include state resets. T o illustrate the pro cedu r e pr esen ted ab o v e we turn to the classic b ounc- ing ball example [15]: w e consider the dynamics of a ball b ou n cing on a flo or with co efficien t of restitution γ < 1 u n der the action of gra vity of uncertain magnitude ( µ ( g ) = 9 . 8 m/s 2 and σ ( g ) = 0 . 2 m/s 2 ). Thus, ev ery time the ball m ak es conta ct with the fl o or the velocit y v − is reset to a new v alue giv en by v + = − γ v − . The guard condition for r esetting the v elo cit y is giv en b y y ( t ) = 0 (where y ( t ) is the heigh t of th e ball ab o ve the flo or at time t ). The equations for the b ouncing b all are giv en by , ˙ y = v , ˙ v = − g , (15) 10 0 2 4 6 −10 0 10 20 30 Time (t) Height (y) Mean Trajectory (a) 0 2 4 6 −10 0 10 20 30 Time (t) Height (y) Boundary Layer Mean Trajectory (b) Figure 3: a) Mont e Carlo sim u lation of the b oun cing ball system with uncer- tain gra vitational acceleration. b ) Nominal b ouncing ball tr a jectory com- pared w ith the m ean tra jectory obtained through the wa v elet-based h yb rid UQ approac h with the b ound ary la yer ap p ro ximation ( ǫ = 0 . 01). with the reset cond ition at y = 0: v + = − γ v − . Thus, if one u ses the b ound ary la ye r approxima tion in Eqn. 14 we get,   ˙ y ˙ v ˙ z   =                             v − g ( v − z ) /ǫ    if : | y | ≥ ǫ    ǫv ( γ z − v ) /ǫ 0    otherwise. (16) W e now use t he hybrid p olynomial c haos exp ansion in Eqn. 9 along with the Wiener-Haar wa v elet basis fun ctions. Th e Mon te Carlo simulatio n s on the b ouncing b all are shown in Fig. 3(a). The av erage or nominal tra jec- tory is also sh own. In Fig. 3(b) we compare the nominal tra j ectory (mean tra jectory fr om Mon te Carlo) w ith the mean pred icted usin g the b oun dary la y er expansion with ǫ = 0 . 01. As sho wn, th e b oundary la yer accurately appro ximates the mean o v er multiple state r esets ev ents (in this case, eac h impact with the flo or). 3.4 Results T o demonstrate the hybrid p olynomial c haos approac h on hybrid dy n ami- cal systems w e c onsider th e simple y et chal lenging example of a switc hing 11 oscillato r give n b y Eqn . 11. ¨ x + c ˙ x + x + λ = 0 if x ≥ 0 ¨ x + c ˙ x + x − λ = 0 if x < 0 . The v alue of c is determin istic and equ al to 0.5. Here w e consider three cases with λ normally distribu ted with: µ ( λ ) = − 10 and σ ( λ ) = 2 (case 1), µ ( λ ) = 10 a nd σ ( λ ) = 2 (case 2), an d µ ( λ ) = 0 and σ ( λ ) = 1 (case 3). In all cases w e a ssume t h at the initial conditions are deterministic a nd g iv en b y [ x (0) , ˙ x (0)] =  10 − 2 , 1 . 0  . 3.4.1 Case 1: µ ( λ ) = − 10 , σ ( λ ) = 2 Let us start with the case when µ ( λ ) = − 10 and σ ( λ ) = 2 in Eqn. 11 . A rep- resen tativ e tra jectory for the dynamics of the system is sho w n in Fig. 4(a) . The corresp onding histogram for x (20 . 0; λ ) is sho wn in Fig. 4(b). Most imp ortant ly , one desires t o compute the mea n and v ariance o f x ( t ; λ ) as a function of time. In the system give n b y Eqn. 11, we expand x ( t ; λ ) using Eqn. 8 and p erform a Galerkin pro jection as sho wn in Eqn . 9. O ne then gets a system of equations for the co efficient s of expansion in Eqn . 8. These co efficien ts, once compu ted, can b e used to calculate the moments of the distribution of x ( t ; λ ). W e compare the results obtained from hybrid p olynomial chao s with those obtained using Mon te C arlo and Qu asi-Mon te Carlo based metho ds. In p articular, we use a W eyl sequence [27] along with in verse transf orm sampling [28 ] to generate the Qu asi-Mont e Carlo s amp les. W e find that the results (in the first tw o moments) from 5000 samp les of Mon te Carlo, 3000 samples of Quasi-Mon te Carlo and the Wiener-Haar h ybrid P C expansion with P = 3 are v isu ally in d istinguishable (see Figs. 5(a) and 5(b)). T reating 5000 Mon te Carlo samples as baseline, we fin d that the hybrid PC expansion has a maxim u m err or of 5 × 10 − 2 in the pr ediction of µ ( x ( t ; λ )). 3.4.2 Case 2: µ ( λ ) = 10 , σ ( λ ) = 2 The case of µ ( λ ) = 10 and σ ( λ ) = 2 is significan tly more c hallenging. A represent ative tra jectory of the system (for λ = 10) is sho wn in Fig. 6(a). When λ > 0 the system switc hes bac k and forth b et wee n mo d es. T he reason for this is th at when the system is in the righ t half-plane, th e equilibriu m of the system is in the left half-plane and vice versa. A h istogram for x (3 . 0; λ ) is d epicted in Fig. 6(b). 12 −15 −10 −5 0 5 10 15 −5 0 5 x dx/dt (a) 0 5 10 15 20 0 500 1000 1500 2000 2500 3000 x Number of occurrences (b) Figure 4: a) A r epresen tativ e tra j ectory for case 1. b) Histogram of x (20; λ ) for case 1. 0 5 10 15 20 0 5 10 15 Time Mean of x MC 5000 samples QMC 3000 samples Hybrid PC (P=3) (a) 0 5 10 15 20 0 2 4 6 8 10 Time Variance of x MC 5000 samples QMC 3000 samples Hybrid PC (P=3) (b) Figure 5: Comparison of a ) predicted mean and b) predicted v ariance of x ( t ; λ ) for v arious UQ metho ds for case 1. 13 −0.1 −0.05 0 0.05 0.1 −2 −1 0 1 2 x dx/dt (a) −0.15 −0.1 −0.05 0 0.05 0.1 0 500 1000 1500 2000 2500 3000 x Number of occurrences (b) Figure 6: a) A represen tativ e tra jectory for case 2. b) Histogram of x (3 . 0; λ ) for case 2. W e again compare h ybrid Wiener-Haar p olynomial chaos to Monte Carlo sampling in Fig. 7. W e find th at hybrid Wiener-Haar p olynomial c h aos ( p = 5) accurately computes the mean and the v ariance of the distribu tion of x . The maxim u m absolute error of h ybrid p olynomial c h aos in mean is µ ( x ( t ; λ ) = 1 . 8 × 10 − 3 and v a r iance is 7 × 10 − 5 . Note that expansions in terms of s tandard b asis fun ctions suc h as Hermite and Legendre p olynomi- als are unable to compute the momen ts of x ( t ; λ ) b ey ond a threshold time that d ep ends weakly on the order of expans ion p (see Fig. 8). The solution in this case is particularly c h allenging b ecause it b ecomes more oscillatory in terms of λ at t increases ( Fig. 9). The Wiener-Haar basis fu n ctions are naturally oscillatory and hence m ore accurate than Hermite p olynomials in captur in g the so lu tion x ( t ; λ ) . Note that, for large time sim ulations the Wiener-Haar expansions will also fail sin ce the solution will ev entually b e- come to o oscillatory for the order of expansion. This p roblem is w ell known in the p olynomial chao s literature [29]. 3.4.3 Case 3: µ ( λ ) = 0 , σ ( λ ) = 1 . 0 W e also consider the case of µ ( λ ) = 0 w ith σ ( λ ) = 1 . 0. T h is case is partic- ularly c hallenging because there is a concentrati on o f probability of x ( t ; λ ) as s ho wn in Fig. 10( b). The reason for this is as follo ws: when µ ( λ ) = 0, the nominal tra jectory con verge s to 0 and so do all tra jectories with λ > 0. Indeed, f or tra jectories with λ > 0 the equilibriu m lies in the opp osite half- plane with resp ect to the cur ren t s tate. This giv es rise to deca ying sw itc hing tra jectories, as case 2 in Fig. 6(a). Note that for λ < 0, the tra jectories are 14 0 1 2 3 −0.04 −0.02 0 0.02 0.04 0.06 Time Mean of x MC 5000 Samples QMC 3000 Samples Hybrid PC (P=5) (a) 0 1 2 3 0 0.5 1 1.5 2 x 10 −3 Time Variance of x MC 5000 Samples QMC 3000 Samples Hybrid PC (P=5) (b) Figure 7: Comp arison of a) p redicted mean and b ) pr edicted v ariance of x ( t ; λ ) b y v arious UQ metho ds for case 2. 0 1 2 3 −0.04 −0.02 0 0.02 0.04 0.06 Time Mean x MC 5000 samples QMC 3000 samples PCM 9 points (a) 0 0.5 1 1.5 2 2.5 3 0 0.5 1 1.5 2 x 10 −3 Time Var x MC 5000 samples QMC 3000 samples PCM 9 points (b) Figure 8: a) Mean and b) V ariance predicted by standard (Hermite) p oly- nomial chao s basis functions for case 2. 15 −0.2 −0.1 0 0.1 0.2 x(t; λ ) t=0 −0.2 −0.1 0 0.1 0.2 x(t; λ ) t=0.5758 −0.2 −0.1 0 0.1 0.2 x(t; λ ) t=1.1818 5 10 15 0 0.5 1 ρ ( λ )/ ρ ( µ ) λ µ =10, σ =2 −0.2 −0.1 0 0.1 0.2 x(t; λ ) t=0 −0.2 −0.1 0 0.1 0.2 x(t; λ ) t=0.5758 −0.2 −0.1 0 0.1 0.2 x(t; λ ) t=1.1818 5 10 15 0 0.5 1 ρ ( λ )/ ρ ( µ ) λ µ =10, σ =2 Figure 9: C ase 2: x ( t ; λ ) b ecomes more oscillatory in λ as t increases. The b ottom plot sh o ws the distribution for λ . 16 −1 −0.5 0 0.5 1 −1 −0.5 0 0.5 1 x dx/dt (a) −0.015 −0.01 −0.005 0 0.005 0.01 0 1000 2000 3000 4000 x Number of Occurrences (b) Figure 10: a) Nominal tra jectory for case 3. b ) Hi s togram of x (2 0; λ ) for case 3. similar to the ones in case 1 (Fig. 4(a)). The Wiener-Haar basis fun ctions along with the hybrid p olynomial c h aos approac h accurately capture the moment s of the distribution for x ( t ; λ ) (see Fig. 11). In f act, an expansion to just P = 3 captures the first t wo m omen ts. The step-fun ction nature of the Wiener-Haar b asis allo ws it to p erform w ell in this scenario. Standard basis f unctions lik e Hermite p olynomials are com- pletely incapable of accurately capturing the momen ts of th e distrib u tion for x ( t ; λ ) sho wn in Fig. 6(b). 4 T ra n sp ort theory approac h for u ncertain t y q u an- tification in h ybrid systems In this section w e presen t a qualitativ ely d ifferen t approac h to UQ in hy- brid systems b ased on transp ort equations. W e wr ite an adv ection equation for the pr obabilit y densit y of the state and expand this equation in an ap- propriate basis, as is done in p olynomial c haos. The resulting equation is equiv alen t t o the F okk er-Planc k equation [30] in the abs en ce of a diffusion term. T hough significant effort has b een pu t into computing solutio n s for the F okk er-Planc k equation in v arious applications [30, 31], our setting is particularly c hallenging due to t h e switc hing dynamics of hybrid systems. W e n ote that advec tion equations f or probabilit y distribution f u nctions ha ve b een used to pr opagate uncertain ty th r ough heterogeneous p orous media with un certain p r op erties [32] an d for hyp erb olic conserv ation la ws with noise [33]. Recen tly , simila r metho d s ha ve b een extended to cumulati ve 17 0 5 10 15 20 0 0.2 0.4 0.6 0.8 Time Mean of x MC 5000 samples Hybrid PC (P=3) (a) 0 5 10 15 20 0 0.2 0.4 0.6 0.8 1 Time Variance of x MC 5000 samples Hybrid PC (P=3) (b) Figure 11: Comp arison of a) pred icted mean and b) predicted v ariance of x ( t ; λ ) for v arious UQ metho ds for case 3. distribution functions in hyp erb olic conserv ation laws [34]. The p olynomial c haos expansion in this setting y ields a system of h yp er- b olic partial differen tial equations for the co efficien ts of the expansion, wh ic h are th en solve d by in tegrating along c h aracteristics. The hybrid nature of the original system is reflected in that th e c h aracteristics exhibit switching. Ev en though we only consider systems without r esets, we can use the results of S ec. 3.3 to treat systems with r esets. As in Sec. 3, let us consider a h yb rid system with ou t resets and uncertain parameters λ with guard cond itions ind ep endent of λ : 2 ˙ x = f i ( x, λ ) w h en G i ( x ) is true. The system has u ncertain initial conditions describ ed by the p robabilit y densit y ρ x 0 ( x ) and the uncertain parameters λ follo w ρ λ ( λ ). W e describ e the system b y the time evo lu tion of the distribution f unction ρ ( x, λ ; t ), whic h has initial condition ρ ( x, λ ; 0) = ρ x 0 ( x ) ρ λ ( λ ) (initial uncertaintie s are in dep end ent) and n ormalizatio n Z ρ ( x, λ ; t ) dxdλ = 1 . Note that, for all time, we ha ve ρ λ ( λ ) = Z ρ ( x, λ ; t ) dx. (17) 2 Note that t h is embo dies the constraint of having no ov erlap in the d omains for different mod es. 18 Our goal is to compute the evol u tion of the den s it y in x (the marginal distribution): ρ x ( x, t ) = Z ρ ( x, λ ; t ) dλ. Ho we ver, without introd ucing assump tions on ρ λ , the equation for ρ x is not closed. W e therefore fo cus on compu ting the ev olution of ρ directly through an expansion. F rom this e v olution, ρ x can then b e calculated at ev ery instant. 4.1 Equation for ρ Let u s defin e the sets S i = { x : G i ( x ) is tru e } and th e ind icator fun ctions 1 i ( x ) = ( 1 if x ∈ S i 0 if not . With this notation, and b ecause λ is constan t along a tra jectory , we ha ve ∂ ρ ∂ t + ∇ · ( ρf ) = 0 ∀ λ ( 18) where f ( x, λ ) = P i 1 i ( x ) f i ( x, λ ) and th e gradien t o p erator acts on ly on x and n ot on λ . 4.2 Boundary conditions at the in terfaces The discon tinuit y in the equation implies that mass ma y accum ulate at the b oundaries b et we en zones where different guard conditions are v alid. In tegrating on a cylinder that crosses one su c h b oun dary we obtain the matc hing condition ∂ σ ∂ t + ∇ s · ( σ f ) = ρ i f i · ˆ n ik − ρ k f k · ˆ n ik , (19) where σ is a surface p robabilit y densit y b et wee n r egions S i and S k , ˆ n ik is the sur f ace norm al fr om S i to S k , ∇ s is the diverge n ce in the space tangent to the surface, and f is th e flo w at the surfac e . 3 This ma y lead to a cascade, with probabilit y condensin g into progressively lo wer dimensional structures: where tw o hypers urfaces meet (the b oundary b et we en three guard condi- tions) th e s ame scenario rep eats, u n til we h a ve mass accum ulating at p oint s . 3 Whether f is f i or f k on the surface will dep end on how the guard conditions are expressed. 19 T o solv e Eqn . 18 the initial condition m ust include initial v alues for σ and the pr obabilit y density on any lo we r dimensional structure wher e mass may accum ulate. The disconti n u it y of f do es not ne c essarily imply accumulat ion. In f act, at an interface we hav e sev eral options: 1. Both f i · ˆ n ik and f k · ˆ n ik are nonzero and hav e the same sign. In th is case, there is n o accum ulation. If we assume that ρ h as a singularit y at the in terface, the flow will mov e the singularit y a wa y from it. 2. f i · ˆ n ik > 0 and f k · ˆ n ik ≤ 0: accumulatio n o ccurs. 3. f i · ˆ n ik ≥ 0 and f k · ˆ n ik < 0: accumulatio n o ccurs. 4. f i · ˆ n ik = f k · ˆ n ik = 0: n o accum u lation. 5. f i · ˆ n ik ≤ 0 and f k · ˆ n ik ≥ 0: n o accumulatio n . If we are in the case without accum ulation and w ithout initial concen tration of dens it y in low er dimensional stru ctures, then σ = 0 and Eqn . 19 b ecomes ρ i f i · ˆ n ik = ρ k f k · ˆ n ik . (20) Theorem 1. Any se c ond or der O D E of the form ¨ X = F i ( X, ˙ X , λ ) when G i ( X ) i s true satisfies the c onditions for no ac cumulation at the interfac e wher e the ODE is disc ontinuous. Pr o o f . Without loss of generalit y we can fo cus on just t w o regions S i and S k and r ewr ite the problem as the first-order ODE ˙ x ≡  ˙ X ˙ Y  = f ( x ) =  Y F i ( X, Y , λ )  when G i ( X ) is tru e. T o find the n ormal ˆ n ik w e consider a C 1 function b ( x ) = b ( X , Y ) = B ( X ) that is positiv e in S k and negativ e S i so that th e in terface is giv en by the lo cus of b ( x ) = 0. The gradient of this fun ction is prop ortional to the normal: ˆ n ik ∝ ∇ b =  ∇ X B 0  and therefore ˆ n ik · f = Y · ∇ X B ||∇ X B || , whic h is con tinuous at the inte r face. 20 4.3 Expansion of the equa t ion for ρ A t every p oint x we expand th e distribu tion in λ : ρ ( x, λ ; t ) = X k a k ( x, t ) w ( λ ) ψ k ( λ ) , (21) where { ψ k } forms an orthogonal basis with r esp ect to w : Z ψ i ( λ ) ψ k ( λ ) w ( λ ) dλ = w k δ ik . W e k eep w k to allo w for a non -n ormalized w eight fun ction w ( λ ). W e replace the expansion in Eqn . 21 in to Eqn . 18 and pro ject onto ψ i to obtain a set of p artial differentia l equations f or the co efficien ts a i ( x, t ): ∂ a i ( x, t ) ∂ t + 1 w i ∇ · X k a k ( x, t ) Z w ( λ ) ψ i ( λ ) ψ k ( λ ) f ( x, λ ) dλ = 0 . (22) Since this equation is lo cal in x there is no question as to whic h f ( x, λ ) m u st b e u sed at an y giv en p oint. 4.4 Example: switc hing oscillator Here we revisit the switc hing oscillator system ¨ x = ( − x − γ ˙ x − λ if x ≥ 0 − x − γ ˙ x + λ otherwise, whic h can b e expressed as the 2–D system  ˙ x ˙ y  =  f x f y  where f x = y and f y = ( − x − γ y − λ if x ≥ 0 − x − γ y + λ otherwise . The transition p oin ts are lo cated at x = 0 and therefore f · ˆ n = f ·  1 0  = f x = y whic h is con tinuous and therefore has the same sign on b oth sides. Therefore there is no mass accumulatio n at the interface for this system. This is a sp ecial case of theorem 1. 21 4.4.1 Case 3 revisited T o connect with the example presented in case 3 ( µ = 0, σ = 1) w e choose w ( λ ) = e − λ 2 / 2 and ψ k ’s as th e pr obabilist’s Hermite polynomials H k , with the f ollo wing prop erties Z H i ( λ ) H k ( λ ) w ( λ ) dλ = k ! √ 2 π δ ik H k +1 ( λ ) = λH k ( λ ) − H ′ k ( λ ) H ′ k ( λ ) = k H k − 1 ( λ ) . Calculating th e terms in Eqn. 22: Z wH i H k f x dλ = y i ! √ 2 π δ ik Z wH i H k f y dλ = − ( x + γ y ) i ! √ 2 π δ ik ∓ i ! √ 2 π ( δ i,k +1 + k δ i,k − 1 ) , where the upp er sign ( − ) is for x ≥ 0 and the lo we r sign (+) is for x < 0. Substituting in to Eqn. 22 we obtain the equations 0 = ∂ t a 0 + ∂ x ( y a 0 ) + ∂ y [ − ( x + γ y ) a 0 ∓ a 1 ] 0 = ∂ t a i + ∂ x ( y a i ) + ∂ y [ − ( x + γ y ) a i ∓ ( i + 1) a i +1 ∓ a i − 1 ] ( i ≥ 1) . (23) Note that, instead of us ing th e Hermite p olynomials, one can u s e the Haar w av elet expansion [22] to r epresen t the solution as wa s done in p revious sections. W e p lan to present this calculation in f uture work. Theorem 2. The system of PDEs for the switching oscil lator is hyp erb olic. Pr o o f . C onsider a system of P DEs of th e form ∂ t u + X ν A ν ∂ ν u = B , where u ( x 1 , . . . , x n , t ) ∈ R m and the A ν are m × m matrices. The system is h yp erb olic if f or any α ν ∈ R the lin ear com bination A = P ν α ν A ν has real eigen v alues. 22 F or the s witc hing oscillator the s y s tem of PDEs can b e written as ∂ t    a 0 a 1 . . .    +    y 0 . . . . . 0 y . . . . . .    ∂ x    a 0 a 1 . . .    +      β ∓ 1 0 . . . . . . . ∓ 1 β ∓ 2 0 ∓ 1 β ∓ 3 . . . . . . . . . . . .      ∂ y    a 0 a 1 . . .    = γ    a 0 a 1 . . .    , where β = − ( x + γ y ). Th u s , any com b ination of th e matrices A ν is going to b e of the trid iagonal form A =         a b 0 . . . 0 b a 2 b 0 . . . 0 b a 3 b . . . 0 b a 0 . . . . . .         . This tridiagonal n on -sym metric m atrix is similar to a tridiagonal symmetric matrix w ith a diagonal similarity matrix: S = D AD − 1 , where D =         √ 0! √ 1! 0 √ 2! 0 √ 3! . . .         . A is similar to S , a symmetric and real matrix, and therefore A has real eigen v alues. The issu e o f h yp erb olicit y is d iscussed in more depth i n [35]. T o solv e the hyp erb olic system from Eqn. 23 we write it in the f orm ∂ t a + y ∂ x a − ( x + γ y ) ∂ y a ∓ A∂ y a = γ a, 23 where A is the tridiagonal m atrix A =         0 1 1 0 2 0 1 0 3 0 1 0 . . .         . (24) W e diagonalize A = P Λ P − 1 and define b = P − 1 a to obtain the set of uncoupled h y p erb olic PDEs ∂ t b i + y ∂ x b i + [ − ( x + γ y ) ∓ λ i ] ∂ y b i = γ b i , (25) where λ i is the i -th eigen v alue. W e will no w prov e that wh en the exp ansion is tru ncated up to a n − 1 , i.e., A is tru ncated to a n × n matrix, the eigen v alues A are the zeros of H n . Theorem 3. Th e eigenvalues of A n , the n × n trunc ate d version of the matrix in E q n. 24, ar e th e zer os of the n -th or der pr ob abilist’s Hermite p olyno mial H n . Pr o o f . W e pro ceed by indu ction to pro ve that d et( A n − λI ) = H n ( − λ ), whic h will then, b y the symmetry of H n , p ro ve our r esu lt. Let B n = A n − λI . I ndeed, d et B 1 = − λ and d et B 2 = λ 2 − 1. F or the general case, B n +1 =        0 B n . . . 0 n 0 · · · 0 1 − λ        Therefore, det B n +1 = − λ det B n + ( − 1) n n ( − 1) n − 1 det B n − 1 = − λ det B n − n det B n − 1 , whic h is the r ecurrence relation satisfied by H n ( − λ ). The charact eristic curve s of Eqn. 25 are giv en b y ˙ x = y ˙ y = − x − γ y ∓ λ i ˙ b i = γ b i . 24 In other wo r ds, the c haracteristics are damp ed oscillators wh ere th e equilib- rium p osition is giv en by the eigenv alues of ∓ A . Th e exp onential growth of b i along a tra jectory is due to the con traction in phase space pro du ced by the d issipation γ . 4.4.2 Results W e no w sh ow results obtained by using transp ort theory approac h on case 3 ( µ = 0, σ = 1) for the switc hing oscillator (Eqn. 11). In Fig. 12, we sh o w a series of pr obabilit y d istribution snapsh ots for Mon te Carlo (5000 samples) and the transp ort op erator metho d (for case 3), gridded in the ( x, y ) plane. Note that we set y = ˙ x , as d efined in the first part of the paper. Th e Mon te Carlo color ma p snapsh ots show that the distri- bution lies on a one-dimensional manifold (in t wo d imensional s p ace). The one-dimensional nature of the distrib ution arises b ecause we chose d eter- ministic initial conditions. As discuss ed pr eviously , all tra j ectories in case 3 with λ ≥ 0 con verge to the origin, resulting in a jump in th e cumulativ e distribution function (CDF). F or λ < 0, how ev er, the tr a jectories con verge asymptotically to either + λ or − λ . This con v ergence to + λ or − λ is highly d ep endent on the individu al tra jectory and results in a fragmen tation of the outpu t distribu tion. Close to con verge n ce ( t = 18), v ery few tra jectories con verge to a p oin t in the range 0 . 2 < x < 0 . 4, as seen in the fl at region of the CDF in Fig. 13. The transp ort-based metho d captures th e sin gularit y at the origin accurately , but is unable to accurately capture the fragmen tation. This is b ecause the method samples the distrib ution sp ars ely (determined by the order of expansion), r esulting in UQ acceleratio n . Ho wev er, this spars ity makes the metho d miss su ch fine details. On using a h igh order of expansion ( n = 70), some samples partially capture th e structure aroun d x = 0. Note that a m uch low er order expans ion accurately captur es th e j ump at the origin and the asymp totic ( x > 1) shap e of th e CDF. As in Fig. 11, w e compare Monte Carlo (5000 s amp les) with the tran s - p ort th eory app roac h (orders 15 and 70 expansion) in Fig. 14. The L ∞ (maxim um) err ors for n = 15 are 9 . 56 × 10 − 2 (mean) and 7 . 95 × 10 − 2 (v ari- ance). F or n = 70 w e get L ∞ errors of 5 . 28 × 10 − 2 and 4 . 92 × 10 − 2 for mean and v ariance resp ectiv ely . As sho wn in Fig. 11, the metho d p erforms rea- sonably well , ho wev er, the r esults are not nearly as go o d as those obtained using h yb rid p olynomial c haos with the Wiener-Haar wa v elet expansion in section 3.1. Ho w ever, with a b etter c h oice of basis fun ctions, one do es exp ect b etter results. The transp ort op erator theory is attractiv e a s it app ears to 25 Transport theory approach n = 70 x −1 0 1 2 3 y Monte Carlo 5,000 samples −1 0 1 y −1 0 1 y −1 0 1 y −1 0 1 x y −1 0 1 2 3 −1 0 1 0 0.5 1 CDF 0 0.5 1 0 0.5 1 0 0.5 1 −1 0 1 2 3 0 0.5 1 x Monte Carlo Transport theory t = 10 t = 8 t = 6 t = 4 t = 2 Figure 12: Snapshots at t = 2 , 4 , 6 , 8 , 10 of th e gridding from 5000 sam- ples of Mon te Carlo and an order 70 expansion using the transp ort theory based metho d. Th e right column compares the one-dimen sional cumulativ e distribution function f or the corresp ond ing times. 26 −1 0 1 2 3 0 0.2 0.4 0.6 0.8 1 x CDF Monte Carlo Transport theory Figure 13: Comp arison at t = 18 of the cumulativ e d istribution fun ction from 5000 samples of Mon te C arlo (solid) and an order 70 expansion using the tr an s p ort theory based metho d (dashed). b e more v ersatile. I n general, the transp ort op erator app roac h is applicable to h yb rid dyn amical systems w ith o verlapping mo des of op eration (by con- structing m u ltiple PDEs f or the o ve r lapping mo d e). In cont r ast, the h ybrid p olynomial chaos m etho d suffers fr om the disadv an tage of b eing inapplicable to such systems. 5 Conclusions As the mo deling of h ybr id dynamical systems b ecomes increasingly imp or- tan t for m o dern day engineering applications such as electrical and biologi- cal net works, air traffic systems, communicati on netw orks, etc., quant ify in g uncertain ty in these systems is going to b ecome a cen tral concern. S ince uncertain ty q u an tification allo ws one to compu te momen ts of output d istri- butions in the presence of parametric uncertain ty , these tec hniqu es will b e used to aid decisions r elated to robu st system d esign and p erf orm ance. In this work, we h av e made th e first attempts to dev elop f ast uncertain ty quan tification metho d s targeted for hybrid d ynamical systems. In particu- lar, w e extended p olynomial c h aos metho ds, a p opular tec h nique for prop- agating uncertain ty through smo oth sy s tems, to h yb rid dynamical systems. W e also deve lop ed metho ds to hand le state resets within the p olynomial 27 0 5 10 15 20 0 0.2 0.4 0.6 0.8 Time Mean of x MC 5000 samples n = 15 n = 70 (a) 0 5 10 15 20 0 0.2 0.4 0.6 0.8 1 Time Variance of x MC 5000 samples n = 15 n = 70 (b) Figure 14: Comp arison of a) pred icted mean and b) predicted v ariance of x ( t ; λ ) for case 3. T he curv es show a 5000-sample Monte Carlo ru n, and 15 & 70 term expansions u sing the transp ort theory approac h. c haos framewo r k by u sing b oun dary lay er appro ximations. W e then applied this new approac h to p erform un certain t y quan tification on switc h ing har- monic oscillators and the b oun cing ball examples. W e also demonstrated the efficacy of u s ing Wiener-Haar expan s ions [9, 22, 23] with our hybrid p olynomial chao s approac h for qu an tifying uncertain t y in hybrid systems that giv e rise to multi-mod al distribu tions or b ecome increasingly oscilla- tory in time. Finally , we sho we d ho w a transp ort theory based app r oac h can capture naturally-e m erging discon tin uities in the distribution. F u tu re efforts in volv e pro vidin g r igorous error b ound s for Wiener-Haar expansions in the hybrid p olynomial c haos setting with b ound ary la yer expansions . W e are also extending our hybrid p olynomial c haos approac h to net wo r ks of hy- brid dynamical systems using our r ecen t w ork on propagating uncertain t y through complex netw orks [36 ]. W e also intend to extend the transp ort op erator based UQ metho d to hybrid systems with o verlapping mo d es of op eration. 6 Ac kno wledgemen ts The authors thank Habib Na jm for p oin ting u s to h is w ork on Wiener-Haar based p olynomial c haos expansion and his insigh tful i n put. W e also thank Alessandro Pin to and George Mathew for v aluable discussions related to h yb rid dyn amical systems. 28 References [1] R. E. Caflisch. Mon te Carlo and qu asi-Mon te C arlo metho ds. A cta Numeric a , 7:1–49, 1998. [2] R. H. My ers , D. C. Montg omery , and C. M. Anderson-Co ok. R esp onse Surfac e Metho dolo g y . Wiley , third ed ition, 2009. [3] T. Hastie, R . Tibsh ir ani, and J. F riedman. The Elements of Statistic al L e arn i ng . S pringer, second edition, 2009. [4] X. W an a n d G. E. Karn iadakis. Re cent adv ances in p olynomial c h aos metho ds and extensions. In Computational U nc ertainty in Military V ehicle Design M e eting Pr o c e e dings . NA TO/OT AN , Pa p er Reference Num b er: R TO -MP-IST-999, 2008. [5] N. Wiener. T he h omogeneous chaos. Americ a n Journal of M athematics , 60(4): 897–9 36, 193 8. [6] M. S . Allen and J. A. C am b eros. Comparison of uncertain ty prop- agatio n / resp onse sur face techniques for t w o aeroelastic systems,. In 50th AIAA Structur es, Structur al Dynamics, and Materials Confer enc e, Palm Springs, California, May 4-7, 2009 , 2009. [7] H. C. Elman, C . W. Miller, E. T. Phip ps, and R. S. T umin aro. Assess- men t of co llo cation and Gal erkin approac hes to lin ear d iffusion equa- tions with random data. Int. J. U nc ertainty Quantific ation , 1:19–2 3, 2011. [8] R. Ghanem. Probabilistic charact er ization of transp ort in heteroge- neous media. Comput. Me tho ds A ppl. Me ch. Engng . , 158:199– 220, 1998. [9] H. N. Na j m, B. J. Debusschere, Y. M. Marzouk, S. Widmer, and O. P . Le Ma ˆ ıtre. Uncertain ty quan tification in c hemical systems. Int. J. Numer. Me th. Engng. , 80:789 –814, 2009. [10] T. Sahai, V. F onob erov, and S. Loire. Unce r tain t y as a stabilizer of the head-tail ordered p hase in carb on-mono xide monola ye r s on graphite. Physic al R eview B , 80(11):11 5413, 2 009. [11] D. Xiu and G. E. Karniadakis. Mo deling uncertain ty in flow simulatio n s via generalized p olynomial chao s. J . Comp. Phys. , 187:137– 167, 20 03. 29 [12] N. Zabaras and B. Ganapathysubramanian. A scalable framew ork for the solution of stochastic in verse p roblems using a sparse grid colloca- tion app roac h. J . Comp. Phys. , 227:4697– 4735, 2 008. [13] X. W an and G. E. Karniadakis. An adaptiv e m u lti-elemen t generalized p olynomial c haos metho d for sto c h astic d ifferen tial equations. Journal of Computational Physics , 209(2):6 17–64 2, 2005 . [14] C. G. Cassandras and J. Lygeros. Sto chast i c Hybrid Systems . CR C T a ylor & F rancis, fir st edition, 1991. [15] K. H. Johansson, M. Egerstedt, J. Lygeros, and S. Sastry . Regulariza- tion of Zeno hybrid automata . Sy stems and Contr o l L etters , 38(3):141– 150, 1999. [16] R. Alur, C . Belta, F. Iv anic, V. Kumar, M. Min tz, G. P app as, H. Rubin, J. S c hug, and G. J. P appas. Hybrid Systems: Computation and Contr ol . Springer-V erlag, fir st edition, 2001. [17] N. Chabrier and F. F ages. Symb olic mo del chec king of b io c hemical net works. I n Computational Metho ds in Systems Biolo gy (CMSB03), volume 2602 of LN CS , p ages 149–162 . Sp r inger-V erlag, 2003. [18] C. T omlin, G. J. P app as, and S. S astry . Conflict resolution for air traffic managemen t: A study in m u ltiagen t h yb rid s y s tems. IEEE T r ansac- tions on Automatic Contr ol , 43:509–5 21, 1998. [19] J. P . He s p anha. S to c hastic h ybrid systems: Application to communi- cation netw orks. In Hybrid Systems: Computation and Contr ol, ser. L e ct. N otes in Comput. Scienc e , pages 387–401. Sprin ger-V erlag, 2004. [20] E. Asarin, O. Maler, and A. Pn ueli. S ym b olic con troller syn thesis for discrete and timed systems. In Hybrid Systems II, LN CS 999 , p ages 1–20. S pringer, 1995. [21] A. Haar. Zur T h eorie d er orthogonalen F u nktionensysteme. Mathema- tische Annalen , 69:331–37 1, 1 910. [22] O. P . Le Ma ˆ ıtre, O. M. Knio, H. N. Na jm, and R. G. Ghanem. Un - certain t y propagation using Wiener-Haar expans ions . J . Comp. Phys. , 197:28 –57, 2004. [23] O. P . Le Ma ˆ ıtre, H. N. Na jm, R. G. Ghanem, and O. M. Kn io. Multi- resolution analysis of Wiener-t yp e u n certain t y propagation sc h emes. J. Comp. P hys. , 197:502–5 31, 2004. 30 [24] H. Niederreiter. Quasi-Mon te Carlo metho ds and pseud o-rand om num- b ers. Bul letin of the A meric an Math ematic al So ciety , 84(6):957 –1041 , 1978. [25] X. W an and G. E. Karniadakis. Be yond Wiener-Ask ey expansions: Handling arbitrary PDFs. Journal of Scientific Computing , 27:455– 464, 2006. [26] H. Ogura. Orthogonal functions of the P oisson pro cesses. IEEE T r ans- actions on Information The ory , 18(4):473 –481, 1 972. [27] H. Niederreiter. R andom Numb er Gener ation and Quasi-Monte Carlo Metho ds . SIAM, 1992. [28] L. De vr o ye . Non-Uniform R andom V ariate Gener ation . Springer- V erlag, firs t edition, 1986. [29] X. W an and G. E. Karn iadakis. Long-term b eha v ior of p olynomial chao s in sto chastic flow sim u lations. Computer Metho ds in Applie d Me chanics and E ngine ering , 195:558 2–5596, 2006. [30] H. Risk en. The F okker- Planck E quation: M e tho ds of Solution and A p- plic ations . Spr inger, second ed ition, 1989. [31] Crispin W Gardiner. H andb o ok of sto c hastic metho ds: for physics, chemistry and the natur al scienc es; 3r d e d. S p ringer series in syn er- getics. Sp ringer, Berlin, 2004. [32] Daniel M. T artako vsky and Svetl ana Bro yda. PDF equations for adv ectiv e-reactiv e transp ort in h eterogeneo u s p orous m edia w ith un- certain pr op erties. J. Contam. Hydr ol. , 120-121 :129–140, 2011. [33] W uan L uo. Wiener Chaos Exp ansion and Numeric al Solutions of Sto chastic Partial D iffer ential Equations . PhD thesis, California In- stitute of T echnolog y , 2006 . [34] P eng W ang and Daniel M. T artak o v s ky . Un certain t y q u an tification in kinematic-w a ve mo dels. J. Comput. Phys. , 231:7868– 7880, 2012. [35] J. T r y o en, O. Le Ma ˆ ıtre, M. Ndjinga, and A. E r n. Intrusiv e pro jection metho ds with up win ding for un certain nonlinear hyp erb olic systems. J. Comp. P hys. , 229:6485– 6511, 2010. 31 [36] Amit Sur ana, T uh in Sahai, and Andrzej Banaszuk. Iterati ve metho ds for scalable uncertaint y qu an tification in complex n et wo r ks. Int. J. Unc ertainty Quantific ation , 2:413–4 39, 2012. 32