key: cord-0047713-aezpamp1
authors: Chen, Xin; Peng, Chao; Lin, Wang; Yang, Zhengfeng; Zhang, Yifang; Li, Xuandong
title: A Novel Approach for Solving the BMI Problem in Barrier Certificates Generation
date: 2020-06-13
journal: Computer Aided Verification
DOI: 10.1007/978-3-030-53288-8_29
sha: 17d41e64e6c07fdbf747dc9c5892ca21a9409f2e
doc_id: 47713
cord_uid: aezpamp1

Barrier certificates generation is widely used in verifying safety properties of hybrid systems because of the relatively low computational complexity it costs. Under sum of squares (SOS) relaxation, the problem of barrier certificate generation is equivalent to that of solving a bilinear matrix inequality (BMI) with a particular type. The paper reveals the special feature of the problem, and adopts it to build a novel computational method. The proposed method introduces a sequential iterative scheme that is able to find analytical solutions, rather than the nonlinear solving procedure to produce numerical solutions used by general BMI solvers and thus is more efficient than them. In addition, different from popular LMI solving based methods, it does not make the verification conditions more conservative, and thus reduces the risk of missing feasible solutions. Benefitting from these two appealing features, it can produce barrier certificates not amenable to existing methods, which is supported by a complexity analysis as well as the experiment on some benchmarks.

Cyber-physical systems (CPS) consists of tightly coupled physical components such as electrical, mechanical, hydraulic, and biological components and software systems. They are deeply involved in many safety-critical systems, for example, high confidence medical devices, traffic control and safety systems, advanced automotive systems and critical infrastructure control systems. Safety verification helps to ensure them not to behave dangerously.

Hybrid systems are popular models used in the verification of Cyber-physical systems, for its ability to describe interacting discrete transitions and continuous dynamics [18] . Safety verification contributes to checking safety properties by determining whether a system can evolve to some states violating desired safety properties when it starts at some initial conditions. A successful verification of a hybrid system can raise our confidence in its corresponding Cyber-physical system.

For Cyber-physical systems with real time constraints, fast verification is a vital requirement. For example, a online verification module in a monitoring system should return the result before the deadline is reached. The paper aims at fast verification of hybrid systems to satisfy the requirement of fast verification of Cyber-physical systems.

Intuitively, safety verification of hybrid systems can be performed by computing the reachable set. Reachable set computation based approaches explicitly computes either exact or approximate reachable sets corresponding to the dynamics in the model, and then compares them with unsafe regions. It has been successfully adopted in verifying behaviors of a system within a finite horizon. However, due to their intrinsic computational difficulty, approaches of this kind can hardly scale up to complex non-linear systems.

Many research efforts have been devoted to barrier certificate generation. A barrier certificate is a function, of which the zero level set separates the unsafe region from all reachable states of a system. It requires all system trajectories starting from some initial conditions fall into one side of the barrier certificate while the unsafe region resides on the other. As the existence of a barrier certificate implies that the unsafe region is not reachable, the safety verification problem can be transformed into the problem of barrier certificate generation. Compared with reachable set computation [31] , barrier certificate generation requires much less computation, since the unsafe region leads to seeking a barrier certificate. Especially, it behaves very well when a safety property concerns infinite time horizon [21, 34] .

Barrier certificate generation is a computation intensive task. A set of verification conditions corresponding to a specific type of barrier certificates is given at first. Then they are encoded into some constraints on state variables and unknown coefficients of barrier certificates of a specific type. Finally, those unknown coefficients are determined by solving the constraints [27] . Thus, how to encode verification conditions and solve them in an effective way is a critical and challenging problem in barrier certificate based verification.

Acting as the barrier between reachable states and the unsafe region, a barrier certificate should always evaluate to be nonnegative or negative accordingly in spite of what type it is. To achieve this, the most popular computational method utilizes the theory of Putinar's Positivstellensatz to derive a sum of squares (SOS) program of the barrier certificate, which results in a bilinear matrix inequality (BMI) solving problem belonging to the class of NP-hard problems [20, 21] . An effective and efficient BMI solver is a prerequisite for success in exploiting SOS relaxation based methods.

The general BMI problem can be solved by the commercial BMI solver PENBMI [14] at the cost of a very high computational complexity, where the (exterior) penalty and (interior) barrier method incorporates with the augmented Lagrangian method. To make it more tractable, the convex SOS relaxation based methods become popular. They transform the BMI problem (non-convex) to a linear matrix inequality (LMI) problem (convex) by fixing some multipliers and then solve it quickly via convex optimization such as semidefinite programming (SDP). Unfortunately, the removal of non-convexity may yield too conservative verification conditions so that the solution to the original BMI problem is invisible to the derived LMI problem.

The paper focuses on quickly solving the BMI problem derived from SOS relaxation by directly attacking the problem without relaxing it to a LMI one. Taking advantage of the special feature of the problem, that is all bilinear terms are cross ones between different parameter vectors, a sequential iterative scheme is proposed. It treats the non-convex BMI problem directly so as to avoid the loss of precision accompanied with non-convexity removing. Meanwhile, it provides much lower computational complexity than the PENBMI solver. Hence, the proposed method spends much less time in computation and has the potential to find solutions beyond the reach of existing methods.

To be specific, a feasible solution to the BMI problem can be found by a dual augmented Lagrangian iterative framework. At each iteration, the minimization over the four sets of primal variables is divided into four sequential minimization problems with respect to one set of primal variables by fixing the other three sets. On the theoretical side, we show that our method returns the feasible solution in cubic time, while the PENBMI solver in quartic time. We have developed a prototyping tool implementing the proposed method and compared it with the PENBMI solver and the LMI solver: SOSTOOLS [22] over a set of benchmarks gathered from the literature. The experiment shows that our tool is more effective than them and provides a much lower computational complexity than the PENBMI solver.

The paper is organized as follows. Section 2 describes the connection between safety verification and barrier certificate generation. Section 3 addresses how to transform the problem of barrier certificate generation into a BMI solving problem. In Sect. 4, a sequential iterative scheme is presented followed by a complexity analysis. Section 5 contains detailed examples illustrating the use of our method as well as the experiment on benchmarks. We compare with related works in Sect. 6 before concluding in Sect. 7.

Notations. Let R be the field of real number. R[x] denotes the polynomial ring with coefficients in R over variables

be the space of SOS polynomials. S n denotes the set of n × n symmetric matrices, and the notation B 0 means that the matrix B ∈ S n is positive semidefinite. A, B denotes the inner product between A and B.

A continuous dynamical system is modeled by a finite number of first-order ordinary differential equationsẋ

whereẋ denotes the derivative of x with respect to the time variable t, and f (x)

We assume that f satisfies the local Lipschitz condition, which ensures that given x = x 0 , there exists a time T > 0 and a unique function τ : [0, T ) → R n such that τ (0) = x 0 . And x(t) is called a solution of (1) that starts at a certain initial state x 0 , that is, x(0) = x 0 . Namely, x(t) is also called a trajectory of (1) from x 0 .

A continuous system over x consists of a tuple S : Θ, f , Ψ , wherein Θ ⊆ R n is a set of initial states, f is a vector field over the domain Ψ ⊆ R n .

A hybrid system is a system which exhibits mixed discrete-continuous behaviors. A popular model for representing hybrid systems is hybrid automata [1] , which combine finite state automata modeling the discrete dynamics, and differential equations modeling the continuous dynamics.

A hybrid automaton is a tuple H : L, X, F , Ψ, E, Ξ, Δ, Θ, 0 , where -L, a finite set of locations (or models); -X ⊆ R n is the continuous state space. The hybrid state space of the system is defined by X = L × X and a state is defined by ( , x) ∈ X ; -F : L → (R n → R n ), assigns to each location ∈ L a locally Lipschitz continuous vector field f ; -Ψ assigns to each location ∈ L a location condition (location invariant) Ψ ( ) ⊆ R n ; -E ⊆ L × L is a finite set of discrete transitions; -Ξ assigns to each transition e ∈ E a switching guard Ξ e ⊆ R n ; -Δ assigns to each transition e ∈ E a reset function Δ e : R n → R n ; -Θ ⊆ R n , an initial continuous state set; -0 ∈ L, the initial location. The initial state space of the system is defined by

Trajectories of hybrid systems combine continuous flows and discrete transitions. Concretely, a trajectory of H is an infinite sequence of states σ = {s 0 , s 1 , s 2 , · · · } such that

Furthermore, for each pair of consecutive state (s i , s i+1 ) ∈ σ with s i = ( i , x i ) and s i+1 = ( i+1 , x i+1 ) satisfies the following one of the two consecution conditions:

If Σ is the set of all possible trajectories of H, the reachable set is defined by R = {s|∃ς ∈ Σ : s ∈ ς}, i.e., R contains all states that are elements of at least one trajectory ς.

In this paper, we focus on semi-algebraic hybrid systems, that is, the corresponding vector fields are polynomials and the sets Θ, Ψ ( ), Ξ e , Δ e in H are semi-algebraic, represented by polynomial equations and inequalities. The semialgebraic sets Θ, Ψ ( ), Ξ e , and Δ e in Definition 2 are represented as follows:

where ∈ L, e ∈ E, θ(x), ψ (x), ρ e (x), and δ e (x ) are vectors of polynomials, and the inequalities are satisfied entry-wise. Suppose that X u assigns to each location ∈ L an unsafe region X u ( ), defined by

where ζ is a vector of polynomials. The safety specification is described over the trace of state ( , x) w.r.t. unsafe regions X u ( ).

. Given a hybrid system H : L, X, F , Ψ, E, Ξ, Δ, Θ, 0 and unsafe regions X u ( ), the safety property holds if there exist no trajectories of H starting from the initial set 0 × Θ, can evolve to any state specified by

For safety verification of hybrid systems, the notion of barrier certificates [21] plays an important role. A barrier certificate maps all the states in the reachable set R to non-negative reals and all the states in the unsafe region to negative reals, thus can be employed to prove safety of hybrid systems. However, the exact reachable set R is usually intractable for most hybrid systems. In [21] , a sufficient inductive condition for barrier certificates is defined as follows.

H for safety w.r.t. unsafe regions X u ( ) is a set of real functions {B (x)} such that, for all ∈ L and e = ( , ) ∈ E, the following conditions hold:

The problem of generating barrier certificates in Definition 4 is an infinitedimensional problem. In order to make it amenable to polynomial optimization, the barrier certificate {B (x)} should be restricted to a set of polynomials with a priori degree bound. Putinar's Positivstellensatz provides a powerful representation for polynomial positivity on semi-algebraic sets, which helps to transform the problem of barrier certificate generation into solving a semidefinite programming via SOS relaxation.

Arising from the second and third conditions of Definition 4, where the parameters of {B (x)} appear on the antecedent sides, the associated SOS representations using Putinar's Positvstellensatz form non-convex BMI constraints, yielded from the polynomial products between the barrier certificate and its polynomial multipliers.

In what follows, the procedure for transforming barrier certificate generation into BMI solving is recapped in detail. Firstly, SOS relaxation is applied to encode the entailment checking in condition (2) as an SOS program. In fact, all the conditions of Definition 4 can be expressed as a unified type, say, a polynomial is nonnegative (positive) on a semi-algebraic set, which can be characterized by Putinar's Positivstellensatz.

Let K be a basic semi-algebraic set defined by:

where g j ∈ R[x], 1 ≤ j ≤ s. Given the finite family g = {g 1 (x), . . . , g s (x)},the polynomial set defined by

is called the quadratic module generated by g.

can be represented as

where

Following Theorem 1, the existence of the representation (4) provides a sufficient and necessary condition of polynomial positivity on a semi-algebraic set K [23] . Although the number of auxiliary polynomials in the representation (4) is only one more than the number of polynomials that define K, the degree bound for σ i (x) is exponential with n and deg(f ). From a computational point of view, the method for finding the above representation has some degree of conservativeness, say, by fixing a priori much smaller degree bound D for σ i (x). Thus, a sufficient condition for the nonnegativity of the given polynomial f (x) on the semi-algebraic set K is provided as (5) ensures that a polynomial is nonnegative on a given semi-algebraic set. At this point, all conditions in Definition 4 can be derived as a unified type, i.e., polynomial nonnegativity on a semi-algebraic set. The representation (5) is used to characterize the conditions of barrier certificate generation, for they are more tractable. 

with the degree bound D, such that the following expressions:

are SOSes for each ∈ L and e ∈ E. Then {B (x)} satisfies the conditions in Definition 4, and therefore guarantees the safety of H.

Remark that a polynomial f (x) with deg(f ) = 2d is a sum of squares if and only if there exists a real symmetric and positive semidefinite matrix Q, called as the Gram matrix, such that

is the vector consisting of all the monomials of degree less than or equal to d. In view of the conditions (6) in Theorem 2, the problem of generating the barrier certificates requires introducing the auxiliary (Gram matrices) variables. In fact, the decision variables in the SOS program (6) are the coefficients of all the unknown polynomials in (6), such as B (x), σ(x), λ e (x) and the associated Gram matrices. The polynomial products, i.e., B (x)η e (x) and B (x)ν (x), derive some quadratic terms of the products of these unknown coefficients, which occur in the second and third constraints of (6). As a consequence, the problem for generating barrier certificates in Theorem 2 derives a non-convex BMI problem. We now show the transformation by a simple example.

Suppose the barrier certificate B(x) with deg(B) = 1, we predetermine its template as B(x) = u 0 + u 1 x with u 0 , u 1 ∈ R and u 1 = 0. For simplicity, here we consider the second condition in Definition 4, that is, to find B(x) which satisfies

Following the SOS relaxation in (6), we need to find B(x) such that

and

Since φ 0 (x) and φ 1 (x) must be SOSes, we have Q 0 and u 2 ≥ 0, which is equivalent to

Therefore, the requirement that φ 0 (x) and φ 1 (x) are SOSes is translated into the BMI constraint of the form

where all B i,j ∈ S 3 are constant matrices.

As illustrated in Example 1, the problem of generating barrier certificates satisfying condition (6) can be transformed into a BMI problem of the form

where all B i,j ∈ S t are constant matrices, u = [u 1 , . . . ,

are parameter coefficients of the unknown polynomials occurring in the original SOS program. Essentially, the BMI problem (9) is NP-hard. To simplify the problem considerably, the canonical approach is to swap v, corresponding to the polynomial multipliers η e (x) and ν (x), with the fixed vector. This strategy can reduce the BMI constraint into the associated LMI one. Unfortunately, the resulting LMI problem is considerably more conservative than the original BMI one. To be specific, the fixed η e (x) and ν (x) may result in too conservative verification conditions that rule out barrier certificates satisfy the non-convex conditions but not the stronger convex conditions. By investigating (9), we can find a crucial feature of B(u, v) , that is, all cross terms between parameters of u and v are of the form u i v j . The feature motivates us to design a more efficient approach for the specific type of BMI problems.

The conventional approaches for solving the BMI problem typically employ the augmented Lagrangian iterative framework, wherein each iteration involves two optimization problems for primal and dual variables. Due to the existence of nonlinear terms (quartic terms) in the associated Lagrangian function, the analytical solutions to the first problem do not exist. The iterative-based nonlinear solving procedure is introduced to obtain the numerical solutions which results in a time-consuming computing process.

Observing the BMI problem (9), we can see that all nonlinear terms are the cross terms between u and v. As a result, the associated dual augmented Lagrangian function is quartic for all variables, but is quadratic with respect to each single variable. Having this crucial feature, if we choose one variable as the independent variable and assign the others with fixed values, we may get the problem of minimizing the quadratic function. According to the firstorder optimality condition, given a quadratic function f (x), the sufficient and necessary condition thatx is a minimizer of f (x) requires that the gradient of f (x) to be zero atx, i.e., ∇f (x) = 0. As a consequence, the analytical solutions to our studied optimization problem can be easily formulated, since the gradient of the associated Lagrangian function is affine.

The analytical optimal solutions can be obtained by calling simple matrix computation, and thus are much more efficient than numerical solutions whose computation relies on complicated nonlinear optimization methods. The computational advantage is further demonstrated by a complexity analysis of our scheme against the existing BMI solving algorithm that combines the (exterior) penalty and (interior) barrier method with the augmented Lagrangian method, presented later in this section.

To utilize the computational advantage of analytical optimal solutions, for the first optimization problem (w.r.t primal variables) involved in each iteration of the augmented Lagrangian iterative framework, rather than using the usual joint minimization for all primal variables, we introduce a sequential minimization scheme, that is, dividing it into four sequential sub-optimization problems over one independent variable while keeping the others fixed. More concretely, the sub-optimization problem with one single primal variable is constructed by replacing the other variables with their optimal solutions obtained from the current iteration (if available) or the last iteration.

This section first introduces an iterative scheme to solve the BMI problem and then illustrates how to derive analytical solutions to the sub-problems in each iteration followed by a complexity analysis against the existing algorithm.

We start by presenting a straightforward reformulation of the BMI problem (9) as follows:

Clearly, there exists a feasible solution (u, v) to the BMI problem (9) if and only if the optimal value of problem (10) is non-positive, i.e., λ * ≤ 0. We try to build an iterative scheme for dealing with the optimization problem (10). The augmented Lagrangian function L associated with (10) is defined as:

where μ > 0, ·, · means the inner product operator, and · F denotes the Frobenius norm of a matrix. Let U ∈ S t be the Lagrangian multiplier associated with the equality constraint, the dual function is defined as

and the Lagrange dual problem associated with (10) is to maximize this dual function g(U ), i.e., max U g(U ). Clearly, the dual function yields lower bounds on the optimal value λ * of the problem (10), that is, g(U ) ≤ λ * for any U .

Applying the dual ascent [17] to the augment Lagrangian function yields the iterative scheme, consisting of the following updates

where the first step is the primal variables update, and the second step is the dual variable update. The first step in (12) consists of quartic terms and is lack of analytical solution. Thus, it requires jointly minimizing L μ (λ, u, v, Z, U k ) with respect to λ, u, v and Z, which can be directly solved by applying the iterative-based nonlinear optimization procedure at the cost of a high computational complexity. Instead of the usual joint minimization solving, we separate the minimization over the primal variables λ, u, v, Z into four steps, that is, λ, u, v and Z are updated in an alternating scheme, that is, minimizing L μ with respect to one primal variable given the others fixed. In detail, the sequential iterative scheme consists of the following new iterations:

The above iterative scheme introduces a sequential minimization that treats the four primal variables one by one. Benefited from the fact that the explicit formulae for the minimizer or maximizer (13) (14) (15) (16) (17) are available, the analytical solutions can be directly derived. Furthermore, as the computation of those analytical solutions involves only simple matrix computation, such as eigenvalue decomposition and matrix inverse, it will be very efficient.

In this subsection, we focus on how to find analytical solutions to problems (13) (14) (15) (16) (17) in terms of the first-order optimality conditions. Theorem 3. The minimizer λ k+1 of (13),i.e.,

has the following analytical formula:

where Tr(U k ) denotes the trace of U k .

Proof. The first-order optimality condition for (13) is

It follows that the specified λ k+1 in (18) is the optimal solution of (13), which concludes the proof.

The first-order optimality condition resembling Theorem 3 can also be invoked to produce the corresponding analytical solutions to (14) and (15), respectively.

Let u k+1 be the minimizer of (14) . Then [j] , and

Proof. The first-order optimality condition for (14) is

Then we have

for i = 1 . . . , p. Thus, ∇ u L μ (λ k+1 , u, v k , Z k , U k ) = 0 yields (19) , which proves the claim.

Let v k+1 be the minimizer of (15) . Then

where [j] , and

Proof. Similar to the proof of Theorem 4.

The theorems below demonstrate the analytical solutions to the Zminimization and U -maximization, respectively. Theorem 6. Let Z k+1 be the minimizer of (16) , and U k+1 be the solution of (17) . Denote by P k+1 the matrix P k+1 := λ k+1 I +B(u k+1 , v k+1 )−μU k . Suppose P k+1 = QΣQ T is a spectral decomposition, namely,

where Σ + and Q † are the nonnegative eigenvalues and the associated orthogonal eigenvectors, while Σ − and Q ‡ are the negative eigenvalues and the associated orthogonal eigenvectors. Then we have

Proof. The first-order optimality condition for (16) is

In view of the terms of (23), the problem (16) is translated to

which reads as

According to the spectral decomposition of P k+1 , the result (21) immediately follows.

From (17), we have

which yields the result (22).

From the above observation in Sect. 4.1 and Sect. 4.2, the detailed procedure for the sequential iterative scheme is summarized in Algorithm 1.

At the beginning of Algorithm 1, u 0 ∈ R p , v 0 ∈ R q are selected randomly, Z 0 = M 0 · M 0 where M 0 ∈ R t is chosen randomly, and heuristically U 0 = δ · I t with δ > 0.

There are several options for the stopping criterion of the loop in Algorithm 1. That is, Algorithm 1 will stop and return the current result when one of the following cases occurs:

where is a given tolerance. A reasonable value for the stopping criterion might be = 10 −6 .

Input: Problem (9); initial values u 0 , v 0 , Z 0 and U 0 . Output: A feasible solution (u * , v * ) of (9). 1 while stopping criterion not met do 2 Compute λ k+1 according (18); 3 Compute u k+1 and v k+1 according to (19) and (20), respectively;

Get the minimal eigenvalue of B k+1 , denoted byλ;

Compute Z k+1 according to (21) ; 10 Compute U k+1 according to (22) .

We analyze the complexity of Algorithm 1 and further compare it with the algorithm in PENBMI solver [14] , which combines the (exterior) penalty and (interior) barrier method with the augmented Lagrangian method. The BMI problem we study corresponds to a nonconvex optimization problem with quartic terms. For the BMI problems of the special form, neither of the two algorithms can guarantee to converge. A complete complexity analysis is not available as the number of iterations is not predictable. Therefore, the computational complexity of one iteration becomes a safe baseline for performance evaluation. In this paper, we follow the same complexity analysis as that in [14] , i.e. analyzing the complexity in one iteration.

Recall that the dimension of the matrix B(u, v) in (9) is t, and the numbers of variables u and v are p and q, respectively. We see that each iteration in Algorithm 1 can be divided into five steps. Firstly, the step of updating λ costs O(t) flops, which is carried out by 3t + 3 adds. In the step of u−update, the complexity is clearly dominated by the computation of the inverse of A u ∈ R p×p , which costs O(p 3 ) flops [5] . Analogously, v−update can be done in O(q 3 ) flops. In the step of Z−update, the critical issue is to compute the eigenvalue decomposition of matrix V k+1 ∈ R t×t , at a cost of about The total cost of each iteration in Algorithm 1 is then O(p 3 + q 3 + t 3 ), while the cost of the algorithm adopted in PENBMI is approximately O((p + q)t 3 + (p + q) 2 t 2 + (p + q) 3 ), as shown in [14] . Assume that p, q and t are bounded by T ∈ Z, i.e., T = max{p, q, t}, the complexity of Algorithm 1 is approximately O(T 3 ), whereas the complexity of PENBMI is approximately O(T 4 ).

In this section, we first show our method by verifying a nonlinear continuous system and then compare our Sequential Iterative Scheme tool: SISBMI solver with the other two solvers: PENBMI and SOSTOOLS.

Example 2. Consider the following nonlinear continuous system [28] ⎡

It is required to verify that all trajectories of the system starting from the initial set Θ = {x ∈ R 3 | (x 1 + 14.5) 2 + (x 2 + 14.5) 2 + (x 3 − 12.5) 2 ≤ 16} will never enter the unsafe region

It suffices to find a barrier certificate B(x), which satisfies all the conditions in Definition 3. Suppose that the degree of B(x) is 4, and the degree bound D = 6. Firstly, we construct a bilinear SOS program (6) , which is further transformed into a BMI problem of the form (9) where the dimension of B(u, v) is 78, and the number of decision variables is 396. By applying our algorithm, we succeed to solve the BMI problem and obtain the following barrier certificate

As shown in Fig. 1 , the zero level set of the barrier certificate B(x) (the steelblue surface) separates X u (the red ball) from all trajectories starting from Θ (the green ball). Therefore, the safety of the above system is verified.

Alternatively, by applying the PENBMI solver to compute the solution of the problem (9), we cannot find barrier certificates with degree less than 6. The system starts in location 1 with the initial set

Our task is to verify that the system will never enter the unsafe set

Applying our SISBMI solver, we obtain the polynomial barrier certificate with degree 4: In Table 1 , n denotes the number of the system variables, and |L| denotes the number of locations; d f denotes the maximal degree of the polynomials in the vector fields; t is the dimension of the matrix B(u, v), and N refers to the number of decision variables appearing in the BMI problem (9), namely, dim(u)+dim(v); d s , d p and d l denote the degrees of the barrier certificates obtained via SISBMI, PENBMI and SOSTOOLS, respectively; I s and I p are the numbers of iterations used by SISBMI and PENBMI, respectively; T s , T p and T l record the time spent by computation in seconds; the symbol-means that the solver was unable to return a feasible solution with the degree bound deg(B) ≤ 6. Table 1 shows that for the 19 examples, our SISBMI solver can successfully handle 17 of them while the numbers of successful examples of PENBMI and SOSTOOLS are 13 and 9, respectively. Our SISBMI solver seems to provide the best solving capability. There are 10 examples that can be treated by BMI solvers (either SISBMI or PENBMI) unable to be solved by the LMI solver SOSTOOLS due to the more conservative conditions in the corresponding LMI problems. To evaluate the best performance of SOSTOOLS, we have tried some widely used multipliers [16, 20] , such as 0, ±1, ±(1+x 2 1 +· · ·+x 2 n ), as well as some polynomial multipliers with random coefficients and the prior degree bound that guarantee the degrees of the polynomials involved in the verification conditions (6) do not increase. Examples C8-C9 show the case where the solver PENBMI performs better than our SISBMI solver as a result of the fact that both SISBMI and PENBMI solvers only find local optimal solutions to the BMI problems.

The above analysis on effectiveness can also be used to support that our SIS-BMI solver is a necessary complement to the existing tools. As shown in Table 1 , PENBMI solver can cover 13 examples. To solve the remaining 6 examples, it has to resort to the SISBMI solver.

Considering the efficiency, the solver SOSTOOLS performs the best for almost all the successful examples because of the much lower computational complexity for solving the relaxed LMI problems. The efficiency comparison between SISBMI and PENBMI solvers can be made by examining the ratio between the execution times of these two solvers in Table 1 . For the 11 examples that are solved by both tools, on average, our SISBMI solver costs 3.4 times than PENBMI solver in the number of iterations while only costs 0.27 times than PENBMI solver in time. That is for all the successful examples, our SISBMI solver takes much less time than PENBMI solver even it spends more iterations, which complies with the complexity analysis of the underlying algorithms. Both the theoretical analysis and the experiments support that our SISBMI solver is more efficient than PENBMI solver.

In theory, the problem of barrier certificate generation is a quantifier elimination problem. The verification conditions corresponding to a barrier certificate can be encoded into a set of constraints on state variables and coefficients where the unknown coefficients are existentially quantified and state variables are universally quantified. Hence, several symbolic computation approaches [11, 19, 29] , such as cylindrical algebraic decomposition (CAD) or Grönber bases computation, have been directly applied to attack the associated quantifier elimination problems. However, due to the high computational complexity, they suffer from the scalability problem.

Due to the relatively low computational complexity, SOS relaxation based methods become popular. Rather than directly handling quantified constraints, they transform them to a non-convex bilinear matrix inequality. Z. Yang et al. [35] relied on the BMI solver PENBMI to compute exact polynomial barrier certificates. O. Bouissou et al. [3] applied interval analysis to handle the BMI problem derived from the dynamical systems whose initial and unsafe regions are restricted to the box form. G. Jessica et al. [10] presented an augmented Lagrangian framework for the special case of bilinear programs that arise from data flow constraints and correspond to the construction of numerical abstract domains aiming at safety verification.

To alleviate its computational intractability, a convex surrogate has been proposed that behaves fairly well. Specifically, once the multipliers are fixed, the BMI problem is further transformed into a LMI problem that can be quickly solved by convex optimization. S. Prajna et al. [20] had first put the idea forward. A. Sogokon et al. [34] employed the comparison principle associated with the convex verification conditions, to generate vector barrier certificates in safety verification.

Inspired by the fact that it is the non-convex feature of verification conditions prevents well-developed convex optimization to be directly applied, many convex but stronger verification conditions are studied. H. Kong et al. [16] proposed an exponential condition for semi-algebraic hybrid systems. Kapinski et al. [12] diagnosed convex verification conditions to Lyapunov-based barrier certificates. C. Sloth et al. [32] considered convex barrier certificates associated with compositional conditions for a group of interconnected hybrid systems. L. Dai et al. [4] studied how to balance the convexity of verification conditions with the expressiveness of barrier certificates. All these convex verification conditions are equivalent forms of LMI problems. They facilitate problem-solving at the risk of losing feasible solutions.

We have presented a sequential iterative scheme for solving the BMI problem derived from the barrier certificate generation of semi-algebraic hybrid systems. Taking advantage of the special feature of the bilinear terms, the proposed approach is more efficient than the existing BMI solver. Furthermore, compared with popular LMI solving based methods, the solving procedure does not make the verification condition more conservative, and thus reduces the risk of missing solutions. In virtue of the two appealing features, our approach can produce barrier certificates not amenable to existing methods, which is evidenced by a theoretical complexity analysis as well as the experiment on some benchmarks.

The algorithmic analysis of hybrid systems

Predicate abstraction for reachability analysis of hybrid systems

Computation of parametric barrier functions for dynamical systems using interval analysis

Barrier certificates revisited

Seeking Darboux polynomials

A semiclosed-loop algorithm for the control of blood glucose levels in diabetics

Quadcopter model

Finding nonpolynomial positive invariants and Lyapunov functions for polynomial systems through Darboux polynomials

Template polyhedra and bilinear optimization

Constraint-based approach for analysis of hybrid systems

Simulationguided Lyapunov analysis for hybrid dynamical systems

Systems Biology in Practice: Concepts, Implementation and Application

PENNON: a code for convex nonlinear and semidefinite programming

Safety verification of nonlinear hybrid systems based on invariant clusters

Exponential-condition-based barrier certificate generation for safety verification of hybrid systems

Numerical Optimization

Virtual substitution & real arithmetic. Logical Foundations of Cyber-Physical Systems

Computing differential invariants of hybrid systems as fixedpoints

Safety verification of hybrid systems using barrier certificates

A framework for worst-case and stochastic safety verification using barrier certificates

SOSTOOLS: sum of squares optimization toolbox for MATLAB

Positive polynomials on compact semi-algebraic sets

Constraints for continuous reachability in the verification of hybrid systems

Safety verification of hybrid systems by constraint propagation-based abstraction refinement

Providing a basin of attraction to a target region of polynomial systems by computation of Lyapunov-like functions

Validating numerical semidefinite programming solvers for polynomial invariants

Lyapunov function synthesis using Handelman representations

Constructing invariants for hybrid systems

Stability and stabilization of polynomial dynamical systems using bernstein polynomials

State estimation of dynamical systems with unknown inputs: entropy and bit rates

Compositional safety analysis using barrier certificates

Non-linear continuous systems for safety verification (benchmark proposal)

Vector barrier certificates and comparison systems

Exact verification of hybrid systems based on bilinear SOS representation

Darboux-type barrier certificates for safety verification of nonlinear hybrid systems