A new constraint programming model and a linear programming-based adaptive large neighborhood search for the vehicle routing problem with synchronization constraints

10/18/2019 ∙ by Minh Hoàng Hà, et al. ∙ 0

We consider a vehicle routing problem which seeks to minimize cost subject to time window and synchronization constraints. In this problem, the fleet of vehicles is categorized into regular and special vehicles. Some customers require both vehicles' services, whose starting service times at the customer are synchronized. Despite its important real-world application, this problem has rarely been studied in the literature. To solve the problem, we propose a Constraint Programming (CP) model and an Adaptive Large Neighborhood Search (ALNS) in which the design of insertion operators is based on solving linear programming (LP) models to check the insertion feasibility. A number of acceleration techniques is also proposed to significantly reduce the computational time. The computational experiments show that our new CP model finds better solutions than an existing CP-based ANLS, when used on small instances with 25 customers and with a much shorter running time. Our LP-based ALNS dominates the cp-ALNS, in terms of solution quality, when it provides solutions with better objective values, on average, for all instance classes. This demonstrates the advantage of using linear programming instead of constraint programming when dealing with a variant of vehicle routing problems with relatively tight constraints, which is often considered to be more favorable for CP-based methods.

READ FULL TEXT VIEW PDF
POST COMMENT

Comments

There are no comments yet.

Authors

This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

1 Introduction

Most of the requirements related to our daily routines are made by a service provider coming to our premises. These types of services can be home care delivery, maintenance operations, public utilities, etc. In such services, efficient delivery and timely service play important roles. This is why the class of the Vehicle Routing Problem, coupled with the Scheduling Problem, comprises a large class of problems with many variations and applications. The main focus of this research is to study the Vehicle Routing Problem with Synchronization Constraints (VRPSC), where both time window and synchronization constraints are present. In the latter constraint, some customers may require the service of two vehicles whose starting service times at the customer must be synchronized.

In a recent industrial project with an Internet Service Provider (ISP) in Vietnam, the authors witnessed several contexts where synchronization constraints arose. The ISP company has to perform installation services for new subscribers and maintenance services for subscribed clients. Both services are performed by technicians who mainly use motorbikes for travelling. In many cases, a customer asks for services by two technicians belonging to two different teams: one being the ”physical” team, which takes charge of the hardware (cable wire, modem, etc.); while the other team manages the signal. To further illustrate the problem, when a customer requires a service from the physical team, two technicians must be mobilized as one helps the other with equipment set-up, such as installing a ladder and other protective gear. In addition, an intern technician, who is in a probationary period, needs to be coupled with an experienced technician at customer locations. When servicing a customer, the company requires that the starting server time of both technicians be as close as possible in order to reduce their waiting time and to limit the disruption to the customer. In the case of our ISP partner, a delay of no more than 15 minutes is permitted. The problem was initially introduced in [8] and can be used to model other real world applications such as home care delivery, aircraft fleet assignment, ground handling, and forest operations (see [3] for more information).

Hojabri et al. [8] proposed a constraint programming-based Adaptive Large Neighborhood Search (cp-ALNS) with insertion operators exploiting constraint propagation capabilities to guarantee the feasibility of a new generated solution. Different from the popular ALNS proposed in [14] to solve VRPs, the cp-ALNS does not try to add all unserved requests one by one but, rather, adds all of them at once to create a new complete solution. Several removal operators were specifically designed for the problem. Numerical results are reported on instances derived from Vehicle Routing Problem with Time Window (VRPTW) benchmark instances, with up to 200 customers and 100 synchronizations.

A dynamic version of the problem was introduced in [3] where customer requests were not foreseen, but arrived one-by-one in real time. The problem was modeled as a CP program from which a metaheuristic was designed. We note that the methods proposed in [3, 8] were based on the same CP model which uses variables representing the successor of a node, and subtour elimination constraints introduced in [18]. However, no result of the pure CP model has been reported.

Apart from the paper of [8], a problem quite similar to the VRPSC is studied in [4]

in the context of home care crew scheduling. The problem is first formulated as a mixed integer programming (MIP) model that constitutes pairwise synchronization and pairwise temporal precedence between customer visits. It is then solved through an optimization-based heuristic. Rasmussen et al.

[2] solves a similar problem with an exact branch-and-price algorithm. Due to the application context, there is no capacity constraint and specific issues about home care crews are taken into account, in addition to the synchronization requirements, like care giver preferences, customer priority and the ability of a particular care giver to serve a given customer. Also, not all customers must be serviced because visits can be rescheduled or canceled.

Afifi et al. [5] propose a simulated annealing algorithm with dedicated local searches (SA-ILS) for the VRPTW with synchronized visits. Their problem differs from the considered problem in this work in two ways; first, the synchronized nodes must be serviced at the same time, i.e., the delay is zero and the fleet of vehicles is not categorized. Second, three objective functions have been considered: minimizing the total travel time; minimizing the sum of negative preferences; and minimizing the maximum difference in service times of the vehicles. Recently, Parragh et al. [7] evaluate several different ways to deal with pairwise synchronization constraints in the context of two problems: the VRPTW with pairwise synchronization, and the service technician routing and scheduling problem. They propose three ways to address the synchronization requirement: individual synchronized timing optimization; global synchronized timing optimization; and adaptive time window. The idea of the first two approaches is to keep the ALNS untouched, while the last one makes some modifications to the insertion scheme in order to identify good starting service times for synchronized visits.

Pillac et al. [6] introduce a Dynamic Technician Routing and Scheduling Problem (D-TRSP) which deals with a limited crew of technicians serving a set of requests. In the D-TRSP, each technician has a set of skills, tools and spare parts required by each request. In addition to designing a route at the beginning of each day, two types of decisions must be managed in real time. First, whenever a new request appears, we must decide whether it is accepted or not. And second, whenever a technician finishes serving a request, we need to find the next request to serve. For a survey on further synchronization issues in the context of vehicle routing problems, we refer the readers to [1].

Our main contributions are as follows: we propose a new CP model and a metaheuristic based on ALNS to address the VRP with time windows and synchronization constraints. Different from the CP model in the literature, our CP uses sequence and interval variables of IBM ILOG CP Optimizier [9] to formulate the problem. The most notable feature of our metaheuristic is we use linear programming and a number of acceleration techniques to quickly check the feasibility of insertion operations integrated in the popular ALNS [14]. This is the first time LP models are used instead of CP to design the ALNS algorithm for the VRP with synchronization constraints. The computational experiments carried on the benchmark data show the performance of our methods. More precisely, our CP model can provide better solutions on small-size instances than the cp-ALNS of [8] in a much shorter running time. Our lp-ALNS dominates the cp-ALNS in terms of solution quality, as it improves 620 over 681 best known solutions. The remainder of the paper is organized as follows: Section 2 introduces the problem definition and our new CP model; the detailed description of the lp-ALNS is provided in Section 3; experimental results are reported in Section 4; and finally, we conclude our work in Section 5.

2 Problem definition and a new constraint programming model

The problem may be formally defined as follows: given , an undirected graph, where is the set of nodes representing customer locations and is the set of edges. is the depot where a set of vehicles is located. Vehicles in set are again divided into two sub-sets: regular vehicles (set ) and special vehicles (set ). All regular vehicles have a capacity of . is the set of regular customers which are visited by regular vehicles only while is the set of special customers, each requires the visits of both types of vehicle. Let be the set of vertices which are the copies of special customers . Let be defined as . It is considered as the set of customer vertices that must be visited by regular vehicles. Each vertex in is associated with a demand . Define and () be the vertices at which the vehicle of type starts from and comes back to, respectively. Here, represents regular vehicles and is the index for special vehicles. Note that, the vertices and () share the same location as . Additionally, we define be the set of vertices appearing on routes of regular vehicles; and be the set of vertices visited by special vehicles. A service time is associated with a vertex , and a service time is with a vertex . Note that the service times at the depot and its copies are set to null. A time window () is imposed on each vertex . Finally, each edge is associated with non-negative values and representing the travel cost and travel time between vertices and of vehicle type .

The problem then consists in constructing routes for the fleet of vehicles such that the total travel cost incurred by the fleet of vehicles is minimized and the following constraints are satisfied:

  • Each vehicle must begin the route at the depot, deliver services to customers and finally return to the depot.

  • Each regular customer is served by exactly one regular vehicle.

  • Each special customer is served by exactly one regular vehicle and one special vehicle.

  • The total demand serviced by a regular vehicle must not exceed its capacity .

  • A regular vehicle must start its service at a vertex within the time window .

  • Starting service time at a special customer visited by a special vehicle must be within a time window . Here, and are given parameters representing a possible delay between regular and special services at vertex . And, is the starting service time at vertex (the copy of in regular vehicle route).

We now present the new CP model for the problem. First, it is worth mentioning that unlike MIP formulations, there is no standard in CP formulation because it strongly depends on each CP package. In this study, we formulate the model using generic keywords and syntaxes of IBM ILOG CP Optimizier [9] adapted from the CP formulations proposed by [10, 11, 12, 13]. Our model uses the following variables:

Variables

interval variable that represents the time interval of size vertex is visited
optional interval variable that represents the time interval vehicle visits vertex ;
sequence variable that represents all working time intervals of vehicle ;
interval variable that represents the time interval of size vertex is visited;
optional interval variable that represents the time interval vehicle visits customer ;
sequence variable that represents all working time intervals of vehicle ;

An interval variable represents the interval of time during which a task can occur. It contains a starting point, an end point, a size, and it can be optional. A decision variable is used to represent whether or not an interval is present. If an interval is marked as optional, it may be absent in the solution. Sequence is a type of variable in IBM ILOG OPL which can be empty or can contain a subset of variables. A sequence represents all intervals that are present in the solution. The constraints in our model are as follows:

Constraints

  1. Function to link each interval to a location.

    (1)
    (2)

    A value is defined for each pair of (, ) and (, ). The type of the interval variable is its last known location.

  2. Each customer must be served by one vehicle of the corresponding type.

    (3)
    (4)

    alternative function ensures that exactly one set of intervals (or ) is present in the solution; and interval variable starts and ends together with the interval variable (or ).

  3. Travel time between two customers must be taken into account. In the following, is the matrix representing travel times between two vertices and in the set .

    (5)
    (6)

    noOverlap constraint on sequence variables and states that the sequence defines a chain of non-overlapping intervals, and any interval in the chain is constrained to end before the start of the next interval in the chain.

  4. All vehicles start their route at the starting depot.

    (7)
    (8)

    first function states that if interval is present, it will be the first interval in the sequence . These two constraints force all regular vehicles to start their routes at and all special vehicles to start at .

  5. All vehicles finish their routes at the corresponding ending depot.

    (9)
    (10)

    Similar to first() function, last() function states that if interval is present, it will be the last interval in the sequence . These two constraints are to force all regular vehicles finishing their routes at ending depot and all special vehicles finishing their routes at ending depot .

  6. Time window constraints.

    (11)

    startOf represents the start of interval whenever the interval variable is present.

  7. Capacity constraints.

    (12)

    presenceOf is equal to 1 if interval variable is present in the solution, 0 otherwise.

  8. Synchronization constraints.

    (13)
    (14)

Objective function

 

We compute the total cost traveled by regular and special vehicles as follows:
= where
     = typeOfNext .
= where
     = typeOfNext .

Then the objective function can be written as:

(15)

3 Linear programming-based adaptive large neighborhood search algorithm

To tackle the VRPSC problem, we design an Adaptive Large Neighborhood Search (ALNS) heuristic, which is based on a Large Neighborhood Search (LNS) introduced by [19]. At each iteration, LNS explores a large neighborhood, which can rearrange a large part of the current solution, therefore allowing the search to move to other promising search spaces.

More precisely, LNS decomposes the original problem by unfixing some decision variables, leading to a partial solution. The unfixed decision variables define a neighborhood of solutions that can be explored by a specific procedure via, possibly, a heuristic or a Mixed Integer Programming (MIP) solver. If the procedure finds an improved solution, it becomes the new current solution and a new large neighborhood is defined around it. This process is repeated until a stopping criterion is reached.

A first key point is the selection of fixed variables to create a partial solution. In fact, the number of fixed variables impacts the size of the neighborhood (the more fixed variables, the narrower the neighborhood). A common strategy is to dynamically vary the number of removed variables. A second key point lies in the selection of fixed/removed variables which can use a random choice or a more sophisticated strategy to guide the search. Finally, the procedure that explores the neighborhood should provide good quality solutions in a short amount of time. Adaptive Large Neighborhood Search is an extension of LNS with a number of different insertion and removal operators. In comparison with LNS, a component that adaptively chooses among a set of removal and insertion heuristics is added to the algorithm. The pseudo-code of a ALNS to solve problems with minimizing objective function is shown in Algorithm 1

. At each iteration, a randomly selected pair of operators (with procedures SelectDestruction and SelectRepair, lines 5 and 6) is applied to the current solution (line 7), with probabilities (respectively, sets

, ) updated by a learning process (line 12). The more an operator has contributed to the solution process, the larger probability it has of being chosen.

1 Create an initial solution s Initialize and while the stop-criterion is not met do
2       := SelectRemoval() := SelectInsertion() solution = GenerateNewSolution() if cost() cost() then
3            
4      if accept(s’, s) then
5            
6      update and
return
Algorithm 1 General ALNS

3.1 Insertion operators

3.1.1 Cheapest insertion heuristic

The purpose of the insertion operation is to reinsert unserviced requests into solution. For this task, one can use the cheapest insertion heuristic which inserts the customer into a route at a feasible position making the objective value increase the least. The process is repeated until all customers are serviced or no more customers can be inserted. The insertion cost of a regular customer into a regular route positioned between two consecutive vertices and (denoted by ) is computed as:

(16)

Whenever a special customer is considered for insertion, it will be added to two positions: one on regular routes and another on special routes. This must be incorporated when computing the insertion cost of special customers. The average value is used to compute the insertion cost of a special customer at positions between vertices and on regular routes and between vertices and on special routes as follows:

(17)

3.1.2 Regret heuristics

As in [14], we also use regret- heuristics as repair operators. Instead of selecting the customer with the least insertion cost in each construction step, the regret heuristics select the customer with the highest regret- value, computed as follows: we denote is the insertion cost when inserting customer into the best position of route . If this insertion is infeasible w.r.t time window and synchronization constraints, the insertion cost is set to infinity, i.e. . Let be the route on which vertex has the -th lowest insertion cost. The regret- value of customer is then calculated as:

(18)

The regret- heuristics choose the unvisited customer with the highest regret- value and insert it into the feasible position leading the least insertion cost. Ties are broken by selecting customers with the lowest insertion cost . Informally speaking, we choose the insertion that leads to the most regret, if it is not done at present. In some situations, if a vertex can be inserted somewhere in the current solution, but cannot be inserted in at least routes, then the vertex that can be inserted in the fewest number of routes is selected. This ensures that the vertex which does not have many insertion options in the current solution will be considered first.

It can be observed that the cheapest insertion heuristic, which we mentioned earlier, is the special case of the regret heuristic with due to the tie-breaking rule. For any , the regret heuristic looks further into future solutions to decide the choice of insertion. In this research, we use regret heuristics with to design insertion operators of our ANLS.

3.1.3 Checking insertion feasibility of regular customers

When inserting a vertex into a position of the current partial solution, it is required to verify if the insertion satisfies the capacity, time window and synchronization constraints. As the insertion operation is repeated multiple times during the search, designing a quick verification procedure is critical to speed up the overall algorithm. As the capacity constraint is easily checked in , we focus on the time window and synchronization constraints only. Verifying the feasibility of an insertion operation w.r.t these constraints are more complex because they delay subsequent visits leading to other violations. As proposed in [15], the time window constraint can be checked in by pre-computing the maximum delay (push forward) that is allowed at each arc of the current solution, without violating time windows. In this research, we also reuse this idea to handle both time window and synchronization constraints when inserting regular vertices.

Given a partial solution, in order to consider all possible positions to insert an unserved regular customer into the routes, we calculate the maximum duration of time (also called maximum delay) that can be spared after a vehicle finishes serving vertex and before it starts to serve next vertex without violating any constraints of other nodes. This value is denoted as , where represents the arc from to . We calculate the maximum delays at all the arcs of the current solution using the following linear programming model:

Let and be the set of all vertices visited by regular and special vehicles in current solution , respectively. Denote is the set of arcs forming the routes in ; and are the sets of arcs on regular and special routes, respectively. We use two types of variables: implying starting service time at customer and representing the maximum delay on arc .

(19)
(20)
(21)
(22)
(23)
(24)
(25)

Objective (19) is to maximize the maximum delay on arc . Constraints (20) and (21) respectively ensure time window and synchronization constraints at all vertices of the current solution. Constraints (22) represent the relationship between the starting times at vertices and when a vehicle travels from to . Constraint (23) has the same meaning as constraints (22) but is written for arc . Finally, constraints (24) and (25) define the domain of variables.

After all the maximum delays are available, we check if an unserved regular node can be inserted at the position between node and node by computing the earliest arrival time () and the waiting time () at as follows:

Finally, we check if feasibility of the insertion satisfies the following constraints:

(26)
(27)

Constraint (26) ensures that the insertion does not lead to a violation of time window and synchronization constraints at all the customers in the current solution. Constraint (27) verifies the time window constraint of vertex .

Although it is fast to solve a LP model of type (F1), the running time of the insertion operators is still expensive due to the large quantity of LP models solved during the search. Through observation, we note that constructing the model to find maximum delay on each arc takes more computational time than solving it. As such, we propose the following model to reduce the running time of the insertion operators:

(28)
(29)
(30)

Objective (28) is to maximize the weighted maximum delay at all arcs of the solution. Here, is a given binary coefficient representing the weight of arc . Constraints (29) indicate the relationship between the variables and . To find the maximum delay on an arc , we just need to set its weight to 1 and weights of other arcs to null. The proposed model (F2) allows us to save a lot of time by constructing a model once and then creating a new one by changing two coefficients in the objective function only. As a result, we can avoid constructing multiple models from scratch. A preliminary experiment shows that using this method helps reduce at least 30% running time of the overall algorithm. After each model calculating the maximum delay at an arc with two extremities to is solved, the value of the variables , , and are also saved for checking insertion feasibility of special customers described in the following section.

3.1.4 Checking insertion feasibility of special customers

Whenever a special customer is selected to be added into the current solution, it will be inserted into two positions: one on regular routes and the other on special routes. Multiple insertion operations make it impossible to use the maximum delay for the feasibility verification purpose. The following LP model, without objective function (F3), is used to check if a special vertex can be added in arc on a regular route and arc on a special route:

(31)
(32)
(33)
(34)
(35)
(36)
(37)
(38)

The meaning of variables , other parameters, and constraints (31)-(33) can be derived from (F1) and (F2). Constraints (34)-(37) are similar to (33), but written for 4 new arcs created by the insertion of vertices and its mirror : (), (), (), and ().

As mentioned above, after solving each model computing the maximum delay on an arc (), the starting service times of nodes and , or in other words, the values of and are saved. can be seen as the earliest time the vertex can be serviced while can be seen as the latest time vertex can be serviced. To avoid misunderstanding these notations, we denote as and as . Using these saved values, we can efficiently check if special node and its copy can be inserted on arc of a special route and arc of a regular route, respectively. First, lower bound (denoted by ) and upper bound (denoted by ) of arrival times at when being inserted on arc and when being inserted on arc are computed as follows:

=

=

=

= .

It can be seen that the synchronization constraints will be violated in the following two cases:

(39)

The insertions of and also need to satisfy the time window constraint at node and the maximum delays computed from the model F2. Thus, we can utilize this property to rapidly verify the feasibility of the insertions. Unlike regular vertices, possible waiting times at and are created by not only time window constraints, but also synchronization constraints. The lower bounds of waiting times created by time window () and synchronization () at a node can be computed as follows:

= max(0, )

= 0

= max(0, )

= max(0, )

Figure 1 illustrates the computation of waiting times in case of and . In Figure 1a, the vehicle arrives at the customer location before the time window and has to wait 30 minutes before starting delivery. In Figure 1b, the special vehicle arrives before the regular vehicle, but it has to wait until the regular vehicle starts to service the customer because the value of is zero. Thus, the waiting time due to the synchronization constraint, in this case, is 45 minutes.

Figure 1: Computing waiting times due to time window and synchronization constraints

Hence, we can calculate the lower bounds of waiting times at node by taking the maximum value between and :

= max(, )

= max(, )

After all the values above are calculated, we can skip the insertions which violate one of the following constraints:

(40)
(41)
(42)

The constraints (40) and (41) verify if the insertions satisfy the maximum delays while constraint (42) checks the time window at the regular vertex . In addition, based on the characteristic of the -regret heuristics, we can only consider the insertions if their cost is smaller than the current -th best insertion cost.

It is worth mentioning that validation procedures (39)-(42), which run in if the values of the variables of each program F2 are available, can detect plenty of infeasible insertions, thus saving a lot of run time of the algorithm. A preliminary experiment on instances with 50 customers and 25 synchronizations shows that our fast validation procedures help to increase the algorithm’s speed by up to 20 times. This ratio increases on larger instances with more synchronizations. However, if the insertion passes all the validations, we cannot ensure that the insertion is indeed feasible w.r.t the time window and synchronization constraints. As a consequence, whenever an insertion of a special request passes the checks above, we will continue to examine if that insertion is feasible by solving a program of type (F3). If that model returns a solution, then the insertion is feasible, and infeasible otherwise.

3.2 Removal operators

The destroy operators remove a fraction of vertices from a complete solution based on different criteria, each guiding the algorithm to another search space. The input of the operators is a complete solution and their outputs are vertices that have been removed from the . We use three destroy operators originally proposed by [14]: random removal, related removal and worst removal.

3.2.1 Random Removal

This is the simplest removal operator. It randomly selects vertices in the solution and removes them. Other vertices remain unchanged. This obviously helps the algorithm diversify the search.

3.2.2 Related Removal

The idea of the related removal, as its name indicates, is to remove similar vertices with the expectation that they could interchange their positions to create a better solution. More specifically, to measure the similarity between two vertices and , we use the relatedness which is calculated as follows:

(43)

To calculate the relatedness of two vertices, we take into account the differences of four characteristics: their starting service times (); the distance between them (); their demand size; and their type. Note that the value of is obtained from solving a program of type F3 whenever a new complete solution is found; and is set to 1 if vertex is special, and to 0 if it is regular. The difference is normalized such that it only takes values from the interval [0, 1]. In this formula, , , and indicate the largest starting service time, the largest distance, and the largest demand of all vertices in the solution, respectively. In addition, each characteristic is associated with a weight to measure its importance.

3.2.3 Worst Removal

The worst removal operator removes the vertices that are very expensive, with the expectation that these vertices might be located in wrong places. Given a request served by some vehicle in a solution , we define the cost of the vertex as the difference between the cost of and the cost of the new solution where vertex is removed completely from . The worst removal heuristic repeatedly chooses a vertex with the largest cost until vertices have been removed.

To add more diversification to our algorithm, the related and worst removal operators are randomized by removing the -th most related (or expensive) request where is the set of vertices in solution and is a random number in [0, 1], and parameter is used to control the randomization. If is small, the most related (in the case of related removal) or expensive (in the case of worst removal) vertex is selected, while less related (or expensive) vertices may be chosen for larger values of with a probability that decreases with the cost . The values of are taken from [14].

Finally, our lp-ALNS also uses acceptance criteria embedded in a simulated annealing framework, adaptive score adjustment to select operators in a dynamic fashion, and adding noise to insertion cost to increase the diversification. All these components and their parameter settings are taken from [14] without any change.

4 Computational results

In this section, the effectiveness of the proposed algorithms is examined. We test our algorithms on the instances proposed in [8] with the number of customer = 25, 50, 100, and 200. These instances are generated from the VRPTW instances of [16, 17] containing three types, depending on the customers’ distribution. The customers are randomly located in the instances of type R, clustered in type C, and mixed between randomly located and clustered in type RC. The instances are also categorized into two classes based on the capacity of vehicles. The first class (including C1, R1, RC1) consists of instances with a relatively small capacity compared to the total customer demand, while in the second class (C2, R2, RC2), the capacity is relatively large. Note that in the instances of types R and RC, the vertices are identically distributed in class 1 and 2, while this is not true for type C. And finally, in these VRPTW instances, the travel time and travel cost between two vertices are set to their Euclidean distance.

The original VRPTW instances are transformed into VRPSC instances as follows: the number of special customers is set to where is the percentage of special customers. There are three values for : 5%, 25%, and 50%. More precisely, the first customer in the VRPTW instances is considered a special customer and the next special customers are selected using a constant interval defined by . In the case of the synchronization constraint, the values of and are set to 0 and 10 for every special customer , respectively. And finally, by private contact, it turns out that the authors in [8] report inexact results for three instances C101, C105 and C106 with 25 customers and 2 synchronizations, so we removed these instances from our experiments.

The CP model is coded in IBM OPL 12.8.0 while the lp-ANLS is implemented in C++ using CPLEX 12.8.0 for the resolution of the linear programs. Both methods are run on a 3.20GHz Xeon(R) E5-2667. Note that our reference algorithm (the cp-ALNS of [8]) was run on a 3.07GHz Xeon(R) X5675, which is a similar generation to our processor. Since different CPU speed conversion techniques can provide very different results, we decided to present the raw running time, letting the readers choose their preferred approach. The parameter setting of the lp-ALNS is chosen empirically. We have tested many settings and the following setting gives the best performance, in terms of both quality and computational time for our algorithm. The number of removed vertices in the removal operators is a random integer between 4 and . In related removal operator, the values of - are set to 4, 2, 1, and 4, respectively. The lp-ALNS stops after 25000 iterations. All the detailed results can be found in http://www.orlab.vn/home/download.

In the first experiment, we compare the results obtained by our CP model and lp-ALNS with those of cp-ALNS proposed by [8]. Because our CP model cannot handle the large instances, we chose the small instances with 25 customers for the experiment. The limited running time of the CP approach for each instance is set to 5 minutes and 3 hours. Table 1 shows the number of times each of our methods finds better solutions (Columns “Better”), equal solutions (Columns “Equal”), and worse solutions (Columns “Worse”) compared to cp-ALNS. The columns “Gap” report the average gaps (in percentage) between solution costs of our methods and those of cp-ALNS. The negative values in these columns indicate that our methods provide better solutions in terms of objective function values. The results obtained show that although our CP model-based algorithms cannot solve any instance to optimality in 3 hours, they do provide quite good solutions. Remarkably, the CP model performs better than cp-ALNS on 56 instances in a much shorter running time (5 minutes vs a couple of hours of cp-ALNS). It can be observed that CP models work better on the instances of the first class (C1, R1, and RC1) and worse on the instances of second class. The instances with shorter routes tend to be easier for our CP model. The results clearly show the performance of our lp-ALNS. It is the most efficient method in terms of solution quality, as it provides 128 better solutions compared to cp-ALNS and is worse on only 7 over 165 instances. Moreover, the gaps, on average, are negative on all instance classes (except RC12). Our lp-ALNS can averagely improve the objective values up to 4.59% (instance class R1).

Data Sync lp-ALNS CP (5 min) CP (3h)
Better Equal Worse Gap Better Equal Worse Gap Better Equal Worse Gap
R1 2 12 0 0 -3.22 9 0 3 -0.82 11 0 1 -2.30
R2 11 0 0 -2.49 1 0 10 10.18 1 0 10 9.69
C1 6 0 0 -2.16 6 0 0 -1.45 5 0 1 -1.33
C2 2 6 0 -0.83 1 3 4 2.72 1 3 4 2.72
RC1 7 1 0 -0.79 3 0 5 5.24 5 0 3 1.05
RC2 4 3 1 0.12 0 0 8 9.77 0 0 8 10.92
R1 7 12 0 0 -4.59 7 0 5 6.36 11 0 1 -2.37
R2 8 0 3 -1.98 0 0 11 13.41 0 0 11 9.57
C1 6 3 0 -2.25 5 0 4 -1.68 5 0 4 -1.95
C2 3 5 0 -0.49 2 3 3 4.00 2 3 3 4.00
RC1 8 0 0 -2.80 5 0 3 3.53 6 0 2 -0.95
RC2 7 0 1 -2.93 1 0 7 9.93 0 0 8 8.95
R1 13 12 0 0 -3.82 6 0 6 7.10 8 0 4 -1.02
R2 9 1 1 -2.74 0 0 11 12.14 0 0 11 8.60
C1 5 3 1 -1.78 3 0 6 1.42 4 0 5 -0.48
C2 2 6 0 -1.84 4 0 4 3.23 4 0 4 2.84
RC1 8 0 0 -1.72 3 0 5 4.58 4 0 4 3.35
RC2 6 2 0 -2.25 0 0 8 15.03 0 0 8 14.60
Table 1: CP and lp-ALNS vs cp-ALNS on 25-customer instances

The second experiment is to investigate the performance of the lp-ALNS on all the instances. The computational results are summarized in Figures 2 and 3, and Table 2 in Appendix. Figure 2 shows that our algorithm clearly dominates the cp-ALNS in terms of solution quality. It provides better than average results on all instance classes. Figure 3 reports the number of new best-known solutions found by the lp-ALNS for each class of instances. In total, 620 best-known solutions have been improved by our lp-ALNS.

Moreover, Table 2 in the Appendix shows that the improvement on the value of the objective function created by the lp-ALNS is significant, especially on large instances (up to 16.75 %). The relatively high gap between final solutions and initial solutions shows the efficiency of construction and deconstruction operators. However, similar to cp-ALNS, the computation time of our algorithm is still high. It depends heavily on the number of customers and that of special customers . More specifically, in these cases, the number of variables and constraints in the LP models increase rapidly, leading to larger programs which are harder to solve.

Figure 2: Comparison between lp-ALNS and cp-ALNS in terms of objective values on average
Figure 3: Number of new best-known solutions found by the lp-ALNS

5 Conclusion

In this research, we study an important variant of the vehicle routing problem with synchronization constraints, which has numerous real-world applications. We propose a new CP model and an ALNS algorithm. The most remarkable feature of our ALNS is that we use linear programming to check the feasibility of insertions. A number of acceleration techniques have been proposed to significantly reduce the computation time of the algorithm. The obtained results on the benchmark instances from the literature show the performance of our method. Our CP model can even provide better solutions than the CP-based metaheuristic for small instances in much shorter running time. Our lp-ALNS dominates the cp-ALNS, in terms of solution quality, when it improves 620 best known solutions over 681 instances, and the improvement gap is relatively high.

The research perspectives are numerous. First, although our CP model provides very good solutions on small instances, it cannot prove any of them optimal. This, combined with the fact that there is no efficient exact method so far to solve the VRP with synchronization constraint, proves the hardness of this class of problem. An efficient and exact method is still an open question. Second, we believe that our lp-ALNS can be a used as a general framework, as it is easy to incorporate other constraints into the LP models. Thus, applying our method to solve other hard variants of VRPs could be an interesting research direction. Finally, our lp-ALNS is still quite time-consuming. Other acceleration techniques exploiting special structures of the LP programs which validate the insertion feasibility are required to make the algorithm become an efficient general solver for VRPs with rich attributes.

References

  • [1] Drexl, M. (2012). Synchronization in vehicle routing—a survey of VRPs with multiple synchronization constraints. Transportation Science 46.3: 297-316.
  • [2] Rasmussen, M.-S., Justesen, T., Dohn, A., & Larsen, J. (2012). The home care crew scheduling problem: Preference-based visit clustering and temporal dependencies. European Journal of Operational Research 219.3: 598-610.
  • [3] Rousseau, L.-M., Gendreau, M., & Pesant, G. (2013). The synchronized dynamic vehicle dispatching problem. INFOR: Information Systems and Operational Research 51.2: 76-83.
  • [4] Bredström, D., & Rönnqvist, M. (2008). Combined vehicle routing and scheduling with temporal precedence and synchronization constraints. European journal of operational research 191.1: 19-31.
  • [5] Afifi, S., Dang, D-C., & Moukrim, A. (2016). Heuristic solutions for the vehicle routing problem with time windows and synchronized visits. Optimization Letters 10.3: 511-525.
  • [6] Pillac, V., Guéret, C., & Medaglia, A. L. (2018). A fast reoptimization approach for the dynamic technician routing and scheduling problem. Recent Developments in Metaheuristics. Springer, Cham, 347-367.
  • [7] Parragh, S-N., & Doerner, K-F. (2018). Solving routing problems with pairwise synchronization constraints. Central European journal of operations research 26.2: 443-464.
  • [8] Hojabri, H., Gendreau, M., Potvin, J-Y., & Rousseau, L.-M. (2018). Large neighborhood search with constraint programming for a vehicle routing problem with synchronization constraints. Computers & Operations Research 92, 87-97.
  • [9] IBM Software, (2015). IBM ILOG CPLEX Optimization Studio V12.6.3.
  • [10] Laborie, P. (2009). IBM ILOG CP Optimizer for detailed scheduling illustrated on three problems

    . In International Conference on AI and OR Techniques in Constraint Programming for Combinatorial Optimization Problems (pp. 148-162). Springer, Berlin, Heidelberg.

  • [11] Ghédira, K. (2013). Constraint satisfaction problems: CSP formalisms and techniques. John Wiley & Sons.
  • [12] Goel, V., Slusky, M., van Hoeve, W.J., Furman, K.C., & Shao, Y. (2015). Constraint programming for LNG ship scheduling and inventory management. European Journal of Operational Research, 241(3), pp. 662-673.
  • [13] Ham, A.M. (2018). Integrated scheduling of m-truck, m-drone, and m-depot constrained by time-window, drop-pickup, and m-visit using constraint programming. Transportation Research Part C: Emerging Technologies 91, 1-14.
  • [14] Pisinger, D. & Røpke, S. (2007). A general heuristic for vehicle routing problems. Computers & Operations Research 34: 2403-2435.
  • [15] Kindervater, G.A.P.. & Savelsbergh M.W.P. (1997). Vehicle routing: handling edge exchanges. in Aarts E.H.L., Lenstra J.K. (eds.), Local Search in Combinatorial Optimization, John Wiley & Sons, pp. 337–360.
  • [16] Solomon, M.M. (1987). Algorithms for the vehicle routing and scheduling problems with time window constraints. Operations Research 35, 254–265.
  • [17] Homberger, J., & Gehring, H. (1999). Two evolutionary metaheuristics for the vehicle routing problem with time windows. INFOR 37, 297–318.
  • [18] Pesant, G., Gendreau, M., Potvin, J.-Y., Rousseau, J.-M. (1998).

    An Exact Constraint Logic Programming Algorithm for the Traveling Salesman Problem with Time Windows

    . Transportation Science, 32: 12–29.
  • [19] Shaw P. (1998). Using Constraint Programming and Local Search Methods to Solve Vehicle Routing Problems. In Principles and Practice of Constraint Programming, Goos, G., Hartmanis, J., van Leeuwen, J., eds., Lecture Notes in Computer Science 1520, Springer, Berlin, pp. 417-431.

Appendix

Table 2 reports the comparison between lp-ALNS and cp-ALNS in terms of objective values on average. The three first columns represent the size, class name, and the number of special customers of each instance class. Columns 4 and 5 report the objective of solutions (Column “FinalObj”) and running time in seconds (Column “RunTime”) on average of the cp-ALNS. Next columns show the results on average for each instance class of our lp-ALNS. More precisely, “InitialObj” is the objective value of the initial solution constructed by regret-2 heuristic. “FinalObj” is the objective value of the final solution after the search is stopped. “RunTime” reports the computational time in seconds. Column “Imp%” shows the improvement of final solutions compared with initial ones. The final column “gap%” presents the gap between lp-ANLS and cp-ALNS solutions. A negative value means lp-ALNS provides a better solution.

width=height=0.5 Data cp-ALNS lp-ALNS gap (%) Size Class Sync FinalObj RunTime InitialObj FinalObj RunTime Imp% 25 R1 2 544.2487 8406.7 667.7745 526.9739 1347.8 20.86% -3.22% R2 458.2109 9978.2 594.1048 446.6154 1319.3 24.42% -2.49% C1 246.9088 9255.3 371.2795 241.5390 1292.7 32.26% -2.16% C2 286.4098 7643.6 342.2520 283.9314 1415.9 16.63% -0.83% RC1 520.1329 6805.4 584.4385 515.7725 1480.5 11.23% -0.79% RC2 485.1764 7672.3 624.6635 485.7430 1650.5 21.82% 0.12% 25 R1 7 692.9658 7846.3 888.4291 661.7204 2697.3 25.68% -4.59% R2 587.3948 10289.6 753.2626 575.4996 3011.7 22.83% -1.98% C1 318.3507 8304.2 485.0793 311.0783 2489.7 34.81% -2.25% C2 325.0070 7308.5 411.1374 323.3936 2090.1 20.95% -0.49% RC1 623.5873 6388.2 900.0518 605.3236 2707.9 31.69% -2.80% RC2 573.4323 8828.2 796.9595 556.5850 2987.9 29.59% -2.93% 25 R1 13 821.1631 7683.6 1113.333 790.7901 3989.3 28.71% -3.82% R2 701.3905 10194.9 930.8132 681.8161 4523.5 26.21% -2.74% C1 356.2599 6594.2 585.0626 349.7486 3481.7 39.33% -1.78% C2 384.0774 6977.4 461.5924 376.7076 3511.1 17.86% -1.84% RC1 657.5279 5892.3 1009.331 644.6960 3844.1 35.93% -1.72% RC2 600.0163 8504.3 905.7041 586.7275 4460.4 33.84% -2.25% 50 R1 3 909.5803 19518.2 1143.763 860.4969 4508.1 24.75% -5.69% R2 747.6849 23888.9 999.7746 710.4967 4496.9 28.77% -5.03% C1 460.1213 39268.2 554.8499 445.4759 4532.7 19.14% -3.13% C2 476.6964 35246.2 633.3015 466.1748 4804.5 25.20% -2.02% RC1 959.1396 17810.0 1165.249 932.1801 4539.8 19.44% -2.98% RC2 788.3569 20046.5 1201.248 775.9619 4857.4 35.07% -1.65% 50 R1 13 1178.034 19893.2 1508.424 1103.480 10834.1 26.81% -6.57% R2 991.7101 23972.1 1352.777 931.9242 13253.3 30.88% -6.03% C1 625.0501 18575.3 951.2536 607.6679 9646.9 35.13% -2.70% C2 610.4964 22449.2 907.3918 603.0134 13477.1 31.94% -1.20% RC1 1231.836 17850.4 1664.564 1167.640 11213.3 29.66% -5.15% RC2 1014.296 22600.4 1629.658 975.5894 12489.1 39.79% -3.96% 50 R1 25 1354.380 19456.9 1786.241 1259.770 18880.4 29.63% -7.15% R2 1134.506 24311.6 1580.263 1057.765 26470.6 32.82% -6.96% C1 689.7339 16392.7 1100.512 674.8498 18270.7 38.23% -2.09% C2 697.0899 20694.2 1069.437 668.8730 23789.4 35.90% -3.85% RC1 1331.455 17536.3 2064.870 1270.645 18259.1 38.31% -4.30% RC2 1087.594 23088.4 1858.651 1041.364 21808.5 43.37% -4.50% 100 R1 5 1429.048 82716.3 1847.463 1349.893 18128.8 26.94% -5.87% R2 1114.927 107910.0 1587.720 1038.766 19082.2 34.41% -6.98% C1 1017.524 41401.6 1493.110 1005.144 21153.7 32.45% -1.20% C2 841.6708 75868.3 1156.707 792.8305 20527.3 30.57% -5.74% RC1 1637.621 73411.3 2187.973 1575.960 19985.8 27.97% -3.86% RC2 1286.366 97137.4 1983.771 1219.061 20101.6 38.63% -5.47% 100 R1 25 1789.477 80930.0 2353.730 1667.065 39167.0 29.28% -7.01% R2 1407.660 106302.7 2108.117 1303.016 55645.6 37.91% -7.66% C1 1441.884 74134.6 2335.434 1396.531 38251.0 39.86% -3.13% C2 1087.373 70059.3 1670.794 1004.834 39562.4 38.45% -7.24% RC1 2146.339 71074.1 2985.869 2045.495 40795.9 31.44% -4.73% RC2 1709.736 99529.7 2594.230 1610.915 51314.6 37.95% -6.06% 100 R1 50 2092.569 72804.5 2831.534 1931.541 59256.5 32.07% -7.91% R2 1693.735 104309.4 2526.176 1513.136 89195.5 40.16% -10.91% C1 1609.269 72829.5 2911.811 1545.038 57106.8 46.57% -3.91% C2 1180.336 54500.8 1927.523 1078.205 69467.0 41.82% -8.42% RC1 2522.189 67991.7 3650.564 2361.475 60008.8 35.39% -6.45% RC2 1945.313 95459.4 3173.071 1817.931 63274.9 42.44% -6.83% 200 R1 10 4144.381 126805.8 5674.585 3893.875 37843.3 31.74% -6.55% R2 3713.917 161933.9 5317.305 3320.119 37547.4 37.72% -11.05% C1 3391.796 100174.6 5108.380 3299.623 39677.3 34.62% -2.67% C2 2578.896 145855.5 3861.907 2364.314 36008.3 38.13% -8.26% RC1 4088.110 138911.3 5642.754 3884.939 38358.8 31.03% -4.95% RC2 3303.763 149639.7 4898.392 2917.336 39697.6 40.55% -11.91% 200 R1 50 5074.646 136440.0 7240.237 4724.183 74670.6 35.08% -7.46% R2 4769.591 163572.0 6608.872 4074.128 93781.0 38.47% -14.99% C1 4238.077 120709.5 7290.512 4126.756 74785.5 43.11% -2.65% C2 3254.876 161222.8 5851.064 3008.252 81829.4 48.25% -7.49% RC1 4944.824 144934.9 7017.894 4573.909 78424.0 34.65% -7.50% RC2 4039.977 155597.6 6108.987 3539.513 92402.6 41.75% -12.49% 200 R1 100 6139.267 137336.0 8718.394 5704.425 116804.2 34.72% -7.55% R2 5616.638 157111.4 7993.287 4699.330 154303.7 41.40% -16.75% C1 4964.759 111991.6 9117.560 4757.915 112905.5 47.57% -4.16% C2 3646.121 132923.9 6969.890 3374.519 146172.7 50.71% -7.31% RC1 5836.127 148733.9 8265.386 5395.134 122680.9 34.72% -7.58% RC2 4882.749 159016.6 7138.363 4204.200 153580.6 41.12% -13.99%

Table 2: Comparison between lp-ALNS and cp-ALNS