Invariant Clusters for Hybrid Systems
Hui Kong∗, Sergiy Bogomolov∗, Christian Schilling† , Yu Jiang‡, Thomas A. Henzinger∗
∗ IST

Austria, Klosterneuburg, Austria
of Freiburg, Freiburg, Germany
‡ University of Illinois at Urbana-Champaign, Illinois, USA

arXiv:1605.01450v1 [math.OC] 4 May 2016

† University

Abstract—In this paper, we propose an approach to automatically compute invariant clusters for semialgebraic hybrid systems.
An invariant cluster for an ordinary differential equation (ODE)
is a multivariate polynomial invariant g(~u, ~x) = 0, parametric
in ~u, which can yield an infinite number of concrete invariants
by assigning different values to ~u so that every trajectory of the
system can be overapproximated precisely by a union of concrete
invariants. For semialgebraic systems, which involve ODEs with
multivariate polynomial vector flow, invariant clusters can be
obtained by first computing the remainder of the Lie derivative
of a template multivariate polynomial w.r.t. its Gröbner basis and
then solving the system of polynomial equations obtained from
the coefficients of the remainder. Based on invariant clusters and
sum-of-squares (SOS) programming, we present a new method
for the safety verification of hybrid systems. Experiments on
nonlinear benchmark systems from biology and control theory
show that our approach is effective and efficient.
Index Terms—ybrid system, nonlinear system, semialgebraic
system, invariant, safety verification, SOS programmingybrid
system, nonlinear system, semialgebraic system, invariant, safety
verification, SOS programmingh

I. I NTRODUCTION
Hybrid systems [1] are models for systems with interacting
discrete and continuous dynamics. Safety verification is among
the most challenging problems in verifying hybrid systems,
asking whether a set of bad states can be reached from a set
of initial states. The safety verification problem for systems
described by nonlinear differential equations is particularly
complicated because computing the exact reachable set is
usually intractable. Existing approaches are mainly based
on approximate reachable set computations [2], [3], [4] and
abstraction [5], [6], [7], [8], [9].
An invariant is a special kind of overapproximation for the
reachable set of a system. Since invariants do not involve
direct computation of the reachable set, they are especially
suitable for dealing with nonlinear hybrid systems. However,
automatically and efficiently generating sufficiently strong
invariants is challenging on its own.
In this work, we propose an approach to automatically compute invariant clusters for a class of nonlinear semialgebraic
systems whose trajectories are algebraic, i.e., every trajectory
of the system is essentially a common zero set (algebraic
variety) of a set of multivariate polynomial equations. An
invariant cluster for a semialgebraic system is a parameterized
multivariate polynomial invariant g(~u, ~x) = 0, with parameter
~u, which can yield an infinite number of concrete invariants
by assigning different values to ~u so that every trajectory of
the system can be overapproximated precisely by a union of
concrete invariants.

The basic idea of computing invariant clusters is as follows.
A sufficient condition for a trajectory of a semialgebraic
system to start from and to always stay in the solution set
of a multivariate polynomial equation g(~x) = 0 is that the
Lie derivative of g(~x) w.r.t. f~ belongs to the ideal generated
by g(~x), i.e., Lf~g ∈ hg(~x)i, where f~ is the vector flow of
the system. Therefore, if some g(~x) satisfies this condition,
g(~x) = 0 is an invariant of the system. Then, according to
Gröbner basis theory, Lf~g ∈ hg(~x)i implies that the remainder
of Lf~g w.r.t. hg(~x)i must be identical to 0. Based on this
theory, we first set up a template polynomial g(~u, ~x) and then
compute the remainder r(~u, ~x) of Lf~g w.r.t. hg(~u, ~x)i. Since
r(~u, ~x) ≡ 0 implies that all coefficients ai (~u) of ~x in r(~u, ~x)
are equal to zero, we can set up a system P of polynomial
equations on ~u from the coefficients ai (~u). By solving P we
get a set C of constraints on ~u. For those elements in C that
are linear in ~u, the corresponding parameterized polynomial
equations g(~u, ~x) = 0 form invariant clusters. Based on
invariant clusters and SOS programming, we propose a new
method for the safety verification of hybrid systems.
The main contributions of this paper are as follows: 1)
We propose to generate invariant clusters for semialgebraic
systems based on computing the remainder of the Lie derivative of a template polynomial w.r.t. its Gröbner basis and
solving the system of polynomial equations obtained from the
coefficients of the remainder. Our approach avoids Gröbner
basis computation and first-order quantifier elimination. 2) We
present a method to overapproximate trajectories precisely by
using invariant clusters. 3) We apply invariant clusters to the
safety verification of semialgebraic hybrid systems based on
SOS programming. 4) We implemented a prototype tool to
perform the aforementioned steps automatically. Experiments
show that our approach is effective and efficient.
The paper is organized as follows. Section II is devoted to
the preliminaries. In Section III, we introduce the approach
to computing invariant clusters and using them to characterize
trajectories. In Section IV, we present a method to verify safety
properties for semialgebraic continuous and hybrid systems
based on invariant clusters. In Section V, we present our
experimental results. In Section VI, we introduce some related
works. Finally, we conclude our paper in Section VII.
II. P RELIMINARIES
In this section, we recall some backgrounds we need
throughout the paper. We first clarify some notation conventions. If not specified otherwise, we decorate vectors ~·, we

use the symbol K for a field, R for the real number field, C
for the complex number field (which is algebraically closed)
and N for the set of natural numbers, and all the polynomials
involved are multivariate polynomials. In addition, for all the
polynomials g(~u, ~x), we denote by ~u the vector composed of
all the ui and denote by ~x the vector composed of all the
remaining variables that occur in the polynomial.
Definition 1 (Ideal): [10] A subset I of K[~x], is called an
ideal if 1) 0 ∈ I; 2) if p, q ∈ I, then p + q ∈ I; and 3) if
p ∈ I and q ∈ K[~x], then pq ∈ I.
Definition 2 (Generated ideal): [10] Let g1 , . . . , gs be
polynomials in K[~x]. The ideal generated by {g1 , . . . , gs }
is
s
X
def 
hg1 , . . . , gs i =
hi gi | h1 , . . . , hs ∈ K[~x] .
i=1

Definition 3 (Algebraic variety): Let K be an algebraically
closed field and I ⊂ K[~x] be an ideal. We define the algebraic
variety of I as
def

V(I) = {~x ∈ Kn | f (~x) = 0 for f ∈ I}.
Next, we present the notation of the Lie derivative, which
is widely used in the discipline of differential geometry.
Definition 4 (Lie derivative): For a given polynomial p ∈
K[~x] over ~x = (x1 , · · · , xn ) and a continuous system ~x˙ = f~,
where f~ = (f1 , . . . , fn ), the Lie derivative of p ∈ K[~x] along
f of order k is defined as follows.
(
p,
k=0
k def
(1)
Lf~p =
p
Pn ∂Lfk−1
~
· fi , k ≥ 1
i=1
∂xi

Essentially, the k-th order Lie derivative of p is the k-th
derivative of p w.r.t. time, i.e., reflects the change of p over
time. We write Lf~p for L1f~p.
We furthermore use the following theorem for deciding
the existence of a real solution of a system of polynomial
constraints.
Theorem 1 (Real Nullstellensatz): [11] The system of multivariate polynomial equations and inequalities p1 (~x) = 0, . . . ,
pm1 (~x) = 0, q1 (~x) ≥ 0, . . . , qm2 (~x) ≥ 0, r1 (~x) > 0, . . . ,
rm3 (~x) > 0 either has a solution in Rn , or there exists the
following polynomial identity
m1
X

βi p i +

v∈{0,1}m2

i=1

+

X

X

v∈{0,1}m3

X
t

ctv

X
t

2

·

m3
m2
Y
2 Y
vj
qj +
rkdk
btv ·
j=1

m3
Y

k=1

rkvk

+

k=1

X

s2w

(2)

=0

w

where dk ∈ N and pi , qj , rk , βj , btv , ctv , sw are polynomials
in R[~x].
Remark 1: Theorem 1 enables us to efficiently decide if
a system of polynomial equations
P and inequalities has a real
solution. By moving the term w s2w in equation (2) to the
right-hand side and denoting
P 2the remaining terms by R(~x, ~y),
x, ~y )
we have −R(~x, ~y) =
w sw , which means that −R(~
is a sum-of-squares. Therefore, finding a set of polynomials

βj , rk , btv , ctv , sw as well as some dk ’s that make −R(~x, ~y )
a sum-of-squares is sufficient to prove that the system has
no real solution, which can be done efficiently by SOS
programming [12].
In this paper, we focus on semialgebraic continuous and hybrid systems which are defined in the following, respectively.
Definition 5 (Semialgebraic system): A semialgebraic sysdef
tem is a tuple M = hX, f~, X0 i, where
1) X is the state space of the system M ,
2) f~ is a Lipschitz continuous vector flow function, and
3) X0 is the initial set described by a semialgebraic set.
The Lipschitz continuity guarantees existence and uniqueness of the differential equation ~x˙ = f~. A trajectory of a
semialgebraic system is defined as follows.
Definition 6 (Trajectory): Let M be a semialgebraic system. A trajectory originating from a point ~x0 ∈ X0 to
time T > 0 is a continuous and differentiable function
~x(t) : [0, T ) → Rn such that i) ~x(0) = ~x0 , and ii)
x
~ x(τ )).
∀τ ∈ [0, T ) : d~
dt t=τ = f (~
Definition 7 (Safety): Given an unsafe set XU , a semialgebraic system M = hX, f~, X0 i is said to be safe if no trajectory
~x(t) of M satisfies both ~x(0) ∈ X0 and ∃τ ∈ R≥0 : ~x(τ ) ∈
XU .
Definition 8 (Hybrid System): A hybrid system is dedef
scribed by a tuple H = hL, X, E, R, G, I, F i, where
• L is a finite set of locations (or modes);
n
• X ⊆ R is the continuous state space. The hybrid state
space of the system is denoted by X = L × X and a state
is denoted by (l, ~x) ∈ X ;
• E ⊆ L × L is a set of discrete transitions;
X
• G : E → 2 is a guard mapping over discrete transitions;
X
• R : E × 2
→ 2X is a reset mapping over discrete
transitions;
X
• I :L →2
is an invariant mapping;
• F : L → (X → X) is a vector field mapping which
assigns to each location l a vector field f~.
The transition and dynamic structure of the hybrid system
defines a set of trajectories. A trajectory is a sequence originating from a state (l0 , ~x0 ) ∈ X0 , where X0 ⊆ X is an initial set,
and consisting of a series of interleaved continuous flows and
discrete transitions. During the continuous flows, the system
evolves following the vector field F (l) at some location l ∈ L
as long as the invariant condition I(l) is not violated. At some
state (l, ~x), if there is a discrete transition (l, l′ ) ∈ E such
that (l, ~x) ∈ G(l, l′ ) (we write G(l, l′ ) for G((l, l′ ))), then
the discrete transition can be taken and the system state can
be reset to R(l, l′ , ~x). The problem of safety verification of a
hybrid system is to prove that the hybrid system cannot reach
an unsafe set XU from an initial set X0 .
III. C OMPUTATION

OF I NVARIANT

C LUSTERS

In this section, we first introduce the notion of invariant
cluster and then show how to compute a set of invariant
clusters and how to use it to represent every trajectory of a
semialgebraic system.

A. Foundations of invariants and invariant clusters
Given a semialgebraic system M , for any trajectory ~x(t)
of M , if we can find a multivariate polynomial g(~x) ∈ R[~x]
such that g(~x(0)) ∼ 0 implies g(~x(t)) ∼ 0 for all t > 0,
where ∼ ∈ {<, ≤, =, ≥, >}, then g(~x) ∼ 0 is an invariant of
the system. For the convenience of presentation, we call g(~x)
an invariant polynomial of M . In addition, a trajectory ~x(t) is
said to be algebraic if there exists an invariant g(~x) = 0 for
~x(t) and g(~x) 6≡ 0. In the following, we present a sufficient
condition for a polynomial g(~x) to be an invariant polynomial.
Proposition 1: Let M = hX, f~, X0 i be a semialgebraic
system and g(~x) ∈ R[~x]. Then g ∼ 0 is an invariant of M for
every ∼ ∈ {<, ≤, =, ≥, >} if g(~x) satisfies
Lf~g ∈ hgi

Fig. 1: Example 2. Curve defined by invariant cluster of
Class(C ∗ , ~x0 ) for ~x0 = (4, 2).

(3)

Proof: See Appendix A.
Proposition 1 states that all the polynomial equations and
inequalities g ∼ 0 are invariants of M if the Lie derivative of
g belongs to the ideal hgi. Note that every invariant satisfying
condition (3) defines a closed set for trajectories, that is, no
trajectory can enter or leave the set defined by the invariant.
For a semialgebraic system whose trajectories are algebraic,
the trajectories can usually be divided into several groups and
in each group all trajectories show similar curves. Essentially,
these similar curves can be described identically by a unique
parameterized polynomial equation, which we characterize
as an invariant cluster. The computation method of invariant
clusters is presented in Subsection III-B.
Definition 9 (Invariant cluster): An invariant cluster C of
a semialgebraic system is a set of invariants which can be uniK
~
formly described
PM as C =i {g(~u, ~x) = 0 | ~u ∈ R \{0}}, where
g(~u, ~x) = i=1 ci (~u)X satisfies Lf~g ∈ hgi and ci (~u) ∈ R[~u]
are fixed linear polynomials on ~u = (u1 , · · · , uK ), X i are
monomials on ~x = (x1 , · · · , xn ), and M, K ∈ N.
Note that by requiring ~u 6= ~0 in Definition 9 and the
following related definitions, we exclude the trivial invariant
0 = 0. Given an invariant cluster, by varying the parameter
~u, we may obtain an infinite set of concrete invariants for
the system. To be intuitive, we present a running example to
demonstrate the related concepts throughout the paper.
Example 1 (running example): Consider

the semialgebraic
system M1 described by [ẋ, ẏ] = y 2 , xy . The set C ∗ =
{u1 − u3 (x2 − y 2 ) = 0 | (u1 , u3 ) ∈ R2 \ {~0}} is an invariant
cluster. It is easy to verify that the polynomial u1 −u3 (x2 −y 2 )
satisfies condition (3) for all (u1 , u3 ) ∈ R2 .

Next we introduce the notion of an invariant class.
Definition 10 (Invariant class): Given a semialgebraic system M with some initial point ~x0 and an invariant cluster
C = {g(~u, ~x) = 0 | ~u ∈ RK \ {~0}} of M , where K ∈ N, an
invariant class of C at ~x0 , denoted by Class(C, ~x0 ), is the
set {g(~u, ~x) = 0 | g(~u, ~x0 ) = 0, ~u ∈ RK \ {~0}}.
Given an invariant cluster C, by substituting a specific point
~x0 for ~x in g(~u, ~x) = 0, we obtain a constraint g(~u, ~x0 ) = 0 on
~u, which yields a subset Class(C, ~x0 ) of C. Apparently, every

Fig. 2: Example 6. X0 : (x + 15)2 + (y − 17)2 ≤ 1, XU :
(x1 − 11)2 + (y1 − 16.5)2 ≤ 1.

member of Class(C, ~x0 ) is an invariant for the trajectory
originating from ~x0 .
Example 2 (running example): For the given invariant cluster C ∗ in Example 1 and a given initial point ~x0 = (4, 2),
we get the invariant class Class(C ∗ , ~x0 ) = {u1 − u3 (x2 −
y 2 ) = 0 | u1 − 12u3 = 0, (u1 , u3 ) ∈ R2 \ {~0}}. Every
member of Class(C ∗ , ~x0 ) is an invariant for the trajectory
of M1 originating from ~x0 . The algebraic variety defined by
Class(C ∗ , ~x0 ) is shown in Figure 1.
An invariant class has the following important properties.
Theorem 2: Given an n-dimensional semialgebraic system
M and an invariant class D = {g(~u, ~x) = 0 | g(~u, ~x0 ) = 0,
~u ∈ RK \ {~0}} of M at a specified point ~x0 , let Dg be the
set of all invariant polynomials occurring in D and π~x0 be the
trajectory of M originating from ~x0 . Then,
1) π~x0 ⊆ V(Dg );
2) there must exist a subset B of Dg consisting of m members such that hDg i = hBi, where m is the dimension
of the hyperplane g(~u, ~x0 ) = 0 in terms of ~u in RK and
hDg i and hBi denote the ideal generated by the members
of Dg and B, respectively.
Proof: See Appendix B.
Remark 2: The first property in Theorem 2 reveals that a
trajectory π~x0 is always contained in the intersection of all the
invariants in the invariant class D of ~x0 . The second property
concludes that the invariant class can be generated by a finite
subset B of D if it consists of an infinite number of invariants.

The algebraic variety V(B) (which is equivalent to V(Dg ))
forms an overapproximation for π~x0 and the quality of the
overapproximation depends largely on the dimension m of
V(B). More specifically, the lower the dimension is, the better
the overapproximation is. The ideal case is that, when m = 1,
V(B) shrinks to an algebraic curve and hence some part of
the algebraic curve matches the trajectory precisely. In the
case of m > 1, V(B) is usually a hypersurface. To make the
overapproximation less conservative, we may take the union
of multiple invariant classes from different invariant clusters
to reduce the dimension of the algebraic variety.

Algorithm 1: Computation of invariant clusters
input : f : n-dimensional polynomial vector field;
N : upper bound for the degree of invariant
polynomials
output: CFamily: a set of invariant clusters
1
2
3

4
5
6

B. Computing invariant clusters

7
8

In the following, we propose a computation method for
invariant clusters. The main idea arises from the following
proposition.
Proposition 2: Let M = hX, f~, X0 i be a semialgebraic
system and suppose for a given multivariate polynomial g ∈
R[~x] that the polynomial Lf~g can be written as Lf~g = pg + r.
Then Lf~g ∈ hgi if r ≡ 0.
According to Proposition 2, if we can find a polynomial
g(~x) such that the remainder of its Lie derivative w.r.t. g(~x)
is identical to 0, then g(~x) = 0 is an invariant of M . The idea
is as follows. We first establish a template g(~u, ~x) for g(~x)
with parameter ~u and then compute the remainder r(~u, ~x) of
Lf~g w.r.t. g(~u, ~x). According to the computing procedure of
PK
a remainder [10], r(~u, ~x) must be of the form i=1 biu(~du) X i ,
j
where d ∈ N, bi (~u) are homogeneous polynomials of degree
d+1 over ~u, uj is the coefficient of the leading term of g(~u, ~x)
by some specified monomial order of ~x, and X i are monomials
on ~x. Since r(~u, ~x) ≡ 0 implies uj 6= 0 and bi (~u) = 0
for all i = 1, . . . , K, we obtain a system C of homogeneous
polynomial equations on ~u plus uj 6= 0 from the coefficients
of r(~u, ~x). Solving C can provide a set of invariant clusters of
M if it exists. Note that all the aforementioned steps can be
performed automatically in a mathematical software such as
Maple. Pseudocode for computing invariant clusters is shown
in Algorithm 1. The principle for the loop in line 5 is that the
remainder may vary from the monomial order of ~x and hence
produce different solutions. Using multiple orders helps to get
more solutions.
Remark 3 (Complexity): In Algorithm 1, the key steps are
computing the remainder in line 4 and solving the system of
equations on ~u in line 8. The former takes only linear time
and hence is very efficient. The latter involves solving a system
of homogeneous polynomial equations, which is known to be
NP-complete [13]. In our implementation in Maple, we use the
command solve for solving the system of equations. In our
experiments on nonlinear (parametric) systems of dimension
ranging from 2 to 8, the solver can return quickly in most cases
whether they have solutions or not. Among the solutions we
obtained, we found some complex solutions, however, we have
not seen nonlinear solutions.
Example 3 (running example): According to Algorithm 1,

9
10
11

CFamily ← ∅;
for i ← 1 to N do
gu~ ,~x ← generate parameterized polynomial over ~x of
degree i;
Lf g ← compute the Lie derivative of gu~ ,~x ;
foreach monomial order O of ~x do
Ru~ ,~x ← compute remainder of Lf g w.r.t. g by O;
Coeffs ← collect coefficients of ~x in Ru~ ,~x ;
Solution ← solve system Coeffs on ~u;
CFamily ← CFamily ∪ Solution;
end
end

the steps for computing the invariant clusters of degree 2 are
as follows:
1) Generate the template polynomial of degree 2:
g2 (~u, ~x) = u6 x2 + u5 xy + u4 x + u3 y 2 + u2 y + u1
2) Compute the Lie derivative Lf~g2 using Definition 4:
Lf~g2 =

∂g2
∂g2
ẋ +
ẏ = u5 x2 y + (2 u3 + 2 u6 ) xy 2
∂x
∂y
+ u2 xy + y 3 u5 + u4 y 2

3) Compute the remainder of Lf~g2 w.r.t. g2 by graded
reverse lexicographic (grevlex) order of (x, y). Using this
order, the leading term of Lf~g2 and g2 is u5 x2 y and
u6 x2 , respectively. Then:

2 u3 u6 − u25 + 2 u26 xy 2
u5 y
r(~u, ~x) = Lf~g2 −
g2 =
u6
u6
(u2 u6 − u4 u5 ) xy
(−u3 u5 + u5 u6 ) y 3
+
+
u6
u6
u1 u5 y
(−u2 u5 + u4 u6 ) y 2
−
+
u6
u6
4) Collect the coefficients of r(~u, ~x):

u2 u6 − u4 u5 −u3 u5 + u5 u6 −u2 u5 + u4 u6
S :=
,
,
,
u6
u6
u6

u1 u5
2 u3 u6 − u25 + 2 u26
,−
u6
u6
5) Solve the system formed by S. To save space, we just
present one of the six solutions we obtained:
C6 = {u6 = −u3 , u2 = u4 = u5 = 0, u3 = u3 , u1 = u1 }
6) Substitute the above solution C6 for ~u in g2 (~u, ~x). We
get the following parameterized invariant polynomial:
g2 (~u, ~x) = −u3 x2 + u3 y 2 + u1

Note that the other five solutions obtained in step 5) are in
fact the products of the invariant polynomials {u2 y, u1 (x +
y), u1 (x − y)}, which have been obtained when initially
computing the invariants of degree 1. Hence they cannot
increase the expressive power of the set of invariant clusters
and should be dropped. The solution presented above is the
one we have given in Example 1.

Algorithm 2: Computation of invariant classes
input : CFamily: set of invariant clusters;
~x0 : an initial point
output: ICls: list of invariant classes
1
2
3
4

C. Overapproximating trajectories by invariant classes
In this section, we address how to overapproximate trajectories precisely by using invariant classes.
Invariant clusters can be divided into two categories according to the number of invariant classes that they can yield by
varying the parameter ~u. 1) finite invariant cluster: This kind
of invariant cluster can yield only one invariant class no matter
how ~u changes. For example, {u1 (x − y) = 0 | u1 ∈ R \ {0}}
is such an invariant cluster for the running example. In this
case, the trajectories that can be covered by the invariant
class is very limited. Moreover, the overapproximation is also
conservative due to the high dimension of the algebraic variety
defined by the invariant class. 2) infinite invariant cluster:
One such invariant cluster C can yield an infinite number of
invariant classes Class(C, ~x0 ) as the initial point ~x0 varies,
e.g., the invariant cluster C ∗ in Example 1. For the trajectory π~x0 , the overapproximating precision of Class(C, ~x0 )
depends largely on the dimension m of the algebraic variety
defined by Class(C, ~x0 ). The best case is m = 1 and then it
gives a curve-to-curve match in part for the trajectory.
Now, we introduce how to identify the invariant classes for
a given point ~x0 from a set of invariant clusters and how to get
a finite representation for it. To be intuitive, we first present a
3-dimensional system and a set of invariant clusters for it.
Example 4: Consider the following semialgebraic system
M2 : [ẋ, ẏ, ż] = [yz, xz, xy]. We obtain a set of invariant
clusters consisting of 7 elements. Here we only present the
infinite invariant cluster (for other clusters see Appendix I).
C7 = {g7 (~u, ~x) = (−u5 − u6 )x2 + u5 y 2 + u6 z 2 + u0
= 0 | ~u ∈ R3 \ {~0}}
The invariant clusters are capable of overapproximating all the
trajectories of the system M2 . For any given initial state, how
can we identify the invariant classes from the set of invariant
clusters? Suppose we want to find the invariant classes that can
overapproximate the trajectory from the state ~x0 = (1, 2, 3).
According to Theorem 2, we have Algorithm 2 for the purpose.
Remark 4: In Algorithm 2, we enumerate the invariant
clusters to find out which one can provide a non-empty
invariant class Class(C, ~x0 ) for ~x0 . For a Class(C, ~x0 ) to
be nonempty, the corresponding hyperplane g(~u, ~x0 ) = 0 must
have at least one solution to ~u ∈ RK \{~0}, which is equivalent
to that its dimension must be at least 1. For a hyperplane in
RK , its dimension is equal to K − 1. Therefore, Class(C, ~x0 )
must be nonempty if K > 1 and the basis of the hyperplane
can be obtained through a basic linear algebraic computation

5

6
7

8

9
10
11
12

ICls ← ∅;
foreach C ∈ CFamily do
D ← Class(C, ~x0 );
if D 6= ∅ then
m ← the dimension of the hyperplane
g(~u, ~x0 ) = 0 defining D;
if m ≥ 1 then
Basis ← basis {u1 , . . . , um } of the
hyperplane g(~u, ~x0 ) = 0;
D ← the polynomials
{g(~u1 , ~x), . . . , g(~um , ~x)}, ui ∈ Basis;
end
ICls ← ICls ∪ D;
end
end

(which will be illustrated in what follows). However, in the
case of g(~u, ~x0 ) being identical to 0, the hyperplane degenerates to the space RK \ {~0} and the dimension will be K.
Therefore, an invariant class with K = 1 is nonempty iff
g(~u, ~x0 ) is identical to 0. For example, given an invariant
cluster C 0 = {u1 (x − y) = 0 | u1 ∈ R \ {0}} and a point
~x0 = (x0 , y0 ), Class(C 0 , ~x0 ) is empty if x0 6= y0 , however,
Class(C 0 , ~x0 ) is equal to C 0 if x0 = y0 .
Example 5: We continue from Example 4. For the given
point ~x0 = (1, 2, 3), according to Algorithm 2, we find that
only Class(C7 , ~x0 ) = {g7 (~u, x, y, z) = 0 | 3u5 + 8u6 + u0 =
0, ~u ∈ R3 \{~0}} is nonempty. The dimension of the hyperplane
H : 3u5 + 8u6 + u0 = 0 is 2. Since u0 = −3u5 − 8u6 , to
get the basis of H, we can write (u0 , u5 , u6 ) = (−3u5 −
8u6 , u5 , u6 ) = u5 (−3, 1, 0) + u6 (−8, 0, 1). Hence, we have
the following basis for H: {(−3, 1, 0), (−8, 0, 1)}. As a result,
we get the finite representation B = {y 2 − x2 − 3 = 0,
z 2 − x2 − 8 = 0} for Class(C7 , ~x0 ). It is easy to check by
the Maple function HilbertDimension that dim(B) = 1.
Therefore, we finally get an algebraic variety V(B) which
provides in part a curve-to-curve match to the trajectory π~x0 .
The 3-D vector field and the algebraic curve V(B) is shown
in Figure 3a and Figure 3b, respectively.
IV. S AFETY V ERIFICATION BASED
C LUSTERS

ON I NVARIANT

A. Safety Verification of Continuous Systems
In this section, we show how to verify a safety property for
a nonlinear system based on invariant clusters.
In Section III, we have demonstrated that a trajectory can
be overapproximated by an invariant class. Since an invariant
class is determined uniquely by a single hyperplane g(~u, ~x0 ) =
0 in RK for an initial point ~x0 , and a hyperplane without

j2 = m1 + 1, . . . , m1 + m2 , k2 = n1 + 1, . . . , n1 + n2 }.
Then we have the following theorem for deciding the safety
of a semialgebraic system.
Theorem 4: Given a semialgebraic system M = hX, f~, X0 i
and invariant cluster C = {g(~u, ~x) = 0 | ~u ∈ RK \ {~0}} of
M with K ≥ 2. Suppose the normal vector of the hyperplane
g(~u, ~x) = 0 over ~u is (1, ψ1 (~x), . . . , ψK (~x)). Then the system
M is safe if there exists the following polynomial identity
(a)

(b)

Fig. 3: (a) 3D vector field of Example 4. (b) The intersection
of the invariants y 2 − x2 − 3 = 0 (blue) and z 2 − x2 − 8 =
0 (orange) overapproximates the trajectory originating from
~x0 = (1, 2, 3) (green ball).

K
X

k=1

γk · (ψk (~x1 ) − ψk (~x2 )) +

+

v∈{0,1}m1 +m2

+
constant term (this is the case for g(~u, ~x0 ) = 0) is uniquely
determined by its normal vector, we can verify if two states
lie on the same trajectory based on the following theorem.
Theorem 3: Given a semialgebraic system M = hX, f~, X0 i
and an invariant cluster C = {g(~u, ~x) = 0 | ~u =
(u1 , · · · , uKP
) ∈ RK \ {~0}} of M with K > 1, where
K
x)ui and ψi (~x) ∈ R[~x], let X0 be the
g(~u, ~x) =
i=1 ψi (~
initial set and XU be the unsafe set. Then, if there exists a
pair of states (~x1 , ~x2 ) ∈ X0 × XU such that ~x1 and ~x2 lie on
the same trajectory, one of the following two formulae must
hold:
(i)
(ii)

∃k ∈ R\{0} : kψi (~x1 ) = ψi (~x2 ), i = 1, . . . , K
ψi (~x1 ) = ψi (~x2 ) = 0, i = 1, . . . , K

(4)
(5)

Moreover, if some ψi (~x) ≡ 1, i.e. g(~u, ~x) contains a
constant term ui , then formula (4) simplifies to
ψi (~x1 )) = ψi (~x2 ), i = 1, . . . , K

(6)

Proof: See Appendix D.
Remark 5: Instead of computing the invariants explicitly,
Theorem 3 provides an alternative way to verify whether two
states ~x1 , ~x2 lie on the same trajectory by checking that the
difference between the normal vectors of g(~u, ~x1 ) = 0 and
g(~u, ~x2 ) = 0. Let us take Example 4 for illustration. We think
of C7 as a hyperplane over ~u ∈ R3 : u0 + (y 2 − x2 )u5 +
(z 2 − x2 )u6 = 0, hence the corresponding normal vector is
~ (~x) = (1, y 2 − x2 , z 2 − x2 ). Given two random points ~x1 =
N
√ √
(1, 2, 3) and ~x2 = (5, 27, 34), it is easy to verify that
~ (~x1 ) 6= N
~ (~x2 ), which means that ~x1 and ~x2 are not on
N
the same trajectory. This can be verified in another way, as
we know that the invariant class of ~x1 is {y 2 − x2 − 3 =
0, z 2 − x2 − 8 = 0} and ~x2 does not belong to its solution set.
Now we demonstrate how to verify a safety property of
semialgebraic systems. Assume X0 and XU can be written
as semialgebraic sets, i.e., X0 = {~x1 ∈ Rn | pi1 (~x1 ) =
0, qj1 (~x1 ) ≥ 0, rk1 (~x1 ) > 0, i1 = 1, . . . , l1 , j1 = 1 . . . m1 ,
k1 = 1, . . . , n1 } and XU = {~x2 ∈ Rn | pi2 (~x2 ) =
0, qj2 (~x2 ) ≥ 0, rk2 (~x2 ) > 0, i2 = l1 + 1, . . . , l1 + l2 ,

X

X

v∈{0,1}n1 +n2

i=1
m1Y
+m2

X
(
btv )2 ·
t

lX
1 +l2

βi p i

v

qj j

(7)

j=1

n1Y
+n2
n1Y
+n2
X
X
(
ctv )2 ·
rkvk +
s2w +
rkdk = 0
t

k=1

w

k=1

where dk ∈ N and βj , γk , btv , ctv , sw are polynomials in
R[~x1 , ~x2 ].
Proof: See Appendix E.
Remark 6: Theorem 4 transforms the safety verification
problem into a decision problem about the existence of a real
solution of a system of polynomial equations and inequalities.
As noted in Remark 1, this decision problem can be solved
by SOS programming. Our implementation uses the efficient
tool SOSTOOLS [12]. Appendix C contains more information
on SOS programming.
In Theorem 4, we deal with a general semialgebraic system
where the initial set and the unsafe set are represented by a
set of polynomial equations and inequalities. However, if the
system is described by much simpler set representations such
as a single polynomial equation or inequality, the programming
problem can be simplified correspondingly. For example, if
both sets can be represented or overapproximated by a single
polynomial equation I(~x1 ) = 0 and U (~x2 ) = 0, respectively,
then the programming problem is simplified to (see [11])
K
X
αj (ψj (~x1 ) − ψj (~x2 )) + βI + θU − 1 is an SOS (8)
j=1

where (ψ1 (~x), . . . , ψK (~x)) is the same as in Theorem 4 and
αj , β, θ ∈ R[~x1 , ~x2 ]. For how to overapproximate a compact
set we refer the reader to Appendix J. The algorithm for
safety verification based on the condition (8) is shown in
Algorithm 3.
The next example demonstrates the application of the verification method.
Example 6 (running example
 2): Given the semialgebraic
system M3 by [ẋ, ẏ] = y 2 , xy and the initial set be
X0 = {(x, y) ∈ R2 | I(x, y) = (x+15)2 +(y −17)2 −1 ≤ 0},
verify that if the unsafe set XU = {(x, y) ∈ R2 | U (x, y) =
(x−11)2 +(y−16.5)2 −1 ≤ 0} can be reached. The parameter
space of C ∗ = {g(~u, ~x) = u1 − u3 (x2 − y 2 ) = 0 | (u1 , u3 ) ∈
R2 \ {~0}} has dimension 1 and hence can provide an invariant
class for every state in X0 and XU . The normal vector of the
hyperplane g(~u, ~x) = 0 is (1, ψ1 (x, y)) = (1, y 2 − x2 ). Let

Algorithm 3: Safety verification
~ the K-dimensional normal vector of an
input : ψ:
invariant cluster;
I(~x1 ): the initial set;
U(~x2 ): the unsafe set;
N : the maximum degree of programming
polynomials α
~ , β, θ
output: IsSafe: whether the system is safe
1
2
3

4
5
6
7
8
9
10
11
12

IsSafe ← False;
for i ← 1 to N do
α
~ ← generate a vector of polynomials of degree i for
ψ;
β ← generate a polynomial of degree i for I(~x1 );
θ ← generate a polynomial of degree i for U(~x2 );
PK
~j (~x1 ) − ψ
~j (~x2 )) + βI+θU −1;
P ← j=1 αj (ψ
Solution ← perform SOS programming on P ;
if Solution is found then
IsSafe ← True;
break;
end
end

ϕ(x1 , y1 , x2 , y2 ) = ψ1 (x1 , y1 ) − ψ1 (x2 , y2 ). By Theorem 4,
to verify the safety property, we only need to verify that the
following system of equations has no real solution.
I(x1 , y1 ) = (x1 + 15)2 + (y1 − 17)2 − 1 = 0

U (x2 , y2 ) = (x2 − 11)2 + (y2 − 16.5)2 − 1 = 0
ϕ(x1 , y1 , x2 , y2 ) = y12 − x21 − (y22 − x22 ) = 0

Note that we substitute (x1 , y1 ), (x2 , y2 ) for (x, y) in I(x, y)
and U (x, y), respectively, to denote the different points in X0
and XU . To prove that the system is safe, we need to find
αi ∈ R[x1 , y1 , x2 , y2 ], i = 1, 2, 3 such that Prog = α1 I +
α2 U + α3 ϕ − 1 is a sum-of-squares. Finally, we found three
polynomials of degree 2 for αi , respectively (see Appendix K
for the expressions of αi and Prog), hence the system is safe.
As shown in Figure 2, although the relative position of XU
to the reachable set of X0 is very close, we can still verify
the safety property using an invariant cluster. However, we
failed to find a barrier certificate for this system by using the
methods in [14], [15].
In Theorem 4, we present a sufficient condition for deciding
if a semialgebraic system is safe. The theory originates from
the fact that the system is safe if there is no invariant class
intersecting both the initial and the unsafe set, which is equivalent to that the formula (7) holds. To verify the latter, we need
to find a set of witness polynomials by SOS programming.
However, as the dimension of the system increases, the number
of parametric polynomials involved increases correspondingly,
which also leads to an increase in computational complexity. In
what follows, we present a new method for safety verification,
which avoids the aforementioned problem. The new method
is based on Proposition 1, that is, for any polynomial g(~x)

satisfying Lf~g ∈ hgi, g(x) ∼ 0 is an invariant for any
∼ ∈ {<, ≤, =, ≥, >}.
Proposition 3: Given a semialgebraic system M =
hX, f~, X0 i and an invariant cluster C = {g(~u, ~x) = 0 | ~u ∈
RK \ {~0}} of M , let X0 and XU be the initial set and the
unsafe set, respectively. Then, the system is safe if there exists
a ~u∗ ∈ RK \ {~0} such that
∀~x ∈ X0 : g(~u∗ , ~x) ≥ 0
∗

∀~x ∈ XU : g(~u , ~x) < 0

(9)
(10)

Proof: See Appendix F.
According to Proposition 3, to verify the safety property,
it suffices to find a ~u∗ ∈ RK \ {~0} which satisfies the
constraints (9) and (10). There are some constraint solving
methods available, e.g, SMT solvers. However, the high complexity of SMT theory limits the applicability of the method.
In the following, we transform the above constraint-solving
problem into an SOS programming problem, which can be
solved efficiently. We write P~ (~x)  ~0 to denote pi (~x) ≥ 0, i =
1, . . . , m for a polynomial vector P~ (x) = (p1 (~x), . . . , pm (~x)).
Proposition 4: Given a semialgebraic system M =
hX, f~, X0 i and an invariant cluster C = {g(~u, ~x) = 0 |,
~u ∈ RK \ {~0}} of M and a constant ǫ ∈ R>0 , let
X0 = {~x ∈ Rn | I~  ~0, I~ ∈ R[~x]m1 } and XU = {~x ∈ Rn |
~  ~0, U
~ ∈ R[~x]m2 }. Then, the system is safe if there exist a
U
∗
K
~u ∈ R \ {~0} and two SOS polynomial vectors ~µ1 ∈ R[~x]m1
and ~µ2 ∈ R[~x]m2 such that the following polynomials are SOS
polynomials.
g(~u∗ , ~x) − µ
~ 1 · I~
~ −ǫ
− g(~u∗ , ~x) − ~µ2 · U

(11)
(12)

Proof: See Appendix G.
Similar to Theorem 4, Proposition 4 also reduces to an
SOS programming problem. However, the ideas behind these
two theories are different in that by Theorem 4 we attempt
to prove no invariant class which overapproximates trajectory
can intersect both X0 and XU , while by Proposition 4 we
mean to find a hypersurface which is able to separate the
reachable set of X0 from the unsafe set XU . Apparently,
there must exist no invariant class intersecting both X0 and
XU if there exists such a hypersurface, but not vice versa.
Hence the latter is more conservative than the former, but it
is also more efficient in theory because it usually involves
less unknown polynomials. For example, for an n-dimensional
system with X0 and XU defined by a single polynomial
inequality, respectively, we usually need n + 1 unknown
polynomials for the former method, however, we need only
2 for the latter. See Algorithm 4 for the pseudocode of the
method based on Proposition 4.
Let us use the following example to demonstrate the application of Algorithm 4.
Example 7: Given the semialgebraic system S4 by
   2

ẋ
y − 2y
= 2
,
ẏ
x + 2x

Algorithm 4: Algorithm for safety verification based on
SOS programming
input : g(~u, ~x): the invariant polynomial defining an
invariant cluster;
~ the polynomial vector describing initial set;
I:
~ : the polynomial vector describing unsafe set;
U
N : the maximum degree of the polynomial
vectors ~µ1 , ~µ2 ;
ǫ: a positive real number
output: IsSafe: whether the system is safe;
g(~u∗ , ~x), µ
~ 1, µ
~ 2 : the feasible solution if IsSafe;
1
2
3

4

5
6
7

8
9
10

11
12
13

IsSafe ← False;
for i ← 0 to N do
~µ1 ← generate a parametric polynomial vector of
~
degree 2i for I;
~µ2 ← generate a parametric polynomial vector of
~;
degree 2i for U
~
p1 ← g(~u, ~x) − ~
µ1 · I;
~ − ǫ;
p2 ← −g(~u, ~x) − ~
µ2 · U
Solution ← perform SOS programming on
{~µ1 , ~µ2 , p1 , p2 };
if Solution is found then
IsSafe ← True;
[g(~u∗ , ~x), ~
µ1 , ~
µ2 ] ← get the feasible solution from
Solution;
break;
end
end

Fig. 4: Example 7. Invariant g(~u∗ , ~x) = 0 (blue curve)
separating reachable set of X0 (green patch) from XU (red
patch).

all locations. The idea is to pick out a polynomial gl (~u∗l , ~x)
from the respective invariant cluster Cl for each location l such
that gl (~u∗l , ~x) ≥ 0 is an invariant for the location l and all the
invariants coupled together by the constraints at the discrete
transitions form a hybrid invariant for the hybrid system. The
aforementioned idea is formalized in Proposition 5.
Proposition 5: Given an n-dimensional hybrid system H =
hL, X, E, R, G, I, F i and a set of invariant clusters {Cl , l =
1, . . . , n}, where Cl = {gl (~ul , ~x) = 0 | ~ul ∈ RKl \ {~0}} with
Kl > 1 is an invariant cluster for location l, the system H is
safe if there exists a set Su~ = {~u∗l ∈ RKl \ {~0}, l = 1, . . . , n}
such that, for all l ∈ L and (l, l′ ) ∈ E, the following formulae
hold:
∀~x ∈ Init(l) : gl (~u∗l , x) ≥ 0
∀~x ∈ G(l, l′ ), ∀~x′ ∈ R((l, l′ ), ~x) : gl (~u∗l , ~x) ≥ 0

2

let the initial set be X0 = {(x, y) ∈ R | I(x, y) = 1 − (x +
6.0)2 − (y + 6.0)2 ≥ 0}, decide whether the unsafe set XU =
{(x, y) ∈ R2 | U (x, y) = 1 − (x − 8.2)2 − (y − 4.0)2 ≥ 0}
can be reached.
By Algorithm 1, we first get an invariant cluster C4 =
{g(~u, ~x) = 13 u3 x3 − 31 u3 y 3 + u3 x2 + u3 y 2 + u1 = 0 | ~u ∈ R2 }
for the system S4 , then by Algorithm 4, we find an invariant
g(~u∗ , ~x) = 0 with ~u∗ = (−3081.9, 7.1798) from the invariant
cluster C4 : g(~u∗ , ~x) = 7.1798x3 − 7.1798y 3 + 21.539x2 +
21.539y 2 − 3081.9 and two polynomials µ1 (~x), µ2 (~x) of
degree 2. The stream plot of S4 and the plot of g(~u∗ , ~x) = 0
are shown in Figure 4. Note that we failed to find a barrier
certificate by using the method in [15] and [14] for this system.
B. Safety Verification of Hybrid Systems
In this section, we extend the safety verification method for
continuous systems based on invariant clusters to semialgebraic hybrid systems.
A hybrid system consists of a set of locations and a set
of discrete transitions between locations. In general, different
locations have different continuous dynamics and hence correspond to different invariant clusters. An invariant for the hybrid
system can be synthesized from the set of invariant clusters of

=⇒ gl′ (~u∗l′ , ~x′ ) ≥ 0
∀~x ∈ I(l) ∩ Uns(l) : gl (~u∗l , ~x) < 0

(13)
(14)
(15)

where Init(l) and Uns(l) denote respectively the initial set
and the unsafe set at location l.
Proof: See Appendix H.
Similar to Proposition 3, we further transform the problem
into an SOS programming problem. Consider a semialgebraic
hybrid system H = hL, X, E, R, G, I, F i, where the mappings
R, G, and I are defined in terms of polynomial inequalities as
follows:
′
~ ll′  0, G
~ ll′ ∈ R[~x]mll′ }
• G : (l, l ) 7→ {~
x ∈ Rn | G
′
~ ll′ ~x  0, R
~ ll′ ~x ∈ R[~x]nll′ }
• R : (l, l , ~
x) 7→ {~x ∈ Rn | R
• I : l 7→ {~
x ∈ Rn | I~l  0, I~l ∈ R[~x]pl }
and the mappings of the initial and the unsafe set are defined
as follows:
~ l  0, Init
~ l ∈ R[~x]rl }
• Init : l 7→ {~
x ∈ Rn | Init
n
~
~ l ∈ R[~x]sl }
• Uns : l 7→ {~
x ∈ R | Unsl  0, Uns
where mll′ , nll′ , rl , pl and sl are the dimensions of the polynomial vector spaces. Then we have the following proposition
for safety verification of the semialgebraic hybrid system H.
Proposition 6: Let the hybrid system H, the initial set
mapping Init, and the unsafe set mapping Uns be defined

as above. Given a set of invariant clusters {Cl , l = 1, . . . , n}
of H, where Cl = {gl (~ul , ~x) = 0 | ~ul ∈ RKl \ {~0}} with
Kl > 1 is an invariant cluster for location l, a set Sγ =
{γll′ ∈ R≥0 , (l, l′ ) ∈ E} of constants, and a constant vector
~ǫ ∈ Rn>0 , the system is safe if there exists a set Su = {~u∗l ∈
RKl \ {~0}, l = 1, . . . , n} and five sets of SOS polynomial
vectors {~θl ∈ R[~x]sl , l ∈ L}, {~κll′ ∈ R[~x]pll′ , (l, l′ ) ∈ E},
ηl ∈ R[~x]tl , l ∈ L}, and
{~σll′ ∈ R[~x]qll′ , (l, l′ ) ∈ E}, {~
wl
{~νl ∈ R[~x] , l ∈ L} such that the following polynomials
are SOS for all l ∈ L and (l, l′ ) ∈ E:
~ l
gl (~u∗l , ~x) − θ~l · Init

gl′ (~u∗l′ , ~x′ )

(16)

γll′ gl (~u∗l , ~x)

− ~κll′ · Gll′ − ~σll′ · Rll′ ~x
~
− ~νl · Il − ~ηl · Unsl − gl (~u∗l , ~x) − ǫl
−

(17)
(18)

Proof: Similar to the proof of Proposition 4, we can
easily derive the formulae (13)–(15) from the SOS’s (16)–(18),
respectively. Then, by Proposition 5, the system is safe.
The algorithm for computing invariants for semialgebraic
hybrid systems based on Proposition 6 is very similar to
Algorithm 4 for semialgebraic continuous systems except that
it involves more SOS constraints on continuous and discrete
transitions.
V. I MPLEMENTATION AND E XPERIMENTS
Based on the approach presented in this paper, we implemented a prototype tool in Maple and Matlab, respectively.
In Maple, we implemented the tool for computing invariant
clusters and identifying invariant classes based on the remainder computation of the Lie derivative of a polynomial w.r.t. its
Gröbner basis and solving the system of polynomial equations
obtained from the coefficients of the remainder. In Matlab,
we implemented the tool for safety verification based on the
SOS programming tool package SOSTOOLS. Currently, we
manually transfer the invariant clusters computed in Maple to
Matlab for safety verification. In the future, we will integrate
the two packages into a single tool.
Now, we present the experimental results on nonlinear
benchmark systems, run on a laptop with an 3.1GHz Intel
Core i7 CPU and 8GB memory.
A. Longitudinal Motion of an Airplane
In this experiment, we study the 6th order longitudinal equations of motion that captures the vertical motion (climbing,
descending) of an airplane ([16], Chapter 5). Let g denote the
gravity acceleration, m the total mass of an airplane, M the
aerodynamic and thrust moment w.r.t. the y axis, (X, Z) the
aerodynamics and thrust forces w.r.t. axis x and z, and Iyy the
second diagonal element of its inertia matrix. Then the motion
of the airplane is described as follows.
X
− g sin(θ) − qw,
m
ẋ = w sin(θ) + v cos(θ),
v̇ =

θ̇ = q,

Z
+ g cos(θ) + qv,
m
ż = −v sin(θ) + w cos(θ),
M
,
q̇ =
Iyy

ẇ =

where the meanings of the variables are as follows: v: axial
velocity, w: vertical velocity, x: range, z: altitude, q: pitch rate,
θ: pitch angle.
To transform the above system into a semialgebraic system,
we first introduce two additional variables d1 , d2 such that
d1 = sin(θ), d2 = cos(θ) and then substitute d1 and d2
respectively for sin(θ) and cos(θ) in the model. In addition, we
get two more constraints d˙1 = qd2 and d˙2 = −qd1 . As a result,
the dimension of the system rises to 8. For this system, using
the method in [17], Ghorbal et al. spent 1 hour finding three
invariant polynomials of degree 3 on a laptop with a 1.7GHz
Intel Core i5 CPU and 4GB memory. Using our method,
we spent only 0.406 seconds obtaining an invariant cluster
g9 (~u, ~x) = 0 of degree 3, which is presented in Appendix L.
By the constraint d21 + d22 = 1, we reduce the normal vector
of the hyperplane g9 (~u, ~x) = 0 in ~u to (1, ψ1 , ψ2 , ψ3 ), where
ψ1 , ψ2 , ψ3 are defined as follows.

M mz gmθ  mqv
ψ1 =
+
+
+ 1 sin (θ)
Iyy Z
Z
Z


X
mqw
+
−
cos (θ)
Z
Z


Xqv qw
gIyy Xθ
Xz
sin(θ)
+x−
− Iyy
+
ψ2 = −
Z
ZM
ZM
M


qv
X2 + Z2
Xqw
cos(θ)
−
−
+ Iyy
ZM
M
ZM m
Mθ
ψ3 = q 2 − 2
Iyy
Given a symbolic initial point ~x0 = (v0 , w0 , x0 , z0 , θ0 , q0 ,
d01 , d02 ), we have verified that our invariant cluster defines the
same algebraic variety as defined by the invariants in [17]
by comparing their Gröbner bases. However, our method is
much more efficient. Moreover, we also obtained the invariant
clusters of higher degrees (4 − 6) quickly. The experimental
result is shown in Table I. The first column is the degree of the
invariants, the second column is the variables to be decided,
the third column is the computing time in seconds, and the
last column is the number of invariant clusters generated. As
can be seen, in the most complicated case, where the number
of the indeterminates reaches up to 3003, we spent only 200.9
seconds to discover an invariant cluster of degree 6. However,
we found that these higher order invariant clusters have the
same expressive power as the invariant cluster of degree 3 in
terms of algebraic variety.
B. Looping particle
Consider a heavy particle on a circular path of radius r. The
motion of the particle is described by the following differential
equation.
 
  
 

˙
rcos(θ)
−r sin(θ)θ̇
ẋ
−yω
˙ 
 ẏ  = 
 rsin(θ)
 =  r cos(θ)θ̇  =  xω 
g cos(θ)
ω̇
− gx
− g(r cos(θ))
−
2
r2
r

r

Note that the above is a parameterized system with gravity
acceleration g and radius r as parameters. Our tool finds the

TABLE I: Benchmark results for the Longitudinal Motion of an Airplane (B1), the Looping Particle system (B2), the Coupled
Spring-Mass system (B3), and the 3D-Lotka-Volterra System (B4).
Degree of
invariants
1
2
3
4
5
6

B1
9
45
165
495
1287
3003

No. of variables
B2
B3
4
6
10
21
20
56
35
126
56
252
84
462

B4
4
10
20
35
56
84

B1
0.016
0.031
0.484
3.844
25.172
200.903

Running time (sec)
B2
B3
0.015
0.047
0.047
0.078
0.049
0.250
0.156
1.109
0.703
6.641
3.000
32.109

following invariant cluster consisting of a parametric polynomial of degree 2: {g(~u, ~x) = 0 | g(~u, ~x) = u5 x2 + u5 y 2 +
u2 ω 2 + 2ur22 g y + u0 , ~u ∈ R3 \ {~0}}. Given an arbitrary point
(x0 , y0 , ω0 ) = (2, 0, ω0 ), we get the invariant class {g(~u, ~x) =
0 | (x20 + y02 )u5 + (ω02 + 2g
u ∈ R3 \ {~0}}.
r 2 y0 )u2 + u0 = 0, ~
According to Algorithm 2, the algebraic variety representing
the trajectory originating from (x0 , y0 , ω0 ) is {(x, y, ω) ∈ R3 |
2g
2
x2 + y 2 − x20 − y02 = 0, ω 2 + 2g
r 2 y − ω0 − r 2 y0 = 0}. The results
in [19] and [18] are special cases of our result when setting
(r, g) = (2, 10) and (r, g, x0 , y0 ) = (2, 10, 2, 0), respectively.
Therefore, our method is more powerful in finding parameterized invariants for parameterized systems. See Table I for
detailed experimental result.
C. Coupled spring-mass system
Consider the system

  
v1
x˙1
k
k
 v˙1  − 1 x1 − 2 (x1 − x2 )
m1

  =  m1

x˙2  
v2
k
2
v˙2
− m2 (x2 − x1 )

The model consists of two springs and two weights w1 , w2 .
One spring, having spring constant k1 , is attached to the ceiling
and the weight w1 of mass m1 is attached to the lower end
of this spring. To the weight w1 , a second spring is attached
having spring constant k2 . To the bottom of this second spring,
a weight w2 of mass m2 is attached. x1 and x2 denote the
displacements of the center of masses of the weights w1 and
w2 from equilibrium, respectively.
In this benchmark experiment, we first tried an instantiated
version of the system by using the same parameters as in [20]:
k2
k1
m1 = m2 = k and m1 = 5m2 . The experimental result
is presented in Table I. We found that the expressive power
of the invariant clusters does not increases any more as the
degree is greater than 3 and it took only 0.250 seconds to
compute the invariant cluster of degree 3. Finally, we perform
the computation directly on the fully parameterized system
and we get the following parameterized invariant cluster.
k2 x1 x2 (m1 u8 − 2m2 u10 )
+ u10 v12 + u1
g(~u, ~x) = u8 v1 v2 +
m1 m2
1 v22 (k1 m2 u8 − k2 m1 u8 + k2 m2 u8 + 2k2 m2 u10 )
+
2
k2 m1
1 (2k1 m2 u10 − k2 m1 u8 + 2k2 m2 u10 )x21
+
2
m1 m2
1 (k1 m2 u8 − k2 m1 u8 + 2k2 m2 u10 )x22
+
2
m1 m2

B4
<0.001
0.031
0.109
0.312
0.750
1.641

B1
0
1
1
1
1
1

No. of invariant clusters
B2
B3
0
0
1
0
0
1
1
1
0
1
1
1

B4
2
3
7
6
6
16

This invariant cluster enables us to analyze the system
properties under different parameter settings.
D. 3D-Lotka-Volterra system
Consider the system [ẋ, ẏ, ż] = [xy − xz, yz − yx, xz − yz].
The experimental result is presented in Table I. Here we
present only two invariant clusters to show their expressive
power: one of degree 1 and one of degree 3. C1 = {g1 (~u, ~x) =
u3 x + u3 y + u3 z + u4 = 0 | ~u ∈ R2 \ {~0}}, C3 =
{g3 (~u, ~x) = 0 | ~u ∈ R4 \ {~0}}, where g3 (~u, ~x) = u5 xyz +
u10 x3 + 3u10 x2 y + 3u10 x2 z + 3xy 2 u10 + 3xz 2 u10 + u10 y 3 +
3u10 y 2 z + 3u10 yz 2 + u10 z 3 + u16 x2 + 2u16 xy + 2u16 xz +
u16 y 2 + 2u16 yz + u16 z 2 + u19 x + u19 y + u19z + u20 . Invariant
cluster C1 shows that every trajectory of the system must
lie in some plane. C3 can overapproximate the trajectories
with a much higher precision because it can provide algebraic
varieties of dimension 1. For example, in [20], the initial
set is given as X0 = {(x, y, z) ∈ R3 | x2 − 1 = 0,
y 2 − 1 = 0, z 2 − 1 = 0} and they discovered four invariants.
In our experiment, we spend 0.109 seconds discovering an
invariant cluster. This invariant cluster can overapproximate
all the trajectories precisely, say, for x0 = (1, −1, 1) ∈ X0 ,
we get the invariant class Class(C3 , x0 ) = {g3 (~u, ~x) = 0 |
−u5 +7u10 +u16 +u19 +u20 = 0, ~u ∈ R5 \{~0}}. By taking the
intersection of its basis xyz + 1 = 0 and x + y + z − 1 = 0, we
immediately obtain a one-dimensional algebraic variety which
overapproximates the trajectory precisely.
E. Hybrid controller
Consider a hybrid controller consisting of two control
modes. The discrete transition diagram of the system is shown
in Fig. 5a and the vector fields describing the continuous
behaviors are as follows:


y 2 + 10y + 25
f1 (~x) =
,
2xy + 10x − 40y − 200


−y 2 − 10y − 25
f2 (~x) =
8xy + 40x − 160y − 800
The system starts from some point in X0 = {(x, y) ∈ R2 |
(x−9)2 +(y−20)2 ≤ 4} and then evolves following the vector
field f1 (~x) at location l1 (Switch-On). The value of x keeps
increasing until it reaches 35, then the system switches immediately to location l2 (Switch-Off) without performing any
reset operation. At location l2 , the system operates following
the vector field f2 (~x) and the value of x keeps decreasing. As
the value of x drops to 5, the system switches immediately

x


 d x d 

x f x

 d x d 

x f x

l

l
x



(a)

(b)

Fig. 5: (a) Hybrid automaton from Subsection V-E. x = 5
and x = 35 are guards for discrete transitions and no reset
operation is performed. (b) Hybrid invariant for the system in
Subsection V-E. Solid patch in green: initial set, curve in blue:
invariant for l1 , curve in purple: invariant for l2 , red shadow
region on the top: unsafe region.

back to location l1 again. Our objective is to verify that the
value of y will never exceed 48 in both locations.
For the convenience of SOS programming, we define the
unsafe set as Uns(l1 ) = Uns(l2 ) = {(x, y) ∈ R2 | 48 < y <
60}, which is sufficient to prove y ≤ 48 in locations l1 and
l2 . According to the theory proposed in Section IV-B, we first
find an invariant cluster for each location, which is composed
of a parameterized polynomial, respectively: g1 (~u1 , ~x) =
1
u12 y 2 + 8u12 x + u12 y + u11 and g2 (~u2 , ~x) =
− 51 u12 x2 + 10
1
4
2
2
5 u22 x + 10 u22 y −32u22 x+u22 y +u21 . In the second phase,
we make use of the constraint condition in Proposition 6 to
compute a pair of vectors ~u∗1 and ~u∗2 . By setting γ12 = γ21 = 1,
our tool found a pair of ~u∗1 = (u11 , u12 ) = (2.9747, 382.14)
and ~u∗2 = (u21 , u22 ) = (2.9747, 138.44). As shown in Fig. 5b,
the curves of g1 (~u∗1 , ~x) = 0 and g2 (~u∗2 , ~x) = 0 form an upper
bound for the reachable set in location l1 and l2 , respectively,
which lie below the unsafe region y ≥ 48. Therefore, the
system is safe.
VI. R ELATED W ORK
In recent years, many efforts have been made toward generating invariants for hybrid systems. Matringe et al. reduce
the invariant generation problem to the computation of the
associated eigenspaces by encoding the invariant constraints
as symbolic matrices [21]. Ghorbal et al. use the invariant
algebraic set formed by a polynomial and a finite set of its successive Lie derivatives to overapproximate vector flows [17].
Both of the aforementioned methods involve minimizing the
rank of a symbolic matrix, which is inefficient in dealing
with parametric systems. In addition, none of these methods
involve how to verify safety properties based on the invariants.
Sankaranarayanan discovers invariants based on invariant ideal
and pseudo ideal iteration [20], but this method is limited to
algebraic systems. Tiwari et al. compute invariants for special
types of linear and nonlinear systems based on Syzygy computation and Gröbner basis theory as well as linear constraint
solving [22]. Platzer et al. use quantifier elimination to find
differential invariants [23]. The methods based on Gröbner
basis and first-order quantifier elimination suffer from the high

complexity significantly. Another approach considers barrier
certificates based on different inductive conditions [15], [14]
which can be solved by SOS programming efficiently but is
limited by the conservative inductive condition. Carbonell et
al. generate invariants for linear systems [24]. Some other
approaches focusing on different features of systems have also
been proposed for constructing inductive invariants [25], [26],
[18], [27], [28].
VII. C ONCLUSION
In this paper, we proposed an approach to automatically
generate invariant clusters for semialgebraic hybrid systems.
The benefit of invariant clusters is that they can overapproximate trajectories of the system precisely. The invariant clusters
can be obtained efficiently by computing the remainder of
the Lie derivative of a template polynomial w.r.t. its Gröbner
basis and then solving a system of homogeneous polynomial
equations obtained from the remainder. Moreover, based on
invariant clusters and SOS programming, we propose a new
method for safety verification of hybrid systems. Experiments
show that our approach is effective and efficient for a large
class of biological and control systems.
ACKNOWLEDGEMENT
This research was supported in part by the European Research Council (ERC) under grant 267989 (QUAREM) and by
the Austrian Science Fund (FWF) under grants S11402-N23
(RiSE) and Z211-N23 (Wittgenstein Award).
R EFERENCES
[1] T. Henzinger, “The theory of hybrid automata,” in Proc. IEEE Symp.
Logic in Computer Science, 1996, pp. 278–292.
[2] E. Asarin, T. Dang, and A. Girard, “Reachability analysis of nonlinear
systems using conservative approximation,” in Hybrid Systems: Computation and Control. Springer, 2003, pp. 20–35.
[3] S. Bogomolov, A. Donzé, G. Frehse, R. Grosu, T. T. Johnson, H. Ladan,
A. Podelski, and M. Wehrle, “Guided search for hybrid systems based
on coarse-grained space abstractions,” International Journal on Software
Tools for Technology Transfer, pp. 1–19, 2015.
[4] S. Bogomolov, G. Frehse, R. Grosu, H. Ladan, A. Podelski, and
M. Wehrle, “A box-based distance between regions for guiding the
reachability analysis of SpaceEx,” in Computer Aided Verification (CAV
2012), ser. LNCS, vol. 7358. Springer, 2012, pp. 479–494.
[5] A. Tiwari, “Abstractions for hybrid systems,” Formal Methods in System
Design, vol. 32, no. 1, pp. 57–83, 2008.
[6] R. Alur, T. Dang, and F. Ivančić, “Progress on reachability analysis
of hybrid systems using predicate abstraction,” in Hybrid Systems:
Computation and Control. Springer, 2003, pp. 4–19.
[7] S. Bogomolov, G. Frehse, M. Greitschus, R. Grosu, C. S. Pasareanu,
A. Podelski, and T. Strump, “Assume-guarantee abstraction refinement
meets hybrid systems,” in HVC 2014. Springer, pp. 116–131.
[8] S. Bogomolov, C. Schilling, E. Bartocci, G. Batt, H. Kong, and R. Grosu,
“Abstraction-based parameter synthesis for multiaffine systems,” in 11th
International Haifa Verification Conference (HVC 2015), ser. LNCS, vol.
9434. Springer, 2015, pp. 19–35.
[9] S. Bogomolov, C. Herrera, M. Muñiz, B. Westphal, and A. Podelski,
“Quasi-dependent variables in hybrid automata,” in Hybrid Systems:
Computation and Control. ACM, 2014, pp. 93–102.
[10] D. A. Cox, J. Little, and D. O’Shea, Ideals, Varieties, and Algorithms:
An Introduction to Computational Algebraic Geometry and Commutative
Algebra. Springer, 2007.
[11] G. Stengle, “A Nullstellensatz and a Positivstellensatz in semialgebraic
geometry,” Mathematische Annalen, vol. 207, no. 2, pp. 87–97, 1974.

[12] S. Prajna, A. Papachristodoulou, P. Seiler, and P. Parrilo, “SOSTOOLS
and its control applications,” Positive polynomials in control, pp. 580–
580, 2005.
[13] A. Ayad, “A survey on the complexity of solving algebraic systems,” in
International Mathematical Forum, vol. 5, no. 7, 2010, pp. 333–353.
[14] H. Kong, F. He, X. Song, W. N. Hung, and M. Gu, “Exponentialcondition-based barrier certificate generation for safety verification of
hybrid systems,” in CAV. Springer, 2013, pp. 242–257.
[15] S. Prajna and A. Jadbabaie, “Safety verification of hybrid systems using
barrier certificates,” HSCC, pp. 271–274, 2004.
[16] R. F. Stengel, “Flight dynamics,” Fluid Dynamics, vol. 1, 2004.
[17] K. Ghorbal and A. Platzer, “Characterizing algebraic invariants by differential radical invariants,” in Tools and Algorithms for the Construction
and Analysis of Systems. Springer, 2014, pp. 279–294.
[18] S. Sankaranarayanan, H. Sipma, and Z. Manna, “Constructing invariants
for hybrid systems,” HSCC, pp. 69–77, 2004.
[19] R. Rebiha, A. V. Moura, and N. Matringe, “Generating invariants for
non-linear hybrid systems,” TCS, vol. 594, pp. 180–200, 2015.
[20] S. Sankaranarayanan, “Automatic invariant generation for hybrid systems using ideal fixed points,” in International Conference on Hybrid
Systems: Computation and Control. ACM, 2010, pp. 221–230.
[21] N. Matringe, A. V. Moura, and R. Rebiha, “Generating invariants
for non-linear hybrid systems by linear algebraic methods,” in Static
Analysis. Springer, 2010, pp. 373–389.
[22] A. Tiwari and G. Khanna, “Nonlinear systems: Approximating reach
sets,” Hybrid Systems: Computation and Control, pp. 171–174, 2004.
[23] A. Platzer and E. Clarke, “Computing differential invariants of hybrid
systems as fixedpoints,” in CAV. Springer, 2008, pp. 176–189.
[24] E. Rodríguez-Carbonell and A. Tiwari, “Generating polynomial invariants for hybrid systems,” HSCC, pp. 590–605, 2005.
[25] T. T. Johnson and S. Mitra, “Invariant synthesis for verification of
parameterized cyber-physical systems with applications to aerospace
systems,” in AIAA Infotech at Aerospace Conference, 2013.
[26] S. Gulwani and A. Tiwari, “Constraint-based approach for analysis of
hybrid systems,” in CAV. Springer, 2008, pp. 190–203.
[27] B. Sassi, M. Amin, A. Girard, and S. Sankaranarayanan, “Iterative
computation of polyhedral invariants sets for polynomial dynamical
systems,” in CDC. IEEE, 2014, pp. 6348–6353.
[28] J. Liu, N. Zhan, and H. Zhao, “Computing semi-algebraic invariants
for polynomial dynamical systems,” in International Conference on
Embedded Software. ACM, 2011, pp. 97–106.
[29] S. Boyd, L. El Ghaoui, E. Feron, and V. Balakrishnan, Linear matrix
inequalities in system and control theory.
Society for Industrial
Mathematics, 1994, vol. 15.
[30] P. Parrilo, “Semidefinite programming relaxations for semialgebraic
problems,” Mathematical Programming, vol. 96, no. 2, pp. 293–320,
2003.

A PPENDIX
A. Proof of Proposition 1
Proof:
Let g(x) ∈ R[~x] be a polynomial satisfying Lf g ∈ hgi, we
first prove that
∀k ≥ 1 : Lkf g ∈ hgi
(19)
The proof is by induction. Since g(x) satisfies Lf g ∈ hgi, we
must have Lf g = p1 g for some p1 ∈ R[x]. Assume ∀k ≤ m :
Lkf g ∈ hgi holds, we only need to prove Lm+1
g ∈ hgi holds.
f
By assumption, we have Lm
g
∈
hgi,
that
is,
Lm
f
f g = pm g
m+1
m
for some pm ∈ R[x]. Then, Lf g = Lf Lf g = Lf pm g =
gLf pm + pm Lf g = gLf pm + pm p1 g = g(Lf pm + pm p1 ).
Therefore, the formula (19) holds.
Second, we prove that g(x) = 0 is an invariant, i.e. for any
trajectory x(t),
g(x(0)) = 0 =⇒ ∀t ≥ 0 : g(x(t)) = 0

(20)

m

g(x(t))
By the formula (19), we have d dt
|t=τ = Lm
m
f g|t=τ =
pm
g(x(τ
))
for
all
m
≥
1,
then
g(x(0))
=
0
=⇒
dm g(x(t))
|
=
0
for
all
m
≥
1.
Consider
the
Taylor
t=0
dtm
expansion of g(x(t)) at t = τ ,

g(x(t)) = g(x(τ )) +

∞
X
1 dm g(x(t))
|t=τ (t − τ )m (21)
m
m!
dt
m=1

Hence, the formula (20) holds.
Next, we prove that for any trajectory x(t),

g(x(0)) 6= 0 =⇒ ∀t ≥ 0 : g(x(t)) 6= 0

(22)

then, according to the continuity of g(x(t)), we can assert that
g(x) ∼ 0 is an invariant for each ∼∈ {<, ≤, =, ≥, >}.
The proof is by contradiction. Given an arbitrary trajectory
x(t), assume g(x(τ )) = 0 for some τ > 0 and
g(x(t)) > 0
m
g(x(t))
(or g(x(t)) < 0) for all t ∈ [0, τ ). Since d dt
|t=τ =
m
Lm
g|
=
p
g(x(τ
))
=
0
for
all
m
≥
1,
then
according
to
t=τ
m
f
the Taylor expansion (21), there must exist a δ ∈ R>0 such
that g(x(t)) = 0 for all t ∈ [τ −δ, τ +δ], which contradicts the
assumption that g(x(t)) > 0 (or g(x(t)) < 0) for all t ∈ [0, τ ).
Thus, the proposition holds.
B. Proof of Theorem 2
Proof: 1. By Definition 10, for all g(~u, ~x) ∈ Dg ,
g(~u, ~x) = 0 is an invariant for the trajectory originating from
~x0 , which means π~x0 ⊆ {~x ∈ Rn | g(~u, ~x) = 0}. Hence,
π~x0 ⊆ V(Dg ) holds.
2. Since the dimension of the hyperplane g(~u, ~x0 ) = 0 is
m, there must exist m vectors {~u1 , . . . , ~um } such that for
all ~u satisfying g(~u, ~x0 ) = 0, P
there exist scalars c1 , . . . , cm ,
m
ui . Therefore,
for evnot all zero, such that ~u =
i=1 ci ~
Pm
ery
g(~
u
,
~
x
)
∈
D
,
we
have
g(~
u
,
~
x
)
=
g(
c
~
u
x) =
g
i=1 i i , ~
Pm
c
g(~
u
,
~
x
),
hence,
g(~
u
,
~
x
)
∈
hg(~
u
,
~
x
),
.
.
.
,
g(~
u
x)i.
i
1
m, ~
i=1 i
For B = {g(~u1 , ~x), . . . , g(~um , ~x)}, we have hDg i ⊆ hBi.
The converse is trivial. Therefore, hDg i = hBi.

C. Principle of SOS programming
The principle of SOS programming is based on the fact that
a polynomial of P
degree 2k can be written as a sum-of-squares
(SOS) P (~x) =
qi (~x)2 for some polynomials qi (~x) of degree k if and only if P (~x) has a positive semidefinite quadratic
form, i.e. P (~x) = ~v (~x)M~v (~x)T , where ~v (~x) is a vector of
monomials with respect to x of degree k or less and M is a real
symmetric positive semidefinite matrix with the coefficients of
P (~x) as its entries. Therefore, the problem of finding a SOS
polynomial P (~x) can be converted to the problem of solving
a linear matrix inequality (LMI) M  0 [29], which can
be solved by semidefinite programming [30]. Currently, there
exists an efficient implementation, named SOSTOOLS [12],
for SOS programming and our implementation is based on
this tool. For details on SOS programming we refer the reader
to the literature.
D. Proof of Theorem 3
Proof: When considering g(~u, ~x) = 0 as a parameterized
hyperplane over ~u, since there is no constant term, the normal
vector of the hyperplane is (ψ1 (~xi ), . . . , ψK (~xi )). Therefore,
if ~x1 and ~x2 are in the same trajectory, they must correspond
to the same hyperplane, namely, the normal vectors of g(~u, ~x1 )
and g(~u, ~x2 ) are either all equal to ~0, which indicates the
formula (5) holds, or are proportional to one another, which
means the formula (4) holds. Moveover, if ψi (~x) ≡ 1, then
formula (6) follows from (4) immediately.
E. Proof of Theorem 4
Proof: By Theorem 1, the identity (7) holds if and only
if there exists no real solution for the following system P of
polynomial equations and inequalities
ψk (~x1 ) − ψk (~x2 ) = 0, k = 1 . . . K

G. Proof of Proposition 4
Proof: Since the polynomial (11) is an SOS, then ∀~x ∈
~
Rn : −g(~u∗ , ~x)−~µ1 · I~ ≥ 0, i.e. ∀~x ∈ Rn : −g(~u∗ , ~x) ≥ ~µ1 · I.
Since ∀~x ∈ X0 : I~  ~0 and µ
~ 1 is an SOS polynomial
vector, which implies ∀~x ∈ Rn : ~µ1  ~0, then ∀~x ∈ X0 :
−g(~u∗ , ~x) ≥ ~µ1 · I~ ≥ 0. Similarly, we can derive from the
SOS polynomial 12 that ∀~x ∈ XU : −g(~u∗ , ~x) < 0. By
Proposition 3, we can conclude that the system is safe.

H. Proof of Proposition 5
Proof: Assume Su is the set that satisfies the formulae (13), (14) and (15), by Proposition 1, every gl (~u∗l , ~x) ≥ 0
is an invariant at location l. To prove this proposition, it
is sufficient to prove that given any trajectory, say π, of
the system H, it cannot reach an unsafe state. Suppose the
infinite time interval R≥0 associated with π is divided into an
infinite
sequence of continuous time subintervals, i.e., R≥0 =
S∞
I
, where In = {t ∈ R≥0 |tn ≤ t ≤ tn+1,tn <tn+1 } is
n
n=0
the time interval that the system spent at location ρ(In ) (where
ρ(In ) returns the location corresponding to In associated to π),
we define the trajectory as π = {(ρ(In ), ~x(t))|t ∈ In , n ∈ N},
where ~x(t0 ) ∈ Init(ρ(I0 )). Then, our objective is to prove the
following assertion:
∀n ∈ N : ∀t ∈ In : gρ(In ) (~u∗ρ(In ) , ~x(t)) ≥ 0.

The basic proof idea is by induction.
Basis: n = 0. Since gρ(I0 ) (~u∗ρ(I0 ) , ~x) ≥ 0 is an invariant, it is
obvious that
∀t ∈ I0 : gρ(I0 ) (~u∗ρ(I0 ) , ~x(t)) ≥ 0.
Induction: n = k. Assume for some k,

pi1 (~x1 ) = 0, pi2 (~x2 ) = 0,
i1 = 1, . . . , l1 , i2 = l1 + 1, . . . l1 + l2

∀j ∈ [0, k] : ∀t ∈ Ij : gρ(Ij ) (~u∗ρ(Ij ) , ~x(t)) ≥ 0.
we mean to prove that

qj1 (~x1 ) ≥ 0, qj2 (~x2 ) ≥ 0,
j1 = 1, . . . , m1 , j2 = m1 + 1 . . . m1 + m2
rk1 (~x1 ) > 0, rk2 (~x2 ) > 0,
k1 = 1, . . . , n1 , k2 = n1 + 1 . . . n1 + n2 .

(24)

∀t ∈ Ik+1 : gρ(Ik+1 ) (~u∗ρ(Ik+1 ) , ~x(t)) ≥ 0.
(23)

That formula (23) has no solution indicates that there exists
no (~x1 , ~x2 ) ∈ X0 × XU such that ψi (~x1 ) = ψi (~x2 ) for all i =
1 . . . K. By the assumption K ≥ 2, we have Class(C, xi ) 6=
∅, i = 1, 2 for any (~x1 , ~x2 ) ∈ X0 × XU . By Theorem 3, the
system is safe.
F. Proof of Proposition 3
Proof: By Proposition 1, g(~u, ~x) ≥ 0 is an invariant of M
for any ~u ∈ RK . Suppose g(~u∗ , ~x) satisfies the constraints (9)
and (10) for some ~u∗ ∈ RK and ~x(t) is an arbitrary trajectory
of M such that ~x(0) ∈ X0 , we must have g(~u∗ , ~x(t)) ≥ 0 for
all t ≥ 0, which means ~x(t) cannot reach XU . Therefore, the
system is safe.

Case 1. (Discrete Transition) By the inductive assumption, we
know that
∀t ∈ Ik : gρ(Ik ) (~u∗ρ(Ik ) , ~x(t)) ≥ 0.
hence
∀t ∈ Ik : ~x(t) ∈ G(ρ(Ik ), ρ(Ik+1 )) =⇒ gρ(Ik ) (~u∗ρ(Ik ) , ~x(t)) ≥ 0
According to condition (14), ∀~x′ (t) ∈ R((l, l′ ), ~x) :
gρ(Ik+1 ) (~u∗ρ(Ik+1 ) , ~x′ (t)) ≥ 0 holds.
Case 2. (Continuous Transition)
Based on Case 1 and the fact that gρ(Ik+1 ) (~u∗ρ(Ik+1 ) , ~x) ≥ 0
is an invariant at ρ(Ik+1 ), we can conclude that ∀t ∈ Ik+1 :
gρ(Ik+1 ) (~u∗ρ(Ik+1 ) , ~x(t)) ≥ 0.
By induction, we know that the assertion (24) holds. Therefore, the system is safe.

I. Invariant clusters of Example 4
C1 = {u2 (x − y) = 0 | u2 ∈ R}, C2 = {u2 (x + y) = 0 | u2 ∈ R},

C3 = {u3 (x + z) = 0 | u3 ∈ R}, C4 = {u3 (x − z) = 0 | u3 ∈ R},
C5 = {u3 (y + z) = 0 | u3 ∈ R}, C6 = {u3 (y − z) = 0 | u3 ∈ R}.
J. Overapproximating a convex compact set
In fact, overapproximating a compact semialgebraic set
by one or multiple ellipsoids is reasonable. If a compact
semialgebraic set is overapproximated by an ellipsoid (~x −
~v )T A(~x−~v ) ≤ 1, where A is a positive definite matrix and x, v
are vectors, then we can use (x−v)T A(x−v) = 1 equivalently
instead of (~x − ~v )T A(~x − ~v ) ≤ 1 for the verification, the
reason is that any trajectory originating from the interior of
the ellipsoid must pass through some point on the surface of
the ellipsoid.
K. Additional information for Example 6
α1 = − 0.027854x21 − 0.048663y12 − 0.014007x22 − 0.031091y22
− 0.019017

α2 = − 0.0227x21 − 0.027808y12 − 0.034068x22 − 0.030096y22
− 0.00016071

α3 = − 1.8255x21 + 1.8251y12 + 1.8211x22 − 1.8282y22 + 0.014587
Prog = 1.8533x41 − 3.5741x21 y12 − 3.6098x21x22 + 3.7075x21 y22
+ 1.8738y14 + 3.688y12x22 − 3.5944y12y22 + 1.8551x42

− 3.5851x22y22 + 1.8583y24 + 0.83561x31 − 0.94703x21y1

− 0.49941x21x2 − 0.74911x21y2 + 1.4599x1y12 + 0.4202x1x22

+ 0.93272x1y22 − 1.6545y13 − 0.61178y12x2 − 0.91767y12y2
− 0.47623y1x22 − 1.0571y1y22 − 0.7495x32 − 1.1242x22y2
− 0.66212x2y22 − 0.99318y23 + 23.1976x21 + 35.9056y12

+ 20.5634x22 + 27.7403y22 + 0.5705x1 − 0.64657y1
− 0.0035356x2 − 0.0053034y2 + 9.8186
L. Invariant Cluster for Subsection V-A


mqwd2
mM z Xd2
mgθ d21
g9 (~u, ~x) = −
+
+
+ d1 +
Z
Iyy Z
Z
Z


2
mgθ d2
Iyy Xqwd2
mqvd1
Iyy qwd1
+
u5 +
+
−
Z
Z
ZM
M
2
2
(Iyy X + Iyy Z )d2
gIyy Xθ d21
Xz
+x−
−
−
Z
ZM m
ZM
gIyy Xθ d22
Iyy qvd2
Iyy Xqvd1
−
u85
−
−
ZM
M
ZM



Mθ
+ d21 + d22 u8 + q 2 − 2
u17 + u1
Iyy

