Temporal constraints and device management for the Skill VRP: Mathematical model and lower bounding techniques *

,


Introduction
Resource constrained routing and scheduling problems are very relevant in theory and practice, and thus they have attracted considerable attention by the scientific literature. In these problems, customers have specific requests that can be met by specialized resources, like technicians, vehicles and devices. In addition, the travel and the delivery of products or services to the customer locations have to be designed. The resources are usually classified as renewable (e.g. vehicles), non-renewable (e.g. spare parts), or doubly constrained (e.g. electric energy).
According to a recent review paper [27], the most popular problems in this field are the Skill VRP and the Technician Routing and Scheduling Problem. The Skill VRP has been introduced in [4] and further investigated in [5]. There, technicians with given skills must perform routes to serve customers, each requiring a set of skills to be served. The aim is to minimize the total routing costs while satisfying the constraints related to the required levels of skill. An extension to the Skill VRP has been studied in [32]. Another generalization has been proposed in [28], where the service of each customer requires one or more resources (e.g. technicians and equipment). There, the resources of each type are limited, and each of them must be assigned to a single route.
The Technician Routing and Scheduling Problem, introduced in [29], can be seen as another generalization of the Skill VRP dealing with multiple resources. Specifically, each technician has a set of skills and may in addition carry a set of tools and spare parts, while each customer requires a subset of them. Tools are renewable resources, whereas spare parts are non-renewable. The goal is to design minimum duration routes for the technicians so that each customer request is served by one technician with the required skills, tools and spare parts. Also in this case the resources of each type are limited, and each of them is carried out by a single technician. Variants to the Technician Routing and Scheduling Problem will be discussed in the section devoted to the literature overview. Interestingly, one of these problems has been the object of a recent challenge facilitated by VeRoLog, the EURO Working Group on Vehicle Routing and Logistics Optimization ( [14]). Specifically, special tools for measuring milk quality have to be delivered to customers, at their request, over a multiple day time horizon. After the measurement, the tools have to be picked up again. The scheduling of these deliveries along the time horizon, and the routing for the planned deliveries and pickups, are the key decisions to address. The studied problem is indeed a simplification of a richer version faced by the service provider. In fact, as outlined in [23], the real problem includes also the scheduling and the routing of some inspectors, who should visit the customers while the measuring tools are present. Such inspectors are different than the technicians in charge of transporting the measuring tools, since they must possess specific qualifications. The need of technicians with different levels of skill, and the simultaneous presence of qualified technicians and devices, with the treatment of the underlying synchronization issues, would make the problem an even greater challenge to address.
The aim of this paper is to investigate in this direction of research, by addressing the management of technicians with different levels of qualification, and the simultaneous presence of a special device and of a properly qualified technician at some customer locations. Such aspects make the problem particularly significant in application contexts such as Home Care and Field Service. The relevance of considering special devices in Home Care is well described in [1], where the authors present a case study to evaluate hospital readmission rates and mortality in elderly patients with chronic obstructive pulmonary disease. The treatment provided at home to such patients requires in fact to consider two types of materials: consumable materials, i.e. nonrenewable (e.g. oxygen replenishment), and specialized equipment, i.e. renewable (e.g. Doppler ultrasonograph and electrocardiograph). See also [31] for the description of a dimensioning problem in Home Care where the importance of managing devices is well outlined.
More in detail, with respect to the kernel Skill VRP problem presented in [4] and in [5], here multiple visits per customer are allowed, i.e. more than one technician may visit a customer in order to operate his required services, although each operation must be performed by a single technician. Furthermore, the generalization considers precedence constraints among the operations, time windows at the customers, and maximum workday duration constraints, which are not addressed in [4] and in [5]. Finally, the handling of a special device, which is required to perform some special operations, and which must be present at the customer location together with a technician qualified to use it, to the best of our knowledge constitutes an original contribution to the literature on Resource constrained routing and scheduling problems. The special device is moved by the technicians during their own tours, independently of their qualifications, and therefore non standard synchronization between the movement of the special device and the tours of the technicians must be imposed, giving rise to a peculiar tour for the special device, defined as a composition of some technician tour fragments. The resulting Skill VRP generalization appears thus to be relevant from an operational perspective, and it may constitute a basic block for interesting multiple day planning extensions.
The plan of the paper is the following. Section 2 proposes an overview of papers dealing with some of the characteristics of the generalized Skill VRP. Section 3 describes the studied problem, while Section 4 presents the Mixed Integer Linear Programming (MILP) model proposed for its mathematical formulation. Some simple, but computationally quite effective, valid inequalities to enhance the Linear Programming (LP) relaxation of the MILP model are presented in Section 5. Finally, Section 6 presents some lower bounding techniques, based on the proposed MILP formulation, and the results of computational experiments, showing that some lower bounding techniques, when enhanced via the proposed valid inequalities, may rapidly produce good lower bounds. The returned percentage optimality gaps, estimated also thanks to simple matheuristics, are in fact quite small for several scenarios of medium to large size, by encouraging the use of the proposed lower bounding techniques both as building blocks for designing exact approaches for the generalized Skill VRP, and also as valuable tools to evaluate the efficacy of more sophisticated heuristic approaches to the problem. Section 7 concludes the paper.

The literature
To the best of our knowledge, in the literature of Resource constrained routing and scheduling problems no paper has dealt with skill, precedence, synchronization and time windows aspects simultaneously with the management of a special device, which must be present at the customer location together with the technician in charge to use it. Rather, skill, precedence, synchronization and time windows constraints have been addressed, with different emphasis also depending on the considered application context. Also, delivery and pickup at the customer locations of renewable resources, such as tools and devices, and of non-renewable resources, such as spare parts, have been addressed, by disregarding however their synchronization with the visits of the technicians who are qualified to use them. Therefore, after listing recent reviews dealing with some of these issues in the Health Care application area, hereafter we shall provide a short overview of papers addressing some of these characteristics in the context of Resource constrained routing and scheduling problems.
Some Health Care reviews.
Cissé et al. [11] detail a comprehensive overview of recent models to Home Health Care routing and scheduling problems (HHCRSP). The authors analyze the existing literature according to the way the constraints and the objective function are formulated, provide an overview of methods developed to solve HHCRSP, and discuss future research directions.
Fikar and Hirsch [19] review relevant problem settings for HHCRSP. They consider singleperiod and multi-period Home Health Care problems and analyze objective functions, constraints and solution methods addressed in the literature.
Paraskevopoulos et al. [27] put in evidence that there exist several interesting variants of HHCRSP, the Skill VRP and the Technician Routing and Scheduling problems being the most prominent. The authors also show that maintenance activities and Home Health Care are the main areas where routing and scheduling of resources are crucial not only in terms of customer satisfaction, but also in terms of operational efficiency.

Time windows and preferences.
Cheng and Rich [10] consider the nurse rostering and routing problem. They propose a MILP formulation for the problem together with a heuristic approach. In particular, time windows are considered. A differentiation is made by considering full-time and part-time nurses, having different costs. Lunch breaks in different times and with different durations are addressed as well.
Berthels and Fahle [2] combine nurse rostering with vehicle routing aspects. They describe the software PARPAP developed for solving the considered problem. In particular, they use a Linear Programming model, constraint programming techniques and metaheuristics. The paper considers time intervals, preferences and other additional constraints as soft constraints.
Nickel et al. [26] consider scheduling and routing in Home Health Care over a weekly time horizon. Hard time windows are present. A two-stage solution approach is proposed which is based on a constraint programming heuristic and an adaptive large neighborhood search in order to improve the initial solution.
In a more general VRP setting, time windows are also addressed in [15]. There, a stochastic variant of the VRP with time windows is considered, where travel times are assumed to be stochastic, and some restrictions are placed on the probability of violation of time window constraints.

Precedence and synchronization.
Bredstrom and Ronnqvist [3] present a set partitioning model for the combined vehicle routing and scheduling problem with time windows and additional temporal constraints. The temporal constraints allow one to impose pairwise synchronization and pairwise temporal precedence between customer visits, independently of the vehicles. The authors propose an optimization based heuristic to solve real size instances.
Korsah et al. [24] address the problem of optimally assigning and scheduling a set of spatially distributed tasks to a fleet of vehicles working together to achieve a high-level goal, in domains where tasks may be related by precedence or synchronization constraints and might have a choice of locations at which they can be performed. Such problems may arise, for example, in disaster preparedness planning, transportation of people, and delivery of supplies. The authors present a novel mathematical model and describe how to solve it to optimality within a Branch-and-Price framework.
Rasmussen et al. [30] address scheduling and routing in Home Health Care, by considering time windows, soft preference constraints and temporal dependencies between the start times of the visits.
In the context of Home Health Care services, synchronization together with time window constraints are also treated in [16], [20] and [21]. In [16], a constructive heuristic is presented together with an iterated local search metaheuristic, whereas in [21] a general variable neighborhood search approach is provided. [20] extends the results in [21] by proposing variable neighborhood search approaches to coordinate multiple Health Care structures and multiple specialties handled by them.
In the same application context, we mention also [12] where some visits may require the presence of two staff members simultaneously, thus requiring visit synchronization. Hard and soft time windows are also taken into account together with multiple visits at some patient locations. A MILP model and a memetic algorithm featuring two specific crossover operators are proposed to solve the problem.
For a survey of vehicle routing problems with multiple synchronization constraints, related to spatial, temporal and load aspects, we refer to Drexl [13]. There, the authors present a classification of different types of synchronization, i.e. task synchronization, operation synchronization, movement synchronization, load synchronization and resource synchronization. Moreover, they comprehensively review the pertinent literature with respect to the application settings as well as with respect to successful solution approaches, both at an exact and at a heuristic level, by emphasizing that this kind of VRP problems are of practical relevance in many application areas, and they considerably complicate the use of standard VRP solution techniques.

Skills of staff.
Eveborn et al. [17] describes the LAPS CARE system, developed and used in Sweden for the daily planning of Home Care activities, by taking into account the skills of the staff members. The problem is formulated in terms of a set partitioning model, and solved via a repeated matching algorithm.
In a more general VRP setting, Cappanera, Gouveia and Scutellà [4] introduce the Skill VRP, where a set of tours, each one operated by a skilled technician, must be designed in such a way that each service requirement, asking for particular skills, is fulfilled by exactly one technician [4,5]. Various MILP models and related valid inequalities are proposed in [4]. Specifically, (i) the models are based on multicommodity flow variables and constraints, (ii) increasing levels of variable disaggregation allow one to strengthen the associated LP bounds, and (iii) the models strongly exploit skill information at the nodes of the logistics network, which is the key aspect differentiating the proposed models from more classical VRP models.
In [5], the authors formulate three MILP models, again characterized by an increasing level of variable disaggregation, and test them on a large suite of randomly generated instances. In particular, it is shown that the behavior of one of such models can be improved further by a cutting plane approach, so allowing to address large scale instances arising in challenging real application contexts.
Skill requirements and related constraints, in the context of assignment, routing and scheduling in Home Care, are also addressed in [6,7,9,35]. See [8] for skill aspects in Home Care under uncertain patient demands.
Pickup and delivery of renewable and non-renewable resources.
In [29], each technician has a set of skills and may carry a set of tools and spare parts to the customers, each requiring a subset of these resources. Tools are renewable resources, whereas spare parts are non-renewable since they are consumed when the technician serves a request. The technicians start their tour from a depot with a set of spare parts and tools, and they also have the option to replenish their tools and spare parts at any time at the depot. The goal is to design minimum duration routes for the technicians so that each customer request is served by one technician, with the required skill, tools, and spare parts, within a given time window. The authors developed a parallel adaptive large neighborhood search algorithm, where a shared pool of promising solutions is maintained. A post-optimization procedure based on a set covering model is then used. Computational experiments on instances with up to 100 service requests are reported.
By considering a Home Health Care context, [25] studies the problem of delivering drugs and medical devices from the company's pharmacy to the patients, and special drugs from the hospital to the patients. Moreover, the pickup of biological samples and unused drugs and medical devices from the patients to the depot or to a medical lab is addressed. Each patient must be visited by one vehicle within a given time window, and vehicles are capacitated. Two MILP models are proposed together with a genetic algorithm and a tabu search method.
Always considering the Home Health Care setting, [33] addresses the case where a set of patient requests have to be accomplished by two types of caregivers in a one-day planning horizon. The first category of caregivers is responsible for primary tasks and to collect biological samples. These caregivers start their routes at a depot and end their working day at the initial depot right after the delivery of the collected samples at a laboratory. The perishability for transferring the biological samples to the laboratory is taken into account. Subsequently, the second category of caregivers performs synchronized services, that can be simultaneous or must respect precedence constraints, within a predetermined time window. In contrast to the first category, such caregivers do not visit the laboratory. The goal is to minimize the transportation cost and the idle time of the caregivers. Simulated annealing and tabu search algorithms are applied in two phases, i.e. to schedule primary tasks in the first phase, and to synchronize services in the subsequent phase.
Finally, the already mentioned problem in [23] concerns a cattle company, that must regularly measure the milk quality at a number of farms, over a multiple day time horizon, by means of special measuring tools. Such tools have therefore to be delivered to the customers at their request. After a given number of days elapsed since delivery, the tools have to be picked up again. Tools are characterized by kind, size and availability. Each request involves only tools of the same kind. Vehicles transporting them are capacitated and the distance they can travel in one day is limited to a given maximum value. During the day, vehicles can visit the depot several times and load/unload tools. On each day, tools unloaded by one vehicle at the depot are not available to other vehicles. But, at the end of a day, the tools still on board of a vehicle are unloaded and made available to vehicles for the next day. Moreover, tools can be picked up at a customer and delivered at another customer during the same day.
It has to be decided on which day to serve each delivery request and, consequently, on which day to pick up the delivered tools. Furthermore, for each day in the planning horizon it has to be decided how to combine deliveries and pickups in the vehicle routes. The goal is to serve all the requests at a minimum cost. [23] provides a description of the solution methods proposed by the winners of the second and of the third place at the VeRoLog Solver Challenge 2016-2017. Specifically, the approach of the second ranked competitor involves sequences of low level heuristics, each of them having a transition probability to move to another low level heuristic. The third ranked method, instead, decomposes the problem into two independent subproblems, which are solved using a genetic algorithm, a deterministic routing algorithm, and a variable neighborhood descent algorithm to further improve the routes found by the deterministic routing algorithm. Table 1 summarizes the different characteristics for the pickup and delivery of renewable and non-renewable resources literature. A set O j is associated with each node j ∈ N , indicating the set of the daily services or operations required by customer j. For example, in a Health Care context the required operations may include assistance with respect to personal needs such as getting out of bed, bathing, dressing, as well as medical services such as drug administration, blood test and ECG test.
A partial ordering relation < is assumed to be associated with O j for each j, which specifies (partial) precedence relationships for the operations at j. Specifically, if o 1 and o 2 belong to O j for some j, o 1 < o 2 means that the operation o 1 must be accomplished before the operation o 2 is started in order to satisfy the requirements of node j. In a Health Care context, for example, dosing medicine must be performed before medication and blood test must be performed before a meal. A time window [a jo , b jo ] is associated with each customer and operation, indicating the time interval in which the operation o at customer j must start, j ∈ N and o ∈ O j . Hereafter O will denote the set of the available services or operations, i.e. O = j∈N O j .
In addition to the normal services, just requiring the presence of a qualified technician, in our study we consider the situation in which a special device, say s, is needed for some of the operations. Let O s denote the subset of the operations needing such a special device. For example, to perform a blood test a centrifuge for blood stabilization is required, while a ECG machine must be available for the ECG test. An operation in O s can be performed by the assigned technician only if the special device is present when the operation is started, and the device must remain at the customer's home for the entire duration of the operation. W.l.o.g. we shall assume that each customer may require at most one operation in O s .
Let T be the set of the available technicians, and let S t denote the skill or ability of technician t, expressed as the set of the operations that t is habilitated to perform. The skills constrain how the customer requirements can be operated: an operation o required by a customer j, i.e. o ∈ O j , can be operated only by a technician t possessing the corresponding skill, i.e. o ∈ S t . Moreover, t can perform multiple operations at j. Concerning the transportation of the special device, we assume that its transportation is in charge of the technicians in T : different technicians may move s, during their own tours, to and from the location of the customers needing the special service; eventually, s can be transported via locations of customers do not requiring s, if this strategy can be useful for improving the overall service. Notice however that a technician visiting a customer has in any case to perform an operation at the customer site. By extending this strategy, the special device could be transported through sites other than the customer locations, like parkings and break centers, if present in the logistics network.
Given nonnegative technician dependent travelling costs c t ij , for each (i, j) ∈ A and technician t ∈ T , we study the problem of defining a set of tours for the technicians, each one starting and ending at the depot d, in such a way that the set of the operations required by each customer j, i.e. O j , is fulfilled by qualified technicians (i.e. the skill constraints are satisfied), by respecting the time window constraints and the partial ordering relation associated with O j . This problem will be referred to as G-Skill VRP. The precedence constraints among the customer operations are imposed as follows: given the nonnegative travelling times τ ij , for each (i, j) ∈ A, and given the operation duration δ t o , for each operation o and technician t, then if o 1 and o 2 belong to O j , and o 1 < o 2 , the completion time of the technician having in care the operation o 1 at node j must be less than or equal to the starting time of the technician having in care the operation o 2 at node j. Indeed, more general precedence constraints are considered in this paper, that include synchronization among operations as their special case. Such constraints will be better specified next. Constraints on the maximum workday duration of the technicians are considered as well. Precisely, the working time of technician t can not exceed his workday duration D t , t ∈ T . The objective function, to be minimized, is the overall travelling cost.
Next we summarize the introduced notation: set of nodes, N is the set of the n customers and d is the depot node; A set of arcs; O set of the available services or operations; O j subset of the daily services or operations required by customer j ∈ N ; O s subset of the operations needing a special device; T set of the available technicians; S t subset of the operations that technician t ∈ T is habilitated to perform; As outlined before, a specificity of the addressed problem, which makes it particularly challenging but also very complex, is the management of a special device, which is required for some of the operations. The special device is moved by the technicians during their own tours, and therefore non standard synchronization between the movement of the special device and the tours of the technicians must be imposed, giving rise to a peculiar tour for the special device, defined as a composition of some technician tour fragments. This is illustrated in   The three technicians have also in charge the movement of the special device, which is needed to perform operation 1. The route of the special device, resulting from the composition of the routes of the technicians moving it, is outlined in red in Figure 2

The MILP model
Hereafter we shall present a MILP formulation to G-Skill VRP. Firstly we shall define the enlarged network used to better capture and model the peculiar problem characteristics. Then the various groups of constraints characterizing the MILP model, which is defined on the enlarged network, will be presented and commented, together with the objective function. The overall model will be reported in the appendix.

The enlarged network
The problem is formulated on an enlarged network EG = (EV d , EA), which is derived from the original logistic network G = (N d , A) as follows. For each customer j ∈ N and each operation o ∈ O j to be performed at node j, we introduce a customer/operation node v = j, o in the enlarged network. For each customer j ∈ N , we also introduce two fictitious nodes, j, s in and j, s out . The two fictitious operations s in and s out associated with each customer j, which are added to O j , are used in order to establish the entrance to and the leaving from the customer location of the special device, as better specified next, and they are also included in the operation set S t of each technician t, since each technician can transport the special device. The aim is to be able to capture, within the model, the entrance and then the leaving of the special device at a customer location, so allowing its synchronization with the movements of the technicians in charge of its transportation and, eventually, with the technician who needs to use it. Based on these modelling issues, we define the set The latter includes the depot node d and a fictitious operation o d ∈ O d such that o d ∈ S t for all technicians t. The fictitious operation o d is introduced in order to habilitate each technician to visit d.
Concerning the set of arcs in the enlarged network, EA contains the arcs Notice that, according to the stated definitions, a complete subgraph is present in the enlarged network for each customer j, linking together all the operations required at j, except for the fictitious operations s in and s out . Hereafter the arcs involving the nodes s in and s out will be also referred to as special arcs. The special arcs (u, v), with v = j, s in , are used only by the (unique) technician transporting the special device to customer j, if this movement is planned. In other words the special device, if required at j or moving along the network through node j, arrives at j exclusively via a special arc entering j, s in , and leaves j via a special arc outgoing j, s out . The technicians do not transporting the special device, instead, do not use the special arcs. In other words, fictitious nodes of type j, s in and j, s out are visited only by a technician transporting the special device.
The enlarged network is represented in Figure 3 (a). There, the little brown nodes represent the fictitious nodes j, s in and j, s out at each customer j. The tours of the technicians for the example depicted in Figure 2 (a) and the tour of the special device, depicted in Figure 2 (b), are represented, in terms of the enlarged network, in Figure 3 (b). The figure shows that the tour of the special device, depicted in red, is defined as the composition of the special arcs visited by the technicians, plus suitable movements inside the customer nodes.  In the enlarged network, the travelling time associated with an arc (u, v), u = i, o 1 , v = j, o 2 , with (i, j) ∈ A (and so i = j), is the travelling time along (i, j) in the original logistic network, i.e. τ ij . Notice that this holds also for the special arcs. On the other hand, the travelling time associated with any arc (u, v) such that u = j, o 1 and v = j, o 2 for some customer j, i.e. with the arcs linking operation nodes related to a given customer, is zero. Moreover, each node v = j, o in the enlarged network inherits the service times associated with the operation o corresponding to the node, which depend on technician t, and were previously defined as δ t o , whereas the service time for the special node d, o d and the fictitious nodes j, s in and j, s out is zero. Concerning the travelling costs in the enlarged network, which are dependent on the technicians, they are set as follows: . This holds true also for the special arcs. Otherwise, i.e. (u, v) is such that u = j, o 1 and v = j, o 2 for some customer j, then c t uv = 0 for all t.
Next we introduce the decision variables and the groups of constraints.

Groups of constraints and objective function
Regarding the formulation, recall that each customer j ∈ N requires a set O j ⊆ O of operations to be performed, where O = ∪ j∈N O j will be used to denote the set of all the required operations. Recall also that an operation o ∈ O j can be performed only by a qualified technician, i.e. by a technician t ∈ T such that o ∈ S t , where S t ⊆ O is the set of operations t can perform. Furthermore, at most one operation o ∈ O j may require the special device, i.e.
Always regarding the formulation, given an arc as the subset of the technicians who can perform both o 1 and o 2 , and so can move along (u, v). Also define EA t as the subset of arcs where technician t ∈ T can move along, that is the subset of arcs (u, v) ∈ EA, with u = i, o 1 and v = j, o 2 , such that both o 1 and o 2 belong to S t .
Assignment of operations to technicians. We start by introducing the group of binary variables modelling the assignment of the operations to the technicians. For each node v = j, o ∈ EV and t ∈ T such that o ∈ O j ∩ S t , let: Any operation of a customer, except the fictitious ones, must be performed by exactly one technician; therefore the assignment variables must satisfy: On the other hand, the special nodes (j, s in ) and (j, s out ) are visited by a technician only if the special device is required at j, or if j has just to be traversed to move the special device along the network. Therefore: Routing of the technicians. Now we introduce the group of the binary variables modelling the routing of the technicians. For each (u, v) ∈ EA and t ∈ T uv : The routing variables must satisfy constraints (4.3) -(4.7). With a little abuse of the notation, hereafter the node d, o d will be simply denoted as d.
Constraints (4.3) guarantee that, if an operation o in O j is assigned to technician t, then t must visit node v = j, o . Notice that this includes also the special nodes j, s in and j, s out . Constraints (4.4) state that, if technician t enters a node v, then t must leave v. Notice that a technician may perform more operations at the same customer j, by moving among nodes of form (j, o), for different operations o ∈ O j , with travelling time zero, according to the definition of the enlarged network. Moreover, a technician may perform more operations at the same customer j by moving to different customers and then coming back to j. Therefore, multiple visits of a customer by the same technician are allowed.
Constraints (4.5) and (4.6) ensure that every technician leaving the depot to perform the assigned operations, then enters the depot after operating. Moreover, to model the case of technicians who are not assigned any operations, and therefore must remain at the depot, we assume that a fictitious arc (d, d), connecting the depot to itself, and having travelling time zero, belongs to the set EA t for each t. In constraints (4.5)-(4.6), the routing variable x t dd will be set to 1 in case of inactivity of technician t, and to 0 otherwise. Constraints (4.7) define the variable domain.
Technician temporal constraints. In order to prevent subtours, that is tours which are disconnected from the depot, the proposed model uses temporal flow variables and standard flow conservation constraints, plus linking constraints to relate the flow variables to the routing variables. The temporal flow variables are used also to impose precedence and synchronization constraints among the customer operations, as better clarified next. Let us discuss the temporal flow variables and the related constraints in more detail.
Let y t uv , with (u, v) ∈ EA and t ∈ T uv , be a flow variable that specifies the cumulative time spent by technician t after traversing the arc (u, v), if (u, v) belongs to the tour of technician t. This cumulative time is defined as the sum of the travelling times associated with the arcs of the subpath from the depot to v, plus the durations of the operations performed by t at the nodes of this subpath, except the one operated at node v. Therefore, y t uv indicates the arrival The temporal flow variables are defined in terms of the technician assignment variables z t v and of the technician routing variables x t uv : Constraints (4.8) take into account that at the depot node d the duration of an operation is zero and therefore the cumulative time of the first visited customer, if any, accounts only for the duration of the travelling from the depot node to the customer node. The flow conservation constraints (4.9) define the total time accumulated by a technician t after the service at a customer node v as the cumulative time from the depot to v plus the duration of the operation at v plus the travelling time along the arc of the subpath leaving node v, if v belongs to the tour of t. Constraints (4.10) and (4.11) state that the cumulative time spent by technician t can be positive only on the arcs belonging to the tour of technician t. Constraints (4.10) play also another role into the model. In fact, they state that the total time of technician t in any point of his tour can not exceed D t , i.e. the workday duration associated with t.
Notice that, in order to allow some waiting time along the technicians route, equalities (4.9) could be replaced by the inequalities Precedences and synchronization. By using the temporal flow variables, the precedence constraints among the customer operations can be modelled as follows: Constraints (4.12) state that, for each customer j, if the operation o 1 must preceed o 2 , and if technicians t and r have been selected, respectively, to operate o 1 and o 2 at node j, then the cumulative time of t at v 1 plus the time spent by t for operating o 1 must be less than or equal to the cumulative time of r at node v 2 . Notice that r = t is allowed, so modelling the case where the precedence relationship involves operations of j assigned to the same technician.  Time windows. Concerning the time windows constraints, stating that the operation o at customer j must be started in the time interval [a jo , b jo ], for j ∈ N and o ∈ O j (o = s in , s out ), they can be expressed as: Routing of the special device. Now, let us model the handling of the special device within the enlarged network. As previously outlined, the modelling idea is that the special device can arrive at customer j only through the node j, s in , thanks to a technician. Then, if the special device is required at j, it is transported directly to the requiring node j, o , with o ∈ O j ∩ O s , and, from there, after the completion of operation o, to the outgoing node j, s out , to leave customer j. Notice that the technician operating o at j may be different than the technician who have transported the special device at j. Also the technician who is in charge of moving the special device out of j may be different. Otherwise, i.e. the special device is not required at j, i.e. j is just used to move the device elsewhere, then the device moves directly from j, s in to j, s out , to leave customer j. Also in this scenario the arrival and the leaving of the special device to/from j are in charge of a technician (possibly two different technicians). In order to model the movement of the special device, let us denote by EA s the set of the arcs along which the special device can move. This set contains: the arcs of form ( i, s out , j, s in ) for each pair of different customers i and j such that an arc (i, j) exists in the original logistic network, modelling the movement of the special device from customer i to customer j; the arcs of form ( j, s in , j, s out ), modelling the case where the special device is not required at j, but just traverses j; the arcs of form ( j, s in , j, o ) and ( j, o , j, s out ), modelling the case where the special device is required at j to perform the special operation o, with o ∈ O j ∩ O s ; and finally the arcs of form (d, j, s in ) and ( j, s out , d), modeling the leaving of the special device from the depot, and its coming back to the depot at the end of the working day. By extending the previous definition, all these arcs will be also referred to as special arcs. Then, let us introduce the following routing variables: Constraints (4.15)-(4.24) are the routing constraints for the special device. In particular, constraints (4.15)-(4.17) guarantee the presence of the special device at each customer requiring it, whereas constraints (4.18)-(4.20) address the possible movement of the special device through the location of a customer does not requiring it. Constraints (4.21) and (4.24) consider the movement of the special device from and to the depot node, while (4.22) and (4.23) are the classical flow conservation constraints for the special device routing.
Constraints (4.25), (4.26) and (4.27), instead, ensure that the special device can move between different customers, and from/to a customer and the depot, if and only if a technician moves along the corresponding special arcs during his tour, by transporting the device. The tour of the special device, at the original logistic network level, is thus the composition of some technician tour fragments, related to those technicians having in charge the transportation of the special device. The enlarged network, instead, allows one a deeper view of the technicians and of the special device tours, by showing their behavior inside the customer nodes.
Precedence and synchronization for the special device. Here we introduce the temporal variables and constraints related to the special device. Hereafter D s will denote the daily availability of the special device.
The temporal variables for the special device, y s uv , are defined only for the special arcs linking different customers in the tour of the special device, and for the special arcs of this tour linking customers to the depot. They are defined as a function of the temporal variables of the technicians according to the following observations, which are based on the previously stated assumption that the special device is moved along the network by means of technicians: (a) the arrival time of the special device at customer j coincides with the arrival time at node j, s in of the (unique) technician visiting j, s in (recall, in fact, that the special device enters customer j only via the node j, s in ); (b) the leaving time of the special device out of customer j coincides with the arrival time at node j, s out of the (unique) technician visiting j, s out (recall, in fact, that the special device leaves customer j only via the node j, s out ). These modelling tools allow one to synchronize the arrival and the leaving time of the special device at customer j, and also to synchronize such times with the time of the technician in charge to use the special device at j, if j requires the special device. Formally: Constraints (4.29) state that the arrival time of the special device at customer j coincides with the arrival time of the (unique) technician transporting the device to j, i.e. visiting j, s in . Constraints (4.30) impose that the special device must arrive at customer j, requiring it, before the arrival of the technician who will use the special device at j, while (4.31) state that the use of the special device at j must be ended before the time when a technician will move the device out of j. Concerning the customers do not requiring the special device, i.e. those customer locations who may be used just to transport the device along the network, (4.32) state that the arrival time of the technician transporting the device at a customer must precede the leaving time of the one who will move out the device. Finally, constraints (4.33) guarantee that the total time of the special device tour can not exceed D s and that the temporal variables for the special device may be positive only along the arcs used to transport it.
Objective function. The objective function, to be minimized, defines the overall routing cost: The overall model is reported in the appendix.

Valid inequalities
In order to enhance the LP relaxation of the proposed formulation, the following valid inequalities, defined on the enlarged network, have been proposed. These inequalities, although quite simple, proved to be very effective from a computational perspective, as shown in the section devoted to the computational results.
Depot cut constraints. These inequalities state that, if a visit is assigned to a technician t, then t must leave the depot node. Equivalently, the routing variable x t dd , modelling the inactivity of t through the fictitious link (d, d), must be set to 0.
See Figure 4 for a graphical explanation.

(j,o)
Customer j d Depot d Technician t Figure 4: Depot cut constraints 2-3 cycle cut constraints. The 2 cycle cut constraints state that, if two operations of a certain customer j, other than the fictitious ones s in and s out , are assigned to a technician t, and thus the corresponding nodes in the enlarged network belong to his tour, then no cycle can exist in the tour of t comprising such nodes. Analogous inequalities can be stated by considering the case of three operations of customer j assigned to the same technician, leading to 3 cycle cut constraints.
These valid inequalities are illustrated in Figure 5.

2-3 aggregated cycle cut constraints.
As only one technician operates each required operation, the 2-3 cycle cut constraints can be aggregated by technicians as follows: At most one arc At most two arcs Aggregated inter-customers cycle cut constraints. They state that, given two customers i and j, if a technician t has to perform both an operation of i and an operation of j, then t can move either from i to j or from j to i. That is, cycles involving the corresponding nodes in the enlarged network are forbidden. Clearly, such valid inequalities can be generalized by considering more than two customers at a time.

Computational experiments
In this section we report computational experiments performed to access the quality of the bounds obtained for G-Skill VRP. All the computational tests were performed using a processor Intel(R) Core(TM) i7-4750HQ CPU @ 2.00 GHz with 8GB of RAM and using the software Xpress 8.5 (Xpress Release 2018 with Xpress-Optimizer 33.01.05 and Xpress-Mosel 4.8.4) [18].
In Section 6.1 we provide details on the generation of the instances. In Section 6.2 we describe the procedures we applied to the models to obtain bounds to the problem. In Section 6.3 we report the computational results obtained for the generated instances, and discuss the main achievements that can be derived based on these results.

Instance description
We experimented two different sets of instances. The first set, comprising 60 instances of medium to large size, has been generated starting from the Euclidean-distance data sets tc that were proposed and used by Gouveia in [22] for the first time, and utilized since then in several computational testings. The name tc of the data sets refers to a central position of the depot in a grid. The second set, instead, is composed of 9 very large instances, and it has been generated starting from the TRSP instances used in Pillac et al.'s experimentation [29], which in turn are based on the Solomon benchmark [34]. We want to emphasize that the Technician Routing and Scheduling Problem, investigated in [29], is different than G-Skill VRP. Therefore, Pillac et al.'s instances have been adapted to face with the specific characteristics of G-Skill VRP.

tc-based instances
We generated a total of 60 instances with |N | = 20 customers, |O| = 6 available operations and |T | = 3 technicians. Two of the technicians are skilled for 4 operations, while the third one is skilled for all the operations. Each technician is able to perform the special operation, i.e. operation 1, that requires the special device, and the working day duration is 480 minutes (8 hours) for all the technicians. Each customer requires 1 or 3 operations and the 20% of the customers require that their operations start in a given time window, either morning (first half of the working day duration) or afternoon (second half of the working day duration). Notice that, since in several application contexts the service delivered at the customer locations may be variable or uncertain, as remarked in [27], it may be reasonable to consider large time windows, such as morning and afternoon, rather than hard time windows, since these could be rarely satisfied in practice.
The 60 instances are divided into four sets of 15 instances each. Sets of instances have different percentages of the number of operations required per customer. Specifically, we have the following sets: (i) set 70-30 in which 70% of the customers (14 customers) require one operation and 30% of the customers (6 customers) require 3 operations, summing up to 32 operations, (ii) set 60-40 in which 60% of the customers require one operation and 40% of the customers require 3 operations (for a total of 36 operations), (iii) set 50-50 in which 50% of the customers require one operation and 50% of the customers require 3 operations (for a total of 40 operations), (iv) set 40-60 in which 40% of the customers require one operation and 60% of the customers require 3 operations (for a total of 44 operations).
We used five Euclidean data sets from [22], namely tc1, tc2, tc3, tc4 and tc5. For each one of the five Euclidean data sets, three different seeds (identified by the numbers 0, 1 and 2) were used to randomly assign operations to customers and technicians. The name of each generated instance reports the name of the original data set used (tc1, tc2, tc3, tc4 or tc5), the set it belongs to (70-30, 60-40, 50-50 or 40-60), and the random seed (0, 1 or 2). As an example, instance tc1-70-30-0 refers to the instance generated from the original instance tc1, belonging to the set 70-30, using the first seed (zero).
The service time δ t o of each one of the |O| = 6 available operations ranges from 10 to 30 minutes. The operation that requires the use of the special device, i.e. operation 1, has a service time of 10 minutes.
The travel time τ ij between customer i and customer j is set to 0.3c ij where c ij is the Euclidean distance value from the data sets tc for nodes i and j. Also the traveling costs c t ij , which depend on the technicians, are defined using the values c ij of the Euclidean data sets. Specifically, two technicians have the same travelling cost, which is equal to c ij , for each pair of customers i and j, whereas the third technician, who is the most skilled among the technicians, has a travelling cost that is two times the cost of the other ones, i.e. 2c ij , for each pair of customers i and j.
As for the precedence constraints, we consider two cases: (i) one type of precedence con-straints in which operation 2 must be performed before operation 3 at each customer requiring them, i.e. 2 < 3, and (ii) two types of precedence constraints in which operation 5 must always be performed before operation 6 in addition to the precedence constraint between operations 2 and 3, i.e. 2 < 3 and 5 < 6.

T RSP -based instances
We selected 9 TRSP instances from [29] which involve 100 customers and 25 technicians. In such instances, each technician has a skill and a set of tools and spare parts he can manage. On the other hand, each customer requires a skill and a set of tools and spare parts.
In order to adapt this different scenario to the G-Skill VRP context, we disregarded spare parts and interpreted tools as operations, i.e. each of the tools required by a customer in the original instance becomes an operation in our instance. Consequently, if a customer i requires a skill k and a tool o in the original TRSP instance, then i becomes a customer requiring skill k and operation o in the corresponding G-Skill VRP instance, and he can be served only by a technician having skill k and possessing o in his tool set. As an example, if i requires skill k and needs tools o 1 and o 2 in the original TRSP instance, in the corresponding G-Skill VRP instance i becomes a customer requiring two operations, i.e. o 1 and o 2 , both with skill k. The total number of operations in the G-Skill VRP instance is thus the total number of tool requirements in the original TRSP instance. Notice that the concept of skill has been generalized with respect to the one presented in Section 3, and the models have been adapted accordingly. In Section 3, skills correspond to operations, and a technician t can perform an operation o at a certain customer if and only if o belongs to set of operations t can perform. In the TRSP instances, instead, skills do not correspond to operations. A technician t can operate operation o requiring skill k at a customer if and only if t has skill k (according to the definition of skill in [29]), and its set of operations (i.e. his set of tools in the original TRSP instance) includes o.
One tool, i.e. tool 1, has then been selected to define the special operation, that is the operation which requires the special device. Recall that no special device is managed in [29].
The size of the T RSP -based instances is summarized in Table 2. The name of each instance is exactly the one as the original instance and it consists of three fields that identify respectively how requests at customers are generated, the length of the planning horizon (1 day in our case), and the instance number. In regards to the first field, requests can be randomly generated (R), organized in clusters (C), or generated combining both (RC). In addition to the instance name, columns in Table 2 report respectively the number of customers (|N |), the total number of operations required by overall customers, and the number of technicians (|T |). Notice that, since we disregarded spare parts, we removed all the customers requiring only spare parts. Similarly, technicians only managing spare parts have been disregarded. Despite of some customer and technician removals, the obtained T RSP -based instances are really very large, with a total number of operations ranging from 86 to 120.
Service times for the operations and workday durations for the technicians, which are not provided in the original instances, are defined as for the tc-based instances. In particular, the working day duration is 480 minutes (8 hours) for all the technicians and service times range in {10, 15, 20}. Technician costs depend on distance and skill. Specifically, we assume that for each technician who is able to perform no more than three operations, the cost of arc (i, j) is exactly the distance between customer i and customer j, whereas for each technician who is able to perform more than three operations the cost of an arc is twice the distance. Concerning time windows, which are not provided in the original instances, we still consider large time windows, i.e. a morning or an afternoon time window. Differently from what happens in the tc-based instances, where only the 20% of the customers has an associated time window, in the T RSP -instance we associated with each customer a morning or an afternoon time window, thus making the instances even more challenging. Finally, no precedence constraints are given.

Solving procedures
In this section, we describe how different variants of the MILP model (4.1)-(4.35) and several procedures have been used to provide lower and upper bounds to the G-Skill VRP. Specifically, we tested four models and three procedures as described in the following.
The first model, referred to as the plain model (p for short) is just the complete model whose feasible set is described by constraints (4.1)-(4.34). The second model, referred to as model with cuts (c for short) is the plain model enhanced with several valid inequalities. Namely, the depot cuts (5.1), the 2-3 aggregated cycle cuts (5.4)-(5.5) and the aggregated inter-customers cycle cuts (5.6). The third model, referred to as the relax model (r for short), is a relaxation of the plain model in which the special device is neglected. The relax model is thus defined by a subset of the constraints that characterize the plain model, namely constraints (4.1)-(4.14). Finally, the fourth model consists in the relax model equipped with the same valid inequalities used in model c, namely (5.1), (5.4)-(5.5), (5.6). This model is referred to as relax model with cuts (x for short). For each of these models, the solver is used to obtain the solution value of the linear programming relaxation and the values of the best lower and upper bounds returned by the Branch and Bound procedure within an a priori fixed maximum amount of computational time. Specifically, for each model m with m in {p, c, r, x}, LP m represents the linear programming relaxation solution value, LB m the best lower bound value and U B m represents the best upper bound solution value. At this regard, notice that U B r and U B x are not necessarily upper bounds to G-Skill VRP. The name of the four alternative models and the constraints describing them are briefly recalled in Table 3. Table 3: Model names and constraints describing them.
Constraints No Cuts Cuts (5.1),(5.4),(5.5),(5.6) (4.1)-(4.34) p c (4.1)-(4.14) r x Preliminary computational results have shown that the relax model is more efficient than the plain model. A possible explanation of why this happens is that by eliminating the complicating constraints related to the special device every single subproblem in the Branch and Bound tree becomes easier to solve. As a result, more nodes of the enumeration tree are explored and this allows to find better lower bounds. Motivated by these preliminary tests, we defined two procedures aiming at improving the lower bound values. In the first procedure, referred to as p 1 , the lower bound value obtained by model x is used as a lower bound value to model c, thus implementing a warm start for model c. As before, LB p 1 and U B p 1 represent, respectively, the best lower and upper bound values returned by procedure p 1 .
In the second procedure, referred to as p 2 , as in the first procedure the lower bound value obtained by model x is used as a lower bound constraint, and in addition we dynamically add violated constraints from the sets (4.15)-(4.34) to the relax model. As before, LB p 2 and U B p 2 represent, respectively, the best lower and upper bound values returned by procedure p 2 . Notice that U B p 2 is not necessarily an upper bound to G-Skill VRP.
While the two procedures above described aim at improving the lower bound value, in some cases returning an upper bound as a collateral effect, the third procedure is a matheuristic, based on the LP solution of model c, that assigns all the nodes requiring the special operation o s to one technician at a time, in an exhaustive way. Its aim is to rapidly estimate the quality of the lower bounding techniques proposed in this paper, thanks to its computational efficiency. Specifically, the LP solution of model c is used to rank the technicians according to the (possibly fractional) number of nodes of type v = j, o s assigned to them in the LP solution, by considering a nonincreasing order. The first technician according to this ordering, say t s , is selected, and all the nodes of type v = j, o s are assigned to him. Therefore, the corresponding assignment variables z ts j,os are set to one within model c, and the plain model so simplified is used to generate a feasible solution with a time limit of 1800 seconds. Then, the next technician according to the considered ordering is selected, all the nodes of type v = j, o s are assigned to him, and again the simplified model c is solved to get a feasible solution, by iterating until all the technicians have been considered, and halving the time limit for the simplified model c solution at each iteration. The resulting matheuristic, referred to as h for short, returns the best feasible solution found among the performed iterations. Its pseudocode is reported in Algorithm 1. By extending the previous notation, U B h will denote the value of the solution returned by h, since it provides an upper bound to G-Skill VRP. Notice that the best lower bound value obtained by solving the simplified versions of model c along the heuristic iterations is not necessarily a lower bound to G-Skill VRP.
Set z ts v = 1 in model c; 16: end for In the next section, detailed results relative to the four models, namely, p, c, r and x, and to the three procedures, namely p 1 , p 2 and h, are presented and discussed separately for the two sets of G-Skill VRP instances.

Computational results for tc-based instances
We start the section with a set of tables that compare pairs of models and/or procedures. Specifically, in Table 4 we compare models p and c, while in Table 5 we compare models r and x. These two tables allow to evaluate the impact of adding valid inequalities on the quality of the lower bounds obtained. Summary information on optimality gaps and computational times for models p, c, r and x is reported in Table 6. Then, after a discussion about the two procedures, p 1 and p 2 , which exploit the relax model with the aim of providing enhanced lower bounds, a lower bound comparison is proposed in Table 7 by reporting for each model, and separately for each set of instances, the average relative gap of the corresponding lower bound with respect to the best lower bound achieved. Information on the best lower and upper bounds obtained via the proposed models and procedures is finally summarized in Table 8, also considering the matheuristic approach h. Additional summary information is given also via a graphical representation.
Computational results obtained when two types of precedence constraints are considered exhibit the same trend of the results obtained with one type of precedence constraint and are thus not reported. Tables 4 and 5 report the same type of information. Specifically, the first column identifies the instance. Then, for each considered model or procedure, identified by a character in the subscript, columns LP display the value of the linear programming relaxation, columns LB display the best lower bound value obtained within the time limit imposed, columns U B display the best upper bound value obtained within the time limit imposed and columns gap display the corresponding optimality gap between the lower and the upper bound found, i.e., gap m = 100 × (U B m − LB m )/LB m where the character m identifies the model or the procedure used. Additionally, in some tables, columns "time" report the computational time in seconds used to obtain the displayed values. A time limit of 10800 seconds was imposed. 1 This is the time used when the column "time" is omitted.
Comparison of the plain model with the model with cuts. Table 4 reports results that allow the comparison of the plain model with the model with cuts. The best results obtained within a time limit of three hours of computations are reported. As indicated before, for each of the two models, either p or c, the gap measures the relative distance between lower and upper bounds obtained with that specific model.
The direct comparison between columns LP p and LP c as well as between columns LB p and LB c highlights that the valid inequalities are largely effective in improving the lower bounds. Furthermore, valid inequalities allow to reduce the gap remarkably and they usually allow to increase the number of instances for which a feasible solution is found. Specifically, from Table 4 we get that, on the set 70-30, the average gap computed on all the instances for which a feasible solution is found decreases from 62.6% obtained with model p to 26.7% obtained with model c.
In addition, the introduction of cuts allows to close the optimality gap for 4 instances over 15 and to determine a feasible solution for every instance. On the set 60-40, the average gap decreases from 119.7% (model p) to 25.8% (model c). In addition, model c allows to find a solution also on the 5 instances on which model p fails, thus determining a solution for all the instances in the set. The gaps decrease remarkably also on instances 50-50 (from 256.2% to 19.6%) and on instances 40-60 (from 255.6% to 66.2). However, on these two sets the introduction of cuts might on some instances prevent the solver from finding a feasible solution. This is probably due to the fact that on the most difficult instances in the test bed the introduction of all the cuts from scratch slows down the solver.
Comparison of the relax model with the relax model with cuts. Table 5 reports results that allow the comparison of the relax model with the relax model with cuts in terms of linear programming relaxation value, best lower bound value and best solution value obtained within a time limit of three hours of computation. Both the models neglect the constraints involving the special device and thus, as observed, their upper bound does not necessarily correspond to a feasible solution to the original problem.
Noticeably, we observe that both the relax model with cuts and the relax model are able to obtain a solution for all the instances. The lower and upper bound values obtained with the relax model with cuts are significantly better than those obtained with the relax model. Further, Table 5 shows that the relax model with cuts is able to close the optimality gap for many instances (36 out of 60 for the one precedence constraint, and indeed the same 36 instances out of 60 for the two precedence constraints). Contrarily, the relax model is able to close the optimality gap only for three instances over 60 (all of them belonging to the 70-30 set). Thus, embedding the model with cuts, also in this case allows to speed up the solver. Overall, the introduction of cuts makes it possible to significantly reduce the optimality gap. On the 70-30, the average gap reduces from 17.2% when model r is used to 3.6% when model x is used. The same figures show that the gap decreases from 26.1% to 2.0% on the 60-40, from 51.3% to 3.1% on the 50-50, and from 111.7% to 5.8% on the 40-60.
When comparing the computational results of Table 5 with the corresponding results displayed in Table 4, we conclude that the presence of the special device makes the problem much harder to solve. In particular, the average computational time of model x halved in comparison with that required by model c. Table 6 summarizes the average results for each of the four sets of instances, by reporting the average optimality gap and the average computational time obtained with models p, c, r and x. The average time of model p is not reported since the time limit was always reached.
Comparison of procedures p 1 and p 2 . We also compared the two procedures using the value of the relax model with cuts as a lower bound cut and then adding the special device constraints, either statically (in procedure p 1 ) or dynamically (in procedure p 2 ). The comparisons are done in terms of the best lower bound value and the best upper bound value obtained within a time limit of three hours of computation. Note that the upper bound values obtained with procedure p 2 are not necessarily upper bound values to G-Skill VRP problem since the maximum computational time can be reached before all the constraints have been added. As a case study, we considered only the instances for which the optimality gap has been closed by the relax model with cuts. As already said, this happens for 36 out of 60 instances.
We report that, noticeably, interesting results were achieved by procedure p 1 , which was able to certify the optimality of the bounds provided by model x for 7 out of 36 instances. In addition, it allowed to reduce the gap obtained by model x on three instances, namely, tc2-70-30-2, tc4-70-30-1, and tc2-50-50-2.
Comparison of the best lower and upper bounds. To better compare the lower bounds computed by means of models p, c, r and x, for each model, and separately for each set of instances, we determined the best lower bound achieved by all the suggested models and procedures, and calculated the average relative gap, over the 15 instances of each set, of the obtained lower bounds with respect to the best lower bound. This is reported in Table 7. Specifically, the quality of both the lower bound given by the linear programming relaxation, i.e. LP, and the best lower bound obtained by each model, i.e. LB, within the given maximum computational time are reported.
Finally, Table 8 reports results that allow the comparison of the best values determined by the proposed models and procedures within the time limit of three hours of computation, also considering the upper bounds obtained by the matheuristic h. Precisely, column bestLB reports the best lower bound achieved by all the proposed models and procedures, while column bestU B is the value of the best feasible solution found, also considering the heuristic h. The specific model or procedure allowing to obtain the best lower bound (bestLBm) and the best upper bound (bestU Bm) is also reported for each instance. To get more information on the behavior of the heuristic, columns U B h report the values of the feasible solutions obtained by h, while columns time h indicate the required computational times. The results displayed in Table 8 show that the proposed models and procedures allow to close the gap for 6 instances in the first two sets, and for one instance in the third set. Furthermore, the gaps are often quite small for the first three sets of instances, whereas they may be quite large for the last set, comprising the more difficult instances of this first test bed.
We conclude the discussion on the tc-based instances with a global overview of the results obtained. Specifically, Figures 6,7,8,and 9, respectively for the instances 70-30, 60-40, 50-50 and 40-60, compare the lower bounds obtained by the four alternative models, namely p, c, r, and x, and the best upper bound, i.e. bestU B, against the best lower bound obtained, i.e. bestLB. That is, the percentage gap between the four lower bounds and bestLB is reported, as well as the percentage gap between bestU B and bestLB. The figures support a series of main observations: the best lower bound is almost always given by the model which disregards the special device constraints, i.e. x, once enhanced with cuts; in several cases, this model is also able to close the optimality gap; the proposed valid inequalities, defined on the enlarged network, are quite effective in improving the quality of the lower bounds; the difficulty of obtaining good bounds seems to increase with the total number of operations considered, thus making the 40-60 instances the most challenging in this data set; on some instances to determine a good upper bound seems to be critical.
Some of the proposed lower bounding techniques appear thus to be a valuable tool to compute good lower bounds to G-Skill VRP. In several cases, in fact, especially when the total number of operations is not very large, the associated optimality gaps, evaluated also via a simple matheuristic approach, are quite small, showing the efficacy of the presented models and techniques.

Computational results for T RSP -based instances
We tested models p, c and x, as well as procedures p 1 and p 2 , to provide lower bounds to the T RSP -based instances. Due to the complexity and the very large size of these instances, no one of such models returned an upper bound as a collateral effect. Moreover, the matheuristic      h failed in computing feasible solutions. We thus generalized h, in an attempt to determine upper bounds also for this second set of instances.
The idea is to handle the huge size of the T RSP -based instances by decomposition, i.e. partitioning the set of the customers of each instance into three subsets: i) a subset named special, which includes all of the customers requiring the special device, thus extending the philosophy of matheuristic h, where all the special operations are dealt together; ii) a subset named morning, which is composed of customers with a morning time window; and iii) a subset called afternoon, comprising customers with an afternoon time window. A generalized bin packing problem is solved for each subset of customers, in order to determine a subset of technicians able to serve them. Three smaller and independent G-Skill VRP subproblems are thus derived, for the special, the morning and the afternoon customer subsets, respectively. The three subproblems are then solved in cascade by model c, along the lines of matheuristic h, by setting a time limit of 2000 seconds for each subproblem. A solution to the overall T RSP -based instance is then obtained by combining the solutions found for the three subproblems.  Figure 10 shows, for each instance, the quality (relative gap) of the lower bounds provided by models p and c with respect to the best lower bound that for this set of instances is always given by procedure p 1 . On this set of instances, procedure p 2 is in fact dominated by procedure p 1 and thus the relative results are not reported. Similarly, we do not report the results of model x, since they are enhanced by p 1 . Figure 10, for each model, reports both the quality of the lower bound given by the linear programming relaxation (LP) and the quality of the best lower bound returned when the time limit is reached (LB). The time limit has been fixed to 3 hours; further computational results obtained with a time limit of six hours demonstrated that there is no substantial improvement when the time limit is extended. The results shown allow to make some interesting observations. First, as for the tc-based instances, the introduction of cuts (model c vs model p) drastically increases the quality of the lower bounds. Specifically, comparing LB p and LB c , we observe that the introduction of valid inequalities makes it possible to reduce the gap of 27.3% on average. Besides, their introduction seems to be even more crucial than for the tc-based instances. In fact, if we consider model p, we observe that there are instances for which the lower bound obtained after three hours of computation is the same as the initial lower bound given by linear programming (see instances C102 and R103 when model p is used), and in general, the improvement between LP p and LB p is still limited. Second, the difficulty in solving the instances does not seem to depend on how the requests are generated (R, C or RC), i.e. quite similar trends are observed on each instance of this set. Unfortunately, the quality of the solutions provided by the matheuristic we proposed is not acceptable. The matheuristic fails in one case over 9 instances and for the remaining ones the average relative gap between the upper bound and the best lower bound is about 164 %. Only a few subproblems (special, morning, afternoon) are solved to optimality and the average computational time on the 8 instances is 4316.35 seconds. This result comes at no surprise since several attempts made to extend the matheuristic used for the tc-based instances revealed that the T RSP -instances are really challenging to solve. Thus, the results obtained stimulate future research in the direction of improving the upper bounds.

Conclusions
We have studied a generalization of the Skill VRP problem, named G-Skill VRP, where multiple visits per customer are allowed in addition to consider precedence constraints among the operations, time windows at the customers, and maximum workday duration constraints. Furthermore, the handling of a special device, which is required to perform some special operations, and which must be present at the customer location together with the technician in charge to use it, constitutes an original and relevant contribution to the literature on VRPs. The resulting Skill VRP generalization appears thus to be very relevant from an operational perspective, and it may constitute a basic block for interesting multi-day planning extensions.
We have proposed and empirically investigated some models and procedures to G-Skill VRP, with the aim of determining good bounds to the problem in an efficient way. The results of a wide computational experimentation show that some of the proposed techniques appear to be a valuable tool to compute good lower bounds to G-Skill VRP. In several cases, in fact, especially when the total number of operations is not so large, the associated optimality gaps, evaluated also via a simple matheuristic approach, are indeed rather small.
We plan to extend the study of G-Skill VRP by proposing more sophisticated heuristic approaches to the problem, to be evaluated also thanks to the presented models and techniques.
Precedence and synchronization constraints.