Nadir Farhi
email: nadir.farhi@ifsttar.fr

Habib Haj-Salem
email: habib.haj-salem@ifsttar.fr

Jean-Patrick Lebacque
email: jean-patrick.lebacque@ifsttar.fr

Algebraic Approach for Performance Bound Calculus on Transportation Networks (Road Network Calculus)

Keywords: traffic flow theory, maximum travel time, network calculus, min-plus algebra

We propose in this article an adaptation of the basic techniques of the deterministic network calculus theory to the road traffic flow theory. Network calculus is a theory based on min-plus algebra. It uses algebraic techniques to compute performance bounds in communication networks, such as maximum end-to-end delays and backlogs. The objective of this article is to investigate the application of such techniques for determining performance bounds on road networks, such as maximum bounds on travel times. The main difficulty to apply the network calculus theory on road networks is the modeling of interaction of cars inside one road, or more precisely the congestion phase. We propose a traffic model for a single-lane road without passing, which is compatible with the network calculus theory. The model permits to derive a maximum bound of the travel time of cars through the road. Then, basing on that model, we explain how to extend the approach to model intersections and large-scale networks.

Notations

𝐹 the set of time functions 𝑓 is non -decreasing and 𝑓 𝑡 = 0, ∀ 𝑡 < 0 .

⊕ element-wise operation (min-plus addition in 𝐹). 𝑓 ⊕ 𝑔 𝑡 ≔ min 𝑓 𝑡 , 𝑔 𝑡 . * min-plus convolution in 𝐹.

𝑓 * 𝑔 𝑡 ≔ inf 0≤𝑠≤𝑡 𝑓 𝑠 + 𝑔 𝑡 -𝑠 .

⊘ min-plus de-convolution in 𝐹.

(𝑓 ⊘ 𝑔)(𝑡) ≔ sup 𝑠≥0 𝑓 𝑡 + 𝑠 -𝑓 𝑠 .

𝜀 the zero element of the dioid (𝐹,⊕, * ). 𝜀 𝑡 = +∞, ∀𝑡 ≥ 0.

𝑒 the unity element of the dioid (𝐹,⊕, * ). 𝑒 𝑡 = +∞, ∀𝑡 > 0, and 𝑒 0 = 0.

𝑓 𝑘 convolution power. 𝑓 𝑘 = 𝑓 * 𝑓 * … * 𝑓 (𝑘 times).

𝐵(𝑡)

the backlog at time 𝑡. 𝐵 𝑡 ≔ 𝑈 𝑡 -𝑌(𝑡).

𝑑(𝑡)

the virtual delay at time 𝑡. 𝑑 𝑡 ≔ Inf ≥ 0, 𝑌 𝑡 + ≥ 𝑈 𝑡 .

𝛾 𝑝 a particular function in 𝐹. 𝛾 𝑝 𝑡 = +∞ ∀𝑡 > 0, and 𝛾 𝑝 0 = 𝑝.

𝛿 𝑇 a particular function in 𝐹. 𝛿 𝑇 𝑡 = 0 ∀𝑡 ≤ 𝑇, and 𝛿 𝑇 𝑡 = +∞ ∀𝑡 > 𝑇.

[expr] + max 0, 𝑒𝑥𝑝𝑟 .

𝑞 car-flow.

𝑞 𝑚𝑎𝑥 the maximum car-flow.

𝑞 𝑖 (𝑡) the car outflow from the 𝑖 th section at time 𝑡.

𝑄 𝑖 (𝑡) the cumulated car outflow from the 𝑖 th section from time zero to time 𝑡.

𝑛 𝑖 the number of cars in the 𝑖 th section at time zero.

𝑛 𝑚𝑎𝑥 the maximum number of cars that a section can contain.

𝜌 the average car-density in the road.

𝜌 𝑐 the critical car-density.

𝜌 𝑗 the jam car-density.

𝑣 the free car speed.

𝑤 the backward wave speed.

∆𝑥 the section length.

𝑞(𝜌)

the car-flow function of the car-density (fundamental traffic diagram).

𝜏(𝜌)

the travel time function of the car-density.

𝜏 𝑚𝑎𝑥 (𝜌) a maximum bound of the travel time, function of the car-density.

𝑓(𝑘)

𝑚 𝑘=1 min 1≤𝑘≤𝑚 𝑓(𝑘).

Introduction

The recent advances in information and communication technologies permit to obtain valuable information on the traffic state by means of probe vehicles [START_REF] Amin | Mobile century-using GPS mobile phones as traffic sensors: a field experiment[END_REF], [START_REF] Herrera | Traffic flow reconstruction using mobile sensors and loop detector data[END_REF]. The information is then either used to derive reliable traffic indicators, or sent (after required analyses, filtrations and reformulations) to connected vehicles via intelligent transportation equipments, in order to improve the traffic conditions. One of the most important traffic indicators that drivers need to receive in order to optimize their trip, is the travel time through the possible paths to their destinations. Even though the average value of the travel time estimation is determinant for drivers, its deviation may be very important in some cases. In order to evaluate the deviation of the travel time, one can determine either the probability distribution of the travel time, or minimum and maximum bounds for it. Several methods exist in the literature to estimate travel times [START_REF] Coifman | Estimating travel times and vehicle trajectories on freeways using dual loop detectors[END_REF], (Claudel & Bayen, 2008), [START_REF] Claudel | Guaranteed bounds on highway travel times using probe and fixed data[END_REF], [START_REF] Ng | Distribution-freetraveltimereliabilityassessment with probabilityinequalities[END_REF].

We present in this article a traffic model that permits to derive a maximum bound of the travel time of cars passing through a single-lane road. The model is deterministic and uses algebraic techniques of the network calculus theory (a theory for performance bound calculus in communication and computer networks) [START_REF] Chang | Performance guarantees in communication networks[END_REF], [START_REF] Boudec | Network calculus -A theory for deterministic queuing systems for the internet[END_REF], [START_REF] Jiang | Stochastic Network Calculus[END_REF]. The objective of our work is to adapt the algebraic approach of the Network Calculus theory to transportation networks.

We propose in this article a first step of applying the deterministic network calculus to determine the maximum bound of the travel time. The main contribution of this article is a traffic model that permits the derivation of that bound with a theoretic formula. Moreover, the maximum bound of the travel time is derived as a function of the average car-density on the road, and is then compared to the formula that gives the average travel time, basing on the same modeling. In section 2 we give a short review in the network calculus theory and in the min-plus algebra [START_REF] Baccelli | Synchronization and Linearity[END_REF], in order to fix the notations and the language. The model is presented in section 3. It is inspired from the cell transmission model [START_REF] Daganzo | The cell transmission model: A dynamic representation of highway traffic consistent with the hydrodynamic theory[END_REF], [START_REF] Daganzo | The cell transmission model, Part II: Network traffic[END_REF], and written on a single-lane road. The cumulative flows are seen as time signals and the car-dynamics is written algebraically as a min-plus linear system [START_REF] Baccelli | Synchronization and Linearity[END_REF]. The impulse response of that system is then interpreted in term of guaranteed minimum service on the road. A maximum travel time of cars through the road is then derived from that guaranteed minimum service, as done basically in deterministic network calculus [START_REF] Chang | Performance guarantees in communication networks[END_REF], [START_REF] Boudec | Network calculus -A theory for deterministic queuing systems for the internet[END_REF]. Thus, a formula giving a maximum bound for the travel time as a function of the average cardensity in the road is obtained. In section 4, we propose a procedure for extending our approach to intersections and large transportation networks.

Review in network calculus

The procedure here is to consider a single-lane road as a server and apply the network calculus theory to determine a maximum bound of the travel time of cars through the road. In this short section we recall some basic results of the network calculus theory. In order to fix notations, a time function 𝑓(𝑡) (written with a small letter) expresses a given flow at time 𝑡, while 𝐹(𝑡) (written with a capital letter) denotes the cumulative flow 𝑓(𝑠)𝑑𝑠 𝑡 0

. For an arrival flow 𝑈 arriving to a server that is empty of data at time zero, a maximum arrival curve 𝛼 is associated in order to upper-bound the arrivals.

Arrival curve. 𝛼 is a (maximum) arrival curve for 𝑈 if 𝑈 𝑡 -𝑈 𝑠 ≤ 𝛼 𝑡 -𝑠 , ∀ 0 ≤ 𝑠 ≤ 𝑡.

In the other side, a minimum service curve 𝛽 is associated to the server, in order to lowerbound the service. If the output from the server is denoted 𝑌, then 𝛽 is defined as follows.

Service curve. 𝛽 is a (minimum) service curve for the server if 𝑌 𝑡 ≥ min 0≤𝑠≤𝑡 𝑈 𝑠 + 𝛽 𝑡 -𝑠 , ∀𝑡 ≥ 0.

Two indicators of the service performance are considered.  The backlog 𝐵(𝑡) of data in the server at a given time 𝑡 is defined by 𝐵 𝑡 = 𝑈 𝑡 -𝑌 𝑡 .

 The virtual delay 𝑑(𝑡) caused by the server at time 𝑡 is defined by 𝑑 𝑡 = Inf ≥ 0, 𝑌 𝑡 + ≥ 𝑈 𝑡 .

In the case where initial data 𝑛 0 is assumed in the server at time zero, the definition of the virtual delay remains correct by replacing the signal 𝑌 by 𝑌 -𝑛 0 .

It is easy to see that arrival and service curves are not unique. In order to obtain good bounds on the backlog and on the virtual delay on a given server, it is necessary to consider "good" arrival and service curves. A good arrival (resp. service) curve is simply the minimal (resp. maximal) one; see (Le [START_REF] Boudec | Network calculus -A theory for deterministic queuing systems for the internet[END_REF] and [START_REF] Chang | Performance guarantees in communication networks[END_REF] for more details.

We recall below a basic result of the deterministic network calculus (Le [START_REF] Boudec | Network calculus -A theory for deterministic queuing systems for the internet[END_REF]) [START_REF] Chang | Performance guarantees in communication networks[END_REF], that gives three bounds for a unique server. If 𝛼 is an arrival curve for an arrival flow 𝑈 to a given server that offers a minimum service curve 𝛽, then we have  The virtual delay is bounded as follows.

𝑑(𝑡) ≤ sup 𝑡≥0 Inf ≥ 0, 𝛼 𝑡 ≤ 𝛽 𝑡 + , ∀ 𝑡 ≥ 0.

 The backlog is bounded as follows. 𝐵(𝑡) ≤ sup 𝑠≥0 𝛼 𝑠 -𝛽 𝑠 , ∀𝑡 ≥ 0.

 The curve 𝑡 ↦ sup 𝑠≥0 (𝛼 𝑡 + 𝑠 -𝛽 𝑠 ) is an arrival curve for the departure flow 𝑌 from the server.

The maximum backlog and delay on a server are then given simply as the maximal vertical and horizontal distances between the arrival and the service curves; see Figure 1.

Figure 1. On the left side: Schema of the server. On the right side: The maximum delay 𝑑 and the maximum backlog 𝑏 determined graphically as the maximum horizontal and vertical distances between the arrival and the service curves, respectively.

Guaranteed service on a single-lane road (the model)

We present in this section an elementary model for the calculus of minimum guaranteed service for a single-lane road seen as a server. The objective of this modeling is to derive the maximum travel time of cars passing through the road, from the guaranteed minimum service of the road. The model is written on the cumulative car-flow variables 𝑄 (see notations below). It is inspired from the cell transmission model [START_REF] Daganzo | The cell transmission model: A dynamic representation of highway traffic consistent with the hydrodynamic theory[END_REF], [START_REF] Daganzo | The cell transmission model, Part II: Network traffic[END_REF]; see also [START_REF] Lebacque | The Godunov scheme and what it means for the first order traffic flow models[END_REF]. The min-plus linear system theory [START_REF] Baccelli | Synchronization and Linearity[END_REF] is then used to derive the minimum guaranteed service.

Let us consider a single-lane road where cars move without passing. In order to be able to fix the car-density, and to simplify the model, we consider a ring road; see Figure 2. The road is divided into 𝑚 sections of length ∆𝑥. The maximum number of cars on one section is denoted by 𝑛 𝑚𝑎𝑥 . The jam density 𝜌 𝑗 on the road is then given by

𝜌 𝑗 = 𝑛 𝑚𝑎𝑥 ∆𝑥 .
(1)

We use the following notations:


Where 𝑞, 𝜌, 𝑣, 𝑤 denote respectively the car-flow, the car-density, the free car-speed and the backward wave speed on the road; see Figure 3. The maximum flow is then given by Let us first recall the well known case [START_REF] Farhi | Modélisation minplus et commande du trafic de villes régulières[END_REF] of an autonomous single-lane ring road, with no entry 𝑈 and no exit 𝑌. In this case the dynamics of the system is written as follows.

𝑞 𝑚𝑎𝑥 = 𝜌 𝑗 1 𝑣 + 1 𝑤 . ( 3
𝑄 1 𝑡 = min 𝑄 𝑚 𝑡 -∆𝑥 𝑣 + 𝑛 𝑚 , 𝑄 2 𝑡 -∆𝑥 𝑤 + 𝑛 1 .

𝑄 𝑖 𝑡 = min 𝑄 𝑖-1 𝑡 -∆𝑥 𝑣 + 𝑛 𝑖-1 , 𝑄 𝑖+1 𝑡 -∆𝑥 𝑤 + 𝑛 𝑖 , ∀ 𝑖 ∈ 2,3, … , 𝑚 -1 . 𝑄 𝑚 𝑡 = min 𝑄 𝑚 -1 𝑡 -∆𝑥 𝑣 + 𝑛 𝑚 -1 , 𝑄 1 𝑡 -∆𝑥 𝑤 + 𝑛 𝑚 , (4)
where 𝑛 𝑘 = 𝜌 𝑗 ∆𝑥 -𝑛 𝑘 = 𝑛 𝑚𝑎𝑥 -𝑛 𝑘 .

The system ( 4) is a discrete time event system for which the dynamics is well understood. Indeed, the system (4) can be represented by an event graph (a class of Petri nets), and its dynamics can be written linearly in min-plus algebra [START_REF] Baccelli | Synchronization and Linearity[END_REF]. Since the event graph representing the system (4) is strongly connected, the asymptotic car flow on the road is then the same at every section, and is given by the minimum over the average weights of the graph circuits [START_REF] Farhi | Modélisation minplus et commande du trafic de villes régulières[END_REF]. Three circuits are distinguished in the event graph of Figure 2.

 The interior circuit with the average growth rate

𝑛 𝑖 𝑚 𝑖=1 𝑚 ∆𝑥/𝑣 = 𝑛 𝑖 𝑚 𝑖=1
𝑚 ∆𝑥 𝑣 = 𝑣𝜌.

 The exterior circuit with the average growth rate 𝑛 𝑖 𝑚 𝑖=1

𝑚 ∆𝑥 𝑤 = 𝑤 𝜌 𝑗 -𝜌 .

 The loops over each section, with the average growth rate

𝑛 𝑚𝑎𝑥 ∆𝑥/𝑣 + ∆𝑥 𝑤 = 𝜌 𝑗 ∆𝑥 ∆𝑥/𝑣 + ∆𝑥 𝑤 = 𝑞 𝑚𝑎𝑥 .
The asymptotic car-flow is then given by 𝑞 𝜌 = 𝑞 𝑖 (𝜌) = min{𝑣𝜌, -𝑤 (𝜌 -𝜌 𝑗 ) } , ∀𝑖.

and retrieves then the fundamental diagram (2) assumed in the model.

The average travel time through the road is then given by

𝜏(𝜌) = 𝑚∆𝑥 𝜌 𝑞(𝜌) = max 𝑚 ∆𝑥 𝑣 , 𝑚 ∆𝑥 𝑤 𝜌 𝜌 𝑗 -𝜌 .
In the following, we are interested on the maximum travel time through the road, rather than the average travel time. For that, we introduce an inflow and outflow of cars on the road. We assume that cars enter into the road from the entry of section 1, and exit it from the same point after passing through the road. The dynamics (4) are then rewritten as follows:

𝑄 1 𝑡 = min 𝑈 𝑡 , 𝑄 𝑚 𝑡 -∆𝑥 𝑣 + 𝑛 𝑚 , 𝑄 2 𝑡 -∆𝑥 𝑤 + 𝑛 1 .

𝑄 𝑖 𝑡 = min 𝑄 𝑖-1 𝑡 -∆𝑥 𝑣 + 𝑛 𝑖-1 , 𝑄 𝑖+1 𝑡 -∆𝑥 𝑤 + 𝑛 𝑖 , ∀ 𝑖 ∈ 2,3, … , 𝑚 -1 .

𝑄 𝑚 𝑡 = min 𝑄 𝑚 -1 𝑡 -∆𝑥 𝑣 + 𝑛 𝑚 -1 , 𝑄 1 𝑡 -∆𝑥 𝑤 + 𝑛 𝑚 .

𝑌 𝑡 = min 𝑄 𝑚 𝑡 -∆𝑥 𝑣 + 𝑛 𝑚 , 𝑄 2 𝑡 -∆𝑥 𝑤 + 𝑛 1 .

(5)

Let us notice the difference between the dynamics of 𝑄 1 which is conditioned by the inflow 𝑈, and the dynamics of 𝑌 which is not conditioned by the inflow 𝑈.

We consider the same modeling approach as above, but here we will see the variables as signals; see [START_REF] Baccelli | Synchronization and Linearity[END_REF]. We first show that the dynamic system ( 5) is linear in a certain algebraic structure, and then use this linearity to derive a guaranteed service of the road seen as a car server.

Let us consider the set of time functions 𝐹 = 𝑓 is non -decreasing and 𝑓 𝑡 = 0, ∀ 𝑡 < 0 , endowed with the following two operations.

 Addition (element-wise minimum): 𝑓 ⊕ 𝑔 𝑡 ≔ min 𝑓 𝑡 , 𝑔 𝑡 , ∀𝑡 ≥ 0.  Product (min-plus convolution): 𝑓 * 𝑔 𝑡 ≔ min 0≤𝑠≤𝑡 {𝑓 𝑠 + 𝑔(𝑡 -𝑠)} , ∀𝑡 ≥ 0.

We then obtain a dioid structure (𝐹, ⊕, * ) [START_REF] Baccelli | Synchronization and Linearity[END_REF]) (Le [START_REF] Boudec | Network calculus -A theory for deterministic queuing systems for the internet[END_REF]. The zero element of that dioid is denoted by 𝜀 and defined by 𝜀 𝑡 = +∞, ∀𝑡 ≥ 0, while the unity element is denoted by 𝑒 and defined by

𝑒 𝑡 = 0 if 𝑡 = 0. +∞ for 𝑡 ≥ 0.
In addition, we consider the following notations:  Power: 𝑓 𝑘 denotes the product (min-plus convolution) of 𝑓 with itself 𝑘 times.  Additive closure:

𝑓 * 𝑡 ≔ 𝑒 ⊕ 𝑓 ⊕ 𝑓 2 ⊕ 𝑓 3 ⊕ …  De-convolution: (𝑓 ⊘ 𝑔)(𝑡) ∶= sup 𝑠≥0 𝑓 𝑡 + 𝑠 -𝑔 𝑠 .
and the two particular signals 𝛾 𝑝 (the gain signal) and 𝛿 𝑇 (the shift signal) in 𝐹 (see Appendix A for more details on those signals).

𝛾 𝑝 𝑡 = 𝑝 if 𝑡 = 0 +∞ for 𝑡 > 0 and 𝛿 𝑇 𝑡 = 0 if 𝑡 ≤ 𝑇 +∞ otherwise
We notice that by using the notation of the min-plus convolution, the definitions of arrival and service curves can be rewritten as follows.

 𝛼 is an arrival curve for 𝑈 if 𝑈 ≤ 𝛼 * 𝑈.  𝛽 is a (minimum) service curve for the server if 𝑌 ≥ 𝛽 * 𝑈.

Now, since the cumulative flows 𝑈, 𝑄 𝑖 , 1 ≤ 𝑖 ≤ 𝑚, and 𝑌 are time functions (or signals) in 𝐹, then by using the notations of addition and product in 𝐹, the system ( 5) is linear according to those operations, and is written as follows1 .

𝑄 1 = 𝛾 𝑛 𝑚 𝛿 ∆𝑥 𝑣 𝑄 𝑚 ⊕ 𝛾 𝑛 1 𝛿 ∆𝑥/𝑤 𝑄 2 ⊕ 𝑈.

𝑄 𝑖 = 𝛾 𝑛 𝑖-1 𝛿 ∆𝑥 𝑣 𝑄 𝑖-1 ⊕ 𝛾 𝑛 𝑖 𝛿 ∆𝑥/𝑤 𝑄 𝑖+1 , ∀ 𝑖 ∈ 2,3, … , 𝑚 -1 .

𝑄 𝑚 = 𝛾 𝑛 𝑚 -1 𝛿 ∆𝑥 𝑣 𝑄 𝑚 -1 ⊕ 𝛾 𝑛 𝑚 𝛿 ∆𝑥/𝑤 𝑄 1 .

𝑌 = 𝛾 𝑛 𝑚 𝛿 ∆𝑥 𝑣 𝑄 𝑚 ⊕ 𝛾 𝑛 1 𝛿 ∆𝑥/𝑤 𝑄 2 .

(6)

Moreover, since we are not only interested in the average quantities, but on the maximum bounds, we include the initial conditions of the system (5) (the signals are null at time zero).

To include those conditions, it is sufficient to add (min-plus addition) to each signal of system (6) the signal unity 𝑒. Then, the dynamics (6) are written as follows.

𝑄 1 = 𝛾 𝑛 𝑚 𝛿 ∆𝑥 𝑣 𝑄 𝑚 ⊕ 𝛾 𝑛 1 𝛿 ∆𝑥/𝑤 𝑄 2 ⊕ 𝑈 ⊕ 𝑒.

𝑄 𝑖 = 𝛾 𝑛 𝑖-1 𝛿 ∆𝑥 𝑣 𝑄 𝑖-1 ⊕ 𝛾 𝑛 𝑖 𝛿 ∆𝑥/𝑤 𝑄 𝑖+1 ⊕ 𝑒, ∀ 𝑖 ∈ 2,3, … , 𝑚 -1 .

𝑄 𝑚 = 𝛾 𝑛 𝑚 -1 𝛿 ∆𝑥 𝑣 𝑄 𝑚 -1 ⊕ 𝛾 𝑛 𝑚 𝛿 ∆𝑥/𝑤 𝑄 1 ⊕ 𝑒. 𝑌 = 𝛾 𝑛 𝑚 𝛿 ∆𝑥 𝑣 𝑄 𝑚 ⊕ 𝛾 𝑛 1 𝛿 ∆𝑥/𝑤 𝑄 2 ⊕ 𝑒.

(7)

The system (7) is then written as follows.

𝑄 = 𝐴 * 𝑄 ⊕ 𝐵 * 𝑈 ⊕ 𝐸 𝑌 = 𝐶 * 𝑄 ⊕ 𝑒 (8)
where

𝑄 = (𝑄 1 𝑄 2 ⋯ 𝑄 𝑚 )′, A = ε 𝛾 𝑛 1 𝛿 ∆𝑥 𝑤 ε ⋯ ε 𝛾 𝑛 𝑚 𝛿 ∆𝑡 𝛾 𝑛 1 𝛿 ∆𝑡 ε 𝛾 𝑛 2 𝛿 ∆𝑥 𝑤 ε ⋯ ε ε 𝛾 𝑛 2 𝛿 ∆𝑡 ε 𝛾 𝑛 3 𝛿 ∆𝑥 𝑤 ε ⋮ ⋮ ⋱ ⋱ ⋱ ⋱ ε ε ⋯ ε 𝛾 𝑛 𝑚 -2 𝛿 ∆𝑡 ε 𝛾 𝑛 𝑚 -1 𝛿 ∆𝑥 𝑤 𝛾 𝑛 𝑚 𝛿 ∆𝑥 𝑤 ε ⋯ ε 𝛾 𝑛 𝑚 -1 𝛿 ∆𝑡 𝜀 , B = e ε ⋮ ⋮ ε ε
, 𝐶 = 𝜀 𝛾 𝑛 1 𝛿 ∆𝑥/𝑤 𝜀 ⋯ 𝜀 𝛾 𝑛 𝑚 𝛿 ∆𝑡 and 𝐸 = (𝑒 𝑒 ⋯ 𝑒)′.

The system (8) is a min-plus linear system. A basic result of the min-plus system theory [START_REF] Baccelli | Synchronization and Linearity[END_REF] gives then the impulse response of system (8):

𝑌 = 𝐶𝐴 * 𝐵𝑈 ⊕ 𝐸 ⊕ 𝑒 = 𝐶𝐴 * 𝐵𝑈 ⊕ 𝐶𝐴 * 𝐸 ⊕ 𝑒. (9)
Note that an impulse response 𝑌 = 𝐶𝐴 * 𝐵𝑈 would be obtained if system (6) is considered instead of system (7) (that is if the initial conditions of system (6) are not taken into account).

In that case, we would conclude directly that 𝐶𝐴 * 𝐵 is a minimum service curve of the road seen as a server (since we have 𝑌 ≥ 𝐶𝐴 * 𝐵 * 𝑈), and derive a maximum bound for the travel time through the road. But here, as mentioned above, it is necessary to take into account the initial conditions of the dynamical system, since we are interested in the maximum bound of the travel time through the road, rather than the average travel time. In order to be able to derive maximum bounds from the formula (9), as done from a service curve, we propose the following extension.

Minimum service couple

We consider here the case where the curve service is given with an additional affine term.

More precisely, we say that 𝛽, 𝜆 is a service couple for a server if

𝑌 ≥ 𝛽 * 𝑈 ⊕ 𝜆.
Then we can easily check (see Appendix B), that to obtain the three bounds given above, for that case, it is sufficient to replace 𝛽 with 𝛽 ⊕ 𝜆. That is to say that:

 The maximum backlog is bounded as follows.

𝐵(𝑡) ≤ sup 𝑠≥0 𝛼 𝑠 -(𝛽 ⊕ 𝜆) 𝑠 , ∀𝑡 ≥ 0.
 The maximum delay is bounded as follows.

𝑑 𝑡 ≤ sup 𝑡≥0 Inf ≥ 0, 𝛼 𝑡 ≤ (𝛽 ⊕ 𝜆) 𝑡 + , ∀ 𝑡 ≥ 0.

 An arrival curve for the departure flow is 𝛼 ⊘ (𝛽 ⊕ 𝜆). That is 𝑌 ≤ (𝛼 ⊘ (𝛽 ⊕ 𝜆)) * 𝑌.

Note that the curve (𝛽 ⊕ 𝜆) is not necessarily a minimu service curve for the server.

The two first bounds are then given by the maximum vertical and horizontal distances between the curves 𝛼 and 𝛽 ⊕ 𝜆.

Theorem 1. A minimum service couple for the single-lane road, seen as a server, is 𝛽, 𝜆 given by

𝛽 = 𝛾 -𝑚𝜌 ∆𝑥 * 𝑎 * + 𝜆 = 𝛽 ⊕ 𝛾 -𝑚𝜌 ∆𝑥 𝛾 𝑚𝜌 -𝑘𝜌 𝑗 + ∆𝑥 𝛿 (𝑚 -𝑘)∆𝑥 𝑣 ⊕ 𝛾 𝑘𝜌 𝑗 -𝑚𝜌 + ∆𝑥 𝛿 𝑘 ∆𝑥 𝑤 𝑚 -1 𝑘=1 𝑚 -1 𝑘=1 ⊕ 𝑒 +
where 𝑎 = 𝛾 𝑚𝜌 ∆𝑥 𝛿 𝑚 ∆𝑥/𝑣 ⊕ 𝛾 𝜌 𝑗 ∆𝑥 𝛿 ∆𝑥 𝑣 +∆𝑥 𝑤 ⊕ 𝛾 𝑚 (𝜌 𝑗 -𝜌)∆𝑥 𝛿 𝑚 ∆𝑥 𝑤 Proof. The system ( 8) is an affine min-plus system. We then have 𝑌 = 𝐶𝐴 * 𝐵𝑈 ⊕ 𝐶𝐴 * 𝐸 ⊕ 𝑒.

We need to compute 𝐴 * . Several methods can be used to compute 𝐴 * . By definition of 𝐴 * , one can simply compute 𝐴 2 , 𝐴 3 , …etc, then deduce 𝐴 * = 𝐴 𝑘 𝑘≥0 by simplifying all the terms. We can easily check that 𝐴 * is given as follows.

(𝐴 * ) 𝑖𝑗 = 𝑎 * if 𝑖 = 𝑗 𝑎 * 𝛾 𝑛 𝑘 𝑖-1 𝑘=𝑗 𝛿 𝑖-𝑗 ∆𝑡 ⊕ 𝛾 𝑛 𝑘 𝑗 -1 𝑘=𝑖 𝛿 𝑗 -𝑖 ∆𝑥 𝑤 if 𝑖 ≠ 𝑗, with k cyclic in 1,2, … , 𝑚 .
In order to compute 𝐶𝐴 * 𝐵, we need only the first column of 𝐴 * , since only the first entry of 𝐵 is not null. Thus, 𝐴 * 𝐵 = (𝐴 * ) •1 and is given by

𝐴 * 𝐵 = 𝑎 * * 𝑒 𝛾 𝑛 1 𝛿 ∆𝑡 ⊕ 𝛾 𝑛 𝑘 𝑚 𝑘=2 𝛿 𝑚-1 ∆𝑥 𝑤 ⋯ 𝛾 𝑛 𝑘 𝑚 -1 𝑘=1 𝛿 𝑚 -1 ∆𝑡 ⊕ 𝛾 𝑛 𝑚 𝛿 ∆𝑥 𝑤 .
Therefore, 𝐶𝐴 * 𝐵 = 𝑎 * * 𝑎.

Similarly, we can check that

𝐶𝐴 * 𝐸 = 𝑎 * * 𝑎 ⊕ 𝛾 𝑛 𝑖 𝑚 𝑖=𝑘 +1 𝛿 (𝑚 -𝑘)∆𝑥 𝑣 ⊕ 𝛾 𝑛 𝑖 𝑘 𝑖=1 𝛿 𝑘 ∆𝑥 𝑤 𝑚 -1 𝑘=1 𝑚 -1 𝑘=1
Then, it is not difficult to check that

𝑛 𝑖 𝑚 𝑖=𝑘+1 ≥ [𝑚𝜌 -𝑘𝜌 𝑗 ] + ∆𝑥. 𝑛 𝑖 𝑘 𝑖=1 ≥ [𝑘𝜌 𝑗 -𝑚𝜌] + ∆𝑥.
Now, since we are interested in the car outflow from the road that comes from the car inflow, without counting the cars being in the road at time zero (and stay there all time because the road is circular), we need to express the variable 𝑍 = [𝑌 -𝑛 𝑖 ] 𝑚 𝑖=1 + = 𝛾 -𝑚𝜌 ∆𝑥 𝑌 + in function of 𝑈. We then conclude that the couple of curves 𝛾 -𝑚𝜌 ∆𝑥 𝐶𝐴 * 𝐵, 𝛾 -𝑚𝜌 ∆𝑥 𝐶𝐴 * 𝐸 ⊕ 𝑒 is a minimum service couple for the ring road (see property (P4) in Appendix A).

Finally, we have 𝛾 -𝑚𝜌 ∆𝑥 * 𝑎 * 𝑎 * + = 𝛾 -𝑚𝜌 ∆𝑥 * 𝑎 * + , since 𝑎 * 𝑎 * 𝑡 = 𝑎 * 𝑡 , ∀𝑡 > 0 and 𝑎 * 𝑎 * 0 = 𝑎(0) < 𝑚𝜌∆𝑥 (see properties (P2) and (P5) in Appendix A).

■

In order to show the shape of the minimum service couple given in Theorem 1, let us take an academic example.

Example 1. Let 𝑚 = 6, ∆𝑥 = 1, 𝑣 = 1, 𝑤 = 1 2 , 𝜌 𝑗 = 1. Then we have 𝑛 𝑚𝑎𝑥 = 1, 𝜌 𝑐 = 1 3 and 𝑞 𝑚𝑎𝑥 = 1 3 . In addition, we take three cases, where we vary the average car-density on the road. Let us first notice that 𝑎 * can also be written as follows (see property (P3) in Appendix A).

𝑎 * = (𝛾 𝑚𝜌 ∆𝑥 𝛿 𝑚∆𝑥/𝑣 ) * * 𝛾 𝜌 𝑗 ∆𝑥 𝛿 ∆𝑥 𝑣 +∆𝑥 𝑤 * * (𝛾 𝑚 𝜌 𝑗 -𝜌 ∆𝑥 𝛿 𝑚∆𝑥 𝑤 ) *

Then we have the three cases:

 𝑛 𝑖 = 1 𝑚 𝑖=1
, that is 𝜌 = 1 6 < 𝜌 𝑐 . Then 𝑎 = 𝛾 1 𝛿 6 ⊕ 𝛾 5 𝛿 12 and 𝑎 * = (𝛾 1 𝛿 6 ) * .

Hence the minimum service couple (𝛽, 𝜆) is given by 𝛽 = 𝛾 -1 𝛾 1 𝛿 6 * + = 𝛾 -1 * 𝛾 1 𝛿 6 * . 𝜆 = 𝛽 ⊕ 𝛾 0 𝛿 5 ⊕ 𝛾 1 𝛿 6 ⊕ 𝛾 2 𝛿 8 ⊕ 𝛾 3 𝛿 10 = 𝛽.

 𝑛 𝑖 = 2 𝑚 𝑖=1
, that is 𝜌 = 1 3 = 𝜌 𝑐 . Then 𝑎 = 𝛾 1 𝛿 3 and 𝑎 * = 𝛾 1 𝛿 3 * . Hence the minimum service couple (𝛽, 𝜆) is given by 𝛽 = 𝛾 -2 𝛾 1 𝛿 3 * + . 𝜆 = 𝛽 ⊕ 𝛾 0 𝛿 8 ⊕ 𝛾 1 𝛿 10 .

 𝑛 𝑖 = 3 𝑚 𝑖=1
, that is 𝜌 = 1 2 > 𝜌 𝑐 . Then 𝑎 = 𝛾 1 𝛿 3 ⊕ 𝛾 3 𝛿 12 , 𝑎 * = 𝛾 1 𝛿 3 * 𝛾 3 𝛿 12 * . Hence the minimum service couple (𝛽, 𝜆) is given by 𝛽 = 𝛾 -3 𝛾 1 𝛿 3 * 𝛾 3 𝛿 12 * + . 𝜆 = 𝛽 ⊕ 𝛾 0 𝛿 10 = 𝛽.

The minimum service couples (given in Theorem 1) corresponding to each of the three cases, are shown in Figure 4.

■

We give below a corollary of Theorem 1, where by relaxing the minimum service couple given in Theorem 1, we obtain practical formulas.

Corollary 1. A minimum service couple (𝛽, 𝜆) for the road is given by

𝛽 𝑡 = q ρ t -τ ρ + 𝜆 𝑡 = 𝛽(𝑡) ⊕ 𝑤𝜌 𝑗 𝑡 - 2𝑚𝜌 𝜌 𝑗 ∆𝑥 𝑤 +
where 𝑞(𝜌) and 𝜏(𝜌) are the average flow and average travel time given by the fundamental diagram (given in section 2). Proof.

1. From the curve 𝛽 given by Theorem 1, we have 2. From the curve λ given in theorem 1, we have

𝛾 𝑚𝜌 -𝑘𝜌 𝑗 + ∆𝑥 𝛿 (𝑚 -𝑘)∆𝑥 𝑣 𝑚 -1 𝑘=1 ≥ max 𝑣𝜌 𝑗 𝑡 -𝑚 - 𝑚𝜌 𝜌 𝑗 ∆𝑥 𝑣 + , 𝛾 𝑚𝜌 -𝜌 𝑗 + ∆𝑥 𝛿 𝑚 -1 ∆𝑥 𝑣 then 𝛾 -𝑚𝜌 ∆𝑥 𝛾 𝑚𝜌 -𝑘𝜌 𝑗 + ∆𝑥 𝛿 (𝑚 -𝑘)∆𝑥 𝑣 𝑚 -1 𝑘=1 ≥ 𝑣𝜌 𝑗 𝑡 -𝑚 ∆𝑥 𝑣 + .
In the other hand, we have

𝛾 𝑘𝜌 𝑗 -𝑚𝜌 + ∆𝑥 𝛿 𝑘 ∆𝑥 𝑤 𝑚 -1 𝑘=1 ≥ max 𝑤𝜌 𝑗 𝑡 - 𝑚𝜌 𝜌 𝑗 ∆𝑥 𝑤 +
, 𝛾 (𝑚 -1)𝜌 𝑗 -𝑚𝜌 + ∆𝑥 𝛿 (𝑚 -1) ∆𝑥 𝑤 then 𝛾 -𝑚𝜌 ∆𝑥 𝛾 𝑘𝜌 𝑗 -𝑚𝜌 + ∆𝑥 𝛿 is more important than the term 𝛽(𝑡) in the service couple given in Corollary 1. Indeed the curve 𝛽(𝑡) gives simply the average service of the road, since 𝑞(𝜌) is the average car-flow and 𝜏 𝜌 is the average travel time on the road.

Corollary 2. If 𝛼 𝑡 = 𝜎 + 𝑟𝑡 is an arrival curve for the inflow 𝑈 to the road, then the following bounds are guaranteed.

 A maximum bound for the travel time of cars through the road is

τ max = max 𝜏 + 𝜎 𝑞 , 2𝑚𝜌 𝜌 𝑗 ∆𝑥 𝑤 + 𝜎 𝑤𝜌 𝑗 .
 A maximum bound for the number of cars waiting at the entry of the road (not confuse with the number of cars queuing at the exit of the road) is

𝑏 𝑚𝑎𝑥 = max 𝜎 + 𝑟𝜏, 𝜎 + 𝑟 2𝑚𝜌 𝜌 𝑗 ∆𝑥 𝑤 .
 An arrival curve for the departure flow from the road is 𝛼 𝑡 = 𝑏 𝑚𝑎𝑥 + 𝑟𝑡.

Proof. It is well known (Le [START_REF] Boudec | Network calculus -A theory for deterministic queuing systems for the internet[END_REF]) if a flow with an arrival curve 𝛼 𝑡 = 𝜎 + 𝑟𝑡 is served in a server with a minimum service curve 𝛽 𝑡 = 𝑅(𝑡 -𝑇) + , then the maximum delay is 𝑇 + 𝜎 𝑅 , the maximum backlog is 𝜎 + 𝑟𝑇, and the curve 𝜎 + 𝑟𝑇 + 𝑟𝑡 is an arrival curve for the departure flow. The result is then obtained by adapting these bounds to the case of couple service instead of minimum service curve (see Appendix B), and by using the couple service given in Corollary 1. ■ In the road traffic, it is probably not interesting to assume arrival curves with no null 𝜎. One may simply estimate the arrival flow rate (linear arrival curve) and the car-density on the road, at a given time, and want to determine the maximum three bounds at the considered time instant. In this case, we simply have from Corollary 2:

τ max (ρ) = max 𝜏 𝜌 , 2𝑚𝜌 𝜌 𝑗 ∆𝑥 𝑤 = max 1 𝑣 , 𝜌 𝜌 𝑗 -𝜌 1 𝑤 , 2𝜌 𝜌 𝑗 1 𝑤 𝑚∆𝑥, (10)
and 𝑏 𝑚𝑎𝑥 = 𝑟 𝜏 𝑚𝑎𝑥 , and 𝛼 𝑡 = 𝑟(𝑡 + 𝜏 𝑚𝑎𝑥 ).

The formula (10) tells that the maximum travel time through the road is greater than the average travel time only in the car-density interval [ 𝜌 𝑗 2 𝑤 𝑣 , 𝜌 𝑗 2 ] when 𝑤 < 𝑣. That is to say that

𝜌 𝑗 2 𝑤 𝑣 < 𝜌 < 𝜌 𝑗 2 ⇔ 𝜏 𝑚𝑎𝑥 𝜌 > 𝜏 𝜌 . (11)
In Figure 5, we show the average and the maximum travel times in function of the car-density in the case where 𝑤 < 𝑣. Let us notice that the choice of the traffic flow model (cell transmission model) and that of the fundamental traffic diagram (triangular one) are very important. Indeed these two choices permitted to write the car-dynamics on the road linearly in the min-plus algebra. This is necessary because the service couple is derived as an impulse response of a min-plus linear system.

Road network calculus

We give in this section the first ideas on how to use the results obtained on a single-lane road to extend the approach to intersections and whole transportation networks. First we notice that for multi-lane roads, one may just use one fundamental diagram for all lanes and apply the same approach as done for a single-lane road. We assume here that the same model and results are used for multi-lane roads. However, other models similar to the one given above can also be developed for multi-lane traffic.

Before presenting the procedure of extending our approach to intersections and networks, we need to recall and extend an important result of the deterministic network calculus on the series composition of servers (it shows in particular the power of the algebraic approach against other approaches). The result tells that a minimum service curve of the series composition of two servers guaranteeing 𝛽 1 and 𝛽 2 as minimum service curve for each of them, is simply the curve 𝛽 1 * 𝛽 2 . The result can easily be proved by using the associative property of the min-plus convolution; see [START_REF] Chang | Performance guarantees in communication networks[END_REF] (Le [START_REF] Boudec | Network calculus -A theory for deterministic queuing systems for the internet[END_REF]) for more details. Moreover, since the min-plus convolution is also commutative, then the order of the composition of the two servers is not important.

For our model, we need to have a similar result for a composition of servers offering minimum service couples (rather than minimum service curves). It is easy to see that the composition of two servers offering two service couples 𝛽 1 , 𝜆 1 and 𝛽 2 , 𝜆 2 guarantees a service couple 𝛽 1 * 𝛽 2 , 𝛽 1 * 𝜆 2 ⊕ 𝜆 1 . Indeed from 𝑌 = 𝛽 1 * 𝑍 ⊕ 𝜆 1 and 𝑍 = 𝛽 2 * 𝑈 ⊕ 𝜆 2 , we get 𝑌 = (𝛽 1 * 𝛽 2 ) * 𝑈 ⊕ (𝛽 1 * 𝜆 2 ⊕ 𝜆 1 ).

In communication and computer networks, one determines first the (residual) guaranteed service for all (input, output) couples through every switching router. Then, it is sufficient to compose (with a min-plus convolution) all the guaranteed services through all the arcs of a given path, to determine the service guaranteed through the whole path. Finally, one determines the maximum end-to-end delays on a communication network by calculating the maximum delay on each (origin-destination) path on the network, by simply using its guaranteed service curve, and considering the whole path as an elementary server. The residual guaranteed service calculus on a given input-output couple of a given router takes into account the control policy set in that router; see [START_REF] Chang | Performance guarantees in communication networks[END_REF] and (Le [START_REF] Boudec | Network calculus -A theory for deterministic queuing systems for the internet[END_REF] for more details.

We think that we can proceed similarly for transportation networks. Indeed, the main difference between data traffic in communication networks and car-traffic in transportation networks is the interaction between particles (drivers observe reaction times contrary to data packets), expressed in the fundamental diagram of the road. Since this difference is already taken into account in the one road model (presented above), then to extend the approach to complicated transportation networks, it remains the adaptation of the residual guaranteed service calculus on input-output couples of routers to apply in intersections, for transportation control policies. Then one only needs to compose elementary road services and residual guaranteed services, to obtain guaranteed services on whole paths. Maximum travel times can then be derived similarly. We notice here that one of the most difficult issues to solve, in order to extend the approach presented in this article to big networks, is the presence of cyclic dependencies of inflows arriving to one intersection, even though some elementary results exist to deal with that issue [START_REF] Chang | Performance guarantees in communication networks[END_REF].

Let us explain the calculus of maximum bounds of travel times in a tree-like network. Let us consider the transportation network of Figure 6. In Figure 6, the notation 𝐴 + indicates some point downstream of intersection 𝐴. The notation 𝐶 𝐴 -indicates some point upstream of intersection 𝐶 in the direction of intersection 𝐴. All other notations are interpreted similarly. We assume that we have the fundamental diagrams of the roads (𝐴 + , 𝐶 𝐴 -), (𝐵 + , 𝐶 𝐵 -), (𝐶 + , 𝐸 𝐶 -), (𝐷 + , 𝐸 𝐷 -) and (𝐸 + , 𝐹 𝐸 -). The approach for calculating a maximum bound for the travel time from 𝐴 to 𝐹 in Figure 6 is the following.

1. Determine the guaranteed service on each of the roads (𝐴 + , 𝐶 𝐴 -), (𝐵 + , 𝐶 𝐵 -), (𝐶 + , 𝐸 𝐶 -), (𝐷 + , 𝐸 𝐷 -) and (𝐸 + , 𝐹 𝐸 -), independent of the intersections. We simply use the fundamental diagrams of the roads, and the model presented above. 2. Calculate the residual services on the merge 𝐶:

a. The service guaranteed on 𝐶 for the aggregate inflows coming from 𝐶 𝐴 -and 𝐶 𝐵 -is assumed to be simply the guaranteed service on the road (𝐶 + , 𝐸 𝐶 -).

b. The control policy on 𝐶 is assumed to be known. c. From these two information (a and b), calculate the residual guaranteed services for the flows 𝐶 𝐴 -, 𝐶 + and (𝐶 𝐵 -, 𝐶 + ). This shall be done by adapting well known results of network calculus theory on the calculus of residual guaranteed service. (This is not yet done). 3. Calculate the residual services on the merge 𝐸. By the same method as in 2. Calculate the residual guaranteed services for the flows 𝐸 𝐶 -, 𝐸 + and (𝐸 𝐷 -, 𝐸 + ).

4. Finally, the guaranteed service for the flow 𝐴 + , 𝐹 𝐸 -, or simply (𝐴, 𝐹), is given by the series composition of the services 𝐴 + , 𝐶 𝐴 -, 𝐶 𝐴 -, 𝐶 + , 𝐶 + , 𝐸 𝐶 -, 𝐸 𝐶 -, 𝐸 + and 𝐸 + , 𝐹 𝐸 -, respectively. The guaranteed services for 𝐵, 𝐹 and 𝐷, 𝐹 are obtained similarly.

The maximum travel times from 𝐴 (resp. 𝐵 and 𝐷) to 𝐹 are then derived from the guaranteed services 𝐴, 𝐹 , 𝐵, 𝐹 and 𝐷, 𝐹 respectively, as done on a single road.

Conclusion

We presented in this article a network calculus traffic model that permits the derivation of a maximum bound for the travel time of cars passing through a single-lane road. An important advantage of this model is its algebraic formulation which is very powerful comparing to other formulations (e.g. the series composition). Even though the model is basic and elementary, its developments and extensions may be promising. Our future work shall be on the realization of the model extension process presented in section 5. We shall also demonstrate the effectiveness of our approach by performing numerical investigations on effective data sets.

 The maximum backlog: That is ∃ 𝑠 0 , 0 ≤ 𝑠 0 ≤ 𝑡 + 𝜏, 𝑌 𝑡 + 𝜏 ≥ min 𝑈 𝑡 + 𝜏 -𝑠 0 + 𝛽 𝑠 0 , 𝜆 𝑡 + 𝜏 .

Then

∃ 𝑠 0 , 0 ≤ 𝑠 0 ≤ 𝑡 + 𝜏, 𝑈(𝑡) ≥ min 𝑈 𝑡 + 𝜏 -𝑠 0 + 𝛽 𝑠 0 , 𝜆 𝑡 + 𝜏 .

  cumulated outflow from time zero to time 𝑡. 𝛼 maximum arrival curve (time function). 𝛽 minimum service curve (time function).

:

  𝑈 𝑡 : The cumulated inflow of cars to the road from time zero up to time 𝑡.  𝑌 𝑡 : The cumulated outflow of cars from the road, from time zero up to time 𝑡.  𝑞 𝑖 (𝑡) : The car outflow from the 𝑖 𝑡 section (the section between positions 𝑖 ∆𝑥 and 𝑖 + 1 ∆𝑥) at time 𝑡.  𝑄 𝑖 (𝑡) = 𝑞 𝑖 𝑠 𝑑𝑠 𝑡 0 The cumulated car outflow from the 𝑖 𝑡 section up to time 𝑡. 𝑛 𝑖 : The number of cars in the 𝑖 𝑡 section (the section between positions 𝑖 ∆𝑥 and 𝑖 + 1 ∆𝑥) at time zero.

Figure 2 .

 2 Figure 2. A single-lane ring road and the corresponding event-graph model.

)Figure 3 .

 3 Figure 3. The fundamental traffic diagram.

Figure 4 .

 4 Figure 4. Calculus of the service couple for the single-lane ring road in three cases of low, critical and high car-density.

Figure 5 .

 5 Figure 5. The average and the maximum travel times through the road in the case where 𝑤 = 𝑣/2.

Figure 6 .

 6 Figure 6. Tree-like network. The car-traffic goes from the left side to the right side.



  The maximum delay: Let 𝑡 ≥ 0. Let 𝜏 ≥ 0 such that 𝜏 < 𝑑 𝑡 ≔ inf{ ≥ 0, 𝑈(𝑡) ≤ 𝑌(𝑡 + )}. Then we have 𝑈(𝑡) > 𝑌(𝑡 + 𝜏). We have 𝑌 𝑡 + 𝜏 ≥ min min 0≤𝑠≤𝑡+𝜏 𝑈 𝑡 + 𝜏 -𝑠 + 𝛽 𝑠 , 𝜆 𝑡 + 𝜏 .

Therefore

  If 𝑈(𝑡) ≥ 𝑈 𝑡 + 𝜏 -𝑠 0 + 𝛽 𝑠 0 then 𝑠 0 > 𝜏 and thus 𝛼 𝑠 0 -𝜏 ≥ 𝑈 𝑡 -𝑈 𝑡 + 𝜏 -𝑠 0 > 𝛽(𝑠 0 ) 𝛼 𝑡 ≤ (𝛽 ⊕ 𝜆) 𝑡 + . If 𝑈(𝑡) ≥ 𝜆(𝑡 + 𝜏) then since 𝑈 0 = 0, we have 𝛼(𝑡) ≥ 𝑈(𝑡) ≥ 𝜆(𝑡 + 𝜏). 𝑠 + 𝑧 -min 𝛽, 𝜆 𝑧 ≤ (𝛼 ⊘ (𝛽 ⊕ 𝜆))(𝑡 -𝑠).

  But since ∀𝜌, min 𝜌𝑣, 𝜌 𝑗 -𝜌 𝑤 ≥ 𝑞 𝑚𝑎𝑥 , then

	Then 𝛾	𝑚 𝑖=1	𝑛 𝑖	𝛿 𝑚 ∆𝑥 𝑣	*	* 𝛾 𝑛 𝑚𝑎𝑥 𝛿	∆𝑥 𝑣	+ ∆𝑥 𝑤	*	* 𝛾	𝑚 𝑖=1	𝑛 𝑖	𝛿 𝑚 ∆𝑥 𝑤	*	≥ min 𝜌𝑣, 𝑞 𝑚𝑎𝑥 , 𝜌 𝑗 -𝜌 𝑤 𝑡.
	𝛾	𝑚 𝑖=1	𝑛 𝑖	𝛿 𝑚 ∆𝑥 𝑣		*		* 𝛾 𝑛 𝑚𝑎𝑥 𝛿	∆𝑥 𝑣	+ ∆𝑥 𝑤	*	* 𝛾	𝑚 𝑖=1	𝑛 𝑖	𝛿 𝑚 ∆𝑥 𝑤	*	≥ min 𝜌𝑣, 𝜌 𝑗 -𝜌 𝑤 𝑡.
	Hence
					𝛾	𝑚 𝑖=1		𝑛 𝑖	𝛿 𝑚 ∆𝑥 𝑣	*	* 𝛾 𝑛 𝑚𝑎𝑥 𝛿	∆𝑥 𝑣	+ ∆𝑥 𝑤	*	* 𝛾	𝑚 𝑖=1	𝑛 𝑖	𝑤 𝛿 𝑚 ∆𝑥	*	-𝜌𝑚∆𝑥	+
																≥ min 𝜌𝑣, 𝜌 𝑗 -𝜌 𝑤 𝑡 -	min 𝜌𝑣, 𝜌 𝑗 -𝜌 𝑤 𝜌𝑚∆𝑥	+
																≥ min 𝜌𝑣, 𝜌 𝑗 -𝜌 𝑤 𝑡 -max 𝑚	∆𝑥 𝑣	,	𝜌 𝜌 𝑗 -𝜌	𝑚	∆𝑥 𝑤
				𝛾	𝑚 𝑖=1	𝑛 𝑖	𝛿 𝑚 ∆𝑥 𝑣	*	≥ 𝜌𝑣 𝑡,
				𝛾 𝑛 𝑚𝑎𝑥 𝛿	∆𝑥 𝑣	+ ∆𝑥 𝑤	*	≥ 𝑞 𝑚𝑎𝑥 𝑡,
				𝛾	𝑚 𝑖=1	𝑛 𝑖	𝛿 𝑚 ∆𝑥 𝑤	*	≥ 𝜌 𝑗 -𝜌 𝑤𝑡.

+

.

Note that, as in the standard algebra, the product operation is sometimes just not symbolized (that is to say that 𝑓 * 𝑔 can simply be written 𝑓𝑔.)

Appendix A (Details on the signals 𝜸 𝒑 and 𝜹 𝑻 )

We give here some particular signals in 𝐹 as well as some properties used in section 3 (traffic model).

 The gain signal 𝛾 𝑝 :

It is then easy to obtain the following signals:

 The signal 𝛾 𝑝 * 𝛿 𝑇 :

Moreover, we can easily check (see Figure 7) that

We explain here the convolution of the signals 𝛾 𝑝 and 𝛿 𝑇 with a signal 𝑓.

 𝛾 𝑝 𝑓 𝑡 ≔ (𝛾 𝑝 * 𝑓) 𝑡 = 𝑓 𝑡 + 𝑝, ∀𝑡.

 𝛿 𝑇 𝑓 𝑡 ≔ 𝛿 𝑇 * 𝑓 𝑡 = 𝑓 𝑡 -𝑇 , ∀𝑡.

 𝛾 𝑝 𝛿 𝑇 𝑓 𝑡 ≔ 𝛾 𝑝 * 𝛿 𝑇 * 𝑓 𝑡 = 𝑓 𝑡 -𝑇 + 𝑝, ∀𝑡. We recall the following additional properties; see [START_REF] Baccelli | Synchronization and Linearity[END_REF] and/or (Le [START_REF] Boudec | Network calculus -A theory for deterministic queuing systems for the internet[END_REF].

Appendix B (minimum service couple)

We clarify here the three bounds of network calculus in the case where we have a minimum service couple 𝛽, 𝜆 instead of a minimum service curve 𝛽. That is to say that we have 𝑌 ≥ 𝛽 * 𝑈 ⊕ 𝜆 instead of 𝑌 ≥ 𝛽 * 𝑈.


Noé Ortega-Quijano

Tahereh Marvdashti

Audrey K Ellerbee Bowden

E L Ginzton

Enhanced depolarization contrast in polarization-sensitive optical coherence tomography

We demonstrate the first application of the Differential Depolarization Index (DDI) for depolarization imaging in polarization-sensitive optical coherence tomography (PS-OCT). Unlike the widely used Degree of Polarization Uniformity (DOPU), the DDI is independent of the incident polarization state and, therefore, more robust to varying system and sample parameters. Moreover, it can be applied to single-input-polarization-state PS-OCT systems, and it also overcomes several limitations of the emerging depolarization index (DI) used in multiple-input-polarization-state systems. Our results on tissue phantoms and human skin prove that DDI yields significant depolarization contrast improvements compared to DOPU, which highlights its potential for depolarization imaging in PS-OCT.

Depolarization of polarized light by biological tissue is linked to tissue heterogeneity and can reveal clinically relevant features [1][2][3][4]. For example, the depolarizing properties of the retinal pigment epithelium have been used to segment layers of the retina in optical coherence tomography images [5,6]. While polarization-sensitive optical coherence tomography (PS-OCT) can measure the reflectivity, linear retardation and optic axis orientation of tissue layers [7,8], these parameters do not directly assess the depolarization properties of a sample.

To date, the Degree of Polarization Uniformity (DOPU) [5] has been the standard metric to characterize depolarization for PS-OCT [9][10][11][12]. Additionally, a closely-related parameter has been demonstrated to detect multiple scattering [13]. Importantly, DOPU can be calculated using data from a single-input-polarization-state (SIPS) system; however, its dependence on the polarization state incident at each layer of the sample [14] leads to an inherent ambiguity in the measurement due to the dependence on illumination conditions [15]. Recently, the Depolarization Index (DI) [16] has been proposed as a more robust metric [15]; however, measurement of the DI requires multiple-input-polarization-state (MIPS) systems, which increases the measurement time and system complexity.

As an alternative to the DI, we propose that the Differential Depolarization Index (DDI) [17], which quantifies depolarization from the variance and covariance of the anisotropic parameters of the sample, is a more robust and physically meaningful method to characterize depolarization. Importantly, DDI overcomes certain inaccuracies of DI [13] and it can be used with SIPS PS-OCT systems. In this work we demonstrate that the contrast in DDI images outperforms that of DOPU in birefringent tissues, which highlights the superior potential of DDI for depolarization imaging in PS-OCT.

Consider a dual-channel polarization-sensitive photodetector that measures the complex signals

H H V V H H V V H V H V H V H V A A A A I A A A A Q A A A A U V i A A A A                                 S                  . ( 1
)
The DOPU [5], which characterizes the degree to which the polarization state of light remains uniform over a given region of tissue, is defined as

      2 2 2 DOPU norm norm norm Q U V       , (2)
with   X  denoting the mean value of X in a tissue region defined by a sliding window, and the subscript norm referring to the normalized Stokes vector that emerges after interaction with the sample. Throughout this work, we use a Gaussian sliding window of 8(x)×8(z) pixels with a standard deviation of 1.3 pixels. The Gaussian window weighting, commonly used as an alternative to rectangular sliding windows [17], is applied to the raw Stokes vector elements before computing the DOPU using Eq. ( 2). The DOPU ranges from 0 for a completely scrambled state of polarization over the evaluation window to 1 for a state that is perfectly defined and constant over the window. Related to the DOPU is its complement-the Degree of Polarization Diversity (DOPD)-which is simply given by

DOPD 1 DOPU   . ( 3
)
Unlike the DOPU, the DOPD yields a value of 1 for a completely depolarized state of polarization over the sliding window, leading to a more intuitive relationship with depolarization (i.e., the more depolarizing the sample is, the greater value the DOPD takes).

Because the measured Stokes vector returning from the sample corresponds to rotation on the Poincaré sphere of the Stokes vector incident on the sample [14,15], the value of DOPU and DOPD are intrinsically linked to the incident polarization state. In contrast with the DOPU and DOPD, the value of the DDI is directly linked to the variance and covariance of the anisotropic properties of the sample [17], which are themselves the inherent source of the depolarization [18][19][20]. These values may be grouped in the following covariance matrix:

                         , ( 4
)
where * is the complex conjugate, 2 i  denotes the variance and , i j  the covariance, with the subscripts xy, 45, and C referring to the complex propagation constant that completely characterizes anisotropy in the x-y, ±45°, and circular directions, respectively [21]. The DDI, denoted as P  , is given by the Frobenius norm of the previous matrix: P    . Importantly, this matrix is unchanged by the incident polarization state.

In the most general case, calculating the DDI in PS-OCT would require measuring the full Jones or Mueller-Jones matrix of the sample using a MIPS system, very much like in Mueller polarimetry [22]. However, SIPS PS-OCT systems exploit the fact that, in most cases, dichroism and circular anisotropy are negligible. Working from this assumption, a clear advantage of SIPS PS-OCT is that the parameters contained in the covariance matrix  can be determined directly from known equations that solve for the linear retardance η and optic axis orientation θ [7,8]. Hence, the covariance matrix takes the simplified form

       
                          , ( 5
)
where xy

 and 45

 are the linear retardance along the x-y and ±45° directions [23]:

    cos 2 sin 2 xy 45         . ( 6
)
The variances and covariances included in Eq. ( 5) are calculated for the sliding window described

P            . ( 7
)
Note that P  is 0 for a uniform sample in which the linear retardance is perfectly defined over a given sliding window and is arbitrarily large for increasingly depolarizing samples. As it is convenient to have a bounded depolarization metric for imaging purposes, we introduce the following normalization:

  DDI tanh P   . ( 8
)
Thus, similar to the DOPD, the DDI varies from 0 for a non-depolarizing sample to 1 for a completely polarization-scrambling one. Notice that the key difference between the two metrics is that the DDI describes the variance in the underlying anisotropic parameters of the sample, while the DOPD merely describes the changes in the polarization state of light that traverses the sample.

In this sense, the DDI is a more accurate metric to characterize the intrinsic depolarizing properties of the sample.

Let us illustrate this by the following a case study. Consider a linearly birefringent tissue layer, whose well-known Mueller matrix is that of a linear retarder

    2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 1 0 0 0 0 1 0 1 0 C S C C C S S S C C S C C S S C S S S C C                                           LR M , (9)
where the compact notation

  cos m m n C n    and   sin m m n S n  
 has been used. In this example, we shall assume that the linear retardance η is deterministic, whereas the axis orientation θ presents spatial stochastic fluctuations. Furthermore, we consider that the probability density function of the axis orientation is a Wrapped-Gaussian (WG) distribution, a convenient angular statistics model that wraps the Gaussian distribution around the unit circle [START_REF] Mardia | Directional Statistics[END_REF]: . This is in contrast to the DDI, for which all three lines for different input polarization states overlap. Finally, recall that every SIPS PS-OCT system is prone to measurement errors when the polarization state incident on a given layer aligns with the optic axis of that layer, which requires an analysis of the "convergence loci" of the system as pointed out in [8]. This inherent limitation would equally affect DOPU and DDI, so the precision of depolarization metrics ultimately depends on the accuracy of the SIPS PS-OCT system itself.

    2 2 2 1 ; , exp 2 2 WG k k f                           , (10)
Note that MIPS PS-OCT systems do not share this limitation and can achieve accurate measurement of anisotropic properties for any sample; importantly, the DDI can also be appliedto MIPS systems and could further characterize depolarization due to diattenuation.

Experiments were carried out using a SIPS PS-OCT system with a configuration similar to that described in previous works [7]. The system comprised a swept-source laser (HSL-2100-ST, Santec) with a center wavelength of 1310 nm and axial and lateral resolutions of 4.7 μm and 10 μm respectively. Measurements were performed at an A-line rate of 20 kHz. The power incident on the sample was 2.5 mW and the SNR of the system was 99 dB. To quantify the difference in contrast between the DOPD and DDI images we relied on the Weber contrast parameter:

strong weak w weak I I C I   , (11)
where strong I and weak I are the average depolarization metric values in a strongly-and weaklydepolarizing tissue region.

To validate the capacity of the DDI to effectively detect depolarization, we measured a tissue phantom composed of several alternating layers of clear PDMS and TiO2+PDMS [26]. We expected the phantom to only show scattering and optical extinction due to the scattering coefficient, with no birefringent effects.  , Eq. ( 6)) are depicted in Figs. 3(b) and 3(c), respectively. given in Fig. 3(e). Note that for visualization purposes, all images were subjected to an intensity threshold determined from the backscattered intensity image by the global image threshold Matlab function using Otsu's method [27], but the depolarization metrics and their corresponding contrast parameters were calculated from raw data prior to applying the threshold.

The strong appearance of a depolarizing structure highlights the significantly better contrast of DDI vs. DOPD. Indeed, the Weber contrast calculation for the boxed regions shown in Figs. 3(d tissue surface, so it may correspond to the basal cell layer located at the epidermis and dermis junction. Melanin is produced by melanocytes in this layer, and melanin granules have been reported to depolarize light [28], so we hypothesize that the depolarization contrast in this case possibly correlates with melanin concentration, which might be increased by the presence of the BCC tumor. However, this possibility would need to be confirmed by a thoughtful study. In summary, we have proposed an alternative method for depolarization imaging in PS-OCT based on the Differential Depolarization Index. Although generally valid for any PS-OCT system, the advantageous fact that DDI is closely related to the variance and covariance of the sample anisotropic properties enables it to be straightforwardly applied to SIPS PS-OCT. The suitability of this metric was validated on a controlled tissue phantom and subsequently applied to image depolarizing structures in human skin and nail fold samples.

Our results show that DDI yields better depolarization contrast compared to DOPU for the same polarization-scrambling structures. Moreover, the inherent independence of the DDI to incident polarization makes this value more convenient to measure and straightforward to interpret. Future works should consider the optimum evaluation window for DDI or investigate the potential of the DDI to aid in segmentation of PS-OCT images.

Noé Ortega-Quijano acknowledges financial support from the France-Stanford Center for Interdisciplinary Studies (Visiting Junior Scholar Fellowship). We also thank Dr. Jean Tang and Dr. Sumaira Aasi for help with collecting the tumor tissue samples.

HA

   and V A  from PS-OCT. The corresponding Stokes vector S  returning from a sample is given by

M[ 25 ]

 25 Fig. 1(c) correspond to 60º 

Fig. 2 (

 2 a) shows the results obtained for the phantom sample, with the reflectance given in dB. Fig. 2(b) includes the normalized DOPU (left) and DOPD (right), while the normalized DDI is presented in Fig. 2(c). All depolarization metrics clearly distinguish the scattering material from the background, which constitute a first confirmation that the DDI is appropriately detecting depolarization in the same way as DOPU and DOPD. The Weber contrast calculated from the mean values in the strongly-and weakly-depolarizing regions shown in Figs. 2(b) and 2(c) is 0.644 and 0.674, respectively. The small relative increase in contrast of DDI over DOPD (4.7%) highlights the similar performance of both metrics in non-birefringent scattering samples. To demonstrate the performance of the DDI in biological tissue, we next measured an ex-vivo human skin sample excised from a region showing a BCC (basal cell carcinoma). The sample was collected under an approved IRB protocol at Stanford University, whereby the subject gave verbal consent. The backscattered intensity image is shown in Fig. 3(a). The linear retardance along the x-y and ±45° directions ( xy  and 45

Fig. 3 (

 3 Fig. 3(d) shows the normalized DOPU (left) and DOPD (right), and the normalized DDI image is

) and 3

 3 (e) yields a relative contrast improvement of roughly 90% (0.491 for DODP vs. 0.933 for DDI). Fig. 3(f) shows the profile of the normalized DOPD and DDI taken along the dashed line shown in the depolarization images. While both DOPU and DDI detect a strongly polarizationscrambling region in the measured tissue, the depolarizing structure is much more sharply localized by DDI. Such a depolarizing tissue region is at a depth of roughly 150 μm below the

Fig. 4

 4 Fig.4shows the results obtained for a human nail fold. A polarization-scrambling layer in the nail plate is observable in both depolarization images, with much greater contrast apparent for DDI compared to DOPD (Figs.4(c) and 4(d)). The contrast parameters calculated for the regions highlighted in Fig.4show a relative contrast increase of 84.7% in the DDI image over DOPD (0.456 for DODP vs. 0.843 for DDI). The normalized DOPD and DDI curves taken from along the dashed line are shown in Fig.4(b). These curves confirm that both metrics are sensitive to tissue depolarization.

Supplemental Information


Haniel Barbosa
email: haniel.barbosa@inria.fr

Efficient Instantiation Techniques in SMT (Work In Progress)

In SMT solving one generally applies heuristic instantiation to handle quantified formulas. This has the side effect of producing many spurious instances and may lead to loss of performance. Therefore deriving both fewer and more meaningful instances as well as eliminating or dismissing, i.e., keeping but ignoring, those not significant for the solving are desirable features for dealing with first-order problems.

This paper presents preliminary work on two approaches: the implementation of an efficient instantiation framework with an incomplete goal-oriented search; and the introduction of dismissing criteria for heuristic instances. Our experiments show that while the former improves performance in general the latter is highly dependent on the problem structure, but its combination with the classic strategy leads to competitive results w.r.t. state-of-the-art SMT solvers in several benchmark libraries.

Introduction

SMT solvers (see [START_REF] Barrett | Satisfiability Modulo Theories[END_REF] for a general presentation of SMT) are extremely efficient at handling large ground formulas with interpreted symbols, but they still struggle to deal with quantified formulas. Quantified first-order logic is best handled with Resolution and Superposition-based theorem proving [START_REF] Bachmair | Rewrite-Based Equational Theorem Proving with Selection and Simplification[END_REF][START_REF] Nieuwenhuis | Paramodulation-Based Theorem Proving[END_REF]. Although there are first attempts to unify such techniques with SMT [START_REF] Moura | Engineering DPLL(T) + Saturation[END_REF], the main approach used in SMT is still instantiation: formulas are freed from quantifiers and refuted with the help of decision procedures for ground formulas.

The most common strategy for finding instances in SMT is the use of triggers [START_REF] Detlefs | Simplify: A Theorem Prover for Program Checking[END_REF]: some terms in a quantified formula are selected to be instantiated and successfully doing so provides a ground instantiation for the formula. These triggers are selected according to various heuristics and instantiated by performing E -matching over candidate terms retrieved from a ground model. The lack of a goal in this technique (such as, e.g., refuting the model) leads to the production of many instances not relevant for the solving. Furthermore, unlike other non-goal-oriented techniques, such as superposition, there are no straightforward redundancy criteria for the elimination of derived instances in SMT solving. Therefore useless instances are kept, potentially hindering the solver's performance.

Our attempt to tackle this issue is two-fold:

• A method for deriving fewer instances by setting the refutation of the current model as a goal, as in [START_REF] Reynolds | Finding Conflicting Instances of Quantified Formulas in SMT[END_REF]. Thus all instances produced by this strategy are relevant.

• Heuristic instantiation is kept under control. Given their speed and the difficulty of the first-order reasoning with interpreted symbols, heuristics are a necessary evil. To reduce side effects, spurious instances are dismissed. The criterion is their activity as reported by the ground solver, in a somehow hybrid approach avoiding both the two-tiered combination of SAT solvers [START_REF] Leino | A Two-tier Technique for Supporting Quantifiers in a Lazily Proof-explicating Theorem Prover[END_REF] and deletion [START_REF] De Moura | Efficient E-Matching for SMT Solvers[END_REF].

We also introduce a lifting of the classic congruence closure procedure to first-order logic and show its suitability as the basis of our instantiation techniques. Moreover, it is shown how techniques common in first-order theorem proving are being implemented in an SMT setting, such as using efficient term indexing and performing Eunification.

Formal preliminaries

Due to space constraints, we refer to the classic notions of many-sorted first-order logic with equality as the basis for the notation in this paper. Only the most relevant are mentioned.

Given a set of ground terms T and a congruence relation on

Congruence Closure with Free Variables

To better handle the quantified formulas during instantiation algorithms we have developed a Congruence Closure with Free Variables (CCFV, for short), which extends the classic congruence closure procedure [START_REF] Nelson | Fast Decision Procedures Based on Congruence Closure[END_REF][START_REF] Nieuwenhuis | Fast congruence closure and extensions[END_REF] into handling conjunctions of equality literals with free variables, performing rigid E-unification: finding solutions to a given set of equations consistent with a set of equations E, assuming that every variable denotes a single term.

L, x y U (RV) L U ∪ {x y} (i) ∈ { , } (ii) x or y is free in U , or E ∪ U |= x y L, x t U (RT) L U ∪ {x t} (i) ∈ { , } (ii) either x is free in U or E ∪ U |= x t L, f (u) f (v) U (Decompose) L, u v U L U (Yield) (i) L = ∅ or E |= L L, f (u) t U (Ematch) L, f (u) f (tn) U . . . L, f (u) f (t 1 ) U (i) ∈ { , } (ii) f (ti) are ground terms from E (iii) E |= t f (ti), for 1 ≤ i ≤ n L, u f (u ) U (Euni) L, u tn,m n , f (u ) f (t n ) U . . . L, u t 1,m 1 , f (u ) f (t 1 ) U . . . L, u t 1,1 , f (u ) f (t 1 ) U (i) ∈ { , } (ii) ti,j, f (t i ) are ground terms from E (iii) E |= ti,j f (t i ), for 1 ≤ i ≤ n, 1 ≤ j ≤ mi L U (Close) ⊥ (i) L is inconsistent modulo E or no other
rule can be applied Our procedure implements the rules 1 shown in Table 1. To simplify presentation, it is assumed, without loss of generality, that function symbols are unary. Rules are shown only for equality literals, as their extension into uninterpreted predicates is straightforward. The calculus operates on conjunctive sets of equality literals containing free variables, annotated with equality literals between these variables and ground terms or themselves. Initially the annotations are empty, being augmented as the rules are applied and the input problem is simplified, embodying its solution.

CCFV algorithm

Given a set of ground equality literals E and a set of non ground equality literals L whose free variables are X, CCFV computes sets of equality literals U 1 , . . . , U n , denoted unifiers. Each unifier associates variables from X to ground terms and allows the derivation of ground substitutions σ 1 , . . . , σ k such that E |= Lσ i , if any:

σ i =
x → t x ∈ X; U |= x t for some ground term t. If x is free in U , t is a ground term selected from its sort class.

Since not necessarily all variables in X are congruent to ground terms in a given unifier U (denoted "free in U "), more than one ground substitution may be obtained by assigning those variables to different ground terms in their sort classes.

A terminating strategy for CCFV is to apply the rules of Table 1 exhaustively over L, except that Ematch may not be applied over the literals it introduces. There must be a backtracking when a given branch results in Close, until being able to apply Yield. In those cases a unifier is provided from which substitutions solving the given E -unification problem can be extracted.

Term Indexing

Performing E -unification requires dealing with many terms, which makes the use of an efficient indexing technique for fast retrieval of candidates paramount.

The Congruence Closure procedure in veriT keeps a signature table, in which terms and predicate atoms are kept modulo their congruence classes. For instance, if a b and both f (a) and f (b) appear in the formula, only f (a) is kept in the signature table. Those are referred to as signatures and are the only relevant terms for indexing, since instantiations into terms with the same signature are logically equivalent modulo the equalities in the current context. The signature table is indexed by top symbol2 , such that each function and predicate symbol points to all their related signatures. Those are kept sorted by congruence classes, to be amenable for binary search. Bitmasks are kept to fast check whether a class contains signatures with a given top symbol, a necessary condition for retrieving candidates from that class.

A side effect of building the term index from the signature table is that all terms are considered, regardless of whether they appear or not in the current SAT solver model. To tackle this issue, an alternative index is built directly from the currently asserted literals while computing on the fly the respective signatures. Dealing directly with the model also has the advantage of allowing its minimization, since the SAT solver generally asserts more literals than necessary. Computing a prime implicant, a minimal partial model, can be done in linear time [START_REF] Déharbe | Computing Prime Implicants[END_REF]. Moreover, the CNF overhead is also cleaned: removing literals introduced by the non-equivalency preserving CNF transformation the formula undergoes, applying the same process described in [START_REF] De Moura | Efficient E-Matching for SMT Solvers[END_REF] for Relevancy.

Implementing E -unification

The main data structure for CCFV is the "unifiers": for a set of variables X, an array with each position representing a valuation for a variable x ∈ X, which consists of:

• a field for the respective variable;

• a flag to whether that variable is the representative of its congruence class;

• a field for, if the variable is a representative, the ground term it is equal to and a set of terms it is disequal to; otherwise a pointer to the variable it is equal to, the default being itself.

Each unifier represents one of the sets U mentioned above. They are handled with a UNION-FIND algorithm with path-compression. The union operation is made modulo the congruence closure on the ground terms and the current assignments to the variables in that unifier, which maintains the invariant of it being a consistent set of equality literals.

The rules in Table 1 are implemented as an adaptation of the recursive descent E-unification algorithm in [START_REF] Baader | Unification theory[END_REF], heavily depending on the term index described shown above for optimizing the search. Currently it does not have a dedicated index for performing unification, rather relying in the DAG structure of the terms. To avoid (usually expensive) re-computations, memoization is used to store the results of E -unifications attempts, which is particularly useful when looking for unifiers for, e.g., f (x) g(y) in which both "f " and "g" have large term indexes. For now these "unification jobs" are indexed by the literal's polarity and participating terms, not taking into account their structure.

Instantiation Framework

Goal-oriented instantiation

In the classic architecture of SMT solving, a SAT solver enumerates boolean satisfiable conjunctions of literals to be checked for ground satisfiability by decision procedures for a given set of theories. If these models are not refuted at the ground level they must be evaluated at the first-order level, which is not a decidable problem in general. Therefore one cannot assume to have an efficient algorithm to analyze the whole model and determine if it can be refuted. This led to the regular heuristic instantiation in SMT solving being not goal-oriented: its search is based solely on pattern matching of selected triggers [START_REF] Detlefs | Simplify: A Theorem Prover for Program Checking[END_REF], without further semantic criteria, which can be performed quickly and then revert the reasoning back to the efficient ground solver.

In [START_REF] Reynolds | Finding Conflicting Instances of Quantified Formulas in SMT[END_REF], Reynolds et al. presented an efficient incomplete goal-oriented instantiation technique that evaluates a quantified formula, independently, in search for conflicting instances: given a satisfiable conjunctive set of ground literals E, a set of quantified formulas Q and some ∀x.ψ ∈ Q it searches for a ground substitution σ such that E |= ¬ψσ. Such substitutions are denoted ground conflicting, with conflicting instances being such that ∀x.ψ → ψσ refutes E ∪ Q.

Since the existence of such substitutions is an NP-complete problem equivalent to Non-simultaneous rigid E-unification [START_REF] Tiwari | Rigid E-Unification Revisited[END_REF], the CCFV procedure is perfectly suited to solve it. Each quantified formula ∀x.ψ ∈ Q is converted into CNF and CCFV is applied for computing sequences of substitutions3 σ 0 , . . . , σ k such that, for

¬ψ = l 1 ∧ • • • ∧ l k , σ 0 = ∅; σ i-1 ⊆ σ i and E |= l i σ i
which guarantees that E |= ¬ψσ k and that the instantiation lemma ∀x.ψ → ψσ k refutes E ∪ Q. If any literal l i+i is not unifiable according to the unifications up to l i , there are no conflicting instances for ∀x.ψ.

Currently our implementation applies a breadth-first search on the conjunction of non-ground literals, computing all unifiers for a given literal l ∈ ¬ψ before considering the next one. Memory consumption is an issue due to the combinatorial explosion that might occur when merging sets of unifiers from different literals. A more general issue is simply the time required for finding the unifiers of a given literal, which can have a huge search space depending on the number of indexed terms. To minimize these problems customizable parameters set thresholds both on the number of potential combinations and of terms to be considered.

Heuristic instantiation with instances dismissal

Although goal-oriented search avoids heuristic instantiation in many cases, triggers and pattern-matching are still the backbone of our framework. A well known side effect of them is the production of many spurious instances which not only interfere with the performances of both the ground and instantiation modules but also may lead to matching loops: triggers generating instances which are used solely to produce new instances in an infinite chain. To avoid this issues, de Moura et al. [START_REF] De Moura | Efficient E-Matching for SMT Solvers[END_REF] mention how they perform clause deletion, during backtracking, of instances which were not part of a conflict. However, this proved to be an engineering challenge in veriT, since its original architecture does not easily allow deletion of terms from the ground solver.

To circumvent this problem, instead of being truly deleted instances are simply dismissed : by labeling them with instantiation levels 4 , at a given level n only instances whose level is at most n -1 are considered. This is done by using the term indexing from the SAT model and eliminating literals whose instantiation level is above the threshold. At each instantiation round, the level is defined as the current level plus one. At the beginning of the solving all clauses are assigned level 0, the starting value for the instantiation level. At the end of an instantiation round, the SAT solver is notified that at that point in the decision tree there was an instantiation, so that whenever there is a backtracking to a point before such a mark the instantiation level is decremented, at the same time that all instances which have participated in a conflict are promoted to "level 0". This ensures that those instances will not be dismissed for instantiation, which somehow emulates clause deletion. With this technique, however, the ground solver will still be burdened by the spurious instances, but they also will not need to be regenerated in future instantiation rounds.

. . .

I 1 . . . . . . I 2 ⊥ 3 I2 ⊥ 2 I1 ⊥ 1 . . . ⊥ Figure 1: Example of instance dismissal
Consider in Figure 1 an example of instance dismissal. I 1 marks an instantiation happening at level 1, in which all clauses from the original formula are considered for instantiation. Those lead to a conflict and a backtrack to a point before I 1 , which decrements the instantiation level to 0. All instances from I 1 which were part of the conflict in ⊥ 1 are promoted to level 0, the rest kept with level 1. At I 1 only terms in clauses with level 0 are indexed. Since subsequently there is no backtracking to a point before I 1 , the instantiation level is increased to 1. At I 2 all clauses of level 1 are considered, thus including those produced both in I 1 and I 1 . After a backtrack, the level is decremented to 1 and the instances participating in ⊥ 2 are promoted. This way at I 2 only the promoted instances from the previous round are considered. Then the ground solver reaches a conflict in ⊥ 3 and cannot produce any more models, concluding unsatisfiability.

Experiments

The above techniques have been implemented in the SMT solver veriT [START_REF] Bouton | veriT: An Open, Trustable and Efficient SMT-Solver[END_REF], which previously offered support for quantified formulas solely through naïve trigger instantiation, without further optimizations 5 . The evaluation was made on the "UF", "UFLIA", "UFLRA" and "UFIDL" categories of SMT-LIB [START_REF] Barrett | The SMT-LIB Standard: Version 2.0[END_REF], which have 10, 495 benchmarks annotated as unsatisfiable. They consist mostly of quantified formulas over uninterpreted functions as well as equality and linear arithmetic. The categories with bit vectors and non-linear arithmetic are currently not supported by veriT and in those in which uninterpreted functions are not predominant the techniques shown here are not quite as effective yet. Our experiments were conducted using machines with 2 CPUs Intel Xeon E5-2630 v3, 8 cores/CPU, 126GB RAM, 2x558GB HDD. The timeout is set for 30 seconds, since our goal is evaluating SMT solvers as backends of verification and ITP platforms, which require fast answers.

The different configurations of veriT are identified in this section according to which techniques they have activated:

• veriT: the solver relying solely on naïve trigger instantiation;

• veriT i: the solver with CCFV and the signature table indexed;

• veriT ig: besides the above, uses the goal-oriented search for conflicting instances;

• veriT igd: does the term indexing from the SAT solver and uses the goal-oriented search and instance dismissal.

Figure 2a shows the big impact of the handling of instantiation by CCFV: veriT i is significantly faster and solves 326 problems exclusively, while the old configuration solves only 32 exclusively. Figure 2b presents a significant improvement in terms of problems solved (474 more against 36 less) by the use of the goal-oriented instantiation, but it also shows a less clear gain of time. Besides the natural chaotic behavior of trigger instantiation, we believe this is due to the more expensive search performed: trying to falsify quantified formulas and handling E -unification, which, in the context of SMT, has a much bigger search space than simply performing E -matching for pattern-matching instantiation. Not always the "better quality" of the conflicting instances offsets the time it took to compute them, which indicates the necessity of trying to identify beforehand such cases and avoid the more expensive search when counter-producent. While veriT ig and veriT igd solve a similar number of problems in the same categories (with a small advantage to the latter), it should be noted that they have quite diverse results depending on the benchmark (a comparison is shown in Figure 4 at Appendix A). Each configuration solves ≈ 150 problems exclusively. This indicates the potential to use both the term indexes, from the signature table and from the SAT solver model with instance dismissal, during the solving.

Regarding overall performance, CVC4 solves the most problems, being the more robust SMT solver for instantiation and also applying a goal-oriented search for conflicting instances. Both configurations of veriT solve approximately the same number of problems as Z3, although mostly because of the better performance on the sledgehammer benchmarks, which have less theory symbols. There are 124 problems solved by veriT igd that neither CVC4 nor Z3 solve, while veriT ig solves 115 that neither of these two do.

Figure 3 shows how the better veriT configuration, with the goal-oriented search and instance dismissal, performs against the other solvers. There are many problems solved exclusively by each system, which indicates the benefit of combining veriT with those systems them in a portfolio when trying to quickly solve a particular problem: while CVC4 alone solves ≈ 92% of the considered benchmarks in 30s, by giving each of the four compared systems 7s is enough to solve ≈ 97% of them.

Conclusion and future work

There is still room for improvement in our instantiation framework. Particularly, a better understanding of the instance dismissal effects is still required. Further analyzing the clauses activity would lead to a more refined promotion strategy and possibly better outcomes. Nevertheless, we believe that our preliminary results are promising.

Regarding the term indexing, besides improving the data structures our main goal is performing it incrementally: by indexing the literals from the SAT model it is not necessary to thoroughly recompute the index at each instantiation round. It is sufficient to simply remove or add terms, as well as update signatures, according to how the model has changed. The same principle may be applied to the memoization of "unification jobs": an incremental term index would allow updating the resulting unifiers accordingly, significantly reducing the instantiation effort over rounds with similar indexes.

Our goal-oriented instantiation has a very limited scope: currently conflicting instances can only be found when a single quantified formula is capable of refuting the model. As it has been shown in [START_REF] Reynolds | Finding Conflicting Instances of Quantified Formulas in SMT[END_REF] and also in our own experiments this is enough to provide large improvements over trigger instantiation, but for many problems it is still insufficient. We intend to combine CCFV with the Connection Calculus [START_REF] Otten | Restricting backtracking in connection calculi[END_REF], a complete goal-oriented proof procedure for first-order logic, in an effort for having a broader approach for deriving conflicting instances. This would present a much more complex search space than the one our strategy currently handles. Therefore the trade-off between expressivity and cost has to be carefully evaluated.

Applying different strategies in a portfolio approach is highly beneficial for solving more problems, but it could be even more so if different configurations were to communicate. Attempting pseudo-concurrent architectures such as described in [START_REF] Reger | Cooperating Proof Attempts[END_REF] in veriT is certainly worth considering.

  Impact of goal-oriented instantiation

Figure

  Figure 2: Comparisons of new term indexing, CCFV and goal-oriented search

  CVC4 vs veriT igd

Figure 3 :

 3 Figure 3: Comparisons with SMT solvers

  T, a congruence C over T is a set C ⊆ {s t | s, t ∈ T} closed under entailment: for all s, t ∈ T, C |= s t iff s t ∈ C. The congruence closure of C is the least congruence on T containing C. Given a consistent set of ground equality literals E, two terms t 1 , t 2 are said congruent iff E |= t 1 t 2 , which amounts to t 1 t 2 being in the congruence closure of the equalities in E, and disequal iff E |= t 1 t 2 . The congruence class of a given term t, represented by [t], is the partition of T induced by E in which all terms are congruent to t.

Table 1 :

 1 CCFV calculus for solving rigid E -unification in equational FOL. Multiple conclusion rules represent branching in the search. x, y refer to variables, t to ground terms, u to non-ground terms and v to terms in general.

Table 2 :

 2 2: Comparisons of new term indexing, CCFV and goal-oriented search Comparison between instantiation based SMT solvers on SMT-LIB benchmarks Our new implementations were also evaluated against the SMT solvers Z3 [8] (version 4.4.2) and CVC4 [3] (version 1.5), both based on instantiation for handling quantified formulas. The results are summarized in Table2, excluding categories whose problems are trivially solved by all systems, which leaves 8, 701 for consideration.

	Logic	Class	CVC4 Z3	veriT igd veriT ig veriT i veriT
	UF	grasshopper sledgehammer 1412 410	418 1249	431 1293	437 1272	418 1134	413 1066
	UFIDL all	61	62	56	58	58	58
		boogie	841	852	722	681	660	661
		sexpr	15	26	15	7	5	5
	UFLIA	grasshopper	320	341	356	367	340	335
		sledgehammer 1892	1581	1781	1778	1620	1569
		simplify	770	831	797	803	735	690
		simplify2	2226	2337 2277	2298	2291	2177
		Total	7947	7697	7727	7701	7203	6916

The calculus still needs to be improved, with a better presentation and the proofs of its properties, which are work in progress.

Since top symbol indexing is not optimal, the next step is to implement fingerprint indexing. The current implementation keeps the indexing as modular as possible to facilitate eventually changing its structure.

Since CCFV is non-proof confluent calculus, as choices may need to be made whenever a matching or unification must be performed, backtracks are usually necessary for exploring different options.

This is done much in the spirit of[START_REF] Ge | Solving Quantified Verification Conditions Using Satisfiability Modulo Theories[END_REF], although their labeling does not take into account the interactions with the SAT solver and is aimed at preventing matching loops, not towards "deletion".

A development version is available at http://www.loria.fr/ ~hbarbosa/veriT-ccfv.tar.gz

Acknowledgments

I would like to thank my supervisors Pascal Fontaine and David Déharbe for all their help throughout the development of this work and comments on the article.

Experiments presented in this paper were carried out using the Grid'5000 testbed, supported by a scientific interest group hosted by Inria and including CNRS, RENATER and several Universities as well as other organizations (see https://www.grid5000.fr).

* This work has been partially supported by the ANR/DFG project STU 483/2-1 SMArT, project ANR-13-IS02-0001 of the Agence Nationale de la Recherche Copyright c by the paper's authors. Copying permitted for private and academic purposes.


Mathieu Cloirec

Nicolas Moës

Gilles Marckmann

P Cartraud

TWO-SCALE ANALYSIS OF CRACK USING THE EXTENDED FINITE ELEMENT METHOD

Keywords: Multi-scale approach, X-FEM, fracture mechanics, lagrange multipliers

may come

Two-scale analysis of crack using the extended finite element method

INTRODUCTION

The aim of this work is to solve problems on a large structure containing a detail with a two-scale analysis. A coarse mesh is used for the structure, while a fine mesh is superimposed in a local zone surrounding the detail. In this region, the stiffness matrices of the coarse mesh elements are modified from the solution of elementary microscopic problems. Therefore, the microscopic and macroscopic problems are treated successively, which differs from other works in this area [START_REF] Cartraud | Analyse mathmatique de la mthode Arlequin mixte: Mathematical analysis of the mixed Arlequin method[END_REF][START_REF] Guidault | Une approche micro-macro pour le suivi de fissure avec enrichissement local[END_REF] . Moreover, the microscopic problems are solved using X-FEM to avoid to mesh complex details or to be able to vary easily the detail geometry (for instance, if the crack grows). The contributions of this work are, first, to preserve at the macroscopic scale a standard finite element method (only the stiffness matrix is affected) and second, to use the X-FEM for the microscopic scale.

The two-scale discrete problem

The finite element space is denoted by V H,h (Ω):

V H,h (Ω) =    v such that v ∈ V H (Ω \ ω) on Ω \ ω v such that v ∈ V h (ω) on ω and v continuous across ∂ω   
where V h (ω) on ω, with a mesh size h, is the microscopic space and V H (Ω \ ω), with a mesh size H, the macroscopic space. The microscopic mesh is a sub-mesh of a set of macroscopic elements (figure 1). Continuous displacements are imposed on the boundary of ω. The variational formulation is to find u H,h ∈ V H,h (Ω) such that:

Ω ε(u H,h ) : E : ε(v) dΩ = Γt F • v ds ∀v ∈ V H,h (Ω) (1)
The discrete problem is considered as a two-scale problem. The solution (u H,h ), defined on is computed like a macroscopic space defined on Ω and improved by a microscopic space (u h,⊥ ) defined on ω.

u H,h = u H + u h,⊥ , u H ∈ V H (Ω), u h,⊥ ∈ V h,⊥ 0 (ω) (2)
with:

V H 0 (ω) = v ∈ V H (ω) and v = 0 on ∂ω , V h 0 (ω) = v ∈ V h (ω) and v = 0 on ∂ω and V h,⊥ 0 (ω) = v ∈ V h 0 (ω) and ω v • λ H dω = 0 ∀λ H ∈ V H 0 (ω) The problem (1) is solved as find u H ∈ V H (Ω), u h ∈ V h 0 (ω), λ H ∈ V H 0 (ω) such that:                a Ω (u H + u h , v H ) = Γ F • v H ds ∀v H ∈ V H (Ω) a Ω (u H + u h , v h ) - ω λ H • v h dω = 0 ∀v h ∈ V h 0 (ω) - ω u h • µ H dω = 0 ∀µ H ∈ V H 0 (ω) (3)
where we use the bilinear form :

a Ω (u, v) = Ω ε(u) : E : ε(v) dΩ.
The system of equation ( 3) is solved by a two-step procedure.

First, with equations (3.2) and (3.3), we eliminate u h and λ H from the system by considering u H as a parameter. The purpose is to obtain elementary solutions as a response to imposed structural loads by activating macroscopic modes.

Then, a macroscopic problem (3.1) for u H is solved.

The use of lagrange multipliers enable to impose an orthogonal condition between the two displacement fields to avoid a linear relation between each other. The microscopic scale field u h incorporates the X-FEM enrichment 3 .

Numerical experiments

A plate with a crack under tension is treated with our approach, and on an overkill mesh taking into account the crack in order to obtain a reference solution (see figure (2.a)).

The figure (2.b) illustrates the microscopic response (u h,⊥ ) of structural load on the fine mesh and the figure (2.c) the total displacement field (u H,h ) for this same zone. The last two figures compare the results, in term of strain, between our approach and the resolution of the same problem on a uniform mesh of size h used over the entire domain.

Figure 1 :

 1 Figure 1: The two-scale problem

  micro displacement (b) Micro displacement and micro mesh.

  Total displacement and micro mesh.

  Total strain for the two-scale approach.

  Strain for the overkill mesh.

Figure 2 :

 2 Figure 2: X-FEM results for the plate with a crack