Adel Djaballah
email: adel.djaballah@ensta-paristech.fr

Alexandre Chapoutot
email: alexandre.chapoutot@ensta-paristech.fr

Michel Kieffer
email: michel.kieffer@lss.supelec.fr

Olivier Bouissou
email: olivier.bouissou@cea.fr

Construction of parametric barrier functions for dynamical systems using interval analysis

Keywords: Formal verification, Dynamic systems, Intervals, Constraint satisfaction problem

ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Formal verification aims at proving that a certain behavior or property is fulfilled by a system. Verifying, e.g., the safety property for a system consists in ensuring that it will never reach a dangerous or an unwanted configuration. Safety verification is usually translated into a reachability analysis problem [START_REF] Asarin | Approximate reachability analysis of piecewise-linear dynamical systems[END_REF][START_REF] Chutinan | Verification of polyhedral-invariant hybrid automata using polygonal flow pipe approximations[END_REF][START_REF] Frehse | SpaceEx: Scalable verification of hybrid systems[END_REF][START_REF] Sun | Controllability and reachability criteria for switched linear systems[END_REF][START_REF] Tiwari | Approximate reachability for linear systems[END_REF]. Starting from an initial region, a system must not reach some unsafe region. Different methods have been considered to address this problem. One may explicitly compute the reachable region and determine whether the system reaches the unsafe region [START_REF] Gulwani | Constraint-based approach for analysis of hybrid systems[END_REF]. An alternative idea is to compute an invariant for the system, i.e., a region in which the system is guaranteed to stay [START_REF] Chutinan | Verification of polyhedral-invariant hybrid automata using polygonal flow pipe approximations[END_REF]. This paper considers a class of invariants determined by barrier functions.

A barrier function [START_REF] Prajna | Barrier certificates for nonlinear model validation[END_REF][START_REF] Prajna | Safety verification of hybrid systems using barrier certificates[END_REF] partitions the state space and isolates an unsafe region from the part of the state space containing the initial region. In [START_REF] Prajna | Safety verification of hybrid systems using barrier certificates[END_REF] polynomial barriers are considered for polynomial systems and semi-definite programming is used to find satisfying barrier functions. Our aim is to extend the class of considered problems to non-polynomial systems and to non-polynomial barriers. This paper focuses on continuous-time systems.

The design of a barrier function is formulated as a quantified constraints satisfaction problem (QCSP) [START_REF] Benhamou | Universally quantified interval constraints[END_REF][START_REF] Ratschan | Efficient solving of quantified inequality constraints over the real numbers[END_REF]. Interval analysis is then used to find the parameters of a barrier function such that the QCSP is satisfied. More specifically, the algorithm presented in [START_REF] Jaulin | Guaranteed tuning, with application to robust control and motion planning[END_REF] for robust controller design is adapted and supplemented with some of the pruning schemes found in [START_REF] Chabert | Contractor programming[END_REF] to solve the QCSP associated to the barrier function design.

The paper is organized as follows. Section 2 introduces some related work. Section 3 defines the notion of barrier functions and formulates the design of barrier functions as a QCSP. Section 4 presents the framework developed to solve the QCSP. Design examples are presented in Section 5. Section 6 concludes the work.

In what follows small italic letters x represent real variables while real vectors x are in bold. Intervals [x] and interval vectors (boxes) [x] are represented between brackets. We denote by IR the set of closed intervals over R, the set of real numbers. Data structures or sets S are in upper-case calligraphic. The derivative of a function x with respect to time t is denoted by ẋ.

Related work

The verification of the safety properties for dynamical systems has been an active field of research in the last years. This related work review focuses on methods involving the computation of invariants for dynamical systems. Alternative methods based on the computation of reachable sets are described in [START_REF] Chen | Taylor model flowpipe construction for non-linear hybrid systems[END_REF][START_REF] Frehse | SpaceEx: Scalable verification of hybrid systems[END_REF] and in the references therein.

An invariant is a part of the state space in which the state of a dynamical system can be proved to remain. Invariants are very useful to prove the safety of dynamical systems. If an invariant does not contain unsafe regions, then the dynamical system is safe. Methods to characterize invariants have been intensively studied for linear and polynomial dynamics but much work has still to be done for non-linear dynamical systems.

Coming from the community interested in hybrid systems, a set of methods has been defined to compute invariants for various classes of systems, for example for linear or affine systems [START_REF] Tiwari | Approximate reachability for linear systems[END_REF] or for polynomial systems [START_REF] Gulwani | Constraint-based approach for analysis of hybrid systems[END_REF][START_REF] Kapinski | Simulation-guided lyapunov analysis for hybrid dynamical systems[END_REF][START_REF] Sankaranarayanan | Constructing invariants for hybrid systems[END_REF][START_REF] Yang | Exact safety verification of hybrid systems based on bilinear SoS representation[END_REF]. These methods introduce a candidate parametric function, which parameter vector has to be adjusted to define an invariant of the considered dynamical system. Various techniques are then employed to determine satisfying parameter vectors. For example, in [START_REF] Sankaranarayanan | Constructing invariants for hybrid systems[END_REF], theory of ideal over polynomials and Gröbner bases are used to define constraints to be satisfied by the parameter vector to be found. These constraints are then solved numerically using tools such as those introduced in [START_REF] Collins | Partial cylindrical algebraic decomposition for quantifier elimination[END_REF]. In [START_REF] Gulwani | Constraint-based approach for analysis of hybrid systems[END_REF], quantified polynomial constraints are introduced for the design of parameters. Then, satisfying parameter vectors are found using Farkas' Lemma and solvers from sat-modulo theory [START_REF] Barrett | Satisfiability modulo theories[END_REF]. Sum-of-Squares (SoS) polynomials are used in [START_REF] Kapinski | Simulation-guided lyapunov analysis for hybrid dynamical systems[END_REF]. The design involves various system simulations and selection of candidate parameter vectors using linear programming. A final validation of the selected parameter vectors is then performed with Mathematica and using interval analysis with dReal [START_REF] Gao | dReal: An SMT solver for nonlinear theories over the reals[END_REF]. Note that our algorithm presented in Section 4.3 could aso be used as validation method for the approach presented in [START_REF] Kapinski | Simulation-guided lyapunov analysis for hybrid dynamical systems[END_REF]. Bilinear SoS programming is considered in [START_REF] Yang | Exact safety verification of hybrid systems based on bilinear SoS representation[END_REF].

An alternative way to find such an invariant is by considering tools such as Lyapunov functions to prove the stability properties of dynamical systems [START_REF] Genesio | On the estimation of asymptotic stability regions: State of the art and new proposals[END_REF]. For example, [START_REF] Parrilo | Semidefinite programming relaxations for semialgebraic problems[END_REF] considers parametric functions to find a Lyapunov function for a system with polynomial dynamics formed by SoS polynomials and employs semidefinite programming (SdP) for the parameter synthesis. In [START_REF] Ratschan | Providing a basin of attraction to a target region by computation of Lyapunov-like functions[END_REF] Lyapunov functions are designed via a branch and relax approach and linear programming to solve the induced constraints. In [START_REF] Goubault | Finding non-polynomial positive invariants and lyapunov functions for polynomial systems through darboux polynomials[END_REF] Darboux polynomials are used to design specific forms of Lyapunov functions involving rational functions, logarithmic and exponential terms. Similar invariants have been also considered in [START_REF] Rebiha | Transcendental inductive invariants generation for non-linear differential and hybrid systems[END_REF].

Safety property verification may also be directly introduced in the design phase, instead of being verified a posteriori, as done in the previous approaches. In [START_REF] Platzer | Differential dynamic logic for verifying parametric hybrid systems[END_REF][START_REF] Platzer | Logic analysis of hybrid systems[END_REF], theorem-proving approaches are employed using symbolic-numeric techniques to synthesize invariants for differential (continuous and hybrid) systems. In particular, quantifier elimination techniques are intensively used and more recently a combination with the approach presented in [START_REF] Kapinski | Simulation-guided lyapunov analysis for hybrid dynamical systems[END_REF] has been defined in [START_REF] Aréchiga | Forward invariant cuts to simplify proofs of safety[END_REF]. Alternatively, techniques searching for barrier certificates aim at determining some parametric function, called barrier, defining an hyper-surface in the state-space which is never crossed by the dynamics of the system, see [START_REF] Dai | Barrier certificates revisited[END_REF][START_REF] Prajna | Optimization-Based Methods for Nonlinear and Hybrid Systems Verification[END_REF][START_REF] Prajna | Safety verification of hybrid systems using barrier certificates[END_REF][START_REF] Prajna | Primal-dual tests for safety and reachability[END_REF][START_REF] Sloth | Compositional safety analysis using barrier certificates[END_REF]. A parameter vector of this barrier has to be found such that the barrier separates the region in which the initial state belongs from the unsafe region. In [START_REF] Prajna | Optimization-Based Methods for Nonlinear and Hybrid Systems Verification[END_REF][START_REF] Prajna | Safety verification of hybrid systems using barrier certificates[END_REF][START_REF] Prajna | Primal-dual tests for safety and reachability[END_REF], polynomial dynamics and barrier function are considered and parameters are designed with SdP, which requires some relaxation to obtain a convex design problem. In [START_REF] Dai | Barrier certificates revisited[END_REF], two candidate functions are combined to define more sophisticated barriers, which parameters are again found via SdP. In [START_REF] Sloth | Compositional safety analysis using barrier certificates[END_REF], linear matrix inequalities and SoS are used to generate the barrier functions for hybrid dynamical systems with polynomial dynamics.

Our work follows this approach for non-linear and possibly non-polynomial continuous-time dynamical systems with bounded perturbations and uses interval analysis for the barrier parameter vector search phase.

Formulation

This section recalls the safety characterization introduced in [START_REF] Prajna | Safety verification of hybrid systems using barrier certificates[END_REF] for continuous-time systems using barrier functions.

Safety for continuous-time systems

Consider the autonomous continuous-time perturbed dynamical system

ẋ = f (x, d), (1)
where x ∈ X ⊆ R n is the state vector and d ∈ D is a constant and bounded disturbance. The set of possible initial states at t = 0 is denoted X 0 ⊂ X . There is some unsafe subset X u ⊆ X that shall not be reached by the system, given any x 0 ∈ X 0 at time t = 0 and any d ∈ D. We assume that classical hypotheses (see, e.g., [START_REF] Bellman | Differential-difference equations[END_REF]) on f are satisfied so that (1) has a unique solution x(t, x 0 , d) ∈ X for any given initial value x 0 ∈ X 0 at time t = 0 and any d ∈ D.

Definition 1 The dynamical system (1) is safe if ∀x 0 ∈ X 0 , ∀d ∈ D and ∀t 0, x(t, x 0 , d) / ∈ X u .

Barrier certificates

A way to prove that (1) is safe is by the barrier certificate approach introduced in [START_REF] Prajna | Safety verification of hybrid systems using barrier certificates[END_REF]. A barrier is a differentiable function B : X → R that partitions the state space X into X -where B(x) 0 and X + where B(x) > 0 such that X 0 ⊆ X -and X u ⊆ X + . Moreover, B has to be such that ∀x 0 ∈ X 0 , ∀d ∈ D, ∀t 0, B(x(t, x 0 , d)) 0.

Proving that B(x(t, x 0 , d)) 0 requires an evaluation of the solution of (1) for all x 0 ∈ X 0 and d ∈ D. Alternatively, [START_REF] Prajna | Safety verification of hybrid systems using barrier certificates[END_REF] provides some sufficient conditions a barrier function has to satisfy to prove the safety of a dynamical system, see Theorem 1.

Theorem 1 Consider the dynamical system (1) and the sets X , D, X 0 and X u . If there exists a function B :

X → R such that ∀x ∈X 0 , B(x) 0, ( 2
)
∀x ∈X u , B(x) > 0, (3)
∀x ∈ X , ∀d ∈ D,

B(x) = 0 =⇒ ∂B(x) ∂x , f (x, d) < 0, (4)
then (1) is safe.

In (4) ., . stands for the dot product in R n . In Theorem 1, (2) and (3) ensure that X 0 ⊆ X -, and X u ⊆ X + , while (4) states that if x is on the border between X -and X + (i.e., B(x) = 0), then the dynamics f pushes the state back in X -whatever the value of the disturbance d.

Parametric barrier functions

The search for a barrier B is challenging since it is over a functional space. As in [START_REF] Prajna | Safety verification of hybrid systems using barrier certificates[END_REF], this paper considers barriers belonging to a family of parametric functions (or templates) B(x, p) depending on a parameter vector p ∈ P ⊆ R m . Then one may search for some parameter value p such that B(x, p) satisfies ( 2)-(4).

If there is no p ∈ P such that B(x, p) satisfies ( 2)-(4), this does not mean that the system is not safe: other structures of functions B(x, p) could provide a barrier certificate.

Characterization using interval analysis

This section presents an approach to find a barrier function that fulfills the constraints of Theorem 1. These constraints are first reformulated to cast the design of a barrier function as a quantified constraint satisfaction problem (QCSP) [START_REF] Ratschan | Efficient solving of quantified inequality constraints over the real numbers[END_REF].

Constraint satisfaction problem

Assume that there exist some functions g 0 : X → R and g u : X → R, such that

X 0 = {x ∈ X | g 0 (x) 0} (5)
and

X u = {x ∈ X | g u (x) 0}. ( 6
)
Theorem 1 may be reformulated as follows.

Proposition 2 If ∃p ∈ P such that ∀x ∈ X , ∀d ∈ D ξ (x, p, d) = (g 0 (x) > 0 ∨ B(x, p) 0) (7) ∧ (g u (x) > 0 ∨ B(x, p) > 0) (8) ∧ B(x, p) = 0 ∨ ∂B ∂x (x, p), f (x, d) < 0 (9)
holds true, then the dynamical system (1) is safe.

PROOF. The first component of ξ (x, p, d),

ξ 0 (x, p) = (g 0 (x) > 0 ∨ B(x, p) 0) (10)
may be rewritten as ξ 0 (x, p) = (g 0 (x) 0 =⇒ B(x, p) 0) , see, e.g., [START_REF] Huth | Logic in Computer Science: Modelling and Reasoning about Systems[END_REF]. If ξ 0 (x, p) holds true for some p ∈ P and x ∈ X , then one has either x ∈ X 0 and B(x, p) 0, or x / ∈ X 0 . In both cases, (2) is satisfied. Similarly, the second component of ξ (x, p, d) representing (3) may be rewritten as

ξ u (x, p) = (g u (x) 0 =⇒ B(x, p) > 0) . (11)
If ξ u (x, p) holds true for some p ∈ P and x ∈ X , then one has either x ∈ X u and B(x, p) 0, or x / ∈ X u . In both cases, (3) is satisfied. Now, one may rewrite the last component of t (x, p, d),

ξ b (x, p, d) = B(x, p) = 0 ∨ ∂B ∂x (x, p), f (x, d) < 0 (12)
as

ξ b (x, p, d) = B(x, p) = 0 =⇒ ∂B ∂x (x, p), f (x, d) < 0 , (13)
which corresponds to (4). If ∃p ∈ P such that ∀x ∈ X , ∀d ∈ D, ξ (x, p, d) holds true, then the conditions of Theorem 1 are satisfied and (1) is safe. 2

In [START_REF] Prajna | Safety verification of hybrid systems using barrier certificates[END_REF], ( 9) is relaxed into

ξ b (x, p, d) = ∂B ∂x (x, p), f (x, d) < 0 , (14)
with the consequence of possible elimination of barrier functions that would satisfy (9) for some p but not [START_REF] Frehse | SpaceEx: Scalable verification of hybrid systems[END_REF]. Our aim in this paper is to design barrier functions without resorting to this relaxation by considering methods from interval analysis [START_REF] Jaulin | Applied Interval Analysis[END_REF] which allow to consider strongly nonlinear dynamics and barrier functions.

Solving the constraints

To find a valid barrier function one needs to find some p ∈ P such that B(x, p) satisfies the conditions of Proposition 2. For that purpose, the Computable Sufficient Conditions-Feasible Point Searcher (CSC-FPS) algorithm [START_REF] Jaulin | Guaranteed tuning, with application to robust control and motion planning[END_REF] is adapted.

In what follows, we assume that X , D, and P are boxes, i.e., X = [x], D = [d], and P = [p]. CSC-FPS may also be applied when X , D, and P consist of a union of non-overlapping boxes.

Consider some function g :

R n × R m → R k and some box [z] ∈ IR k . CSC-FPS is designed to determine whether ∃p ∈ [p] , ∀x ∈ [x] , g (x, p) ∈ [z] (15)
and to provide some satisfying p. We extend CSC-FPS to handle conjunctions and disjunctions of constraints and supplement it with efficient pruning techniques involving contractors provided by interval analysis [START_REF] Jaulin | Applied Interval Analysis[END_REF].

FPS branches over the parameter search box [p]

. Branching is performed based on the results provided by CSC. For a given box [p] 0 ⊆ [p], CSC returns true when it manages to prove that ( 15) is satisfied for some p ∈ [p] 0 . CSC returns false when it is able to show that there is no p ∈ [p] 0 satisfying [START_REF] Gao | dReal: An SMT solver for nonlinear theories over the reals[END_REF]. CSC returns unknown in the other cases.

In Proposition 2, ξ (x, p, d) consists of the conjunction of three terms of the form

τ (x, p, d) = (u(x, p) ∈ A) ∨ (v(x, p, d) ∈ B) . (16)
For ξ 0 (x, p), defined in (

To illustrate the main ideas of CSC-FPS combined with contractors, one focuses on the generic QCSP

∃p ∈ [p] , ∀x ∈ [x] , ∀d ∈ [d] , τ (x, p, d) holds true. ( 20
)
Finding a solution for such QCSP involves three steps: validation, reduction of the parameter and state spaces, and bisection.

Validation

In the validation step, one tries to prove that some vector p

∈ [p] is such that ∀x ∈ [x], ∀d ∈ [d], τ (x, p, d) holds true. By definition of τ (x, p, d), one has to prove that ∃p ∈ [p] , ∀x ∈ [x] , ∀d ∈ [d] , (u(x, p) ∈ A) ∨ (v(x, p, d) ∈ B) . ( 21
)
For that purpose, one chooses some arbitrary p ∈ [p] and evaluates the set of values u (

[x] , p) = {u (x, p) | x ∈ [x]} and v ([x] , p, [d]) = {v (x, p, d) | x ∈ [x], d ∈ [d]}
IR n → IR k for a function f : R n → R k satisfies for all [x] ∈ IR n , f ([x]) = {f (x) | x ∈ [x]} ⊆ [f ] ([x]) . (22)
The natural inclusion function is the simplest to obtain: all occurrences of the real variables are replaced by their interval counterpart and all arithmetic operations are evaluated using interval arithmetic. More sophisticated inclusion functions such as the centered form, or the Taylor inclusion function may also be used, see [START_REF] Jaulin | Applied Interval Analysis[END_REF].

Using inclusion functions, one may evaluate whether

[u]([x], p) ⊆ A or [v]([x], p, [d]) ⊆ B holds true
for the various A and B defined in ( 17), [START_REF] Gulwani | Constraint-based approach for analysis of hybrid systems[END_REF], and [START_REF] Huth | Logic in Computer Science: Modelling and Reasoning about Systems[END_REF].

Different choices can be considered for p: one can take a random point in [p], the middle, or one of the edges of [p].

Here, only the middle of [p] is considered.

Reduction of the parameter and state spaces

To facilitate the search for p ∈ [p] one may previously eliminate parts of [p] which may be proved not to contain any p satisfying [START_REF] Jaulin | Guaranteed tuning, with application to robust control and motion planning[END_REF]. The elimination process can be done by evaluation or by using contractors [START_REF] Chabert | Contractor programming[END_REF].

Evaluation

Considering the negation of ( 21), one deduces that a box [p] can be eliminated, i.e., shown not to contain any p satisfying ( 16), if

∀p ∈ [p], ∃x ∈ [x] , ∃d ∈ [d] , u(x, p) / ∈ A ∧ v(x, p, d) / ∈ B, (23)
If, using again inclusion functions, one is able to prove that

∃x ∈ [x] , ∃d ∈ [d] , [u](x, [p]) ∩ A = ∅ ∧ [v](x, [p] , d) ∩ B = ∅ (24)
then ( 23) holds true and one can thus eliminate [p].

In general, [u](x, [p]) and [v](x, [p] , d) are not degenerated intervals, i.e., are not reduced to a real value. When A and B are half-lines, as is the case for ( 7) and ( 8), one may be able to show that (24) holds true. Nevertheless, this will be impossible to show for [START_REF] Chabert | Contractor programming[END_REF] since in this case

A = R \ {0} and [u](x, [p]) is not a degenerated interval. To show that [u](x, [p]) ∩ A = ∅, one needs to have [u](x, [p]) = [0, 0], which is impossible in general.

Contractors

Consider some function g : R n → R k and some set Z ⊂ R k .

Definition 4 A contractor C c : IR n → IR n associated to the generic constraint c ≡ g(x) ∈ Z ( 25
)
is a function taking a box [x] as input and returning a box

C c ([x]) satisfying C c ([x]) ⊆ [x] ( 26
)
and g ([x]) ∩ Z = g (C c ([x])) ∩ Z. (27)
C c provides a box containing the set {x ∈ [x] | g(x) ∈ Z} of solutions of g(x) ∈ Z included in [x]: [START_REF] Platzer | Differential dynamic logic for verifying parametric hybrid systems[END_REF] ensures that the returned box is included in [x] and [START_REF] Platzer | Logic analysis of hybrid systems[END_REF] ensures that no solution of g(x) ∈ Z in [x] is lost.

Example 1 Consider the constraint c ≡ x 1 exp (-x 2 ) ∈ [4, ∞[ (28)
and the initial box 28), one deduces that

[x] = ([x 1 ] , [x 2 ]) T = ([0.1, 10] , [0.1, 10]) T . From (
x 1 ∈ ([4, ∞[ / exp (-[x 2 ])) (29)
Moreover, since one also has

x 1 ∈ [x 1 ], one deduces that x 1 ∈ C 1 ([x]), with C 1 ([x]) = ([4, ∞[ / exp (-[x 2 ])) ∩ [x 1 ] (30)
Similarly one has

x 2 ∈ C 2 ([x]), with C 2 ([x]) = -log ([4, ∞[ / [x 1 ]) ∩ [x 2 ] (31)
A contractor C c associated to c is then

C c ([x]) = (C 1 ([x]) , C 2 ([x])) T . ( 32
)
This contractor for

[x] = ([0.1, 10] , [0.1, 10]) T provides C 1 ([x]) = ([0.1165, 10] , [0.1, 0.9163])
T . This contractor corresponds to the generic forward-backward contractor, also called HC4-revise in [START_REF] Chabert | Contractor programming[END_REF]. It may be simply implemented, e.g., using IBEX as follows. Various contractors have been proposed in the literature, e.g., the forward-backward contractor, the contractor by parallel linearization, the Newton contractor, the Krawczyk contractor, etc. [START_REF] Jaulin | Applied Interval Analysis[END_REF]. Most of them are available, e.g., in IBEX. Example 1 shows that a contractor may be easily implemented from the expression of a generic constraint.

Consider now two functions g 1 : R n → R k1 and g 2 : R n → R k2 , two sets Z 1 ⊂ R k1 and Z 2 ⊂ R k2 , and the associated constraints c 1 ≡ g 1 (x) ∈ Z 1 and c 2 ≡ g 2 (x) ∈ Z 2 . Assume that two contractors C c1 and C c2 are available for c 1 and c 2 . A contractor C c1∧c2 for the conjunction c 1 ∧ c 2 of c 1 and c 2 may be obtained as

C c1∧c2 ([x]) = C c1 ([x]) ∩ C c2 ([x]), (33)
or by composition of contractors

C c1∧c2 ([x]) = C c2 (C c1 ([x])). ( 34
) A contractor C c1∨c2 for the disjunction c 1 ∨ c 2 of c 1 and c 2 may be obtained as follows

C c1∨c2 ([x]) = 2{C c1 ([x]) ∪ C c2 ([x])}, (35)
see [START_REF] Chabert | Contractor programming[END_REF], with 2{•} the interval hull of a set.

Using a contractor C c for [START_REF] Parrilo | Semidefinite programming relaxations for semialgebraic problems[END_REF], one is able to characterize some

[x] ⊂ [x] such that ∀x ∈ [x], g(x) / ∈ Z.
Proposition 5 Consider a box [x], the elementary constraint [START_REF] Parrilo | Semidefinite programming relaxations for semialgebraic problems[END_REF], and the contracted box

C c ([x]) ⊆ [x]. Then, ∀x ∈ [x] \C c ([x]) , one has g(x) / ∈ Z, (36)
where

[x] \C c ([x]) denotes the box [x] deprived of C c ([x]
), which is not necessarily a box.

PROOF. Consider x ∈ [x] \C c ([x]
) and assume that g(x) ∈ Z. Since g(x) ∈ Z and x ∈ [x], one should have x ∈ C c ([x]), according to [START_REF] Platzer | Logic analysis of hybrid systems[END_REF], which contradicts the fact that

x ∈ [x] \C c ([x]). 2
Proposition 5 can be used to eliminate [p] or a part of [p] for which it is not possible to find any p satisfying [START_REF] Jaulin | Guaranteed tuning, with application to robust control and motion planning[END_REF].

Consider the constraint τ ≡ (u(x, p) ∈ A) ∨ (v(x, p, d) ∈ B) (37)
and a contractor C τ for this constraint. It involves elementary contractors for the components of the disjunction in [START_REF] Sloth | Compositional safety analysis using barrier certificates[END_REF], combined as in [START_REF] Rebiha | Transcendental inductive invariants generation for non-linear differential and hybrid systems[END_REF]. For the boxes [x], [p], and [d], one gets

[x] ′ , [p] ′ , [d] ′ = C τ [x] , [p] , [d] , (38)
where [x] ′ , [p] ′ , and [d] ′ are the contracted boxes. Three cases may then be considered.

(

) If [p] \ [p] ′ = ∅, then ∀p ∈ [p] \[p] ′ , ∀x ∈ [x], ∀d ∈ [d], u(x, p) / ∈ A ∧ v(x, p, d) / ∈ B, 1
and there is no p

∈ [p] \ [p]
′ such that (37) holds true for all x ∈ [x] and for all d ∈ [d]. Consequently, the search space for p can be reduced to [p] ′ , see Figure 1

(a). (2) If [x] \ [x] ′ = ∅ then, from Proposition 5, one has ∀p ∈[p], ∀x ∈[x] \[x] ′ , ∀d ∈[d], u(x, p) / ∈ A ∧ v(x, p, d) / ∈ B, (40)
and there is no p ∈ [p] such that (37) holds true for all x ∈ [x], see Figure 1

(b). (3) If [d] \ [d] ′ = ∅, then ∀p ∈ [p], ∀x ∈ [x], ∀d ∈ [d] \ [d] ′ , u(x, p) / ∈ A ∧ v(x, p, d) / ∈ B, (41)
and there is no p ∈ [p] such that [START_REF] Sloth | Compositional safety analysis using barrier certificates[END_REF] holds true for all d ∈ [d], see Figure 1 (b).

One can reduce the size of the sets for the state x and the disturbance d on which [START_REF] Sloth | Compositional safety analysis using barrier certificates[END_REF] has to be verified using the contraction on the negation of this constraint. Consider the negation τ of τ

τ ≡ u(x, p) ∈ A ∧ v(x, p, d) ∈ B (42)
and a contractor C τ for this constraint. Assume that after applying C τ for the boxes [x], [p], and [d], one gets

[x] ", [p] ", [d]" = C τ [x] , [p] , [d] . (43)
From Proposition 5, one knows that [p] [p]

∀ (x, p, d) ∈ ([x] × [p] × [d]) \ ([x] " × [p] " × [d] ") , u(x, p) ∈ A ∨ v(x,
[x] d] , [ [p]' [x] d] , [ ' ' [p] [x] d] , [ [p]' [x] d] , [ ' ' ( ) a
[x] d] , [ [p]'' [x] d] , [ '' '' ( ) b p p [p] [x] d] , [ [p]'' [x] d] , [ '' '' ( ) a p p

Bisection

One is unable to decide whether some p ∈ [p] satisfies [START_REF] Jaulin | Guaranteed tuning, with application to robust control and motion planning[END_REF] when

[u]([x], p) ∩ A = ∅ and [u]([x], p) A (45)
or when

[v]([x], p, [d]) ∩ B = ∅ and [v]([x], p, [d]) B. (46)
This situation occurs in two cases. First, when the selected p does not satisfy [START_REF] Jaulin | Guaranteed tuning, with application to robust control and motion planning[END_REF] for all x ∈ [x] and for all d ∈ [d].

Second, when inclusion functions introduce some pessimism, i.e., they provide an over-approximation of the range of functions over intervals.

To address both cases, one may perform bisections of [x] × [d] and try to verify [START_REF] Jaulin | Guaranteed tuning, with application to robust control and motion planning[END_REF] on the resulting sub-boxes for the same p. Bisection allows to isolate subsets of [x] × [d] on which one may show that (23) holds true. Bisections also reduce pessimism, and may thus facilitate the verification of [START_REF] Jaulin | Guaranteed tuning, with application to robust control and motion planning[END_REF]. The bisection of [x] × [d] is performed within CSC as long as the width of the bisected boxes are larger than some ε x > 0. When CSC is unable to prove [START_REF] Jaulin | Guaranteed tuning, with application to robust control and motion planning[END_REF] or [START_REF] Papachristodoulou | On the construction of Lyapunov functions using the sum of squares decomposition[END_REF] and when all bisected boxes are smaller than ε x , CSC returns unknown.

FPS performs similar bisections on [p] and stops when the width of all bisected boxes are smaller than ε p > 0.

Composition of constraints

To prove the safety of the dynamical system (1), Proposition 2 shows that one has to find some p ∈ [p] such that ∀x ∈ [x], ∀d ∈ [d], ξ (x, p, d) holds true. Since ξ (x, p, d) consists of the conjunction of three elementary constraints of the form [START_REF] Jaulin | Applied Interval Analysis[END_REF], validation requires the verification of ( 21) for the same p considering ( 7), [START_REF] Benhamou | Universally quantified interval constraints[END_REF], and ( 9) simultaneously.

Invalidation may be performed as soon as one is able to prove that one of the constraints ( 7), ( 8), or (9) does not hold true using [START_REF] Papachristodoulou | On the construction of Lyapunov functions using the sum of squares decomposition[END_REF]. Contraction may benefit from the conjunction or disjunction of these constraints, as introduced in ( 34) and ( 35).

CSC-FPS algorithms with contractors

The CSC-FPS algorithm, presented in [START_REF] Jaulin | Guaranteed tuning, with application to robust control and motion planning[END_REF] is supplemented with the contractors introduced in Section 4.2 to improve its efficiency. No specific contractor is indicated in this section. Nevertheless, the forward-backward contractor used in the experimental part provides good results and is provided directly from the explicit expression of the constraint using tools such as IBEX. ′ 0 is too small, even if one was not able to decide whether it contains a satisfying p, it is not further considered to ensure termination of FPS in finite time. The price to be paid in such situation is the impossibility to state whether the initial box [p] contains some satisfying p. This is done by setting flag to unknown at Line 31. Finally, when Q = ∅, no satisfying p has been found. Whether For that purpose, due to the pessimism of inclusion functions, it may be necessary to bisect

FPS, described in

[x] × [d] in subboxes stored in a stack S initialized with [x] × [d]. Algorithm 1 FPS 1: procedure FPS(ξ0, ξu, ξ b , [p],[x], [d]) 2:
⊲ ξ0, ξu, and ξ b associated to ( 7), [START_REF] Benhamou | Universally quantified interval constraints[END_REF], and ( 9 ⊲ Call CSC for each constraint ( 7), [START_REF] Benhamou | Universally quantified interval constraints[END_REF], and (9) 13:

r0:= CSC(ξ0,[p] ′ 0 , [x],[d]) 14: ru:= CSC(ξu,[p] ′ 0 , [x],[d]) 15: r b := CSC(ξ b ,[p] ′ 0 , [x],[d]) 16: ⊲ m([p] ′ 0 ) is
For each subbox [x] 0 × [d] 0 ⊆ [x] × [d]
, CSC determines at Line 7 whether m([p] 0 ) is satisfying. Alternatively, CSC tries to prove that [p] 0 does not contain any satisfying p for all (x, d) ∈ [x] 0 × [d] 0 . This is done in two steps. First, at Line 11, one applies (24) considering some (x, d) ∈ [x] 0 × [d] 0 , here taken as the midpoint of [x] 0 × [d] 0 . Second, one applies Proposition 5 at Lines 15-17 using the result of a contractor for τ , as described in [START_REF] Vaněček | Control systems: from linear analysis to synthesis of chaos[END_REF] and [START_REF] Yang | Exact safety verification of hybrid systems based on bilinear SoS representation[END_REF].

When one is not able to conclude and provided that w(

[x] 0 × [d] 0 ) is larger than ε x , some parts of [x] 0 × [d] 0 for which m([p] 0 )
is satisfying are removed at Line 23, e.g., using the forward-backward contractor, before performing a bisection and storing the resulting subboxes in S. When w(

[x] 0 × [d] 0 ) is less than ε x , to ensure completion of CSC in a finite time, [x] 0 × [d] 0 is

Examples

This section presents experiments for the characterization of barrier functions. The considered dynamical systems are described first before providing numerical results, comparison of different approaches and discussions.

Algorithm 2 CSC 1: procedure CSC(τ ,[p] 0 ,[x],[d]) ⊲ τ is of the form (16) 2: stack S := [x] × [d] 3: flag:=true 4: while S = ∅ do 5: [x] 0 × [d] 0 :=pop(S) 6: ⊲ Validation of m([p] 0 ) using (21) 7: if [u]([x] 0 , m([p] 0 )) ⊆ A∨[v]([x] 0 , m([p] 0 ), [d] 0 ) ⊆ B then 8: continue 9: end if 10: ⊲ Elimination of [p] 0 using (24) 11: if [u] m [x] 0 , [p] 0 ∩ A = ∅ ∧ [v] m [x] 0 , [p] 0 , m [d] 0 ∩ B = ∅ then 12:
return(false) 13: end if 14:

⊲ Reduction of [p] 0 using (38) 15: ([x] ′ 0 , [p] ′ 0 , [d] ′ 0 ) := Cτ ([x] 0 , [p] 0 , [d] 0 ) 16:
⊲ Prove with Proposition 5 that [p] 0 does not contain any satisfying p 17

: if [x] ′ 0 = [x] 0 ∨ [d] ′ 0 = [d] 0 then 18: return(false) 19: end if 20: ⊲ [x] 0 × [d] 0 undetermined and large enough to be bisected 21: if (w([x] 0 × [d] 0 ) εx) then 22:
⊲ Reduction of the state space using (43) and (44) 23:

(

[x] ′′ 0 , [d] ′′ 0 ) := C τ ([x] 0 , m([p] 0 ), [d] 0 ) 24:
⊲ Bisection 25:

(

[x] 1 × [d] 1 , [x] 2 × [d] 2 ):=bisect([x] ′′ 0 × [d] ′′ 0 ) 26: push([x] 1 × [d] 1 ) in S 27: push([x] 2 × [d] 2 ) in S 28: else 29: ⊲ [x] 0 × [d]
0 too small to be further explored 30:

flag:=unknown 31: end if 32: end while 33: return(flag) 34: end procedure

Dynamical system descriptions

For the following examples, one provides the dynamics of the system, the constraints g 0 and g u for the definition of the sets X 0 and X u , the state space [x], and the parametric expression of the barrier function. In all cases, the parameter space is chosen as [p] = [-10, 10] m where m is the number of parameters.

Example 2 Consider the system Example 5 Consider the disturbed system from [START_REF] Prajna | Safety verification of hybrid systems using barrier certificates[END_REF] Example 6 Consider the system with a limit cycle

ẋ1 ẋ2 = x 1 + x 2 x 1 x 2 -0.
ẋ1 ẋ2 = x 2 -x 1 + d 3 x 3 1 -x 2 with g 0 (x) = (x 1 -1.5) 2 + x 2 2 -0.
ẋ1 ẋ2 = x 2 + (1 -x 2 1 -x 2 2 )x 1 + ln(x 2 1 + 1) -x 1 + (1 -x 2 1 -x 2 2 )x 2 + ln(x 2 2 + 1)
with g 0 (x) = (x 1 -1) 2 + (x 2 + 1.5) 2 -0.05, g u (x) = (x 1 + 0.

    ẋ1 ẋ2 ẋ3     =     10(x 2 -x 1 )
x 1 (28x 3 )x 2

            ẋ1 ẋ2 ẋ3 ẋ4 ẋ5 ẋ6             =             -x 1 + 4x 2 -6x 3 x 4 -x 1 -x 2 + x 3 5 x 1 x 4 -x 3 + x 4 x 6 x 1 x 3 + x 3 x 6 -x 3 4 -2x 3 2 -x 5 + x 6 -3x 3 x 4 -x 3 5 -x 6             with g 0 (x) = (x 1 -3.05) 2 + (x 2 -3.05) 2 + (x 3 -3.05)

Experimental conditions and results

CSC-FPS, presented in Section 4, has been implemented using the IBEX library [START_REF] Araya | Upper bounding in inner regions for global optimization under inequality constraints[END_REF][START_REF] Chabert | Contractor programming[END_REF]. The selection of candidate barrier functions is performed choosing polynomials with increasing degree, except for Examples 2, 3, and 5, where parametric functions taken from [1,[START_REF] Amir | A globally asymptotically stable polynomial vector field with no polynomial Lyapunov function[END_REF] are considered. For each example, the computing time to get a valid barrier function and the number of bisections of the search box [p] are provided. Table 1 summarizes the results for the versions of CSC-FPS with and without contractors.

As in [START_REF] Jaulin | Guaranteed tuning, with application to robust control and motion planning[END_REF], we choose ε x = 10 -1 and ε p = 10 -5 . Computations were done on an Intel core 1.7 GHz processor with 8 GB of RAM. If after 30 minutes of computations no valid barrier function has been found, the search is stopped. This is denoted by T.O. (time out) in Table 1. Moreover, for Examples 2, 5, and 6, graphical representation of the computed barrier functions are provided. In Figures 3,4, and 5, X 0 is in green, X u is in red, the bold line is the barrier function and some trajectories starting from X 0 are also represented.

The results in Table 1 show the importance of contractors which are beneficial in all cases. Thanks to contractors, valid barrier functions were obtained for all examples, which is not the case employing the original version of CSC-FPS proposed in [START_REF] Jaulin | Guaranteed tuning, with application to robust control and motion planning[END_REF]. In theory, both FPS and CSC are of exponential complexity in the dimension m of the parameter space and n of the state space. In practice, contractors allow, on the considered examples, to get solutions in a reasonable amount of time, even for systems with a larger number of states and parameters.

The search for barrier functions using the relaxed version ( 14) of the constraint [START_REF] Collins | Partial cylindrical algebraic decomposition for quantifier elimination[END_REF] as in [START_REF] Prajna | Safety verification of hybrid systems using barrier certificates[END_REF], has also been performed using the version of our algorithms with contractors considering the same parametric barrier functions.

We were not able to find a valid barrier function in less than 30 minutes. This shows the detrimental effect of the relaxation ( 14) on the search technique.

Some examples were also addressed using RSolver, which is a tool to solve some classes of QCSP [START_REF] Ratschan | Efficient solving of quantified inequality constraints over the real numbers[END_REF]. Nevertheless, this requires some modifications of the examples, since RSolver does not address dynamics or constraints involving divisions [START_REF] Ratschan | Rsolver user manuel[END_REF]. Moreover, the type of constraints processed by RSolver, for problems such as [START_REF] Jaulin | Applied Interval Analysis[END_REF], allows only parametric barrier functions that are linear in the parameters. Outer-approximating boxes for the initial and the unsafe regions X 0 , X u defined by g 0 and g u are given to Rsolver. As a consequence, Examples 2, 3 and 4 could not be considered. Only Examples 5, 6, and 7 were tested by Rsolver, which was able to find a satisfying barrier function only for Lorenz in less than 1s. RSolver was unable to find a solution for Examples 5 and 6. Rsolver is well designed for problems with parametric barriers linear in the parameters, but has difficulties for non-linear problem and non-convex constraints such as [START_REF] Collins | Partial cylindrical algebraic decomposition for quantifier elimination[END_REF].

Table 2 provides results obtained using dReal [START_REF] Gao | dReal: An SMT solver for nonlinear theories over the reals[END_REF], which is sat-modulo-theory solver [START_REF] Barrett | Satisfiability modulo theories[END_REF], based on interval analysis tools to handle first-order logic formulas over the real numbers. In the tested version (3.15.11.03) of dReal, the support for QCSP such as [START_REF] Jaulin | Applied Interval Analysis[END_REF] does not benefit from contractors. A randomization heuristic is used in the search of satisfying parameter vectors. This leads to results obtained in a amount of time varying with the runs. A precision parameter of 0.001 (a tentative to have a fair comparison with respect to values of ε X and ε P used in our algorithm) was employed. Due to the non-deterministic behavior of dReal, each example has be run 10 times and we report minimal and maximal values for the execution time and the number of bisections. When no solution is found in less than 30 minutes, the value T.O. (time out) is indicated.

Table 2 shows that dReal can handle all considered problems despite the absence of contractors. In particular, one observes that randomization speed up the search for valid parameters, as for example, in case of examples 2, 5, 6 and 7. Nevertheless, as a consequence, the overall algorithm has a varying execution time. The absence of contractors may increase the parameter search time, as for example, in the case of examples 1, 3, and 4. Combining randomized methods and contractor programming should be very beneficial to efficiently solve QCSP.

Conclusion

This paper presents a new method to find parametric barrier functions for nonlinear continuous-time perturbed dynamical systems. The proposed technique addresses the design of quite general barrier functions for continuoustime systems described by nonlinear dynamics that are also general. The search for barrier functions is formulated as an interval quantified constraint satisfaction problem. A branch-and-prune algorithm proposed in [START_REF] Jaulin | Guaranteed tuning, with application to robust control and motion planning[END_REF] has been supplemented with contractors to address this problem. Contractors are instrumental in solving problems with large number of parameters. The proposed approach can thus find barrier functions for a large class of possibly perturbed dynamical systems.

Alternative techniques based on RSolver or dReal may be significantly more efficient for some specific classes of problems, e.g., where the parameters appear linearly in the parametric barrier functions. Combining ideas from RSolver, dReal and our approach may be very useful to further improve the global efficiency of barrier function characterization.

Future work will be dedicated to the search for the class of parametric barrier functions to consider. This may be done by exploring a library of candidate barrier functions. In our approach rejection of a candidate function occurs mainly after a timeout. Even if contractors aiming at eliminating some parts of the parameter space were defined, their efficiency is limited. Better contractors for that purpose may be very helpful.

An other future research direction is to extend the proposed method to hybrid systems as done in [START_REF] Prajna | Safety verification of hybrid systems using barrier certificates[END_REF], i.e., to consider a set of quantified constraints for each location of an hybrid automaton and the constraints associated to the transitions.

  , A = ]0, +∞[ and B = ]-∞, 0] ; (17) for ξ u (x, p), defined in (8), A = ]0, +∞[ and B = ]0, +∞[ ; (18) for ξ b (x, p, d), defined in (9), A = ]-∞, 0[ ∪ ]0, +∞[ and B = ]-∞, 0[ .

Definition 3

 3 taken by u(x, p) and v(x, p, d) for all x ∈ [x] and for all d ∈ [d]. Outer-approximations of u ([x] , p) and v ([x] , p, [d]) are easily obtained using inclusion functions provided by interval analysis. An inclusion function [f ] :

  Variable x1,x2; NumConstraint c(x1, x2, x1*exp(-x2) >= 4.0); CtcFwdBwd contractor(c); IntervalVector x(2); x[0]=Interval(0.1,10); x[1]=Interval(0.1,10); contractor.contract(x);

  p, d) ∈ B. (44) Indeed, if [p] " = [p], one can focus on the search for some p ∈ [p] satisfying (37) by considering only [x]" × [d] ", since for all (x, d) ∈ ([x] × [d]) \ ([x] " × [d] "), (37) is satisfied for all p ∈ [p], see Figure 2 (a). Now, if [p] " = [p], then any p ∈ [p] \ [p] " will satisfy (37) for all (x, d) ∈ ([x] × [d]), see Figure 2 (b).

Fig. 1 .

 1 Fig. 1. Contractions using Cτ ; the set for which (37) is satisfied is in gray; (a) [p] \ [p] ′ = ∅ and the search space for satisfying p can be reduced to [p] ′ ; (b) [x] ′ = [x] or [d] ′ = [d], it is thus not possible to find some p ∈ [p] such that (37) is satisfied for all x ∈ [x] and all d ∈ [d].

Fig. 2 .

 2 Fig. 2. Contractions using C τ ; the set for which (42) is satisfied is in white; (a) [x]" = [x] and/or [d]" = [d] and one has only to find some suitable p ∈ [p] such that (37) is satisfied for all x ∈ [x] " and all d ∈ [d] "; (b) [p]" = [p], one may choose any p ∈ [p] \ [p] " (for example the value of p indicated in red) and (37) will hold true for all (x, d) ∈ ([x] × [d]).

′ 0 ,

 0 the midpoint of [p] ′ 0 , is satisfying. When all calls of CSC return true at Line 17, a barrier function with parameter m [p] ′ 0 is found. When one of the calls of CSC returns false at Line 21, [p] ′ 0 is proved not to contain any satisfying p. In all other cases, when w([p] ′ 0 ), the width of [p] ′ 0 , is larger than ε p , [p] ′ 0 is bisected and the resulting subboxes are stored in Q, see Lines 24-28. When [p]

  [p] may however contain some satisfying p in an unexplored subbox [p] ′ 0 depends on the value of flag. CSC, described in Algorithm 2, determines whether m([p]) satisfies (16) by showing that τ (x, m([p]), d) holds true for all x ∈ [x] and all d ∈ [d]. Alternatively, CSC may prove that there is no p ∈ [p] satisfying (16) for all x ∈ [x] and all d ∈ [d].

) 3 :

 3 queue Q := [p] 4: flag := true 5: while Q = ∅ do 6: [p] 0 := dequeue(Q) 7: ⊲ Reduction of [p] 0 using (34), (38), and (39) 8: [p] ′ 0 := C ξ 0 ∧ξu (m([x]), [p] 0 ) 9: if [p] ′ 0 = ∅ then 10: continue ⊲ There is no satisfying p ∈ [p] 0 11: end if 12:

  not further explored. The price to be paid is an impossibility to determine whether m([p] 0 ) is satisfying for all (x, d) ∈ [x] 0 × [d] 0 . This is indicated by setting flag to unknown at Line 30. Nevertheless, one may still prove that [p] 0 does not contain any satisfying p considering other subboxes of [x] × [d].

5x 2 2 with g 0 (x0+p3) 2 +p 2 2 + x 1 + p 4 . 3 1 + x 1 x 2 -x 2 with g 0 2 + p 3 . 4 x 2 2 - 1 + p 2 x 2 2 + p 3 x 1 x 2 + p 4 x 1 + p 5 x 2 +

 2021431220234212212 (x) = (x 1 +1.25) 2 +(x 2 -1.25) 2 -0.05, g u (x) = (x 1 +2.5) 2 +(x 2 -0.8) 2 -0.05, and [x] = [-10 3 , 0]×[-10 3 , 10 3 ]. The parametric barrier function is B(x, p) = p1p2(x0+p3) Example Consider the system from [2] (x) = (x 1 -1.125) 2 + (x 2 -0.625) 2 -0.0125, g u (x) = (x 1 -0.875) 2 + (x 2 -0.125) 2 -0.0125, and [x] = [-100, 100] × [-100, 100]. The parametric barrier function used is B(x, p) = ln(p 1 x 1 )ln(x 2 ) + p 2 x Example Consider the system from [24] 0.5, g u (x) = (x 1 -3.5) 2 + (x 2 -0.5) 2 -0.5 , and [x] = [-10 3 , 10 3 ] × [-100, 100]. A quadratic parametric barrier function is chosen B(x, p) = p 1 x 2 p 6 .

2 + x2+p3 p4 2 - 1 . 7

 2217 6) + (x 2 -1) 2 -0.05, and [x] = [-10 3 , 10 3 ] × [-10 3 , 10 3 ]. The parametric barrier function used is B(x, p) = x1+p1 p2Example Consider the Lorenz system from[START_REF] Vaněček | Control systems: from linear analysis to synthesis of chaos[END_REF]

x 1 x 2 1 + p 2 x 1 + p 3 x 3 + p 4 . 8

 211348 g 0 (x) = (x 1 +14.5) 2 +(x 2 +14.5) 2 +(x 3 -12.5) 2 -0.25, g u (x) = (x 1 +16.5) 2 +(x 2 +14.5) 2 +(x 3 -2.5) 2 -0.25, and [x] = [-20, 20]× [-20, 0]× [-20, 20]. The considered parametric barrier function is B(x, p) = p 1 x 2Example Consider the system from[START_REF] Papachristodoulou | On the construction of Lyapunov functions using the sum of squares decomposition[END_REF]

Fig. 3 .

 3 Fig. 3. Results for Example 2.

Fig. 4 .

 4 Fig. 4. Results for Example 5 with various values of the disturbance d.

2 Fig. 5 .

 25 Fig. 5. Results for Example 6.

  25, g u (x) = (x 1 + 0.8) 2 + (x 2 + 1) 2 -0.25, [x] = [-100, 100] × [-10, 10], and d ∈ [0.9, 1.1]. The parametric barrier function B(x, p) = p 1 x 2 1 + p 2 x 2 2 + p 3 x 1 x 2 + p 4 x 1 + p 5 x 2 + p 6 is considered.

  2 + (x 4 -3.05) 2 + (x 5 -3.05) 2 + (x 6 -3.05) 2 -0.0001, g u (x) = (x 1 -7.05) 2 + (x 2 -3.05) 2 + (x 3 -7.05) 2 + (x 4 -7.05) 2 + (x 5 -7.05) 2 + (x 6 -7.05) 2 -0.0001, and [x] = [0, 10] × [0, 10] × [2, 10] × [0, 10] × [0, 10] × [0, 10]. The considered parametric barrier function is B(x, p) =

	p 1 x 2 1 + p 2 x 4 2 + p 3 x 2 3 + p 4 x 2 4 + p 5 x 4 5 + p 6 x 2 6 + p 7 .

Table 1

 1

	Results of CSC-FPS without and with contractors
				Without contr.	With contr.
	Example n m	time	bisect. time bisect.
	1	2	4	36s	4520	16s	4553
	2	2	3	T.O.	/	1s	159
	3	2	6	1133s	20388	1s	6
	4	2	6	253s	14733	7s	435
	5	2	4	T.O.	/	98s	4072
	6	3	4	167s	1753	21s	47
	7	6	7	697s	67600	1s	261

Table 2

 2 Results of dReal for 10 executions with precision of 0.001

		Time (in s)	Bisections
	Example	Min	Max	Min	Max
	1	285.1 T.O. 31337	/
	2	0.06	0.19	35	94
	3	113.3 T.O.	61	/
	4	533.3 T.O.	51	/
	5	0.05	0.06	34	39
	6	0.08	0.08	54	55
	7	0.4	0.4	107	109

Acknowledgements

Authors are thankful to Sicun Gao and Soonho Kong for their valuable help on dReal.

⋆ This research was partially supported by Labex DigiCosme (project ANR-11-LABEX-0045-DIGICOSME) operated by ANR as part of the program "Investissement d'Avenir" Idex Paris-Saclay (ANR-11-IDEX-0003-02), by the ANR INS Project CAFEIN (ANR-12-INSE-0007), and by the iCODE Institute, research project of the IDEX Paris-Saclay.