MODELING AND SOLVING A RICH VEHICLE ROUTING PROBLEM FOR THE DELIVERY OF GOODS IN URBAN AREAS

This work addresses a vehicle routing problem that aims at representing delivery operations of large volumes of products in dense urban areas. Inspired by a case study in a drinks producer and distributor, we propose a mathematical programming model and solution approaches that take into account costs with own and chartered vehicles, multiple deliverymen, time windows in customers, compatibility of vehicles and customers, time limitations for the circulation of large vehicles in city centers and multiple daily trips. Results with instances based on real data provided by the company highlight the potential of applicability of some of the proposed methods.


INTRODUCTION
Vehicle routing problems arise in many practical situations, including the delivery of goods to customers, the pick-up and the transportation of urban waste to regulated landfills, the elaboration of itineraries for electrical vehicles with stops in recharging stations, the transportation of people with reduced mobility, and equipment repair in different houses.The literature presents several surveys on problem variants, solution methods and real applications, such as Bodin et al. (1983), Assad (1988), Ronen (1988), Osman (1993), Desroisiers et al. (1995), Cunha (2000), Breedam (2001), Bräysy & Gendreau (2005a, 2005b), Parragh et al. (2008), Laporte (2009), Baldacci et al. (2010), Belfiore & Yoshizaki (2013) and Schneider et al. (2014).The large number of contributions shows that a great deal of effort has been made for over 50 years by the scientific community when tackling these problems.This is due not only to their practical relevance but also to the difficulties faced when solving them.Nevertheless, as each studied case can present distinct features from what has already been discussed in the literature, new variations of the general problem continue to appear and require differentiated treatments.*Corresponding author.Departamento de Engenharia de Produc ¸ão, Universidade Federal de São Carlos -UFSCar, Rodovia Washington Luís, km 235, SP-310, 13565-905 São Carlos, SP, Brazil.E-mails: josefneto@yahoo.com.br;vpureza@dep.ufscar.br In this paper, we consider a variant that represents the situation of many companies that deliver large quantities of goods in highly dense urban areas in a daily basis.Given the difficulties in driving and parking, customers close to each other are clustered, and for each cluster, a single parking slot is used by the vehicle assigned to serve those customers.The deliveries are then performed on foot by the vehicles' crew from the parking slot.The relatively long service times in each cluster when compared to the vehicles traveling times justify the use of multiple deliverymen as it allows the reduction of the time required by a single deliveryman to serve a cluster and consequently increases the number of customers visited during working hours.Avoiding the violation of working hours is a primary objective of such companies given its direct impact on overtime costs.
Situations such as the one described are common in Brazilian drinks and tobacco companies, for which a large number of customers consists of small and medium sized retailers located in urban areas.In addition to the traditional routing and scheduling decisions, the route plan must elect the crew size in each vehicle as service times depend on the number of deliverymen employed.Time windows for the deliveries are often present since they avoid conflicts in critical hours (for example, lunch time in restaurants) or with the deliveries from other suppliers.In addition, due to large demands or constraints that impede some customers from being visited by large vehicles, it is often necessary that one or more vehicles make a second tour (trip) from the depot (warehouse) to serve additional customers.The problem can therefore be referred as a Multi-trip Vehicle Routing Problem with Time Windows and Multiple Deliverymen (MTVRPTWMD).
The concept of routes with multiple trips for each vehicle was first addressed for the VRP in the working paper by Fleischmann (1990).The author propose an extension of the Clark & Wright (1964) savings heuristic to generate trips and use a heuristic for the bin packing problem to assign the trips to a homogeneous fleet.As pointed out in Petch & Salhi (2004), enabling multiple trips can be essential for tactical and strategic planning, allowing savings in all transportation costs.Since the introduction of the concept, several papers explore heuristic solution approaches such as tabu search (Taillard et al., 1996;Brandão & Mercer, 1997;Olivera & Viera, 2007;Alonso et al., 2008;Seixas & Mendes, 2013), genetic algorithms (Salhi & Petch, 2007), and exact algorithms (Azi et al., 2010).
Routing with multiple deliverymen, in turn, was introduced in Ferreira & Pureza (2012) with the variant Vehicle Routing Problem with Multiple Deliverymen (VRPMD) and in Pureza et al. (2012) with the Vehicle Routing Problem with Time Windows and Multiple Deliverymen (VRPTWMD).Ferreira & Pureza present a mathematical model for the case of a limited fleet of identical vehicles, which is then solved by an extension of the Clark and Wright's savings heuristic and a tabu search algorithm.Pureza et al. formulate the problem with an unlimited fleet of vehicles and propose an ant colony optimization algorithm and a tabu search algorithm for solving the model.In both papers, and also in Grancy & Reimann (2014), Alvarez Diaz & Munari (2016) and Munari & Morabito (2016), the mathematical models and solution approaches were motivated by the practical relevance of the issue; however, no real-life cases are discussed or solved.
Our research, on the other hand, is inspired by the urban delivery operations of a drinks manufacturer and distributor located in the state of São Paulo, Brazil.The delivery is carried out by multiple deliverymen in each vehicle (truck) and the logistic operation considers different characteristics and additional constraints from the models of the previous works, such as a limited heterogeneous fleet, the use of chartered trucks, multiple trips for the same truck, varying number of deliverymen in each trip, the existence of dangerous routes, limitations on circulation times of truck types in city centers and incompatibility of truck types and customers.
Delivery operations in urban areas are often difficult to implement due to several real-life constraints.The consideration of multiple deliverymen and multiple trips adds more complexity to the problem, therefore, well-structured procedures capable to provide efficient solutions is of great interest to distributors.As far as we know, the decision on the number of deliverymen has not been explored in current computational systems designed to support vehicle routing and scheduling decisions.Nevertheless, its practical and theoretical relevance is highlighted in the taxonomy of routing problems, recently proposed in Braekers et al. (2015).
In addition to this motivation, in the case of the company studied, the number of deliverymen in each truck is fixed.Therefore, it is interesting to verify whether setting the number of deliverymen based on the specificities of each trip can bring benefits, such as cost reduction, satisfaction of constraints formerly violated, and the increase of the number of demand clusters visited within the deliverymen working hours.The last advantage is particularly relevant and directly related to the first; in our study case, we observe that even with multiple deliverymen and multiple trips, meeting the daily demand usually require the use of overtime.As the company is typical of the drinks sector, the model and solution approaches proposed can be useful for other suppliers, as well as for companies from other sectors with similar distribution operations.Besides the practical motivation, vehicle routing problems with multiple deliverymen remain little explored, which means that its continuing study can represent contributions to the body of knowledge in combinatorial optimization and distribution logistics.
The remainder of this paper is organized as follows.Section 2 presents a brief description of the company and the relevant aspects of its distribution operations.The mathematical model is presented in Section 3, while Section 4 describes a GRASP approach for solving instances based on real data provided by the company.Section 5 presents the results and analysis of computational experiments with the model, the heuristic and a hybrid approach, followed by conclusions and perspectives of future research in Section 6.

THE COMPANY'S DISTRIBUTION LOGISTICS
The drinks company addressed in this research produces and delivers a large variety of bottled products such as different types of water, juices, beer, tea, soft, energy and sports drinks to approximately 25,000 customers per week.The customers are distributed in 131 cities of the states of São Paulo and Minas Gerais and consist of retailers of different market segments and demand volumes, ranging from snack bars to hypermarkets.
In 2014, the company's fleet accounted a total of 222 trucks for product delivery.The fleet encompasses six types of trucks, characterized by different numbers of doors, weight and space capacities and functionalities.In particular, some truck types have mobile platforms that allow lifting the cargo to a higher ground.In addition to the own fleet, the company hires trucks for periods of the year during which the requirements of the demand exceeds the transportation capacity (e.g., school vacations and national holidays).
Nevertheless, even with the use of additional trucks, a second trip from the warehouse may be necessary for each available truck, not only to meet high demands but also as the result of regulations that restrain the use of certain types of trucks.One example is a municipal law, valid in most large Brazilian cities, that allows the circulation of large trucks in city centers only in early working hours (usually up to 9:00 or 10:00 a.m.).The company also prohibits the use of some of its largest trucks in city centers due to their riding difficulties when the traffic is intense.Other constraints include the compatibility between the truck type and the customer.For instance, most supermarkets and hypermarkets require trucks equipped with mobile platforms while retailers whose location can only be accessed through narrow streets must be served by smaller trucks.However, in order to perform a second trip, the truck capacity utilization must be equal or greater than a minimum percentage (usually 83%).
In addition to the driver, the standard crew in each truck comprises a second man and they both deliver the goods to the customers.However, if the cargo is addressed to areas controlled by drug gangs (called "dangerous routes" by the company), a security guard is added to the standard crew.It should be noted that the guard does not act as a deliveryman; his/her role is solely to ensure the security of the crew, the cargo and the truck.For the two-man crew, service times in each customer are estimated by the company as 5 minutes plus 20 seconds per product pack (in case of hypermarkets or supermarkets) or 30 seconds per product pack (other market segments).
The trucks routes are elaborated by the logistic operators with the aid of a commercial software.The routes follow a geographical sectorization of the customers, which simplifies the planning process.Nevertheless, the crews often review the prescribed itineraries for several reasons.First, they have more updated information on the traffic and street conditions.Second, it is sometimes more efficient to serve a customer from a given sector in a route of an adjacent sector if the customer is located close to the border between these sectors.Third, given the difficulty in parking the trucks due to scarcity of parking slots in city centers, the truck does not necessarily travel from one customer to the next; instead, if there are customers sufficiently close to each other, these customers are seen as a demand cluster, and a single parking slot within a maximum distance from each customer (approximately 150 meters) is elected for the truck serving the cluster.The goods are then delivered at each customer by the two-man crew, who visits the customers on foot (Fig. 1).This means that in practice, the distribution can be seen as a two-echelon routing problem: in the first level, the trucks travel from the warehouse to the parking slot elected to the cluster, and in the second level the truck's crew delivers the goods to the customers' locations on foot.The crews also change the routes sequencing by taking into account the customers' preferences of delivery time and location, which are not included in the data base available to the routing software.In fact, the only recorded scheduling constraints on customers are time windows in restaurants.For this market segment, deliveries must occur at least two hours prior the opening hours (mostly 9:00 a.m.) in order to cool the products.At the time of the study case, the major challenge faced by the company regarding their delivery operations was cost reduction.Since the Brazilian drinks market is highly competitive, the company's policy is to meet all daily requests on the same day, which often exceeds the deliverymen working hours.In these cases, instead of overtime payment, the company offers compensatory time off.However, if the banked time is not taken by the end of the current fiscal year, the employees must be paid for the overtime hours in the beginning of the following year.In order to avoid these costs, the company wants to eliminate overtime while maintaining the service level provided to the customers.

PROBLEM DESCRIPTION AND MODELING
Formally, the problem tackled in this paper consists of planning daily routes for a fleet of heterogeneous vehicles in order to deliver goods to customers in urban areas.The fleet departs (and returns) to the company's warehouse and encompasses own and chartered vehicles.For each of these categories, distinct costs are applied; fuel and maintenance (variable) costs incur in own vehicles while the affreightment contract between the company and the autonomous truck owners prescribe only hiring (fixed) costs to the chartered vehicles.
A route here denotes the itinerary followed by a vehicle between its departure from the warehouse at the beginning of the day and its return to the warehouse when all of its assigned deliveries are completed.A route can comprise two tours (trips) at most, as long as the total route time does not violate the working hours decremented by the time spent on administrative tasks and the deliverymen lunch time (i.e., no overtime is allowed).If a second trip is performed, the vehicle is reloaded at the warehouse.
Each node in a given route represents the warehouse or a demand node.A demand node is a single customer or a cluster of customers, and it is represented by the associated parking site.We assume that clusters were previously defined (either manually or by a clustering algorithm) in order to ensure a maximum radial distance between each costumer and the parking site, and that parking sites are in the exact location of the single customer or one of the customers in case of clusters.
In each trip, the total cargo must not exceed the vehicle capacity, and the crew size is limited to 3, the capacity of the vehicle's cabin.The number of deliverymen in a vehicle can vary from the first to the second trip, since changes in the workload may require a different number of deliverymen.In case of dangerous routes, the maximum number of deliverymen reduces to 2 in order to accommodate the security guard.Customers served in a given route must be compatible to the type of vehicle.In addition, large vehicles are limited to ride in city centers in determined periods of day.
Given the aforementioned constraints, the route plan should otimize three objectives according to the following order of priority: (1) maximize the number of demand nodes served; (2) minimize the cost with own and chartered vehicles; and (3) minimize the number of deliverymen employed.Consider the following notation: • Indexes i, j, h, m Demand nodes (represented by associated parking sites) and the warehouse copies.Pmin Minimum capacity utilization of a vehicle in the second trip (%); • Variables For i ∈ D, it is also the arrival time in i.
The problem is formulated as a linear mixed integer program, as follows.We also consider that some variables are previously fixed.For example, x i jklr = {0, 1}, The model's objective ( 1) is to minimize costs with own and chartered vehicles as well as the number of deliverymen decremented by the aggregate "priority" of the costumers associated to each visited parking site.The values of parameters P j , Co, Cs and G are set to guarantee the lexicographic order <number of served clusters, own and chartered vehicles cost, number of deliverymen>.In our experiments, all costumers have the same priority, therefore P j is replaced by P.
Constraints (2) impose that at most one vehicle k in mode l in a trip r arrives at a demand node j from other demand node i or from the origin node of r, while constraints (3) force that at most one vehicle k in mode l in a trip r departs from a demand node i to another demand node j or to the destination node of the trip.Constraints (4) guarantee that the same vehicle k that enters a demand node h in mode l leaves h to another demand node or the node destination of the trip in the same mode l.
Constraints (5) guarantee that each vehicle k in trip r arrives at the destination node of r in mode l from a single node i (demand type or origin node of trip r).Constraints (6) ensure that each vehicle k in trip r leaves the origin node of r in mode l to a single node j (demand type or the destination node of r).Constraints (7) ensure that each vehicle k in the first trip arrives at its destination node in a mode g and departs in the second trip to a node j (demand type or destination node of the second trip) in a mode l.
Constraints (8) impose that the load onboard each vehicle k in the origin node of each trip r as the total demand of the nodes served by k in r.Constraints (9) compute the load onboard vehicle k in a mode l in trip r when it arrives at node j after visiting node i. Constraints (10) define the starting time of service in each node j visited after node i in route r with vehicle k in mode l.
Constraints (11) prescribe that the total route time must not exceed the maximum time that guarantees the return to the warehouse within the deliverymen working hours.Constraints (12) ensure that vehicles with circulation limitations do not ride in the city center after the maximum circulation time F.This is done by imposing that if the vehicle departs from a node i located in the city center to visit a node j outside that area (among which, the warehouse) then the vehicle must cross the border of the city center at time F at most.Therefore, for each pair of nodes (i, j ) (i ∈ B; j / ∈ B) the most probable border crossing is identified in a pre-processing phase by, for example, computing the shortest path between i and j .As the distance between i and the border is a function of the destination j , it can be refereed as DC i j , i.e., in terms of i and j .
Constraints (13) express the relationship between the arrival time of vehicle k in mode l at the destination node of the first trip and the departure time of k in mode g from the origin node of the second trip.As the origin node of the second trip is equal to the destination node of the first trip, the constraints link these time instants considering that a change of mode can occur.The departure time of k from the origin node of the second trip (in case it occurs) is equal to the arrival time of k at the destination node of the first trip added by the vehicle reloading time.
Constraints (14) guarantee the precedence between the first and second trips, eliminating symmetric solutions.Constraints (15) ensure the minimum capacity utilization of the vehicles in the second trip.Finally, constraints (16) to (18) impose the decision variables domain.It should be noted that the model can be easily extended to cope with situations with more than two trips.

A HEURISTIC APPROACH
This section describes a GRASP algorithm (Feo & Resende, 1995) for solving model ( 1)- (18), presented in the previous section.GRASP is an iterative procedure for which each iteration consists of the construction of an initial solution followed by an improvement stage.The initial solution is constructed by randomly selecting solution components from a restricted list of candidates, which seeks to ensure a good quality selection.The best generated (and improved) solution obtained along the iterations is, therefore, the solution returned by the method.
The general steps of our algorithm are presented in Figure 3, and a detailed description of each step is shown in the following sections.Note that the construction stage (step 2.1) is applied with set FD, corresponding to the original fleet (own and chartered vehicles) duplicated.The fleet duplication allows us to easily address the first and the second trip of a given vehicle.Specifically, each vehicle k has an exact copy k , so that in case k is used in the solution and a second trip is required, the first trip is performed by k while k serves the customers of the second trip.The coupling of the two trips is ensured by imposing that the departure time of k from the warehouse is equal to the arrival of k to the warehouse plus the vehicle loading time.

Construction stage:
For the set of demand nodes C and the set of duplicated vehicles FD, obtain a starting solution with the insertion heuristic described in Section 4.1.If the solution comprises one or more second trips with violated minimum load constraints, apply the feasibility procedure described in Section 4.2.

Improvement stage:
Apply the trip reduction procedure to the starting solution, followed by the travelled distance reduction procedure, followed by the deliverymen reduction procedure (Section 4.3).Let S be the resulting solution.If f (S) ≤ f (S * ), make S * = S.

3.
Return S * .The construction stage comprises the insertion heuristic described in Section 4.1, which guarantees feasible solutions except for the minimum load constraints for the second trip.In case these constraints are violated, the heuristic uses a feasibility procedure (described in Section 4.2) to transfer demand nodes from feasible trips to infeasible second trips.If there are still infeasible second trips, these are deconstructed and the demand nodes are reconsidered for insertion.
The algorithm improvement stage (step 2.2), in turn, consists of the three-phase local search described in Section 4.3.The first phase aims at reducing the number of trips, followed by the reduction of the travelled distance and finally, by the decrease of the number of deliverymen.
Note that the order of the phases reflects the order of importance of the model's objectives: by reducing the number of trips, the number of routes may also decrease, making the resulting idle vehicles available for serving yet unrouted customers and perhaps eliminating the need for chartered vehicles; by reducing the travelled distance, the cost with own vehicles lessens, while reducing the number of deliverymen completes the effort to make the solution more efficient.

The insertion heuristic
Starting solutions are generated as described in Figure 4. Due to the problem highly constrained nature, the limited resources (vehicles) must be well allocated, therefore, in step 2 the set of the duplicated vehicles FD is partitioned into subsets (i) own vehicles, (ii) own vehicles copies, (iii) chartered vehicles, and (iv) chartered vehicles copies.In each of these subsets, vehicles are sorted according to the following criteria: (1) smaller number of compatible customers, (2) subjected to circulation limitations, and (3) smaller capacity, so that criterion i decides the sorting if criterion i − 1 declares the vehicles tied.A list LV is then constructed by sorting the subsets from (i) to (iv), which reflects the vehicles priority for selection when a new trip is initialized in step 3.2.1.Note that this ranking favors the selection of lower cost, more constrained vehicles, and complies with the trips sequencing as well.
Also due to the constrained nature of the problem, the restricted list of candidates LRC prioritizes demand nodes with smaller chances of insertion.First, the list exclusively comprises all nodes with time windows (step 3.1), which are randomly selected, one by one, and inserted in the least costly feasible position among all current trips (step 3.2.1).In case there is no feasible insertion, the selected demand node p starts a new trip with the first compatible vehicle k from the sorted list LV.The largest possible crew is assigned to the trip, i.e., two deliverymen if p is located in a drug gang area or three deliverymen, otherwise.The departure of k from the warehouse occurs at time 0 if k is a first trip vehicle; otherwise, its departure is equal to the return time to the warehouse of the associated first trip vehicle added by the reloading time.
Nodes that cannot be inserted in any partial route or start a new trip are removed from LRC and stored in set NCI (step 3.2.1).After concluding the insertion and scheduling of nodes with time windows, the solution construction resumes with LRC redefined with α% of the best evaluated nodes (without time windows) according to the following criteria: (4) smaller number of compatible vehicles, and (5) larger demand.As before, criterion i decides the sorting if criterion i − 1 declares the nodes tied.

1.
Let NCR = set of unserved demand nodes = C, FA = set of idle vehicles = FD, and S = ∅.

2.
Let NCI = set of unserved demand nodes whose insertion failed = ∅, CT = set of unserved demand nodes with time windows, and new = false.Build the sorted list LV of available vehicles.

3.
For j = 1 to 2: Select the first node p from list LRC that can be served in one of the existent trips, and its feasible insertion position (disregarded the minimum load constraints) that result in the smallest cost increase.If there is no such p, make new = true and go to step 3.2; otherwise, insert p, update the vehicle capacity and reschedule the trip.

3.2.2
Let S the current partial solution.Make NCR = NCR − {p} and return to step 3.1.

The feasibility procedure
In case the solution produced by the insertion heuristic (Fig. 4) violates the minimum load constraints in one or more second trips, a feasibility procedure is called to relocate demand nodes from feasible second trips, served by other vehicles, to the infeasible second trips and, in case the latter do not become feasible, to relocate demand nodes from first trips.Figure 5 presents the procedure steps.Note that node relocations are allowed as long as they do not violate the feasibility of the trip of origin and do not produce further infeasibilities in the trip of destination.

1.
Let LI be a list of second trip vehicles with violated minimum load constraints and LF a list of feasible second trip vehicles, both sorted in non-increasing order of capacity utilization.

2.
Following the lists, for each pair of vehicles k and k (k in LI and k in LF): 2.1.Identify the demand nodes served by k that can relocated to k, maintaining the trip of k feasible and not producing additional infeasibilities in the trip of k.Create a list LP comprising the node with the smallest demand that makes the trip of k feasible, or in case no such a node exists, of the nodes sorted in non-increasing order of demand.

2.2.
While LP = ∅ and the trip of k is infeasible, insert the first node ( p) listed in LP in the insertion least costly position identified in step 2.1 and make LP = LP − {p}.After each node relocation, update LP by removing demand nodes that can no longer be relocated to the trip of k and the capacity utilization of vehicles k and k .

3.
Update LI by removing vehicles whose trips are no longer infeasible.If LI = ∅, redefine list LF with first trip vehicles sorted in non-increasing order of capacity utilization and redo steps 2 and 3.

4.
If there are still infeasible second trips, deconstruct these trips by returning their demand nodes to set NCR and their vehicles to set FA, and apply the insertion heuristic (described in Fig. 4) from step 2.

5.
Return the resulting solution S. In step 1, unfeasible second trip vehicles regarding minimum load constraints are sorted in list LI in non-increasing order of capacity utilization, and feasible second trip vehicles are sorted in list LF, also in non-increasing order of capacity utilization.The vehicles sorting in both lists is justified by the fact that unfeasible trips carried out by vehicles with larger capacity utilization require less additional load to become feasible, while feasible trips have a larger number of demand nodes that can be relocated to other trips without the former become infeasible.
A list LP is then built with nodes that can be relocated from a given vehicle k in LF to a vehicle k in LI.LP either comprises the node with the smallest demand that makes the trip of k feasible (unitary list), or in case no such a node exists, it includes all nodes sorted in non-increasing order of demand (step 2.1).The first node in LP is then inserted in the least costly position of vehicle k's trip.This procedure is repeated until LP = ∅ or all trips in LI become feasible (step 2.2).
If there are still infeasible second trips (step 3), list LF is redefined with first trip vehicles.The prioritization of feasible second trips over feasible first trips for node relocation is due to the fact that second trips mostly comprise less constrained demand nodes, which facilitates their insertion in other trips.
Remaining infeasible trips are deconstructed in step 4, with their demand nodes returned to the set of unserved nodes NCR and their vehicles returned to set FA.The procedure then takes advantage of possible route slack resulting from the node relocations in step 2.2 and idle vehicles to include unserved nodes in the solution.

The local search procedure
The local search procedure (LS) proposed in this work is inspired by the tabu search algorithm by Pureza et al. ( 2012) and aims at improving the construction heuristic solution by applying three phases in sequence.Firstly, it seeks to reduce the number of vehicles/routes (with the possibility of creating new routes for serving unrouted nodes), followed by the minimization of the total distance travelled by own vehicles, and finally, it attempts to decrement the crew size in each trip.Figure 6 illustrates the steps of the procedure.Note that the phases sequencing follows the hierarchy of the problem's objectives.

1.
From solution S obtained by the construction stage described in Figure 3: 1.1.Apply the route reduction phase by relocating nodes, one by one, from the trip r with the least number of nodes to other routes.The nodes are inserted in the least costly feasible position.Let S be the resulting solution.
If r is completely emptied, make S = S .

1.2.
Apply the distance reduction phase to S by performing an intra-route 2opt, an inter-route two-node exchange, and an inter-route one-node insertion.Each neighborhood is searched until a neighbor solution with shorter own vehicle travelled distance is found or all moves are investigated.Repeat the step until no improved solution is found for all three neighborhoods.Let S be the resulting solution.

1.3.
Apply the crew size reduction phase to S by decrementing the number of deliverymen in each trip, one by one, until any additional decrement makes the trip infeasible or the crew size is equal to 1. Let S be the resulting solution.While improvements are obtained in the first and third phases (steps 1.1 and 1.3, respectively) by using a single type of move each, the second phase (step 1.2) consists of a sequence of three move types: 2-opt (intra-route), exchange of two nodes in different trips (inter-route), and insertion of a node from trip i to a distinct trip j (inter-route).According to the results of extensive preliminary experiments, this sequence of moves shows the best trade-off between time and solution quality.Only improving moves are accepted, and the selection of the neighbor solution follows the first-improve policy, i.e., each neighborhood is searched until a solution with shorter own vehicle travelled distance is found or all neighbors are investigated.The types of moves are applied on a cyclical basis until no improved solution is found for all three neighborhoods.

COMPUTATIONAL EXPERIMENTS
The purpose of the experiments discussed in this section is to verify the applicability of model ( 1)-( 18) and our solution approach and to compare the results obtained with the company's estimated solutions.For this, we collected data on some daily deliveries in a major city served by the company, and present the solutions of 14 instances.As depicted in Table 1, the instances range from 9 to 209 customers and are divided in 5 sets according to the city area type (Central, Peripheral and Dangerous) or market segment (Supermarkets and All segments).The real data provided by company consists of the customers' location and demand, the available fleet, as well as the routes sequencing in the company's solutions.On the other hand, the travelled distance and scheduling of the company's routes had to be estimated (thus, also the routes costs) because no information regarding customer clusters was available.As discussed in Section 2, the routes planned by the logistic operators assume that the trucks travel from one customer to another, even though in practice, clusters take shape when customers close to the truck parking site are served on foot by the crew.For this reason, the travelled distance, route times and costs of the company's solutions were computed considering the most probable clusters, defined a posteriori by examining the actual route sequencing, truck capacity and the maximum radial distance between a customer and the truck parking site adopted by the deliverymen (150 meters).Distances between any two nodes were computed considering the city real network.Once the clusters were defined, route times were calculated based on the demand of each cluster, average traveling speed, reloading time for a second trip (60 minutes) and service carried out by the standard crew (two deliverymen).As discussed in Section 2, the company provided us estimated service times for a two-man crew, hence we set the service time in a given cluster as the sum of the estimated service time of the customers in the cluster.
Instances with distances and times computed according to previous paragraph are denoted by RC01, RC02 and RC03 (set C); RP01, RP02 and RP03 (set P); RS01 and RS02 (set S); RD01, RD02 and RD03 (set D); and RA01, RA02 and RA03 (set A).Instance RA01 comprises the data of RC02, RP01 and RD01; instance RA02 comprises the data of RC01, RP02 and RD02 and instance RA03 comprises the data of RC03, RP03 and RD03.Their solutions correspond therefore to what we call the company's estimated solutions.
The corresponding instances used to solve model ( 1)-( 18) are, in turn, denoted by MC01, MC02 and MC03 (set C); MP01, MP02 and MP03 (set P), MS01 and MS02 (set S); MD01, MD02 and MD03 (set D); and MA01, MA02 and MA03 (set A).Note that the only difference between RC01..RA03 and the corresponding MC01..MA03 is the clusters formation, which in the latter case was done a priori, according to the following requirements: (i) clusters are formed exclusively by either customers with time window or customers without time windows; (ii) the radial distance between each customer in a cluster and the truck parking site cannot exceed 150 meters; (iii) the total demand in a given cluster must be less than or equal to the capacity of the smallest truck that can serve the corresponding region (central, peripheral, dangerous).Single customer clusters are formed if the customer's demand is greater than the capacity of the smallest truck.
The modeling language GAMS with solver CPLEX 12.5.0.1 was used to implement and solve model ( 1)- (18).Based on real data, we set model parametersP i , Co, Cs and G to 900, 1.1, 500 and 1, respectively.CPLEX's branch-and-cut ran with standard parameters, except for options fpheur = 1, heurfreq = 100, lbheur = 1, threads = 4, and maximum processing time (τ max ) equal to 18000 seconds (5 hours).GRASP was implemented in Oracle 11g with the programming language PL/SQL and run with parameter α = 30%.According to the company, each route is planned in 3 minutes, so as an attempt to reproduce the conditions under which the route planning is carried out, the maximum processing time employed in GRASP was set to the size of the available fleet multiplied by 180 seconds.For the sake of simplicity, for the cases with one and three deliverymen (modes 1 and 3), we simply divided the estimated service times for a two-man crew (provided by the company) by 1/2 and 3/2, respectively.From the perspective of the company's deliverymen, these estimates are close enough to their expectations of service time.
All experiments were conducted in a Dell computer, model Optiplex 9010 with processor Intel i7 with 3.4 GHz, 16 GB RAM and operational system Windows 7 Professional with 64 bits.

Results
We first compare the solutions obtained with CPLEX with the corresponding company's estimated solutions.Regarding the latter, Table 2 presents for each instance of the sets, the solution cost, the fleet size employed, the distance travelled by the vehicles, the total crew size (excluding the security guard, and equal to 2 per vehicle), and if the solution is feasible according to the model constraints.For the CPLEX's solutions, Table 2 presents the number of clusters defined according to the requirements discussed in the beginning of Section 5, the solution cost, the fleet size employed, the distance travelled by the vehicles, the total crew size (variable and equal to the sum of the largest delivery crews in the two possible trips of each vehicle), the number of unserved clusters, and the solutions gap between the best feasible MIP solution found and the LP relaxation.Table 3 summarizes the comparisons by depicting the average percent difference between GAMS/CPLEX's solutions and the company's estimated solutions.
As shown in the sixth column of Table 2, the company's solutions are infeasible for all instances but one, either for violating route time, time windows, truck capacity or minimum load in the second trip.Even though route scheduling had to be estimated, the company's logistics operators confirm that such violations are, in fact, very common.As an example, the company's estimated solution for instance RC02 exhibits violation of the time window upper limit of 3 clusters of restaurants, ranging from 47 to 110 minutes beyond the maximum limit of 120.We note that if we increase the time window upper limit of the only unserved cluster in the corresponding instance MC02 from 120 to 180 minutes (i.e., a 60 minute violation), CPLEX succesfully incorporates the cluster in a route, resulting in a less infeasible option than the company's solution in less than 100 seconds of computing time.
Table 3 shows that CPLEX produced average reductions in cost, number of vehicles, route time, and travelled distance of 41%, 2%, 30%, 31%, respectively.The improvements were achieved at the expense of only 6% average increase in the number of deliverymen.We note that the large crew size increments in sets C and P compensate the overtime hours present in the company's solutions (set P) or the clusters' time windows violation (set C).In fact, overtime accounts for an average of 2:45 h per vehicle in set P and time windows violations accounts for an average of 18 min per cluster with time window in set C. Similarly, the large decrement in set S crew size indicates that some of the vehicles in the company's solutions are overmanned.This can be clearly observed in the results of instances RS01 (the only case with a feasible company's estimated route plan) and MS01.Both solutions have equal or very similar costs, number of vehicles and total distance, whereas the number of deliverymen is four men larger in the company's solution.That is, setting the crew size based on the specificities of each trip can indeed bring benefits, either by impeding overtime/time windows violations or overmanning.
It should be noted though that the cost reductions obtained in instance set A are the result of unserved clusters, some of which could certainly be served, since the smaller instances that make them up have feasible solutions.For example, instance MA01 comprises the data of instances MC01, MP01 and MD01; when these are solved individually, all clusters are served while MA01 has 3 unserved clusters.This observation indicates that the computational time provided to GAMS/CPLEX (18000 seconds) is insufficient for tackling the largest instances.Moreover, despite the general good results, the incumbents of 10 of the 14 instances were found very close to the maximum computational time.Given the practical motivation for this research,    we felt necessary to verify the quality of the solutions obtained within the route planning time used in practice (3 minutes per vehicle).Table 4 presents the results.
We note in Table 4 that limiting the computation time to the company route planning time results in solutions substantianly inferior, with increases of the travelled distance cost of 58%, 11% and 6% in instances MP01, MD02 and MS02, respectively.For instances MC01, MC02, MP03, MD01, MD03 and MA01, the gaps are substantially higher since the average served clusters goes down to 32%.In other four instances (MC03, MP02, MA02 and MA03) the computational time is not even sufficient to obtain a feasible solution.We therefore conclude that the time required by CPLEX's branch-and-cut (and possibly, other exact methods) to generate good quality solutions, makes its use unviable in practical contexts, as one might expect, and points to the investigation of heuristic or hybrid methods.Table 5 presents the results from the application of our GRASP implementation, limited to the route planning time practiced by the company.We also include the number of unserved clusters and the gaps of the solutions found with GAMS/CPLEX for a quick comparison and Table 6 summarizes the differences.
Table 5 shows that GRASP produced solutions with smaller gaps than GAMS/CPLEX in 43% of the instances as the result of a larger number of served clusters.For example, among the 372 clusters of set A, GRASP fails to serve 2 clusters while GAMS/CPLEX's solutions exhibit 27 unserved clusters.Average gaps with GAMS/CPLEX and GRASP are equal to 2.79% and 1.17%, respectively, indicating that the heuristic approach is capable of generating better quality solutions in much shorter computational times.However, regarding instances MD02 and MD03 (dangerous areas), GRASP's solutions present an average of 4.3% of unserved clusters.CPLEX's routes, on the other hand, incorporated all clusters in these two instances.

A simple hybrid approach
As means to investigate the performance of a hybrid method, we solve the 14 instances with GAMS/CPLEX from the best solution produced by GRASP.Since these best solutions are usually obtained within 40% of the computational time, we set the time for the application of the metaheuristic and GAMS/CPLEX as 0.4 (Fleet size ×180) seconds and 0.6 (Fleet size ×180), respectively.That is, the time provided to the hybrid method is equal to the route planning time practiced by the company.Table 7 shows the results, along with the number of unserved clusters and the gaps of the solutions found by GAMS/CPLEX without the starting solution and τ max = 18000 seconds.Table 8 summarizes the comparisons.
The application of the hybrid method (within the company route planning time) improved 6 of the 14 solutions of GAMS/CPLEX (obtained within 18000 seconds of processing time) regarding the number of served clusters.In particular, Table 7 shows that service is provided to all clusters of instances MC02, MC03, MP03 and MA02, despite the fact that the travelled distance is 6.5% shorter.When compared to the solutions provided by GRASP, the average percentage of unserved clusters in instances MD02, MD03, MA01 and MA03 decreases from 2.4% to 0.9%, while solutions for set D exhibit an increased number of served clusters in 2 instances and reductions of 16.5% on the travelled distance.For set S, the improvements are even larger, disclosing an average of 22.2% decrease on the number of vehicles employed.In addition, reductions of 1, 7 and 12 deliverymen are observed in the solutions of instances MP02, MS01 and MS02, respectively, representing a 38.4% savings on the crew size prescribed by GRASP.
We conclude that the performance of the hybrid method is superior to the other proposed methods, as it combines mostly high quality solutions to low computational time.The method can be considered satisfactory for all instances but MD03, MA01 and MA03; in these cases, an average of 1.2% of the clusters remained unserved (note in Table 2 that all clusters in instance MD03 can actually be served).Table 9 summarizes the comparisons between the hybrid method solutions and the company's solution.

CONCLUSIONS AND PERSPECTIVES OF FUTURE RESEARCH
This paper addressed the distribution operation of large volumes of products in dense urban areas.
We proposed an optimization model that incorporates most of the relevant conditions imposed to  the fleet that carry out the distribution in a large Brazilian drinks producer and distributor, disclosing the complexity and some of the difficulties pertinent to this important logistic activity.In addition to merging three variants of the Vehicle Routing Problem (multi-trip, multi-deliverymen and with time windows), inherent to the studied application, a highlight of the proposed model is the diversity of constraints, characterizing the delivery operation as a rich vehicle routing problem.As the company is typical in the drinks sector, this study may be helpful to others with similar delivery or pickup operations.Sets of instances based on real data provided by the company and comprising 1 to 3 market areas were solved with the branch-and-cut algorithm of the software GAMS/CPLEX, a GRASP metaheuristic and a simple hybrid method that applies GAMS/CPLEX from the incumbent solution provided by the metaheuristic.In particular, the hybrid method was capable to generate high quality solutions for most of the instances within the route planning time required by the company, in these cases presenting an impressive average cost reduction of 37% relative to the company's estimated solutions.It should be emphasized that the resulting solutions were analyzed and validated by the company's logistic operators, indicating the potential of the application of operations research techniques in real environments.The cases for which the methods failed to route all demand clusters motivate the development of more sophisticated metaheuristics or hybrid methods.
Since limited experiments during this research revealed a close relation between cluster formation and the routing quality, another interesting avenue of research is to study the problem that combines clustering to the routing/scheduling/crew size decisions.Finally, considering uncertainties in travel times and/or demand with the application of robust optimization techniques also follows as a suggestion of future research.For instance, it is clear that urban centers are particularly susceptible to traffic jams, which may not justify the trucks travelling times based on an average speed adopted in the present study.
The problem combines characteristics of the VRPTWMD addressed in Pureza et al. (2012) and the MTVRPTW presented in Seixas & Mendes (2013), with several additional constraints.The flow network (Fig. 2) is represented by a graph G(n + 2, A) with three types of nodes: n − 1 demand nodes (represented by the associated parking sites) and three origins and/or destinations of the trips.Note that the warehouse is both origin and destination of trips, therefore, it is represented by the following three copies: nodes O 1 (origin of the first trip), O 2 (destination of the first trip and origin of the second trip) and O 3 (destination of the second trip).

Figure 3 -
Figure 3 -Steps of the GRASP algorithm.

Figure 4 -
Figure 4 -Steps of the insertion heuristic.

Figure 5 -
Figure 5 -Steps of the feasibility procedure.

Figure 6 -
Figure 6 -Steps of the local search heuristic.

(
the two possible trips of each vehicle. build LRC with demand nodes ∈ (NCR ∩ CT) − NCI (i.e., unserved nodes with time windows).Otherwise ( j = 2), build LRC with α% of the best evaluated nodes with no time windows from set NCR − CT − NCI.
3.2.While LRC = ∅: 3.2.1If FA = FD or new = true, start a new trip.Randomly select a node p from list LRC and the first vehicle k from list LV that can serve p without violating any of the problem constraints, disregarded the minimum load in the second trip.If there is no such ( p, k), make NCI = NCI + {p} and go to step 3.1.Otherwise, start the trip associated to vehicle k with node p and the maximum possible crew.Make LV = LV − {k}, update the capacity of vehicle k and schedule the trip.Otherwise, if FA = FD and new = false, insert nodes in the current solution.

Table 1 -
Sets of instances.

Table 3 -
Average percentage difference between GAMS/ CPLEX's solutions and the company's estimated solutions.

Table 4 -
GAMS/CPLEX performance with two computational times.
max = Fleet size ×180 seconds.Largest crew in the two possible trips of each vehicle.

Table 6 -
Average percentage difference between GRASP's solutions and GAMS/CPLEX's solutions.

Table 8 -
Average percentage difference between the hybrid method's solutions and GAMS/CPLEX's solutions.

Table 9 -
Average percentage difference between the hybrid method's solutions and the company's estimated solutions.