Wenjuan Gu

Diego Cattaruzza
email: diego.cattaruzza@centralelille.fr

Maxime Ogier
email: maxime.ogier@centralelille.fr

Frédéric Semet
email: frederic.semet@centralelille.fr

Adaptive large neighborhood search for the commodity constrained split delivery VRP

Keywords: vehicle routing problem, multiple commodities, adaptive large neighborhood search, local search

come   L'archive ouverte pluridisciplinaire

Introduction

The vehicle routing problem (VRP) and its variants have been widely studied in the literature [START_REF] Toth | Vehicle Routing: Problems, Methods, and Applications[END_REF]. In the VRP with multiple commodities, capacitated vehicles are used to deliver a set of commodities and meet the demands of customers. Four different strategies to deliver a set of commodities to customers were presented by [START_REF] Archetti | Multicommodity vs. single-commodity routing[END_REF]: separate routing, mixed routing, split delivery mixed routing and commodity constrained split delivery mixed routing.

In the separate routing strategy, a specific set of vehicles is dedicated to each commodity, and any commodity has to be delivered to any customer by a single vehicle visit. In this case, considering each commodity separately is a classical capacitated VRP (CVRP). A customer is visited as many times as the number of commodities that he/she requires. In other words, if a customer needs multiple commodities, this customer needs to be served several times, even if all the commodities could be delivered by only one vehicle.

In the mixed routing strategy any set of commodities can be mixed in the same vehicle, and all customers must be delivered at once. If a customer requires one or more commodities, all of them are delivered by a single vehicle in a single visit. Here again, the problem that arises corresponds to a single CVRP. In this delivery strategy, when the remaining capacity of a vehicle is not sufficient to deliver all the demand of a given customer, it is wasted with a possible increase in transportation costs.

In the split delivery mixed routing strategy, any set of commodities can be mixed in the same vehicle. The commodities can be split in any possible way. Moreover, a commodity can be delivered to a customer by several vehicles. The problem that arises corresponds to the split delivery VRP (SDVRP) that was introduced by [START_REF] Dror | Savings by split delivery routing[END_REF] and [START_REF] Dror | Split delivery routing[END_REF]. For an extensive review of SDVRP, the interested reader is referred to the survey by [START_REF] Archetti | Vehicle routing problems with split deliveries[END_REF]. In the split delivery mixed routing strategy, a customer can be visited several times if this is beneficial, even if this customer requires only one commodity. This strategy minimizes transportation costs but may cause inconvenience to customers [START_REF] Archetti | To split or not to split: That is the question[END_REF]. For example, if one commodity is delivered by several vehicles, the handling time for the customer may increase significantly.

The commodity constrained split delivery mixed routing (C-SDVRP) is a strategy recently proposed by [START_REF] Archetti | Multicommodity vs. single-commodity routing[END_REF] to deliver multiple commodities. In the C-SDVRP, any set of commodities can be mixed in the same vehicle. The demand of a customer can be split, but only as many times as the required commodities and when a commodity is delivered to a customer, the entire required amount is handed over. As a consequence, one commodity can be delivered by only one vehicle.

The C-SDVRP shares similarities with the SDVRP. However, a significant difference is that the demand of a customer cannot be arbitrarily split since a commodity has to be delivered by the same vehicle. Thus, only splitting a request of a customer according to different commodities is allowed. Considering the convenience of customers and transportation costs, the C-SDVRP is more interesting than the other three strategies to deliver multiple commodities. This problem arises in several real-life situations, for instance in the delivery of groups of products: dairy products, fresh fruits, or vegetables to supermarkets, catering services or restaurants. Each group of products represents a commodity. These commodities can be mixed in the same vehicle. For the customer, it is acceptable to have more than one delivery, but splitting the delivery of a specific commodity (a group of products) is not practical at all. Few commodities are considered. The number of deliveries for a customer is therefore acceptable. Moreover, considering to split deliveries is beneficial for the entire logistic system since it reduces transportation costs.

Despite its practical relevance, the C-SDVRP has received very little attention. The C-SDVRP was first introduced by [START_REF] Archetti | Multicommodity vs. single-commodity routing[END_REF]. A branch-and-cut algorithm was proposed by the authors and is able to solve to optimality 25 out of 64 small instances (15 customers) within 30 minutes. [START_REF] Archetti | Multicommodity vs. single-commodity routing[END_REF] also proposed a heuristic method to tackle this problem. The heuristic consists in making copies of each customer (one for each required commodity) and uses an open-access injection-ejection algorithm for the CVRP. This solving method is simple but does not seem to be very efficient. Indeed, customer replicas share the same location, so the resulting capacitated VRP has many equivalent solutions. We use an example to highlight this point.

Figure 1 shows an example where customers require multiple commodities. A square indicates the location of the depot and circles represent the locations of the three clients. The number of commodities is three. The number on each edge corresponds to the associated travel cost, and the numbers in the dotted ellipses are the demands for each commodity. To increase readability, a different colored ellipse represents each commodity. For instance, for each customer, the number in the red ellipse is the demand of commodity one required by this customer. Each vehicle has a capacity of Q = 10 units.

According to the heuristic used by [START_REF] Archetti | Multicommodity vs. single-commodity routing[END_REF], the C-SDVRP is solved as a CVRP where the set of customers contains all the customer replicas. Figure 2 shows two solutions of this example.

Solution 1 is composed of three routes to serve all of the demand. Route 1 serves the whole demand for customer 1 (6 units) and commodity 1 of customer 2 (3 units) for a total of 9 units. Route 2 delivers the remaining of the demand of customer 2 (5 units) and commodities 1 and 2 (5 units) of customer 3. The remaining demand for customer 3 (1 unit) is delivered in route 3. The three routes cost 11.5, 11.5 and 8 respectively. The total cost is then 31. Solution 2 is also composed of three routes. Route 1 is the same as in solution 1. Route 2 delivers the remaining of the demand of customer 2 (5 units) and commodity 1 (3 units) of customer 3. The remaining demand for customer 3 (3 units) is delivered in route 3. The three routes cost 11.5, 11.5 and 8 respectively, for a total cost of 31.

These two solutions are different. However, in essence, the two solutions are equivalent in that customers' replicas share the same location. Taking this specific feature of the problem into account to avoid exploring many equivalent solutions is a real challenge.

In a more recent work [START_REF] Archetti | A branch-price-and-cut algorithm for the commodity constrained split delivery vehicle routing problem[END_REF] proposed an extended formulation for the C-SDVRP and developed a branch-price-and-cut algorithm. This algorithm has limitations for solving large-scale instances: optimal solutions were obtained with up to 40 customers and 3 commodities per customer within 2 hours. We are not aware of any other study in the literature solving the C-SDVRP.

This paper aims at proposing an efficient heuristic to tackle medium and large sized C-SDVRP instances. To this end, we propose an adaptive large neighborhood search (ALNS) heuristic taking into account the specialty of the C-SDVRP. ALNS is an efficient metaheuristic proposed by [START_REF] Ropke | An adaptive large neighborhood search heuristic for the pickup and delivery problem with time windows[END_REF] and extends the large neighborhood search (LNS) heuristic [START_REF] Shaw | A new local search algorithm providing high quality solutions to vehicle routing problems[END_REF][START_REF] Shaw | Using constraint programming and local search methods to solve vehicle routing problems[END_REF] by allowing multiple destroy and repair methods to be used within the same search. Recently, ALNS has been successfully applied to the capacitated VRP [START_REF] Sze | A hybridisation of adaptive variable neighbourhood search and large neighbourhood search: Application to the vehicle routing problem[END_REF] and to many variants of the VRP [START_REF] Azi | An adaptive large neighborhood search for a vehicle routing problem with multiple routes[END_REF][START_REF] François | Large neighborhood search for multi-trip vehicle routing[END_REF][START_REF] Masson | An adaptive large neighborhood search for the pickup and delivery problem with transfers[END_REF][START_REF] Sze | The cumulative capacitated vehicle routing problem with min-sum and min-max objectives: An effective hybridisation of adaptive variable neighbourhood search and large neighbourhood search[END_REF].

The contributions of this paper are as follows. First, a new heuristic method for the C-SDVRP is proposed. As mentioned earlier, if a customer is replicated as many time as the required commodities, many equivalent solutions exist. Hence, to avoid this pitfall, the proposed method explicitly takes this feature into account. To improve a solution, we adapt existing local search (LS) operators to deal with a customer as a whole (i.e., with the whole demand he/she requires) or only as a part (i.e., with a single commodity he/she requires). Second, in order to further improve the quality of a new and better encountered solution, a mixed integer programming (MIP) based operator is developed. Finally, we provide a large number of new best-known solutions for medium and large sized C-SDVRP instances up to 100 customers and 3 commodities per customer within about 10 minutes of computing time.

The rest of this paper is organized as follows. Section 2 defines the C-SDVRP. The proposed algorithm for the C-SDVRP is described in Section 3. Section 4 reports the computational results. Section 5 concludes the paper and suggests new directions for future research.

Problem definition

The C-SDVRP is defined on a directed graph G = (V, A) in which V = {0} V C is the set of vertices, and A is the set of arcs. More precisely, V C = {1, ..., n} represents the set of customer vertices, and 0 is the depot. A cost c ij is associated with each arc (i, j) ∈ A and represents the non-negative cost of travel from i to j. Let M = {1, ..., M } be the set of available commodities. Each customer i ∈ V C requires a quantity d im ≥ 0 of commodity m ∈ M. Note that a customer i ∈ V C may request any subset of commodities, namely d im may be equal to zero for some m ∈ M.

A fleet of identical vehicles with capacity Q is based at the depot and is able to deliver any subset of commodities.

The problem is to find a solution that minimizes the total transportation cost, and that involves two related decisions such as finding a set of vehicle routes that serve all customers and selecting the commodities that are delivered by a vehicle route to each customer. Moreover, each solution must be such that:

1. each route starts and ends at the depot; 2. the total quantity of commodities delivered by each vehicle does not exceed the vehicle capacity; 3. each commodity requested by each customer must be delivered by a single vehicle; 4. the demands of all customers need to be satisfied.

We use the example presented in Figure 1 to illustrate an optimal solution of a C-SDVRP instance. The solution is provided in Figure 3. In the C-SDVRP case, two vehicle routes are required to deliver all the commodities required by the customers. One route (black line) delivers all the commodities of customer 1 (6 units) and delivers commodities 1 and 3 of customer 3. The cost of this route is 13. The other route (purple line) delivers all the commodities of customer 2 (8 units) and commodity 2 of customer 3 (2 units). The cost of this route is 11.5. The total cost of the solution is 24.5. For the solutions obtained using the other three delivery strategies (separate routing, mixed routing and split delivery mixed routing), the interested reader is referred to [START_REF] Archetti | Multicommodity vs. single-commodity routing[END_REF].

Figure 3: An optimal solution for the C-SDVRP instance proposed in Figure 1.

Adaptive large neighborhood search

In order to tackle the C-SDVRP for medium and large instances, we propose a heuristic method based on the ALNS framework of [START_REF] Ropke | An adaptive large neighborhood search heuristic for the pickup and delivery problem with time windows[END_REF]. We make use of local search, and we develop a mathematical programming based operator to improve the quality of solutions.

Due to the characteristics of the problem under study, one node can be duplicated as many times as the number of commodities required by the associated customer [START_REF] Archetti | Multicommodity vs. single-commodity routing[END_REF]. With each duplicated node, we then associate the demand of the customer for the corresponding commodity. However, the simple duplication of customers without further consideration of the customer location can produce several equivalent solutions as illustrated previously in the paper. To enhance the performance of the algorithm, we explicitly consider customer replications for each commodity and the customer as a single entity associated with its total demand. To this intent and for the sake of clarity, in the following, we will call the duplicated nodes customer-commodity, and we will use the term customer to refer to the customer associated with the total demand.

We represent a solution of C-SDVRP as the set of routes needed to serve all customers. In order to take into account the specific features of the C-SDVRP, a route can be represented by (1) a sequence of customers, or by (2) a sequence of customer-commodities. Note that, because a customer can be delivered by several vehicles, in the representation with customers, it is possible that a customer appears in several routes (with different commodities).

The second representation gives more flexibility but increases in complexity. To clarify this, we consider the case of removing a customer from a solution. In the first case, when removing a customer from a route, the customer associated with all the commodities delivered by the route is removed. In the second case, it is possible to remove only one commodity. We illustrate the two representations in Figures 4 and5. We use the example presented in Figure 1. Figure 4 shows a set of routes represented as two sequences of customers. If customer 1 is removed from route 1, then customer 1 with all the three commodities is removed. If customer 3 in route 2 is removed, then customer 3 with commodity 2 is removed. Once a customer is removed, the remaining capacity of this route increases and the cost decreases. Note that even if customer 3 has been removed from route 2, he is still present in route 1.

Figure 5 shows a set of routes represented as two sequences of customer-commodities. In order to better understand the feature of C-SDVRP, we hide the circle which represents the customer. In fact, the two routes imply the same solution as shown in Figure 4. If commodity 2 of customer 3 in route 2 is removed, then the remaining capacity of this route increases and the cost decreases. However, if one commodity (like commodity 1) of customer 1 is removed, then the remaining capacity of this route increases but the cost of this route does not change.

(a) Initial solution.

(b) Solution after removing customer 1 from route 1 and customer 3 from route 2. As mentioned above, in this paper we present operators for customers and customer-commodities that work based on the specific need. It is important to note that to deal with both customers and customer-commodities, each route in the solution has two concurrent representations. In the first (a) Initial solution.

(b) Solution after removing commodity 1 of customer 1 from route 1 and commodity 2 of customer 3 from route 2. representation, each route contains a sequence of customers, and a set of commodities is associated with each customer (see Figure 4 for an example). In the second representation, each route contains a sequence of customer-commodities (see Figure 5 for an example). When using an operator, the corresponding representation (customers or customer-commodities) is used. The following considerations are then taken into account when translating one representation to the other:

• In the second representation with customer-commodities, when a route contains several commodities of the same customer, it is always optimal to group them (since there is zero cost to travel between these customer-commodities). Thus, in the first representation, a customer will not appear twice on the same route.

• When dealing with the second representation with customer-commodities, it is possible that a customer appears in different routes (with different commodities). Hence, it is possible in the first representation with customers that a customer appears in several routes (e.g., customer 3 in Figure 4(a)).

• When dealing with the first representation with customers, moving a customer means to move that customer with the associated commodities in the current route. These commodities may be a subset of the commodities required by this customer.

General framework

The basic idea of ALNS is to improve the current solution by destroying it and rebuilding it. It relies on a set of removal and insertion heuristics which iteratively destroy and repair solutions. The removal and insertion heuristics are selected using a roulette wheel mechanism. The probability of selecting a heuristic is dynamically influenced by its performance in past iterations [START_REF] Pisinger | Large neighborhood search[END_REF]. A sketch of the method is outlined in Algorithm 1.

In Algorithm 1, s best represents the best solution found during the search, while s is the current solution at the beginning of an iteration. The cost of a solution s is denoted by f (s).

A removal heuristic h rem and an insertion heuristic h ins are applied to the current solution s. We indicate by s rem and s ins the intermediate solutions obtained after applying h rem and h ins respectively.

h rem removes and h ins inserts ρ customers or customer-commodities, where ρ is a parameter that varies between ρ min and ρ max . We adapt the strategy proposed by [START_REF] François | Large neighborhood search for multi-trip vehicle routing[END_REF] to set the value of ρ: we slightly destroy the current solution when a new solution has just been accepted (small values of ρ), and we increase the value of ρ proportionally to the number of iterations without improvements. The probabilities of selecting a removal or an insertion heuristic are dynamically adjusted during the search (Section 3.8).

Once the heuristics h rem and h ins have been applied, the solution obtained s ins is possibly improved by applying a local search. The resulting solution is denoted by s ′ (Section 3.3).

We allow the insertion heuristics to propose solutions that violate the capacity constraint. This is done with the aim of reducing the number of routes in the solution. Infeasibility is then penalized in the objective function by adding a factor proportional to the violation. Details on the penalization of the load violation are given in Section 3.9. We then try to recover feasibility by applying the local search.

Whenever a new best solution is obtained, a mathematical programming based operator (MPO) is applied to further improve the new best solution (Section 3.7). This can be seen as an intensification phase of the algorithm.

The new solution s ′ is then subject to an acceptance rule. If accepted, the new solution becomes the current solution. Otherwise, the current solution does not change. This is repeated until a stopping criterion is met and the best solution found s best is returned.

Algorithm 1 Adaptive large neighborhood search framework.

1: s init ←generate an initial feasible solution using split procedure and LS 2: ρ ← ρ min , s ← s init , s best ← s init 3: repeat 4:

Roulette wheel: select a removal heuristic h rem and an insertion heuristic h ins 5:

Destroy: s rem ← remove ρ customers (or customer-commodities) from s applying h rem 6:

Repair: s ins ← insert removed customer-commodities into s rem applying h ins 7:

Improve: s ′ ← improve solution s ins with local search (Algorithm 2)

8: if f (s ′ ) < f (s best ) then 9:
s best ← improve solution s ′ with MPO (Section 3.7) In the following, each component of this algorithm is explained in detail.

Initial solution

Let n cc be the number of customer-commodities. A feasible initial solution is constructed as follows. First, we randomly determine a sequence of customer-commodities and we construct a giant tour S = (S 0 , S 1 , . . . , S ncc ), where S 0 represents the depot and S i is the i th customer-commodity in the sequence. Then, we apply a split procedure to get a feasible solution. This procedure is inspired by the works of [START_REF] Beasley | Route-first cluster-second methods for vehicle routing[END_REF] and [START_REF] Prins | A simple and effective evolutionary algorithm for the vehicle routing problem[END_REF]. It works on an auxiliary graph H = (X , A cc , Z), where X contains n cc + 1 nodes indexed from 0 to n cc , where 0 is a dummy node and node i, i > 0, represents customer-commodity S i . A cc contains one arc (i, j), i < j, if a route serving customer-commodities S i+1 to S j is feasible with respect to the capacity Q of the vehicle. The weight z ij of arc (i, j) is equal to the cost of the trip that serves (S i+1 , S i+2 , . . . , S j ) in this same order. The optimal splitting of the giant tour S corresponds to a minimum cost path from 0 to n cc in H. Finally, this feasible solution is improved by applying local search (Algorithm 2).

Local search

In order to improve a solution, a local search procedure (LS) is applied. LS is based on a set of classical operators that work on customers and on customer-commodities.

Let us first introduce some notation that will be useful in the remaining part of the section. Let u and v be two different nodes. They are associated with a customer or a customer-commodity, or one of them may be the depot depending on the operator. These nodes may belong to the same route or different routes. Let x and y be the successors of u and v in their respective routes. R(u) denotes the route that visits node u.

Operators on customers

Here we present the operators that are defined for customers. These operators consider a customer together with all the commodities delivered to this customer in a given route. The different operators are illustrated in Figure 6 (where we only represent the intra-route cases). Insert customer: this operator removes a customer u and inserts it after customer v. Swap customers: this operator swaps the positions of customer u and customer v.

2-opt on customers: if R(u) = R(v), this operator replaces (u, x) and (v, y) by (u, v) and (x, y).

Operators on customer-commodities

Here, we present the operators that are defined for customer-commodities. The different operators are illustrated in Figure 7.

Insert customer-commodity : this operator removes a customer-commodity u then inserts it after customer-commodity v.

Swap customer-commodities: this operator swaps customer-commodity u and customer-commodity v.

We can notice that the insert and swap operators for customer-commodities allow some moves that are not feasible when we consider only customers. In most cases, the insertion of a customer or of only one of the commodities to be delivered has the same cost. However, during the search, we allow the violation of the vehicle capacity, and the infeasibility is penalized (see Section 3.9). In an infeasible solution with some overcapacity in a route, it may be possible that inserting customers in other routes also leads to an infeasibility; while inserting only a customer-commodity decreases the infeasibility or leads to a feasible solution. Similarly, a swap of customers may not be feasible (or lead to an increase in the cost) due to the vehicle capacity restriction. At the opposite, when the commodities of some customers are delivered in several routes, swapping customer-commodities may be feasible with respect to the capacity and decreases the solution cost. As an example, in Figure 7, case (b), v and w are two commodities for the same customer, while u is a commodity for another customer who requires two commodities. Swapping u with both v and w may not be possible because of vehicle capacity while swapping customer-commodities u and v decreases the cost of the solution.

Note that we do not consider a 2-opt operator based on customer-commodities. Since this operator works on elements of the same route, it is never beneficial to split apart customer-commodities associated to the same customer. It follows that its behavior would be the same as the operator 2-opt on customers.

LS is applied when the split procedure has generated an initial solution or when the removal and insertion heuristics have modified the solution. In the first case, we do not allow the solution to be infeasible with respect to the vehicle capacity. In the second case, the insertion heuristic can produce infeasible solutions. As a consequence, LS (and thus the operators) may have to address the infeasibility of the current solution. LS is depicted in Algorithm 2. Given a solution, first, four local search operators are applied: insert customer, insert customer-commodity, swap customers and swap customer-commodities. They are invoked iteratively until there is no further improvement. Then the operator 2-opt of customers is applied. If it improves the solution, we reiterate with the first four operators. After exploring the neighbor defined by each operator, the move that improves the most is implemented. When all the local search operators fail, routes of the current solution are concatenated and split algorithm is applied again. This strategy is inspired by [START_REF] Prins | A GRASP × evolutionary local search hybrid for the vehicle routing problem. Bio-inspired algorithms for the vehicle routing problem[END_REF]. If this provides a better solution, the whole procedure is repeated.

Removal heuristics

This section describes the set of removal heuristics we propose to destroy the current solution. Heuristics Shaw removal and worst removal use the representation of a solution with customers, while random removal can be applied to both representations with customers and customer-commodities. Another operator that randomly removes one route is also considered. Obtain s ′ by applying the insert customer operator; s = s ′

6:

Obtain s ′ by applying the insert customer-commodity operator; s = s ′

7:

Obtain s ′ by applying the swap customers operator; s = s ′ 8:

Obtain s ′ by applying the swap customer-commodities operator; s = s ′ 9:

until No improvement has been obtained 10:

Obtain s ′ by applying the 2-opt of customers operator; s = s ′ 11:

until No improvement has been obtained

12:

Obtain a giant tour by concatenation of trips 13:

Apply the split procedure 14: until No improvement has been obtained specified criterion (e.g., location or demand). When customers with really different characteristics are removed, it is likely that each customer is then reinserted at the same position in the solution. Hence, by removing similar customers, Shaw removal aims to provide a different solution once an insertion heuristic has been applied.

Here, we define similarity between two customers as the distance between these two customers. The heuristic works as follows: a first customer i is randomly selected and removed. We then compute the similarity between customer i and the other customers (here, it is the distance), and sort the customers in a list L, according to the similarity with customer i. A determinism parameter p d (p d ≥ 1) is used to have some randomness in the customer selection to be removed in L, with a higher probability for the firsts customers. The removed customer then plays the role of customer i and the procedure is repeated until ρ customers have been removed. The interested reader can find a detailed Shaw removal pseudocode in [START_REF] Ropke | An adaptive large neighborhood search heuristic for the pickup and delivery problem with time windows[END_REF].

Worst removal: this heuristic aims at removing the customers who induce a high cost in the solution. More precisely, at each iteration, we first calculate for each customer the cost decrease if it is removed from the solution. Then customers are sorted in decreasing order according to these values. As in Shaw removal, a determinism parameter p d controls the randomization in the choice of the worst customer to remove.

Random removal: this removal heuristic randomly chooses ρ customers and removes them from the current solution. It can also be applied with customer-commodities, by randomly removing ρ customer-commodities.

Route removal: in this removal heuristic, an entire route from the current solution is randomly selected, and all the customer-commodities on this route are removed.

Insertion heuristics

In this section, we describe the insertion heuristics implemented in the proposed ALNS algorithm. In this work, all insertion heuristics consider customer-commodities.

Greedy insertion: in this insertion heuristic, at each iteration, for each removed customercommodity i, we first compute the best insertion cost ∆f 1 i : cost of inserting the customer-commodity at its best position into the solution (i.e., the insertion that minimizes the increase of the cost of the solution). Then, the customer-commodity with the minimum insertion cost is selected to be inserted at its best position. After each iteration, the insertion costs of the remaining removed customercommodities are recomputed. This process stops when all customer-commodities have been inserted.

Regret insertion: this insertion heuristic chooses, at each iteration, the removed customer-commodity which produces the biggest regret if it is not inserted at its best position at the current iteration. In regret-k heuristic, at each iteration, we first calculate, for each removed customercommodity i, ∆f 1 i the cost of inserting i at its best position, and ∆f η i (η ∈ {2; ...; k} the cost of inserting i at its η th best position. Then, for each customer-commodity i the regret value is computed as:

reg i = k η=2 (∆f η i -∆f 1 i ).
This represents the extra cost if i is not inserted at the current iteration in its best position. Finally, the customer-commodity with highest regret value reg i is inserted at its best position into the solution. The heuristic continues until all customer-commodities have been inserted. For regret-k insertion heuristics in this work, we have considered the values of k to be two and three.

Random insertion: this insertion heuristic randomly chooses a removed customer-commodity, and randomly chooses the insertion position in the solution.

Note that when inserting customer-commodities, violations of vehicle capacity are allowed and penalized in the cost function (see Section 3.9). However, we also impose a maximum capacity violation on each route. Hence, it is possible that a customer-commodity cannot be inserted in any route of the current solution. In this case, we create one additional route which includes this customer-commodity.

Acceptance and stopping criterion

When the removal, insertion and LS steps have been applied, we use a simulated annealing criterion to determine if the new solution obtained s ′ is accepted. However, a deterministic decision rule is applied in two cases. At each iteration of the algorithm, if s ′ has a lower cost than the current solution s (f (s ′ ) < f (s)), then s ′ is accepted. The solution s ′ is rejected if the costs f (s ′ ) and f (s) are equal. We reject solutions with the same cost in order to avoid working with equivalent solutions where some customer-commodities belonging to the same customer have been exchanged.

When

f (s ′ ) > f (s), then s ′ is accepted with probability e -(f (s ′ )-f (s))/T
, where T > 0 is the temperature. As proposed in [START_REF] Ropke | An adaptive large neighborhood search heuristic for the pickup and delivery problem with time windows[END_REF], the initial temperature is set such that a solution which is w% worst than the initial solution s init is accepted with a probability p accept . More formally, T is chosen such that e -(w•f (s init ))/T = p accept . Then, at each iteration of the ALNS, the temperature T is decreased using the formula T = T • γ, where γ ∈ [0, 1] is the cooling factor.

The stopping criterion for the whole procedure is a fixed number of ALNS iterations.

Mathematical programming based operator to reassign commodities

When a new best solution is identified, we intensify the search by applying a mathematical programming based operator (MPO). The main purpose is to assign the visits to a specific customer among the solution routes in a different way by solving a capacitated facility location problem.

We use Figures 8 (a) and (b) to explain the idea behind MPO. In the example, we assume the vehicle capacity is 10 units. The number d im in the ellipse reflects the demand of commodity m required by customer i. We focus on the commodities of customer 2: the ellipses with a solid line in Figure 8. We assume that Figure 8 (a) is a solution obtained after LS. Customer 2 has two commodities: the first is delivered on route 2, and the second is delivered on route 3. Inserting or swapping one of these two customer-commodities does not provide a better solution. However, we can notice that the deliveries to customer 2 can be reassigned to the first route and the total cost decreases, as shown in Figure 8 (b). The reader can notice that the operators implemented in LS do not consider this kind of movements.

Let us introduce the notation that we need to formally present MPO. First, we assume that customer i is considered. We indicate by: • M i the set of commodities required by customer i;

• s i the solution obtained from the current solution by removing all the visits to customer i; • R i the set of routes in s i ;

• c i r the cost to insert (best insertion) customer i in route r ∈ R i ;

• Q i r the remaining capacity in route r ∈ R i .
Then, we introduce the following binary decision variables:

x i mr =
1 if the delivery of commodity m of customer i is assigned to route r ∈ R i ; 0 otherwise.

x i r =
1 if at least the delivery of one commodity required by customer i is assigned to route r ∈ R i ; 0 otherwise.

The mathematical program that we solve in order to apply MPO is the following:

(IP M P O ) min r∈R i c i r x i r (1) s.t. r∈R i x i mr = 1, ∀ m ∈ M i (2) m∈M i d im x i mr ≤ Q i r x i r , ∀ r ∈ R i (3) x i mr ∈ {0, 1}, ∀ m ∈ M i , r ∈ R i (4) x i r ∈ {0, 1}, ∀ r ∈ R i (5)
This mathematical program corresponds to a capacitated facility location problem, where only the costs related to the inclusion of a new route in the solution (the fixed cost) are taken into account. The objective function (1) aims to minimize the total insertion cost. Constraints (2) require that the delivery of each commodity (i.e., the delivery of each customer-commodity) must be assigned to one route. Constraints (3) impose that the total quantity of commodities assigned to a selected vehicle does not exceed its remaining capacity. Constraints (4)-( 5) define the decision variables.

(IP M P O ) is solved for each i ∈ V C , but only the reassignment of visits associated with the highest cost reduction is effectively implemented.

Adaptive weight adjustment

A roulette wheel is used to give more or less importance to the removal and insertion heuristics to be used. The procedure implemented is based on the principles described in [START_REF] Ropke | An adaptive large neighborhood search heuristic for the pickup and delivery problem with time windows[END_REF] and [START_REF] François | Large neighborhood search for multi-trip vehicle routing[END_REF]. This last variant includes a normalization process for the score. The main difference in our approach is that removal and insertion heuristics are not selected independently. A pair of removal and insertion heuristics h p = {h rem ; h ins } is chosen in each iteration. We denote by H p the set of pair of heuristics to be used. Each pair h p is associated with a weight ω hp , a score π hp , and θ hp the number of times that pair h p has been used. Initially, all pairs of heuristics have the same weight and p∈Hp ω hp = 1.

We define a segment as a fixed number of ALNS iterations in the proposed algorithm. During a segment, the weights of all pairs are kept constant. Before starting a new segment, for each pair h p ∈ H p , the score π hp and the number of times the pair is used θ hp are reset to 0. During a segment, each time a pair h p is used, θ hp is increased by 1, and if the new solution s ′ is accepted after using pair h p , the score π hp is updated according to: π hp ← π hp + σ µ , where σ µ (µ ∈ {1; 2; 3}) reflects different cases regarding the score change π hp . That is, the score π hp is increased by σ 1 when s ′ is a new best solution, or σ 2 when s ′ is a new improved solution (f (s ′ ) < f (s)) but not a new best solution, or σ 3 when s ′ is not an improved solution but has been accepted according to the simulated annealing criterion. The values σ µ (σ µ ∈ [0, 1]) are normalized to satisfy

σ 1 + σ 2 + σ 3 = 1.
At the end of each segment, we update all the weights of the pairs of heuristics based on the recorded scores. First, the score π hp is updated as: π hp ← π hp θ hp , where θ hp is the number of times that pair h p was used in this segment. If θ hp = 0, we set π hp to the same value as in the previous segment. Then, the scores of all pairs of heuristics are normalized:

πhp = π hp h∈Hp π h . (6)
Let ω hp,q be the weight of pair h p used in segment q, and λ ∈ [0, 1] a factor which determines the change rate in the performance of the pair of heuristics. At the end of segment q, the weight of all pairs of heuristics h p to be applied in segment q + 1 is updated as:

ω hp,q+1 = (1 -λ)ω hp,q + λπ hp .
(7)

Infeasibility penalization scheme

In our implementation of the ALNS, we allow some violations of the vehicle capacity in order to reduce the number of routes. Let K init be the number of vehicles used in the initial solution s init . Then, the vehicle capacity Q can be extended by an amount of Q extra = Q K init . If a vehicle delivers more than Q units of products, we penalize the infeasibility by adding to the solution cost a term proportional to the load violation which is βl(s), where l(s) is the total load violation of the solution s and β is the penalty rate.

The penalty rate β is related to capacity violations. Initially, β is set equal to a minimum value β min computed as

β min = 10 • f (s init ) i∈V C ,m∈M d im
, where f (s init ) is the cost of the initial solution obtained after applying the split procedure and LS. This ratio approximates an average transportation cost per unit of demand.

The penalty rate β is dynamically modified during the search since a high rate may eliminate infeasibility quickly, but may also prevent searching other promising regions. We keep track of the number of infeasible solutions obtained during the last consecutive I iterations of the ALNS algorithm. If I inf infeasible solutions are obtained consecutively, the value of β is increased to 2β. Similarly, if I f eas feasible solutions are generated consecutively, the value of β is decreased to max {β min ; β/2}.

Computational experiments

In this section, we present the results obtained by the proposed ALNS heuristic. The algorithm was implemented in C++ and ran on an Intel (R) Core(TM) i7-4600U, 2.10GHz, and 16GB of RAM.

We first describe in Section 4.1 the test instances on which we evaluate our algorithm. In Section 4.2 we report the values of the parameters that we set for the ALNS. Then, in Section 4.3 we perform a sensitivity analysis to determine which insertion and removal heuristics lead to the best performance of the proposed algorithm. A set of settings is tested on a subset of instances, and the best is used to perform all other computational tests. In Section 4.4, we analyze the effect of the number of iterations on the computational time and the solution quality. Then, Section 4.5 reports the results on the testbed. The effectiveness of the MPO is studied in Section 4.6, while in Section 4.7 we validate the effectiveness of integrating a local search in the large neighborhood search. We examine how the computational time of our method varies according to the instance size in Section 4.8. Finally, Section 4.9 provides managerial insights about the split customers that are delivered by more than one vehicle.

Instances

To assess the efficiency of our algorithm, we perform computational experiments on the benchmark instances proposed by [START_REF] Archetti | Multicommodity vs. single-commodity routing[END_REF]. These instances are built from the R101 and C101 Solomon instances for the VRP with Time Windows [START_REF] Solomon | Algorithms for the vehicle routing and scheduling problems with time window constraints[END_REF]), from which the customer locations are kept. From instances R101 and C101, several instances for the C-SDVRP were generated by considering the n first customers, with n ∈ {15, 20, 40, 60, 80, 100}.

In the testbed, up to 3 commodities are taken into account. The number of commodities is indicated by M (M ∈ {2; 3}). Each commodity is required by a customer with a probability p of 0.6 (on average each customer needs 2 out of 3 commodities) or with probability p equal to 1 (each customer requires all commodities). The quantity of each commodity required by a customer varies within the intervals ∆ = [1, 100] or ∆ = [40,60]. A last parameter α ∈ {1.1, 1.5, 2, 2.5} is used to determine the vehicle capacity from the original one in Solomon instances. We indicate by P = (I, M, p, ∆, α) (where I ∈ {R101, C101}) the set of parameters used by [START_REF] Archetti | Multicommodity vs. single-commodity routing[END_REF] to generate instances.

When n = 15, 64 instances are available, one for each combination of parameters in set P. These are referred as small instances. When n ∈ {20, 40, 60, 80}, 80 mid-size instances are available: 20 instances for each value of n. For the mid-size instances, the combination of parameters in P is restricted to M = 3, α = 1.5 and ∆ = [1, 100]. Hence I and p can take two values each, leading four combinations of parameters. For each combination, five instances have been generated, leading to 20 instances in total. When n = 100, 320 large instances are available, that means five instances for each combination of parameters in set P.

In the following sections, we present the parameter setting for our algorithm and the results obtained. We provide detailed results on each of the 464 instances in the Appendix A. Due to the high number of instances, we only report average results. In particular, we present results for each group of instances, where a group is defined by a quartet (n, I, M, p). Results for a group are averaged on the values of α and ∆.

We summarize all the notations used to present the results in Table 1.

ALNS parameters

In this section, we present the values of the parameters set in our ALNS algorithm. The probabilities to select the pairs of removal and insertion heuristics are updated after a number of iterations called a segment. In our implementation, we define a segment as 100 iterations. To set the values of σ 1 , σ 2 , σ 3 , preliminary experiments were carried out on instances with n = 80 customers for some combinations. The best results were obtained when σ 1 = 0.7, σ 2 = 0.1, σ 3 = 0.2. We set the reaction factor λ that appears in Equation (7) to 0.5 as proposed in [START_REF] Masson | An adaptive large neighborhood search for the pickup and delivery problem with transfers[END_REF] (instead of λ = 0.1 as in [START_REF] Ropke | An adaptive large neighborhood search heuristic for the pickup and delivery problem with time windows[END_REF]) to ensure higher reactivity when performing fewer iterations. The determinism parameter p d in Shaw removal and worst removal heuristics is equal to 6 as in [START_REF] Ropke | An adaptive large neighborhood search heuristic for the pickup and delivery problem with time windows[END_REF].

The acceptance of a new current solution is based on a simulated annealing criterion. The initial value of the temperature T is set so that a solution that is w% worse than the current solution is accepted with probability p accept with w = 0.35 and p accept = 0.7. In addition, we set the cooling factor γ to 0.999.

In order to destroy a solution, we need to determine the number ρ of customers (or customercommodities) to be removed. We follow the scheme proposed by [START_REF] François | Large neighborhood search for multi-trip vehicle routing[END_REF]: small moves are applied when a new solution has just been accepted, while large moves are applied when no new solution has been accepted in the most recent iterations. The value of ρ evolves in the interval [ρ min , ρ max ]. For small instances (n = 15), we set ρ min = N/2, ρ max = N , while for the other instances, we set ρ min = N/10 and ρ max = N/4. When the removal heuristic is customer-based, then N represents the number of customers (N = n). When the removal heuristics are defined for customer-commodities, N represents the number of customer-commodities (N = n cc ).

To update the penalization rate for capacity violations, we consider the number of infeasible and feasible solutions that are obtained consecutively. These values are respectively set to I inf = 50 and I f eas = 5 after performing preliminary experiments.

Efficiency assessment for the removal and insertion heuristics

In Sections 3.4 and 3.5, we propose several removal and insertion heuristics to be used inside the ALNS framework. Before testing the proposed algorithm on the set of instances, we use a subset of instances in order to determine the removal and insertion heuristics that are used inside the ALNS framework. We also consider configurations with only one removal and one insertion heuristic, which corresponds to designing a large neighborhood search (LNS). This permits to emphasize the efficiency of individual insertion and removal heuristics, and also to point out the benefit of the adaptive approach included in the design of the ALNS when choosing the removal and insertion heuristics.

This analysis is performed to tune the proposed algorithm. We only use the 20 instances with n = 80 customers. A summary of these experiments is reported in Table 2. The first column shows which configuration is considered (LNS or ALNS). The second column enumerates the configurations that we tested. Columns 3 to 7 indicate which removal heuristics are used in a specific configuration. Among them, 'RandC' and 'RandCC' represent the random removal heuristic considering customers and customer-commodities, respectively. Columns 8 to 11 indicate which insertion heuristics are used. In columns 12 to 14 we report the following average statistics for each configuration: the gap with the best-known solutions, the computational time in seconds, the number of new best-known solutions (see Table 1). In the first five configurations, we study the influence of the removal heuristic. To this intent, the insertion heuristic is kept fixed. Using the Shaw removal heuristic, the LNS is able to improve all the best-known solutions.

Configurations 6 to 8 point out the impact of the insertion heuristic. We fixed the removal heuristic to the Shaw removal since it provided the best results in the previous experiments. These tests show that all insertion heuristics perform well except the random insertion. To further evaluate the behavior of the random insertion heuristic we consider configuration 9 where we modify the removal heuristic. It can be observed that the results do not improve. Among the LNS configurations, combining the Shaw removal with the regret-3 insertion provides the best results.

In the second part of Table 2, we consider different configurations for the ALNS framework, with several removal and insertion heuristics. Configuration 10 shows the performance of the ALNS algorithm using all proposed heuristics. In other configurations, only a subset of them is invoked in the ALNS. As we can see, using all the proposed removal and insertion heuristics does not provide the best results. The best configuration is configuration number 13: an ALNS with two removal heuristics: Shaw and random removal of customers, and three insertion heuristics: greedy, regret-2 and regret-3. In the following sections, all reported computational experiments have been conducted using this configuration.

Analysis with respect to the number of iterations

We examine the impact of the number of iterations for different instance sizes. Preliminary computational experiments showed that the algorithm converges after a certain number of iterations. Thus, we aim to determine the number of iterations that would be a good compromise between the solution quality and the computational time.

For small instances and medium instances with 20 customers (mid-20), we run our algorithm with the number of iteration iter limited to 100, 1000, 3000 and 5000. We compare the results with those reported in [START_REF] Archetti | Multicommodity vs. single-commodity routing[END_REF] and [START_REF] Archetti | A branch-price-and-cut algorithm for the commodity constrained split delivery vehicle routing problem[END_REF]. In Table 3 and Table 4 we report average statistics for different values of iter: the gap with respect to the optimal values, the computational time in seconds, the number of optimal solutions obtained (see Table 1).

Results over the testbed are indicated in bold. Detailed results are provided in Appendix A.

According to Table 3, the proposed algorithm solves to optimality 57 out of 64 small instances when 100 iterations are allowed, within an average computational time of 2 seconds. When we increase the number of iterations to 3000, the number of optimal solutions reaches 63 with an average computational time of 12 seconds. The average gap with optimal solutions obtained by [START_REF] Archetti | A branch-price-and-cut algorithm for the commodity constrained split delivery vehicle routing problem[END_REF] varies from 0.10% to 0.03%. Increasing the number of iterations to 5000 does not improve the quality of the results.

When n = 20, results reported in Table 4 indicate that our method identifies an optimal solution in 5000 iterations for 18 out of the 19 instances for which [START_REF] Archetti | A branch-price-and-cut algorithm for the commodity constrained split delivery vehicle routing problem[END_REF] provided optimal values. The average gap with the best-known values is less than 0.01%. Over the whole set, the average CPU time is less than 1 minute.

Since large instances usually require more iterations to obtain high-quality results, we compare the ALNS algorithm behavior with 5000 and 10000 iterations on mid-80 instances. In Table 5, we report average results and the gap with respect to the best-known values. Detailed results are provided in Appendix A. ALNS can provide best-known solutions for the 20 mid-80 instances in the benchmark. Within 5000 iterations, the best-known solutions can be improved by about 0.79% on average and the CPU times is less than 9 minutes. When 10000 iterations are executed, the solutions are averagely 0.92% better than the best-known values. On the other side, the average computational times are almost doubled.

Computational experiments on the whole testbed

In this section, we consider the results obtained on the whole set of instances for the C-SDVRP with the designed ALNS algorithm. As determined in the previous sections, the removal and insertion heuristics are the ones of configuration 13, and the number of iterations is equal to 3000 (resp. 5000) on small (resp. medium and large) instances. The algorithm is run once on each instance.

The comparison is made with the results reported in [START_REF] Archetti | Multicommodity vs. single-commodity routing[END_REF] and [START_REF] Archetti | A branch-price-and-cut algorithm for the commodity constrained split delivery vehicle routing problem[END_REF]. [START_REF] Archetti | A branch-price-and-cut algorithm for the commodity constrained split delivery vehicle routing problem[END_REF] propose a branch-and-price algorithm that can solve to optimality instances with up to 40 customers but that can systematically close instances with up to 20 customers. When the optimal value is not available we compare to the best known solution (BKS) given either by [START_REF] Archetti | Multicommodity vs. single-commodity routing[END_REF] or [START_REF] Archetti | A branch-price-and-cut algorithm for the commodity constrained split delivery vehicle routing problem[END_REF].

For the 64 small instances, an optimal solution is known. Results reported in Table 3 indicate that our algorithm finds an optimal solution for 63 out of 64 instances. The average optimality gap is 0.03%. For the mid-20 instances (n = 20), 19 out of 20 instances have a known optimal value. The proposed algorithm can find 18 out of 19 optimal values. On average, the gap with respect to BKS (either optimal or feasible) is lower than 0.01%.

For the mid-40 instances (n = 40), average results are reported in Table 6. Only 5 optimal values are known. On average, the proposed ALNS finds solutions that are 0.26% better than the best-known solutions (either optimal or feasible). Our algorithm finds 4 out of the 5 known optimal solution while 9 new best-known solutions have been identified. The algorithm runs in less than 2 minutes on average. For instances with n ≥ 60, no optimal solution is available. For the mid-60 instances (n = 60) (see Table 7), the ALNS finds solutions that are on average 0.49% better than the BKS. Among them, 15 new best-known values are obtained. The CPU times are less than 5 minutes on average. As presented in Table 5 on mid-80 instances (n = 80), the ALNS identifies solutions that are 0.79% better than the BKS, and new best-known values are found for all instances.

Table 8 reports the results on the 320 large instances (n = 100). The ALNS finds 300 new bestknown values with an average improvement of 0.70%. The computational time is around 10 minutes which is very reasonable when considering the instance size.

Effectiveness of MPO operator in the ALNS algorithm

As described in Section 3.7 we developed a mathematical programming based operator (MPO) to re-assign commodities for one customer to the routes of a solution. Due to the increase observed in computational time consumption, we decided to use it only to improve a new global best solution further. In this section, we analyze the results on medium and large instances to prove the effectiveness of MPO in our algorithm. In Table 9, we indicate the average number of times that MPO is called for each group of instances as well as the average number of times that MPO improves the new best solution for each group of instances. Detailed results are reported in Appendix A (Table 13 and Table 14). For medium instances (namely for n = 20, 40, 60, 80), the impact of MPO on the quality of the solution increases as the instance size increases. Especially, when p = 1, MPO has a greater impact than when p = 0.6. The main reason is that, when p = 1, all customers require all the commodities while a smaller number of commodities are required when p = 0.6. As a consequence, for the same number of customers in the instance, the number of customer-commodities is larger. Therefore, when p = 1, the difficulty of solving the problem increases as well as the possibilities to reassign commodities with MPO. For large instances (n = 100), the same conclusions on the performance of MPO can be drawn.

Evaluation of the LS in the ALNS algorithm

In this section, we evaluate the importance of the local search in our ALNS algorithm. We run, on mid-80 instances, the ALNS algorithm without LS and MPO with a limit of 400000 iterations. We then compare the results obtained by the ALNS with LS and MPO on 5000 iterations. This choice is made to have a comparison fair: we remove two optimization components, but we allow a large number of iterations.

We compare the results obtained by the proposed algorithm, indicated by ALNS+LS+MPO, with the same algorithm when LS and MPO are deactivated. This last version is indicated as ALNS-LS-MPO. The performance comparison of both variants is reported in Table 10 for mid-80 and large instances.

It is clear from Table 10 that after 400000 iterations the results obtained by the ALNS without LS and MPO are of lower quality than those obtained by the original ALNS. We only obtain 15 new best-known values for the 20 mid-80 instances. Moreover, the average improvement (0.39%) does not compete with the improvement obtained with the proposed algorithm (ALNS+LS+MPO) while the computational times are similar. The same observations can be made from the results for large instances. We then conclude on the importance of the LS and MPO in the proposed ALNS algorithm.

Trend between instance size and computational time

Last, we examine how the CPU time required by the ALNS varies according to the instance size. We consider the results obtained with 5000 iterations (even for the small instances) to perform the analysis. We sort the 464 (small, medium and large) instances according to the number of customercommodities instead of the number of customers. The variation of the computational time avg.t(s) according to the number of customer-commodities avg.n cc is depicted on Figure 9. When the size of the instances increases, the average computational time significantly increases (not in a linear fashion). This behavior is not surprising since when the size of the instances increases, it takes more time to operate the LS operators, which computational complexity is O(n 2 cc ). Nevertheless, the average computational time of our ALNS algorithm for the large instances (n cc = 300) remains reasonable with less than 19 minutes.

Characteristics of split customers

As mentioned above, considering customer-commodity makes the problem more complex with equivalent solutions since several commodities are related to the same customer and then have the same location. Identifying the customers who are good candidates for being split into customer-commodities can be beneficial for decision-makers. In this section, we propose to study the characteristics of customers who are delivered by more than one vehicle. We name these customers split customers. We provide detailed results on the 20 mid-80 instances with 5000 iterations for the ALNS algorithm. Table 11 reports these detailed results for several characteristics about the split customers in the best solution obtained on each instance. The first three columns report the characteristics of the instance, and the fourth column indicates the identification number of the instance. Column nbSplit reports the number of split customers (out of 80), and columns nb2-split and nb3-split indicate the number of customers delivered by respectively 2 and 3 vehicles. Columns nbN earDepot and nbLargeDemand report the number of split customers located near the depot and the number of split customers with a large demand respectively. As [START_REF] Nagy | The Vehicle Routing Problem with Divisible Deliveries and Pickups[END_REF], we consider that a customer is located near the depot if he or she is one of the 25% of customers closest to the depot; and we consider that a customer has a large demand if he or she is one of the 25% of customers with the largest demand. Note that the demand of a customer is the sum of the demands for the commodities they need. Column nbCluster reports the number of split customers located inside a cluster of customers. As [START_REF] Nagy | The Vehicle Routing Problem with Divisible Deliveries and Pickups[END_REF], we consider a customer to be in a cluster if at least five other customers are inside its neighborhood. Two customers are neighbors if the distance between these two customers is less than 30% of the average distance between the depot and all the customers in the instance. Using this definition, in clustered instances (C101) more than 70% of customers are in a cluster, while in random instance (R101) less than 30% of customers are in a cluster. From the results in Table 11, it is clear that there are more split customers when customers require more commodities. p denotes the probability that a customer requires a commodity. When p = 0.6, the average number of split customers is 6.7, while when p = 1, there are 15.5 split customers on average. Moreover, very few split customers are delivered by three vehicles, i.e., one vehicle for each of the required commodities. In the 20 instances, there are only three split customers in this case. In addition, one important feature of split customers is to be inside a cluster. Indeed, for C101 instances, 86% of split customers are inside a cluster. For R101 instances, 46% of split customers are inside a cluster, while less than 30% of customers are in a cluster in these instances. Proximity to the depot is also an important feature for split customers since 36% of split customers are near the depot. It is less obvious to link large demands to split customers since 30% of split customers have a large demand.

Conclusions

In this paper, we presented a dedicated heuristic algorithm for the C-SDVRP. The proposed algorithm is based on an adaptive large neighborhood search framework introduced by [START_REF] Ropke | An adaptive large neighborhood search heuristic for the pickup and delivery problem with time windows[END_REF]. This is the first heuristic specifically designed with the aim to provide high-quality solutions for the medium and large size instances. According to the main feature of the C-SDVRP, i.e. the requirement of different commodities, we adapt some classical local search moves to consider either with a customer (i.e., a customer and all its commodities) or a customer-commodity (namely, a single commodity required by a customer). We developed a mathematical programming based operator to intensify the search and further improve the best solutions. The results show that our ALNS algorithm is very effective in finding high-quality solutions on large size instances. In particular, our method outperforms the algorithms proposed in [START_REF] Archetti | Multicommodity vs. single-commodity routing[END_REF] and in [START_REF] Archetti | A branch-price-and-cut algorithm for the commodity constrained split delivery vehicle routing problem[END_REF].

The proposed ALNS algorithm could then be adapted to tackle other variants of routing problems with split deliveries. One of these variants is the case with multiple depots and available quantities at each depot. In this case, because of the limited quantities available for each commodity at each depot, it is worth considering splitting deliveries to find feasible solutions. Another interesting variant of the problem is the VRP with divisible deliveries and pickup [START_REF] Gribkovskaia | General solutions to the single vehicle routing problem with pickups and deliveries[END_REF], [START_REF] Nagy | The Vehicle Routing Problem with Divisible Deliveries and Pickups[END_REF]): in this case, delivery and pickup naturally represent two different commodities, but a pickup operation increases the use of the vehicle capacity, and it could then be optimal to visit twice the same customer in one route.

Figure 1 :

 1 Figure 1: An example from Archetti et al. (2014).

Figure 2 :

 2 Figure 2: Two solutions of the instance in Figure 1.

Figure 4 :

 4 Figure 4: Two sequences of customers.

Figure 5 :

 5 Figure 5: Two sequences of customer-commodities.

Figure 6 :

 6 Figure 6: Local search operators with customers.

  (a) Insert customer-commodity. (b) Swap customer-commodities.

Figure 7 :

 7 Figure 7: Local search operators with customer-commodities.

  (a) Solution before using MPO.(b) Solution after using MPO.

Figure 8 :

 8 Figure 8: MPO for reassigning commodities.

Figure 9 :

 9 Figure 9: The avg.t(s) of the ALNS with respect to avg.n cc

Table 1 :

 1 Notations for computational results Percentage of improvement between Cost and OPT (∆ O = 100 * (Cost -OP T )/OP T ) ∆ B Percentage of improvement between Cost and BKS (∆ B = 100 * (Cost -BKS)/BKS) ∆ O/B To indicate that the value is ∆ O if an optimal value if available, and ∆ B otherwise avg.∆ O Average percentage of improvement between Cost and OPT (∆ O = 100 * (Cost -OP T )/OP T ) avg.∆ B Average percentage of improvement between Cost and BKS (∆ B = 100 * (Cost -BKS)/BKS) min.∆ O Minimum percentage of improvement between Cost and OPT (∆ O = 100 * (Cost -OP T )/OP T ) max.∆ B Maximum percentage of improvement between Cost and BKS (∆ B = 100 * (Cost -BKS)/BKS) Number of times that MPO is called for a instance nbM P Oimp Number of times that MPO improves the new best solution during the search avg.nbM P O Average number of times that the MPO is called for each group of instances avg.nbM P Oimp Average number of times that the MPO improves the current best solution for each group of instances * Indicates the solution value is optimal nbSplit Number of split customers (that are delivered by more than one vehicle)

	Symbol	Meaning
	C101/R101	Name of the original Solomon instance
	n	Number of customers
	M	Number of commodities
	p	Probability that a customer requires a commodity
	id	Instance id (5 mid and large sized instances are generated for each combination of parameters in P)
	∆	Demand range of each customer for each commodity
	α	Value for generating the vehicle capacity (see Archetti et al. (2015))
	CCm	Total number of customers that require commodity m
	ncc	Total number of customer commodities
	avg.ncc	Average total number of customer commodities in each group of instances
	nbIns	Number of instances in each group
	OP T	Optimal solution value from literature (when available)
	BKS	Best known solution value from literature (when available)
	Cost	Best solution cost found by the proposed algorithm
	∆ O
	t(s)	CPU time in second
	avg.t(s)	Average CPU time in seconds
	nbOP T	Number of optimal solutions obtained
	nbE	Number of solution values obtained equal to the best known
	nbN BK	Number of new best known solutions obtained
	nbR	Number of routes
	nbM P O
	nb2-split	Number of split customers delivered by exactly two vehicles
	nb3-split	Number of split customers delivered by exactly three vehicles
	nbN earDepot	Number of split customers that are close to the depot
	nbLargeDemand	Number of split customers with a large demand
	nbCluster	Number of split customers located inside a cluster of customers

Table 2 :

 2 LNS configurations compared to ALNS configurations.

			Removal heuristics				Insertion heuristics			Results
		conf.	Shaw	RandC	RandCC	Worst Route	Greedy	Regret-2	Regret-3	Rand	avg.∆ B	avg.t(s)	nbN BK
	LNS	1	×						×			-0.77	521.18
		2		×					×			-0.72	438.34
		3			×				×			-0.54	276.76
		4				×			×			-0.23	714.33
		5					×		×			0.01	333.60
		6	×					×				-0.61	533.30
		7	×							×		-0.78	541.81
		8	×								×	-0.07	1106.18
		9		×							×	-0.10	1146.84
	ALNS	10	×	×	×	×	×	×	×	×	×	-0.71	562.29
		11	×	×				×	×			-0.72	513.49
		12	×	×					×	×		-0.77	500.97
		13	×	×				×	×	×		-0.79	511.00
		14	×	×	×			×	×			-0.68	467.99
		15	×	×		×		×	×			-0.63	555.34
		16	×	×		×		×	×	×		-0.66	530.26
		17	×	×	×			×	×	×		-0.79	477.49

Table 3 :

 3 Impact of the number of iterations on small instances.

			instances		ALNS (100 iterations)	ALNS (1000 iterations)	ALNS (3000 iterations)	ALNS (5000 iterations)
		n	M	p	avg.ncc nbIns	avg.∆ O	avg.t(s)	nbOP T	avg.∆ O	avg.t(s)	nbOP T	avg.∆ O	avg.t(s)	nbOP T	avg.∆ O	avg.t(s)	nbOP T
	C101	15	2	0.6	22	8	0.00	0.93	8	0.00	2.82	8	0.00	7.24	8	0.00	11.65
	C101	15	2	1	30	8	0.01	1.64	7	0.00	4.45	8	0.00	10.63	8	0.00	16.86
	C101	15	3	0.6	28	8	0.00	1.13	8	0.00	3.92	8	0.00	10.40	8	0.00	17.12
	C101	15	3	1	45	8	0.47	3.72	4	0.33	9.81	7	0.21	22.26	7	0.21	34.62
	R101	15	2	0.6	22	8	0.00	1.07	8	0.00	2.89	8	0.00	7.02	8	0.00	11.05
	R101	15	2	1	30	8	0.00	2.10	8	0.00	4.75	8	0.00	10.76	8	0.00	16.99
	R101	15	3	0.6	28	8	0.00	1.26	8	0.00	3.88	8	0.00	9.65	8	0.00	15.59
	R101	15	3	1	45	8	0.35	3.36	6	0.00	8.63	8	0.00	20.21	8	0.00	31.57
			total			64	0.10	1.90	57	0.04	5.14	63	0.03	12.27	63	0.03	19.43	63

Table 4 :

 4 Impact of the number of iterations on mid-20 instances.

		instances			ALNS (100 iterations)		ALNS (1000 iterations)		ALNS (3000 iterations)		ALNS (5000 iterations)
	n	M	p	avg.ncc	nbIns avg.∆ O/B	avg.t(s)	nbOP T	avg.∆ O/B	avg.t(s)	nbOP T	avg.∆ O/B	avg.t(s)	nbOP T	avg.∆ O/B	avg.t(s)	nbOP T
	C101	3	0.6	37.4	5	0.26	2.41	3	0.10	8.11	4	0.10	20.66	4	0.00	33.38	5
	C101	3	1	60	5	0.66	5.75	1	0.11	21.13	2	0.04	49.66	4	0.00	77.14	5
	R101	3	0.6	37.4	5	0.01	3.87	4	0.00	9.60	5	0.00	22.09	5	0.00	34.87	5
	R101	3	1	60	5	0.82	10.68	1	0.16	25.66	2	0.01	54.17	3	0.01	81.72	3
		total			20	0.44	5.68	9	0.09	16.12	13	0.04	36.65	16	0.00	56.78	18

Table 5 :

 5 Impact of the number of iterations on mid-80 instances.

	instances					ALNS (5000 iterations)	ALNS (10000 iterations)
		n	M	p	avg.ncc	nbIns avg.∆ B	avg.t(s)	nbN BK	avg.∆ B	avg.t(s)	nbN BK
	C101	80	3	0.6	150.4	5	-0.77	373.36	5	-0.88	690.87	5
	C101	80	3	1	240	5	-0.75	704.20	5	-0.83	1339.55	5
	R101	80	3	0.6	150.4	5	-0.84	319.83	5	-1.06	605.74	5
	R101	80	3	1	240	5	-0.81	646.61	5	-0.93	1196.69	5
		total				20	-0.79	511.00	20	-0.92	958.21	20

Table 6 :

 6 Summary of results on mid-40 sized instances.

	instances					ALNS results
		n	M	p	avg.ncc	nbIns avg.∆	min∆	max∆	avg.t(s)	nbOP T	nbE	nbN BK
	C101	40	3	0.6	76.2	5	-0.12	-0.42	0.15	78.14	2	-	2
	C101	40	3	1	120	5	-0.76	-1.91	0.00	161.39	-	-	4
	R101	40	3	0.6	76.2	5	0.15	0.00	0.29	80.70	1	-	-
	R101	40	3	1	120	5	-0.32	-0.78	0.00	148.33	1	1	3
			total			20	-0.26	-1.91	0.29	117.14	4	1	9

Table 7 :

 7 Summary of results on mid-60 sized instances.

	instances					ALNS results
		n	M	p	avg.ncc	nbIns avg.∆ B	min∆ B	max∆ B	avg.t(s)	nbE	nbN BK
	C101	60	3	0.6	110	5	-0.56	-1.40	0.00	193.12	1	4
	C101	60	3	1	180	5	-0.38	-1.13	0.25	403.78	-	3
	R101	60	3	0.6	110	5	-0.41	-0.83	0.00	177.20	1	3
	R101	60	3	1	180	5	-0.61	-1.20	-0.14	353.24	-	5
		total				20	-0.49	-1.40	0.25	281.83	2	15

Table 8 :

 8 Summary of results on large sized instances.

	instances					ALNS results
		n	M	p	avg.ncc	nbIns avg.∆ B	min∆ B	max∆ B	avg.t(s)	nbN BK
	C101	100	2	0.6	136.4	40	-0.66	-1.69	0.24	330.02	38
	C101	100	2	1	200	40	-0.77	-1.64	0.03	569.11	39
	C101	100	3	0.6	188.4	40	-0.77	-2.04	0.73	557.46	39
	C101	100	3	1	300	40	-1.17	-2.88	0.32	1106.23	37
	R101	100	2	0.6	136.4	40	-0.53	-1.50	1.05	323.48	36
	R101	100	2	1	200	40	-0.55	-1.28	0.29	557.74	36
	R101	100	3	0.6	188.4	40	-0.53	-1.23	0.02	548.02	39
	R101	100	3	1	300	40	-0.62	-1.72	0.25	1078.70	36
			total			320	-0.70	-2.88	1.05	633.85	300

Table 9 :

 9 Effectiveness of MPO in the ALNS algorithm.

	instances			results
		n	M	p	avg.ncc	nbIns avg.nbM P O	avg.nbM P Oimp
	C101	20	3	0.6	37.4	5	2.60	0.00
	C101	20	3	1	60	5	7.00	0.20
	R101	20	3	0.6	37.4	5	4.00	0.20
	R101	20	3	1	60	5	9.80	0.20
	C101	40	3	0.6	76.2	5	8.00	0.20
	C101	40	3	1	120	5	16.40	2.00
	R101	40	3	0.6	76.2	5	12.40	0.00
	R101	40	3	1	120	5	18.60	1.80
	C101	60	3	0.6	110	5	20.80	0.20
	C101	60	3	1	180	5	30.80	1.60
	R101	60	3	0.6	110	5	18.00	0.80
	R101	60	3	1	180	5	24.60	1.60
	C101	80	3	0.6	150.4	5	29.00	1.60
	C101	80	3	1	240	5	32.00	6.40
	R101	80	3	0.6	150.4	5	23.80	1.00
	R101	80	3	1	240	5	29.60	2.40
	C101	100	2	0.6	136.4	40	27.05	0.85
	C101	100	2	1	200	40	29.65	2.00
	C101	100	3	0.6	188.4	40	31.60	1.33
	C101	100	3	1	300	40	37.28	3.33
	R101	100	2	0.6	136.4	40	27.15	0.53
	R101	100	2	1	200	40	30.15	2.10
	R101	100	3	0.6	188.4	40	32.40	2.23
	R101	100	3	1	300	40	34.48	3.58

Table 10 :

 10 Comparison between two ALNS variants on mid-80 and large instances.

			instances				ALNS+LS+MPO (5000 iterations)			ALNS-LS-MPO (400000 iterations)
		n	M	p	avg.ncc	nbIns avg.∆ B	min∆ B	max∆ B	avg.t(s)	nbN BK	avg.∆ B	min∆ B	max∆ B	avg.t(s)	nbN BK
	C101	80	3	0.6	150.4	5	-0.77	-1.48	-0.01	373.36	5	-0.70	-1.50	0.11	330.77	4
	C101	80	3	1	240	5	-0.75	-1.36	-0.24	704.20	5	-0.49	-1.08	0.02	686.68	4
	R101	80	3	0.6	150.4	5	-0.84	-1.24	-0.05	319.83	5	-0.49	-1.03	0.57	348.14	4
	R101	80	3	1	240	5	-0.81	-1.39	-0.50	646.61	5	0.11	-0.83	1.71	712.20	3
			total			20	-0.79	-1.48	-0.01	511.00	20	-0.39	-1.50	1.71	519.45
	C101	100	2	0.6	136.4	40	-0.66	-1.69	0.24	330.02	38	-0.23	-1.52	1.21	360.30
	C101	100	2	1	200	40	-0.77	-1.64	0.03	569.11	39	-0.30	-1.58	1.30	606.03
	C101	100	3	0.6	188.4	40	-0.77	-2.04	0.73	557.46	39	-0.28	-1.66	1.66	576.18
	C101	100	3	1	300	40	-1.17	-2.88	0.32	1106.23	37	-0.57	-2.47	1.20	1256.61
	R101	100	2	0.6	136.4	40	-0.53	-1.50	1.05	323.48	36	0.00	-0.75	1.12	366.70
	R101	100	2	1	200	40	-0.55	-1.28	0.29	557.74	36	-0.08	-1.24	1.83	614.35
	R101	100	3	0.6	188.4	40	-0.53	-1.23	0.02	548.02	39	0.21	-0.94	1.61	576.27
	R101	100	3	1	300	40	-0.62	-1.72	0.25	1078.70	36	0.10	-1.78	2.67	1211.77
			total			320	-0.70	-2.88	1.05	633.85	300	-0.14	-2.47	2.67	696.03	187

Table 11 :

 11 Characteristics of split customers in the best solutions of mid-80 instances (5000 iterations).

	instances
		n	p	id	nbSplit nb2-split nb3-split nbN earDepot nbLargeDemand	nbCluster
	C101	80	0.6	1	7	7	0	2	3	7
				2	5	5	0	3	1	5
				3	7	7	0	4	2	5
				4	9	8	1	3	5	7
				5	10	10	0	1	7	6
			1	1	20	19	1	5	5	18
				2	14	14	0	4	4	14
				3	18	18	0	6	4	13
				4	14	14	0	4	4	14
				5	18	18	0	6	2	16
	R101	80	0.6	1	6	6	0	4	6	4
				2	5	5	0	3	2	2
				3	7	7	0	2	3	4
				4	5	5	0	3	1	3
				5	6	6	0	6	3	3
			1	1	13	13	0	3	3	7
				2	14	13	1	7	1	4
				3	13	13	0	4	2	6
				4	14	14	0	4	4	6
				5	17	17	0	6	5	7

Table 12 :

 12 Detailed computational results for the small sized instances (n = 15).

	instances

Table 13 :

 13 Detailed computational results for the mid sized instances (n ∈ {20; 40; 60; 80}).

	instances

Table 13 (

 13 continued).

	instances

Table 14 :

 14 Detailed computational results for the large sized instances (n = 100).

	instances

Table 14 (

 14 continued).

	instances

Table 14 (

 14 continued).

	instances

Table 14 (

 14 continued).

	instances

Table 14 (

 14 continued).

	instances

Table 14 (

 14 continued).

	instances

Acknowledgments

This work is partially supported by the CSC (China Scholarship Council) and by the ELSAT 2020 project. The authors thank C. Archetti and N. Bianchessi for providing the benchmark instances. Thanks are also due to the referees for their valuable comments.

A Detailed results on the benchmark instances

Table 12, Table 13 and Table 14 report the detailed results for the small, medium and large instances respectively. We report values in bold whenever we improve (or optimize or equal to) the respective instances.


Wilfried Raballand

Pierre-Antoine Laharotte

SIMPLIFIED MODELLING OF TRAVEL IN LOW DENSITY AREAS

INTRODUCTION AND OBJECTIVES

At a time when institutions are involved with the issue of energy transition, communities with low density areas must equip themselves with tools to deal with the issue of travel needs and transport supply. In addition to environmental considerations, [START_REF] Schutz | Transportation Planning Needs for Small and Medium-Sized Communities[END_REF] has identified a wide range of transportation planning needs in small communities, including safety, development, diversification and changing modes of transport. Owing to lack of sufficient financial resources, these areas cannot collect the expensive data needed to develop a travel model according to the usual tested methods. It is therefore necessary to develop tools to propose "simplified" modelling. Such basic models have been developed since the 1970s. [START_REF] Ortúzar | Modelling Transport. 4th Edition[END_REF] identify the usual methods that can be applied in such cases. More specifically, [START_REF] Bera | Estimation of origin-destination matrix from traffic counts: the state of the art[END_REF] have established a state of the art for models that use counting data.

As part of the study of a new road infrastructure begun in 2017 in the French town of Vesoul (35,000 inhabitants), a significant amount of data was collected, consisting of 70 automatic vehicle counting points, and 36 number plate reading (NPR) stations. This set of data enables to build an Origin-Destination (OD) matrix, which is assigned to the road network of the area. This "conventional model" based on observed flow data is used as a ground truth resource in the introduced results. Strong challenges related to the new infrastructure made this financial investment possible in terms of data and studies. The present work uses the Vesoul case study to develop and test a simplification process to model the road trips in the area. The objective of the study lies in reducing the amount of data and assessing the performance impact. The models developed require a smaller amount of expensive data. This means that it is possible to propose a soberer, more systematic modelling approach to the local authority of the low density areas sheltering between 20000 and 50000 inhabitants. Furthermore this approach enables to spare time compared to usual studies, which is an additional saving on the cost of the assessment tool.

METHOD

From the conventional model developed during the study begun in 2017, we keep the geographical features of the transport supply and demand:

• the road network on which the demand is assigned,

• the division of the area into zones which are either the destination or the starting point for trips. We propose a simplified model that uses the usual 4-step scheme: generation, distribution, modal choice, assignment. The study makes a focus on demand estimation, and on the category of light vehicles. The modal choice is considered as near-homogeneous for all OD relationships. This is because in this type of urban area, private car usage is largely predominant among all modes of transport, due to inefficient public transport, low traffic congestion and ease of parking. In addition, in order to focus on the first steps of the model, we will use the Vesoul road network which was already tuned during the path allocation stage, as well as the same commercial software to assess the simplified model. In order to reduce the use of expensive data, we have included in the model the free, public statistical data published by the French National Institute of Statistics and Economic Studies (INSEE), including population and employment data from the general census. These data help to explain people's demand for mobility. They make it possible to construct regression models and assumptions concerning the generation and the attractiveness of trips for each zone of the area, distinguishing between morning rush hour (MRH) and evening rush hour (ERH).

A 3 stages process to reproduce OD matrices The simplified model is built in three stages. Each corresponds to the treatment of one of the three types of trip that make up the demand: internal flows, exchange flows, and transit flows. At each stage, the work involves testing and validating the assumptions of the simplified generation-distribution model, using field data collected as part of the initial study. For the two stages concerning internal flows and exchange flows, we follow the usual method, which involves aggregating and standardizing the mobility behaviour of the users.

Stage 1: reproducing the internal flows

The first step, generation, is based on a simple expression of the emission of trips Ei as a linear combination of population Popi and employment Empi data from the internal area of origin i (see Eq. 1). This principle is validated by means of a linear regression. The second step, distribution, is based partly on a trip attractiveness function, and partly on a function of the location and/or distance of the destination zone. We begin by defining attractiveness Aj as a linear combination of population and employment data from the internal area of destination j (see Eq. 2). Next, we analyse the residual function fI of Eq. 3 to estimate the plausible impact of the impedance zij between zones i and j, on the distributed volume of trips Dij, as in the case of gravity models. 𝑓 𝐼 (𝑧 𝑖𝑗 ) is defined empirically, and K is the constant of proportionality of a doubly-constrained gravity model.

(1)

𝐸 𝑖 = 𝛾 𝑒 𝑃𝑜𝑝 ⁄ 𝑃𝑜𝑝 𝑖 + 𝛾 𝑒 𝐸𝑚𝑝 ⁄ 𝐸𝑚𝑝 𝑖 (2) 𝐴 𝑗 = 𝛾 𝑎 𝑃𝑜𝑝 ⁄ 𝑃𝑜𝑝 𝑗 + 𝛾 𝑎 𝐸𝑚𝑝 ⁄ 𝐸𝑚𝑝 𝑗 (3) 𝐷 𝑖→𝑗 = 𝐾 𝐸 𝑖 𝐴 𝑗 𝑓 𝐼 (𝑧 𝑖𝑗 )
Stage 2: reproducing the transit flows For the third stage concerning transit flows, the process is different, because it is not possible to solve generation and attractiveness on the fictitious external zones simply and jointly. In order to replace the NPR data, it is necessary to define a model specific to transit that can be based on automatic counts only. We propose three assumptions (see Eq. 4,5 and 6), making it possible, respectively, to estimate:

• the isotropic fictive relationship between incoming exchange flows and outgoing exchange flows at rush hours,

(4) ∑ 𝐷 𝑟 𝑘→ 𝑗 𝑖𝑛𝑐𝑜𝑚𝑖𝑛𝑔 𝑒𝑥𝑐ℎ𝑎𝑛𝑔𝑒 𝑗 = 𝜏 𝑀𝑅𝐻 (𝑜𝑟 𝐸𝑅𝐻) ∑ 𝐷 𝑖→𝑟 𝑘 𝑜𝑢𝑡𝑔𝑜𝑖𝑛𝑔 𝑒𝑥𝑐ℎ𝑎𝑛𝑔𝑒 𝑖 , ∀𝑘 ∈ [1, 𝑛]
For all n entry and exit points (as associated to the external fictive areas rk, 𝑘 ∈ [1, 𝑛]), we have a relation of proportionality between the number of incoming exchange trips and the number of outgoing exchange trips, with a given coefficient  that depends of the rush hour.  can be expressed as:

𝜏 𝑀𝑅𝐻 (𝑜𝑟 𝐸𝑅𝐻) = min 𝑘∈[1,𝑛] { ∑ 𝐷 𝑟 𝑘→ 𝑗 𝑗 ∑ 𝐷 𝑖→𝑟 𝑘 𝑖 ⁄ }
where ∑ 𝐷 𝑟 𝑘→ 𝑗 𝑗 is equal to the automatic count of incoming vehicles, and ∑ 𝐷 𝑖→𝑟 𝑘 𝑖 is equal to the automatic count of outgoing vehicles at the entry and exit point k.

• the total volume of transit flows (passing through the urban area), is defined by the combination of the residual flows resulting from the difference between total flows and the fictive exchange flows defined in Eq. 4:

(5)

∑ 𝐷 𝑟 1 →𝑟 2 𝑡𝑟𝑎𝑛𝑠𝑖𝑡 𝑟 1 ,𝑟 2 = ∑ | ∑ 𝐷 𝑖→𝑟 𝑘 𝑖 - ∑ 𝐷 𝑟 𝑘 →𝑗 𝑗 | ,  = { 𝜏 𝑀𝑅𝐻 , 𝑎𝑡 𝑀𝑅𝐻 𝜏 𝐸𝑅𝐻 -1 , 𝑎𝑡 𝐸𝑅𝐻 𝑘∈[1,𝑛]
The volume of transit flows is expressed as the anisotropic result of the difference between the total incoming flows and the incoming exchange flows. • the distribution of transit flows between the various entry and exit points in the area.

(6)

𝐷 𝑟 1 →𝑟 2 𝑡𝑟𝑎𝑛𝑠𝑖𝑡 = 𝑓 𝑡𝑟𝑎𝑛𝑠𝑖𝑡 () ∑ 𝐷 𝑖→𝑟 2 𝑖 ∑ 𝐷 𝑟 1 →𝑗 𝑗 If  (also noted 𝑟 1 𝑂𝑟 2 ̂)
is the angle formed by the entry point r1, the center of the urban area O, and the exit point r2, ftransit() is a strictly increasing function from =0 to =.

Stage 3: reproducing the exchange flow

The exchange flow is defined in 2 steps:

• assessment of the total incoming (respectively outcoming) exchange flow for any crossing point 𝑟 𝑘 : this volume results from the difference between the counting volume and the combination of transit flows incoming (resp. outcoming) by 𝑟 𝑘 : ∑ 𝐷 𝑟 𝑘→ 𝑗 𝑗 -∑ 𝐷 𝑟 𝑘→ 𝑟 𝑖 𝑖 ;

• distribution of the total incoming / outcoming exchange flow on the OD relationship according to the emission and attractiveness of the internal areas.

At the end of each stage, we obtain a new travel demand. This new demand is then substituted within the global matrix that has been built from the field data. The resulting matrix is assigned with the same tool and the same settings used during the initial study.

The results of the simplified model in terms of traffic on the road network are evaluated on the basis of a comparison with the results of the conventional model, notably concerning:

• the reproduction of actual traffic observed on the entire road network,

• the estimation of simulated traffic shifts onto the new infrastructure project. In addition, an inventory of the expensive data used and the corresponding budget required is drawn up in each phase. This assessment serves to quantify the savings made in technical and financial means compared to the initial data collection. It also illustrates the correlation between the amount of field data used and the quality of modelling the actual situation.

RESULTS

In the following, the results of the simplified model are compared either with data collected directly in the field, or with the demand tuned during use in the 2017 study.

Emission and attractiveness assessment

In the case of simplified modelling of internal flows and exchange flows, we apply the usual method to solve Eqs. 1 and 2. The variables explaining generation and attractiveness by zone that give the best results are as follows:

• the population for generation, and the sum (population + jobs) for attractiveness for the MRH, • the sum (population + jobs) for generation, and also for attractiveness for the ERH. The generation and attractiveness models obtained are validated by linear regression (see Table 1). A satisfactory result overall is to be noted, with a lower quality for exchange traffic, and at the MRH. Concerning any possible impact of the impedance between zone of origin and zone of destination on distribution, we propose testing a phenomenon of resistance to the distance travelled (see Eq. 7). Fig. 2 illustrates the resolution of coefficient β by regression. The results show that there is no significant resistance to distance, such that β is strictly negative. After an alternative analysis of the resistance to travel time, we obtain a comparable result.

(7) 𝑓 𝐼 (𝑧 𝑖𝑗 ) = 𝑑 𝑖𝑗 𝛽 Therefore, as we are unable to resolve the dependence of distribution on remoteness, we propose to assume that fI(zij)=1.

Evaluation of the assumptions about the transit flows

In the case of transit flows, we first examine the assumption that the ratio between incoming and outgoing exchange flows is isotropic (Eq. 4) for a given rush hour.

Regression analysis gives us R²=0.94 for the MRH, and R²=0.96 for the ERH, and so validates this assumption.

Using the relation (Eq. 5), we estimate the total volume of transit flows on the basis of the anisotropic part of the traffic entering and leaving the area. Comparison between the survey data and the results of our model (Table 2) shows that the assumptions about volume estimation are acceptable. The distribution of the transit movements between entries into and exits from the urban area is calculated using Eq. 6 and the following definition (Eq. 8) of the function ftransit(ψ). A regression analysis of the real transit flows gives a satisfactory validation of this distribution function as it is defined: R²=0.76 for the MRH, and R²=0.72 for the ERH.

(8)

𝑓 𝑡𝑟𝑎𝑛𝑠𝑖𝑡 () = | 1 𝑖𝑓  > 𝜋 2 ⁄ √2 2 ⁄ 𝑖𝑓  = 𝜋 2 ⁄ 𝜖 ≪ 1 𝑖𝑓  < 𝜋 2 ⁄

Global comparison

For each development phase of the simplified model we can assess the performance of the simplified model by comparing the results obtained with the data collected in the field. We consider as comparison criterion:

• aggregate indicators describing the correlation between observed traffic and modelled traffic: R², and more specifically for this type of data, the Geoffrey E. Havers (GEH) indicator (more accurately the value of the third quartile of all GEHs). • the estimation of the volume of traffic during a working day on the new road infrastructure. In order to more clearly situate one model in relation to the others, we plot any criteria according to the overall cost of the expensive data needed to develop it. These overall costs are shown in Table 3, and they are the X axis of the following graphs. Fig. 3 gives the changes in the comparison criteria chosen according to the cost of the data used.

In these representations, moving from right to left, we reduce the use of expensive data and increase the use of the simplified model. We observe:

• for R²: a moderate but homogeneous decline in data correlation,

• for the GEH: a significant, heterogeneous increase according to the rush hour. Simplified modelling significantly affects the individual reproduction of traffic volumes as a whole. However, the divergence found is acceptable and does not presuppose by itself the capabilities of the tool developed. In particular, we can target the new infrastructure and analyse traffic estimates calculated over a working day (see Fig. 4). To within the margin of error, the estimated volume on the infrastructure does not vary significantly and retains the same order of magnitude. We can therefore conclude that the simplified model provides a comparable assessment of the project on this criterion. We can also provide a detailed inventory of the data needed for each modelling phase in terms of the quantity and types of collections used. Table 4 gives the comparison for all the models tested mentioned above. In addition to the savings achieved, the field survey is observed to be simplified and the system required less cumbersome.

CONCLUSION

The results show that modelling of travel in this type of area is not totally dependent on field data collections, and therefore on the budget available. In particular, there is a tested and validated alternative to demand estimation, which merges traffic data with public statistical data. This work provides a development basis for modelling solutions in less urbanized areas. The Vesoul case study explored the potential and the limits of such a "simplified" modelling approach. Despite significant differences with results from a conventional model, the simplified model is a robust tool, simple to develop, and well suited for assistance with decision-making in low density areas. In particular, it can be applied in the study of the development of green and shared mobility in rural areas.

Fig. 1

 1 Fig.1represents, on the left, a map of the urban area, and, on the right, the layout of the road network and the zoning of the model. We propose a simplified model that uses the usual 4-step scheme: generation, distribution, modal choice, assignment. The study makes a focus on demand estimation, and on the category of light vehicles. The modal choice is considered as near-homogeneous for all OD relationships. This is because in this type of urban area, private car usage is largely predominant among all modes of transport, due to inefficient public transport, low traffic congestion and ease of parking. In addition, in order to focus on the first steps of the model, we will use the Vesoul road network which was already tuned during the path allocation stage, as well as the same commercial software to assess the simplified model.

Fig. 1 :

 1 Fig. 1: Map of the urban area of Vesoul (left). Network and zoning of the Vesoul model (right).

Fig. 2 :

 2 Fig. 2: Determination of the coefficient β by linear regression on internal flows. Estimated value of β at MRH

Fig. 3 :

 3 Fig. 3: Aggregate indicators of discrepancies between observed traffic and modelled traffic.R² at MRH in blue, at ERH in red

Table 1 :

 1 Validation by regression (coefficient R²) of the generation and attractiveness laws

		Ei -Internal flows	Aj -Internal flows	Ei -Exchange flows	Aj -Exchange flows
	MRH	0.90	0.91	0.64	0.73
	ERH	0.91	0.90	0.71	0.84

Table 2 :

 2 Validation of the estimation of the volume of transit flows

		Field survey	Simplified model	Absolute difference	Relative difference
	MRH	477	526	+49	+10 %
	ERH	543	549	+6	+1 %

Table 3 :

 3 Overall cost in expensive data of the different models

		Conventional model	Simplified internal	Simplified internal	Overall simplified
			mode	and exchange model	model
	Amount	€41,400	€30,600	€12,300	€5,100
	Proportion simp/conv	1	0.74	0.30	0.12

Table 4 :

 4 Number and types of data used according to the model examined

	Type of data	Conventional model	Simplified internal	Simplified internal	Overall simplified
		"Ground Truth Data"	mode	and exchange model	model
	Vehicle counts within	58	22	<10	<10
	the perimeter
	Cordon survey	12	12	12	12
	vehicle counts
	NPR internally	22	22	-	-
	Cordon survey NPR	12	12	12	-
	INSEE population	-	X	X	X
	and employment