A New Approach to Nonlinear Invariants for Hybrid Systems Based on the Citing Instances Method

In generating invariants for hybrid systems, a main source of intractability is that transition relations are first-order assertions over current-state variables and next-state variables, which doubles the number of system variables and introduces many more free variables. The more variables, the less tractability and, hence, solving the algebraic constraints on complete inductive conditions by a comprehensive Gröbner basis is very expensive. To address this issue, this paper presents a new, complete method, called the Citing Instances Method (CIM), which can eliminate the free variables and directly solve for the complete inductive conditions. An instance means the verification of a proposition after instantiating free variables to numbers. A lattice array is a key notion in this paper, which is essentially a finite set of instances. Verifying that a proposition holds over a Lattice Array suffices to prove that the proposition holds in general; this interesting feature inspires us to present CIM. On one hand, instead of computing a comprehensive Gröbner basis, CIM uses a Lattice Array to generate the constraints in parallel. On the other hand, we can make a clever use of the parallelism of CIM to start with some constraint equations which can be solved easily, in order to determine some parameters in an early state. These solved parameters benefit the solution of the rest of the constraint equations; this process is similar to the domino effect. Therefore, the constraint-solving tractability of the proposed method is strong. We show that some existing approaches are only special cases of our method. Moreover, it turns out CIM is more efficient than existing approaches under parallel circumstances. Some examples are presented to illustrate the practicality of our method.


Introduction
In the real world, there exist many systems exhibiting mixed discrete-continuous behavior, which cannot be described in a proper way by using either a discrete or continuous model. The notion of a hybrid automaton has been introduced for modeling such systems [1]. For example, embedded systems are often modeled as hybrid systems due to their involvement in both digital control software and analog plants, the physical process of which is often specified in the form of differential equations.
Safety verification is among the most challenging problems in verifying hybrid systems, which consists of asking whether a set of bad (unsafe) states can be reached from a set of initial states. The safety verification problem for systems described by non-linear differential equations is particularly complicated, as computing the exact reachable set is usually intractable [2][3][4], with the exception of some severely restricted sub-classes, such as timed automata [5] and initialized rectangular automata [6].
However, for purposes of verification of safety properties, it often suffices to compute an over-approximation of the reachable set of states-if the over-approximation does not intersect the set of bad states, then the original system will never reach a bad state. So far, the existing approaches are mainly based on approximate reachable set computations [7][8][9] and abstraction [10][11][12][13][14].
An over-approximation of the reachable states is also called an invariant of the system. The most precise invariant of a system is its exact reach set. The standard technique for proving a safety property ϕ is to generate an inductive invariant ψ that implies ϕ. Moreover, invariants have the benefit of avoiding computing the exact reachable set of hybrid systems and are very useful for hybrid systems which are described by non-linear differential equations which cannot be solved analytically.
Some typical techniques for invariants are based on templates, which are used to search for inductive invariant assertions with standard computational techniques in algebraic geometry involving Gröbner basis [3,15]. The intuitive idea behind these techniques is that of first fixing a template assertion (i.e., a parametric polynomial with unknown coefficients and of bounded degree in the system variables) and encoding inductive conditions into constraints on unknown coefficients, such that any solution to these constraints is an inductive assertion. This technique guarantees that the invariant must be found if the system has it in that form. However, the key challenges for these techniques is how to define an inductive condition with completeness and how to efficiently compute an inductive invariant that satisfies the inductive condition. Usually, these two aspects contradict each other; that is, an inductive condition with completeness often encounters the computability or complexity problem. Sankaranarayanan et al.'s approach [3], for example, makes significant strides in generating non-linear polynomial equations as invariants of a hybrid system. Their approach consists of the following four main steps: 1. Guessing a template polynomial of fixed degree as a candidate inductive invariant. 2. Defining complete inductive conditions. 3. Utilizing a comprehensive Gröbner basis to encode complete inductive conditions into constraints on parameters and, therefore, reducing the invariant generation problem to a non-linear constraint solving problem. 4. Solving constraints on parameters. Any solution with respect to the constraints is an inductive invariant.
In order to encode the complete inductive conditions, a comprehensive Gröbner basis is involved, which is an exact but impractical encoding as the construction of the comprehensive Gröbner basis is very expensive. Hence, Sankaranarayanan et al. resorted to more tractable (but incomplete) inductive conditions (as shown in Tables 7 and 8) to avoid the explicit construction of a comprehensive Gröbner basis. However, in doing so, their approach sacrifices completeness [3]. The exact steps of [3] go as follows.
1. Guessing a template polynomial of fixed degree as candidate inductive invariant. 2. Defining alternative inductive conditions-which are incomplete but more tractable-to take the place of complete inductive conditions. 3. Utilizing a Gröbner basis to encode the alternative inductive conditions into constraints on parameters. 4. Solving the constraints. Any solution with respect to the constraints is an inductive invariant.
To address this issue, instead of computing a comprehensive Gröbner basis, this paper presents a new technique, called Citing Instances Method (CIM). By substituting real numbers for the variables in the complete inductive conditions, CIM can easily derive a set of constraints on the parameters, which can be solved more efficiently by the elimination of free variables. This idea is inspired by the Parallel Numerical Method, which was first presented in [16,17]. The paper showed that a theorem can be proved by carrying out a series of numerical verifications of one certain instance. Their method was firstly used to prove geometric theorems; furthermore, many classic theorems, such as Simson's Theorem and Ptolemy's Theorem, were proved using this technique. Differing from the parallel numerical method, CIM is used to generate constraints on parameters. For example, it is well-known that there are many methods to generate a constraint on {a, b} which guarantees that the identity below holds, 1. Zero Polynomial Theorem [3]: By the Zero Polynomial Theorem, three constraint equations are obtained synchronously, {a x in the left side of (1), respectively. Then, equate the left side of (1) to zero (e.g., when x = 0, then 1 + b = 0 is obtained). Hence, four constraint equations can be derived independently: Here is an intuitive explanation of CIM: when {b = −1, a = −3} and the degree of the polynomial in the left side of (1) is 3, the equation must have at most three roots. However, the left side of (1) vanishes over {−1, 0, 1, 2}, which means that the equation (1) has at least four roots. Therefore, the equation (1) must be identically equal to zero .
In this paper, by an instance, we mean a verification of the proposition by substituting numbers for the variables in the proposition. CIM exhibits two interesting features: (1) the constraint equations can be obtained in parallel. As a result, we can make hla clever use of this feature to speed up our computation. In fact, some special assignments to variables tend to generate constraints easily and some parameters may be determined in the early stage. For example, (1): this solution can be used to simplify the rest of the constraint equations on the parameters. Thus, solving constraint equations becomes easier. (2) By substituting real values for the variables, all or some of the variables are, in fact, eliminated by an instance. The main source of the intractability of hybrid systems is that transition relations are first-order assertions on current-state variables and next-state variables, which doubles the number of system variables and may introduce many more free variables and encoding inductive conditions by a comprehensive Gröbner basis will be very expensive. In contrast, CIM can be used to eliminate these free variables. Generally, to automatically discover invariants, the more variables, the less tractability.
Inspired by [16,17], this paper presents a new and complete algorithm for invariant generation in hybrid systems. A lattice array, which is a key notion introduced in this paper and essentially a finite set of instances, will be directly applied to the complete inductive conditions to generate the constraint on parameters, and any solution to the constraint produces an inductive assertion. The main idea of our approach is sketched as follows: 1. Guessing a template of a fixed degree as an invariant template. 2. Defining the complete inductive conditions. 3. Applying a lattice array to the complete inductive conditions to generate the constraint on the parameters. 4. Solving the constraint.
Any solution to the constraint guarantees the template an inductive invariant.
The main contributions of this paper are as follows: Firstly, we propose a complete method for constructing invariants of hybrid systems which can solve the complete inductive conditions rather than alternative incomplete inductive conditions. Secondly, CIM takes the place of a comprehensive Gröbner basis in the invariant generation process, and the tractability of constraint-solving is stronger. On one hand, CIM utilizes instances to eliminate the free variables (fewer variables leads to higher tractability). On the other hand, the generation of constraint equations by CIM is parallelized by their independence from each other. Therefore, we can make a clever use of the parallelizability of CIM to start with constraint equations generated by special instances and spread the solutions to some parameters determined in the early stage to other constraint equations, which hence simplifies the constraint equations obtained; this whole process works a bit like the Domino Effect. Therefore, the tractability of the constraint-solving is stronger than comprehensive Gröbner basis computation. Thirdly, CIM involves less symbolic computation and that is why it requires fewer computational resources.

Related Work
Recently, many researchers have devoted effort towards finding the non-linear invariants of hybrid systems. Based on the theory of ideals and Gröbner bases, Sankaranarayanan et al. [3] presented an approach for generating polynomial equation invariants for hybrid systems with more general (non-linear) polynomial dynamics. To control the complexity of the constraint solving, however, this method has to make a trade-off between the complexity of the invariant generation process and the strength of the resulting invariants and several stronger conditions replace the complete conditions. Differing from Sankaranarayanan et al.'s approach, the paper [18] presented a complete method that was not based on guessing a template. However, it is only complete for linear systems. Without resorting to Gröbner bases, the paper [19] implemented the promising algorithm Fastind. Although its use is limited to a discrete system, Fastind executes significantly faster than implementations using Gröbner bases. Fastind is based on remainder computations over parameterized polynomials, and is still an incomplete method. Not coming singly but in pairs, Kong et al. [15] proposed an approach to automatically generate invariant clusters for semi-algebraic hybrid systems by computing the remainder of the Lie derivative of a template polynomial with respect to its Gröbner basis. The benefit of invariant clusters is that they can precisely overapproximate trajectories of the system. Another approach considered barrier certificates based on different inductive conditions [4,[20][21][22] which can be solved efficiently by sum-of-squares (SOS) programming. The zero level set of barrier certificates forms the boundary of the reach set of a hybrid system and, hence, is an invariant. However, this approach is limited by the conservative inductive condition. On the whole, the verification problem of hybrid systems is undecidable: it is doomed to be impossible to find a universal approach for all hybrid systems. This implies that various inductive invariants and computational methods can be proposed for different classes of hybrid systems with some simplification or restriction. Some other approaches, which focus on different features of systems, have also been proposed for the construction of inductive invariants [23][24][25][26][27].
The paper is organized as follows. Section 2 is devoted to details about the modeling framework and some elementary lemmas on which CIM is based. We introduce the theory of the proposed CIM and describe the algorithms in Section 3. Constraint generation by CIM and techniques for solving these constraint equations are discussed in Section 4. In Section 5, we show that the approach of [3] is a special case of CIM. In Section 6, the technique is illustrated with a few examples. Finally, Section 8 concludes our work and discusses the future work in this direction.

Preliminaries
In this section, we introduce the lemmas on which CIM is based, as well as presenting our computational model of algebraic hybrid systems. First, we clarify some notation used throughout the paper.
We denote by K[x 1 , · · · , x n ] the polynomial ring in n indeterminates {x 1 , · · · , x n } over the field K. For conciseness, we also use boldface lowercase letters to denote vectors throughout the paper (e.g., x = x = (x 1 , · · · , x n )). If n = 1, K[x] is called a univariate polynomial ring; it is called a multivariate polynomial ring when n > 1. Let the variables be ordered as x 1 ≺ x 2 ≺ · · · ≺ x n . Moreover, we assume K = R throughout the paper unless otherwise specified.

Basic Lemmas
Lemma 1 ([16]). Let f (x) and g(x) be univariate polynomials of degree less than n. If there are n + 1 distinct numbers b 0 , b 1 , · · · , b n ∈ K such that for every b k one has f (b k ) = g(b k ) with k = 0, 1, · · · , n. Then, It is very easy to understand Lemma 1. If h(x) = f (x) − g(x) = 0, then h(x) has n + 1 roots at least. Meanwhile, h(x) is a univariate equation, the degree of which is less than n and which has n roots at most. As h(x) is not identical to zero, it contradicts the fact that h(x) has at least n + 1 roots. In this paper, we call the n + 1 distinct numbers {b 0 , b 1 , · · · , b n } n + 1 instances. Lemma 1 shows that, in order to prove that two univariate polynomials are identical, it suffices to test n + 1 instances. Moreover, the n + 1 instances can be chosen arbitrarily. Lemma 1 can be extended to multivariate polynomials. For this purpose, we need the following definition.
Definition 1 (Lattice Array [16]). Suppose S 1 , S 2 , . . . , S m ⊂ K such that |S j | = t j with j ∈ {1, · · · , m}, where |S j | denotes the size of S j . We call the Cartesian product of the above m subsets, the m-dimensional lattice array on K. Clearly, S has t 1 t 2 . . . t m elements.
In order to safely conclude whether a proposition is true or false, Multi-Instances Numerical Verification needs a finite set of instances to test it, which has requirements for both the numbers of instances and the relations between instances. The lattice array is defined for this purpose. We first assign variables freely; once variables are assigned, the instances are determined by Definition 1.
In our algorithm, the size of the lattice array depends on the degree of the pseudo-remainder. The following lemma is given for the purpose of determining an upper bound on the degree.

Lemma 4 ( [16]
). Let f and g be the polynomials in the variables (u 1 , u 2 , . . . , u n , x) over a field K: where m ≥ 1, ≥ 1, a k , b k are polynomials in the variables (u 1 , u 2 , . . . , u n ) over K, and a m , b are not zero polynomials. Then there exists a mechanical method to determine the polynomials P, Q, and R in the variables (u 1 , u 2 . . . u n ) over K such that P f + Qg = R and deg(P, The bound on deg(R, u j ) is estimated as follows: In particular, if 1 ≤ m ≤ 2, deg(R, u j ) ≤ A m can be simply expressed as follows:

Characteristic Set
The mathematical concept of a characteristic set was discovered in the late forties by J.F. Ritt [29]. In the late seventies, the Chinese mathematician Wen-Tsün Wu specialized it with modifications to commutative algebra and demonstrated its power for mechanical theorem proving [30,31].
We call x c , denoted by lv( f ), the leading variable of f . Considering f as a polynomial in x c , we can write f as a n x n c + a n−1 x n−1 where a n , · · · , a 0 are in K[x 1 , · · · , x c−1 ], and a n = 0. We call a n the initial or leading coe f f icient of f and n the leading degree of f , denoting them as lc( f ) and ld( f ), respectively.
x n ] is said to be an ascending set (or chain), if one of the following two conditions holds: 1. r = 1 and F 1 s not identically zero; 2. r > 1 and 0 < class(F 1 ) < class(F 2 ), · · · , < class(F r ) ≤ n, and each F i is reduced with respect to the preceding polynomials, F j s(1 ≤ j < i). [32]). Let P ⊆ K[x] is nonempty set of polynomials, and I = P be the ideal generated by P, ascending set F = [F 1 , · · · , F r ] is the Characteristic Set of P, if F ⊆ P , and prem(P, F) = {0}.

Definition 3 (Characteristic Set
The Wu-Ritt Process [32] described how to obtain such the ascending set. Besides Gröbner basis method, characteristic set provides an alternative algorithmic way for solving multivariate polynomial equations or differential equations.

Hybrid System
In this paper, we adopt the hybrid automata proposed in [1] as our modeling framework. Many other models for hybrid systems can be found in [33][34][35][36].
A hybrid system can be defined as follows: L is a finite set of locations (or modes); • V is a set of real-valued system variables. The hybrid system state space is denoted by Σ = L × R |V| , a state is denoted by , s ∈ Σ, and s ∈ R |V| is a continuous state of the variables over the real numbers; is a map that maps each location ∈ L to a differential rule f ; that is, F( ) = f , of the formv i = f i (x). The differential rule F( ) specifies how the system variables evolve at the location , which is also known as a vector field or a flow field; • F : L → (R |V| → R |V| ) is a map that maps each location ∈ L to a differential rule f ; that is, F( ) = f , of the formv i = f i (x). The differential rule F( ) specifies how the system variables evolve at the location , which is also known as a vector field or a flow field; • I : L → 2 R |V| is a map that maps each location ∈ L to a location condition (location invariant), which is an assertion over V and defines all possible continuous states that the system is allowed to move to while at location ; • Θ is an assertion specifying the initial condition; and • 0 ∈ L is the initial location. We assume that the initial condition satisfies location invariance at the initial location; that is, Θ |= I( 0 ).
The transition and dynamic structures of the hybrid system define a set of trajectories. A trajectory is a sequence of states starting from a state 0 ,x 0 ∈ Θ, where Θ ∈ Σ is an initial state set, consisting of a series of interleaved continuous flows and discrete transitions. During the continuous flows, the system evolves following the vector fields at some location until the invariant condition I( ) is violated. At some state ,x , if there is a discrete transition τ from location to such that ρ τ (x) = true, then the discrete transition can be taken.

such that
Initiation:x 0 |= Θ specifies an initial state. Furthermore, for each consecutive state pair i ,x i , i+1 ,x i+1 , one of the two consecution conditions below is satisfied.
Discrete Consecution: there exists a transition τ : 1 , 2 , ρ τ , α τ ∈ T such that τ is enabled; that is, Continuous Consecution : i = i+1 = ; in other words, the location does not change and there exists a time interval, δ > 0, along which a smooth (continuous and differentiable to all orders) function Ψ : [0, δ] → R |V| exists, such that Ψ evolves fromx i tox i+1 according to the differential rule at location while satisfying the location condition I( ). Formally, A hybrid state ,x is called a reachable state if it appears in some trajectory of H. A linear constraint over V is an inequality of the form a 1 x 1 + . . . + a n x n + b ≤ 0, and a linear assertion over a set of variables V is a conjunction of linear constraints over V. The set of points satisfying a linear assertion forms a polyhedron.
A non-linear constraint is an inequality of the form P ≤ 0, where P is a polynomial in {x 1 , . . . , x n }. The constraint is said to be algebraic if P = 0. A non-linear assertion is a conjunction of non-linear constraints. The set of points satisfying a non-linear assertion is called a semi-algebraic set .
Throughout the paper, given an assertion ψ over the variables V, ψ denotes the assertion obtained by replacing each variable v ∈ V by v ∈ V .
General hybrid systems can be specialized into algebraic hybrid systems such that, for each transition τ, the transition relation ρ τ is an algebraic assertion over V ∪ V and the initial condition Θ is an algebraic assertion over V. Definition 6 (Algebraic Hybrid Systems ). An algebraic hybrid system is a hybrid system V, L, T , 0 , Θ , where: 1. For each transition τ : 1 , 2 , ρ τ , α τ , α τ is an algebraic assertion. 2. The initial condition Θ and the location conditions I( ) are also algebraic assertions.
This paper focuses on algebraic hybrid systems.
Definition 7 (Invariant). An invariant of a hybrid system H at a location is an assertion ψ such that, for any reachable state ,x at ,x |= ψ.
Template polynomials play an important role in Sankaranarayanan et al.'s approach. Given a degree d, a template polynomial is essentially a generic degree d polynomial; that is, p = ∑ i 1 +···+i n ≤d L (a 1 , · · · , a ) x i 1 1 ...x i n n , where {a 1 , · · · , a } are coefficients to be decided. Sankaranarayanan et al. treat these coefficients {a 1 , · · · , a } as unknowns and encode inductive conditions by a Gröbner basis to generate constraints on the coefficients such that any solution corresponds to an inductive assertion. Formally, Definition 8 (Template ). Let a = (a 1 , · · · , a l ) be template variables and L(a) be the polynomial of the form Hybrid systems generally consist of many locations and, hence, an invariant can be seen as a mapping to map each location to an assertion which is true under any system state reaching the location. Thus, the following two definitions come very naturally.
Definition 9 (Algebraic Assertion Map ). Given a domain of algebraic assertions D, an algebraic assertion map for an algebraic hybrid system is a map η : L → D that associates each location of the system with an algebraic assertion, where each algebraic assertion η( ) is of the form p = 0.
For simplicity, we shall use η( ) to denote both the algebraic assertion p = 0 and the polynomial p, as long as it will not cause any ambiguity.
Definition 10 (Template Map). Let H ≡ V, L, T , 0 , Θ be an algebraic hybrid system. Assuming a set of template variables a = (a 1 , · · · , a l ), a template map over H is a map η : L → L(a)[V] that maps each location in L to a template over {a, V}.
It is a well-known fact, from the pioneering work of Floyd and Hoare [37,38], that if η is an inductive assertion map, then η( ) is invariant at . In fact, all known invariant generation methods are inductive assertion generation methods.
An important concept used in this paper is the Lie derivative. In our context, the Lie derivative evaluates the change of a scalar function ϕ(x) along the flow of a vector field of the formẋ = f (x), where f (x) = ( f 1 (x), . . . , f n (x)). Formally,
However, for the convenience of discussion, this paper divides the {l + d + k} indeterminates into three categories: x = (x 1 , · · · , x k ), u = (u 1 , · · · , u n ), and a = (a 1 , · · · , a l ), which are called dependent variables, independent variables, and template variables (parameters), respectively. The ring of polynomials with coefficients in K can also be written as K[a, u, x].
In Sankaranarayanan et al.'s approach, encoding complete inductive conditions into the constraint by a comprehensive Gröbner basis will lead to computability problems. In our approach, CIM is used to take the place of a comprehensive Gröbner basis in encoding the complete inductive conditions. Concretely, given an algebraic hybrid system H, we first define a d-degree polynomial as template (i.e., ∑ i 1 +···+i n ≤d L(a)x i 1 1 ...x i n n , where a = (a 1 , · · · , a l ) are the template variables to be decided). CIM is an algorithm for finding all the assignmentsâ = (â 1 , · · · ,â ) to the template variables a that guarantee the truth of a formula of the form: where a = (a 1 , · · · , a l ), u = (u 1 , · · · , u n ), and x = (x 1 , · · · , x k ), f m , G ∈ L(a)[u, x] are all algebraic assertions. The indeterminates in formula (5) are divided into three groups: a denotes the template variables (parameters), u denotes independent variables, and x denotes dependent variables, as they are constrained by k polynomial equations in the hypothesis of formula (5). Note: For an invariant generation, the implication (5) will contain more system variables, as transition relations are algebraic constraints on both current-state and next-state variables, which doubles the number of system variables and introduces many more free variables.

Basic Lemma and Theorem
We first prove a lemma. Lemma 5. Let R ∈ L(a)[u] with u = (u 1 , · · · , u n ), such that deg(R, u i ) = d i ; S R = S 1 × S 2 × · · · × S n be an n-dimensional lattice array of size (d 1 + 1, · · · , d n + 1); and Ω be the following system of equations over a: whereû i ∈ S R . Then, for any vectorâ,â is a solution to Ω iffâ satisfies the following identity: Proof. =⇒. Suppose thatâ is a solution to the system of equations (6); that is, In other words, R(â, u) vanishes over S R . By Lemma 2, R(â, u) ≡ 0. ⇐=. Suppose that there existsâ which satisfies R(â, u) ≡ 0, then R(â, u) must vanish over S R ; that is, for everyû i ∈ S R , R(â,û i ) = 0.
In other words, (8) holds; that is,â is a solution to (6).  Table 1. Table 1. The six constraint equations generated in parallel.
Problem: Given a formula in the form of (5), how do we compute all the possible valuesâ for a for which (5) holds?
The traditional method is to construct a Gröbner basis [3] of { f 1 , · · · , f k }. Let I = Ideal({ f 1 , · · · , f k }) be the ideal generated by { f 1 , · · · , f k } and G be a polynomial. By Hilbert's Nullstellensatz [40], the formula (5) holds is equivalent to that there exists an integer m ≥ 1 such that G m belongs to I. Therefore, according to Hilbert's Nullstellensatz, to compute all the possibleâ that make the formula (5) hold, one has to enumerate all the m ≥ 1 to find all theâ that make G m ∈ I based on Gröbner basis, which is apparently not possible because there is an infinite number of m.
For the above reason, Sankaranarayanan et al. chose to set m = 1 to find all theâ that make G ∈ I in [3]. There are two drawbacks to this approach. Firstly, since G ∈ I is just a sufficient condition for the formula (5) to hold, the approach can only find part of the solutions to (5) which is not complete.
Secondly, computingâ in this way involves constructing a comprehensive Gröbner basis [41], a variant of the Gröbner basis, which is very expensive (double exponential in the dimension of the variables). To avoid this issue, the basic idea of CIM is "turning the difficulty of quality into the complexity of quantity". We now introduce our Theorem 1: Theorem 1. Given x 1 ≺ x 2 ≺ · · · ≺ x k ≺ u 1 , · · · , ≺ u n ≺ a 1 , · · · , ≺ a l , let F = { f i , i = 1, · · · , k} be defined as in formula (5) and F i (a, u, x 1 , · · · , x i ), i = 1, · · · , k be the ascending set of F obtained by Wu-Ritt's Algorithm [39]. Then, there must exist a sequence of polynomials P i (a, u, x), i = 1, · · · , k, Q (a, u, x), and R(a, u), such that the following equation holds: Moreover, for any assignmentâ to a, the formula (5) holds if R(â, u) ≡ 0 and the following system of equations has no solution Proof. For the convenience of proof, let formula (11) denote the hypothesis of (5) and formula (12) denote the conclusion of the formula (5): G(a, u, x) = 0.
Now, we show how to construct such a series of polynomials P i , Q, and R ∈ L(a)[u] satisfying equation (9).
By applying Lemma 4 to F k and G, we can obtain P k F k + Q k G = R k (a, u, x 1 , . . . , x k−1 ). Similarly, by applying Lemma 4 to F k−i and R k−i repeatedly, we get the following sequence of equations, .............
Next, by a sequence of substitutions of R k−i in (13), in a top-down order, we can easily derive the following equation: Now, let R = R 1 , P 1 = P 1 , P j = ∏ j−1 i=1 Q i P j , j = 2, . . . , k, and Q = ∏ k j=1 Q j . Then, we can obtain the following equation: Suppose that there exists an assignmentâ to a that satisfies R(â, u) ≡ 0 and the system of equations (10) has no solution; we prove the formula (5) holds. As R ≡ 0, we get In addition, as the system of equations (10) has no solution, for any (û,x), F i (â,û,x) = 0 we can deduce that Q(â,û,x) = 0. Then, by the equation (15), we deduce that G(â,û,x) = 0 as well. Therefore, formula (5) holds. Note that, in Wu-Ritt's algorithm, Q = ∏ k j=1 Q j = 0 is called the non-degenerate condition [39].

Remark 1.
According to Theorem 1, we know that the solutionsâ to a that satisfy R(â, u) ≡ 0 and the system of equations (10) have no solution satisyfing formula (5) as well. Therefore, the idea for findingâ to make formula (5) hold is that we must first solve the constraint on a derived from R(a, u) ≡ 0 to find all the solutionsâ and decide whether the system of equations (10) has a solution, given a =â. For the first step, we use Lemma 5. Regarding whether the non-degenerate condition is true or not, we can make use of the criterion presented in [39] (i.e., I 1 I 2 · · · I k = 0, I i is the initial of F i ). From the above analysis, we can see that the solution to R(a, u) ≡ 0 cannot ensure that formula (5) holds. If the non-degenerate condition is false, formula (5) may or may not hold, which is crucial to the theorem's proof by CIM or by Wu-Ritt's Algorithm [39]. Fortunately, for invariant generation, it only means we need another step to verify the obtained invariant. Generally, verifying invariants is less expensive than generating invariants, which also means that our invariant generation algorithm is a two-phase algorithm.

Generate the Constraint for Implication by CIM
By Lemma 5, in order to generate a constraint over the template variables a that guarantees R(a, u) to be identical to zero, we only need test whether R(a, u) vanishes over a lattice array S R of size (deg(R, u 1 ) + 1, deg(R, u 2 ) + 1, · · · , deg(R, u n ) + 1), where the upper bound of deg(R, u i ) can be estimated by Lemma 4.
More specifically, takeû ∈ S R , substitute for u in (11) and (12), compute R(a) by the division algorithm (Lemma 4), and then equate R(a) = 0 to obtain the constraint equations. This process is similar to Wu's division method, but is simpler as the u = (u 1 , · · · , u n ) have been replaced by the numbers (instances). However, the aforementioned benefits come at a cost. The benefits are that the free variables are eliminated, the cost is that the implication (5) turns into |S R | simpler implications without free variables. This is the main idea of CIM: "turning the difficulty of quality into the complexity of quantity". If (5) is too difficult to deal with, CIM is meaningful. We outline our algorithm for generating constraints by computing R in Algorithm 1.

Algorithm 1: Generating constraint for implication by R.
input : The hypothesis of implication f m (a 1 , · · · , a , u 1 , · · · , u n , x 1 , · · · x k ), the conclusion of implication G(a 1 , · · · , a , u 1 , · · · , u n , x 1 , · · · , x k ) = 0, and A i the upper bound of deg(R, u i ), U is the set of independent variables, MaxA is the maximum ofA i , inst is a an n-dimension vector of the form û 1 , · · · ,û n , LA is the set of n-dimension vector. output : C, the set of constraints on template variables {a 1 , · · · , a }, a log file to record the process details. 1 /*by default C and LA are ∅ */; 2 C = ∅ ; 3 LA = ∅ ; 4 /*The nested loop structure aims to generate a lattice array for independent variables of {u 1 , · · · , u n } of size (MaxA n ); for simplicity, u i = 0, · · · , MaxA */ 5 for i = 0 to MaxA do 6 /*generate an instance*/  1. u = {u 1 , · · · , u n } and x = {x 1 , · · · , x k } are independent variables and dependent variables, respectively. 2. Lines 5-13 aim to generate Lattice Array. 3. Lines 16-22 aim to substitute real values for {u 1 , · · · , u n } in f i and G and, hence, eliminate the independent variables. It is obvious that instances will simplify the computation of proposition (5). 4. Lines 14-33 are a ForEach loop, hence, the algorithm is easy to parallelize. 5. To make Algorithm 1 simpler, we can use Max(deg( f m , u i ), deg(G, u i )) to take the place of A i , the upper bound of deg(R, u i ), which is not allowed for theorem proving, as |S R | is crucial to theorem proving by CIM. However, for invariant generation, it only means we need to verify the obtained constraint in the second phase (Section 4).

Complexity Analysis
Assuming that the degree of each variable in f i is no more than d (i.e., deg( f i ) ≤ d, 1 ≤ i ≤ k), and no more than δ in g (i.e., deg(g) ≤ δ), our method consists of three main steps. In the first step, we transform the equations (11) into triangular sets (i.e., equations in triangular form) by Wu-Ritt's method. By [32], the complexity of computing the characteristic set of (11) is O(k O( +n+k) (d +

1) O ( +n+k) 3
); however, CIM eliminates the independent variables {u 1 , · · · , u n }, such that the complexity of computing the characteristic set for every instance is O(k O( +k) (d + 1) O (( +k) 3 ) ). As the size of the lattice array is n (d+1) at most, the total complexity is O(n (d+1) k O( +k) (d + 1) O ( +k) 3 ). The second step is the process of constructing R by applying the pseudo-division algorithm to f and g, which can be decomposed into a series of steps to eliminate the highest power in x i . Assuming deg(F i ) ≤ λ, by [32], the complexity of computing R is O(δ O( +n+k) (λ + 1) O (( +n+k)k) ) at most, and CIM eliminates the independent variables {u 1 , · · · , u n } and, so, the complexity of computing R is O(n (d+1) δ O( +k) (λ + 1) O (( +k)k) ). Thus, complexities of Wu-Ritt's method and CIM in steps 1 and 2 are ), respectively. Therefore, according to the above analysis, we can see that our approach is exponential in the degree of the involved polynomials, the number of variables as well as the number of the involved polynomials. However, according to [42], the complexity of computing Gröbner basis is d 2 ( +n+k+O(1)) , which is double exponential. Therefore, the approach in [3] is obviously more expensive than CIM approach.

Illustration of CIM
In this section, we illustrate CIM by an example taken from [3], and discuss some skills in computing R by different kinds of implications generated in every step. Note that the main idea of CIM is "turning the difficulty of quality into the complexity of quantity" and, hence, can solve the complete inductive conditions, under which the approach by a comprehensive Gröbner basis in [3] is intractable.
F(l) = (ẏ = v y ∧δ = 1 ∧v y = −10), Example 3 (Bouncing Ball, from [3]). Figure 2 shows a graphical representation of a ball bouncing on a soft floor, which can be modeled as a hybrid system. The variable y represents the position of the ball (obviously, y = 0 represents the ball being at floor level), v y represents its velocity, and δ denotes the time elapsed since its last bounce. A bounce is modeled by the transition τ, in which the velocity v y of the ball is halved and the ball reverses direction. Figure 1. The hybrid automaton for a bouncing ball Furthermore, for each consecutive state pair l i , x i , l i+1 , x i+1 , one of the two consecution conditions below is satisfied.
Discrete Consecution: there exists a transition τ : 1 , 2 , ρ τ ∈ T such that l i = 1 , l i+1 = 2 , and x i , x i+1 |= ρ τ , where the unprimed variables refer to x i and the primed variables to x i+1 , or Continuous Consecution: l i = l i+1 = , and there exists a time interval δ > 0, along with a smooth (continuous and differentiable to all orders) function f : [0, δ] → R n , such that f evolves from x i to x i+1 according to the differential rule at location , while satisfying the location condition I( ). Formally, Our approach consists of the following steps: 1. Guessing a template of fixed degree as an invariant template. 2. Defining the complete inductive conditions. 3. Applying a lattice array to the complete inductive conditions to generate the constraint on the parameters 4. Solving the constraint. Any solution to the constraint guarantees the template an inductive invariant Different from the existing approaches, we use CIM to encode the complete inductive conditions.
Step 1. Predefine the template map Just as in [3], we set the degree of the invariant at 2, as there is only one location l. We set η(l) to be a generic quadratic form on y, v y , δ , as follows: η (l) : a 1 y 2 + a 2 v y 2 + a 3 δ 2 + a 4 yv y + a 5 v y δ +a 6 yδ + a 7 y + a 8 v y + a 9 δ + a 10 .
Step 2. Generate constraint for the initial condition In (16), there are three variables, which are constrained by three equations in Θ. Thus, there are no independent variables. By Algorithm 1, R = a 10 + 256a 2 + 16a 8 and, then, we obtain the constraint a 10 + 256a 2 + 16a 8 = 0; the same result as in [3].
Step 3. Encode the discrete consecution By Definition 11, the discrete consecution can be expressed by the following implication: 1. In (17), there are six variables that are constrained by five equations. We might as well assume that {y, y , v y , δ, δ } are dependent variables, and {v y } is the independent variable. The degree of {v y } is 2 and the size of the lattice array is (2+1). and, so, three implications are obtained (Table 2). 2. Applying Algorithm 1. We get three simpler implications and three corresponding R by Lemma 3 ( Table 2). Table 2. Implications and R generated by Citing Instances Method (CIM) on the encoding initial condition.
Step 4. Encode the Continuous Consecution By Definition 11, the continuous consecution can be expressed by the implication whereη(l) is the Lie Derivative of η(l) with respect to F(l), that is,η ( ) = a 4 v 2 y + 2a 1 yv y + a 6 δv y + (−20a 2 + a 5 + a 7 )v y + (2a 3 − 10a 5 )δ + (−10a 4 + a 6 )y + (a 9 − 10a 8 ). (19) In (18), there are three variables which are constrained by one equations (we ignore I(l) = y ≥ 0). We might as well assume that {v y } is the dependent variable and that {y, δ} are independent variables. The degree of y and δ is 2 and the size of the lattice array is (2+1,2+1). Applying Algorithm 1 to (18). Thus, nine implications are obtained (Table 3 ). We analyze some implications to illustrate the domino effect: y + (a 5 + a 7 )v y + a 9 = 0 By Algorithm 1,F i is a constant, the consequent of implication is R, i.e., R = a 4 v 2 y + (a 5 + a 7 )v y + a 9 , and the constraint equations are {a 5 = −a 7 , a 4 = 0, a 9 = 0}. We will use this solution to simplify the remaining eight implications and the same below. Thus, the resulting implications become more and more simple (Table 4), we call it the domino effect, which is the benefit of parallelism in CIM.

A Special Case
In this section, we show that Sankaranarayanan et al.'s approach is a special case of our method. In Sankaranarayanan et al.'s approach, the antecedent of the discrete consecution implication {η ( i ) = 0 ∧ ρ τ ∧ α τ |= η j = 0} and continuous consecution implication {I( ) ∧ η( ) = 0 |=η( ) = 0} contain template variables, which requires the construction of a comprehensive Gröbner basis, a variant of the Gröbner basis [41]. Unfortunately, encoding inductive conditions by comprehensive Gröbner bases is an exact but impractical approach, as the non-linear constraints produced make the constraint-solving problem intractable, even for a simple hybrid system. Tables 7 and 8) to eliminate template variables in the antecedent of implication. In doing this, (17), (18), and (5) are transformed into (20), (21), and (22), respectively. Obviously, (22) is a special case of (5), the antecedent of which does not contain any more parameters. In particular, if the transition relations are separable [3]-that is, each variable in V is expressed as a polynomial expression over the variables in V-we can generate constraints more easily by CIM, as the antecedent of implication has been the ascending set.

Name Alternative Continue Consecution Condition
We apply Algorithm 1 to (20) and (21). 1. For the initial condition case, it is same as Step 1.

Remark 2.
The main shortcoming of our approach lies in that CIM tends to produce many redundant instances, which results in more time consumption than Sankaranarayanan et al.'s approach (4.9 s vs. 3.7 s in Maple, respectively). However, CIM is intrinsic to be parallelized; under parallel circumstances, our method is less time-consuming (0.9 s vs. 3.7 s in Maple). In addition, the generation procedure for constraint equations in our approach is simpler (which can even be done by hand), as numerical computation takes the place of symbolic computation.

Experiments
To show the practicality of CIM, we present another two application examples. One is a train system, the other is a charged particle in a magnetic field.

Experiment 1
Consider a train system (from [3]). Figure 3 shows a hybrid automaton modeling a train accelerating (location l 0 ), traversing at constant speed (l 1 ), and decelerating (l 2 ). There exist three continuous variables in the train system: x, the position of the train; v, the train's velocity; and t, a master clock. The system has one discrete variable s, representing the number of stops made so far. The initial condition is given by {x = s = v = t = 0}. There are three discrete transitions τ 1 , τ 2 , and τ 3 ; the transition relations are as follows: Train System: Figure 6 shows a hybrid automaton modeling a train accelerating (location l 0 ), traversing at constant speed (l 1 ), and decelerating (l 2 ). The system has three continuous variables: x, the position of the train, v, the train's velocity, and t, a master clock. The system has one discrete variable s, representing the number of stops made so far. The initial condition is given by x = s = t = v = 0. There are three discrete transitions, τ 1 , τ 2 , and τ 3 , with transition relations wherein id(X) denotes x∈X (x = x). Application of our technique resulted in the following assertion map: With v = 5 at l 1 the assertion η(l 1 ) can be simplified to 4x + 115s − 20t + 25 = 0. An analytic argument for the assertion η(l 0 ) is as follows. Consider the system at the state l 0 , s, v, x, t . Each stop s consists of accelerating from 0 to 5 in l 0 and decelerating from 5 to 0 in l 2 . The distance covered in these two modes is 25 4 and 25 2 respectively. Furthermore, accelerating from 0 to v in l 0 advances the position another v 2 4 . Hence the total distance traveled is given by x = s( 25 4 + 25 2 ) + v 2 4 + 5t l1 , where t l1 , the time spent in location l 1 is given by t − t l0 − t l2 = t − (v/2 + s( 5 2 )) − ( 5 1 s + 2s). Substituting, we obtain the inductive differential.tex; 18/02/2005; 9:15; p.31 η (l i ) , i = 0, 1, 2 is the template map, a 1 x 2 + a 2 xs + a 3 xt + a 4 xv +a 5 s 2 + a 6 st + a 7 sv +a 8 t 2 + a 10 v 2 + a 11 x +a 12 s + a 13 t + a 14 v + a 15 where λc i denotes continuous consecution λ at location l i , and λd ij denotes discrete consecution λ for translation from l i to l j . When b 11 = 1, λd 12 = λd 20 = 1, the invariant is obtained as follows: which the same as that derived in [3].

Experiment 2
In this experiment, we consider a charged particle in a magnetic field (from [3]). Figure 4 shows a charged particle in a 2D-plane with a reflecting barrier at x = 0 and a magnetic field at x ≤ d ≤ 0. There are eight system variables in total: the particle's position x, y; its velocity v x , v y ; a bounce counter b, which is incremented every time the particle collides against the barrier at x = 0; along with the parameters a, d, and time t. There also exist three locations: magnetic, right, and left. 729 instances are generated by Algorithm 1 and the following invariants were obtained at last: which are the same as those obtained in [3].

A Clever Use of CIM
The other method is that, in order to generate constraints without computing R, we can make a clever use of the left side of (9); that is, we only need to make G vanish over a lattice array S R . Once G vanishes over S R , R naturally vanishes over S R , which means that R ≡ 0, according to Lemma 2. More specifically, for eachû ∈ S R , first substitute it for u in (11) (this also means that free variables are eliminated). Then, compute the roots r(a 1 , · · · , a n ) of x in a and, finally, substitute r(a 1 , · · · , a n ) for x in (12). Now, we have a polynomial constraint G(a 1 , · · · , a n ) = 0 in a. Solving G(a 1 , · · · , a n ) = 0 generates a set of solutions to a which satisfy formula (5). However, there is an issue with this approach: r i is not always easy to compute even though the free variables are replaced by numbers. In that case, we can experimentally obtain the roots of dependent variables. Here is an example.
We illustrate, in this example, how to experimentally obtain the roots of dependent variables. Figure 5 shows a graphical representation of an accumulator with varying accumulation every second, which can be modeled as a hybrid system where t is a continuous-time variable, s is the accumulator, and i is the varying accumulation; the system will terminate when i increases to 100. Obviously, Figure 5 has an invariant s = i(i+1) 2 . The target template invariant is a generic degree-two template polynomial: η = a 0 s 2 + a 1 si + a 2 st + a 3 s + a 4 i 2 + a 5 it + a 6 i + a 7 t 2 + a 8 t + a 9 . For the discrete consecution condition, we have assuming {s, i , s , t, t } are the dependent variables and {i} is the independent variable. By CIM, we need compute the roots of {s, i , s , t, t } under the lattice array {i = 0, 1, 2, 3}.
Generally speaking, it is difficult to collect accurate continuous states during the dynamics of a system. When discrete translations are taken, however, we can easily collect continuous states and discrete states (as discrete translations are often controlled by computer software), such that we can collect the roots of {s, s , i , t, t } in the log file by adding two statements to the discrete translation ( Figure 6). The starting point is to run the system with i = 0, s = 0, t = 0 and collect the roots of {s, i , s , t, t } in the log file (Table 9) during every discrete transition.  Finally, we obtain the invariant s = i(i+1) 2 . Without directly computing the roots, CIM can dramatically reduce the complexity. Similar ideas appeared in [43,44].

Conclusions
In this paper, we presented a new approach (called CIM) for invariant generation in hybrid systems. Some examples are given to illustrate how CIM works. The cornerstones of our technique are theories based on the solution number of polynomial equations and Wu-Ritt's Well-Ordering Theorem.
CIM can take the place of a comprehensive Gröbner basis for invariant generation. Furthermore,