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

Abstract

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.

1. 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 $\phi$ is to generate an inductive invariant $\psi$ that implies $\phi$. 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:

• Guessing a template polynomial of fixed degree as a candidate inductive invariant.
• Defining complete inductive conditions.
• 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.
• 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.

• Guessing a template polynomial of fixed degree as candidate inductive invariant.
• Defining alternative inductive conditions—which are incomplete but more tractable—to take the place of complete inductive conditions.
• Utilizing a Gröbner basis to encode the alternative inductive conditions into constraints on parameters.
• 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,

$$x^3-4x^2+2x+1-(x-1)(x^2+ax+b)\equiv 0.$$

(1)

• Zero Polynomial Theorem [3]: $x^3-4x^2+2x+1-(x-1)(x^2+ax+b)\iff x^3-4x^2+2x+1-(x^3+(a-1)x^2+(b-a)x-b)\iff -(a+3)x^2+(2-b+a)x+1+b$. By the Zero Polynomial Theorem, three constraint equations are obtained synchronously, $\{a+3=0,2-b+a=0,1+b=0\}$, solving $\{b=-1,a=-3\}$.
• CIM (this paper): Take four arbitrary distinct values $\{-1,0,1,2\}$ of x and substitute them for 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: $\{-4-2a+2b=0,1+b=0,0=0,-7-2a-b=0\}$. Solving, we get $\{b=-1,a=-3\}$. 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, $x=0\Rightarrow b=-1$ in (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:

• Guessing a template of a fixed degree as an invariant template.
• Defining the complete inductive conditions.
• Applying a lattice array to the complete inductive conditions to generate the constraint on the parameters.
• 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.

2. 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 $\mathbb{K}[x_1,\cdots ,x_n]$ the polynomial ring in n indeterminates $\{x_1,\cdots ,x_n\}$ over the field $\mathbb{K}$. For conciseness, we also use boldface lowercase letters to denote vectors throughout the paper (e.g., $\overrightarrow{x}=\mathit{x}=(x_1,\cdots ,x_n)$). If $n=1,\mathbb{K}\left[\mathit{x}\right]$ is called a univariate polynomial ring; it is called a multivariate polynomial ring when $n>1$. Let the variables be ordered as $x_1\prec x_2\prec \cdots \prec x_n$. Moreover, we assume $\mathbb{K}=\mathbb{R}$ throughout the paper unless otherwise specified.

For a multivariate polynomial ring, $\mathbb{K}\left[\mathit{x}\right]=\mathbb{K}[x_1,\cdots ,x_{k-1},$ $x_k,x_{k+1},\cdots ,x_n]=\mathbb{K}[x_1,\cdots ,x_{k-1},x_{k+1},\cdots ,x_n]\left[x_k\right]$. This means that a multivariate polynomial $p\in \mathbb{K}\left[\mathit{x}\right]$ can be written as a univariate polynomial; that is, $p=\sum a_{i_1{i}_2\dots i_n}{x}_1^{i_1}\dots x_n^{i_n}=c_m{x}_k^m+c_{m-1}{x}_k^{m-1}+\cdots +c_0={\displaystyle \sum _{i=0}^m}{c}_i{x}_k^i$, where $c_i\in \mathbb{K}[x_1,\cdots ,x_{k-1},x_{k+1},\cdots ,x_n],m=deg(p,x_k)$. For example, $p(x_1,x_2,x_3,x_4)=x_2^5+x_3^4{x}_4^2+(2x_2+x_1)x_4^3$ can be considered as a univariate polynomial in $x_2$, namely $p(x_1,x_2,x_3,x_4)=x_2^5+2x_4^3{x}_2+x_1{x}_4^3+x_3^4{x}_4^2.$

An atomic algebraic assertion $\varphi$ over $\mathbb{K}\left[\mathit{x}\right]$ has the form of $p(x_1,x_2,\cdots ,x_n)=0.$ An algebraic assertion is a conjunction of atomic algebraic assertions (i.e., ${\bigwedge }_i{p}_i(x_1,x_2,\cdots ,x_n)=0,$ where $p_i\in \mathbb{K}\left[\mathit{x}\right]$). We denote by $\widehat{\mathit{x}}$ an assignment of $\mathit{x}$.

2.1. Basic Lemmas

Lemma 1

([16]). Let $f\left(x\right)$ and $g\left(x\right)$ be univariate polynomials of degree less than n. If there are $n+1$ distinct numbers $b_0,b_1,\cdots ,b_n\in \mathbb{K}$ such that for every $b_k$ one has $f\left(b_k\right)=g\left(b_k\right)$ with $k=0,1,\cdots ,n$. Then,

$$f\left(x\right)\equiv g\left(x\right).$$

(2)

It is very easy to understand Lemma 1. If $h\left(x\right)=f\left(x\right)-g\left(x\right)=0$, then $h\left(x\right)$ has $n+1$ roots at least. Meanwhile, $h\left(x\right)$ is a univariate equation, the degree of which is less than n and which has n roots at most. As $h\left(x\right)$ is not identical to zero, it contradicts the fact that $h\left(x\right)$ has at least $n+1$ roots. In this paper, we call the $n+1$ distinct numbers $\{b_0,b_1,\cdots ,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,\dots ,S_m\subset \mathbb{K}$ such that $|S_j|=t_j$ with $j\in \{1,\cdots ,m\}$, where $|S_j|$ denotes the size of $S_j$. We call the Cartesian product of the above m subsets,

$$S=S_1\times S_2\times \dots \times S_m,$$

the m-dimensional lattice array on $\mathbb{K}$. Clearly, S has $t_1{t}_2\dots t_m$ elements.

For the convenience of discussion in what follows, we denote $\left|S\right|$ by $(t_1,t_2,\dots ,t_m)$.

Lemma 2

(Multi-Instances Numerical Verification [16]). Let $f\left(x_1,\cdots ,x_m\right)$$\in \mathbb{K}\left[x_1,\cdots ,x_m\right]$, $deg(f,x_k)=n_k,n_k\in {\mathbb{Z}}^{+},k=1,2,\cdots ,m$. If there exists an m-dimensional lattice array S of size $(n_1+1,n_2+1,\cdots ,n_m+1)$ such that $f\left(x_1,\cdots ,x_m\right)$ vanishes along S, i.e., for every $\left({\widehat{x}}_1,{\widehat{x}}_2,\cdots ,{\widehat{x}}_m\right)\in S,f({\widehat{x}}_1,{\widehat{x}}_2,\cdots ,{\widehat{x}}_m)=0$, then f is identically equal to zero.

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.

Example 1.

Prove (3) by Lemma 2:

$$(x+y)\left(x-y\right)\equiv x^2-y^2.$$

(3)

1. 

Determining the size of the lattice array. By observation, the degrees of (3) in x and y are both less than or equal to 2, such that the lattice array size is $\left(2+1,2+1\right)$.

2. 

Determining the members of the lattice array. Essentially, the variables $x,y$ can be assigned freely, so we can assign variables which are easy to test. For our convenience, let $x=0,1,2$ and $y=0,1,2$. Then, the nine instances in the lattice array are $(0,0),(0,1),(0,2),(1,0),(1,1),(1,2),(2,0),$$(2,1),$ and $(2,2)$.

3. 

Substituting instances for x and y one by one, test whether the left side of (3) is equal to the right side of (3). If one of them leads to the left side of (3) differing from the right side of (3), it is a counterexample for (3) not being an identity. Otherwise, we can conclude that (3) is an identity.

Note: Arbitrarily taking nine points, such as (0,0), (1,1), (3,2), (4,0), (7,1), (2,2), (9,10), (2,7), and (8,9), can these nine points be used to safely prove (3)? Of course, the answer is no. The reason is that these nine points do not conform to the definition of a lattice array. It is easy to find a counterexample (see Figure 1): $x+2y+3z-7=14x+9y-z-15$ is an equation, rather than an equality; however, it is easy to find (1+1,1+1,1+1) points (lying on the intersection) which satisfy $x+2y+3z-7$ (the left side) = $14x+9y-z-15$ (the right side).

Lemma 3

(Pseudo-Remainder Formula for Multivariate Polynomials [28]). Let $F,G\in \mathbb{K}\left[x_1,\dots ,x_n\right],$ and $x_k$ be a fixed variable, such that

$$\begin{array}{ccc}F& =& a_m{x}_k^m+a_{m-1}{x}_k^{m-1}+\dots +a_0,\hfill \\ G& =& b_{\ell }{x}_k^{\ell }+b_{\ell -1}{x}_k^{\ell -1}+\dots +b_0,\hfill \end{array}$$

where $\ell =deg(G,x_k),m=deg(F,x_k)$, and $a_i,b_i\in \mathbb{K}\left[x_1,\dots ,x_{k-1},x_{k+1}\dots ,x_n\right]$. Let $G\ne 0$ and $m\ge \ell$. To pseudo-divide F by G, there exists an algorithm to obtain two polynomials, Q and R, such that $lc{(G,x_k)}^{s}F=QG+R$, where $lc(G,x_k)$ is the leading coefficient of G with respect to $x_k$, $0\le s\le m-\ell +1$, $deg(R,x_k)<\ell$. Q and R are called the pseudo-quotient and the pseudo-remainder, respectively.

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,\dots ,u_n,x)$ over a field $\mathbb{K}$:

$$\begin{array}{cc}f=& a_m{x}^m+a_{m-1}{x}^{m-1}+\dots +a_0,\hfill \\ g=& b_{\ell }{x}^{\ell }+b_{\ell -1}{x}^{\ell -1}+\dots +b_0,\hfill \end{array}$$

where $m\ge 1,\ell \ge 1$, $a_k,b_k$ are polynomials in the variables $(u_1,u_2,\dots ,u_n)$ over $\mathbb{K}$, and $a_m,b_{\ell }$ 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\dots u_n)$ over $\mathbb{K}$ such that

$$Pf+Qg=R$$

and $deg(P,x)\le \ell -1$, $deg(Q,x)\le m-1$, $deg(R,x)=0$. The bound on $deg(R,u_j)$ is estimated as follows:

Suppose $\ell \ge m$, and

$$\begin{array}{ccc}{A}_0& =& deg\left(f,u_j\right)\hfill \\ A_1& =& deg\left(g,u_j\right)+\left(\ell -m+1\right)A_0\hfill \\ A_{k+1}& =& 2A_k+A_{k-1}.\hfill \end{array}$$

Then $deg\left(R,u_j\right)\le A_m$

In particular, if $1\le m\le 2$, $deg(R,u_j)\le A_m$ can be simply expressed as follows:

$$deg(R,u_j)\le mdeg(g,u_j)+(m\ell -m+1)deg(f,u_j).$$

2.2. 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].

Let $f\in \mathbb{K}[x_1,\cdots ,x_n]$, the class of f, denoted by $class\left(f\right)$, is the largest i such that $x_i$ occurs in f. If $f\in \mathbb{K}$, then $class\left(f\right)=0$. Let $c=class\left(f\right)>0$. We call $x_c$, denoted by $lv\left(f\right)$, the $leadingvariable$ of f. Considering f as a polynomial in $x_c$, we can write f as

$$a_n{x}_c^n+a_{n-1}{x}_c^{n-1}+\cdots +a_0$$

where $a_n,\cdots ,a_0$ are in $\mathbb{K}[x_1,\cdots ,x_{c-1}]$, and $a_n\ne 0$. We call $a_n$ the $initial$ or $leadingcoefficient$ of f and n the $leadingdegree$ of f, denoting them as $lc\left(f\right)$ and $ld\left(f\right)$, respectively.

Definition 2

(Ascending Set [32]). A sequence of polynomials $F=[F_1,\cdots ,F_r]\subseteq \mathbb{K}[x_1,\cdots ,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\left(F_1\right)<class\left(F_2\right),\cdots ,<class\left(F_r\right)\le n$, and each $F_i$ is reduced with respect to the preceding polynomials, $F_j^{\prime }s(1\le j<i).$

Definition 3

(Characteristic Set [32]). Let $P\subseteq K\left[x\right]$ is nonempty set of polynomials, and $I=\langle P\rangle$ be the ideal generated by P, ascending set $F=[F_1,\cdots ,F_r]$ is the Characteristic Set of P, if $F\subseteq \langle P\rangle$, and $prem(P,F)=\left\{0\right\}$.

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.

2.3. 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:

Definition 4

(Hybrid System ). A hybrid system is a tuple $\mathcal{H}:〈V,L,\mathcal{T},{\ell }_0,\mathsf{\Theta }〉$, where

• 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 $\sum =L\times {\mathbb{R}}^{\left|V\right|}$, a state is denoted by $〈\ell ,\mathit{s}〉\in \sum$, and $\mathit{s}\in {\mathbb{R}}^{\left|V\right|}$ is a continuous state of the variables over the real numbers;
• $\mathcal{T}\subseteq L\times L\times 2^{{\mathbb{R}}^{\left|V\right|}}\times ({\mathbb{R}}^{\left|V\right|}\mapsto {\mathbb{R}}^{\left|V\right|})$ is a set of discrete transitions. A discrete transition $\tau =〈{\ell }_1,{\ell }_2,{\rho }_{\tau },{\alpha }_{\tau }〉$ consists of the pre- and post-locations ${\ell }_1,{\ell }_2$, a guard ${\rho }_{\tau }$ (which is a boolean function of the variables V), and an action ${\alpha }_{\tau }$, which is a first-order assertion over $V^{\prime }\cup V$, where V denotes the current-state variables and $V^{\prime }$ denotes the next-state variables;
• $F:L\mapsto ({\mathbb{R}}^{\left|V\right|}\mapsto {\mathbb{R}}^{\left|V\right|})$ is a map that maps each location $\ell \in L$ to a differential rule $\mathit{f}$; that is, $F\left(\ell \right)=\mathit{f}$, of the form ${\dot{v}}_i=f_i\left(\mathit{x}\right)$. The differential rule $F\left(\ell \right)$ specifies how the system variables evolve at the location ℓ, which is also known as a vector field or a flow field;
• $F:L\mapsto ({\mathbb{R}}^{\left|V\right|}\mapsto {\mathbb{R}}^{\left|V\right|})$ is a map that maps each location $\ell \in L$ to a differential rule $\overrightarrow{f}$; that is, $F\left(\ell \right)=\overrightarrow{f}$, of the form ${\dot{v}}_i=f_i\left(\mathit{x}\right)$. The differential rule $F\left(\ell \right)$ specifies how the system variables evolve at the location ℓ, which is also known as a vector field or a flow field;
• $I:L\mapsto 2^{{\mathbb{R}}^{\left|V\right|}}$ is a map that maps each location $\ell \in 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
• ${\ell }_0\in L$ is the initial location. We assume that the initial condition satisfies location invariance at the initial location; that is, $\mathsf{\Theta }\vDash I\left({\ell }_0\right)$.

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 $〈{\ell }_0,{\widehat{\mathit{x}}}_0〉\in \mathsf{\Theta }$, where $\mathsf{\Theta }\in \sum$ 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\left(\ell \right)$ is violated. At some state $〈\ell ,\widehat{\mathit{x}}〉$, if there is a discrete transition $\tau$ from location ℓ to ${\ell }^{\prime }$ such that ${\rho }_{\tau }\left(\widehat{\mathit{x}}\right)=true$, then the discrete transition can be taken.

Definition 5

(Trajectory). A trajectory of a hybrid system $\mathcal{H}$ is an infinite sequence of states $〈\ell ,\widehat{\mathit{x}}〉\in L\times {\mathbb{R}}^{\left|V\right|}$ of the form

$$〈{\ell }_0,{\widehat{\mathit{x}}}_0〉,〈{\ell }_1,{\widehat{\mathit{x}}}_1〉,〈{\ell }_2,{\widehat{\mathit{x}}}_2〉,\dots ,$$

such that

Initiation: ${\widehat{\mathit{x}}}_0\vDash \mathsf{\Theta }$ specifies an initial state.

Furthermore, for each consecutive state pair $〈{\ell }_i,{\widehat{\mathit{x}}}_i〉$, $〈{\ell }_{i+1},{\widehat{\mathit{x}}}_{i+1}〉$, one of the two consecution conditions below is satisfied.

Discrete Consecution: there exists a transition $\tau :〈{\ell }_1,{\ell }_2,{\rho }_{\tau },{\alpha }_{\tau }〉\in \mathcal{T}$ such that τ is enabled; that is, ${\rho }_{\tau }\left({\widehat{\mathit{x}}}_i\right)$= true and ${\widehat{\mathit{x}}}_{i+1}={\alpha }_{\tau }\left({\widehat{\mathit{x}}}_i\right)$. Or,

Continuous Consecution: ${\ell }_i={\ell }_{i+1}=\ell$; in other words, the location does not change and there exists a time interval, $\delta >0$, along which a smooth (continuous and differentiable to all orders) function $\mathsf{\Psi }:[0,\delta ]\mapsto {\mathbb{R}}^{\left|V\right|}$ exists, such that Ψ evolves from ${\widehat{\mathit{x}}}_i$ to ${\widehat{\mathit{x}}}_{i+1}$ according to the differential rule at location ℓ while satisfying the location condition $I\left(\ell \right)$. Formally,

1. 

$\mathsf{\Psi }\left(0\right)={\widehat{\mathit{x}}}_1,\mathsf{\Psi }\left(\delta \right)={\widehat{\mathit{x}}}_2$ and $\forall t\in [0,\delta ],\mathsf{\Psi }\left(t\right)\vDash I\left(\ell \right),$

2. 

$\forall t\in [0,\delta ],\langle \mathsf{\Psi }\left(t\right),\dot{\mathsf{\Psi }}\left(t\right)\rangle \vDash F\left(\ell \right).$

A hybrid state $〈\ell ,\widehat{\mathit{x}}〉$ is called a reachable state if it appears in some trajectory of $\mathcal{H}$.

A linear constraint over V is an inequality of the form $a_1{x}_1+\dots +a_n{x}_n+b\le 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\le 0$, where P is a polynomial in $\{x_1,\dots ,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 $\psi$ over the variables V, ${\psi }^{\prime }$ denotes the assertion obtained by replacing each variable $v\in V$ by $v^{\prime }\in V^{\prime }$.

General hybrid systems can be specialized into algebraic hybrid systems such that, for each transition $\tau$, the transition relation ${\rho }_{\tau }$ is an algebraic assertion over $V^{\prime }\cup V$ and the initial condition $\mathsf{\Theta }$ is an algebraic assertion over V.

Definition 6

(Algebraic Hybrid Systems ). An algebraic hybrid system is a hybrid system $〈V,L,\mathcal{T},{\ell }_0,\mathsf{\Theta }〉$, where:

1. 

For each transition $\tau :〈{\ell }_1,{\ell }_2,{\rho }_{\tau },{\alpha }_{\tau }〉$, ${\alpha }_{\tau }$ is an algebraic assertion.

2. 

The initial condition Θ and the location conditions $I\left(\ell \right)$ are also algebraic assertions.

3. 

Each rule $F\left(\ell \right)$ is of the form $\bigwedge \dot{v_i}=f_i(v_1,\dots ,v_n)$, where $f_i\in \mathbb{R}[v_1,\dots ,v_n].$

This paper focuses on algebraic hybrid systems.

Definition 7

(Invariant). An invariant of a hybrid system $\mathcal{H}$ at a location ℓ is an assertion ψ such that, for any reachable state $〈\ell ,\widehat{\mathit{x}}〉$ at ℓ, $\widehat{\mathit{x}}\vDash \psi$.

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={\sum }_{i_1+\cdots +i_n\le d}\mathcal{L}\left(a_1,\cdots ,a_{\ell }\right)x_1^{i_1}\dots x_n^{i_n}$, where $\{a_1,\cdots ,a_{\ell }\}$ are coefficients to be decided. Sankaranarayanan et al. treat these coefficients $\{a_1,\cdots ,a_{\ell }\}$ 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 $\mathit{a}=(a_1,\cdots ,a_l)$ be template variables and $\mathcal{L}\left(\mathit{a}\right)$ be the polynomial of the form ${\sum }_{\mathit{\beta }}{c}_{\mathit{\beta }}{\mathit{a}}^{\mathit{\beta }}$, where each $c_{\beta }$ is a real-valued coefficient $\mathit{\beta }=({\beta }_1,\cdots ,{\beta }_n),{\beta }_n\ge 0$. A d-degree template over $\{\mathit{a},\mathit{x}\}$ is a polynomial in the variables $\mathit{x}$ with coefficients $\mathcal{L}\left(\mathit{a}\right)$ (i.e., $p={\sum }_{i_1+\cdots +i_n\le d}\mathcal{L}\left(a_1,\cdots ,a_{\ell }\right)x_1^{i_1}\dots x_n^{i_n}$).

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 $\mathcal{D}$, an algebraic assertion map for an algebraic hybrid system is a map $\eta :L\to \mathcal{D}$ that associates each location of the system with an algebraic assertion, where each algebraic assertion $\eta \left(\ell \right)$ is of the form $p=0$.

For simplicity, we shall use $\eta \left(\ell \right)$ 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 $\mathcal{H}\equiv 〈V,L,\mathcal{T},{\ell }_0,\mathsf{\Theta }〉$ be an algebraic hybrid system. Assuming a set of template variables $\mathit{a}=(a_1,\cdots ,a_l)$, a template map over $\mathcal{H}$ is a map $\eta :L\to \mathcal{L}\left(\mathit{a}\right)\left[V\right]$ that maps each location in L to a template over $\{\mathit{a},V\}$.

It is a well-known fact, from the pioneering work of Floyd and Hoare [37,38], that if $\eta$ is an inductive assertion map, then $\eta \left(\ell \right)$ is invariant at ℓ. In fact, all known invariant generation methods are inductive assertion generation methods.

Definition 11

(Inductive Algebraic Assertion Map ). An inductive algebraic assertion map $\eta \left(\ell \right)$ is a map that associates with each location $\ell \in L$ an assertion $\eta \left(\ell \right)$ that holds initially and is preserved by all discrete transitions and continuous flows. More formally, an inductive assertion map satisfies the following requirements:

Initiation: The algebraic assertion at ${\ell }_0$ subsumes the initial condition; that is, $\mathsf{\Theta }\vDash \eta \left({\ell }_0\right)$.

Discrete Consecution: For each transition $\tau :〈{\ell }_i,{\ell }_j,{\rho }_{\tau },{\alpha }_{\tau }〉$ starting from a state satisfying $\eta \left({\ell }_i\right)$, taking τ leads to a state satisfying $\eta \left({\ell }_j\right)$. Formally, $\eta \left({\ell }_i\right)\wedge {\rho }_{\tau }\wedge {\alpha }_{\tau }\vDash \eta {\left({\ell }_j\right)}^{\prime }$.

Continuous Consecution: For every location $\ell \in L$ and states $〈\ell ,{\widehat{\mathit{x}}}_1〉$, $〈\ell ,{\widehat{\mathit{x}}}_2〉$ such that ${\widehat{\mathit{x}}}_1$ evolves from ${\widehat{\mathit{x}}}_2$ according to the differential rule $F\left(\ell \right)$ at ℓ, if ${\widehat{\mathit{x}}}_1\vDash \eta \left(\ell \right)$, then ${\widehat{\mathit{x}}}_2\vDash \eta \left(\ell \right)$. If $\eta \left(\ell \right)$ is an assertion of the form $f_{\ell }\left(\mathit{x}\right)=0$ and $f_{\ell }\in \mathbb{R}[v_1,\dots ,v_{\left|V\right|}]$ is a real-valued smooth function, we can express continuous consecution by the following condition:

$$I\left(\ell \right)\left(\mathit{x}\right)\wedge (f_{\ell }\left(\mathit{x}\right)=0)\vDash \dot{f_{\ell }}\left(\mathit{x}\right)=0.$$

Note that $\dot{f_{\ell }}$ denotes the Lie derivative of $f\left(\ell \right)$ along the vector field $F\left(\ell \right)$.

An important concept used in this paper is the Lie derivative. In our context, the Lie derivative evaluates the change of a scalar function $\phi \left(\mathit{x}\right)$ along the flow of a vector field of the form $\dot{\mathit{x}}=\mathit{f}\left(\mathit{x}\right)$, where $\mathit{f}\left(\mathit{x}\right)=(f_1\left(\mathit{x}\right),\dots ,f_n\left(\mathit{x}\right))$. Formally,

Definition 12

(Lie Derivative).

$${\mathcal{L}}_{\mathit{f}}\phi \left(\mathit{x}\right)\triangleq \frac{\partial \phi }{\partial \mathit{x}}\mathit{f}\left(\mathit{x}\right)=\sum _{i=1}^n\frac{\partial \phi }{\partial x_i}\dot{x_i}=\sum _{i=1}^n\frac{\partial \phi }{\partial x_i}{f}_i\left(\mathit{x}\right).$$

3. Theory of CIM

Assuming $\{a_1,\cdots ,a_l\}$ are the parameters to be decided, there exists a system of k polynomial equations that involve $\{l+d+k\}$ indeterminates in formula (4); that is, $\{a_1,\cdots ,a_l,u_1,\cdots ,u_n,x_1,\cdots ,x_k\}$. Let the $\{l+d+k\}$ indeterminates be ordered as $x_1\prec x_2\prec \cdots \prec x_k\prec u_1,\cdots ,\prec u_n\prec a_1,\cdots ,\prec a_l$. By Wu–Ritt’s Algorithm [39], the indeterminates are divided into two categories: dependent variables ($x_1,\cdots ,x_k$) and independent variables ($u_1,\cdots ,u_n,a_1,\cdots ,a_l$).

$$\begin{array}{c}\underset{m=1}{\overset{k}{\bigwedge }}{f}_m(a_1,\cdots ,a_l,u_1,\cdots ,u_n,x_1,\cdots ,x_k)=0.\hfill \end{array}$$

(4)

However, for the convenience of discussion, this paper divides the $\{l+d+k\}$ indeterminates into three categories: $\mathit{x}=(x_1,\cdots ,x_k)$, $\mathit{u}=(u_1,\cdots ,u_n)$, and $\mathit{a}=(a_1,\cdots ,a_l)$, which are called dependent variables, independent variables, and template variables (parameters), respectively. The ring of polynomials with coefficients in $\mathbb{K}$ can also be written as $\mathbb{K}[\mathit{a},\mathit{u},\mathit{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 $\mathcal{H}$, we first define a d-degree polynomial as template (i.e., ${\sum }_{i_1+\cdots +i_n\le d}\mathcal{L}\left(\mathit{a}\right)x_1^{i_1}\dots x_n^{i_n}$, where $\mathit{a}=(a_1,\cdots ,a_l)$ are the template variables to be decided). CIM is an algorithm for finding all the assignments $\widehat{\mathit{a}}=({\widehat{a}}_1,\cdots ,{\widehat{a}}_{\ell })$ to the template variables $\mathit{a}$ that guarantee the truth of a formula of the form:

$$\begin{array}{c}\hfill \begin{array}{c}\forall (\mathit{u},\mathit{x}):\underset{m=1}{\overset{k}{\bigwedge }}{f}_m(\widehat{\mathit{a}},\mathit{u},\mathit{x})=0\vDash G(\widehat{\mathit{a}},\mathit{u},\mathit{x})=0,\hfill \end{array}\end{array}$$

(5)

where $\mathit{a}=(a_1,\cdots ,a_l)$, $\mathit{u}=(u_1,\cdots ,u_n)$, and $\mathit{x}=(x_1,\cdots ,x_k)$, $f_m,G\in \mathcal{L}\left(\mathit{a}\right)[\mathit{u},\mathit{x}]$ are all algebraic assertions. The indeterminates in formula (5) are divided into three groups: $\mathit{a}$ denotes the template variables (parameters), $\mathit{u}$ denotes independent variables, and $\mathit{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.

3.1. Basic Lemma and Theorem

We first prove a lemma.

Lemma 5.

Let $R\in \mathcal{L}\left(\mathit{a}\right)\left[\mathit{u}\right]$ with $\mathit{u}=(u_1,\cdots ,u_n)$, such that $deg(R,u_i)=d_i$; $S_R=S_1\times S_2\times \cdots \times S_n$ be an n-dimensional lattice array of size $(d_1+1,\cdots ,d_n+1)$; and Ω be the following system of equations over $\mathit{a}$:

$$\mathsf{\Omega }:R(\mathit{a},{\widehat{\mathit{u}}}_i)=0,i=1,\cdots ,\left|S_R\right|,$$

(6)

where ${\widehat{\mathit{u}}}_i\in S_R$. Then, for any vector $\widehat{\mathit{a}}$, $\widehat{\mathit{a}}$ is a solution to Ω iff $\widehat{\mathit{a}}$ satisfies the following identity:

$$R(\widehat{\mathit{a}},\mathit{u})\equiv 0.$$

(7)

Proof. 

⟹. Suppose that $\widehat{\mathit{a}}$ is a solution to the system of equations (6); that is,

$$\left\{\begin{array}{c}R(\widehat{\mathit{a}},{\widehat{\mathit{u}}}_1)=0\hfill \\ R(\widehat{\mathit{a}},{\widehat{\mathit{u}}}_2)=0\hfill \\ \cdots \hfill \\ R(\widehat{\mathit{a}},{\widehat{\mathit{u}}}_{|S_R|})=0\hfill \end{array}\right..$$

(8)

In other words, $R(\widehat{\mathit{a}},\mathit{u})$ vanishes over $S_R$. By Lemma 2, $R(\widehat{\mathit{a}},\mathit{u})\equiv 0$.

⟸. Suppose that there exists $\widehat{\mathit{a}}$ which satisfies $R(\widehat{\mathit{a}},\mathit{u})\equiv 0$, then $R(\widehat{\mathit{a}},\mathit{u})$ must vanish over $S_R$; that is, for every ${\widehat{\mathit{u}}}_i\in S_R,$

$$R(\widehat{\mathit{a}},{\widehat{\mathit{u}}}_i)=0.$$

In other words, (8) holds; that is, $\widehat{\mathit{a}}$ is a solution to (6). □

Example 2.

(from [3]) Compute the constraint on $\{a_1,a_2,a_3\}$ such that $p=(2a_2+3)x_1{x}_2^2+3a_3{x}_2+4(a_3+a_1+10)\equiv 0$.

1. 

By observation, the lattice array size of p is (1+1,2+1). Let $x_1=0,1$ and $x_2=-1,0,1$. Then, six instances in the lattice array are $\left(0,-1\right),\left(0,0\right),\left(0,1\right),\left(1,-1\right)$, $\left(1,0\right),\left(1,1\right)$.

2. 

For each instance, the constraint is obtained by substituting, which consists of the six equations shown in

Table 1. The six constraint equations generated in parallel.

${\mathit{x}}_1$	${\mathit{x}}_2$	${\mathit{a}}_1,{\mathit{a}}_2,{\mathit{a}}_3$
0	0	$4a_3+a_1+10=0$
0	1	$7a_3+a_1+10=0$
0	$-1$	$a_3+a_1+10=0$
1	0	$4a_3+a_1+10=0$
1	1	$7a_3+2a_2+a_1+13=0$
1	$-1$	$a_3+2a_2+a_1+13=0$

.

3. 

Solving the above six equations: $a_1=-10,a_2=-\frac{3}{2},a_3=0$, as in [3].

Problem: Given a formula in the form of (5), how do we compute all the possible values $\widehat{\mathit{a}}$ for $\mathit{a}$ for which (5) holds?

The traditional method is to construct a Gröbner basis [3] of $\{f_1,\cdots ,f_k\}$. Let $I=Ideal\left(\{f_1,\cdots ,f_k\}\right)$ be the ideal generated by $\{f_1,\cdots ,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\ge 1$ such that $G^m$ belongs to I. Therefore, according to Hilbert’s Nullstellensatz, to compute all the possible $\widehat{a}$ that make the formula (5) hold, one has to enumerate all the $m\ge 1$ to find all the $\widehat{a}$ that make $G^m\in 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 $\widehat{a}$ that make $G\in I$ in [3]. There are two drawbacks to this approach. Firstly, since $G\in 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 $\widehat{a}$ 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\prec x_2\prec \cdots \prec x_k\prec u_1,\cdots ,\prec u_n\prec a_1,\cdots ,\prec a_l$, let $F=\{f_i,i=1,\cdots ,k\}$ be defined as in formula (5) and $F_i(\mathit{a},\mathit{u},x_1,\cdots ,x_i),i=1,\cdots ,k$ be the ascending set of F obtained by Wu–Ritt’s Algorithm [39]. Then, there must exist a sequence of polynomials $P_i(\mathit{a},\mathit{u},\mathit{x}),i=1,\cdots ,k$, $Q(\mathit{a},\mathit{u},\mathit{x})$, and $R(\mathit{a},\mathit{u})$, such that the following equation holds:

$$\sum _{i=1}^k{P}_i{F}_i+QG=R.$$

(9)

Moreover, for any assignment $\widehat{\mathit{a}}$ to $\mathit{a}$, the formula (5) holds if $R(\widehat{\mathit{a}},\mathit{u})\equiv 0$ and the following system of equations has no solution

$$\begin{array}{c}\hfill \left\{\begin{array}{ccc}\hfill F_i(\widehat{\mathit{a}},\mathit{u},x_1,\cdots ,x_i)& =& 0,i=1,\cdots ,k,\hfill \\ \hfill Q(\widehat{\mathit{a}},\mathit{u},\mathit{x})& =& 0.\hfill \end{array}\right.\end{array}$$

(10)

Proof. 

For the convenience of proof, let formula (11) denote the hypothesis of (5) and formula (12) denote the conclusion of the formula (5):

$$\underset{m=1}{\overset{k}{\bigwedge }}{f}_m(\mathit{a},\mathit{u},\mathit{x})=0,$$

(11)

$$G(\mathit{a},\mathit{u},\mathit{x})=0.$$

(12)

Now, we show how to construct such a series of polynomials $P_i,Q$, and $R\in \mathcal{L}\left(\mathit{a}\right)\left[\mathit{u}\right]$ satisfying Equation (9).

By applying Lemma 4 to $F_k$ and G, we can obtain ${\tilde{P}}_k{F}_k+Q_{k}G=R_k(\mathit{a},\mathit{u},x_1,\cdots ,x_{k-1})$. Similarly, by applying Lemma 4 to $F_{k-i}$ and $R_{k-i}$ repeatedly, we get the following sequence of equations,

$$\begin{array}{c}{\tilde{P}}_k{F}_k+Q_{k}G=R_k,\hfill \\ {\tilde{P}}_{k-1}{F}_{k-1}+Q_{k-1}{R}_k=R_{k-1},\hfill \\ {\tilde{P}}_{k-2}{F}_{k-2}+Q_{k-2}{R}_{k-1}=R_{k-2},\hfill \\ \dots \dots \dots \dots .\hfill \\ {\tilde{P}}_2{F}_2+Q_2{R}_3=R_2,\hfill \\ {\tilde{P}}_1{F}_1+Q_1{R}_2=R_1,\hfill \end{array}$$

(13)

where $R_1\in \mathcal{L}\left(\mathit{a}\right)\left[\mathit{u}\right]$, $R_i\in \mathcal{L}\left(\mathit{a}\right)[\mathit{u},$$x_1,\cdots ,x_{i-1}],i=2,\cdots ,k$ and ${\tilde{P}}_i$ and $Q_i\in \mathcal{L}\left(\mathit{a}\right)[\mathit{u},$$x_1,\cdots ,x_i],i=1,\cdots ,k$. Note that $R_1$ does not contain the variable $x_i$ any more.

Next, by a sequence of substitutions of $R_{k-i}$ in (13), in a top-down order, we can easily derive the following equation:

$$\begin{array}{c}\hfill {\tilde{P}}_1{F}_1+\sum _{j=2}^k\left(\prod _{i=1}^{j-1}{Q}_i\right){\tilde{P}}_j{F}_j+\left(\prod _{j=1}^k{Q}_j\right)G=R_1.\end{array}$$

(14)

Now, let $R=R_1$, $P_1={\tilde{P}}_1$, $P_j=\left({\prod }_{i=1}^{j-1}{Q}_i\right){\tilde{P}}_j,j=2,\dots ,k$, and $Q=\left({\prod }_{j=1}^k{Q}_j\right)$. Then, we can obtain the following equation:

$$\sum _{i=1}^k{P}_i{F}_i+QG=R.$$

Suppose that there exists an assignment $\widehat{\mathit{a}}$ to $\mathit{a}$ that satisfies $R(\widehat{\mathit{a}},\mathit{u})\equiv 0$ and the system of Equations (10) has no solution; we prove the formula (5) holds. As $R\equiv 0$, we get

$$\sum _{i=1}^k{P}_i{F}_i=-QG.$$

(15)

In addition, as the system of equations (10) has no solution, for any $(\widehat{\mathit{u}},\widehat{\mathit{x}})$, $F_i(\widehat{\mathit{a}},\widehat{\mathit{u}},\widehat{\mathit{x}})=0$ we can deduce that $Q(\widehat{\mathit{a}},\widehat{\mathit{u}},\widehat{\mathit{x}})\ne 0$. Then, by the Equation (15), we deduce that $G(\widehat{\mathit{a}},\widehat{\mathit{u}},\widehat{\mathit{x}})=0$ as well. Therefore, formula (5) holds. Note that, in Wu–Ritt’s algorithm, $Q={\prod }_{j=1}^k{Q}_j\ne 0$ is called the non-degenerate condition [39]. □

Remark 1.

According to Theorem 1, we know that the solutions $\widehat{\mathit{a}}$ to $\mathit{a}$ that satisfy $R(\widehat{\mathit{a}},\mathit{u})\equiv 0$ and the system of Equations (10) have no solution satisyfing formula (5) as well. Therefore, the idea for finding $\widehat{\mathit{a}}$ to make formula (5) hold is that we must first solve the constraint on $\mathit{a}$ derived from $R(\mathit{a},\mathit{u})\equiv 0$ to find all the solutions $\widehat{\mathit{a}}$ and decide whether the system of Equations (10) has a solution, given $\mathit{a}=\widehat{\mathit{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\cdots I_k\ne 0$, $I_i$ is the initial of $F_i$). From the above analysis, we can see that the solution to $R(\mathit{a},\mathit{u})\equiv 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.

3.2. Generate the Constraint for Implication by CIM

By Lemma 5, in order to generate a constraint over the template variables $\mathit{a}$ that guarantees $R(\mathit{a},\mathit{u})$ to be identical to zero, we only need test whether $R(\mathit{a},\mathit{u})$ vanishes over a lattice array $S_R$ of size $(deg(R,u_1)+1,deg(R,u_2)+1,\cdots ,deg(R,u_n)+1)$, where the upper bound of $deg(R,u_i)$ can be estimated by Lemma 4.

More specifically, take $\widehat{\mathit{u}}\in S_R$, substitute for $\mathit{u}$ in (11) and (12), compute $R\left(\mathit{a}\right)$ by the division algorithm (Lemma 4), and then equate $R\left(\mathit{a}\right)=0$ to obtain the constraint equations. This process is similar to Wu’s division method, but is simpler as the $\mathit{u}=(u_1,\cdots ,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,\cdots ,a_{\ell },u_1,\cdots ,u_n,x_1,\cdots x_k)$, the conclusion of implication     $G(a_1,\cdots ,a_{\ell },u_1,\cdots ,u_n,x_1,\cdots ,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 of$A_i$, $inst$ is a an n-dimension vector of the form $\langle {\widehat{u}}_1,\cdots ,{\widehat{u}}_n\rangle$, $LA$ is the set of n-dimension vector. output: C, the set of constraints on template variables $\{a_1,\cdots ,a_{\ell }\}$, a log file to record the process details. 1 /*by default C and $LA$ are ∅ */; 2 $C=\varnothing$; 3 $LA=\varnothing$; 4 /*The nested loop structure aims to generate a lattice array for independent variables of $\{u_1,\cdots ,u_n\}$ of size (${MaxA}^n$); for simplicity, $u_i=0,\cdots ,MaxA$ */ 5 for$i=0$ to $MaxA$do 13 end 14 foreachinst in LAdo 33 end 34 returnC;

• $\mathit{u}=\{u_1,\cdots ,u_n\}$ and $\mathit{x}=\{x_1,\cdots ,x_k\}$ are independent variables and dependent variables, respectively.
• Lines 5–13 aim to generate Lattice Array.
• Lines 16–22 aim to substitute real values for $\{u_1,\cdots ,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).
• Lines 14–33 are a ForEach loop, hence, the algorithm is easy to parallelize.
• 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\left(f_i\right)\le d,1\le i\le k$), and no more than $\delta$ in g (i.e., $deg\left(g\right)\le \delta$), 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 $\mathcal{O}(k^{\mathcal{O}(\ell +n+k)}{(d+1)}^{{\mathcal{O}}^{{(\ell +n+k)}^3}})$; however, CIM eliminates the independent variables $($, such that the complexity of computing the characteristic set for every instance is $\mathcal{O}(k^{\mathcal{O}(\ell +k)}{(d+1)}^{{\mathcal{O}}^{({(\ell +k)}^3)}})$. As the size of the lattice array is $n^{(d+1)}$ at most, the total complexity is $\mathcal{O}(n^{(d+1)}{k}^{\mathcal{O}(\ell +k)}{(d+1)}^{{\mathcal{O}}^{{(\ell +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\left(F_i\right)\le \lambda$, by [32], the complexity of computing R is $\mathcal{O}({\delta }^{\mathcal{O}(\ell +n+k)}{(\lambda +1)}^{{\mathcal{O}}^{((\ell +n+k)k)}})$ at most, and CIM eliminates the independent variables $(u_1,\cdots ,u_n)$ and, so, the complexity of computing R is $\mathcal{O}(n^{(d+1)}{\delta }^{\mathcal{O}(\ell +k)}{(\lambda +1)}^{{\mathcal{O}}^{((\ell +k)k)}})$. Thus, complexities of Wu–Ritt’s method and CIM in steps 1 and 2 are $\mathcal{O}(k^{\mathcal{O}{(\ell +n+k)}^2}{(d+1)}^{{\mathcal{O}}^{{(\ell +n+k)}^4}})$ and $\mathcal{O}(n^{(d+1)}{k}^{\mathcal{O}{(\ell +k)}^2}{(d+1)}^{{\mathcal{O}}^{{(\ell +k)}^4}})$, 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^{(\ell +n+k+\mathcal{O}(1\left)\right)}}$, which is double exponential. Therefore, the approach in [3] is obviously more expensive than CIM approach.

4. 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.

$$\begin{array}{c}V=\{y,v_y,\delta \},\hfill \\ L=\left\{l\right\},\hfill \\ \mathcal{T}=\left\{\tau \right\},where\hfill \\ \tau :\langle l,l,\left[\begin{array}{c}\delta >0\wedge y=0\wedge y^{\prime }=y\\ \wedge v_y^{\prime }=\frac{v_y}{2}\wedge {\delta }^{\prime }=0\end{array}\right]\rangle ,\hfill \\ \mathsf{\Theta }=(y=0\wedge v_y=16\wedge \delta =0),\hfill \\ F\left(l\right)=(\dot{y}=v_y\wedge \dot{\delta }=1\wedge \dot{v_y}=-10),\hfill \\ I\left(l\right)=y\le 0,\hfill \\ {\ell }_0=l.\hfill \end{array}$$

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.

Our approach consists of the following steps:

• Guessing a template of fixed degree as an invariant template.
• Defining the complete inductive conditions.
• Applying a lattice array to the complete inductive conditions to generate the constraint on the parameters
• 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 $\eta \left(l\right)$ to be a generic quadratic form on $\left\{y,v_y,\delta \right\}$, as follows:

$$\eta \left(l\right):\left[\begin{array}{c}{a}_1{y}^2+a_2{v_y}^2+a_3{\delta }^2+a_{4}y{v}_y+a_5{v}_y\delta \hfill \\ +a_{6}y\delta +a_{7}y+a_8{v}_y+a_9\delta +a_{10}\hfill \end{array}\right].$$

Step 2. Generate constraint for the initial condition

$$(y=0\wedge v_y=16\wedge \delta =0)\vDash \eta \left(l\right)$$

(16)

In (16), there are three variables, which are constrained by three equations in $\mathsf{\Theta }$. 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:

$$\begin{array}{cc}\hfill & \eta \left(l\right)=0\wedge {\rho }_{\tau }\wedge {\alpha }_{\tau }\vDash \eta {\left(l\right)}^{\prime }=0,\hfill \\ \hfill & {\rho }_{\tau }:y=0,\hfill \\ \hfill & {\alpha }_{\tau }:\left\{\begin{array}{c}{y}^{\prime }=y\hfill \\ {\delta }^{\prime }=0\hfill \\ v_y^{\prime }=-\frac{v_y}{2}\hfill \end{array}\right..\hfill \end{array}$$

(17)

• In (17), there are six variables that are constrained by five equations. We might as well assume that $\{y,y^{\prime },v_y^{\prime },\delta ,{\delta }^{\prime }\}$ are dependent variables, and $\left\{v_y\right\}$ is the independent variable. The degree of $\left\{v_y\right\}$ is 2 and the size of the lattice array is (2+1). and, so, three implications are obtained (

  Table 2. Implications and R generated by Citing Instances Method (CIM) on the encoding initial condition.

  ${\mathit{v}}_{\mathit{y}}$	$\mathbf{Implication}$	R
  0	$a_3{\delta }^2+a_9\delta +a_{10}=0\vDash a_{10}=0$	$a_{10}$
  1	$a_3{\delta }^2+a_5\delta +a_9\delta +a_{10}+a_2+a_8=0\vDash a_{10}+\frac{1}{4}{a}_2-\frac{1}{2}{a}_8=0$	$a_{10}+\frac{1}{4}{a}_2-\frac{1}{2}{a}_8$
  2	$a_3{\delta }^2+2a_5\delta +a_9\delta +a_{10}+4a_2+2a_8=0\vDash a_{10}+a_2-a_8=0$	$a_{10}+a_2-a_8$

  ).
• Applying Algorithm 1. We get three simpler implications and three corresponding R by Lemma 3 (Table 2).

At last, we get three constraint equations $\{a_{10}=0,a_{10}+\frac{1}{4}{a}_2-\frac{1}{2}{a}_8=0,a_{10}+a_2-a_8=0\}$. Solving, we obtain $\{a_{10}=0,a_2=0,a_8=0\}$.

Step 4. Encode the Continuous Consecution

By Definition 11, the continuous consecution can be expressed by the implication

$$I\left(l\right)\wedge (\eta \left(l\right)=0)\vDash (\dot{\eta \left(l\right)}=0),$$

(18)

where $\dot{\eta \left(l\right)}$ is the Lie Derivative of $\eta \left(l\right)$ with respect to $F\left(l\right)$,

$$\dot{\eta \left(l\right)}=\frac{\partial \eta \left(l\right)}{\partial y}\dot{y}+\frac{\partial \eta \left(l\right)}{\partial \delta }\dot{\delta }+\frac{\partial \eta \left(l\right)}{\partial v_y}\dot{v_y};$$

that is,

$$\begin{array}{cc}\hfill \dot{\eta \left(\ell \right)}& =a_4{v}_y^2+2a_{1}y{v}_y+a_6\delta v_y+(-20a_2+a_5+a_7)v_y\hfill \\ & +(2a_3-10a_5)\delta +(-10a_4+a_6)y+(a_9-10a_8).\hfill \end{array}$$

(19)

In (18), there are three variables which are constrained by one equations (we ignore $I\left(l\right)=y\ge 0$). We might as well assume that $\left\{v_y\right\}$ is the dependent variable and that $\{y,\delta \}$ are independent variables. The degree of y and $\delta$ 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. Implications generated by CIM on the encoding continuous consecution.

$\mathit{\delta }$	y	Implication
0	0	$0=0$$\vDash a_4{v}_y^2+(a_5+a_7)v_y+a_9=0$
0	1	$\begin{array}{l}{a}_4{v}_y+a_1+a_7=0\hfill \\ \vDash a_4{v}_y^2+2a_1{v}_y+(a_5+a_7)v_y-10a_4+a_6+a_9=0\hfill \end{array}$
0	2	$\begin{array}{l}2a_4{v}_y+4a_1+2a_7=0\hfill \\ \vDash a_4{v}_y^2+4a_1{v}_y+(a_5+a_7)v_y-20a_4+2a6+a_9=0\hfill \end{array}$
1	0	$\begin{array}{l}{a}_5{v}_y+a_3+a_9=0\hfill \\ \vDash a_4{v}_y^2+a_6{v}_y+(a_5+a_7)v_y+2a_3-10a_5+a_9=0\hfill \end{array}$
1	1	$\begin{array}{l}{a}_4{v}_y+a_5{v}_y+a_1+a_3+a_6+a_7+a_9=0\hfill \\ \vDash a_4{v}_y^2+2a_1{v}_y+a_6{v}_y+(a_5+a_7)v_y+2a_3-10a_5-10a_4+a_6+a_9=0\hfill \end{array}$
1	2	$\begin{array}{l}2a_4{v}_y+a_5{v}_y+4a_1+a_3+2a_6+2a_7+a_9=0\hfill \\ \vDash a_4{v}_y^2+4a_1{v}_y+a_6{v}_y+(a_5+a_7)v_y+2a_3-10a_5-20a_4+2a_6+a_9=0\hfill \end{array}$
2	0	$\begin{array}{l}2a_5{v}_y+4a_3+2a_9=0\hfill \\ \vDash a_4{v}_y^2+2a_6{v}_y+(a_5+a_7)v_y+4a_3-20a_5+a_9=0\hfill \end{array}$
2	1	$\begin{array}{l}{a}_4{v}_y+2a_5{v}_y+a_1+4a_3+2a_6+a_7+2a_9=0\hfill \\ \vDash a_4{v}_y^2+2a_1{v}_y+2a_6{v}_y+(a_5+a_7)v_y+4a_3-20a_5-10a_4+a_6+a_9=0\hfill \end{array}$
2	2	$\begin{array}{l}2a_4{v}_y+2a_5{v}_y+4a_1+4a_3+4a_6+2a_7+2a_9=0\hfill \\ \vDash a_4{v}_y^2+4a_1{v}_y+2a_6{v}_y+(a_5+a_7)v_y+4a_3-20a_5-20a_4+2a_6+a_9=0\hfill \end{array}$

). We analyze some implications to illustrate the domino effect:

• $0=0\vDash a_4{v}_y^2+(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}_y^2+(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. Domino effect (1).

  $\mathit{\delta }$	y	$\mathbf{Implication}$
  0	0	$0=0$$\vDash 0=0$
  0	1	$a_1+a_7=0$$\vDash 2a_1{v}_y+a_6=0$
  0	2	$4a_1+2a_7=0$$\vDash 4a_1{v}_y+2a6=0$
  1	0	$-a_7{v}_y+a_3=0$$\vDash a_6{v}_y+2a_3+10a_7=0$
  1	1	$-a_7{v}_y+a_1+a_3+a_6+a_7=0$$\vDash 2a_1{v}_y+a_6{v}_y+2a_3+10a_7+a_6=0$
  1	2	$-a_7{v}_y+4a_1+a_3+2a_6+2a_7=0$$\vDash 4a_1{v}_y+a_6{v}_y+2a_3+10a_7+2a_6=0$
  2	0	$-2a_7{v}_y+4a_3=0$$\vDash 2a_6{v}_y+4a_3+20a_7=0$
  2	1	$-2a_7{v}_y+a_1+4a_3+2a_6+a_7=0$$\vDash 2a_1{v}_y+2a_6{v}_y+4a_3+20a_7+a_6=0$
  2	2	$-2a_7{v}_y+4a_1+4a_3+4a_6+2a_7=0$$\vDash 4a_1{v}_y+2a_6{v}_y+4a_3+20a_7+2a_6=0$

  ), we call it the domino effect, which is the benefit of parallelism in CIM.
• $-a_7{v}_y+2a_3=0\vDash a_6{v}_y+2a_3+10a_7=0$
  By Algorithm 1, $R=2a_3{a}_7+10a_7^2+2a_3{a}_6$ and $2a_3{a}_7+10a_7^2+2a_3{a}_6=0$ is the constraint equation. However, by comprehensive Gröbner basis, we need to compute the remainder under two situations: $a_7=0$ and $a_7\ne 0$. If $a_7$ is polynomial, situation becomes very complicated. This is why Sankaranarayanan et al. defined Alternative Consecution Relations to eliminate template variables in the antecedent of implication to avoid construction of comprehensive Gröbner basis. However, every coin has two sides. In Wu–Ritt’s Algorithm, $a_7=0$ means non-degenerate condition is false, the solution to $2a_3{a}_7+10a_7^2+2a_3{a}_6=0$ may or may not satisfy the formula (5), which is crucial to the theorem’s proof by CIM. 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.
• $a_1+a_7=0\vDash 2a_1{v}_y+a_6=0$
  By Algorithm 1, $a_1+a_7$ is a constant, so $R=2a_1{v}_y+a_6$ and $a_1=0,a_6=0$. We continue to simplify the implications(

  Table 5. Domino effect (2).

  $\mathit{\delta }$	y	$\mathbf{Implication}$
  0	0	$0=0$$\vDash 0=0$
  0	1	$a_7=0$$\vDash 0=0$
  0	2	$2a_7=0$$\vDash 0=0$
  1	0	$-a_7{v}_y+a_3=0$$\vDash 2a_3+10a_7=0$
  1	1	$-a_7{v}_y+a_3+a_7=0$$\vDash 2a_3+10a_7=0$
  1	2	$-a_7{v}_y+a_3+2a_7=0$$\vDash 2a_3+10a_7=0$
  2	0	$-2a_7{v}_y+4a_3=0$$\vDash 4a_3+20a_7=0$
  2	1	$-2a_7{v}_y+4a_3+a_7=0$$\vDash 4a_3+20a_7=0$
  2	2	$-2a_7{v}_y+4a_3+2a_7=0$$\vDash 4a_3+20a_7=0$

  ).
• $-a_7{v}_y+a_3=0\vDash 2a_3+10a_7=0$
  By Algorithm 1, $R=2a_3+10a_7$ and the constraint is $a_3=-5a_7$, continue to simplify the implications(

  Table 6. Domino effect (3).

  $\mathit{\delta }$	y	$\mathbf{Implication}$
  0	0	$0=0$$\vDash 0=0$
  0	1	$a_7=0$$\vDash 0=0$
  0	2	$2a_7=0$$\vDash 0=0$
  1	0	$-a_7{v}_y-5a_7=0$$\vDash 0=0$
  1	1	$-2a_7{v}_y-4a_7=0$$\vDash 0=0$
  1	2	$-a_7{v}_y-3a_7=0$$\vDash 0=0$
  2	0	$-2a_7{v}_y-20a_7=0$$\vDash 0=0$
  2	1	$-2a_7{v}_y-19a_7=0$$\vDash 0=0$
  2	2	$-2a_7{v}_y-18a_7=0$$\vDash 0=0$

  ).

At last, the final constraint equations in the template variables are

$$a_{\{1,2,4,6,8,9,10\}}=0,$$

$$a_5+a_7=0,$$

$$5a_7+a_3=0,$$

and the corresponding invariant is $y=v_y\delta +5{\delta }^2$ [3].

5. 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 $\{\eta \left({\ell }_i\right)=0\wedge {\rho }_{\tau }\wedge {\alpha }_{\tau }\vDash \eta {\left({\ell }_j\right)}^{\prime }=0\}$ and continuous consecution implication $\{I\left(\ell \right)\wedge \eta \left(\ell \right)=0\vDash \dot{\eta \left(\ell \right)}=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. Sankaranarayanan et al. defined as Alternative Consecution Relations (as shown in

Table 7. Alternative discrete consecution conditions. LC: local consecution; CV: constant-value consecution; CS: constant-scale consecution; PS: polynomial-scale consecution.

Name	Alternative Discrete Consecution Condition
LC	${\alpha }_{\tau }\wedge {\rho }_{\tau }\vDash \eta {\left(l_j\right)}^{\prime }=0$
CV	${\alpha }_{\tau }\wedge {\rho }_{\tau }\vDash (\eta {\left(l_j\right)}^{\prime }-\eta \left(l_i\right))=0$
CS	${\alpha }_{\tau }\wedge {\rho }_{\tau }\vDash (\eta {\left(l_j\right)}^{\prime }-\lambda \eta \left(l_i\right))=0$
PS	${\alpha }_{\tau }\wedge {\rho }_{\tau }\vDash (\eta {\left(l_j\right)}^{\prime }-p\eta \left(l_i\right))=0$

and

Table 8. Alternative Continuous Consecution Conditions. LC: local consecution; CV: constant-value consecution; CS: constant-scale consecution; PS: polynomial-scale consecution.

Name	Alternative Continue Consecution Condition
LC	$I\left(l\right)\vDash \dot{\eta \left(l\right)}=0$
CV	$I\left(l\right)\vDash (\dot{\eta \left(l\right)}-\eta \left(l\right))=0$
CS	$I\left(l\right)\vDash (\dot{\eta \left(l\right)}-\lambda \eta \left(l\right))=0$
PS	$I\left(l\right)\vDash (\dot{\eta \left(l\right)}-p\eta \left(l\right))=0$

) 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^{\prime }$ 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.

$${\rho }_{\tau }\wedge {\alpha }_{\tau }\vDash \eta {\left(l\right)}^{\prime }-\lambda \eta \left(l\right)=0,$$

(20)

$$I\left(l\right)\vDash \dot{\eta \left(l\right)}-\lambda \eta \left(l\right)=0,$$

(21)

$$\begin{array}{cc}\hfill & \underset{m=1}{\overset{k}{\bigwedge }}{f}_m(u_1,\dots ,u_n,x_1,\dots ,x_k)=0\vDash \hfill \\ \hfill & G(a_1,\dots ,a_{\ell },u_1,\dots ,u_n,x_1,\dots ,x_k)=0.\hfill \end{array}$$

(22)

We apply Algorithm 1 to (20) and (21).

• For the initial condition case, it is same as Step 1.
• For the discrete consecution case (with $\lambda d$), assume that $\{y,y^{\prime },v_y^{\prime },{\delta }^{\prime }\}$ are dependent variables, $\{v_y,\delta \}$ are independent variables, and the size of lattice array is (2+1, 2+1).
• For the continuous consecution case (with $\lambda c$), there are no dependent variables as $I\left(l\right)\ge 0$ is ignored, $\{v_y,\delta ,y\}$ are independent variables, and the size of lattice array is (2+1, 2+1,2+1).
• Solving these 37 implications, the following two groups of constraints are obtained:

$$\begin{array}{cc}\hfill & a_1=0,a_2=0,a_3=0,a_4=0,a_5=0,\hfill \\ \hfill & a_6=0,a_7=0,a_8=0,a_9=0,a_{10}=0,\hfill \\ \hfill & \lambda d=\lambda d,\lambda c=\lambda c;\hfill \end{array}$$

(23)

$$\begin{array}{cc}\hfill & a_1=0,a_2=0,a_3=-5a_7,a_4=0,a_5=-a_7,\hfill \\ \hfill & a_6=0,a_7=a_7,a_8=0,a_9=0,a_{10}=0,\hfill \\ \hfill & \lambda d=0,\lambda c=0.\hfill \end{array}$$

(24)

Obviously, (23) is a trivial solution, (24) is the same as in [3], and the corresponding invariant is $y=v_y\delta +5{\delta }^2$.

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.

6. 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.

6.1. 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 ${\tau }_1,{\tau }_2,$ and ${\tau }_3$; the transition relations are as follows:

$$\begin{array}{cc}{\rho }_{{\tau }_1}\wedge {\alpha }_{{\tau }_1}:\hfill & v=5\wedge id(s,x,v,t),\hfill \\ {\rho }_{{\tau }_2}\wedge {\alpha }_{{\tau }_2}:\hfill & id(s,x,v,t),\hfill \\ {\rho }_{{\tau }_3}\wedge {\alpha }_{{\tau }_3}:\hfill & \left[v=0\wedge s^{\prime }=s+1\wedge t^{\prime }=t+2\wedge id(x,v)\right],\hfill \end{array}$$

(25)

where $id\left(x\right)$ denotes $x^{\prime }=x$.

$\eta \left(l_i\right),i=0,1,2$ is the template map,

$$\eta \left(l_0\right):\left[\begin{array}{c}{a}_1{x}^2+a_{2}xs+a_{3}xt+a_{4}xv\hfill \\ +a_5{s}^2+a_{6}st+a_{7}sv\hfill \\ +a_8{t}^2+a_{10}{v}^2+a_{11}x\hfill \\ +a_{12}s+a_{13}t+a_{14}v+a_{15}\hfill \end{array}\right],$$

$$\eta \left(l_1\right):\left[\begin{array}{c}{b}_1{x}^2+b_{2}xs+b_{3}xt+b_{4}xv\hfill \\ +b_5{s}^2+b_{6}st+b_{7}sv\hfill \\ +b_8{t}^2+b_{10}{v}^2+b_{11}x\hfill \\ +b_{12}s+b_{13}t+b_{14}v+b_{15}\hfill \end{array}\right],$$

$$\eta \left(l_2\right):\left[\begin{array}{c}{c}_1{x}^2+c_{2}xs+c_{3}xt+c_{4}xv\hfill \\ +c_5{s}^2+c_{6}st+c_{7}sv\hfill \\ +c_8{t}^2+c_{10}{v}^2+c_{11}x\hfill \\ +c_{12}s+c_{13}t+c_{14}v+c_{15}\hfill \end{array}\right],$$

Applying Algorithm 1, 2250 instances in total and the following constraints are obtained:

$$\begin{array}{c}{a}_{\{1,\dots ,9,15\}}=0,\hfill \\ a_{10}=-\frac{1}{4}{b}_{11}\lambda d_{12}\lambda d_{20}\hfill \\ a_{11}=b_{11}\lambda d_{12}\lambda d_{20}\hfill \\ a_{12}=\frac{115}{4}{b}_{11}\lambda d_{12}\lambda d_{20}\hfill \\ b_{\{1,\dots ,9\}}=0\hfill \\ b_{10}=\frac{1}{2}{b}_{11},b_{11}=b_{11}\hfill \\ b_{12}=\frac{115}{4}{b}_{11},b_{13}=-5b_{11}\hfill \\ b_{14}=-5b_{11},b_{15}=-\frac{75}{4}{b}_{11},\hfill \\ c_{\{1,\dots ,9\}}=0\hfill \\ c_{10}=\frac{1}{2}{b}_{11}\lambda d_{12}\hfill \\ c_{11}=b_{11}\lambda d_{12}\hfill \\ c_{12}=\frac{115}{4}{b}_{11}\lambda d_{12},\hfill \\ c_{13}=-5b_{11}\lambda d_{12},\hfill \\ c_{14}=-5b_{11}\lambda d_{12},\hfill \\ c_{15}=-\frac{75}{4}{b}_{11}\lambda d_{12},\hfill \\ \lambda c_0=0,\lambda c_1=0,\lambda c_2=0,\hfill \\ \lambda d_{01}=\frac{1}{\lambda d_{12}\lambda d_{20}},\hfill \\ \lambda d_{12}=\lambda d_{12},\hfill \\ \lambda d_{20}=\lambda d_{20},\hfill \end{array}$$

(26)

where $\lambda c_i$ denotes continuous consecution $\lambda$ at location $l_i$, and $\lambda d_{ij}$ denotes discrete consecution $\lambda$ for translation from $l_i$ to $l_j$.

When $b_{11}=1,\lambda d_{12}=\lambda d_{20}=1$, the invariant is obtained as follows:

$$\begin{array}{cc}\eta \left(l_0\right):\hfill & v^2-4x-10v+115s+20t=0,\hfill \\ \eta \left(l_1\right):\hfill & 5v^2+4xv+115vs-20vt=0,\hfill \\ \eta \left(l_2\right):\hfill & 2v^2+4x-20v+115s-20t+75=0,\hfill \end{array}$$

(27)

which the same as that derived in [3].

6.2. 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\le d\le 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:

$$\begin{array}{c}left:\left\{\begin{array}{c}{v}_y+2=0\hfill \\ v_x-2=0\hfill \\ a(x+y)=4b-4ab\hfill \end{array}\right.\hfill \\ magnetic:\left\{\begin{array}{c}{v}_y+2=a(x-2)\hfill \\ v_x^2+v_y^2=8\hfill \\ v_x-2=-a(y+2)+4b(1-a)\hfill \end{array}\right.\hfill \\ right:\left\{\begin{array}{c}{v}_y+2=0\hfill \\ v_x+2=0\hfill \\ a(x-y)=-4(b+1)(1-a)\hfill \end{array}\right.\hfill \end{array}$$

which are the same as those obtained in [3].

7. 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\equiv 0$, according to Lemma 2. More specifically, for each $\widehat{\mathit{u}}\in S_R$, first substitute it for $\mathit{u}$ in (11) (this also means that free variables are eliminated). Then, compute the roots $\mathit{r}(a_1,\cdots ,a_n)$ of $\mathit{x}$ in $\mathit{a}$ and, finally, substitute $\mathit{r}(a_1,\cdots ,a_n)$ for $\mathit{x}$ in (12). Now, we have a polynomial constraint $G(a_1,\cdots ,a_n)=0$ in $\mathit{a}$. Solving $G(a_1,\cdots ,a_n)=0$ generates a set of solutions to $\mathit{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=\frac{i(i+1)}{2}$.

The target template invariant is a generic degree-two template polynomial: $\eta =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

$$\eta =0\wedge i^{\prime }-i-1=0\wedge s^{\prime }-s-i^{\prime }=0\wedge t=1\wedge t^{\prime }=0\vDash {\eta }^{\prime }=0,$$

assuming $\{s,i^{\prime },s^{\prime },t,t^{\prime }\}$ are the dependent variables and $\left\{i\right\}$ is the independent variable. By CIM, we need compute the roots of $\{s,i^{\prime },s^{\prime },t,t^{\prime }\}$ 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^{\prime },i^{\prime },t,t^{\prime }\}$ 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^{\prime },s^{\prime },t,t^{\prime }\}$ in the log file (

Table 9. The experimentally collected hybrid states of $\{i,s,i^{\prime },s^{\prime },t,t^{\prime }\}$.

i	s	${\mathbf{i}}^{\prime }$	${\mathbf{s}}^{\prime }$	t	${\mathbf{t}}^{\prime }$
0	0	1	1	1	0
1	1	2	3	1	0
2	3	3	6	1	0
3	6	4	10	1	0

) during every discrete transition.

Finally, we obtain the invariant $s=\frac{i(i+1)}{2}$. Without directly computing the roots, CIM can dramatically reduce the complexity. Similar ideas appeared in [43,44].

8. 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, the authors are confident that most of the symbolic algebra used in existing methods can be replaced using this technique. Comparing to the well-established approaches in this field, CIM exhibits the following features which make it interesting:

• Applying instances by essentially instantiating the free variables to real numbers. Hence, the free variables are removed.
• Comparing to the existing approaches, CIM can solve the complete inductive conditions directly and is a complete approach.
• According to our analysis, Wu–Ritt’s method has exponential complexity while Gröbner basis method has double exponential complexity. Therefore, our approach is more efficient.
• The main idea of CIM is “turning the difficulty of quality into the complexity of quantity”. If Formula (5) is too difficult to deal with, CIM turns (5) into $|S_R|$ simpler implications without free variables, then CIM is meaningful.
• Parallelization is another advantage of CIM. CIM was created to be a parallel method. To raise the computer’s calculation capacity, one important method is to develop both parallel machines and parallel algorithms.

We also wish to extend CIM to the inequality field. In fact, the famous algorithm of cylindrical algebraic decomposition (CAD) [45] is a typical algorithm of proving by instances.

Many algorithms rely on the computation of (comprehensive) Gröbner basis [19]. If the whole task is too difficult to be accomplished using a Gröbner basis in one stroke, CIM can solve it by using many instances. From another perspective, CIM can be regarded as a Divide and Conquer (DAC)-type method. In Example 3, as it is too difficult to encode continuous consecution using a comprehensive Gröbner basis, we can first use CIM to obtain nine simpler implications and then apply (comprehensive) Gröbner basis, one-by-one. The combination of the CIM and (comprehensive) Gröbner basis is another interesting research goal.