Sonia Cafieri
email: sonia.cafieri@enac.fr

Andrew R Conn

Marcel Mongeau
email: marcel.mongeau@enac.fr

Andy Conn

An MINLP and a continuous optimization formulations for aircraft conflict avoidance via heading and speed deviations

Keywords: mixed-integer nonlinear optimization, continuous optimization, global optimization, aircraft conflicts

We introduce two new optimization models for the aircraft conflict avoidance problem that aims at issuing decisions on both speed and headingangle deviations to keep aircraft pairwise separated by a given separation distance. The first model is a new mixed-integer nonlinear formulation. The second model is a continuous optimization formulation, less typical in aircraft conflict avoidance. The advantages of the two models are combined within a three-phase methodology that we propose to solve the problem to global optimality. Computational experiments on various instances from the literature yield very promising results, and show the effectiveness of the proposed models and of the three-phase solution approach.

Introduction

We address a real-world application arising in transportation engineering, and more specifically a traffic conflict avoidance problem, where it is essential to maintain vehicles pairwise separated by a threshold distance. The pairwise separation between moving objects arises as a fundamental component in a number of applications spanning from transportation engineering, to robotics, computational geometry, computer graphics, and computer-aided design [START_REF] Basch | Kinetic collision detection between two simple polygons[END_REF]. In transportation engineering, traffic conflict avoidance applies to various modes of transportation; recent applications concern autonomous vehicles [START_REF] Ahn | Safety verification and control for collision avoidance at road intersections[END_REF]. In robotics, motion planning optimization [START_REF] Ecker | Performance of several optimization methods on robot trajectory planning problems[END_REF] must encompass collision avoidance amongs robots. This can be modeled considering a safety zone around each robot, to be avoided [START_REF] Bezzo | A disjunctive programming approach for motion planning of mobile router network[END_REF]. In computational geometry, the problem we consider could be seen as a special case of collision detection of simple 2D geometrical objects with rectilinear motion, as in [START_REF] Cameron | Collision detection by four-dimensional intersection testing[END_REF].

In this paper we consider in particular the context of air traffic management, and we address the aircraft conflict avoidance problem. The objective is to keep vehicles (in our case, aircraft) pairwise separated by a given standard distance all along their trajectories; the loss of such a separation leading to a conflict. Aircraft separation is achieved in practice by modifying the aircraft headings, speeds or altitudes, or sometimes by combining some of these possible maneuvers. Decision variables of the optimization models from the literature are then related to the chosen maneuvers (see [START_REF] Cafieri | MINLP in Air Traffic Management: Aircraft conflict avoidance[END_REF]).

Separation constraints depend continuously on time; the simplest way to deal with them is based on discretizing the time horizon and to impose pairwise separation at each time step. Some approaches discretize the decision variables as well (see [START_REF] Lehouillier | Two decomposition algorithms for solving a minimum weight maximum clique model for the air conflict resolution problem[END_REF] and references therein), while others rely on optimal control [START_REF] Cafieri | Combination of optimal control approaches for aircraft conflict avoidance via velocity regulation[END_REF]. Space discretization (where separation is built upon a set of relevant points of the 2D space) is also proposed in the literature, see e.g., [START_REF] Omer | A space-discretized mixed-integer linear model for air-conflict resolution with speed and heading maneuvers[END_REF], where moreover a discrete set of turning angles is used. Such discretizations potentially yield very large combinatorial optimization problems. Other ways to model separation constraints rely on geometrical constructions or on analytical descriptions of aircraft motion. One such geometrical construction, proposed in [START_REF] Pallottino | Conflict resolution problems for air traffic management systems solved with mixed integer programming[END_REF] and extended later (e.g., [START_REF] Alonso-Ayuso | Collision avoidance in air traffic management: A mixed-integer linear optimization approach[END_REF][START_REF] Alonso-Ayuso | A mixed 0-1 nonlinear optimization model and algorithmic approach for the collision avoidance in ATM: Velocity changes through a time horizon[END_REF]), is based on disks representing safety areas around aircraft (shadows), and on their interception. In this paper, we build on the analytical description of aircraft motion presented in [START_REF] Cafieri | Aircraft deconfliction with speed regulation: New models from mixed-integer optimization[END_REF], and exploited subsequently e.g., in [START_REF] Cafieri | Feasibility pump for aircraft deconfliction with speed regulation[END_REF][START_REF] Cafieri | Mixed-integer nonlinear programming for aircraft conflict avoidance by sequentially applying velocity and heading angle changes[END_REF][START_REF] Cafieri | Maximizing the number of conflict-free aircraft using mixed-integer nonlinear programming[END_REF]. It is tailored to optimization models that do not rely on discretizing either time, or decision variables. The resulting formulation of the separation constraints is however nonlinear.

Most mathematical optimization approaches of the literature involve mixedinteger (either linear or nonlinear) formulations [START_REF] Cafieri | Feasibility pump for aircraft deconfliction with speed regulation[END_REF][START_REF] Cafieri | Aircraft deconfliction with speed regulation: New models from mixed-integer optimization[END_REF][START_REF] Cafieri | Mixed-integer nonlinear programming for aircraft conflict avoidance by sequentially applying velocity and heading angle changes[END_REF][START_REF] Cafieri | Maximizing the number of conflict-free aircraft using mixed-integer nonlinear programming[END_REF][START_REF] Pallottino | Conflict resolution problems for air traffic management systems solved with mixed integer programming[END_REF][START_REF] Rey | Subliminal speed control in air traffic management: Optimization & simulation[END_REF][START_REF] Vela | Near real-time fuel-optimal en route conflict resolution[END_REF]. The only fully-continuous optimization approach, which is based on B-splines to build aircraft conflict-free trajectories, is proposed in [START_REF] Peyronne | Solving air traffic conflict problems via local continuous optimization[END_REF]. It relies on transforming a semi-infinite programming formulation of the pairwise separation constraints into a single equality constraint. Only locally-optimal solutions are obtained; moreover, the authors could not achieve computing explicit derivatives, therefore resorting to finite-difference approximations.

As regards the possible types of separation maneuvers studied in the airtraffic conflict avoidance literature, most works consider either only speed regulation [START_REF] Cafieri | Feasibility pump for aircraft deconfliction with speed regulation[END_REF][START_REF] Cafieri | Aircraft deconfliction with speed regulation: New models from mixed-integer optimization[END_REF][START_REF] Cafieri | Mixed-integer nonlinear programming for aircraft conflict avoidance by sequentially applying velocity and heading angle changes[END_REF][START_REF] Cafieri | Maximizing the number of conflict-free aircraft using mixed-integer nonlinear programming[END_REF][START_REF] Pallottino | Conflict resolution problems for air traffic management systems solved with mixed integer programming[END_REF][START_REF] Rey | Subliminal speed control in air traffic management: Optimization & simulation[END_REF], or only heading angle changes [START_REF] Alonso-Ayuso | Collision avoidance in air traffic management: A mixed-integer linear optimization approach[END_REF][START_REF] Alonso-Ayuso | Exact and approximate solving of the aircraft collision resolution problem via turn changes[END_REF][START_REF] Pallottino | Conflict resolution problems for air traffic management systems solved with mixed integer programming[END_REF][START_REF] Peyronne | Solving air traffic conflict problems via local continuous optimization[END_REF]. Combination of maneuvers is rarely considered, despite the interest of resorting to complementary maneuvers in an operational context. A single type of maneuver may indeed be insufficient to resolve all conflicts. Simply consider the (extreme) situation of face-to-face aircraft: speed regulation cannot evidently be enough. Moreover, it is possible to exhibit problem instances that sole headingangle deviations will not resolve. Combining speed and heading angle deviations is considered in [START_REF] Alonso-Ayuso | Exact and approximate solving of the aircraft collision resolution problem via turn changes[END_REF], leading to a formulation characterized by highly-nonlinear trigonometric constraints. Multiple types of maneuvers are considered in [START_REF] Lehouillier | Two decomposition algorithms for solving a minimum weight maximum clique model for the air conflict resolution problem[END_REF], where the authors propose a maneuver-discretized approach, including altitude (flight level) changes. However, such altitude changes are rarely employed in practice, as they are associated with several operational drawbacks. Another approach based on discretization, combining speed and heading changes, is that of [START_REF] Omer | A space-discretized mixed-integer linear model for air-conflict resolution with speed and heading maneuvers[END_REF]. Speed and heading deviations are applied sequentially in a two-step solution approach in [START_REF] Cafieri | Mixed-integer nonlinear programming for aircraft conflict avoidance by sequentially applying velocity and heading angle changes[END_REF]. In [START_REF] Rey | Complex number formulation and convex relaxations for aircraft conflict resolution[END_REF], the authors propose an approach combining speed and heading angle deviations, that is extended in [START_REF] Dias | Disjunctive linear separation conditions and mixed-integer formulations for aircraft conflict resolution[END_REF] to allow also flight level changes. This approach leads to mixed-integer quadratically-constrained problems, solved through a sequence of convex relaxations. As a result, the authors obtain global optima for several problem instances. In this paper, we consider the same maneuver combination as in [START_REF] Rey | Complex number formulation and convex relaxations for aircraft conflict resolution[END_REF] (speed and heading angle deviations); we propose two new models, allowing us to obtain globally-optimal solutions for instances not solved in [START_REF] Dias | Disjunctive linear separation conditions and mixed-integer formulations for aircraft conflict resolution[END_REF][START_REF] Rey | Complex number formulation and convex relaxations for aircraft conflict resolution[END_REF].

Note finally that, in an operational context, each aircraft needs to recover its original trajectory after being deviated. This subsequent problem is out of the scope of the present paper; let us simply refer the reader to [START_REF] Alonso-Ayuso | Exact and approximate solving of the aircraft collision resolution problem via turn changes[END_REF] for a simple approach to trajectory recovery (it involves solving a convex unconstrained quadratic problem for each pair of aircraft).

The main contribution of this paper is the introduction of two new models for aircraft conflict resolution relying on both aircraft heading and speed deviations, and stemming from two different ways to address the conditional separation constraints characterizing the problem. The inherent combinatorics in conflict avoidance problems leads naturally to conditional constraints: if two vehicles are converging spatially, then their smallest inter-distance must be bounded below at all time. Moreover, the separation constraint that must be imposed to avoid conflicts often involves nonconvex, nonlinear constraints. We introduce a mixed-integer nonlinear programming (MINLP) formulation, and a continuousoptimization model that is based on an original penalty-function strategy introduced in [START_REF] Cafieri | The continuous quadrant penalty formulation of logical constraints[END_REF] to deal with conditional constraints. The MINLP model has a continuous relaxation that is not convex, and the continuous-optimization formulation is a constrained non-convex problem.

As a second contribution, we propose an effective three-phase solution approach whose first phase is the solution of the continuous optimization problem, and the other two phases rely on the MINLP, thus enabling us to obtain exact global solutions.

The remainder of this paper is structured as follows. Section 2 presents the aircraft conflict avoidance problem, and proposes a reformulation of the anglerelated functions involved in the two models we are introducing. Sections 3 and 4 introduce respectively these two models: a mixed-integer nonlinear optimiza-tion formulation, and a continuous optimization model based on a penalty function. Section 5 reports numerical experiments on the aircraft conflict avoidance problem, performed using each of the two models. Section 6 proposes an effective global-optimization approach combining the advantages of the two models proposed above. Finally, conclusions are drawn in Section 7.

Aircraft conflict avoidance: Towards new models

This section first briefly defines the problem to be solved, and explains the important separation constraints, then it introduces some linearizations that will be used in the two new optimization models that we are proposing.

Problem statement and separation constraints

This subsection defines the air traffic conflict avoidance problem, and presents how we model the aircraft separation constraints in a continuous time interval.

In the sequel, for any given vector, x, the notation x i represents its i th component. We assume that we are given an index set, A = {1, 2, . . . , n}, corresponding to aircraft flying at the same constant altitude (flight level), on straight-line segment trajectories. The aim is to keep aircraft pairwise separated during the time horizon these aircraft are monitored. Two aircraft are separated if their relative distance is greater than or equal to a given standard separation distance (commonly, 5 nautical miles (NM) for commercial en-route flights). A conflict corresponds to the loss of such a separation. For each aircraft i ∈ A, we know its initial position, (x 0 i , y 0 i ) ∈ R 2 , and its velocity (i.e., its heading angle, ϕ i , and its speed, v i ). Given this configuration, aircraft conflict avoidance consists in issuing separation maneuvers (simultaneously, at time t = 0) so as to ensure that the aircraft will remain separated.

The main difficulty in aircraft conflict avoidance problems comes from the separation constraint that must hold for each pair i, j ∈ A (i < j) of aircraft:

||x ij (t)|| ≥ d ∀ t ≥ 0, ( 1
)
where d is the minimum required separation distance, x ij (t) = x i (t)-x j (t) ∈ R 2 is a vector representing the position of aircraft i relative to that of aircraft j, and we use the usual, Euclidean norm. This constraint represents in general a non-convex set. A further difficulty is that it must hold for all time t. As mentioned in Section 1, the latter difficulty has been addressed in various ways in the literature. Here we choose the simple "parabola approach" of Cafieri et al. [START_REF] Cafieri | Aircraft deconfliction with speed regulation: New models from mixed-integer optimization[END_REF] to address this continuum of constraints. Roughly speaking, one remarks that, for each pair of aircraft i and j, satisfying constraint (1) amounts to ensuring that a simple univariate quadratic separation function, f ij (t), is always non-negative. More precisely: where v ij is the velocity vector of aircraft i relative to that of aircraft j, and x 0 ij is the vector of their relative initial positions. These vectors will in turn be straightforwardly computable explicitly in terms of the decision variables and of the input data of the problem. We assume that aircraft are separated at t = 0 (verified in a pre-processing stage). Without loss of generality, we suppose v ij ̸ = 01 . Since the quadratic function f ij is strictly convex, one can easily verify that the condition (1) is equivalent to:

f ij (t) = ∥v ij ∥ 2 t 2 + 2x 0 ij • v ij t + ∥x 0 ij ∥ 2 -d 2 ≥ 0 ∀ t ≥ 0, (2)
if t m ij > 0, then f ij (t m ij ) ≥ 0, (3)
where

t m ij is the minimizer of f ij . As t m ij = -x 0 ij • v ij ∥v ij ∥ 2 ,
finally the separation constraint (1) is equivalent to the following logical condition:

if t m ij > 0, then ∥v ij ∥ 2 (∥x 0 ij ∥ 2 -d 2 ) -(x 0 ij • v ij ) 2 ≥ 0. ( 4
)
Figure 1 provides an intuitive interpretation of the separation condition (3): the parabola gives, at every time t, the value of the separation function f ij (which is positive if and only if aircraft i and j are separated by a distance of at least d) i.e., it gives the square of the distance between aircraft i and j plus d 2 , at any time t. One observes that, in the displayed example, initially (at time t = 0) the aircraft are separated (the curve goes through the green vertical half-line) but are converging (their inter-distance is decreasing). The separation constraint therefore requires the "lowest" point of the parabola, p m := (t m ij , f ij (t m ij )), to lie outside the fourth quadrant.

Towards new optimization models

We aim at proposing mathematical models involving decision variables controlling both speed and heading angle of each aircraft. We first present our main decision variables and the (nonlinear) way they are related to express the separation conditions. Then, we propose some reformulations to avoid explicitly handling trigonometric functions.

For each aircraft i ∈ A, let q i be its speed variation, and let θ i be its headingangle variation. More precisely, q i is the proportion of variation of the aircraft original speed, v i , so that the actual aircraft speed is q i v i , and θ i is added (its value may be negative) to the original angle, ϕ i . Both speed and angle variations are assumed to be bounded. If we suppose that the speed can be decreased by no more than 6% and increased by no more than 3%, (common operational and cost constraints are even stricter), we have that the actual aircraft speed is

q i v i ∈ [0.94v i , 1.03v i ].
We shall assume that the heading angles, θ i , are bounded in the interval [θ i , θ i ], where θ i = -π/6 and θ i = π/6, as it is common practice in air traffic control for operational reasons.

Recall from Subsection 2.1, that in order to express the separation condition for a pair of aircraft i, j ∈ A (i < j), we suppose that their relative initial position,

x 0 ij ∈ R 2 , is given: x 0 ij := x 0 i -x 0 j y 0 i -y 0 j ,
and, taking into account possible speed and heading angle variations, their relative velocity, v ij ∈ R 2 , is:

v ij := cos(ϕ i + θ i )(q i v i ) -cos(ϕ j + θ j )(q j v j ) sin(ϕ i + θ i )(q i v i ) -sin(ϕ j + θ j )(q j v j ) . (5)
To avoid explicit trigonometric functions in our optimization models, we define for each aircraft i ∈ A:

c i := cos(ϕ i + θ i ),
and

s i := sin(ϕ i + θ i ). (6)
In order to linearize (5), we further define ω i and π i for each i ∈ A as follows:

ω i := c i q i v i (7)
π i := s i q i v i , (8)
so that

v ij = ω i -ω j π i -π j , i, j ∈ A : i < j. (9)
In both optimization models to be introduced in Sections 3 and 4, our main decision-variable vectors will be: ω, π and q, (and v remains as an auxiliary decision-variable vector). Let us consider a particular aircraft i ∈ A. The variable q i controls the variation of the speed of aircraft i. The other two variables, ω i and π i , although they both depend on q i , account for variations in the heading angle of aircraft i. Note that, as c 2 i + s 2 i = 1, it follows that:

ω 2 i + π 2 i = (q i v i ) 2 (c 2 i + s 2 i ) = (q i v i ) 2 . ( 10
)
Thus, in fact we only have two degrees of freedom for each aircraft i ∈ A. Equation ( 10) is a nonlinear equality constraint that is included in our optimization models for each i ∈ A.

Remark: Going from the pair of variables (q i , θ i ) to the triple (q i , ω i , π i ) is crucial. As already mentioned, it still involves only two degrees of freedom for each aircraft i ∈ A, taking into account equation [START_REF] Cafieri | Combination of optimal control approaches for aircraft conflict avoidance via velocity regulation[END_REF]. Expressing our optimization problem in terms of the above triple has the strong advantage of avoiding trigonometric functions, and of removing the bilinear product c i q i and s i q i , thus linearizing ( 5) into ( 9), thereby enabling a more efficient use of optimization solvers.

We introduce in the following two sections two new models for aircraft conflict avoidance based on the above setting: a mixed-integer non-linear formulation, and a continuous penalty-based nonlinear model. Both formulations are constructed using the main optimization variables, q i , ω i and π i , i ∈ A, and the quadratic equality constraints [START_REF] Cafieri | Combination of optimal control approaches for aircraft conflict avoidance via velocity regulation[END_REF] linking the above variables (together with the auxiliary variables v ij , i, j ∈ A (i < j), and their defining constraints). The two formulations address differently the crucial issue of the logical separation condition:

if t m ij > 0, then ∥v ij ∥ 2 (∥x 0 ij ∥ 2 -d 2 ) -(x 0 ij • v ij ) 2 ≥ 0 ( 11
)
that one must consider for every possible pair of aircraft, i, j ∈ A, i < j.

MINLP formulation

The first model that we introduce for aircraft conflict avoidance is an MINLP, in which the separation constraint ( 11) is modelled via the complementary formulation (see, e.g., [START_REF] Belotti | On handling indicator constraints in mixed integer programming[END_REF]) of logical constraints:

t m ij (y ij -1) ≥ 0, ( 12
)
y ij ∥v ij ∥ 2 ∥x 0 ij ∥ 2 -d 2 -(x 0 ij • v ij ) 2 ≥ 0, ( 13
)
y ij ∈ {0, 1}, (14)
where y ij is an extra auxiliary binary decision variable. Remark by the way that constraint [START_REF] Cafieri | The continuous quadrant penalty formulation of logical constraints[END_REF] is equivalent to the following implication (without a strict inequality in the condition on the sign of t m ij ):

if t m ij ≥ 0, then ∥v ij ∥ 2 (∥x 0 ij ∥ 2 -d 2 ) -(x 0 ij • v ij ) 2 ≥ 0, (15)
since by assumption, aircraft i and j are separated when t m ij = 0 (it is rather the form (15) that is seen in the literature on aircraft conflict avoidance).

In order to formulate our optimization model, it remains to define an appropriate objective function. A natural objective is to minimize the angle and speed variations with respect to the given initial conditions ϕ i and v i (i ∈ A). However, this requires expliciting each angle change in the objective function through trigonometric functions, losing thereby the computational benefit of our model for angle change. Note first that, as in air-traffic practice the allowed angle deviation, θ i , is a bounded variable, one has for each i ∈ A:

a ≤ cos(θ i ) ≤ 1 and -b ≤ sin(θ i ) ≤ b, (16)
for some values of a ≥ 0 and b ≤ 1. In practice, each |θ i | is typically bounded by π/6; so that one generally sets a = √ 3

2 and b = 1 2 . Let us now introduce, for each i ∈ A, an auxiliary variable, b i , defined as an upper bound on | sin(θ i )|:

-b ≤ -b i ≤ sin(θ i ) ≤ b i ≤ b. ( 17
)
In order to minimize angle deviations, one then aims at minimizing the b i 's. Thus, in order to minimize both speed and angle deviations, one may consider the following objective function:

(1 -λ) i∈A (1 -q i ) 2 + λ i∈A b i , (18)
where λ ∈ [0, 1] is a user-defined weighting parameter that can be adjusted to reach a good compromise between speed and angle deviations.

Linking constraints

Introducing the auxiliary-variable vector b requires defining constraints to link these new variables to the original decision-variable vectors q, ω and π. To that aim, let i ∈ A, and let us consider the trigonometric identity:

cos(ϕ i + θ i ) = cos ϕ i cos θ i -sin ϕ i sin θ i , together with a ≤ cos(θ i ) ≤ 1 and -b i ≤ sin(θ i ) ≤ b i .
Considering the four cases that can arise, depending on the signs of sin ϕ i and cos ϕ i , and recalling that each (given) heading angle, ϕ i , satisfies -π ≤ ϕ i < π, one obtains the following relationship:

min (a cos ϕ i , cos ϕ i )-b i | sin ϕ i | ≤ c i ≤ max (a cos ϕ i , cos ϕ i )+b i | sin ϕ i | (19)
Similarly, the trigonometric identity: [START_REF] Fourer | AMPL: A Modeling Language for Mathematical Programming[END_REF] where c i and s i are the intermediate cosine and sine variables defined by [START_REF] Belotti | On handling indicator constraints in mixed integer programming[END_REF]. From the definitions ( 7) and ( 8) of ω i and π i , we then obtain the following linking constraints:

sin(ϕ i + θ i ) = sin ϕ i cos θ i + cos ϕ i sin θ i yields: min (a sin ϕ i , sin ϕ i )-b i | cos ϕ i | ≤ s i ≤ max (a sin ϕ i , sin ϕ i )+b i | cos ϕ i |,
ω i ≥ (min (a cos ϕ i , cos ϕ i ) -b i | sin ϕ i |) q i v i ( 21
)
ω i ≤ (max (a cos ϕ i , cos ϕ i ) + b i | sin ϕ i |) q i v i ( 22
)
π i ≥ (min (a sin ϕ i , sin ϕ i ) -b i | cos ϕ i |) q i v i ( 23
)
π i ≤ (max (a sin ϕ i , sin ϕ i ) + b i | cos ϕ i |) q i v i . (24)
Valid (bound-constraint) inequalities Let us now determine lower and upper bounds for the decision variables ω i and π i , i ∈ A, in order to render a branch-and-bound solution procedure more efficient. Let us introduce the following notation, related to inequalities [START_REF] Ecker | Performance of several optimization methods on robot trajectory planning problems[END_REF] and ( 20):

c i (b i ) := min (a cos ϕ i , cos ϕ i ) -b i | sin ϕ i | (25)
c i (b i ) := max (a cos ϕ i , cos ϕ i ) + b i | sin ϕ i | (26)
s i (b i ) := min (a sin ϕ i , sin ϕ i ) -b i | cos ϕ i | (27)
s i (b i ) := max (a sin ϕ i , sin ϕ i ) + b i | cos ϕ i |. ( 28
)
Using the bounds on the variable q i , q i ≤ q i ≤ q i (where, as mentioned above, commonly q i = 0.94 and q i = 1.03), and following the possible cases for the signs of (the given data) cos ϕ i and sin ϕ i in the expressions of c i (b i ) and c i (b i ) above, we obtain, for each aircraft i ∈ A, the bound-constraint valid inequality:

ω i ≤ ω i ≤ ω i , (29)
where

ω i = c i ( b) q i v i , if c i ( b) ≥ 0 c i ( b) q i v i , otherwise,
and

ω i = c i ( b) q i v i , if c i ( b) ≤ 0 c i ( b) q i v i , otherwise.
Recall that v i , c i ( b) and c i ( b) are constants (the latter will be denoted simply c i and c i respectively in the sequel), and the speed bounds, q i and q i , are always positive.

In an analogous manner, we obtain for π i :

π i ≤ π i ≤ π i , (30)
where The MINLP formulation can be summarized as follows: min q,ω,π,t m ,v,b,y

π i = s i q i v i , if s i ≥ 0 s i q i v i , otherwise, and
π i = s i q i v i , if s i ≤ 0 s i q i v i ,
(1 -λ) i∈A (1 -q i ) 2 + λ i∈A b i (31)
s.t.

y ij ∥v ij ∥ 2 ∥x 0 ij ∥ 2 -d 2 -(x 0 ij • v ij ) 2 ≥ 0 i, j ∈ A : i < j (32) t m ij ∥v ij ∥ 2 = -x 0 ij • v ij i, j ∈ A : i < j (33) t m ij (y ij -1) ≥ 0 i, j ∈ A : i < j (34) v ij = ω i -ω j π i -π j i, j ∈ A : i < j (35) ω 2 i + π 2 i = (q i v i ) 2 i ∈ A ( 36
)
ω i ≥ (min (a cos ϕ i , cos ϕ i ) -b i | sin ϕ i |) q i v i i ∈ A ( 37
)
ω i ≤ (max (a cos ϕ i , cos ϕ i ) + b i | sin ϕ i |) q i v i i ∈ A ( 38
)
π i ≥ (min (a sin ϕ i , sin ϕ i ) -b i | cos ϕ i |) q i v i i ∈ A ( 39
)
π i ≤ (max (a sin ϕ i , sin ϕ i ) + b i | cos ϕ i |) q i v i i ∈ A ( 40
)
q i ≤ q i ≤ q i i ∈ A (41)
ω i ≤ ω i ≤ ω i i ∈ A (42)
π i ≤ π i ≤ π i i ∈ A (43) 0 ≤ b i ≤ b i ∈ A ( 44
)
y ij ∈ {0, 1} i, j ∈ A : i < j. (45)
Remarks:

1. Constraints (37) to (40) are bilinear constraints since, as cos(ϕ i ) and sin(ϕ i ) are known input data, their sign, and thereby the maxima and minima involved, are known in advance.

2. Replacing constraints (33) and (34) with the equivalent constraints:

(x 0 ij • v ij )(y ij -1) ≤ 0 i, j ∈ A : i < j, (46)
one can eliminate the optimization variables t m ij 's and the highly non-linear constraints (33) from the formulation.

3. In our implementation, we do not use the decision variables v ij , but rather auxiliary variables representing ∥v ij ∥ 2 and x 0 ij • v ij , i, j ∈ A : i < j, and we adapt thereby constraints (35).

Penalty continuous-optimization formulation

This section introduces our second optimization formulation, named NLP-Penalty in the sequel, of the aircraft conflict avoidance problem. The construction of the NLP-Penalty model differs from that of our MINLP formulation in the way the logical separation constraint [START_REF] Cafieri | The continuous quadrant penalty formulation of logical constraints[END_REF] is addressed. Whereas the MINLP formulation models this logical constraint via the complementary formulation, the NLP-Penalty model addresses it via the continuous quadrant penalty formulation of logical constraints introduced in [START_REF] Cafieri | The continuous quadrant penalty formulation of logical constraints[END_REF].

The continuous optimization formulation we are proposing involves the optimization variables:

q i , ω i , π i , i ∈ A and (t m ij , f m ij ), v ij , i, j ∈ A : i < j, (47)
where t m ij = argminf ij (t), using the notation (2), and letting

f m ij := f ij (t m ij ). The corresponding feasible domain is a subset of R n × R n × R n × (R 2 \ S) ( n 2 ) × R ( n 2 ) . ( 48
)
Letting i, j be a pair of aircraft (i, j ∈ A, i < j), remark that the separation condition [START_REF] Cafieri | The continuous quadrant penalty formulation of logical constraints[END_REF] then simply reads:

if t m ij > 0, then f m ij ≥ 0. ( 49
)
We reformulate the logical constraint (49) using the continuous quadrant penalty function g β : (t, f ) ∈ R 2 → R introduced in [START_REF] Cafieri | The continuous quadrant penalty formulation of logical constraints[END_REF]:

g β (t, f ) =                    0 if t ≤ 0 or f ≥ 0 t 2 if 0 < t ≤ - f 3 -(t 2 + 6tf + f 2 ) 8 if - f 3 < t < -3f f 2 if - t 3 ≤ f < 0 (50)
(setting β = 3 as proposed by the authors), and illustrated on Figure 2. Summing up over all pairs i, j ∈ A : i < j, this yields the penalty function:

i,j∈A:i<j

g β (t m ij , f m ij ). (51)
Each of the term functions g β is shown in [START_REF] Cafieri | The continuous quadrant penalty formulation of logical constraints[END_REF] to be smooth. To give a geometrical interpretation of the objective function (51), remark that the logical constraint (49) is equivalent to requiring: (t m ij , f m ij ) ∈ R 2 \ S (where S is the "forbidden" open fourth quadrant). Since v ij , t m ij , and f m ij depend on q i , ω i , π i and q j , ω j , π j , these six optimization variables control the location of the point p m = (t m ij , f m ij ) on the graph of f ij (Figure 1). Each term of the objective function (51) therefore guides a continuous optimization descent method towards the subset of the feasible set where the corresponding separation constraint is satisfied. To summarize, the (continuous) NLP-Penalty formulation we are proposing penalizes constraint (49), i, j ∈ A : i < j, and preserves the other constraints:

t f 0 0 0 f 2 -1 8 (t 2 + 6 t f + f 2 ) t 2
min q,ω,π,t m ,f m ,v i,j∈A:i<j g β (t m ij , f m ij ) (52) s.t.
f m ij ∥v ij ∥ 2 = ∥v ij ∥ 2 (∥x 0 ij ∥ 2 -d 2 ) -(x 0 ij • v ij ) 2 i, j ∈ A : i < j (53) t m ij ∥v ij ∥ 2 = -x 0 ij • v ij i, j ∈ A : i < j (54) v ij = ω i -ω j π i -π j i, j ∈ A : i < j (55) ω 2 i + π 2 i = (q i v i ) 2 i ∈ A ( 56
)
q i ≤ q i ≤ q i i ∈ A ω i ≤ ω i ≤ ω i i ∈ A (57) π i ≤ π i ≤ π i i ∈ A. (58)
Remarks:

1. In our implementation, the f m ij 's and the t m ij 's do not in fact appear as auxiliary variables: they are explicitly replaced in the model using constraints (53) and (54). As for the MINLP model, we also do not use the decision variables v ij (but rather auxiliary variables representing ∥v ij ∥ 2 and x 0 ij • v ij , i, j ∈ A : i < j), and we adapt thereby constraints (55).

2. Since the objective function of the NLP-penalty formulation does not involve angle-related deviations, the vector, b, of extra decision variables introduced for the MINLP does not intervene in the calculation of the bound constraints (57) and ( 58), but only its extreme values through b, the upper bound on the | sin(θ i )|'s.

As demonstrated in [START_REF] Cafieri | The continuous quadrant penalty formulation of logical constraints[END_REF], the penalty function, g β , features several desirable properties in view of its use with a state-of-the-art local optimization solver. However, the fact that each term function, g β : R 2 → R, leans outwards the (open fourth quadrant) set S (which corresponds to a non-separated pair of aircraft) does not imply that the objective function i,j∈A:i<j

g β (t m ij , f m ij ) : R 2( n 2 ) → R leans outwards the Cartesian-product set S ( n 2 ) (⊆ R 2( n 2
) ). Indeed, remark first that this objective function is not separable: each of the speed and anglerelated variables q i , ω i , π i corresponding to one aircraft i ∈ A intervenes in several terms of the objective function (each aircraft i is involved in several pairs of aircraft). Moreover, this objective function is not convex, and may thereby feature local optima that do not satisfy the constraint(s) the penalty function aims at modelling. Finally, one must also take into account the remaining constraints (all the ones that are not penalized), which include nonlinear equality constraints such as (56), yielding thereby a nonconvex feasible domain.

The next two sections show that good results are obtained despite the abovementionned potential difficulties.

Numerical results

This section presents the performance of the MINLP and the NLP-Penalty formulations of Sections 3 and 4, through numerical experiments.

We test our models on two types of benchmark instances: the Circle Problem (CP), and the Random Circle Problem (RCP). In both cases, instances are constructed as follows. Aircraft are uniformly positioned on the circumference of a circle, representing the observed airspace portion, and are given an initial speed and a heading to fly towards the opposite side of the circle. This kind of instances are widely used in the literature [START_REF] Cafieri | Aircraft deconfliction with speed regulation: New models from mixed-integer optimization[END_REF][START_REF] Cafieri | Mixed-integer nonlinear programming for aircraft conflict avoidance by sequentially applying velocity and heading angle changes[END_REF][START_REF] Cafieri | Maximizing the number of conflict-free aircraft using mixed-integer nonlinear programming[END_REF][START_REF] Dias | Disjunctive linear separation conditions and mixed-integer formulations for aircraft conflict resolution[END_REF][START_REF] Rey | Complex number formulation and convex relaxations for aircraft conflict resolution[END_REF]. We use the instances proposed by (and available via) [START_REF] Dias | Disjunctive linear separation conditions and mixed-integer formulations for aircraft conflict resolution[END_REF][START_REF] Rey | Complex number formulation and convex relaxations for aircraft conflict resolution[END_REF]. The radius of the circle is 200 nautical miles (NM). The CP instances are generated with initial identical speeds (v i = 500 NM/h, i ∈ A), and with all aircraft heading towards the center of the circle, thus resulting in a highly symmetric problem. The RCP instances are such that aircraft initial speeds and headings are randomly deviated within specified ranges (486-594 NM/h for the speed, and ±π/6 with respect to the heading towards the center, for the angle), resulting thereby in more realistic and less structured instances. The standard separation distance, d, to be respected between aircraft, is equal to 5 NM. Aircraft speeds and heading angles are bounded due to operational constraints. For each aircraft i ∈ A, we set the bounds on variables q i to: q i = 0.94 and qi = 1.03, that correspond to a socalled subliminal speed change (the original aircraft speed is allowed to vary only between -6% and +3%). We allow heading angles, θ i , to vary between -π/6 and +π/6, and compute bounds on variables ω i and π i accordingly (i ∈ A). The β parameter of the continuous quadrant penalty function, g β , is set to β = 3.

The proposed models are implemented using the AMPL [START_REF] Fourer | AMPL: A Modeling Language for Mathematical Programming[END_REF] modeling language. Numerical tests are performed on a 2.66 GHz Intel Xeon (octo core) processor with 32 GB of RAM and Linux operating system.

We first present numerical results obtained by considering our MINLP model (Section 3). The problem is solved by running the MINLP solver COUENNE [START_REF] Belotti | Branching and bounds tightening techniques for non-convex MINLP[END_REF], version 0.5, with all its default values, except for its feasibility tolerance which we increase to 10 -5 . As this is recommended only for small problems, we disable both the optimality-based bound tightening as well as the "aggressive" feasibility-based bound tightening options, that are generally computationally expensive. Tests are run with a time limit of 600 seconds.

Note that, in contrast with [START_REF] Rey | Complex number formulation and convex relaxations for aircraft conflict resolution[END_REF], whose objective function involves a fixed, arbitrary nonlinear combination of the speed and angle deviations, our model allows the user to decide a particular speed-angle trade off, expressed through the weighting parameter, λ in (31), or to compensate any unfortunate choice of units, for instance to scale the angle deviation (expressed in radiants, degrees, etc.) with the speed deviation (which can be in percentage as here, in mile per hour, or in knots). In our experiments, we set λ = 10 -6 . The recourse to such a small value of λ is not surprising when one considers the difference in the scales of the two criteria, speed and angle deviations, that constitute the objective function [START_REF] Duff | MA57-A code for the solution of sparse symmetric definite and indefinite systems[END_REF]. Indeed, a typical extreme speed deviation for an aircraft i corresponds to accelerating it up to a fraction q i = 1.03 of its initial speed, which yields a corresponding penalty term value of (1 -q i ) 2 = 9 10,000 ; whereas an extreme angle deviation involves a heading angle change of |θ i | = π 6 which yields a corresponding penalty term value of b i = 1 2 : the latter is 3 to 4 orders of magnitude larger than the former.

Table 1 reports the results obtained solving our MINLP problem. The headings of the table are as follows: Name, name of the instance; n, number of aircraft; n c , number of initial conflicts (conflicts resulting from no speed change and no heading deviation); n hthc , number of initial head-to-head conflicts; time, computing time in seconds, and speed dev., value of the speed deviation in the computed optimal solution (this later value will be useful in the tests of Section 6). One observes that our MINLP model allows us to find global optima for all the classical CP instances except for the largest one (n=20), including the ones involving 11 ≤ n ≤ 19 aircraft for which only local optima could be found in [START_REF] Dias | Disjunctive linear separation conditions and mixed-integer formulations for aircraft conflict resolution[END_REF][START_REF] Rey | Complex number formulation and convex relaxations for aircraft conflict resolution[END_REF]. The computing time is always less than 4 seconds for these academic and highly symmetric instances, even when a large number of conflicts are to be solved. Moreover, all the 10-and 20-aircraft RCP instances are solved to global optimality, the computing time varying from 0 to a little more than 5 seconds (note that among the randomly-generated instances of [START_REF] Rey | Complex number formulation and convex relaxations for aircraft conflict resolution[END_REF], RCP 10 10 features no conflicts to be solved). For the, larger, 30-aircraft RCP instances, 7 out of 15 instances are solved to global optimality, while our MINLP ends with a time limit (tlim in Table 1) for the remaining cases. In such cases, a feasible solution is however always computed (no guarantee of global optimality).

Let us now concentrate on the solution of the asymmetrical, and more realistic, RCP instances using our NLP-Penalty model. We choose the IPOPT solver for NLPs [START_REF] Wächter | On the implementation of an interior-point filter line-search algorithm for large-scale nonlinear programming[END_REF], version 3.12, to solve the optimization problem. The inner linear algebra solver is MA57 [START_REF] Duff | MA57-A code for the solution of sparse symmetric definite and indefinite systems[END_REF], and we set the feasibility tolerance to 10 -5 , exactly as for the MINLP solver. Note that if this local optimization approach attains the lower bound zero, then global optimality is achieved with respect to the penalty objective function (51). Such solutions correspond however only to feasible solutions for the MINLP problem (which aims at deciding angle and speed for each aircraft so as to minimize their deviations from original values). Recall moreover that our NLP-Penalty problem is a nonconvex constrained NLP; it can thereby feature multiple local minima. Further, the penalization approach is likely to yield local minima that violate the (penalized) separation constraints (referred to as local infeasibility). The strategy we propose to address these issues is to solve the problem (to local optimality), using randomly-generated starting points through a simple multistart heuristic. We first run the solver from the starting guess corresponding to no deviation in speed and no deviation in heading angle (i.e., setting for all i ∈ A: q i = 1, and values of ω i and π i corresponding to θ i = 0). Then, within a maximum of N = 5 trials, we generate new starting guesses choosing randomly, for each i ∈ A: q i uniformly in [0.94, 1.03], and ω i , π i such that the value of θ i varies uniformly in [-π/6, π/6]. Table 2 reports the numerical results obtained using our NLP-Penalty model. The first four columns report the characteristics of the test instances as for Table 1; the following four columns report respectively the computing time, in seconds, for the successful run (time success); the total time, in seconds, obtained by summing up the computing times of all trials (total time); the number of trials (nb trials); and the speed deviation that corresponds to the computed solution (speed dev.). The later is computed a posteriori. One first observes that this local-optimization approach always yields, on all the considered test cases, global optimality with respect to the penalty objective function (51) (since the lower bound zero is attained). Moreover, for 32 instances out of 34 the solution is computed in one run only, while for the other ones no more than 2 runs are needed. Computing times are always below 6 seconds, except for one instance whose multiple starting-point runs yielded a total of 20 seconds. These results show the interest of the proposed continuous optimization approach.

Three-phase solution method

The proposed NLP-Penalty and MINLP formulations yield solutions that are not comparable: the former yields a locally-optimal solution while the latter ends with a global minimizer and this, for objective functions that are different. In this subsection, we propose an effective global optimization approach that attempts at exploiting the advantages of both models. It is a three-phase methodology whose first phase solves the NLP-Penalty problem; the other two phases rely on the MINLP to obtain an exact global minimizer.

The solution method we are proposing applies to the case of the minimization of aircraft speed changes only (heading deviation changes remain however as optimization variables). This choice is based on the observation that, in practice, air traffic controllers (ATC) prefer handling aircraft heading deviation changes rather than (subliminal) speed changes which cannot be easily visualized on air situation displays. Thus, we consider λ = 0 in our MINLP model. This yields an optimization problem where decisions are taken on both speed and heading angle changes, and the objective is to minimize i∈A (1

-q i ) 2 .
Here is a summary of our approach : First, we solve the NLP-Penalty problem by a state-of-the-art local optimization solver in order to find a conflict-free solution. In an attempt to avoid convergence towards a point of local infeasibility, the same simple multistart heuristic described above is applied. From the obtained solution, we then compute the corresponding (angle-related) values of b i , for all aircraft i ∈ A, and of the binary variables y ij , for all pairs i, j, i < j (recall that the b i 's and y ij 's are not optimization variables of the NLP-Penalty problem). We also compute, a posteriori, the value of the speed deviation i∈A (1 -q i ) 2 . A speed deviation that is zero up to a small tolerance, denoted tol and set to tol = 10 -6 , allows us to stop the process: the computed feasible solution is globally optimal.

Second, we fix the variables y ij for all i, j, i < j, to the values computed in the first phase, and solve the problem (referred to, in the following, as MINLP f ix ) starting from the initial solution corresponding to the values of q i , ω i , π i , b i and y ij computed in the first phase. The value of the speed deviation computed at the end of the first phase is used as an upper bound (cutoff value) in the branchand-bound addressing MINLP f ix , in order to remove sub-optimal regions from the feasible set. This second phase consists therefore in a global optimization on a continuous nonconvex problem. If the computed solution corresponds to a null (up to the tolerance tol) speed deviation (which is the objective function to be minimized, since we set λ = 0), then it is a globally optimal solution, and we stop.

Third, if the speed deviation is positive (larger than tol), the value of the binary variables y ij are unfixed, and the MINLP is solved, starting from the last computed solution, and using the corresponding speed deviation as an upper bound.

This three-phase method is described in Algorithm 1.

We apply this three-phase solution method on the 30-aircraft RCP instances. We use again IPOPT to solve the NLP-Penalty problem, and COUENNE to solve the problems MINLP f ix and MINLP. The solvers' setting is the same as in the previous tests, with again a time limit of 600 seconds. The values of the other parameters involved in our models are also the same as in the previous tests.

Table 3 reports the results. The first four columns are the same as in previous tables, followed by two columns for each of the second and the third phases, reporting the computing time in seconds, and the corresponding speed deviations. The last column displays the total time in seconds, i.e., the sum of the computing time for solving the problem through the three phases of computation. The computing time for the first phase is the one reported in Table 2 and is therefore not repeated in Table 3. A dash ("-") in the 3rd-phase columns indicates that global optimality was already proved when exiting the 2nd phase.

One observes that all the instances, except RCP 30 5, are solved to global optimality. At the exception of three instances, one does not need to go beyond Algorithm 1 Aircraft conflict avoidance: Three-phase method Require: n, number of aircraft;

for each aircraft i: v i , initial speed; ϕ i initial heading; (x 0 i , y 0 i ) initial position tol, speed-deviation stopping criterion; N , number of local-opt. starting points -If (no initial pairwise conflicts), then stop -Current iterate: q c = 1, and ω c , π c such that θ = 0; upper bound := +∞ Phase 1:

-Solve NLP-Penalty by local optimization (multistart N times) -If (zero penalty value), then update (q c , ω c , π c ) {conflict-free solu-tion} else {local infeasibility} re-initialize q c = 1, and ω c , π c such that θ = 0; compute b c , y c , and go to Phase 3 -Compute b c , y c , and q dev := i∈A (1 -q c i ) 2

-If (q dev ≤ tol) then stop, else upper bound := q dev Phase 2:

-Solve MINLP f ix by continuous global optimization, using upper bound starting from (q c , ω c , π c , b c , y c ), with λ = 0 {minimizing the speed deviation} to get new (q c , ω c , π c , b c , y c ). Compute q dev := i∈A (1 -q c i ) 2

-If (q dev ≤ tol) then stop, else upper bound := q dev Phase 3:

-Solve MINLP by mixed-integer global optimization, using upper bound starting from (q c , ω c , π c , b c , y c ), with λ = 0 {minimizing the speed deviation} the second phase to obtain an optimal solution. For all tests, the speed deviation computed through the second phase is at most 3.3 ×10 -6 . The third phase is run for the three instances for which the speed deviation is not sufficiently small (larger than tol). In two cases out of these three instances, the third phase then reduces the speed deviation to 10 -17 , whereas for the remaining instance the time limit is attained. The total computing time to obtain a globally-optimal solution is lower than one minute for 8 out of the 15 instances. This is very promising in view of the aircraft conflict avoidance application.

Conclusion

This paper proposes two new mathematical optimization models for the airtraffic conflict avoidance problem involving finding aircraft speed and heading (angle) deviations. Both models benefit from reformulations to remove nonlinearities that render large instances difficult to tackle. One model is a mixedinteger nonlinear optimization formulation. The other model constitutes one of the only two continuous-optimization formulations of the aircraft conflict avoidance problem; it relies on the continuous quadrant penalty formulation of logical constraints introduced in [START_REF] Cafieri | The continuous quadrant penalty formulation of logical constraints[END_REF].

We propose moreover a three-phase global-optimization methodology combining the advantages of both models, that allows us to address efficiently aircraft conflict avoidance problems.

Future tracks of research include traffic conflict problems dealing with autonomous vehicles involving for instance separation maneuvers applied not simultaneously, in a dynamic environment.

Figure 1 :

 1 Figure 1: The separation function for two initially-separated, converging aircraft

  otherwise, with the notation: s i := s i ( b), and s i := s i ( b).

Figure 2 :

 2 Figure 2: The piecewise quadratic penalty function g β (left), and its 3D graph (right), when β = 3 (from [11]).

Table 1 :

 1 Results on CP and RCP instances: MINLP, with λ = 10 -6 . Global exact solutions.

								MINLP (λ = 10 -6 )
					Name	n nc n hth time (s) speed dev.
					RCP 10 1 10 2	1	0.368	6.1e-12
					RCP 10 2 10 3	1	0.276	4.6e-12
					RCP 10 3 10 2	1	0.276	7.4e-13
					RCP 10 4 10 1	1	1.172	1.2e-12
					RCP 10 5 10 5	1	0.268	3.6e-11
					RCP 10 6 10 4	1	0.488	5.6e-12
	MINLP (λ = 10 -6 ) Name n nc n hth time (s) speed dev.	RCP 10 7 10 4 RCP 10 8 10 4 RCP 10 9 10 3	1 1 1	0.324 0.512 0.624	4.2e-12 8.4e-12 4.7e-12
	CP 4 4 6	2	0.056	3.4e-18	RCP 10 10 10 0	0	0	0
	CP 5 5 10	0	0.080	4.5e-15	RCP 20 1 20 8	2	2.320	1.1e-09
	CP 6 6 15	3	0.232	4.5e-15	RCP 20 2 20 9	2	2.244	2.0e-10
	CP 7 7 21	0	0.256	4.5e-11	RCP 20 3 20 13	2	3.860	3.5e-10
	CP 8 8 28	4	0.248	8.5e-12	RCP 20 4 20 9	2	5.712	8.9e-06
	CP 9 9 36	0	0.224	4.1e-11	RCP 20 5 20 12	3	2.264	5.5e-11
	CP 10 10 45	5	0.408	2.5e-10	RCP 20 6 20 13	2	2.840	1.0e-11
	CP 11 11 55	0	0.420	9.2e-13	RCP 20 7 20 9	2	2.476	2.8e-10
	CP 12 12 66	6	0.372	4.0e-10	RCP 20 8 20 9	2	2.552	2.9e-11
	CP 13 13 78	0	1.096	2.8e-11	RCP 20 9 20 19	4	2.100	7.4e-11
	CP 14 14 91	7	1.012	1.9e-11	RCP 20 10 20 15	3	2.088	1.3e-10
	CP 15 15 105 0	2.764	1.4e-09	RCP 30 1 30 35	1	25.02	5.1e-09
	CP 16 16 120 8	1.700	5.1e-10	RCP 30 2 30 38	1	11.59	3.1e-09
	CP 17 17 136 0	2.304	8.1e-10	RCP 30 3 30 46	1	47.13	8.5e-06
	CP 18 18 153 9	2.788	4.7e-10	RCP 30 4 30 39	1	301.6	0.000214
	CP 19 19 171 0	3.672	2.7e-10	RCP 30 5 30 36	2	10.31	3.6e-10
	CP 20 20 190 10	tlim	0.00084	RCP 30 6 30 32	2	tlim	2.1e-05
					RCP 30 7 30 18	1	13.89	8.4e-10
					RCP 30 8 30 40	1	21.55	2.5e-10
					RCP 30 9 30 41	2	24.42	6.5e-09
					RCP 30 10 30 46	1	tlim	0.000776
					RCP 30 11 30 34	2	tlim	0.000154
					RCP 30 12 30 36	1	tlim	0.000914
					RCP 30 13 30 30	1	tlim	7.7e-05
					RCP 30 14 30 39	2	tlim	6.8e-05
					RCP 30 15 30 30	1	tlim	0.000281

Table 2 :

 2 Results on RCP instances: NLP-Penalty problem.

	NLP-Penalty

Table 3 :

 3 Results for large RCP instances: three-phase method. Global exact solutions (except for RCP30 5)

					2nd phase: MINLP f ix	3rd phase: MINLP	Total
	Name	n	nc	n hth	time (s)	speed dev.	time (s) speed dev.	time (s)
	RCP 30 1	30	35	1	34.93	1.4e-06	7.900	4.22e-17	45.94
	RCP 30 2	30	38	1	59.69	6.5e-15	-	-	60.76
	RCP 30 3	30	46	1	2.920	1.7e-16	-	-	3.816
	RCP 30 4	30	39	1	62.98	2.6e-16	-	-	65.60
	RCP 30 5	30	36	2	43.98	2.2e-06	tlim	2.24e-06	tlim
	RCP 30 6	30	32	2	43.99	1.5e-16	-	-	46.20
	RCP 30 7	30	18	1	105.2	1.2e-16	-	-	106.8
	RCP 30 8	30	40	1	2.700	1.3e-17	-	-	4.124
	RCP 30 9	30	41	2	51.23	3.3e-06	20.38	1.89e-17	76.56
	RCP 30 10 30	46	1	67.59	7.3e-07	-	-	73.46
	RCP 30 11 30	34	2	51.79	4.4e-16	-	-	52.79
	RCP 30 12 30	36	1	86.90	4.8e-15	-	-	88.45
	RCP 30 13 30	30	1	4.092	3.6e-16	-	-	6.632
	RCP 30 14 30	39	2	3.752	8.8e-18	-	-	23.81
	RCP 30 15 30	30	1	79.32	9.9e-16	-	-	84.05

Otherwise, the fact that i and j are separated at t = 0 is equivalent to having i and j separated for all t ≥ 0; equivalently, setting v ij =0 in (2) boils down to the initial separation ∥x 0 ij ∥ ≥ d.

Acknowledgements S. Cafieri and M. Mongeau thank Barbara Conn for proofreading the manuscript.


domain: spi sde

text:
Maxime Jean
email: maxime.jean@ifpen.fr

Pascal Granier
email: pascal.granier@ifpen.fr

Thomas Leroy
email: thomas.leroy@ifpen.fr

Combustion Stability Control Based on Cylinder Pressure for High Efficiency Gasoline Engines

Keywords: hybrid gasoline engine, exhaust gas recirculation, combustion stability, in-cylinder pressure sensor

Introduction

The transport sector is today responsible for a large part of greenhouse gas emissions and poor air quality in urban areas. Thus, reducing greenhouse gas emissions and improving air quality constitute the main challenges for the transport sector today. To meet this ambition, automotive industry must improve vehicle fuel consumption and pollutant emissions via the development of new technologies.

Electrification is a key solution to this challenge, particularly in hybridized powertrains. The penetration rate of hybrid vehicles is expected to grow in the 2030 horizon and is encouraged to maintain a massive use of combustion engines [START_REF] Kapus | Passenger car powertrain 4. x-from vehicle level to a cost optimized powertrain system[END_REF]. According to these trends, gasoline engines will be the most common technology for the combustion engine. Yet efforts can still be made to improve its efficiency. Indeed, it is known that combustion engine efficiency is directly related to the overall hybrid vehicle consumption [START_REF] Gautrot | Swumble In-Cylinder Fluid Motion for High Efficiency Gasoline SI Engines: Development of the second generation[END_REF].

A widespread axis of research is to increase the compression ratio. To avoid limitations of the consequent knocking phenomenon, this high compression ratio combustion system can be associated with complex air loop systems including Miller cycle operation and Exhaust Gas Recirculation (EGR). On this kind of engine, a precise control of the engine air filling is essential to manage pollutant emissions (through air-fuel ratio control), and to manage the overall fuel consumption (through air/EGR/ignition control).

EGR control is a major challenge of new high efficient gasoline engines. A high level of EGR is required to minimize the consumption, but crossing a certain limit would compromise the stability of the combustion. The EGR feedback control strategy uses an estimation of the EGR quantity as an input. Nevertheless, EGR estimation is highly scattered due to part to part dispersion and sensor accuracy [START_REF] Kiwan | Effects of differential pressure measurement characteristics on low pressure-EGR estimation error in SI-engines[END_REF][START_REF] Nyerges | Verification and Comparison of Nine Exhaust Gas Recirculation Mass Flow Rate Estimation Methods[END_REF]. To prevent excessive EGR filling that would jeopardize combustion stability, a safety margin is added to the EGR target calibration map. Stability is thus ensured but it leads to a fuel consumption increase.

The challenge is to operate the engine near the stability limits to obtain the lowest consumption. A real-time estimation of the combustion stability would be precious information for improving the robustness of the engine control. For a few years, in-cylinder pressure sensors have been spreading among production engines [START_REF] Booven | Potential of Series-Compatible In-Cylinder Pressure Sensors for Gasoline Engines Using the Example of Ignition Angle Control[END_REF][START_REF] Booven | Real-time assessment of in-cylinder pressure for fresh air mass estimation using the example of modern spark-ignited combustion engines[END_REF]. This kind of sensor could be leveraged to obtain the stability information and thus allow feedback control.

Some reference articles deal with the problem of controlling engine combustion. In [START_REF] Maldonado | Artificial-intelligence-based prediction and control of combustion instabilities in spark-ignition engines[END_REF] and [START_REF] Singh | Neural network control of spark ignition engines with high EGR levels[END_REF], the authors model the combustion process thanks to Artificial Neural Networks to control the cycle-to-cycle injected fuel mass. In [START_REF] Norouzi | Model predictive control of internal combustion engines: A review and future directions[END_REF], the Model Predictive Control technique is used as a multi-variable controller to control IMEP and combustion phasing. Those papers deal with combustion control using IMEP or combustion phasing measurements. They do not use any direct estimation of the stability level as feedback information.

In [START_REF] Yao | Combustion Stability Control of Dieseline PPCI Based on In-Cylinder Pressure Signals[END_REF], the authors propose an online estimation of combustion stability of a Diesel PPCI (Partially Premixed Compression Ignition) engine. They use the indicator to control injection timing and EGR rate. In [START_REF] Haertl | Experimental Investigation of a Control Strategy Based on Combustion Stability and Combustion Phasing for a Multi-Cylinder Engine with Fueled Pre-Chambers and Cylinder Pressure Transducers[END_REF], a control strategy estimating combustion stability and phasing for a multicylinder engine with fueled prechambers is explored. However, for both references, the proposed estimators rely on a class of Coefficients of Variance (cov) that are highly biased during transient operating conditions. This makes it irrelevant to use these estimators in a real-time control strategy.

The contribution of this paper is a robust estimation of combustion instabilities, based on the in-cylinder pressure sensor signals, even during transient operating conditions. A feedback loop controller can maintain the engine as close as possible to the limit of stability, thus improving the fuel consumption. The strategy is implemented on an ECU and experimental results obtained on a high efficiency gasoline engine stress the relevance of the proposed approach.

This paper is organised as follows. Section 2 presents the compromise between fuel consumption minimization and combustion instabilities. Section 3 presents the proposed combustion instability estimator combined with its feedback loop controller. Section 4 gives some details of the experimental setup and control implementation. Finally, Section 5 shows the experimental results of the proposed controller on a high efficiency gasoline engine performing hybrid real-life typical operations.

Engine Stability and Fuel Consumption Compromise

Stationary Regime

It is well known [START_REF] Wei | Gasoline engine exhaust gas recirculation-A review[END_REF][START_REF] Tornatore | Experimental and numerical study on the influence of cooled EGR on knock tendency, performance and emissions of a downsized spark-ignition engine[END_REF] that the engine combustion efficiency increases with EGR rate. Up to a certain limit, excessive EGR rate generates combustion instabilities [START_REF] Randolph | The Effects of EGR Composition on Combustion Performance and Efficiency[END_REF]: Indicated Mean Effective Pressure (I MEP) fluctuations and misfires. The I MEP is derived cycle by cycle from in-cylinder pressure sensor signal processing. Combustion stability can also be derived from an ionization measurement signal. This paper uses a cylinder pressure measurement but interested readers can refer to [START_REF] Byttner | Estimation of Combustion Variability Using in-Cylinder Ionization Measurements[END_REF][START_REF] Kaiadi | Closed-Loop Combustion Control Using Ion-Current Signals in a 6-Cylinder Port-Injected Natural-Gas Engine[END_REF][START_REF] Byttner | Closed-loop control of EGR using ion currents[END_REF] for other techniques (the proposed approach can also be used with this alternative measurement). It is well established that the Coefficient of Variation of the I MEP defines the cyclic variability of the combustion. Standard estimation of combustion stability is the Coefficient of Variance (cov) of the I MEP, given by:

cov(I MEP) = σ(I MEP) µ(I MEP) , (1)
where σ(I MEP) and µ(I MEP) are, respectively, the standard deviation and the mean value of the random variable I MEP. The estimator cov of cov is generally computed with the conventional estimators of the mean value µ and of the standard deviation σ. They are estimated with the I MEP observations of a window of n consecutive cycles. A typical value for n is several hundreds, such as 300. It is generally estimated on stationary operating conditions on an engine test bench for mapping purposes.

Figure 1 shows the experimental evolution of the corresponding Indicated Specific Fuel Consumption (ISFC) (Gross IMEP is considered in the experimental results, the impact of pumping losses on the system is out of the scope of the article), and of cov(I MEP), both estimated on windows of 300 cycles. It shows a progressive decrease of the fuel consumption (until frequent misfires increase the consumption), while cov(I MEP) stays constant until it reaches the limit of stability and rapidly increases. The corresponding experimental setup is described in Section 4.

Standard Engine Calibration

The usual calibration of EGR maps dedicated to open loop controllers involves experimental acquisitions pictured in Figure 1. For a given operating point defined by the engine speed and load, the limit of stability is identified with a gradual increase of EGR rate. An added protection margin prevents overshooting the limit in the case of bad estimation of the EGR rate or unexpected changes due to environmental conditions, as shown in Figure 2. In function of the performance of the EGR estimation, this margin can rise to more than 5% of EGR, and thus limits the efficiency benefit brought by EGR technology.

Benefit from a Closed Loop Controller

The margin described in Section 2.2 prevents unacceptable instabilities at the expense of combustion efficiency. An online estimation of cov(I MEP) based on in-cylinder pressure sensor signals could describe the combustion instabilities. A feedback loop controlling the EGR rate could then maintain cov(I MEP) in the vicinity of its target value (usually set around 3%). It would both remove the conservative margin and avoid unstable operating conditions.

Estimation Bias in Transient States

The traditional estimation cov(I MEP) is performed on an engine test bench during several hundred consecutive cycles. However, it would not correspond to the desired dynamics for real-time control applications. Shortening the window size is a first answer to circumventing the problem. In the following sections, the window size n will be set to 30. However, decreasing the window size increases the variance of the estimator cov. Also, how short the window can be, the proposed estimator is naturally biased during transient states. Indeed, σ(I MEP) sums up the undesired variability of the I MEP around its set-point value with the variability in time of this set-point value. The engine running frequently in rapid transient states, this estimator can not be used in a feedback loop controller. An example of such bias can be seen in Figure 3. It shows a portion of WLTP cycle. The top plot shows the I MEP and the bottom plot (blue) shows the estimator cov(I MEP). It is particularly biased on the load transients at t = 23 s, 28 s, 45 s and 58 s.

Proposition

Combustion Stability Indicator

The estimator cov(I MEP) does not adequately reflect the combustion instabilities in transient regimes, even if it stays valid in the stationary regime of the protocol described in Section 2.1.

A straightforward approach would consist of a high-pass filtering of the IMEP signal. A high-pass approach would assume that the discrete time IMEP is the combination of low frequency changes of engine torque (desired), and of high frequency combustion instabilities (undesired), with no overlapping between these two sets. However, such a separation cannot be made in practice: some very fast variations cannot be considered noise since they correspond to desired engine load transient. The existence of an appropriate cutoff frequency being questionable, its calibration in practice is not an easy task. To avoid these drawbacks, we propose another approach.

The measurement of interest is not so much the normalized spread of I MEP; it is the undesired variations of I MEP around its moving setpoint target. The proposed solution to overcome the natural bias of cov(I MEP) lies in the variable X, measurable cycle by cycle, whose cov appropriately reflects combustion instabilities, while invariant under load changes.

X = I MEP Q f uel η , ( 2
)
where Q f uel is the injected mass of fuel and η is the engine efficiency given by a map.

The implementation proposed in this paper uses an efficiency map function of engine speed and IMEP. However, it could also be function of any other observed or estimated variable.

Another writing of X is: X = k η meas η , with η meas being the instantaneous efficiency computed online from the I MEP measurement and from the injected mass of fuel, and k being a constant given by the engine characteristics (number of cylinders, engine displacement and fuel Lower Heating Value). For the engine configuration and fuel properties described in Section 4.2, k = 1.03.

If the efficiency map had no bias and if there were no combustion instabilities, X would stay constant, during stationary and transient states. If there are some instabilities, cov(X) would then properly quantify them. Constant multiplicative biases of the map does impact X, but not cov(X) (since the cov is invariant under multiplication by a constant). During transient regimes, efficiency map biases could have an impact, but only if the bias significantly changes between operating points within the same estimation window. In this study, the efficiency map is only a function of the two main variables (engine speed and IMEP). Taking into account additional variables would slightly improve the prediction.

The variable X can be seen in the example in Figure 3 (middle plot), and its estimated cov(X) (bottom plot, red). It shows that X is far less subject to transient variations than I MEP, and consequently that cov(X) is less impacted by the transient regimes.

Combustion Stability Control

Given that the engine has several cylinders (three in the configuration given Section 4.2), but one air loop, the largest estimated cov(X) among the different cylinders is used by the controller, corresponding to the most unstable cylinder. This feedback can then be compared to the target value θ and used in an integral feedback loop to control an EGR offset ∆EGR from the default calibration EGR map . The gain K of the controller must be calibrated to avoid entering the unstable regime. The control law is then:

EGR sp = ∆EGR + EGR map , with ∆EGR = K ∑ t max i∈cyclinders [ cov(X)] i -θ . ( 3
)
This implementation applies the same offset on all engine operating points. This global correction is the simplest implementation of this control algorithm, tested in this paper. However, a more suitable algorithm would finely tune different EGR offset ∆EGR for the different operating points: it would consist in implementing an adaptive map. This work is out of scope of this paper, but could have a straightforward implementation.

If a misfire is encountered, the stability indicator computation is reset, leading to an adjustable decrease of the offset over the 30 following cycles. The instability estimation features some variability because of the short estimation window (n = 30). However, a small integral gain makes its impact on the controlled EGR setpoint insignificant. The stability estimator target must be carefully calibrated as well for the same reasons. We started from the usual value of 3% and decreased it to 2% according to our experimental results.

Experimental Setup 4.1. Engine Setup

The control strategy is validated on a cutting edge prototype engine. A three cylinders internal combustion engine was designed and manufactured [START_REF] Leroy | New Combustion Process for Higher Efficiency of Gasoline Engines[END_REF], based on a PSA EB2ADTS 96 kW engine. The cylinder head was replaced by an IFPEN cylinder head featuring Swumble TM technology. It was equipped with intake and exhaust Variable Valve Timing systems, featuring a Miller intake lift duration of 140°. Injection pressure was boosted to 350 bar to improve air-fuel mixing and thus reduce particulate emissions. The engine characteristics are given in Table 1. The air loop was also redesigned. It implements Exhaust Gas Recirculation. EGR is an important lever to reduce the engine knock limitations. These modifications induce strong constraints on the intake pressure and reduce exhaust enthalpy. Therefore, the choice of a variable geometry turbocharger has been made for its versatility and efficiency. Overall, the air loop has been changed for cutting edge performances, Figure 4.

Test-Bench Setup

In-cylinder pressure is measured for each cylinder with a Kistler 6041B sensor. The fuel consumption is measured with a Coriolis mass flow meter AVL KMA 4000.

A fully open engine control unit was developed using a McLaren TAG400i unit. Cylinder pressure sensors are sent to the ECU via CAN thanks to a rapid-prototyping platform [START_REF] Grondin | Rapid-prototyping multi-sensors processing platform for real time engine control and diagnosis[END_REF].

The fuel used for these tests is standard E10 RON 95 gasoline. Humidity of the air is regulated to 50% of relative humidity and 20 °C at sea level pressure. The measurements were taken in hot engine conditions, at 90 °C for oil and water temperatures. Water Charge Air Coolant coolant temperature is regulated at 25 °C.

Control Strategy Implementation

The engine control strategy detailed in Section 3 is implemented in the ECU. Figure 5 gives some details of the implementation. I MEP computation is performed in the rapid prototyping platform and sent via CAN to the ECU. The variable X given in (2) uses the IMEP, the fuel mass flow Q f uel coming from the fuel controller and the efficiency η computed through a look-up table function of the engine speed and I MEP. A buffer is then filled with 30 values of X and its cov is computed. The difference between the cov and its target is integrated and is added to the EGR target value EGR map provided by a calibration map. This EGR setpoint is then provided to the EGR control that manages the engine actuators.

Usage of any cylinder pressure sensor able to provide a high frequency signal allowing an IMEP computation could be suitable for the control strategy, interested readers can refer to [START_REF] Vollberg | Cylinder pressure sensors for smart combustion control[END_REF][START_REF] Hellemans | In-cylinder pressure sensor[END_REF], for example.

Experimental Results

Torque Trajectory on WTLP Cycle

In this Section, the controller is tested in typical driving conditions. The I MEP trajectory is obtained thanks to hybrid vehicle simulation [START_REF] Leroy | Swumble 3-Cylinder High Efficiency Gasoline Engine for Future Electrified Powertrains[END_REF] implementing energy management strategy to split the IMEP request between the engine and the electric motor. Figure 6 presents engine speed and IMEP as a function of time for a part of the WTLP cycle. The right hand side figure presents the engine efficiency as a function of engine speed and I MEP. The bubbles highlight the most frequently used operating points. As one can notice, the engine is used in its high efficiency area, taking full advantage of the hybrid architecture. In this area, highest EGR rates are required to optimize the fuel consumption. This is then a very interesting situation to validate the presented strategy.

The engine speed trajectory tracking is ensured by the generator load controller of the test bench. In parallel, the IMEP trajectory setpoint is sent to the ECU.

Strategy Validation

The estimation window of 30 cycles offers a compromise between precision and estimation speed, required to operate in transient regimes. The integral gain of the controller has to be small enough to filter out the estimation variability (because the size of the window is small, the value of cov(X) highly depends on the realization of the 30 random variables, the estimator has then a large variance) and to encompass several operating point changes within its filter characteristic time. In fact, the purpose of the proposed controller is not to quickly react to engine transients, but to slowly adapt the EGR prepositioning map EGR map . Given all these aspects, the gain K was set at a value of 0.001 in the experiments. The gain is the same for positive and negative actions but could be differentiated if the application requires it.

Since a small integral gain is selected, the controller needs to be tested on a long duration. A specific pattern of 86 s of the WLTP cycle presented in Section 5.1 has been repeated 21 times. The strategy is deactivated during engine stops. The feedforward EGR map map has been deliberately lowered (towards the safety direction) to highlight the behaviour of the controller. The response of the controller during one pattern is shown Figure 7. Every curve is given as a function of the time (one second corresponds to around 70 cycles at 2700 rpm). The controller increases the EGR correction ∆EGR of 0.6% during this pattern. Despite the IMEP dynamics (top plot), cov(X) (second plot) still provides a relevant estimation of the engine stability. It is not much affected by load transients in themselves, but correctly reflects the I MEP variability around its set trajectory such as times 40 s and 72 s. cov(I MEP) is also plotted in the second figure to stress that cov(X) is much more robust to IMEP transient variations. Figure 8 shows the behaviour of the controller during the 21 patterns. At the beginning of the test, EGR setpoint is set to its feedforward mapped value (as discussed previously, this value is derated to take into account estimation uncertainties and variation of environmental conditions). The objective of this experiment is to show the benefit of the feedback loop over consumption. From the first to the 16th cycle, the average cov(X) along each pattern has increased of 0.1%. After 16 repeated patterns, the controller achieves the stabilisation of the offset ∆EGR at +3.5% (third plot). When the controller converges, there is a 1.7% gain in IFSC in regard to the initial conditions. The correction maintains its average value during the next 5 patterns. The experiment shows that the controller can shift a conservative calibration map towards the limit of stability until cov(X) reaches the target (2%). No misfires were observed during the cycle.

Conclusions

This paper presents a first experimental implementation of a feedback loop control of the EGR rate ensuring higher efficiency while avoiding unstable combustion. An online estimation of the Coefficient of Variance of the engine I MEP for control purposes has been proposed. It has been used to add a correction to the EGR propositioning map. It allows us to maintain the engine near the limit of stability and alleviates the efficiency losses due to the conservative prepositioning EGR map. The simple implementation proposed in the paper applies a global correction on the EGR map, given by an integral controller and based on the error between the maximum (among the cylinders) estimated Coefficient of Variance and a target value. Further studies could focus on adaptive map implementations.

Experimental results showed, on a repetition of a portion of hybrid WLTP cycle, that the controller moves the engine towards the limit of stability, and that it converges before reaching the unstable domain. The test cycle is relevant since it explores limited areas of high efficiency of the engine operating map. This area is small because it corresponds to a hybridized architecture which limits the transient rates. This algorithm could be easily implemented with no extra costs in engines already equipped with in-cylinder pressure sensors and could improve the engine efficiency by a few percent.

Figure 1 .

 1 Figure 1. Evolution of the Indicated Specific Fuel Consumption (ISFC) and of the estimated engine stability cov(I MEP), function of EGR rate.

Figure 2 .

 2 Figure 2. EGR calibration. 1 Calibration set point. 2 Optimal value.

Figure 3 .

 3 Figure 3. Portion of WLTP cycle. Top plot: I MEP, middle plot: X defined Equation (2), bottom plot: cov estimators. A clear bias of the first stability estimator can be seen during transient states.

Figure 4 .

 4 Figure 4. Prototype engine.

Figure 5 .

 5 Figure 5. Engine control architecture.

ersion 11 Figure 5 . 5 . 1 .Figure 6 Figure 6 .

 1155166 Figure 5. Engine control architecture.5. Experimental results206

5. 2 .

 2 Strategy validation 219 The estimation window of 30 cycles offers a compromise between precision and 220 estimation speed, required to operate in transient regimes. The integral gain of the con-221 troller has to be small enough to filter out the estimation variability 2 and to encompass 222 several operating point changes within its filter characteristic time. In fact, the purpose 223 of the proposed controller is not to quickly react to engine transients, but to slowly adapt 224 the EGR prepostioning map EGR map . Given all these aspects, the gain K was set at a 225 value of 0.001 in the experiments. The gain is the same for positive and negative actions

Figure 6 .

 6 Figure 6. Hybrid vehicle typical transient trajectory for strategy validation.

Figure 7 .

 7 Figure 7. Control of the EGR offset during one pattern of WLTP cycle. Top plot: I MEP of the third cylinder I MEP cyl3 and I MEP set point I MEP sp . Middle plot: maximum estimated cov over the 3 cylinders and its target θ. Bottom plot: EGR correction ∆EGR.

Figure 8 .

 8 Figure 8. Control of the EGR offset subject to repeated pattern of a portion of WLTP cycle.

Table 1 .

 1 Multi-cylinder engine main features.

	Engine displacement [L] Vol. compression ratio [-] Bore × Stroke [mm] Number of intake/exhaust valves Valve lift CA duration, at 1 mm [°C A] Injection system Injector Injection pressure [bar]	1.2 13.65:1 75 × 90.5 2/2 Intake 140, Exhaust 210 Central direct injection Bosch HDEV6, six holes 350

Funding: This research received no external funding.

Acknowledgments:

The authors would like to acknowledge Adrien Nail, Gérard Martin and Ludovic Nowak for the fine work performed in operating the multicylinder engine.

Author Contributions: Conceptualization, M.J. and T.L.; methodology, M.J.; software, P.G.; validation, P.G., M.J. and T.L.; formal analysis, M.J. and P.G.; investigation, M.J. and P.G.; writing-review and editing, M.J., P.G. and T.L.; visualization, M.J. and P.G.; supervision, T.L.; project administration, T.L. All authors have read and agreed to the published version of the manuscript.

Conflicts of Interest:

The authors declare no conflict of interest.