Nonlinear Analysis 65 (2006) 1150–1167 www.elsevier.com/locate/na
Hybrid abstractions of affine systems Marie-Anne Lefebvre, Herv´e Gu´eguen ∗ Sup´elec - IETR, avenue de la Boulaie, BP 81127, F35511 Cesson-S´evign´e Cedex, France
Abstract This paper considers the problem of building a set of hybrid abstractions for affine systems in order to compute over approximations of the reachable space. Each abstraction is based on a decomposition of the continuous state space that is defined by hyperplanes generated by linear combinations of two vectors. The choice of these vectors is based on consideration of the dynamics of the system and uses, for example, the left eigenvectors of the matrix that defines these dynamics. We show that the reachability calculus can then be performed on a combination of such abstractions and how its accuracy depends on the choice of hyperplanes that define the decomposition. c 2005 Elsevier Ltd. All rights reserved. Keywords: Hybrid systems; Reachability; Abstractions
1. Introduction Recent advances in control and the increasing need for integration and validation of control applications have led to the development of research of automated verification for continuous and hybrid systems. Safety properties are especially important and are generally expressed as regions of the state space that must not be reached by the system. Verification may then be based on logical characterizations of the reachable space [1,2]. Other approaches aim at abstracting, from the existing initial model, a discrete event model on which discrete verification can be performed [3,4]. Finally, some approaches are based on an explicit reachability calculus in the hybrid state space [5,6]. The last two approaches require computation of the reachable continuous state space. This task is complex when the continuous dynamics are nontrivial, and remains one of the main obstacles in safety verification of hybrid systems [7]. Most solutions that have been ∗ Corresponding author. Tel.: +33 2 99 84 45 45 04; fax: +33 2 99 84 45 45 99.
E-mail address:
[email protected] (H. Gu´eguen). c 2005 Elsevier Ltd. All rights reserved. 0362-546X/$ - see front matter doi:10.1016/j.na.2005.12.016
M.-A. Lefebvre, H. Gu´eguen / Nonlinear Analysis 65 (2006) 1150–1167
1151
proposed are based on the computation of an over-approximation of the reachable space, as this is sufficient to decide whether the property is satisfied. To compute the reachable space, it is possible to explicitly calculate an over-approximation of the continuous reachable space from a time-sampled calculation, expanded to include all the trajectories. The various propositions then differ in the management of these over-approximations and in the representation (polyhedral, ellipsoidal, etc.) of the manipulated regions (see for example [8–10]). Other approaches, called ‘hybridization-based’ methods in [11], aim at developing a hybrid automaton with simpler dynamics from the continuous system [6]. Asarin, Dang and Girard [12], for example, proposed the use of linear approximations with bounded input, but generally the goal is to produce a linear hybrid automaton with differential inclusions as an abstract model to be able to compute reachability, for example, with a tool such as HyTech [13]. In these ‘hybridization-based’ methods, to approximate conservatively the continuous dynamics, the continuous state space is decomposed using linear constraints and a location is created in the abstract model for each cell. A polyhedron that includes the vector field for all points of the cell is then computed and it defines the differential inclusion associated to that location. Finally, a transition from one location to another location is introduced when their relative regions in the state space have a common boundary and there exists at least one continuous trajectory that crosses this frontier from the first to the second region. Each trajectory of the concrete system can be mapped to a trajectory of the hybrid automaton, thereby defining an abstraction. Naturally, some trajectories of this later system, called spurious trajectories, are not real trajectories of the system. The reachable space of the abstract model is then an over-approximation of the reachable space of the concrete system. The problem with these approaches is then to find a trade-off between the accuracy of the abstraction and its simplicity, i.e., chiefly the number of locations that are introduced. In a previous paper [14], we proposed a method to guide the decomposition of the state space to abstract affine planar regular systems. This method relies on the properties of the continuous dynamics to obtain an abstraction that is relatively simple to analyze but does not introduce too many spurious trajectories. This paper proposes a generalization of this approach to higher-dimension systems that may also be singular. Therefore, it presents a method of computing an abstraction of an affine continuous system by a linear hybrid automaton [6], i.e., which is defined by linear predicates over the state variables and their derivatives. This method differs from other works [6,15] by an intelligent choice of the hyperplanes used to decompose the state space. In Section 2, the approach for planar systems is presented to illustrate the basic concepts of the approach. The properties of affine systems that are used to guide the abstraction are then presented in Section 3, and the global abstraction and reachability procedure is detailed in Section 4. Finally, Section 5 considers the question of accuracy of the abstraction and gives some heuristics to produce more accurate approximations. 2. Planar regular systems In this section, the system under study is a planar affine system whose dynamics are defined by Eq. (1), where the matrix A is regular so that there exists an equilibrium point denoted xe (2). x˙ = Ax + b
(1)
−1
(2)
xe = −A
b.
1152
M.-A. Lefebvre, H. Gu´eguen / Nonlinear Analysis 65 (2006) 1150–1167
Fig. 1. Cells and associated differential inclusions for planar systems.
The problem is then to define a set of lines in order to determine a decomposition of the state space leading to an abstract model with properties that simplify the reachability computation. Dynamics properties of the system described by (1) can be used to select these lines because it is possible to prove (see [14] for details) that, for all points of a line passing through the equilibrium point, the derivative vectors are collinear. For a line defined by Eq. (3), it is then possible to characterize the vector field by Eq. (4). qT (x − xe ) = 0
(3)
γT x˙ = 0
(4)
−1 q. with γ = AT
Therefore, for all points of a cell defined by (5) and delimited by two such lines, the vector field belongs to the region defined by (6) whose boundaries are collinear with the two derivative vectors obtained on each boundary of the cell, as illustrated in Fig. 1. T qiT (x − xe ) ≥ 0 ∧ qi+1 (5) (x − xe ) ≥ 0 T x˙ ≤ 0 . (6) γiT x˙ ≥ 0 ∧ γi+1 Moreover, if the vector used to specify the boundary (3) is a real left eigenvector of the matrix A, this boundary cannot be crossed by any continuous trajectory; otherwise it is always crossed in the same direction. Any set of vectors {qi } then specifies an automaton such that each location is associated with an invariant given by the definition of the cell (5) and the continuous dynamics specified by the inclusion (6) where the two vectors γi and γi+1 are given by (4). As each boundary is only crossed in one direction, each location is the source and the target of at most one transition. This procedure is illustrated in Fig. 2 for a system with complex eigenvalues. From the decomposition of the state space (Fig. 2(a)) and the direction of the vector field on the boundaries (symbolized by arrows), the automaton of Fig. 2(b) is realized. One situation si is associated to the cell Si , and its associated differential inclusion, given by (6) and symbolized by the couples of arrows, is specified by the direction of the vector field on the boundaries. The guard of a transition from a situation to its successor is defined by the equation of the common boundary of the two adjacent cells. Note that (6) imposes constraints only on the direction and not on the norm of the vector field. This does not pose any problem, as the abstraction is used to calculate the reachable space without temporal information. The reachable space from one point x0 of a cell defined by (5)
M.-A. Lefebvre, H. Gu´eguen / Nonlinear Analysis 65 (2006) 1150–1167
1153
Fig. 2. Example of decomposition of the state space (a) and the associated automaton (b).
Fig. 3. Reachable space from one point (a), a line segment (b), and propagation within two adjacent cells (c).
within this cell can be simply computed and it is specified by the conjunction of (5) and (7), as illustrated in Fig. 3(a). T γiT (x − x0 ) ≥ 0 ∧ γi+1 (7) (x − x0 ) ≤ 0 . This result is proved by considering the evolution of the state vector from an initial value given by: t t x˙ (τ ) dτ = x0 + A (x (τ ) − xe ) dτ. x (t) = x0 + t0
t0
When the state vector is within the cell defined by (5), Eq. (6) can be used to write: t t γiT (x (t) − x0 ) = γiT A (x (τ ) − xe ) dτ = γiT A (x (τ ) − xe ) dτ =
t0
t t0
and
T γi+1
(x (t) − x0 ) =
t0
qiT (x (τ ) − xe ) dτ ≥ 0
t t0
T qi+1 (x (τ ) − xe ) dτ ≤ 0.
1154
M.-A. Lefebvre, H. Gu´eguen / Nonlinear Analysis 65 (2006) 1150–1167
Fig. 4. Reachable space for example ‘saddle’.
Because all the regions of the state space considered are convex, it is possible to compute the reachable space, from a convex polytope within a cell, as the convex hull of the reachable spaces from its vertices, as illustrated in Fig. 3(b). From the computation of the intersection of this reachable space with the boundary of the cell, specified by its vertices, it is possible to obtain the reachable region in the next cell (Fig. 3(c)) and to iteratively compute the global reachable space. The application of this method to planar affine systems is illustrated by two examples taken from Dang and Maler [9]. For the first example, named ‘saddle’, the dynamics and initial region are defined by (8) and the reachable space is represented in Fig. 4, whereas the second, named ‘sink’, is defined by (9) and illustrated in Fig. 5. In both cases the results are very similar to those obtained with existing tools, such as d/dt [16] or Coho [17]. 0 ≤ x 1 ≤ 0.4 −5 0 0 (8) A= b= X0 : 0 ≤ x 2 ≤ 0.4 0 4 0 −2 −3 0.1 ≤ x 1 ≤ 0.3 0 A= (9) b= X0 : 0.1 ≤ x 2 ≤ 0.3. 3 −2 0 3. General affine systems To abstract an affine continuous system by a linear hybrid automaton it is necessary to decompose the state space and to associate a differential inclusion to each region. The study of planar regular systems has shown that some lines have particular properties to define this decomposition and that the resulting abstract model then has ‘good’ properties for reachability calculus. The aim is then to generalize this to general affine dynamics specified by Eq. (10), without any assumptions on the properties of matrix A, and to find hyperplanes to split the state space, so that the abstract hybrid systems might have the same ‘good’ properties for reachability calculus. When defining the hyperplanes the objectives are twofold: – for all points of a plane, the vector field must be normal to the same vector, – whenever possible, all the continuous trajectories must cross the boundary of a region in the same direction.
M.-A. Lefebvre, H. Gu´eguen / Nonlinear Analysis 65 (2006) 1150–1167
1155
Fig. 5. Reachable space for example ‘sink’.
The first property ensures that the differential inclusion can be easily computed and the second property ensures that each boundary will be associated with at most one transition. x˙ = Ax + b.
(10)
3.1. First criterion The first criterion that a hyperplane H (q, k) defined by Eq. (11) must fulfill can be formalized by Eq. (12). (11) H (q, k) = x/qT x = k ∃γ s.a ∀x ∈ H (q, k) ,
γT x˙ = 0.
(12)
Property 1: The conditions that must hold on q, γ and k so that Eq. (12) is satisfied are given by Eq. (13). q = AT γ k = −γT b.
(13)
This can be proved by considering that if x1 ∈ H (q, k) and x2 ∈ H (q, k) then qT (x1 − x2 ) = 0, and if γ is such as γT x˙ 1 = 0, then: γT x˙ 2 = γT (Ax2 + b) = γT (Ax1 + b + A (x2 − x1 )) = γT A (x2 − x1 ) . Therefore, if γT x˙ 2 = 0, q and AT γ are collinear, and as the direction of these vectors is the only important property, the first component of Eq. (13) holds without loss of generality. Then,
1156
M.-A. Lefebvre, H. Gu´eguen / Nonlinear Analysis 65 (2006) 1150–1167
γT x˙ 1 = γT (Ax1 + b) = q T x1 + γ T b = k + γT b which proves the second component of (13). The choice of one vector q, in the image of AT , then implies some constraints on the constant k to satisfy Eq. (13), to obtain a relevant hyperplane. This leads to either one or an infinite number of values for k, according to the characteristics of matrix A: – If matrix A is invertible, from (13) we deduce: γ = A−T q, and therefore k = −γT b = −qT A−1 b. That is, the value of k is unique. – If matrix A is singular, then A has some eigenvectors, denoted w0i , associated with the null eigenvalue. The general solution of the first equation of (13) is then:
αi w0i (14) γ = γ0 + where γ0 is a particular solution of (13), and therefore the constant k may be written as:
T k = −γ0T b − αi w0i b. (15) T b = 0, then k is fixed to a unique value, otherwise k is free. If for all w0i , w0i This property is of significance when considering hyperplanes, orthogonal to one vector q, verifying Eq. (12), and may be summarized as follows:
Property 2.a: If q does not belong to the image of AT then no such hyperplane exists. Property 2.b: If the matrix A is either regular or such as all eigenvectors, denoted w0i , associated T b = 0, then a vector q in the image of AT defines only one with a null eigenvalue satisfy w0i hyperplane to satisfy Eq. (13). Property 2.c: If the matrix A is singular and at least one of the eigenvectors w0i associated with T b = 0, then a vector q in the range of AT defines an infinite the null eigenvalue verifies w0i number of hyperplanes to satisfy Eq. (13). 3.2. Differential inclusion and reachable space This first criterion is valuable in computing the differential inclusions and the reachable space for the abstract model, as the following properties are satisfied. Property 3: For all points of a region defined by (16), where the boundaries are hyperplanes, H (qi , ki ), such that (13) is satisfied, the predicate (17) is true. T qT (16) 1 x ≥ k 1 ∧ q2 x ≤ k 2 ˙ ≥ 0 ∧ γT ˙≤0 . (17) γT 1x 2x Property 3 is directly proved by γiT x˙ = γiT (Ax + b) = qiT x − ki ≺i 0, where ≺i denotes ≥ when i = 1 and denotes ≤ when i = 2.
M.-A. Lefebvre, H. Gu´eguen / Nonlinear Analysis 65 (2006) 1150–1167
1157
Property 4: An over-approximation of the reachable space from a point x0 , within a region defined by (16), is then given by the conjunction of (16) and (18). T γT (18) 1 (x − x0 ) ≥ 0 ∧ γ2 (x − x0 ) ≤ 0 . This relation generalizes the result obtained previously for planar systems in the regular case (Eq. (7)) and results from: t t γiT x (t) = γiT x0 + x˙ (τ ) dτ = γiT x0 + γiT Ax (τ ) + γiT b dτ t0
t T = γi x0 + qiT x (τ ) − ki dτ. t0
That is, using (16), γiT (x(t) − x0 ) = previously for ≺i .
t
T t0 (qi x(τ )
t0
− ki )dτ ≺i 0 with the same convention as
3.3. Second criterion The second criterion used to choose the hyperplanes H (q, k) that define the decomposition is the uniqueness of the direction in which continuous trajectories cross the boundaries of a cell. This direction changes according to another hyperplane H S (q) [18] defined by: qT (Ax + b) = 0. In general, nothing can be said about this hyperplane, but some special cases, where the vector q is generated by left eigenvectors of A, are worth consideration. The use of left eigenvectors allows, first, a formal expression for the vector γ associated to a hyperplane H (q, k) to be derived, and, second, a study of how this hyperplane is crossed by continuous trajectories. It is closely akin to the principles used to characterize reachable space by [1] or [19], but is used here to meet the requirement of producing the most deterministic abstraction, which is very close to that studied in [20] using another approach. The first interesting case occurs for the hyperplanes defined by a left eigenvector w of A associated with a real eigenvalue (λ = 0). Eqs. (14) and (15) then give:
1 αi w0i γ= w+ λ (19)
1 T k = k0 − αi w0i b with k0 = − wT b. λ The hyperplane H (w, k0 ) is such as: wT x˙ = wT (Ax + b) = λwT x + wT b 1 1 = λ wT x + wT b = λ k 0 + wT b λ λ = 0. This hyperplane is a separating hyperplane that cannot be crossed by the continuous trajectories. Under the assumptions of Property 2.b, this hyperplane is also the only hyperplane defined by w to be useful in relation to criterion 1. Under the assumptions of Property 2.c, k is free and w defines a family of hyperplanes H (w, k). All the hyperplanes located on the
1158
M.-A. Lefebvre, H. Gu´eguen / Nonlinear Analysis 65 (2006) 1150–1167
same side of the separating plane are crossed in the same direction as wT x˙ = λk + wT b is constant. The hyperplane H (q, k) with q = α1 w1 + α2 w2 , where w1 and w2 are left eigenvectors of A associated with real eigenvalues (λ1 = 0, λ2 = 0), is also an interesting case. Eqs. (14) and (15) then give:
α2 α1 w1 + w2 + α0i w0i γ= λ1 λ
2 k = k0 − α0i w0i T b (20) α1 α2 with k0 = − w1T b − w2T b. λ1 λ2 The condition on a change of the direction of crossing is: qT x˙ = α1 λ1 w1T x + α2 λ2 w2T x + α1 w1T b + α2 w2T b α1 α2 = λ1 α1 w1T x + α2 w2T x + w1T b + w2T b λ1 λ2 λ 1 − α2 λ1 w2T x + α2 λ2 w2T x − α2 w2T b + α2 w2T b λ2 1 T T = λ1 (k − k0 ) + α2 (λ2 − λ1 ) w2 x + w2 b . λ2 As above, if Property 2.b is verified then q defines only one hyperplane H (q, k0 ) according to Eq. (20) and the sign of the expression qT x˙ depends on the position of x with respect to the separating hyperplane defined by w2 . Therefore, it is constant when the decomposition of the half space defined by this hyperplane is considered. Conversely, if the conditions of Property 2.c are satisfied, then k is free and criterion 1 specifies a family of hyperplanes for each (α1 , α2 ). However, consistent with the second criterion, the only member of this family that is interesting is defined by α0i = 0. With equivalent considerations, not detailed here but which can be found in [21], this result is still valid if q is defined as a real combination of two conjugate complex left eigenvectors. Note that these results are consistent with the previous results for regular planar systems, because the hyperplanes that satisfy (13) are then the lines to which the equilibrium point belongs. Finally, when q is a linear combination of a real left eigenvector with a combination of complex eigenvectors, the plane H S that defines the change of crossing direction is unique and only depends on the eigenvectors considered. Therefore, it has to be computed once only for all values of the coefficients of combination [21]. Note that when the assumption of Property 2.c is satisfied, all the hyperplanes H (q, k) (with k varying) have the same hyperplane H S for the change of direction of trajectories. Any other choice of the vector q requires consideration of different hyperplanes H S that have to be computed for each H (q, k). 4. Abstraction and reachability The global procedure to abstract an affine continuous system and to compute an overapproximation of its reachable space is summarized in Algorithm 1. Each step of this procedure is detailed below. The basic concept of the proposed approach is that it is easier to consider a set of decompositions that defines a set of abstractions and then to compose them using an automata product than to build a global abstraction.
M.-A. Lefebvre, H. Gu´eguen / Nonlinear Analysis 65 (2006) 1150–1167
1159
Algorithm 1 (Global Procedure). Input: (A, b) continuous dynamics, Inv considered region of the state space, X 0 initial region, Result: over-approximation of Rc (X 0 ), the region of Inv reachable from X 0 with respect to the dynamics constraint defined by (A, b) step 1: determination of basic vectors step 2: choice of the set of generating elements for decomposition step 3: computation of the abstract automaton for each generating element step 4: reachability computation for the composition of all abstract automata. 4.1. Determination of basic vectors This first step aims at determining a set of linearly independent vectors to be used in the next steps to generate the decomposition of the state space. On the one hand, this choice is constrained by the first part of Eq. (13), that is to say that, the vectors must be in the range of AT , and on the other hand, it is guided by the consideration of Section 3.3. This step is summarized in Algorithm 2. Algorithm 2 (Choice of Basic Vectors). (A, b) continuous dynamics basis for the range of AT {wi } with the associated eigenvalues. computation of the left eigenvectors and generalized vectors of A computation of independent real combinations of complex conjugate vectors selection of real left eigenvectors and generalized vectors associated with non-null eigenvalue and real combinations computed at step 2 step 4: selection of the vectors associated with the null eigenvalue.
Input: Result: step 1: step 2: step 3:
The left eigenvectors and the generalized vectors that are computed in the first step are deduced from the Jordan decomposition of the matrix A and the recurrent relation between these vectors denoted ui,m , and an eigenvector ui,1 associated to an eigenvalue λi is given by: T T A λi ui,1 = ui,1 T T T λi ui,m + ui,m−1 = ui,m A
m = 2..n i .
These vectors are linearly independent and when their associated eigenvalue is not null they are all in the range of AT . When the eigenvalue is null, the first n i − 1 vectors are in the range of AT but not the last. This is why the two cases are distinguished in the algorithm in steps 3 and 4. 4.2. Choice of the set of generating elements A generating element will be used in the next step to generate a family of hyperplanes that define one abstraction of the continuous dynamics. A generating element can be a single vector when the conditions of Property 2.c are satisfied or a pair of vectors otherwise. These vectors are selected in the set {wi } computed in the previous step. From consideration of Section 3.3, the single-vector generating elements are chosen first from the real left eigenvectors, referred to below as category G 1 , and second from the other vectors of the set referred to as category G 2 . The other generating elements are: – first, pairs of real left eigenvectors, referred to as category G 3 – second, pairs of the two real vectors associated with complex conjugate eigenvectors, referred to as category G 4
1160
M.-A. Lefebvre, H. Gu´eguen / Nonlinear Analysis 65 (2006) 1150–1167
– third, pairs made up of one real eigenvector and one real combination of complex conjugate eigenvectors, referred to as category G 5 – last, pairs made up of other choices of two vectors in {wi }, referred to as category G 6 . Unfortunately, this step is not yet completely formalized and the choice of the vectors and the pairs that are selected to define generating elements is still based on prior knowledge and heuristics. For all hyperplanes considered below, Property 1 will be satisfied, therefore the category of the generating element has no effect on the differential inclusions but only on the transition structure of the abstract automaton. Therefore, to simplify computations in the following steps, categories G 1 , G 3 and G 4 are favored and categories G 5 and G 2 are preferred to G 6 . The number of generating elements that have to be selected is also an open question that is considered below. The only constraint on this number is that the regions that will be considered in the reachability computation have to be polytopes in order to be specified by their vertices. When the region of interest in the state space is bounded, one generating element is sufficient to fulfill this constraint. As a result of this step of the procedure, a set of generating elements {gi } is defined. Each element of this set belongs to one of the six specified categories and consists of one vector (categories G 1 and G 2 ) or a pair of vectors (categories G 3 , G 4 , G 5 and G 6 ), and all the vectors belonging to the set are computed in step 1. 4.3. Computation of abstract automaton Each generating element is used to generate a decomposition of the region of the state space under consideration and its associated abstract automaton. The procedure is linked to the category of the generating element and is summarized in Algorithm 3 for categories G 1 and G 2 . Algorithm 3 (Abstraction for Generating Elements of Categories G 1 and G 2 ). Input: gi = {w j } generating element and Inv the region of the state space considered. Result: HAi associated linear hybrid automaton. step 1: select a finite number of hyperplanes H (w j , kl ) specified by the value kl and compute the vector γil associated with each hyperplane from Eq. (13) step 2: if gi is in category G 2 , compute the hyperplane H S that specifies for all hyperplanes the change of crossing direction step 3: create a location for each cell that is the intersection of (wTj x ≥ kl ) ∧ (wTj x ≤ kl+1 ) and Inv. The differential inclusion of this location is given by (γTjl x˙ ≥ 0) ∧ (γTj(l+1)x˙ ≤ 0)
step 4: create the transitions associated with the crossing of hyperplanes H w j , kl , taking into account H S when required. The procedure for generating elements of categories G 3 , G 4 , G 5 and G 6 is summarized in Algorithm 4. Algorithm 4 (Abstraction for Generating Element of Categories G 3 , G 4 , G 5 and G 6 ). Input: gi = {wm , wn } generating element and Inv the region of the state space considered. Result: HAi associated linear hybrid automaton. step 1: compute a finite number of vectors ql = α1l wm + α2l wn , the associated hyperplane H (ql , kl ), and the associated vector γil according to Eq. (13)
M.-A. Lefebvre, H. Gu´eguen / Nonlinear Analysis 65 (2006) 1150–1167
1161
Fig. 6. Forms of abstract automaton.
step 2: if gi is in category G 5 , compute the hyperplane H S that specifies for all hyperplanes H (ql , kl ) the change of crossing direction; if gi is in category G 6 , compute for each hyperplane H (ql , kl ) the hyperplane H Sl that specifies the change of crossing direction T x≤k step 3: create a location for each cell that is the intersection of (qlT x ≥ kl )∧(ql+1 l+1 ) and T T Inv. The differential inclusion of this location is given by (γil x˙ ≥ 0) ∧ (γi(l+1) x˙ ≤ 0) step 4: create the transitions associated with the crossing of hyperplanes H (ql , kl ) taking into account H S or H Sl when needed. In these algorithms, an important consideration is the number and the values of the parameters k (Algorithm 3) and (α1 , α2 ) (Algorithm 4), and Section 5 gives some relationships between these choices and the accuracy of the over-approximation. Naturally, it is not of interest to consider values of these parameters if the cell, considered at step 3 of the two algorithms, is empty. Finally, when one of the vectors of the generating element is an eigenvector it is also important that the separating hyperplane that it specifies belongs to the set of hyperplanes so that the decomposition respects the invariant subspace of the continuous dynamics. Note that the basic principles of these algorithms are not new, but the choice of the hyperplanes that specify the decomposition simplifies their computation. In particular, step 1 of Algorithm 4 is simple to implement as it is mainly based on linear combinations because, if H (q1 , k1 ) and H (q2 , k2 ) are two hyperplanes that satisfy Property 1, then for all (α1 , α2 ), H (α1q1 + α2 q2 , α1 k1 + α2 k2 ) is also a hyperplane that satisfies this property. Furthermore, the vector γ = α1 γ1 +α2 γ2 (where γi is associated with H (qi , ki )) is the vector normal to the vector field. The six categories of generating elements lead to four general forms for the abstract automaton. Categories G 1 and G 3 are associated with automata as shown in Fig. 6(a) with a sink location and only one outgoing transition for each location. Category G 4 is associated with automata as shown in Fig. 2(b) with only one outgoing transition for each location. Categories G 2 and G 5 are associated with automata as shown in Fig. 6(b) where each location may be left by two transitions according to the two possible positions of the trajectories with respect to the unique hyperplane H S ; these are represented by ‘m’ and ‘p’ in the figure. Finally, category G 6 is associated with automata as shown in Fig. 6(c) where each location may be left by two transitions according to the positions of the trajectory with respect to the hyperplane H Si , which is specific to each hyperplane used to decompose the state space. Naturally, these are
1162
M.-A. Lefebvre, H. Gu´eguen / Nonlinear Analysis 65 (2006) 1150–1167
the general forms that may occasionally be simplified according to the region of the state space considered. 4.4. Reachability computation This step in the global procedure is a classic procedure (see for example [5]) to compute the reachable space of a set of linear hybrid automata and could be performed using, for example, HyTech [13]. However, the automata considered are specific linear hybrid automata as they are all the result of the abstraction of the same dynamics with the same principles. It is then possible to adapt a specific procedure to their characteristics. The global model is the product of the abstract automata HAi computed in the previous step. Each location of this global model is the composition of one location of each of the automata. Its invariant is then the intersection of the invariants of the composed locations and its differential inclusion is the intersection of their differential inclusions. The reachable space, from an initial point x0 , within this invariant is then specified by the conjunction of all the constraints (16) and (18) for the relevant locations of each basic automaton HAi . From the computation of the intersection of this reachable space, with the boundary of the invariant and with the guard of the transitions, it is possible to determine the new reachable locations of the global model, and their initial region, and then, the reachable space within their relative invariants. To avoid pointless computation the product automaton is not computed and the composition is performed for the locations that are reachable from the initial region. Naturally, the complexity of the computation is linked to the number of basic automata HAi that are combined. As stated in Section 4.2, the constraint on the number of automata is that the regions that are considered have to be characterized by their vertices in the reachability computation. It is then often useful to first undertake computations with this minimum number as they can lead to relevant conclusions. This point is illustrated in Fig. 7, which shows the reachable space for abstractions calculated with one generating element for a system whose dynamics are defined by (21), and for the initial domain and the invariant defined by (22). ⎛ ⎞ ⎛ ⎞ −1 2 0 0 0 ⎜−2 −1 0 ⎜0⎟ 0⎟ ⎜ ⎟ ⎜ A=⎝ b=⎝ ⎟ (21) 0 0 −3 0 ⎠ 0⎠ 0 0 0 −4 0 R0 = [4.2; 4.8] × [3.4; 4] × [4.2; 4.8] × [3.4; 4] Inv = [−10, 10] × [−10, 10] × [0, 10] × [0, 10] .
(22)
The generating element is first, g1 (category G 4 ), the pair of combinations of the two complex conjugate eigenvectors, that constrains the coordinates x 1 , x 2 and leads to the reachable space in Fig. 7(a); second, g2 (category G 3 ), composed of the two real eigenvectors, that constrains the coordinates x 3 , x 4 and leads to the reachable space in Fig. 7(b); and third, g3 (category G 5 ), composed of one eigenvector and one combination of complex eigenvectors, that leads to the reachable space in Fig. 7(c). With these three simple computations it is possible to conclude that some regions of the state space are unreachable even if the result of the composition of two automata is more accurate, as can be seen in Fig. 7(d) where the reachable space for the composition of automata based on g1 and g3 is drawn.
M.-A. Lefebvre, H. Gu´eguen / Nonlinear Analysis 65 (2006) 1150–1167
1163
Fig. 7. Examples of computed reachable spaces for abstraction based on: (a) g1 , (b) g2 , (c) g3 , (d) g1 and g3 .
5. Accuracy As stated in the introduction, the problem with these approaches of abstraction is to find a trade-off between the accuracy of the abstraction and its simplicity, i.e., mainly the number of locations that are introduced. The accuracy of an abstraction generated by a pair of given vectors is clearly linked to the number of hyperplanes (and locations): if more hyperplanes are considered the regions are smaller and the differential inclusions are closer, leading to a more accurate approximation but requiring more intersection computations. However, the position of the hyperplanes is also important. This point can be illustrated for the example defined by (23), where the region of interest is the cube [0, 4] × [0, 4] × [0, 4]. In this example, the generating elements are the pair of left eigenvectors of A that are the basis vectors of the space. For one of these generating elements, a set of vectors qi that define the decomposition hyperplanes and are generated by regular variation of the coefficients of linear combination in step 1 of Algorithm 4 is shown in Fig. 8(a). The relative vectors γi that define the differential inclusions are shown in Fig. 8(b). As can be seen, the regular disposition of the vectors in the state space corresponds to a very irregular disposition in the derivative space, and the importance of all regions, with respect to the final result, is not the same. ⎛ ⎞ ⎛ ⎞ −10 0 0 0 −5 0 ⎠ b = ⎝0⎠ X 0 = (4; 4; 4)T . (23) A=⎝ 0 0 0 −0.1 0 To qualify the precision of the abstraction, it is possible to consider the ‘width’ of the intersection of the reachable subregion from a point with the boundary of the region. This
1164
M.-A. Lefebvre, H. Gu´eguen / Nonlinear Analysis 65 (2006) 1150–1167
Fig. 8. Set of vectors defining the decomposition of the state space (a), and relative vectors defining the constraints for the differential inclusions (b).
Fig. 9. Accuracy of reachable space by abstraction.
precision is intuitively linked, for decomposition generated by combinations of two vectors (Fig. 9), to the spread of the normal vectors in the state space on the one hand and in the derivative space on the other. For example, for decompositions based on generating elements of categories G 3 , G 4 , G 5 and G 6 , the value defined by (24) can then be used when designing the family of hyperplanes in step 1 of Algorithm 4 to select a new hyperplane when it becomes smaller than a given threshold. qT q γT γ i i+1 · i i+1 . qi qi+1 γ γ i i+1
(24)
Results of the computation of the reachable space from one point, for the example defined by (23), are given in Fig. 10 and illustrate the usefulness of this criterion. The result in Fig. 10(a) is based on a regular decomposition, whereas the result in Fig. 10(b) is based on a decomposition with the same number of regions but where a hyperplane is added when the value of the measure defined by (24) becomes less than 0.85. It can be seen that the result is superior for the same computation load for the reachability calculus, and a little more for the abstraction. Similar results are obtained for the example in [16], the dynamics of which are defined by Eq. (25) and the initial domain by Eq. (26). The result of the reachability calculus based on a regular decomposition is given in Fig. 11(a), and the result based on decomposition with the same number of hyperplanes but where a hyperplane is added when the value of the measure becomes less than 0.95 is in Fig. 11(b). It can be seen that the second result is more accurate,
M.-A. Lefebvre, H. Gu´eguen / Nonlinear Analysis 65 (2006) 1150–1167
1165
Fig. 10. Influence of the position of the decomposition hyperplanes.
Fig. 11. Influence of the position of the decomposition hyperplanes.
especially in the divergence region. ⎛ ⎞ ⎛ ⎞ −1 −4 0 0 A = ⎝ 4 −1 0 ⎠ b = ⎝0⎠ 0 0 0.5 0 R0 = [0.025; 0.05] × [0.1; 0.15] × [0.05; 0.1] .
(25) (26)
6. Conclusion This paper presents an approach to abstract affine continuous systems using a set of hybrid automata to perform reachability calculus. Each abstraction is based on a decomposition of the state space by hyperplanes such that they also define hyperplanes for the vector fields to infer the differential inclusions of the abstraction. The choice of the vectors that define these hyperplanes is based on crossing considerations and favors left eigenvectors. This abstracting approach, which is based on the decomposition of the state space using dynamics properties, has proved to be more interesting than an a priori decomposition, as it naturally generates less spurious trajectories. Moreover, the calculus of the elementary abstractions and the reachable space are quite simple and straightforward. However, it may be pointless to decompose a priori regions of the state space that are not reachable. The modification of Algorithm 1, aimed at linking more closely the decomposition–abstraction phase with the reachability calculus to decompose regions as the computation progresses, as applied with other principles in [15], is under investigation.
1166
M.-A. Lefebvre, H. Gu´eguen / Nonlinear Analysis 65 (2006) 1150–1167
Fig. 12. Reachability for two generating elements (a and b), and three generating elements (c).
All the abstractions of the examples have been computed using specific Matlab macros, and the reachability computations have been performed using a specific program based on the Polylib library [22]. However, a problem that occurs during the implementation is the complexity of the polytopes, and that increases with the number of composed automata and the number of steps in the reachability computation. To enhance the applicability of the approach it should then be useful to study the simplification of the representation. The number and the position of the hyperplanes influence the precision of the reachability, but a second factor that is also important in this approach is the number of abstractions that are composed to perform the calculus, as stated above. An important open problem that will be the subject of further work is the definition of a strategy to choose the initial set of generating elements, and for each element the set of decomposition hyperplanes, and to modify these when the reachability calculus is inconclusive. Should one decomposition be refined (step 1 of Algorithms 3 or 4) and if so which one, or should a new generating element be chosen (step 2 of Algorithm 1) and if so how? For example, according to the region of interest it is not obvious that the result in Fig. 12(c) computed with three generating elements is really better than the result in Fig. 12(a) computed with two generating elements. On the other hand, the latter seems better than the result in Fig. 12(b), which is computed with two other generating elements. References [1] A. Tiwari, Approximate reachability for linear systems, in: HSCC2003: Hybrid Systems Computation and Control, Prague, The Czech Republic, 2003, pp. 514–525. [2] S. Prajna, A. Jadbabaie, Safety verification of hybrid systems using barrier certificates, in: HSCC2004: Hybrid Systems: Computation and Control, Philadelphia, USA, 2004, pp. 477–492.
M.-A. Lefebvre, H. Gu´eguen / Nonlinear Analysis 65 (2006) 1150–1167
1167
[3] A. Chutinan, B. Krogh, Verification of infinite-state dynamics systems using approximate quotient transition systems, IEEE Transactions on Automatic Control 46 (2001) 1401–1410. [4] R. Alur, T. Dang, F. Ivancic, Reachability of hybrid systems via predicate abstractions, in: HSCC2003: Hybrid Systems Computation and Control, Prague, 2003, pp. 35–48. [5] H. Gu´eguen, J. Zaytoon, On the formal verification of hybrid systems, Control Engineering Practice 12 (2004) 1253–1267. [6] T. Henzinger, P. Ho, H. Wong-Toi, Algorithmic analysis of non-linear hybrid systems, IEEE Transactions on Automatic Control 43 (1998) 540–554. [7] H. Yazarel, G. Pappas, Geometric programming relaxations for linear systems reachability, in: American Control Conference, Boston, USA, 2004, pp. 553–559. [8] A. Kurzhanski, P. Varaiya, Ellipsoidal technique for reachability analysis, in: HSCC2000: Hybrid Systems Computation and Control, Pittsburgh, USA, 2000, pp. 202–214. [9] T. Dang, O. Maler, Reachability analysis via face lifting, in: HSCC1998: Hybrid Systems Computation and Control, Berkeley, USA, 1998, pp. 96–109. [10] T. Henzinger, B. Horowitz, R. Majumbar, H. Wong -Toy, Beyond HyTech: Hybrid systems analysis using interval numerical methods, in: HSCC2000: Hybrid Systems Computation and Control, Pittsburgh, USA, 2000, pp. 130–144. [11] E. Asarin, T. Dang, Abstraction by projection and application to multi affine systems, in: HSCC2004: Hybrid Systems: Computation and Control, Philadelphia, USA, 2004, pp. 16–31. [12] E. Asarin, T. Dang, A. Girard, Reachability analysis of nonlinear systems using conservative approximation, in: HSCC2003: Hybrid Systems Computation and Control, Prague, Czech Republic, 2003, pp. 20–35. [13] T. Henzinger, P.-H. Ho, H. Wong-Toi, HyTech: A model checker for hybrid systems, International Journal on Software Tools for Technology Transfer 1 (1997) 110–122. [14] M.-A. Lefebvre, H. Gu´eguen, J. Buisson, Structured hybrid abstractions of continuous systems, in: IFAC World Congress B02, Barcelona, Spain, 2002. [15] G. Frehse, Algorithmic verification of hybrid systems past HyTech, in: HSCC2005: Hybrid Systems Computation and Control, Zurich, Switzerland, 2005, pp. 258–273. [16] T. Dang, Verification and synthesis of hybrid systems, Ph.D. Thesis, Verimag, INPG, Grenoble, France, 2000. [17] M. Greenstreet, I. Mitchell, Reachability analysis using polygonal projections, in: HSCC1999: Hybrid Systems Computation and Control, Berg an Dal, The Netherlands, 1999, pp. 103–116. [18] N. Pettit, Analysis of Piecewise Linear Dynamical Systems, Research Studio Press, Taunton, Somerset, England, 1995. [19] G. Nenninger, V. Krebs, Hybrid regions of abstraction of piecewise affine hybrid systems, in: ADPM2000: Automation of Mixed Processes, Dortmund, Germany, 2000, pp. 87–97. [20] J. Lunze, B. Nixdorf, J. Schr¨oder, Deterministic discrete-event representation of linear continuous-variable systems, Automatica 35 (1999) 395–406. [21] M.-A. Lefebvre, Abstractions pour la v´erification de sˆuret´e des syst`emes hybrides, Ph.D. Thesis, IETR-Sup´elec, Univ Rennes 1, Rennes, France, 2004. [22] D.K. Wilde, A library for doing polyhedral operations, IRISA, Rennes, France, 785, 1993.