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 setup, 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 programmingbased Adaptive Large Neighborhood Search (cpALNS) 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 cpALNS 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 onebyone 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 optimizationbased heuristic. Rasmussen et al.
[2] solves a similar problem with an exact branchandprice 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 (SAILS) 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 (DTRSP) which deals with a limited crew of technicians serving a set of requests. In the DTRSP, 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 smallsize instances than the cpALNS of [8] in a much shorter running time. Our lpALNS dominates the cpALNS 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 lpALNS 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 subsets: 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 nonnegative 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

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.

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 ).

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 nonoverlapping intervals, and any interval in the chain is constrained to end before the start of the next interval in the chain.

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 .

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 .

Time window constraints.
(11) startOf represents the start of interval whenever the interval variable is present.

Capacity constraints.
(12) presenceOf is equal to 1 if interval variable is present in the solution, 0 otherwise.

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 programmingbased 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 pseudocode 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.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 tiebreaking 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 precomputing 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.
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 lpALNS 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 lpANLS 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) E52667. Note that our reference algorithm (the cpALNS 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 lpALNS 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 lpALNS 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 lpALNS with those of cpALNS 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 cpALNS. The columns “Gap” report the average gaps (in percentage) between solution costs of our methods and those of cpALNS. 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 modelbased algorithms cannot solve any instance to optimality in 3 hours, they do provide quite good solutions. Remarkably, the CP model performs better than cpALNS on 56 instances in a much shorter running time (5 minutes vs a couple of hours of cpALNS). 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 lpALNS. It is the most efficient method in terms of solution quality, as it provides 128 better solutions compared to cpALNS and is worse on only 7 over 165 instances. Moreover, the gaps, on average, are negative on all instance classes (except RC12). Our lpALNS can averagely improve the objective values up to 4.59% (instance class R1).
Data  Sync  lpALNS  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 
The second experiment is to investigate the performance of the lpALNS 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 cpALNS in terms of solution quality. It provides better than average results on all instance classes. Figure 3 reports the number of new bestknown solutions found by the lpALNS for each class of instances. In total, 620 bestknown solutions have been improved by our lpALNS.
Moreover, Table 2 in the Appendix shows that the improvement on the value of the objective function created by the lpALNS 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 cpALNS, 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.
5 Conclusion
In this research, we study an important variant of the vehicle routing problem with synchronization constraints, which has numerous realworld 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 CPbased metaheuristic for small instances in much shorter running time. Our lpALNS dominates the cpALNS, 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 lpALNS 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 lpALNS is still quite timeconsuming. 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: 297316.
 [2] Rasmussen, M.S., Justesen, T., Dohn, A., & Larsen, J. (2012). The home care crew scheduling problem: Preferencebased visit clustering and temporal dependencies. European Journal of Operational Research 219.3: 598610.
 [3] Rousseau, L.M., Gendreau, M., & Pesant, G. (2013). The synchronized dynamic vehicle dispatching problem. INFOR: Information Systems and Operational Research 51.2: 7683.
 [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: 1931.
 [5] Afifi, S., Dang, DC., & Moukrim, A. (2016). Heuristic solutions for the vehicle routing problem with time windows and synchronized visits. Optimization Letters 10.3: 511525.
 [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, 347367.
 [7] Parragh, SN., & Doerner, KF. (2018). Solving routing problems with pairwise synchronization constraints. Central European journal of operations research 26.2: 443464.
 [8] Hojabri, H., Gendreau, M., Potvin, JY., & Rousseau, L.M. (2018). Large neighborhood search with constraint programming for a vehicle routing problem with synchronization constraints. Computers & Operations Research 92, 8797.
 [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. 148162). 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. 662673.
 [13] Ham, A.M. (2018). Integrated scheduling of mtruck, mdrone, and mdepot constrained by timewindow, droppickup, and mvisit using constraint programming. Transportation Research Part C: Emerging Technologies 91, 114.
 [14] Pisinger, D. & Røpke, S. (2007). A general heuristic for vehicle routing problems. Computers & Operations Research 34: 24032435.
 [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. 417431.
Appendix
Table 2 reports the comparison between lpALNS and cpALNS 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 cpALNS. Next columns show the results on average for each instance class of our lpALNS. More precisely, “InitialObj” is the objective value of the initial solution constructed by regret2 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 lpANLS and cpALNS solutions. A negative value means lpALNS provides a better solution.
Comments
There are no comments yet.