Hybrid models for future Event prediction

hybrid models in artificial intelligence and hybrid models for future event prediction and hybrid choice models progress and challenges
Dr.MohitBansal Profile Pic
Published Date:26-10-2017
Your Website URL(Optional)
HYBRID SYSTEMS Hybrid systems comprise continuous and discrete-event subsystems that interact. The continuous parts, unlike the discrete-event parts, cannot be simulated in their origi- nal form. Instead, continuous models are approximated with discrete-event models amenable to computer simulation. Any number of approximating systems can be used, each answering in a different way three fundamental questions: 1. In any finite interval of time, the continuous system traverses an infinity of states, but estimates can be computed for only a finite number of them (i.e., the discrete system must be legitimate); which points are picked? 2. To calculate these points requires knowledge of the trajectory between them; what is assumed? 3. The system interacts with its environment; how are inputs and outputs, contin- uous and discrete, handled? For the moment, consider only questions 1 and 2 in relation to a single, ordinary differential equation x ˙ = f (x) This continuous model can be approximated by (1) discretizing time into points separated by intervals of duration h and (2) assuming that x follows a line in thoseHYBRID SYSTEMS 183 intervals. This yields the discrete-time system x(t + (n+ 1)h)= x(t + nh)+ hf (x(t + nh)) 0 0 0 which estimates x at times t + kh for k = 0, 1, 2,... . Indeed, this is Euler’s method 0 for the numerical solution of an ordinary differential equation. Another approach is to discretize x into points separated by intervals of length q, thus yielding a set of discrete states x + kq for k =...,−2,−1, 0, 1, 2,....The 0 goal now is to determine when x takes on these discrete values. By again assuming that x is a line between calculated points, the discrete-event system  q/ f (x + kq) if f (x + kq)= 0 0 0 ta(k)= ∞ otherwise δ (k)= k+ f (x + kq)/ f (x + kq) int 0 0 λ(k)= x +δ (k)q 0 int generates an output trajectory which approximates the continuous x. This is a first- order quantized state system (see, e.g., Refs. 70, 100, and 158). To be concrete, consider a rock dropped from a bridge. Its downward velocity v is described by the ordinary differential equation D v ˙ = g− v (5.1) m where m is the rock’s mass, g is acceleration due to gravity, and D is the coefficient of aerodynamic drag. It is not necessary to simulate this system; the rock’s speed can be written down once and for all if the rock’s initial velocity v is known. With this 0 information, we have      g D D v(t)= 1− exp − t + v exp − t (5.2) 0 D/m m m and that this is the solution to Equation 5.1 can be verified by taking its derivative with respect to time. Applying Euler’s method, we select an interval of duration h and approximate this continuous model with the discrete system   D v = v + h g− v k+1 k k m184 HYBRID SYSTEMS 2 Exact QSS 1.8 Euler 1.6 1.4 1.2 1 0.8 0.6 0.4 0.2 0 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 time (s) FIGURE 5.1 The Euler method and quantized state method applied to the model of a falling 2 rock with g= 9.8m/s , m= 1 kg, D= 5, and v = 0 m/s. The discrete models use h= 0.1 0 and q = 0.2 respectively. where v is the speed of the rock at time kh. A quantized state approximation using k q yields the discrete-event model  q/g− (D/m)(v + kq) if g− (D/m)(v + kq)= 0 0 0 ta(k)= ∞ otherwise δ (k)= k+ (g− (D/m)(v + kq))/g− (D/m)(v + kq) int 0 0 λ(k)= v +δ (k)q 0 int Figure 5.1 compares the exact trajectory and the event trajectories produced by these approximations. Euler’s method and its kin, which discretize time, are the mainstays of numerical analysis and will be the focus of this chapter. Quantized state methods are nonetheless attractive for many applications (see, e.g., Refs. 70, 93, 94, 103, and 138) and can be incorporated with ease into the framework developed here; this intriguing topic is discussed in several textbooks (see, e.g., those by Zeigler et al. 157 and Cellier and 1 Kofman 20). 1 Notation in this chapter differs from the previous parts of the book and may be confusing. The set Q stands for a set of discrete states rather than set of total states; q is a member of Q or, in the code, a continuous state variable. So, too, x is a continuous state variable or, in the code, a discrete input. Usage of any symbol is consistent within a section, but beware of changes between sections. This unfortunate situation is a compromise between conflicting uses of x, q,and Q in the literature on discrete-event and hybrid systems. velocity (m/s)AN ELEMENTARY HYBRID SYSTEM 185 5.1 AN ELEMENTARY HYBRID SYSTEM Consider a system with a state that evolves continuously but an input that is discrete; the lopsided square wave that drives the tank’s motors is one example. The input signal u takes values u , u ,..., u at times t , t ,..., t . The equation that describes 0 1 n 0 1 n this system is ˙ x = f (x, u) (5.3) Begun in state x = x(t ) and fed the input trajectory ut , t), with t t , the state 0 0 0 n trajectory of the system is the solution to Equation 5.3: t j+1 t   n−1  x(t)= x + f (x(τ), u ) dτ + f (x(τ), u ) dτ (5.4) 0 j n j=0 t t j n Euler’s method is adapted to approximate this system in the following way. A fixed stepsize h is used to compute x(t) through intervals in which u is constant. Changes in u are inputs to the discrete-event system. The external transition function responds to them by using the elapsed time for the integration step and then storing the new value of u. This ensures that the endpoints of the integrals are handled precisely. This method is realized by the (output-free) discrete-event model X =R S=R×R δ ((x, u))= (x+ hf (x, u), u) int ′ ′ δ ((x, u), e, u )= (x+ ef (x, u), u ) ext ′ ′ δ ((x, u), u )= (x+ hf (x, u), u ) con ta((x, u))= h which has a state trajectory that approximates the continuous x. A simplified model of an unloaded electric motor will illustrate the method. The current i through the motor is related to the voltage v at its terminals by ˙ i = (v− iR)/L (5.5) where R is the resistance of the motor windings and L is their inductance. To have a simple calculation, take L = R= 1 so that Euler’s method, applied to Equation 5.5, gives i(t+ h)= i(t)+ h(v(t)− i(t))= (1− h)i(t)+ hv(t) (5.6) where i and v are known at time t. Observe that h 2 is necessary for the calculated i to remain bounded, and h 2 is needed to compute something similar to the186 HYBRID SYSTEMS a TABLE 5.1 Simulation of Equation 5.5 Using Euler’s Method ti,uv Event (0, 0) 0, 0 — init,0 (0.5, 1) 0, 0 — int (0.75, 0) — 1 in (0.75, 1) 0,1— ext (1.25, 1) 0.5, 1 — int (1.5, 0) — 0 in (1.5, 1) 0.625,0— ext (2, 1) 0.3125, 0 — int (2.25, 0) — 1 in (2.25, 1) 0.234375, 1 ext,final,0 a Where h= 0.5, v is a square wave with period 0.75 and amplitude 1, v((0, 0))= 0, and i = u= 0 initially. Output events have been omitted. correct continuous trajectory. Arbitrarily choosing h= 0.5 and v to be a square wave with period 0.75 and amplitude of 1, the discrete-event model calculates i at times 0.5, 0.75, 1.25, 1.5, 2, 2.25,... . These discrete instants capture precisely the moments when v changes and approximate i in the intervening spans. Table 5.1 shows the state and input trajectory for this simulation. 5.2 NETWORKS OF CONTINUOUS SYSTEMS The network structure used for discrete systems is also applicable to continuous systems (see, e.g., Refs. 149 and 157), but it is most valuable as a modeling aide and rarely used in simulations. Because continuous systems are concerned with the real or, more generally, complex numbers, algebraic operations can reduce every network to a set of ordinary differential equations, or if not, then this indicates a special difficulty with the model. If a hard case is not ill-defined, then it is a differential algebraic model; Cellier and Kofman 20 give a good introduction to simulation algorithms for these systems. The formulation in Section 5.1 can be applied to any network that is reducible to a set of ordinary differential equations. A continuous input u is incorporated into the state space of the model as follows. c If u is a function of x alone, replace f (x, u (x)) with z(x)= f (x, u (x)) and solve c c c the differential equation x ˙ = z(x). If u is a function of t as well, then add a new c state variable τ that satisfies τ ˙ = 1, define a new function z(x,τ)= f (x, u (x,τ)), c and solve the pair of differential equations x ˙ = z(x,τ) τ ˙ = 1 which is an input free system. Armed with this trick, continuous inputs do not require special handling by the simulation program.HYBRID MODELS AS DISCRETE-EVENT SYSTEMS 187 5.3 HYBRID MODELS AS DISCRETE-EVENT SYSTEMS A model of a sampled data system must provide output at fixed intervals, regardless of how the continuous equations are solved; similarly, a model of a threshold sensor must provide output at the moment of a threshold crossing, regardless of the timestep used for integration. Simulation methods for a hybrid system must therefore coordinate discrete-events and the timestep of the numerical integration scheme. A continuous model in its state space form can be approximated by an atomic, discrete-event model that encapsulates methods for event detection and numerical integration. If the continuous model is connected to discrete components, then inter- actions with those occur by discrete-events: discrete sensor readings, square waves, and so forth. This arrangement of a hybrid model, which splits continuous and discrete-event elements, is illustrated in Figure 5.2. m The hybrid model has a vector x ∈R of continuous state variables and a set Q of discrete states that describe switch settings, operating modes, and so forth. The differential equation that governs the model’s continuous evolution depends, in general, on the discrete state q and so has the form x ˙ = f (x, q) (5.7) The occurrence of autonomous events in this system is conditional on both q and x. Just as in a pure discrete-event system, the real time that will elapse prior to the next internal event is a function of the present state. The hybrid analog of the time advance function has the form m ∞ G :(R × Q)→ R 0 and events happen when G((x, q))= 0. Referring again to Figure 5.2, the continuous variables are exposed only through the event trajectory that is produced as output. If the last event was at time t and the k value of x at that time is x , then the value of x at the next event, some time h later, is k t +h k  x = x + f (x(t), q) dt (5.8) k+1 k t k FIGURE 5.2 Separation of a model into discrete-event and continuous components which interact by discrete-events.188 HYBRID SYSTEMS If the next event is autonomous, then h= G((x , q)) in the formulation given above. k Otherwise, the system was interrupted by an external event and h≤ G((x , q)). k Regardless, the value of x at the event can be obtained with any suitable numerical method and Equation 5.8 therefore poses no special difficulty. The effect of the event is a separate matter: x and q can change discontinuously in response to discrete input, the occurrence of an internal event, or both. Just as with a discrete-event system, three transition functions account for the three possible cases: m m ˆ δ :(R × Q)→ (R × Q) int m ∞ b m ˆ δ :(R × Q)× R × X → (R × Q) ext 0 m b m ˆ δ :(R × Q)× X → (R × Q) con define the response of the system to internal, external, and confluent events respec- tively. Similarly, the output function m b ˆ λ :(R × Q)→ Y generates discrete output when an autonomous event occurs. These “hatted” functions are an abbreviated description of the discrete-event system ˆ δ ((x , q))=δ ((x , q)) int k int k+1 b b ˆ δ ((x , q), e, x )=δ ((x , q), e, x ) ext k ext k+1 b b ˆ δ ((x , q), x )=δ ((x , q), x ) (5.9) con k con k+1 ta((x , q))= G((x , q)) k k ˆ λ((x , q))=λ((x , q)) k k+1 where x is calculated by Equation 5.8 with h= ta((x , q)) for internal and con- k+1 k fluent events and in the output function, and with h= e for external events. A simple model will illustrate the essential concepts. Consider the system shown in Figure 5.3, which has a single continuous variable constrained by surfaces at 0 and 1. It follows a line from 1 to 0, bounces off that surface, returns to 1, bounces again, and this repeats forever. This system has a single discrete variable d, with the range1,−1, that gives the continuous variable x its direction. The system can be interrupted at any time by an input u∈1,−1, which immediately changes the value of d. The model produces an output each time it encounters a constraining surface. Beginning with x ∈ 0, 1, the continuous trajectory between events is described by x ˙ = d (5.10)HYBRID MODELS AS DISCRETE-EVENT SYSTEMS 189 x 1 d=−1 d=1 d= −1 0 t (0,0) (1,0) (1,1) (2,0) (2,1) FIGURE 5.3 State trajectory of the bouncing hybrid system. the time remaining to the next internal event is  x if d=−1 G((x, d))= (5.11) 1− x if d = 1 and the system’s response to discrete-events is ˆ δ ((x, d))= (x,−d) int ˆ δ ((x, d), e, u)= (x, u) ext ˆ δ ((x, d), u)= (x, u) con ˆ λ((x, d))= x It is trivial to solve Equation 5.8 for this system; doing so and then substituting the 2 preceding equations into Equation 5.9 gives the discrete-event model ˆ δ ((x, d))=δ ((x , d))= ((d+ 1)/(2d),−d) int int k+1 ˆ δ ((x, d), e, u)=δ ((x , d), e, u)= (x+ ed, u) ext ext k+1 ˆ δ ((x, d), u)=δ ((x , d), u)= ((d+ 1)/(2d), u) con con k+1  x if d=−1 ta((x, d))= G((x, d))= 1− x if d = 1 ˆ λ((x, d))=λ((x , d))= (d+ 1)/(2d) k+1 2 The equation (d+ 1)/(2d) is zero when d=−1and x is descending. It is one when d =1and x is climbing.190 HYBRID SYSTEMS The output of this model when begun in state (0, 1) and fed, for instance, the input trajectory u(0, 0), (∞, 0)) with  1if t = (1.5, 0) u(t)=  otherwise is y(0, 0), (∞, 0)) with ⎧ ⎪ 1if t = (1, 0) or t = (2m, 0), m∈N ⎨ y(t)= 0if t = (2m+ 1, 0), m∈N−0 ⎪ ⎩  otherwise which can be confirmed by a table-driven simulation or logical argument (both are excellent exercises). The continuous trajectory, of course, is calculated only at discrete points. If a more detailed sampling of x is required, then G must be reformulated to invoke calculations at smaller intervals. This is accomplished by (1) introducing an upper ˆ limit for G equal to the largest sampling period h and (2) adding checks toλ and max ˆ δ for satisfaction of the conditions x = 1 and x = 0. The reformulated equations int are  minh , x if d=−1 max G((x, d))= minh , 1− x if d = 1 max  (x,−d)if x = 1∨ x = 0 ˆ δ ((x, d))= int (x, d) otherwise ˆ δ ((x, d), e, u)= (x, u) ext ˆ ˆ ˆ δ ((x, d), u)=δ (δ ((x, d)), 0, u) con ext int  x if x = 1∨ x = 0 ˆ λ((x, d))=  otherwise That the output is unchanged is easily verified. The state trajectory, however, is sampled more often: in the absence of input, with a frequency 1/h . max 5.4 NUMERICAL SIMULATION OF HYBRID SYSTEMS The formulation above assumes that G can be calculated without future knowledge of x. Lacking a direct expression for the location of events in time, they must be definedNUMERICAL SIMULATION OF HYBRID SYSTEMS 191 instead by locations in the model’s state space. This is reminiscent of the activity- scanning approach to discrete-event simulation (see, e.g., Ref. 157), but the event condition depends simultaneously on continuous and discrete variables. Specifically, autonomous events are described by a set of threshold functions g (x(t + h), q)= 0 1 k g (x(t + h), q)= 0 2 k . . . g (x(t + h), q)= 0 p k and G(x(t ), q) is defined implicitly as the smallest, nonnegative h that satisfies at k least one of these equalities or∞ if no such h exists. To illustrate this construction, consider the time advance of the bouncing sys- tem described in the previous section. Equation 5.11 is defined implicitly by the surface g((x(t + h), d))= x(t + h)− (d+ 1)/2= 0 k k and, because we know that x(t + h)= x + hd, this equation can be solved directly k k to obtain G. Two methods are required for simulation of the hybrid system: a numerical in- tegration method for advancing the continuous solution and a root-finding method to locate events between integration steps. These topics are covered by most in- troductory textbooks on numerical analysis (see, e.g., Refs. 71, 114, and 151) and in the specific context of hybrid simulation by a growing body of literature (see, e.g., Refs. 20 and 36). Almost any method is suitable, but those that readily ac- commodate changes in the step size and discrete adjustments of the state vari- ables are most desirable because discrete-events will frequently occasion the need for both. To demonstrate how a simulator is constructed, two particular methods are con- sidered here. Corrected Euler with error control is used for numerical integration. This method is part of the Runge–Kutta family. Given a step size h and present value x at time t , the next value at time t + h is computed with the two-stage k k k rule k = hf (x ) 1 k k = hf (x + k /2) 2 k 1 x = x + k k+1 k 2192 HYBRID SYSTEMS The magnitude of the errorǫ incurred by a step is approximately ǫ≈k − k 2 1 If after a step the error is too large, the step size is reduced and the integration retried. Conversely, if the error is very small, then the step size is increased to speed up the simulation. Having obtained the estimate ǫ using the step size h, the adjusted step size h that satisfies an error toleranceǫ is estimated by tol tol hǫ tol h ≈ tol ǫ For our purposes, the step size will begin at some maximum value h after each max event and then be reduced if necessary; in this case the new trial step is 0.8minh , h. tol Having selected a step size h that satisfies the error criteria and then tentatively advancing the trajectory, we now look to see if an event threshold lies between the old state x and the new, tentative state x . A simple, but often effective, method is to k k+1 interpolate each g in the interval t , t + h and look for a root of the interpolating j k k function. The two endpoints for the interpolation are g = g (x , q) j,k j k g = g (x , q) j,k+1 j k+1 There are three possibilities for this pair of points: 1. If g ≈ 0or g ≈ 0, then an event has been found and an internal transition j,k j,k+1 occurs at time t or t + h, respectively. k k 2. If both have the same sign, so that g g 0, then g does not cross zero j,k j,k+1 j and there is no event in the interval. 3. If the points have opposite signs, so that g g 0, then there is an event j,k j,k+1 in the interval. In possibilities 1 and 2, the time of the next autonomous event is known: it occurs now or at the next integration step. In possibility 3, the existence of an intermediate event is known, but its time of occurrence must still be found. A procedure for doing so is illustrated in Figure 5.4. To start, we take t = 0, which is possible because the system is time-invariant, and k approximate the function g (t) in the interval 0, h with the line j g (h− t)+ g t j,k j,k+1 g (t)≈ j hNUMERICAL SIMULATION OF HYBRID SYSTEMS 193 g(x(t),q) 0 g(x(t ),q) k h’ h t g(x(t ),q) k+1 FIGURE 5.4 Using linear interpolation to locate an event in time. Setting this equal to zero and solving for t gives the location of the event in time at hg j,k t ≈ g − g j,k j,k+1 Now the integration step is reduced to this estimate of the event time and a new x k+1 is calculated. If this state satisfies condition 1 or 2 in the preceding paragraph, then the search ends; the event is found or the interval has shrunk so as not to include it. Otherwise, these steps are repeated. To implement this procedure, the discrete-event model has, in addition to its the discrete state q, two copies of the continuous state: the present value x and tentative ′ next value x . For each g , there is a Boolean variable e that is true if g (x, q)≈ 0or k k k ′ ′ g (x , q)≈ 0 and false otherwise. The variable h keeps the time separating x and x . k ′ The model first looks ahead by a single integration step to calculate x . Next, it looks for events at the ends of the interval 0, h and for events within it that might have been missed. There are four cases. The simplest is where the g neither are zero at k the ends of the interval nor change sign. In this case, there are no events and, barring ′ external input, the state h units of time later is x . The other three cases involve discrete-events. If at the left end of the interval any of the g are withinǫ of zero, then the e are set accordingly and an internal event k tol k occurs immediately. Conversely, if this is true at the right end of the interval and none ′ of the g have otherwise changed sign, then an event occurs at h, with state x , and the k e are again set accordingly. Otherwise, an event occurs inside the interval and must k ′ be located with greater precision—a new h is estimated, a new x calculated, and the process repeated. Algorithm 5.1 shows this procedure in detail, using the corrected Euler and linear interpolation methods described above. With Algorithm 5.1, the discrete-event model that simulates a hybrid system can ′ be written compactly as follows. Let s= ((x, q), x , e ¯, h) be the state of the model where e ¯ is the vector of event flags e ,..., e . The error toleranceǫ and maximum 1 p tol194 HYBRID SYSTEMS Input: x,q,h ,ǫ max tol ′ Output: x ,h,e ,..., e 1 p 1 Tentatively advance x to satisfy the error criteria 2 h← h max 3 repeat 1 4 k ← hf (x, q), k ← hf (x+ k , q) 1 2 1 2 ′ 5 x ← x+ k ,ǫ←k − k 2 2 1 6 if ǫ ≤ǫ then h← 0.8minh, hǫ /ǫ tol tol 7 untilǫ≤ǫ tol 8 Look for events in 0, h 9 foreach k∈ 1, p do g ← g (x, q), e ←g≤ǫ k k k k tol 10 Found an event at the start of the interval ′ 11 if (∃e )(e = true) then h← 0, x ← x, return k k 12 repeat 13 foreach k∈ 1, p do ′ ′ ′ 14 g ← g (x , q), e ←g≤ǫ k k tol k k ′ ′ 15 if g g 0 then h← minh, hg /(g − g ) k k k k k 16 end 1 17 k ← hf (x, q), k ← hf (x+ k , q) 1 2 1 2 ′ 18 x ← x+ k 2 ′ 19 until (∀k)(0 g g ∨ e ) k k k Algorithm 5.1 Locate the internal event and tentatively advance x to it. integration step h are fixed parameters. Let F be a function that applies Algorithm max ′ 5.1 to the state; F(s) is equal to s except for x , e ¯, and h which the algorithm computes. ′ ′ Define also the function I ((x, q),τ)= (x , q), where x is x advanced by the step sizeτ. The discrete-event model is  ′ ′ ˆ F((δ ((x , q)), x , e ¯, h)) if (∃e )(e = true) int k k ′ δ (((x, q), x , e ¯, h))= int ′ ′ F(((x , q), x , e ¯, h)) otherwise ′ b b ′ ˆ δ (((x, q), x , e ¯, h), e, u )= F((δ (I ((x, q), e), e, u ), x , e ¯, h)) ext ext ′ b ′ b ′ ˆ δ (((x, q), x , e ¯, h), u )= F((δ ((x , q), u ), x , e ¯, h)) (5.12) con con ta(s)= h  ′ ˆ λ((x , q)) if (∃e )(e = true) k k λ(s)=  otherwise To summarize, at each event the continuous variable x is advanced to the present time. If there is an input or if an event threshold is reached, then the discrete transition function of the hybrid system is applied and, in the latter case, an output is also produced. The time advance at each event is the smaller of the integration step size and the time to the next internal event of the hybrid system.NUMERICAL SIMULATION OF HYBRID SYSTEMS 195 v c R s v + C s − R l FIGURE 5.5 A model of a circuit with two discrete elements. The electric circuit shown in Figure 5.5 will illustrate the simulation procedure, and the outcome of the simulation is listed in Table 5.2. This circuit has a load that operates when the voltage v is greater than the diode’s threshold voltage. A capacitor stores c energy to drive the load for short periods of time when the main power is disconnected. This capacitor charges rapidly when the load is disconnected, and it trickle-charges while the load is operating. The model has two discrete variables; s∈0, 1 for the switch with s= 0 open (not conducting) and s= 1 closed (conducting) and likewise for the diode; its state is indicated by d. The switch is an input to the system, and the diode state is determined by the capacitor voltage. The voltage v is the model’s only continuous state variable. It evolves in each c mode as follows: ⎧ ⎪ 0if s= 0, d = 0 ⎪ ⎪ ⎪ ⎪ v ⎪ c ⎪ ⎪ − if s= 0, d = 1 ⎪ ⎨ CR l v ˙ = v − v (5.13) c s c ⎪ if s= 1, d = 0 ⎪ ⎪ CR ⎪ s ⎪   ⎪ ⎪ 1 v − v v ⎪ s c c ⎪ ⎩ − if s= 1, d = 1 C R R s l The diode conducts when v exceeds the threshold voltage v and stops conducting c op when the voltage drops below v , where v v . A single function defines both cl op cl cases:  v − v if d = 0 c op g(v , d)= c v − v if d = 1 c cl Significantly, the diode does not open and close at the same voltage level. If this were the case, then the model would be illegitimate, switching endlessly between its two states when v reached a threshold value. The model’s output is the state of the diode. c196 a TABLE 5.2 Discrete-Events in the Circuit Simulation ′ tv (s, d) eh ¯ x uy Event c (0, 0) 0 (1, 0) False 0.20.18  init,0 (0.2, 1) 0.18 (1, 0) False 0.20.3276  int (0.4, 1) 0.3276 (1, 0) False 0.20.448632  int (0.6, 1) 0.448632 (1, 0) True 0.1035162642 0.5027534314  int (0.7035162642, 0) 0.448632 (1, 0) True 0.1035162642 0.5027534314  1 out (0.7035162642, 1) 0.5027534314 (1, 1) False 0.20.5018723334  int (0.9035162642, 1) 0.5018723334 (1, 1) False 0.20.5012731867  int (1, 0) 0.5018723334 (1, 1) False 0.20.5012731867 1  in (1, 1) 0.5015458935 (0, 1) False 0.20.4112676327  ext (1.2, 1) 0.4112676327 (0, 1) False 0.20.3372394588  int (1.4, 1) 0.3372394588 (0, 1) False 0.20.2765363562  int (1.6, 1) 0.2765363562 (0, 1) True 0.1013635251 0.2499263016  int (1.7013635251, 0) 0.2765363562 (0, 1) True 0.1013635251 0.2499263016  0 out (1.7013635251, 1) 0.2499263016 (0, 0) False 0.20.2499263016  int (1.9013635251, 1) 0.2499263016 (0, 0) False 0.20.2499263016  int,final,0 a The location of internal events in time can be calculated directly or with Algorithm 5.1. Output events for which y= have been omitted.NUMERICAL SIMULATION OF HYBRID SYSTEMS 197 The discrete transition functions for this circuit are ˆ δ ((v , s, d))= (v , s,¬d) int c c ˆ δ ((v , s, d), e, u)= (v , u, d) ext c c ˆ ˆ ˆ δ ((v , s, d), u)=δ (δ ((v , s, d)), 0, u) con c ext int c ˆ λ((v , s, d))=¬d c 1 To keep the numbers simple, use R = R = 1, C = 1, v = 1, v = , and v = l s s cl op 4 1 and begin with d = 0, s= 1, and v = 0. The discrete-event model uses h = 0.2 c max 2 andǫ= 0.05. The capacitor begins to charge, and until the diode closes, the discrete equations (i.e., corrected Euler equations) that simulate v are c 2 2 v = v (1− h+ h /2)+ h− h /2 c,k+1 c,k the error at each step is 2 ǫ=(1− v )h /2 c,k and the interpolating function for g, solved for time, is (v − 0.5)h c,k t = v − v c,k c,k+1 causing the discrete-event model to undergoe the first three internal events listed in Table 5.2. The diode closes at the fourth event when the real time is approximately 0.7. ˆ ˆ Because this internal event coincides with an event surface,λ andδ determine the int output and next state of the model. In its new mode, v is simulated with the difference c equation v = v + h(h− 1)(2v − 1) c,k+1 c,k c,k the error at each step is 2 ǫ=(1− 2v )h c,k the interpolated g solved for time is (v − 0.25)h c,k t = v − v c,k c,k+1 and the discrete-event model evolves in this way until the switch is opened. If the simulation closed the diode at exactly v = 0.5, then both the actual system and its c198 HYBRID SYSTEMS approximation would remain at that value. Instead, the simulated voltage overshoots the mark and descends slowly back to the correct solution. At t = (1, 0) the switch is opened. This external event advances v by the elapsed ˆ time and invokes theδ function of the hybrid system to change its state. In this new ext configuration, the discrete-event system simulates v by c 2 v = (1− h+ h /2)v c,k+1 c,k the error at each step is ′ 2 ǫ =h v /2 c,k and the interpolation function for g, solved for time, is as before. The capacitor is drained to the cutoff point after another three events, and the diode opens at the real time of t ≈ 1.7. Now the switch and diode are open and the equations governing the simulation are v = v c,k+1 c ǫ= 0 t=∞ and so time advances with h= h = 0.2 but nothing actually changes until the max switch is closed again. Working through this example, it quickly becomes apparent that an explicit formula for G saves a great deal of time and effort. When G can be calculated directly from the system’s state, and so the location of the next event in time is known explicitly, the resulting change is commonly called a time event and G a time event function. If the next event time must be inferred from event surfaces in the state space, then the resulting change is called a state event and the g are state event functions. State k events require much more computational effort on the part of the simulator, and so result in slower-running simulations. Time events are therefore preferable when they can be constructed. 5.5 A SIMULATOR FOR HYBRID SYSTEMS The Hybrid class extends the Atomic class of the discrete-event simulator from Chap- ter 4 to implement the approximating system constructed above. Algorithms for finding event surfaces in time and for solving the differential equation are encapsu- lated in two classes: the event locator and the ode solver. The differential function, state and time event functions, and discrete transition functions of the hybrid systemA SIMULATOR FOR HYBRID SYSTEMS 199 linear_event_locator linear_event_locator linear_event_locator FIGURE 5.6 UML diagram of the classes that constitute the hybrid simulator. are contained in a third class: the ode system.The Hybrid class uses these three supporting classes to simulate a hybrid model. Figure 5.6 shows their relationships to each other and the simulation engine. The event locator uses the single method find events, which looks for state events in a prescribed interval. Its arguments are the continuous state q of the hybrid model,200 HYBRID SYSTEMS the width h of the interval to search, the state q at the end of the interval, and an end integration method for calculating states within the interval. The method returns true if an event is found and false otherwise. In the positive case, the time of the event is written to the supplied h, the elements of the array events are set to indicate which threshold function caused the event, and the state of the model at the event time is written to q . In the negative case, the elements of the events array are set to false end but h and q retain the values supplied by the caller. end The ode solver has two methods for advancing the continuous state variables. The ′ integrate method takes the system from the state q at time t to the state q at time t+ h where h≤ h . Integration methods with error control may use a value smaller max than h to satisfy their tolerance for error. The method returns the step size that max ′ is actually used, and the final state q is copied to q before returning. The advance method is identical to the integrate method except that it is guaranteed to advance the continuous trajectory to the end of the requested interval. The hybrid model itself is implemented by extending the ode system class. It has methods for computing the state event functions, the time event function, and the differential function. These are used by the event locator and ode solver classes. There are also four methods for implementing the discrete dynamics: external event, internal event, confluent event, and output func. These correspond to the hatted tran- sition and output functions of the hybrid system. All of these methods accept an array q that holds the model’s continuous state. This is maintained by the Hybrid class, which supplies it to the ode system as needed. A model’s discrete state, however, is captured by the member variables of its implementing class. Also observe that the discrete transition functions and output functions are supplied with an array of Boolean values; this array holds the e flags plus one extra in the last array position k to indicate a time event. The gc output method is used to reclaim objects created by the output func. Equation 5.12 is implemented by the Hybrid class. Its private tentative step method implements the F function, which represents Algorithm 5.1 in the def- inition of the discrete-event model. The integration function I is realized with the ode solver, which is also used by the event locator to find state events (i.e., roots of the g ) in the course of taking a tentative step. Time events are scheduled k with the ode system’s time event func method, which returns the time remaining to the nearest, explicitly known autonomous event. The time advance selected by the tentative step method is the smallest timestep used by the ode solver, the near- est state event found by the event locator, and the next time event reported by the ode system. The member variables of the Hybrid class are exactly the state variables of the discrete-event model, except for the extra Boolean variable event exists, which saves a search through the e s by being true if any of them is true and false k otherwise. The Hybrid Class 1 ifndef _adevs_hybrid_h_ 2 define _adevs_hybrid_h_ 3 include algorithmA SIMULATOR FOR HYBRID SYSTEMS 201 4 include "adevs_models.h" 5 6 namespace adevs 7 8 9 template typename X class ode_system 10 11 public: 12 // Make a system with N state variables and M state event functions 13 ode_system(int N_vars, int M_event_funcs): 14 N(N_vars),M(M_event_funcs) 15 // Get the number of state variables 16 int numVars() const return N; 17 // Get the number of state events 18 int numEvents() const return M; 19 // Copy the initial state of the model to q 20 virtual void init(double q) = 0; 21 // Compute the derivative for state q and put it in dq 22 virtual void der_func(const double q, double dq) = 0; 23 // Compute the state event functions for state q and put them in z 24 virtual void state_event_func(const double q, double z) = 0; 25 // Compute the time event function using state q 26 virtual double time_event_func(const double q) = 0; 27 // The internal transition function 28 virtual void internal_event(double q, 29 const bool state_event) = 0; 30 // The external transition function 31 virtual void external_event(double q, double e, 32 const BagX& xb) = 0; 33 // The confluent transition function 34 virtual void confluent_event(double q, const bool state_event, 35 const BagX& xb) = 0; 36 // The output function 37 virtual void output_func(const double q, const bool state_event, 38 BagX& yb) = 0; 39 // Garbage collection function 40 virtual void gc_output(BagX& gb) = 0; 41 virtual ode_system() 42 private: 43 const int N, M; 44 ; 45 46 template typename X class ode_solver 47 48 public: 49 ode_solver(ode_systemX sys):sys(sys) 50 // Take an integration step from state q of at most size h and 51 // return the step size that was actually used. Copy the result of 52 // the integration step to q.