Jean-Baptiste Jeannin
email: jeannin@cs.cmu.edu

Khalil Ghorbal
email: kghorbal@cs.cmu.edu

Yanni Kouskoulas
email: yanni.kouskoulas@jhuapl.edu

Ryan Gardner
email: ryan.gardner@jhuapl.edu

Aurora Schmidt
email: aurora.schmidt@jhuapl.edu

Erik Zawadzki

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

Formal Verification of ACAS X, an Industrial Airborne Collision Avoidance System

Formal verification of industrial systems is very challenging, due to reasons ranging from scalability issues to communication difficulties with engineering-focused teams. More importantly, industrial systems are rarely designed for verification, but rather for operational needs. In this paper we present an overview of our experience using hybrid systems theorem proving to formally verify ACAS X, an airborne collision avoidance system for airliners scheduled to be operational around 2020. The methods and proof techniques presented here are an overview of the work already presented in [8], while the evaluation of ACAS X has been significantly expanded and updated to the most recent version of the system, run 13. The effort presented in this paper is an integral part of the ACAS X development and was performed in tight collaboration with the ACAS X development team.

INTRODUCTION

In recent years, aircraft collision avoidance has been a case study of choice for formal verification of cyber-physical systems, in particular thanks to its very clear specification -absence of collision [START_REF] Dowek | Provably safe coordinated strategy for distributed conflict resolution[END_REF][START_REF] Galdino | Formal verification of an optimal air traffic conflict resolution and recovery algorithm[END_REF][START_REF] Ghorbal | Hybrid theorem proving of aerospace systems: Applications and challenges[END_REF][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][START_REF] Tomlin | Conflict resolution for air traffic management: A study in multiagent hybrid systems[END_REF]. However few of those studies were done on industrial systems -i.e., systems which will ultimately equip commercial aircraft -or as an integral part of the system's development, with potential impact on the system's final design.

On the industrial side, two main systems exist, both developed under the lead of the Federal Aviation Administration (FAA): the Traffic Collision Avoidance System (TCAS) [START_REF]Introduction to TCAS II[END_REF], which is currently mandated on all large passenger aircraft and whose design started in the late 1970s; and the Next-Generation Airborne Collision Avoidance System (ACAS X) [START_REF]Introduction to TCAS II[END_REF][START_REF] Holland | Optimizing the next generation collision avoidance system for safe, suitable, and acceptable operational performance[END_REF][START_REF] Kochenderfer | Next generation airborne collision avoidance system[END_REF], a complete redesign of the TCAS system initiated by the FAA to accommodate recent challenges such as more crowded airspace and the arrival of Unmanned Aerial Vehicles. This paper describes our work using hybrid systems modeling and theorem proving to formally assess the safety of ACAS X. This work is part of ACAS X's development and has been done in collaboration with the ACAS X development team. This paper focuses on what is called run 13 of the XA version of ACAS X (simply denoted ACAS X throughout this paper) -the most recent version used for normal operations of manned aircraft.

A typical collision avoidance encounter involves two aircraft: the equipped ownship aircraft actively trying to avoid colliding with a non-equipped intruder aircraft. Both TCAS and ACAS X avoid such collisions by giving exclusively vertical advisories to the pilot of the ownship -such as "Climb at a rate of 1,500 ft/min" (CL1500) or "Do not descend" (DND). If both aircraft are equipped with a collision avoidance system, both pilots can get coordinated advisories. A collision -or Near Mid-Air Collision (NMAC) -is formally defined as the two aircraft being within rp = 500 ft horizontally and hp = 100 ft vertically -thus forming a safety puck around the ownship where the intruder should not enter to guarantee the safety of both aircraft. While the design of TCAS was based on geometric considerations on paths taken by aircraft, the design of ACAS X is radically different. ACAS X is designed around a Markov Decision Process (MDP) estimating the probabilities that two aircraft may climb, descend or turn, in the presence or absence of an advisory [START_REF] Kochenderfer | Robust airborne collision avoidance through dynamic programming[END_REF]. Based on those probabilities, the MDP is optimized to minimize the probability of an NMAC, while also minimizing (with lesser priority) advisories, disruptive to the pilot. This optimization results in a score table giving advisories for different possible configurations. The table is then queried and interpolated in flight; the result is sometimes corrected with online scores to take into account factors not considered by the MDP such as proximity to the ground or coordination with other ACAS X-equipped aircraft, before being communicated to the pilot. The score table combines two sub-tables interpolated in succession, one with 698,819 interpolation points and the other with 38,053,125 interpolation points, resulting in a large number of combinatorial cases.

Given this design, what does it mean to formally and effectively verify ACAS X? Observe that absence of NMAC cannot be guaran-teed in all cases: if the two aircraft are flying head-on at the same altitude and are too close to maneuver, no flyable advisory can avoid an NMAC. Moreover, directly verifying each point in the ACAS X score table is infeasible because of the size of the table. Instead, we seek to formally quantify the behavior of ACAS X. Our approach focuses on the score table and is composed of three major steps. First, based on a model of the typical behavior of the aircraft, we symbolically characterize safe regions, geometric configurations of the ownship and intruder's positions and speeds such that a given advisory, if followed correctly, will not result in an NMAC. Second, we model both aircraft as a hybrid system combining differential equations representing continuous motion of the aircraft, and discrete actions modeling advisories -issued according to the safe zones -and compliance of the pilot. Using this model, we formally prove that, if the two aircraft are in a safe region, following the advisory associated with that safe region will not result in an NMAC. Third and finally, we exhaustively compare the safe regions to the advisories given by the score table, thus transferring our formal argument of safety to the operational ACAS X. We identify parts of the table that are provably safe, but also counterexamples of situations where following the table could lead to an NMAC. We feed back those results to the ACAS X development team. The general formal verification approach used in this work -modeling, proving and comparing -is not limited to ACAS X, but can potentially be reused for other collision avoidance systems, or even other systems whose design is based on optimization or machine learning.

Related Work. Our approach is different from previous complementary work in various ways. We analyze an independent industrial system instead of hypothetical safe-by-design systems as in [START_REF] Dowek | Provably safe coordinated strategy for distributed conflict resolution[END_REF][START_REF] Galdino | Formal verification of an optimal air traffic conflict resolution and recovery algorithm[END_REF][START_REF] Ghorbal | Hybrid theorem proving of aerospace systems: Applications and challenges[END_REF][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][START_REF] Tomlin | Conflict resolution for air traffic management: A study in multiagent hybrid systems[END_REF]. Our hybrid model uses realistic, continuous dynamics, and does not discretize the trajectories of aircraft as in [START_REF] Dowek | Provably safe coordinated strategy for distributed conflict resolution[END_REF][START_REF] Von Essen | Analyzing the next generation airborne collision avoidance system[END_REF][START_REF] Galdino | Formal verification of an optimal air traffic conflict resolution and recovery algorithm[END_REF]. We provide mechanized proofs of correctness of our model instead of simulations or stress testing approaches as in [START_REF] Chludzinski | Evaluation of TCAS II version 7.1 using the FAA fast-time encounter generator model[END_REF][START_REF] Holland | Optimizing the next generation collision avoidance system for safe, suitable, and acceptable operational performance[END_REF][START_REF] Kochenderfer | Correlated encounter model for cooperative aircraft in the national airspace system version 1[END_REF][START_REF] Lygeros | On the formal verification of the TCAS conflict resolution algorithms[END_REF][START_REF] Tomlin | Conflict resolution for air traffic management: A study in multiagent hybrid systems[END_REF]. Finally, unlike [START_REF] Von Essen | Analyzing the next generation airborne collision avoidance system[END_REF][START_REF] Kochenderfer | Robust airborne collision avoidance through dynamic programming[END_REF], our hybrid model and proofs are independent from the different versions of the score table used by ACAS X, thus allowing us to provide universal safe regions that can be reused for ACAS X or even new systems. A more detailed discussion of related work can be found in [START_REF] Jeannin | A formally verified hybrid system for the next-generation airborne collision avoidance system[END_REF].

This paper gives a brief overview of the modeling approach in Sect. 2 and focuses on the comparison and its results in Sect. 3. Our ongoing work and future avenues are reported in Sect. 4.

OVERVIEW OF THE APPROACH

We present and motivate the hybrid model and construction of the safe regions. We sketch the main steps of our formal proof of safety. For more details, we refer the reader to [START_REF] Jeannin | A formally verified hybrid system for the next-generation airborne collision avoidance system[END_REF] and its companion technical report [START_REF] Jeannin | A formally verified hybrid system for the next-generation airborne collision avoidance system[END_REF].

Modelling ACAS X

Let us consider a 3-dimensional encounter involving an ownship O and an intruder I. For the purposes of this study, we make a number of assumptions on the capabilities and behavior of both aircraft during an encounter. We assume that the ownship has the capability to perfectly estimate the geometrical configuration of the encounter, i.e., the relative position of the intruder, as well as the speed and direction of both aircraft; in practice such information is provided by different sensors as well as communication data exchanged between aircraft. In particular, we discard all sources of noise and uncertainty related to the state estimation. We further assume that the intruder has a constant velocity vector, both horizontal and vertical, with no turn, acceleration or change in its vertical speed, thus describing a straight-line uniform trajectory. As for the ownship, we assume that it has a straight-line uniform horizontal trajectory, with no turn or horizontal acceleration. We, however, allow nondeterministic changes in the ownship vertical acceleration that are used to come into compliance with and follow vertical advisories. These assumptions limit the generality of the result, but the proof presented in this paper is still powerful because most planes, most of the time, fly in straight lines; therefore we expect ACAS X to be able to handle this sort of steady-state pilot model. We make these assumptions throughout the technical development of the paper, including in Sect. 3.

ACAS X can issue one of 16 vertical advisories, summarized in Table 1. For example, if provided with the (corrective) advisory "Climb 1,500" (CL1500), the pilot is instructed to climb at a rate of at least 1, 500 ft/min. Similarly, the (preventive) advisory "Do Not Climb" (DNC) instructs the pilot to not climb, while the advisory "Maintain Climb" (MCL) instructs the pilot to maintain their current climb rate. Our model, based on [START_REF] Kochenderfer | Robust airborne collision avoidance through dynamic programming[END_REF], supposes that the pilot complies with the advisory with acceleration at least ar = g/4 starting after a delay modeling reaction time of at most dp = 5 seconds; then the pilot continues following the advisory forever -or until a new advisory is issued. Four advisories (SDES1500, SCL1500, SDES2500 and SCL2500) can only be used to strengthen a previous advisory; as such the compliance of the pilot is typically better, and our model supposes a compliance with an acceleration of ar = g/3 and within a delay of dp = 3 seconds. ACAS X can also issue the advisory "Clear of Conflict" (COC), meaning that no specific action is required from the pilot. All advisories, except COC, request the pilot to stay either above or below a given vertical velocity, as shown in Table 1. We model such advisories with two variables: the target vertical velocity given to the pilot ḣf and a binary variable w, taking value -1 if the pilot is asked to stay below that velocity, or +1 if she is asked to stay above.

Reduction to a 2-dimensional Encounter

Under our assumptions, the 3-dimensional encounter can be reduced to a 2-dimensional encounter without loss of generality. This observation makes the formal model much easier to reason about. This section gives a geometric intuition of the reduction. A formal proof, based on a similar argument, is presented in [START_REF] Jeannin | A formally verified hybrid system for the next-generation airborne collision avoidance system[END_REF]. For this reduction let us put ourselves in the reference frame of the ownship. A 3-dimensional encounter is represented in Fig. 1 -top view at the top and side view at the bottom. Let r be the distance between the aircraft, rv the relative velocity of the intruder with respect to the ownship, and θv the angle formed by (OI) and rv.

DNC -∞ 0 g/4 5 -1 0 DND 0 +∞ g/4 5 +1 0 MDES -∞ current g/4 5 -1 current MCL current +∞ g/4 5 +1 current DES1500 -∞ -1500 g/4 5 -1 -1500 CL1500 +1500 +∞ g/4 5 +1 +1500 SDES1500 -∞ -1500 g/3 3 -1 -1500 SCL1500 +1500 +∞ g/3 3 +1 +1500 SDES2500 -∞ -2500 g/3 3 -1 -2500 SCL2500 +2500 +∞ g/3 3 +1 +2500 COC -∞ +∞ Not applicable
The green line (DC), along the vector rv represents the horizontal trajectory of the intruder I relative to the ownship O. The gray disk surrounding the ownship represents the projected safety puck.

Recall that to avoid an NMAC, the trajectory of the intruder must not intersect the safety puck. If r| sin θv| > rp, then the green projected trajectory of the intruder does not intersect the projected safety puck and there is no possible NMAC. Otherwise, let us look at a side view of the vertical plane along the trajectory (DC), represented at the bottom of Fig. 1. In this plane, the intruder avoids an NMAC if and only if it does not enter a virtual puck centered on a virtual ownship M middle of [DC], with a modified puck radius √ (r 2 -r 2 p sin 2 θv) and a modified distance to M, r| cos θv|. This reduces the 3D encounter to a virtual planar head-on encounter.

Construction of the Safe Regions

Let us now change the reference frame to put ourselves in the reference frame fixed to the intruder and centered on the original position of the ownship. Let h be the relative height of the intruder with respect to the ownship, ḣ0 the relative vertical velocity of the ownship with respect to the intruder, and a the relative vertical acceleration of the ownship with respect to the intruder.

Let us consider an ownship initially descending and receiving a CL1500 advisory, meaning that the pilot should start climbing at 1,500 ft/min. Let us further assume that the pilot reacts immediately. The nominal trajectory of the ownship is illustrated in red on Fig. 2: it consists of a parabola at acceleration g/4 until reaching the advised 1,500 ft/min, followed by a straight line climbing at 1,500 ft/min. If the pilot follows the advisory, the ownship will always be at or above this nominal trajectory. This nominal trajectory gives rise to a first construction of what we call the safe region: if the intruder is in the green area to the right, left or below all pucks centered on a point of the red trajectory, then the intruder is in a safe region, and no NMAC can occur as long as the straight line encounter assumptions are satisfied. Based on the nominal trajectory, we can mathematically express a symbolic implicit formulation of the safe region in full generality, given in Fig. 3.

The implicit formulation involves three universal quantifiers over t, rt and ht. Eliminating those quantifiers is crucial when it comes to comparing the actual advisories of ACAS X. Quantifier elimination algorithms could be used but often lead to formulations that are very large and hard to interpret. Instead, we go back to Fig. 2 and manually construct the boundaries of the safe region. They consist of a vertical line on the left side created by the left side of the puck, followed by a half parabola created by the bottom left-hand corner of the puck, then a horizontal segment of length 2rp created by the bottom of the puck, and finally another half parabola and a straight-line at 1,500 ft/min, both created by the bottom right-hand corner of the puck. This observation leads to an explicit formulation of the safe regions which can be expressed mathematically in full generality and is given in Fig. 3.

Both formulations are equivalent, but each comes with an advantage: the implicit formulation makes the formal proof of the hybrid model easier, while the explicit formulation enables an easy comparison with the ACAS X table. We formally prove in the KeYmaera theorem prover [START_REF] Platzer | KeYmaera: A hybrid theorem prover for hybrid systems[END_REF] that the implicit formulation of the safe region C impl (r, h, ḣ0) and its explicit formulation C expl (r, h, ḣ0) are equivalent. This formal proof, along with all the KeYmaera models and proofs of this paper are available online at http://www.ls.cs.cmu.edu/pub/acasx.zip.

Proof of Safety

In order to formally prove the correctness of the safe regions, we build in Eq. (1) a hybrid model described by the differential dynamic logic [START_REF] Platzer | Differential dynamic logic for hybrid systems[END_REF]. We then formally prove, using the hybrid systems theorem prover KeYmaera [START_REF] Platzer | KeYmaera: A hybrid theorem prover for hybrid systems[END_REF], that no NMAC can occur, provided all advisories satisfy the condition C impl (r, h, ḣ0) (or equivalently C expl (r, h, ḣ0)). Implicit formulation

1 rp ≥ 0 ∧ hp > 0 ∧ rv ≥ 0 ∧ ar > 0∧ 2 (w = -1 ∨ w = 1) ∧ C impl (r, h, ḣ0) →
7 {r = -rv, h = -ḣ0, ḣ 0 = a & w ḣ0 ≥ w ḣf ∨ wa ≥ ar} 8 ) * ] (|r| > rp ∨ |h| > hp) (1)
A(t, ht, ḣ0) ≡ 0 ≤ t < max(0, w( ḣf -ḣ0)) ar ∧ ht = war 2 t 2 + ḣ0t ∨ t ≥ max(0, w( ḣf -ḣ0)) ar ∧ ht = ḣf t - w max(0, w( ḣf -ḣ0)) 2 2ar C impl (r, h, ḣ0) ≡ ∀t.∀rt.∀ht. rt = rvt ∧ A(t, ht, ḣ0) → (|r -rt| > rp ∨ w(h -ht) < -hp) Explicit formulation case1(r, ḣ0) ≡ -rp ≤ r < -rp - rv min(0, w ḣ0) ar case2(r, ḣ0) ≡ -rp - rv min(0, w ḣ0) ar ≤ r ≤ rp - rv min(0, w ḣ0) ar bound1(r, h, ḣ0) ≡ wr 2 v h < ar 2 (r + rp) 2 + wrv ḣ0(r + rp) -r 2 v hp bound2(r, h, ḣ0) ≡ wh < - min(0, w ḣ0) 2 2ar -hp case3(r, ḣ0) ≡ rp - rv min(0, w ḣ0) ar < r ≤ rp + rv max(0, w( ḣf -ḣ0)) ar case4(r, ḣ0) ≡ rp + rv max(0, w( ḣf -ḣ0)) ar < r bound3(r, h, ḣ0) ≡ wr 2 v h < ar 2 (r -rp) 2 + wrv ḣ0(r -rp) -r 2 v hp bound4(r, h, ḣ0) ≡ (rv = 0) ∨ wrvh < w ḣf (r -rp) - rv max(0, w( ḣf -ḣ0)) 2 2ar -rvhp case5(r, ḣ0) ≡ -rp ≤ r < -rp + rv max(0, w( ḣf -ḣ0)) ar case6(r, ḣ0) ≡ -rp + rv max(0, w( ḣf -ḣ0)) ar ≤ r bound5(r, h, ḣ0) ≡ wr 2 v h < ar 2 (r + rp) 2 + wrv ḣ0(r + rp) -r 2 v hp bound6(r, h, ḣ0) ≡ (rv = 0 ∧ r > rp) ∨ wrvh < w ḣf (r + rp) - rv max(0, w( ḣf -ḣ0)) 2 2ar -rvhp C expl (r, h, ḣ0) ≡ w ḣf ≥ 0 → 4 i=1 (casei(r, ḣ0) → boundi(r, h, ḣ0)) ∧ w ḣf < 0 → 6 i=5 (casei(r, ḣ0) → boundi(r, h, ḣ0))
Figure 3: Implicit and explicit formulations of the safe region for an immediate response [START_REF] Jeannin | A formally verified hybrid system for the next-generation airborne collision avoidance system[END_REF] Eq. ( 1) is of the form p → [α]q, which says all executions of program α starting in a state satisfying logical formula p end up in a state satisfying q. It is akin to the Hoare triple {p}α{q} with precondition p and postcondition q. The precondition in Eq. ( 1) imposes constraints on several constants, as well as the formula C impl (r, h, ḣ0) that forces the intruder to be in a safe region for an initial advisory (w, ḣf ). We cannot guarantee safety if the intruder starts in an unsafe region. The postcondition encodes absence of NMAC. Lines 3 to 5 express the action of ACAS X. The nondeterministic choice operator ∪ expresses that the system can either continue with the same advisory by doing nothing -testing the trivial condition ?true -which ensures it always has a valid choice and cannot get stuck. Otherwise it can choose a new advisory (w, ḣf ) that passes the safety condition C impl (r, h, ḣ0) -the variable "advisory" will be the next message to the pilot. Lines 6 and 7 express the action of the ownship, first nondeterministically choosing an arbitrary acceleration (a := * ) then following the continuous dynamics. The evolution of the variables r, h and ḣ0 is expressed by a differential equation, and requires (using the operator &) that the ownship evolves towards its target vertical velocity ḣf at acceleration ar (condition wa ≥ ar), unless it has already reached vertical velocity ḣf (condition w ḣ0 ≥ w ḣf ). Finally, the star * on line 8 indicates that the program can be repeated any number of times, allowing the system to go through several safe advisories.

Delayed Pilot Response and "Clear of Conflict" Advisory

A similar process can be followed to add a delay parameter that models the reaction time of the pilot before she starts following the advisory. During that delay we consider a worst case scenario where the pilot can climb or descend at an acceleration of g/3. Details of the approach, including a formally proved model of the system and equations of the safe regions, can be found in [START_REF] Jeannin | A formally verified hybrid system for the next-generation airborne collision avoidance system[END_REF].

This approach can be adapted to model a "Clear of Conflict" (COC) advisory. When a COC advisory is issued, our model allows the pilot to climb or descend with an acceleration up to g/3. A COC advisory can be considered safe if and only if it can be followed by another safe advisory 1 second later -since ACAS X can issue advisories every second. Thus, COC is safe if and only if there exists an advisory that is safe for a delay augmented by 1 second.

COC has the same representation as other advisories in the model -a period of delay followed by compliance with an advisory -but has a different interpretation. The COC advisory itself only advises the pilot that the g/3 action for the delay period in the model is safe.

The vertical velocity limit attached to COC represents a future safe advisory that may be issued after the COC to ensure safety after the COC period is over. As such, the value of its vertical velocity limit is variable, and depends on what could be safely issued by the system after a period of pilot delay and free vertical acceleration.

Limitations

Our approach and model have some limitations. Our model of a straight-line uniform trajectory for the intruder disallows any turn or reaction on its part; in particular it does not handle an intruder equipped with a collision avoidance system. Similarly, our model does not allow any turn on the part of the ownship. As a result, our safe regions do not guarantee NMAC avoidance for different, nonstraight-line pilot behavior, and not satisfying our safety conditions does not necessarily lead to an NMAC. Furthermore, our model permits but does not enforce follow-on advisories. Our analysis thus conservatively rejects advisories that are not safe for all future, even those that ACAS X would ultimately strengthen or reverse. But follow-on advisories can change the outcome of the encounter if there is enough time for the pilot to react to the modified advisory.

In addition of the advisories presented in Table 1, ACAS X can issue the recently-added Multi-Threat Level-Off (MTLO) advisory that we do not handle yet in our model or safe regions. Finally, our comparison with the ACAS X system focuses on its score table exclusively, and ignores online scores that account for example for ground proximity or coordination.

COMPARISON

For each configuration of the aircraft, the safety conditions presented in Fig. 3 give us a means of evaluating the safety of the ACAS X table, by performing a comparison between the provably safe set of advisories allowed by the safety conditions, and the advisory recommended by the ACAS X table. Any mismatch gives us a counterexample of a potentially unsafe behavior. This section presents the method and results of the comparison, then examines a few automatically generated counterexamples.

Method

The ACAS X score table used to determine advisories is stored as a (high-dimensional) grid of points, often called cutpoints. When the ACAS X table is queried online, it performs a multilinear interpolation to combine the table data from the surrounding cutpoints. Thus, in a sense, all of the information in the table is determined by its cutpoints. As a first step in analyzing the safety of the ACAS X table, we compare the advisory given by the system for each cutpoint with the set of advisories that are proved safe in our alternate model. Throughout this section we call a cutpoint safe if its corresponding ACAS X advisory respects the safe regions.

The comparison we perform covers all of the cutpoints in the table exhaustively (except for the recently-added advisory MTLO).

For the run 13 of ACAS X, cutpoints are defined by: relative horizontal motion (angle and magnitude, with 187 and 37 discrete values); horizontal separation (101 discrete distances); vertical velocities (25 discrete vertical velocities each for ownship and intruder); relative vertical separation (45 discrete distances); the current state of the advisory, and whether the system estimates it is being followed (33 distinct states). Some continuous dimensions have equally spaced cutpoints (e.g. θv is equally spaced in radians), but others have variable spacing, with points closer in areas that may be more dangerous (e.g. vertical separation is discretized in steps of 50 ft between -400 ft and 400 ft, but the discretization step then becomes increasingly larger -100 ft, 200 ft, 400 ft, and finally 4,000 ft at the last step -as the separation increases). All in all, there are 648,591,384,375 (6.49e11) unique cutpoints to analyze. Because the table may have arbitrarily different entries at each cutpoint, we cannot rely on any sort of symmetry arguments, but have to explore the full state space exhaustively. Fortunately however, we do not have to predict all of the infinitely many possible future trajectories from a cutpoint to check whether an advisory is safe, but we can exploit the safe regions identified in Sect. 2 to reduce this to a simple arithmetic check of the formula in Fig. 3.

Improving upon our previous paper [START_REF] Jeannin | A formally verified hybrid system for the next-generation airborne collision avoidance system[END_REF], we present new results and counterexamples featuring a delayed pilot response within the limits assumed in Sect. 2, and counterexamples where the system issued COC. Because the safe regions are symbolic formulas, the same comparison can be repeated with any other instantiation of the parameters such as the worst case vertical acceleration (assumed to be g/3) or the pilot delay (3 seconds). The proofs for the safe regions identify them as safe for any instantiation within the model, which shows a clear advantage of the formal symbolic approach whenever achievable.

When the system issues a COC, the comparison computes a twosided envelope, since the system isn't asking the pilot to go in any particular direction. Each side of the envelope represents the most extreme acceleration in that direction during the delay (g/3), followed by the most extreme advisory after 1 second in the other direction (either DES1500 or CL1500). The upper and lower bounds of the envelope eventually cross to create a closed area, whose inside is considered unsafe. The area outside this curve is safe for any intruder, and guarantees we can avoid collision if we issue COC now, by issuing the appropriate follow-on advisory the next second.

Results

We ran two different analyses, one comparison assuming that the pilot reacts immediately (no delay), and another one assuming that the pilot reacts after a 3-second delay. We present the results in Table 2. For the 3-second-delay comparison, 97.7% of the cutpoints we examined were either safe, or had no alternate advisory for which they could be made safe. Only about 2.3% of the cutpoints had alternate advisories that could be issued to improve the table's behavior in straight-line encounters. For the no-delay comparison, we found similar statistics of 97.6% and 2.4%, respectively.

Counterexamples

We discuss four counterexamples found by the comparison, using either no delay for pilot response, or a 3-second delay prior to compliance. The counterexamples are visualized by plotting the predicted ownship and intruder paths in altitude versus time, afer horizontal reduction if necessary (Figures 4 and5). When visualizing the counterexample, a particular acceleration and vertical rate target is chosen during the delay period that when combined with later compliance to an advisory, results in an NMAC.

The two panels of Fig. 4 show examples where the ACAS X advisory appears to increase the risk of an NMAC over the nominal risk had the pilot continued on the same trajectory. Panel (a) shows a counterexample found using no delay in pilot response. In this case, the ownship was originally planning to cross in altitude with the intruder. The table issues a DNC advisory, causing the ownship to reduce its upward velocity to zero at a rate of g/4. This results in a collision course with the intruder and an NMAC approximately 20 seconds later. After inspecting the ACAS X table, it was determined that a major cause for this behavior was inaccurate horizontal timing estimation in the table. As a result of this counterexample, improvements are currently being made to the system's timing estimation components, thus demonstrating the benefit of applying formal verification techniques early on in the development.

Panel (b) of Fig. 4 shows a counterexample where the advisory issued induces a risk for a pilot that responds after a delay. In the example plotted, during the delay of 3 seconds, the ownship pilot increases climbing with acceleration g/3 to reach 5,000 ft/min, before complying with the DNC advisory with acceleration g/4. This possible behavior puts the ownship on a near collision course with the intruder. Based on the anticipated miss distance of the nominal course, advisories CL1500 or MCL would have lowered the risk here. These counterexamples where ACAS X advisories appear to induce a problem are collected and communicated to the system designers for studying and improving the system. Fig. 5 shows counterexamples where the ACAS X table does not issue an advisory (i.e., issues COC), but the scenario would be significantly safer with an advisory for straight-line encounters. Panel (a) shows an identified counterexample where, due to the very close horizontal range of only 150 feet, the intruder is already within horizontal conflict range. However, the intruder begins 1,000 feet above the ownship and is descending. Without an advisory, the ownship could decide to climb for 1 second before an advisory can be issued by the system. Even if the most extreme available advisory (DES1500) were issued 1 second later, the intruder would enter the back edge of the puck. If, however, a descend advisory were issued in the first place and the pilot responds with no delay, then the observed NMAC be avoided.

In panel (b) of Fig. 5, due to the high horizontal closure rate the time to horizontal conflict is quite short (about 7 seconds). By not issuing an advisory at this state, the pilot may reduce climb to only 8, 000 ft/min during the 1 second where no advisory is issued and keep the same vertical velocity during the assumed delay of 3 seconds after a corrective advisory, such as SDES1500 (descend at rate 1,500 ft/min) is issued. By the time the pilot starts complying with the advisory it is already too late to avoid NMAC. Due to the high vertical rates of this example, the descend advisories are the only type of advisories available to ACAS X, because the system is designed to not increase the climbing rate beyond 2,500 ft/min. These two counterexamples show identified states where not issuing an advisory may actually augment the risk of NMAC.

The comparison we are performing generates billions of counterexamples. By categorizing and using them to create simulation tests, we provide an important body of stress testing scenarios. Analysis and discussion of counterexamples with system designers has already led to design improvements of the ACAS X system.

CURRENT AND FUTURE WORK

Aircraft collision avoidance systems such as ACAS X have a huge potential to make air travel safer and more efficient. They can help utilize the increasingly crowded airspace more efficiently. And when aircraft get too close to one another, they can help the pilot prevent a mid-air collision. Aircraft collision avoidance systems are expected to give good advisories under all circumstances while simultaneously being minimally invasive for the pilot. Among all possible advisories, the system is supposed to prefer choices that minimize the impact on the original flight trajectories while also avoiding the potential of NMAC. The advisories have to be issued fast enough so that they still resolve the NMAC before it is too late, respecting that the pilot may take a moment to react to the advisory and have flexibility in how exactly they follow it. And the system should not give advisories that could confuse the pilot under stress; for example, it should avoid changing its mind too quickly about whether to climb or descend. Furthermore, the system should work well whether the intruder is equipped with a collision avoidance system or not. These are a lot of responsibilities and safety requirements for a system like ACAS X, which make its design very challenging. How could it ever be possible to design a system respecting those constraints? And, once designed, how could one argue that it will really work as intended under all circumstances?

The safety analysis of aircraft collision avoidance systems is challenging, because of the many possible geometric configurations to consider, as well as complicated predictions of how the positions of the aircraft may change over time. Additionally, ACAS X uses a large machine-optimized score table for its decisions. This setting leads to a disciplined way of reaching optimal compromises between safety and operational suitability, but makes verification more difficult. through the table would be impossible and not even necessarily conclusive, even if performed at approximate grid resolutions.

This paper illustrates that the verification problem for ACAS X is not completely insurmountable, however. An appropriate separation of concerns in the verification design made it possible to verify collision freedom of controls respecting safe regions, under a broad class of behaviors in a precisely specified hybrid system model, and then subsequently compare those safe regions to the decisions of the score table. This is a big step forward in the safety analysis of ACAS X and contributed to a deeper understanding of ACAS X itself. At the same time, the formal verification of ACAS X presented in this paper can be improved in several different ways.

More Accurate Model. The hybrid model presented in Sect. 2 has a number of limitations. In many of our counterexamples, ACAS X would issue an advisory later on, correcting its initial assessment. In contrast, our current model assumes that an advisory is followed forever. We are working on an improved model handling such subsequent advisories. Our model also assumes that both aircraft follow a straight-line trajectory, with no possible turn.

In order to relax this limitation, we are working on extending the reduction presented in Sect. 2.2 to possibly curved trajectories of both aircraft.

Another limitation of our approach is that we do not handle intruders equipped with a collision avoidance system, or encounters with multiple intruders. This is made more challenging because those cases are typically handled by online scores after interpolation in the score table. But currently, our comparison looks only at the score table and ignores online scores -partly because considering them would add several orders of magnitude to an already large comparison calculation. We would like to extend our model to handle equipped and multiple intruders, and extend our comparison to handle some amount of online scores.

Fewer Counterexamples. The safe regions we identified may be overly conservative -in particular for the COC advisory -which leads to false alarms in the comparison with the actual advisories issued by the ACAS X table. Thus, we may flag an advisory as unsafe, although it is not, simply because of our worst case assumptions. One approach is to improve the precision of our safe regions, to make them less conservative -while keeping them provably safe. Another complementary approach is to sort and categorize the automatically generated counterexamples. This will make it easier for the development team to navigate through the millions of counterexamples we are generating.

Continuous Region Comparisons. So far, the table compar-

ison has been performed at the discrete grid points of the ACAS X table . A natural question one may ask is what happens at the points off this grid? Are the advisories safe for those points? Recall that the ACAS X decisions for off-grid points result from a multi-linear interpolation of the values on the neighboring grid points (vertices).

An important fact about the safe regions we identified is that they characterize which actions are safe at any point in the state space regardless of the ACAS X grid. We can thus reuse the same model, proof and safe regions to assess the safety of any ACAS X decision.

Although it is conceptually straight-forward to check if one particu-Figure 6: A 2-dimensional example of how interpolation can induce an advisory a to be optimal in the middle of a cell even if a is not optimal on any vertex. The regions in the diagram correspond to areas where different actions are optimal. The color of the circles represents the optimal action at the vertices.

lar point and advisory respect the safe region, checking the safety of an entire region -typically a hyperrectangle -requires first understanding the distribution of optimal ACAS X advisories within the hyperrectangle before checking a geometrical inclusion. Given the large number of cells in the score table, we crucially need a good understanding of the multi-linear interpolation as well as a scalable method to tackle the inclusion problem. For example, whenever all the vertices of the considered hyperrectangle have the same advisory, it is sufficient to only check the safety of the vertices. What is more challenging, however, is when a hyperrectangle exhibits different advisories on its vertices. Figure 6 shows a 2-dimensional projection where the vertices have different optimal advisories. Observe how the multi-linear interpolation computes a different optimal advisory in a non-obvious zone in the middle of the cell.

As a first attempt, our strategy to tackle this problem will be as follows. If a cell lies, even partially, in an area where a particular advisory a is unsafe while there exist some points in the cell where a is an optimal advisory, then we may flag these cell as potentially unsafe. This is conservative -the points where a is unsafe may be disjoint from the points where a is optimal. These potentially unsafe cells can be investigated more closely using sampling methods focused on the cell, or further decomposed into smaller cells and recursively analyzed. If, on the other hand, we can prove that no unsafe action can ever be an optimal action within the cell, then we have conclusively proven that the cell is safe under our model. We have a method of answering this cell-checking question by converting it into a sequence of simple convex optimization problems that soundly relax the original question.

Optimality vs. Safety. The ACAS X table is optimized to minimize the probability of a future NMAC while trying to maintain a good operational suitability. In contrast, the regions we identified focus solely on (local) safety aspects for straight-line encounters with no optimality objectives.

How can one get the best of both worlds, i.e., an optimal ACAS X table that is provably safe -at least for the most common encounters? One way of achieving such objective is to integrate the symbolic safe regions into the decisions from the ACAS X table, such that optimal advisories can only be issued if they are deemed safe with respect to the safe regions. One may also prefer provably safe advisories over ACAS X advisories whenever the safety conditions are violated. We are currently discussing with the development team the following three possible integration schemes:

1. The safe regions could be incorporated into the MDP optimization itself, e.g. by changing the objective function to have a prohibitive score for advisories outside the provable safe regions.

2. As in ModelPlex [START_REF] Mitsch | ModelPlex: Verified runtime validation of verified cyber-physical system models[END_REF], the final optimal advisory from the ACAS X table could be checked for compatibility with the safe regions and replaced by another safe action if needed.

3. Recall that ACAS X uses online scores which adapt the score from the ACAS X table, for example for ground proximity or coordination. The safe regions can be added as a new online score to ACAS X such that it avoids issuing advisories that are not known to be safe.

The above integration propositions feed the knowledge gained from the safety analysis back into the design of ACAS X to help make it an even safer system while preserving its main objective. Approach 1 gives the most direct feedback into the original system optimization but leaves the prohibitive score to use unspecified. Approach 2 is conceptually very clean and comes with strong theoretical guarantees [START_REF] Mitsch | ModelPlex: Verified runtime validation of verified cyber-physical system models[END_REF] about the final system whenever the assumptions of the hybrid model are met. It requires, however, an additional decision component into the system architecture. Finally, approach 3 leads to a particularly smooth integration, because the ACAS X infrastructure already uses online scores for various purposes on a routine basis. Further evaluation is needed to determine which score to use and how it would interact with other online scores.

Figure 1 :

 1 Figure 1: Reduction to a 2-dimensional encounter

3 [

 3 ( ( ?true ∪ 4 ḣf := * ; (w := -1 ∪ w := 1); 5 ?C impl (r, h, ḣ0); advisory := (w, ḣf ) );

6 a

 6 := * ;

Figure 2 :

 2 Figure2: Trajectory of ownship (red) and safe region for the intruder (green), immediate response[START_REF] Jeannin | A formally verified hybrid system for the next-generation airborne collision avoidance system[END_REF]

  (a) Pilot responds immediately to the advisory DNC (Do Not Climb). Starting state is r = 4,000 ft, rv = 180 ft/s , θv = 180 • , h = 1,200 ft, ḣ0 = 3,000 ft/min, ḣ1 = -3,000 ft/min, previous state = 'None-None'. There is an NMAC at time 0.00 s with 66.66 ft vertical separation and 0.40 ft horizontal separation. (b) Pilot responds to the advisory DNC after 3 seconds. Starting state is r = 650 ft, rv = 10 ft/s , θv = 130 • , h = 150 ft, ḣ0 = 2,000 ft/min, ḣ1 = 500 ft/min, previous state = 'None-None'. There is an NMAC at time 0.00 s with 81.79 ft vertical separation and 0.01 ft horizontal separation.

Figure 4 :

 4 Figure 4: Counterexamples with Induced Risk: Original ownship path (dashed blue) and intruder path (red) vs. potential ownship path (blue) responding to a do-not-climb (DNC) advisory with no delay (panel a) and with a 3-second delay (panel b).

  (a) Pilot responds immediately to COC (no advisory) followed 1 s later by Descend 1, 500 ft/min. Starting state is r = 150 ft, rv = 20 ft/s , θv = 170 • , h = 1,000 ft, ḣ0 = 0 ft/min, ḣ1 = -3,000 ft/min, previous state = 'None-None'. There is an NMAC at time 24.20 s with 99.84 ft vertical separation and 484 ft horizontal separation. (b) Pilot responds to COC (no advisory) after 3 s followed 1 s later by Subsequent-Descend 1, 500 ft/min. Starting state is r = 17,500 ft, rv = 2,350 ft/s θv = 180 • , h = 800 ft, ḣ0 = 10,000 ft/min, ḣ1 = 1,500 ft/min, previous state = 'DND-DND'. There is an NMAC at time 0.00 s with 10.99 ft vertical separation and 7.50 ft horizontal separation.

Figure 5 :

 5 Figure 5: Counterexamples of Insufficient Advice: Original ownship path (dashed blue) and intruder path (red) vs. potential ownship path (blue) responding to a COC advisory, i.e., continuing with no advisory for 1 second followed with the most extreme advisory available with no delay (panel a) and with 3-second delay (panel b).

Table 1 :

 1

			Advisories and their modeling variables [8, 7]
		ACAS X Run 13 Specification [9]			Our model
		Vertical Rate Range	Strength Delay (s)
	Advisory	Min (ft/min) Max (ft/min)	ar	dp	w	ḣf (ft/min)
	DNC2000	-∞	+2000	g/4	5	-1	+2000
	DND2000	-2000	+∞	g/4	5	+1	-2000
	DNC1000	-∞	+1000	g/4	5	-1	+1000
	DND1000	-1000	+∞	g/4	5	+1	-1000
	DNC500	-∞	+500	g/4	5	-1	+500
	DND500	-500	+∞	g/4	5	+1	-500

Table 2 :

 2 Results of the comparison

		With 3-second delay	With no delay
		Number of cutpoints Percentage Number of cutpoints Percentage
	Points examined	648,591,384,375	100.0%	648,591,384,375	100.0%
	Safe or unresolvable points	633,430,949,641	97.7%	632,951,653,646	97.6%
	Counterexamples	15,160,434,734	2.3%	15,639,730,729	2.4%

  For one thing, classical testing techniques that try to reach some kind of coverage (e.g. one test case per path through a program) are hopelessly insufficient for ACAS X, because most of its crucial decisions come from interpolating the individual data entries in the table. Testing what happens when following all paths

Acknowledgments

This research was conducted under the sponsorship of the Federal Aviation Administration Traffic Alert & Collision Avoidance System (TCAS) Program Office (PO) AJM-233 under contract number DTFAWA-11-C-00074. Additionally, support for the basic verification technology used as a foundation for this research was provided by the National Science Foundation under NSF CAREER Award CNS-1054246.

The authors would like to warmly thank Neal Suchy for his lead of the ACAS X project and support of this work, as well as Stefan Mitsch and Jan-David Quesel for their support of the KeYmaera tool. The authors would also like to thank Jeff Brush, Jessica Holland, Robert Klaus, Barbara Kobzik-Juul, Mykel Kochenderfer, Ted Londner, Sarah Loos, Ed Morehouse, Wes Olson, Michael Owen, Joshua Silbermann, and the ACAS X development team for interesting remarks.


Khalil Ghorbal

André Platzer

Characterizing Algebraic Invariants by Differential Radical Invariants

Keywords: Characterizing Algebraic Invariants by Differential Radical Invariants invariant algebraic sets, polynomial vector fields, real algebraic geometry, Zariski topology, higher-order Lie derivation, automated generation and checking, symbolic linear algebra, rank minimization, formal verification

HAL is

Introduction

Reasoning about the solutions of differential equations by means of their conserved functions and expressions is ubiquitous all over science studying dynamical processes. It is even crucial in many scientific fields (e.g. control theory or experimental physics), where a guarantee that the behavior of the system will remain within a certain predictable region is required. In computer science, the interest of the automated generation of these conserved expressions, so-called invariants, was essentially driven and motivated by the formal verification of different aspects of hybrid systems, i.e. systems combining discrete dynamics with differential equations for the continuous dynamics.

The verification of hybrid systems requires ways of handling both the discrete and continuous dynamics, e.g., by proofs [START_REF] Platzer | Differential dynamic logic for hybrid systems[END_REF], abstraction [START_REF] Rodríguez-Carbonell | An abstract interpretation approach for automatic generation of polynomial invariants[END_REF][START_REF] Tiwari | Abstractions for hybrid systems[END_REF], or approximation [START_REF] Lanotte | Taylor approximation for hybrid systems[END_REF]. Fundamentally, however, the study of the safety of hybrid systems can be shown to reduce constructively to the problem of generating invariants for their differential equations [START_REF] Platzer | Logics of dynamical systems[END_REF]. We focus on this core problem in this paper. We study the case of algebraic invariant equation, i.e. invariants described by a polynomial equation of the form h = 0 for a polynomial h. We also only consider algebraic differential equations (or algebraic vector fields), i.e. systems of ordinary differential equations in (vectorial) explicit form dx dt = p(x), with a polynomial right-hand side, p. The class of algebraic vector fields is far from restrictive and many analytic nonalgebraic functions, such as the square root, the inverse, the exponential or trigonometric functions, can be exactly modeled as solutions of ordinary differential equations with a polynomial vector field (a concrete example will be given in Section 6.2).

While algebraic invariant equations are not the only invariants of interest for hybrid systems [START_REF] Platzer | The structure of differential invariants and differential cut elimination[END_REF][START_REF] Platzer | A differential operator approach to equational differential invariants -(invited paper)[END_REF], they are still intimately related to all other algebraic invariants, such as semialgebraic invariants. We, thus, believe that the characterization we achieve in this paper to be an important step forward in understanding the invariance problem of polynomial vector fields, and hence the hybrid systems with polynomial vector fields.

Our results indicate that algebraic geometry is well suited to reason about and effectively compute algebraic invariant equations. Relevant concepts and results from algebraic geometry will be introduced and discussed as needed. The proofs of all presented results are available in [START_REF] Ghorbal | Characterizing algebraic invariants by differential radical invariants[END_REF].

Content In Section 2, we introduce a precise algebraic abstraction of the reachable set of the solution of a generic algebraic initial value problem. This abstraction is used to give a necessary and sufficient condition for a polynomial h to have the reachable set of the solution as a subset of the set of its roots. Section 3 builds on top of this characterization to, firstly, check the invariance of a variety candidate (Section 3.1) and, secondly, give an algebraic characterization for a variety to be an invariant for a polynomial vector field (Section 3.2). The characterization of invariant varieties is exploited in Section 4 where the generation of invariant varieties is reduced to symbolic linear algebra computation. The contributions of this work are summarized in Section 5. Finally, Section 6 presents three case studies to highlight the importance of our approach through concrete and rather challenging examples.

Sound and Precise Algebraic Abstraction by Zariski Closure

We consider autonomous1 algebraic initial value problems (see Def. 1 below). A nonautonomous system with polynomial time dependency can be reformulated as an autonomous system by adding a clock variable that reflects the progress of time. In this section, we investigate algebraic invariant equations for the considered initial value problems. This study is novel and will turn out to be fruitful from both the theoretical and practical points of view. The usual approach which assumes the initial value to be in a region of the space, often an algebraic set, will be discussed in Section 3.

Let x = (x 1 , . . . , x n ) ∈ R n , and x(t) = (x 1 (t), . . . , x n (t)), where x i : R → R, t → x i (t). The initial value x(t ι ) = (x 1 (t ι ), . . . , x n (t ι )) ∈ R n , for some t ι ∈ R, will be denoted by x ι . We do not consider any additional constraint on the dynamics, that is the evolution domain corresponds to the domain of definition.

Definition 1 (Algebraic Initial Value Problem). Let p i , 1 ≤ i ≤ n, be multivariate polynomials of the polynomial ring R[x]. An algebraic initial value problem is a pair of an explicit algebraic ordinary differential equations system (or polynomial vector field), p, and an initial value, x ι ∈ R n :

dx i dt = ẋi = p i (x), 1 ≤ i ≤ n, x(t ι ) = x ι . (1)
Since polynomial functions are smooth (C ∞ , i.e. they have derivatives of any order), they are locally Lipschitz-continuous. By Cauchy-Lipschitz theorem (a.k.a. Picard-Lindelöf theorem), there exists a unique maximal solution to the initial value problem (1) defined on some nonempty open set U t ⊆ R. A global solution defined for all t ∈ R may not exist in general. For instance, the maximal solution x(t) of the

1-dimensional system ẋ = x 2 , x(t ι ) = x ι = 0 is defined on R \ {t ι + x -1
ι }. Algebraic invariant equations for initial value problems are defined as follows.

Definition 2 (Algebraic Invariant Equation (Initial Value Problem)).

An algebraic invariant equation for the initial value problem (1) is an expression of the form h(x(t)) = 0 that holds true for all t ∈ U t , where h ∈ R[x] and x : U t → R n , is the (unique) maximal solution of (1).

Notice that any (finite) disjunction of conjunctions of algebraic invariant equations over the reals is also an algebraic invariant equation (w.r.t. Def. 2) using the following equivalence (R[x] is an integral domain):

i j f i,j = 0 ←→ i j f 2 i,j = 0 . (2)
In Def. 2, the function h(x(t)), and hence the polynomial h(x), depend on the fixed but unknown initial value x ι . We implicitly assume this dependency for a clearer notation and will emphasize it whenever needed. Also, observe that h(x(t)), seen as a real valued function of time t, is only defined over the open set U t ⊆ R since the solution x(t) is itself only defined over U t . The polynomial function h : R n → R; x → h(x) is, however, defined for all R n . Definition 3 (Orbit). The reachable set, or orbit, of the solution of Eq. (1),

x(t) is defined as O(x ι ) def = {x(t) | t ∈ U t } ⊆ R n .
The complete geometrical characterization of the orbit requires the exact solution of Eq. ( 1). Very few initial value problems admit an analytic solution, although a local approximation can be always given using Taylor series approximations (such approximation is for instance used in [START_REF] Lanotte | Taylor approximation for hybrid systems[END_REF] for the verification of hybrid systems). In this work, we introduce a sound abstraction of the orbit, O(x ι ), using (affine) varieties 2 . The idea is to embed the orbit (which is not a variety in general) in a variety to be defined. The embedding we will be using is a well-known topological closure operation in algebraic geometry called the Zariski closure ([6, Chapter 1]). Varieties, which are sets of points, can be represented and computed efficiently using their algebraic counterpart: ideals of polynomials. Therefore, we first recall three useful definitions: an ideal of the ring R[x], the variety of a subset of R[x], and finally the vanishing ideal of a subset of R n . Definition 4 (Ideal). An ideal I is a subset of R[x] that contains the polynomial zero (0), is stable under addition, and external multiplication. That is, for all h 1 , h 2 ∈ I, the sum h 1 + h 2 ∈ I; and if h ∈ I, then, qh ∈ I for all q ∈ R[x].

For a finite natural number r, we denote by h 1 , . . . , h r the subset of R[x] generated by the polynomials {h 1 , . . . , h r }, i.e. the set of linear combinations of the polynomials h i (where the coefficients are themselves polynomials):

h 1 , . . . , h r def = r i=1 g i h i | g 1 , . . . , g r ∈ R[x] .
By Def. 4, the set h 1 , . . . , h r is an ideal. More interestingly, by Hilbert's Basis Theorem [START_REF] Hilbert | Über die Theorie der algebraischen Formen[END_REF], any ideal I of the Noetherian ring R[x] can be finitely generated by, say {h 1 , . . . , h r }, so that I = h 1 , . . . , h r .

Given Y ⊆ R[x], the variety (over the reals), V (Y ), is a subset of R n defined by the common roots of all polynomials in Y . That is,

V (Y ) def = x ∈ R n | ∀h ∈ Y, h(x) = 0 .
The vanishing ideal (over the reals), I(S), of S ⊆ R n is the set of all polynomials that evaluates to zero for all x ∈ S:

I(S) def = h ∈ R[x] | ∀x ∈ S, h(x) = 0 .
The Zariski closure Ō(x ι ) of the orbit O(x ι ) is defined as the variety of the vanishing ideal of O(x ι ):

Ō(x ι ) def = V (I(O(x ι ))) . (3)
That is, Ō(x ι ) is defined as the set of all points that are common roots of all polynomials that are zero everywhere on the orbit O(x ι ). The variety Ō(x ι ) soundly overapproximates all reachable states x(t) in the orbit of O(x ι ), including the initial value x ι :

Proposition 1 (Soundness of Zariski Closure). O(x ι ) ⊆ Ō(x ι ) .
Therefore, all safety properties that hold true for Ō(x ι ), are also true for O(x ι ). The soundness in Proposition 1 corresponds to the reflexivity property of the Zariski closure: for any subset S of R n , S ⊆ V (I(S)). Besides, the algebraic geometrical fact that the Zariski closure Ō(x ι ) is the smallest3 variety containing O(x ι ) corresponds to the fact that Ō(x ι ) is the most precise algebraic abstraction of O(x ι ).

Observe that if the set of generators of I(O(x ι )) is only the zero polynomial, I(O(x ι )) = 0 , then Ō(x ι ) = R n (the whole space) and the Zariski closure fails to be informative. For instance, for (non-degenerated) one dimensional vector fields (n = 1) that evolve over time, the only univariate polynomial that has infinitely many roots is the zero polynomial. This points out the limitation of the closure operation used in this work and raises interesting question about how to deal with such cases (this will be left as future work).

The closure operation abstracts time. This means that Ō(x ι ) defines a subset of R n within which the solution always evolves without saying anything about where the system will be at what time (which is what a solution would describe and which is exactly what the abstraction we are defining here gets rid off). In particular, Ō(x ι ) is independent of whether the system evolves forward or backward in time.

Although, we know that I(O(x ι )) is finitely generated, computing all its generators may be intractable. By the real Nullstellensatz, vanishing ideals over R are in fact exactly the real radical ideals [1, Section 4.1]. In real algebraic geometry, real radical ideals are notoriously hard 4 to compute. However, we shall see in the sequel that Lie derivation will give us a powerful computational handle that permits to tightly approximate (and even compute in some cases) I(O(x ι )). The Lie derivative of a polynomial along a vector field is defined as follows.

Definition 5 (Lie Derivative). The Lie derivative of h ∈ R[x] along the vector field p = (p 1 , . . . , p n ) is defined by:

L p (h) def = n i=1 ∂h ∂x i p i . (4)
Higher-order Lie derivatives are defined recursively:

L (k+1) p (h) def = L p (L (k) p (h)), with L (0) p (h) def = h.
We state an important property of the ideal I(O(x ι )). Similar result is known under different formulations ([23, Theorem 3.1] and [START_REF] Platzer | Logical Analysis of Hybrid Systems: Proving Theorems for Complex Dynamics[END_REF]Lemma 3.7]).

Proposition 2. I(O(x ι )) is a differential ideal for L p , i.e. it is stable under the action of the L p operator. That is, for all h ∈ I(O(x ι )), L p (h) ∈ I(O(x ι )).

In the next section, we give a necessary and sufficient condition for a polynomial h to be in I(O(x ι )), that is for the expression h = 0 to be an algebraic invariant equation for the initial value problem (1), i.e. h evaluates to 0 all along the orbit of x ι .

Differential Radical Characterization

In this section, we study the algebraic properties of the Zariski closure Ō(x ι ) defined in the previous section. We then define and characterize invariant algebraic sets of polynomial vector fields.

For h ∈ R[x], we recursively construct an ascending chain of ideals of R[x] by appending successive Lie derivatives of h to the list of generators:

h ⊂ h, L (1) p (h) ⊂ • • • ⊂ h, . . . , L (N -1) p (h) = h, . . . , L (N ) p (h) .
Since the ring R[x] is Noetherian, the chain above has necessarily a finite length: the maximal ideal (in the sense of inclusion), so-called the differential radical ideal 5 of h, will be denoted by Lp h . Its order is the smallest N such that:

L (N ) p (h) ∈ L (0) p (h), . . . , L (N -1) p (h) . (5)
The following theorem, an important contribution of this work, states a necessary and sufficient condition for a polynomial h to be in I(O(x ι )).

Theorem 1 (Differential Radical Characterization). Let h ∈ R[x], and let N denote the order of

Lp h . Then, h ∈ I(O(x ι )) if and only if 0≤i≤N -1 L (i) p (h)(x ι ) = 0 . (6)
The statement of Theorem 1 is general and assumes nothing about x ι ∈ R n . A natural question to ask is how differential radical characterization can be used to reason about invariant regions of a given polynomial vector field. By invariant (or stable) regions, we mean, regions S ⊂ R n from which the trajectory of the solution of the initial value problem (1), with x ι ∈ S, can never escape. In particular, we focus on invariant algebraic sets where S is variety. Definition 6 (Invariant Variety). The variety S is an invariant variety for the vector field p if and only if ∀x ι ∈ S, O(x ι ) ⊆ S.

Dual to the geometrical point of view in Def. 6, the algebraic point of view is given by extending the definition of algebraic invariant equation for initial value problems (Def. 2), to algebraic invariant equation for polynomial vector fields. Definition 7 (Algebraic Invariant Equation (Vector Field)). The expression h = 0 is an algebraic invariant equation for the vector field p if and only if V ( h ) is an invariant variety for p. Unlike Def. 2, Def. 7, or its geometrical counterpart, Def. 6, corresponds to the typical object of studies in hybrid system verification as they permit the abstraction of the continuous part by means of algebraic equations. In the two following sections, we show how differential radical characterization (Theorem 1) can be used to address two particular questions: checking the invariance of a variety candidate (Section 3.1) and characterizing invariant varieties (Section 3.2).

We will say that the polynomial h is a differential radical invariant (for p) if and only if V Lp h is an invariant variety for p. 5 The construction of Lp h is very similar to the construction of the radical of an ideal except with higher-order Lie derivatives in place of higher powers of polynomials.

Checking Invariant Varieties by Differential Radical Invariants

The problem we solve in this section is as follows: given a polynomial vector field p, can we decide whether the equation h = 0 is an algebraic invariant equation for the vector field p ? Dually, we want to check whether the variety V ( h ) is invariant for p. Theorem 2 solves the problem.

Theorem 2. Let h ∈ R[x]
, and let N denote the order of Lp h . Then, V ( h ) is an invariant variety for the vector field p (or equivalently h = 0 is an algebraic invariant equation for p) if and only if

h = 0 → 1≤i≤N -1 L (i) p (h) = 0 . ( 7
)
Corollary 1 (Decidability). It is decidable whether the expression h = 0 is an algebraic invariant equation for the vector field p assuming real algebraic coefficients for h and p.

The sound and complete related proof rule from Theorem 2 can be written as follows (N denotes the order of Lp h ):

(DRI) h = 0 → 1≤i≤N -1 L (i) p (h) = 0 (h = 0) → [ ẋ = p](h = 0) . (8)
Using the naive trick in Eq. ( 2), theoretically, the proof rule can be easily extended to check for the invariance of any finite disjunction of conjunctions of algebraic invariant equations for p. This means that we can check for the invariance of any variety for p, given its algebraic representation. However, in practice, other techniques, outside the scope of this paper, should be considered to try to keep the degree of the involved polynomials as low as possible.

Differential Radical Characterization of Invariant Varieties

In the previous section, we were given a variety candidate of the form V ( h ) and asked whether we can decide for its invariance. In this section, we characterize all invariant varieties of a vector field p using a differential radical criterion. The following theorem fully characterizes invariant varieties of polynomial vector fields.

Theorem 3 (Characterization of Invariant Varieties).

A variety S is an invariant variety for the vector field p if and only if there exists a polynomial h such that S = V Lp h . As a consequence, every invariant variety corresponds to an algebraic invariant equation involving a polynomial and its higher-order Lie derivatives (N denotes the order of Lp h ):

0≤i≤N -1 L (i) p (h) = 0 . (9)
Observe how Theorem 3 proves, from the differential radical characterization point of view, the well-known fact about invariant polynomial functions [17, Theorem 3]: If L p (h(x)) = 0, then, for any c ∈ R, Lp h(x) -c = h(x) -c , and so S = V ( h(x) -c ) is an invariant variety for p.

An algebraic invariant equation for p is defined semantically (Def. 7) as a polynomial that evaluates to zero if it is zero initially (admits x ι as a root). Differential radical invariants are, on the other hand, defined as a structured, syntactically computable, conjunction of polynomial equations involving one polynomial and its successive Lie derivatives. By Theorem 3, both coincide.

The explicit formulation of Eq. ( 5), namely

L (N ) p (h) = N -1 i=0 g i L (i) p (h), (10)
for some

g i ∈ R[x],
is computationally attractive as it only involves polynomial arithmetic on higher-order Lie derivatives of one polynomial, h, which in turn can be computed automatically by symbolic differentiation. Section 4 exploits this fact to automatically generate differential radical invariants and consequently invariant varieties.

Effective Generation of Invariant Varieties

In the previous section, we have seen (Theorem 3) that differential radical ideals characterize invariant varieties. Based on Eq. ( 10), we explain in this section how we automatically construct differential radical ideals given a polynomial vector field p by deriving a set of constraints that the coefficients of a parametrized polynomial have to satisfy. The degree of a polynomial in R[x] is defined as the maximum degree among the (finite) set of degrees of its monomials 6 . When the degrees of all nonzero monomials of a polynomial h are equal, we say that h is homogeneous, or a form, of degree d.

By introducing an extra variable x 0 and multiplying all monomials by a suitable power of x 0 , any polynomial of R[x] can be homogenized to a form in R[x 0 ][x]. The additional variable x 0 is considered as a time-independent function: its time derivative is zero ( ẋ0 = p 0 = 0). "De-homogenizing" the vector field corresponds to instantiating x 0 with 1, which gives back the original vector field. Geometrically, the homogenization of polynomials corresponds to the notion of projective varieties in projective geometry, where the homogenized polynomial is the algebraic representative of the original variety in the projective plane [START_REF] Cox | Ideals, Varieties, and Algorithms: An Introduction to Computational Algebraic Geometry and Commutative Algebra[END_REF]Chapter 8].

From a computational prospective, working in the projective plane offers a more symmetric representation: all monomials of a form have the same degree. The arithmetic of degrees over forms is also simplified: the degree of a product is the sum of the degrees of the operands. In the reminder of this section, we only consider forms of R[x 0 , . . . , x n ]. The symbol x will denote the vector of all involved variables.

If h denotes a form of degree d, and d the maximum degree among the degrees of p i , then the degree of the polynomial L (k) p (h) is given by:

deg(L (k) p (h)) = d + k(-1 + d ) . (11)
Recall that a form of degree d in R[x 0 , . . . , x n ] has

m d def = n+d d (12)
monomials (the binomial coefficient of n + d and d). A parametrized form h α of degree d can therefore be represented by its symbolic coefficients' vector α : R m d . For this representation to be canonical, one needs to fix an order over monomials of the same degree. We will use the usual lexicographical order, except for x 0 :

x 1 > x 2 > • • • > x n > x 0 .
We first compare the degrees of x 1 , if equal, we compare the degrees of x 2 and so on till reaching x n and then x 0 . For instance, for n = 2, a parametrized form

h α of degree d = 1 is equal to α 1 x 1 + α 2 x 2 + α 3 x 0 . Its related coefficients' vector is α = (α 1 , α 2 , α 3 ).
Let h α be a parametrized form of degree d and let α = (α 1 , . . . , α m d ) denote the coefficients' vector with respect to the monomial order defined above. Since all polynomials in Eq. ( 10) are forms in projective coordinates, the degree of each term g i L (i) p (h α ) matches exactly the degree of L (N ) p (h α ). Hence, by Eq. ( 11

): deg(g i ) = (N -i)(-1 + d ).
The coefficients' vector of each form g i is then a vector, β i , of size m (N -i)(-1+d ) (see Eq. ( 12)). So that we obtain m d+N (-1+d ) biaffine equations: linear in α i , 1 ≤ i ≤ m d , and affine

β i,j , 0 ≤ i ≤ N -1, 1 ≤ j ≤ m (N -i)(-1+d ) . A concrete example is as follows. Example 1. Suppose we have n = 2, d = 1, p 1 = a 1 x 1 + a 2 x 2 and p 2 = b 1 x 1 + b 2 x 2 .
For d = 1, the parametrized form h α is equal to α 1 x 1 + α 2 x 2 + α 3 x 0 . Let N = 1. The first-order Lie derivative, L p (h α ), has the same degree, 1, and is equal to α 1 (a 1 x 1 + a 2 x 2 ) + α 2 (b 1 x 1 + b 2 x 2 ). In this case, g is a form of degree 0, that is a real number. So it has one coefficient β ∈ R. We, therefore, obtain m 1 = 3 1 = 3 constraints:

(-a 1 + β)α 1 + (-b 1 )α 2 = 0 (-a 2 )α 1 + (-b 2 + β)α 2 = 0 (β)α 3 = 0 ↔   -a 1 + β -b 1 0 -a 2 -b 2 + β 0 0 0 β   .   α 1 α 2 α 3   = 0 .
As suggested in Example 1, for a given d and N , and if we concatenate all vectors β i into one vector β, the equational constraints can be rewritten as a symbolic linear algebra problem of the following form:

M d,N (β)α = 0, (13)
where α and β are decoupled. The matrix M d,N (β) is called the matrix representation of the ideal membership problem L (N ) p (h α ) ∈? h α , . . . , L (N -1) p (h α ) . Recall that the kernel (or null-space) of a matrix M ∈ R r×c , with r rows and c columns is the subspace of R c defined as the preimage of the vector 0 ∈ R c :

ker(M ) def = {x ∈ R c | M x = 0} . Let s = dim(ker(M d,N (β))) ≤ m d .
If, for all β, s = 0, then the kernel is {0}. Hence, α = 0 and, for the chosen N , we have h α = 0: the only differential radical ideal generated by a form of degree d is the trivial ideal 0 . If, however, s ≥ 1, then, by Theorem 3, we generate an invariant (projective) variety for p. In this case, dehomogenizing is not always possible. In fact, the constraint on the initial value could involve x 0 , which prevents the de-homogenization (see Example 2). Otherwise, we recover an invariant (affine) variety for the original vector field. This is formally stated in the following theorem.

Theorem 4 (Effective Generation of Projective Invariant Varieties).

Let h α denote a parametrized form of degree d. There exists a real vector β such that dim(ker

(M d,N (β))) ≥ 1 if and only if for α ∈ ker(M d,N (β)), V Lp h α ⊂
R n+1 is a projective invariant variety for the homogenized vector field.

When s = dim(ker(M d,N (β))) ≥ 1, the subspace ker(M d,N (β)) is spanned by s vectors, e 1 , . . . , e s ∈ R m d , and for α = γ 1 e 1 +• • •+γ s e s , for arbitrarily (γ 1 , . . . , γ s ) ∈ R s , the variety V Lp h α is a family of invariant varieties of p (parametrized with γ). In the sequel, we give a sufficient condition, so that, for any given initial value, one gets a variety (different from the trivial whole space) that embeds the reachable set of the trajectory, O(x ι ). For instance, for conservative Hamiltonian system, if the total energy function, h, is polynomial (such as the energy function of the perfect pendulum), then, for any initial value 6), then, by Theorem 1, h α ∈ I(O(x ι )), and Ō(x ι ) ⊆ V Lp h α ([5, Corollary 1]). However, for x ι to satisfy Eq. ( 6), α must be in the intersection of N hyperplanes, H 0 , . . . , H N -1 , each defined explicitly by the condition L (i) p (h α )(x ι ) = 0:

x ι , O(x ι ) ⊆ V Lp h(x) -h(x ι ) = V ( h(x) -h(x ι ) ). For a generic x ι ∈ R n , if x ι satisfies Eq. (
H i def = α ∈ R m d | L (i) p (h α )(x ι ) = 0 . ( 14
)
Proposition 3 (Effective Sound Approximation of O(x ι )). Let h α be a parametrized form of degree d, and M d,N (β) the matrix representation of Eq. [START_REF] Lanotte | Taylor approximation for hybrid systems[END_REF].

Let H i ⊆ R m d , 0 ≤ i ≤ N -1
, be the hyperplanes defined in Eq. [START_REF] Neuhaus | Computation of real radicals of polynomial ideals II[END_REF]. Then, O(x ι ) ⊆ V Lp h α , if there exists β such that:

dim(ker(M d,N (β))) > m d -dim N -1 i=0 H i . ( 15
)
The remainder of this section discusses our approach to maximize the dimension of the kernel of M d,N (h), as well as the complexity of the underlying computation.

Gaussian Elimination Let β = (β 1 , . . . , β s ) : R s . By Theorem 4, we want to find an instance, β * , of β that maximizes dim ker(M d,N (β)), where all the elements of M d,N (β) are affine in β. At each iteration, our algorithm [5, Algorithm 1] assigns new values to the remaining coefficients in β for the matrix M d,N (β) to maximize the dimension of its kernel. A set, M , gathers all the instantiations of M d,N (β). The procedure ends when no further assignment can be done. The algorithm is in fact a typical MapReduce procedure which can be parallelized. A naive approach would be to first extract a basis for the matrix M d,N (β) (which requires symbolic computation capabilities for linear algebra), then, solves for βs that zero the determinant. In practice, however, row-reducing speeds up the computation: we row-reduce M d,n (β), and record any divisions by the pivot element: we then branch with any β that zero the denominator.

Example 2. We apply the algorithm sketched above to Example 1. The determinant of the matrix

M 1,1 (β) is β β 2 -(a 1 + b 2 )β -a 2 b 1 + a 1 b 2 .
Since we do not have any constraints on the parameters a 1 , a 2 , b 1 , b 2 , the only generic solution for the determinant is β = 0. The kernel of M 1,1 (0), of dimension 1, is generated by (0, 0, 1). The only candidates in this case are h α (x) = γx 0 , γ ∈ R. If we de-homogenize (set x 0 to 1), then, γ = 0 and we find the trivial invariant variety, R n .

The result of Example 2 is expected as it studies a generic linear vector field without any a priori constraints on the parameters. This triggers, naturally, an interesting feature of the differential radical characterization: its ability to synthesize vector fields to enforce an invariant variety. For instance, in Example 2, let

δ def = (a 1 -b 2 ) 2 + 4a 2 b 1 .
If δ ≥ 0, and a 2 = 0, then the kernel of M 1,1 (β) is generated by the vector a 1 -b 2 ± √ δ, 2a 2 , 0 (which is an eigenvector of the matrix M 1,1 (β)). By Theorem 4, we have an invariant variety given by: a 1 -b 2 ± √ δ x 1 + 2a 2 x 2 = 0. This is also expected for linear systems as eigenvectors span stable subspaces.

Complexity By Theorem 4, the generation of invariant varieties is equivalent to maximizing the dimension of the kernel of the matrix M d,N (β) over unconstrained β, which is in turn equivalent to the following unconstrained minimal rank problem:

min β rank(M d,N (β)), (16)
where the elements of the vector β are in R. If the vector field p has no parameters, then the entries of the matrix M d,N (β) are either elements of β or real numbers. Under these assumptions, the problem ( 16) is in PSPACE [2, Corollary 20] over the field of real numbers 7 , and is at least NP-hard (see [START_REF] Buss | The computational complexity of some problems of linear algebra[END_REF]Corollary 12] and [8, Theorem 8.2]) independently from the underlying field. In fact, deciding whether the rank of M d,N (β) is less than or equal to a given fixed bound is no harder than deciding the corresponding existential first-order theory.

On the other hand, there is an NP-hard lower bound for the feasibility of the original set of (biaffine) equations in β and α given in Eq. ( 13). In the simpler bilinear case and, assuming, as above, that the vector field has no parameters, finding a nontrivial solution (α = 0 is trivial) is also NP-hard [8, Theorems 3.7 and 3.8].

Related Work and Contributions

The contribution of this work is fourfold.

Sound and Precise Algebraic Abstraction of Reachable Sets (Section 2) Unlike previous work [START_REF] Tiwari | Nonlinear systems: Approximating reach sets[END_REF][START_REF] Sankaranarayanan | Automatic invariant generation for hybrid systems using ideal fixed points[END_REF][START_REF] Matringe | Generating invariants for non-linear hybrid systems by linear algebraic methods[END_REF][START_REF] Liu | Computing semi-algebraic invariants for polynomial dynamical systems[END_REF], we start by studying algebraic initial value problems. We propose a sound abstraction (Proposition 1) to embed (overapproximate) the reachable set. Our abstraction relies on the Zariski closure operator over affine varieties (closed sets of the Zariski topology), which allows a clean and sound geometrical abstraction. From there, we define the vanishing ideal of the closure, and give a necessary and sufficient condition (Theorem 1) for a polynomial equation to be an invariant for algebraic initial value problems.

Checking Invariant Varieties by Differential Radical Invariants (Section 3.1) The differential radical characterization allows to check for and falsify the invariance of a variety candidate. Unlike already existing proof rules [START_REF] Tiwari | Nonlinear systems: Approximating reach sets[END_REF][START_REF] Matringe | Generating invariants for non-linear hybrid systems by linear algebraic methods[END_REF][START_REF] Platzer | A differential operator approach to equational differential invariants -(invited paper)[END_REF], which are sound but can only prove a restrictive class of invariants. From Theorem 2, we derive a sound and complete proof rule (Eq. ( 8)) and prove that the problem is decidable (Corollary 1) over the real-closed algebraic fields.

Differential Radical Characterization of Invariant Varieties (Section 3.

2) The differential radical criterion completely characterizes all invariant varieties of polynomial vector fields. This new characterization (Theorem 3) permits to relate invariant varieties to a purely algebraic, well-behaved, conjunction of polynomial equations involving one polynomial and its successive Lie derivatives (Eq. ( 9)). It naturally generalizes [START_REF] Lafferriere | Symbolic reachability computation for families of linear vector fields[END_REF][START_REF] Tiwari | Approximate reachability for linear systems[END_REF] where linear vector fields are handled and [START_REF] Sankaranarayanan | Constructing invariants for hybrid systems[END_REF][START_REF] Matringe | Generating invariants for non-linear hybrid systems by linear algebraic methods[END_REF] where only a restrictive class of invariant varieties is considered.

Effective Generation of Invariant Varieties (Section 4) Unlike [START_REF] Tiwari | Nonlinear systems: Approximating reach sets[END_REF][START_REF] Sankaranarayanan | Automatic invariant generation for hybrid systems using ideal fixed points[END_REF][START_REF] Liu | Computing semi-algebraic invariants for polynomial dynamical systems[END_REF][START_REF] Rodríguez-Carbonell | Generating polynomial invariants for hybrid systems[END_REF], we do not use quantifier elimination procedures nor Gröbner Bases algorithms for the generation of invariant varieties. We have developed and generalized the use of symbolic linear algebra tools to effectively generate families of invariant varieties (Theorem 4) and to soundly overapproximate reachable sets (Proposition 3). In both cases, the problem requires maximizing the dimension of the kernel of a symbolic matrix. The complexity is shown to be NP-hard, but in PSPACE, for polynomial vector fields without parameters. We also generalize the previous related work on polynomial-consecution. In particular, Theorems 2 and 4 in [START_REF] Matringe | Generating invariants for non-linear hybrid systems by linear algebraic methods[END_REF] are special cases of, respectively, Theorem 4 and Proposition 3, when the order of differential radical ideals is exactly 1.

Case Studies

The following challenging example comes up as a subsystem we encountered when studying aircraft dynamics:

p 1 = -x 2 , p 2 = x 1 , p 3 = x 2
4 , p 4 = x 3 x 4 . It appears frequently whenever Euler angles and the three dimensional rotational matrix is used to describe the dynamics of rigid body motions. For some chosen initial value, such as x ι = (1, 0, 0, 1), it is an exact algebraic encoding of the trigonometric functions : x 1 (t) = cos(t), x 2 (t) = sin(t), x 3 (t) = tan(t), x 4 (t) = sec(t). When d = 2 and N = 1, the matrix M 2,1 (β) is 35 × 15, with 90 (out of 525) nonzero elements, and |β| = 5. The maximum dimension of ker(M 2,1 (β)) is 3 attained for β = 0. The condition of Proposition 3 is satisfied and, for any x ι , we find the following algebraic invariant equations for the corresponding initial value problem:

h 1 = x 2 1 + x 2 2 -x ι 2 1 -x ι 2 2 = 0, h 2 = -x 2 3 + x 2 4 + x ι 2 3 -x ι 2 4 = 0 .
In particular, for the initial value x ι = (1, 0, 0, 1), one recovers two trigonometric identities, namely cos(t) 2 + sin(t) 2 -1 = 0 for h 1 and -tan(t) 2 + sec(t) 2 -1 = 0 for h 2 .

For N = 3, the matrix M 2,3 (β) is 126 × 15, with 693 (out of 1890) nonzero elements, and |β| = 55. We found a β for which the dimension of ker(M 2,3 (β)) is 5. By Theorem 4, we have a family of invariant varieties for p encoded by the following differential radical invariant: h = γ 1 -x 2 3 γ 2 +x 2 4 γ 2 +x 2 x 4 γ 3 +x 2 1 γ 4 +x 2 2 γ 4 +x 1 x 4 γ 5 , where γ i , 1 ≤ i ≤ 5, are real numbers. In particular, when (γ 1 , γ 2 , γ 3 , γ 4 , γ 5 ) = (1, 0, 0, 0, 1), we have the following algebraic invariant equation for p:

-1 + x 1 x 4 = 0 ∧ -x 2 x 4 + x 3 = 0 ∧ -1 -x 2 3 + x 2 4 = 0 . (17)
Interestingly, since x ι = (1, 0, 0, 1) satisfies the above equations, we recover, respectively, the following trigonometric identities:

-1 + cos(t) sec(t) = 0 ∧ -sin(t) sec(t) + tan(t) = 0 ∧ -1 -tan(t) 2 + sec(t) 2 = 0 .
We stress the fact that Eq. ( 17) is one algebraic invariant equation for p. In fact, any conjunct alone, a part from -1 -x 2 3 + x 2 4 = 0, of Eq. ( 17) is not an algebraic invariant equation for p. Indeed, we can falsify the candidate -1 + x 1 x 4 = 0 using Theorem 2: the implication -1 + x 1 x 4 = 0 -→ -x 2 x 4 + x 3 = 0 is obviously false in general.

Notice that h 1 and h 2 can be found separately by splitting the original vector field into two separate vector fields since the pairs (p 1 , p 2 ) and (p 3 , p 4 ) can be decoupled. However, by decoupling, algebraic invariant equation such as Eq. ( 17) cannot be found. This clearly shows that in practice, splitting the vector field into independent ones should be done carefully when it comes to generating invariant varieties. This is somehow counter-intuitive as decoupling for the purpose of solving is always desirable. In fact, any decoupling breaks an essential link between all involved variables: time.

We proceed to discuss collision avoidance of two airplanes (Section 6.1) and then the use of invariant varieties to tightly capture the vertical motion of an airplane (Section 6.2).

Collision Avoidance

We revisit the linear vector field encoding Dubin's vehicle model for aircrafts [START_REF] Dubins | On curves of minimal length with a constraint on average curvature, and with prescribed initial and terminal positions and tangents[END_REF]. Although the system was discussed in many recent papers [START_REF] Platzer | Computing differential invariants of hybrid systems as fixedpoints[END_REF][START_REF] Sankaranarayanan | Automatic invariant generation for hybrid systems using ideal fixed points[END_REF][START_REF] Liu | Computing semi-algebraic invariants for polynomial dynamical systems[END_REF], we want to highlight an additional algebraic invariant equation that links both airplanes when turning with the same angular velocity. The differential equation system is given by:

p 1 = ẋ1 = d 1 , p 2 = ẋ2 = d 2 , p 3 = ḋ1 = -ω 1 d 2 , p 4 = ḋ2 = ω 1 d 1 , p 5 = ẏ1 = e 1 , p 6 = ẏ2 = e 2 , p 7 = ė1 = -ω 2 e 2 , p 8 = ė2 = ω 2 e 1 .
The angular velocities ω 1 and ω 2 can be either zero (straight line flight) or equal to a constant ω which denotes the standard rate turn (typically 180 • /2mn for usual commercial airplanes). When the two airplanes are manoeuvring with the same standard rate turn ω, apart from the already known invariants, we discovered the following differential radical invariant (which corresponds to a family of invariant varieties):

h 1 = γ 1 d 1 + γ 2 d 2 + γ 3 e 1 + γ 4 e 2 = 0 ∧ h 2 = γ 2 d 1 -γ 1 d 2 + γ 4 e 1 -γ 3 e 2 = 0,
for an arbitrarily (γ 1 , . . . , γ 4 ) ∈ R 4 . We have Lp h 1 = Lp h 2 = h 1 , h 2 . Observe also that V ( h 1 ) and V ( h 2 ) are not invariant varieties for p.

Longitudinal Motion of an Airplane

The full dynamics of an aircraft are often separated (decoupled) into different modes where the differential equations take a simpler form by either fixing or neglecting the rate of change of some configuration variables [START_REF] Stengel | Flight Dynamics[END_REF]. The first standard separation used in stability analysis gives two main modes: longitudinal and lateral-directional. We study the 6th order longitudinal equations of motion as it captures the vertical motion (climbing, descending) of an airplane. We believe that a better understanding of the envelope that soundly contains the trajectories of the aircraft will help tightening the surrounding safety envelope and hence help trajectory management systems to safely allow more dense traffic around airports. The current safety envelope is essentially a rough cylinder that doesn't account for the real capabilities allowed by the dynamics of the airplane. We use our automated invariant generation techniques to characterize such an envelope. The theoretical improvement and the effective underlying computation techniques described earlier in this work allow us to push further the limits of automated invariant generation. We first describe the differential equations (vector field) then show the nontrivial energy functions (invariant functions for the considered vector field) we were able to generate. Let g denote the gravity acceleration, m the total mass of an airplane, M the aerodynamic and thrust moment w.r.t. the y axis, (X, Z) the aerodynamics and thrust forces w.r.t. axis x and z, and I yy the second diagonal element of its inertia matrix. The restriction of the nominal flight path of an aircraft to the vertical plane reduces the full dynamics to the following 6 differential equations [25, Chapter 5] (u: axial velocity, w: vertical velocity, x: range, z: altitude, q: pitch rate, θ: pitch angle): u = X m -g sin(θ) -qw ẋ = cos(θ)u + sin(θ)w θ = q ẇ = Z m + g cos(θ) + qu ż = -sin(θ)u + cos(θ)w q = M I yy .

We encode the trigonometric functions using two additional variables for cos(θ) and sin(θ), making the total number of variables equal to 8. The parameters are considered unconstrained. Unlike [START_REF] Sankaranarayanan | Automatic invariant generation for hybrid systems using ideal fixed points[END_REF], we do not consider them as new time independent variables. So that the total number of state variables (n) and hence the degree of the vector field are unchanged. Instead, they are carried along the symbolic row-reduction computation as symbols in M d,N (β). For the algebraic encoding of the above vector field (n = 8), the matrix M 3,1 (β) is 495 × 165, with 2115 (out of 81675) nonzero elements, and |β| = 9.

We were able to automatically generate the following three invariant functions: We substituted the intermediate variables that encode sin and cos back to emphasize the fact that algebraic invariants and algebraic differential systems are suitable to encode many real complex dynamical systems. Using our Mathematica implementation, the computation took 1 hour on a recent laptop with 4GB and 1.7GHz Intel Core i5.

Conclusion

For polynomial vector fields, we give an algebraic characterization of invariant varieties. This so-called differential radical characterization makes it possible to decide for the invariance of a given variety candidate. It is, in addition, computationally attractive: generating invariant varieties requires minimizing the rank of a symbolic matrix and is hence at least NP-hard. The case studies show how the technique applies successfully to rather complex systems. We also revisited some known problems in the literature to exemplify the benefits of having a necessary and sufficient condition: all other known sound approaches generate a special class of invariant varieties (i.e. miss others).

In the future, we plan to investigate upper bounds for the order of the differential radical ideal of a given polynomial. Also, invariant varieties are not the only invariant of interest for polynomial vector fields, we want to consider semialgebraic sets as they play a prominent role in both hybrid systems and control theory. Finally, the effective use of algebraic invariants in general in the context of hybrid systems is still a challenging problem that we want to explore in more depth.

Autonomous means that the rate of change of the system over time depends only on the system's state, not on time.

In the literature, some authors use the terminology algebraic sets so that varieties is reserved for irreducible algebraic sets. Here we will use both terms equally.

Smallest here is to be understood w.r.t. to the usual geometrical sense, that is, any other variety containing O(xι), contains also its closure Ō(xι).

Given an ideal I ⊆ R[x], the degree of the polynomials that generate its real radical is bounded by the degree of polynomials that generate I to the power of 2 O(n 2 ) [14, Theorem

5.9].

The degree of the zero polynomial (0) is undefined. We assume in this work that all finite degrees are acceptable for the zero polynomial.

The complexity class depends on the underlying field and is worse for fields with nonzero characteristic.

Acknowledgments We thank the anonymous reviewers for their careful reading and detailed comments. We also would like to very much thank JEAN-BAPTISTE JEANNIN and ANDREW SOGOKON for the multiple questions, various comments and fruitful objections they both had on an early version of this work. We are finally grateful to ERIC GOUBAULT and SYLVIE PUTOT for the relevant references they pointed out to us on the integrability theory of nonlinear systems.

This material is based upon work supported by the National Science Foundation by NSF CA-REER Award CNS-1054246, NSF EXPEDITION CNS-0926181 and grant no. CNS-0931985. This research is also partially supported by the Defense Advanced Research Agency under contract no. DARPA FA8750-12-2-0291.