Linear Phase-portrait Approximations for ? ?? Nonlinear Hybrid Systems Thomas A. Henzinger
Howard Wong-Toi
Department of Computer Science, Cornell University, Ithaca, NY 14853 (tahjhoward)@cs.cornell.edu
Abstract. We use linear hybrid automata to de ne linear approxima-
tions of the phase portraits of nonlinear hybrid systems. The approximating automata can be analyzed automatically using the symbolic model checker HyTech. We demonstrate the technique through the study of predator-prey systems, where we compute population bounds for both species. We also identify a class of nonlinear hybrid automata for which linear phase-portrait approximations can be generated automatically.
1 Introduction Hybrid systems combine discrete and continuous dynamics, and their analysis requires techniques from both computer science and control theory. Computer scientists typically model hybrid systems as discrete programs augmented with continuous environment variables undergoing simple dynamics, whereas control theorists typically study complex behaviors of continuous parameters within a simple discrete structure of control modes. This paper attempts to narrow the gap between the two disciplines by pushing the computer science end of the frontier. Algorithmic program analysis techniques can be extended for analyzing certain properties, such as reachability, of linear hybrid automata [ACHH93, ACH+ 95]. These automata model hybrid systems with linearity restrictions on discrete jumps (linear inequalities between sources and targets of jumps) and continuous ows (dierential inequalities of the form Ax_ b). Algorithmic analysis techniques for linear hybrid automata have been implemented in HyTech [AHH93, HHWT95b] and used to check properties of distributed real-time protocols [HH95b, HW95]. The de nition of linearity for hybrid automata is dierent from that for linear systems in systems theory. For instance, linear hybrid automata cannot model directly continuous ows of the form x_ = x. Since most control applications ? This research was supported in part by the ONR YIP award N00014-95-1-0520,
by the NSF CAREER award CCR-9501708, by the NSF grants CCR-9200794 and CCR-9504469, by the AFOSR contract F49620-93-1-0056, and by the ARPA grant NAG2-892. ?? A preliminary version of this paper appeared in Hybrid Systems III: Veri cation and Control (R. Alur, T.A. Henzinger, E.D. Sontag, eds.), Lecture Notes in Computer Science 1066, Springer-Verlag, 1996, pp. 377{388.
cannot be modeled directly using linear hybrid automata, we employ the idea of abstract interpretation [CC77] to establich properties of hybrid systems using linear hybrid automata as approximations. Suppose that we are given a complex hybrid system A and a property . Typically, checking if A satis es is dicult. Our goal is to perform the simpler task of checking, automatically, if a linear hybrid automaton B that approximates A satis es . We choose B such that if holds for B , then holds also for A. If, on the other hand, is not satis ed by B , then we cannot infer anything about A, and instead must try to nd a closer approximation B 0 of A. Our main thesis is that linear automata approximations are adequate for checking certain interesting properties of complex hybrid systems. The approximation of hybrid systems using linear hybrid automata has been advocated previously. Two translation techniques from nonlinear to linear hybrid automata are suggested in [HH95a]. The rst translation, called clock translation, requires that the continuous ows are described by solvable dierential equations, and replaces continuous variables by clock variables with constraints that are equivalent to the constraints on the original variables. The second translation, called rate translation, approximates continuous behavior with piecewiseconstant bounds on the rst derivatives. The rate translation conservatively approximates a set of continuous trajectories using piecewise-linear envelopes over a state space that is partitioned into linear regions. The method is applied in two steps. First, every control mode of a hybrid system is split into several copies, each of which corresponds to a region of the partitioned state space. Second, within every new control mode, the nonlinear dynamics is replaced by a linear dynamics: the rst derivative of each variable is bounded by two constants that correspond to the extremal rates (with respect to time) of the variable in the corresponding region of the state space. In this paper, we improve on the rate translation technique by suggesting linear phase-portrait approximations that are strictly more accurate. These approximations are motivated by a time-invariant view of hybrid automata [Hen95]. They dier from the rate translation by approximating directly over the vector eld of ow tangents, thus taking into account the relative rates of variables. This leads to more accurate approximations that still allow algorithmic analysis using HyTech. (The rate translation, which approximates the rate of each variable independently, was designed for automatic analysis using the restricted version of HyTech available at the time. Recently HyTech has been extended to permit arbitrary linear ow elds [HHWT95a], thereby enabling us to approximate with a more expressive form of linear hybrid automata.) Since rate translations are linear phase-portrait approximations, the linear phase-portrait approximation technique inherits the asymptotic completeness property of rate translations, which states that well-behaved nonlinear hybrid automata can be approximated arbitrarily closely [HH95a]. Outline We rst de ne the model of nonlinear hybrid automata (Section 2) and formalize the approximation of hybrid automata (Section 3). We then de ne linear hybrid automata, the class that can be analyzed algorithmically, and 2
y_
y
9
x=1 ^ y=2 q0 1 x 3 y = 7 12 y_ = xx_ ^ x_ = x
q1 y_ 3x
1
(a) A, nonlinear hybrid automaton Fig. 1.
1 3 x (b) Flow eld for location q0 of A
The hybrid automaton A
1 3 x_ (c) Endpoints of
ow vectors for q0
linear phase-portrait approximations (Section 4). Section 5 applies linear phaseportrait approximation to a simple predator-prey ecology and veri es a control strategy for maintaining population bounds. We do not claim that linear phaseportrait approximation can replace standard analytic techniques for nonlinear behavior. On the contrary, our examples show that some understanding of the solutions of given dierential inequalities is often necessary for obtaining useful approximations, and better understanding leads to better approximations. Section 6 presents as simple exception to this rule the class of pseudo-linear hybrid automata, for which linear phase-portrait approximations can be generated automatically. Related work Phase portraits have been studied extensively in the literature on dynamical systems [HS74, Arn83]. Typically, researchers concentrate on the continuous dynamics of a system and analyze extremely nontrivial properties, such as stability and convergence. Our work diers in two respects. First, we consider products of nondeterministic dynamical systems with discrete transition structures. Second, our goal is to analyze and derive simple properties of such systems automatically. The theory of hybrid automata is presented in [Hen95, HKPV95]. In [HRP94, HH95c], abstract interpretation techniques are used to provide linear approximations of linear hybrid systems, whereas here we approximate nonlinear hybrid systems.
2 Hybrid Automata We de ne hybrid automata, which are used to model systems consisting of mixed discrete and continuous components. Informally, a hybrid automaton consists of a nite set X of real-valued variables and a labeled multigraph (V; E ). The edges E represent discrete jumps and are labeled with guarded assignments to variables in X . The vertices V represent continuous ows and are labeled with dierential inequalities over the rst derivatives of the variables in X . The state of the automaton changes either instantaneously when a discrete jump occurs (\transition step") or, while time elapses, through a continuous ow (\time step"). 3
De nition A hybrid automaton A consists of the following components: Variables. A nite ordered set X = fx ; : : :; xng of real-valued variables. For 1
example, the automaton A in Figure 1, part (a), has two variables, x and y. A valuation is a point (s1 ; : : :; sn) in the n-dimensional real space Rn, or equivalently, a function that maps each variable xi to a value ai . Each predicate over X de nes a set [ ] Rn of valuations, such that s 2 [ ] if [X := s] is true. Locations. A nite set V of vertices called locations, used to model control modes. For example, the automaton A has two locations, q0 and q1. A state (v; s) consists of a location v 2 V and a valuation s 2 Rn. A region is a set of states. Invariants. A labeling function inv that assigns to each location v 2 V a predicate inv (v) over X , the invariant of v. For example, location q0 of the automaton A has the invariant 1 x 3. Invariants of the form true are omitted from gures, as in location q1. A state (v; s) is admissible if s 2 [ inv (v)]]. The automaton control may reside in location v only while the invariant inv (v) is true. Therefore invariants can be used to force transitions. Initial conditions. A labeling function init that assigns to each location v 2 V a predicate init (v) over X , the initial condition of v. For example, location q0 of the automaton A has the initial condition x = 1 ^ y = 2. Initial conditions of the form false are omitted from gures, as in location q1. A state (v; s) is initial if it is admissible and s 2 [ init (v)]]. The automaton control may start in location v when the initial condition init (v) and the invariant inv (v) are true. Flow conditions. A labeling function that assigns a ow condition to each location v 2 V . The ow condition (v) is a predicate over X [ X_ , where X_ = fx_ 1; : : :; x_ ng. The dotted symbol x_ refers to the rst derivative of the variable x with respect to time, i.e. x_ = dx=dt. For example, location q1 of the automaton A has the ow condition y_ 3x. While the automaton control resides in location v, the variable values change along dierentiable trajectories whose rst derivatives satisfy the ow condition. Formally, for each real 0, we de ne the binary time-step relation ! such that (v; s)! (v0 ; s0 ) if v0 = v, and there is a dierentiable function : [0; ] ! Rn such that { the endpoints of the ow match those of , i.e. (0) = s and () = s0 ; { the invariant is satis ed, i.e. for all reals t 2 [0; ], (t) 2 [ inv (v)]]; and { the ow condition is satis ed, i.e. for all reals t 2 (0; ), (v)[X; X_ := (t); _(t)] is true, where _(t) = (1 (t)=dt; : : :; n (t)=dt). Transitions. A nite multiset E of edges called transitions. Each transition (v; v0 )
identi es a source location v 2 V and a target location v0 2 V . For example, the automaton A has a single transition, (q0 ; q1). Jump conditions. A labeling function jp that assigns a jump condition to each transition e 2 E . The jump condition jp (e) is a predicate over X [ X 0 , where X 0 = fx01; : : :; x0ng. While the unprimed symbol x refers to the value of the variable x before the jump, the primed symbol x0 refers to the value of x after 4
the jump. Formally, for the transition e = (v; v0 ), we de ne the binary transitione step relation ! on the admissible states such that (v; s)!e (v0 ; s0 ) if jp (e)[X; X 0 := s; s0 ] is true. When V writing jump conditions in the gures of this paper, we omit the conjunct x2X x0 = x. For example, the transition (q0; q1) of the automaton A has the jump condition y = 7 21 ^ x0 = x ^ y0 = y. Veri cation Let A be a hybrid automaton with n variables. The state space V Rn of A is the set of admissible We de ne the binary successor S states. [S ! e . For a region R , relation !A on the state space as 2 0 ! e2E we de ne the successor region post (R) to be the set of states that can be reached from some state in R via a single time step or transition step, i.e. post (R) = f0 j 9 2 R such that !A 0g. The region post (R) reachable from R is the set of states S that can be reached from R by a nite number of steps, i.e. post (R) = i0 post i (R). In practice, many veri cation tasks about hybrid systems can be posed in a natural way as reachability problems. Often, the system can be composed with a monitor process that \watches" the system and enters a violation state whenever a system trajectory violates a given safety or timing requirement. For a hybrid automaton A, the reachable region reach (A) is the region that is reachable from the set of initial states. The automaton A is correct with respect to a set S of violation states if the reachable region contains no violation states, i.e. reach (A) \ S is empty. Example 1. Consider the hybrid automaton A of Figure 1. Throughout the paper, we consider the task of verifying that location q1 is not reachable, i.e. reach (A) contains no state of the form (q1; v). The ow condition for location q0 is x_ = x ^ y_ = xx_ . We may solve these simple equations to see the solutions are of the form y = 21 x2 + k. From the initial valuation (1; 2), the valuation (3; 6) is reached, and we see that location q1 is not reached.
R
3 Approximation A hybrid automaton may be approximated by another automaton by relaxing the invariants, the initial conditions, the ow conditions, or the jump conditions. Flow elds and jump elds Given a set X of n real-valued variables, a
ow vector for X is a npoint in Rn. A (nondeterministic) ow eld for X is a map ow : Rn ! 2 that assigns to each valuation a set of ow vectors. Intuitively, if the vector (a1 ; : : :; an) is in ow (s), then the rate of change with respect to time of each variable xi at the valuation s 2 Rn can be ai . We de ne an order on ow elds such that ow 1 ow 2 if for all valuations s 2 Rn,
ow 1 (s) ow 2(s). With each location of a hybrid automaton, we associate a ow eld in the natural way: the ow eld for location v assigns to each valuation s 2 [ inv (v)]] the set fa j (v)[X; X_ := s; a] is trueg of ow vectors. For example, for the automaton A of Figure 1, the ow eld for location q0 is given by ow (x; y) = f(x; x2) j 1 x 3g. The automaton ow eld ow (A) of a hybrid automaton A is the function mapping each location to its ow eld.
R
5
R
Given a set X of n real-valued variables, a jump vector for X is a point in Rn. A (nondeterministic) jump eld for X is a map jump : Rn ! 2 n that assigns to each valuation a set of jump vectors. Intuitively, if the vector b is in jump (s), then the valuation s 2 Rn can be instantaneously incremented by b. We de ne an order on jump elds such that jump 1 jump 2 if for all valuations s 2 Rn, jump 1 (s) jump 2 (s). With each location of a hybrid automaton, we associate a jump eld in the natural way: the jump eld for transition e = (v; v0 ) assigns to each valuation s 2 [ inv (v)]] the set fs0 ? s j jp (e)[X; X 0 := s; s0 ] is true and s0 2 [ inv (v0 )]]g of jump vectors. The automaton jump eld jump (A) of a hybrid automaton A is the function mapping each transition to its jump eld. The automaton ow eld and the automaton jump eld are collectively referred to as the phase portrait of a hybrid automaton [Hen95]. Phase-portrait approximation A hybrid automaton B approximates the hybrid automaton A if reach(A) reach (B ). If the approximating automaton B is correct with respect to a set of violation states, then so is the approximated automaton A. The following lemma states a sucient condition for B to approximate A. Let two hybrid automata be compatible if they have the same set of variables, locations, and transitions. For compatible hybrid automata A and B , we de ne inv (A) inv (B ) if for all locations v 2 V , inv (A)(v) inv (B )(v); init (A) init (B ) if for all locations v 2 V , init (A)(v) init (B )(v); ow (A) ow (B ) if for all locations v 2 V , ow (A)(v) ow (B )(v); and jump (A) jump (B ) if for all transitions e 2 E , jump (A)(e) jump (B )(e). We de ne A B if (1) inv (A) inv (B ), (2) init (A) init (B ), (3) ow (A) ow (B ), and (4) jump (A) jump (B ). The following lemma follows from the fact that is a simulation over compatible hybrid automata [HH95a].
Lemma 1. Let A and B be compatible hybrid automata. If A B, then reach (A) reach (B ).
For approximating a hybrid automaton A, it may be bene cial to divide the state space of A into subsets, and then approximate over each subset of the state space. A partitioning P for A is a function that maps each location v of A to a nite set finv 1; : : :; inv k(v)g of invariants, such that (1) the disjunction of the invariants is equivalent to the invariant of v, and (2) each valuation in the interior of [ inv (v)]] either lies in the interior of [ inv i ] for some 1 i k(v), or lies in the intersection [ inv i ^ inv j ] for some 1 i < j k(v). The partitioning P leads to the partitioned automaton P (A). For each location v of A, P (A) has the locations v1; : : :; vk(v), with inv (vi ) = P (v)i , init (vi ) = init (v), and
(vi ) = (v) for all 1 i k(v). For each transition e = (v; v0 ) of A, P (A) has the transitions eij = (vi ; vj0 ) for 1 i k(v) and 1 j k(v0 ), with jp (eij ) = jp (e). In addition, P (A) has transitions that let control pass freely between locations of P (A) that correspond to a common location of A, i.e. for each location v of A and Vall 1 i 6= j k(v), P (A) has the transition (vi ; vj ) with the jump condition x2XA x0 = x. 6
We de ne the function homP : VP (A) ! VA such that for all locations vi of P (A), homP (vi ) = v. The function extends to states, and sets of states, in the natural way, e.g. homP (vi ; s) = (v; s). Lemma 2. Let A be a hybrid automaton A, and let P be a partitioning for A. Then homP (reach (P (A))) = reach (A). A hybrid automaton B is a phase-portrait approximation of the hybrid automaton A with partitioning P if B is compatible with P (A) and P (A) B . It follows from the previous two lemmas that reach (A) homP (reach (B )). This gives an upper bound on the reachable region reach (A). Iterating approximations Bounds on the reachable region can be used to obtain tighter approximations of ow elds. Suppose that B is a phase-portrait approximation of A with the trivial partitioning. Reachability analysis using B may show that only some proper subset of each invariant is reachable. A second phase-portrait approximation may then be generated, where, in essence, the approximating ow elds need only include ow vectors for the valuations that are reachable in the rst approximation. Iterating this procedure may lead to more and more accurate approximations. Similar reasoning may be used to generate tighter jump elds in the approximating automata.
4 Linear Approximation Linear hybrid automata Linear hybrid automata are a subclass of hybrid automata whose reachability problem can be analyzed algorithmically. A linear expression over a set X of real-valued variables is a linear combination of variables with rational coecients. A linear inequality is an inequality between linear expressions; a convex linear predicate is a nite conjunction of linear inequalities; and a linear predicate is a nite disjunction of convex linear predicates. A hybrid automaton with the set X of variables is linear if (1) for all locations v, the invariant inv (v) is a convex linear predicate over X , (2) for all locations v, the initial condition init (v) is a convex linear predicate over X , (3) for all locations v, the ow condition (v) is a convex linear predicate over X_ , and (4) for all transitions e, the jump condition jump (e) is a convex linear predicate over X [ X 0 (the convexity conditions can be relaxed by splitting locations and transitions). In ow conditions, linear dependencies between the rates of variables can be expressed, although the ow eld must be independent of the values of the variables. For example, the automaton A00 in part (b) of Figure 2 is a linear hybrid automaton. Algorithmic analysis techniques for linear hybrid automata are implemented in tools such as HyTech [AHH93, HH95b, HHWT95a]. In particular, HyTech includes a semidecision procedure for solving the reachability problem of linear hybrid automata. Rate translation The rate translation [HH95a] approximates nonlinear hybrid automata by linear hybrid automata. It consists of two steps: partitioning the state space within each location, and replacing nonlinear dynamics within each 7
x=1 ^ y=2 q00 1x2 y_ 2 [1; 4] ^ x_ 2 [1; 2]
x=1 ^ y=2
y = 7 21 q01 2 x 3 y = 7 21 y_ 2 [4; 9] ^ x_ 2 [2; 3]
q1
q1
(b) A00 , a linear phase-portrait approximation for A
(a) A0 , a rate translation for A Fig. 2.
q0 1x3 y = 7 12 y_ 6(x_ ? 3=2) ^ y_ 4(x_ ? 3=4) ^ y_ x_
Linear phase-portrait approximations for A
region of the partitioned state space by piecewise-constant bounds on derivatives. For simplicity, we assume that all invariants, initial conditions, and jump conditions of the given nonlinear hybrid automaton A are convex linear predicates. A linear hybrid automaton B is a rate approximation of A if (1) B is a phase-portrait approximation of A, and (2) for all locations v of B , (v) is of the V form xi 2XB li x_ i ui, where li and ui are rational or in nite. Thus the ow conditions for each location are independent of the valuations, and rectangular, i.e. they consist of independent upper and lower bounds on the derivatives of all variables. Given a partitioning P , if there exists a minimal rate approximation of A with P with respect to , it is referred to as the rate translation rate P (A) of A with P . Example 2. Consider automaton A of Figure 1. The invariant at q0 requires the value of x to remain between 1 and 3. The ow vectors are of the form (x; x2), for all 1 x 3. Their endpoints, with respect to the origin, are shown in part (c) of the gure. Figure 2, part (a), shows the rate translation A0 for the hybrid automaton A with the partitioning that splits the invariant along the line x = 2 in location q0 . Since x_ may vary only from 1 to 2 while the automaton control remains in location q00 , the ow condition in A0 for q00 is x_ 2 [1; 2] ^ y_ 2 [1; 4]. Computing the time-step successors from the valuation (1; 2) shows that all valuations with x = 2 and 2 21 y 6 are reachable. Similarly, reachability for location q01 shows that all valuations with x = 3 and 3 65 y 10 21 are reachable, and thus location q1 is reachable. Indeed, no partitioning that divides q0 into two locations avoids trajectories to q1.
Linear phase-portrait approximation We generalize rate approximation in two respects. First, rather than individually approximating the rate of each variable, we approximate the phase portrait of the joint evolution of all variables. This extension, which enables approximations with fewer trajectories than rate translation, has been made viable by our recent reimplementation of HyTech [HHWT95a], which now allows the modeling of nonrectangular linear ow elds. Second, rather than partitioning with rectangular grids over Rn [HH95a], we partition using arbitrary linear regions. This, again, enables approximations with fewer trajectories than rate translation. A hybrid automaton B is a linear phaseportrait approximation of A if (1) B is a phase-portrait approximation of A, and 8
x0^y 0 x_ = (A ? By ? x)x ^ y_ = (Cx ? D ? y)y
x = x0 ^ y = y 0
Fig. 3.
Predator-prey hybrid automaton
(2) B is linear. Given a partitioning P , if there exists a minimal linear phaseportrait approximation of A with P with respect to , it is referred to as the linear phase-portrait translation phase P (A) of A with P . In many cases, linear phase-portrait approximations can be obtained by approximating, for each location, the ow eld using a convex linear predicate containing the convex hull of the set of ow vectors. Example 3. The automaton A00 in Figure 2, part (b), is a linear phase-portrait
approximation of A (with the trivial partitioning). It has the ow condition y_ x_ ^ y_ 6(x_ ? 3=2) ^ y_ 4(x_ ? 3=4), de ning the convex hull of the ow vectors f(x; x2) j 1 x 3g at q0 . Reachability analysis reveals that the maximal value of y obtained in location q0 of A00 at x = 3 is 8. The location q1 is reachable. While this approximation has fewer trajectories than the rate approximation A0 , it is still insucient to establish that location q1 is not reachable. Both the rate approximation A0 and the linear phase-portrait approximation A00 above can be improved through use of a ner partitioning. For linear-phase portrait approximation, using the partitioning of the rate translation A0 shows that y never exceeds 7 in location q0, and hence establishes that q1 is unreachable. For rate approximation, dividing the invariant into three equal parts | x 2 [1; 5=3], x 2 [5=3; 7=3], and x 2 [7=3; 3]| instead of two is still insucient to bound y suciently, but dividing into four equal parts suces to show that q1 is unreachable.
Asymptotic completeness A hybrid automaton can be approximated arbitrarily closely by choosing a suciently ne partitioning of the invariants. Given a real 0, the -extension of a predicate over Rn is the predicate 0 such that [ 0] = fs0 j 9s 2 [ ] such that dist (s; s0 ) g, where dist (s; s0 ) is the maximal value of the absolute dierence between corresponding coordinates of s and s0 . For a hybrid automaton A, let A be the hybrid automaton obtained from A where all invariants, initial conditions, ow conditions, and jump conditions are replaced by their -extensions. Intuitively, A models a system that diers from A by at most in its measurements of the values for continuous variables. A partitioning P is of width if, for each invariant occurring in P and for all valuations s1 ; s2 2 [ ] , dist (s1 ; s2 ) . A hybrid automaton is monotonic if, in each location, each variable is either nondecreasing or nonincreasing. A hybrid automaton is bounded if for all locations v, inv (v) is a bounded set (this requirement can be relaxed without aecting the following two results nor Theorem 5; see [HH95a]). Rate approximations, and therefore linear phase-portrait approximations, are asymptotically complete for monotonic bounded automata. 9
Theorem3 [HH95a]. Let A be a monotonic bounded hybrid automaton. For all reals > 0, (1) there exists a rate approximation B with some partitioning P such that P (A) B P (A ), and (2) there exists a real > 0 such that for all partitionings P of width , if rate P (A) exists, then P (A) rate P (A) P (A ). Corollary4. Let A be a monotonic bounded hybrid automaton. For all reals > 0, (1) there exists a linear phase-portrait approximation B with some partitioning P such that P (A) B P (A ), and (2) there exists a real > 0 such that for all partitionings P of width , if phase P (A) exists, then P (A) phase P (A) P (A ). While splitting locations and nding rate approximations can be used for arbitrary accuracy, it may be easier to obtain suciently accurate linear phaseportrait approximations with fewer locations, as suggested above for the example automaton A.
5 Example: Predator-Prey Systems We demonstrate the use of linear phase-portrait approximations on nonlinear systems modeling the population growth of two interacting species [Lot20]. We show that some interesting properties of such systems can be discovered automatically through algorithmic analysis. A predator-prey ecology with limited growth Much of our exposition de ning predator-prey systems is derived from Chapter 12 of [HS74]. One species is the predator, whose population is modeled by the variable y, and the other species is the prey, modeled by x. The prey forms the entire food supply for the predator, and we assume that the per capita food supply for the predator at any instant of time is proportional to the number of prey. The growth of the predator population is proportional to the dierence between its actual per capita food supply and a basic per capita food supply required to maintain its population. The population of the prey is subject to two competing forces. First, the population may grow because there is a constant food supply available: the prey's population would increase without bound in the absence of predators. Furthermore, we assume this rate of increase would be proportional to the number of prey. Second, the predators consume the prey at a rate that is proportional to the number of predators and to the number of prey. No population has the potential to increase without bound. There are social phenomena, such as overcrowding, spread of disease, and pollution, that imply that populations experience negative growth once they exceed a threshold limiting population. Assuming these negative growth factors are proportional to the species population and its difference from the threshold population leads to the Volterra-Lotka predator-prey equations. x_ = (A ? By ? x)x and y_ = (Cx ? D ? y)y, where A, B , C , D, , and are all positive real-valued constants. The hybrid automaton for the system appears in Figure 3. By examining the
ow eld determined by the ow condition above, we partition the state space 10
y
x_ 0 y_ 0
L x_ = 0
M y_ = 0
y R
0 x_ 0 xy__ y_ 0 0
x_ 0 y_ 0
Fig. 4.
x=0^ y A=B x_ = 0 ^ y_ 0
.
S3
(x0 ; y0 )
A= D=C
x
(b) Reachability using phase-portrait approx.
Predator-prey populations
x=0 L x_ 0 ^ y_ 0
L ^ M x_ 0 ^ y_ 0
L
y = A=B x=0 x=0^ y=0 y=0 y A=B x_ = 0 y = 0 ^ y=0^ ^ y_ 0 y = 0 x A= x = A= x A= ^ x D=C y_ = 0 y_ = 0 ^ x_ 0 ^ x_ 0 Fig. 5.
1
S2
A=
x D=C (a) Phase portrait: L, M non-intersecting
M y_ = 0 R y_ ?C x=B _ S y_ 0
x_ 0 y_ 0
L x_ = 0
M
x = D=C
M y_ ?C x=B _ ^ y_ 0 y=0^ x D=C y_ = 0 ^ x_ 0
Linear phase-portrait approximation for the predator-prey system
along the lines where either x_ or y_ have value 0, i.e. along the coordinate axes, and along the lines L : A ? By ? x = 0 and M : Cx ? D ? y = 0. For the case that the lines L and M do not intersect in the northeastern quadrant of R2, the phase portrait of the system is shown in Figure 4, part (a). Linear phase-portrait approximation In the region R, to the right of the line M in Figure 4, we infer tighter constraints on x_ and y_ than their signs. The values taken by the function f (x; y) = y= _ x_ in the region R determine the ow vectors in R, because x_ is nonpositive and y_ nonnegative. The absolute value of f (x; y) is bounded above by any Max =Min , where Max is an upper bound on the value of y_ in R, and Min is a lower bound on the absolute value of x_ in R. We can take Cxy for Max , because D + y is always positive. Since the lines L and M do not intersect in the northeastern quadrant, we know that A= < D=C , and hence that A ? D=C < 0. Since x is no less than D=C in R, we infer that A ? x < 0, and hence that A ? x ? By < ?By. We may therefore take Byx for Min . We conclude that f (x; y) is bounded below by ?Cxy=(Bxy) = ?C=B , and thus that all ow vectors in R have directions between (?B; C ) and (?1; 0), i.e. they all satisfy the ow condition y_ 0 ^ y_ ?C x=B _ . The automaton for a linear phase-portrait approximation appears in Figure 5. The layout of the locations matches the partitioning of the state space as 11
shown in Figure 4. The implicit invariant constraint x 0 ^ y 0 has been omitted from all invariants. The constraint M refers to all valuations on the line M , i.e. all valuations where Cx ? D ? y = 0. The constraint M refers to all valuations at, or to the right of, the line M , i.e. M is Cx D + y. Similarly M is Cx D + y, L is A ? By ? x = 0, L is x A ? By, and L is x A ? By. Computing population bounds The phase-portrait approximation above can be used to compute, for given starting populations, bounds on the populations of both species. For simplicity, we consider only the case where the initial populations x0 and y0 lie in the region R of the state space. The set of time-step successors of the state (x0 ; y0) is obtained by following all ow vectors in the cone indicated in Figure 4, part (b). First, the states in region S1 are reached. Control may then pass to the location corresponding to the central region in the partition, where both x_ and y_ are nonpositive. After adding the states in region S2 , and then S3 , reachability analysis terminates. The maximum value of y among the reachable states is (By0 + Cx0 ? D)=(B + ). For example, using the equations x_ = (2000 ? y ? 5x)x and y_ = (4x ? 2600 ? 4y)y, and the initial population vector (900; 150), HyTech automatically computes a bound of 230 for y. The given phase-portrait approximation yields a strictly tighter bound than the rate translation when applied to the same partitioning. The rate translation for each of the quadrants above provides no further information than the sign of the derivatives x_ and y_ . Therefore, given initial starting populations x0 and y0 in R, the upper bound we obtain for y is (Cx0 ? D)= by following the initial point vertically up to the line M , as opposed to following the diagonal to M for the linear phase-portrait approximation. Naturally, either approximation method can yield tighter bounds at the expense of a ner partitioning. (If the lines L and M do intersect in the northeastern quadrant, reachability analysis on a linear phase-portrait approximation can be used to obtain bounds on both populations in this case, too.) Iterating approximations As suggested in Section 3, bounds on the reachable region can be used to generate better phase-portrait approximations. For example, it can be shown that in the rightmost region R of Figure 4, the absolute value of f (x; y), the ow tangent at (x; y), is bounded above by Cy=(x + By ? A). Let Z be a bounded subset of R. Let ymax be an upper bound for y over all valuations in Z , and xmin (resp. ymin ) be a lower bound for x (resp. y) over Z . It follows that jf (x; y)j Cymax =(xmin + Bymin ? A), provided that (xmin + Bymin ? A) 0. Previously, we showed how reachability analysis of the automaton in Figure 5 leads to the region S1 in R, from which we infer bounds of ymax = 230, ymin = 150, and xmin = 800. We can therefore replace the ow condition y_ ?C x=B _ in the region S1 with y_ ?92x= _ 215. Recomputation now shows that only a subset of the region S1 is reachable. In particular, we obtain the tighter bound of ymax = 55250=307 180 for ymax . We can iterate this process, obtaining successively lower values of ymax , and more restrictive ow conditions. 12
Controlling the ecology Standard analysis techniques can be used to show that in our setup, the predator population always tends toward 0, while the prey population tends to A=. Suppose, however, that we wish to keep the predator population above a nontrivial minimal value, or more generally, that the populations need to be controlled so that they remain within given lower and upper bounds. Assume that the prey population can be measured accurately, but that the predator population is unobservable. Our control strategy consists of monitoring the prey population, and releasing a xed number k of additional prey into the system whenever it reaches its minimal allowable value. In general, it is unwise to increase the prey population to its maximal allowable value, because the abundance of prey may cause the predator population to grow too large. In our example, we require the predator population to lie within the range [100; 350], and the prey population within [800; 1100]. HyTech con rms automatically that both population bounds are successfully maintained when k 200. For larger values of k, our linear phase-portrait approximation admits trajectories where the predator population exceeds the upper bound of 350. (This does not imply that all values of k greater than 200 lead to excessively large predator populations, because the approximation yields more reachable states than the original system.)
6 Example: Pseudo-linear Hybrid Automata It is nontrivial to construct useful linear phase-portrait approximations for arbitrary hybrid automata, because this may involve understanding a set of dierential inequalities. Here we identify a restricted class of hybrid automata for which linear phase-portrait approximations may be constructed automatically. A hybrid automaton with the set X of variables is pseudo-linear if (1) for every ow condition, replacing all occurrences of undotted variables with arbitrary rational constants results in a convex linear predicate over the dotted variables in X_ , and (2) for every location v, every operator (function symbol) in the ow condition
(v) is monotonic in each argument over the domain inv (v). For instance, the automaton A of Figure 1 is pseudo-linear. Given a pseudo-linear hybrid automaton A, the following transformation yields a linear phase-portrait approximation (A). For each location v, every occurrence of a variable x in the ow condition (v) is replaced by either xmax or xmin , depending on the force of the occurrence of x in (v), where xmax is a rational or in nite upper bound on the values that x may take along any trajectory satisfying the invariant inv (v) and the ow condition (v), and xmin is a rational or in nite lower bound. The replacement of each variable occurrence by a constant is chosen so that the transformed ow condition is implied by the conjunction of the invariant with the original ow condition. For simplicity, we assume that (v) is a nite conjunction of inequalities. Then, an occurrence of x at the right-hand (left-hand) side of an inequality has positive force, and is replaced with xmax , if the occurrence lies within the scope of an even (resp. odd) 13
number of operators that are monotonically nonincreasing in the occurrence of x over inv (v). Otherwise, the occurrence of x has negative force and is replaced with xmin. Theorem5. Given a pseudo-linear hybrid automaton A, the automaton (A) is a linear phase-portrait approximation of A (with the trivial partitioning). Example 4. The automaton A of Figure 1 is pseudo-linear. The transformation above leads to the linear hybrid automaton (A) with the ow condition y_ x_ ^ y_ 3x_ ^ x_ 3 ^ x_ 1 for q0. The rst two conjuncts are derived from y_ = xx_ , and the latter two from x_ = x. In this case, the derived ow conditions exactly match those of the linear phase-portrait translation with the trivial partitioning. Example 5. The predator-prey system of the previous section is also pseudolinear. However, the phase-portrait approximation we derive by the transformation above is the rate translation. This demonstrates the limitation of automatic methods: often a more careful analysis can provide tighter approximations. Acknowledgement We thank Pei-Hsin Ho for interesting discussions about approximating nonlinear hybrid systems and the cumulative errors involved.
References [ACH+95] R. Alur, C. Courcoubetis, N. Halbwachs, T.A. Henzinger, P.-H. Ho, X. Nicollin, A. Olivero, J. Sifakis, and S. Yovine. The algorithmic analysis of hybrid systems. Theoretical Computer Science, 138:3{34, 1995. [ACHH93] R. Alur, C. Courcoubetis, T.A. Henzinger, and P.-H. Ho. Hybrid automata: an algorithmic approach to the speci cation and veri cation of hybrid systems. Hybrid Systems I, LNCS 736, pp. 209{229. Springer-Verlag, 1993. [AHH93] R. Alur, T.A. Henzinger, and P.-H. Ho. Automatic symbolic veri cation of embedded systems. 14th Annual Real-time Systems Symposium, pp. 2{11. IEEE Computer Society Press, 1993. [Arn83] V. I. Arnol'd. Geometric Methods in the Theory of Ordinary Dierential Equations. Springer-Verlag, New York, 1983. [CC77] P. Cousot and R. Cousot. Abstract interpretation: a uni ed lattice model for the static analysis of programs by construction or approximation of xpoints. 4th Annual Symposium on Principles of Programming Languages. ACM Press, 1977. [Hen95] T.A. Henzinger. Hybrid automata with nite bisimulations. ICALP 95: Automata, Languages, and Programming, LNCS 944, pp. 324{335. Springer-Verlag, 1995. [HH95a] T.A. Henzinger and P.-H. Ho. Algorithmic analysis of nonlinear hybrid systems. CAV 95: Computer-aided Veri cation, LNCS 939, pp. 225{238. SpringerVerlag, 1995. [HH95b] T.A. Henzinger and P.-H. Ho. HyTech: The Cornell Hybrid Technology Tool. Hybrid Systems II, LNCS 999, pp. 265{293. Springer-Verlag, 1995. [HH95c] T.A. Henzinger and P.-H. Ho. A note on abstract-interpretation strategies for hybrid automata. Hybrid Systems II, LNCS 999, pp. 252{264. Springer-Verlag, 1995.
14
[HHWT95a] T.A. Henzinger, P.-H. Ho, and H. Wong-Toi. HyTech: the next generation. 16th Annual Real-time Systems Symposium, pp. 56{65. IEEE Computer Society Press, 1995. [HHWT95b] T.A. Henzinger, P.-H. Ho, and H. Wong-Toi. A user guide to HyTech. TACAS 95: Tools and Algorithms for the Construction and Analysis of Systems, LNCS 1019, pp. 41{71. Springer-Verlag, 1995. Full version available as Technical Report CSD-TR-95-1532, Cornell University (write to
[email protected], or check http://www.cs.cornell.edu/Info/People/tah/hytech.html). [HKPV95] T.A. Henzinger, P.W. Kopke, A. Puri, and P. Varaiya. What's decidable about hybrid automata? 27th Annual Symposium on Theory of Computing, pp. 373{ 382. ACM Press, 1995. [HRP94] N. Halbwachs, P. Raymond, and Y.-E. Proy. Veri cation of linear hybrid systems by means of convex approximation. SAS 94: Static Analysis Symposium, LNCS 864, pp. 223{237. Springer-Verlag, 1994. [HS74] M.W. Hirsch and S. Smale. Dierential Equations, Dynamical Systems, and Linear Algebra. Academic Press, 1974. [HW95] P.-H. Ho and H. Wong-Toi. Automated analysis of an audio control protocol. CAV 95: Computer-aided Veri cation,LNCS 939, pp. 381{394. Springer-Verlag, 1995. [Lot20] A.J. Lotka. Analytical note on certain rhythmic relations in organic systems. Proceedings of the National Academy of Sciences of the United States of America, 6:410{415, 1920.
15