Synthesizing Invariant Barrier Certificates via Difference-of-Convex Programming

A barrier certificate often serves as an inductive invariant that isolates an unsafe region from the reachable set of states, and hence is widely used in proving safety of hybrid systems possibly over the infinite time horizon. We present a novel condition on barrier certificates, termed the invariant barrier-certificate condition, that witnesses unbounded-time safety of differential dynamical systems. The proposed condition is by far the least conservative one on barrier certificates, and can be shown as the weakest possible one to attain inductive invariance. We show that discharging the invariant barrier-certificate condition -- thereby synthesizing invariant barrier certificates -- can be encoded as solving an optimization problem subject to bilinear matrix inequalities (BMIs). We further propose a synthesis algorithm based on difference-of-convex programming, which approaches a local optimum of the BMI problem via solving a series of convex optimization problems. This algorithm is incorporated in a branch-and-bound framework that searches for the global optimum in a divide-and-conquer fashion. We present a weak completeness result of our method, in the sense that a barrier certificate is guaranteed to be found (under some mild assumptions) whenever there exists an inductive invariant (in the form of a given template) that suffices to certify safety of the system. Experimental results on benchmark examples demonstrate the effectiveness and efficiency of our approach.


Introduction
Hybrid systems are mathematical models that capture the interaction between continuous physical dynamics and discrete switching behaviors, and hence are widely used in modelling cyber-physical systems (CPS).These CPS may be complex and safety-critical, with sensitive variables of the environment in its sphere of control.Everyday examples include process control at all scales, ranging from household appliances to nuclear power plants, or embedded systems in transportation domain, such as autonomous driving maneuvers in automotive, aircraft collision-avoidance protocols in avionics, or automatic train control applications, as well as a broad range of devices in health technologies, such as cardiac pacemakers.
The safety-critical feature of these CPS, with increasingly complex behaviors, has initiated automatic safety or, dually, reachability verification of hybrid systems [1,15].The problem of reachability verification is undecidable in general [1], albeit with decidable families of sub-classes (see, e.g., [2,[16][17][18]32]) identified in the literature.The hard core of the verification problem lies in reasoning about the continuous dynamics, which are often characterized by ordinary differential equations (ODEs).In particular, when nonlinearity arises in the ODEs, the explicit computation of the exact reachable set is usually intractable even for purely continuous dynamics [51].
Therefore in the literature, a plethora of approximation schemes, as surveyed in [15], for reachability analysis of hybrid systems has been developed, including an invariant-style reasoning scheme known as barrier certificate [42].A barrier certificate often serves as an inductive invariant that isolates an unsafe region from the reachable set, thereby witnessing safety of hybrid systems possibly over the infinite time horizon.A common way to synthesize barrier certificates is to reduce the condition defining barrier certificates to a numerical optimization or constraint solving problem.There is, however, a trade-off between the expressiveness of the barrier-certificate condition and the efficiency in discharging the reduced constraints.Hence, to enable efficient algorithmic synthesis of barrier certificates via, e.g., linear programming (LP), second-order cone programming (SOCP), semidefinite programming (SDP) and interval analysis [11,31,63], the general condition on inductive invariance (that a barrier certificate defines an invariant, see [8,53]) has been strengthened into a spectrum of different shapes, e.g., [8,30,53,62,64].It has been, nevertheless, a long-standing challenge to find a barrier-certificate condition that is as weak as possible while admitting efficient synthesis algorithms.
In this paper, we present a new condition on barrier certificates, termed the invariant barrier-certificate condition, based on the sufficient and necessary condition on being an inductive invariant [37].Our invariant barrier-certificate condition is by far, to the best of our knowledge, the least conservative one on barrier certificates, and can be shown as the weakest possible one to attain inductive invariance.We show, by leveraging Putinar's Positivstellensatz [33], that discharging the invariant barrier-certificate condition -thereby synthesizing invariant barrier certificates-can be encoded as solving an optimization problem subject to bilinear matrix inequalities (BMIs).We further show that general bilinear matrix-valued functions can be decomposed as a difference of two psd-convex (extension of convexity to matrix-valued functions) functions using eigendecomposition, thus resulting in a synthesis algorithm as per difference-ofconvex programming (DCP) [34,54], which solves a series of convex sub-problems (in the form of linear matrix inequalities (LMIs)) that approaches (arbitrarily close to) a local optimum of the BMI problem.This algorithm is incorporated in a branch-and-bound framework that searches for the global optimum in a divideand-conquer fashion.We present a weak completeness result of our method, in the sense that a barrier certificate is guaranteed to be found (under some mild assumptions) whenever there exists an inductive invariant (in the form of a given template) that suffices to certify the system's safety.A similar result on completeness is previously provided only by symbolic approaches, yet to the best of our knowledge, not by methods base on numerical constraint solving, e.g., [4,61,62].Experiments on a collection of examples suggested that our invariant barriercertificate condition recognizes more barrier certificates than existing conditions, and that our DCP-based algorithm is more efficient than directly solving the BMIs via off-the-shelf solvers.We defer all proofs to Appendix A.

A Bird's-Eye Perspective
We use the following example to give a bird's-eye view of our approach.
Example 1 ( overview [11]).Consider the following continuous-time dynamical system modelled by an ordinary differential equation: .
The verification obligation is to show that the system trajectory originating from any state in the initial set X 0 = {x | I(x) ≤ 0} with I(x) = x 2 1 + (x 2 − 2) 2 − 1 will never enter the unsafe set X u = {x | U(x) ≤ 0} with U(x) = x 2 + 1.
A barrier certificate satisfying our condition in Definition 4 serves as an inductive invariant that suffices to isolate the unsafe region X u from the set of reachable states from X 0 , thereby proving safety of the system over the infinite time horizon.To this end, we proceed in the following steps.
1) Encode as Sum-of-Squares (SOS) Constraints.We set a (polynomial) barrier-certificate template B(a, x) = ax 2 with unknown coefficient a ∈ R. According to Theorem 1, we only need to consider Lie derivatives up to order . By Theorem 5, B(a, x) is an invariant barrier certificate if there exists a polynomial v(x), SOS polynomials σ(x), σ (x) and a constant > 0 such that are SOS polynomials.We set = 0.01 in this example.
2) Reduce to a BMI Optimization Problem.Observe that the above SOS constraints can be formulated as BMI constraints.For instance, let us assume that (1.2) is an SOS polynomial of degree at most 2 and v(s, x) = s 0 +s 1 x 1 +s 2 x 2 is a template polynomial with unknown coefficients s.Then constraint (1.2) is equivalent to the BMI constraint meaning that the bilinear matrix (LHS of ) is negative semidefinite.Note that the bilinearity arises due to the coupling of the unknown coefficients a and s.
Constraints (1.1) and (1.3) can be reduced to BMI constraints in an analogous way 1 , yielding F 1 and F 3 .It then follows that, to solve the SOS constraints, we need to find a feasible solution (a, s) such that 2 To exploit well-developed optimization techniques, the feasibility problem ( 2) is transformed to an optimization problem subject to BMI constraints: maximize λ, a, s λ where I is the identity matrix with compatible dimensions.Note that problem (2) has a feasible solution if and only if the optimal value λ * in (3) is non-negative.
3) Decompose as Difference-of-Convex Problems.The problem (3) contains non-convex constraints and hence does not admit efficient (polynomialtime) algorithms tailored for convex optimizations.However, by our technique presented in Section 5, a non-convex function B i (λ, a, s) can be decomposed as the difference of two psd-convex (defined later) matrix-valued functions: The decomposition of B 2 (λ, a, s), for instance, gives 1 Despite that no bilinearity is involved in constraints (1.1) and (1.3), they can be processed in the same way as (1.2), yielding LMI constraints. 2Extra constraints on σ(x) and σ (x) being SOS polynomials can be encoded analogously in the feasibility problem, yet are omitted here for the sake of simplicity.
4) Solve a Series of Convex Sub-problems.Now, we apply a standard iterative procedure in difference-of-convex programming [10] as follows.Given a feasible solution z k = (λ k , a k , s k ) to the BMI optimization problem (3), the concave part −B − i (λ, a, s) in ( 4) is linearized around z k , thus yielding a series of convex programs (k = 0, 1, . ..): where DB − i denotes the derivative of the matrix-valued function B − i .The soundness of our approach asserts that the feasible set of the linearized program (5) under-approximates the feasible set of the original BMI program (3).Therefore, if λ k ≥ 0 after iteration k, we can safely claim that (a k , s k ) is a feasible solution to (2).A barrier certificate B(x) is then obtained by substituting a k in B(a, x).Moreover, if we take the optimum z * ,k of ( 5) to be the next linearization point z k+1 , the solution sequence {z k } k∈N converges to a local optimum of (3).We show that the linearized program ( 5) is equivalent to an LMI optimization problem admitting polynomial-time algorithms, say the well-known interiorpoint methods supported by most offthe-shelf SDP solvers.Our iterative procedure starts with a strictly feasible initial solution z 0 to program (3) and terminates with λ 2 ≥ 0 (subject to numerical round-off) and a 2 = −0.00363421,yielding the barrier certificate Fig. 1 depicts the system dynamics and the synthesized barrier certificate.
We remark that the aforementioned iterative procedure on solving a series of convex optimizations converges only to a local optimum of the BMI problem (3).This means that, in some cases, it may miss the global optimum that induces a non-negative λ * .We will present in Section 6 a solution to this problem by incorporating our iterative procedure into a branch-and-bound framework that searches for the global optimum in a divide-and-conquer fashion.

Mathematical Foundations
Notations.Let N, N + , R, R + and R + 0 be respectively the set of natural, positive natural, real, positive real and non-negative real numbers.For a vector x ∈ R n , x i refers to its i-th component and x denotes the 2 -norm; for a matrix A ∈ R n×m , A(i, j) refers to its (i, j)-th element.Let R[x] be the polynomial ring in We denote by Σ[x] ⊂ R[x] the set of SOS polynomials over x. S n denotes the space of n × n real, symmetric matrices.For A ∈ S n , A 0 means that A is positive semidefinite (psd, for short) 3 , i.e., ∀x ∈ R n : Differential Dynamical Systems.We consider a class of continuous dynamical systems modelled by ordinary differential equations of the autonomous type: where x ∈ R n is the state vector, ẋ denotes its temporal derivative dx/dt, with t ∈ R + 0 modelling time, and f : R n → R n is a polynomial flow field (or vector field ) that governs the evolution of the system.A polynomial vector field is local Lipschitz, and hence for some T ∈ R + ∪ {∞}, there exists a unique solution (or trajectory) ζ x0 : [0, T ) → R n originating from any initial state x 0 ∈ R n such that (1) ζ x0 (0) = x 0 , and (2) ∀τ ∈ [0, T ) : We assume in the sequel that T is the maximal instant up to which ζ x0 exists for all x 0 .Remark 1.Our techniques on synthesizing barrier certificates in this paper focus on differential dynamics of the form (6).However, we foresee no substantial difficulties in extending the results to multi-mode hybrid systems where extra constraints on the system evolution, e.g., guards, are present.Safety Verification Problem.Given a domain set X ⊆ R n , an initial set X 0 ⊆ X and an unsafe set X u ⊆ X , the reachable set of a dynamical system of the form (6) at time instant t ∈ [0, T ) is defined as R X0 (t) = {ζ x0 (t) | x 0 ∈ X 0 }.We denote by R X0 the aggregated reachable set, i.e., the union of R X0 (t) over t ∈ [0, T )4 .The system is said to be safe iff R X0 ∩ X u = ∅, and unsafe otherwise.For simplicity, we consider X = R n throughout this paper.
To avoid the explicit computation of the exact reachable set, which is usually intractable for nonlinear hybrid systems (cf., e.g., [15]), barrier-certificate methods make use of a partial differential operator, termed the Lie derivative, to capture the evolution of a barrier function along the vector field: Definition 1 (Lie derivative [29]).Given a vector field f : R n → R n over x, the Lie derivative of a polynomial function The Lie derivative L k f B(x) is essentially the k-th temporal derivative of the (barrier) function B(x), and thus captures the change of B(x) over time.
An inductive invariant Ψ ⊆ R n of a dynamical system is a set of states such that all the trajectories starting from within Ψ remain in Ψ : Definition 2 (Inductive invariant [41]).Given a system (6), a set Ψ ⊆ R n is an inductive invariant of system (6) if and only if In the sequel, we refer to inductive invariants simply as invariants.In [37], a sufficient and necessary condition on being a polynomial invariant is proposed: Theorem 1 (Invariant condition [37]).Given a polynomial B ∈ R[x], its zero sub-level set {x | B(x) ≤ 0} is an invariant of system (6) if and only if5 where N B,f ∈ N + is a completeness threshold, i.e., a finite positive integer that bounds the order of Lie derivatives, which can be computed using Gröbner bases 6 .
In contrast, a barrier certificate is a function whose zero sub-level set isolates an unsafe region X u from the reachable set R X0 w.r.t.some initial set X 0 : Definition 3 (Semantic barrier certificate [53]).Given a system (6), an initial set X 0 and an unsafe set X u , a barrier certificate of (6) is a differentiable function B : R n → R satisfying The existence of such a barrier certificate trivially implies safety of the system.Moreover, one may readily verify that if some set As observed in [53], however, the semantic statement in Definition 3 encodes merely the general principle of barrier certificates [8], yet in itself is not that useful for safety verification because it explicitly involves the system solutions.Therefore, in order to enable efficient synthesis, the semantic condition on barrier certificates has been strengthened into a handful of different shapes (see, e.g., [8,30,42,62], which all imply inductive invariance).It has been yet a long-standing challenge to find a barrier-certificate condition that is as weak as possible while admitting efficient synthesis algorithms.
Our BMI encoding of the invariant barrier-certificate condition (cf.Section 4) roots in Putinar's Positivstellensatz, which characterizes positivity of polynomials on a semi-algebraic set defined by a system of polynomial inequalities: Assume the Archimedean condition holds7 , i.e., there exists We now recall a key technique used in our reduction to semidefinite optimizations.Given a symmetric matrix X ∈ S n partitioned as An important property of the Schur complement X/A is that it characterizes the positive semidefiniteness of the block matrix X: We apply the Schur complement in Section 5 to transform nonlinear convex constraints into linear constraints.

Invariant Barrier-Certificate Condition as BMIs
In this section, we present our invariant barrier-certificate condition (see Definition 4) based on the necessary and sufficient condition on being an inductive invariant (cf.Theorem 1), and show how to encode it as BMI constraints.

Invariant Barrier-Certificate Condition
Definition 4 (Invariant barrier certificate).Given a system (6), an initial set X 0 and an unsafe set X u , a polynomial function B : R n → R is an invariant barrier certificate of system (6) if and only if Notice that the consecution constraint in Definition 4 involves Lie derivatives of orders up to N B,f ∈ N + , as is the case in Theorem 1.Our invariant barriercertificate condition hence generalizes existing conditions on barrier certificates, e.g., [4,62,65], which consider Lie derivatives only up to the first order.
The consecution condition in Definition 4 is in fact equivalent to the invariant condition (8) in Theorem 1 (cf.Lemma 2 in Appendix A), thereby revealing the relation between an inductive invariant and an invariant barrier certificate: Theorem 4 (Inductive invariance).Given a system (6), an initial set X 0 and an unsafe set It follows from Theorem 4 that our invariant barrier-certificate condition is the least conservative one on barrier certificates to attain inductive invariance.
Remark 2. We do not employ the invariant condition (8) in Theorem 1 as the constraint on the consecution of Lie derivatives.This is because our consecution condition in Definition 4 is simpler, and in particular, amenable to more straightforward transformations to SOS constraints via Putinar's Positivstellensatz, as shown later in Subsection 4.2.
Remark 3.For a fixed 0 < N < N B,f , the consecution condition in Definition 4 can be strengthened in the following way while preserving inductive invariance: where for the N-th Lie derivative, one needs L N f B(x) < 0 (rather than L N f B(x) ≤ 0).In practice, using such a strengthened consecution condition -with less subconstraints to solve-may yield more efficient synthesis.

Encoding as BMI Optimizations
Next, we show how to encode synthesizing an invariant barrier certificate (cf.Definition 4) as an optimization problem subject to BMIs.To this end, we first recast the invariant barrier-certificate condition into a collection of SOS constraints8 .Theorem 5 (Sufficient condition for invariant barrier certificate).Given a system (6), an initial set X 0 = {x | I(x) ≤ 0} and an unsafe set By enforcing the Archimedean condition and applying Putinar's Positivstellensatz, we further derive a necessary condition of invariant barrier certificate: Theorem 6 (Necessary condition for invariant barrier certificate).Given a system (6), an initial set X 0 = {x | I(x) ≤ 0} and an unsafe set are SOS polynomials.
Notice that a polynomial B(x) satisfying the sufficient condition in Theorem 5 suffices as an invariant barrier certificate that witnesses safety of the system.In contrast, a polynomial B(x) satisfying the necessary condition in Theorem 6 may serve as a candidate invariant barrier certificate, and safety of the system can be concluded via a posterior check9 of B(x) per Definition 4.
Next we show how to encode an SOS constraint of the shape "h(x) ∈ Σ[x]" in Theorems 5 and 6 as a BMI constraint.To this end, we first set a template polynomial10 B(a, x) parameterized by unknown real coefficients a as the barrier certificate.We then proceed by setting templates for the remaining unknown polynomials (e.g., v i,j (x)) and SOS polynomials (e.g., σ(x) and ρ(x)) in h(x), with all the parameters in these templates grouped into s.Observe that the parameterized SOS polynomial h(a, s, x) is of a bilinear form on the parameter spaces, i.e., h(a, s, x) is linear in a and s separately.However, nonlinearity arises in the combined parameter space (a, s) due to the product couplings of a and s, i.e., v i,j (s i,j , x)L j f B(a, x) in the consecution constraint.Now the problem of synthesizing an invariant barrier certificate boils down to searching for an instantiation of the parameters a and s such that the sufficient condition in Theorem 5 holds (or alternatively, the necessary condition in Theorem 6 holds and the posterior check passed).Such an instantiation of a (making B(a, x) an invariant barrier certificate) will be called valid in the sequel.
Suppose that a parameterized SOS polynomial h(a, s, x) is of degree at most 2d, with user-specified d ∈ N. Then h(a, s, x) can always be written in quadratic form as h(a, s, x) = b T Q(a, s)b, where b = (1, x 1 , x 2 , x 1 x 2 , . . ., x d n ) is the basis vector of size p = n+d n containing all monomials of degree up to d, and Q(a, s) ∈ S p is a parameterized real symmetric matrix known as the Gram matrix [6] 11 .An important fact states that h(a, s, x) is SOS if and only if Q(a, s) 0.
Let F(a, s) = −Q(a, s).As per h(a, s, x), the matrix-valued function F(a, s) is bilinear in (a, s).Observe that h(a, s, x) is SOS if and only if the BMI constraint F(a, s) 0 holds.See Example 1 for an illustration of this BMI encoding.
In general, F(a, s) can be flattened in an expanded bilinear form as where m and n are the size of a and s, respectively; F, H i , G j , F i,j ∈ S p are constant matrices.Discharging the conditions of invariant barrier certificates hence amounts to solving the BMI feasibility problem of finding a and s s.t.
Here F(a, s) is indexed by ι and l is the number of SOS constraints involved.
To exploit well-developed techniques in optimization, the feasibility problem ( 10) is transformed to an optimization problem subject to BMI constraints: maximize λ, a, s λ A solution (λ, a, s) to ( 11) is feasible if it satisfies the BMIs in (11), and strictly feasible if all the BMIs are satisfied with strict inequalities.We sometimes drop the λ component in the solution when it is clear from the context.Notice that problem (10) has a feasible solution if and only if the optimal value λ * in the BMI optimization problem (11) is non-negative.
To achieve (weak) completeness of our method in subsequent sections on solving the BMI optimization problem, we make the following assumption on the boundedness of the search space (a, s) of the optimization.
Assumption 1 (Boundedness on the parameters) Every feasible solution (a, s) to the BMI problem (11) is in a compact set with non-empty interior, i.e., Remark 4. The boundedness on a in Assumption 1 makes sense in practice since we usually prefer barrier certificates with bounded coefficients.Moreover, when the bilinear functions F ι (a, s) in (11) are affine in a and s, i.e., with a zero constant matrix F , the parameters a and s can be scaled independently by any positive factor.Therefore in this case, w.l.o.g, one may simply set L a = L s = 1.

Solving BMI Optimizations via DCP
The BMI optimization problem (11), derived from the synthesis problem, is known to be NP-hard and contains non-convex constraints [57], and hence is not amenable to efficient (polynomial-time) algorithms committed to solving convex optimizations.In this section, we present an algorithm for solving general BMI optimizations via difference-of-convex programming [34,54], which solves a series of convex sub-problems that approaches a local optimum of (11).
For brevity, we consider optimization problems with a single BMI constraint 12 : where the objective function g : R m+n → R is linear in z = (x, y); F, H i , G j , F i,j ∈ S p are constant symmetric matrices.

Difference-of-Convex Decomposition
The key challenge in solving the BMI problem ( 12) is its non-convexity, that is, the matrix-valued function B(x, y) is, in general, not psd-convex.
There have been attempts, most pertinently in [10], to decompose a bilinear function as a difference between two psd-convex functions, known as the difference-of-convex (DC) decomposition, such that the optimization in its decomposed form enjoys well-established techniques in difference-of-convex programming [34,54].The DC decomposition in [10], however, is confined to BMIs of a specific structure, namely, X T Y + Y T X 0, where X and Y are matrix variables containing variables x i and y j , respectively.The more general bilinear function B(x, y) in (12) does unfortunately not admit straightforward forms of decomposition such as those in [10, Lemma 3.1].
In what follows, we present a difference-of-convex decomposition of the matrixvalued function B(x, y), inspired by [59], using eigendecomposition.
First, observe that the function B(x, y) can be written as where ⊗ denotes the Kronecker product: for two matrices represents the zero matrices with compatible dimensions, and The form of (13) implies that B(x, y) is psd-convex if the matrix M = 0 Γ Γ T 0 is positive semidefinite.Unfortunately, as [59, Theorem 1] points out, for a non-trivial bilinear function B(x, y), M may not be positive semidefinite.
Nevertheless, the matrix M can always be decomposed as M = M 1 − M 2 with M 1 , M 2 0, i.e., a difference between two psd-matrices.One way to do so is to use the eigendecomposition of the (real symmetric 13 ) matrix M ∈ S (m+n)p .That is, M = V T DV , where the orthogonal matrix V contains the eigenvectors of M ; D is a diagonal matrix whose diagonal elements are the eigenvalues of M .
Let D + be the matrix obtained by setting all negative elements of D to zero and D − = D + − D. We have It follows that M 1 , M 2 0 and therefore we find a DC decomposition of B(x, y): Theorem 7 (Difference-of-convex decomposition).The following form where is a difference-of-convex decomposition of B(x, y).Namely, the matrix-valued functions B + (x, y) and B − (x, y) are psd-convex on R m+n .
Remark 5.In practice, the aforementioned matrices M , M 1 and M 2 induced by eigendecomposition are often highly sparse.One can hence exploit the sparsity to improve the algorithmic performance of the DCP-based synthesis approach.

Reduction to LMIs
On top of the DC decomposition (cf.Theorem 7), we can now apply a standard iterative procedure in difference-of-convex programming [10] to solve the BMIs.The core idea of the procedure is to iteratively solve a series of convex subproblems.More specifically, given a feasible solution z k = (x k , y k ) to the BMI optimization problem (12), the "concave part" −B − (x, y) in ( 14) is linearized around z k , thereby yielding a series of convex programs (k = 0, 1, . ..): where DB − (z) : R m+n → S p is the derivative of the matrix-valued function B − at z, i.e., a linear mapping from a vector u ∈ R m+n to a matrix in S p : 13 M thus only has real eigenvalues.S keeps track of visited points An extra regularization term 1  2 δ z − z k 2 with δ < 0 is added in (15) to enforce that g(z) strictly increases after each iteration until it stabilizes, which can be encoded as a second-order cone constraint and embedded in SDP solving.
Note that the linearized problem ( 15) is convex and therefore can be solved efficiently 14 via methods including, among others, augmented Lagrangian methods [36] and gradient descent methods [3].Furthermore, the Schur complement in Theorem 3 implies that (15) can be reformulated as an LMI problem: Theorem 8.The quadratic matrix inequality (QMI) constraint where N is the square root matrix of M 1 , i.e., M 1 = N T N , and Theorem 8 entails that the series of linearized convex sub-problems of the form (15) can be solved alternatively by most off-the-shelf SDP solvers designated for discharging LMIs via polynomial-time algorithms, say the interiorpoint methods.Furthermore, by taking the optimum of the k-th sub-problem to be the next linearization point z k+1 , we obtain an iterative procedure for solving general BMIs, as depicted in Algorithm 1.
Algorithm 1 falls into the DCP framework [10] and thus enjoys useful properties, e.g., soundness, termination and convergence as follows.
Theorem 9 (Soundness).Every solution z i = (x i , y i ) ∈ S with i = 0, . . ., k returned by Algorithm 1 is a feasible solution to the original BMI problem (12).
The result below states termination and convergence of Algorithm 1 in terms of KKT points of (12), i.e., solutions fulfilling the KKT conditions [3] of (12) 16 .
We remark that, under some sufficient KKT conditions and regularity conditions [3], a KKT point suffices as a local optimum.In this case, the infinite sequence {z i } i∈N of points visited by Algorithm 1 (for ε = 0) converges to a local optimum of (12).

Finding the Initial Solution
The iterative procedure in Algorithm 1 starts with a fed-by-oracle strictly feasible initial solution z 0 to the BMI problem (12).Finding such an initial solution, however, is non-trivial in general due to the non-convexity of (12).We argue though, that a strictly feasible initial solution can be obtained for the BMI problem of the form (11) induced by the barrier-certificate synthesis problem.
Recall that in the BMI problem (11), bilinearity arises from the multiplication of B(a, x) with some unknown multiplier polynomials parameterized by s.One way to reduce the BMI constraints to LMIs is to fix every multiplier polynomial to be a non-negative constant, thereby yielding a linear program: maximize λ, a λ subject to F ι (a, s) s=(cι,0,...,0) + λI 0, ι = 1, 2, . . ., l where s in F ι (a, s) is substituted by (c ι , 0, . . ., 0) with c ι ∈ R + 0 , which encodes a non-negative constant multiplier polynomial.Observe that no s-variable is involved in (16) and the constraints therein are linear in a.
As a consequence, a strictly feasible solution to the BMI problem (11) can be obtained by solving the LMI problem (16).In fact, when considering Lie derivatives only up to the first order, solving (the feasibility counterpart of) ( 16) is exactly the procedure to synthesize either an exponential barrier certificate [30] (with c ι ∈ R + ) or a convex barrier certificate [42] (with c ι = 0).Algorithm 1 therefore subsumes existing synthesis techniques in the sense that any valid barrier certificate synthesized by methods in [30,42] can also be discovered by Algorithm 1.Moreover, an alternative way to reduce the BMI constraints to LMIs is to fix the multipliers to be some given non-trivial (SOS) polynomials [64].
Algorithm 2: Branch-and-Bound: Searching for a valid parameter ā input: A BMI optimization problem of the form (11)  6 Incorporating in a Branch-and-Bound Framework The aforementioned iterative procedure on solving a series of convex optimizations converges only to a local optimum of the BMI problem (11) (or more generally, (12)).This means that, in some cases, it may miss the global optimum that induces a non-negative λ * .We present in this section a solution to this problem by incorporating the iterative procedure into a branch-and-bound framework that searches for the global optimum in a divide-and-conquer fashion, as is a common technique in non-convex optimizations.The basic idea is as follows.We first try to solve the BMI problem (11) by Algorithm 1 over the compact parameter space (C a , C s ).If a valid solution, (i.e, a solution that contains a valid parameter ā ∈ C a such that B(ā, x) is an invariant barrier certificate) is found, then the corresponding barrier certificate can be obtained.Otherwise, we keep bisecting C a and apply Algorithm 1 over each bisection 17 .The procedure, as depicted in Algorithm 2 in a recursive manner, terminates when a valid parameter is found or the partition is fine enough.
Algorithm 2 takes as input a BMI problem of the form (11) that encodes either the sufficient condition in Theorem 5 or the necessary condition in Theorem 6 for invariant barrier certificates.In the former case, a sample-and-check process (Line 2-3) is necessary to attain (weak) completeness (see Theorem 11).The conditional statement in Line 4 rules out parameter (sub-)spaces that have already been explored, which is the case when the projection of some visited point in S glb (a global set that keeps track of visited points by Algorithm 1, initialized as ∅) onto a is in the current parameter space.
The following theorem claims a weak completeness result: our method guarantees to find a barrier certificate when there exists an inductive invariant (in the form of a given template) that suffices to certify safety of the system.
Theorem 11 (Weak Completeness).Algorithm 2 returns a valid parameter ā ∈ C a , if (1) the partition granularity is fine enough (i.e., small enough η ∈ R + ), (2) the degrees of multiplier polynomials and SOS polynomials used to form (11) are large enough, and (3) there exists, for the given template B(a, x), a strictly valid parameter â ∈ C a (i.e., any parameter in some neighborhood of â is valid).Remark 7. The bisection operation in Algorithm 2 induces -in the worst casean exponential blow-up in the number of branches.In practice, one can prune branches inducing only negative objective values, via, e.g., convex relaxation [27].

Experimental Results
We have carried out a prototypical implementation 18 of our synthesis techniques in Wolfram Mathematica, which was selected due to its built-in primitives for SDP, polynomial algebra and matrix operations.Given a safety verification problem as input, our implementation works toward discovering an invariant barrier certificate (in the form of a given template) that witnesses unbounded-time safety of the system.A collection of benchmark examples (detailed in Appendix B) has been evaluated on a 2.10GHz Intel Xeon processor with 376GB RAM running 64-bit CentOS Linux 7.
Table 1 reports the empirical results.BMI-DC concerns our locally-convergent Algorithm 1 for solving BMIs (encoding the sufficient condition in Theorem 5) based on DC decomposition.We compare our approach with PENLAB [14] -an off-the-shelf solver in Matlab for directly discharging the same BMI problems (with no guarantee on convergence)-and SOSTOOLS [40] -for solving LMIs derived from Prajna and Jadbabaie's original barrier-certificate condition [42].The comparison is performed under the same problem configurations 19 .Due to numerical errors caused by floating-point computations and the fact that reaching the local/global optimum does not necessarily yield a valid barrier certificate, we additionally perform a posterior check, via both the quantifierelimination procedure in Mathematica and the SMT solver Z3 [38], of the synthesized candidate barrier certificate per Definition 4.
Table 1 shows that BMI-DC suffices to synthesize valid barrier certificates in most of the examples within a reasonable number of iterations (i.e., the number of convex sub-problems solved by SDP).This however does not cover all the cases: for the focus example, the solution is close enough to a local optimum (after 100 iterations) but yields still an invalid barrier certificate.This problem can be solved (if there exists an invariant barrier certificate as specified) by enforcing the branch-and-bound framework as presented in Section 6.The phase portraits of a selected set of examples and the synthesized invariant barrier certificates are depicted in Fig. 2 (see more in Appendix B).
The comparison in Table 1 suggests that (1) Our invariant barrier-certificate condition recognizes more barrier certificates than the original (more conservative) condition as implemented in SOSTOOLS.In particular, the lie-high-order example does admit an inductive invariant in the form of the given template, but none of the existing barrier-certificate conditions [4,62,65] -concerning Lie derivatives only up to the first order-recognizes it, since we have L 1 f B(x) = 0 for some x on the boundary of B and hence it requires to exploit the second-order Lie derivative L 2 f B; (2) Our DCP-based synthesis algorithm finds more barrier certificates in less time than directly solving the BMI problems via non-convex optimization techniques as implemented in PENLAB.
We remark that symbolic methods based on, e.g., quantifier elimination [37], can hardly deal with any of the examples listed in Table 1 due to the prohibitively high computation complexity.Moreover, it would be desirable to pursue a comparison with the augmented Lagrangian method for solving BMIs as proposed in [4], which unfortunately is not yet possible due to the unavailability of the implementation thereof.We will discuss crucial differences to [4] in Section 8.

Related Work
As surveyed in [15], the research community has, over the past three decades, extensively addressed the automatic verification of safety-critical hybrid systems.The almost universal undecidability of the unbounded-time reachability problem [1], however, confines the sound key-press routines to either semi-decision procedures or approximation schemes, most of which address bounded-time verification by, e.g., computing the finite-time image of a set of initial states.
Invariant generation [26], amongst others, is a well-established approximation scheme that provides a reliable witness for safety (or equivalently, unreachability) of dynamical systems over the infinite time horizon.Invariants can be constructed in various forms, e.g., barrier certificates [42,53] and differential invariants [37,41].With a priori specified templates, the invariant synthesis problem can be reduced to numerical optimizations or constraint solving, as in, e.g., [22,25,47,56].
Most pertinently, Prajna and Jadbabaie proposed in their seminal work [42] a concept coined barrier certificate to encode invariants.To enable efficient synthesis via semidefinite programming, the barrier-certificate condition in [42] strengthens the general condition encoding inductive invariance.Since then, significant efforts have been investigated in developing more relaxed (i.e., weaker) forms of barrier-certificate condition that still admit efficient synthesis, thereby leading to, e.g., exponential-type barrier certificates [30], Darboux-type barrier certificates [64], general barrier certificates [8] and vector barrier certificates [53].To attain efficient synthesis, these barrier-certificate conditions share a common property on convexity.That is, if for some a 1 , a 2 ∈ R m , B(a 1 , x) and B(a 2 , x) both satisfy the barrier-certificate condition, then for any 0 < µ < 1, B(µa 1 + (1 − µ)a 2 , x) must also satisfy the barrier-certificate condition.
However, neither the semantic barrier-certificate condition (9) encoding the general principle of barrier certificates [8,53] nor the inductive invariant condition ( 8) is convex.This means, when resorting to convex barrier-certificate conditions, one may miss some potential barrier certificates that suffice as inductive invariants witnessing safety.Therefore, non-convex conditions were suggested [62], for which the synthesis problem can be reduced to BMI problems solvable via customized schemes, e.g., the augmented Lagrangian method [4] and the alternating minimization algorithm [65].Our synthesis techniques also exploit a BMI reduction, with three crucial differences: (1) our invariant barriercertificate condition is equivalent to the inductive invariant condition in the sense of Theorem 4, and thus is less conservative than all the aforementioned conditions which consider Lie derivatives only up to the first order; (2) our DCP-based techniques for solving BMIs naturally inherit appealing results on convergence and (weak) completeness, which are not (and can hardly be) provided by the approaches in [4,62,65]; (3) our DCP-based iterative procedure visits only feasible solutions to the original BMI problem, and hence whenever a solution that induces a non-negative objective value is found, we can safely terminate the algorithm and claim a feasible solution to the original BMI problem, which may yield a valid barrier certificate.This is not the case for the approaches in [4,62,65].
Beyond barrier certificates, Wang and Rajamani [59] investigated the feasibility problem of general BMI problems with an application to multi-objective nonlinear observer design.The technique of eigendecomposition was also used therein to conduct the DC decomposition.The decomposed concave part, however, is simply ignored and no iterative procedure that exhibits convergence to a local optimum can be provided.
The idea of augmenting a locally-convergent algorithm with a branch-andbound framework to find the global optimum has been exploited in the realm of optimization [20] and control [58].In contrast, our method is designed for the specific problem of barrier-certificate synthesis, and hence our branch-andbound algorithm concerns only the parameter space of a, i.e., coefficients of the template barrier certificate.

Conclusion
Barrier certificates are powerful tools to prove time-unbounded safety of hybrid systems.We have presented a new condition on barrier certificates -the invariant barrier-certificate condition.This condition is by far the least conservative one on barrier certificates, and can be shown as the weakest possible one to attain inductive invariance.We showed that our invariant barrier-certificate condition can be reformulated as an optimization problem subject to bilinear matrix inequalities, which can be solved by our locally-convergent algorithm based on difference-of-convex programming.By incorporating this algorithm into a branch-and-bound framework, we obtained a weak completeness result.Experiments on benchmark examples suggested that our invariant barrier-certificate condition recognizes more barrier certificates than existing conditions, and that our DCP-based algorithm is more efficient than directly solving the BMIs via off-the-shelf solvers.
We stress that our techniques for solving BMIs are of a general nature rather than being confined to barrier-certificate synthesis.Interesting future directions include to extend our method to other synthesis problems, e.g., discovering invariants and/or termination proofs of deterministic/probabilistic programs.

Appendix A Proofs of Lemmas and Theorems
Lemma 2 (Equivalence of Lie consecution).The consecution condition in Definition 4 holds if and only if the invariant condition (8) in Theorem 1 holds.
Proof.We prove both the "if" and the "only if" part by contradiction.
For the "if" part, suppose that the invariant condition (8) holds but the consecution condition is invalid.The latter implies that for some x 0 ∈ R n and 1 Note that (17) implies B(x 0 ) = 0. From (8), it follows that either holds, or there exists 0 holds.However, - cannot hold as (17).
For the "only if" direction, suppose that the consecution condition in Definition 4 holds but the invariant condition ( 8) is invalid.The latter implies that there exists x 0 such that B(x 0 ) ≤ 0 and holds for any 0 ≤ i ≤ N B,f .For i = 0, (20) yields that B(x 0 ) ≥ 0. Together with the premise B(x 0 ) ≤ 0, we have B(x 0 ) = L 0 f B(x 0 ) = 0. Now, by taking the case i = 1 in the consecution condition, we deduce L 1 f B(x 0 ) ≤ 0. Meanwhile, for i = 1, (20) yields L 1 f B(x 0 ) ≥ 0. It thus follows that L 1 f B(x 0 ) = 0. Analogously, by taking i = 2, . . ., N B,f , we conclude L i f B(x 0 ) = 0 for all 0 ≤ i ≤ N B,f .This is exactly encoded in (8) (the rightmost conjunctive clause) and hence contradicts the assumption that (8) is invalid.Therefore, the consecution condition implies (8).
Proof (of Theorem 4).The claim is an immediate consequence of Lemma 2.
Proof (of Theorem 5).It can be shown that the k-th condition in Theorem 5 implies the k-th condition in Definition 4, for k = 1, 2, 3.For instance, the first condition in Theorem 5 requires that −B(x) + σ(x)I(x) is an SOS polynomial (and thus non-negative), we therefore have B(x) ≤ σ(x)I(x).It follows that ∀x ∈ {x | I(x) ≤ 0} : B(x) ≤ 0. A similar argument applies to the other two conditions.
Proof (of Theorem 6).The invariant barrier-certificate condition in Definition 4 characterizes positivity of polynomials over certain sets.By adding a "ball" constraint x 2 −L ≤ 0 to those sets (thus achieving the Archimedean condition), we can apply Putinar's Positivstellensatz to rewrite those polynomials into SOS forms.
Proof (of Theorem 8).Note that the square root matrix N of M 1 exists since M 1 020 .The claim then follows immediately by applying the Schur complement in Theorem 3.
Proof (of Theorem 9).We prove by induction on i.The base case holds as z 0 is assumed to be a feasible solution to (12).For the induction step, we show that z i+1 is a feasible solution to (12) if z i is a feasible solution to (12).Since z i+1 is a feasible solution to (15) linearized at z i , it suffices to show that the feasible set of ( 15) is a subset of the feasible set of (12).
Theorem 7 shows that B − (z) is psd-convex, then by [10, Lemma 2.2 (b)], we have In the meantime, z i is a feasible solution to (15) and thus fulfils Combining ( 21) and ( 22), we have B(x, y) = B + (z) − B − (z) 0 which is exactly the BMI constraint of ( 12).This completes the proof.
Proof (of Theorem 10).Let S = {z i } i∈N be the infinite sequence of visited points for ε = 0.
We first show that (2) implies (1).Assume that (2) holds, i.e., S converges (to a KKT point of ( 12)), then by Cauchy's criterion for convergence, we have ∀ε It then remains to show that S converges to a KKT point of ( 12) if the set of KKT points of ( 12) is finite.This is in fact a straightforward corollary of [10,Theorem 4.3], by noticing that the assumptions thereof can be readily verified.For simplicity, we highlight the validity of only a few of these assumptions: Since z 0 in Algorithm 1 is a strictly feasible solution to (12), the relative interior of the feasible set of ( 12) is non-empty and thus Assumption A1 in [10] holds; Our Assumption 1 on the boundedness of the search space ensures that g(z) in ( 12) is bounded from above over a bounded feasible set, and therefore the boundedness assumptions in [10,Theorem 4.3] holds.
Furthermore, under Assumption 1 on the boundedness of parameter a ∈ C a , Λ a can be shown to be bounded by the well-known Gershgorin circle theorem.
Therefore, by taking an interior point of C a as ã, and λ = Λ ã − for some ∈ R + , we obtain a strictly feasible solution ( λ, ã) to program (16).
Proof (of Theorem 11).When the assumptions (1)-(3) hold, Algorithm 2 will eventually visit a branch wherein any parameter is valid (in case a valid parameter has not been found yet).If the necessary condition for invariant barrier certificates in Theorem 6 is used to form the BMI problem (11), Line 7 ensures to return a valid parameter ā ∈ C a ; Otherwise if the BMI problem (11) encodes the sufficient condition in Theorem 5 which strengthens the invariant barrier-certificate condition in Definition 4, a valid parameter ā may not induce a non-negative objective value of (11).In this case, however, any parameter sampled and returned by Line 2-3 in the branch is valid, as it contains only valid parameters.Example 2 ( contrived).The vector flow field is:

Appendix B Benchmark Examples
-B(a, x) includes all monomials up to degree 2.
-B(a, x) includes all monomials up to degree 1.

Fig. 1 :
Fig. 1: Phase portrait of the system in Example 1.The arrows indicate the vector field and the solid curves are randomly sampled trajectories.

Fig. 2 :
Fig. 2: Phase portraits of a selected set of examples with the synthesized invariant barrier certificates.The arrows indicate the vector field (hidden in 3D-graphics for a presentation) and the solid curves are randomly sampled trajectories.

Fig. 3 :
Fig. 3: Phase portraits of another selected set of examples with the synthesized invariant barrier certificates.The arrows indicate the vector field (hidden in 3D-graphics for a clear presentation) and the solid curves are randomly sampled trajectories.
with Ca = {a | a 2 ≤ La}.output: A valid parameter ā, or otherwise ⊥ indicating a failure. 1 if La < η then return ⊥; abort on fine-enough partitions (η ∈ R + ) /* sample-and-check is not necessary if Theorem 6 is used */ 2 ā ← a randomly-sampled point in Ca; 3 if ā is valid then return ā; check validity (inductive invariance) 4 if proj a (S glb ) ∩ Ca = ∅ then S glb contains a global set of visited points Remark 6. Different choices of the multiplier constants c ι in (16) may lead to different initial solutions fed to Algorithm 1, thereby considerably different number of iterations until termination.In practice, techniques like randomization are worth exploring when choosing these multiplier constants.

Table 1 :
Empirical results on benchmark examples (time in seconds) nsys: system dimension; d flow : maximal flow-field degree; d BC : degree of the template barrier certificate.#iter.: number of iterations.0 means that the initial solution (cf.Subsection 5.3) is valid.verified: the synthesized barrier certificate is valid (), invalid () or inconclusive (?, beyond the capability of quantifier elimination in Mathematica and nonlinear reasoning in Z3).time: CPU-time, excluding that for casting the BMIs/LMIs.Boldface marks the winner among 's.