Pascal Lefèvre
email: pascal.lefevre@univ-poitiers.fr

Philippe Carré
email: philippe.carre@univ-poitiers.fr

Philippe Gaborit
email: philippe.gaborit@xlim.fr

WATERMARKING AND RANK METRIC CODES

Keywords: watermarking, rank metric codes, gabidulin, luminance, structured errors

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

INTRODUCTION

Image watermarking is an important area of research. An example of motivation is the strong need to protect online multimedia contents. To insure copyright and intellectual property over massive online distribution, we need efficient protection to control the distribution, stop manipulations and duplications by pirates or unaware normal users.

To be efficient, a watermark needs to be imperceptible, needs to embed high capacity payloads and has be to robust [START_REF] Miller | A review of watermarking principles and practices[END_REF] against the most common image processings (malicious or not) while ensuring a secure transmission of the payload [START_REF] Cayre | Watermarking security: theory and practice[END_REF].

A very useful tool to enchance the robustness of a mark is the use of error-correcting codes which permit to correct errors induced by a given attack. Now, depending on the attack and on the structure of the induced errors, the type of codes used can be more or less efficient. For instance, if the error induced on the mark is random, the best results are obtained with binary codes, like BCH codes for instance (in the case of small lengths).

For other attacks, it may happen that the error comes in packet. In that case, it is better to use more structured codes over a larger alphabet (say GF (2 m )), like Reed-Solomon codes where decoding is done by packets [START_REF] Abdul | Error correcting codes for robust color wavelet watermarking[END_REF]. One does not decode error independently on each bit, but on packet of m Thanks to Délégation générale de l'Armement for funding. bits, so that an error on each bit of the packet or only one error on one bit of the packet is corrected the same way. Therefore, based on the attack, i.e. based on the error type, we can choose an adapted error correcting code. This point of view on the error type is rather well known and led to numerous industrial applications of these Hamming codes.

In this paper, we consider a new type of metric called rank metric. Error correcting codes using this metric are already often used in network coding [START_REF] Silva | On metrics for error correction in network coding[END_REF] and cryptography [START_REF] Gaborit | New results for rank-based cryptography[END_REF]. It permits to correct errors with a specific structure. If one considers a code over GF (2 m ) of length m. Each coordinate of a codeword over GF (2 m ) is encoded by m bits, and since the code has length m, any codeword can be seen as a m × m matrix. Now, it is possible to correct errors for m × m error matrix of low rank. For instance, consider an attack on the mark which flips every bits of the mark, if one consider the usual Hamming metric, it is not possible to correct this error since all bits are false. Meanwhile, in terms of rank metric, the associated error matrix has rank 1 (because it is filled with ones only) and hence, the received modified matrix can be decoded and the original message retrieved.

Our contribution: we introduce the rank metric (section 2) and Lattice QIM method (section 3). Then, we propose a watermarking process that uses both concepts able to deal with structured errors produced by a luminance modification that are not handled by conventional Hamming codes. We explain why such structure exists in this case and provide an enhanced Lattice QIM decoder for theoretical error free decoding against this attack.

RANK METRIC CODES 2.1. Definitions and properties

We refer to [START_REF] Mukhamedovich | Theory of codes with maximum rank distance[END_REF] for more details on rank metric codes. Let us consider B = (β 1 , ..., β m ) a basis of GF (q m ) over GF (q). Consider now x = (x 1 , ..., x n ) ∈ GF (q m ) n . The rank of x over GF (q) is the rank of the matrix X = (x ij ), where

x ij = m i=1 x ij β i .
It is denoted by Rank(x|GF (q)), or by Rank(x) when there is no ambiguity. Let x and y in GF (q m ) n , one defines the rank distance between x and y as d R (x, y) = Rank(x -y). As the rank of a vector is independant of the basis, it is obvious that the rank metric is less precise than the Hamming metric as two vectors with different Hamming distance could have the same rank. The linear rank code C of length n and dimension k over GF (q m ) is a subspace of GF (q m ) n . The minimum rank distance of C is defined as

d = min(d R (c 1 , c 2 ), c 1 , c 2 ∈ C|c 1 = c 2 ).

Decoding rank metric codes

Classical bounds for Hamming distance can be adapted in rank metric context and the same type of decoding results hold. If a rank code C has minimum distance d and one receives a vector y = c + e, for c ∈ C and e an error vector of rank less than (d -1)/2, it may be possible to uniquely decode y in c. Unlike classical Hamming codes, only few families of codes with easy rank metric decoding are known. Gabidulin codes is one code family among them and has parameters [n, k, n -k + 1] over GF (q n ). These codes can decode up to (n -k)/2 errors and can be seen as an analogous family in rank metric of the well-known Reed-Solomon codes family. Different approaches to decode Gabidulin codes have been proposed in the literature such as [START_REF] Ernst | A fast matrix decoding algorithm for rank-error-correcting codes[END_REF][START_REF] Loidreau | A Welch-Berlekamp Like Algorithm for Decoding Gabidulin Codes[END_REF].

Rank metric codes in practice

In practice, we use codes in an extension GF (q m ) of GF (2), and one associates a binary vector of length m to any coordinate of the codeword, so that a codeword c can be seen as a m×m binary matrix. After an attack on the watermark, codeword c is modified with error e, which also is a m × m binary matrix. To evaluate if the rank metric is better than the classical Hamming metric, we compare the embedded watermark (a codeword c) with the modified watermark (y = c+e). Suppose m = 4. Let c be a codeword and y = c + e a modified codeword such that :

c =     1 0 1 0 1 1 0 1 1 0 0 1 1 1 0 1     , y =     1 0 1 0 1 0 0 0 1 0 0 1 1 1 1 0    
Then, the error matrix is:

e =     0 0 0 0 0 1 0 1 0 0 0 0 0 0 1 1    
The error matrix e has rank 2, hence, if the code can correct up to 2 rank errors, then it is possible to decode y into c. In terms of Hamming metric, if we had started from a length 16 binary code, it would correspond to an error of weight 4. In that particular case, it is possible to find both Hamming or rank metric codes which can decode this type of errors, for reasonable dimensions k. Suppose we now have an error matrix such that:

e =     1 1 0 1 1 1 0 1 1 0 0 0 1 0 1 1    
The Hamming weight of e is 9 and the rank of e is 4. We see that it is not possible to decode with both metrics with this error matrix. In fact, rank metric is more interesting when the error has a particular structure such as:

e =     1 1 1 0 1 1 1 0 1 1 1 0 0 0 0 0    
With Hamming metric, there are still 9 errors out of 16 bits transmitted, no binary code of length 16 is able to decode properly while with the rank metric, the rank of e is only 1. We can easily decode with such error (for instance, with a Gabidulin code of parameters [START_REF] Silva | On metrics for error correction in network coding[END_REF][START_REF] Cayre | Watermarking security: theory and practice[END_REF][START_REF] Abdul | Error correcting codes for robust color wavelet watermarking[END_REF]). Of course, such a structured error does not happen necessarily often. In the following, we show that this type of error may happen in a certain case where the rank metric behaves better than the classical Hamming metric in watermarking applications. In the next section, we introduce an embedding method we used to combine with rank metric codes.

LATTICE QIM (LQIM)

Vector quantization was introduced by B. Chen and Gregory W. Wornell ( [START_REF] Chen | Quantization index modulation: A class of provably good methods for digital watermarking and information embedding[END_REF][START_REF] Moulin | Data-hiding codes[END_REF]). To embed binary information, we have two cosets of the lattice ∆Z L cosets of dimension L and a quantizer Q m :

Λ 0 = - ∆ 4 + ∆Z L , Λ 1 = ∆ 4 + ∆Z L y = Q m (x, ∆) = x ∆ ∆ + (-1) m+1 ∆ 4 , (1)
with x a host sample, y the quantized sample and m = 0, 1.

• • • • × × × × ∆ Fig. 1.
Representation of the lattice or quantization space with L = 2 and quantization step ∆. Crosses represents quantization vectors carrying bit 1 and bit 0 are carried by circles. Figure 1 illustrates an example of the lattice space in dimension L = 2. For any circle or cross (say a quantized vector y), the diamond delimited by the dotted lines will be denoted by "quantization cell". To embed information in x, one can quantize x to the nearest quantization cell center y.

For detection, we compute which coset is closer to the received vector z:

m = arg min m∈{0,1} dist(z, Λ m ), dist(z, Λ) = min y∈Λ z -y 2 (2)
As an illustration, z will be decoded into y the quantization cell center where z is located. In the next section, we explain how the experiment results allow us to see the error structure produced by a luminance modification on a LQIM watermark coupled with a rank metric code. Results were made using the Corel image database where 1000 images where randomly chosen for our tests among the 10000 images available. Embedded messages are randomly generated binary sequences of 49 bits with L = 6 and ∆ = 28 so that image quality is maintained at DWR = 35db in the spatial domain. Pixel values were randomly chosen in an image to embed the message (in total, 49×6 = 294 coefficients).

STUDY OF THE LUMINANCE ATTACK

Detection analysis after attack

The first experiment shows binary error rates (BERs) between the original payload and the decoded payload in function of the luminance constant β. The following equation illustrates how the luminance additive constant change of parameter β acts on a quantized vector y:

z = y + β × u (3)
with u = (1, ..., 1) ∈ R L and z a corrupted version of y. When β increases, z saturates. From a geometrical point of view (in 2D for figure 1), every z travel from one quantization cell to another flipping the embedding bit at every quantization cell change. In figure 2, one can notice that the red curve looks like a square waveform and we can distinguish three cases: β = 0, 0.5, 1. In the first case, there is no error at detection at regular intervals (such as β ∈ [14, 21[). Then, the third case is similar to the first case; BER = 1 also happens at regular intervals (such as β ∈ [5, 14[). This third case represents situations when every bits of the payload are flipped. The second case (BER = 0.5), detected payloads are random sequences. Unlike the previous cases, β = 0.5 only happens for specific values instead of intervals (β = 4, 21, 40, ...). This curve clearly shows the existence of a partially structured error form. Moreover, rank metric codes are well suited against this type of error.

Luminance error structure

When we apply a luminance modification on an image, every quantized vectors will suffer the same distortion of the form z -y = βu according to equation 3. In other words, if we imagine the quantization space as in figure 1, every corrupted vectors z will travel through the quantization cells. Indeed, BER = 0 means that every z are located in a cell carrying the original bit, not necessarily the one it was quantized at embedding. When BER = 1, every z are located in a cell carrying the opposite version of the embedded bit. In the next subsection, we show how a rank metric code can remove the majority of errors.

Rank metric codes application

As a second experiment, we used a rank metric code of parameters [START_REF] Ernst | A fast matrix decoding algorithm for rank-error-correcting codes[END_REF][START_REF] Abdul | Error correcting codes for robust color wavelet watermarking[END_REF][START_REF] Gaborit | New results for rank-based cryptography[END_REF] (corrects at most errors of rank 2) and measured Image Error Rates (IERs) and embedded a codeword as the watermark payload. IERs are the ratio of images where the message was not decoded by the rank metric, i.e., the error rank Rank(e) ≥ 2. In figure 2, this code is very efficient because IERs are 0 for every β except at four values. When z makes the transition from one cell to another, LQIM decoder decodes with BER = 0.5 (which means we have full rank errors). On the other hand, a similar parameters Hamming code is [47, 23, 11] and corrects at most 5 binary errors over 47 bits. In that case, the IER curve obtained with this code is identical to the red BER curve because either we have no binary error either every bits of the payload is flipped, i.e., the Hamming metric is inefficient against this attack. In practice, the decoder cannot guess β and the probability to find β such that errors are not structured for rank metric codes depends on ∆ : a small value means a higher probability to decode with errors. Moreover, the red curve in figure 2 does not look like a square waveform for some images due to the random nature of the pixel values they contain. Some BER values might be slightly under 1 or slightly above 0, this is why we choose, for a start, a rank metric code of parameters (7, 3, 5) correcting at most errors of rank 2 to correct more errors. Theoretically, a code correcting errors of rank at most 1 is enough. In the context of a luminance image processing, these codes provide almost perfect error correction. In the next subsection, we show how to improve the LQIM rank metric decoder.

Enhanced LQIM rank metric decoder

The luminance channel is parametrized by the additive constant β. Suppose a watermarked image is damaged by this channel. At decoding, we have the modified versions z = y + β. Equation 3 shows how to improve the decoder performances by adding a controlled luminance modifications. Periodically, one can notice that we cannot properly decode for β = √ 2∆/4 with k ∈ Z. This case represents transition states traveling from one quantization cell to another, i.e., vectors z are located at the boundaries of quantization cells. Let

δ i ≤ √ 2∆/4, δ i ∈ N, 1 ≤ i ≤ n.
Then, from figure 3, we deduce the following property: there is a unique i such that the corrupted image with z + δ i cannot be well decoded with LQIM rank metric decoder and for every j = i, corrupted image with z + δ j is perfectly well decoded with LQIM rank metric decoder. By modifying z with δ i , we guarantee to have the majority of z + δ j perfectly decoded. A majority vote strategy on the decoding of multiple attacked image can get rid of the spikes at the only cost of time decoding. Taking n = 3 suffices to have good results with this decoding strategy. We have d = √ 2∆/4, δ 1 = 0, δ 2 = d/3 and δ 3 = 2d/3. In the experiments, d 6 and we used δ 1 = 0, δ 2 = 2 and δ = 4 and they represent modified versions of transmitted z. Then, we extract 3 estimations of the original payload. Using the proposed property, two out of the three payloads are correct given fixed β. In figure 4, error rates are 0 for almost every β. For non-zero values, some images may have lost some information because of the image processing, so it is impossible for the proposed method to retrieve those information. As a summary, we proposed a strategy to improve the LQIM decoder combined with a rank metric code. Cases where BER = 0.5 can be avoided by taking an average estimation of multiple decoded payloads where the luminance parameter was modified on purpose. With this enhanced decoder, it is possible to use smaller quantization steps for more invisible watermark (see figure 4 with ∆ = 8). Indeed, the correction ability of rank metric codes eliminates all errors except at some particular β values. Using the proposed property, we can shift curves i.e. IER spikes and compute the majority vote of payloads. Compared to other approaches against luminance image processing (such as embedding in the low frequency coefficients in the DCT domain), we innovate with a theoretically perfect resistance with weak quantization noise at the cost of some capacity (code rate k/n). Due to space constraints, we chose to focus on the luminance image processing to fully understand how rank metric codes work. Results are very interesting for luminance because of the error structure but less interesting for others attacks (such as JPEG compression and additive white gaussian noise) where the error type is less structured.

CONCLUSION

We introduced a new kind of error correcting codes which uses the rank distance instead of the usual Hamming distance. Then, we studied the robustness of the association of LQIM method and rank metric codes against luminance image processing. We proposed an analysis of the luminance attack and showed that rank metric codes provide great correction power because of the errors structured nature. A first application of these codes with Lattice QIM method provided good results but some errors still remain (spikes). Finally, we gave a method to enhance the LQIM rank metric decoder to obtain a theoretically error free decoding. Rank metric codes can be of great interest in digital watermarking. Of course, this result is limited to luminance modifications and similar attacks. As a perspective, we plan to study more about how rank metric codes can be applied in digital watermarking.

Fig. 2 .

 2 Fig.2. Red curve : Binary error rate of LQIM method in function of β. This curve is similar when β is negative. Blue curve : Image error rate of LQIM method combined with a rank metric code in function of β. Each point of the blue curve represents the ratio of images where the error rank Rank(e) ≥ 2 (failed decoding).

Fig.

  Fig.Image error rates of LQIM method combined with a rank metric code in function of β with controlled distortions using δ = 0, 2, 4 to shift error rates curves from one to another (spike shift).

Fig. 4 .

 4 Fig. 4. Red and blue curves respectively represents BER and IER of LQIM embedding with enhanced LQIM rank metric decoder in function of β. Experiments with ∆ = 28, 8 are showed.


Augustin Lombard

Simon Tamayo

Frédéric Fontane

Frédéric Fontane Modelling

Modelling the Time-dependent VRP through Open Data

Keywords: GRASP, open data, time-dependency modeling, vehicle routing problem

published or not. The documents may come   L'archive ouverte pluridisciplinaire

I. INTRODUCTION

HIS paper presents an approach to model the static and deterministic time-dependent vehicle routing problem by using online cartography services information.

The Vehicle Routing Problem (VRP) aims at finding the optimal set of routes to be performed by a fleet of vehicles to serve a given set of customers [1].

The VRP was first introduced in 1959 by Dantzig and Ramser [START_REF] Dantzig | The Truck Dispatching Problem[END_REF] when addressing gasoline delivery optimization problems. The VRP can be solved using exact methods (Fisher et al. [START_REF] Fisher | Vehicle routing with time windows: Two optimization algorithms[END_REF]; Bard et al. [4]; Toth and Vigo [1]; Chabrier [START_REF] Chabrier | Vehicle Routing Problem with elementary shortest path based column generation[END_REF]; Qureshi et al. [START_REF] Qureshi | An exact solution approach for vehicle routing and scheduling problems with soft time windows[END_REF]; Azi et al. [START_REF] Azi | An exact algorithm for a single-vehicle routing problem with time windows and multiple routes[END_REF]) or heuristic algorithms (Badeau et al. [START_REF] Badeau | A parallel tabu search heuristic for the vehicle routing problem with time windows[END_REF]; Cordeau et Maischberger [START_REF] Cordeau | A Parallel Iterated Tabu Search Heuristic for Vehicle Routing Problems[END_REF]; Laporte et al. [START_REF] Laporte | Classical and modern heuristics for the vehicle routing problem[END_REF]; Jones et al. [START_REF] Jones | Multi-objective metaheuristics: An overview of the current state-of-the-art[END_REF]; Montemanni et al. [START_REF] Montemanni | Ant colony system for a dynamic vehicle routing problem[END_REF]; Polimeni and Vitettaa [START_REF] Polimeni | An Approach for Solving Vehicle Routing Problem with Link Cost Variability in the Time[END_REF]; Adulyasak and Jaillet [START_REF] Adulyasak | Models and Algorithms for Stochastic and Robust Vehicle Routing with Deadlines[END_REF]). Nonetheless, the exact approaches have limitations related to the computing time, which makes them less adapted to large problems.

The VRP has many important variants that take into account different delivery configurations and constraints. A quite exhaustive taxonomy of the VRP variants can be found in the work of Eksioglu et al. [START_REF] Eksioglu | The vehicle routing problem: A taxonomic review[END_REF]. One of these variants is the time-dependent VRP also known as TDVRP, which Manuscript received December 22, 2017; revised January 30, 2018. This work was supported by ADEME (French environment and energy management agency), La Poste, Marie de Paris (City of Paris), Pomona Group and Renault, within the context of the Urban Logistics Chair at Mines ParisTech.

A. Lombard. (phone: +33663013516 e-mail: augustin.lombard@minesparistech.fr).

S. Tamayo. (e-mail: simon.tamayo@mines-paristech.fr). F. Fontane. (e-mail: frederic.fontane@mines-paristech.fr).

All authors are with MINES ParisTech, PSL Research University, 60 boulevard Saint-Michel, Paris, France. allows introducing a notion of "variable cost" in the objective function that is measured in time instead of distance. The TDVRP allows integrating the notion of travel time variability due to road congestion.

The TDVRP was first introduced in 1992 by Malandraki and Daskin [START_REF] Malandraki | Time Dependent Vehicle Routing Problems: Formulations, Properties and Heuristic Algorithms[END_REF] and by Hill and Benton [START_REF] Hill | Modelling Intra-City Time-Dependent Travel Speeds for Vehicle Scheduling Problems[END_REF]. In both cases the problem was formulated as a Mixed Integer Linear Program (MILP) problem and solved through heuristic approaches. The TDVRP has been largely explored in the literature. Gendreau et al. [START_REF] Gendreau | Time-dependent routing problems: A review[END_REF] proposed a survey on the problem in which they classify the problem in 3 categories depending on the quality and the evolution of the information related to travel times, namely: 1) Static and deterministic time-dependent VRP; 2) Static and stochastic time-dependent VRP; and 3) Dynamic time-dependent VRP.

One of the limits highlighted in the existing literature is that the time-dependent function (i.e. the evolution of the speed in time) is difficult to obtain, yet it is the basis for a pertinent solution. This paper aims at proposing a model to solve the TDVRP using open-access cartography data to model the time-dependent travel times. More precisely, a static and deterministic time function is represented using a multi-layer travel times matrix, which is built by querying cartography services -such as Google Maps. Thus, with a list of addresses or coordinates, one can find a pseudooptimum itinerary, taking traffic into account.

This paper is organized as follows. Section 2 states the problem and the research questions. Section 3 presents the proposed approach, the mathematical model and the data collection procedure. Section 4 proposes an application problem in the urban area of Paris and its results. Finally, section 5 concludes with remarks and openings.

II. PROBLEM STATEMENT

A. Model the time-dependent VRP with up to date, realistic and open access data

Professional software presently used to create delivery itineraries often avoid taking traffic into account. However, today, traffic is responsible for significant difficulties, such as late deliveries, loss of time for drivers or increased fuel consumption, to name a few. Open access cartography tools such as Google Maps exist to improve a given itinerary and reduce driving times. This paper aims at modeling and optimizing the VRP based on open web services to realistically model traffic during a given period of time.

B. Research questions

This work aims at formalizing a simple and straight forward model for solving the TDVRP. The proposed

III. PROPOSED APPROACH

A. Multi-layer travel times matrix

In order to model time-dependent travel times, this paper proposes a matrix representing travel times at different moments of a planning horizon; a multi-layer matrix, with one layer for each time-step, with each layer containing the travel times between any points in the problem.

The matrix focuses only in the time needed to go from one point to another (which will depend on traffic). As a result, it is assumed that each layer is not necessarily symmetrical.

Fig. 1 shows an example of the proposed multi-layer matrix.

As shown in Fig. 1 each layer is a matrix containing the average travel times between the nodes of the problem for a time-step. For instance, 𝑡 !" ! ! in the first layer indicates the average travel time from A to B in the first hour of the day.

B. Mathematical model

The time-dependent VRP can be modelled as an oriented graph with a node set V, and arc set A. The first node corresponds to the depot and the rest to the clients. Each arc (i, j) has an associated time-dependent cost C ij (k i ), equal to the travel time between nodes i and j at a departure time k i . In any given time-step, cost respects the triangular inequality. A given route is modelled as a list R of nodes; R a is the a th node served in the route. P i is the position of node i in a given route R. T represents the length of the planning horizon; and u i is a dummy variable used to express condition (1d). x ij ∈ {0,1} indicates if in the solution, the vehicle travels from node i to node j.

In the model, drop-off and handling times are not considered. Thus, arrival time to a node is also the departure time from the same node. The vehicle must start its tour from the depot at the beginning of the planning horizon (k 0 = 0).

With the previous definitions, the problem can be modelled as follows:

𝑚𝑖𝑛 𝑥 !" 𝐶 !" 𝑘 ! !∈! !!! !∈! (1a) 𝑈𝑛𝑑𝑒𝑟: ∀𝑖 ∈ 𝑉, 𝑥 !" = 1 !∈! (1b) ∀𝑗 ∈ 𝑉, 𝑥 !" = 1 !∈! (1c) 𝑢 ! -𝑢 ! + 𝑛 𝑥 !" ≤ 𝑛 -1 2 ≤ 𝑖, 𝑗 ≤ 𝑛, 𝑖 ≠ 𝑗 (1d) ∀𝑖, 𝑗 ∈ 𝑉, 𝑥 !" ∈ 0,1 (1e)
∀𝑖 ∈ 𝑉, 𝑢 ! ∈ ℤ (1f) ∀𝑖 ∈ 𝑉, 𝑘 ! ∈ 0, 𝑇 (1g)
Constraints (1b) and (1c) ensure flow conservation: there is exactly one arrival and one departure from every node. Constraint (1d) ensures that every route passes through node 1. Combined with the flow conservation, it proves that all nodes are served in the same route (i.e. there is no subroute). Constraint (1d) comes from the integer formulation of the travelling salesman problem of Miller et al. [START_REF] Miller | Integer Programming Formulation of Traveling Salesman Problems[END_REF]. It can be proven that P i are values for the dummy variables u i that verify constraint (1d).

In a route R, node i is at position P i . The nodes visited before node i are nodes R a with a = 1, …, P i -1. Then, departure time from i, k i , is the sum of the time needed to go from node R a to R a+1 , with a = 1, …, P i -1. Departure time k i from a node i can be computed recursively as follows:

𝑘 ! = 0 (2a) 𝑘 ! = 𝑥 ! ! ! !!! 𝐶 ! ! ! !!! (𝑘 ! ! ) ! ! !! !!! (2b)

C. Open data fetching approach

To create the multi-layer time matrix, we propose using open data. This paper uses Google Maps, and its Distance Matrix API. In order to obtain a time-dependent matrix that represents the average behavior of the modelled zone (i.e. to limit the influence of rare incidents, such as car crashes, etc.), the query must be done in the future. In the present example the query is done in the near future, two weeks after present day.

It is important to note that the time matrix is not symmetrical because of one-way roads, and different levels of congestion.

In terms of distance, the triangular inequality is not necessarily verified because of traffic and the different types of roads (e.g. highway, residential areas). The fastest road is not necessary the shortest. However, the cost in the objective function, the driving time, verifies the triangular inequality. The driving time between two points will always be the smallest.

The matrix is built as follows:

• Select N points on a map representing the clients and the depot.

0 ≤ t < 1 A B C A B C … … … … … … … … … 1 ≤ t < 2 23 ≤ t < 24 … A B C A B C 0 0 0 t AB S 2 t AC S 2 t BA S 2 t BC S 2 t CA S 2 t CB S 2 A B C A B C 0 0 0 t AB S 1 t AC S 1 t BA S 1 t BC S 1 t CA S 1 t CB S 1
• Retrieve longitude and latitude coordinates for the previously pinned points. • Discretize the time horizon in a given number of time steps. • Query online cartography services and obtain the travel times between all points, for all time steps. • Assemble the different layers into a single matrix.

D. Optimization procedure

This paper uses a greedy randomized adaptive search procedure (GRASP) and an enhancement heuristic to optimize the model. This method builds on the work of Reyes et al. [START_REF] Reyes | Vehicle routing with roaming delivery locations[END_REF], that solved a version of the VRP, the VRP with roaming delivery locations (VRPRDL).

In the first step of the optimization, a GRASP is conducted. N grasp different routes will be created to serve all of the nodes and the best will be kept. Routes are created as follows. For each node not yet inserted on the road, cost of its insertion at any place is calculated. The insertions are sorted in ascending order and the K grasp smallest are kept. One insertion is randomly chosen and the corresponding node is inserted at the corresponding location of the route. When all the nodes are placed, a total cost is calculated. At the end of the GRASP, among the N grasp routes, the one with the lowest cost is kept. The output of the GRASP is a single solution (the best one amongst the N grasp trials).

After the GRASP, an "insertion-deletion" heuristic is used to improve the current solution. L nodes are deleted and then reinserted, following a first-deleted, first-reinserted rule. First, we build a list containing the nodes that would lead to the greatest cost saving in the solution if they were to be deleted. The size of the list is stated by a given parameter K del , then the deleted nodes are randomly chosen one at a time and the list is updated after each deletion. Nodes are reinserted in the order they were deleted, at a place chosen randomly in the list of K ins places leading to the lowest cost. This heuristic is repeated N' times.

If parameters K grasp , K del and K ins are large numbers, the stochastic search of the algorithm increases, which could help avoid local optima. Indeed, higher stochasticity would need more search iterations to ensure convergence. However, if the algorithm is intended to run during a given computation time, high values for these parameters lead to less robust solutions because of the uncertainty of the insertions.

IV. APPLICATION

The approach is tested on the urban area of Paris. In the recent years, through-roads along the river Seine were closed for traffic and opened for pedestrians. This resulted in an increase of traffic on the ring road and other roads to cross the city [START_REF] Awada | Fermeture des voies sur berges rive droite à Paris[END_REF]. Managing deliveries with traffic is now a key issue for delivery services.

The example is composed of 30 clients, 1 depot, all in Paris and its urban area. The multi-layer time matrix is queried for the December 12 th 2017, from 8 a.m. to 8 p.m. with time steps of 2 hours. As a result, the multi-layer matrix is composed of 16 31x31 matrices. Fig. 2 shows clients and depot on a map of Paris urban area and table I shows the latitude and longitude of all the chosen clients. The matrix was generated by querying Google Maps API [START_REF] Google | Google Maps Distance Matrix API[END_REF] using the parameters mode = "driving" and traffic_model = "best_guess".

A. Parameters

In the first phase of the optimization, the GRASP, N = 30 routes were created with each the 30 nodes. The best route was kept. The insertion was chosen amongst the K = 3 best insertions. In the second phase of the optimization, the heuristic is run N' = 20 times, removing L = 6 nodes each time. These parameters were chosen as a trade off between efficiency and computation time. It is important to note that improvement in the choice of these parameters will lead to a better solution.

B. Results

To assess the quality of the optimization process, the best route was plotted on Google Maps to show the real itinerary taken by the vehicle to go from a node to another. Time-dependent VRP is also compared to the classical VRP, where the travel times matrix is the average of the multi-layer matrix along the time steps. To do so, total cost was compared in both cases.

On Fig. 3, the itinerary of the delivery vehicle can be seen. The vehicle starts with customer 1, does his loop, and finishes with customer 30. The itinerary seems rational, even if there are a few crossings. Crossings in the best itineraries can be explained by sub-optimality of the solution or by traffic and road types (e.g. highway versus city streets). It can be faster to exit a highway after the delivery location is passed by, which will imply crossing roads. The general pattern is the same between the two scenarios, with an 8-shape versus an "almost" 8-shape, and slight modifications in the north of the city. The modifications are explained by the traffic. The order the clients are served can change to minimize the time the vehicle is used.

In Fig. 3, the classical VRP, the pseudo-optimum route lasts 7 hours and 28 minutes. In Fig. 4, the pseudo-optimum route lasts 7 hours and 35 minutes. The relative gap is 1.7%, which seems rather small since the driving time can double during the day on some arcs. In the model, service time is neglected. When a vehicle arrives at a node, it leaves immediately. Since the vehicle starts at 8 a.m. from the depot, it has to go through morning congestion. However, in the solution, the vehicle is back before the afternoon congestion.

Fig. 3 and Fig. 4 indicate some of the travel times between nodes at the beginning and at the end of the route. For the first three nodes, the classical VRP is fastest by 8 minutes, whereas at the end of the itinerary it is slower by 15 minutes. This explains why the TDVRP can lead to faster solutions compared to the classical VRP. That is to say, that although the chosen roads are congested, it is faster to use them in the morning because the return roads are even more congested at the same time. Waiting for the return roads to be less congested leads to time savings.

The solution returned by the TDVRP is not necessarily better in terms of cost, because the cost in the classical VRP is computed with the average travel times matrix. Nonetheless it must be highlighted that the classical VRP will indeed encounter surprises related the congestion, that the TDVRP will be certainly avoiding.

We computed the cost (i.e. travel time) of the itinerary presented in Fig. 4 taken backwards. To go from the depot to node 28 (passing through node 30 and 29), takes 74 minutes at 6 p.m. whereas it takes 109 minutes at 8 a.m.

C. Sensibility of the model

In order to verify the robustness of the algorithm, a set of 20 instances of optimization was run.

Table II presents the result of these 20 simulations with the set of customers and depot proposed above. C ML indicates the cost of the instance using the multi-layer matrix, C 2D indicates the cost of the instance where the travel times matrix is the average along the time steps of the multi-layer matrix. Instances differ by the seed given to the pseudo-random number generator.

Results range from 6h59 to 7h53. This is a difference of just less than an hour. This variability comes from the parameters K grasp , K del and K ins of the optimization method, which allows a non-optimal insertion or deletion of a node in a route, to get out of local minima. There is a great sensibility to these parameters. If K grasp , K del and K ins were fixed to 1, then the optimization would be fully deterministic, and the results would all be the same. The variability would be reduced by performing more steps of the GRASP or "insertion-deletion" method. If a very large number of solutions was generated in the GRASP method, the solution would be expected to be the same in each column of the table. To be able to have 20 instances tested, each with the multi-layer and the 2D matrixes, in a reasonable time, 30 GRASP and 20 insertions-deletions are done. This is not a search of the global-optimum.

Sometimes, the multi-layer matrix leads to better solutions than the 2D matrix. On average, the multi-layer matrix is better by 4.3%. Taking traffic variability into account helps building better solutions.

D. Technical considerations in the proposed approach

The numerical method and optimizer were programmed using Python. The GRASP was based on the pseudo-code presented by Reyes et al. [START_REF] Reyes | Vehicle routing with roaming delivery locations[END_REF]. No optimizer module was used.

Retrieving data Google Maps limits the size of the test problem, since there are quotas for the number of queries per day. For free, one can make 2,500 queries per day per project. The limit is 100,000 per project when paying. If traffic information is required, no more than 100 queries at a time are allowed. A multi-layer matrix of 24 time-steps and N nodes will require 24*N 2 queries. The biggest problem that can be queried in a day has 64 nodes. Because of the 100 queries at a time limit, the method has to be well structured.

Concerning calculation time, the machine has a 2.5 GHz Intel Core i5. One core is used. The 30 GRASP iterations take 80 seconds on average, and the 20 insertions-deletions 10 seconds.

V. CONCLUSION AND PERSPECTIVES

This paper proposed a new approach to model and optimize the time-dependent vehicle routing problem (TDVRP). A discretization of time was used to create a multi-layer matrix, in which each layer contains the travel times matrix for the nodes of the problem at a given time of the day. One of the originalities of this work yields in the proposed mathematical model, in which the costs of arcs depend on the time at which the vehicle passes through them. As a result, when minimizing overall cost, optimization algorithms will favor solutions in which the "most expensive arcs" are taken in their "less expensive" times. In other words, if a given delivery point is located in a road that is heavily congested between 8 a.m. and 9 a.m., the optimization logic will avoid it until 9 a.m., creating solutions that perform other deliveries before going to the aforementioned point, positively avoiding the congested road at the most congested time.

The proposed model presents a concrete example that can be applied to any city, of planning deliveries taking traffic into account. Yet there are several improvement perspectives, especially to integrate more realistic traffic considerations. Possible improvements include parallelizing the code to compute more difficult problems more efficiently or to add complexity to the problem, such as limited capacity, stochasticity, time windows or roaming delivery locations.

Open data cartography services, such as Google maps, offer accessible data to build realistic models. The proposed approach is based in a discretization of time. It is important to consider that this choice leads to faster computation times but it also implies a rough approach from the trafficmodeling standpoint. An interesting perspective would be to create a multi-layer matrix with a denser discretization, which could be used to interpolate the time horizon and thus, obtain a continuous time-dependency function.

Fig. 2 .

 2 Fig. 2. Map of Paris urban area with 30 clients and one depot (D).

Fig. 3

 3 Fig. 3 Optimization results with the classical VRP.

Fig. 4

 4 Fig. 4 Optimization results with the TDVRP.

TABLE I

 I

		COORDINATES OF THE CLIENTS
	Client index	Latitude (°)	Longitude (°)
	0 (Depot)	48.845761	2.339546
	1	48.847397	2.348344
	2	48.758064	2.169373
	3	48.841295	2.588015
	4	48.841575	2.347142
	5	48.858862	2.294403
	6	48.886944	2.343126
	7	48.745642	2.369067
	8	48.854819	2.306285
	9	48.870349	2.33312
	10	48.892779	2.243434
	11	48.890971	2.252033
	12	48.893764	2.178884
	13	48.834553	2.237409
	14	48.788931	2.327215
	15	48.797628	2.309662
	16	48.809578	2.375374
	17	48.805202	2.388838
	18	48.827828	2.350659
	19	48.82934	2.34343
	20	48.834854	2.367539
	21	48.841673	2.351875
	22	48.852078	2.329325
	23	48.867007	2.349496
	24	48.876527	2.316532
	25	48.879352	2.348189
	26	48.909403	2.339389
	27	48.875115	2.393429
	28	48.919156	2.388138
	29	48.917793	2.243537
	30	48.718288	2.210583