Amanda Matson
email: amatson2@kent.edu

Claude-Michel Brauner
email: claude-michel.brauner@u-bordeaux.fr

Peter V Gordon
email: gordon@math.kent.edu

Uniqueness of traveling fronts in premixed flames with stepwise ignition-temperature kinetics and fractional reaction order

Keywords: Reaction-diffusion systems, traveling front solution, uniqueness of solution, qualitative dependency of solution on parameters AMS subject classifications: 35K57, 35C07, 34B08, 34E05, 80A25

Introduction

The canonical constant density approximation model of flame propagation in one dimensional formulation reads [START_REF] Berestycki | Traveling wave solutions to combustion models and their singular limits[END_REF][START_REF] Law | Combustion Physics[END_REF][START_REF] Williams | Combustion theory[END_REF][START_REF] Ya | The mathematical theory of combustion and explosions[END_REF]:

Tt = T xx + Ω(T, C), Ct = ΛC xx -Ω(T, C), (1.1)
where T and C are appropriately normalized temperature and concentration of the deficient reactant, x ∈ R, t > 0 are normalized spatiotemporal coordinates, Λ > 0 is an inverse of the Lewis number and Ω is the reaction rate. The reaction rate is typically specified as: Ω(T, C) := C α F (T ) T ≥ θ and C > 0, 0 T < θ and/or C = 0, (

where α ≥ 0 is the reaction order, θ ∈ (0, 1) is the ignition temperature and F (T ) is a positive non-decreasing function that characterizes the enhancement of chemical reaction with the increase of the temperature. The studies of system (1.1) trace back to pioneering works of Frank-Kamenetskii, Semenov and Zeldovich in 1930's and 1940's [START_REF] Frank-Kamenetskii | Diffusion and heat transfer in chemical kinetics[END_REF][START_REF] Semenov | Thermal theory of combustion and explosion[END_REF][START_REF] Ya | The mathematical theory of combustion and explosions[END_REF]. This system was then analyzed by mathematicians and physicists alike. The substantial body of literature concerning system (1.1) is dedicated to the analysis of traveling front solutions, that is special solutions of the form: T (x, t) := u(ξ), C(x, t) := v(ξ), ξ = x + c t, (1.3) where c > 0 is an a-priori unknown speed of propagation. In a context of combustion, such solutions represent flame fronts propagating with a constant speed from burned state far to the right to unburned state far to the left. When α ≥ 1 system (1.1), after substitution of an ansatz (1.3), reduces to the following system of ODE's on a real line:

u ξξ -cu ξ + Ω(u, v) = 0, Λv ξξ -cv ξ -Ω(u, v) = 0, ξ ∈ R, (1.4)
complemented with boundary like conditions:

u → 1, v → 0 as ξ → ∞, and u → 0, v → 1 as ξ → -∞.

(1.5)

Conditions (1.5) prescribe the steady temperature and reactant concentration far ahead (ξ → -∞) and far behind (ξ → ∞) the flame front. Since solutions of (1.4), (1.5) are translationally invariant, we fix translations by imposing a constraint:

u(ξ ign ) = θ, (1.6)
where ξ ign is an arbitrary fixed number. The constraint (1.6) fixes the position of an ignition interface, the unique position where the temperature is equal to the ignition temperature θ. Hence, when crossing the ignition interface, the reaction rate jumps from some positive value to zero while preserving continuity of the temperature and concentration of reactant as well as their fluxes. Consequently, the mixture ahead of the ignition interface (ξ < ξ ign ) is in a non-reactive state, whereas at and behind the ignition interface (ξ ≥ ξ ign ), the chemical reaction takes place. We note that uniqueness of the ignition interface follows from the monotonicity of any solution of problem (1.4), (1.5) that results from the fact that Ω(u, v) ≥ 0 and can be directly verified.

In view of the discussion above, system (1.4) is equivalent to the following one:

u ξξ -cu ξ = 0, Λv ξξ -cv ξ = 0, ξ ∈ (-∞, ξ ign ), u ξξ -cu ξ + v α F (u) = 0, Λv ξξ -cv ξ -v α F (u) = 0, ξ ∈ (ξ ign , ∞), (1.7)
complemented with conditions of continuity of a solution and its first derivatives when crossing the interface:

[u] = [v] = [u ξ ] = [v ξ ] = 0 at ξ = ξ ign . (1.8)
Here and below [•] stands for a jump of the quantity. A sketch of a traveling front solution for system (1.1) with the reaction order α ≥ 1 is depicted on Figure 1. System (1.4), (1.5), (1.6) (equivalently system (1.5), (1.6), (1.7), (1.8)) was extensively studied in the literature. In the special case when Λ = 1, this system reduces to a single equation. This equation is well understood by now, details can be found in many books and review articles on the subject, for example [START_REF] Gilding | Travelling waves in nonlinear diffusion-convection reaction[END_REF][START_REF] Volpert | Traveling wave solutions of parabolic systems[END_REF][START_REF] Xin | Front propagation in heterogeneous media[END_REF][START_REF] Ya | The mathematical theory of combustion and explosions[END_REF]. In particular, it is well known that in this case (1.4), (1.5), (1.6) admits a unique traveling front solution for any Λ ∈ (0, ∞) fixed. The general case when Λ = 1 is substantially more complex, and the complete understanding of this case is still lacking. Review of many relevant results concerning this general case can be found in [START_REF] Volpert | Traveling wave solutions of parabolic systems[END_REF][START_REF] Volpert | Elliptic partial differential equations[END_REF]. In particular, it was shown in [START_REF] Berestycki | Traveling wave solutions to combustion models and their singular limits[END_REF] that when the reaction rate α is a positive integer then system (1.4), (1.5), (1.6) admits a traveling front solution. Moreover, when α = 1 and Λ ∈ (0, 1] this solution is unique [START_REF] Kanel | A stationary solution to a system of equations in the theory of burning[END_REF]. The question whether uniqueness of a solution holds for α = 1 and Λ > 1 is still open.

When 0 < α < 1, non-linear term (1.2) becomes non-Lipschitz at C = 0. Consequently, the system that describes traveling fronts for (1.1) changes substantially. After substitution of an ansatz (1.3) into (1.1) this system reads:

   u ξξ -cu ξ = 0, Λv ξξ -cv ξ = 0, ξ ∈ (-∞, ξ ign ), u ξξ -cu ξ + v α F (u) = 0, Λv ξξ -cv ξ -v α F (u) = 0, ξ ∈ (ξ ign , ξ tr ), u ξξ -cu ξ = 0, v = 0, ξ ∈ (ξ tr , ∞), (1.9)
where ξ tr > ξ ign is an a-priori unknown constant which we refer to as a position of trailing interface. The trailing interface, in the context of combustion, indicates the leftmost point where entire reactant available for the reaction is fully consumed. As a result of it, as can be easily checked, the temperature at the trailing interface as well as at all points to the right of it is equal to the temperature of the fully reacted mixture (u = 1, ξ ≥ ξ tr ). Hence, the distance between the ignition and trailing interfaces R = ξ tr -ξ ign is the length of the reaction zone.

In view of the discussion above, the last equation in (1.9) can be solved explicitly and the solution reads:

u = 1, v = 0, ξ ∈ (ξ tr , ∞). (1.10)
Hence, system (1.9) reduces to

u ξξ -cu ξ = 0, Λv ξξ -cv ξ = 0, ξ ∈ (-∞, ξ ign ), u ξξ -cu ξ + v α F (u) = 0, Λv ξξ -cv ξ -v α F (u) = 0, ξ ∈ (ξ ign , ξ tr ), (1.11)
complemented with boundary conditions at the trailing interface:

u = 1, u ξ = 0, v = v ξ = 0 at ξ = ξ tr , (1.12)
that manifest continuity of the temperate and reactant and their fluxes when crossing the trailing interface. Problem (1.11), (1.12) should be considered with constraint (1.6), continuity conditions at the ignition interface (1.8) and boundary like conditions far ahead of the ignition interface:

u → 0, v → 1 as ξ → -∞. (1.13)
We note that the main principle difference of a system describing traveling front solution for system (1.1) with α ≥ 1 and 0 < α < 1 is that the latter involves a trailing interface which position is a-priori unknown. We also note that any solution of (1.6), (1.8), (1.11),(1.12), (1.13) is monotone so for any such solution positions of both ignition and trailing interfaces are uniquely defined. The sketch of the traveling front solution for this system is depicted in Figure 2. Despite the fact that reactions of fractional order (0 < α < 1) are common in various combustion applications, for example in high pressure combustion systems (see e.g [START_REF] Law | Combustion Physics[END_REF]), problem (1.6), (1.8), (1.11),(1.12), (1.13) has not received appropriate attention in mathematical literature partially due to the fact that the non-Lipschitz reaction rate creates multiple technical difficulties for the analysis. To the best of our knowledge, a recent paper [START_REF] Brauner | Existence of a traveling wave solution in a free interface problem with fractional order kinetics[END_REF] is the only mathematical paper where this system was analyzed. In [START_REF] Brauner | Existence of a traveling wave solution in a free interface problem with fractional order kinetics[END_REF] the authors considered (1.6), (1.8), (1.11), (1.12), (1.13) with a stepwise ignition-temperature kinetics, that is with F (T ) being a Heaviside step-function The position of the trailing interface, the point to the right of which temperature and concentration of deficient are identically one and zero (u = 1, v = 0), respectively is indicated by ξ tr .

F (T ) := 1 T ≥ θ, 0 T < θ. (1.14)
The particular choice of stepwise ignition-temperature kinetics, from perspective of applications, is justified by the fact that such kinetic is apparently the most appropriate approximation of overall reaction kinetics for certain hydrogen-oxygen/air and ethylene-oxygen mixtures as evident from multiple theoretical and numerical studies of detailed reaction mechanism in such mixtures see [START_REF] Sánchez | Recent advances in understanding of flammability characteristics of hydrogen[END_REF] and references therein. Under the assumption of stepwise ignition-temperature kinetics, system (1.11) simplifies to :

u ξξ -cu ξ = 0, Λv ξξ -cv ξ = 0, ξ ∈ (-∞, ξ ign ), u ξξ -cu ξ + v α = 0, Λv ξξ -cv ξ -v α = 0, ξ ∈ (ξ ign , ξ tr ).
(1.15)

In [START_REF] Brauner | Existence of a traveling wave solution in a free interface problem with fractional order kinetics[END_REF] it was shown that problem (1.6), (1.8), (1.15), (1.12), (1.13) admits a solution, however, the question of multiplicity of solution was not addressed there. In this paper, we prove that a solution constructed in [START_REF] Brauner | Existence of a traveling wave solution in a free interface problem with fractional order kinetics[END_REF] is unique and discuss some qualitative properties of the solution.

The paper is organized as follows. In section 2, we set up the stage by giving some preliminaries and stating the main result. In section 3, we give a proof of necessary lemmas and propositions used in the proof of the main result. In the last section, we discuss certain properties of traveling front solutions and present results of numerical simulations. In particular, we give some formal arguments based on asymptotic and numerical studies of the problem that the velocity of propagation of the flame front is a decreasing function of the ignition temperature θ, inverse of the Lewis number Λ and the reaction order α.

2 Preliminaries and the statement of the main result.

In this section, we give an alternative formulation of problem (1.6), (1.8), (1.15), (1.12), (1.13), and state the main result of this paper. Let us first note that in a view of translational invariants of traveling fronts, we can fix the position of trailing interface ξ tr and treat ξ ign as a-priori unknown. Hence, we set ξ tr = 0 and ξ ign = -R where R > 0 is a-priori unknown. With this alteration, the problem describing traveling front solutions for (1.1) with α ∈ (0, 1) and step-wise reaction kinetics discussed in the introduction reads:

u ξξ -cu ξ = 0, Λv ξξ -cv ξ = 0, v > 0, ξ ∈ (-∞, -R), u ξξ -cu ξ = -v α , Λv ξξ -cv ξ = v α , v > 0, ξ ∈ (-R, 0), (2.1)
This system is complemented with boundary like conditions far ahead of the combustion front:

u → 0 v → 1 as ξ → -∞, (2.2)
continuity conditions for a solution and its first derivatives at the ignition interface:

[u] = [v] = [u ξ ] = [v ξ ] = 0 at ξ = -R, (2.3)
boundary condition at the ignition interface:

u = θ, at ξ = -R, (2.4)
and boundary conditions at the trailing interface

u = 1, u ξ = 0, v = v ξ = 0 at ξ = 0 (2.5)
Straightforward computations show that for any solution of problem (2.1) satisfying (2.3) and (2.4) has:

u(ξ) = θ exp(c(ξ + R)), ξ ≤ -R, (2.6)
and thus

u = θ, u ξ = cθ at ξ = -R. (2.7)
Combining elementary computations above, we conclude that any solution of (2.1)-(2.5) on an interval [0, R] must verify the following boundary value problem:

   u ξξ -cu ξ = -v α , Λv ξξ -cv ξ = v α v > 0, ξ ∈ (-R, 0), u = θ, u ξ = cθ ξ = -R, u = 1 u ξ = 0, v = v ξ = 0, ξ = 0. (2.8)
Therefore, the necessary condition for problem (2.1)-(2.5) to have a solution is the existence of solution for problem (2.8). We note that the existence of solution of (2.8) does not guarantee existence of solution for (2.1)-(2.5). Indeed, one can easily verify that for any solution of (2.1)-(2.5) we must have

v(ξ) = 1 -a exp c Λ (ξ + R) ξ ∈ (-∞, -R), (2.9)
for some constant a ∈ (0, 1). This constant, as follows from continuity conditions (2.3), has to be chosen from the overdetermined system:

1 -a = lim ξ→-R + v(ξ), -ca = Λ lim ξ→-R + v ξ (ξ), (2.10)
which is not necessarily consistent. However, [3, Theorem 1.1] guarantees that for any set of relevant parameters ( θ ∈ (0, 1), α ∈ (0, 1), Λ ∈ (0, ∞)) system (2.1)-(2.5) admits a solution. In view of this result, problem (2.8) admits a solution and uniqueness of the solution for problem (2.8) implies uniqueness of the solution for (2.1)-(2.5). The main result of this section is as follows:

Theorem 2.1. Let θ ∈ (0, 1), α ∈ (0, 1), Λ ∈ (0, ∞).
Then, problem (2.8) admits a unique solution.

Let us now show that the proof of Theorem 2.1 reduces to the analysis of the solution of a single second order ODE. First, let us observe that integration of the equation for temperature u in (2.8) and integration of the same equation multiplied by exp(-cξ) taking into account boundary conditions give:

c = ˆ0 -R v α (ξ)dξ, (2.11) cθ = ˆ0 -R v α (ξ) exp(-c(ξ + R))dξ.
(2.12)

Next we introduce the following rescaling:

v(ξ) = A w(x), (2.13)
where

x := cξ, A := c -2 1-α . (2.14)
Substituting (2.13), (2.14) into the equation for concentration of reactant v in (2.8) and taking into account boundary conditions, we have:

Λ w -w -wα = 0, w > 0, x ∈ (-cR, 0), w = 0, w = 0, x = 0, (2.15)
where the prime denotes the derivative with respect to the variable x.

Consider now the first equation of (2.15) extended on the half line, that is, the following problem

Λw -w -w α = 0, x < 0, w(0) = 0, w (0) = 0. (2.16)
Properties of (2.16) were studied in great detail in [3, Section 4] via a topological approach, see also appendix below. In particular, the following result comes directly from Lemmas 4.7 and 4.8 of [START_REF] Brauner | Existence of a traveling wave solution in a free interface problem with fractional order kinetics[END_REF], see appendix for more details.

Lemma 2.1. There exists a unique function w positive on (-∞, 0) that verifies (2.16).

It is clear that w is smooth as long as it does not vanish. The function w and its derivative can be extended by 0 up to 0 -. More precisely, it follows from [3, Lemmas 6.5 & 6.6] that the Hölder regularity of the function

w near x = 0 reads: w ∈ C ∞ (-∞, 0)) ∩ C 2+[β],β-[β] (-∞, 0]), β = 2α 1-α , 0 < α < 1. Moreover, w is decreasing on (-∞, 0].
As a consequence, for an arbitrary c > 0, R > 0 there is a unique decreasing positive solution to (2.15) and w(x) = w(x) on (-cR, 0].

In view of observations above we note that in terms of rescaled variables (2.11), (2.12) read:

c 2 1-α = ˆ0 -cR w α (s)ds, (2.17) θc 2 1-α = ˆ0 -cR w α (s) exp(-(s + cR))ds, (2.18)
which in particular imply

θ = ´0 -cR w α (s) exp(-(s + cR))ds ´0 -cR w α (s)ds . (2.19)
Introduce the following function:

φ(x) := ´0 -x w α (s) exp(-(s + x))ds ´0 -x w α (s)ds , x ∈ (0, ∞). (2.20)
The key observation that allows us to prove the main Theorem 2.1 is the following.

Proposition 2.1. Let w be the solution of (2.16) and set

φ(x) := ´0 -x w α (s) exp(-(s + x))ds ´0 -x w α (s)ds , x ∈ (0, ∞). ( 2

.21)

The function φ defined by (2.21) is strictly decreasing on (0, ∞) and

lim x→0 φ(x) = 1, lim x→∞ φ(x) = 0. (2.22)
Proof of Proposition 2.1 is given in the next section. This proposition makes the proof of Theorem 2.1 given below elementary.

Proof of Theorem 2.1. From (2.17) and (2.18) we have

c = ˆ0 -cR w α (s)ds 1-α 2 , (2.23) θ = ´0 -cR w α (s) exp(-(s + cR))ds ´0 -cR w α (s)ds . (2.24)
By strict monotonicity of φ (Proposition 2.1), there is a unique value σ * ∈ (0, ∞) such that φ(σ * ) = θ. Thus, in (2.23) we must have cR = σ * . This observation and uniqueness of the solution of (2.16) imply uniqueness of the speed c = c * with

c * = ˆ0 -σ * w α (s)ds 1-α 2
.

(2.25)

The latter immediately implies uniqueness of position of the ignition interface

R = R * = σ * /c * .
3 Proof of Proposition 2.1

In this section we present a proof of Proposition 2.1 which is based on the following three lemmas.

Lemma 3.1. Function φ(x) defined by (2.21) has the following properties:

lim x→0 φ(x) = 1, lim x→∞ φ(x) = 0. ( 3.1)
Proof. Let us first prove the first equality in (3.1). Observe that

e -x ˆ0 -x w α (s)ds < ˆ0 -x w α (s)e -(s+x) ds < ˆ0 -x w α (s)ds, x > 0. (3.2)
Hence,

e -x < φ(x) < 1, x > 0. (3.3)
Sending x to zero in the inequality above we have the first equality in (3.1). Now let us prove the second equality in (3.1). Integrating the equation in (2.16) and taking into account boundary conditions we have:

-Λw (-x) + w(-x) = ˆ0 -x w α (s)ds. (3.4)
Since w (-x) < 0 for x > 0 we have:

w(-x) < ˆ0 -x w α (s)ds. (3.5) Let ρ(x) = ´0 -x w α (s)ds w α (-x) . (3.6)
In view of the inequality above we have:

ρ(x) > w 1-α (-x), (3.7)
and thus

ρ(x) → ∞ as x → ∞, (3.8) as w(-x) → ∞ as x → ∞. Next observe that ˆ0 -x w α (s)e -(s+x) ds = ˆ0 -x+ √ ρ(x) w α (s)e -(s+x) ds + ˆ-x+ √ ρ(x) -x w α (s)e -(s+x) ds ≤ e - √ ρ(x) ˆ0 -x w α (s)ds + w α (-x) ρ(x) (3.9)
Dividing the expression above by ´0 -x w α (s)ds and using the definitions of φ(x) and ρ(x) we have: Assume there is

φ(x) ≤ e - √ ρ(x) + 1 ρ(x) . ( 3
x 0 ∈ R such that Y (x 0 ) > 0, Y (x 0 ) < 0, Y (x 0 )Y (x 0 ) > (Y (x 0 )) 2 .
(3.12)

Then, inf x∈[x0,∞) Y (x) > 0. (3.13)
Proof. Let

Y (x 0 ) = Y 0 , Y (x 0 ) = Y 0 , (3.14)
and assume that (3.12) hold. Then, as follows from (3.11) and the last inequality in (3.12)

Λ (Y 0 ) 2 Y 0 -Y 0 -Y α 0 < ΛY (x 0 ) -Y 0 -Y α 0 = 0, (3.15)
and hence

Λ Y 0 Y 0 2 - Y 0 Y 0 -Y α-1 0 < 0. (3.16) Therefore, µ -< Y 0 Y 0 < 0 < µ + , (3.17)
where

µ ± = 1 ± 1 + 4ΛY α-1 0 2Λ , (3.18)
are roots of the quadratic equation

Λµ 2 -µ -Y α-1 0 = 0. (3.19)
Next observe that a linear equation

ΛZ -Z -Y α-1 0 Z = 0, (3.20) satisfying Z(x 0 ) = Y 0 , Z (x 0 ) = -γY 0 , (3.21)
has the following solution

Z(x) =   2ΛY 0 1 + 4ΛY α-1 0   exp x -x 0 2Λ cosh   1 + 4ΛY α-1 0 2Λ (x -x 0 )   ×    |µ -| + 1 2Λ - 1 2Λ + γ tanh   1 + 4ΛY α-1 0 2Λ (x -x 0 )      . (3.22)
In view that for x ≥ 0 cosh(x) ≥ 1, 0 ≤ tanh(x) < 1 we have from (3.22) that

Z(x) > Z 0 = 2ΛY 0 1 + 4ΛY α-1 0 (|µ -| -γ) , x ∈ [x 0 , ∞). (3.23) Hence, Z(x) is uniformly bounded away from zero on [x 0 , ∞) provided γ < |µ -|. (3.24) Let x M be the largest x > x 0 such that 0 < Y (x) < Y 0 x ∈ I = (x 0 , x M ). ( 3
γ ∈ |Y 0 | Y 0 , |µ -| , (3.27)
fixed. Indeed, in view of the facts that Y (x 0 ) = Z(x 0 ) and Z (x 0 ) < Y (x 0 ) < 0 (which is guaranteed by our choice of parameter γ) we first observe that (3.26) holds for x ∈ (x 0 , x 0 + δ) with δ > 0 sufficiently small. Assume now x * is a smallest value of x ∈ I at which (3.26) is violated. Consider two possibilities of how it can happen. In the first one we have

Y (x * ) = Z(x * ) and Y (x * ) > Z (x * ), (3.28)
but then

Y (x * -ε) < Z(x * -ε), (3.29)
for ε > 0 sufficiently small. This contradicts the definition of x * and hence this situation is impossible. The second possibility is

Y (x * ) ≥ Z(x * ) and Y (x * ) = Z (x * ), (3.30)
In this case we have

ΛY (x * ) = Y (x * ) + Y α (x * ) > Y (x * ) + Y α-1 0 Y (x * ) ≥ Z (x * ) + Y α-1 0 Z(x * ) = ΛZ (x * ), (3.31)
and thus

Y (x * ) > Z (x * ). (3.32) But then Y (x * -ε) < Z (x * -ε). (3.33)
for ε > 0 sufficiently small, which again contradicts the definition of x * . Hence for x ∈ I we have Y 0 > Y (x) ≥ Z 0 > 0 and thus by continuity

Z 0 ≤ Y (x) ≤ Y 0 x ∈ [x 0 , x M ]. (3.34)
Clearly at x = x M we necessarily have

Y (x M ) = Y 0 , Y (x M ) ≥ 0. ( 3
Y (x) = Y (x M )e x-x M Λ + Λ -1 ˆx x M Y α (s)e
x-s Λ ds.

(3.36)

We now claim that Y (x) > Y 0 for x > x M . Indeed, by continuity Y (x) > Y 0 /2 for x ∈ [x M , x M + δ] for some possibly small δ > 0. Then by (3.36) we have

Y (x) ≥ Y 0 2 α e x-x M Λ -1 x ∈ (x M , x M + δ]. (3.37)
This implies that

Y (x) > 0 on (x M , x M + δ], (3.38
Y (x * M ) ≥ Y (x M + δ)e x * M -(x M +δ) Λ + Y α 0 e x * M -(x M +δ) Λ -1 > 0, (3.41)
and hence Y (x) is strictly increasing on (x M , ∞). Thus

Y (x) > Y 0 x ∈ (x M , ∞). ( 3
Y (x) ≥ Z 0 > 0 x ∈ [x 0 , ∞), (3.43)
that imply (3.13) and thus completes the proof.

Corollary 3.1. It follows immediately from Lemma 3.2 that the solution of (2.16) must satisfy

w (x)w(x) ≤ (w (x)) 2 , ∀x ∈ (-∞, 0). (3.44)
Indeed, if there was a point x 0 < 0 where (3.44) didn't hold, then it would imply that w(x) > 0 for x ∈ [x 0 , 0]. Thus such w would not satisfy boundary condition w(0) = 0 and hence would not solve (2.16) which gives a contradiction.

We note that formula (3.44) has a clear geometric interpretation. Consider problem (2.16) as a dynamical system on the plane. Introducing change of variables t = x, q(t) = w(x), p(t) = w (x), (3.45) we obtain the following system of two first order ODE's: q(t) = p(t), Λ ṗ(t) = p(t) + q α (t). (3.46) where dot stands for the derivative with respect to t. This problem is considered for t < 0 in the quadrant Q = {q 0, p 0}. Rewriting p and q in polar coordinates we have q(t) := r(t) cos(ϑ(t)), p(t) := r(t) sin(ϑ(t)), (3.47) where, t ∈ (-∞, 0), r > 0, θ ∈ (-π/2, 0). In terms of the new variables condition (3.44) is equivalent to

θ(t) ≤ 0. (3.48)
We present derivation of this formula in the appendix.

Lemma 3.3. Let w be the solution of (2.16) and set

f (x) := w α (x), x ≤ 0, 0, x > 0. (3.49)
For any fixed y > 0, the function

ψ(x, y) = ´∞ -x+y f (s)ds ´∞ -x f (s)ds , (3.50)
is non-decreasing function of x on (0, ∞) and strictly increasing on (y, ∞).

Proof. First we observe that ψ(x, y) ≡ 0 for x ≤ y.

(3.51) Hence, ψ(x, y) is non-decreasing in x for x ∈ (0, y]. For x > y we have,

∂ ∂x ψ(x, y) = f (y -x) ´∞ -x f (s)ds -f (-x) ´∞ y-x f (s)ds ´∞ -x f (s)ds 2 . (3.52)
Thus to show that ψ(x, y) is strictly increasing in x for x ∈ (y, ∞), we need to prove the following inequality ´∞ y-x f (s)ds This in particular implies that for x > y > 0 we have g(s + y -x) ≤ g(s -x) s ∈ (0, x -y).

f (y -x) < ´∞ -x f (s)ds f (-x) x ∈ (y, ∞). (3.53) Define g(x) := f (x) f (x) = α w (x) w(x) x ∈ (-∞, 0). (3.54) Observe that g (x) = α w 2 (x) w (x)w(x) -(w (x)) 2 , ( 3
(3.57)

Consequently, ˆs 0 g(τ + y -x)dτ ≤ ˆs 0 g(τ -x)dτ.

(3.58)

By the definition of the function g (see Eq. (3.54)) we have ˆs 0 g(τ + y -x)dτ = ˆs 0

d dτ f (τ + y -x) f (τ + y -x) dτ = ˆs 0 d dτ log f (τ + y -x) dτ = log f (s + y -x) f (y -x) , ˆs 0 g(τ -x)dτ = ˆs 0 d dτ f (τ -x) f (τ -x) dτ = ˆs 0 d dτ log f (τ -x) dτ = log f (s -x) f (-x) . (3.59)
Therefore, by (3.58) and (3.59) we have that for x > y > 0 and s ∈ (0, x -y)

log f (s + y -x) f (y -x) ≤ log f (s -x) f (-x) , (3.60)
that is

f (s + y -x) f (y -x) ≤ f (s -x) f (-x) . (3.61)
Using the facts that f (x) ≡ 0 for x > 0 and f (x) > 0 for s ∈ (-y, 0) for y > 0 we observe that ˆx-y In computation above we used the definition of f (see Eq. (3.49) ) and ψ (see Eq.(3.50)).

0 f (s + y -x)ds = ˆ0 y-x f (s)ds = ˆ∞ y-x f (s)ds, ˆx-y 0 f (s -x)ds = ˆ-y -x f (s)ds < ˆ0 -x f (s)ds = ˆ∞ -x f (s)ds. ( 3
In view of the statement of Lemma 3.3, χ(x, y) is non-increasing function of x on (0, ∞) and is strictly decreasing function of x on [y, ∞). Consequently, ´∞ 0 χ(x, y)e -y dy is a strictly decreasing function, and therefore, φ(x) is strictly decreasing as follows from (3.65).

4 Dependency of the speed of propagation on parameters: asymptotics and numerics.

In the previous sections, we established the uniqueness of a solution for problem (2.8). Our main result states that for an arbitrary fixed θ ∈ (0, 1), Λ ∈ (0, ∞) and α ∈ (0, 1) there exists a unique pair (c * , R * ) for which this problem admits a solution. This pair represents the velocity of propagation and width of the reaction zone for this set of parameters. The natural question is then to investigate the dependency of these quantities on the parameters of the problem, that is, dependencies c * (θ, Λ, α) and R * (θ, Λ, α). Of a particular interest is how the velocity of propagation changes with the parameters of the problem as this quantity is the main characteristic of the flame front.

As we showed in the previous sections, (c * , R * ) are uniquely determined from the solution of problem (2.16) that represents scaled distribution of the deficient reactant over the reaction zone. Hence, we start with a discussion of qualitative properties of the solution of (2.16) . While the solution of this problem can not be obtained in the closed form, it can be computed numerically. Figure 3 depicts function w for different values of α and Λ = 1. We note that solutions of (2.16) for Λ = 1 can be obtained from solutions of this equation with Observe that solutions for (2.16) for different values of α are not ordered. Indeed, for fixed |x| > 0 sufficiently small, w(x|Λ, α) is a decreasing function of α. Whereas for |x| sufficiently large, w(x|Λ, α) is increasing function of α. Formal considerations (which can be made rigorous) show that asymptotic behavior of the solution of (2.16) is as follows:

w(x) = (1 -α) 2 2Λ(1 + α) 1 1-α (-x) 2 1-α (1 + o(1)) |x| 1, w(x) = (1 -α) 1 1-α (-x) 1 1-α (1 + o(1)) |x| 1. (4.2)
That is, as α increases, w gets flatter and flatter near the origin and grows faster and faster for large values of |x|. Another important observation is that the solution (2.16) near the origin is heavily influenced by the specific value of the Lewis number as in this region the solution is obtained (in the first approximation) by balancing diffusion and reaction whereas the transport term is negligible. In contrast, far from the origin, the key ingredients are transport and reaction while diffusion is negligible in this region. Consequently, the solution of (2.16) is essentially independent of the Lewis number when |x| 1.

Let us also discuss the limiting behavior of the function w as α → 0 and α → 1. It is straightforward to verify that as α → 0 the solution of (2.16) approaches, on a compact sets, to the limiting profile

w 0 (x) = -x + Λ exp x Λ -1 , (4.3)
that verifies the limiting problem Λw 0 -w 0 -1 = 0, x < 0, w 0 (0) = 0, w 0 (0) = 0. (

We now claim that as α → 1 the function w approaches zero on compact sets. Indeed, let η = w -Λw . In view that w(x) > 0 and w (x) < 0 on (-∞, 0) and w(0) = w (0) = 0, we have η(x) > 0 on (-∞, 0) as η(0) = 0. By (2.16) we then have

-η = w α , x < 0, η(0) = 0. (4.5)
Taking into account that η, w are positive on (-∞, 0) we the have

-η = w α ≤ η α (4.6)
Integrating this inequality and taking into account the initial condition we have

η(x) ≤ [(1 -α)(-x)] 1 1-α , x ≤ 0 (4.7)
Since w < η we conclude

w(x) ≤ [(1 -α)(-x)] 1 1-α , x ≤ 0. (4.8)
Fixing x < 0 and taking a limit α → 1 in the inequality above, we conclude that w → 0 as α → 0 on compact sets. Now we turn to the evaluation of the pair c * (θ, Λ, α), R * (θ, Λ, α). An algorithm for finding this pair is rather straightforward. First, we fix Λ and α and solve problem (2.16). Then, we follow the procedure outlined in the proof of the main Theorem 2.1. Namely, using the solution of (2.16), we evaluate φ(x) given by (2.21) and a function given as:

ζ(x) = ˆ0 -x w α (s)ds 1-α 2 . ( 4.9)
We next fix θ and find a unique number σ * > 0 such that

φ(σ * ) = θ. (4.10)
The existence and uniqueness of such σ * follows from Proposition 2.1. Hence as follows from (2.24)

c * = ζ(σ * ), (4.11)
and then by (2.25)

R * = σ * /c * . (4.12)
The function φ(x) is decreasing, and the function ζ(x) is increasing on x ∈ (0, ∞). Moreover, by (4.1) we have x as α → 0 which follows directly from the fact that w approaches to w 0 in this limit. We now claim that φ(x) approaches unity and ζ(x) approaches zero on any fixed compact subset of x ∈ (0, ∞) as α → 1. To see that let us observe first that by (2.16) and (2.21) after integration by parts we have:

ζ(x|Λ, α) = √ Λζ x Λ Λ = 1, α .
φ(x) = 1 -φ(x), (4.14) with φ
(x) = Λ + (1 -Λ) ´0 -x w(s) exp(-(s + x))ds/w(-x) 1 + Λ|w (-x)/w(-x)| ≤ Λ + |1 -Λ| 1 + Λ|w (-x)/w(-x)| , x > 0, (4.15)
where the last inequality follows from the monotonicity of w. Next by (2.16) and (3.44) we obtain w (x)w(x) = 1 Λ (w (x) + w α (x)) w(x) ≤ (w (x)) 2 , x < 0. (4.16) Dividing the expression above by w 2 and taking into account that w < 0 we have

Λ w (x) w(x) 2 + w (x) w(x) ≥ 1 w 1-α (x)
, x < 0. (4.17)

Combining this observation with (4.8) we have

Λ w (x) w(x) 2 + w (x) w(x) ≥ 1 (1 -α)(-x) , x < 0. (4.18)
Fixing x < 0 and taking a limit as α → 1 in the expression above, we conclude that |w /w| → ∞ as α → 1 on compact sets. This observation implies that φ(x) → 0 as α → 1 and hence φ(x) → 1 in this limit. Finally by (4.8) and (4.9) we have

ζ(x) ≤ (1 -α) √ x. (4.19)
Taking a limit as α → 1 for x > 0 fixed we obtain ζ → 0 as α → 1 on compact sets. We now will discuss the dependency of φ, ζ of α. Figure 5 depicts functions φ and ζ for several values of Λ and α = 1/2. For a fixed x ∈ (0, ∞), the function φ is increasing in Λ, and the function ζ is decreasing in Λ.

One obvious observation that follows from the monotonicity of φ is that σ * is a decreasing function of θ which together with the monotonicity of ζ immediately implies that c * is a decreasing function of θ. This result is very natural in a physical context of the problem as an increase in ignition temperature decreases the speed of propagation.

Let us now discuss the behavior of φ(x), ζ(x) near the origin which will allow us to obtain asymptotic expressions for (c * , R * ) for ignition temperatures near unity. This regime is of particular interest as flame fronts are known to become unstable when ignition temperature approaches one. Discussion of flame front instabilities in this regime for the cases of zero and first order kinetics can be found in [START_REF] Brailovsky | Diffusive-thermal instabilities in premixed flames: Stepwise ignition-temperature kinetics[END_REF].

The behavior of functions φ(x) and ζ(x) for small values of x can be reconstructed from asymptotic formula (4.2). After rather tedious but straightforward computations, we have:

φ(x|Λ, α) ≈ 1 - 1 -α 2 x, ζ(x|Λ, α) = (2Λ) -α 2 (1 -α) 1+α 2 √ 1 + α x 1+α 2 for x 1. (4.20)
This observation immediately implies that for ignition temperatures close to unity we have:

c * (θ, Λ, α) 2 1 + α Λ -α 2 (1 -θ) 1+α 2 , R * (θ, Λ, α) 2(1 + α) 1 -α Λ α 2 (1 -θ) 1-α 2 , |1 -θ| 1.(4.21)
Direct verification shows that, in this regime, c * is a decreasing function of both α and Λ, whereas R * is an increasing function of both of these parameters. Moreover, R * is decreasing function of θ. Figure 6 depicts dependency of c * and R * on α for several values of Λ with θ = 0.98, and Figure 7 shows the dependency of the velocity c * on α and Λ. These figures were generated using asymptotic formulas (4.21) which are extremely close to their numerical counterparts . Moreover, one can verify that c * and R * given by (4.21) fully reproduce the limiting behavior of these functions in the limits α → 0 and α → 1. Indeed, when α = 1, the velocity of propagation and width of the reaction zone are given by [

2] c * (θ, Λ, α = 1) = θ 1 -θ + Λ θ 1 -θ 2 -1 2 , R * (θ, Λ, α = 1) = ∞. ( 4
c * (θ, Λ, α = 1 -θ √ Λ , c * (θ, Λ, α = 0) 2(1 -θ) for |θ -1| 1. (4.25) Hence c * (θ, Λ, α → 1) → c * (θ, Λ, α = 1), R * (θ, Λ, α → 1) → ∞, c * (θ, Λ, α → 0) → c * (θ, Λ, α = 0), R * (θ, Λ, α → 0) → R * (θ, Λ, α = 0). ( 4

.26)

Now consider the regime of intermediate values of θ, we study this regime numerically. Figures 8,9, 10 depict dependency of the velocity of propagation and width of the reaction zone on the reaction order α for Λ = 0.2, 1, 5 and θ = 0.75, 0.5, 0.25, and Figure 11 shows dependency of the velocity of propagation on α and Λ for θ = 0.5. According to numerics for intermediate values of θ, the character of the dependency of c * and R * on α and Λ remains similar to the one for θ near unity. However, dependency of c * on both α and Λ becomes weaker as θ decreases. The dependency of R * on α in this regime is still very strong, but dependency on Λ becomes weaker as θ decreases. When θ becomes sufficiently small, the asymptotic behavior of φ(x) and ζ(x) can again be recovered from the asymptotic formula (4.2) and reads:

φ(x|Λ, α) ≈ 1 1 -α 1 x , ζ(x|Λ, α) = √ 1 -α √ x for x 1. (4.27)
Therefore, in this regime we have:

c * (θ) ≈ 1 √ θ , R * (θ) ≈ 1 (1 -α) 1 √ θ , θ 1.
(4.28)

Consequently in this regime, the velocity of propagation (in the first approximation) depends exclusively on the ignition temperature θ, and the reaction width is independent of Λ and is an increasing function of α. As in the regime θ near unity, the velocity of propagation and width of the reaction zone approaches to formal limits as α → 0 and α → 1. The discussion above strongly suggests that the velocity of propagation decreases with the increase of the reaction order. Consequently, the velocity of propagation with the reaction order α ∈ (0, 1) is bounded from below by the velocity of propagation with the reaction order unity and from above by the velocity of propagation with zero reaction order. This, far from obvious observation, is quite in line with the physical intuition as an increase of the reaction order decreases the reaction rate which, in turn, slows the flame front. Another observation which is less surprising is that the increase of the molecular diffusivity with regard to the thermal diffusivity decreases the speed of the flame front. This is clearly the case for the reactions of first order but remains true for reaction order α ∈ (0, 1).

We hence formulate the following:

Conjecture 4.1. The velocity of propagation c * (θ, Λ, α) is a decreasing function of all of its arguments whereas the reaction width R * (θ, Λ, α) is a decreasing function of θ and an increasing function of Λ and α.

Appendix

In this appendix we show that inequality (3.44) is equivalent to (3.47) and briefly discuss solution of problem (2.16) from dynamical systems point of view.

In what follows we denote the vector field in the right hand side of (3.46) by X c . The unique critical point of this vector field in Q is the origin O = (0, 0). System (3.46) was thoroughly investigated in [3, Section 4] using a weak version of Poincaré-Bendixson theorem (see [START_REF] Roussarie | Some Applications of the Poincaré-Bendixson Theorem[END_REF]). The following proposition summarizes the results of Lemmas 4.7 and 4.8 of [START_REF] Brauner | Existence of a traveling wave solution in a free interface problem with fractional order kinetics[END_REF]. Proposition 5.1. There exists a unique global stable manifold at the origin given by a trajectory of X c in Q converging toward the origin when t → 0 -such that the orbit defined by this trajectory is the graph of a function p = p (q) with q ∈ (0, +∞). Moreover, p is analytic for q > 0 and extends continuously at 0 with the value p (0) = 0. The trajectory for system (3.46) can be obtained numerically. Figures 12 and13 depict such a trajectory on the phase plane for Λ = 1, α = 1/2 in cartesian and polar coordinates respectively. Qualitative behavior of the trajectory is similar for other values of parameters Λ, α.

Figure 1 :

 1 Figure 1: Sketch of a traveling front solution for system (1.1) with α ≥ 1. Flame (red) and deficient reactant (blue) fronts propagating from burned ((u, v) = (1, 0)) to unburned ((u, v) = (0, 1)) states at ξ → ∞ and ξ → -∞ respectively with the constant speed c > 0. The arrow indicates the direction of propagation.The position of an ignition interface where the temperature is equal to an ignition one (u = θ) is indicated by ξ ign .

Figure 2 :

 2 Figure2: Sketch of a traveling front solution for system (1.1) with 0 < α < 1. Flame (red) and deficient reactant (blue) fronts propagating from burned ((u, v) = (1, 0)) to an unburned ((u, v) = (0, 1)) states at ξ → ∞ and ξ → -∞ respectively with the constant speed c > 0. The arrow indicates the direction of propagation.The position of an ignition interface where the temperature is equal to an ignition one (u = θ) is indicated by ξ ign . The position of the trailing interface, the point to the right of which temperature and concentration of deficient are identically one and zero (u = 1, v = 0), respectively is indicated by ξ tr .

. 10 )

 10 Sending x → ∞ in this inequality and taking into account (3.8) we obtain the second equation in (3.1). Lemma 3.2. Let Y be a non-negative solution of ΛY -Y -Y α = 0. (3.11)

  .25) Note that the existence of such x M is guaranteed since Y (x 0 ) = Y 0 and Y (x 0 ) < 0.We claim that for any solution of (3.11) satisfying (3.12) the following inequalities holdY (x) > Z(x), Y (x) > Z (x),(3.26) for x ∈ I. Here Z is the solution of (3.20), (3.21) with arbitrary

  .35) By (3.11) for x > x M we have

YY

  (x) > Y 0 on (x M , x M + δ]. (x M + δ) > Y 0 , Y (x M + δ) > 0. (3.40) Now assume that x *M is the first point greater than x M + δ such that Y (x * M ) = 0, but by (3.36) and(3.40)

  .55) and thus by Corollary 3.1 g (x) ≤ 0 on (-∞, 0). (3.56)

f

  .62) Hence integrating(3.61) in s from 0 to x -y and taking into account (3.62) we conclude that (3.53) holds for arbitrary x > y > 0 and hence ψ(x, y) is indeed strictly increasing in x for x > y > 0.Now we are ready to give a proof of the propositionProof of Proposition 2.1. Letχ(x, y) = 1 -ψ(x, y) = ´y-x -x f (s)ds ´∞ -x f (s)ds . (s)ds e -y dy = ˆ∞ 0 ˆy 0 f (s -x)e -y ds dy = ˆ∞ 0 ˆ∞ s f (s -x)e -y dy ds = ˆ∞ 0 ˆ∞ s e -y dy f (s -x)ds = ˆ∞ 0 f (s -x)e -s ds = ˆ∞ -x f (s)e -(s+x) ds. y)e -y dy = ´∞ -x f (s)e -(s+x) ds ´∞ -x f (s)ds = ´0 -x w α (s)e -(s+x) ds ´0 -x w α (s)ds = φ(x).(3.65)

Figure 3 :

 3 Figure 3: Numerical solutions of problem (2.16) for α = 1/4 (blue), α = 1/2 (orange) and α = 3/4 (green).

(4. 13 )Figure 4 :

 134 Figure 4: Functions φ and ζ for α = 0.75 (blue), α = 0.5 (orange) and α = 0.25 (green) and Λ = 1. The arrow indicates direction of increase of α.

Figure 5 :

 5 Figure 5: Functions φ and ζ for Λ = 5 (blue), Λ = 1 (orange) and Λ = 0.2 (green) and α = 1/2. The arrow indicates direction of increase of Λ.

  .22) whereas when α = 0, we have[START_REF] Brailovsky | Diffusive-thermal instabilities in premixed flames: Stepwise ignition-temperature kinetics[END_REF] c * (θ, Λ, α = 0) = R * (θ, Λ, α = 0

Figure 6 :

 6 Figure 6: Dependency of the velocity of propagation c * and width of the reaction zone R * on α for Λ = 0.2 (blue), Λ = 1 (orange) and Λ = 5 (green) for the ignition temperature θ = 0.98. The arrow indicates direction of increase of Λ.

Figure 7 :

 7 Figure 7: Dependency of the velocity of propagation c * on α and Λ for θ = 0.98. Left panel represents a three dimensional plot of c * (θ = 0.98, Λ, α) and the right plane depicts level sets of this function

Figure 8 :

 8 Figure 8: Dependency of the velocity of propagation c * and width of the reaction zone R * on α for Λ = 0.2 (blue), Λ = 1 (orange) and Λ = 5 (green) for the ignition temperature θ = 0.75. The arrow indicates direction of increase of Λ.

Figure 9 :

 9 Figure 9: Dependency of the velocity of propagation c * and width of the reaction zone R * on α for Λ = 0.2 (blue), Λ = 1 (orange) and Λ = 5 (green) for the ignition temperature θ = 0.5. The arrow indicates direction of increase of Λ.

Figure 10 :

 10 Figure 10: Dependency of the velocity of propagation c * and width of the reaction zone R * on α for Λ = 0.2 (blue), Λ = 1 (orange) and Λ = 5 (green) for the ignition temperature θ = 0.25. The arrow indicates direction of increase of Λ.

Figure 11 :

 11 Figure 11: Dependency of the velocity of propagation c * on α and Λ for θ = 0.5. Left panel represents a three dimensional plot of c * (θ = 0.5, Λ, α) and the right plane depicts level sets of this function.

Figure 12 :

 12 Figure 12: Phase portrait for system (3.46) with Λ = 1 and α = 1/2 in cartesian coordinates (q, p).

Figure 13 :

 13 Figure 13: Phase portrait for system (3.46) with Λ = 1 and α = 1/2 in polar coordinates (ϑ, r).

Acknowledgments. The work of AM and PVG was supported in a part by US-Israel BSF grant 2020005. PVG would like to thank Fedor Nazarov for multiple valuable discussions and substantial help with proving the main result of this paper.

Next observe that by (3.47) we have tan(θ(t)) = p(t) q(t) .

(5.1)

In view of Proposition 5.1, the trajectory (q, p)(t) is smooth and hence (differentiating the expression above) we have θ(t) = q(t) ṗ(t) -p(t) q(t) r 2 (t) .

(5.2)

Since

Declaration of Competing Interest:

The authors declare that they have no competing interests.

Data availability: No data was used for the research described in the article.


Noudéhouénou L J Houssou

Jean-Loup Guillaume

Armelle Prigent
email: armelle.prigent@univ-lr.fr

Edit distance with Quasi Real Penalties: a hybrid distance for network-constrained trajectories

Keywords: Trajectory analysis, edit distance, road network, clustering, graphs

In this paper, we propose a new distance for network-constrained trajectories named Edit distance with Quasi Real Penalties (EQRP). Depending on the case, it can compare trajectories as non-ordered sets and as sequences while other distances only compare trajectories as non-ordered sets or as sequences. Moreover, it is parameter-free, manages local time shifting, and respects triangle inequality; three properties expected from a trajectory distance that are not satisfied simultaneously by any other distance to the best of our knowledge. To demonstrate the pertinence of our idea, we benchmark our distance against some state-of-the-art distances for networkconstrained trajectories. Specifically, for each distance, we determine its capability to identify precisely similar trajectories. We also determine their respective performance for trajectory clustering. Our results show the predominance of EQRP over the existing edit distances and in some cases a more precise ability to evaluate the dissimilarity between network-constrained trajectories compared to other measures.

I. INTRODUCTION

Location-based services (LBS) are widely used in our modern societies to provide assistance in many contexts such as social interactions, navigation, entertainment, or fitness activities for example. This results in massive production and storage of spatio-temporal traces also called trajectories, which mainly consist of ordered sequences of locations and describe their authors' movements in time and space. We distinguish two kinds of trajectories: free-space or nearly free-space ones which can take any form and occur in an environment with no or few limitations (e.g. on the sea) and networkconstrained trajectories corresponding to travel that only take place along the edges of a given network such as a road network (see Figure 1). Compared to free-space trajectories, network-constrained trajectories are more related to human activities and have broader socio-economic applications. They are for example essential to understand mobility patterns in cities, predict traffic congestion, select the best advertisement and business locations or recommend optimal travel routes.

Initially, network-constrained trajectories were processed without taking into account the underlying network structure as their free-space counterpart. This is because most methods and measures dedicated to trajectory analysis have been originally designed for free-space trajectories. However, such an approximation can induce significant bias in trajectories comparison. For example, trajectories can be mistakenly declared similar because some of their points appear spatially close while the network distance between them is actually large. To overcome these biases, a bunch of measures specially developed for network-constrained trajectories have been recently proposed in the literature. We can classify them into four main families: simple distances, node-to-trajectory distances, set distances, and warping distances. Simple distances compare trajectories by evaluating the distance between nodes of the same index (i.e. comparing first node of each trajectory, then second node, and so on) or using a set of points of interest. Node-to-trajectory distances determine the minimum spacing between each node of a first trajectory and all the nodes of a second one. Set distances compare two trajectories based on the size of their intersection, i.e. the number of nodes or edges they have in common. Finally, warping distances take into account the order of the nodes and try to find the optimal alignment between trajectories when comparing them. This brings an additional dimension to the trajectories comparison process which is not present in the three other families. While warping distances consider trajectories as sequences of nodes and the others mostly as non-ordered sets of nodes, trajectories actually bear both characteristics and should be compared along these two dimensions. For example, a trajectory and its copy in the opposite direction have exactly the same nodes but are not identical, without being completely different. However, to the best of our knowledge, none of the existing distances are able to process trajectories both as non-ordered sets and sequences. In addition, we compared state-ofthe-art distances by checking three major properties: triangle inequality, local time shifting handling, and parameter dependence. Triangle inequality ensures that the distance measure discriminates correctly the compared trajectories. It is also used to reduce the number of distance computations through pruning strategies when comparing large sets of trajectories. Local time shifting corresponds to the presence of very similar sub-sequences, offset from each other, in different trajectories. A measure able to handle it can optimally align trajectories during their comparison. Concerning parameter dependence, we are actually evaluating if the measure is parameter-free or not as it is more suitable to have no parameters to calibrate before applying a measure. None of the existing distances satisfy simultaneously these three properties even though the shortcomings of some distances are compensated by other distances and vice-versa. These observations led us to formulate the following hypothesis: is it possible to hybridize two distances to obtain a measure that respects the three forementioned properties but also treats trajectories as non-ordered sets and sequences concomitantly? We respond positively to this question by proposing Edit distance with Quasi Real Penalties (EQRP) a new network-constrained trajectory distance that combines warping and node-to-trajectory distances.

The remaining of the paper is structured as follows: to support our motivations, we start by explaining some notions related to network-constrained trajectories. Then we present the related works. After, we define our new distance EQRP and describe the methodology adopted for the comparison of some state-of-art network-constrained distances to the newly proposed one. Later, we present the implementation details of our experiments, the results obtained and a discussion of these results before concluding.

II. DEFINITIONS

A road network is a set of interconnected road segments that support mobility in a given geographical space. Usually, it is modeled as a graph G(V, E) where V represents the set of nodes which are often road intersections and E corresponds to the set of edges which are road segments linking the nodes.

A trajectory is a time-ordered sequence of geographic positions, usually identified by latitudes and longitudes, which describe a movement. In the context of a network-constrained environment such as an urban road network, the geographic positions can be replaced by the nodes of the graph representing the road network. In the rest of the paper, we consider trajectories as sequences of road network nodes. We formally define a trajectory T as:

T = ((v 1 , t 1 ), ..., (v i , t i ), ..., (v n , t n )) (1)
Where v i ∈ V are nodes and t 1 < ... < t i < ... < t n are timestamps associated to nodes. Timestamps can be omitted in various situations.

III. RELATED WORK

Although trajectory data analysis has been an active research field for decades, similarity/distance measures specifically designed for network-constrained trajectories emerged only recently. One of the first distance integrating road network structure in its definition was proposed in 2005 by Hwang et al. [START_REF] Hwang | Spatio-temporal similarity analysis between trajectories on road networks[END_REF]. It compares trajectories according to a set of Points Of Interest (POIs) which can be road intersections or places. Two trajectories are considered similar if all the POIs lie on both of them otherwise there are not. Thus, this distance measure, which does not respect triangle inequality, can only take two values 0 or 1, and is not able to compare trajectories accurately. Moreover, POIs must be predefined by users who may not have extended knowledge of the road network. Following this first distance, Tiakas et al. [START_REF] Tiakas | Searching for similar trajectories in spatial networks[END_REF] presented a distance for network-constrained trajectories, similar to the Euclidean distance for free-space trajectories. This distance corresponds to the average shortest path distances between nodes of same index (i.e nodes that appear at the same position in compared trajectories). It respects triangle inequality, but cannot handle local time shifting or trajectories with different lengths.

To ensure that, no matter the size of trajectories, it is always possible to compute a network distance between them, Evans et al. [START_REF] Evans | Summarizing trajectories into k-primary corridors: a summary of results[END_REF] opted for a node-to-trajectory approach. Here, the distance between two trajectories is equal to the average of minimum distances between each node of one trajectory and all the nodes of the second. In the same vein, [START_REF] Shang | Trajectory similarity join in spatial networks[END_REF] defined a spatio-temporal similarity measure by converting the node-totrajectory distance to a similarity measure using an exponential function. This second measure has been labeled as Spatiotemporal Linear Combine STLC in [START_REF] Su | A survey of trajectory distance measures and performance evaluation[END_REF]. Notice that both nodeto-trajectory measures do not respect triangle inequality.

Another node-to-trajectory distance is the Hausdorff distance. It measures the mutual proximity of two sets, by indicating the maximal distance from any node in one set to the other set. [START_REF] Roh | Nncluster: An efficient clustering algorithm for road network trajectories[END_REF] and [START_REF] Evans | Fast and exact network trajectory similarity computation: a case-study on bicycle corridor planning[END_REF] proposed versions of Hausdorff distance adapted to network constrained trajectories by replacing Euclidean distance with network distance. The version in [START_REF] Roh | Nncluster: An efficient clustering algorithm for road network trajectories[END_REF] respect triangle inequality contrary to the one in [START_REF] Evans | Fast and exact network trajectory similarity computation: a case-study on bicycle corridor planning[END_REF].

Besides, we also have distances based on the intersection of sets. They are used for applications like carpooling or network flow analysis where similar trajectories are those which share common road segments. Usually, these distances are computed edge-wise to ease the detection of common road segments and to integrate the length of road segments. As examples, we mention SIM TRAJ [START_REF] Kharrat | Clustering algorithm for network constraint trajectories[END_REF], a similarity that equals to the ratio between the length of two trajectories' common parts and the length of one of them selected as reference, DSL [START_REF] Won | Trajectory clustering in road network environment[END_REF] a dissimilarity which corresponds to the length of two trajectories uncommon parts over the total length of all road segments forming them and EBD [START_REF] Wang | Fast large-scale trajectory clustering[END_REF] equivalent to the length of the longest trajectory minus the length of the compared trajectories' intersection. Among the three, only EBD respect triangle inequality.

The previously presented distances are all sequence-independent, i.e. they don't consider the order of trajectories' nodes, contrary to warping measures. Warping measures for network-constrained trajectories are mainly extensions of measures for free-space trajectories such as: Longest Common Subsequence (LCSS) [START_REF] Vlachos | Discovering similar multidimensional trajectories[END_REF], Dynamic Time Warping (DTW) [START_REF] Berndt | Finding patterns in time series: a dynamic programming approach[END_REF], Edit Distance on Real sequence (EDR) [START_REF] Chen | Robust and fast similarity search for moving object trajectories[END_REF] and Edit distance with Real Penalty (ERP) [START_REF] Chen | On the marriage of lp-norms and edit distance[END_REF]. They are sequence-dependent and able to handle local time shifting and trajectories with different lengths. Shang et al. [START_REF] Shang | Personalized trajectory matching in spatial networks[END_REF] proposed the first warping measure adapted to network-constrained trajectories. Their measure extends LCSS which compares two trajectories by determining the size of their longest common subsequence. The points constituting this sub-sequence must not necessarily occupy consecutive positions but must appear in the same order within the two trajectories. The particularity of LCSS is to match two points, even if they do not coincide, as long as the distance between them is below a threshold. Later, [START_REF] Wang | Torch: A search engine for trajectory data[END_REF] defined LORS (Longest Overlapping Road Segments) which is also a variation of LCSS using the sequence of edges. In addition to the trajectory modeling difference, LORS matches edges only when they are identical while the former measure depends on a threshold value to loosely match nodes.

Then, [START_REF] Yuan | Distributed in-memory trajectory similarity search and join on road network[END_REF] built LCRS (Longest Common Road Segments) on top of LORS using Jaccard index. Note that the idea of using the Jaccard index to define trajectory similarity was first presented by Xia et al. in [START_REF] Xia | Spatiotemporal similarity measure for network constrained trajectory data[END_REF]. Network versions of DTW [START_REF] Shang | Dita: Distributed in-memory trajectory analytics[END_REF], EDR and ERP [START_REF] Koide | Fast subtrajectory similarity search in road networks under weighted edit distance constraints[END_REF] were also proposed. These last three distances have the advantage to be complete match warping measures, i.e. all the nodes of compared trajectories are used in the calculation of the distance, contrary to LCSS based distances which are partial match distances. A match corresponds to a pairing between two nodes, each belonging to one of the compared trajectories. The distance between these matching nodes is used to compute the global distance between the trajectories. Concerning DTW, its particularity is to allow sequence stretching or shrinking in time by replicating previously processed points when there is no match. This ensures that there is always a matching pair for any trajectory point. EDR and ERP are edit distances that evaluate the cost of edit operations namely, insertion, deletion, and substitution necessary to transform a trajectory into another. Substitution corresponds to a match while insertion and deletion occur when there is no match which is also called a gap. The cost of each edit operation is designated as a subcost.

Despite their qualities, edit distances also bear shortcomings. EDR, for example, depends on a matching threshold which is often difficult to fix. It also uses fixed subcosts that are not correlated to the actual distances between nodes. This negatively impacts the ability to accurately assess the distances between trajectories. ERP, for its part, is defined according to real distances between nodes but requires the specification of a reference node in the evaluation of the distances between trajectories to handle gaps. The impact of this reference node on ERP accuracy in graph space had not been studied. Moreover, the use of real distances between nodes makes ERP sensitive to the presence of deformations in the trajectories. Among warping distances, only ERP respect triangle inequality.

As a wrap-up, we retain that some network-constrained trajectory distances consider trajectories only as non-ordered sets while others treat them only as sequences although trajectories bear both properties. Moreover, despite their respective qualities, none of these distances concomitantly respect triangle inequality, is parameter-free, and able to manage local time shifting, even though each of these characteristics is important to guarantee a good comparison of trajectories. These are the reasons that led us to define a new measure that meets all these needs.

IV. EDIT DISTANCE WITH QUASI REAL PENALTIES

The measure we propose differs from existing edit distances by how it handles gaps during the matching process. More precisely, when there is a match, the classical network distance between pair of nodes is used. However, if there is a gap, the insertion or deletion cost is fixed at the Hausdorff distance between the compared trajectories. Let's recall that Hausdorff distance measures the mutual proximity of two sets, by indicating the maximal distance from any node in one set to the other set. Doing so, we replace the fixed cost of EDR (always equal to 1) by a semi-fixed cost, i.e. which is fixed only for a pair of trajectories. This also allows us to avoid the use of a reference node contrary to ERP. The result is an edit distance that considers trajectories both as sequences and non-ordered sets of nodes (thanks to Hausdorff distance). This definition also makes it possible to indirectly set the matching threshold between nodes to the Hausdorff distance between trajectories. The new edit distance we introduce is named Edit distance with Quasi Real Penalties (EQRP). It is formally defined between two trajectories T 1 and T 2 of respective sizes n and m by the recurrence formula:

EQRP (T1, T2) =                0, If n = 0 and m = 0 ∞, If n = 0 or m = 0 min      EQRP (Rest(T1), Rest(T2)) + dG(Head(T1), Head(T2)) + D EQRP (Rest(T1), T2) + Hausdorf f (T1, T2) + D EQRP (T1, Rest(T2)) + Hausdorf f (T1, T2) + D Otherwise (2)
Head designates the first node of the considered sequence, Rest refers to the whole sequence except the first node and d G is the network distance between two nodes, and D is the diameter of the network. The Hausdorff distance formula is :

Hausdorf f (T 1 , T 2 ) = max max vi∈T1 min vj ∈T2 d G (v i , v j ) max vj ∈T2 min vi∈T1 d G (v i , v j )
(3) Adding the diameter to the different subcosts ensures that EQRP respects the triangle inequality: for any triplet of distinct trajectories R, S, and T we have

EQRP (R, S) ≤ EQRP (R, T ) + EQRP (T, S) (4)
Indeed, for an edit distance to respect this inequality, it is sufficient that the distance between elements on which it is based also respects it even if there is a gap [START_REF] Waterman | Some biological sequence metrics[END_REF]. When there is no gap, we use the graph distance (the shortest path between nodes) which respects the triangle inequality. Conversely, if there is a gap, the inequality remains respected because the Hausdorff distance is necessarily smaller than or equal to D.

As other warping distances, EQRP can be calculated using a dynamic programming algorithm to solve its recurrence formula. This method of calculation has a complexity of O(nm(|V |+|E|)) with n and m, the respective sizes of the compared trajectories and (|V |+|E|) the cost of computing the shortest path between two nodes. This complexity includes the cost of pre-calculating the Hausdorff distance whose cost is also O(nm(|V |+|E|)).

V. COMPARISON TO OTHER DISTANCES

After having formally defined our new distance, it is important to ensure its relevance by comparing it to existing network-constrained distances. However, comparing distances is not trivial because the values returned by different distances for the same pair of trajectories are not always consistent as each distance has its own semantics. Therefore, it is not appropriate to rely directly on the values and we need to proceed through indirect approaches. In this section, we describe two alternative approaches: resistance to transformations and clustering ability.

A. Comparison by transformations

The method of comparison by transformations was initially proposed by Han et al in their review [START_REF] Su | A survey of trajectory distance measures and performance evaluation[END_REF]. It aims to compare distances between free-space trajectories by evaluating their ability to resist different transformations: addition or deletion of sampling points; trajectory resampling; time dilation and compression; spatial expansion and compression; noise addition. Here, we do not explicitly take into account the temporal data of the trajectories, because most of the compared distances do not integrate them, so the transformations relating exclusively to the time domain are discarded (time dilation and compression, trajectory resampling). Additionally, the trajectories are modeled as sequences of nodes, it is subsequently not possible to add sample points without modifying the road network. The constraints of the road network also make it very difficult to expand or compress the trajectories without changing their shape. Due to these limitations, we retain only two of the transformations initially proposed and adapt them to network-constrained trajectories. We renamed the deletion of sample points as nodes deletion and we decline the noise addition in two versions: loops and detours addition.

1) Nodes deletion: Before any node deletion, we identify first the pivot nodes which are nodes where the trajectory noticeably changes shape. This is done using the Ramer-Douglas-Peuchker (RDP) algorithm [START_REF] Heckbert | Survey of polygonal surface simplification algorithms[END_REF] whose principle is to replace a polyline (a broken line) made up of several points by a direct line, i.e. a shortcut, between the endpoints. This modification is only applied if the distance between the direct line and the furthest point of the polyline is less than a predefined threshold. To adapt this idea to network-constrained trajectories, we add the GPS coordinates as attributes to the nodes. Once pivot nodes are identified, we randomly delete some nodes among those remaining. The number of nodes removed depends on a ratio r whose values are between 0 and 1. For a trajectory T , of size n, with k pivot nodes, the number of nodes to remove is N B pts = min{n * r, n -k}. The min ensures that only the non-pivot nodes can be deleted.

2) Loops and detour addition: In free-space, noise addition to a trajectory results in adding Gaussian noise to points' coordinates, while for network-constrained trajectories, it corresponds to the addition of deviations in the route of the trajectories. The reason is that the space of the road network is discrete, which prevents the addition of points other than those present in the set of nodes forming the network. Two types of deviations can be generated: a loop or a detour. A loop is a circuit in the road network, a sequence of nodes starting and ending with the same node. A detour is an alternative route between two distinct nodes of the same trajectory. Usually, these two disturbances occur during the map-matching process when the raw trajectories contain noise or due to an error in the matching process between points and road segments. The objective here is to try to reproduce them artificially.

To add a loop, we randomly choose a node u in the trajectory and a node v outside the trajectory but within a predefined radius around u. We calculate the shortest paths from u to v and vice-versa. The resulting path, which forms a loop, is then inserted in the trajectory: once arrived at u, we go through the loop before continuing the initial trajectory. The number of loops to add to a trajectory is also managed by a ratio between 0 and 1.

The procedure to add a detour is similar. First, we choose two distinct nodes u and v of the same trajectory separated by a given number of nodes p (for projection). Then, we select a random node w outside the trajectory but accessible from u within a predefined radius. Then we determine the shortest path going from u to v via w and we replace the sub-sequence of the original trajectory between u and v by the one passing through w. Note that the detour is not necessarily longer than the path it replaces in the original trajectory, because the sub-trajectory from u to v had no reason to be the shortest path between these two nodes. The detour is therefore just an alternative path passing through the vertex w.

3) Evaluation of transformations impact: To quantify the robustness of distances against transformations, we rely on the method described in [START_REF] Su | A survey of trajectory distance measures and performance evaluation[END_REF]. First, we randomly select a reference trajectory from the original dataset. Second, we extract the ordered list of trajectories closest to the reference trajectory, using the distance we want to evaluate, sorted by increasing value. This list is the original list. Third, we apply a given transformation to all trajectories in the original dataset apart from the reference trajectory, to obtain a derived dataset. We then order the trajectories of the derived dataset by increasing distance to the reference (and untouched) trajectory to obtain a transformed list. Fourth, we evaluate the similarity between the original list and the transformed list. This is repeated for each transformation. The intuition behind this approach is that a distance that resists well to transformations will tend to similarly order the closest trajectories to a reference trajectory, both in the original and in derived datasets. The similarity between the ordered lists is computed using Spearman's rank correlation coefficient [START_REF] Gauthier | Detecting trends using spearman's rank correlation coefficient[END_REF] that measures the correlation between two rankings X and Y of the same size as:

r s = cov(X, Y ) σ X .σ Y (5)
With cov the covariance and σ the standard deviation. Spearman's measure takes values from -1 to 1. A value close to 1 means that the rankings are similar which also indicates a good robustness of the distances against transformations. A value close to 0 indicates an absence of link between the two rankings and therefore a low, or non-existent, resistance of distances against transformations. A value close to -1 indicates a decreasing trend: close trajectories in the original data are distant in the derived data and vice-versa.

Comparing distances through correlations can be a source of bias in the case of measures such as EDR and LCSS because of the way they are calculated. Indeed, these measures consider any pair of trajectories having no matching pair of nodes as being at a maximum distance from each other, i.e. at distance 1 when distances are normalized. This poses a problem of resolution in the trajectories comparison and biases the nearest neighbor ordering. More specifically, we note that from a certain index, the order of the nearest neighbors does not vary anymore since all the trajectories are found at distance 1 from the reference trajectory regardless of the transformations. We then obtain high correlation values that do not reflect the ability of the distance to manage a transformation but rather its inability to do so. We could possibly choose to exclude the trajectories at a maximum distance but this problem also arises in a more general way when distances can only take a limited set of values, small compared to the number of trajectories, and that many trajectories are therefore at the same distance from the reference trajectory. To overcome this bias, we break the equality of the trajectories located at the same distance from the reference trajectory by a random choice. For a better understanding, consider six trajectories T 2; T 3; T 4; T 5; T 6; T 7 whose respective distances to the reference trajectory T 1 are: 0.3; 1; 0.5; 1; 1; 0.2. The natural ordering of the nearest neighbors of T 1 would be T 7; T 2; T 4; T 3; T 5; T 6 if we used the order of the indices in case of equality. If the distance used has a limited resolution as explained above, the order of the last three trajectories T 3; T 5; T 6 has very little chance to change because of transformations. We break this constant ordering by randomly ranking T 3; T 5; T 6 each time they are all at the maximum distance from the reference trajectory. In an extreme case, if all the trajectories are at distance 1 from the reference trajectory in the original dataset and in a derived dataset, this solution will lead to a total lack of correlation between the two rankings. For the sake of impartiality, we apply this random ranking technique to all distances.

Comparison of distances according to their robustness to transformations does not assess their resolution, i.e. their ability to accurately account for the relative differences be-tween trajectories. In order to better assess this resolution, we propose to also use a comparison by clustering.

B. Comparison by Clustering

A way to compare network-constrained trajectory distances is to measure how well they can cluster trajectories. The underlying idea is that if we have groups of similar trajectories (visually close for example) then the more discriminating a distance is, the better it should distribute the trajectories in their original clusters. To compare distances by clustering, we select a sample of trajectories that we divided into distinct groups. Each group consists of trajectories that are geographically close to each other and the clusters do not overlap. Trajectories are clustered using Hierarchical Cluster Analysis (HCA) algorithm. The comparison of the clustering results is done using the Adjusted Rand Index (ARI) measure [START_REF] Hubert | Comparing clusterings[END_REF] and the Silhouette coefficient [START_REF] Rousseeuw | Silhouettes: a graphical aid to the interpretation and validation of cluster analysis[END_REF].

ARI estimates the quality of clustering by comparing the result obtained with the ground truth while ignoring permutations. Indeed, it may happen that the trajectories are perfectly distributed in clusters but the order of these clusters differs from that provided by the ground truth. Such a result must be counted as correct, hence the need to ignore permutations between clusters. The other specificity of ARI measure is that random clustering returns a value close to 0 regardless of the dataset size and the number of clusters. The formula of the ARI measure is:

ARI = RI -E[RI] max(RI) -E[RI] ( 6
)
where RI is the Rand Index measure from which ARI is derived, RI = (a + b)/ n 2 ; a is the number of pairs of elements that belong to the same cluster in the ground truth and in the result of the clustering; b is the number of pairs of elements that belong to the same cluster in the ground truth and to different clusters in the clustering result; E[RI] is the expectation of the Rand Index and is obtained by calculating this measurement on random distributions of elements in clusters.

The Silhouette coefficient measures the similarity of an object with its own group compared to other groups. It does not require a ground truth and is defined as:

S = b -a max(a, b) (7)
With a the average distance between an element and all the other members of its cluster and b the average distance between an element and all the elements of the closest cluster. The Silhouette coefficient for all the elements of a dataset corresponds to the average of the Silhouette coefficients of the different elements. The values of the ARI measure and the Silhouette coefficient range from -1 (worst case) to 1 (best case). Negative values usually indicate the misassignment of elements to clusters. Values close to 0 indicate clusters overlapping for the Silhouette coefficient and random assignments for the ARI measure. The advantage of the Silhouette coefficient compared to the ARI measure is that it does not require ground truth.

VI. EXPERIMENTATION

Experiments were performed on a server running UBUNTU 16.04. Geographic data are processed with Quantum GIS 3.4 Madeira (QGIS 3.4).

A. Road network and trajectories

The road network is extracted from the map of the administrative region of Porto (Portugal) downloaded on Open-StreetMap 1 and using the open-source tool Osm2Po2 . The directed graph obtained is strongly connected with 350189 nodes representing intersections or road ends and 811583 edges.

Trajectories were extracted from Porto trajectory dataset [START_REF] Moreira-Matias | Predicting taxi-passenger demand using streaming data[END_REF] which contain 1710670 raw trajectories corresponding to 442 taxis' journeys in the city of Porto and its surroundings over a whole year, from 01/07/2013 to 30/06/2014. Each trip is described by an identifier, a set of GPS coordinates measured every 15 seconds, and various additional information about the taxi, the origin of the trip request, etc. The trajectories selected have been map-matched using barefoot 3 , an open-source library written in Java, and implementing the Hidden Markov Model (HMM). Each trajectory is then transformed into a sequence of road network nodes.

B. Distances

We compare EQRP to network versions of STLC, Hausdorff distance, DTW, LCSS, EDR and ERP in terms of clustering performances and transformations handling. Table I presents the formulas of the competing distances network versions which all have a complexity of O(nm(|V |+|E|)). Sizes of trajectories T 1 and T 2 are respectively n and m. Head designates the first node of the considered sequence, Rest refers to the whole sequence except the first node, v g is the reference node, and d G is the network distance between two nodes. To implement the network version of ERP, we had to define a reference node for the case when there is a gap. Among diverse possibilities, we choose the reference node based on the criterion of eccentricity. Recall that the eccentricity of a vertex is the maximum distance between this vertex and all the other vertices of the graph. According to this choice, two variants of the graph version of ERP were implemented. In the first variant (ERP G1 ), the reference node is one of those having the minimum eccentricity, and in the second (ERP G2 ) it is a node of maximum eccentricity which is considered. A node of minimum eccentricity naturally constitutes a reference in a graph since it belongs to its center. Likewise, a maximum eccentricity node can act as a reference node because it is the most eccentric of all. Concerning LCSS and EDR, the value of the threshold is set at 0 because the imprecision associated with the geographic coordinates had been removed by the map-matching operation. All the distances are implemented in python. We also use python implementation of hierarchical clustering provided by SciPy4 library.

C. Transformations parameters

We choose to delete 20% of vertices in the nodes deletion transformation. Concerning the loops addition transformation, we select 5% of nodes in each trajectory and set the neighborhood radius to 500 meters. For the detours adding transformation, we set the projection parameter to 5 (nodes) in order to simulate medium detours. The neighborhood radius parameter remains fixed at 500 meters. These parameters make it possible to limit the number of generated loops or detours on the trajectories in order to avoid strong distortions. This choice of neighborhood radius aims to produce loops of intermediate size (neither too small nor too large) and also ensure the presence of neighboring nodes in the graph representing the road network. Note that we conducted various experiments with different parameters and the results obtained follow the same tendency as those corresponding to the parameters mentioned above.

D. Clustering

The ground truth for the clustering consists of a sample of 100 trajectories split into five non-overlapping clusters of 20 trajectories each. Figure 2 illustrates the disjoint zones as well as the trajectories that are extracted there. We compute the pairwise distance matrices for each distance and apply Hierarchical Cluster Analysis. We use the ward linkage measure for HCA to minimize the variance of clusters being merged.

VII. RESULTS

A. Impact of nodes deletion

Figure 3 presents the results for nodes deletion. Among the distances compared, Hausdorff appears to be the most robust, with correlations regularly between 80% and more than 90%. STLC is the second distance that best handles the effects of nodes deletion transformation with correlation

Distances

Formulas

ST LC(T 1 , T 2 ) = v i ∈T 1 e -d G (v i ,T 2 ) n + v j ∈T 2 e -d G (v j ,T 1 ) m d G (v i , T 2 ) = min v j ∈T 2 d G (v i , v j ) DT W (T 1 , T 2 ) =              0, If n = 0 and m = 0 ∞, If n = 0 or m = 0 d G (Head(T 1 ), Head(T 2 )) + min      DT W (T 1 , Rest(T 2 )) DT W (Rest(T 1 ), T 2 ) DT W (Rest(T 1 ), Rest(T 2 ))
Otherwise

LCSS(T 1 , T 2 ) =          0, If n = 0 or m = 0 1 + LCSS(Rest(T 1 ), Rest(T 2 )), if d G (Head(T 1 ), Head(T 2 )) ≤ , max LCSS(Rest(T 1 ), T 2 ) LCSS(T 1 , Rest(T 2 )) Otherwise EDR(T 1 , T 2 ) =              n, If m = 0 m, If n = 0 min      EDR(Rest(T 1 ), Rest(T 2 )) + subcost EDR(Rest(T 1 ), T 2 ) + 1 EDR(T 1 , Rest(T 2 )) + 1 Otherwise subcost = 0, If d G (Head(T 1 ), Head(T 2 )) ≤ 1, Otherwise ERP (T 1 , T 2 ) =              n l d G (v i , v g ), If m = 0 m l d G (v j , v g ), If n = 0 min      ERP (

B. Impact of adding loops

The curves presented in Figure 4 show that in the case of loop addition, the three best distances are in order Hausdorff, STLC and EQRP. The Hausdorff distance obtains the best performance, with correlations reaching 80% while STLC correlation peak is between 70% and 80%. The curves describing the evolution of the EQRP correlations indicate a good capability to resist the addition of loops with correlation values that reach 70%. The graph versions of ERP, EDR, LCSS and DTW remain less efficient than the leading trio. Especially ERP correlations are close to 0 as for the nodes deletion transformation. However DTW stands out with better correlation values that reach almost 50%. This gives it a medium resistance to transformation by adding loops.

C. Impact of adding detours

From Figure 5 we note that adding detours has quite a similar impact on distances as loops addition. That could be expected as loops are kind of detours. ERP variants are again the least resilient to this transformation with correlations close to 0. EDR and LCSS follow respectively with correlation also close to 0, which also show a weak capacity to manage this transformation. DTW, presents an average resistance with correlations reaching at best 50%. The three best distances remain in order: Hausdorff, STLC and EQRP. Hausdorff distance and STLC present performances that reach 80%, while EQRP reach 70%. This shows their good ability to manage the addition of detours.

D. Clustering performances

The results from the clustering comparison, presented in table II, corroborate the conclusions resulting from experiences related to trajectory transformations. Hausdorff distance, STLC and EQRP manage to distribute perfectly (or almost perfectly) trajectories into distinct clusters. The same goes for DTW which stands out here, with substantially equal performance unlike what was observed for the comparison by transformations. LCSS, EDR and ERP variants produce low quality partitioning. ERP remains the worst-performing distance.

VIII. DISCUSSION

The Hausdorff distance resists well to various transformations because of its formulation based on the evaluation of the minimum distances between nodes and trajectories. Actually, these distances tend to remain stable even when the trajectories are deformed. Indeed, the minimum distance between a node u and a trajectory varies only if the trajectory noticeably changes its shape so that it modifies its node closest to u. STLC good abilities can also be explained by the node-to-trajectory distance computation scheme. The case of DTW is particular in the sense that it is a relatively sensitive distance against transformations but it has a good discrimination capability. This is explained by the fact that DTW computes the real distances between the nodes of the compared trajectories. It subsequently returns a relatively precise distance between trajectories but any modification in the position of the nodes is automatically reflected in the distance calculated. Concerning warping distances LCSS and EDR, their weak performance is due to their calculation method that does not include the actual distance between nodes and uses matching thresholds. In the case of ERP we suspect a negative impact of the reference nodes even if it remains to be proved. Our choice of a double implementation of ERP network version, aimed to highlight the impact of the reference node choice on its robustness and accuracy but the results do not lead to a clear conclusion. Although the results seem to show a slight predominance of Hausdorff distance and STLC over EQRP regarding transformations, it should be remembered that these two distances do not take into account the order of the nodes and are therefore limited in their evaluation of the distances between trajectories. This limitation effect is mainly perceptible during the comparison to similar trajectories. Take the example of three car trajectories T 1 , T 2 and T 3 such that the first two are disjoint and T 3 is identical to T 2 in terms of nodes but their sequences are reversed, i.e. T 2 is from node A to node B while T 3 is in the opposite direction. When comparing T 1 to T 2 and T 3 with Hausdorff distance and STLC, it is not possible to distinguish T 2 from T 3 which will appear perfectly identical because they are at the same distance of T 1 . However EQRP will highlight the difference due to the direction of each trajectory because it aligns trajectories according to their nodes' sequences.

IX. CONCLUSION

We propose a new distance for network-constrained trajectories combining the qualities of node-to-trajectory distances with those of the edit distances. This distance also considers trajectories both as non-ordered sets and as sequences. We compare our distance to some of the most popular distances for network-constrained trajectories by assessing their resistance to transformations and their clustering performances. The result shows that our distance outperforms all the existing edit distances and is in some cases more effective than node-totrajectory distances.

As future work, we are interested in the use of EQRP to query large network-constrained trajectories databases. We also project to create a network-constrained trajectory analysis platform integrating a search engine. This is a perspective that would lead to the development of parallel or distributed architectures for trajectories storage and distances computation and will include indexes specifically defined to the EQRP distance, based for example on the triangle inequality.

Fig. 1 :

 1 Fig. 1: Network-constrained trajectories.

Fig. 2 :

 2 Fig. 2: Distinct clusters with their trajectories.

Fig. 3 :

 3 Fig. 3: Spearman's correlations by distance as a function of the number of nearest neighbors for the nodes deletion transformation (Ratio = 20%).

Fig. 4 :

 4 Fig. 4: Spearman's correlations by distance as a function of the number of nearest neighbors for the adding loop transformation (Ratio = 5% & Radius = 500 meters).

Fig. 5 :

 5 Fig. 5: Spearman's correlations by distance as a function of the number of nearest neighbors for the adding detours transformation (Projection = 5 & Radius = 500 meters).

TABLE I :

 I Competing distances formulas

  Rest(T 1 ), Rest(T 2 )) + d G (Head(T 1 ), Head(T 2 )) ERP (Rest(T 1 ), T 2 ) + d EQRP proves its robustness by generating correlations around 70 % in the best case. Moreover, it is much more efficient than the graph versions of the other warping distances as ERP, EDR, DTW and LCSS. Only DTW produce average correlation values while the other three have correlations close to 0. The performance of ERP, regardless of the variant considered, stands out as the worst.

G (Head(T 1 ), v g ) ERP (T 1 , Rest(T 2 )) + d G (Head(T 2 ), v g )

Otherwise

values regularly above 80%.

TABLE II :

 II Clustering performances by distance

	Distances	ARI	Silhouette
	ST LC	1.0	0.887
	DT W	1.0	0.557
	EQRP	1.0	0.538
	Hausdorf f	0.975	0.488
	LCSS	0.283	0.147
	EDR	0.173	0.106
	ERP G1	0.022	0.367
	ERP G2	0.006	0.434

https://www.openstreetmap.org

https://osm2po.de

https://github.com/bmwcarit/barefoot

https://docs.scipy.org