Stefan Mitsch
email: smitsch@cs.cmu.edu

Khalil Ghorbal
email: kghorbal@cs.cmu.edu

André Platzer
email: aplatzer@cs.cmu.edu

On Provably Safe Obstacle Avoidance for Autonomous Robotic Ground Vehicles

Nowadays, robots interact more frequently with a dynamic environment outside limited manufacturing sites and in close proximity with humans. Thus, safety of motion and obstacle avoidance are vital safety features of such robots. We formally study two safety properties of avoiding both stationary and moving obstacles: (i) passive safety, which ensures that no collisions can happen while the robot moves, and (ii) the stronger passive friendly safety in which the robot further maintains sufficient maneuvering distance for obstacles to avoid collision as well. We use hybrid system models and theorem proving techniques that describe and formally verify the robot's discrete control decisions along with its continuous, physical motion. Moreover, we formally prove that safety can still be guaranteed despite location and actuator uncertainty.

I. INTRODUCTION

With the increased introduction of autonomous robotic ground vehicles as consumer products-such as autonomous household appliances [START_REF] Fiorini | Cleaning and Household Robots: A Technology Survey[END_REF] or driverless cars on regular Californian roads 1 -we face an increased need for ensuring safety not only for the sake of the consumer, but also the manufacturer. Since those robots are designed for environments that occupy stationary as well as moving obstacles, motion safety and obstacle avoidance are vital safety features for such robots [START_REF] Bouraine | Provably safe navigation for mobile robots with limited field-of-views in dynamic environments[END_REF][START_REF] Platzer | Formal Verification of Curved Flight Collision Avoidance Maneuvers: A Case Study[END_REF][START_REF] Täubig | Guaranteeing functional safety: design for provability and computer-aided verification[END_REF][START_REF] Wu | Guaranteed infinite horizon avoidance of unpredictable, dynamically constrained obstacles[END_REF].

In this paper, we prove safety for the obstacle avoidance control algorithm of a robot. One of the conceptual difficulties is what safety means for an autonomous robot. We would want it to be collision-free, but that usually requires other vehicles to be sensible, e. g., not actively try to run into the robot when it is just stopped in a corner. One way of doing that is to assume stringent constraints on the behavior of obstacles (see, e. g., [START_REF] Bouraine | Provably safe navigation for mobile robots with limited field-of-views in dynamic environments[END_REF][START_REF] Loos | Adaptive Cruise Control: Hybrid, Distributed, and Now Formally Verified[END_REF]). In this paper, we want to refrain from doing so and allow arbitrary obstacles with an arbitrary continuous motion with a known upper bound on their velocity. Then our robot is safe, intuitively, if no collision can ever happen where the robot is to blame. The first notion we consider is passive safety [START_REF] Maček | Towards Safe Vehicle Navigation in Dynamic Urban Scenarios[END_REF]. Passive safety can be guaranteed with minimal assumptions about obstacles. It requires that the robot does not actively collide, i. e., if a collision occurs at all then only while the robot is stopped and the (moving) obstacle ran into the robot. The difficulty with that notion is that it still allows the robot to go kamikaze and stop in unsafe places, creating unavoidable collision situations in which an obstacle has no control choices left that would prevent a collision. The second notion we consider is passive friendly safety [START_REF] Maček | Towards Safe Vehicle Navigation in Dynamic Urban Scenarios[END_REF], which aims 1 http://www.nytimes.com/2010/10/10/science/10google.html? r=0 for more careful robot decisions that respect the features of moving obstacles (e. g., their braking capabilities). Thus, a passive friendly robot not only ensures that it is itself able to stop before a collision occurs, but it also maintains sufficient maneuvering room for obstacles to avoid a collision as well.

Motion safety and obstacle avoidance lead to interesting cognitive robotics questions: how much does the robot have to know about the goals and constraints of other vehicles so as not to be considered to blame? In this paper, we successively construct models and proofs that increase the level of assumed knowledge and explicitly expressed uncertainty. We start with (i) passive safety, which assumes a known upper bound on the velocity of obstacles. Then we extend to (ii) passive friendly safety for a known lower bound of an obstacle's braking power and consider an upper bound on the reaction time that we allow an obstacle to start collision avoidance attempts. Note that all our models use symbolic bounds so they hold for all choices of the bounds. As a result, we can account for uncertainty in several places (e. g., by instantiating upper bounds on acceleration or time with values including uncertainty). We additionally show how further uncertainty that cannot be attributed to such bounds (in particular location and actuator uncertainty) can be modeled and verified explicitly.

We use the well-known dynamic window algorithm [START_REF] Fox | The dynamic window approach to collision avoidance[END_REF] and verify passive safety and passive friendly safety. Unlike existing work on obstacle avoidance (e. g., [START_REF] Althoff | Safety assessment of robot trajectories for navigation in uncertain and dynamic environments[END_REF][START_REF] Jia Pan | Collision-free and smooth trajectory computation in cluttered environments[END_REF][START_REF] Seward | Safe and effective navigation of autonomous robots in hazardous environments[END_REF][START_REF] Täubig | Guaranteeing functional safety: design for provability and computer-aided verification[END_REF][START_REF] Van Den | LQG-MP: Optimized Path Planning for Robots with Motion Uncertainty and Imperfect State Information[END_REF]), we use hybrid models and verification techniques that describe and verify the robot's discrete control choices along with its continuous, physical motion. In summary, our contributions are (i) hybrid models of navigation and obstacle avoidance control algorithms of robots, and (ii) proofs that they guarantee passive and passive friendly safety in the presence of stationary and moving obstacles despite sensor and actuator uncertainty. The models and proofs of this paper are bundled with KeYmaera. 2The paper is organized as follows. In the next section, we discuss related work on navigation and obstacle avoidance of robots that focuses on verification. Section III recalls differential dynamic logic that we use as a modeling formalism for the hybrid dynamics of a robot and the safety constraints. In Sect. IV, we introduce models of obstacle avoidance control and physical motion, and prove that they guarantee passive safety and passive friendly safety with stationary as well as moving obstacles. We then model uncertainty explicitly and prove that motion is still safe. Sect. V concludes the paper.

II. RELATED WORK

Isabelle has recently been used to formally verify that a C program implements the specification of the dynamic window algorithm [START_REF] Täubig | Guaranteeing functional safety: design for provability and computer-aided verification[END_REF]. This is interesting, but the algorithm itself and its impact on motion of the robot was considered in an informal pen-and-paper argument only. We, instead, formally verify the correctness of the dynamic window algorithm using a hybrid verification technique. Thus, we complement the work in [START_REF] Täubig | Guaranteeing functional safety: design for provability and computer-aided verification[END_REF] in a twofold manner: First, we create hybrid models of the control and the motion dynamics of the robot and formally verify correctness of the dynamic window algorithm control for the combined dynamics. Second, we model stationary as well as moving obstacles and prove two versions of safety. These complementary results together present a strong safety argument from concept (this paper) to implementation [START_REF] Täubig | Guaranteeing functional safety: design for provability and computer-aided verification[END_REF].

PASSAVOID [START_REF] Bouraine | Provably safe navigation for mobile robots with limited field-of-views in dynamic environments[END_REF] is a navigation scheme, which avoids braking inevitable collision states (i. e., states that regardless of the robot's trajectory lead to a collision) to achieve safety in the presence of moving obstacles. Since PASSAVOID is designed to operate in completely unknown environments, it ensures that the robot is at rest when a collision occurs (passive safety). The motion dynamics of the robot have only been considered in simulation. We prove the stronger passive friendly safety using a hybrid verification technique (i. e., algorithm and motion dynamics), which ensures that the robot does not create unavoidable collision situations by stopping in unsafe places.

For the purpose of guaranteeing infinite horizon safety, velocity obstacle sets [START_REF] Wu | Guaranteed infinite horizon avoidance of unpredictable, dynamically constrained obstacles[END_REF] assume unpredictable behavior for obstacles with known forward speed and maximum turn rate (i. e., Dubin's cars). The authors focus only on the obstacle behavior; the robot's motion is explicitly excluded from their work. We complement their work and show that a robot, which has a known upper bound on its reaction time and considers discs as velocity obstacle sets (i. e., known forward speed and unknown turn rate, as allowed by [START_REF] Wu | Guaranteed infinite horizon avoidance of unpredictable, dynamically constrained obstacles[END_REF]), moves safely.

Hybrid models of driver support systems in cars [START_REF] Loos | Adaptive Cruise Control: Hybrid, Distributed, and Now Formally Verified[END_REF][START_REF] Mitsch | Towards formal verification of freeway traffic control[END_REF] have been verified with a model of the continuous dynamics of cars. That points out interesting safety conditions for vehicles on straight lines, but not in the general two-dimensional case that we consider in this work.

Safety of aircraft collision avoidance maneuvers in the two-dimensional plane was verified for constant translational velocity and a rotational velocity that stays constant during the maneuver [START_REF] Loos | Formal verification of distributed aircraft controllers[END_REF][START_REF] Platzer | Formal Verification of Curved Flight Collision Avoidance Maneuvers: A Case Study[END_REF]. Our models include acceleration for both translational and rotational velocity.

LTLMoP contains an approach [START_REF] Sarid | Guaranteeing High-Level Behaviors while Exploring Partially Known Maps[END_REF] to guarantee highlevel behavior for map exploration when the environment is continuously updated. The approach synthesizes and resynthesizes plans, expressed in linear temporal logic, of a hybrid controller, when new map information is discovered. This work focuses on preserving the state and task completion history, and thus on guaranteeing that the robot will follow a high-level behavior (e. g., visit all rooms) even when the controller is re-synthesized, not on safe obstacle avoidance.

Pan et al. [START_REF] Jia Pan | Collision-free and smooth trajectory computation in cluttered environments[END_REF] propose a method to smooth the trajectories produced by sampling-based planners in a collision-free manner. Our paper proves that such trajectories are indeed safe when considering the control choices of a robot and its continuous dynamics. [START_REF] Van Den | LQG-MP: Optimized Path Planning for Robots with Motion Uncertainty and Imperfect State Information[END_REF] is a motion planning approach that takes into account the sensors, controllers, and motion dynamics of a robot while working with uncertain information about the environment. The approach assesses randomly generated paths by the approximated probability of a collision with an obstacle. One goal is to select paths with low collision probability; however, guaranteeing collision-free motion is not their focus, since a collision-free path may not even have been generated.

LQG-MP

Althoff et al. [START_REF] Althoff | Safety assessment of robot trajectories for navigation in uncertain and dynamic environments[END_REF] use a probabilistic approach to rank trajectories according to their collision probability. To further refine such a ranking, a collision cost metric is proposed, which derives the cost of a potential collision by considering the relative speeds and masses of the collision objects.

Seward et al. [START_REF] Seward | Safe and effective navigation of autonomous robots in hazardous environments[END_REF] try to avoid potentially hazardous situations by using Partially Observable Markov Decision Processes. Their focus, however, is on a user-definable tradeoff between safety and progress towards a goal. Hence, safety is not guaranteed under all circumstances.

In summary, this paper addresses safety of robot obstacle avoidance in the following manner.

• Unlike [START_REF] Bouraine | Provably safe navigation for mobile robots with limited field-of-views in dynamic environments[END_REF][START_REF] Sarid | Guaranteeing High-Level Behaviors while Exploring Partially Known Maps[END_REF][START_REF] Seward | Safe and effective navigation of autonomous robots in hazardous environments[END_REF][START_REF] Täubig | Guaranteeing functional safety: design for provability and computer-aided verification[END_REF][START_REF] Wu | Guaranteed infinite horizon avoidance of unpredictable, dynamically constrained obstacles[END_REF], we study combined models of the hybrid dynamics in terms of discrete control and differential equations for continuous physical motion of the robot as well as the obstacles, not discrete control alone or only the behavior of obstacles. • Unlike [START_REF] Loos | Adaptive Cruise Control: Hybrid, Distributed, and Now Formally Verified[END_REF][START_REF] Loos | Formal verification of distributed aircraft controllers[END_REF][START_REF] Mitsch | Towards formal verification of freeway traffic control[END_REF][START_REF] Platzer | Formal Verification of Curved Flight Collision Avoidance Maneuvers: A Case Study[END_REF][START_REF] Täubig | Guaranteeing functional safety: design for provability and computer-aided verification[END_REF], we verify safety in the two-dimensional plane not one-dimensional space and do not assume constant translational and rotational velocity, but include accelerations for both, as needed for ground vehicles.

• Unlike [START_REF] Althoff | Safety assessment of robot trajectories for navigation in uncertain and dynamic environments[END_REF][START_REF] Jia Pan | Collision-free and smooth trajectory computation in cluttered environments[END_REF][START_REF] Sarid | Guaranteeing High-Level Behaviors while Exploring Partially Known Maps[END_REF][START_REF] Seward | Safe and effective navigation of autonomous robots in hazardous environments[END_REF][START_REF] Van Den | LQG-MP: Optimized Path Planning for Robots with Motion Uncertainty and Imperfect State Information[END_REF], we produce formal, deductive proofs in a formal verification tool. Note, that in [START_REF] Täubig | Guaranteeing functional safety: design for provability and computer-aided verification[END_REF] the correctness of the implementation w.r.t. the algorithm's specification was formally verified, but not the motion of the robot. • Unlike [START_REF] Jia Pan | Collision-free and smooth trajectory computation in cluttered environments[END_REF][START_REF] Sarid | Guaranteeing High-Level Behaviors while Exploring Partially Known Maps[END_REF][START_REF] Täubig | Guaranteeing functional safety: design for provability and computer-aided verification[END_REF], we verify safety even in the presence of moving obstacles. • Unlike [START_REF] Bouraine | Provably safe navigation for mobile robots with limited field-of-views in dynamic environments[END_REF][START_REF] Jia Pan | Collision-free and smooth trajectory computation in cluttered environments[END_REF][START_REF] Sarid | Guaranteeing High-Level Behaviors while Exploring Partially Known Maps[END_REF][START_REF] Täubig | Guaranteeing functional safety: design for provability and computer-aided verification[END_REF], we verify passive friendly safety, which is important because passive (non-friendly) safe robots may cause unavoidable collisions by stopping in unsafe places so that obstacles will collide with them. • Unlike [START_REF] Bouraine | Provably safe navigation for mobile robots with limited field-of-views in dynamic environments[END_REF][START_REF] Loos | Adaptive Cruise Control: Hybrid, Distributed, and Now Formally Verified[END_REF][START_REF] Platzer | Formal Verification of Curved Flight Collision Avoidance Maneuvers: A Case Study[END_REF][START_REF] Sarid | Guaranteeing High-Level Behaviors while Exploring Partially Known Maps[END_REF], we consider sensor and actuator uncertainty in our verification results. • Unlike [START_REF] Althoff | Safety assessment of robot trajectories for navigation in uncertain and dynamic environments[END_REF][START_REF] Seward | Safe and effective navigation of autonomous robots in hazardous environments[END_REF][START_REF] Van Den | LQG-MP: Optimized Path Planning for Robots with Motion Uncertainty and Imperfect State Information[END_REF], we do not minimize or probabilistically minimize collisions, but prove that collisions can never occur (as long as the robot system fits to the model).

III. PRELIMINARIES: DIFFERENTIAL DYNAMIC LOGIC

A robot and the moving obstacles in its environment form a hybrid system: they make discrete control choices (e. g., compute the actuator set values for acceleration, braking, or steering), which in turn influence their actual physical behavior (e. g., slow down to a stop, move along a curve). Hybrid systems have been considered as joint models for both components, since verification of either component alone does not capture the full behavior of a robot and its environment. nondeterministic repetition, repeats α n ≥ 0 times x := θ assign value of term θ to variable x (discrete jump) x := * assign arbitrary real number to variable x x 1 = θ1, . . . , evolve xi along differential equation system x i = θi x n = θn & F restricted to maximum evolution domain F

In order to verify safe obstacle avoidance, we use differential dynamic logic dL [START_REF] Platzer | Differential Dynamic Logic for Hybrid Systems[END_REF][START_REF] Platzer | Logical Analysis of Hybrid Systems: Proving Theorems for Complex Dynamics[END_REF][START_REF] Platzer | The Complete Proof Theory of Hybrid Systems[END_REF], which has a notation for hybrid systems as hybrid programs. We use hybrid programs for modeling a robot that follows the dynamic window algorithm as well as for modeling the behavior of moving obstacles. dL allows us to make statements that we want to be true for all runs of the program (safety) or for at least one run (liveness). Both constructs are necessary to verify safety: for all possible control choices and entailed physical motion, our robot must be able to stop, while at the same time there must be at least one possible execution in which the obstacle is able to stop without collision as well. Table I summarizes the syntax of hybrid programs together with an informal semantics. Sequential composition α; β says that β starts after α finishes (e. g., first let the robot choose acceleration, then steering angle). The nondeterministic choice α ∪ β follows either α or β (e. g., let the robot decide nondeterministically between remaining stopped or accelerating). The nondeterministic repetition operator α * repeats α zero or more times (e. g., let the robot repeatedly control and drive). Assignment x := θ instantaneously assigns the value of the term θ to the variable x (e. g., let the robot choose maximum braking), while x := * assigns an arbitrary value to x (e. g., let the robot choose any acceleration). x = θ & F describes a continuous evolution of x within the evolution domain F (e. g., let the velocity of the robot change according to its acceleration, but never lower than zero). The test ?F checks that a particular condition F holds, and aborts if it does not (e. g., test whether or not the distance to an obstacle is large enough). A typical pattern that involves assignment and tests is to limit the assignment of arbitrary values to known bounds (e. g., limit an arbitrarily chosen acceleration to the physical limits of the robot, as in x := * ; ?x ≤ A, which says x is any value less or equal A).

The set of dL formulas is generated by the following EBNF grammar (where ∼ ∈ {<, ≤, =, ≥, >} and θ 1 , θ 2 are arithmetic expressions in +, -, •, / over the reals):

φ ::= θ 1 ∼ θ 2 | ¬φ | φ∧ψ | φ∨ψ | φ → ψ | ∀xφ | [α]φ | α φ
Further operations, such as Euclidian norm θ and infinity norm θ ∞ of a vector θ, are definable.

To specify the desired correctness properties of hybrid programs, dL formulas of the form F → [α]G mean that all executions of the hybrid program α, which start at a state in which formula F is true, lead to states in which formula G is true. Dually, formula F → α G expresses that there is a state reachable by the hybrid program α that satisfies formula G.

Differential dynamic logic comes with a verification technique to prove those correctness properties. We did all our proofs in the verification tool KeYmaera, which implements this verification technique [START_REF] Platzer | Differential Dynamic Logic for Hybrid Systems[END_REF][START_REF] Platzer | Differential-algebraic Dynamic Logic for Differential-algebraic Programs[END_REF][START_REF] Platzer | Logical Analysis of Hybrid Systems: Proving Theorems for Complex Dynamics[END_REF]. KeYmaera supports hybrid systems with nonlinear discrete jumps, nonlinear differential equations, differential-algebraic equations, differential inequalities, and systems with nondeterministic discrete or continuous input. This makes KeYmaera more readily applicable to robotic verification than other hybrid system verification tools, such as SpaceEx [START_REF] Frehse | SpaceEx: Scalable verification of hybrid systems[END_REF], which focuses on piecewise linear systems. KeYmaera implements automatic proof strategies that decompose hybrid systems symbolically [START_REF] Platzer | Logical Analysis of Hybrid Systems: Proving Theorems for Complex Dynamics[END_REF]. This compositional verification principle helps scaling up verification, because KeYmaera verifies a big system by verifying properties of subsystems. Strong theoretical properties, including relative completeness results, have been shown about dL [START_REF] Platzer | Differential Dynamic Logic for Hybrid Systems[END_REF][START_REF] Platzer | The Complete Proof Theory of Hybrid Systems[END_REF] indicating how this composition principle can be successful.

IV. ROBOTIC GROUND VEHICLE NAVIGATION

The robotics community has come up with a large variety of robot designs, which differ not only in their tool equipment, but also (and more importantly for the discussion in this paper) in their kinematic capabilities. We focus on wheel-based vehicles. In order to make our models applicable to a large variety of robots, we use only limited control options (e. g., do not move sideways to avoid collisions). We consider robots that drive forward (non-negative translational velocity) in sequences of arcs in two-dimensional space. 3 Such trajectories can be realized by robots with single wheel drive, differential drive, Ackermann drive, synchro drive, or omni drive [START_REF] Bräunl | Driving robots[END_REF].

Many different navigation and obstacle avoidance algorithms have been proposed for such robots, e. g. dynamic window [START_REF] Fox | The dynamic window approach to collision avoidance[END_REF], potential fields [START_REF] Khatib | Real-time obstacle avoidance for manipulators and mobile robots[END_REF], or velocity obstacles [START_REF] Fiorini | Motion Planning in Dynamic Environments Using Velocity Obstacles[END_REF]. For an introduction to various navigation approaches for mobile robots, see [START_REF] Bonin-Font | Visual Navigation for Mobile Robots: A Survey[END_REF][START_REF] Choset | Principles Of Robot Motion[END_REF]. In this paper, we focus on the dynamic window algorithm [START_REF] Fox | The dynamic window approach to collision avoidance[END_REF], which is derived from the motion dynamics of the robot and thus discusses all aspects of a hybrid system (models of discrete and continuous dynamics). Other approaches, such as velocity obstacles [START_REF] Fiorini | Motion Planning in Dynamic Environments Using Velocity Obstacles[END_REF], are interesting for further work on provably safe moving obstacle avoidance.

We want to prove motion safety of a robot that avoids obstacles by dynamic window navigation. Under passive safety [START_REF] Maček | Towards Safe Vehicle Navigation in Dynamic Urban Scenarios[END_REF], the vehicle is in a safe state if it is able to come to a full stop before making contact with an obstacle (i. e., the vehicle does not itself collide with obstacles, so if a collision occurs at all then while the vehicle was stopped). Passive safety, however, puts the burden of avoiding collisions mainly on other objects. The safety condition that we want to prove additionally is the stronger passive friendly safety [START_REF] Maček | Towards Safe Vehicle Navigation in Dynamic Urban Scenarios[END_REF]: we want to guarantee that our robot will come to a full stop safely under all circumstances, but will also leave sufficient maneuvering room for moving obstacles to avoid a collision. 4 We consider both models and both safety properties to show the differences between the assumed knowledge and the safety guarantees that can be made. The verification effort and complexity difference is quite instructive: passive safety can be guaranteed by proving safety of all robot choices, whereas passive friendly safety additionally requires liveness proofs for 3 If the radius of such a circle is infinite, the robot drives (forward) on a straight line. 4 The robot ensures that there is enough room for the obstacle to stop before a collision occurs. If the obstacle decides not to, the obstacle is to blame and our robot is still considered safe. the obstacle. In the following sections, we discuss models and verification of the dynamic window algorithm in detail.

A. Passive Safety of Obstacle Avoidance

The dynamic window algorithm is an obstacle avoidance approach for mobile robots equipped with synchro drive [START_REF] Fox | The dynamic window approach to collision avoidance[END_REF] but can be used for other drives too [START_REF] Brock | High-speed navigation using the global dynamic window approach[END_REF]. It uses circular trajectories that are uniquely determined by pairs of translational and rotational velocity (v, ω). The algorithm is roughly organized into two steps: (i) The range of all possible pairs of translational and rotational velocities is reduced to admissible ones that result in safe trajectories (i. e., avoid collisions since those trajectories allow the robot to stop before it reaches the nearest obstacle). The admissible pairs are further restricted to those that can be realized by the robot within a short time frame (the so-called dynamic window). If the set of admissible and realizable velocities is empty, the algorithm stays on the previous safe curve (such curve exists unless the robot started in an unsafe state). (ii) Progress towards the goal is optimized by maximizing a goal function. For safety verification, we can omit the second step and verify the stronger property that all choices that are fed into the optimization are safe. a) Modeling: We develop a model of the principles in the dynamic window algorithm as a hybrid program in dL. The dynamic window algorithm safely abstracts the robot's shape to a single point, since other shapes reduce to adjusting the (virtual) shapes of the obstacles (cf. [START_REF] Minguez | Abstracting Vehicle Shape and Kinematic Constraints from Obstacle Avoidance Methods[END_REF] for an approach to attribute robot shape to obstacles). We also use this abstraction to reduce the verification complexity. The robot has state variables describing its current position p r = (p x r , p y r ), translational velocity v r ≥ 0, translational acceleration a r , a direction vector d r = (cos θ, sin θ), 5 and angular velocity θ = ω r . 6The translational and rotational velocities are linked w.r.t. the rigid body planar motion by the formula p r -p c ω r = v r , where the radius p r -p c is the distance between the robot and the center of its current curve p c = (p x c , p y c ). Following [START_REF] Platzer | Differential-algebraic Dynamic Logic for Differential-algebraic Programs[END_REF], we encode sine and cosine functions in the dynamics using the extra variables d x r = cos θ and d y r = sin θ to avoid undecidable arithmetic. The continuous dynamics for the dynamic window algorithm [START_REF] Fox | The dynamic window approach to collision avoidance[END_REF] can be described by the differential equation system of ideal-world dynamics of the planar rigid body motion (p

x r = v r d x r , p y r = v r d y r , v r = a r , ( p r -p c ω r ) = a r , d x r = -ω r d y r , d y r = ω r d x r
) where the condition ( p r -p c ω r ) = a r encodes the rigid body planar motion p r -p c ω r = v r that we consider. The dynamic window algorithm assumes direct control of the translational velocity v r . We, instead, control acceleration a r and do not perform instant changes of the velocity. Our model is closer to the actual dynamics of a robot, which cannot really change its velocity instantly. The realizable velocities, then, follow from the differential equation system. The robot starts at p r = (0, 0) as initial position, with v r = 2 as initial translational velocity and ω r = 1 as initial rotational velocity; Figure 1(d) shows the resulting circular trajectory. Figure 1(b) and Figure 1(e) show the resulting curve when braking (the robot brakes along the curve and comes to a complete stop before completing the circle). If the rotational velocity is constant (ω r = 0), the robot drives an Archimedean spiral with the translational and rotational accelerations controlling the spiral's separation distance (a r /ω 2 r ). The corresponding plots are shown in Figures 1(c) and 1(f).

As in the dynamic window algorithm, we assume bounds for the acceleration a r in terms of a maximum acceleration A ≥ 0 and braking power b > 0, as well as a bound Ω on the rotational velocity ω r . We use ε to denote the upper bound for the control loop time interval (e. g., sensor and actuator delays, sampling rate, and computation time). That is, the robot may react as quickly as it wants, but it can take no longer than ε. Notice that, without such a time bound, the robot would not be safe, because its control might never run. In our model, all these bounds will be used as symbolic parameters and not concrete numbers. Therefore, our results apply to all values of these parameters and can be enlarged to include uncertainty.

An obstacle has vector state variables describing its current position p

o = (p x o , p y o ) and velocity v o = (v x o , v y o )
. The obstacle model is very liberal. The only restriction about the dynamics is, that within ε time units, they move continuously with bounded velocity v o ≤ V . Note, that the dynamic window algorithm considers a special case V = 0 (obstacles are stationary). Depending on the relation of V to ε, moving obstacles can make quite a difference, e. g., when fast obstacles meet communication-based virtual sensors as in RoboCup. 7 In order to determine admissible velocities, the dynamic window algorithm requires the distance to the nearest obstacle for every possible curve. In the presence of moving obstacles, however, all obstacles must be considered and tested for safety (e. g., in a loop). To capture this requirement, our model nondeterministically picks any obstacle p o := ( * , * ) and tests its safety, which includes safety for the worst-case behavior V of the nearest obstacle (ties are included) and is thus safe for all possible obstacles. In the case of non-point obstacles, p o denotes the obstacle perimeter point that is closest to the robot (this fits naturally to obstacle point sets delivered by radar sensors, from which the closest point will be chosen). At each active step of the robot, the position p o is updated (again nondeterministically to summarize all obstacles). In this process, the robot may or may not find another safe trajectory. If it does, the robot can follow that safe trajectory w.r.t. any nondeterministically chosen obstacle (again, V that all other obstacles will stay more distant than the worst case of the nearest). If not, the robot can still brake on the previous trajectory, which is known to be safe. Model 1 represents the common controller-plant model: it repeatedly executes control choices followed by dynamics, cf. [START_REF] Althoff | Safety assessment of robot trajectories for navigation in uncertain and dynamic environments[END_REF]. For the sake of clarity we restrict the study to circular trajectories with constant positive radius, that is p r -p c > 0, which also covers straight line trajectories if the radius is infinite. Thus, the condition ( p r -p c ω r ) = a r can be rewritten as ω r = ar pr-pc . We have bundled adjusted dynamics for spinning behavior ( p r -p c = 0, ω r = 0) and Archimedean spiral (ω r = 0, a r = 0) with KeYmaera. The control of the robot is executed in parallel to that of the obstacle, cf. [START_REF] Bonin-Font | Visual Navigation for Mobile Robots: A Survey[END_REF]. The obstacle itself may choose any velocity

?curve ∧ safe) (9) curve ≡ p r -p c > 0 ∧ ω r p r -p c = v r ( 10
)
∧ d r = (p r -p c ) ⊥ p r -p c safe ≡ p r -p o ∞ > v 2 r 2b + A b + 1 A 2 ε 2 + εv r (11)
+ V ε + v r + Aε b ( 12
)
dyn ≡ (t := 0;

p x r = v r d x r , p y r = v r d y r , (13)
d x r = -ω r d y r , d y r = ω r d x r , (14)
p x o = v x o , p y o = v y o , (15)
v r = a r , ω r = a r p r -p c , t = 1 (16) & v r ≥ 0 ∧ t ≤ ε) (17)
in any direction up to the worst-case velocity V assumed about obstacles ( v o ≤ V ), cf. (3). This uses the modeling pattern from Section III: we assign an arbitrary value to the obstacle's velocity (v o := ( * , * )), which is then restricted to any value up to the worst-case velocity using a subsequent test (?

v o ≤ V ).
The robot is allowed to brake at all times since (4) has no test. If the robot is stopped, it may choose to stay in its current spot, cf. [START_REF] Brock | High-speed navigation using the global dynamic window approach[END_REF]. Finally, the robot may choose a new safe curve in its dynamic window: it chooses any acceleration within the bounds of its braking power and acceleration [START_REF] Choset | Principles Of Robot Motion[END_REF], and any rotational velocity in the bounds, cf. [START_REF] Fiorini | Cleaning and Household Robots: A Technology Survey[END_REF]. This corresponds to testing all possible acceleration and rotational velocity values at the same time. An actual implementation would use loops to enumerate all possible values and all obstacles and test each pair (v r , ω r ) separately w.r.t. every obstacle. The admissible pairs would be stored in a data structure (as e. g., in [START_REF] Täubig | Guaranteeing functional safety: design for provability and computer-aided verification[END_REF]).

The curve is determined by the robot following a circular trajectory around the rotation center p c , starting in direction d r with angular velocity ω r , cf. [START_REF] Fiorini | Motion Planning in Dynamic Environments Using Velocity Obstacles[END_REF]. The distance to the nearest obstacle on that curve is measured in [START_REF] Fox | The dynamic window approach to collision avoidance[END_REF]. The trajectory starts at p r tangential to the circle centered at p c -which has positive radius-with translational velocity v r and rotational velocity ω r as defined by ω r p r -p c = v r in (10). 8 A circular trajectory around the rotation center p c ensures passive safety if it allows the robot to stop before it collides with the nearest obstacle. Consider the worst case where the center p c is at infinity and the robot travels on a straight line. In this case, the distance between the robot's current position p r and the nearest obstacle p o must account for the following components: First, the robot needs to be able to brake from its current velocity v r to a complete stand still, cf. [START_REF] Khatib | Real-time obstacle avoidance for manipulators and mobile robots[END_REF]. Second, the robot may not react before ε time units. Thus, we must additionally take into account the distance that the robot may travel w.r.t. the worst-case acceleration A and the distance needed for compensating its potential acceleration of A during that time with braking power b, cf.

v 2 r 2b = vr/b 0 (v r -bt)dt in
A b + 1 A 2 ε 2 + εv r = ε 0 (v r +At)dt+
Aε/b 0 (v r +Aε-bt)dt in [START_REF] Khatib | Real-time obstacle avoidance for manipulators and mobile robots[END_REF]. Third, during the worst-case time (ε + vr+Aε b ) entailed by [START_REF] Khatib | Real-time obstacle avoidance for manipulators and mobile robots[END_REF], the obstacle may approach the robot in the worst case on a straight line with maximum velocity V , cf. [START_REF] Loos | Adaptive Cruise Control: Hybrid, Distributed, and Now Formally Verified[END_REF]. To simplify the proof, we measure the distance between the robot's position p r and the obstacle's position p o in the infinitynorm p r -p o ∞ (i. e., both |p x r -p x o | and |p y r -p y o | must be safe) and thus over-approximate the Euclidean norm distance

p r -p o 2 = (p x r -p x o ) 2 + (p y r -p y o ) 2 by a factor of at most √ 2.
Finally, the continuous dynamics of the robot and the obstacle are defined in ( 13)-( 17). b) Verification: We verify the safety of the control algorithm modeled as a hybrid program in Model 1, using a formal proof calculus for dL [START_REF] Platzer | Differential Dynamic Logic for Hybrid Systems[END_REF][START_REF] Platzer | Differential-algebraic Dynamic Logic for Differential-algebraic Programs[END_REF][START_REF] Platzer | Logical Analysis of Hybrid Systems: Proving Theorems for Complex Dynamics[END_REF]. The robot is at a safe distance from the obstacle, if it is able to brake to a complete stop at all times before the approaching obstacle reaches the robot. The following condition captures this requirement as a property that we prove to hold for all executions:

φ ps ≡ (v r = 0) ∨ p r -p o > v 2 r 2b + V v r b . (18)
The formula [START_REF] Platzer | Differential Dynamic Logic for Hybrid Systems[END_REF] states that the robot is stopped or the robot and the obstacle have different positions safely apart. This accounts for the robot's braking distance v 2 r 2b while the obstacle is allowed to approach the robot with its worst case travel distance V vr b . We prove that the property [START_REF] Platzer | Differential Dynamic Logic for Hybrid Systems[END_REF] holds for all executions of Model 1 under the assumption that we start in a state satisfying the following conditions: 9

ψ ps ≡ φ ps ∧ ( p r -p c > 0) ∧ ( d r = 1) . (19)
The first condition of the conjunction formalizes the fact that we are starting in a position that is passively safe, that is φ ps holds. The second conjunct states that the robot is not initially spinning. The last conjunct d r = 1 says that the direction d r is a unit vector.

Theorem 1 (Passive safety): If the robot starts in a state where ψ ps [START_REF] Platzer | Differential-algebraic Dynamic Logic for Differential-algebraic Programs[END_REF] holds, then the control model dw ps (Model 1) always guarantees the passive safety condition (φ ps ), as expressed by the provable dL formula: ψ ps → [dw ps ]φ ps .

We proved Theorem 1 for circular trajectories, spinning, and spiral trajectories using KeYmaera, a theorem prover for hybrid systems. In the next section, we explore the stronger requirements of passive friendly safety, where the robot will not only safely stop itself, but also allow for the obstacle to stop before a collision occurs. 9 The formal proof uses the additional constraints stated earlier, vr ≥ 0, A ≥ 0, V ≥ 0, Ω ≥ 0, b > 0, and ε > 0, which we leave out for simplicity.

B. Passive Friendly Safety of Obstacle Avoidance

Passive friendly safety, as introduced above, requires the robot to take careful decisions that respect the features of moving obstacles. The definition of Maček et al. [START_REF] Maček | Towards Safe Vehicle Navigation in Dynamic Urban Scenarios[END_REF] requires that the robot respects the worst-case braking time of the obstacle. In our model, the worst-case braking time is represented as follows. We assume an upper bound τ on the reaction time of the obstacle and a lower bound b o on its braking capabilities. Then, τ V is the worst-case distance that the obstacle can travel before actually reacting and V 2 2bo is the distance for the obstacle to stop from the worst-case velocity V with an assumed minimum braking capability b o . c) Modeling: Model 2 uses the same obstacle avoidance algorithm as Model 1. The essential difference reflects what the robot considers to be a safe distance to an obstacle. As shown in [START_REF] Platzer | The Complete Proof Theory of Hybrid Systems[END_REF] the distance not only accounts for the robot's own braking distance, but also for the braking distance V 2 2bo and reaction time τ of the obstacle. The verification of passive friendly safety is more complicated than passive safety as it accounts for the behavior of the obstacle discussed below.

Model 2 Dynamic window with passive friendly safety dw pfs ≡ (ctrl; dyn) * (see Model 1 for details) [START_REF] Platzer | Logical Analysis of Hybrid Systems: Proving Theorems for Complex Dynamics[END_REF] safe

≡ p r -p o ∞ > v 2 r 2b + V 2 2b o + τ V (21)
+ A b + 1 A 2 ε 2 + εv r + V ε + v r + Aε b d) Verification:
We verify the safety of the robot's control choices as modeled in Model 2. Unlike the passive safety case, the passive friendly safety property φ pfs should guarantee that if the robot stops, moving obstacles (cf. Model 3) still have enough time and space to avoid a collision. This requirement can be captured by the following dL formula:

η pfs ≡ (v r = 0) ∧ p r -p o > V 2 2b o + τ V ∧ (0 ≤ v o ≤ V ) → obstacle (( p r -p o > 0) ∧ (v o = 0)) . ( 22
)
Equation ( 22) says that, once the robot stops (v r = 0), there exists an execution of the hybrid program obstacle, (existence of a run is formalized by the diamond operator obstacle ), that allows the obstacle to stop (v o = 0) before a collision occurs ( p r -p o > 0). Passive friendly safety is now stated as φ pfs ≡ (ϕ ps ∨ v r = 0) ∧ η pfs , where the property ϕ ps is analogue to the passive safety property φ ps w.r.t. Model 2:

ϕ ps ≡ p r -p o > v 2 r 2b + V 2 2b o + V v r b + τ .
We study the passive friendly safety with respect to the initial feasible states satisfying the following property:

ψ pfs ≡ ϕ ps ∧ ( p r -p c > 0) ∧ ( d r = 1) . (23)
Observe that, in addition to the condition η pfs , the difference with the passive safety condition is reflected in the special treatment of the case v r = 0. Indeed, if the robot would start with a null translational velocity (which is passive safety) while not satisfying ϕ ps , then we cannot prove the passive friendly safety as the obstacle may be unsafely close already initially. Besides, we can see in φ pfs that we are required to prove η pfs even when the robot comes to a full stop.

In Model 2 the hybrid program ctrl o is a coarse model given by equation ( 3), which only constrains its non-negative velocity to be less than or equal to V . Such an obstacle could trivially prevent a collision by stopping instantaneously. In η pfs , we consider a more interesting refined obstacle behavior modeled by the hybrid program given in Model 3. (25)

a o := * ; ?v o + a o ε o ≤ V (26) dyn or ≡ (t := 0; p x o = v o d x o , p y o = v o d y o , (27)
v o = a o , t = 1 & t ≤ ε o ∧ v o ≥ 0)
The refined obstacle may choose any direction as described by the unit vector d o in [START_REF] Täubig | Guaranteeing functional safety: design for provability and computer-aided verification[END_REF] and any acceleration a o , as long as it does not exceed the velocity bound V , cf. [START_REF] Van Den | LQG-MP: Optimized Path Planning for Robots with Motion Uncertainty and Imperfect State Information[END_REF]. The dynamics of the obstacle are straight ideal-world translational motion in the two-dimensional plane, see [START_REF] Wu | Guaranteed infinite horizon avoidance of unpredictable, dynamically constrained obstacles[END_REF]. We verified Theorem 2 in KeYmaera. The symbolic bounds on velocity, acceleration, braking, and time in the above models represent uncertainty implicitly (e. g., one may instantiate the braking power b with the minimum specification of the robot's brakes, or with the actual braking power achievable w.r.t. the current terrain). Dually, whenever knowledge about the current state is available, the bounds can be instantiated more aggressively to allow efficient robot behavior (e. g., in a rare worst case we may face a particularly fast obstacle, but right now there are only slow-moving obstacles around). Theorems 1 and 2 are verified for all those values. Other aspects of uncertainty need explicit changes in the models and proofs, as discussed in the next section.

C. Safety of Obstacle Avoidance Despite Uncertainty

Robots have to deal with uncertainty in almost every aspect of their interaction with the environment, ranging from sensor inputs (e. g., inaccurate localization, distance measurement) to actuator effects (e. g., wheel slip depending on the terrain). In this section, we show two examples of uncertainty explicitly handled by the models. First, we allow localization uncertainty (the robot knows its position only approximately). We then consider imperfect actuator commands (braking and acceleration) considered to be within an interval. Such intervals are instantiated, e. g., according to sensor or actuator specification (e. g., GPS error), or w.r.t. experimental measurements.

(29) curve ≡ pr -p c > 0 ∧ ω r pr -p c = v r (30) ∧ d r = (p r -p c ) ⊥ pr -p c safe ≡ pr -p o ∞ > v 2 r 2b + A b + 1 A 2 ε 2 + εv r (31)
+ V ε + v r + Aε b + U p
Model 4 introduces location uncertainty. It adds a location measurement pr before the control decisions are made. This location measurement may deviate from the real position p r by no more than symbolic parameter U p , cf. (29). The measured location pr is used in all control decisions of the robot (e. g., in (30) to compute the curve's center point and rotational velocity), while the robot's physical motion is still computed on the symbolic real position p r . We prove in KeYmaera that the robot is still safe (passive safety), if it accounts for the location uncertainty as stated in the safety constraint (31). dyn with a r replaced by ãr Model 5 introduces uncertainty in actuation with an actuator perturbation between control and dynamics, cf. (32). Actuator perturbation affects the acceleration by a damping factor u m , known to be within [U m , 1], cf. (33). Note, that this damping factor can change arbitrarily often, but is assumed to be constant for ε time units. In the worst case, the robot has full acceleration (u m = 1) but fully reduced braking (u m = U m ). For instance, the robot accelerates on perfect terrain, but is on slippery terrain whenever it brakes. The robot considers this worst case scenario during control in its safety constraint (34). Using KeYmaera, we prove that the robot is still safe (passive safety) when it follows Model 5.

V. CONCLUSION AND FUTURE WORK

Robots are hybrid systems, because they share continuous physical motion with complicated computer algorithms controlling their behavior. We demonstrate that this understanding also helps proving robots safe. We develop hybrid system models of the dynamic window algorithm for autonomous ground vehicles and prove that the algorithm guarantees both passive safety and passive friendly safety in the presence of moving obstacles. All the proofs were achieved automatically with some manual guidance using the verification tool KeYmaera. We manually provided inductive and differential invariants (crucial), as well as reduced the number of terms and variables with simple interactions to reduce arithmetic complexity. 85% of the proof steps were automatic. Most interactive steps were simple arithmetic simplifications that KeYmaera's current proof strategies do not yet automate. We augment these models and safety proofs with robustness for localization uncertainty and imperfect actuation. We observe that incremental revision of models and proofs helps reducing the verification complexity and understanding the impact of uncertainty on the behavior of the robot. Further, note that all symbolic bounds in our models-such as those for maximum obstacle velocity and sensor/actuator uncertainty-can be instantiated with probabilistic models (e. g., assume the 2σ confidence interval of the distribution of obstacle velocities as maximum obstacle velocity). In this case, our verified safety guarantees translate into a safety probability.

Future work includes refined kinematic capabilities (e. g., going sideways with omni-drive) and additional sources of uncertainties, such as distance measurements. Also, to derive less conservative safety bounds for specific environments, one could introduce specific obstacle models (e. g., pedestrians).

  Position (p x r , p y r ), translational velocity vr and rotational velocity ωr for positive acceleration on a circle.

  Position (p x r , p y r ), translational velocity vr and rotational velocity ωr for braking to a complete stop on a circle.

  Position (p x r , p y r ), translational velocity vr and rotational velocity ωr for translational acceleration on a spiral.

  (p x r , p y r ) motion plot for acceleration (a).

5 (

 5 e) (p x r , p y r ) motion plot for braking (b).

  (p x r , p y r ) motion plot for (c).

Fig. 1 :

 1 Fig. 1: Trajectories of the robot.

Figure 1 (

 1 Figure 1(a) depicts the position and velocity changes of a robot accelerating on a circle around a center point p c = (2, 0). The robot starts at p r = (0, 0) as initial position, with v r = 2 as initial translational velocity and ω r = 1 as initial rotational velocity; Figure 1(d) shows the resulting circular trajectory. Figure 1(b) and Figure 1(e) show the resulting

Model 1

 1 7 http://www.robocup.org/ Dynamic window with passive safetydw ps ≡ (ctrl; dyn) * (1) ctrl ≡ ctrl o || ctrl r (2) ctrl o ≡ v o := ( * , * ); ? v o ≤ V ctrl r ≡ (a r := -b)(4) ∪ (?v r = 0; a r := 0; ω r := 0)(5) ∪ (a r := * ; ?-b ≤ a r ≤ A;(6) ω r := * ; ? -Ω ≤ ω r ≤ Ω; (7) p c := ( * , * ); d r := ( * , * ); (8) p o := ( * , * );

  8 dr = (pr -pc) ⊥ pr -pc means component-wise equality of the orthonormal vector, normalized by the radius, as in d x r = -

Model 3

 3 Refined obstacle controls acceleration obstacle ≡ (ctrl or ; dyn or ) * (24) ctrl or ≡ d o := ( * , * ); ? d o = 1;

Theorem 2 (

 2 Passive friendly safety): If the property ψ pfs (23) holds initially, then the control model dw pfs (Model 2), guarantees the passive friendly safety condition φ pfs in presence of obstacles per Model 3, as expressed by the provable dL formula: ψ pfs → [dw pfs ]φ pfs .

Model 5

 5 Safety despite actuator uncertainty dw uas ≡ (ctrl; act; dyn) * (see Model 1 for details) (32) act ≡ u m := * ; ãr := u m a r ; ?0 < U m ≤ u m ≤ 1 (33) safe ≡ p r -p o ∞

TABLE I :

 I Hybrid program representations of hybrid systems.

	Statement	Effect
	α; β	sequential composition, first run α, then β
	α ∪ β	nondeterministic choice, following either α or β
	α *

  10

  10 Instantiation with probabilistic bounds means that the symbolically guaranteed safety is traded for a safety probability. Safety despite uncertain position measurements dw ups ≡ (locate; ctrl; dyn) * (see Model 1 for details) (28) locate ≡ pr := ( * , * ); ? pr -p r ≤ U p

	Model 4

http://symbolaris.com/info/KeYmaera.html

As stated earlier, we study unidirectional motion: the robot moves along its direction, that is the vector dr gives the direction of the velocity vector.

The derivative with respect to time is denoted by superscript prime ( ).

ACKNOWLEDGMENTS

This material is based upon work supported by NSF CAREER Award CNS-1054246 and NSF EXPEDITION CNS-0926181, by DARPA FA8750-12-2-0291, and by US DOT UTC TSET award # DTRT12GUTC11. This work was also supported by the Austrian BMVIT under grant FIT-IT 829598, FFG BRIDGE 838526, and FFG Basisprogramm 838181.


Khalil Ghorbal

Franjo Ivancic

Gogul Balakrishnan

Naoto Maeda

Aarti Gupta

Franjo Ivančić

Donut Domains: Efficient Non-convex Domains for Abstract Interpretation

à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Efficient program analysis using abstract interpretation [START_REF] Cousot | Abstract Interpretation: A unified lattice model for static analysis of programs by construction or approximation of fixpoints[END_REF] typically uses convex domains such as intervals, octagons, zonotopes or polyhedra [START_REF] Cousot | Static determination of dynamic properties of programs[END_REF][START_REF] Cousot | Automatic discovery of linear restraints among the variables of a program[END_REF][START_REF] Ghorbal | The zonotope abstract domain taylor1+[END_REF][START_REF] Goubault | Static analysis of numerical algorithms[END_REF][START_REF] Miné | The octagon abstract domain[END_REF]. However, certain properties of interest require reasoning about non-convex structures. One approach to non-convex reasoning is to utilize powerset domains of elementary convex domains [START_REF] Bagnara | Widening operators for powerset domains[END_REF][START_REF] Gurfinkel | Boxes: A symbolic abstract domain of boxes[END_REF][START_REF] Halbwachs | Verification of linear hybrid systems by means of convex approximations[END_REF]. In general, it has proved to be difficult to provide satisfactory improvements over elementary convex domains with powerset domains while maintaining small enough performance degradation. Furthermore, it would be difficult to maintain enough disjunctions in the powerset depending on the particular non-convex shape being approximated. Note, however, that the recently proposed Boxes domain by Gurfinkel and Chaki [START_REF] Gurfinkel | Boxes: A symbolic abstract domain of boxes[END_REF] can potentially represent exponentially many interval constraints compactly. It utilizes a BDD-like extension to elementary range constraints called LDD [START_REF] Chaki | Decision diagrams for linear arithmetic[END_REF]. However, we are considering relational domains such as octagons, zonotopes or polyhedra as well.

Additional non-convex domains based on congruence (either linear [START_REF] Granger | Static analysis of linear congruence equalities among variables of a program[END_REF] or trapezoid [START_REF] Masdupuy | Array abstractions using semantic analysis of trapezoid congruences[END_REF]) analysis have been developed. Such domains catch a congruence relation that variables satisfy and are suitable for the analysis of indexes of arrays for instance. Recent work by Chen et al. considered a polyhedral abstract domain with interval coefficients [START_REF] Chen | Interval polyhedra: An abstract domain to infer interval linear relationships[END_REF]. This abstract domain has the ability to express certain non-convex invariants. For example, in this domain some multiplications can be evaluated precisely (see Figure 6 for instance). Other interesting nonconvex abstract domains were introduced to capture specific invariants such as min-max invariants [START_REF] Allamigeon | Inferring Min and Max Invariants Using Max-plus Polyhedra[END_REF] and quadratic templates [START_REF] Adjé | Coupling policy iteration with semi-definite relaxation to compute accurate numerical invariants in static analysis[END_REF].

We address a different type of non-convexity commonly occurring in software, which relates to small sub-regions of instability within a normal operating (convex) region of interest. The non-convex region of values that may cause the bug is (under-)approximated using a convex inner region (or hole) that is subtracted from a convex outer region. We call this representation donut domains. Our approach relies on the usual operations defined on (convex) sub-domains, except for the need to compute under-approximations in domain A 2 . The donut domains can be considered extensions to the work on signed types domain introduced by the authors in [START_REF] Prabhu | Interprocedural exception analysis for C++[END_REF]. There, we start with a finite set of types, and allow a set-minus operation only from the universal set.

Under-Approximations of Polyhedra. Under-approximations have been utilized for applications such as test vector generation and counterexample generation, by providing must-reach sets.

Bemporad et al. introduced the notion of inner approximations of polyhedra using intervals in [START_REF] Bemporad | Inner and outer approximations of polytopes using boxes[END_REF]. In [START_REF] Kanade | Generating and analyzing symbolic traces of Simulink/Stateflow models[END_REF], polyhedra are under-approximated for test vector generation of Simulink/Stateflow models using bounded vertex representation (BVR). Goubault and Putot describe a method to compute an underapproximating zonotope [START_REF] Goubault | Under-approximations of computations in real numbers based on generalized affine arithmetic[END_REF] using modal intervals [START_REF] Goldsztejn | Modal intervals revisited: a mean-value extension to generalized intervals[END_REF] for non-linear operations.

In this work, we propose a novel technique to find under-approximations of polyhedra based on a fixed template. We first re-formulate the problem by introducing an auxiliary matrix. This matrix represents the fact that we are looking for an inner polyhedral object of a particular shape. Using this auxiliary matrix re-formulation, we can then use classical convex analysis techniques to characterize an under-approximations of polyhedra.

Motivating Example. Figure 1 highlights a code snippet taken from XTide 3 . The XTide package provides accurate tide and current predictions in a number of formats based on algorithms. Similar patterns may exist in controller-related software to avoid regions of controller or numerical instability.

After the step marked initializations, (dx, dy) could be any point in R 2 except the origin (0, 0), this particular point is kept and propagated forward as a "hole". After the if-statement, the set of reachable values is: (dy > dx ∧ dy > -dx) ∨ (-dy > dx ∧ -dy > -dx). The above region is non-convex; therefore, any classical abstract domain will end up at this control point with for both variables. Moreover, here, the interpretation of the strict inequality of the test is required to prove that dx = 0. The else case is even harder: in addition to the non-convexity of the set of possible values, one needs to consider the full-zero-test together with the negation of |dy| > |dx|, to prove that the division by dy is safe. All known non-convex domains fail to prove the safety of such example. Typically, the powerset domain cannot encode the full-zero-test which is mandatory here (for the else branch).

Content. The rest of this paper is organized as follows: we define, Section 2, a new set of domains called donut domains. Section 3 proposes a novel method to compute polyhedral under-approximations for arbitrary linear templates. Finally, in Section 4, first experiments and promising results are discussed.

Donut Abstract Domains

In this section we introduce donut domains, and define the operation on donut domains based on operations in the component domains.

Lattice Structure

Let (A 1 , ≤ 1 , ∪ 1 , ∩ 1 , ⊥ 1 , 1 , γ 1 ) and (A 2 , ≤ 2 , ∪ 2 , ∩ 2 , ⊥ 2 , 2 , γ 2
) denote two classical numerical abstract domains, where ≤ , ∪ ,∩ ,⊥ , ,γ denote the partial order, the join and meet operations, the bottom and top elements and the concretization function of the classical abstract domain for ∈ {1, 2}, respectively. In this work, we extend a given abstract domain with an under-approximation operator α, such that for any concrete object X, we have γ • α(X) ⊆ X. An abstract object X 1\2 of the domain A 1 \ A 2 is defined by a pair of objects (X 1 , X 2 ), such that X 1 ∈ A 1 and X 2 ∈ A 2 . The object X 1\2 abstracts the set of possible values reached by the variables as follows:

-The object X 1 ∈ A 1 represents an over-approximation of the set of reachable values.

-The object X 2 ∈ A 2 represents an under-approximation of the set of unreachable values (usually within γ 1 (X 1 )).

The concretisation function is defined as follows.

γ 1\2 def = γ 1\2 (X 1 , X 2 ) def = γ 1 (X 1 ) \ γ 2 (X 2 ) .
Figure 2 depicts a concretization of a typical donut object where the domain A 1 is the Affine Sets domain [START_REF] Ghorbal | A logical product approach to zonotope intersection[END_REF] and A 2 is the octagons domain. One should keep in mind the implicit set of unreachable values implied by γ 1 (X 1 ) -namely R p \ γ 1 (X 1 ) denoted in the sequel by γ1 (X 1 ). Indeed, the set of unreachable values is actually γ1 (X 1 ) ∪ γ 2 (X 2 ). As said earlier, γ 2 (X 2 ) is a (convex) under-approximation of the set of unreachable values. The fact that the intersection γ 1 (X 1 ) ∩ γ 2 (X 2 ) is not empty permits to encode a hole inside γ 1 (X 1 ) (see Figure 2).

Interval Concretisation. The interval concretization of the variable x

k , 1 ≤ k ≤ p, denoted by [x k ], is defined by π k (γ 1 (X 1 ) \ γ 2 (X 2 ))
, where π k denotes the orthogonal projection of a given set onto dimension k. Note that

[x k ] ⊇ π k (γ 1 (X 1 )) \ π k (γ 2 (X 2 )). For instance in ([-2, 2] × [-2, 2], [-1, 1] × [-∞, +∞]), we have [x 2 ] = [-2, 2], whereas [-2, 2] \ [-∞, +∞] = ∅.
We embed A 1 \ A 2 with a binary relation and prove that it is a pre-order.

Definition 1. Given X 1 , Y 1 ∈ A 1 and X 2 , Y 2 ∈ A 2 , we say that (X 1 , X 2 ) is less than or equal to (Y 1 , Y 2 ) denoted by (X 1 , X 2 ) ≤ 1\2 (Y 1 , Y 2 ) if and only if X 1 ≤ 1 Y 1 and γ1 (X 1 ) ∪ γ 2 (X 2 ) ⊇ γ1 (Y 1 ) ∪ γ 2 (Y 2 ) . (1)
Proposition 1. The binary relation ≤ 1\2 is a pre-order over

A 1 \ A 2 . It defines an equivalence relation ∼ defined by (X 1 , X 2 ) ≤ 1\2 (Y 1 , Y 2 ) and (Y 1 , Y 2 ) ≤ 1\2 (X 1 , X 2 ) and characterized by X 1 = Y 1 (X 1 ≤ 1 Y 1 and Y 1 ≤ 1 X 1 ), γ 2 (X 2 ) ⊆ γ 2 (Y 2 ) ∪ γ1 (Y 1 ) and γ 2 (Y 2 ) ⊆ γ 2 (X 2 ) ∪ γ1 (X 1
). We reuse the symbol ≤ 1\2 to also denote the partial order quotiented by the equivalence relation ∼.

With respect to ≤ 1\2 , we have

(⊥ 1 , ⊥ 2 ) ∼ (⊥ 1 , 2 ) ≤ 1\2 ( 1 , 2 ) ≤ 1\2 ( 1 , ⊥ 2 );
therefore, we define the bottom and top elements of

A 1 \ A 2 by ⊥ 1\2 def = (⊥ 1 , -) 1\2 def = ( 1 , ⊥ 2 ) .

Decidability of the Order

Despite the non-convexity of γ, the equivalence class introduced in Proposition 1 suggests particular representatives of objects (X 1 , X 2 ) which are easily comparable. Indeed, γ is no longer involved when the concretization of the hole X 2 is included in the concretization of X 1 . Observe moreover, that the definition of the order relation ≤ 1\2 allows comparing two abstract objects having their holes in two different abstract domains, since only the concretization functions are involved in (1).

Proposition 2. Let (X 1 , X 2 ) and (Y 1 , Y 2 ) be two elements of A 1 \ A 2 such that γ 2 (X 2 ) ⊆ γ 1 (X 1 ), and γ 2 (Y 2 ) ⊆ γ 1 (Y 1 ). Therefore, (X 1 , X 2 ) ≤ 1\2 (Y 1 , Y 2 ) if and only if X 1 ≤ 1 Y 1 and γ 1 (X 1 ) ∩ γ 2 (Y 2 ) ⊆ γ 2 (X 2 ). The condition γ 1 (X 1 ) ∩ γ 2 (Y 2 ) ⊆ γ 2 (X 2 )
, can be checked in the abstract world rather than in the concrete domain up to the use of an expressive enough domain for both A 2 and A 1 : for instance a box and an octagon can be seen as special polyhedra and the meet operation of Polyhedra abstract domain can be used.

Let X P 1 denote the abstract representation in Polyhedra domain of the abstract object X 1 , that is α P (γ 1 (X 1 )). The general procedure to decide whether (X 1 , X 2 ) is less than or equal to (Y 1 , Y 2 ), is as follows:

1. First, we "upgrade" X 2 and Y 2 to Polyhedra domain, which gives (X 1 , X P

2 ) and (Y 1 , Y P 2 ). 2. Then, we derive our particular representatives, namely (X 1 , X

P 1 ∩ P X P 2 ) for (X 1 , X P 2 ) and (Y 1 , Y P 1 ∩ P Y P 2 ) for (Y 1 , Y P
2 ) (∩ P being the meet operation in Polyhedra domain). 3. Finally, we use Proposition 2 by checking for the inequalities X 1 ≤ 1 Y 1 and

X P 1 ∩ P Y P 1 ∩ P Y P 2 ≤ P X P 1 ∩ P X P 2 .

Meet and join operations

We start with a simple example to clarify the intuition behind the formal definition given later. What the example suggests is that when computing a join of two elements (X 1 , X 2 ) and (Y 1 , Y 2 ), we often end up with multiple (not necessarily convex nor connex) holes defined by (γ

2 (X 2 ) ∪ γ1 (X 1 )) ∩ (γ 2 (Y 2 ) ∪ γ1 (Y 1 )
), which gives by distributing the meet over the join:

(γ 2 (X 2 ) ∩ γ 2 (Y 2 )) ∪ (γ 2 (X 2 ) ∩ γ1 (Y 1 )) ∪ (γ 2 (Y 2 ) ∩ γ1 (X 1 )) ∪ (γ 1 (X 1 ) ∩ γ1 (Y 1 )
) .

An under-approximation of the final element γ1 (X 1 ) ∩ γ1 (Y 1 ) is implicit since the over-approximation of reachable values is given by X 1 ∪ 1 Y 1 . Thus, only the intersection of the first three sets will be considered (which is sound). In our example, γ(

[1, 6]) = [-∞, 1) ∪ (6, +∞], and γ([0, 3]) = [-∞, 0) ∪ (3, +∞], this gives [1, 2] ∩ [2, 5] = [2, 2] and [1, 2] ∩ ([-∞, 1) ∪ (6, +∞]) = ∅ [2, 5] ∩ ([-∞, 0) ∪ (3, +∞]) = (3, 5] .
As said earlier, the intersection ([-∞, 1)

∪ (6, +∞]) ∩ ([-∞, 0) ∪ (3, +∞]) is implicit since it is covered by γ1 ([0, 3] ∪ [1, 6]).
We now formalize the join operator:

(X 1 , X 2 ) ∪ 1\2 (Y 1 , Y 2 ) def = (X 1 ∪ 1 Y 1 , (X 1 , X 2 ) ∩(Y 1 , Y 2 )),
where ∩ is defined by:

(X 1 , X 2 ) ∩(Y 1 , Y 2 ) def = α((γ 2 (X 2 )∩γ 2 (Y 2 ))∪(γ 2 (X 2 )∩γ 1 (Y 1 ))∪(γ 2 (Y 2 )∩γ 1 (X 1 ))).
We may perform heuristic checks to prioritize which hole (if many) to keep, which may also depend on the under-approximation abstraction function α. For instance we may choose an inner approximation (if working with closed domains) of the hole [START_REF] Bagnara | Not necessarily closed convex polyhedra and the double description method[END_REF][START_REF] Bagnara | Exact join detection for convex polyhedra and other numerical abstractions[END_REF] instead of choosing the hole [START_REF] Allamigeon | Inferring Min and Max Invariants Using Max-plus Polyhedra[END_REF][START_REF] Allamigeon | Inferring Min and Max Invariants Using Max-plus Polyhedra[END_REF].

Notice also that we have a straightforward fallback operator ∩fb , that involves only X 2 and Y 2 :

X 2 ∩fb Y 2 def = α(γ 2 (X 2 ) ∩ γ 2 (Y 2 )) .
The operator is sound with respect to under-approximation. It focuses only on a particular hole, namely γ 2 (X 2 ) ∩ γ 2 (Y 2 ), instead of considering all possibilities.

In our current implementation, we use this fallback operator in a smart manner: before computing the meet of both holes, we relax, whenever possible, in a convex way, these holes. This relaxation is performed by removing all constraints that could be removed while preserving γ 1 (X 1 ). For instance, if the hole is the point (0, 0), and the abstraction of X 1 is given by the conjunction y ≥ x ∧ -y ≥ x, then the hole (0, 0) is relaxed to x ≥ 0 (see Figure 3). For the meet operation, we proceed in a similar manner. If the domain A 2 is closed under the meet operation (almost all polyhedra-like abstract domains), it is possible to replace α by α, and ∩fb by ∩ 2 . In our example, the fallback operator gives the box [START_REF] Allamigeon | Inferring Min and Max Invariants Using Max-plus Polyhedra[END_REF][START_REF] Allamigeon | Inferring Min and Max Invariants Using Max-plus Polyhedra[END_REF].

The meet operator ∩ 1\2 is defined in a similar manner:

(X 1 , X 2 ) ∩ 1\2 (Y 1 , Y 2 ) def = (X 1 ∩ 1 Y 1 , X 2 ∪Y 2 )
where

X 2 ∪Y 2 def = α2 (γ 2 (X 2 ) ∪ γ 2 (Y 2 )) .
We deliberately omit γ1 (X 1 ) ∪ γ1 (Y 1 ) in the above definition of ∪ because it is implicit from X 1 ∩ 1 Y 1 . If the domain A 2 is closed under the join operation, then ∪ is exactly equal to ∪ 2 . Very often, however, the join operation leads to an over-approximation. Therefore the detection of an exact join as in [START_REF] Bemporad | Convexity recognition of the union of polyhedra[END_REF][START_REF] Bagnara | Exact join detection for convex polyhedra and other numerical abstractions[END_REF] is of particular interest. In our current implementation, if X 2 and Y 2 overlap, we soundly extend, in a convex way, the non empty intersection. For instance, if Example 2. Consider a 2-dim simple abstract objects. Figure 4 shows a graphical representation of two overlapping objects. The remaining sub-figures highlight some of the pertinent steps with respect to the computation of ∪ 1\2 and ∩ 1\2 for such overlapping objects.

X 2 = [-2, 1] × [-1, 1] and Y 2 = [-1, 2] × [-2,

Loop widening

When processing loop elements in abstract interpretation, we may require widening to guarantee termination of the analysis. For donut domains, we extend the widening operations defined on the component abstract domains. We use the pair-wise definition of widening operators ∇. We thus define widening of donut domains as:

(X 1 , X 2 )∇ 1\2 (Y 1 , Y 2 ) = (X 1 ∇ 1 Y 1 , X 2 ∩ 2 Y 2 ) .
We use the standard widening operator ∇ 1 for abstract domain A 1 . Similarly, we use standard meet operator ∩ 2 of abstract domain A 2 for the inner region, which ensures the soundness of ∇ 1\2 . The convergence of the first component is guaranteed by the widening operator ∇ 1 . The convergence of the second component needs however more attention. Note that the simple use of narrowing operator of A 2 is unsound as it may give a donut object which is not an upper bound. To ensure the termination we add a parameter k which will encode the maximal number of allowed iterations. If the donut object does not converge within those k iterations, the hole component is reduced to ⊥ 2 . Note that the use of the narrowing operator of A 2 instead of ∩ 2 does not give in general an upper bound of (X 1 , X 2 ) and (Y 1 , Y 2 ).

Interpretation of Tests

The ability to express holes allows us to better handle a wide range of non-convex tests such as the = test or the strict inequality test. We start with classical tests. For ∈ {=, ≤} :

x k 0 (X 1 , X 2 ) def = ( x k 0 1 (X 1 ), x k 0 1 (X 2 )), where • 2 def = α2 • • 2 .
Such under-approximation is required so that the newly computed (exact) hole can be encoded in A 2 . Therefore, if the exact hole fits naturally in A 2 (say we have a linear constraint and A 2 is Polyhedra domain), there is no need to under-approximate ( • 2 = • 2 ). In Section 3, we detail how we compute such an under-approximation, whenever needed. If no algorithm is available for the under-approximation, we keep the object X 2 unchanged, which is sound.

The non-equality test = is defined as follows:

x k = 0 (X 1 , X 2 ) def = ( x k = 0 (X 1 ), α(γ 2 (X 2 ) ∪ x k = 0 2 )) .
Although x k = 0 (X 1 ) is interpreted as the identity function in standard implementations, nothing prevents the use of any available enhancement proposed by the used analyzer. For the hole, we compute the join of the new hole implied by the constraint x k = 0 together with the already existing hole X 2 . If holes γ 2 (X 2 ) and x k = 0 2 do not overlap, we discard X 2 . In fact, very often (as will be seen in experiments), the hole induced by the constraint x k = 0 is mandatory in order to prove the safety of subsequent computations. Finally, our approach offers, for free, an interesting abstraction of the strict inequality tests. A comparison with Not Necessarily Closed domains [START_REF] Bagnara | Not necessarily closed convex polyhedra and the double description method[END_REF] is planned as future work.

x k < 0 (X 1 , X 2 ) def = x k = 0 • x k ≤ 0 (X 1 , X 2 ) .

Abstract Assignment

We define in this section the abstraction of the assignment transfer function in A 1 \ A 2 . We first give an abstraction of the forget transfer function (nondeterministic assignment) :

x k ←? 1\2 (X 1 , X 2 ) def = (Y 1 , Y 2 ),
where

Y 1 def = x k ←? 1 (X 1 ) Y 2 def = x k ←? 2 (X 2 ) if γ 1 (X 1 ) ∩ γ 2 ( x k ←? 2 (X 2 )) ⊆ γ 2 (X 2 ) ⊥ 2 otherwise .
For Y 2 , we basically check whether applying the forget operator to X 2 intersects γ 1\2 (X 1 , X 2 ), by checking if this newly computed hole is included on the original hole, that is

γ 2 (X 2 ). If yes, Y 2 is set to ⊥ 2 . For instance, forget- ting x 2 in (X 1 , X 2 ) def = ([-2, 2] × [-2, 2], [-1, 1] × [-∞, +∞]) gives ([-2, 2] × [-∞, +∞], [-1, 1] × [-∞, +∞]): since x 2 ←? 2 (X 2 ) = [-1, 1] × [-∞, +∞], γ 1 (X 1 ) ∩ γ 2 ( x 2 ←? 2 (X 2 )) = [-1, 1] × [-2, 2] which is included in γ 2 (X 2 ). For- getting x 1 , however, makes Y 2 = ⊥ 2 .
The assignment could be seen as a sequence of multiple basic, already defined, operations. We distinguish two kind of assignments x ← e, where e is an arithmetic expression: (ı) non-invertible assignments, where the old values of x are lost, such as x ← c, c ∈ R, and (ıı) invertible assignments, such as x ← x + y. For non-invertible assignment, we have:

x k ← e 1\2 def = x k = e 1\2 • x k ←? 1\2 .
Invertible assignments are defined in a similar manner. It augments first the set of variables by a new fresh variable, say v, then enforces the test v = e, and finally forgets x and (syntactically) renames v to x. Notice that augmenting the set of variables in A 1 \ A 2 makes the newly added variable, v, unconstrained in both components, X 1 and X 2 . We can suppose that such a variable v already exists, and used whenever we have an invertible assignment; hence, we obtain:

x k ← e 1\2 def = swap(x k , v) in x k ←? 1\2 • v = e 1\2 .

Template-based Under-Approximations of Polyhedra

In this section we develop a new technique to under-approximate holes obtained after linear tests. Holes obtained after non-linear tests are so far reduced to ⊥ 2 , which is sound. We plan to improve this as a future work. Consider for instance the object ( 5 depicts the exact evaluation of a linear assignment. If we use boxes to encode holes, we need to compute a box inside the white polytope. In Figure 6, an under-approximation is needed for all convex domains, whereas a non-convex domain such as Interval Polyhedra [START_REF] Chen | Interval polyhedra: An abstract domain to infer interval linear relationships[END_REF] can express exactly this kind of pattern. The problem can be seen as follows: given a polyhedron P, we seek to compute a maximal (in a sense to define) inner polyhedron T (could be boxes, zones, octagons, linear-templates, etc. depending on A 2 ), which obeys the template pattern matrix T .

[-2, 3] × [-2, 2], [-1, 1] × [0, 1]). Figure
Let P = {x ∈ R p |Ax ≤ b} be a non-empty polyhedron, where A is a known m×p matrix, b a known vector of R m , and x a vector of R p . The inner polyhedron T is expressed in a similar manner: T = {x ∈ R p |T x ≤ c}, where T is a known n × p matrix, and c and x are unknown vectors within R n and R p , respectively. The inclusion T ⊆ P holds if and only if

∃c ∈ R n , such that T is consistent, and ∀x ∈ R p : T x ≤ c =⇒ Ax ≤ b .
The consistency of T (that is the system admits a solution in R p ) discards the trivial (and unwanted) cases where the polyhedron T is empty. For the non trivial cases, the existence of the vector c and the characterization of the set of its possible values are given by Proposition 3. Proposition 3. Let C be the set of c such that T is consistent. There exists a vector c ∈ C such that T ⊆ P if and only if there exists an n × m matrix Λ, such that λ i,j , the elements of the matrix Λ, are non-negative and ΛT = A. For a given possible Λ, the set c Λ ⊆ C is characterized by

{c ∈ R n | Λc ≤ b} .
Proof. Let x denote a vector of R p , and b denote a known vector of R m . Let A and T be two known matrices with p columns and m and n lines, respectively. Suppose that c is such that T is consistent. Therefore, we can assume that

t i , x ≤ c i , 1 ≤ i ≤ n,
where t i denotes the ith line of the matrix T , is consistent. For a fixed j, 1 ≤ j ≤ m, the inequality a j , x ≤ b j , is then a consequence of the system T x ≤ c if and only if there exist non-negative real numbers λ i,j , 1

≤ i ≤ n, such that n i=1 λ i,j t i = a j and n i=1 λ i,j c i ≤ b j .
The previous claim of the existence of the non-negative λ i,j is a generalization of the classical Farkas' Lemma (see for instance [START_REF] Rockafellar | Convex Analysis[END_REF]Section 22,Theorem 22.3] for a detailed proof). The matrix Λ is then constructed column by column using the elements λ i,j , 1 ≤ i ≤ n for the jth column. Of course, by construction, such a Λ has non-negative elements, and satisfies ΛT = A, and Λc ≤ b.

On the other hand, if such a matrix exists, and the set {c ∈ R n | Λc ≤ b} is not empty, we have by the fact that Λ has non-negative elements

T x ≤ c =⇒ ΛT x ≤ Λc .
Therefore, ΛT = A and Λc ≤ b, gives Ax ≤ b.

On the Consistency of T x ≤ c. It not obvious in general, given a matrix T , to characterize the set of c such that T is consistent. However, given a vector c, we can efficiently check whether the system is consistent or not using its dual form and a LP solver. Indeed, the system T x ≤ c is inconsistent if and only if there exists a non negative vector λ ∈ R n such that T t λ = 0 and λ, c < 0, where T t denotes the transpose of T . Therefore, given a vector c, if the objective value of the following problem min λ, c s.t. T t λ = 0 .

(

) 2
is non negative, the system is consistent. Observe that, for simple patterns such as boxes, the characterization of the set of c that makes the system consistent is immediate.

Computing Λ. The matrix Λ is built column by column. Let us denote by λ -,j ∈ R n the jth column of Λ, by a j ∈ R p , 1 ≤ j ≤ m, the jth line of A, by b j ∈ R the jth component of b, and by t i ∈ R p , 1 ≤ i ≤ m, the ith line of T . The vector λ -,j satisfies n i=1 λ i,j t i = a j . To each feasible λ -,j corresponds a pattern

P λ-,j def = {x ∈ R p | λi,j >0 t i , x ≤ 0},
which is included in the affine subspace P j def = {x ∈ R p | a j , x ≤ 0}. The maximal pattern (with respect to set inclusion) corresponds to λ defined as the solution of the following linear program.

min n i=1 λ i,j t i s.t. n i=1 λ i,j t i = a j ∀0 ≤ i ≤ n, λ i,j ≥ 0 . (3)
Therefore, computing Λ needs solving p instances of the LP (3).

Computing c. We have already established (Proposition 3) that the vector c verifies Λc ≤ b. Since Λ is known, any feasible c (that is such that Λc ≤ b) that makes the system T x ≤ c consistent (the objective value of the LP (2) is non negative) gives an under-approximation of P that respects our initial template T . Of course, it is immediate to see that the set of c that lies on the boundaries of the feasible region (that is by making Λc = b) gives, in general, a "better" under-approximation than the strict feasible solutions since the saturation makes some of the facets of the inner pattern (T ) included in those of the underapproximated polyhedron P. Moreover, in some cases, the saturation gives a unique consistent solution for c. For instance, when we under-approximate a shape P which respects already the pattern T , c is uniquely determined and gives actually b using our technique. In other words, under-approximating an octagon (for instance) with an octagonal pattern gives exactly the first octagon.

Implementation

We have implemented donut domains on the top of Apron library [START_REF] Jeannet | Apron: A library of numerical abstract domains for static analysis[END_REF]. The domains A 1 and A 2 are parameters of the analysis and can be specified by the user among already existing Apron domains. The current version uses an enhanced implementation of the set-theoretic operators, mainly based on already existing routines of the underlying abstract domains, as described earlier and relies on ∪fb and ∩fb as fallback operators. This very simple approach allows to build the donut domain without an additional effort on the top of already existing domains. As Table 1 shows, our approach permits to catch almost all division-byzero false positives that classical domains (even non-convex) fail to prove. The WCfS column indicates the weakest condition that we need to infer to prove the [START_REF] Adjé | Coupling policy iteration with semi-definite relaxation to compute accurate numerical invariants in static analysis[END_REF]. Here the safety property can not be proved with this hole and relies on an earlier (disjoint) hole created by a former test, namely usemax ∈ {0}. We could also choose systematically (as a heuristic) the hole that contains "zero", which is sufficient here to discard the remaining false alarm. Such a property-driven hole behaviour would be an interesting direction for future research. The proof of the motivating example is really challenging as it requires to handle both the hole that comes from the full-zero-test, together with strict inequalities tests and the over-approximation that comes from the join operation. Our technique that consists of relaxing the hole in a convex way before using the fallback operator works here and is able to prove that in both branches the division is safe. In goc example, we can see one interesting ability of donuts domain: when we compute a convex join of two non-overlapping objects, the hole in between is directly captured which permits a better precision. Finally, example x2 needs a precise interpretation of strict inequalities.

Under-Approximation. We have implemented our technique of Section 3 using GLPK [START_REF] Makhorin | The GNU Linear Programming Kit (GLPK)[END_REF] solver. Some experiments, obtained for randomly generated polyhedra with octagonal template, are shown in Figure 7. Although all shown polyhedra are bounded, our technique works perfectly well for unbounded shapes. The rate of volume, volT volP , is used as a metric for the quality of the underapproximation (shown near each pattern in Figure 7). All obtained octagons are maximal with respect to set inclusion. It is not clear which choice among many (see the left graph), is the best. Indeed, such a choice depends on the future computations and the properties one would like to prove. The donut domains can be viewed as an effort to make some Boolean structure in the underlying concrete space visible at the level of abstract domains as a "setminus" operator. This allows optimization of the related abstract operators (such as meet and join) to take full advantage of its semantics in terms of excluded states. While the powerset domains allow an arbitrary Boolean combination, and automata-based abstractions [START_REF] Dams | Automata as abstractions[END_REF] using may-and must-transitions offer full generality in verification, this comes at significant cost. In practice, the full expressiveness may not be needed. We exploit the set-minus operator, which is quite versatile in capturing many problems of interest -division by zero, instability regions in numeric computations, sets excluded by contracts in a modular setting, etc.

In the future, we wish to expand the experiments performed using donut domains. Furthermore, other non-convexity issues may be addressed by trying to combine the work on LDDs with insights gained here to allow handling many holes in an efficient manner.

Fig. 2 .

 2 Fig. 2. The concretization of a typical non-convex abstract object

Example 1 .

 1 Consider a one-dimensional donut domain where A 1 and A 2 are Intervals domain. Assume we are interested in computing ([0, 3], [1, 2]) ∪ ([1, 6], [2, 5]) . The above join yields the following union of four intervals: [0, 1) ∪ (2, 3] ∪ [1, 2) ∪ (5, 6], which can be combined without loss of precision into [0, 2) ∪ (2, 3] ∪ (5, 6], or equivalently [0, 6] \ ([2] ∪ (3, 5]) .

Fig. 3 .

 3 Fig. 3. Relaxing the hole (0, 0) (red circle in the left hand side figure) to x ≥ 0

  0], the intersection gives the box [-1, 1] × [-1, 0], and the extension we compute gives the box [-2, 2] × [-1, 0]. If, however, the holes are disjoint, we randomly pick up one of them.

Fig. 4 .

 4 Fig. 4. Illustrating the join and meet operators using interval component domains. The donut holes are highlighted using dashed lines. (a) Two initial abstract objects. (b) The concrete union of the objects. (c) The abstract object representing ∪ 1\2 . (d) The concrete intersection of the objects. (e) The abstract object representing ∩ 1\2 .

Fig. 5 .Fig. 6 .

 56 Fig. 5. Evaluation of a linear expression x2 ← x1 + x2 1\2

Fig. 7 .

 7 Fig. 7. Under-approximation of randomly generated polyhedra with octagons

Table 1 .

 1 Division-by-zero analysis resultssafety of the program. Whenever the negation of this condition is verified by (included in) the donut hole, the program is proved to be safe. The third column shows the inferred donut holes when using a non-relational domain (boxes) to encode holes. The analyzed examples 4 use mainly the absolute value function to avoid the division by zero (widely used technique).The motiv example is the motivating example with its two branches. The gpc code is extracted from the Generic Polygon Clipper project. The examples xcor, goc and x2 are extracted from a geometric object contact detection library. Observe that the use of boxes is sufficient to eliminate almost all false alarms here. In the last example, among the two possible holes, namely usemax ∈ [1, 10] and usemax ∈ {0}, we choose by default the one created immediately after the test (usemax > 10 or usemax <

www.nec-labs.com/research/system/systems SAV-website/benchmarks.php. The C files are the real source code, while the SPL files extracts the hard piece of code that leads to false alarms, and with which we feed our proof of concept implementation.

Acknowledgment. The authors would like to thank Sriram Sankaranarayanan, Enea Zaffanella and anonymous reviewers for their valuable comments on an earlier draft of this work.